Mathematical Engineering - Fondamenti di Automatica

Lab 2 - Soluzione 2

Laboratory

Matlab - Laboratorio 2 Esercizio 1 (file esercizio1.mdl) Per studiare la stabilità asintotica del sistema utilizzando Simulink occorre innanzitutto realizzare il sistema in modo che sia simulabile. Dopodiché, per ingresso nullo, se ne osserverà il movimento libero a partire da differenti stati iniziali. Dopo avere digitato simulink al prompt di Matlab, aprire un nuovo file su cui verranno “trascinati” gli elementi necessari per la definizione e simulazione del sistema. Avendo a disposizione le matrici è possibile utlizzare il blocco State-Space di Simulink (libreria Continuous): Per poter simulare il sistema, occorre definire come ingresso del sistema una costante. Si utilizzi il blocco Constant, all’interno della libreria Sources: P er rappresentare l’uscita in un grafico, utilizziamo il blocco X-Y graph, all’interno della libreria Sinks; questo blocco necessita di due ingressi: il primo è il vettore che verrà rappresentato lungo l’asse x, il secondo quello che verrà rappresentato sull’asse y. La variabile indipendente, in questo caso è il tempo. Prendiamo dunque il blocco clock (libreria Sources) per generare il tempo e colleghiamolo alla porta “x” del blocco X-Y graph. Connettendo i blocchi come descritto si trova il seguente schema: Si noti che il valore dell’ingresso è stato già messo a 0. Ora occorre definire le matrici del blocco State-Space. Cliccando due volte sul blocco ne appaiono le proprietà; qui si possono definire le matrici come indicate nel testo dell’esercizio; dunque, per p = 2, si avrà: Si noti che è stato impostato uno stato iniziale non nullo (in (1,1)) per poter osservare il moto libero del sistema. Ora è possibile simulare il comportamento; infine, una volta che si ha un riscontro grafico è possibile cambiare la scala del blocco X-Y graph in modo da ottenere una visualizzazione più chiara, agendo sulle proprietà del blocco: Simulando ora per 10 secondi (Simulation, Run) si ottiene il seguente risultato: L ’uscita del sistema tende a zero. Dato che tale comportamento sussiste variando a piacere la condizione iniziale del sistema (verificare!!), il sistema è asintoticamente stabile (NOTA: l’asintotica stabilità viene dedotta via simulazione, in modo ovviamente né esaustivo né rigoroso). Ciò poteva essere osservato calcolando gli autovalori della matrice A: >> A=[0 -1;-2 -2]; >> e ig(A) ans = -1.0000 +1.0000i -1.0000 -1.0000i Pertanto, avendo il sistema due autovalori complessi coniugati a parte reale negativa, si spiega il comportamento oscillatorio smorzato del sistema (fuoco stabile). Se vogliamo rappresentare x 1 occorre prendere come uscita l’equazione: y(t) = x 1( t) da cui c = [1 0]. Analogamente, per rappresentare x 2 deve essere: y(t) = x 2( t) da cui c = [0 1]. Quindi, basta copiare il sistema precedentemente costruito cambiando i parametri del vettotre c, ottenendo lo schema: Per rappresentare le traiettorie nel piano ( x 1, x 2) occorre aggiungere un altro X-Y graph alimentato da x1 e x 2:Simulando il sistema si ottengono i seguenti andamenti ( y è già stato riportato): Il piano delle fasi risulta essere: Infine, per simulare il sistema sottoposto a un ingresso forzato u = 1, è sufficiente cambiare il valore della costante in ingresso, modificando il parametro del blocco Constant: Le uscite di simulazione e il quadro delle traiettorie saranno analoghi a quanto trovato in precedenza a parte il valore di equilibrio raggiunto dal sistema. Per studiare la stabilità del sistema per gli altri valori di p proposti occorre modificare di volta in volta tutti i blocchi State-Space (matrice A) e ripetere quanto visto in precedenza per p = 2. I risultati sono i seguenti: p = 0.5 autovalori = 0.2929,  1.7071 nodo stabile, asintotica stabilità p =  autovalori = 0, 2 semplice stabilità (il movimento libero è limitato e gli stati con x 2 nullo sono tutti stati di equilibrio semplicemente stabili) p =  autovalori = 2.4142, 0.4142 sella, instabilità (le variabili divergono, nello spazio di stato si osserva l’effetto sella in prossimità dell’equilibrio nullo corrispondente a ingresso nullo) Esercizio 2 (file esercizio2.mdl) In questo caso viene data la funzione di trasferimento del sistema (che è a tempo discreto). Pertanto, in Simulink, si può utilizzare direttamente il blocco transfer fcn, che nel caso discreto si chiama Discrete Transfer Fcn (libreria Discrete): Come ingresso diamo una costante u = 1 e come uscita rappresentiamo il tutto, come fatto nell’esercizio precedente, in un grafico X-Y; lo schema è dunque: Per inserire i parametri nel blocco della funzione di trasferimento basta aprire le proprietà e scrivere numeratore e denominatore. Ciò va fatto sotto forma di vettori riga: Una volta inseriti i parametri, si può simulare il sistema e successivamente cambiare, se necessario, la scala del grafico e il tempo di simulazione. In questo caso si è selezionato un tempo di simulazione di 15 istanti di tempo; i limiti dati al grafico sono i seguenti: Il grafico che ne risulta è il seguente: Si noti che il comportamento lascia presupporre che il sistema sia stabile asintoticamente; calcolando i poli della funzione di trasferimento si ottiene: >> poli=roots([1 0.7 0.12 -0.09]) poli = -0.4755 + 0.3641i -0.4755 - 0.3641i 0.2509 >> abs(poli) ans = 0.5989 0.5989 0.2509 Il che conferma la stabilità asintotica del sistema (tutti gli auto valori hanno modulo inferiore a 1) e la presenza di oscillazioni (vi sono auto valori complessi). Per rappresentare le variabili di stato è possibile riscrivere il sistema in forma di stato e utilizzare un metodo analogo a quanto già visto nell’esercizio 1; per trovare il sistema in forma di stato basta utilizzare il comando di matlab tf2ss; successivamente si utilizzi il blocco state-space discreto inserendo le matrici (A, B,C,D) trovate e proseguire come per l’esercizio 1 (scegliendo ora Discrete State-Space nella libreria Discrete), duplicando il sistema per selezionare, una per volta, le variabili di stato. >> [A,B,C,D]=tf2ss([0 3 -0.1 -0.43],[1 0.7 0.12 -0.09]) A = -0.7000 -0.1200 0.0900 1.0000 0 0 0 1.0000 0 B = 1 0 0 C = 3.0000 -0.1000 -0.4300 D = 0 Esercizio 3 (file esercizio3.mdl) Per risolvere questo esercizio occorre tradurre lo schema a blocchi dato nel testo in uno schema Simulink: Supponendo che i blocchi 1, 2 e 3 abbiano funzione di trasferimento per rappresentare lo schema precedente in Simulink si può utilizzare il blocco Transfer Fcn per tre volte in serie: Per sommare o sottrarre due segnali esiste il nodo sommatore (denominato appunto Sum, libreria Commonly Used Blocks); per cambiare i segni del nodo sommatore occorre entrare nelle proprietà dello stesso e specificarli nel campo “list of signs”. Infine, per moltiplicare un segnale (guadagno “k”) si utilizza un blocco denominato gain (libreria Commonly Used Blocks). L’ingresso, come richiesto dal testo, sarà una costante con valore 1. Si può pertanto costruire il seguente schema: + \ + d  1 1   s sG Si noti che è stato inserito il solito blocco Clock per generare il tempo da inserire nel blocco X-Y Graph. La richiesta dell’esercizio è di confermare i calcoli svolti a lezione fossero corretti, ovvero se il limite di stabilità del sistema fosse raggiunto per K = 8. Per fare ciò svolgiamo 2 simulazioni con 2 diversi valori di K (minore e maggiore di 8); annullando il disturbo (la stabilità interna di un sistema non dipende dagli ingressi) otteniamo: - K = 4 : Come si nota il sistema perde la stabilità per valori di K elevati. Per diagrammare l’andamento dell’uscita all’equilibrio al variare di K (nel caso di asintotica stabilità), fissiamo il parametro di ingresso pari a 5 e il parametro di disturbo pari a 1; ripetiamo diverse simulazioni per differenti valori di K riportando il valore dell’uscita corrispondente in una tabella. Il risultato ottenuto è il seguente: K 0 0.3 0.6 1 2 3 4 5 6 7 uscita 6 5.769 5.625 5.5 5.33 5.25 5.2 5.166 5.143 5.125 NOTA : per visualizzare l’uscita nel workspace di Matlab, collegare il blocco To Workspace (libreria Sinks) all’uscita y selezionando array alla voce Save format. I valori ottenuti (asterischi blu) ben approssimano il risultato ottenuto in teoria (linea rossa) - K = 10: Pertanto, per ridurre l’effetto del disturbo sull’uscita occorre portare K a valori elevati, avvicinandosi però così al limite di stabilità per il sistema. K d uy   1 01234567855.25.45.65.866.2 Kuscita Esercizio 4 (file rosmca.mdl) Per la definizione del sistema, si può utilizzare il blocco Simulink User Defined Function; infatti il sistema si può vedere come: Pertanto basta definire due funzioni, che abbiano ognuna i due ingressi x e y e siano fatte come f e g. Lo schema si semplifica molto, ma ora occorrono dei blocchi multiplexer (mux in Commonly Used Blocks) per generare gli ingressi delle funzioni; infatti il blocco User Defined Function accetta in ingresso un vettore colonna del tipo: Nel nostro caso u deve contenere x e y. Quindi lo schema, completo anche dei grafici per osservare i risultati della simulazione, diventa:     yxgyyxfx ,,                     tu tu tu n 1 e = 0.5 m = 1 Tenuto conto che i parametri dati nel testo del problema sono r = 1 K = 10 a = 6 b = 2 si deve porre Il piano delle fasi risulta essere: Si noti come nasca un comportamento ciclico. Per valutare come cambia il comportamento al cambiare del parametro K, simuliamo dapprima il comportamento del sistema per i seguenti valori di K: f(u)=1*u(1)*(1-u(1)/10)-6*u(1)*u(2)/(2+u(1)) g(u)=0.5*6*u(1)*u(2)/(2+u(1))-1*u(2) Le condizioni iniziali u(1) = 4 e u(2) = 2 vanno fissate nel blocco integratore 1/ s.Simulando si ottiene: K = 5 K = 3 Come si può notare, al diminuire di K, il comportamento passa da coesistenza ciclica a coesistenza stazionaria di prede e predatori. Inoltre, il ciclo diminuisce in ampiezza. Pertanto il sistema subisce una biforcazione chiamata biforcazione di Hopf. Diminuendo ulteriormente K si ottiene, fissando ad esempio K al valore 0.5, si ottiene Pertanto, al diminuire di K, il comportamento passa da coesistenza stazionaria di prede e predatori a presenza di sole prede alla capacità portante K. Il sistema subisce una biforcazione chiamata biforcazione transcritica o scambio di stabilità (l’equilibrio (K,0), al diminuire di K, da instabile diventa stabile). Esercizio 5 (file rosmca3.mdl) Per risolvere questo esercizio bisogna schematizzare in Simulink il sistema dinamico non lineare espandendo e modificando il modello dell’esercizio precedente ottenendo il seguente schema Simulink: Fissando i parametri ai valori forniti nel testo e simulando il sistema si ottengono i seguenti grafici: Il comportamento è pertanto quello tipico di un attrattore caotico visualizzabile digitando al prompt di Matlab i seguenti comandi >> plot3(x1,x2,x3) >> grid >> axis([0 1.1 0 1.5 4 14]) >> xlabel('x'); ylabel('y'); zlabel('z'); Esercizio 6 (file competizione.mdl) Lo schema Simulink che caratterizza il modello di competizione è dato da: Le funzioni F(u) e G(u) sono definite come specificato nel testo: F(u) = 1*u(1)*(1-u (1)/3)-1*u(1)*u (2) G(u) = 4*u(2)*(1-u (2)/2)-2*u (1)*u (2) Per verificare che per K 1 = 3 e K 2 = 2 il sistema ammette due equilibri stabili alternativi ( K 1, 0) e (0, K 2), occorre simulare il sistema per differenti condizioni iniziali e osservare che vi sono condizioni iniziali che fanno tendere il sistema verso l’equilibrio di dominanza della specie x 1 ed altre condizioni iniziali che fanno tendere il sistema verso l’equilibrio di dominanza della specie x 2. A tale scopo fissiamo in (2, 1) le condizioni iniziali del sistema (scegliendo 2 (1) come Initial condition del blocco integratore che restituisce la variabile x 1 ( x 2) ). Simulando si osserva che il00.20.40.60.81 00.511.5468101214 xyz sistema tende verso l’equilibrio ( K 1, 0). La condizioni iniziale (2, 1.4) porta invece il sistema verso l’equilibrio alternativo (0, K 2). Provare con le condizioni iniziali (1, 0.2) e (1, 0.1 ). Il quadro delle traiettorie risultante dalle 4 simulazioni è dato da che mostra quanto si voleva verificare. Verifichiamo ora che al diminuire del parametro K 1 il sistema tende sempre verso la dominanza della specie x 2. Fissiamo, ad esempio, K 1 = 1.5 andando a modificare l’espressione di F(u ): F(u) = 1*u(1)*(1-u (1)/1.5)-1*u(1)*u (2) e simuliamo il sistema. Come si può notare, anche partendo da condizioni iniziali prossime a (K 1, 0), il sistema tende sempre verso l’equilibrio (0, K 2) (si provi a simulare con condizioni iniziali (1.5, 0.1), (1, 0.1) e (2, 0.1)). Riportiamo ora il parametro K 1 al suo valore iniziale F(u) = 1*u(1)*(1-u (1)/3)-1*u(1)*u (2) Verifichiamo ora che al diminuire del parametro K 2 il sistema tende sempre verso la dominanza della specie x 1. Fissiamo, ad esempio, K 2 = 0.5 andando a modificare l’espressione di G(u ): G(u) = 4*u(2)*(1-u (2)/0.5)-2*u(1)*u (2) 00.511.522.533.500.511.522.5 x 1x2 K1 = 3 K 2 = 2 00.511.522.533.500.511.522.5 x 1x2 K1 = 1.5 K 2 = 2 e simuliamo il sistema. Come si può notare, anche partendo da condizioni iniziali prossime a (0, K 2), il sistema tende sempre verso l’equilibrio (K 1, 0) (si provi a simulare con condizioni iniziali (0.1, 0.5), (0.1, 0.1) e (0.1, 1)). Verifichiamo ora che al diminuire dei parametri K 1 e K 2 il sistema tende verso la coesistenza stazionaria della due specie. Fissiamo, ad esempio, K 1 = 1.5 e K 2 = 0.5 andando a modificare l’espressione di F(u) (G(u) presenta già il valore desiderato per K 2): F(u) = 1*u(1)*(1-u (1)/1.5)-1*u(1)*u (2) e simuliamo il sistema. Come si può notare, partendo da condizioni iniziali qualsiasi (anche prossime a (K 1, 0) e (0, K 2), il sistema tende sempre verso la coesistenza stazionaria di x 1 e x 2. Per giustificare i cambiamenti di comportamento al variare di K 1 e K 2 dal punto di vista grafico, si faccia riferimento al file addendum_ese6.pdf. 00.511.522.533.500.511.5 x 1x2 K1 = 3 K 2 = 0.5 00.20.40.60.811.21.41.61.8200.10.20.30.40.50.60.70.80.91 x 1x2 K1 = 1.5 K 2 = 0.5