fix proposal for bug 13831: ss2ss does not transform initial state
[scilab.git] / scilab / modules / cacsd / macros / ss2ss.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA -
3 //
4 // Copyright (C) 2012 - 2016 - Scilab Enterprises
5 //
6 // This file is hereby licensed under the terms of the GNU GPL v2.0,
7 // pursuant to article 5.3.4 of the CeCILL v.2.1.
8 // This file was originally licensed under the terms of the CeCILL v2.1,
9 // and continues to be available under such terms.
10 // For more information, see the COPYING file which you should have received
11 // along with this program.
12
13 function [Sl1,right,left]=ss2ss(Sl,T,F,G,flag)
14     // State-space to state-space conversion
15     // Returns the linear system Sl1=[A1,B1,C1,D1]
16     // where A1=inv(T)*A*T,B1=inv(T)*B,C1=C*T,D1=D.
17     // Optional parameters F and G are state feedback
18     // and output injection respectively. For example,
19     // Sl1=ss2ss(Sl,T,F) returns Sl1=[A1,B1,C1,D1] with
20     // A1=inv(T)*(A+B*F)*T;B1=inv(T)*B;C1=(C+D*F)*T;D1=D;
21     // If F is given as input then right is a non singular
22     // linear system such that Sl1=Sl*right.
23     // Sl1*invsyslin(right) is a factorization of Sl.
24     // Idem for left: if F and G are given, Sl1=left*Sl*right.
25     // Example: Sl=ssrand(2,2,5); trzeros(Sl);
26     // Sl1=ss2ss(Sl,rand(5,5),rand(2,5),rand(5,2));
27     // trzeros(Sl1), trzeros(rand(2,2)*Sl1*rand(2,2))
28     // See also : projsl
29
30     [A,B,C,D]=abcd(Sl);
31     [LHS,RHS]=argn(0);
32     select RHS
33     case 2 then
34         Ti=inv(T)
35         Sl1=syslin(Sl.dt,Ti*A*T,Ti*B,C*T,D,Ti*Sl.X0);
36         right=eye(A);left=right;
37     case 3 then
38         A1=A+B*F;C1=C+D*F;
39         A1=inv(T)*A1*T;B1=inv(T)*B;C1=C1*T;D1=D
40         Sl1=syslin(Sl.dt,A1,B1,C1,D1,Sl.X0);
41         right=syslin(Sl.dt,A+B*F,B,F,eye(F*B));
42         left=eye(size(C,1),size(C,1));
43     case 4 then
44         A1=A+B*F+G*C+G*D*F;C1=C+D*F;B1=B+G*D
45         A1=inv(T)*A1*T;B1=inv(T)*B1;C1=C1*T;D1=D
46         Sl1=syslin(Sl.dt,A1,B1,C1,D1,Sl.X0);
47         right=syslin(Sl.dt,A+B*F,B,F,eye(F*B));
48         // Warning left is computed as [ At + G*Ct,G;Ct,I]
49         // where [At Bt; Ct Dt] is Sl1*right
50         At=A+B*F; Ct=C+D*F
51         left=syslin(Sl.dt,At+G*Ct,G,Ct,eye(Ct*G));
52     case 5 then
53         if flag==1 then
54             // x in R^n , y in R^p, u in R^m
55             // output injection [ A + GC, (B+GD,-G)]
56             //                  [   C   , (D   , 0)]
57             // then feeback ( since output injection increase the
58             // size of the feedback the F matrix must be of size
59             // (m+p,n) --> F=[F1;F2] with F1 (m,n) and F2 (p,n)
60             //
61             // Sl1= [ A+GC +BF1+G*D*F1 -GF2, (B+GD,-G)]
62             //  [ C+D*F1               , (D   , 0)]
63             //
64             // We have then the following property
65             // Sl1 equiv  left*sysdiag(sys,eye(p,p))*right
66             //
67             //
68             n=size(A,"r");p=size(C,"r");
69             A1=A+G*C+[B+G*D,-G]*F;B1=[B+G*D,-G];C1=C+[D,zeros(p,p)]*F;
70             D1=[D,zeros(p,p)];
71             A1=inv(T)*A1*T;B1=inv(T)*B1;C1=C1*T;D1=D1
72             Sl1=syslin(Sl.dt,A1,B1,C1,D1,Sl.X0);
73             left=syslin(Sl.dt,A+G*C,[G,-G],C,[eye(p,p),zeros(p,p)]);
74             // Now we compute the right associated to left*Sl1
75             A1=A+G*C;B1=[B+G*D,-G];C1=C;D1=[D,zeros(p,p)];
76             right=syslin(Sl.dt,A1+B1*F,B1,F,eye(F*B1));
77             return
78         end
79         if flag==2 then
80             // x in R^n , y in R^p, u in R^m
81             // feedback first F of size(m,n)
82             //                      [ A+BF,B]
83             //                  [ C+DF,D]
84             // then output injection
85             // Sl1= [ A+GC +BF+G*D*F, (B+GD,-G)]
86             //  [ C+D*F         , (D   , 0)]
87             // this is a generalization of the case 4
88             // We have then the following property
89             // Sl1 equiv left*sysdiag(sys*right,eye(p,p)))
90             //
91             p = size(C,"r");
92             A1=A+B*F+G*C+G*D*F;
93             B1=[B+G*D,-G];
94             C1=C+ D*F;
95             D1=[D,zeros(p,p)];
96             A1=inv(T)*A1*T;B1=inv(T)*B1;C1=C1*T;D1=D1
97             Sl1=syslin(Sl.dt,A1,B1,C1,D1,Sl.X0);
98             right=syslin(Sl.dt,A+B*F,B,F,eye(F*B));
99             // Warning left is computed as [ At + G*Ct,(G,-G);
100             //                             [    Ct    ,(I, 0)]
101             // where [At Bt; Ct Dt] is Sl1*right
102             At=A+B*F; Ct=C+D*F
103             left=syslin(Sl.dt,At+G*Ct,[G,-G],Ct,[eye(Ct*G),zeros(Ct*G)]);
104         end
105     end
106 endfunction