// Commande par placement de pôles, partie offline
n=3
Te=0.1
// Matrices du système
A = eye(3,3)+Te*[ 0,1,0 ; 0,0,1 ; -0.1,-0.2,-0.3 ]
B = Te*[ 0 ; 0 ; 1 ]
C = [ 1 , zeros(1,n-1) ]
D = zeros(1,1);
x0=[10;0;0]
// Contrôleur
G=ppol(A,B,0.5*[1,1,1])
Gy=G(1)
Gz=G(2:n)
GI=0
// alternative avec intégrateur
q=1
p=1
Ae = [A,zeros(n,q);Te*eye(q,q),zeros(q,n-q),eye(q,q)]
Be = [B;zeros(q,p)]
//G=ppol(Ae,Be,[0.5*ones(p,n),0.5*ones(p,q)])
Q=diag([100,10,1,1000]);
R=1;
Big=sysdiag(Q,R);               //Now we calculate C1 and D12
[w,wp]=fullrf(Big);
C1=wp(:,1:(n+q));
D12=wp(:,(n+q+1):$);            //[C1,D12]'*[C1,D12]=Big
P=syslin('d',Ae,Be,C1,D12);     // Discrete time plant
[G,X]=lqr(P)     
G=-G;
//G=ppol(Ae,Be,[0.5*ones(p,n),0.5*ones(p,q)])
disp('Spectre du contrôleur')
disp(abs(spec(Ae-Be*G)))
Gy=G(:,1)
GI=G(:,(n+1):(n+q))
G=G(:,1:n)
// Observateur LQ
Q=diag([1,10,100]);
R=1;                            //Usual notations x'Qx + u'Ru
Big=sysdiag(Q,R);               //Now we calculate C1 and D12
[W,Wt]=fullrf(Big);
B1=W(1:n,:);
D21=W(($+1-q):$,:);
P21=syslin('d',A,B1,C,D21);     // Discrete time plant
[K,X]=lqe(P21);
K=-K;
//K=ppol(A',C',0.2*[1,1,1])'
disp('Spectre de l''observateur')
disp(abs(spec(A-K*C)))
Aobs=A-K*C
Bobs=[B,K]