Résolution de systèmes d'équations différentielles

Bonjour à tous . Je veux résoudre ce système d'équations différentielles où il y a une dépendance entre les deux variables d'état.
Voici le système en question:

dy(1)= (Achap+Bchap*Ft*Kt)*y(1) + Bchap*lda*dy(2);
dy(2)= Attild*y(2) + Bttild*u0

et on veut tracer la loi de commande u de la forme :
u = +Ft*Kt*y(1)+lda*dy(2);

Sachant que :
- Les variables d'état sont y(1) et y(2) , tout le reste est connu.
- dy(1) est de dimension 20*1 , dy(2) est de dimension 10*1 .
- Les autres matrices sont constantes (j'ai bien vérifié les dimensions de ces matrices).

Voila ça fait qqs semaines que je ne trouve pas de solution alors j'espère que quelqu'un ici pourra m'aider!

Merci d'avance.

code matlab de ce programme :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
function  dy = ilmi4tankdelay(t,y)
%Conditions initiales des variables d'état 
 
%dy = zeros(30,1) ;  
dy=[zeros(20,1);zeros(10,1)];
 
 
%%%%%%%%%%%%%%% quadruple tank process with delays %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Définition des paramètres du système
S1=28; S2=32 ; S3=28 ; S4=32 ;
a1=0.071; a2=0.057; a3=0.071; a4=0.057 ;
g= 981 ;  kc =0.5;
 
 
% Définition de paramètres pour un système à minimum de phase (zéros stables)
h10 = 12.4 ; h20 = 12.7; 
h30 = 1.8 ; h40 = 1.4; 
v10 = 3 ; v20 = 3 ;
k1=3.33 ; k2 = 3.35 ; 
gamma1 = 0.7 ; gamma2 = 0.6  ;
 
 
%%%%%%%%%%%%%
% système original : Quadruple tank with delay
A0= [-a1/S1*sqrt(g/2*h10) 0 0 0;0 -a2/S2*sqrt(g/2*h20) 0 0;0 0 -a3/S3*sqrt(g/2*h30) 0;0 0 0 -a4/S4*sqrt(g/2*h40)];
A1=[0 0 a3/S1*sqrt(g/2*h30) 0;0 0 0 a4/S2*sqrt(g/2*h40);0 0 0 0;0 0 0 0];
B0=[(k1*gamma1)/S1 0; 0 (k2*gamma2)/S2; 0 0; 0 0];
B1=[0 0; 0 0;0 (1-gamma2)*k2/S3;(1-gamma1)*k1/S4 0];
C=[kc 0 0 0;0 kc 0 0];
 
tau1=5 ;
tau2=2 ;
tau3=4 ;
 
% tau1=5 ;
% tau2=2 ;
% tau3=4 ;
 
%Système augmenté 1
Atild=[A0 A1 zeros(4,2);zeros(4,10);C zeros(2,4) zeros(2,2)];
Btild=[B0 B1; zeros(4,2) zeros(4,2);zeros(2,2) zeros(2,2)];
Ctild=[C zeros(2,4) zeros(2,2)]; 
 
 
%Paramètres de la méthode de collocation orthogonale
A1t = [6.6640    2.8405   -1.2325    1.1613   -0.6550
     -5.1848    0.7688    2.9413   -2.2497    1.2291
     2.2497   -2.9413   -0.7688    5.1848   -2.4952
    -1.1613    1.2325   -2.8405   -6.6640    8.7783
     1.7631   -1.8126    3.6799  -23.6304   21.0000];
B1t = [-8.7783 -8.7783
     2.4952 2.4952
     -1.2291 -1.2291 
     0.6550 0.6550
     -1.0000 -1.0000];
C1t=[0 0 0 0 1;0 0 0 0 1];
 
A2t = [6.6640    2.8405   -1.2325    1.1613   -0.6550
   -5.1848    0.7688    2.9413   -2.2497    1.2291
    2.2497   -2.9413   -0.7688    5.1848   -2.4952
   -1.1613    1.2325   -2.8405   -6.6640    8.7783
    1.7631   -1.8126    3.6799  -23.6304   21.0000];
B2t = [-8.7783 -8.7783
     2.4952 2.4952
     -1.2291 -1.2291 
     0.6550 0.6550
     -1.0000 -1.0000];
C2t=[0 0 0 0 1;0 0 0 0 1];
 
%u11=C1t*w1;
%u12=C2t*w2;
 
At=[A1t zeros(5,5);zeros(5,5) A2t];
Bt=[B1t;B2t];
Ct=[C1t zeros(2,5);zeros(2,5) C2t];
 
%Système augmenté 2
Attild=[(1/tau2)*A1t zeros(5,5);zeros(5,5) (1/tau3)*A2t];
Bttild=[(1/tau2)*B1t;(1/tau3)*B2t];
Achap=[Atild Btild*Ct; zeros(10,10) Attild];
Bchap=[zeros(10,2) ; Bttild];
%Cchap=[Ctild zeros(2,10)];
 
 
 
 
%%%%%%%%%%%%
% % Définition de paramètres pour un système à non-minimum de phase (zéros instables sur le demi-plan gauche)
% h10 = 12.6 ;h20 = 13; 
% h30 = 4.8 ;h40 = 4.9; 
% v10 = 3.15 ; v20 = 3.15 ;
% k1=3.14 ; k2 = 3.29 ; 
% gamma1 = 0.43 ; gamma2 = 0.34 ;
 
K1t=[Ctild zeros(2,10)];
K2t=[zeros(2,8) eye(2,2) zeros(2,10)];
K3t=[Ctild*Atild zeros(2,10)];
 
% At= [A zeros(4,2); C zeros(2,2)];
% Bt = [B ;zeros(2,2)];
Kt = [K1t' K2t' K3t']';
 
Ft = 1.0e+009 * [-5.7619   -5.8320   -3.6609   -3.7854   -1.4857   -0.1871 ; -5.7619   -5.8320   -3.6609   -3.7854   -1.4857   -0.1871];
 
 
F1 = [3.1559    0.3655
    0.2936    3.1942];
 
F2 = [0.2196    0.0838
    0.0708    0.1541];
 
 
F3 = [24.3199    0.0646
    0.0208   32.2304];
 
 
%u = +Ft*Kt*y(1) +lambda*dy(2); % Retour de sortie dynamique 
% 
% dksi = (Achap+Bchap*Fi*Kt)*ksi + Bchap*lambda*domega;
% domega = Attild*omega + Bttild*u
u0=[3;3];
lda = -F3*Ctild*Btild*Ct*inv(Attild) ;
 
%dy=[dy(1);dy(2)];
%[dy]= [(Achap+Bchap*Ft*Kt)*y(1) + Bchap*lda*dy(2); Attild*y(2) + Bttild*u0 ];
dy(1)= (Achap+Bchap*Ft*Kt)*y(1) + Bchap*lda*dy(2); 
dy(2)= Attild*y(2) + Bttild*u0 ;

%%%%%%%%%%%%%%%%%%Simulation hors ode45%%%%%%%%%%%%%%%%%%%%%

Code:

1
2
3
4
5
6
7
8
9
10
11
%commande workspace de la fonction ode 45
options = odeset('RelTol',1e-5,'AbsTol',[1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7 1e-7]);
[t,y] = ode45(@ilmi4tankdelay,[0 50],[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0],options);
 
figure(1)
plot(t,+Ft*Kt*y(1))
xlabel ('temps (s) '), ylabel('Loi de commande (V) '), grid ;
 
figure(1)
plot(t,+Ft*Kt*y(1)+lda*dy(2))
xlabel ('temps (s) '), ylabel('Loi de commande (V) '), grid ;