feat: Add first three sections of third lab
This commit is contained in:
76
Lab03.m
Normal file
76
Lab03.m
Normal file
@@ -0,0 +1,76 @@
|
|||||||
|
%% Setup
|
||||||
|
% Inizializzazione costanti
|
||||||
|
p1 = 0.003;
|
||||||
|
p2 = 0.025;
|
||||||
|
p3 = 0.000013;
|
||||||
|
V1 = 12;
|
||||||
|
VG = 126;
|
||||||
|
n = 5/54;
|
||||||
|
Gb = 81;
|
||||||
|
Ib = 15;
|
||||||
|
|
||||||
|
% Inizializzazione solver simbolico con funzioni del sistema
|
||||||
|
syms G beta I gamma r
|
||||||
|
f = [-p1*(G-Gb)-beta*G+gamma/VG, -n*I + r/V1, -p2*beta+ p3*(I-Ib)];
|
||||||
|
g = G;
|
||||||
|
x = [G, beta, I];
|
||||||
|
u = [r gamma];
|
||||||
|
|
||||||
|
%% Calcolo punto di equilibrio
|
||||||
|
% L'ingresso del sistema è dato dal testo
|
||||||
|
r_eq = 16.66667;
|
||||||
|
gamma_eq = 0;
|
||||||
|
u_eq = [r_eq gamma_eq];
|
||||||
|
|
||||||
|
%%% Soluzione a mano
|
||||||
|
% Ponendo il sistema di equazioni a zero e risolvendo a mano risulta
|
||||||
|
I_eq = r_eq/(n*V1);
|
||||||
|
beta_eq = p3/p2 .* (I_eq - Ib);
|
||||||
|
G_eq = (gamma_eq/VG + p1*Gb)/(p1+beta_eq);
|
||||||
|
|
||||||
|
X_eq = [ G_eq I_eq beta_eq ];
|
||||||
|
% Osserviamo che G_eq corrisponde al valore corretto di glicemia (81)
|
||||||
|
|
||||||
|
%%% Soluzione con solver simbolico
|
||||||
|
% Essendo il sistema non-lineare, è necessario utilizzare il motore
|
||||||
|
% simbolico per risolvere il sistema di equazioni
|
||||||
|
[G_eq, I_eq, beta_eq] = solve(subs(f,u,u_eq)==[0 0 0]);
|
||||||
|
% Il comando double calcola il valore delle frazioni
|
||||||
|
X_eq = double([ G_eq I_eq beta_eq ])
|
||||||
|
|
||||||
|
% Il risultato calcolato da MATLAB è identico a quello trovato a mano
|
||||||
|
|
||||||
|
%% Linearizzazione del sistema
|
||||||
|
% Lo jacobiano si può calcolare sia a mano che utilizzando il motore
|
||||||
|
% simbolico di MATLAB. La seguente soluzione implementa il secondo metodo
|
||||||
|
|
||||||
|
% Per ogni matrice si calcola prima lo jacobiano corrispondente e poi si
|
||||||
|
% effettua la sostituzione utilizzando il comando subs, che ci permette di
|
||||||
|
% sostituire le variabili di stato / l'input del sistema con il
|
||||||
|
% corrispondente valore nel punto di equilibrio
|
||||||
|
|
||||||
|
A = jacobian(f,x);
|
||||||
|
A = double(subs(A,x,X_eq));
|
||||||
|
|
||||||
|
B = jacobian(f,u);
|
||||||
|
B = double(subs(B,u,u_eq));
|
||||||
|
|
||||||
|
C = jacobian(g,x);
|
||||||
|
C = double(subs(C,x,X_eq));
|
||||||
|
|
||||||
|
D = jacobian(g,u);
|
||||||
|
D = double(subs(D,u,u_eq));
|
||||||
|
|
||||||
|
%% Studio stabilità
|
||||||
|
s = tf('s');
|
||||||
|
W = minreal(zpk(inv(s*eye(size(A))-A)));
|
||||||
|
|
||||||
|
%%% Stabilità interna
|
||||||
|
E = real(eig(A))
|
||||||
|
% Il sistema è instabile in quanto è presente un autovalore maggiore di zero
|
||||||
|
|
||||||
|
%%% Stabilità BIBO
|
||||||
|
H = C*W*B;
|
||||||
|
p = pole(H)
|
||||||
|
|
||||||
|
% Il sistema non è stabile BIBO in quanto è presente un polo positivo
|
||||||
Reference in New Issue
Block a user