bug 13751 fix +help revisited 80/18180/4
Serge Steer [Tue, 24 May 2016 10:29:37 +0000 (12:29 +0200)]
Change-Id: I87c206212a5de4baa37d928150f7af0514d5712a

scilab/CHANGES
scilab/modules/cacsd/help/en_US/lqg2stan.xml
scilab/modules/cacsd/macros/lqg2stan.sci
scilab/modules/cacsd/tests/unit_tests/lqg2stan.dia.ref [new file with mode: 0644]
scilab/modules/cacsd/tests/unit_tests/lqg2stan.tst [new file with mode: 0644]

index cac5d0e..35fe140 100644 (file)
@@ -376,6 +376,8 @@ In 6.0.0:
 
 * Bug #13490 fixed - histc help page fixed to match the macro (by default, normalize the result).
 
+* Bug #13751 fixed - lqg2stan returned wrong (inverted) values.
+
 * Bug #13769 fixed - t = "abc..//ghi" was parsed as a continued + comment
 
 * Bug #13810 fixed - householder(v, k*v) returned column of %nan. Input parameters were not checked. The Householder matrix could not be returned. Help pages were inaccurate and without examples. There was no householder() demo.
index 0f5358c..fe1c5c9 100644 (file)
@@ -20,7 +20,7 @@
     </refnamediv>
     <refsynopsisdiv>
         <title>Syntax</title>
-        <synopsis>[P,r]=lqg2stan(P22,bigQ,bigR)</synopsis>
+        <synopsis>[P_aug,r]=lqg2stan(P,Qxu,Qwv)</synopsis>
     </refsynopsisdiv>
     <refsection>
         <title>Arguments</title>
                 <term>P22</term>
                 <listitem>
                     <para>
-                        <literal>syslin</literal> list (nominal plant) in state-space form
+                        State space representation of the nominal plant (<literal>nu</literal> inputs, <literal>ny</literal> outputs, <literal>nx</literal> states).
                     </para>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>bigQ</term>
+                <term>Qxu</term>
                 <listitem>
                     <para>
-                        <literal>[Q,S;S',N]</literal> (symmetric) weighting matrix
+                        <literal>[Q,S;S',N]</literal> symmetric <literal>nx+nu</literal> by <literal>nx+nu</literal> weighting matrix.
                     </para>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>bigR</term>
+                <term>Qwv</term>
                 <listitem>
                     <para>
-                        <literal>[R,T;T',V]</literal> (symmetric) covariance matrix
+                        <literal>[R,T;T',V]</literal> symmetric <literal>nx+ny</literal> by <literal>nx+ny</literal> covariance matrix.
                     </para>
                 </listitem>
             </varlistentry>
                 <term>r</term>
                 <listitem>
                     <para>
-                        <literal>1</literal>x<literal>2</literal> row vector = (number of measurements, number of inputs)  (dimension of  the 2,2 part of <literal>P</literal>)
+                        Row vector <literal>[ny nu]</literal>.
                     </para>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>P</term>
+                <term>P_aug</term>
                 <listitem>
                     <para>
-                        <literal>syslin</literal> list (augmented plant)
+                        Augmented plant state space representation (see: <link linkend="syslin">syslin</link>)
                     </para>
                 </listitem>
             </varlistentry>
     <refsection>
         <title>Description</title>
         <para>
-            <literal>lqg2stan</literal>  returns the augmented plant for linear LQG (H2) controller 
-            design.
-        </para>
-        <para>
-            <literal>P22=syslin(dom,A,B2,C2)</literal> is the nominal plant; it can be in continuous 
-            time (<literal>dom='c'</literal>) or discrete time (<literal>dom='d'</literal>).
-        </para>
-        <programlisting role=""><![CDATA[ 
-. 
-x = Ax + w1 + B2u
-y = C2x + w2
- ]]></programlisting>
-        <para>
-            for continuous time plant.
-        </para>
-        <programlisting role=""><![CDATA[ 
-x[n+1]= Ax[n] + w1 + B2u
-    y = C2x + w2
- ]]></programlisting>
-        <para>
-            for discrete time plant.
-        </para>
-        <para>
-            The (instantaneous) cost function is <literal>[x' u'] bigQ [x;u]</literal>.
-        </para>
-        <para>
-            The covariance of <literal>[w1;w2]</literal> is <literal>E[w1;w2] [w1',w2'] = bigR</literal>
+            <literal>lqg2stan</literal> returns the augmented plant for linear LQG (H2) controller 
+            design problem defined by:
         </para>
+        <itemizedlist>
+            <listitem>
+                <para>
+                    The nominal plant <literal>P22</literal>:described by
+                </para>
+                <latex style="display"
+                   align="left">\left\{\begin{array}{l}\dot{x}=A x+B
+                    u +w \\y=C x+D u +v\end{array}\right. \text{ for continuous time or }
+                    \left\{\begin{array}{l}x^+=A x+B u+w\\y=C
+                    x+D u +v\end{array}\right. \text{ for discrete time.}
+                </latex>
+            </listitem>
+            <listitem>
+                <para>
+                    The (instantaneous) cost function
+                    <latex><![CDATA[\left[\begin{array}{ll} x'&
+              u'\end{array}\right] Q_{xu}\left[\begin{array}{l} x\\
+              u\end{array}\right]]]></latex>.
+            </para>
+          </listitem>
+          <listitem>
+            <para>
+              The noises covariance matrix
+              <latex><![CDATA[\mathbb{E}(\left[\begin{array}{l}w\\v\end{array}\right]
+              \left[\begin{array}{ll}w'&v'\end{array}\right])=Q_{wv}
+              ]]></latex>
+                </para>
+            </listitem>
+        </itemizedlist>
+        
+        <latex style="display"
+        align="left">\text{P_aug }=\left\{\begin{array}{l}\dot{x}=A x+B_1 W+B
+            u  \\y_1=C_1 x+D_{12}W\\y_2=-C x+D_{21} W+D u
+            \end{array}\right. \text{ for continuous time or }
+            \left\{\begin{array}{l}x^+=A x+B_1 W+B
+            u  \\y_1=C_1 x+D_{12} W \\y_2=-C x+D_{21} W+D u
+            \end{array}\right. \text{ for discrete time.}
+        </latex>
+        
+        <caution>
+            <para>
+                Up to Scilab-5.5.2 lqg2stan returns wrong inverted values
+                (see <ulink url="http://bugzilla.scilab.org/show_bug.cgi?id=13751"> bug 13751</ulink>) 
+                to obtain the good result one had to use <code>[P,r]=lqg2stan(-P,Qxu,Qwv)</code>.
+            </para>
+            <para>
+                This bug is fixed since Scilab-6.0.0, old codes must be modified accordingly.
+            </para>
+        </caution>
+    </refsection>
+    <refsection>
+        <title>Algorithm</title>
         <para>
-            If <literal>[B1;D21]</literal> is a factor of <literal>bigQ</literal>, <literal>[C1,D12]</literal>
-            is a factor of <literal>bigR</literal> and <literal>[A,B2,C2,D22]</literal> is
-            a realization of P22, then <literal>P</literal> is a realization of
-            <literal>[A,[B1,B2],[C1,-C2],[0,D12;D21,D22]</literal>.
-            The (negative) feedback computed by <literal>lqg</literal> stabilizes <literal>P22</literal>,
-            i.e. the poles of <literal>cl=P22/.K</literal> are stable.
+            If <literal>[B1;D21]</literal> is a factor of
+            <literal>Qxu</literal>, <literal>[C1,D12]</literal> is a
+            factor of <literal>Qwv</literal> (see: <link
+            linkend="fullrf">fullrf</link>) then
+            <code>P_aug=syslin(P.dt,P.A,[B1,P.B],[C1;-P.C],[0,D12;D21,P.D])</code> 
         </para>
+        
     </refsection>
+    
     <refsection>
         <title>Examples</title>
         <programlisting role="example"><![CDATA[ 
 ny=2;nu=3;nx=4;
 P22=ssrand(ny,nu,nx);
-bigQ=rand(nx+nu,nx+nu);bigQ=bigQ*bigQ';
-bigR=rand(nx+ny,nx+ny);bigR=bigR*bigR';
-[P,r]=lqg2stan(P22,bigQ,bigR);K=lqg(P,r);  //K=LQG-controller
-spec(h_cl(P,r,K))      //Closed loop should be stable
+Qxu=rand(nx+nu,nx+nu);Qxu=Qxu*Qxu';
+Qwv=rand(nx+ny,nx+ny);Qwv=Qwv*Qwv';
+[P_aug,r]=lqg2stan(P,Qxu,Qwv);
+K=lqg(P_aug,r);  //K=LQG-controller
+spec(h_cl(P_aug,r,K))      //Closed loop should be stable
 //Same as Cl=P22/.K; spec(Cl('A'))
+
 s=poly(0,'s')
 lqg2stan(1/(s+2),eye(2,2),eye(2,2))
  ]]></programlisting>
index 5fdfbce..0115dcb 100644 (file)
 // along with this program.
 
 function [P,m]=lqg2stan(P22,Q,R)
-    //P = standard plant for LQG control problem
-    //described by the triple (A,B,C)
-    //  .
-    //  x = Ax + w1 + Bu   (resp x[n+1]= ... if dom='d')
-    //
-    //  y = Cx + w2
-    //
-    //  cov(w1,w2)=R;
-    //
-    //  mininize (x,u)'Q(x,u)
-    //
+//P = standard plant for LQG control problem
+//  described by the triple (A,B,C)
+//  .
+//  x = Ax + w1 + Bu   (resp x[n+1]= ... if dom='d')
+//  y = Cx + w2
+//
+//  w1 and w2 white noises with unit variance and cov(w1,w2)=R;
+//
+//  mininize (x,u)'Q(x,u)
+//
 
     flag=0;
-    P221=P22(1);
-    if P221(1)=="r" then
-    P22=tf2ss(P22);flag=1;end
-    P22=-P22;
-    [A,B,C,D22]=P22(2:5);
-    [nx,nu]=size(B);
-    [ny,nx]=size(C);
-    Qhalf=real(sqrtm(Q));
-    Rhalf=real(sqrtm(R));
-    B1=Rhalf(1:nx,:);D21=Rhalf(nx+1:nx+ny,:);
-    B2=B;
-    C1=Qhalf(:,1:nx);D12=Qhalf(:,nx+1:nx+nu);
-    C2=C;
-    D11=0*C1*B1;
-    dom=P22(7);
+    if typeof(P22)=="rational" then
+      P22=tf2ss(P22);
+      flag=1;
+    elseif typeof(P22)=="constant" then
+      P22=syslin("c",[],[],[],1)
+    end
+    if typeof(P22)<>"state-space" then 
+      error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space expected.\n"),"lqg2stan",1))
+    end
+    //P22=-P22;//bug 13751 fix
+    [ny,nu,nx]=size(P22);
+    if or(size(Q)<>nx+nu) then
+      error(msprintf(_("%s: Wrong size for input argument #%d: %d-by-%d matrix expected.\n"),...
+                       "lqg2stan",2,nx+nu,nx+nu))
+    end
+    if or(size(R)<>nx+ny) then
+      error(msprintf(_("%s: Wrong size for input argument #%d: %d-by-%d matrix expected.\n"),...
+                       "lqg2stan",3,nx+ny,nx+ny))
+    end
+    if norm(Q.'-Q,1)>100*%eps*norm(Q,1) then
+      error(msprintf(_("%s: Wrong value for input argument #%d: Must be symmetric.\n"),"lqg2stan",2))
+    end
+    if norm(R.'-R,1)>100*%eps*norm(R,1) then
+      error(msprintf(_("%s: Wrong value for input argument #%d: Must be symmetric.\n"),"lqg2stan",3))
+    end
+    dom=P22.dt;
     if dom==[] then
         warning(msprintf(gettext("%s: Input argument #%d is assumed continuous time.\n"),"lqg2stan",1));
+        dom="c";
+    end
+    Qhalf=sqrtm(Q);
+    if norm(imag(Qhalf),1)>100*%eps then
+       error(msprintf(_("%s: Wrong value for input argument #%d: Must be a symmetric non negative matrix.\n"),"lqg2stan",2))
     end
-    P=syslin(dom,A,real([B1,B2]),real([C1;C2]),real([D11,D12;D21,D22]));
-    m=size(C2*B2);
+    Qhalf=real(Qhalf);
+    Rhalf=sqrtm(R);
+    if norm(imag(Qhalf),1)>100*%eps then
+       error(msprintf(_("%s: Wrong value for input argument #%d: Must be a symmetric non negative matrix.\n"),"lqg2stan",3))
+    end
+    Rhalf=real(Rhalf);
+    
+    B1=Rhalf(1:nx,:);D21=Rhalf(nx+1:nx+ny,:);
+    C1=Qhalf(:,1:nx);D12=Qhalf(:,nx+1:nx+nu);
+    D11=0*C1*B1;
+
+    P=syslin(dom,P22.A,real([B1,P22.B]),real([C1;P22.C]),real([D11,D12;D21,P22.D]));
+    m=[ny,nu];
     if flag==1 then
-    P=ss2tf(P);end
+      P=ss2tf(P);
+    end
 endfunction
diff --git a/scilab/modules/cacsd/tests/unit_tests/lqg2stan.dia.ref b/scilab/modules/cacsd/tests/unit_tests/lqg2stan.dia.ref
new file mode 100644 (file)
index 0000000..6908063
--- /dev/null
@@ -0,0 +1,45 @@
+//<-- CLI SHELL MODE -->
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2016 - INRIA - Serge Steer
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+s=poly(0,'s');
+P22=syslin("c",1/(s+2));
+[P,r]=lqg2stan(P22,eye(2,2),eye(2,2));
+assert_checktrue(and(P.num==[1,0,1;0,0,1;1,1,1]));
+assert_checktrue(and(P.den==[2+s,1,2+s;1,1,1;2+s,1,2+s]));
+assert_checktrue(and(r==[1 1]));
+P22=syslin("c",[1;1]/(s+2));
+[P,r]=lqg2stan(P22,eye(2,2),eye(3,3));
+Nref=[1,0,0,sqrt(2);0,0,0,1;sqrt(2)/2,1,0,1;sqrt(2)/2,0,1,1];
+Dref=[2,1,1,2,1,0,0,1;1,1,1,1,0,0,0,0;2,1,1,2,1,0,0,1;2,1,1,2,1,0,0,1];
+assert_checkalmostequal(coeff(P.num),Nref,1e-10,1e-10);
+assert_checkalmostequal(coeff(P.den),Dref,1e-10,1e-10);
+assert_checktrue(and(r==[2 1]));
+P22=syslin("c",[1 1]/s);
+[P,r]=lqg2stan(P22,eye(3,3),eye(2,2));
+Nref=[1,0,-1,-1;0,0,1,0;0,0,0,1;-1,1,1,1];
+Dref=[0,1,0,0,1,0,1,1;1,1,1,1,0,0,0,0;1,1,1,1,0,0,0,0;0,1,0,0,1,0,1,1];
+assert_checkalmostequal(coeff(P.num),Nref,1e-10,1e-10);
+assert_checkalmostequal(coeff(P.den),Dref,1e-10,1e-10);
+assert_checktrue(and(r==[1 2]));
+A=[1 2;-1 0];B=[1;0];C=[0 1];
+P22=syslin("c",A,B,C);
+[P,r]=lqg2stan(P22,eye(3,3),eye(3,3));
+assert_checkequal(P.A,P22.A);
+assert_checkequal([1,0,0,1;0,1,0,0],P.B);
+assert_checkequal([1,0;0,1;0,0;0,1],P.C);
+assert_checkequal([0,0,0,0;0,0,0,0;0,0,0,1;0,0,1,0],P.D);
+assert_checkequal(P22,P($,$));
+assert_checktrue(and(r==[1 1]));
+A=[1 2;-1 0];B=[1 0;1 1];C=[0 1];
+P22=syslin("c",A,B,C);
+[P,r]=lqg2stan(P22,eye(4,4),eye(3,3));
+assert_checkequal(P.A,P22.A);
+assert_checkequal([1,0,0,1,0;0,1,0,1,1],P.B);
+assert_checkequal([1,0;0,1;0,0;0,0;0,1],P.C);
+assert_checkequal([0,0,0,0,0;0,0,0,0,0;0,0,0,1,0;0,0,0,0,1;0,0,1,0,0],P.D);
+assert_checkequal(P22,P($,$-1:$));
+assert_checktrue(and(r==[1 2]));
diff --git a/scilab/modules/cacsd/tests/unit_tests/lqg2stan.tst b/scilab/modules/cacsd/tests/unit_tests/lqg2stan.tst
new file mode 100644 (file)
index 0000000..a54e661
--- /dev/null
@@ -0,0 +1,50 @@
+//<-- CLI SHELL MODE -->
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2016 - INRIA - Serge Steer
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+s=poly(0,'s');
+P22=syslin("c",1/(s+2));
+[P,r]=lqg2stan(P22,eye(2,2),eye(2,2));
+assert_checktrue(and(P.num==[1,0,1;0,0,1;1,1,1]));
+assert_checktrue(and(P.den==[2+s,1,2+s;1,1,1;2+s,1,2+s]));
+assert_checktrue(and(r==[1 1]));
+
+P22=syslin("c",[1;1]/(s+2));
+[P,r]=lqg2stan(P22,eye(2,2),eye(3,3));
+
+Nref=[1,0,0,sqrt(2);0,0,0,1;sqrt(2)/2,1,0,1;sqrt(2)/2,0,1,1];
+Dref=[2,1,1,2,1,0,0,1;1,1,1,1,0,0,0,0;2,1,1,2,1,0,0,1;2,1,1,2,1,0,0,1];
+assert_checkalmostequal(coeff(P.num),Nref,1e-10,1e-10);
+assert_checkalmostequal(coeff(P.den),Dref,1e-10,1e-10);
+assert_checktrue(and(r==[2 1]));
+
+P22=syslin("c",[1 1]/s);
+[P,r]=lqg2stan(P22,eye(3,3),eye(2,2));
+Nref=[1,0,-1,-1;0,0,1,0;0,0,0,1;-1,1,1,1];
+Dref=[0,1,0,0,1,0,1,1;1,1,1,1,0,0,0,0;1,1,1,1,0,0,0,0;0,1,0,0,1,0,1,1];
+assert_checkalmostequal(coeff(P.num),Nref,1e-10,1e-10);
+assert_checkalmostequal(coeff(P.den),Dref,1e-10,1e-10);
+assert_checktrue(and(r==[1 2]));
+
+A=[1 2;-1 0];B=[1;0];C=[0 1];
+P22=syslin("c",A,B,C);
+[P,r]=lqg2stan(P22,eye(3,3),eye(3,3));
+assert_checkequal(P.A,P22.A);
+assert_checkequal([1,0,0,1;0,1,0,0],P.B);
+assert_checkequal([1,0;0,1;0,0;0,1],P.C);
+assert_checkequal([0,0,0,0;0,0,0,0;0,0,0,1;0,0,1,0],P.D);
+assert_checkequal(P22,P($,$));
+assert_checktrue(and(r==[1 1]));
+
+A=[1 2;-1 0];B=[1 0;1 1];C=[0 1];
+P22=syslin("c",A,B,C);
+[P,r]=lqg2stan(P22,eye(4,4),eye(3,3));
+assert_checkequal(P.A,P22.A);
+assert_checkequal([1,0,0,1,0;0,1,0,1,1],P.B);
+assert_checkequal([1,0;0,1;0,0;0,0;0,1],P.C);
+assert_checkequal([0,0,0,0,0;0,0,0,0,0;0,0,0,1,0;0,0,0,0,1;0,0,1,0,0],P.D);
+assert_checkequal(P22,P($,$-1:$));
+assert_checktrue(and(r==[1 2]));