%*****************************************************************************
% COMM.m :
% Synthse et simulation de lois de commande et d'observateurs,
% monovariables et multivariables.
%
% Ce programme permet la rsolution des exercices de l'ouvrage :
%
% COMMANDE ET ESTIMATION MULTIVARIABLES
% Mthodes linaires et optimisation quadratique
% Eric Ostertag
%
% Ed. ELLIPSES, Collection TECHNOSUP, fvrier 2006
%
%*****************************************************************************
% COMM permet le calcul de lois de commande par retour d'tat ou
% d'observateurs (identit, rduit, de perturbation, filtre de Kalman), pour
% cinq procds figurant parmi les exercices rsolus du livre, sous forme
% de modle temps continu ou discret. Trois d'entre eux sont des procds
% multivariables.
%
% Parmi les mthodes proposes figurent le placement de ples, la synthse
% modale complte, les mthodes de dcouplage, les mthodes d'optimisation
% quadratique (LQC) et un algorithme gnral prsent dans le livre.
%
% Le programme propose ensuite de tracer des courbes de rponses diverses,
% puis de passer en simulation pour tester les algorithmes, sans sortir
% dfinitivement, de telle sorte que l'utilisateur peut continuer essayer
% d'autres mthodes tant qu'il travaille sur le mme procd.
%
% Il utilise :
% sous-programmes : Roppenecker, Falb_Wolovich_Roppenecker, Becker_Ostertag,
% saisievalpro, ChoixVpropVparam, arrondi_matrice, responses, scopes_simul
% modles de simulation pour Matlab R2007b (Matlab 7.5 / Simulink 7.0) :
% Simul_univ_continu.mdl, Simul_univ_discret.mdl,
% pendule_et_frottements.mdl
%
% N.B.: en raison de la compatibilit ascendante mais non descendante
% des modles Simulink, il est recommand aux utilisateurs de versions
% antrieures d'utiliser, en les renommant comme ci-dessus, les fichiers
% de simulation suivants:
% modles de simulation pour Matlab 7.1 / Simulink 6.3 (R14SP3) :
% Simul_univ_continu_71.mdl, Simul_univ_discret_71.mdl,
% pendule_et_frottements_71.mdl,
% modles de simulation pour Matlab 6.5.1 / Simulink 5.1 (R13SP1) :
% Simul_univ_continu_65.mdl, Simul_univ_discret_65.mdl,
% pendule_et_frottements_65.mdl,
%
%
% Auteur : Eric Ostertag, 27 juin 2005
% Dernire mise jour : 19 dcembre 2007
%
%*****************************************************************************
echo on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% Commande et Observation Monovariable et Mutlivariable (COMM) %%%%%%%%%
% (Version du 16 avril 2009) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off;
clear; close all; format short g; format compact;
disp(['Date et heure : ' num2str(datestr(now)) ])
echo on;
% ***
% Choix du procd
% ***
echo off;
disp('*');
disp('Procd tudi :');
disp(' 1) pendule invers LIP 100 (Amira)');
disp(' 2) systme trois cuves');
disp(' 3) arotherme PT326 (Feedback)');
disp(' 4) mouvement latral du Boeing 747');
disp(' 5) lecteur de bande magntique');
proc = 0;
while proc < 1 || proc > 5,
proc = input('Votre choix (1 5) --> ');
if isempty(proc), proc = 0; beep; end
end
nom_y = cell(3,1);
if proc == 1, % Pendule invers LIP 100 (Amira)
% Variables d'tat exprimes en units lectriques fournies par les capteurs
% x1 = position du chariot (xc) ; x2 = angle du pendule (theta)
% x3 = vitesse du chariot (xc_dot) ; x4 = vitesse angulaire du pendule (theta_dot)
% temps : s
tscale = 't (s)';
Ac = [0 0 -1.9503 0; 0 0 0 1; 0 -0.1289 -1.9148 0.0008185; 0 21.473 26.3389 -0.1362];
Bc = [0; 0; -6.1344; 84.3027];
Em = [1 0 0 0]'; % E = Em = Matrice de pondration du bruit d'quation : xdot = Ax + Bu + Ev
Cmm = [1 0 0 0; 0 1 0 0]; % Position du chariot et angle du pendule seuls mesurs
nom_y{1} = 'y(1) = position du chariot (unit du capteur, V)';
nom_y{2} = 'y(2) = angle du pendule (unit du capteur, V)';
troismesures = input('Deux mesures : position et angle (par dfaut); changer en trois mesures (O/N) ? [N] --> ','s');
if isempty(troismesures), troismesures = 'N'; end
if troismesures == 'o' || troismesures == 'O'
Cmm(3,:) = [0 0 1 0]; % Vitesse du chariot mesure en plus
nom_y{3} = 'y(3) = vitesse du chariot (unit du capteur, V)';
end
Te = 0.03; % Priode d'chantillonnage du modle discret: Te = 0,03 s
elseif proc == 2, % Systme trois cuves
% x1, x2, x3 = carts des 3 hauteurs de liquide (m)
% u1 ,u2 = dbits d'alimentation (m3/mn)
% temps : mn
tscale = 't (mn)';
Ac = [-0.332 0.332 0; 0.332 -0.664 0.332; 0 0.332 -0.524];
Bc = [0.764 0; 0 0; 0 0.764];
Cmm = [1 0 0; 0 0 1]; % Seuls x1 et x3 sont mesurs
nom_y{1} = 'y(1) = cart de la hauteur de liquide de la 1re cuve (m)';
nom_y{2} = 'y(2) = cart de la hauteur de liquide de la 3me cuve (m)';
Em = [0 1 0]';
Te = 0.1; % Echantillonnage du modle discret Te = 0,1 mn
elseif proc == 3, % Arotherme PT326 (Feedback
s = tf('s');
Gs = 1/(1+0.25*s)^2;
Gs.inputdelay = 0.2;
[num den] = tfdata(Gs, 'v');
[Ar Br Cr Dr] = tf2ss(num, den);
Ac = Ar'; Bc = Cr'; Cmm = Br'; D = Dr'; % Passage la forme canonique d'observabilit
% x1 = temprature de sortie (units du capteur)
% x2 = grandeur non physique
% temps : s
tscale = 't (s)';
nom_y{1} = 'y(1) = temprature la sortie de l''arotherme (units du capteur)';
Em = [1 1]';
Te = 0.2; % Echantillonnage du modle discret Te = 0,2 s
elseif proc == 4, % Mouvement latral du Boeing 747
% x1 = angle de drapage latral, bta (rad); x2 = vitesse angulaire de lacet, r (rad/s)
% x3 = vitesse angulaire de roulis, p(rad/s); x4 = angle de roulis, phi (rad)
% u1 = angle de braquage du gouvernail, delta_g (rad)
% u2 = angle de braquage des ailerons, delta_a (rad)
% temps : s
tscale = 't (s)';
Ac = [-0.0558 -0.9968 0.0802 0.0415; 0.598 -0.115 -0.0318 0; -3.05 0.388 -0.4650 0; 0 0.0805 1 0];
Bc = [0.00729 0.0583; -0.475 -2.009; 0.153 0.0241; 0 0];
Cmm = [0 1 0 0; 0 0 1 0]; % Seuls x2 et x3 sont mesurs
nom_y{1} = 'y(1) = vitesse angulaire de lacet, r (rad/s)';
nom_y{2} = 'y(2) = vitesse angulaire de roulis, p(rad/s)';
Em = [1 0 0 1]';
Te = 0.1; % Echantillonnage du modle discret Te = 0,1 s
elseif proc == 5, % Lecteur de bande magntique
% x1, x2 = position horizontale de la bande la priphrie du cabestan (mm)
% x3, x4 = vitesse angulaire du cabestan (rad/s);
% u1, u2 = courant d'induit des deux moteurs (A)
% temps : s
tscale = 't (s)';
Ac = [0 0 -100 0; 0 0 0 100; 1.25 -1.25 -0.2 -0.2; 1.25 -1.25 -0.2 -0.2];
Bc = [0 0; 0 0; 0.4 0; 0 0.4];
Cmm = [0.5 0.5 0 0; -2 2 0.32 0.32];
% position de la bande au regard de la tte, (x1+x2)/2, et tension de bande
% sont mesures
nom_y{1} = 'y(1) = position de la bande au regard de la tte (mm)';
nom_y{2} = 'y(2) = tension de la bande (N)';
Em = [0 0; 0 0; 1 0; 0 1];
Te = 0.01; % Echantillonnage du modle discret Te = 0,01 s
end
n = size(Ac,1); p = size(Bc,2); q = size(Cmm,1);
D = zeros(q,p); re = size(Em,2);
disp('Ples et zros du procd temps continu, en boucle ouverte : ')
[Poles Zeros] = pzmap(Ac, Bc, Cmm(1:p,:), D(1:p,:)),
[Phi,Gamma,Cmm,D] = c2dm(Ac,Bc,Cmm,D,Te,'zoh'); % Phi et Gamma, mme dim. que A et B
algo = 0; icre = 0; iobs = 0; integ = 0; old_plant = ' '; isimul = 0;
while algo ~= 100,
echo off;
disp('*');
disp('Choix du modle :');
disp(' C) continu (par dfaut)');
if proc == 2
disp(' D) discret (priode Te = 0,1 mn)');
else
disp([' D) discret (priode Te = ' num2str(Te) ' s)']);
end
disp(' Q) quitter le programme');
plant = input('Votre choix --> ','s');
if isempty(plant)
plant = 'C';
end
if plant == 'Q' || plant == 'q',
algo = 100;
elseif plant == 'D' || plant == 'd',
plant = 'D'; A = Phi; B = Gamma; E = Em; T = Te; VpVp = ' temps discret ';
disp('Ples et zros du modle temps discret, en boucle ouverte : ');
[Poles_z Zeros_z] = pzmap(A, B, Cmm(1:p,:), D(1:p,:))
else
plant = 'C'; A = Ac; B = Bc; E = Em; T = 0; VpVp = ' ';
end
if algo ~= 100,
n = size(A,1); Cm = Cmm; integ = 0; qw = 0;
if plant ~= old_plant,
icre=0; iobs=0; old_plant=plant;
end
disp('*');
disp('Menu principal : synthse d''un ...');
disp(' A) retour d''tat (par dfaut)');
disp(' B) observateur identit');
disp(' C) observateur rduit');
if proc == 1,
disp(' D) observateur de "perturbation" (pendule invers)');
end
disp(' E) filtre de Kalman');
syn = input('Votre choix --> ','s');
if syn == 'a',
syn = 'A';
elseif syn == 'b',
syn = 'B';
elseif syn == 'c',
syn = 'C';
elseif syn == 'd',
syn = 'D';
elseif syn == 'e',
syn = 'E';
else
syn = 'A';
end
if syn == 'A', algo = 0;
disp(' Composantes du vecteur de mesure (sortie) :');
for i=1:q
disp([' ' nom_y{i}]);
end
if p == q,
disp('Nombre d''entres = nombre de sorties ==>');
disp([' toutes les ' num2str(q) ' composantes de y seront asservies leur valeur de consigne']);
Sc = eye(p);
elseif p < q,
if p == 1,
l_c = 0;
while l_c < 1 || l_c > q,
l_c = input(['Commande u scalaire ==> \n' ...
' numro de l''unique composante de y asservir sa valeur de consigne --> ']);
if isempty(l_c), l_c = 0; end
end
qc = 1;
else
qc = p+1;
while qc > p,
l_c = input(['Vecteur-ligne avec les numros des composantes de y asservir \n' ...
' une valeur de consigne (' num2str(p) ' valeurs max.) --> ']);
if isempty(l_c),
l_c = 1;
end
qc = length(l_c); l_c = sort(l_c);
end
end
Sc = zeros(qc, q);
for i=1:qc,
Sc(i,l_c(i)) = 1;
end
else
disp(['Systme ayant p > q : supprimez ' num2str(p-q) 'entres de commande']);
beep
break
end
integ = input('Action intgrale : 0 = sans (par dfaut), 1 = avec; Votre choix --> ');
if isempty(integ) || integ ~= 1,
integ = 0; qw = 0;
end
if integ == 1,
qi = q+1;
while qi > q,
l_int = input(['Vecteur-ligne avec les numros des composantes de y reboucler ' ...
'par intgrateur(s) (' num2str(q) ' valeurs max.) --> ']);
qi = length(l_int); l_int = sort(l_int);
end
Si = zeros(qi, q);
for i=1:qi,
Si(i,l_int(i)) = 1;
end
Ci = Si*Cm;
B = [B; zeros(qi,p)]; Cm = [Cm zeros(q,qi)]; B_cons = [zeros(n,q); Si];
if plant == 'C',
a22 = zeros(qi);
else
a22 = eye(qi);
end
A = [A zeros(n,qi); -Ci a22], B, n = n+qi; qw = qi;
else
Si = zeros(q,q); qi = 0;
end
disp(['Rang de la matrice de commandabilit Qs : ' num2str(rank(ctrb(A,B)))])
if rank(ctrb(A,B)) ~= n,
disp(' *** Procd non commandable ***'); beep;
return;
end
disp('Algorithme utiliser :');
disp(' 1) placement de ples (modules MATLAB "acker.m" ou "place.m")');
if integ == 0,
disp(' 2) commande modale simple');
end
disp(' 3) synthse modale complte (formule de Roppenecker)');
if integ == 0,
disp(' 4) mthode du dcouplage (Falb-Wolovich ou Roppenecker)');
end
disp(' 5) formule gnrale (Becker-Ostertag)');
disp(' 6) critre quadratique (LQC), avec pondration de l''tat : xT*Q*x');
disp(' 7) critre quadratique (LQC), avec pondration de la sortie : xT*CT*Q*C*x');
elseif syn == 'B' || syn == 'C',
algo = 0;
if syn == 'C',
c_shape = input('Forme de C : 1 = (0|Iq), 2 = (Iq|0); Votre choix --> ');
if c_shape == 1, % Cas C=(0|Iq)
A11 = A(1:n-q, 1:n-q); A12 = A(1:n-q, n-q+1:end);
A21 = A(n-q+1:end, 1:n-q); A22 = A(n-q+1:end, n-q+1:end);
B1 = B(1:n-q, :); B2 = B(n-q+1:end, :);
else % Cas C=(Iq|0) : permutation des indices 1 et 2
A22 = A(1:q, 1:q); A21 = A(1:q, q+1:end);
A12 = A(q+1:end, 1:q); A11 = A(q+1:end, q+1:end);
B2 = B(1:q, :); B1 = B(q+1:end, :);
end
end
if (syn == 'B' && rank(obsv(A,Cm)) ~= n) || (syn == 'C' && rank(obsv(A11,A21)) ~= n-q),
disp(' *** Procd non observable ***'); beep;
return;
end
disp('Algorithme utiliser :');
disp(' 1) placement de ples multivariable (module "place.m")');
disp(' 3) formule de Roppenecker');
disp(' 5) formule gnrale (Becker-Ostertag)');
disp(' 6) critre quadratique "LQC"');
elseif syn == 'D',
algo = 10;
elseif syn == 'E',
algo = 11;
end
if algo ~= 10 && algo ~= 11,
while algo<1 || algo>7 && algo~=99,
algo = input('Votre choix --> ');
if isempty(algo),
algo = 99;
end
end
end
disp (['Valeurs propres et vecteurs propres du procd' VpVp ':']);
[Vp, Lambda] = eig(A,'nobalance')
end
switch algo
case 1,
echo on;
%
% *** 1.- Placement de ples mono- ou multivariable ("acker.m" ou "place.m")
%
echo off;
if syn == 'A',
disp(['N.B.: Pour pouvoir utiliser "place.m", ordre de multiplicit des valeurs propres '...
'<= ' num2str(p) ' (nombre d''entres)']);
else
disp(['N.B.: Pour pouvoir utiliser "place.m", ordre de multiplicit des valeurs propres '...
'<= ' num2str(q) ' (nombre de sorties)']);
end
if syn == 'A' || syn == 'B',
lambda = saisievalpro(n,T); rB = rank(B); nb = n;
end
if syn == 'C',
lambda = saisievalpro(n-q,T); rB = rank(A21); nb = n-q;
end
ps = sort(lambda); pmult = 0;
for i=1:nb-rB
imax = min(nb,i+rB);
if all(ps(i:imax) == ps(i)),
pmult = 1;
end
end
if (syn == 'A' && p == 1) || ((syn == 'B' || syn == 'C') && q == 1) ,
if pmult == 0,
method = input('Choix du module : 1 = acker.m, 2 = place.m; Votre choix --> ');
else
method = 1;
end
else
if pmult == 0,
method = 2;
else
method = 0;
end
end
if method == 0,
if syn == 'A',
disp('Il y a des ples de multiplicit > p ==> faites un autre choix');
else
disp('Il y a des ples de multiplicit > q ==> faites un autre choix');
end
iobs = 0;
algo = 99;
beep;
else
if method == 1,
if syn == 'A',
Lpp = acker(A,B,lambda); L = Lpp, icre = 1;
end
if syn == 'B',
Gpp = acker(A',Cm',lambda); G = Gpp', iobs = 1;
end
if syn == 'C',
Gpp = acker(A11',A21',lambda); Gr = Gpp', iobs = 1;
end
elseif method == 2,
if syn == 'A',
Lpp = place(A,B,lambda); L = Lpp, icre = 1;
end
if syn == 'B',
Gpp = place(A',Cm',lambda); G = Gpp', iobs = 1;
end
if syn == 'C',
Gpp = place(A11',A21',lambda); Gr = Gpp', iobs = 1;
end
end
end
case 2,
echo on;
%
% *** 2.- Calcul d'une commande modale simple
%
echo off;
if integ == 1,
disp('*** Synthse uniquement sans intgrateur ***'); icre=0; beep;
else
Wp = inv(Vp); % les lignes de Wp sont les vecteurs propres gauches de A
disp(['Seule(s) ' num2str(p) ' valeur(s) propre(s) pourra(ont) tre dplace(s)']);
j = 0;
for i=1:n
dep = input(['Dplacer la valeur propre lambda' num2str(i) ' = ' num2str(Lambda(i,i)) ' (O/N) ? --> '],'s');
if dep == 'O' || dep == 'o',
j = j+1; lambda = input(' Nouvelle valeur (continue) --> ');
if plant == 'D',
lambda = exp(lambda*Te);
end
Ti_p(j,:) = Wp(i,:); lambdadiff(j) = Lambda(i,i)-lambda;
end
if j == p,
break;
end
end
if j>0,
Bpchap = Ti_p*B; Lmodale = inv(Bpchap)*diag(lambdadiff)*Ti_p;
L = Lmodale, icre = 1;
end
end
case 3,
echo on;
%
% *** 3.- Synthse modale complte (Formule de Roppenecker)
%
echo off;
if syn == 'A',
[Lrop,lambda_cre] = Roppenecker(A,B,Lambda,Vp,[],T); L = Lrop, icre = 1;
end
if syn == 'B',
[Grop,lambda_obs] = Roppenecker(A',Cm',Lambda,Vp,[],T); G = Grop',iobs = 1;
end
if syn == 'C',
[Grop,lambda_obs] = Roppenecker(A11',A21',Lambda,Vp,[],T);Gr=Grop', iobs = 1;
end
case 4,
echo on;
%
% *** 4.- Mthode du dcouplage, selon Falb-Wolovich ou Roppenecker
%
echo off;
if integ == 1,
disp('*** Synthse uniquement sans intgrateur ***'); icre=0; beep;
else
[L, M] = Falb_Wolovich_Roppenecker(A,B,Cm,T);
if L == 0,
icre = 0;
else
L, M,
% Vrification de la position des zros ("instables" ou non ?)
A_BF = arrondi_matrice(A-B*L); G_BF = arrondi_matrice(Cm(1:p,:)*inv(-A_BF)*B);
BM = B*M; sys_BF = ss(A_BF, BM, Cm, D, T);
disp('Matrice de transfert du systme boucl :');
Gw = tf(sys_BF)
disp('Zro(s) du systme boucl : ');
Z = zero(sys_BF)
% Affichage de la matrice de transfert dcouple
Msys_BF = minreal(sys_BF); % supprimer modes non commandables et non observables
icre = 2;
end
end
case 5,
echo on;
%
% *** 5.- Formule gnrale (Becker-Ostertag)
%
echo off;
if syn == 'A',
[Lbo,lambda_cre] = Becker_Ostertag(A,B,Lambda,[],T); L = Lbo, icre = 1;
end
if syn == 'B',
[Gbo,lambda_obs] = Becker_Ostertag(A',Cm',Lambda,[],T); G = Gbo', iobs = 1;
end
if syn == 'C',
[Gbo,lambda_obs] = Becker_Ostertag(A11',A21',Lambda,[],T); Gr = Gbo', iobs = 1;
end
case {6,7},
echo on;
%
% *** 6-7.- Calcul par critre quadratique (LQC), avec pondration de l'tat ou de la sortie
%
echo off;
Q = []; diag_Q = []; R = []; diag_R = [];
if algo == 6,
disp('*** Pondration de l''tat : xT*Q*x ***');
Cw = eye(n);
else
disp('*** Pondration de la sortie : yT*Q*y ***');
Cw = [Cm; zeros(qw,n-qw) eye(qw)];
end
nw = size(Cw,1);
diagoui = input('Matrice Q diagonale (O/N) ? [O] --> ','s');
if isempty(diagoui)
diagoui = 'O';
else
diagoui = diagoui(1);
end
while size(Q,1) ~= nw,
if diagoui == 'O' || diagoui == 'o'
diag_Q = input(['Vecteur-ligne contenant les ' num2str(nw) ' lments diagonaux qii de Q --> ']);
Q = diag(diag_Q);
else
Q = input(['Matrice Q (' num2str(nw) 'x' num2str(nw) ') = ']);
end
end
Qw = Cw'*Q*Cw;
if syn == 'A',
nR = p;
elseif syn == 'B',
nR = q;
end
diagoui = input('Matrice R diagonale (O/N) ? [O] --> ','s');
if isempty(diagoui)
diagoui = 'O';
else
diagoui = diagoui(1);
end
while size(R,1) ~= nR,
if diagoui == 'O' || diagoui == 'o'
diag_R = input(['Vecteur-ligne contenant les ' num2str(nR) ' lments diagonaux rii de R --> ']);
R = diag(diag_R);
else
R = input(['Matrice R (' num2str(nR) 'x' num2str(nR) ') = ']);
end
end
if plant == 'C',
if syn == 'A',
[Lopt, P, Lambda_BF] = lqr(A,B,Qw,R); L = Lopt, icre = 1;
end
if syn == 'B',
[Gopt, Sigmac, Lambda_obs] = lqr(A',Cm',Qw,R); G = Gopt', iobs = 1;
end
elseif plant == 'D',
if syn == 'A',
[Lopt, P, Lambda_BF] = dlqr(A,B,Qw,R); L = Lopt, icre = 1;
end
if syn == 'B',
[Gopt, Sigmad, Lambda_obs] = dlqr(A',Cm',Qw,R); G = Gopt', iobs = 1;
end
end
case 10,
echo on;
%
% *** 10.- Calcul d'un observateur du systme discret augment de la perturbation constante
% (frottement sec), pour le pendule invers, uniquement.
%
echo off;
if plant == 'C',
disp('*** Synthse uniquement en discret ***'); algo = 0; beep;
elseif plant == 'D',
if proc ~= 1,
disp('*** Pendule invers uniquement ***'); algo = 0; beep;
else
% Modle continu augment de la perturbation constante (frottement sec)
A_tild = [Ac Bc; zeros(1,n) zeros(1,p)], B_tild = [Bc; zeros(1,p)],
C_tild = [Cm zeros(q,p)], D_tild = D,
% Modle discret augment correspondant
[A,B,Cm,D] = c2dm(A_tild,B_tild,C_tild,D_tild,Te,'zoh'); Phi_tild = A, Gamma_tild = B,
n = size(A,1); q = size(Cm,1);
[Vp_aug,Lambda_aug] = eig(A,'nobalance');
typeobs = input('Type d''observateur : 1 = identit, 2 = rduit; Votre choix --> ');
disp('Mthode de synthse choisie :');
disp(' 1 = placement de ples (module "place.m") (dfaut)');
disp(' 3 = formule de Roppenecker');
disp(' 5 = formule de Becker-Ostertag');
method = input('Votre choix --> ');
if isempty(method)
method = 1;
end
if typeobs == 1
switch method;
case 3
[Lobs1,lambda_obs] = Roppenecker(A',Cm',Lambda_aug,Vp_aug,[],Te);
G = Lobs1'
case 5
[Lobs2,lambda_obs] = Becker_Ostertag(A',Cm',Lambda_aug,[],Te);
G = Lobs2'
otherwise
disp(['Introduisez n=' num2str(n) ' valeurs propres (continues)' ...
' pour l''observateur identit ===> ...']);
disp(['N.B.: ordre de multiplicit <= ' num2str(q) ' (nombre de sorties)']);
lambda_obs = saisievalpro(n,Te);
G = (place(A',Cm',lambda_obs))'
end
elseif typeobs == 2
% Partitionnement des matrices entre parties mesures : x1 xq,
% et parties estimes : x(q+1) xn, plus la perturbation
c_shape = 2;
A22 = A(1:q, 1:q); A21 = A(1:q, q+1:end);
A12 = A(q+1:end, 1:q); A11 = A(q+1:end, q+1:end);
B2 = B(1:q, :); B1 = B(q+1:end, :);
switch method;
case 3
[Grop,lambda_obs] = Roppenecker(A11',A21',Lambda_aug,Vp_aug,[],Te); Gr = Grop'
case 5
[Gbo,lambda_obs] = Becker_Ostertag(A11',A21',Lambda_aug,[],Te); Gr = Gbo'
otherwise
disp(['Introduisez r=' num2str(n-q) ' valeurs propres (continues)' ...
' pour l''observateur rduit ===> ...']);
disp(['N.B.: ordre de multiplicit <= ' num2str(q) ' (nombre de sorties)']);
lambda_obs = saisievalpro(n-q,Te);
Gr = (place(A11',A21',lambda_obs))'
end
end
iobs = 2;
end
end
case 11,
echo on;
%
% *** 11.- Calcul d'un filtre de Kalman ("kalman.m")
%
echo off;
diag_Qbar = []; diag_Rbar = [];
% a) Adaptation du modle LTI du procd la structure de "kalman.m" (cf. aide en ligne)
ext_mod = ss(Phi,[Gamma Gamma],Cm,[D zeros(q,p)],Te);
% b) Calcul du filtre de Kalman
disp('Matrice de variance du bruit d''quation (bruit de procd) :');
while length(diag_Qbar) ~= p,
diag_Qbar = input(['Vecteur-ligne contenant les ' num2str(p) ' lments diagonaux qii de Q_barre --> ']);
end
disp('Matrice de variance du bruit de mesure :');
while length(diag_Rbar) ~= q,
diag_Rbar = input(['Vecteur-ligne contenant les ' num2str(q) ' lments diagonaux rii de R_barre --> ']);
end
Qbar = diag(diag_Qbar); Rbar = diag(diag_Rbar);
[KEST,PhiK,Sigma_moins,Kk] = kalman(ext_mod,Qbar,Rbar); Kk,
E = Gamma; G = PhiK, iobs = 1;
case 100,
break
otherwise
disp('*** Refaites un choix ***'); beep;
end
%
% Achvement des synthses (observateur et systme boucl)
%
if algo ~= 99,
if iobs == 1 || iobs == 2,
if syn == 'B' || (syn == 'D' && typeobs == 1) || syn == 'E'
F = A-G*Cm, J = B; H = eye(n); K = zeros(n,q); % Terminer synthse des observateurs
elseif syn == 'C' || (syn == 'D' && typeobs == 2),
F = A11-Gr*A21, G = F*Gr+A12-Gr*A22, J = B1-Gr*B2
if c_shape == 1
H = [eye(n-q); zeros(q,n-q)], K = [Gr; eye(q)];
else
H = [zeros(q,n-q); eye(n-q)], K = [eye(q); Gr]
end
end
end
if icre == 1 || icre == 2, % Terminer synthse des retours d'tat et de la matrice M
if integ == 1
L1 = L(:,1:n-qi), L2 = L(:,n-qi+1:n),
else
L1 = L; L2 = zeros(p,q);
end
if plant == 'C',
A_BF1 = Ac - Bc*L1; G_BF = Sc*Cmm*inv(-A_BF1)*Bc;
else
A_BF1 = Phi - Gamma*L1; G_BF = Sc*Cmm*inv(eye(size(A_BF1))-A_BF1)*Gamma;
end
if icre ~= 2,
if rank(G_BF) < p,
M = eye(p),
else
M = inv(G_BF),
end
end
disp('Ples et zros du systme boucl : ')
if integ == 1,
A_BF = A-B*L;
B_BF = B*M*Sc + B_cons;
C_BF = Cm;
D_BF = zeros(q,q);
else
A_BF = A_BF1;
B_BF = B*M*Sc;
C_BF = Cmm;
D_BF = zeros(q,q);
end
sys_BF = ss(A_BF, B_BF, C_BF, D_BF, T);
[Poles Zeros] = pzmap(sys_BF)
%
% Paramtres des divers procds pour le trac de courbes et les simulations
%
if proc == 1 % Pendule invers LIP 100 (Amira)
x0 = [0 -2 0.5 0]; yc = [1 0]; start_consigne = [0 0]; tfinal = 9;
p_charge = [0.5]; start_p_charge = [5];
p_mesures = [0.05 0.1]; start_p_mesures = [7]; horizon_simul = 300*Te;
eqn_pwr = 0.01*Te; eqn_seed = [23341];
meas_pwr = [0.02 0.01]*Te ; meas_seed = [14223 42313];
if troismesures == 'o' || troismesures == 'O'
yc = [yc 0]; start_consigne = [start_consigne 0];
p_mesures = [p_mesures 0];
meas_pwr = [meas_pwr 0.01*Te] ; meas_seed = [meas_seed 31424];
end
end
if proc == 2, % Systme trois cuves
x0 = [-0.8 -0.5 -0.2]; yc = [0.4 0.8]; start_consigne = [0 0]; tfinal = 20;
p_charge = [0.5]; start_p_charge = [15];
p_mesures = [1 1]; start_p_mesures = [20 25]; horizon_simul = 300*Te;
eqn_pwr = 0.01*Te; eqn_seed = [23341];
meas_pwr = [0.02 0.01]*Te ; meas_seed = [14223 42313];
end
if proc == 3, % Arotherme PT326 (Feedback
x0 = [10 0]; yc = [1]; start_consigne = [0]; tfinal = 30;
p_charge = [0.5]; start_p_charge = [15];
p_mesures = [1]; start_p_mesures = [20]; horizon_simul = 100*Te;
eqn_pwr = 0.01*Te; eqn_seed = [23341]; meas_pwr = [0.01]*Te ; meas_seed = [14223];
end
if proc == 4, % Mouvement latral du Boeing 747
x0 = [0.15 0 0 0.5]; yc = [1 1]; start_consigne = [0 0]; tfinal = 30;
p_charge = [0.5]; start_p_charge = [15];
p_mesures = [1]; start_p_mesures = [20]; horizon_simul = 100*Te;
eqn_pwr = 0.01*Te; eqn_seed = [23341];
meas_pwr = [0.02 0.01]*Te ; meas_seed = [14223 42313];
end
if proc == 5, % Lecteur de bande magntique
x0 = [-1 2 0 0]; yc = [1 2]; start_consigne = [0 0]; tfinal = 2;
p_charge = [1 2]; start_p_charge = [50*Te 60*Te];
p_mesures = [0.2]; start_p_mesures = [75*Te]; horizon_simul = 100*Te;
eqn_pwr = 0.01*Te; eqn_seed = [23341 34412];
meas_pwr = [0.02^2 0.01^2]*Te ; meas_seed = [14223 42313];
end
xx0 = x0;
if T == 0,
t = 0:tfinal/200:tfinal;
else
t = 0:T:200*T;
end
%
% Trac de rponses indicielles du systme boucl
%
icre = 0; isimul = 1;
responses
end
%
% Simulations, si dsir
%
if isimul == 1,
sim = input('Dsirez-vous simuler la rgulation (O/N) ? [N] --> ','s');
if isempty(sim)
sim = 'N';
else
sim = sim(1);
end
if sim == 'O' || sim == 'o',
if integ == 1 || iobs == 2
if plant == 'C'
A = Ac; B = Bc;
else
A = Phi; B = Gamma;
end
end
n = size(A,1); Cm = Cmm; q = size(Cm,1); z0 = 0;
if iobs == 0,
F = zeros(n); G=zeros(n,q); J=zeros(n,p); H=zeros(n); K=zeros(n,q);
end
x0_nul = input('Etat initial du procd nul (O/N) ? [O] --> ','s');
if isempty(x0_nul)
x0_nul = 'O';
else
x0_nul = x0_nul(1);
end
if x0_nul == 'O' || x0_nul == 'o',
x0 = zeros(n,1);
else
x0 = xx0;
end
if plant == 'C',
disp('Simulez avec "Simul_univ_continu.mdl" ===>');
elseif plant == 'D',
if iobs == 2,
disp('Simulez avec "pendule_et_frottements.mdl" ===>');
else
disp('Simulez avec "Simul_univ_discret.mdl" ===>');
end
end
disp('(Aprs la simulation, tapez "return" pour retourner au programme)')
keyboard
scopes_simul % Trac d'oscillogrammes relevs dans la simulation
end
end
end
end