Mathematical Engineering - Fondamenti di Automatica

Lab 4 - Soluzioni 4

Laboratory

Matlab - Laboratorio 4 Esercizio 1: Diagrammi di Bode e teorema della risposta in frequenza a) Per tracciare il diagramma di Bode reale è sufficiente utilizzare la funzione Matlab bode (help bode); la chiamata base richiede soltanto la definizione del sistema dinamico di cui si vuole tracciare il diagramma. Per prima cosa definiamo il sistema % Definizione del sistema num = 10; den = poly([-1 -1 -1]); sistema = tf(num,den) Da notare come per definire il denominatore sia stata usata la funzione poly: questa funzione restituisce in un vettore i coefficienti ordinati per potenze decrescenti che definiscono il polinomio e le cui radici sono quelle fornite in ingresso alla funzione. In questo caso è immediato riconoscere come le radici del denominatore siano 3, tutte pari a -1. Un metodo alternativo prevede l’utilizzo della funzione conv che restituisce i coefficienti del polinomio prodotto di due polinomi in ingresso. Nel caso in esame vi sono 3 polinomi da moltiplicare, perciò sarà opportuno usare la funzione 2 volte come segue: den = conv(conv([1 1],[1 1]),[1 1]) A questo punto è possibile utilizzare la funzione bode per il tracciamento del relativo diagramma. % Visualizzazione diagrammi con la funzione bode figure; wmin=0.01; wmax=10; bode(sistema,{wmin,wmax}) grid on dove il sec ondo ingresso corrisponde alle pulsazioni in rad/s in cui si vuole visualizzare il grafico. Il risultato è riportato in figura. Notare il guadagno pari a 20db (10!) e le frequenze di visualizzazione del grafico, come richieste nella chiamata. Sovrapponiamo ora il diagramma approssimato con il comando bodeasin (NOTA: il comando non fa parte della comune libreria Control System Too lbox ma è stato sviluppato in modo autonomo; help bodeasin per il suo uso) -60-50-40-30-20-1001020Magnitude (dB) 10 -210 -110 010 1-270-180-900 Phase (deg) Bode Diagram Frequency (rad/sec) % Sovrapposizione diagrammi approssimati con la funzione bodeasin >> [omega,bms,omegasf,bfsf]=bodeasin(num,den,wmin,wmax); >> [ampiezza,fase,puls]=bode(sistema,{wmin,wmax}); >> figure; >> subplot(2,1,1); semilogx(omega,bms,'r'); >> subplot(2,1,2); semilogx(omegasf,bfsf,'r'); >> subplot(2,1,1); hold on; semilogx(puls,20*log10(ampiezza(1,:)')); grid on; >> subplot(2,1,2); hold on; semilogx(puls,fase(1,:)'); grid on; L’ errore massimo si trova in corrispondenza della pulsazione 1 (il sistema ha tre poli in 1); dato che in ogni singolarità l’errore è pari a circa 3 dB, l’errore massimo è pari a circa 3*3 dB = 9 dB = = 2.8184. Numericamente (in dB): >> [ampiezza1,fase1]=bode(sistema,1); >> erroremax=20-20*log10(ampiezza1) erroremax = 9.0309 b) Per il teorema della risposta in frequenza , l’uscita a transitorio esaurito del sistema è pari a: ( ) | ( )| ( ( ( )) )dove | ( )| è il modulo e ( ( )) è la fase della funzione di trasferimento, valutata in ( = pulsazione e U = ampiezza della sinusoide in ingresso) Pertanto per rispondere al quesito è necessario valutare la funzione di trasferimento in , e calcolarne modulo e fase. Per fare ciò viene in aiuto ancora il comando bode, che restituisce le grandezze desiderate. Oppure grazie alle funzioni abs e phase si possono calcolare modulo e fase di del numero complesso ricavato sostituendo . Da notare che essendo un sistema lineare, vale la sovrapposizione degli effetti e dunque si calcoleranno separatamente le risposte all’ingresso sinusoidale a frequenza 0.5 e a frequenza 0 (una costante è un coseno a frequenza 0!): l’uscita totale sarà la somma delle rispettive uscite. Il codice matlab è il seguente: % Calcolo della risposta in frequenza nel punto desiderato [G0modulo,G0fase] = bode(sistema,0)% Per l'ingresso costante [G05modulo,G05fase] = bode(sistema,0.5) % Per l'ingresso a 0.5 rad/s % Controllo con metodo alternativo per omega = 0.5 10 -110 010 1-60-40-20020magnitude (dB) 10 -110 010 1-300-200-1000 frequency (rad/sec) phase (deg) G05 = 10/(0.5j + 1)^3 % Valutazione di G(s) in s = j*omega G05ma = abs(G05)% Calcolo del modulo G05fa = phase(G05)*180/pi% Calcolo della fase Il risultato è il seguente: G05modulo = G05ma = 7.1554 7.1554 G05fase = G05fa = -79.6952 -79.6952 Ovviamente i risultati sono identici. Concludendo, l’uscita a transitorio esaurito sarà calcolabile grazie al teorema della risposta in frequenza: % Uscita secondo teorema di risposta in frequenza t = 0:0.001:100; % Definizione del vettore tempo y_rf = 15*G0modulo - 3*G05modulo*sin(0.5*t + G05fase*pi/180); % Uscita % grafici u = 15 - 3*sin(0.5*t);% Definizione ingresso y_sim = lsim(sistema,u,t); % Simulazione uscita da condizione iniziale nulla figure% Grafico dei risultati plot(t,y_sim,t,y_rf,'g','LineWidth',2) grid on legend('Simulata','Risposta in frequenza') xlabel('Tempo [s]') ylabel('Uscita') Il risultato è la seguente figura: Come si può notare le due uscite sono equivalenti A REGIME, ovvero quando il transitorio iniziale si è esaurito (tempo di risposta pari a 5). 0102030405060708090100020406080100120140160180 Tempo [s]Uscita [EU] Simulata Risposta in frequenza Esercizio 2 : Diagrammi di Bode – analisi sistema del secondo ordine La funzione di trasferimento del sistema meccanico è data da ( ) I poli sono in √ e sono complessi coniugati se . Tale condizione è verificata per i valori proposti nel testo dell’esercizio, così che la funzione di trasferimento è del tipo ( ) con √ , √ e . Sostituendo i valori per m e k si ha , , e ( ) Tracciamo ora i diagrammi di Bode per h [ ] che corrispondono a valori di smorzamento pari a [ ] ; % Definizione di sistemi num = 1; den1 = [2 7 8]; G1 = tf(num,den1); den2 = [2 4 8]; G2 = tf(num,den2); den3 = [2 0.4 8]; G3 = tf(num,den3); den4 = [2 0.04 8]; G4 = tf(num,den4); den5 = [2 0.01 8]; G5 = tf(num,den5); % Diagrammi di Bode bode(G1,G2,G3,G4,G5) legend('7' ,'4','0.4','0.04','0.01') -100-50050Magnitude (dB) 10 -210 -110 010 110 2-180-135-90-450 Phase (deg) Bode Diagram Frequency (rad/sec)7 4 0.4 0.04 0.01 Si può notare come , al diminuire dello smorzamento, cresca l’ampiezza del picco di risonanza e cresca la derivata della fase in corrispondenza della pulsazione di risonanza che è pari a 2 rad/s. Per simulare la risposta a ingresso sinusoidale di ampiezza unitaria e pulsazione 2 utilizziamo i seguenti comandi % risposte a ingresso sinusoidale [G2modulo,G2fase] = bode(G2,2); [G3modulo,G3fase] = bode(G3,2); [G4modulo,G4fase] = bode(G4,2); % Uscita secondo teorema di risposta in frequenza t = 0:0.001:100;% Definizione del vettore tempo y_rf_2 = G2modulo*sin(2*t + G2fase*pi/180); y_rf_3 = G3modulo*sin(2*t + G3fase*pi/180); y_rf_4 = G4modulo*sin(2*t + G4fase*pi/180); figure; plot(t,y_rf_2,t,y_rf_3,t,y_rf_4); legend('4','0.4','0.04') xlabel('tempo'); ylabel('risposta'); Al diminuire dello smorzamento si nota come le oscillazioni aumentino in ampiezza. Tale comportamento si spiega con la presenza del picco di risonanza nella pulsazione 2. Esercizio 3 : Sistemi di controllo – Analisi di stabilità a) La funzione d’anello L= ( ) ( ) rispetta i criteri per l’applicazione del criterio di Bode per l’analisi di stabilità del sistema di controllo dato (L non ha poli a parte reale positiva e il diagramma di Bode del modulo di L(s) interseca una e una sola volta l’asse a 0 dB). % Definizione L num = 10*[10 1]; den = poly([- 1 -1 -1 -1]); 0102030405060708090100-15-10-5051015 temporisposta 4 0.4 0.04 L = tf(num,den); bode(L); grid; Grazie al comando margin è possibile calcolare margine di fase e margine di guadagno. % Calcolo del margine di fase e di guadagno [GainM,PhaseM,Wcg,Wcp] = margin(L) I valori restituiti da Matlab sono i seguenti GainM = PhaseM = 0.1796 -41.1034 Da cui è possibile notare che, nonostante il guadagno sia positivo, essendo margine di fase (PhaseM) negativo il sistema in anello chiuso è instabile. Utilizzando il semplice comando margin senza richiedere alcuna uscita, viene visualizzato anche il diagramma di Bode \ di L: figure margin(L) -40-2002040Magnitude (dB) 10 -310 -210 -110 010 1-270-180-90090 Phase (deg) Bode Diagram Frequency (rad/sec) La risposta a gradino conferma l’instabilità del sistema: % Simulazione al gradino del sistema sotto controllo step(feedback(L,1),10) -40-30-20-10010203040Magnitude (dB) 10 -310 -210 -110 010 1-270-180-90090 Phase (deg) Bode Diagram Gm = -14.9 dB (at 2.34 rad/sec) , Pm = -41.1 deg (at 4.49 rad/sec) Frequency (rad/sec) 012345678910-1-0.500.511.52 x 10 4Step Response Time (sec)Amplitude b) Modificare il guadagno di L(s) =( ) ( ) significa traslare il modulo della L(s) valutata al punto precedente verso l’alto ( >10) o verso il basso ( 