% Sintesi per tentativi
% Esempio di uso combinato di funzioni anticipatrici e
% attenuatrici
% cfr Lanari,Oriolo "Controlli Automatici: Esercizi di Sintesi"
% es. 2.3 p. 79
% versione 19/4/99
 
close all;clear all;

% processo modificato Fhat(s) ottenuto dal processo P(s)
% con l'aggiunta di un polo nell'origine e opportuno guadagno
% la pulsazione di attraversamento desiderata � w*_t=4
% il margine di fase desiderato � m*_fi=50

numFhat=2500*[1 1];
denFhat=conv([1 5 0],[1 12 20]);

% 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
close all

% calcolo modulo e fase di Fhat in corrispondenza alla 
% pulsazione di attraversamento desiderata 

w=4;
[magw,phasew]=bode(numFhat,denFhat,w);
magwdB=20*log10(magw)
phasew
pause

% in corrispondenza a w*_t=4 e' necessario ottenere simultaneamente   
% un anticipo di fase di circa 8 gradi e un'attenuazione di 18 dB
% allo scopo, si puo' usare una rete a sella (attenuatrice in
% bassa frequenza + anticipatrice in alta frequenza)

% primo passo: scelta rete anticipatrice
% in previsione della successiva introduzione della rete 
% attenuatrice, conviene ottenere un anticipo di fase maggiore
% pari a 8+10=18 gradi

% per ottenere un anticipo di fase di circa 18 gradi aumentando
% il meno possibile la pulsazione di attraversamento si sceglie m_a=2
% e w tau_a =1 (con w=w*_t=4) dai diagrammi universali

tau_a=1/4;
m_a=2;
numR_a=[tau_a, 1];
denR_a=[tau_a/m_a, 1];

% calcolo modulo e fase del processo intermedio Fint(s)= R_a(s) Fhat(s) 
% in corrispondenza alla pulsazione di attraversamento desiderata 

numFint=conv(numFhat,numR_a);
denFint=conv(denFhat,denR_a);

w=4;
[magw,phasew]=bode(numFint,denFint,w);
magwdB=20*log10(magw)
phasew
pause

% secondo passo: scelta rete attenuatrice
% per ottenere una attenuazione di circa magwdB=20 dB contenendo entro
% 10 gradi il ritardo di fase si sceglie m_i=10
% e w tau_i=100 (con w=w*_t=4) dai diagrammi universali

tau_i=100/4;
m_i=10;
numR_i=[tau_i/m_i, 1];
denR_i=[tau_i, 1];

% funzione compensatrice complessiva

figure(1)
bode(conv(numR_a,numR_i),conv(denR_a,denR_i))
pause

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

numF=conv(numFint,numR_i);
denF=conv(denFint,denR_i);

w1=logspace(-1,3,1000);
[magFhat,phaseFhat]=bode(numFhat,denFhat,w1);
[magFint,phaseFint]=bode(numFint,denFint,w1);
[magF,phaseF]=bode(numF,denF,w1);
magFhatdB=20*log10(magFhat);
magFintdB=20*log10(magFint);
magFdB=20*log10(magF);
figure(1)
subplot(211)
semilogx(w1,magFhatdB,'y-',w1,magFintdB,'g-.',w1,magFdB,'c--'),grid
xlabel('Frequency (rad/sec)')
ylabel('Magnitude (dB)')
title('processo compensato (blu), intermedio (verde) e modificato (giallo)')
subplot(212)
semilogx(w1,phaseFhat,'y-',w1,phaseFint,'g-.',w1,phaseF,'c--',w1,-180*(1+0*w1),'w-'),grid
xlabel('Frequency (rad/sec)')
ylabel('Phase (deg)')
pause

% zoom diagrammi precedenti nella zona di interesse

w2=logspace(0,2,100);
[magFhat,phaseFhat]=bode(numFhat,denFhat,w2);
[magFint,phaseFint]=bode(numFint,denFint,w2);
[magF,phaseF]=bode(numF,denF,w2);
magFhatdB=20*log10(magFhat);
magFintdB=20*log10(magFint);
magFdB=20*log10(magF);
figure(1)
subplot(211)
semilogx(w2,magFhatdB,'y-',w2,magFintdB,'g-.',w2,magFdB,'c--'),grid
xlabel('Frequency (rad/sec)')
ylabel('Magnitude (dB)')
axis([1 100 -50 40])
title('processo compensato (blu), intermedio (verde) e modificato (giallo)')
subplot(212)
semilogx(w2,phaseFhat,'y-',w2,phaseFint,'g-.',w2,phaseF,'c--',w2,-180*(1+0*w2),'w-'),grid
xlabel('Frequency (rad/sec)')
ylabel('Phase (deg)')
axis([1 100 -270 -90])
pause

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

figure(2)
margin(numF,denF)
pause
close all

% se le specifiche su margine di fase e pulsazione di attraversamento
% provenivano da specifiche su modulo alla risonanza e banda passante
% ad anello chiuso si possono osservare i 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(1)
semilogx(w3,magWdB,'-',w3,-3+0*w3),grid
xlabel('Frequency (rad/sec)')
ylabel('Magnitude (dB)')
title('anello chiuso')
pause

% ...oppure si pu˜ usare la carta di Nichols 

figure(2)
deg=-360:0;
ngrid('new')
nichols(numF,denF,w1)
plot(deg,0*deg,'w')
hold off



