clc, clear all, close all

%% Quesito 1 e 2
% Analisi FdT a c.a. per il loop interno di portata:
G=tf(1,[10 1]);
% essendo un sistema del I° ordine, ed in forma di Bode, abbiamo:
% polo in p=-0.1=> tau=1/|-0.1|=10, K=1, B3=(1/tau)/(2*pi)=0.0159 cioè la pulsazioni 
% corrispondente al punto di rottura divisa per (2*pi), e smorzamento xi=1 in quanto un % modo aperiodico essendo sull’asse reale ha smorzamento xi=asin(pi/2)=1, vedi slide 122. 
% A conferma:
damp(G) % vedi “help damp”
% f.d.t in forma zpk
zpk(G)
% guadagno in alta frequenza: K'=bm/an=10

% Analisi FdT a c.c. per il loop interno di portata:
% C1(s) f.d.t. Regolatore Loop di Portata
C1=tf(100);
W1=feedback(C1*G,1)
% K=0.99, p=-1/10.1, tau=1/|-10.1|=0.099, B3=(1/tau)/(2*pi)=1.6, K'=bm/an=100/10=10,
% A conferma:
W1=zpk(W1);
W1.DisplayFormat = 'frequency' % forma di bode
damp(W1)

% Nota il SDC a c.c. è decisamente più pronto ma ha perso un pelo di accuratezza in 
% termini di risposta al gradino, infatti K=0.99.
% Questo significa: 1% di errore a regime in risposta al gradino. Detto ciò, visto che % il vero obbiettivo del controllo in cascata è controllare la variabile d’uscita 
% dell’anello più esterno questo potrebbe essere un non problema. Vedi soluzione punti % 5 e 6 per dettagli riguardo le prestazioni dell’anello esterno.
figure(1)
subplot(211)
bode(G,W1), grid
subplot(212)
step(G,W1)

%% Quesito 3 e 4

% P(s) f.d.t. Ambiente da condizionare
kp=0.1;
P=tf(kp,[1000 70 1]);
% H2(s) f.d.t. Sensore di Temperatura 
H2=tf(1,[1 1]);
G2=W1*P*H2
zpk(G2)
% La precedente istruzione mostra che il polo di W1 è 10 volte del più veloce tra tutti % i poli e zeri di G2. Per cui in accordo con quanto riportato in slide 133 può essere % trascurato senza compromettere l’analisi.
% Sulla base di considerazioni analoghe si potrebbe pensare sia legittimo trascurare 
% anche il polo di H2. Ciò nonostante, questa approssimazione all’interno della f.d.t. 
% ad anello aperto risulta pericolosa in quanto porterebbe ad una f.d.t. con margini di % stabilità infiniti in presenza di una sola correzione statica, i.e. C2(s)=kc2, cosa 
% questa non verrà. D’altro canto, come verra mostrato nei punti 5 e 6, trascurare W1, % non modifica particolarmente i margini di stabilità del sistema a c.c.
figure(2)
subplot(211)
bode(P,W1*P*H2, P*H2), grid
subplot(212)
step(P,W1*P*H2, P*H2)
% Le figure precedenti confermano che i modi di P sono dominanti, detto ciò con 
% riferimento al progetto dell’anello di controllo esterno, trascurare W1 è legittimo 
% in quanto trascurare tale polo non modifica i margini di stabilità della f.d.t. ad 
% anello aperto, discorso diverso sarebbe se trascurassimo anche H2.
[Gm,Pm,Wcg,Wcp] = margin(W1*P*H2)
[Gm1,Pm1,Wcg1,Wcp2] = margin(P*H2) % Approssimazione 1, legittima
[Gm2,Pm2,Wcg2,Wcp2] = margin(P) 	   % Approssimazione 2, pericolosa 

%% Quesito 5 e 6

% Analisi Loop esterno di Temperatura
% C2(s) f.d.t. Regolatore Loop di Temperatura
kc2=1;
C2=tf(kc2*[100 1],[1 0]);

Wr=feedback(C2*W1*P,H2);
% Wr(s) nell’ipotesi di considerare istantaneo il loop di portata
Wr_approx1=feedback(C2*P,H2); 
% Confronto diagrammi di bode
figure(1)
bode(Wr,Wr_approx1), grid
legend('Wr(s)','Wr_{approx1}(s)');
                              
zeriWr_approx1=zero(Wr_approx1)
poliWr_approx1=pole(Wr_approx1)
% Da ciò è facile che la posizione dei due zeri a c.c. è molto prossima alla posizione % dei due poli, rendendo tali cancellazioni legittime.
% Wr(s) che tiene conto delle cancellazione polo-zero e di W1(s)=1
denWr_approx2=conv([-1/(-0.025-0.098*j) 1],[-1/(-0.025+0.098*j) 1])
% Il comando conv permette di calcolare (1+s/p1)*(1+s/p2)
Wr_approx2=tf(1,denWr_approx2)
 
% Confronto diagrammi di bode
figure(2)
bode(Wr,Wr_approx1,Wr_approx2), grid
legend('Wr(s)','Wr_{approx1}(s)','Wr_{approx2}(s)');
% Confronto risposte indiciali
figure(3)
step(Wr,Wr_approx1,Wr_approx2), grid
legend('Wr(s)','Wr_{approx1}(s)','Wr_{approx2}(s)');
axis([0 250 0 1.5])
 
% Analisi f.d.t. Disturbo/Uscita
Wdt=feedback(1,C2*W1*P*H2);
Wdt_approx1=feedback(1,C2*P*H2);
 
Wdq=1/(C1*G)*feedback(W1*P,C2*P*H2);
Wdq_approx1=1/(C1*G)*feedback(W1*P,C2*P*H2);
 
figure(4)
subplot(1,2,1)
    step(Wdt,Wdt_approx1), grid
    legend('W \Delta \theta(s)','W \Delta \theta_{approx1}(s)')
subplot(1,2,2)
    step(Wdq,Wdq_approx1), grid
    legend('W \Delta q(s)','W \Delta q_{approx1}(s)')

figure(5)
subplot(1,2,1)
    bode(Wdt,Wdt_approx1), grid
    legend('W \Delta \theta(s)','W \Delta \theta_{approx1}(s)')
subplot(1,2,2)
    bode(Wdq,Wdq_approx1), grid
    legend('W \Delta q(s)','W \Delta q_{approx1}(s)')

% Le figure mostrano che il legame setpoint-uscita è di tipo passa basso con 
% B3=0.15/(2*pi) e Bp=[0, B3]. La f.d.t. disturbo uscita Wdt è invece di tipo passa 
% alto con B3=0.06/(2*pi) e Bp=[B3, +infinito]. Mentre Wdq è di tipo passa banda, ma 
% con una attenuazione davvero importante e Bp=[0.01/(2*pi), 0.02/(2*pi)].

%% Quesito 7 e 8

% f.d.t. ad anello aperto per il loop di temperatura 
F=C2*W1*P*H2;
Fapprox1=C2*P*H2;
 
figure(6)
subplot(1,2,1), 
    margin(F), legend('F(j\omega)')
subplot(1,2,2), 
    margin(Fapprox1), legend('F_{approx1}(j\omega)')  

kcr=10^(16.1/20);

figure(7)
subplot(1,2,1),
    nyquist(Fapprox1,kcr*Fapprox1,10*Fapprox1)
    legend('F_{approx1}(j\omega)| k_{C2}=1',...
        'F_{approx1}(j\omega)| k_{C2}=6.38','F_{approx1}(j\omega)| k_{C2}=10')
subplot(1,2,2), 
    nyquist(Fapprox1,6.38*Fapprox1,10*Fapprox1)
    axis([-2 1 -0.3 0.3])    

figure(8)
    nyquist(F,Fapprox1)
    axis([-1e-4 5e-5 -1e-3 1e-3])
    legend('F(j\omega)','F_{approx1}(j\omega)')

[Gm,Pm,Wcg,Wcp] = margin(F)
[Gm_approx1,Pm_approx1,Wcg_approx1,Wcp_approx1] = margin(Fapprox1)


%% Quesito 10   
G11=tf(-1,[-10 1]);
Fp=G11;
figure(4)
    nyquist(0.5*Fp,Fp,2*Fp)
    legend('F_p(j\omega)| k_{C1}=0.5',...
           'F_p(j\omega)| k_{C1}=1','F_p(j\omega)| k_{C1}=2')
% Questo sistima è un tipico esempio di sistema a stabilità paradossale.
