% Sintesi per tentativi
% Esempio di uso della funzione anticipatrice in congiunzione
% con un guadagno
% cfr Lanari,Oriolo "Controlli Automatici: Esercizi di Sintesi"
% es. 2.2, p. 71
% 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 e' w*_t=12
% il margine di fase desiderato e' m*_fi=40

numFhat=[10];
denFhat=conv([1,0],conv([0.2,1],[0.05,1]));

% 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

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

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

% per ottenere un anticipo di fase di circa 50 gradi si sceglie m_a=12
% e w tau_a =1.5 (con w=w*_t=12) dai diagrammi universali
% la corrispondente amplificazione e' benefica perche' si deve
% aumentare comunque la pulsazione di attraversamento 

tau_a=1.5/12;
m_a=12;
numR=[tau_a, 1];
denR=[tau_a/m_a, 1];

figure(2)
bode(numR,denR)
pause

close all

% calcolo valore del modulo del processo intermedio Fint(s)=Fhat(s)*R(s) 
% alla pulsazione di attraversamento desiderata e determinazione del 
% guadagno necessario per conseguire quest'ultima

numFint=conv(numFhat,numR);
denFint=conv(denFhat,denR);

w=12;
[magw,phasew]=bode(numFint,denFint,w);

K=1/magw;

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

numF=K*numFint;
denF=denFint;

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,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)')
title('processo compensato (blu), intermedio (verde) e modificato (giallo)')
subplot(212)
semilogx(w2,phaseFhat,'y-',w2,phaseF,'c--',w2,-180*(1+0*w2),'w-'),grid
xlabel('Frequency (rad/sec)')
ylabel('Phase (deg)')
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



