feat: Add problem 3

.gitattributes

@@ -0,0 +1,28 @@
+* text=auto
+
+*.fig binary
+*.mat binary
+*.mdl binary diff merge=mlAutoMerge
+*.mdlp binary
+*.mex* binary
+*.mlapp binary
+*.mldatx binary
+*.mlproj binary
+*.mlx binary
+*.p binary
+*.sfx binary
+*.sldd binary
+*.slreqx binary merge=mlAutoMerge
+*.slmx binary merge=mlAutoMerge
+*.sltx binary
+*.slxc binary
+*.slx binary merge=mlAutoMerge
+*.slxp binary
+
+## Other common binary file types
+*.docx binary
+*.exe binary
+*.jpg binary
+*.pdf binary
+*.png binary
+*.xlsx binary

Lab03.m

@@ -11,9 +11,9 @@ 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)];
+f = [-p1*(G-Gb)-beta*G+gamma/VG, -n*I + r/V1, -p2*beta + p3*(I-Ib)];
 g = G;
-x = [G, beta, I];
+x = [G, I, beta];
 u = [r gamma];
 
 %% Calcolo punto di equilibrio
@@ -28,7 +28,7 @@ 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 ];
+x_eq = [ G_eq I_eq beta_eq ];
 % Osserviamo che G_eq corrisponde al valore corretto di glicemia (81)
 
 %%% Soluzione con solver simbolico
@@ -36,7 +36,7 @@ X_eq = [ G_eq I_eq beta_eq ];
 % 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 ])
+x_eq = double([ G_eq I_eq beta_eq ])
 
 % Il risultato calcolato da MATLAB è identico a quello trovato a mano
 
@@ -50,13 +50,13 @@ X_eq = double([ G_eq I_eq beta_eq ])
 % corrispondente valore nel punto di equilibrio
 
 A = jacobian(f,x);
-A = double(subs(A,x,X_eq));
+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));
+C = double(subs(C,x,x_eq));
 
 D = jacobian(g,u);
 D = double(subs(D,u,u_eq));
@@ -67,10 +67,24 @@ 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
+% Il sistema è stabile perché tutti gli autovalori sono negativi
 
 %%% Stabilità BIBO
 H = C*W*B;
 p = pole(H)
 
-% Il sistema non è stabile BIBO in quanto è presente un polo positivo
+% Il sistema + stabile BIBO in quanto tutti i poli sono negativi
+
+%% Studio raggiungibilità e controllabilità
+M_R = ctrb(A,B);
+rank(M_R)
+% La matrice è di rango 3, quindi il sistema è completamente raggiungibile
+
+%% Progettazione della legge di controllo
+lambda_k = [-1 -2 -3];
+
+% Calcolo K
+K = place(A,B(:,1),lambda_k);
+
+% Non serve calcolare alfa in quanto il controllore serve solo a far si che
+% il sistema non si sposti dallo stato di equilibrio

Lab04.m

@@ -0,0 +1,84 @@
+%% Setup
+% Inizializzazione costanti
+m = 0.02;
+G = 9.81;
+Kt = 708.27;
+Km = 1.52*10^-4;
+
+% Inizializzazione solver simbolico con funzioni del sistema
+syms x1 x2 Im
+f = [x2, G-(Km.*Im.^2)/(m.*x1.^2)];
+x = [x1, x2];
+u = Im;
+
+%% Calcolo punto di equilibrio
+% L'ingresso del sistema è dato dal testo
+u_eq = 0.8;
+
+% La soluzione si può sia calcolare a mano che usando il solver simbolico
+% di matlab
+[x1_eq, x2_eq] = solve(subs(f,u,u_eq)==[0 0]);
+% Il comando double calcola il valore delle frazioni
+x_eq = double([ x1_eq, x2_eq ]);
+
+% Si osserva che sono presenti due punti di equilibrio. Prendiamo in
+% considerazione solo quello con variabili positive
+x_eq = x_eq(2,:)
+
+%% 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
+
+% Definisco la funzione uscita attorno al punto di equilibrio
+g = Kt.*(x1-x_eq(1));
+
+A = jacobian(f,x);
+A = double(subs(A,[x u],[x_eq u_eq]));
+
+B = jacobian(f,u);
+B = double(subs(B,[x u],[x_eq 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 perché è presente un autovalore positivo
+
+%%% Stabilità BIBO
+H = C*W*B;
+p = pole(H)
+
+% Il sistema è instabile BIBO è presente un polo positivo
+
+%% Studio raggiungibilità e osservabilità
+% Raggiungibilità
+M_R = ctrb(A,B);
+rank(M_R)
+% Il rango è due, quindi il sistema è completamente raggiungibile
+
+% Osservabilità
+M_O = obsv(A,C);
+rank(M_O)
+% Il rango è due, quindi il sistema è completamente osservabile
+
+%% Progettazione sistema di controllo
+% Scelta autovalori
+lambda_k = [-0.7 -0.8];
+lambda_o = [-10 -11];
+
+% Calcolo K
+K = place(A,B,lambda_k);
+
+% Calcolo L
+L = place(A',C',lambda_o)';
+
+% Calcolo alpha
+alpha = inv(-C*((A-B*K)\B));

Lab05.m

@@ -0,0 +1,57 @@
+%% Setup variabili
+A = [0 1; -3 -4];
+B = [0 ; 1];
+C = [2 1];
+D = 0;
+
+%% Punto 1-2
+lambda_k = [-0.7 -0.8];
+lambda_o = [-10 -11];
+
+% Controllo raggiungibilità
+M_R = ctrb(A,B);
+rank(M_R)
+% Rango 2: completamente raggiungibile
+
+% Controllo osservabilità
+M_O = obsv(A,C);
+rank(M_O)
+% Rango 2: completamente osservabile
+
+% Calcolo K
+K = place(A,B,lambda_k);
+
+% Calcolo L
+L = place(A',C',lambda_o)';
+
+% Calcolo alpha
+alpha = inv(-C*((A-B*K)\B));
+
+%% Punto 3
+% (Su simulink) Si osserva che l'uscita del sistema tende ad 1
+
+%% Punto 4
+lambda_k_4 = [-10 -12];
+K_4 = place(A,B,lambda_k_4);
+alpha_4 = inv(-C*((A-B*K_4)\B));
+% Si osserva che con autovalori troppo elevati il sistema converge lo stesso, ma
+% all'inizio si osserva un transitorio con valori troppo grandi (più sono grandi gli autovalori, più
+% l'errore cresce)
+
+% Con questa tecnica non riesco a controllare il sistema velocemente per
+% via dello zero nella H(s)
+minreal(zpk(ss(A,B,C,D)))
+
+%% Punto 5
+% Nota: ripristinare gli autovalori del punto 1.
+epsilon = 10.^(-2/20);
+B_reale = B.*epsilon;
+% Osserviamo che il sistema "reale" (con la perturbazione) è instabile.
+
+%% Punto 6
+lambda_o_6 = [-1000 -1100];
+L_6 = place(A',C',lambda_o_6)';
+
+%% Punto 7
+lambda_o_7 = [-1000 -2];
+L_7 = place(A',C',lambda_o_7)';

Prob1.m

@@ -0,0 +1,67 @@
+%% Dati
+s = tf('s');
+Gp = 25/(s^3 + 3.3*s^2 + 2*s);
+Gs = 1;
+Ga = 0.095;
+Gr = 1;
+Gd = 1;
+
+% Specifiche
+Kd = 1;
+
+%% Calcolo Kc da errore di inseguimento
+% L'input di riferimento è una rampa
+% Le formule sono state prese dalle tabelle sulle slide
+rho_r = 0.15;
+Kp = 25/2;
+Kc_r = 1/(rho_r*Kp*Ga) % Kc deve essere maggiore di questo
+
+%% Calcolo Kc da disturbo su attuatore
+% L'input di riferimento è un gradino
+% La formula si ottiene calcolando il limite manualmente
+rho_a = 0.015;
+Da0 = 0.0055;
+Kc_a = Da0/(Ga*rho_a)
+
+%% Calcolo MFLS e omega_c da disturbo dell'impianto
+% Il disturbo è di tipo sinusoidale
+rho_p = 0.0005;
+ap = 0.02;
+omega_p = 0.02;
+MLFS = mag2db(rho_p/ap)
+omega_l = omega_p*10.^(-MLFS/40);
+omega_c_inf = 2*omega_l
+
+%% Calcolo MHFT e omega_c da disturbo del sensore
+% Il disturbo è di tipo sinusoidale
+rho_s = 0.0005;
+as = 0.1;
+omega_s = 40;
+MHFT = mag2db((rho_s*Gs)/as)
+omega_h = omega_s*10.^(MHFT/40);
+omega_c_inf = omega_h/2
+
+%% Calcolo picco di risonanza e omega_c da tempo di salita e di assestamento
+scap = 0.1;
+tr = 3;
+ts = 12;
+alpha = 0.05;
+xi = abs(log(scap))/sqrt(pi.^2+log(scap).^2)
+Tp = 1/(2*xi*sqrt(1-xi.^2))
+Sp = (2*xi*sqrt(2+4*xi^2+2*sqrt(1+8*xi^2)))/(sqrt(1+8*xi^2)+4*xi^2-1)
+omega_c_rise = (((pi-acos(xi))*sqrt(sqrt(1+4*xi^4)-2*xi^2))/sqrt(1-xi^2))/tr
+omega_c_settle = (-log(alpha)*sqrt(sqrt(1+4*xi^4)-2*xi^2))/(xi*ts)
+
+% Diagramma di Nichols e Nyquist
+myngridst(Tp,Sp)
+
+%% Definizione di Gc
+% Determinazione del modulo di Kc
+Kc = max(Kc_a, Kc_r);
+
+% Analisi del diagramma di Nyquist per la determinazione del segno di Kc
+L = -Kc*Gp*Ga;
+[num, den] = tfdata(L, "v");
+figure
+nyquist1(num, den)
+grid on

Prob2.m

@@ -0,0 +1,71 @@
+%% Dati
+s = tf('s');
+Gp = 40/(s^2 + 3*s + 4.5);
+Gs = 1;
+Ga = -0.09;
+Gr = 1;
+Gd = 1;
+
+% Specifiche
+Kd = 1;
+p = 0; % L'impianto non ha poli nell'origine
+
+%% Calcolo Kc da errore di inseguimento
+% L'input di riferimento è una rampa
+% Le formule sono state prese dalle tabelle sulle slide
+
+h_r = 1; % input è una rampa
+nu_r = h_r - p % limite inferiore per nu_r
+
+rho_r = 0.35;
+Kp = evalfr(Gp, 0);
+Kc_r = -1/(rho_r*Kp*Ga) % Kc deve essere maggiore di questo
+
+%% Calcolo Kc da disturbo su attuatore
+% Il disturbo è un gradino
+h_a = 0;
+nu_a = h_a - p
+% da questo calcolo risulta nu>=0, ma considerando che il punto precedente 
+% impone nu >= 1, allora risulta calcolando il limite che il contributo del
+% disturbo tende a zero, pertanto non possiamo ricavare nessuna condizione
+% su Kc da questo requisito
+
+%% Calcolo Kc da disturbo dell'impianto
+% Il disturbo è una rampa
+rho_p = 0.001;
+Dp0 = 0.003;
+Kc_p = -Dp0/(rho_p*Kp*Ga)
+
+%% Calcolo MHFT e omega_c da disturbo del sensore
+% Il disturbo è di tipo sinusoidale
+rho_s = 0.0002;
+as = 0.01;
+omega_s = 50;
+MHFT = mag2db((rho_s*Gs)/as)
+omega_h = omega_s*10.^(MHFT/40);
+omega_c_inf = omega_h/2
+
+%% Calcolo picco di risonanza e omega_c da tempo di salita e di assestamento
+scap = 0.08;
+tr = 2.5;
+ts = 10;
+alpha = 0.05;
+xi = abs(log(scap))/sqrt(pi.^2+log(scap).^2)
+Tp = 1/(2*xi*sqrt(1-xi.^2))
+Sp = (2*xi*sqrt(2+4*xi^2+2*sqrt(1+8*xi^2)))/(sqrt(1+8*xi^2)+4*xi^2-1)
+omega_c_rise = (((pi-acos(xi))*sqrt(sqrt(1+4*xi^4)-2*xi^2))/sqrt(1-xi^2))/tr
+omega_c_settle = (-log(alpha)*sqrt(sqrt(1+4*xi^4)-2*xi^2))/(xi*ts)
+
+% Diagramma di Nichols e Nyquist
+myngridst(Tp,Sp)
+
+%% Definizione di Gc
+% Determinazione del modulo di Kc
+Kc = max(Kc_a, Kc_p);
+
+% Analisi del diagramma di Nyquist per la determinazione del segno di Kc
+L = -Kc*Gp*Ga;
+[num, den] = tfdata(L, "v");
+figure
+nyquist1(num, den)
+grid on

Prob3.m

@@ -0,0 +1,74 @@
+%% Dati
+s = tf('s');
+Gp = 100/(s^2 + 5.5*s + 4.5);
+Gs = 1;
+Ga = 0.014;
+Gr = 1;
+Gd = 1;
+
+% Specifiche
+Kd = 1;
+p = 0; % L'impianto non ha poli nell'origine
+
+%% Calcolo Kc da errore di inseguimento
+% L'input di riferimento è una rampa
+% Le formule sono state prese dalle tabelle sulle slide
+
+h_r = 1; % input è una rampa
+nu_r = h_r - p % limite inferiore per nu_r
+
+rho_r = 0.15;
+Kp = evalfr(Gp, 0);
+Kc_r = 1/(rho_r*Kp*Ga) % Kc deve essere maggiore di questo
+
+%% Calcolo Kc da disturbo su attuatore
+% Il disturbo è un gradino
+h_a = 0;
+nu_a = h_a - p
+% da questo calcolo risulta nu>=0, ma considerando che il punto precedente 
+% impone nu >= 1, allora risulta calcolando il limite che il contributo del
+% disturbo tende a zero, pertanto non possiamo ricavare nessuna condizione
+% su Kc da questo requisito
+
+%% Calcolo MFLS e omega_c da disturbo dell'impianto
+% Il disturbo è di tipo sinusoidale
+rho_p = 0.002;
+ap = 0.16;
+omega_p = 0.03;
+MLFS = mag2db(rho_p/ap)
+omega_l = omega_p*10.^(-MLFS/40);
+omega_c_inf = 2*omega_l
+
+%% Calcolo MHFT e omega_c da disturbo del sensore
+% Il disturbo è di tipo sinusoidale
+rho_s = 0.0008;
+as = 0.2;
+omega_s = 60;
+MHFT = mag2db((rho_s*Gs)/as)
+omega_h = omega_s*10.^(MHFT/40);
+omega_c_inf = omega_h/2
+
+%% Calcolo picco di risonanza e omega_c da tempo di salita e di assestamento
+scap = 0.12;
+tr = 2;
+ts = 8;
+alpha = 0.05;
+xi = abs(log(scap))/sqrt(pi.^2+log(scap).^2)
+Tp = 1/(2*xi*sqrt(1-xi.^2))
+Sp = (2*xi*sqrt(2+4*xi^2+2*sqrt(1+8*xi^2)))/(sqrt(1+8*xi^2)+4*xi^2-1)
+omega_c_rise = (((pi-acos(xi))*sqrt(sqrt(1+4*xi^4)-2*xi^2))/sqrt(1-xi^2))/tr
+omega_c_settle = (-log(alpha)*sqrt(sqrt(1+4*xi^4)-2*xi^2))/(xi*ts)
+
+% Diagramma di Nichols e Nyquist
+myngridst(Tp,Sp)
+
+%% Definizione di Gc
+% Determinazione del modulo di Kc
+Kc = max(Kc_r);
+
+% Analisi del diagramma di Nyquist per la determinazione del segno di Kc
+L = -Kc*Gp*Ga;
+[num, den] = tfdata(L, "v");
+figure
+nyquist1(num, den)
+grid on