feat: Add problem 1

Lab03.m

@@ -73,4 +73,18 @@ E = real(eig(A))
 H = C*W*B;
 p = pole(H)
 
-% Il sistema + stabile BIBO in quanto tutti i poli sono negativi
+% 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