Mathematical Engineering - Fondamenti di Automatica

Lab 1 - Soluzione Exercizi

Laboratory

Matlab - Laboratorio 1Tracce di possibili soluzioni: Esercizio 1: La stabilità del sistema dipende dagli autovalori della matrice dinamica A: in particolare il sistema è asintoticamente stabile se e solo se tutti gli autovalori di A hanno modulo inferiore a 1. Quindi, occorre innanzitutto definire la matrice A come definita dal testo dell’esercizio, per poi calcolarne gli autovalori tramite la funzione eig(). Successivamente si può calcolarne il modulo tramite la funzione abs(). Pertanto, si può scrivere: A=[-0.8 1 -1; -0.5 0.6 -1; -1 1 -0.5]; % matrice A autoval=eig(A) % autovalori di A autoval = -0.4755 + 0.3641i -0.4755 - 0.3641i 0.2509 modautoval=abs(autoval) % modulo degli autovalori di A modautoval = 0.5989 0.5989 0.2509 Come si vede, il modulo dei tre autovalori è strettamente minore di 1 quindi il sistema è asintoticamente stabile. Inoltre, dato che il sistema ha due autovalori complessi coniugati, ci si può aspettare una risposta allo scalino che presenta oscillazioni smorzate. Il tempo di risposta del sistema sarà pari a () ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛− = D R T λ ln 1 5 dove Dλ è l’autovalore dominate del sistema cioè quello/i di modulo massimo. Pertanto, () = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛− = 5989.0ln 1 5 RT 9.753 circa. Per calcolare l’equilibrio del sistema sottoposto a ingresso costante 2 = u occorre porre: ubxAx += da cui si ottiene ( ) ubAIx 1 − −= e udxcy +=. Per invertire la matrice (con il comando inv() ) occorre innanzitutto costruire una matrice identità delle stesse dimensioni di A; per questo sono utili sia il commando eye() che il commando size() : ( AI −) b=[1; 1; 1]; % Vettore b c=[1 1 1]; % Vettore c d=0; % Vettore d u=2; % Vettore degli ingressi % Calcolo dell’equilibrio: mi occorre una ma trice identità –eye()- delle di mensioni di A –s ize()-; pertanto mi occo rre eye( size(A)); quindi l’equilibrio si può calcolare come: xeq=inv(eye(size(A))-A)*b*u xeq = 0.8092 0.7514 1.2948 yeq=c*xeq+d*u yeq = 2.8555 Per sim ulare il com portam ento del sistem a occorre innanzitutto defi nire il m odello dinam ico a tem po discreto in matlab e poi utilizzare il comando lsi m() . Attenzione, essendo un sistem a a tem po discreto, le traiettorie saranno sequenze di punti, n on linee continue. I com andi ch e occorre conoscere sono ss() per l a definizione del sistem a dinam ico e, appunto , lsim() ; in quest’ult im o co mando occorre specificare, o ltre al sistem a da simu lare, l’ingre sso al siste ma (nel nostro caso 2=u ) unitam ente al vettore dei tem pi che lo ha generato. Nel com ando ss() , invece, bisogna specificar e le matrici del si stem a e, nel caso di sistema a tempo discreto, il tem po di c ampionam ento, in quest o esercizio unitario (il valore -1 i ndica un te mpo d i cam pionam ento generico). Dunque si p uò scrivere: sistema=ss(A,b,c,d,1) %sistema a tempo discreto – tsampling = 1 a = x1 x2 x3 x1 -0.8 1 -1 x2 -0.5 0.6 -1 x3 -1 1 -0.5 b = u1 x1 1 x2 1 x3 1 c = x1 x2 x3 y1 1 1 1 d = u1 y1 0 Sampling time: 1 Discrete-time model. %creo un ingresso costante che duri un certo numero di campioni t=0:20; %vettore di tempi a passo unitario u=2*ones(size(t)); %ingresso: vale sempre due %simulo il sistema, salvando uscita, tempo e traiettorie degli stati: [out,time,state]=lsim(sistema,u,t); %NOTA: la condizione iniziale sullo stato è nulla; in alternativa va specificata come X0 in lsim(sistema,u,t,X0) %rappresento i risultati: figure %apro una nuova figura plot(time,out, '-*' ) %uscita hold all plot(time,u) %ingresso xlim([-1,20]) %visualizzo l’asse x da -1 a 20 ylim([0,7]) %visualizzo l’asse y da 0 a 7 grid on title( 'Uscita del sistema a un ingresso costante u=2' ) xlabel( 'tempo [campioni]' ) legend( 'output' ,'input' ) %aggiungo la legenda 0 2 4 6 8 10 12 14 16 18 20 0 1 2 3 4 5 6 7 Uscita del sistema a un ingresso costante u=2 tempo [ campioni] output input %rappresento l’andamento dei tre stati: figure plot(time,state(:,1), '-*' ) %prendo la prima colonna del ve ttor e degli stati hold all plot(time,state(:,2), '-*' ) plot(time,state(:,3), '-*' ) xlim([-1,20]) ylim([0,3]) grid on title( 'Movimento per ingresso costante u=2' ) xlabel( 'tempo [campioni]' ) legend( 'x_1' ,'x_2' ,'x_3' ) 0 2 4 6 8 10 12 14 16 18 20 0 0.5 1 1.5 2 2.5 3 Movimento pe r ingresso costante u=2 tempo [ cam pioni] x1 x2 x3 Infine, volendo tracciare l e traiettorie nello spazio di stato, si utilizza il co mand o plot3() co me segue plot3(state(:,1),state(:,2),state(:,3),'*') Eser cizio 2 Analogam ente a quanto vi sto nell’esercizio 1, per st udiare la stabilità del sistema e per val utarne il co mportam ento qualitativo, occorre conoscere gli autovalori della matrice A. Differentemente dal caso precedente, il sistema è a tem po continuo; la st abilità asintotica è quindi garantita se e solo se tutti gli aut ovalori di A han no parte reale negativa. Dunque, analogamente a quanto visto nell’ esercizio precedente si può scrivere: A=[-1 -1 2; 0 -2 10; 0 -10 -2]; % matrice A autoval=eig(A) % autovalori di A autoval = -1.0000 -2.0000 +10.0000i -2.0000 -10.0000i REval= real(autoval) REval = -1 -2 -2 Pertanto si conclude che il sistema è asintoticam ente st abile. Inoltre la m atrice A ha due autovalori c omplessi coniugati; pert anto, ci si può aspettare una risposta allo scalino che presenta oscillazioni s morzate. La costante d i te mpo dom inante del sist ema è pari a () ( )max 1 λ ℜ −= DT m entre il tempo di risposta è pari a . Pertanto e D R T T 5= 1= DT 5=RT . Per calcolare l’equilibri o de l sistem a, occorre i mporre ub xA + =0 da cui si ottiene ub A x 1− −= e ud xc y + = . Co me nell’ esercizio 1, occorre prim a definire le matrici del sistema b=[1; 1; 1]; % Vettore b c=[1 1 1]; % Vettore c d=0; % Vettore d u=1; % Vettore degli ingressi % Calcolo dell’equilibrio: xeq=-inv(A)*b*u xeq = 0.7308 0.1154 -0.0769 yeq=c*xeq+d*u yeq = 0.7692 Per simulare il comportamento occorre definire (con il comando ss() ) il modello dinamico. Successivamente, grazie al comando lsim() , si potrà osservare l’andamento delle variabili di interesse. Attenzione, considerando ora un sistema a tempo continuo , le traiettorie saranno linee continue e non sequenze di punti. sistema=ss(A,b,c,d) % sistema a tempo continuo a = x1 x2 x3 x1 -1 -1 2 x2 0 -2 10 x3 0 -10 -2 b = u1 x1 1 x2 1 x3 1 c = x1 x2 x3 y1 1 1 1 d = u1 y1 0 Continuous-time model. % Definisco il tempo di durata della simulazione (NOTA: è maggiore del tempo di risposta TR del sistema) t=0:0.001:7; % Definisco l’ingresso u, costante e pari a 1 u=ones(size(t)); % Simulo il sistema: [uscita,tempo,stato]=lsim(sistema,u,t); %NOTA: la condizione iniziale sullo stato è nulla; in alternativa va specificata come X0 in lsim(sistema,u,t,X0) %rappresento i risultati: figure %apro una nuova figura plot(tempo,uscita) %uscita hold all plot(tempo,u) %ingresso ylim([0,1.2]) %visualizzo l’asse y da 0 a 1.2 grid on title( 'Uscita del sistema a ingresso unitario’ ) xlabel( 'tempo [s]' ) % rappresento anche il guadagno statico del sistema: guadagno=dcgain(sistema)*ones(size(t)); % è circa pari a 0.7692 plot(tempo,guadagno) legend( 'uscita' ,'ingresso' ,'guadagno statico del sistema' ) 0 1 2 3 4 5 6 7 0 0.2 0.4 0.6 0.8 1 Uscita del sistema a un ingresso unitario tempo [s] uscita ingresso guadagno statico del sistema % rappresento gli stati: figure plot(tempo,stato(:,1)) %prendo la prima colonna del ve ttor e degli stati hold all plot(tempo,stato(:,2)) plot(tempo,stato(:,3)) grid on title('Movimento per ingresso costante u=1') xlabel('tempo [s]') legend('x_1','x_2','x_3') 0 1 2 3 4 5 6 7 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 M ovim ento per ingres s o c os tante u= 1 tempo [s] x1x2x3 Infine, volendo tracciare le traiettorie nello spazio di stato, si utilizza il comando plot3() come segue plot3(stato(:,1),stato(:,2),stato(:,3)). I comandi hold on e plot3(xeq(1),xeq(2),xeq(3),'*') permettono di sovrapporre lo stato di equilibrio che viene asintoticamente raggiunto. Si definisca semplicemente ora un nuovo ingresso come 5 volte l’ingresso precedente; dunque: u2=5*u; [uscita2,tempo2,stato2]=lsim(sistema,u2,t); %rappresento i risultati: figure %apro una nuova figura plot(tempo2,uscita2) %uscita hold all plot(tempo2,u2) %ingresso ylim([0,5.2]) %visualizzo l’asse y da 0 a 5.2 grid on title('Uscita del sistema a un ingresso u=5’) xlabel('tempo [s]') % rappresento anche l’uscita di regime: si noti che l’uscita di regime, essendo l’ingresso costante, non è altro che l’ingresso\ per il guadagno della funzione di trasferimento del sistema regime=5*guadagno; plot(tempo2,regime) legend('uscita','ingresso','regime del sistema per u=5') 0 1 2 3 4 5 6 7 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Uscita del sistema a un ingr esso u=5 tempo [s] uscita ingresso regime del sistema per u=5 % rappresento gli stati: figure plot(tempo2,stato2(:,1)) %prendo la prima colonna del ve ttor e degli stati hold all plot(tempo2,stato2(:,2)) plot(tempo2,stato2(:,3)) grid on title( 'Movimento per ingresso costante u=5' ) xlabel( 'tempo [s]' ) legend( 'x_1' ,'x_2' ,'x_3' ) 0 1 2 3 4 5 6 7 -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 Movimento per ingresso costante u=5 tempo [s] x1x2x3 Co me si può notare la risposta è del tutto analoga a quella trovata al punto 3, fatto salvo per un fattore moltiplicativo 5; ciò era prevedibile, vi ste le proprietà di line arità del sistema. Inoltre l’uscit a di regim e è pari al nuovo ingre sso 5 m oltiplicato per il guadagno statico del sistema (circa pari a 0 .769 2). Eserciz io 3 Innanzitutto occorre scrivere il sistem a dinam ico nella form a: () () () () () ()t du t cx ty t bu t Ax tx + = + = & Per scrivere il sistem a dinam ico in form a matriciale si può porre () () () ()tx t x tx t x &= = 2 1 A questo pun to si può scrivere: () () () () () () () ()t x ty t x k t x a M t x t x t x s s s 1 1 2 2 2 1 1 = − − = = & & da cui ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − − = s s s s M a M k A 1 0 ⎥⎦ ⎤ ⎢⎣ ⎡= 0 0 b [ ]0 1=c d = 0 Per studiare la stabilità del sistema, occo rre conoscere gli autovalori della matrice A. Il polinom io caratteri stico del sistema è: () s s s s s s s s M k M a M a M k A I + + =⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ + − = − λ λ λ λ λ 2 1 det det da cui si possono calcolare gli auto valori: s s s s M k a a 2 2 2 2,1 − ± −= λ . Essendo sia as che ks positivi, gli autova lori avranno sem pre parte reale negativa. Pertanto il si stema è sem pre asin toticam ente st abile (co me d oveva essere rappresentan do u n sistem a dissipativo) Controlliam o num ericamente che la soluzione analitica sia corretta. Assegniam o i valori dati all e variabili del problem a: M_s = 120; % Kg k_s = 160; % N/m a_s = 80; % N*s/m Ora scriviamo la matrice A e valutiamone gli autovalori: A=[0 1;   -k_s/M_s -a_s/M_s]; autoval=eig(A) -0.3333 + 1.1055i -0.3333 - 1.1055i Quindi, essendo entrambi gli autovalori a parte reale negativa, il sistema risulta asintoticamente stabile. Per risolvere il secondo punto occorre riscrivere l'equazione del sistema. Nell’equazione precedente, infatti, la variabile x descriveva la dinamica della massa sedile+passeggero quando l’auto era ferma. Per studiare i casi in cui l’auto si muove verticalmente a causa del manto stradale irregolare occorre riscrivere l'equazione tenendo conto che il moto del sedile sia relativo al moto del manto stradale, ovvero scrivendo: relsrels sedilesxaxkxM& && −−= dove: strada sedile relxxx −= Quindi, l'equazione dinamica diventa: stradasstradassedilessedilessedilesxaxkxkxaxM& &&& + = + + Considerando come ingresso la forma della strada si ottiene un sistema dinamico della forma: () () () () () () () () txtyM xkxa txktxa M txtxtx s stradasstradas s s s 11 2 221 1 =+ +−−= = & & & dunque si scriverà l'ingresso come: ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = strada stradax x u& B sarà dunque una matrice della forma: ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ =s s s s M a Mk B00 Si può ora scrivere il nuovo sistema, aggi ornando i valori delle matrici B e d: B=[0 0; k_s/M_s a_s/M_s]; c=[1 0]; d=[0 0]; sys_sedile_2=ss(A,B,c,d) a = x1 x2 x1 0 1 x2 -1.333 -0.6667 b = u1 u2 x1 0 0 x2 1.333 0.6667 c = x1 x2 y1 1 0 d = u1 u2 y1 0 0 Continuous-time model. A questo punto, p er simulare cor rettamente il sistema, o ccor re costruire l'ingresso e la sua derivata. La funzione genera_dosso(vel) (digita help genera_dosso) modellizza il do sso come u na rampa a pendenza costante e di lunghezza pari a 30 cm e altezza 20 cm (hd osso), seguita da u n piano di 30 cm (ld osso) e da una rampa simmetr ica alla precedente. Chiamand o S la p endenza della rampa si può modellizzare il dosso (cioè l’altezza del manto stradale) in funzione d ella velocità di percorrenza del veicolo ( vca r): () () () ⎪ ⎩ ⎪ ⎨ ⎧