Previous Next Contents

4  General System and Control macros

4.1   abcd state-space matrices




CALLING SEQUENCE :

[A,B,C,D]=abcd(sl)



PARAMETERS :




DESCRIPTION :

returns the A,B,C,D matrices from a linear system Sl.

Utility function. For transfer matrices Sl is converted into state-space form by tf2ss.

The matrices A,B,C,D are the elements 2 to 5 of the syslin list Sl, i.e. [A,B,C,D] = Sl(2:5) .


EXAMPLE :

A=diag([1,2,3]);B=[1;1;1];C=[2,2,2];
sys=syslin('c',A,B,C);
sys("A")
sys("C")
[A1,B1,C1,D1]=abcd(sys);
A1
systf=ss2tf(sys);
[a,b,c,d]=abcd(systf)
spec(a)
c*b-C*B
c*a*b-C*A*B



See Also : syslin X, ssrand X

4.2   abinv AB invariant subspace




CALLING SEQUENCE :

[X,dims,F,U,k,Z]=abinv(Sl,alfa,beta)



PARAMETERS :




DESCRIPTION :

Output nulling subspace (maximal unobservable subspace) for Sl = linear system defined by a syslin list containing the matrices [A,B,C,D] of Sl. The vector dims=[dimR,dimVg,dimV,noc,nos] gives the dimensions of subspaces defined as columns od X. The dimV first columns of X i.e V=X(:,1:dimV), span the AB-invariant subspace of Sl i.e the unobservable subspace of (A+B*F,C+D*F). (dimV=nx iff C^(-1)(D)=X).

The dimR first columns of X i.e. R=X(:,1:dimR) spans the controllable part of Sl in V, (dimR<=dimV). (dimR=0 for a left invertible system). R is the maximal controllability subspace of Sl in kernel(C).

The dimVg first columns of X spans Vg=maximal AB-stabilizable subspace of Sl. (dimR<=dimVg<=dimV).

F is a decoupling feedback: for X=[V,X2] (X2=X(:,dimV+1:nx)) one has X2'*(A+B*F)*V=0 and (C+D*F)*V=0.

The zeros od Sl are given by : X0=X(:,dimR+1:dimV); spec(X0'*(A+B*F)*X0) i.e. there are dimV-dimR closed-loop fixed modes.

If the optional parameter alfa is given as input, the dimR controllable modes of (A+BF) in V are set to alfa (or to [alfa(1), alfa(2), ...]. (alfa can be a vector (real or complex pairs) or a (real) number). Default value alfa=-1.

If the optional real parameter beta is given as input, the noc-dimV controllable modes of (A+BF) "outside" V are set to beta (or [beta(1),beta(2),...]). Default value beta=-1.

In the X,U bases, the matrices [X'*(A+B*F)*X,X'*B*U;(C+D*F)*X,D*U] are displayed as follows: 

[A11,*,*,*,*,*]  [B11 * ]
[0,A22,*,*,*,*]  [0   * ]
[0,0,A33,*,*,*]  [0   * ]
[0,0,0,A44,*,*]  [0  B42]
[0,0,0,0,A55,*]  [0   0 ]
[0,0,0,0,0,A66]  [0   0 ]

[0,0,0,*,*,*]    [0  D2]
where the X-partitioning is defined by dims and the U-partitioning is defined by k.

A11 is (dimR x dimR) and has its eigenvalues set to alfa(i)'s. The pair (A11,B11) is controllable and B11 has nu-k columns. A22 is a stable (dimVg-dimR x dimVg-dimR) matrix. A33 is an unstable (dimV-dimVg x dimV-dimVg) matrix (see st_ility).

A44 is (noc-dimV x noc-dimV) and has its eigenvalues set to beta(i)'s. The pair (A44,B42) is controllable. A55 is a stable (nos-noc x nos-noc) matrix. A66 is an unstable (nx-nos x nx-nos) matrix (see st_ility).

Z is a column compression of Sl and k is the normal rank of Sl i.e Sl*Z is a column-compressed linear system. k is the column dimensions of B42,B52,B62 and D2. [B42;B52;B62;D2] is full column rank and has rank k.

This function can be used for the (exact) disturbance decoupling problem. 
DDPS:
   Find u=Fx+Rd which reject Q*d and stabilizes the plant:

    xdot=Ax+Bu+Qd
    y=Cx+Du
DDPS has a solution iff Im(Q) is included in Vg + Im(B). 
Let G=(X(:,dimVg+1:nx))'= left annihilator of Vg i.e. G*Vg=0;
B2=G*B; Q2=G*Q; DDPS solvable iff B2(:,1:k)*R1 + Q2 =0 has a solution.
R=[R1;*] is the solution  (with F=output of abinv).
Im(Q2) is in Im(B2) means row-compression of B2=>row-compression of Q2
Then C*[(sI-A-B*F)^(-1)+D]*(Q+B*R) =0   (<=>G*(Q+B*R)=0)



EXAMPLE :

nu=3;ny=4;nx=7;
nrt=2;ngt=3;ng0=3;nvt=5;rk=2;
flag=list('on',nrt,ngt,ng0,nvt,rk);
Sl=ssrand(ny,nu,nx,flag);alfa=-1;beta=-2;
[X,dims,F,U,k,Z]=abinv(Sl,alfa,beta);
[A,B,C,D]=abcd(Sl);dimV=dims(3);dimR=dims(1);
V=X(:,1:dimV);X2=X(:,dimV+1:nx);
X2'*(A+B*F)*V
(C+D*F)*V
X0=X(:,dimR+1:dimV); spec(X0'*(A+B*F)*X0)
trzeros(Sl)
spec(A+B*F)   //nr=2 evals at -1 and noc-dimV=2 evals at -2.
clean(ss2tf(Sl*Z))
A=diag(1:6);A(2,2)=-7;A(5,5)=-9;B=[1,2;0,3;0,4;0,5;0,0;0,0];C=[zeros(3,3),eye(3,3)];
sl=syslin('c',A,B,C);sl=ss2ss(sl,rand(6,6));
[X,dims,F,U,k,Z]=abinv(sl,alfa,beta);
[A,B,C,D]=abcd(sl);clean(X'*(A+B*F)*X)
clean(X'*B*U)



Author : F.D.


See Also : cainv X, st_ility X, ssrand X, ss2ss X

4.3   arhnk Hankel norm approximant




CALLING SEQUENCE :

[slm]=arhnk(sl,ord,[tol])



PARAMETERS :




DESCRIPTION :

computes slm, the optimal Hankel norm approximant of the stable continuous-time linear system sl with matrices [A,B,C,D].


EXAMPLE :

A=diag([-1,-2,-3,-4,-5]);B=rand(5,1);C=rand(1,5);
sl=syslin('c',A,B,C);
slapprox=arhnk(sl,2);
[nk,W]=hankelsv(sl);nk
[nkred,Wred]=hankelsv(slapprox);nkred



See Also : equil X, equil1 X, hankelsv X

4.4   arl2 SISO model realization by L2 transfer approximation




CALLING SEQUENCE :

h=arl2(y,den0,n [,imp])
h=arl2(y,den0,n [,imp],'all')
[den,num,err]=arl2(y,den0,n [,imp])
[den,num,err]=arl2(y,den0,n [,imp],'all')



PARAMETERS :




DESCRIPTION :

[den,num,err]=arl2(y,den0,n [,imp]) finds a pair of polynomials num and den such that the transfer function num/den is stable and its impulse response approximates (with a minimal l2 norm) the vector y assumed to be completed by an infinite number of zeros.

If y(z) = y(1)(1/z)+y(2)(1/z^2)+ ...+ y(ny)(1/z^ny)

then l2-norm of num/den - y(z) is err.

n is the degree of the polynomial den.

The num/den transfer function is a L2 approximant of the Fourier's series of the rational system.

Various intermediate results are printed according to imp.

[den,num,err]=arl2(y,den0,n [,imp],'all') returns in the vectors of polynomials num and den a set of local optimums for the problem. The solutions are sorted with increasing errors err. In this case den0 is already assumed to be poly(1,'z','c')


EXAMPLE :

v=ones(1,20);
xbasc();
plot2d1('enn',0,[v';zeros(80,1)],2,'051',' ',[1,-0.5,100,1.5])

[d,n,e]=arl2(v,poly(1,'z','c'),1)
plot2d1('enn',0,ldiv(n,d,100),2,'000')
[d,n,e]=arl2(v,d,3)
plot2d1('enn',0,ldiv(n,d,100),3,'000')
[d,n,e]=arl2(v,d,8)
plot2d1('enn',0,ldiv(n,d,100),5,'000')

[d,n,e]=arl2(v,poly(1,'z','c'),4,'all')
plot2d1('enn',0,ldiv(n(1),d(1),100),10,'000')



See Also : ldiv X, imrep2ss X, time_id X, armax X, frep2tf X

4.5   balreal balanced realization




CALLING SEQUENCE :

[slb [,U] ] = balreal(sl)



PARAMETERS :




DESCRIPTION :

Balanced realization of linear system sl=[A,B,C,D]. sl can be a continuous-time or discrete-time state-space system. sl is assumed stable. 
slb=[inv(U)*A*U ,inv(U)*B , C*U , D]
is the balanced realization.

slb is returned as a syslin list.


EXAMPLE :

A=diag([-1,-2,-3,-4,-5]);B=rand(5,2);C=rand(1,5);
sl=syslin('c',A,B,C);
[slb,U]=balreal(sl);
Wc=clean(ctr_gram(slb))
W0=clean(obs_gram(slb))



See Also : ctr_gram X, obs_gram X, hankelsv X, equil X, equil1 X

4.6   bilin general bilinear transform




CALLING SEQUENCE :

[sl1]=bilin(sl,v)



PARAMETERS :




DESCRIPTION :

Given a linear system in state space form, sl=syslin(dom,A,B,C,D) (syslin list), sl1=bilin(sl,v) returns in sl1 a linear system with matrices [A1,B1,C1,D1] such that the transfer function H1(s)=C1*inv(s*eye()-A1)*B1+D1 is obtained from H(z)=C*inv(z*eye()-A)*B+D by replacing z by z=(a*s+b)/(c*s+d). One has w=bilin(bilin(w,[a,b,c,d]),[d,-b,-c,a])


EXAMPLE :

s=poly(0,'s');z=poly(0,'z');
w=ssrand(1,1,3);
wtf=ss2tf(w);v=[2,3,-1,4];a=v(1);b=v(2);c=v(3);d=v(4);
[horner(wtf,(a*z+b)/(c*z+d)),ss2tf(bilin(w,[a,b,c,d]))]
clean(ss2tf(bilin(bilin(w,[a,b,c,d]),[d,-b,-c,a]))-wtf)



See Also : horner X, cls2dls X

4.7   cainv Dual of abinv




CALLING SEQUENCE :

[X,dims,J,Y,k,Z]=cainv(Sl,alfa,beta)



PARAMETERS :




DESCRIPTION :

cainv finds a bases (X,Y) (of state space and output space resp.) and output injection matrix J such that the matrices of Sl in bases (X,Y) are displayed as: 

                  [A11,*,*,*,*,*]                [*]
                  [0,A22,*,*,*,*]                [*]
   X'*(A+J*C)*X = [0,0,A33,*,*,*]   X'*(B+J*D) = [*]
                  [0,0,0,A44,*,*]                [0]
                  [0,0,0,0,A55,*]                [0]
                  [0,0,0,0,0,A66]                [0]

          Y*C*X = [0,0,C13,*,*,*]          Y*D = [*]
                  [0,0,0,0,0,C26]                [0]
The partition of X is defined by the vector dims=[nd1,nu1,dimS,dimSg,dimN] and the partition of Y is determined by k.

Eigenvalues of A11 (nd1 x nd1) are unstable. Eigenvalues of A22 (nu1-nd1 x nu1-nd1) are stable.

The pair (A33, C13) (dimS-nu1 x dimS-nu1, k x dimS-nu1) is observable, and eigenvalues of A33 are set to alfa.

Matrix A44 (dimSg-dimS x dimSg-dimS) is unstable. Matrix A55 (dimN-dimSg,dimN-dimSg) is stable

The pair (A66,C26) (nx-dimN x nx-dimN) is observable, and eigenvalues of A66 set to beta.

The dimS first columns of X span S= smallest (C,A) invariant subspace which contains Im(B), dimSg first columns of X span Sg the maximal "complementary detectability subspace" of Sl

The dimN first columns of X span the maximal "complementary observability subspace" of Sl. (dimS=0 iff B(ker(D))=0).

This function can be used to calculate an unknown input observer: 
// DDEP: dot(x)=A x + Bu + Gd
//           y= Cx   (observation)
//           z= Hx    (z=variable to be estimated, d=disturbance)
//  Find: dot(w) = Fw + Ey + Ru such that
//          zhat = Mw + Ny
//           z-Hx goes to zero at infinity
//  Solution exists iff Ker H contains Sg(A,C,G) inter KerC
//i.e. H is such that:
// For any W which makes a column compression of [Xp(1:dimSg,:);C]
// with Xp=X' and [X,dims,J,Y,k,Z]=cainv(syslin('c',A,G,C));
// [Xp(1:dimSg,:);C]*W = [0 | *] one has
// H*W = [0 | *]  (with at least as many aero columns as above).



See Also : abinv X, dt_ility X

4.8   calfrq frequency response discretization




CALLING SEQUENCE :

[frq,split]=calfrq(h,[fmin,fmax])



PARAMETERS :




DESCRIPTION :

frequency response discretization ; frq is the discretization of [fmin,fmax] such that the peaks in the frequency response are well represented.

Default values for fmin and fmax are 1.d-3, 1.d+3 if h is continuous-time or 1.d-3, 1/(2*h('dt')) if h is discrete-time.

Singularities are located between frq(split(k)) and frq(split(k)+1) for k>1.




EXAMPLE :

s=poly(0,'s')
h=syslin('c',(s^2+2*0.9*10*s+100)/(s^2+2*0.3*10.1*s+102.01))
h1=h*syslin('c',(s^2+2*0.1*15.1*s+228.01)/(s^2+2*0.9*15*s+225)) 
[f1,spl]=calfrq(h1,0.01,1000);
rf=repfreq(h1,f1);
plot2d(real(rf)',imag(rf)')



See Also : bode X, black X, nyquist X, freq X, repfreq X, logspace X

4.9   canon canonical controllable form




CALLING SEQUENCE :

[Ac,Bc,U,ind]=canon(A,B)



PARAMETERS :




DESCRIPTION :

gives the canonical controllable form of the pair (A,B).

Ac=inv(U)*A*U, Bc=inv(U)*B

The vector ind is made of the epsilon_i's indices of the pencil [sI - A , B] (decreasing order). For example with ind=[3,2], Ac and Bc are as follows: 
      [*,*,*,*,*]           [*]
      [1,0,0,0,0]           [0]
Ac=   [0,1,0,0,0]        Bc=[0]
      [*,*,*,*,*]           [*]
      [0,0,0,1,0]           [0]
If (A,B) is controllable, by an appropriate choice of F the * entries of Ac+Bc*F can be arbitrarily set to desired values (pole placement).


EXAMPLE :

A=[1,2,3,4,5;
   1,0,0,0,0;
   0,1,0,0,0;
   6,7,8,9,0;
   0,0,0,1,0];
B=[1,2;
   0,0;
   0,0;
   2,1;
   0,0];
X=rand(5,5);A=X*A*inv(X);B=X*B;    //Controllable pair 
[Ac,Bc,U,ind]=canon(A,B);  //Two indices --> ind=[3.2];
index=1;for k=1:size(ind,'*')-1,index=[index,1+sum(ind(1:k))];end
Acstar=Ac(index,:);Bcstar=Bc(index,:);
s=poly(0,'s');
p1=s^3+2*s^2-5*s+3;p2=(s-5)*(s-3);   
//p1 and p2 are desired closed-loop polynomials with degrees 3,2
c1=coeff(p1);c1=c1($-1:-1:1);c2=coeff(p2);c2=c2($-1:-1:1);
Acstardesired=[-c1,0,0;0,0,0,-c2];  
//Acstardesired(index,:) is companion matrix with char. pol=p1*p2
F=Bcstar\\(Acstardesired-Acstar);   //Feedbak gain
Ac+Bc*F         // Companion form 
spec(A+B*F/U)   // F/U is the gain matrix in original basis.



See Also : obsv_mat X, cont_mat X, ctr_gram X, contrss X, ppol X, contr X, stabil X


Author : F.D.

4.10   cls2dls bilinear transform




CALLING SEQUENCE :

[sl1]=cls2dls(sl,T [,fp])



PARAMETERS :




DESCRIPTION :

given sl=[A,B,C,D] (syslin list),a continuous time system cls2dls returns the sampled system obtained by the bilinear transform s=(2/T)*(z-1)/(z+1).


EXAMPLE :

s=poly(0,'s');z=poly(0,'z');
sl=syslin('c',(s+1)/(s^2-5*s+2));  //Continuous-time system in transfer form
slss=tf2ss(sl);  //Now in state-space form
sl1=cls2dls(slss,0.2);  //sl1= output of cls2dls
sl1t=ss2tf(sl1) // Converts in transfer form
sl2=horner(sl,(2/0.2)*(z-1)/(z+1))   //Compare sl2 and sl1



See Also : horner X

4.11   colregul removing poles and zeros at infinity




CALLING SEQUENCE :

[Stmp,Ws]=colregul(Sl,alfa,beta)



PARAMETERS :




DESCRIPTION :

computes a prefilter Ws such that Stmp=Sl*Ws is proper and with full rank D matrix.

Poles at infinity of Sl are moved to alfa;

Zeros at infinity of Sl are moved to beta;

Sl is a assumed to be a left invertible linear system (syslin list) in state-space representation with possibly a polynomial D matrix.


See Also : invsyslin X, inv X, rowregul X, rowshuff X


Author : F. D. , R. N.



4.12   cont_frm transfer to controllable state-space




CALLING SEQUENCE :

[sl]=cont_frm(NUM,den)



PARAMETERS :




DESCRIPTION :

controllable state-space form of the transfer NUM/den.


EXAMPLE :

s=poly(0,'s');NUM=[1+s,s];den=s^2-5*s+1;
sl=cont_frm(NUM,den); 
slss=ss2tf(sl);       //Compare with NUM/den



See Also : tf2ss X, canon X, contr X

4.13   cont_mat controllability matrix




CALLING SEQUENCE :

Cc=cont_mat(A,B) 
Cc=cont_mat(sl)



PARAMETERS :




DESCRIPTION :

cont_mat returns the controllability matrix of the pair A,B (resp. of the system sl=[A,B,C,D]). 
Cc=[B, AB, A^2 B,..., A^(n-1) B]



See Also : ctr_gram X, contr X, canon X, st_ility X

4.14   contr controllability, controllable subspace




CALLING SEQUENCE :

[n [,U]]=contr(A,B [,tol]) 
[A1,B1,U,ind]=contr(A,B [,tol])



PARAMETERS :




DESCRIPTION :

[n,[U]]=contr(A,B,[tol]) gives the controllable form of an (A,B) pair.(dx/dt = A x + B u or x(n+1) = A x(n) +b u(n)). The n first columns of U make a basis for the controllable subspace.

If V=U(:,1:n), then V'*A*V and V'*B give the controllable part of the (A,B) pair.

[A1,B1,U,ind]=contr(A,B) returns the Hessenberg controllable form of (A,B).


EXAMPLE :

W=ssrand(2,3,5,list('co',3));  //cont. subspace has dim 3.
A=W("A");B=W("B");
[n,U]=contr(A,B);n
A1=U'*A*U;
spec(A1(n+1:$,n+1:$))  //uncontrollable modes
spec(A+B*rand(3,5))



See Also : canon X, cont_mat X, unobs X, stabil X

4.15   contrss controllable part




CALLING SEQUENCE :

[slc]=contrss(sl [,tol])



PARAMETERS :




DESCRIPTION :

returns the controllable part of the linear system sl = (A,B,C,D) in state-space form.


EXAMPLE :

A=[1,1;0,2];B=[1;0];C=[1,1];sl=syslin('c',A,B,C);  //Non minimal
slc=contrss(sl);
sl1=ss2tf(sl);sl2=ss2tf(slc);      //Compare sl1 and sl2



See Also : cont_mat X, ctr_gram X, cont_frm X, contr X

4.16   csim simulation (time response) of linear system




CALLING SEQUENCE :

[y [,x]]=csim(u,t,sl,[x0])



PARAMETERS :




DESCRIPTION :

simulation of the controlled linear system sl. sl is assumed to be a continuous-time system represented by a syslin list.

u is the control and x0 the initial state.

y is the output and x the state.

The control can be:

1. a function : [inputs]=u(t)

2. a list : list(ut,parameter1,....,parametern) such that: inputs=ut(t,parameter1,....,parametern) (ut is a function)

3. the string "impuls" for impulse response calculation (here sl is assumed SISO without direct feed through and x0=0)

4. the string "step" for step response calculation (here sl is assumed SISO without direct feed-through and x0=0)


EXAMPLE :

s=poly(0,'s');rand('seed',0);w=ssrand(1,1,3);w('A')=w('A')-2*eye();
t=0:0.05:5;
//impulse(w) = step (s * w)
xbasc(0);xset("window",0);xselect();
plot2d([t',t'],[(csim('step',t,tf2ss(s)*w))',0*t'])
xbasc(1);xset("window",1);xselect();
plot2d([t',t'],[(csim('impulse',t,w))',0*t'])
//step(w) = impulse (s^-1 * w)
xbasc(3);xset("window",3);xselect();
plot2d([t',t'],[(csim('step',t,w))',0*t'])
xbasc(4);xset("window",4);xselect();
plot2d([t',t'],[(csim('impulse',t,tf2ss(1/s)*w))',0*t'])



See Also : syslin X, dsimul X, flts X, ltitr X, rtitr X, ode X, impl X

4.17   ctr_gram controllability gramian




CALLING SEQUENCE :

[Gc]=ctr_gram(A,B [,dom])
[Gc]=ctr_gram(sl)



PARAMETERS :




DESCRIPTION :

Controllability gramian of (A,B) or sl (a syslin linear system).

dom character string giving the time domain : "d" for a discrete time system and "c" for continuous time (default case). 

                  /+inf                             ---+inf
                  [  At  ' A't                      \\   k      k
             Gc = | e  BB e   dt               Gc =  > A BB' A'
                  ]                                 /
                 /0                                 ---0



EXAMPLE :

A=diag([-1,-2,-3]);B=rand(3,2);
Wc=ctr_gram(A,B)
U=rand(3,3);A1=U*A/U;B1=U*B;
Wc1=ctr_gram(A1,B1)    //Not invariant!



See Also : equil1 X, obs_gram X, contr X, cont_mat X, cont_frm X, contrss X


Author : S. Steer INRIA 1988

4.18   dbphi frequency response to phase and magnitude representation




CALLING SEQUENCE :

[db,phi] =dbphi(repf)



PARAMETERS :




DESCRIPTION :



db(k) is the magnitude of repf(k) expressed in dB i.e. db(k)=20*log(abs(repf(k)))/log(10) and phi(k) is the phase of repf(k) expressed in degrees.


See Also : repfreq X, bode X

4.19   des2tf descriptor to transfer function conversion




CALLING SEQUENCE :

[S]=des2tf(sl)
[Bfs,Bis,chis]=des2tf(sl)



PARAMETERS :




DESCRIPTION :

Given the linear system in descriptor form i.e. Sl=list('des',A,B,C,D,E), des2tf converts sl into its transfer function representation: 
S=C*(s*E-A)^(-1)*B+D
Called with 3 outputs arguments des2tf returns Bfs and Bis two polynomial matrices, and chis polynomial such that: 
S=Bfs/chis - Bis
chis is the determinant of (s*E-A) (up to a xcative constant);


EXAMPLE :

s=poly(0,'s');
G=[1/(s+1),s;1+s^2,3*s^3];
Descrip=tf2des(G);Tf1=des2tf(Descrip)
Descrip2=tf2des(G,"withD");Tf2=des2tf(Descrip2)
[A,B,C,D,E]=Descrip2(2:6);Tf3=C*inv(s*E-A)*B+D



See Also : glever X, pol2des X, tf2des X, ss2tf X, des2ss X, rowshuff X


Author : F. D.

4.20   dscr discretization of linear system




CALLING SEQUENCE :

[sld [,r]]=dscr(sl,dt [,m])



PARAMETERS :




DESCRIPTION :

Discretization of linear system. sl is a continuous-time system:

dx/dt=A*x+B*u (+ noise).

sld is the discrete-time system obtained by sampling sl with the sampling period dt.


EXAMPLE :

s=poly(0,'s');
Sys=syslin('c',[1,1/(s+1);2*s/(s^2+2),1/s])
ss2tf(dscr(tf2ss(Sys),0.1))



See Also : syslin X, flts X, dsimul X

4.21   dsimul state space discrete time simulation




CALLING SEQUENCE :

y=dsimul(sl,u)



PARAMETERS :




DESCRIPTION :

Utility function. If [A,B,C,D]=abcd(sl) and x0=sl('X0'), dsimul returns y=C*ltitr(A,B,u,x0)+D*u i.e. the time response of sl to the input u. sl is assumed to be in state space form (syslin list).


EXAMPLE :

z=poly(0,'z');
h=(1-2*z)/(z^2-0.2*z+1);
sl=tf2ss(h);
u=zeros(1,20);u(1)=1;
x1=dsimul(sl,u)   //Impulse response
u=ones(20,1);
x2=dsimul(sl,u);  //Step response



See Also : syslin X, flts X, ltitr X

4.22   dt_ility detectability test




CALLING SEQUENCE :

[k, [n [,U [,Sld ] ] ]]=dt_ility(Sl [,tol])



PARAMETERS :




DESCRIPTION :

Detectability test for sl, a linear system in state-space representation. U is a basis whose k first columns span the unstable, unobservable subspace of Sl (intersection of unobservable subspace of (A,C) and unstable subspace of A). Detectability means k=0.

Sld = (U'*A*U,U'*B,C*U,D) displays the "detectable part" of Sl=(A,B,C,D), i.e. 
         [*,*,*]
U'*A*U = [0,*,*]
         [0,0,*]
 
C*U = [0,0,*]
with (A33,C3) observable (dimension nx-n), A22 stable (dimension n-k) and A11 unstable (dimension k).


EXAMPLE :

A=[2,1,1;0,-2,1;0,0,3];
C=[0,0,1];
X=rand(3,3);A=inv(X)*A*X;C=C*X;
W=syslin('c',A,[],C);
[k,n,U,W1]=dt_ility(W);
W1("A")
W1("C")



See Also : contr X, st_ility X, unobs X, stabil X

4.23   equil balancing of pair of symmetric matrices




CALLING SEQUENCE :

T=equil(P,Q)



PARAMETERS :




DESCRIPTION :

equil returns t such that:

T*P*T' and inv(T)'*Q*inv(T) are both equal to a same diagonal and positive matrix.


EXAMPLE :

P=rand(4,4);P=P*P';
Q=rand(4,4);Q=Q*Q';
T=equil(P,Q)
clean(T*P*T')
clean(inv(T)'*Q*inv(T))



See Also : equil1 X, balanc X, ctr_gram X

4.24   equil1 balancing (nonnegative) pair of matrices




CALLING SEQUENCE :

[T [,siz]]=equil1(P,Q [,tol])



PARAMETERS :




DESCRIPTION :

equil1 computes t such that:

P1=T*P*T' and Q1=inv(T)'*Q*inv(T) are as follows:

P1 = diag(S1,S2,0,0) and Q1 = diag(S1,0,S3,0) with S1,S2,S3 positive and diagonal matrices with respective dimensions siz=[n1,n2,n3]

tol is a threshold for rank determination in SVD


EXAMPLE :

S1=rand(2,2);S1=S1*S1';
S2=rand(2,2);S2=S2*S2';
S3=rand(2,2);S3=S3*S3';
P=sysdiag(S1,S2,zeros(4,4));
Q=sysdiag(S1,zeros(2,2),S3,zeros(2,2));
X=rand(8,8);
P=X*P*X';Q=inv(X)'*Q*inv(X);
[T,siz]=equil1(P,Q);
P1=clean(T*P*T')
Q1=clean(inv(T)'*Q*inv(T))



See Also : balreal X, minreal X, equil X, hankelsv X


Author : S. Steer 1987

4.25   feedback feedback operation




CALLING SEQUENCE :

Sl=Sl1/.Sl2



PARAMETERS :




DESCRIPTION :

The feedback operation is denoted by /. (slashdot). This command returns Sl=Sl1*(I+Sl2*Sl1)^-1, i.e the (negative) feedback of Sl1 and Sl2. Sl is the transfer v -> y for y = Sl1 u , u = v - Sl2 y.

The result is the same as Sl=LFT([0,I;I,-Sl2],Sl1).

Caution: do not use with decimal point (e.g. 1/.1 is ambiguous!)


EXAMPLE :

S1=ssrand(2,2,3);S2=ssrand(2,2,2);
W=S1/.S2;
ss2tf(S1/.S2)
//Same operation by LFT:
ss2tf(lft([zeros(2,2),eye(2,2);eye(2,2),-S2],S1))
//Other approach: with constant feedback
BigS=sysdiag(S1,S2); F=[zeros(2,2),eye(2,2);-eye(2,2),zeros(2,2)];
Bigclosed=BigS/.F;
W1=Bigclosed(1:2,1:2);   //W1=W (in state-space).
ss2tf(W1)
//Inverting
ss2tf(S1*inv(eye()+S2*S1))



See Also : lft X, sysdiag X, augment X, obscont X

4.26   flts time response (discrete time, sampled system)




CALLING SEQUENCE :

[y [,x]]=flts(u,sl [,x0])
[y]=flts(u,sl [,past])



PARAMETERS :




DESCRIPTION :



State-space form:

sl is a syslin list containing the matrices of the following linear system

sl=syslin('d',A,B,C,D) (see syslin): 
  x[t+1] = A x[t] + B u[t]
  y[t] = C x[t] + D u[t]
or, more generally, if D is a polynomial matrix (p = degree(D(z))) : 
  D(z)=D_0 + z D_1 + z^2 D_2 +..+ z^p D_p
  y[t] = C x[t] + D_0 u[t] + D_1 u[t+1] +..+ D_[p] u[t+p]

  u=[u0,u1,... un]  (input)
  y=[y0,y1,... yn-p]  (output)
  x=x[n-p+1]  (final state, used as \fVx0\fR at next call to flts)
Transfer form:

y=flts(u,sl[,past]). Here sl is a linear system in transfer matrix representation i.e

sl=syslin('d',transfer_matrix) (see syslin). 
    past = [u     ,...,  u   ]
           [ -nd           -1]
           [y     ,...,  y   ]
           [ -nd           -1]
is the matrix of past values of u and y.

nd is the maximum of degrees of lcm's of each row of the denominator matrix of sl. 
u=[u0 u1 ... un]  (input)
y=[y0 y1 ... yn]  (output)
p is the difference between maximum degree of numerator and maximum degree of denominator


EXAMPLE :

sl=syslin('d',1,1,1);u=1:10;
y=flts(u,sl); 
plot2d2("onn",(1:size(u,'c'))',y')
[y1,x1]=flts(u(1:5),sl);y2=flts(u(6:10),sl,x1);
y-[y1,y2]

//With polynomial D:
z=poly(0,'z');
D=1+z+z^2; p =degree(D);
sl=syslin('d',1,1,1,D);
y=flts(u,sl);[y1,x1]=flts(u(1:5),sl);
y2=flts(u(5-p+1:10),sl,x1);  // (update)
y-[y1,y2]

//Delay (transfer form): flts(u,1/z)
// Usual responses
z=poly(0,'z');
h=(1-2*z)/(z^2+0.3*z+1)
u=zeros(1,20);u(1)=1;
imprep=flts(u,tf2ss(h));   //Impulse response
plot2d2("onn",(1:size(u,'c'))',imprep')
u=ones(1,20);
stprep=flts(u,tf2ss(h));   //Step response
plot2d2("onn",(1:size(u,'c'))',stprep')
//
// Other examples
A=[1 2 3;0 2 4;0 0 1];B=[1 0;0 0;0 1];C=eye(3,3);Sys=syslin('d',A,B,C);
H=ss2tf(Sys); u=[1;-1]*(1:10);
//
yh=flts(u,H); ys=flts(u,Sys);
norm(yh-ys,1)    
//hot restart
[ys1,x]=flts(u(:,1:4),Sys);ys2=flts(u(:,5:10),Sys,x);
norm([ys1,ys2]-ys,1)
//
yh1=flts(u(:,1:4),H);yh2=flts(u(:,5:10),H,[u(:,2:4);yh(:,2:4)]);
norm([yh1,yh2]-yh,1)
//with D<>0
D=[-3 8;4 -0.5;2.2 0.9];
Sys=syslin('d',A,B,C,D);
H=ss2tf(Sys); u=[1;-1]*(1:10);
rh=flts(u,H); rs=flts(u,Sys);
norm(rh-rs,1)
//hot restart
[ys1,x]=flts(u(:,1:4),Sys);ys2=flts(u(:,5:10),Sys,x);
norm([ys1,ys2]-rs,1)
//With H:
yh1=flts(u(:,1:4),H);yh2=flts(u(:,5:10),H,[u(:,2:4); yh1(:,2:4)]);
norm([yh1,yh2]-rh)



See Also : ltitr X, dsimul X, rtitr X

4.27   frep2tf transfer function realization from frequency response




CALLING SEQUENCE :

[h [,err]]=frep2tf(frq,repf,dg [,dom,tols,weight])



PARAMETERS :




DESCRIPTION :

Frequency response to transfer function conversion. The order of h is a priori given in dg which must be provided. The following linear system is solved in the least square sense. 
 
        weight(k)*(n( phi_k) - d(phi_k)*rep_k)=0, k=1,..,n
where phi_k= 2*%i*%pi*frq when dom='c' and phi_k=exp(2*%i*%pi*dom*frq if not. If the weight vector is not given a default penalization is used (when dom='c').

A stable and minimum phase system can be obtained by using function factors.


EXAMPLE :

s=poly(0,'s');
h=syslin('c',(s-1)/(s^3+5*s+20))
frq=0:0.05:3;repf=repfreq(h,frq);
clean(frep2tf(frq,repf,3))

Sys=ssrand(1,1,10); 
frq=logspace(-3,2,200);
[frq,rep]=repfreq(Sys,frq);  //Frequency response of Sys
[Sys2,err]=frep2tf(frq,rep,10);Sys2=clean(Sys2)//Sys2 obtained from freq. resp of Sys
[frq,rep2]=repfreq(Sys2,frq); //Frequency response of Sys2
xbasc();bode(frq,[rep;rep2])   //Responses of Sys and Sys2
[sort(trzeros(Sys)),sort(roots(Sys2('num')))]  //zeros
[sort(spec(Sys('A'))),sort(roots(Sys2('den')))] //poles

dom=1/1000; // Sampling time 
z=poly(0,'z');
h=syslin(dom,(z^2+0.5)/(z^3+0.1*z^2-0.5*z+0.08))
frq=(0:0.01:0.5)/dom;repf=repfreq(h,frq);
[Sys2,err]=frep2tf(frq,repf,3,dom);
[frq,rep2]=repfreq(Sys2,frq); //Frequency response of Sys2
xbasc();plot2d1("onn",frq',abs([repf;rep2])');



See Also : imrep2ss X, arl2 X, time_id X, armax X, frfit X

4.28   freq frequency response




CALLING SEQUENCE :

[x]=freq(A,B,C [,D],f)
[x]=freq(NUM,DEN,f)



PARAMETERS :




DESCRIPTION :

x=freq(A,B,C [,D],f) returns a real or complex mxp*t matrix such that:

x(:,k*p:(k+1)*p)= C*inv(f(k)*eye()-A)*B + D.

Thus, for f taking values along the imaginary axis or on the unit circle x is the continuous or discrete time frequency response of (A,B,C,D).

x=freq(NUM,DEN,f) returns a real or complex matrix x such that columns k*(p-1)+1 to k*p of x contain the matrix NUM(f(k))./DEN(f(k))


EXAMPLE :

s=poly(0,'s');
sys=(s+1)/(s^3-5*s+4)
rep=freq(sys("num"),sys("den"),[0,0.9,1.1,2,3,10,20])
[horner(sys,0),horner(sys,20)]
//
Sys=tf2ss(sys);
[A,B,C,D]=abcd(Sys);
freq(A,B,C,[0,0.9,1.1,2,3,10,20])



See Also : repfreq X, horner X

4.29   freson peak frequencies




CALLING SEQUENCE :

fr=freson(h)



PARAMETERS :




DESCRIPTION :

returns the vector of peak frequencies in Hz for the SISO plant h


EXAMPLE :

h=syslin('c',-1+%s,(3+2*%s+%s^2)*(50+0.1*%s+%s^2))
fr=freson(h)
bode(h)
g=20*log(abs(repfreq(h,fr)))/log(10)



See Also : frep2tf X, zgrid X, h_norm X

4.30   g_margin gain margin




CALLING SEQUENCE :

[gm [,fr]]=g_margin(h)



PARAMETERS :




DESCRIPTION :

returns gm, the gain margin in dB of h (SISO plant) and fr, the achieved corresponding frequency in hz. The gain margin is values of the system gain at points where the nyquist plot crosses the negative real axis.


EXAMPLE :

h=syslin('c',-1+%s,3+2*%s+%s^2)
[g,fr]=g_margin(h)
[g,fr]=g_margin(h-10)
nyquist(h-10)



See Also : p_margin X, black X, chart X, nyquist X

4.31   gfrancis Francis equations for tracking




CALLING SEQUENCE :

[L,M,T]=gfrancis(Plant,Model)



PARAMETERS :




DESCRIPTION :

Given the the linear plant: 
    x'= F*x + G*u 
    y = H*x + J*u
and the linear model 
    xm'= A*xm + B*um
    ym = C*xm + D*um
the goal is for the plant to track the model i.e. e = y - ym ---> 0 while keeping stable the state x(t) of the plant. u is given by feedforward and feedback 
 u = L*xm + M*um + K*(x-T*xm) = [K , L-K*T] *(x,xm) + M*um
The matrices T,L,M satisfy generalized Francis equations 
    F*T + G*L = T*A
    H*T + J*L = C
          G*M = T*B
          J*M = D
The matrix K must be chosen as stabilizing the pair (F,G) See example of use in directory demos/tracking.


EXAMPLE :

Plant=ssrand(1,3,5);
[F,G,H,J]=abcd(Plant);
nw=4;nuu=2;A=rand(nw,nw);
st=maxi(real(spec(A)));A=A-st*eye(A);
B=rand(nw,nuu);C=2*rand(1,nw);D=0*rand(C*B);
Model=syslin('c',A,B,C,D);
[L,M,T]=gfrancis(Plant,Model);
norm(F*T+G*L-T*A,1)
norm(H*T+J*L-C,1)
norm(G*M-T*B,1)
norm(J*M-D,1)



See Also : lqg X, ppol X

4.32   imrep2ss state-space realization of an impulse response




CALLING SEQUENCE :

[sl]=imrep2ss(v [,deg])



PARAMETERS :




DESCRIPTION :

Impulse response to linear system conversion (one input). v must have an even number of columns.


EXAMPLE :

s=poly(0,'s');
H=[1/(s+0.5);2/(s-0.4)]   //strictly proper
np=20;w=ldiv(H('num'),H('den'),np);
rep=[w(1:np)';w(np+1:2*np)'];   //The impulse response
H1=ss2tf(imrep2ss(rep))
z=poly(0,'z');
H=(2*z^2-3.4*z+1.5)/(z^2-1.6*z+0.8)     //Proper transfer function
u=zeros(1,20);u(1)=1;
rep=rtitr(H('num'),H('den'),u);   //Impulse rep. 
//   <=> rep=ldiv(H('num'),H('den'),20)
w=z*imrep2ss(rep)   //Realization with shifted impulse response 
// i.e strictly proper to proper
H2=ss2tf(w);



See Also : frep2tf X, arl2 X, time_id X, armax X, markp2ss X, ldiv X

4.33   invsyslin system inversion




CALLING SEQUENCE :

[sl2]=invsyslin(sl1)



PARAMETERS :




DESCRIPTION :

Utility function. Computes the state form of the inverse sl2 of the linear system sl1 (which is also given in state form).

The D-matrix is supposed to be full rank. Old stuff used by inv(S) when S is a syslin list.


See Also : rowregul X, inv X

4.34   kpure continuous SISO system limit feedback gain




CALLING SEQUENCE :

g=kpure(sys)



PARAMETERS :




DESCRIPTION :

kpure(sys) computes the gains g such that the system sys fedback by g (sys/.g) has poles on imaginary axis.


EXAMPLE :

s=poly(0,'s');
h=syslin('c',(s-1)/(1+5*s+s^2+s^3))
xbasc();evans(h)
g=kpure(h)
hf=h/.g(1)
roots(denom(hf))



See Also : evans X, krac2 X

4.35   krac2 continuous SISO system limit feedback gain




CALLING SEQUENCE :

g=krac2(sys)



PARAMETERS :




DESCRIPTION :

krac2(sys) computes the gains g such that the system sys fedback by g (sys/.g) has 2 real equal poles.


EXAMPLE :

h=syslin('c',352*poly(-5,'s')/poly([0,0,2000,200,25,1],'s','c'));
xbasc();evans(h,100)
g=krac2(h)
hf1=h/.g(1);roots(denom(hf1))
hf2=h/.g(2);roots(denom(hf2))



See Also : evans X, kpure X

4.36   lin linearization




CALLING SEQUENCE :

[A,B,C,D]=lin(sim,x0,u0)
[sl]=lin(sim,x0,u0)



PARAMETERS :




DESCRIPTION :

linearization of the non-linear system [y,xdot]=sim(x,u) around x0,u0.

sim is a function which computes y and xdot.

The output is a linear system (syslin list) sl or the four matrices (A,B,C,D) For example, if ftz is the function passed to ode e.g. 
[zd]=ftz(t,z,u)
and if we assume that y=x

[z]=ode(x0,t0,tf,list(ftz,u) compute x(tf).

If simula is the following function: 
deff('[y,xd]=simula(x,u)','xd=ftz(tf,x,u); y=x;');
the tangent linear system sl can be obtained by: 
[A,B,C,D]=lin(simula,z,u)
sl = syslin('c',A,B,C,D,x0)



EXAMPLE :

deff('[y,xdot]=sim(x,u)','xdot=[u*sin(x);-u*x^2];y=xdot(1)+xdot(2)')
sl=lin(sim,1,2);



See Also : external X, derivat X

4.37   lqe linear quadratic estimator (Kalman Filter)




CALLING SEQUENCE :

[K,X]=lqe(P21)



PARAMETERS :




DESCRIPTION :

lqe returns the Kalman gain for the filtering problem in continuous or discrete time.

P21 is a syslin list representing the system P21=[A,B1,C2,D21]

The input to P21 is a white noise with variance: 
     [B1 ]               [Q  S]
BigV=[   ] [ B1' D21'] = [    ]
     [D21]               [S' R]
X is the solution of the stabilizing Riccati equation and A+K*C2 is stable.

In continuous time: 
(A-S*inv(R)*C2)*X+X*(A-S*inv(R)*C2)'-X*C2'*inv(R)*C2*X+Q-S*inv(R)*S'=0

K=-(X*C2'+S)*inv(R)
In discrete time: 
X=A*X*A'-(A*X*C2'+B1*D21')*pinv(C2*X*C2'+D21*D21')*(C2*X*A'+D21*B1')+B1*B1'
K=-(A*X*C2'+B1*D21')*pinv(C2*X*C2'+D21*D21')

xhat(t+1)= E(x(t+1)| y(0),...,y(t)) (one-step predicted x) satisfies the recursion: 
xhat(t+1)=(A+K*C2)*xhat(t) - K*y(t).



See Also : lqr X


Author : F. D.



4.38   lqg LQG compensator




CALLING SEQUENCE :

[K]=lqg(P,r)



PARAMETERS :




DESCRIPTION :

lqg computes the linear optimal LQG (H2) controller for the "augmented" plant P=syslin('c',A,B,C,D) (continuous time) or P=syslin('d',A,B,C,D) (discrete time).

The function lqg2stan returns P and r given the nominal plant, weighting terms and variances of noises.

K is given by the following ABCD matrices: [A+B*Kc+Kf*C+Kf*D*Kc,-Kf,Kc,0] where Kc=lqr(P12) is the controller gain and Kf=lqe(P21) is the filter gain. See example in lqg2stan.


See Also : lqg2stan X, lqr X, lqe X, h_inf X, obscont X


Author : F.D.

4.39   lqg2stan LQG to standard problem




CALLING SEQUENCE :

[P,r]=lqg2stan(P22,bigQ,bigR)



PARAMETERS :




DESCRIPTION :

lqg2stan returns the augmented plant for linear LQG (H2) controller design.

P22=syslin(dom,A,B2,C2) is the nominal plant; it can be in continuous time (dom='c') or discrete time (dom='d'). 
  . 
  x = Ax + w1 + B2u
  y = C2x + w2
for continuous time plant. 
  x[n+1]= Ax[n] + w1 + B2u
      y = C2x + w2
for discrete time plant.

The (instantaneous) cost function is [x' u'] bigQ [x;u].

The covariance of [w1;w2] is E[w1;w2] [w1',w2'] = bigR

If [B1;D21] is a factor of bigQ, [C1,D12] is a factor of bigR and [A,B2,C2,D22] is a realization of P22, then P is a realization of [A,[B1,B2],[C1,-C2],[0,D12;D21,D22]. The (negative) feedback computed by lqg stabilizes P22, i.e. the poles of cl=P22/.K are stable.


EXAMPLE :

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
//Same as Cl=P22/.K; spec(Cl('A'))
s=poly(0,'s')
lqg2stan(1/(s+2),eye(2,2),eye(2,2))



See Also : lqg X, lqr X, lqe X, obscont X, h_inf X, augment X, fstabst X, feedback X


Author : F.D.

4.40   lqr LQ compensator (full state)




CALLING SEQUENCE :

[K,X]=lqr(P12)



PARAMETERS :




DESCRIPTION :

lqr computes the linear optimal LQ full-state gain for the plant P12=[A,B2,C1,D12] in continuous or discrete time.

P12 is a syslin list (e.g. P12=syslin('c',A,B2,C1,D12)).

The cost function is l2-norm of z'*z with z=C1 x + D12 u i.e. [x,u]' * BigQ * [x;u] where 
      [C1' ]               [Q  S]
BigQ= [    ]  * [C1 D12] = [    ]
      [D12']               [S' R]
The gain K is such that A + B2*K is stable.

X is the stabilizing solution of the Riccati equation.

For a continuous plant: 
(A-B2*inv(R)*S')'*X+X*(A-B2*inv(R)*S')-X*B2*inv(R)*B2'*X+Q-S*inv(R)*S'=0

K=-inv(R)*(B2'*X+S)
For a discrete plant: 
X=A'*X*A-(A'*X*B2+C1'*D12)*pinv(B2'*X*B2+D12'*D12)*(B2'*X*A+D12'*C1)+C1'*C1;

K=-pinv(B2'*X*B2+D12'*D12)*(B2'*X*A+D12'*C1)
An equivalent form for X is 
X=Abar'*inv(inv(X)+B2*inv(r)*B2')*Abar+Qbar
with Abar=A-B2*inv(R)*S' and Qbar=Q-S*inv(R)*S'

The 3-blocks matrix pencils associated with these Riccati equations are: 
               discrete                           continuous
   |I   0    0|   | A    0    B2|         |I   0   0|   | A    0    B2|
  z|0   A'   0| - |-Q    I    -S|        s|0   I   0| - |-Q   -A'   -S|
   |0   B2'  0|   | S'   0     R|         |0   0   0|   | S'  -B2'   R|
Caution: It is assumed that matrix R is non singular. In particular, the plant must be tall (number of outputs >= number of inputs).


EXAMPLE :

A=rand(2,2);B=rand(2,1);   //two states, one input
Q=diag([2,5]);R=2;     //Usual notations x'Qx + u'Ru
Big=sysdiag(Q,R);    //Now we calculate C1 and D12
[w,wp]=fullrf(Big);C1=w(:,1:2);D12=w(:,3);   //[C1,D12]'*[C1,D12]=Big
P=syslin('c',A,B,C1,D12);    //The plant (continuous-time)
[K,X]=lqr(P)
spec(A+B*K)    //check stability
norm(A'*X+X*A-X*B*inv(R)*B'*X+Q,1)  //Riccati check
P=syslin('d',A,B,C1,D12);    // Discrete time plant
[K,X]=lqr(P)     
spec(A+B*K)   //check stability
norm(A'*X*A-(A'*X*B)*pinv(B'*X*B+R)*(B'*X*A)+Q-X,1) //Riccati check



See Also : lqe X, gcare X, leqr X


Author : F.D.

4.41   ltitr discrete time response (state space)




CALLING SEQUENCE :

[X]=ltitr(A,B,U,[x0]) 
[xf,X]=ltitr(A,B,U,[x0])



PARAMETERS :




DESCRIPTION :

calculates the time response of the discrete time system 
         x[t+1] = Ax[t] + Bu[t].
The inputs ui's are the columns of the U matrix 
        U=[u0,u1,...,un];
x0 is the vector of initial state (default value : 0) ;

X is the matrix of outputs (same number of columns as U). 

        X=[x0,x1,x2,...,xn]
xf is the vector of final state xf=X[n+1]


EXAMPLE :

A=eye(2,2);B=[1;1];
x0=[-1;-2];
u=[1,2,3,4,5];
x=ltitr(A,B,u,x0)
x1=A*x0+B*u(1)
x2=A*x1+B*u(2)
x3=A*x2+B*u(3) //....



See Also : rtitr X, flts X

4.42   markp2ss Markov parameters to state-space




CALLING SEQUENCE :

[sl]=markp2ss(markpar,n,nout,nin)



PARAMETERS :




DESCRIPTION :

given a set of n Markov parameters stacked in the (row)-matrix markpar of size noutX(n*nin) markp2ss returns a state-space linear system sl (syslin list) such that with [A,B,C,D]=abcd(sl): 
   C*B = markpar(1:nout,1:nin),
   C*A*B =markpar(1:nout,nin+1:2*nin),....



EXAMPLE :

W=ssrand(2,3,4);   //random system with 2 outputs and 3 inputs
[a,b,c,d]=abcd(W);
markpar=[c*b,c*a*b,c*a^2*b,c*a^3*b,c*a^4*b];
S=markp2ss(markpar,5,2,3);
[A,B,C,D]=abcd(S);
Markpar=[C*B,C*A*B,C*A^2*B,C*A^3*B,C*A^4*B];
norm(markpar-Markpar,1)
//Caution... c*a^5*b is not C*A^5*B !



See Also : frep2tf X, tf2ss X, imrep2ss X

4.43   minreal minimal balanced realization




CALLING SEQUENCE :

slb=minreal(sl [,tol])



PARAMETERS :




DESCRIPTION :

[ae,be,ce]=minreal(a,b,c,domain [,tol]) returns the balanced realization of linear system sl (syslin list).

sl is assumed stable.

tol threshold used in equil1.


EXAMPLE :

A=[-eye(2,2),rand(2,2);zeros(2,2),-2*eye(2,2)];
B=[rand(2,2);zeros(2,2)];C=rand(2,4);
sl=syslin('c',A,B,C);
slb=minreal(sl);
ss2tf(sl)
ss2tf(slb)
ctr_gram(sl)
clean(ctr_gram(slb))
clean(obs_gram(slb))



See Also : minss X, balreal X, arhnk X, equil X, equil1 X


Author : S. Steer INRIA 1987

4.44   minss minimal realization




CALLING SEQUENCE :

[slc]=minss( sl [,tol])



PARAMETERS :




DESCRIPTION :

minss returns in slc a minimal realization of sl.


EXAMPLE :

sl=syslin('c',[1 0;0 2],[1;0],[2 1]);
ssprint(sl);
ssprint(minss(sl))



See Also : contr X, minreal X, arhnk X, contrss X, obsvss X, balreal X

4.45   obs_gram observability gramian




CALLING SEQUENCE :

Go=obs_gram(A,C [,dom])
Go=obs_gram(sl)



PARAMETERS :




DESCRIPTION :

Observability gramian of the pair (A,C) or linear system sl (syslin list). dom is the domain which can be 
                  /+inf                             ---+inf
                  [  A't    At                      \\   k      k
             Go = | e  C'C e   dt       or      Go = > A' C'C A
                  ]                                 /
                 /0                                 ---0



EXAMPLE :

A=-diag(1:3);C=rand(2,3);
Go=obs_gram(A,C,'c');     // <=> w=syslin('c',A,[],C); Go=obs_gram(w);
norm(Go*A+A'*Go+C'*C,1)
norm(lyap(A,-C'*C,'c')-Go,1)
A=A/4; Go=obs_gram(A,C,'d');    //discrete time case
norm(lyap(A,-C'*C,'d')-Go,1)



See Also : ctr_gram X, obsvss X, obsv_mat X, lyap X

4.46   obscont observer based controller




CALLING SEQUENCE :

[K]=obscont(P,Kc,Kf)
[J,r]=obscont(P,Kc,Kf)



PARAMETERS :




DESCRIPTION :

obscont returns the observer-based controller associated with a nominal plant P with matrices [A,B,C,D] (syslin list).

The full-state control gain is Kc and the filter gain is Kf. These gains can be computed, for example, by pole placement.

A+B*Kc and A+Kf*C are (usually) assumed stable.

K is a state-space representation of the compensator K: y->u in:

xdot = A x + B u, y=C x + D u, zdot= (A + Kf C)z -Kf y +B u, u=Kc z

K is a linear system (syslin list) with matrices given by: K=[A+B*Kc+Kf*C+Kf*D*Kc,Kf,-Kc].

The closed loop feedback system Cl: v ->y with (negative) feedback K (i.e. y = P u, u = v - K y, or xdot = A x + B u, y = C x + D u, zdot = (A + Kf C) z - Kf y + B u, u = v -F z) is given by Cl = P/.(-K)

The poles of Cl ( spec(cl('A')) ) are located at the eigenvalues of A+B*Kc and A+Kf*C.

Invoked with two output arguments obscont returns a (square) linear system K which parametrizes all the stabilizing feedbacks via a LFT.

Let Q an arbitrary stable linear system of dimension r(2)xr(1) i.e. number of inputs x number of outputs in P. Then any stabilizing controller K for P can be expressed as K=lft(J,r,Q). The controller which corresponds to Q=0 is K=J(1:nu,1:ny) (this K is returned by K=obscont(P,Kc,Kf)). r is size(P) i.e the vector [number of outputs, number of inputs];


EXAMPLE :

ny=2;nu=3;nx=4;P=ssrand(ny,nu,nx);[A,B,C,D]=abcd(P);
Kc=-ppol(A,B,[-1,-1,-1,-1]);  //Controller gain
Kf=-ppol(A',C',[-2,-2,-2,-2]);Kf=Kf';    //Observer gain
cl=P/.(-obscont(P,Kc,Kf));spec(cl('A'))   //closed loop system
[J,r]=obscont(P,Kc,Kf);
Q=ssrand(nu,ny,3);Q('A')=Q('A')-(maxi(real(spec(Q('A'))))+0.5)*eye(Q('A')) 
//Q is a stable parameter
K=lft(J,r,Q);
spec(h_cl(P,K))  // closed-loop A matrix (should be stable);



See Also : ppol X, lqg X, lqr X, lqe X, h_inf X, lft X, syslin X, feedback X, observer X


Author : F.D.



4.47   observer observer design




CALLING SEQUENCE :

Obs=observer(Sys,J)
[Obs,U,m]=observer(Sys [,flag,alfa])



PARAMETERS :




DESCRIPTION :

Obs=observer(Sys,J) returns the observer Obs=syslin(td,A+J*C,[B+J*D,-J],eye(A)) obtained from Sys by a J output injection. (td is the time domain of Sys). More generally, observer returns in Obs an observer for the observable part of linear system Sys: dotx=A x + Bu, y=Cx + Du represented by a syslin list. Sys has nx state variables, nu inputs and ny outputs. Obs is a linear system with matrices [Ao,Bo,Identity], where Ao is no x no, Bo is no x (nu+ny), Co is no x no and no=nx-m.

Input to Obs is [u,y] and output of Obs is:

xhat=estimate of x modulo unobservable subsp. (case flag='pp') or

xhat=estimate of x modulo unstable unobservable subsp. (case flag='st')

case flag='st': z=H*x can be estimated with stable observer iff H*U(:,1:m)=0 and 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 H2*Obs with H2=H*U(:,m+1:nx) i.e. Co, the C-matrix of the observer for H*x, is Co=H2.

In the particular case where the pair (A,C) of Sys is observable, one has m=0 and the linear system U*Obs (resp. H*U*Obs) is an observer for x (resp. Hx). The error spectrum is alpha(1),alpha(2),...,alpha(nx).


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);Sys3=ss2tf(H*U(:,m+1:$)*Obs*[Idu;Sys])  
//Transfer u-->[u;y=Sys*u]-->Obs-->xhat-->HUxhat=zhat  i.e. u-->output of Obs
//this transfer must equal Sys2, the u-->z transfer  (H2=eye).



See Also : dt_ility X, unobs X, stabil X


Author : F.D.

4.48   obsv_mat observability matrix




CALLING SEQUENCE :

[O]=obsv_mat(A,C)
[O]=obsv_mat(sl)



PARAMETERS :




DESCRIPTION :

obsv_mat returns the observability matrix: 
O=[C; CA; CA^2;...; CA^(n-1) ]



See Also : contrss X, obsvss X, obs_gram X

4.49   obsvss observable part




CALLING SEQUENCE :

[Ao,Bo,Co]=obsvss(A,B,C [,tol])
[slo]=obsvss(sl [,tol])



PARAMETERS :




DESCRIPTION :

slo=(Ao,Bo,Co) is the observable part of linear system sl=(A,B,C) (syslin list)

tol threshold to test controllability (see contr); default value = 100*%eps


See Also : contr X, contrss X, obsv_mat X, obs_gram X

4.50   p_margin phase margin




CALLING SEQUENCE :

[phm,fr]=p_margin(h)
phm=p_margin(h)



PARAMETERS :




DESCRIPTION :

The phase margin is the values of the phase at points where the nyquist plot of h crosses the unit circle.


EXAMPLE :

h=syslin('c',-1+%s,3+2*%s+%s^2)
[p,fr]=p_margin(h)  
[p,fr]=p_margin(h+0.7)  
nyquist(h+0.7)
t=(0:0.1:2*%pi)';plot2d(sin(t),cos(t),-3,'000')



See Also : chart X, black X, g_margin X, nyquist X


Author : S. S.

4.51   pfss partial fraction decomposition




CALLING SEQUENCE :

elts=pfss(Sl)
elts=pfss(Sl,rmax)
elts=pfss(Sl,'cord')
elts=pfss(Sl,rmax,'cord')



PARAMETERS :

Sl : syslin list (state-space or transfer linear system) rmax : real number controlling the conditioning of block diagoanalization cord : character string 'c' or 'd'.


DESCRIPTION :

Partial fraction decomposition of the linear system Sl (in state-space form, transfer matrices are automatically converted to state-space form by tf2ss):

elts is the list of linear systems which add up to Sl i.e. elts=list(S1,S2,S3,...,Sn) with:

Sl = S1 + S2 +... +Sn.

Each Si contains some poles of S according to the block-diagonalization of the A matrix of S.

For non proper systems the polynomial part of Sl is put in the last entry of elts.

If Sl is given in transfer form, it is first converted into state-space and each subsystem Si is then converted in transfer form.

The A matrix is of the state-space is put into block diagonal form by function bdiag. The optional parameter rmax is sent to bdiag. If rmax should be set to a large number to enforce block-diagonalization.

If the optional flag cord='c' is given the elements in elts are sorted according to the real part (resp. magnitude if cord='d') of the eigenvalues of A matrices.


EXAMPLE :

W=ssrand(1,1,6);
elts=pfss(W); 
W1=0;for k=1:size(elts), W1=W1+ss2tf(elts(k));end
clean(ss2tf(W)-W1)



See Also : pbig X, bdiag X, coffg X, dtsi X


Author : F.D.

4.52   phasemag phase and magnitude computation




CALLING SEQUENCE :

[phi,db]=phasemag(z [,mod])



PARAMETERS :




DESCRIPTION :

phasemag computes the phases and magnitudes of the entries of a complex matrix. For mod='c' phasemag computes phi(:,i+1) to minimize the distance with phi(:,i), i.e. it tries to obtain a "continuous representation" of the phase.

To obtain the phase between -%pi and %pi use phi=atan(imag(z),real(z))


EXAMPLE :

s=poly(0,'s');
h=syslin('c',1/((s+5)*(s+10)*(100+6*s+s*s)*(s+.3)));
[frq,rf]=repfreq(h,0.1,20,0.005);
xbasc(0);
plot2d(frq',phasemag(rf,'c')');
xbasc(1);
plot2d(frq',phasemag(rf,'m')');



See Also : repfreq X, gainplot X, atan X, bode X

4.53   ppol pole placement




CALLING SEQUENCE :

[K]=ppol(A,B,poles)



PARAMETERS :




DESCRIPTION :

K=ppol(A,B,poles) returns a mxn gain matrix K such that the eigenvalues of A-B*K are poles. The pair (A,B) must be controllable. Complex number in poles must appear in conjugate pairs.

An output-injection gain F for (A,C) is obtained as follows:

Ft=ppol(A',C',poles); F=Ft'

The algorithm is by P.H. Petkov.


EXAMPLE :

A=rand(3,3);B=rand(3,2);
F=ppol(A,B,[-1,-2,-3]);
spec(A-B*F)



See Also : canon X, stabil X

4.54   projsl linear system projection




CALLING SEQUENCE :

[slp]=projsl(sl,Q,M)



PARAMETERS :




DESCRIPTION :

slp= projected model of sl where Q*M is the full rank factorization of the projection.

If (A,B,C,D) is the representation of sl, the projected model is given by (M*A*Q,M*B,C*Q,D).

Usually, the projection Q*M is obtained as the spectral projection of an appropriate auxiliary matrix W e.g. W = product of (weighted) gramians or product of Riccati equations.


EXAMPLE :

rand('seed',0);sl=ssrand(2,2,5);[A,B,C,D]=abcd(sl);poles=spec(A)
[Q,M]=pbig(A,0,'c');  //keeping unstable poles
slred=projsl(sl,Q,M);spec(slred('A'))
sl('D')=rand(2,2);  //making proper system
trzeros(sl)  //zeros of sl
wi=inv(sl);  //wi=inverse in state-space
[q,m]=psmall(wi('A'),2,'d');  //keeping small zeros (poles of wi) i.e. abs(z)<2
slred2=projsl(sl,q,m);
trzeros(slred2)  //zeros of slred2 = small zeros of sl
//  Example keeping second order modes
A=diag([-1,-2,-3]);
sl=syslin('c',A,rand(3,2),rand(2,3));[nk2,W]=hankelsv(sl)
[Q,M]=pbig(W,nk2(2)-%eps,'c');    //keeping 2 eigenvalues of W
slr=projsl(sl,Q,M);  //reduced model
hankelsv(slr)



See Also : pbig X


Author : F. D.

4.55   repfreq frequency response




CALLING SEQUENCE :

[ [frq,] repf]=repfreq(sys,fmin,fmax [,step])
[ [frq,] repf]=repfreq(sys [,frq])
[ frq,repf,splitf]=repfreq(sys,fmin,fmax [,step])
[ frq,repf,splitf]=repfreq(sys [,frq])



PARAMETERS :




DESCRIPTION :

repfreq returns the frequency response calculation of a linear system. If sys(s) is the transfer function of Sys, repf(k) equals sys(s) evaluated at s= %i*frq(k)*2*%pi for continuous time systems and at exp(2*%i*%pi*dt*frq(k)) for discrete time systems (dt is the sampling period).

db(k) is the magnitude of repf(k) expressed in dB i.e. db(k)=20*log10(abs(repf(k))) and phi(k) is the phase of repf(k) expressed in degrees.

If fmin,fmax,step are input parameters, the response is calculated for the vector of frequencies frq given by: frq=[10.^((log10(fmin)):step:(log10(fmax))) fmax];

If step is not given, the output parameter frq is calculated by frq=calfrq(sys,fmin,fmax).

Vector frq is splitted into regular parts with the split vector. frq(splitf(k):splitf(k+1)-1) has no critical frequency. sys has a pole in the range [frq(splitf(k)),frq(splitf(k)+1)] and no poles outside.


EXAMPLE :

A=diag([-1,-2]);B=[1;1];C=[1,1];
Sys=syslin('c',A,B,C);
frq=0:0.02:5;w=frq*2*%pi; //frq=frequencies in Hz ;w=frequencies in rad/sec;
[frq1,rep] =repfreq(Sys,frq);
[db,phi]=dbphi(rep);
Systf=ss2tf(Sys)    //Transfer function of Sys
x=horner(Systf,w(2)*sqrt(-1))    // x is Systf(s) evaluated at s = i w(2)
rep=20*log(abs(x))/log(10)   //magnitude of x in dB
db(2)    // same as rep
ang=atan(imag(x),real(x));   //in rad.
ang=ang*180/%pi              //in degrees
phi(2)
repf=repfreq(Sys,frq);
repf(2)-x



See Also : bode X, freq X, calfrq X, horner X, nyquist X, dbphi X


Author : S. S.

4.56   ricc Riccati equation




CALLING SEQUENCE :

[X]=ricc(A,B,C,"cont")
[X]=ricc(F,G,H,"disc")



PARAMETERS :




DESCRIPTION :

Riccati solver.

Continuous time: 
  X=ricc(A,B,C,'cont')
gives a solution to the continuous time ARE 
  A'*X+X*A-X*B*X+C=0 .
B and C are assumed to be nonnegative definite. (A,G) is assumed to be stabilizable with G*G' a full rank factorization of B.

(A,H) is assumed to be detectable with H*H' a full rank factorization of C.

Discrete time: 
  X=ricc(F,G,H,'disc')
gives a solution to the discrete time ARE 
  X=F'*X*F-F'*X*G1*((G2+G1'*X*G1)^-1)*G1'*X*F+H
F is assumed invertible and G = G1*inv(G2)*G1'.

One assumes (F,G1) stabilizable and (C,F) detectable with C'*C full rank factorization of H. Use preferably ric_desc.


EXAMPLE :

//Standard formulas to compute Riccati solutions
A=rand(3,3);B=rand(3,2);C=rand(3,3);C=C*C';R=rand(2,2);R=R*R'+eye();
B=B*inv(R)*B';
X=ricc(A,B,C,'cont');
norm(A'*X+X*A-X*B*X+C,1)
H=[A -B;-C -A'];
[T,d]=gschur(eye(H),H,'cont');T=T(:,1:d);
X1=T(4:6,:)/T(1:3,:);
norm(X1-X,1)
[T,d]=schur(H,'cont');T=T(:,1:d);
X2=T(4:6,:)/T(1:3,:);
norm(X2-X,1)
//       Discrete time case
F=A;B=rand(3,2);G1=B;G2=R;G=G1/G2*G1';H=C;
X=ricc(F,G,H,'disc');
norm(F'*X*F-(F'*X*G1/(G2+G1'*X*G1))*(G1'*X*F)+H-X)
H1=[eye(3,3) G;zeros(3,3) F'];
H2=[F zeros(3,3);-H eye(3,3)];
[T,d]=gschur(H2,H1,'disc');T=T(:,1:d);X1=T(4:6,:)/T(1:3,:);
norm(X1-X,1)
Fi=inv(F);
Hami=[Fi Fi*G;H*Fi F'+H*Fi*G];
[T,d]=schur(Hami,'d');T=T(:,1:d);
Fit=inv(F');
Ham=[F+G*Fit*H -G*Fit;-Fit*H Fit];
[T,d]=schur(Ham,'d');T=T(:,1:d);X2=T(4:6,:)/T(1:3,:);
norm(X2-X,1)



See Also : riccati X, ric_desc X, schur X, gschur X

4.57   rowregul removing poles and zeros at infinity




CALLING SEQUENCE :

[Stmp,Ws]=rowregul(Sl,alfa,beta)



PARAMETERS :




DESCRIPTION :

computes a postfilter Ws such that Stmp=Ws*Sl is proper and with full rank D matrix.

Poles at infinity of Sl are moved to alfa;

Zeros at infinity of Sl are moved to beta;

Sl is a assumed to be a right invertible linear system (syslin list) in state-space representation with possibly a polynomial D matrix.

This function is the dual of colregul (see function code).


EXAMPLE :

s=%s;
w=[1/s,0;s/(s^3+2),2/s];
Sl=tf2ss(w);
[Stmp,Ws]=rowregul(Sl,-1,-2);
Stmp('D')     // D matrix of Stmp
clean(ss2tf(Stmp))



See Also : invsyslin X, colregul X


Author : F. D. , R. N.



4.58   rtitr discrete time response (transfer matrix)




CALLING SEQUENCE :

[y]=rtitr(Num,Den,u [,up,yp])



PARAMETERS :




DESCRIPTION :

y=rtitr(Num,Den,u [,up,yp]) returns the time response of the discrete time linear system with transfer matrix Den^-1 Num for the input u, i.e y and u are such that Den y = Num u at t=0,1,...

If d1=maxi(degree(Den)), and d2=maxi(degree(Num)) the polynomial matrices Den(z) and Num(z) may be written respectively as: 
  D(z)= D_0  + D_1  z + ... + D_d1   z^d1
  N(z)= N_0  + N_1  z + ... + N_d2   z^d2
and Den y = Num u is interpreted as the recursion: 
  D(0)y(t)+D(1)y(t+1)+...+ D(d1)y(t+d1)= N(0) u(t) +....+ N(d2) u(t+d2)
It is assumed that D(d1) is non singular.

The columns of u are the inputs of the system at t=0,1,...,T: 
  u=[u(0) , u(1),...,u(T)]
The outputs at t=0,1,...,T+d1-d2 are the columns of the matrix y: 
  y=[y(0), y(1),  .... y(T+d1-d2)]
up and yp define the initial conditions for t < 0 i.e 
  up=[u(-d1), ..., u(-1)  ]
  yp=[y(-d1), ...  y(-1)  ]
Depending on the relative values of d1 and d2, some of the leftmost components of up, yp are ignored. The default values of up and yp are zero: up = 0*ones(m,d1), yp=0*ones(n,d1)


EXAMPLE :

z=poly(0,'z');
Num=1+z;Den=1+z;u=[1,2,3,4,5];
rtitr(Num,Den,u)-u
//Other examples
//siso
//causal
n1=1;d1=poly([1 1],'z','coeff');       // y(j)=-y(j-1)+u(j-1)
r1=[0 1 0 1 0 1 0 1 0 1 0];
r=rtitr(n1,d1,ones(1,10));norm(r1-r,1)
//hot restart
r=rtitr(n1,d1,ones(1,9),1,0);norm(r1(2:11)-r)
//non causal
n2=poly([1 1 1],'z','coeff');d2=d1;    // y(j)=-y(j-1)+u(j-1)+u(j)+u(j+1)
r2=[2 1 2 1 2 1 2 1 2];
r=rtitr(n2,d2,ones(1,10));norm(r-r2,1)
//hot restart
r=rtitr(n2,d2,ones(1,9),1,2);norm(r2(2:9)-r,1)
//
//MIMO example
//causal
d1=d1*diag([1 0.5]);n1=[1 3 1;2 4 1];r1=[5;14]*r1;
r=rtitr(n1,d1,ones(3,10));norm(r1-r,1)
//
r=rtitr(n1,d1,ones(3,9),[1;1;1],[0;0]);
norm(r1(:,2:11)-r,1)
//polynomial n1  (same ex.)
n1(1,1)=poly(1,'z','c');r=rtitr(n1,d1,ones(3,10));norm(r1-r,1)
//
r=rtitr(n1,d1,ones(3,9),[1;1;1],[0;0]);
norm(r1(:,2:11)-r,1)
//non causal
d2=d1;n2=n2*n1;r2=[5;14]*r2;
r=rtitr(n2,d2,ones(3,10));norm(r2-r)
//
r=rtitr(n2,d2,ones(3,9),[1;1;1],[10;28]);
norm(r2(:,2:9)-r,1)
//
//  State-space or transfer
a = [0.21 , 0.63 , 0.56 , 0.23 , 0.31
     0.76 , 0.85 , 0.66 , 0.23 , 0.93
     0 , 0.69 , 0.73 , 0.22 , 0.21
     0.33 , 0.88 , 0.2 , 0.88 , 0.31
     0.67 , 0.07 , 0.54 , 0.65 , 0.36];
b = [0.29 , 0.5 , 0.92
     0.57 , 0.44 , 0.04
     0.48 , 0.27 , 0.48
     0.33 , 0.63 , 0.26
     0.59 , 0.41 , 0.41];
c = [0.28 , 0.78 , 0.11 , 0.15 , 0.84
     0.13 , 0.21 , 0.69 , 0.7 , 0.41];
d = [0.41 , 0.11 , 0.56
     0.88 , 0.2 , 0.59];
s=syslin('d',a,b,c,d);
h=ss2tf(s);num=h('num');den=h('den');den=den(1,1)*eye(2,2);
u=1;u(3,10)=0;r3=flts(u,s);
r=rtitr(num,den,u);norm(r3-r,1)



See Also : ltitr X, exp X, flts X

4.59   sm2des system matrix to descriptor




CALLING SEQUENCE :

[Des]=sm2des(Sm);



PARAMETERS :




DESCRIPTION :

Utility function: converts the system matrix: 
 Sm = [-sE + A   B;
      [    C     D]
to descriptor system Des=list('des',A,B,C,D,E)).


See Also : ss2des X, sm2ss X

4.60   sm2ss system matrix to state-space




CALLING SEQUENCE :

[Sl]=sm2ss(Sm);



PARAMETERS :




DESCRIPTION :

Utility function: converts the system matrix: 
 Sm = [-sI + A   B;
      [    C     D]
to linear system in state-space representation (syslin) list.


See Also : ss2des X

4.61   specfact spectral factor




CALLING SEQUENCE :

[W0,L]=specfact(A,B,C,D)



DESCRIPTION :

Given a spectral density matrix phi(s): 
                 -1                   -1
     R + C*(s*I-A) * B  +  B'*(-s*I-A') * C'  with R=D+D' > 0
specfact computes W0 and L such that W(s)=W0+L*(s*I-A)^-1*B is a spectral factor of of PHI(s), i.e.

phi(s)=W'(-s)*W(s)


EXAMPLE :

A=diag([-1,-2]);B=[1;1];C=[1,1];D=1;s=poly(0,'s');
W1=syslin('c',A,B,C,D);
phi=gtild(W1,'c')+W1;
phis=clean(ss2tf(phi))
clean(phis-horner(phis,-s)');   //check this is 0...
[A,B,C,D]=abcd(W1);
[W0,L]=specfact(A,B,C,D);
W=syslin('c',A,B,L,W0)
Ws=ss2tf(W);
horner(Ws,-s)*Ws



See Also : gtild X, sfact X, fspecg X


Author : F. D.



4.62   ss2des (polynomial) state-space to descriptor form




CALLING SEQUENCE :

S=ss2des(Sl)
S=ss2des(Sl,flag)



PARAMETERS :




DESCRIPTION :

Given the linear system in state-space representation Sl (syslin list), with a D matrix which is either polynomial or constant, but not zero ss2des returns a descriptor system as list('des',A,B,C,0,E) such that: 
Sl=C*(s*E-A)^(-1)*B
If the flag "withD" is given, S=list('des',A,B,C,D,E) with a D matrix of maximal rank.


EXAMPLE :

s=poly(0,'s');
G=[1/(s+1),s;1+s^2,3*s^3];Sl=tf2ss(G);
S=ss2des(Sl)
S1=ss2des(Sl,"withD")
Des=des2ss(S);Des(5)=clean(Des(5))
Des1=des2ss(S1)



See Also : pol2des X, tf2des X, des2ss X


Author : F. D.

4.63   ss2ss state-space to state-space conversion, feedback, injection




CALLING SEQUENCE :

[Sl1,right,left]=ss2ss(Sl,T, [F, [G , [flag]]])



PARAMETERS :




DESCRIPTION :

Returns the linear system Sl1=[A1,B1,C1,D1] where A1=inv(T)*A*T, B1=inv(T)*B, C1=C*T, D1=D.

Optional parameters F and G are state feedback and output injection respectively.

For example, Sl1=ss2ss(Sl,T,F) returns Sl1 with: 

 
              | inv(T)*(A+B*F)*T  ,    inv(T)*B |
        Sl1 = |                                 |
              |     (C+D*F)*T     ,      D      |
and right is a non singular linear system such that Sl1=Sl*right.

Sl1*inv(right) is a factorization of Sl.

Sl1=ss2ss(Sl,T,0*F,G) returns Sl1 with: 
 
              | inv(T)*(A+G*C)*T  , inv(T)*(B+G*D) |
        Sl1 = |                                    |
              |      C*T          ,      D         |
and left is a non singular linear system such that Sl1=left*Sl (right=Id if F=0).

When both F and G are given, Sl1=left*Sl*right.


EXAMPLE :

Sl=ssrand(2,2,5); trzeros(Sl)       // zeros are invariant:
Sl1=ss2ss(Sl,rand(5,5),rand(2,5),rand(5,2)); 
trzeros(Sl1), trzeros(rand(2,2)*Sl1*rand(2,2))
// output injection [ A + GC, (B+GD,-G)]
//                  [   C   , (D   , 0)]
p=1,m=2,n=2; sys=ssrand(p,m,n);

// feedback (m,n)  first and then output injection.

F1=rand(m,n);
G=rand(n,p);
[sys1,right,left]=ss2ss(sys,rand(n,n),F1,G,2);

// Sl1 equiv left*sysdiag(sys*right,eye(p,p)))

res=clean(ss2tf(sys1) - ss2tf(left*sysdiag(sys*right,eye(p,p))))

// output injection then feedback (m+p,n) 
F2=rand(p,n); F=[F1;F2];
[sys2,right,left]=ss2ss(sys,rand(n,n),F,G,1);

// Sl1 equiv left*sysdiag(sys,eye(p,p))*right 

res=clean(ss2tf(sys2)-ss2tf(left*sysdiag(sys,eye(p,p))*right))

// when F2= 0; sys1 and sys2 are the same 
F2=0*rand(p,n);F=[F1;F2];
[sys2,right,left]=ss2ss(sys,rand(n,n),F,G,1);

res=clean(ss2tf(sys2)-ss2tf(sys1))



See Also : projsl X, feedback X

4.64   ss2tf conversion from state-space to transfer function




CALLING SEQUENCE :

[h]=ss2tf(sl) 
[Ds,NUM,chi]=ss2tf(sl)



PARAMETERS :




DESCRIPTION :

Called with three outputs [Ds,NUM,chi]=ss2tf(sl) returns the numerator polynomial matrix NUM, the characteristic polynomial chi and the polynomial part Ds separately i.e.: 
h=NUM/chi + Ds
Method:

One uses the characteristic polynomial and det(A+Eij)=det(A)+C(i,j) where C is the adjugate matrix of A.


EXAMPLE :

s=poly(0,'s');
h=[1,1/s;1/(s^2+1),s/(s^2-2)]
sl=tf2ss(h);
h=clean(ss2tf(sl))
[Ds,NUM,chi]=ss2tf(sl)



See Also : tf2ss X, syslin X, nlev X, glever X

4.65   st_ility stabilizability test




CALLING SEQUENCE :

[ns, [nc, [,U [,Slo] ]]]=st_ility(Sl [,tol])



PARAMETERS :




DESCRIPTION :

Slo=( U'*A*U, U'*B, C*U, D, U'*x0 ) (syslin list) displays the stabilizable form of Sl. Stabilizability means ns=nx (dim. of A matrix). 
         [*,*,*]            [*]
U'*A*U = [0,*,*]     U'*B = [0]
         [0,0,*]            [0]
where (A11,B1) (dim(A11)= nc) is controllable and A22 (dim(A22)=ns-nc) is stable. "Stable" means real part of eigenvalues negative for a continuous linear system, and magnitude of eigenvalues lower than one for a discrete-time system (as defined by syslin).


EXAMPLE :

A=diag([0.9,-2,3]);B=[0;0;1];Sl=syslin('c',A,B,[]);
[ns,nc,U]=st_ility(Sl);
U'*A*U
U'*B
[ns,nc,U]=st_ility(syslin('d',A,B,[]));
U'*A*U
U'*B



See Also : dt_ility X, contr X, stabil X, ssrand X


Author : S. Steer INRIA 1988

4.66   stabil stabilization




CALLING SEQUENCE :

F=stabil(A,B,alfa)
K=stabil(Sys,alfa,beta)



PARAMETERS :




DESCRIPTION :

F=stabil(A,B,alfa) returns a gain matrix F such that A+B*F is stable if pair (A,B) is stabilizable. Assignable poles are set to alfa(1),alfa(2),.... If (A,B) is not stabilizable a warning is given and assignable poles are set to alfa(1),alfa(2),.... If alfa is a number all eigenvalues are set to this alfa (default value is alfa=-1).

K=stabil(Sys,alfa,beta) returns K, a compensator for Sys such that (A,B)-controllable eigenvalues are set to alfa and (C,A)-observable eigenvalues are set to beta.

All assignable closed loop poles (which are given by the eigenvalues of Aclosed=h_cl(Sys,K) are set to alfa(i)'s and beta(j)'s.


EXAMPLE :

// Gain:
Sys=ssrand(0,2,5,list('st',2,3,3));
A=Sys('A');B=Sys('B');F=stabil(A,B);
spec(A) //2 controllable modes 2 unstable uncontrollable modes
//and one stable uncontrollable mode
spec(A+B*F) //the two controllable modes are set to -1.
// Compensator:
Sys=ssrand(3,2,5,list('st',2,3,3)); //3 outputs, 2 inputs, 5 states
//2 controllables modes, 3 controllable or stabilizable modes.
K=stabil(Sys,-2,-3);  //Compensator for Sys.
spec(Sys('A'))
spec(h_cl(Sys,K))   //K Stabilizes what can be stabilized.



See Also : st_ility X, contr X, ppol X

4.67   svplot singular-value sigma-plot




CALLING SEQUENCE :

[SVM]=svplot(sl,[w])



PARAMETERS :




DESCRIPTION :

computes for the system sl=(A,B,C,D) the singular values of its transfer function matrix: 
             G(jw) = C(jw*I-A)B^-1+D
or
             G(exp(jw)) = C(exp(jw)*I-A)B^-1+D
or
             G(exp(jwT)) = C(exp(jw*T)*I-A)B^-1+D
evaluated over the frequency range specified by w. (T is the sampling period, T=sl('dt') for sampled systems).

sl is a syslin list representing the system [A,B,C,D] in state-space form. sl can be continous or discrete time or sampled system.

The i-th column of the output matrix SVM contains the singular values of G for the i-th frequency value w(i). 
SVM = svplot(sl)
is equivalent to 
SVM = svplot(sl,logspace(-3,3))  (continuous)

SVM = svplot(sl,logspace(-3,%pi)) (discrete)



EXAMPLE :

x=logspace(-3,3);
y=svplot(ssrand(2,2,4));
xbasc();plot2d1("oln",x',20*log(y')/log(10));
xgrid(12)
xtitle("Singular values plot","(Rd/sec)", "Db");



Author : F.D

4.68   syssize size of state-space system




CALLING SEQUENCE :

[r,nx]=syssize(Sl)



PARAMETERS :




DESCRIPTION :

returns in r the vector [number of outputs, number of inputs] of the linear system Sl. nx is the number of states of Sl.


See Also : size X

4.69   tf2ss transfer to state-space




CALLING SEQUENCE :

sl=tf2ss(h [,tol])



PARAMETERS :




DESCRIPTION :

transfer to state-space conversion:

h=C*(s*eye()-A)^-1*B+D(s)


EXAMPLE :

s=poly(0,'s');
H=[2/s,(s+1)/(s^2-5)];
Sys=tf2ss(H)
clean(ss2tf(Sys))



See Also : ss2tf X, tf2des X, des2tf X

4.70   time_id SISO least square identification




CALLING SEQUENCE :

[H [,err]]=time_id(n,u,y)



PARAMETERS :




DESCRIPTION :

Identification of discrete time response. If y is strictly proper (y(1)=0) then time_id computes the least square solution of the linear equation: Den*y-Num*u=0 with the constraint coeff(Den,n):=1. if y(1)~=0 then the algorithm first computes the proper part solution and then add y(1) to the solution


EXAMPLE :

z=poly(0,'z');
h=(1-2*z)/(z^2-0.5*z+5)
rep=[0;ldiv(h('num'),h('den'),20)]; //impulse response
H=time_id(2,'impuls',rep)
//  Same example with flts and u
u=zeros(1,20);u(1)=1;
rep=flts(u,tf2ss(h));        //impulse response
H=time_id(2,u,rep)
//  step response
u=ones(1,20);
rep=flts(u,tf2ss(h));     //step response.
H=time_id(2,'step',rep)
H=time_id(3,u,rep)    //with u as input and too high order required



Author : Serge Steer INRIA


See Also : imrep2ss X, arl2 X, armax X, frep2tf X

4.71   trzeros transmission zeros and normal rank




CALLING SEQUENCE :

[tr]=trzeros(Sl)
[nt,dt,rk]=trzeros(Sl)



PARAMETERS :




DESCRIPTION :

Called with one output argument, trzeros(Sl) returns the transmission zeros of the linear system Sl.

Sl may have a polynomial (but square) D matrix.

Called with 2 output arguments, trzeros returns the transmission zeros of the linear system Sl as tr=nt./dt;

(Note that some components of dt may be zeros)

Called with 3 output arguments, rk is the normal rank of Sl

Transfer matrices are converted to state-space.

If Sl is a (square) polynomial matrix trzeros returns the roots of its determinant.

For usual state-space system trzeros uses the state-space algorithm of Emami-Naeni & Van Dooren.

If D is invertible the transmission zeros are the eigenvalues of the "A matrix" of the inverse system : A - B*inv(D)*C;

If C*B is invertible the transmission zeros are the eigenvalues of N*A*M where M*N is a full rank factorization of eye(A)-B*inv(C*B)*C;

For systems with a polynomial D matrix zeros are calculated as the roots of the determinant of the system matrix.

Caution: the computed zeros are not always reliable, in particular in case of repeated zeros.


EXAMPLE :

W1=ssrand(2,2,5);trzeros(W1)    //call trzeros
roots(det(systmat(W1)))         //roots of det(system matrix)
s=poly(0,'s');W=[1/(s+1);1/(s-2)];W2=(s-3)*W*W';[nt,dt,rk]=trzeros(W2);
St=systmat(tf2ss(W2));[Q,Z,Qd,Zd,numbeps,numbeta]=kroneck(St);
St1=Q*St*Z;rowf=(Qd(1)+Qd(2)+1):(Qd(1)+Qd(2)+Qd(3));
colf=(Zd(1)+Zd(2)+1):(Zd(1)+Zd(2)+Zd(3));
roots(St1(rowf,colf)), nt./dt     //By Kronecker form



See Also : gspec X, kroneck X

4.72   unobs unobservable subspace




CALLING SEQUENCE :

[n,[U]]=unobs(A,C,[tol])



PARAMETERS :




DESCRIPTION :

[n,[U]]=unobs(A,C,[tol]) gives the unobservable form of an (A,C) pair. The n first columns of U make a basis for the controllable subspace.

The (2,1) block (made of last nx-n rows and n first columns) of U'*A*U is zero and and the n first columns of C*U are zero.


EXAMPLE :

A=diag([1,2,3]);C=[1,0,0];
unobs(A,C)



See Also : contr X, contrss X, canon X, cont_mat X, spantwo X, dt_ility X

4.73   zeropen zero pencil




CALLING SEQUENCE :

[Z,U]=zeropen(Sl)



PARAMETERS :




DESCRIPTION :

Z = sE - F is the zero pencil of the linear system Sl with matrices [A,B,C,D]. Utility function.

With U row compression of [B;D] i.e, U*[B;D]=[0;*]; one has: 
U*[-sI+A |B;   [ Z |0;      
      C  |D] = [ * |*]
The zeros of Z are the zeros of Sl.


See Also : systmat X, kroneck X

Previous Next Contents