diff --git a/Lab03.m b/Lab03.m new file mode 100644 index 0000000..aa03913 --- /dev/null +++ b/Lab03.m @@ -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 \ No newline at end of file