- return distinct entities in their initial order (rather than sorted), with the `"keepOrder"` option.
- consider all `Nan` values as the same one, with the `"uniqueNan"` option.
* `ellipj()` function introduced, to compute the `sn`, `cn`, `dn`, `ns`, `nc` and `nd` Jacobi elliptic functions.
+* `riccati()` is upgraded:
+ - `riccati(H)` and `riccati(H,E)` syntaxes added, to describe the Riccati equation through its Hamiltonian H or (E,H) pencil.
+ - The residual is returned as new 3rd output argument.
+ - When no solution is found, `X=[]` | `X1=[]` is returned instead of yielding an error.
Help pages:
* `setPreferencesValue` will be removed from Scilab 6.1.x. Please use `xmlSetValues` instead.
* `%sn()` is obsolete. Please use `ellipj()` instead.
* `sysdiag()` is obsolete. Please use `blockdiag()` instead.
+* `ric_desc` is obsolete and will be removed from Scilab 6.1.x. Please use `riccati` instead.
Removed Functions
* [#14604](http://bugzilla.scilab.org/show_bug.cgi?id=14604): `emptystr()` is 40x slower with 6.0.0 wrt 5.5.2
* [#14605](http://bugzilla.scilab.org/show_bug.cgi?id=14605): fixed - `bench_run` was too strict about the specification of tests names.
* [#14606](http://bugzilla.scilab.org/show_bug.cgi?id=14606): Memory used by variables returned by `[names,mem]=who()` was always zero.
+* [#14610](http://bugzilla.scilab.org/show_bug.cgi?id=14610): `x = ric_desc(H,E)` always yielded an error. [x1,x2,residual] = ric_desc(..) returned a wrong `residual` value.
* [#14642](http://bugzilla.scilab.org/show_bug.cgi?id=14642): No more "\r" carriage return with printf.
* [#14741](http://bugzilla.scilab.org/show_bug.cgi?id=14741): The syntax `[m,e]=log2(x)` was not documented. As public function `frexp()` was in duplicate with `[m,e]=log2(x)`.
* [#14746](http://bugzilla.scilab.org/show_bug.cgi?id=14746): Tiny numbers were sometimes displayed as 0.
* [#15269](http://bugzilla.scilab.org/show_bug.cgi?id=15269): `xgetech` was poor and stiff compared to any combination of `gca()` properties `.axes_bounds`, `.data_bounds`, `.log_flags`, and `.margins`. It is removed.
* [#15271](http://bugzilla.scilab.org/show_bug.cgi?id=15271): `bitget` needed to be upgraded.
* [#15321](http://bugzilla.scilab.org/show_bug.cgi?id=15321): `lu()` was leaking memory.
+* [#15350](http://bugzilla.scilab.org/show_bug.cgi?id=15350): `ric_desc()` should be merged into `riccati()`.
* [#15368](http://bugzilla.scilab.org/show_bug.cgi?id=15368): `freson()` silently returned frequencies not corresponding to a maximum, or returned [] instead of some still computable maxima frequencies.
* [#15392](http://bugzilla.scilab.org/show_bug.cgi?id=15392): `comet` and `comet3d` did not allow specifying colors with colors names.
* [#15425](http://bugzilla.scilab.org/show_bug.cgi?id=15425): The Kronecker product `a.*.b` failed when `a` or `b` or both are hypermatrices, with one or both being polynomials or rationals.
<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:lang="en" xml:id="ric_desc">
<refnamediv>
<refname>ric_desc</refname>
- <refpurpose>Riccati equation</refpurpose>
+ <refpurpose>Riccati equation <emphasis role="bold">(obsolete)</emphasis></refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Syntax</title>
- <synopsis>X=ric_desc(H [,E))
+ <synopsis>
+ X=ric_desc(H [,E))
[X1,X2,zero]=ric_desc(H [,E])
</synopsis>
</refsynopsisdiv>
</refsection>
<refsection>
<title>Description</title>
+ <warning>
+ <para>
+ This function is obsolete and will be removed from Scilab 6.1.x. Please use
+ <literal>riccati(H)</literal> or <literal>riccati(H,E)</literal> instead.
+ </para>
+ </warning>
<para>
Riccati solver with hamiltonian matrices as inputs.
</para>
</member>
</simplelist>
</refsection>
+ <refsection role="history">
+ <title>History</title>
+ <revhistory>
+ <revision>
+ <revnumber>6.1.0</revnumber>
+ <revdescription>
+ <literal>ric_desc()</literal> is declared obsolete. <literal>riccati()</literal>
+ replaces it.
+ </revdescription>
+ </revision>
+ </revhistory>
+ </refsection>
</refentry>
<para>
One assumes <literal>(F,G1)</literal> stabilizable and <literal>(C,F)</literal> detectable
with <literal>C'*C</literal> full rank factorization of <literal>H</literal>. Use preferably
- <literal>ric_desc</literal>.
+ <literal>riccati()</literal>.
</para>
<para>
C, D are symmetric .It is assumed that the matrices A, C and D are such that the corresponding matrix pencil has N eigenvalues with moduli less than one.
<link linkend="riccati">riccati</link>
</member>
<member>
- <link linkend="ric_desc">ric_desc</link>
- </member>
- <member>
<link linkend="schur">schur</link>
</member>
</simplelist>
<!--
* Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
* Copyright (C) INRIA -
- *
* Copyright (C) 2012 - 2016 - Scilab Enterprises
+ * Copyright (C) 2020 - Samuel GOUGEON
*
* This file is hereby licensed under the terms of the GNU GPL v2.0,
* pursuant to article 5.3.4 of the CeCILL v.2.1.
xml:lang="en" xml:id="riccati">
<refnamediv>
<refname>riccati</refname>
- <refpurpose>Riccati equation</refpurpose>
+ <refpurpose>Solves the matricial Riccati equation (continuous | discrete time domain)</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Syntax</title>
- <synopsis>X=riccati(A,B,C,dom,[typ])
- [X1,X2]=riccati(A,B,C,dom,[typ])
+ <synopsis>
+ X = riccati(H)
+ X = riccati(H, E)
+ X = riccati(A, B, C, dom)
+ X = riccati(A, B, C, dom, method)
+ [X1, X2, residual] = riccati(...)
</synopsis>
</refsynopsisdiv>
<refsection>
- <title>Arguments</title>
+ <title>Arguments</title>
<variablelist>
<varlistentry>
- <term>A,B,C</term>
+ <term>A, B, C</term>
+ <listitem>
+ Square matrices of real numbers, of size n x n: Matricial coefficients
+ of the equation.
+ <para/>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>H, E</term>
<listitem>
- <para>
- real matrices nxn, <literal>B</literal> and <literal>C</literal> symmetric.
- </para>
+ Square matrices of real numbers, of size 2n x 2n: Hamiltonian matrices | pencil
+ of the equation.
+ <para/>
</listitem>
</varlistentry>
<varlistentry>
<term>dom</term>
<listitem>
- <para>
- <literal>'c'</literal> or <literal>'d'</literal> for the time domain (continuous or discrete)
- </para>
+ Type / time domain of the Riccati equation:
+ <literal>"c"</literal> or <literal>"continuous"</literal>, or
+ <literal>"d"</literal> or <literal>"discrete"</literal>.
+ <para/>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>method</term>
+ <listitem>
+ string : <literal>"eigen"</literal> for block diagonalization (default), or
+ <literal>"schur"</literal> for Schur method.
+ <para/>
</listitem>
</varlistentry>
<varlistentry>
- <term>typ</term>
+ <term>X1, X2, X</term>
<listitem>
- <para>
- string : <literal>'eigen'</literal> for block diagonalization or <literal>schur'</literal> for Schur method.
- </para>
+ square matrices of real numbers (X2 invertible). X is symetric when B and C
+ are so.
+ <para/>
</listitem>
</varlistentry>
<varlistentry>
- <term>X1,X2,X</term>
+ <term>residual</term>
<listitem>
- <para>square real matrices (X2 invertible), X symmetric</para>
+ real number: norm<subscript>1</subscript> of the actual Equation(X) residue.
+ <para/>
</listitem>
</varlistentry>
</variablelist>
+ <para/>
</refsection>
+
<refsection>
<title>Description</title>
- <para>
- <literal>X=riccati(A,B,C,dom,[typ])</literal> solves the Riccati equation:
- </para>
- <programlisting role=""><![CDATA[
-A'*X+X*A-X*B*X+C=0
- ]]></programlisting>
- <para>
- in continuous time case, or:
- </para>
- <programlisting role=""><![CDATA[
-A'*X*A-(A'*X*B1/(B2+B1'*X*B1))*(B1'*X*A)+C-X
- ]]></programlisting>
- <para>
- with <literal>B=B1/B2*B1'</literal> in the discrete time case.
- If called with two output arguments, <literal>riccati</literal> returns <literal>X1,X2</literal>
- such that <literal>X=X1/X2</literal>.
- </para>
+ <refsect3>
+ <title>Continuous time domain</title>
+ <para>
+ <emphasis role="bold">X = riccati(A, B, C, "c")</emphasis> solves the matricial
+ Riccati equation:
+ </para>
+ <para>
+ <literal>A'*X + X*A - X*B*X + C = 0</literal>
+ </para>
+ <para>
+ <emphasis role="bold">X = riccati(H)</emphasis>
+ where <literal>H = [A, -B ; -C, -A']</literal> is the Hamiltonian matrix,
+ does the same, using the 'eigen' method, but is more stable.
+ </para>
+ </refsect3>
+ <refsect3>
+ <title>Discrete time domain</title>
+ <para>
+ <emphasis role="bold">X = riccati(A, B, C, "d")</emphasis> solves the Riccati equation:
+ </para>
+ <para>
+ <literal>A'*X*A - (A'*X*B<subscript>1</subscript>) / (B<subscript>2</subscript> +
+ B<subscript>1</subscript>'*X*B<subscript>1</subscript>) *
+ (B<subscript>1</subscript>'*X*A) + C - X = 0</literal>
+ </para>
+ <para>
+ with <literal>B = B<subscript>1</subscript> / B<subscript>2</subscript> *
+ B<subscript>1</subscript>'</literal>.
+ </para>
+ <para>
+ <emphasis role="bold">X = riccati(H, E)</emphasis> does the same, where
+ <literal>H = [A, zeros(n,n) ; -C, eye(n,n)]</literal> and
+ <literal>E = [eye(n,n), B ; zeros(n,n), A']</literal>
+ define the Hamiltonian pencil <literal>(E,H)</literal>.
+ </para>
+ </refsect3>
+ <refsect3>
+ <title>Output options</title>
+ <para>
+ <emphasis role="bold">[X1, X2] = riccati(..)</emphasis> provides
+ <literal>X1</literal> and <literal>X2</literal>, with <literal>X2</literal> invertible,
+ such that <literal>X = X1 / X2</literal>.
+ </para>
+ <para>
+ <literal>residual</literal> is the L1-norm of the actual equation's result. If
+ <literal>X</literal> is an actual solution, <literal>residual</literal> should be 0.
+ Most often its value is close to <literal>%eps*norm(X)</literal>.
+ </para>
+ <warning>
+ In the <literal>discrete</literal> case, sometimes
+ <literal>B<subscript>1</subscript></literal> and
+ <literal>B<subscript>2</subscript></literal> can't be retrieved from <literal>B</literal>.
+ Then <literal>residual</literal> can't be assessed and is set to <literal>%nan</literal>.
+ </warning>
+ </refsect3>
+ <warning>
+ When no solution is found, <literal>X=[]</literal> or <literal>X1=[]</literal> is
+ returned, and a warning is yielded.
+ </warning>
</refsection>
+
<refsection>
<title>Examples</title>
+ <emphasis role="bold">Continuous time domain</emphasis>
+ <para>
<programlisting role="example"><![CDATA[
-// Continuous
-n = 10;
-A = rand(n,n);
-B = rand(n,n);
-C = rand(n,n);
-C = C*C';
-R = rand(n,n);
-R = R*R'+eye();
-B = B*inv(R)*B';
+n = 3;
+// [A, B, C] = (grand(n,n,"uin",-9,9), grand(n,n,"uin",-9,9), grand(n,n,"uin",-9,9))
+A = [
+ -62. 91. 57.
+ -43. -45. -19.
+ 58. 83. -62.
+ ];
+B = [
+ 75. -31. -10.
+ -79. 70. 68.
+ -72. -5. 32.
+ ];
+C = [
+ -56. 70. 58.
+ -41. 54. 50.
+ 90. 2. -40.
+ ];
+// With A, B, C
+// ------------
+X = riccati(A, B, C, "c")
+clean(A'*X + X*A - X*B*X + C)
+[x1, x2, res] = riccati(A,B,C, "c");
+x = x1 / x2;
+and(x==X)
+res
-X = riccati(A,B,C,'c','eigen')
- ]]></programlisting>
- <programlisting role="example"><![CDATA[
-// Discrete
+// With the Hamiltonian
+// --------------------
+H = [A, -B; -C, -A'];
+X = riccati(H)
+clean(A'*X + X*A - X*B*X + C)
+[x1, x2] = riccati(H);
+x = x1 / x2;
+and(x==X)
+ ]]></programlisting>
+ <screen>
+--> X = riccati(A, B, C, "c")
+ X =
+ -0.1790367 0.4166084 0.2319152
+ -0.4977093 0.7993495 0.3086213
+ 0.5595286 0.3202094 -0.103394
+
+--> clean(A'*X + X*A - X*B*X + C)
+ ans =
+ 0. 0. 0.
+ 0. 0. 0.
+ 0. 0. 0.
+
+--> [x1, x2, res] = riccati(A,B,C, "c");
+--> x = x1 / x2;
+--> and(x==X)
+ ans =
+ T
+
+--> res
+ res =
+ 3.340D-13
+
+
+--> // With the Hamiltonian
+--> // --------------------
+--> H = [A, -B; -C, -A'];
+--> X = riccati(H);
+ X =
+ -0.1790367 0.4166084 0.2319152
+ -0.4977093 0.7993495 0.3086213
+ 0.5595286 0.3202094 -0.103394
-n = 10;
-F = rand(n,n);
-G1 = rand(n,n);
-G2 = rand(n,n);
-G2 = G2*G2'+eye();
-G = G1/G2*G1';
-H = rand(n,n);
-H = H*H';
+--> clean(A'*X + X*A - X*B*X + C)
+ ans =
+ 0. 0. 0.
+ 0. 0. 0.
+ 0. 0. 0.
-[X1,X2]= riccati(F,G,H,'d','schur')
- ]]></programlisting>
+--> [x1, x2] = riccati(H);
+--> x = x1 / x2;
+--> and(x==X)
+ ans =
+ T
+</screen>
+ </para>
+
+ <emphasis role="bold">Discrete time domain</emphasis>
+ <para/>
+ <programlisting role="example"><![CDATA[
+// Eq: A'*X*A - (A'*X*B1) / (B2 + B1'*X*B1) * (B1'*X*A) + C - X = 0
+n = 4;
+[A, B1, B2, C] = (rand(n,n), rand(n,n), rand(n,n), rand(n,n));
+B = B1 / B2 * B1';
+
+X = riccati(A, B, C, 'd', 'schur')
+if X <> [] then
+ clean(A'*X*A - (A'*X*B1) / (B2 + B1'*X*B1) * (B1'*X*A) + C - X)
+ [x1, x2, res] = riccati(A, B, C, 'd', 'schur');
+ x = x1 / x2;
+ and(x==X)
+ res
+else
+ disp("Riccati: No solution found")
+end
+ ]]></programlisting>
</refsection>
<refsection role="see also">
<link linkend="ricc">ricc</link>
</member>
<member>
- <link linkend="ric_desc">ric_desc</link>
- </member>
- <member>
<link linkend="srfaur">srfaur</link>
</member>
<member>
</member>
</simplelist>
</refsection>
+
+ <refsection role="history">
+ <title>History</title>
+ <revhistory>
+ <revision>
+ <revnumber>6.1.0</revnumber>
+ <revdescription>
+ <itemizedlist>
+ <listitem>
+ <literal>riccati(H)</literal> and <literal>riccati(H,E)</literal> added.
+ </listitem>
+ <listitem>
+ <literal>residual</literal> output added.
+ </listitem>
+ <listitem>
+ When no solution is found, <literal>X=[]</literal> |
+ <literal>X1=[]</literal> is now returned, without error.
+ </listitem>
+ </itemizedlist>
+ </revdescription>
+ </revision>
+ </revhistory>
+ </refsection>
</refentry>
<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:lang="ja" xml:id="ric_desc">
<refnamediv>
<refname>ric_desc</refname>
- <refpurpose>リカッチ方程式</refpurpose>
+ <refpurpose>リカッチ方程式 <emphasis role="bold">(obsolete)</emphasis></refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>呼び出し手順</title>
</refsection>
<refsection>
<title>説明</title>
+ <warning>
+ <para>
+ This function is obsolete and will be removed from Scilab 6.1.x. Please use
+ <literal>riccati(H)</literal> or <literal>riccati(H,E)</literal> instead.
+ </para>
+ </warning>
<para>
ハミルトン行列を入力とするリカッチソルバ.
</para>
</member>
</simplelist>
</refsection>
+ <refsection role="history">
+ <title>履歴</title>
+ <revhistory>
+ <revision>
+ <revnumber>6.1.0</revnumber>
+ <revdescription>
+ <literal>ric_desc()</literal> is declared obsolete. <literal>riccati()</literal>
+ replaces it.
+ </revdescription>
+ </revision>
+ </revhistory>
+ </refsection>
</refentry>
<literal>(F,G1)</literal>は可安定, <literal>(C,F)</literal> は
<literal>H</literal>のフルランク分解 <literal>C'*C</literal>を
用いて可検出であることとする.
- より適する場合, <literal>ric_desc</literal> を使用すること.
+ より適する場合, <literal>riccati()</literal> を使用すること.
</para>
<para>
C, D は対象行列である.
<link linkend="riccati">riccati</link>
</member>
<member>
- <link linkend="ric_desc">ric_desc</link>
- </member>
- <member>
<link linkend="schur">schur</link>
</member>
</simplelist>
+++ /dev/null
-<?xml version="1.0" encoding="UTF-8"?>
-<!--
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) INRIA -
- *
- * Copyright (C) 2012 - 2016 - Scilab Enterprises
- *
- * This file is hereby licensed under the terms of the GNU GPL v2.0,
- * pursuant to article 5.3.4 of the CeCILL v.2.1.
- * This file was originally licensed under the terms of the CeCILL v2.1,
- * and continues to be available under such terms.
- * For more information, see the COPYING file which you should have received
- * along with this program.
- *
- -->
-<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink"
- xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML"
- xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org"
- xml:lang="ja" xml:id="riccati">
- <refnamediv>
- <refname>riccati</refname>
- <refpurpose>リカッチ方程式</refpurpose>
- </refnamediv>
- <refsynopsisdiv>
- <title>呼び出し手順</title>
- <synopsis>X=riccati(A,B,C,dom,[typ])
- [X1,X2]=riccati(A,B,C,dom,[typ])
- </synopsis>
- </refsynopsisdiv>
- <refsection>
- <title>パラメータ</title>
- <variablelist>
- <varlistentry>
- <term>A,B,C</term>
- <listitem>
- <para>
- 実数行列 nxn, <literal>B</literal> および対称行列 <literal>C</literal> .
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>dom</term>
- <listitem>
- <para>
- 時間領域(連続または離散)に関する <literal>'c'</literal> または <literal>'d'</literal>
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>typ</term>
- <listitem>
- <para>
- 文字列 : ブロック対角化の場合は<literal>'eigen'</literal>,
- または Schur 法の場合は<literal>schur'</literal>.
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>X1,X2,X</term>
- <listitem>
- <para>正方実数行列 (X2 は可逆), 対称行列 X </para>
- </listitem>
- </varlistentry>
- </variablelist>
- </refsection>
- <refsection>
- <title>説明</title>
- <para>
- <literal>X=riccati(A,B,C,dom,[typ])</literal> は次のリカッチ方程式を解きます:
- </para>
- <programlisting role=""><![CDATA[
-A'*X+X*A-X*B*X+C=0
- ]]></programlisting>
- <para>
- (連続系の場合),または
- </para>
- <programlisting role=""><![CDATA[
-A'*X*A-(A'*X*B1/(B2+B1'*X*B1))*(B1'*X*A)+C-X
- ]]></programlisting>
- <para>
- (離散時間系の場合), ただし <literal>B=B1/B2*B1'</literal>.
- 出力引数2個でコールされた場合,
- <literal>riccati</literal> は
- <literal>X=X1/X2</literal>となるような<literal>X1,X2</literal>を返します.
- </para>
- </refsection>
- <refsection role="see also">
- <title>参照</title>
- <simplelist type="inline">
- <member>
- <link linkend="ricc">ricc</link>
- </member>
- <member>
- <link linkend="ric_desc">ric_desc</link>
- </member>
- <member>
- <link linkend="srfaur">srfaur</link>
- </member>
- <member>
- <link linkend="lindquist">lindquist</link>
- </member>
- </simplelist>
- </refsection>
-</refentry>
b=-b;
h=[-a',c'*c;
0*eye(a),a];
- x=ric_desc(h);h=[]
+ x=riccati(h);h=[]
gm=syslin("c",-a'+c'*c*x,-c',b'-d'*c*x,d')';
endfunction
Ar=A-B*Si*D'*C;
H=[Ar,-B*Si*B';
-C'*inv(R)*C,-Ar'];
- X=ric_desc(H);
+ X=riccati(H);
F=-Si*(D'*C+B'*X)
endfunction
Ar=A-B*Si*D'*C;
H=[Ar',-C'*Ri*C;
-B*Si*B',-Ar];
- Z=ric_desc(H);
+ Z=riccati(H);
H=-(B*D'+Z*C')*Ri
endfunction
indic=1;test=1;
end
if indic ==0 then
- [X1,X2,errx]=ric_desc(H);
+ [X1,X2,errx]=riccati(H);
if errx > 1.d-4 then
mprintf(gettext("%s: Riccati solution inaccurate: equation error = %g.\n"),"h_inf",errx);
end
indic=1 ;test=1;
end
if indic==0 then
- [Y1,Y2,erry]=ric_desc(J);
+ [Y1,Y2,erry]=riccati(J);
if erry > 1.d-4 then
mprintf(gettext("%s: Riccati solution inaccurate: equation error = %g.\n"),"h_inf",erry);
end
mprintf(gettext("%s: An eigenvalue of %s (controller) is close to Imaginary axis.\n"),"h_inf","H");
end
- [X1,X2,errx]=ric_desc(H);
+ [X1,X2,errx]=riccati(H);
if errx > 1.d-4 then
mprintf(gettext("%s: Riccati solution inaccurate: equation error = %g.\n"),"h_inf",errx);
end
if dy < 1.d-6 then
mprintf(gettext("%s: An eigenvalue of %s (observer) is close to Imaginary axis.\n"),"h_inf","J");
end
- [Y1,Y2,erry]=ric_desc(J);
+ [Y1,Y2,erry]=riccati(J);
if erry > 1.d-4 then
mprintf(gettext("%s: Riccati solution inaccurate: equation error = %g.\n"),"h_inf",erry);
end
//
// (solution X is also given by X=riccati(A,G,C,'d') with G=B/R*B')
//!
+ warnobsolete("riccati()", "6.1.x")
[LHS,RHS]=argn(0);
if RHS==1 then
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA -
-//
// Copyright (C) 2012 - 2016 - Scilab Enterprises
+// Copyright (C) 2020 - Samuel GOUGEON
//
// This file is hereby licensed under the terms of the GNU GPL v2.0,
// pursuant to article 5.3.4 of the CeCILL v.2.1.
// For more information, see the COPYING file which you should have received
// along with this program.
-function [x1,x2]=riccati(a,b,c,dom,typ)
- //[x1,[x2]]=riccati(a,b,c,dom,[typ]) is a Riccati solver
- // x=x1/x2 solves:
- // a'*x+x*a-x*b*x+c=0 (continuous time case)
- // a'*x*a-(a'*x*b1/(b2+b1'*x*b1))*(b1'*x*a)+c-x with b=b1/b2*b1'
- // (discrete time case)
- // If called with LHS=1 (one output argument) riccati returns x.
+function [x1, x2, residual] = riccati(a, b, c, dom, method)
+ // x = riccati(h) // continuous time domain
+ // x = riccati(h, e) // discrete time domain
+ // x = riccati(a, b, c, dom)
+ // x = riccati(a, b, c, dom, method)
+ //[x1,x2,residual] = riccati(..) such that x = x1/x2 solve :
+ // a'*x + x*a - x*b*x + c = 0 (continuous time case)
+ // a'*x*a - (a'*x*b1/(b2+b1'*x*b1)) * (b1'*x*a) + c - x = 0 (discrete time case)
+ // with b = b1 / b2 * b1'
//
- // -- a,b,c real matrices nxn, b and c symmetric.
+ // -- a,b,c real matrices of size n x n
+ // For the discrete time domain: b and c must be symmetric.
+ // -- h, e real matrices of size 2n x 2n
// -- dom = 'c' or 'd' for the time domain (continuous or discrete)
- // -- typ = 'eigen' --->block diagonalization
- // typ = 'schur' --->schur method
- // See also ric_desc
- //!
+ // -- method = 'eigen' --->block diagonalization
+ // method = 'schur' --->schur method
+
+ [x1, x2, residual] = ([], [], [])
+ [lhs,rhs] = argn(0)
+ if ~or(rhs==[1 2 4 5]) then
+ msg = _("%s: Wrong number of input arguments: %s expected.\n")
+ error(msprintf(msg, "riccati", "1,2,4 or 5"))
+ end
+ msg = gettext("%s: No solution found. Dimension(stable subspace)=%d. %d expected.\n")
+
+ if rhs > 2 then
+ if ~isdef("method","l") then
+ method = "eigen"
+ end
+ ham = [a, -b ; -c, -a']
+ [n,n] = size(a)
+ if part(dom,1)=="c" then // riccati(a, b, c, "continuous")
+ select method,
+ case "schur" then
+ [s,d] = schur(ham, "c")
+ if d <> n then
+ warning(msprintf(msg, "riccati", d, n))
+ return
+ end
+ s = s(:,1:n),
+ case "eigen" then
+ [hb,u1] = bdiag(ham),
+ [u2,d] = schur(hb,"c")
+ u = u1*u2,
+ if d <> n then
+ warning(msprintf(msg, "riccati", d, n))
+ return
+ end
+ s = u(:,1:n)
+ end
+
+ else // riccati(a, b, c, "discrete")
+ E = [eye(n,n), b ; zeros(n,n), a']
+ H = [a, zeros(n,n) ; -c, eye(n,n)]
+ [bs,as,s,d] = schur(H, E, "d");
+ if d <> n then
+ warning(msprintf(msg, "riccati", d, n))
+ return
+ end
+ s = s(:,1:n);
+ end
+ x1 = s(n+1:2*n,:)
+ x2 = s(1:n,:),
+ if lhs > 2
+ x = x1 / x2
+ if part(dom,1)=="c"
+ residual = norm(a'*x + x*a -x*b*x + c,1)
+ else
+ b1 = real(sqrtm(b))
+ b2 = eye(b)
+ if norm(b1/b2*b1'-b) / norm(b) > 1e3*%eps
+ // This can occur when b is not symetric
+ // Then sqrtm(b) may be highly complex
+ b1 = eye(b)
+ try
+ b2 = inv(b) // b may not be invertible
+ catch
+ b2 = %nan
+ end
+ end
+ // A'*X*A - (A'*X*B1) / (B2 + B1'*X*B1) * (B1'*X*A) + C - X = 0
+ residual = norm(a'*x*a - (a'*x*b1 / (b2+b1'*x*b1)) * (b1'*x*a)+c-x,1)
+ end
+ end
- [lhs,rhs]=argn(0),
- if rhs==4 then typ="eigen",end,
- ham=[a -b;-c -a'],
- [n,n]=size(a),
- if part(dom,1)=="c" then
- select typ,
- case "schur" then
- [s,d]=schur(ham,"c"),
- if d<>n then
- error(msprintf(gettext("%s: Wrong dimension (%d) of stable subspace: %d expected.\n"),"riccati",d, n))
+ else // from the Hamiltonian H and the Hamiltonian pencil (E,H)
+ H = a
+ n2 = size(H,2)
+ n1 = n2/2
+ if rhs==1 then // riccati(H): Continuous
+ A = H(1:n1,1:n1);
+ [Hb,W1] = bdiag(H);
+ if cond(W1) > 1.d10*norm(H,1) then
+ // write(%io(2),'Warning : Bad conditioning => balancing');
+ [Hb,W1] = balanc(H);
+ end
+ if cond(W1) > 1.d+10*norm(H,1) then
+ Hb = H
+ W1 = eye(W1)
end
- s=s(:,1:n),
- case "eigen" then
- [hb,u1]=bdiag(ham),
- [u2,d]=schur(hb,"c"),
- u=u1*u2,
- if d<>n then
- error(msprintf(gettext("%s: Wrong dimension (%d) of stable subspace: %d expected.\n"),"riccati",d, n))
+ [W2,n] = schur(Hb,"c")
+ if n <> n1 then
+ warning(msprintf(msg, "riccati", n, n1))
+ else
+ W1 = W1*W2
+ UV = W1(:,1:n1)
+ x2 = UV(1:n1,:)
+ x1 = UV(n1+1:n2,:)
+ if lhs > 2
+ x = x1 / x2
+ // H = [A, -b ; -c, -A']
+ b = -H(1:n1,n1+1:$)
+ c = -H(n1+1:$,1:n1)
+ residual = norm(A'*x + x*A - x*b*x + c, 1)
+ end
+ end
+
+ else // riccati(H, E): Discrete
+ E = b
+ [UV,n] = schur(H,E,"d")
+ if n <> n1 then
+ warning(msprintf(msg, "riccati", n, n1))
+ else
+ x2 = UV(1:n,1:n)
+ x1 = UV(n+1:2*n,1:n)
+ if lhs > 2 then
+ A = H(1:n,1:n)
+ B = E(1:n,n+1:2*n)
+ C = -H(n+1:2*n,1:n)
+ B1 = real(sqrtm(B))
+ B2 = eye(A)
+ if norm(B1/B2*B1'-B) / norm(B) > 1e3*%eps
+ // This can occur when B is not symetric
+ // Then sqrtm(B) may be complex with non residual imaginary part
+ B1 = eye(B)
+ try
+ B2 = inv(B)
+ // B may not be invertible => no way to assess residual
+ catch
+ B2 = %nan
+ end
+ end
+ X = x1 / x2
+ // A'*X*A - (A'*X*B1) / (B2 + B1'*X*B1) * (B1'*X*A) + C - X = 0
+ residual = norm(A'*X*A - (A'*X*B1/(B2+B1'*X*B1))*(B1'*X*A) + C - X, 1)
+ end
end
- s=u(:,1:n),
- end,
- else
- aa=[eye(n,n), b; zeros(n,n), a'],bb=[a, zeros(n,n); -c, eye(n,n)],
- [bs,as,s,n1]=schur(bb,aa,"d");
- if n1<>n then
- error(msprintf(gettext("%s: Wrong dimension (%d) of stable subspace: %d expected.\n"),"riccati",n1, n))
end
- s=s(:,1:n1);
- end,
- x1=s(n+1:2*n,:),x2=s(1:n,:),
- if lhs==1 then x1=x1/x2,end
+ end
+ if lhs==1 & x1 <> []
+ x1 = x1 / x2
+ end
+
endfunction
+++ /dev/null
-// =============================================================================
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2012 - INRIA - Serge Steer
-//
-// This file is distributed under the same license as the Scilab package.
-// =============================================================================
-// <-- CLI SHELL MODE -->
-// <-- ENGLISH IMPOSED -->
-//
-//
-// <-- Non-regression test for bug 11092 -->
-//
-// <-- Bugzilla URL -->
-// http://bugzilla.scilab.org/show_bug.cgi?id=11092
-//
-// <-- Short Description -->
-//Incorrect argument check in h_inf
-G=syslin("c",1/%s^3);
-[P,r]=macglov(G);
-assert_checktrue(execstr("[K,ro]=h_inf(P,r,0,1,30)","errcatch")==0);h_inf: P22 is stabilizable.
-h_inf: P22 is detectable.
- gama = 1.4142135624 Unfeasible (Hx hamiltonian) test = 0.30000E+01
- gama = 2.0000000000 Unfeasible (Hx hamiltonian) test = 0.15613E+01
- gama = 2.8284271247 Unfeasible (Hx hamiltonian) test = 0.13582E+01
- gama = 4.0000000000 Unfeasible (Hx hamiltonian) test = 0.20258E+01
- gama = 5.6568542495 Unfeasible (Hx hamiltonian) test = 0.16526E+02
- gama = 5.8423739467 Unfeasible (Hx hamiltonian) test = 0.60242E+02
- gama = 5.8916775545 Unfeasible (Hx hamiltonian) test = 0.19635E+03
- gama = 5.9041997314 Unfeasible (Hx hamiltonian) test = 0.45849E+03
- gama = 5.9104908357 Unfeasible (Hx hamiltonian) test = 0.13893E+04
- gama = 5.9120667565 Unfeasible (Hx hamiltonian) test = 0.28256E+04
- gama = 5.9128551898 Unfeasible (Hx hamiltonian) test = 0.58522E+04
- gama = 5.9132495247 Unfeasible (Hx hamiltonian) test = 0.12603E+05
- gama = 5.9134467218 Unfeasible (Hx hamiltonian) test = 0.29787E+05
- gama = 5.9135453277 Unfeasible (Hx hamiltonian) test = 0.93598E+05
- gama = 5.9135699799 Unfeasible (Hx hamiltonian) test = 0.20153E+06
- gama = 5.9135823062 Unfeasible (Hx hamiltonian) test = 0.47597E+06
- gama = 5.9135884693 Unfeasible (Hx hamiltonian) test = 0.14915E+07
- gama = 5.9135900101 Unfeasible (Hx hamiltonian) test = 0.31966E+07
- gama = 5.9135907805 Unfeasible (Hx hamiltonian) test = 0.74615E+07
- gama = 5.9135911657 Unfeasible (Hx hamiltonian) test = 0.22414E+08
- gama = 5.9135912620 Unfeasible (Hx hamiltonian) test = 0.44919E+08
// <-- CLI SHELL MODE -->
// <-- ENGLISH IMPOSED -->
+// <-- NO CHECK REF -->
//
//
// <-- Non-regression test for bug 11092 -->
--- /dev/null
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2020 - Samuel GOUGEON
+//
+// This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- CLI SHELL MODE -->
+// <-- NO CHECK REF -->
+//
+// <-- Non-regression test for bug 14610 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/14610
+//
+// <-- Short Description -->
+// x = riccati(h,e) yielded an error
+
+h = rand(2,2);
+e = rand(2,2);
+assert_checkequal(execstr("x = riccati(h, e)", "errcatch"), 0);
-//<-- CLI SHELL MODE -->
// =============================================================================
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) 2009 - INRIA - Serge Steer
// <-- Non-regression test for bug 4596 -->
//
+//<-- CLI SHELL MODE -->
+//<-- NO CHECK REF -->
+//
// <-- Bugzilla URL -->
// http://bugzilla.scilab.org/show_bug.cgi?id=4596
//
-//<-- CLI SHELL MODE -->
// =============================================================================
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) 2007-2008 - INRIA
+// Copyright (C) 2020 - Samuel GOUGEON
//
// This file is distributed under the same license as the Scilab package.
// =============================================================================
// Added : 25/06/2007
// Thanks to Sabine GAUZERE
// ============================================================================
+// <-- CLI SHELL MODE -->
+// <-- ENGLISH IMPOSED -->
n = 10;
A = rand(n,n);
B = rand(n,n);
C = rand(n,n);
C = C*C';
R = rand(n,n);
-R = R*R'+eye();
-B = B*inv(R)*B';
+R = R*R' + eye();
+B = B * inv(R) * B';
// Test de l'équation en temps continu
+// -----------------------------------
+// (Eq) A'*X + X*A - X*B*X + C = 0
+// From the matricial coefficients of the equation
X = riccati(A,B,C,'c','eigen');
-// Vérification
C_cont = A'*X+X*A-X*B*X;
+assert_checktrue(norm(C_cont+C,1) < (10000*%eps));
+[X1,X2,residual] = riccati(A,B,C,'c','eigen');
+assert_checkequal(X, X1/X2);
+assert_checktrue(residual < 1e-12);
+// From the Hamiltonian of the equation
+H = [A, -B ; -C, -A'];
+x = riccati(H);
+assert_checkequal(x,X);
+[x1, x2, residual] = riccati(H);
+assert_checkequal(x, x1/x2);
+assert_checktrue(residual < 1e-12);
+// No solution
+a = zeros(3,3);
+b = a;
+c = eye(3,3);
+X = riccati(a, b, c, 'c','eigen');
+WARNING: riccati: No solution found. Dimension(stable subspace)=0. 3 expected.
+assert_checkequal(X, []);
+[X1, X2, residual] = riccati(a, b, c, 'c', 'eigen');
+WARNING: riccati: No solution found. Dimension(stable subspace)=0. 3 expected.
+assert_checkequal(X1, []);
+assert_checkequal(residual, []);
// Test de l'équation en temps discret
+// -----------------------------------
+// From the matricial coefficients of the equation
F = A;
B = rand(n,n);
G1 = B;
G2 = R;
G = G1/G2*G1';
-H = C;
-[X1,X2]= riccati(F,G,H,'d','schur');
-// Vérification
-X = X1/X2;
+X = riccati(F,G,C,'d','schur');
C_disc = F'*X*F-(F'*X*G1/(G2+G1'*X*G1))*(G1'*X*F)-X;
-// Comparaison des différents résultats obtenus
-if norm(C_cont+C,1) > (10000*%eps) then bugmes();quit;end
-if norm(C_disc+H,1) > (10000*%eps) then bugmes();quit;end
+assert_checktrue(norm(C_disc+C,1) < 10000*%eps)
+ ans =
+ T
+[X1, X2, residual] = riccati(F,G,C,'d','schur');
+assert_checkequal(X, X1/X2);
+assert_checktrue(residual < 1e-12);
+// From the Hamiltonian pencil (E,H) of the equation
+H = [F, zeros(n,n) ; -C, eye(n,n)];
+E = [eye(n,n), G ; zeros(n,n), F'];
+x = riccati(H, E);
+assert_checkequal(x, X);
+[x1, x2, residual] = riccati(H, E);
+assert_checkequal(x, x1/x2);
+assert_checktrue(residual < 1e-12);
-//<-- CLI SHELL MODE -->
// =============================================================================
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) 2007-2008 - INRIA
+// Copyright (C) 2020 - Samuel GOUGEON
//
// This file is distributed under the same license as the Scilab package.
// =============================================================================
// Thanks to Sabine GAUZERE
// ============================================================================
+// <-- CLI SHELL MODE -->
+// <-- ENGLISH IMPOSED -->
+
n = 10;
A = rand(n,n);
B = rand(n,n);
C = rand(n,n);
C = C*C';
R = rand(n,n);
-R = R*R'+eye();
-B = B*inv(R)*B';
+R = R*R' + eye();
+B = B * inv(R) * B';
// Test de l'équation en temps continu
+// -----------------------------------
+// (Eq) A'*X + X*A - X*B*X + C = 0
+// From the matricial coefficients of the equation
X = riccati(A,B,C,'c','eigen');
-// Vérification
C_cont = A'*X+X*A-X*B*X;
+assert_checktrue(norm(C_cont+C,1) < (10000*%eps));
+
+[X1,X2,residual] = riccati(A,B,C,'c','eigen');
+assert_checkequal(X, X1/X2);
+assert_checktrue(residual < 1e-12);
+
+// From the Hamiltonian of the equation
+H = [A, -B ; -C, -A'];
+x = riccati(H);
+assert_checkequal(x,X);
+[x1, x2, residual] = riccati(H);
+assert_checkequal(x, x1/x2);
+assert_checktrue(residual < 1e-12);
+
+// No solution
+a = zeros(3,3);
+b = a;
+c = eye(3,3);
+X = riccati(a, b, c, 'c','eigen');
+assert_checkequal(X, []);
+[X1, X2, residual] = riccati(a, b, c, 'c', 'eigen');
+assert_checkequal(X1, []);
+assert_checkequal(residual, []);
+
// Test de l'équation en temps discret
+// -----------------------------------
+// From the matricial coefficients of the equation
F = A;
B = rand(n,n);
G1 = B;
G2 = R;
G = G1/G2*G1';
-H = C;
-[X1,X2]= riccati(F,G,H,'d','schur');
-// Vérification
-X = X1/X2;
+X = riccati(F,G,C,'d','schur');
C_disc = F'*X*F-(F'*X*G1/(G2+G1'*X*G1))*(G1'*X*F)-X;
+assert_checktrue(norm(C_disc+C,1) < 10000*%eps)
+
+[X1, X2, residual] = riccati(F,G,C,'d','schur');
+assert_checkequal(X, X1/X2);
+assert_checktrue(residual < 1e-12);
+
+// From the Hamiltonian pencil (E,H) of the equation
+H = [F, zeros(n,n) ; -C, eye(n,n)];
+E = [eye(n,n), G ; zeros(n,n), F'];
+x = riccati(H, E);
+assert_checkequal(x, X);
+[x1, x2, residual] = riccati(H, E);
+assert_checkequal(x, x1/x2);
+assert_checktrue(residual < 1e-12);
-// Comparaison des différents résultats obtenus
-if norm(C_cont+C,1) > (10000*%eps) then pause,end
-if norm(C_disc+H,1) > (10000*%eps) then pause,end
map.emplace(L"nanmin", L"min");
map.emplace(L"nanmax", L"max");
map.emplace(L"numer", L".num");
+ map.emplace(L"ric_desc", L"riccati");
map.emplace(L"square", L"replot");
map.emplace(L"sysdiag", L"blockdiag");
map.emplace(L"with_tk", L"with_module('tclsci')");
projects / scilab.git / commitdiff