% rif. "Progetto in frequenza: esercizio riassuntivo"
% specifiche a), b) e c)

close all;clear all;

% processo P = 10/(1+0.1s) (1+0.01s)
 
numP=10;
denP=conv([0.1 1],[0.01 1]);

% processo modificato Fhat(s) 
% Fhat(s) = 100/s(1+0.1s)(1+0.01s)
% ottenuto dal processo P(s)
% con l'aggiunta di un polo nell'origine e opportuno guadagno
% la pulsazione di attraversamento desiderata e' w*_t=30
% il margine di fase desiderato e' m*_fi=20

numFhat=conv([10],numP);
denFhat=conv([1,0],denP);

% diagrammi di Bode Fhat(s)

figure(1)
bode(numFhat,denFhat)
pause

% calcolo valori attuali pulsazione di attraversamento e margine di fase 

figure(2)
margin(numFhat,denFhat)
pause

% anello chiuso senza compensazione
% sistema "quasi" al limite di stabilita'
 
[numWhat,denWhat] = cloop(numFhat,denFhat,-1);
roots(denWhat)
pause

w1=logspace(0,2,100);
[magWhat,phaseWhat]=bode(numWhat,denWhat,w1);
magWdBhat=20*log10(magWhat);

figure(3)
semilogx(w1,magWdBhat,'-',w1,-3+0*w1),grid
xlabel('pulsazione (rad/sec)')
ylabel('modulo (dB)')
title('risposta armonica sistema retroazionato senza compensazione')
pause

% risposta indiciale: notare sovraelongazione e oscillazioni 

t=0:0.01:2;
figure(4)
step(numWhat,denWhat,t)
title('risposta indiciale sistema retroazionato senza compensazione')
pause

% per ottenere un anticipo di fase di circa 20 gradi lasciando
% quasi inalterata la pulsazione di attraversamento si sceglie m_a=10
% e w tau_a =0.4 (con w=w*_t=30) dai diagrammi universali

tau_a=0.4/30;
m_a=10;
numR=[tau_a, 1];
denR=[tau_a/m_a, 1];

figure(5)
bode(numR,denR)
pause

% verifica soddisfacimento specifiche sui 
% diagrammi di Bode di F(s)=R(s)Fhat(s)=G(s)P(s)

numF=conv(numFhat,numR);
denF=conv(denFhat,denR);

w1=logspace(0,4,1000);
[magFhat,phaseFhat]=bode(numFhat,denFhat,w1);
[magF,phaseF]=bode(numF,denF,w1);
magFhatdB=20*log10(magFhat);
magFdB=20*log10(magF);
figure(6)
subplot(211)
semilogx(w1,magFhatdB,'r-',w1,magFdB,'b-','LineWidth',2),grid
xlabel('Frequency (rad/sec)')
ylabel('Magnitude (dB)')
title({'risposta armonica ad anello aperto','prima della compensazione (rosso) e dopo la compensazione (blu)'})
subplot(212)
semilogx(w1,phaseFhat,'r-',w1,phaseF,'b-',w1,-180*(1+0*w1),'w-','LineWidth',2),grid
xlabel('Frequency (rad/sec)')
ylabel('Phase (deg)')
pause

% zoom diagrammi precedenti nella zona di interesse

w2=logspace(1,2,100);
[magFhat,phaseFhat]=bode(numFhat,denFhat,w2);
[magF,phaseF]=bode(numF,denF,w2);
magFhatdB=20*log10(magFhat);
magFdB=20*log10(magF);
figure(7)
subplot(211)
semilogx(w2,magFhatdB,'r-',w2,magFdB,'b-','LineWidth',2),grid
xlabel('Frequency (rad/sec)')
ylabel('Magnitude (dB)')
title({'risposta armonica ad anello aperto','prima della compensazione (rosso) e dopo la compensazione (blu)'})
subplot(212)
semilogx(w2,phaseFhat,'r-',w2,phaseF,'b-',w2,-180*(1+0*w2),'w-','LineWidth',2),grid
xlabel('Frequency (rad/sec)')
ylabel('Phase (deg)')
pause

% calcolo valori pulsazione di attraversamento e margine di fase 
% dopo la compensazione

figure(8)
margin(numF,denF)
pause

% diagrammi di Bode della funzione di trasferimento 
% ad anello chiuso W(s)=F(s)/1+F(s)...

[numW,denW] = cloop(numF,denF,-1);

w3=logspace(0,2,100);
[magW,phaseW]=bode(numW,denW,w3);
magWdB=20*log10(magW);

figure(9)
semilogx(w3,magWdB,'-',w3,-3+0*w3),grid
xlabel('pulsazione (rad/sec)')
ylabel('modulo (dB)')
title('risposta armonica sistema retroazionato con compensazione')
pause

clsys=tf(numW,denW);

bandwidth(clsys)                    % banda passante finale
20*log10(getPeakGain(clsys))        % modulo alla risonanza finale

% risposta indiciale finale 
figure(10)
step(clsys,t)
title('risposta indiciale sistema retroazionato con compensazione')

