1 files changed, 68 insertions, 0 deletions
diff --git a/modules/cacsd/macros/fstabst.sci b/modules/cacsd/macros/fstabst.sci
new file mode 100755
index 000000000..cf66f7221
--- /dev/null
+++ b/modules/cacsd/macros/fstabst.sci
@@ -0,0 +1,68 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA -
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [J]=fstabst(Stplant,r)
+    //[J]=fstabst(Stplant,r)
+    // Parametrization of all stabilizing feedbacks.
+    //-- Stplant=[ Stplant11 Stplant12;
+    //         Stplant21 Stplant22]
+    //  (in state-space or transfer form: conversion in state-space is
+    //   done for the computations)
+    //-- r = size of  Stplant22 subsystem (2,2) block of Stplant
+    //
+    //          j =[ j11 j12;
+    //               j21 j22]
+    // K is a stablizing controller iff K=LFT(J,r,Q) with Q stable
+    //!
+    if and(typeof(Stplant)<>["rational","state-space"]) then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n"),"fstabst",1))
+    end
+    if Stplant.dt<>"c" then
+        error(msprintf(gettext("%s: Wrong type for argument #%d: In continuous time expected.\n"),"fstabst",1))
+    end
+    flag=0;
+    if typeof(Stplant)=="rational" then Stplant=tf2ss(Stplant),flag=1;end
+
+
+    [LHS,RHS]=argn(0);
+    [a,b1,b2,c1,c2,d11,d12,d21,d22]=smga(Stplant,r),
+    Rd1=d12'*d12,
+    R12=sqrtm(Rd1);
+    Rd2=d21*d21',
+    R22=sqrtm(Rd2);
+    p=r(2),r=r(1);
+    //-------------
+    [s1,s2,t]=size(Stplant);
+    [w1,k1]=rowcomp(d12),do1=w1(k1+1:s1-r,:)',
+    if do1==[] then do1=0,end,
+    [w2,k2]=rowcomp(d21'),do2=w2(k2+1:s2-p,:),
+    if do2==[] then do2=0,end,
+    //gains f and h
+    //-------------
+    ar1=a-b2*inv(Rd1)*d12'*c1,br1=b2*inv(Rd1)*b2',
+    cr1=c1'*do1*do1'*c1,
+    xc=riccati(ar1,br1,cr1,"cont"),
+    f=inv(Rd1)*(b2'*xc+d12'*c1),
+    ar2=(a-b1*d21'*inv(Rd2)*c2)',br2=c2'*inv(Rd2)*c2,
+    cr2=b1*do2'*do2*b1',
+    xo=riccati(ar2,br2,cr2,"cont"),
+    h=(b1*d21'+xo*c2')*inv(Rd2),
+    // J:
+    //---
+    aj=a-b2*f-h*c2+h*d22*f;
+    bj=[h ,(b2-h*d22)*inv(R12)],
+    cj=[-f;
+    -inv(R22)*(c2-d22*f)],
+    dj=[0*ones(p,r)            eye(p,p)*inv(R12);
+    inv(R22)*eye(r,r)         d22],
+    //Normalization to get inner and co-inner factors
+    J=syslin("c",aj,bj,cj,dj),
+    if flag==1 then J=ss2tf(J);end
+endfunction