1 files changed, 96 insertions, 0 deletions
diff --git a/modules/cacsd/macros/observer.sci b/modules/cacsd/macros/observer.sci
new file mode 100755
index 000000000..65d884866
--- /dev/null
+++ b/modules/cacsd/macros/observer.sci
@@ -0,0 +1,96 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - F. Delebecque
+//
+// 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 [Obs,U,m]=observer(Sys,flag,alfa)
+    //Obs=observer for (observable part of) linear system Sys
+    //Obs is a linear system with matrices [Ao,Bo,Identity].
+    //where Ao is no x no, Bo is no x (nu+ny) and Co is no x no
+    //and no=nx-m;
+    //input to Obs is [u,y] (assuming Sys: dotx=A x + Bu, y=Cx + Du)
+    //output of Obs is:
+    // xhat=estimate of x modulo unobservable subsp. (case 'pp')
+    // or
+    // xhat=estimate of x modulo unstable unobservable subsp. (case 'st')
+    //
+    //case flag='st':
+    // z=H*x can be estimated with stable observer iff H*U(:,1:m) = 0
+    // assignable poles of the observer are set to alfa(1),alfa(2),...
+    //
+    //case flag='pp':
+    // z=H*x can be estimated with given error spectrum iff H*U(:,1:m)=0
+    //all poles of the observer are assigned and set to alfa(1),alfa(2),...
+    //
+    //If H satifies the constraint: H*U(:,1:m)=0  (ker(H) contains unobs-subsp
+    //of Sys) one has H*U=[0,H2] and the observer for
+    // z=H*x is is H2*Obs with H2=H*U(:,m+1:nx) i.e. Co, the C-matrix of the
+    // observer for H*x, is Co=H2.
+    //
+    //EXAMPLE:
+    // nx=5;nu=1;ny=1;un=3;us=2;Sys=ssrand(ny,nu,nx,list('dt',us,us,un));
+    // nx=5 states, nu=1 input, ny=1 output,
+    // un=3 unobservable states, us=2 of them unstable.
+    // [Obs,U,m]=observer(Sys);  Stable observer (default)
+    // W=U';H=W(m+1:nx,:);[A,B,C,D]=abcd(Sys);  //H*U=[0,eye(no,no)];
+    // Sys2=ss2tf(syslin('c',A,B,H))  //Transfer u-->z
+    // Idu=eye(nu,nu);ss2tf(Obs*sysdiag(Idu,Sys)*[Idu;Idu])
+    // Transfer u-->[u;u]-->w=[u;y=Sys*u]-->Obs*w  i.e. u-->output of Obs
+    // this transfer must equal Sys2, the u-->z transfer  (H2=eye).
+
+    [LHS,RHS]=argn(0);
+    if typeof(Sys)<>"state-space" then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space expected.\n"),"observer",1))
+    end
+
+    [nx,nx]=size(Sys.A);
+    td=Sys.dt;x0=Sys.X0;
+
+    if RHS<>2 then [m1,m2,U,sl2]=dt_ility(Sys);end
+    if RHS==1 then
+        flag="st";alfa=-ones(1,nx);
+    end
+    if RHS==2 then
+        //poles are not given-->set to -ones
+        alfa=-ones(1,nx);
+        [A,B,C,D]=abcd(Sys);  //
+        J=flag;
+        //  F=A+J*C;G=[B+J*D,-J]; //
+        Obs=syslin(td,A+J*C,[B+J*D,-J],eye(A));U=[];m=[];return;  //Ao
+    end
+    if RHS==3 then
+        if size(alfa,"*")==1 then alfa=alfa*ones(1,nx);end
+    end
+    select flag
+    case "pp"
+        m=m2;
+        no=nx-m;
+        alfa=alfa(1:no);
+        [A,B,C,D]=abcd(sl2);
+        Ao=A(m+1:nx,m+1:nx);
+        Bo=B(m+1:nx,:);
+        Co=C(:,m+1:nx);
+        Do=D;
+        J=-ppol(Ao',Co',alfa);J=J';
+        F=Ao+J*Co;G=[Bo+J*Do,-J];
+        Obs=syslin(td,F,G,eye(Ao));
+        return;
+    case "st"
+        m=m1;
+        no=nx-m;
+        alfa=alfa(1:no);
+        [A,B,C,D]=abcd(sl2);
+        Ao=A(m+1:nx,m+1:nx);
+        Bo=B(m+1:nx,:);
+        Co=C(:,m+1:nx);
+        Do=D;
+        J=stabil(Ao',Co',alfa);J=J';
+        F=Ao+J*Co;G=[Bo+J*Do,-J];
+        Obs=syslin(td,F,G,eye(Ao));
+        return;
+    end
+endfunction