École Nationale des Sciences Appliquées d'Agadir

Génie Électrique  Option GSEMD

20252026

TP 02  Commande Prédictive (MPC)

Approche 2 : Modèle augmenté
avec état étendu (x, u(k − 1))
Élément 2 : Commandes Intelligentes

Réalisé par :

BENJIMA Mohammed

May 2, 2026
ENSA Agadir  2025/2026 TP 02  MPC Approche 2

Contents
1 Code MATLAB commenté et interprété 2
1.1 Section 1  Initialisation et paramètres . . . . . . . . . . . . . . . . . . . . 2
1.2 Section 2  Réponse indicielle en boucle ouverte . . . . . . . . . . . . . . . 3
1.3 Section 3  Calcul de Σi , Φp et Gf . . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Section 4  Calcul de Ψp , Gc et gc⊤ . . . . . . . . . . . . . . . . . . . . . . 5
1.5 Section 5  Simulation en boucle fermée . . . . . . . . . . . . . . . . . . . 6
1.6 Section 6  Inuence des paramètres Np , Nc , λ . . . . . . . . . . . . . . . . 8

2 Conclusion 11

1
ENSA Agadir  2025/2026 TP 02  MPC Approche 2

1 Code MATLAB commenté et interprété

1.1 Section 1  Initialisation et paramètres
Rappel mathématique
On charge les matrices A, B , C du modèle discret et les scalaires Nc , Np , λ qui
pilotent tout le calcul MPC.

1 clear ; clc ; close all ;

2

3 % Matrices du modele discret

4 A = [0 , 1 ;
5 -1.015 , 1.4]; % matrice d ' etat (2 x2 )
6

7 B = [0; 1]; % matrice d ' entree (2 x1 )

8 C = [ -0.7 , 1]; % matrice de sortie (1 x2 )
9

10 % Parametres MPC
11 Nc = 2; % Horizon de commande ( nombre d ' increments
libres )
12 Np = 4; % Horizon de prediction
13 lambda = 0.5; % Poids sur Delta u ( regularisation )
14

15 fprintf ( 'Nc = %d , Np = %d , lambda = %.2 f\ n ' , Nc , Np , lambda );

Listing 1: Initialisation générale et paramètres du modèle

Interprétation
ˆ A possède des valeurs propres λ1,2 telles que |λi | > 1 ⇒ le système est instable
en boucle ouverte.
ˆ Np = 4 : le contrôleur anticipe sur 4 pas futurs.

ˆ Nc = 2 < Np : seuls 2 incréments ∆u sont optimisés, ce qui réduit la dimension

du problème QP et améliore la robustesse.

ˆ λ = 0,5 : pénalise les variations trop brusques de la commande ; un λ plus

grand ⇒ commande plus lente mais plus douce.

2
ENSA Agadir  2025/2026 TP 02  MPC Approche 2

1.2 Section 2  Réponse indicielle en boucle ouverte

Rappel mathématique
On simule x(k +1) = Ax(k)+Bu(k), y(k) = Cx(k) pendant 1000 pas avec u(k) = 1
(échelon unité) et conditions initiales nulles.

1 N_sim = 1000;
2 x_bo = zeros (2 , N_sim +1) ; % etat : conditions initiales nulles
3 y_bo = zeros (1 , N_sim ); % sortie
4 u_bo = ones (1 , N_sim ); % echelon unite u =1
5

6 for k = 1: N_sim
7 y_bo ( k) = C * x_bo (: , k ); % y (k) = C x(k )
8 x_bo (: , k +1) = A * x_bo (: , k) + B * u_bo ( k); % x (k +1) = Ax (k)+
Bu ( k)
9 end
10

11 figure ( ' Color ' ,' white ');

12 plot (1: N_sim , y_bo , 'b ', ' LineWidth ', 1.5) ;
13 xlabel ( ' Pas k ') ; ylabel ( 'y(k ) ') ;
14 title ( ' Reponse indicielle en boucle ouverte (1000 pts ) ');
15 grid on ;

Listing 2: Simulation de la réponse indicielle en boucle ouverte

Interprétation
La sortie y(k) diverge (amplitude croissante) car les valeurs propres de A sont en
dehors du cercle unité. Cela conrme la nécessité d'une commande en boucle fermée
(MPC).

1.3 Section 3  Calcul de Σi , Φp et Gf

Rappel mathématique
i−1
X
Σi = CAj B, i = 1, . . . , Np .
j=0

Σi représente la réponse impulsionnelle cumulée du système sur i pas.

1 % --- 2 a) Calcul des Sigma_i (i = 1 : Np ) -----------------------

2 Sigma = zeros (1 , Np );

3
ENSA Agadir  2025/2026 TP 02  MPC Approche 2

3 for i = 1: Np
4 s = 0;
5 for j = 0:i -1
6 s = s + C * (A^j ) * B; % CA ^j B ( scalaire ici )
7 end
8 Sigma (i ) = s;
9 end
10

11 fprintf ( ' Sigma_i = ') ; disp ( Sigma );

12

13 % --- 2 b) Formation de Phi_p ( Np x 1) ---------------------------

14 Phi_p = Sigma (:) ; % vecteur colonne
15

16 % --- Formation de G_f ( Np x Nc ) : matrice de Toeplitz ----------

17 % G_f (i ,j ) = Sigma (i -j +1) si i >= j , 0 sinon
18 G_f = zeros ( Np , Nc ) ;
19 for i = 1: Np
20 for j = 1: Nc
21 if i >= j
22 G_f (i ,j) = Sigma (i - j + 1) ;
23 end
24 end
25 end
26

27 fprintf ( ' Phi_p =\ n '); disp ( Phi_p );

28 fprintf ( ' G_f =\ n ') ; disp ( G_f );

Listing 3: Calcul des Σi , de Φp et de Gf

Interprétation
ˆ Φp (Np × 1) : traduit l'eet de la commande passée u(k − 1) sur les sorties
futures. Sa i-ème composante est Σi .

ˆ Gf (Np ×Nc ) : matrice de Toeplitz inférieure. La colonne j contient les Σi−j+1

décalés ; chaque colonne représente l'impact d'un incrément ∆u à l'instant
k + j − 1 sur les sorties futures.

4
ENSA Agadir  2025/2026 TP 02  MPC Approche 2

1.4 Section 4  Calcul de Ψp , Gc et gc⊤

Rappel mathématique
 
CA
 CA2
 

−1
Gc = λI + G⊤ G⊤ gc⊤ = 1ère ligne de Gc .

Ψp = 
 ..
, f Gf f,
 .


CANp

1 n = size (A , 1) ; % dimension de l ' etat ( n =2)

2

3 % --- Psi_p ( Np x n) : [ CA ; CA ^2 ; ... ; CA ^ Np ] ---------------

4 Psi_p = zeros ( Np , n);
5 for i = 1: Np
6 Psi_p (i ,:) = C * (A ^i) ;
7 end
8

9 % --- G_c = ( lambda *I + G_f '* G_f ) ^( -1) * G_f ' ------------------
10 % Resolution du systeme lineaire ( plus stable que inv () )
11 G_c = ( lambda * eye ( Nc ) + G_f ' * G_f ) \ ( G_f ') ;
12

13 % --- g_c ^ T = premiere ligne de G_c (1 x Np ) -------------------

14 gc_T = G_c (1 ,:) ;
15

16 fprintf ( ' Psi_p =\ n '); disp ( Psi_p );

17 fprintf ( ' G_c =\ n ') ; disp ( G_c );
18 fprintf ( ' gc_T =\ n '); disp ( gc_T );

Listing 4: Calcul de Ψp , Gc et gc⊤

Interprétation
ˆ Ψp (Np × n) : encode la propagation libre de l'état x(k) sur l'horizon de
prédiction, pondérée par la matrice C .

ˆ Gc : solution de la minimisation du critère quadratique J = ∥ŷf − yR ∥2 +

λ∥∆uf ∥2 . L'opérateur (λI + G⊤
f Gf ) Gf est la
−1 ⊤
pseudo-inverse régularisée de
Gf  régularisation de Tikhonov.

ˆ gc⊤ : seule la première ligne est utilisée car on n'applique que ∆u∗ (k) (stratégie
receding horizon ). Les autres lignes ∆u∗ (k + 1), . . . sont recalculées au pas
suivant.

5
ENSA Agadir  2025/2026 TP 02  MPC Approche 2

1.5 Section 5  Simulation en boucle fermée

Rappel mathématique
Loi de commande appliquée à chaque instant k :

u(k) = u(k − 1) + gc⊤ yR − Ψp x(k) − Φp u(k − 1) .

Signal de référence :

yR = [50, . . . , 50, −50, . . . , −50, 50, . . . , 50, −50, . . . , −50] ∈ R1000 .

| {z } | {z } | {z } | {z }
300 200 250 250

1 % --- Signal de reference (1000 points ) -------------------------

2 yr_signal = [50* ones (1 ,300) , -50* ones (1 ,200) , ...
3 50* ones (1 ,250) , -50* ones (1 ,250) ];
4 N_bf = length ( yr_signal );
5

6 % --- Initialisation --------------------------------------------

7 x_bf = zeros (2 , N_bf +1) ; % etat
8 y_bf = zeros (1 , N_bf ); % sortie
9 u_bf = zeros (1 , N_bf +1) ; % commande : u_bf (1) = u ( -1) = 0
10 e_bf = zeros (1 , N_bf ); % erreur
11

12 for k = 1: N_bf
13 y_bf ( k) = C * x_bf (: , k); % sortie courante y( k) = C x (k)
14

15 % Vecteur de reference future yR ( Np x 1)

16 yR = zeros (Np , 1) ;
17 for j = 1: Np
18 idx = min (k+j -1 , N_bf ) ; % saturation en fin de
signal
19 yR ( j) = yr_signal ( idx );
20 end
21

22 % Loi de commande MPC Approche 2

23 u_prev = u_bf ( k); % u (k -1)
24 delta_u = gc_T * ( yR - Psi_p * x_bf (: , k) ...
25 - Phi_p * u_prev ); % Delta u
*( k )
26 u_bf ( k +1) = u_prev + delta_u ; % u (k)
27

28 % Erreur de poursuite

6
ENSA Agadir  2025/2026 TP 02  MPC Approche 2

29 e_bf ( k) = y_bf (k ) - yr_signal ( k);

30

31 % Evolution de l ' etat

32 x_bf (: , k +1) = A* x_bf (: , k) + B* u_bf (k +1) ;
33 end
34

35 % --- Trace des resultats (3 sous - figures ) ----------------------

36 figure ( ' Color ' ,' white ', ' Position ' ,[100 100 900 700]) ;
37

38 subplot (3 ,1 ,1) ;
39 plot (1: N_bf , u_bf (2: end ) , 'k ', ' LineWidth ', 1.4) ;
40 xlabel ( 'k ') ; ylabel ( 'u(k ) ');
41 title ( ' Subplot 1 : Signal de commande u (k) ');
42 grid on ;
43

44 subplot (3 ,1 ,2) ;
45 plot (1: N_bf , yr_signal , 'r -- ', ' LineWidth ', 1.5) ; hold on ;
46 plot (1: N_bf , y_bf , 'b ', ' LineWidth ' , 1.4) ;
47 xlabel ( 'k ') ; ylabel ( ' Amplitude ');
48 title ( ' Subplot 2 : Reference y_r (k ) vs Sortie y(k ) ') ;
49 legend ( ' y_r (k) ','y(k) ', ' Location ',' best '); grid on ;
50

51 subplot (3 ,1 ,3) ;
52 plot (1: N_bf , e_bf , 'm ', ' LineWidth ' , 1.2) ;
53 xlabel ( 'k ') ; ylabel ( 'e(k ) ');
54 title ( ' Subplot 3 : Erreur e( k) = y(k ) - y_r (k) ');
55 grid on ;
56

57 sgtitle ( sprintf ( ' MPC Approche 2 -- Nominal : Np =%d , Nc =%d ,

lambda =%.2 f ' ,...
58 Np , Nc , lambda ) , ' FontWeight ',' bold ',' FontSize ' ,13) ;

Listing 5: Simulation en boucle fermée avec la loi MPC

7
ENSA Agadir  2025/2026 TP 02  MPC Approche 2

Interprétation
ˆ Subplot 1 : la commande u(k) eectue de grands sauts lors des changements
de référence, puis se stabilise  comportement attendu d'un MPC.

ˆ Subplot 2 : y(k) suit dèlement yr (k) après un transitoire. L'erreur statique

doit être nulle en régime (grâce à l'intégration implicite via ∆u).

ˆ Subplot 3 : l'erreur e(k) converge vers zéro à chaque échelon. La vitesse de

convergence dépend de Np , Nc et λ.

1.6 Section 6  Inuence des paramètres Np , Nc , λ

1 cases = {
2 % Nom Np Nc lambda
3 ' Nominal ', 4 , 2, 0.50;
4 ' Cas 1 ', 12 , 2, 2.00;
5 ' Cas 2 ', 8, 4 , 0.50;
6 ' Cas 3 ', 5, 4 , 0.05;
7 ' Cas 4 ', 8, 2 , 0.50;
8 ' Cas 5 ', 4, 1 , 5.00;
9 };
10

11 fprintf ( '\n % -10 s % -5 s % -5 s % -8 s % -15 s % -15 s\n ' ,...

12 ' Cas ' ,'Np ', 'Nc ',' lambda ' ,' Depassement (%) ', ' Temps reponse ');
13 fprintf ( '%s \n ', repmat ( '-' ,1 ,65) ) ;
14

15 for c = 1: size ( cases ,1)

16 cas_name = cases {c ,1};
17 Np_c = cases {c ,2}; Nc_c = cases {c ,3}; lam_c = cases {c ,4};
18

19 % Recalcul complet pour ce jeu de parametres

20 Sigma_c = zeros (1 , Np_c );
21 for i = 1: Np_c
22 s = 0;
23 for j = 0:i -1 , s = s + C *( A^j )*B ; end
24 Sigma_c ( i) = s;
25 end
26

27 Phi_c = Sigma_c (:) ;

28 Gf_c = zeros ( Np_c , Nc_c );
29 for i = 1: Np_c

8
ENSA Agadir  2025/2026 TP 02  MPC Approche 2

30 for j = 1: Nc_c
31 if i >= j , Gf_c (i , j) = Sigma_c (i -j +1) ; end
32 end
33 end
34

35 Psi_c = zeros ( Np_c , n );

36 for i = 1: Np_c , Psi_c (i ,:) = C *( A^ i); end
37

38 Gc_c = ( lam_c * eye ( Nc_c ) + Gf_c '* Gf_c ) \ Gf_c ';

39 gcT_c = Gc_c (1 ,:) ;
40

41 % Simulation boucle fermee

42 xc = zeros (2 , N_bf +1) ;
43 yc = zeros (1 , N_bf );
44 uc = zeros (1 , N_bf +1) ;
45 for k = 1: N_bf
46 yc ( k) = C * xc (: , k );
47 yR_c = zeros ( Np_c , 1) ;
48 for j = 1: Np_c
49 yR_c ( j) = yr_signal ( min (k +j -1 , N_bf ));
50 end
51 u_prev_c = uc (k) ;
52 uc ( k +1) = u_prev_c + gcT_c *( yR_c - Psi_c * xc (: , k ) -
Phi_c * u_prev_c );
53 xc (: , k +1) = A* xc (: , k ) + B * uc (k +1) ;
54 end
55

56 % Indicateurs de performance ( premier palier : ref = 50)

57 seg = yc (1:300) ;
58 overshoot = max (( max ( seg ) - 50) /50 * 100 , 0) ;
59 tr_idx = find ( abs ( seg - 50) <= 2.5 , 1 , ' first ') ; % 5%
de 50
60 if isempty ( tr_idx ) , tr_str = ' >300 '; else , tr_str = num2str
( tr_idx ) ; end
61

62 fprintf ( '% -10 s % -5 d % -5 d % -8.2 f % -15.1 f % -15 s \n ' ,...

63 cas_name , Np_c , Nc_c , lam_c , overshoot , tr_str );
64

65 % Figure par cas

66 figure ( ' Name ' ,[ ' Cas : ' cas_name ], ' Color ' ,' white ', ' Position '
,[150 150 800 450]) ;

9
ENSA Agadir  2025/2026 TP 02  MPC Approche 2

67 subplot (2 ,1 ,1) ;
68 plot (1: N_bf , yr_signal , 'r -- ',' LineWidth ' ,1.5) ; hold on ;
69 plot (1: N_bf , yc , 'b ',' LineWidth ' ,1.4) ;
70 title ( sprintf ( '% s : Np =% d , Nc =% d , lambda =%.2 f ' , cas_name ,
Np_c , Nc_c , lam_c ));
71 legend ( ' y_r (k) ','y(k ) ', ' Location ',' best '); grid on ;
72 ylabel ( ' Amplitude '); xlabel ( 'k ');
73 subplot (2 ,1 ,2) ;
74 plot (1: N_bf , yc - yr_signal , 'm ' ,' LineWidth ' ,1.2) ;
75 ylabel ( 'e(k ) ') ; xlabel ( 'k ') ; title ( ' Erreur e (k) '); grid on ;
76 end
77 fprintf ( '\n === FIN DU TP 02 ===\ n ');

Listing 6: Boucle de simulation pour les 6 cas paramétriques

Interprétation
Tableau synthétique des eets paramétriques attendus :
Cas Np Nc λ Eet attendu
Nominal 4 2 0.50 Référence de base

Cas 1 12 2 2.00 Np ↑ + λ ↑ ⇒ réponse plus lente, moins de dépassement

Cas 2 8 4 0.50 Np ↑ + Nc ↑ ⇒ plus agile, risque de dépassement

Cas 3 5 4 0.05 λ ↓ ⇒ commande plus agressive, dépassement possible

Cas 4 8 2 0.50 Np ↑ seul ⇒ meilleure anticipation

Cas 5 4 1 5.00 Nc ↓ + λ ↑↑ ⇒ réponse très amortie

Règles qualitatives :
ˆ Np ↑ : l'horizon de prédiction plus long améliore la qualité de suivi mais
augmente le coût de calcul.

ˆ Nc ↑ (avec Nc ≤ Np ) : plus de degrés de liberté ⇒ réponse plus rapide mais

potentiellement plus oscillatoire.

ˆ λ ↑ : régularisation forte ⇒ commande douce, réponse lente ; λ ↓ ⇒ commande

agressive, risque d'instabilité.

ˆ Le meilleur compromis (rapidité, dépassement, stabilité) est typiquement

obtenu pour des valeurs intermédiaires de λ avec Nc < Np .

10
ENSA Agadir  2025/2026 TP 02  MPC Approche 2

2 Conclusion
Ce TP a permis de concevoir et simuler une commande prédictive MPC selon l'Approche
2 (modèle augmenté avec état étendu). Les points clés sont :

1. Le modèle augmenté intègre u(k−1) dans l'état, ce qui permet de formuler directement
u(k) sans passer par les incréments d'état.

2. Les matrices Σi , Φp , Gf , Ψp encodent complètement la dynamique prédite sur

l'horizon Np .

3. La loi de commande est une correction proportionnelle à l'erreur de prédiction,

pondérée par le gain gc⊤ issu de la résolution QP hors-ligne.

4. L'étude paramétrique montre que le réglage Np = 8, Nc = 2, λ = 0,5 (Cas 4) ore

généralement un bon compromis entre rapidité et stabilité.

11