diff options
Diffstat (limited to 'modules/cacsd/macros/ssrand.sci')
| -rwxr-xr-x | modules/cacsd/macros/ssrand.sci | 275 |
1 files changed, 275 insertions, 0 deletions
diff --git a/modules/cacsd/macros/ssrand.sci b/modules/cacsd/macros/ssrand.sci new file mode 100755 index 000000000..344e3465d --- /dev/null +++ b/modules/cacsd/macros/ssrand.sci @@ -0,0 +1,275 @@ +// 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 [sl,U]=ssrand(nout,nin,nstate,flag) + //flag=list('co',dim_cont_subs) + //flag=list('uo',dim_unobs_subs) + //*******************controllable-unobservable*********** + //display a general Kalman form + //flag=list('ncno',dim_cno,dim_ncno,dim_co,dim_nco) + // + // (A,B,C) U-similar to + // + // [*,*,*,*; [*; + // [0,*,*,*; [0; + // A= [0,0,*,*; B= [*; + // [0,0,0,*] [0]; + // + // C= [0,0,*,*] D= [0] + // + // (A11,B1) --> controllable, unobservable modes , dimension flag(2) + // (A22,B2) --> uncontrolable,unobservables modes, dimension flag(3) + // (A33,B3) --> controllable,observable modes , dimension flag(4) + // (A44,B4) --> uncontrollable,observable modes, dimension flag(5) + //*******************stabilizability********************* + //flag=list('st',dim_cont_subs,dim_stab_subs,dim_stab0) + // dim_cont_subs<=dim_stab_subs<=dim_stab0 + //pair (A,B) U-similar to: + // [*,*,*,*; [*; + // [0,s,*,*; [0; + //A= [0,0,i,*; B=[0; + // [0,0,0,u] [0] + // + // (A11,B1) controllable s=stable matrix i=neutral matrix u=unstsble matrix + //********************detectability*********************** + //flag=list('dt',dim_inst_unob,dim_instb0,dim_unobs) + // dim_inst_unob<=dim_instb0<=dim_unobs + // + //pair (A,C) U-similar to: + // + // [u,*,*,*; + // [0,i,*,*; + //A= [0,0,s,*; + // [0,0,0,*] + // + //C= [0,0,0,*] + // + // (A44,C4) observable s=stable matrix i=neutral matrix u=unstsble matrix + // + //********************output-nulling*********************** + //flag=list('on',nr,ng,ng0,nv,rk) + // nr<=ng<=ng0<=nv<=nstate and rk<=nu (rk<nu if nr>0) + // nstate-nv>=rk ny>=rk + // + //system (A,B,C,D) U-equivalent to + // + // [*,*,*,*,*; [0,*; + // [0,s,*,*,*; [0,0; + //A= [0,0,n,*,*; B=[0,0; + // [0,0,0,u,*; [0,0; + // [*,*,*,*,*] [*,0] + // + //C= [*,*,*,*,*] D=[*,0] + // + // with (A11,B21) controllable + // A22 stable + // A33 neutral + // A44 unstable + // [B51;D11] full column rank + // range([A51,A52,A53,A54;C1,C2,C3,C4]) included in range([B51;D11]) + // dimension A11=nr + // A22=ng-nr + // A33=ng0-ng + // A44=nv-ng0 + // A55=nstate-nv + // column dimension of B1. = column dimension of D1. =rk + //**********************'unknown-input'*********************** + //flag=list('ui',nw,nwu,nwui,nwuis,rk) + // nw<=nwu<=nwui<=nwuis<=nstate rk<=nout + // + //system (A,B,C,D) U-equivalent to + // + // [*,*,*,*,*; [*; + // [*,u,*,*,*; [*; + //A= [*,0,i,*,*; B=[*; + // [*,0,0,s,*; [*; + // [*,0,0,0,*] [*] + // + // [0,0,0,0,*; [0; + //C= [*,0,0,0,0] D=[*; + // + // with (A55,C15) observable + // [C21,d2] full row-rank + // row-image([[A21;A31;A41;A51]],[B2;B3;B4;B5]) + // included in row-image([C21,D2]); + // A22 unstable + // A33 neutral + // A44 stable + // dimension A11=nw + // dimension A22=nwu-nw + // dimension A33=nwui-nwu + // dimension A44=nwuis-nwui + // row dimension C2.= row dimension of D2. =rk + //*************************************************************** + + deff("[w]=st_able(w,margin)","if w~=[] then w=w-(max(real(spec(w)))+margin)*eye();end") + margin=0.5; //M "stable" will mean real-part(M) < -margin + [lhs,rhs]=argn(0) + //rand('seed',0) + rand("normal") + if rhs==3 then flag=[];end + select flag(1) + case [] + sl=syslin("c",rand(nstate,nstate),rand(nstate,nin),rand(nout,nstate),.. + zeros(nout,nin),zeros(nstate,1));U=eye(nstate,nstate); + case "fullV" + //C^(-1)(D) = X + rd=flag(2) //rank(D) needs rd<nin + sl=syslin("c",rand(nstate,nstate),rand(nstate,nin),... + [rand(rd,nstate);zeros(nout-rd,nstate)],[rand(rd,nin);zeros(nout-rd,nin)]); + case "noW" + // B ker(D)=0 + rd=flag(2) //rank(D) needs nout>rd or C=0 + sl=syslin("c",rand(nstate,nstate),[zeros(nstate,nin-rd),rand(nstate,rd)],... + rand(nout,nstate),[zeros(nout,nin-rd),rand(nout,rd)]); + case "fullW" + //admits polynomial inverse + iv=rand(nin,nout) + for k=1:nstate + iv=iv+(%s^k)*rand(nin,nout); + end + sl=minss(inv(tf2ss(iv))); + case "co" + nc=flag(2); + iii=nstate-nc; + a=[rand(nc,nc),rand(nc,iii); + zeros(iii,nc),rand(iii,iii)]; + b=[rand(nc,nin); + zeros(iii,nin)];c=rand(nout,nstate); + U=rand(nstate,nstate); + sl=ss2ss(syslin("c",a,b,c),U); + case "uo" + no=flag(2) + iii=nstate-no; + a=[rand(no,no),rand(no,iii);zeros(iii,no),rand(iii,iii)]; + c=[zeros(nout,no),rand(nout,iii)];b=rand(nstate,nin); + U=rand(nstate,nstate); + sl=ss2ss(syslin("c",a,b,c),U); + case "ncno" + cno=flag(2);ncno=flag(3);co=flag(4);nco=flag(5); + a=sysdiag(rand(cno,cno),rand(ncno,ncno),rand(co,co),rand(nco,nco)); + b=[rand(cno,nin); + zeros(ncno,nin); + rand(co,nin); + zeros(nco,nin)]; + c=[zeros(nout,cno+ncno),rand(nout,nstate-cno-ncno)]; + U=rand(nstate,nstate); + sl=ss2ss(syslin("c",a,b,c),U); + case "st" + [nc,ns,ns0]=flag(2:4);nsc=ns-nc;nsc0=ns0-ns;nxs=nstate-ns; + a=sysdiag(rand(nc,nc),imag_axis(ns-nc,ns0-ns,nstate-ns0)), + b=[rand(nc,nin);zeros(nstate-nc,nin)]; + c=rand(nout,nstate); + U=rand(nstate,nstate); + sl=ss2ss(syslin("c",a,b,c),U); + case "dt" + [n1,n2,n3]=flag(2:4); + w=imag_axis(n3-n2,n2-n1,n1,"uis"); + a=sysdiag(w,rand(nstate-n3,nstate-n3)); + c=[zeros(nout,n3),rand(nout,nstate-n3)];b=rand(nstate,nin); + U=rand(nstate,nstate); + sl=ss2ss(syslin("c",a,b,c),U); + case "on" + [nr,ng,ng0,nv,rk]=flag(2:6); + p1=1;p2=1; + a11=sysdiag(rand(nr,nr),imag_axis(ng-nr,ng0-ng,nv-ng0)); + a12=rand(nv,nstate-nv); + a21=rand(nstate-nv,nv);a22=rand(nstate-nv,nstate-nv); + b11=zeros(nv,rk);b12=[rand(nr,nin-rk);zeros(nv-nr,nin-rk)]; + b22=zeros(nstate-nv,nin-rk);d2=zeros(nout,nin-rk); + c1=rand(nout,nv);c2=rand(nout,nstate-nv); + b21d1=rand(nstate-nv+nout,rk); + if rk<>0 then a21c1=b21d1*rand(rk,nv);else a21c1=c1;end + b21=b21d1(1:nstate-nv,:);d1=b21d1(nstate-nv+1:nstate-nv+nout,:); + a21=a21c1(1:nstate-nv,:);c1=a21c1(nstate-nv+1:nstate-nv+nout,:); + + if rk>=nout then + error(msprintf(gettext("%s: Wrong values for input argument #%d: you must choose rk<nout.\n"),"ssrand",4)) + end + // if rk>=nout then j=-b21*pinv(d1);a22=-j*c2;end + [pp,qq]=size(a21);a21(p1+1:pp,1:nv)=0*a21(p1+1:pp,1:nv); + b21(p1+1:pp,:)=0*b21(p1+1:pp,:); + [pp,qq]=size(c1);c1(p2+1:pp,:)=0*c1(p2+1:pp,:); + d1(p2+1:pp,:)=0*d1(p2+1:pp,:); + w=syslin("c",[a11,a12;a21,a22],[b11,b12;b21,b22],[c1,c2],[d1,d2]); + U=eye(nstate,nstate);sl=ss2ss(w,U); + case "ri" + [nv,nr,nw,rk]=flag(2:5) + n1=nv-nr;n2=nr;n3=nstate-nv; + a11=rand(n1,n1);a12=zeros(n1,n2);a21=rand(n2,n1); + a22=rand(n2,n2);a23=rand(n2,n3); + a33=rand(n3,n3); + b2=rand(n2,nin);b31=rand(n3,rk);b32=zeros(n3,nin-rk); + d1=rand(nout,rk);d2=zeros(nout,nin-rk); + c3=rand(nout,n3); + a31=b31*rand(rk,n1);a32=b31*rand(rk,n2); + c1=d1*rand(rk,n1);c2=d1*rand(rk,n2); + a13=rand(n1,nout)*c3;b1=rand(n1,nout)*[d1,d2]; + b3=[b31,b32];d=[d1,d2] + sl=syslin("c",[a11,a12,a13;a21,a22,a23;a31,a32,a33],[b1;b2;b3],[c1,c2,c3],d); + case "ui" + [nw,nwu,nwui,nwuis,rk]=flag(2:6) + a11=rand(nw,nw); + a22=sysdiag(imag_axis(nwu-nw,nwui-nwu,nwuis-nwui,"uis"),... + rand(nstate-nwuis,nstate-nwuis)) + a12=rand(nw,nstate-nw); + b1=rand(nw,nin); + c21d2=rand(rk,nw+nin); + a21b2=rand(nstate-nw,rk)*c21d2; + a21=a21b2(:,1:nw);b2=a21b2(:,nw+1:nw+nin); + c21=[c21d2(:,1:nw),zeros(rk,nwuis-nw)]; + d2=c21d2(:,nw+1:nw+nin); + c11=zeros(nout-rk,nwuis);c12=rand(nout-rk,nstate-nwuis); + d1=zeros(nout-rk,nin);c22=zeros(rk,nstate-nwuis); + w=syslin("c",[a11,a12;a21,a22],[b1;b2],[c11,c12;c21,c22],[d1;d2]); + U=rand(nstate,nstate); + sl=ss2ss(w,U); + end + rand("uniform"); + //rand('seed',0) + +endfunction +function w=imag_axis(ns,nn,nu,flag); + //w=random block diagonal matrix with ns stable evals + //nn imaginary-axis evals and nu unstable evals. + //flag='uis' blocks along main diagonal appear in the + // order "unstable", "imaginary axis", "stable". + [LHS,RHS]=argn(0); + if RHS==3 then flag="siu";end + if flag=="siu" then + deff("[w]=st_able(w,margin)","if w~=[] then w=w-(max(real(spec(w)))+margin)*eye();end") + margin=0.5; //M "stable" will mean real-part(M) < -margin + w=[];k=int(nn/2); + rand("normal"); + //rand('seed',0); + w=sysdiag(w,st_able(rand(ns,ns),margin)); + x=rand(1,k); + for u=x + w=sysdiag(w,[u,-2*u;u,-u]);end + if abs(2*k-nn)>1.d-5 then w=sysdiag(w,0);end + w=sysdiag(w,-st_able(rand(nu,nu),margin)); + rand("uniform"); + //rand('seed',0); + return + end + if flag=="uis" then + deff("[w]=st_able(w,margin)","if w~=[] then w=w-(max(real(spec(w)))+margin)*eye();end") + w=[];k=int(nn/2); + rand("normal"); + //rand('seed',0); + w=sysdiag(w,-st_able(rand(nu,nu),margin)); + x=rand(1,k); + for u=x + w=sysdiag(w,[u,-2*u;u,-u]);end + if abs(2*k-nn)>1.d-5 then w=sysdiag(w,0);end + w=sysdiag(w,st_able(rand(ns,ns),margin)); + rand("uniform"); + //rand('seed',0); + end +endfunction |