* Bugs 14610 15350 fixed: ric_desc() fixed & merged in riccati() 09/21209/9
Samuel GOUGEON [Wed, 8 Jan 2020 20:03:10 +0000 (21:03 +0100)]
  http://bugzilla.scilab.org/14610
  http://bugzilla.scilab.org/15350

  riccati() page updated (PDF): http://bugzilla.scilab.org/attachment.cgi?id=5052

Change-Id: I45de428c6a32dacf7979d7f6ed6342130ba03d5c

20 files changed:
scilab/CHANGES.md
scilab/modules/cacsd/help/en_US/matrix_computation/ric_desc.xml
scilab/modules/cacsd/help/en_US/matrix_computation/ricc.xml
scilab/modules/cacsd/help/en_US/matrix_computation/riccati.xml
scilab/modules/cacsd/help/ja_JP/matrix_computation/ric_desc.xml
scilab/modules/cacsd/help/ja_JP/matrix_computation/ricc.xml
scilab/modules/cacsd/help/ja_JP/matrix_computation/riccati.xml [deleted file]
scilab/modules/cacsd/macros/fspecg.sci
scilab/modules/cacsd/macros/gcare.sci
scilab/modules/cacsd/macros/gfare.sci
scilab/modules/cacsd/macros/h_inf.sci
scilab/modules/cacsd/macros/ric_desc.sci
scilab/modules/cacsd/macros/riccati.sci
scilab/modules/cacsd/tests/nonreg_tests/bug_11092.dia.ref [deleted file]
scilab/modules/cacsd/tests/nonreg_tests/bug_11092.tst
scilab/modules/cacsd/tests/nonreg_tests/bug_14610.tst [new file with mode: 0644]
scilab/modules/cacsd/tests/nonreg_tests/bug_4596.tst
scilab/modules/cacsd/tests/unit_tests/riccati.dia.ref
scilab/modules/cacsd/tests/unit_tests/riccati.tst
scilab/modules/slint/src/cpp/DeprecatedChecker.cpp

index 976b77e..1662f05 100644 (file)
@@ -125,6 +125,10 @@ Feature changes and additions
   - 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:
@@ -159,6 +163,7 @@ Obsolete functions or features
 * `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
@@ -200,6 +205,7 @@ Bug Fixes
 * [#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.
@@ -213,6 +219,7 @@ Bug Fixes
 * [#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.
index 51996d9..7ba76e0 100644 (file)
 <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>
@@ -123,4 +130,16 @@ x=ric_desc(h)
             </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>
index 66348ef..9357b8f 100644 (file)
@@ -103,7 +103,7 @@ X=F'*X*F-F'*X*G1*((G2+G1'*X*G1)^-1)*G1'*X*F+H
         <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.
@@ -170,9 +170,6 @@ norm(X2-X,1)
                 <link linkend="riccati">riccati</link>
             </member>
             <member>
-                <link linkend="ric_desc">ric_desc</link>
-            </member>
-            <member>
                 <link linkend="schur">schur</link>
             </member>
         </simplelist>
index f8692f0..972f4e2 100644 (file)
@@ -2,8 +2,8 @@
 <!--
  * 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">
@@ -121,9 +265,6 @@ H = H*H';
                 <link linkend="ricc">ricc</link>
             </member>
             <member>
-                <link linkend="ric_desc">ric_desc</link>
-            </member>
-            <member>
                 <link linkend="srfaur">srfaur</link>
             </member>
             <member>
@@ -131,4 +272,27 @@ H = H*H';
             </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>
index 20e8787..34e4846 100644 (file)
@@ -16,7 +16,7 @@
 <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>
@@ -117,4 +123,16 @@ E=[eye(n,n),G;               H=[A, 0*ones(n,n);
             </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>
index cb838d4..997a26d 100644 (file)
@@ -109,7 +109,7 @@ X=F'*X*F-F'*X*G1*((G2+G1'*X*G1)^-1)*G1'*X*F+H
             <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 は対象行列である.
@@ -161,9 +161,6 @@ norm(X2-X,1)
                 <link linkend="riccati">riccati</link>
             </member>
             <member>
-                <link linkend="ric_desc">ric_desc</link>
-            </member>
-            <member>
                 <link linkend="schur">schur</link>
             </member>
         </simplelist>
diff --git a/scilab/modules/cacsd/help/ja_JP/matrix_computation/riccati.xml b/scilab/modules/cacsd/help/ja_JP/matrix_computation/riccati.xml
deleted file mode 100644 (file)
index 90826c3..0000000
+++ /dev/null
@@ -1,104 +0,0 @@
-<?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>
index cee1214..9911643 100644 (file)
@@ -32,6 +32,6 @@ function [gm]=fspecg(g)
     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
index 547c765..3dbf613 100644 (file)
@@ -28,6 +28,6 @@ function [X,F]=gcare(Sl)
     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
index eb3cd57..db95977 100644 (file)
@@ -33,6 +33,6 @@ function [Z,H]=gfare(Sl)
     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
index 8e140d6..1c2ab5e 100644 (file)
@@ -395,7 +395,7 @@ function [P6,Kinf,tv,Uc#i,Yc#i]=h_test(P2,r,mu)
         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
@@ -420,7 +420,7 @@ function [P6,Kinf,tv,Uc#i,Yc#i]=h_test(P2,r,mu)
             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
@@ -520,7 +520,7 @@ function [Sk,polesH,polesJ]=h_contr(P,r,mu,U2i,Y2i)
         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
@@ -537,7 +537,7 @@ function [Sk,polesH,polesJ]=h_contr(P,r,mu,U2i,Y2i)
     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
index 68d5fc1..cdb1649 100644 (file)
@@ -41,6 +41,7 @@ function [X1,X2,zero]=ric_desc(H,E)
     //
     //  (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
index 3e7ba29..f74b009 100644 (file)
@@ -1,7 +1,7 @@
 // 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
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_11092.dia.ref b/scilab/modules/cacsd/tests/nonreg_tests/bug_11092.dia.ref
deleted file mode 100644 (file)
index 9cda1db..0000000
+++ /dev/null
@@ -1,42 +0,0 @@
-// =============================================================================
-// 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
index 4fcfc18..b2af9a6 100644 (file)
@@ -7,6 +7,7 @@
 
 // <-- CLI SHELL MODE -->
 // <-- ENGLISH IMPOSED -->
+// <-- NO CHECK REF -->
 //
 //
 // <-- Non-regression test for bug 11092 -->
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_14610.tst b/scilab/modules/cacsd/tests/nonreg_tests/bug_14610.tst
new file mode 100644 (file)
index 0000000..961b12e
--- /dev/null
@@ -0,0 +1,20 @@
+// =============================================================================
+// 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);
index 7b5b75b..7dca14b 100644 (file)
@@ -1,4 +1,3 @@
-//<-- CLI SHELL MODE -->
 // =============================================================================
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) 2009 - INRIA - Serge Steer
@@ -8,6 +7,9 @@
 
 // <-- Non-regression test for bug 4596 -->
 //
+//<-- CLI SHELL MODE -->
+//<-- NO CHECK REF -->
+//
 // <-- Bugzilla URL -->
 // http://bugzilla.scilab.org/show_bug.cgi?id=4596
 //
index 2b7b0d6..7db68d7 100644 (file)
@@ -1,7 +1,7 @@
-//<-- 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);
index e2df6ae..c36c2a4 100644 (file)
@@ -1,7 +1,7 @@
-//<-- 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
index a44380d..aa2d589 100644 (file)
@@ -111,6 +111,7 @@ std::unordered_map<std::wstring, std::wstring> DeprecatedChecker::initDep()
     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')");