"COMM" : Commande et Observation Mono- et Multivariables

by

 

27 Jun 2005 (Updated )

Companion software to french book "Commande et estimation multivariables"

COMM.m
%*****************************************************************************
%    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

Contact us