Università degli Studi di Roma Tor Vergata FACOLTÀ DI INGNEGNERIA Corso di Laurea Magistrale in Ingegneria dell’automazione Progetto per il corso di controllo dei processi Analisi e controllo di uno scambiatore di calore Candidato: Relatore: Francesco Ferrante Prof Riccardo Marino Matricola 0157498 Anno Accademico 2011–2012 francesco ferrante A N A L I S I E C O N T R O L LO D I U N O S C A M B I ATO R E D I C A LO R E Il matematico gioca un gioco in cui egli stesso inventa le regole. Il fisico gioca un gioco in cui le regole sono fornite dalla Natura. Ma, con il passare del tempo, diventa sempre più evidente che le regole che il matematico trova interessanti sono quelle che la Natura ha scelto. —Paul Adrien Maurice Dirac INDICE Notazione I 1 2 3 II 4 5 6 12 Analisi dello scambiatore di calore descrizione dell’impianto 14 15 modello matematico dello scambiatore di calore 2.1 Modellazione degli organi di misura e di attuazione 2.2 Analisi F.E.M 24 17 23 analisi numerica dello scambiatore di calore 29 3.1 Caratterizzazione input-output dello scambiatore di calore 29 m 3.1.1 Analisi del legame v − ϑlu 29 m 3.1.2 Analisi del legame ϑli − ϑlu 36 m 3.1.3 Analisi del legame ϑv − ϑlu 38 3.2 Analisi dell’impianto basata sul modello 42 3.2.1 Analisi del comportamento dinamico del processo attorno al punto di lavoro nominale 44 3.2.2 Analisi del comportamento dinamico del processo attorno ad alcuni punti di lavoro ammissibili 49 Controllo dello scambiatore di calore controllo dello scambiatore di calore 60 4.1 Descrizione del problema di controllo 60 4.2 Problema di regolazione asintotica 61 4.2.1 Primi passi verso la soluzione del problema di regolazione e assunzioni sul problema 62 4.2.2 Regolazione con azione feed-back+feed-forward 68 4.2.3 Controllo Observer-Based 70 4.3 Metodi basati sulla risposta armonica per la regolazione 70 4.3.1 Controllo PID based 71 4.4 Valutazione delle prestazioni del sistema di controllo in piccolo 73 sintesi dei regolatori per lo scambiatore di calore 5.1 Verifica delle condizioni per la sintesi del controllore 5.2 Controllo PID dello scambiatore di calore 78 5.3 Controllo PID con rete compensatrice 84 5.4 Controllo LQG 89 5.5 Scelta del controllore più soddisfacente 97 alcuni risultati numerici 101 6.1 Controllo PID con azione in avanti 4 102 76 77 59 6.2 6.3 7 Controllo PID con rete compensatrice con azione in avanti Controllo LQG con azione in avanti 115 conclusioni 121 a strumenti necessari per la sintesi dei regolatori bibliografia 132 122 indice 109 5 ELENCO DELLE FIGURE Figura 1 Figura 2 Figura 3 Figura 4 Figura 5 Figura 6 Figura 7 Figura 8 Figura 9 Figura 10 Figura 11 Figura 12 Figura 13 Figura 14 Figura 15 Figura 16 Scambiatore di calore 15 Tratto di condotta di lunghezza infinitesima considerato per il bilancio energetico 18 Profili di temperatura in condizioni nominali 21 Grafico della mappa della velocità all’equilibrio 22 Piano (θli θv ) dei punti di equilibrio. In blu vengono indicati i punti al di fuori della regione ammissibile mentre in rosso i punti appartenenti alla regione ammissibile 23 Caratteristica statica delle servovalvola 23 Andamento dell’indice J ( N ) 27 Profili di temperatura del liquido e del rame in condizioni nominali confrontati con le relative approssimazione ottenute dal modello ad elementi finiti 28 Risposta ad un gradino di velocità pari a 0,31 m s−1 a partire dalla condizione indicata in 3.1.1.1 30 Sollecitazione dell’impianto con un ingresso di velocità v(t) = v̄ + 0.1 sin(0.6 t) [m s−1 ] 32 Sollecitazione dell’impianto con un ingresso di velocità v(t) = v̄ + 0.2 sin(0.6 t) [m s−1 ] 33 m ( t ) a fronAndamento della variazione di temperatura ∆ ϑlu te della variazione della velocità del fluido rispetto al valore nominale: ∆ v(t) = −0.0203 δ−1 (t − 20) [m s−1 ] 33 Andamento reale ed approssimato con il modello (F.O.P.D.T.) m ( t ) rispetto al valore della variazione di temperatura ∆ ϑlu nominale a fronte della variazione della velocità del fluido rispetto al valore nominale ∆ v(t) = −0.0203 δ−1 (t − 20) [m s−1 ] 35 Andamento reale ed approssimato con il modello (S.O.R.P.P.D.T.) m ( t ) rispetto al valore della variazione di temperatura ∆ ϑlu nominale a fronte della variazione della velocità del fluido rispetto al valore nominale ∆ v(t) = −0.0203 δ−1 (t − 20) [m s−1 ] 36 Andamento reale ed approssimato con i modelli (S.O.R.P.P.D.T.) m ( t ) rie (F.O.P.D.T. della variazione di temperatura ∆ ϑlu spetto al valore nominale a fronte della variazione della velocità del fluido rispetto al valore nominale ∆ v(t) = −0.0203 δ−1 (t − 20) [m s−1 ] 36 m ( t ) a fronAndamento della variazione di temperatura ∆ ϑlu te della variazione della temperatura del fluido in ingresso rispetto al valore nominale ∆ θli = 5 δ−1 (t − 20) [◦ C ] 37 6 Figura 17 Figura 18 Figura 19 Figura 20 Figura 21 Figura 22 Figura 23 Figura 24 Figura 25 Figura 26 Figura 27 Figura 28 Figura 29 Figura 30 Figura 31 Figura 32 Figura 33 Elenco delle figure Andamento reale ed approssimato con il modello (F.O.P.D.T.) m ( t ) rispetto al valore della variazione di temperatura ∆ ϑlu nominale a fronte della variazione della temperatura del fluido in ingresso rispetto al valore nominale ∆ θli (t) = 5 δ−1 (t − 20) [◦ C ] 38 m ( t ) a fronAndamento della variazione di temperatura ∆ ϑlu te della variazione della temperatura del fluido in ingresso rispetto al valore nominale ∆ θv = 5 δ−1 (t − 20) [◦ C ] 39 Andamento reale ed approssimato con il modello (F.O.P.D.T.) m ( t ) rispetto al valore della variazione di temperatura ∆ ϑlu nominale a fronte della variazione della temperatura del fluido in ingresso rispetto al valore nominale ∆ θv (t) = 5 δ−1 (t − 20) [◦ C ] 39 Andamento reale ed approssimato con il modello (S.O.R.P.D.T.) m ( t ) rispetto al valodella variazione di temperatura ∆ ϑlu re nominale a fronte della variazione della temperatura del vapore rispetto al valore nominale ∆ θv (t) = 5 δ−1 (t − 20) [◦ C ] 40 Andamento reale ed approssimato con il modello (F.O.P.D.T.) e con il modello (S.O.R.P.P.D.T.) della variazione di temm ( t ) rispetto al valore nominale a fronte della peratura ∆ ϑlu variazione della temperatura del vapore rispetto al valore nominale ∆ θv (t) = 5 δ−1 (t − 20) [◦ C ] 41 Diagramma di Bode della funzione di trasferimento Fy v (s) 46 Diagramma di Bode della funzione di trasferimento Fy w1 (s) 46 Diagramma di Bode della funzione di trasferimento Fy w2 (s) 47 Diagramma di Bode della funzione di trasferimento Fy v (s) e della sua approssimante (3.3) 47 Diagramma di Bode della funzione di trasferimento Fy w1 (s) e della sua approssimante (3.5) 48 Diagramma di Bode della funzione di trasferimento Fy w2 (s) e della sua approssimante (3.7) 48 Diagramma di Bode della funzione di trasferimento Fy v (s) per il punto di lavoro θli = 75 ◦C θv = 140 ◦C 49 Diagramma di Bode della funzione di trasferimento Fy w1 (s) per il punto di lavoro θli = 75 ◦C θv = 140 ◦C 50 Diagramma di Bode della funzione di trasferimento Fy w2 (s) per il punto di lavoro θli = 75 ◦C θv = 140 ◦C 50 Diagramma di Bode della funzione di trasferimento Fy v (s) per il punto di lavoro θli = 50 ◦C θv = 140 ◦C 51 Diagramma di Bode della funzione di trasferimento Fy w1 (s) per il punto di lavoro θli = 50 ◦C θv = 140 ◦C 52 Diagramma di Bode della funzione di trasferimento Fy w2 (s) per il punto di lavoro θli = 50 ◦C θv = 140 ◦C 52 7 Figura 34 Figura 35 Figura 36 Figura 37 Figura 38 Figura 39 Figura 40 Figura 41 Figura 42 Figura 43 Figura 44 Figura 45 Figura 46 Figura 47 Figura 48 Figura 49 Figura 50 Figura 51 Figura 52 Figura 53 Elenco delle figure Diagramma di Bode della funzione di trasferimento Fy v (s) per il punto di lavoro θli = 75 ◦C θv = 105 ◦C 53 Diagramma di Bode della funzione di trasferimento Fy w1 (s) per il punto di lavoro θli = 75 ◦C θv = 105 ◦C 54 Diagramma di Bode della funzione di trasferimento Fy w2 (s) per il punto di lavoro θli = 75 ◦C θv = 105 ◦C 54 Diagramma di Bode della funzione di trasferimento Fy v (s) per il punto di lavoro θli = 50 ◦C θv = 105 ◦C 55 Diagramma di Bode della funzione di trasferimento Fy w1 (s) per il punto di lavoro θli = 50 ◦C θv = 105 ◦C 56 Diagramma di Bode della funzione di trasferimento Fy w2 (s) per il punto di lavoro θli = 50 ◦C θv = 105 ◦C 56 Diagrammi di Bode della funzione di trasferimento Fy v (s) per i diversi punti di lavoro. Le diverse funzioni di trasferimento nel diagramma sono indicate con Fi (s) con i ∈ {0, 1, 2, 3, 4} a seconda del punto di lavoro considerato, in maniera consistente alla numerazione riportata in (Tabella 1) 57 Diagrammi di Bode della funzione di trasferimento Fy w1 (s) per i diversi punti di lavoro. Le diverse funzioni di trasferimento nel diagramma sono indicate con Fi (s) con i ∈ {0, 1, 2, 3, 4} a seconda del punto di lavoro considerato, in maniera consistente alla numerazione riportata in (Tabella 1) 57 Diagrammi di Bode della funzione di trasferimento Fy w2 (s) per i diversi punti di lavoro. Le diverse funzioni di trasferimento nel diagramma sono indicate con Fi (s) con i ∈ {0, 1, 2, 3, 4} a seconda del punto di lavoro considerato, in maniera consistente alla numerazione riportata in (Tabella 1) 58 Schema di controllo adottato 65 Schema di controllo con azione integrale 68 Schema di controllo con azione avanti 68 Diagramma di Bode della risposta armonica del controllore PID 80 Diagramma di Bode della funzione d’anello L(s) 80 Diagramma di Nyquist della funzione d’anello L(s) 81 Diagrammi di Bode della funzione di sensitività complementare T (s) e della funzione H(s) T (s) 82 Diagramma di Bode del modulo della funzione di sensitività S(s) 82 Diagrammi di Bode delle funzioni anello ottenute per i diversi punti di lavoro. Le funzioni di trasferimento nel diagramma sono indicate con Li (s) con i ∈ {0, 1, 2, 3, 4} a seconda del punto di lavoro considerato, in maniera consistente alla numerazione riportata in (Tabella 1) 83 Diagrammi di Bode della funzione di compensazione G (s) Diagrammi di Bode del controllore C (s) G (s) 85 8 85 Figura 54 Figura 55 Figura 56 Figura 57 Figura 58 Figura 59 Figura 60 Figura 61 Figura 62 Figura 63 Figura 64 Figura 65 Figura 66 Figura 67 Figura 68 Figura 69 Figura 70 Figura 71 Figura 72 Figura 73 Figura 74 Figura 75 Figura 76 Figura 77 Elenco delle figure 9 Diagrammi di Bode della funzione di trasferimento P(s) a confronto con la funzione d’anello L(s) = P(s) C (s) G (s) 86 Diagramma di Nyquist della funzione di anello L(s) nel caso nominale 86 Diagramma di Bode della funzione sensitività nel caso nominale 87 Diagramma di Bode della funzione T (s) H(s) nel caso nominale 87 Diagrammi di Bode delle funzioni anello ottenute per i diversi punti di lavoro. Le funzioni di trasferimento nel diagramma sono indicate con Li (s) con i ∈ {0, 1, 2, 3, 4} a seconda del punto di lavoro considerato, in maniera consistente alla numerazione riportata in (Tabella 1) 88 Valori singolari della funzione di trasferimento del processo linearizzato attorno al punto di lavoro nominale 90 Diagrammi di Bode della funzione di trasferimento del processo nominale e della sua approssimazione di ordine 6 91 Diagrammi di Bode della funzione di trasferimento del compensatore LQG con azione integrale 92 Valori singolari associati ai modi naturali presenti nel legame ingresso-uscita del controllore LQG 93 Diagramma di Bode di Cr (s) e di C (s) a confronto 93 Diagramma di Bode della funzione di trasferimento Cd (s) 94 Diagramma di Bode della funzione di trasferimento L(s) 94 Diagramma di Bode del modulo della funzione di sensitività S(s) 95 Diagramma di Bode della funzione di trasferimento H(s) T (s) 95 Diagramma di Nyquist della funzione d’anello L(s) 96 Diagrammi di Bode delle funzioni anello ottenute per i diversi punti di lavoro. Le funzioni di trasferimento nel diagramma sono indicate con Li (s) con i ∈ {0, 1, 2, 3, 4} a seconda del punto di lavoro considerato, in maniera consistente alla numerazione riportata in (Tabella 1) 96 Diagrammi di Bode del modulo della funzione di sensitività per i tre controllori nel caso nominale 98 Diagrammi di Bode della funzione di trasferimento tra la variazione del riferimento di temperatura e l’errore di regolazione H(s) T (s) per i tre controllori nel caso nominale 99 Diagrammi di Bode delle funzioni di trasferimento C (s) e Cd (z) a confronto 100 Andamento del riferimento y∗ (t) 101 Andamento delle variabili ϑli (t) e ϑv (t) 101 Andamento dell’uscita controllata y(t) 102 Andamento del controllo u(t) 102 Variazione delle temperatura di riferimento da 98 ◦ a 94 ◦ 103 Figura 78 Figura 79 Figura 80 Figura 81 Figura 82 Figura 83 Figura 84 Figura 85 Figura 86 Figura 87 Figura 88 Figura 89 Figura 90 Figura 91 Figura 92 Figura 93 Figura 94 Variazione delle temperatura di riferimento da 94 ◦ a 98 ◦ 104 Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (65◦ , 120◦ ) a (75◦ , 140◦ ) 105 Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (75◦ , 140◦ ) a (50◦ , 140◦ ) 106 Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (50◦ , 140◦ ) a (75◦ , 105◦ ) 107 Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (75◦ , 105◦ ) a (50◦ , 105◦ ) 108 ◦ Variazione delle temperatura di riferimento da 98 a 94 ◦ 109 Variazione delle temperatura di riferimento da 94 ◦ a 98 ◦ 110 Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (65◦ , 120◦ ) a (75◦ , 140◦ ) 111 Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (75◦ , 140◦ ) a (50◦ , 140◦ ) 112 Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (50◦ , 140◦ ) a (75◦ , 105◦ ) 113 Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (75◦ , 105◦ ) a (50◦ , 105◦ ) 114 ◦ Variazione delle temperatura di riferimento da 98 a 94 ◦ 115 Variazione delle temperatura di riferimento da 94 ◦ a 98 ◦ 116 Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (65◦ , 120◦ ) a (75◦ , 140◦ ) 117 Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (75◦ , 140◦ ) a (50◦ , 140◦ ) 118 Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (50◦ , 140◦ ) a (75◦ , 105◦ ) 119 Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (75◦ , 105◦ ) a (50◦ , 105◦ ) 120 E L E N C O D E L L E TA B E L L E Tabella 1 Tabella 2 Valore della temperatura del liquido in ingresso e del vapore per i punti di lavoro scelti. Il punto di lavoro nominale viene indicato con il numero 0. 44 Valore dei margini di guadagno e di fase ottenuti per i punti di lavoro scelti con il controllo PID. Il punto di lavoro nominale viene indicato con il numero 0. 83 10 Tabella 3 Tabella 4 Tabella 5 Elenco delle tabelle Valore dei margini di guadagno e di fase ottenuti per i punti di lavoro scelti con il controllo PID+ Rete compensatrice. Il punto di lavoro nominale viene indicato con il numero 0. 88 Valore dei margini di guadagno e di fase ottenuti per i punti di lavoro scelti con il controllo LQG. Il punto di lavoro nominale viene indicato con il numero 0. 92 Valore dei margini di guadagno e di fase ottenuti per i punti di lavoro scelti con il controllo LQG. Il punto di lavoro nominale viene indicato con il numero 0. 98 11 N OTA Z I O N E ( δ−1 (t − t0 ) := 1 0 t ≥ t0 t < t0 Funzione gradino C Insieme dei numeri complessi C+ Insieme dei numeri complessi con parte reale non negativa R Insieme dei numeri reali R+ Insieme dei numeri reali non negativi Rn Vettori reali di dimensione n Rn × m Matrici reali di dimensione n × m In || u|| s Matrice identità di dimensione n √ u T u Norma euclidea Variabile complessa del dominio di Laplace T.H.D.% Distorsione armonica percentuale (Total Harmonic Distortion) x T , ( A T ) Vettore x trasposto (matrice A trasposta) C k ( D ) Insieme delle funzioni continue e con le prime k derivate continue nell’insieme D L∞ Insieme delle funzioni reali limitate 12 INTRODUZIONE I n questa trattazione si vuole affrontare il problema della analisi e del controllo di uno scambiatore di calore. La trattazione verrà affrontata dividendo il testo in due parti: La prima relativa all’analisi del processo considerato mentre la seconda dedicata alla sintesi del sistema di controllo. In particolare nella prima parte verrà derivato un modello matematico del processo basato sulle equazioni che governano il fenomeno fisico portante, successivamente su tale modello verrano condotte delle analisi atte a caratterizzare la dinamica del processo fisico e il comportamento ingresso uscita del medesimo. Nella seconda parte invece verrà formalizzato il problema di controllo dello scambiatore di calore e verranno presentati dei risultati che condurranno alla soluzione del problema. Infine si mostreranno i risultati ottenuti dal controllo dello scambiatore di calore. 13 Parte I. Analisi dello scambiatore di calore 14 1 I D E S C R I Z I O N E D E L L’ I M P I A N T O n termotecnica con il termine scambiatore di calore si fa riferimento ad un componente atto a realizzare uno scambio di energia termica tra due fluidi a temperature diverse. In generale gli scambiatori sono sistemi aperti che operano senza scambio di lavoro, ovvero presentano un flusso costante di fluido e una distribuzione di temperatura a regime costante. Nella realtà industriale esistono molteplici metodologie per effettuare lo scambio termico tra due fluidi ossia esistono diverse tecnologie per realizzare uno scambiatore di calore. Nella trattazione che ci appresta a presentare verrà presa in considerazione una particolare tipologia di scambiatore di calore detto Scambiatore a superficie. In questo tipo di scambiatore le correnti assorbono il calore dalle superfici con le quali sono separati, inoltre in tali dispositivi si riconoscono sempre due lati contenenti i fluidi. In particolare la struttura che verrà presa in considerazione in questo ambito è mostrata in (Figura 1). In questa struttura il fluido ad Figura 1.: Scambiatore di calore un temperatura ϑlim , acqua nel caso in esame, scorre con velocità v all’interno di una serpentina in rame la quale è immersa in una camera contenente del vapore saturo a temperatura ϑv ed infine il fluido sbocca dalla serpentina con una temperatura ϑlu . In particolare nell’impianto in esame la velocità del fluido v può essere variata, agendo sulla apertura di una servovalvola calettata sulla serpentina stessa, mentre la temperatura del fluido in uscita ϑlu può essere misurata attraverso una termocoppia installata allo sbocco della condotta dello scambiatore di calore. E’ evidente che a seconda del valore delle variabili appena menzionate l’impianto si troverà in una condizione operativa piuttosto che in un altra. Ebbene tra tutte le possibili condizioni operative in cui l’impianto può trovarsi si può prendere in considerazione il punto di lavoro riferito ai 15 valori nominali delle variabili in gioco, ossia: descrizione dell’impianto • ϑlu = 98◦ C • ϑv = 120◦ C • ϑli = 65◦ C tale punto viene denominato punto di lavoro nominale e verrà considerato come il punto di lavoro di riferimento per condurre la maggior parte delle analisi sullo scambiatore di calore. 16 2 A M O D E L LO M AT E M AT I C O D E L LO S C A M B I ATO R E D I C A LO R E l fine di poter effettuare un analisi numerica dello scambiatore di calore è necessario derivare per tale sistema un modello matematico sufficientemente affidabile. Prima di passare alla scrittura del modello matematico del processo in esame è bene fare alcune assunzioni semplificative [3]. • La temperatura del vapore saturo è uguale in ogni punto della camera • Il fluido è assunto essere incompressibile • Il fluido è assunto essere in regime laminare • La serpentina è considerata come un tubo rettilineo Si noti che l’aver assunto il fluido in regime laminare impone un vincolo sulla velocità v dello stesso, la quale ovviamente non può assumere valori qualunque. In particolare tale assunzione implica che: v ∈ [0.1, 0.7] m s−1 Inoltre si supporrà nel seguito della trattazione che: ϑli ≤ 80◦ C (2.1) ◦ ϑli ≥ 100 C (2.2) A fronte delle ipotesi semplificative appena esposte, al fine di formulare un modello matematico del processo in esame, si può effettuare un bilancio delle energie scambiate da una sezione di lunghezza infinitesima della serpentina in un intervallo temporale infinitesimo. A tal scopo è conveniente dapprima fissare una ascissa di riferimento, che indicheremo con ξ, in modo da poter individuare in maniera univoca ogni punto della serpentina. Detto ciò si indichino ora con ϑl (ξ, t) e con ϑr (ξ, t) rispettivamente, la temperatura del generico punto del liquido e del rame nell’istante di tempo t e con ϑv (t) la temperatura del vapore nella camera nell’istante di tempo t, dunque il bilancio energetico può essere effettuato considerando dapprima i flussi di energia entranti ed uscenti, dalla porzione di liquido contenuta nel tratto di condotta di lunghezza infinitesima in esame e successivamente prendendo in considerazione i medesimi flussi relativamente alle pareti della condotta (Figura 2). In particolare relativamente ai flussi di energia riguardanti il liquido, si possono mettere in evidenza principalmente tre forme di scambio termico: • Energia che il liquido scambia al suo interno • Energia scambiata con la parete della serpentina • Energia che il liquido trasporta a causa del suo moto nella serpentina 17 modello matematico dello scambiatore di calore Figura 2.: Tratto di condotta di lunghezza infinitesima considerato per il bilancio energetico Ebbene indicando con: • αlm coefficiente di adduzione acqua-rame [ J m−2 s−1 C −1 ] • cl calore specifico dell’acqua [ J Kg−1 C −1 ] • Dl diametro interno della condotta [m] • µl densità lineare dell’acqua [Kg m−1 ] • v velocità del fluido nella condotta [m s−1 ] gli scambi termici sopra elencati in maniera qualitativa possono essere espressi in modo quantitativo come segue: • Energia che il liquido scambia al suo interno −cl µl dξ (ϑl (ξ, t + dt) − ϑl (ξ, t)) • Energia scambiata con la parete della serpentina αlm dξ (ϑr (ξ, t) − ϑl (ξ, t))π Dl dξ dt • Energia che il liquido trasporta a causa del suo moto nella serpentina µl v dt cl (ϑl (t, ξ ) − ϑl (t, ξ + dξ )) Viceversa relativamente ai flussi di energia riguardanti le pareti della condotta si possono mettere in evidenza, anche in questo caso, tre forme di scambio termico: • Energia che la parete scambia con se stessa • Energia che la parete scambia con il liquido • Energia che la parete scambia con il vapore In questo caso indicando con: 18 modello matematico dello scambiatore di calore • αv coefficiente di adduzione vapore-rame [ J m−2 s−1 C −1 ] • cr calore specifico del rame [ J Kg−1 C −1 ] • Dv diametro interno della condotta [m] • µr densità lineare del rame [Kg m−1 ] tali termini possono essere quantificati come segue: • Energia che la parete scambia con se stessa −µr dξ cr (ϑr (ξ, t + dt) − ϑr (ξ, t)) • Energia che la parete scambia con il liquido −αlm π Dl dξ dt (ϑr (ξ, t) − ϑl (ξ, t)) • Energia che la parete scambia con il vapore αv π Dv ξ dt (ϑr (ξ, t) − ϑv (t)) A questo punto imponendo la conservazione dell’energia si ottengono due bilanci energetici uno per il liquido ed uno per il rame della condotta. Bilancio energetico del liquido: − cl µl dξ (ϑl (ξ, t + dt) − ϑl (ξ, t)) + αlm dξ (ϑr (ξ, t) − ϑl (ξ, t))π Dl dξ dt+ + µl v dt c1 (ϑl (t, ξ ) − ϑl (t, ξ + dξ )) = 0 (2.3) Bilancio energetico del rame: µr dξ cr (ϑr (ξ, t + dt) − ϑr (ξ, t)) + αlm π Dl dξ dt (ϑr (ξ, t) − ϑl (ξ, t))+ + αv π Dv ξ dt (ϑr (ξ, t) − ϑv (t)) = 0 (2.4) Pertanto al fine di ottenere delle equazioni di validità locale si può effettuare un passaggio al limite e confondere gli incrementi parziali delle varie temperature con le relative derivate prime. In tal modo si ottiene il seguente sistema di equazioni alle derivate parziali (E.D.P.): ∂ϑr (ξ, t) αv (ϑv (t) − ϑr (ξ, t)) π Dv − αlm π Dl (ϑr (ξ, t) − ϑl (ξ, t)) − µr cr =0 dt ∂(ϑl (ξ, t)) ∂ϑ (ξ, t) − cl µl + αlm π Dl (ϑr (ξ, t) − ϑl (ξ, t)) + µl v c1 l =0 dt dξ ϑl (0, t) = ϑli (t) ∀t ≥ 0 ϑr (ξ, 0) = ϑr ( x ) ∀ ξ ∈ [0 , L] ϑl (ξ, 0) = ϑl ( x ) ∀ ξ ∈ (0 , L] 19 modello matematico dello scambiatore di calore (2.5) Il significato delle boundary-conditions alle quali è soggetto il sistema di (E.D.P.) è presto detto: La temperatura del liquido in ingresso deve poter essere specificata istante per istante mentre i profili spaziali delle temperature del rame e del liquido nell’istante iniziale devono essere specificate al fine di ottenere un problema ben posto. Fatte le necessarie precisazioni si può riscrivere il sistema di (E.D.P.) (2.5) in maniera più compatta come segue: ∂ϑr (ξ, t) =0 (ϑv (t) − ϑr (ξ, t)) − β (ϑr (ξ, t) − ϑl (ξ, t)) − τ2 dt ∂ϑ (ξ, t) ∂(ϑl (ξ, t)) + (ϑr (ξ, t) − ϑl (ξ, t)) + v τ1 l =0 − τ1 dt dξ ϑl (0, t) = ϑli (t) ∀t ≥ 0 ϑr (ξ, 0) = ϑr ( x ) ∀ ξ ∈ [0 , L] ϑl (ξ, 0) = ϑl ( x ) ∀ ξ ∈ (0 , L] (2.6) Dove le variabili τ1 , τ2 , β sono espresse come segue: β := Dl αlm Dv α v (2.7) (2.8) τ1 := cl µl [s] π Dl αlm (2.9) (2.10) cr µr τ2 := [s] π Dv α v (2.11) Ottenuta una scrittura più snella per il sistema (E.D.P.) (2.5) si può passare ora ad ottenere una soluzione analitica di tipo stazionario di tale sistema al fine di poter valutare il profilo di temperatura che si stabilisce lungo la condotta a transitorio esaurito a fronte di una una temperatura costante per il liquido in ingresso ϑli (t) ≡ ϑli0 Ebbene indicando con (ϑr (ξ ) , ϑl (ξ ) la soluzione cercata, allora essa certamente deve soddisfare il seguente sistema di equazioni misto Algebrico-E.D.O: 20 modello matematico dello scambiatore di calore (ϑv (t) − ϑr (ξ )) − β (ϑr (ξ ) − ϑl (ξ )) = 0 (ϑ (ξ ) − ϑ (ξ )) + v τ ∂ϑl (ξ ) = 0 r 1 l dξ 0 ϑl (0, t) = ϑli ∀t ≥ 0 (2.12) A fronte di opportune manipolazione e di una successiva integrazione rispetto alla variabile spaziale ξ si ottengono le seguenti espressioni per la soluzione cercata: − τ v (xβ+1) 1 ϑ ( ξ ) = e (ϑli0 − ϑv ) + ϑv l ϑr (ξ ) = β β+1 e x 1 v ( β +1) −τ (2.13) ! (ϑli0 − ϑv ) + ϑv +ϑv Attraverso l’uso delle relazioni ottenute è possibile graficare i profili di temperatura del liquido e del rame che ci stabiliscono in condizioni stazionarie, mostrati in (Figura 3) 120 115 110 105 100 θ(ξ ) ◦ C 95 90 85 80 75 θ r (ξ ) 70 θ l (ξ ) 65 60 0 0.5 1 1.5 2 2.5 ξ [m] Figura 3.: Profili di temperatura in condizioni nominali In particolare essendo la temperatura del fluido in uscita la grandezza di maggior interesse è conveniente disporre di un espressione di tale valore dipendente dagli altri parametri di processo, tale espressione può essere ottenuta a partire dalla (2.13) valutando ϑl (ξ ) nel punto ξ = L, l’espressione che si ottiene è la seguente: ϑ ( L) = ϑlu = e L 1 v ( β +1) −τ (ϑli0 − ϑv ) + ϑv 21 modello matematico dello scambiatore di calore Dal legame appena mostrato è possibile derivare il valore della velocità del fluido necessario ad ottenere una data temperatura per il liquido in uscita dalla condotta, in particolare la relazione suddetta è la seguente: v= L (2.14) ϑv τ1 ( β + 1) ln( ϑϑluli − − ϑv ) il legame appena pervenuto è di fondamentale importanza in quanto permette di stabilire, a partire dalle condizioni operative del processo, se sia possibile o meno condurre l’impianto in un equilibrio assegnato, ossia permette di stabilire se il valore di velocità da imporre al fluido in una data condizione operativa dell’impianto si ammissibile o meno. A titolo di esempio si consideri il caso in cui si voglia mantenere la temperatura del fluido in uscita al suo valore nominale ossia 98◦ C, anche a fronte di variazioni della temperatura del fluido in ingresso e della temperatura del vapore, ebbene è evidente che esisteranno delle condizioni operative per le quali il valore della velocità richiesto non rientra nei vincoli imposti dalla assunzione di regime laminare. Per ragioni che saranno chiare più avanti è conveniente definire l’insieme ammissibile dei punti di equilibrio ammissibili in maniera precisa. Definizione 2.1: Si definisce insieme ammissibile delle condizioni operative Θϑlu l’insieme definito come segue: ( ) L Θϑlu : (ϑli , ϑv ) ∈ R2 : 0.1 ≤ ≤ 0.7 ϑv τ1 ( β + 1) ln( ϑϑluli − − ϑv ) A supporto di quanto appena illustrato, in (Figura 4), viene riportato il grafico della relazione (2.14) ottenuto facendo variare ϑv e ϑli in un intervallo opportuno e fissando ϑlu al suo valore nominale ossia 98 [◦ C ] Inoltre per rendersi con- 0.7 0.8 0.7 0.6 0.6 v [m s−1] 0.5 0.5 0.4 0.4 0.3 0.2 0.3 0.1 0 120 120 100 100 60 60 θli[°C] 0.2 80 80 40 40 20 20 0 0 θv[°C] 0.1 Figura 4.: Grafico della mappa della velocità all’equilibrio to di quale sia regione ammissibile, nel caso di temperatura del fluido di uscita nominale, in (Figura 5) vengono discriminati, tramite un diversa colorazione, i punti appartenenti alla regione ammissibile Θ98 dai restanti. Prima di concludere questa sezione, si osservi che la relazione appena discussa può essere utilizzata per derivare il valore della velocità del fluido in condizioni nominali v̄ (velocità nominale del fluido), il quale a priori risulta non noto, 22 2.1 modellazione degli organi di misura e di attuazione 140 135 130 θv [ ◦C ] 125 120 115 110 105 100 50 55 60 65 70 75 80 θl i [ ◦C ] Figura 5.: Piano (θli θv ) dei punti di equilibrio. In blu vengono indicati i punti al di fuori della regione ammissibile mentre in rosso i punti appartenenti alla regione ammissibile ebbene a riguardo la ( 2.14 nella pagina precedente) porge il seguente valore per la velocità nominale: v̄ ≈ 0,31 m s−1 2.1 modellazione degli organi di misura e di attuazione Come già indicato espressamente nella situazione in esame sono presenti degli organi di misura e degli organi di attuazione, in particolare sono presenti una termocoppia per la misura della temperatura del fluido in uscita e una servovalvola per la regolazione della velocità del fluido nella serpentina. Tali dispositivi debbono essere modellati oculatamente al fine di ottenere un modello accurato del processo sotto controllo. In particolare per quanto riguarda il dispositivo di attuazione della valvola è possibile, data la sua elevata rapidità, trascurare la sua presenza e tenere in conto solamente dell’intervallo delle velocità ammissibili per il fluido modellando la servovalvola come un dispositivo non lineare statico descritto dalla caratteristica in (Figura 6). Relativamente alla termocop- Figura 6.: Caratteristica statica delle servovalvola pia, è possibile modellare tale dispositivo come un sistema a singola costante m la misura della temperatura del liquido di tempo. Con ciò indicando con ϑlu 23 2.2 analisi f.e.m in ingresso fornita dalla termocoppia, con z lo stato associato a tale dispositivo e scegliendo ragionevolmente per la costante di tempo del dispositivo il valore di 1 [s], si ottiene il modello del dispositivo suddetto: ż = −z + ϑlu m ϑlu 2.2 (2.15) =z (2.16) analisi f.e.m Il modello ottenuto, essendo descritto da un set di (E.D.P.), è di fatto un modello a dimensione infinita, dunque per poter effettuare una simulazione numerica del comportamento dello scambiatore di calore si deve realizzare una analisi di tipo (F.E.M), ossia si effettua una discretizzazione spaziale (mesh) sufficientemente fitta della condotta, ottenendo in tale modo un set di equazioni di (E.D.O.) ognuna delle quali rappresenta l’evoluzione della temperatura di ognuno dei campioni spaziali (elementi finiti) ottenuti a fronte dell’operazione di discretizzazione. In letteratura esistono molte metodologie per effettuare l’operazione di discretizzazione del dominio spaziale per il quale si sta cercando la soluzione del sistema di (E.D.P.). Tuttavia in questa sede si utilizzerà il metodo più semplice di essi il quale, nel caso in esame, prevede di suddividere il dominio spaziale in N elementi uguali fra loro di lunghezza l= L N e di considerare come stato la temperatura del liquido e del rame per ognuno dei segmenti di condotta considerati. In tal modo si è implicitamente assunto che la temperatura all’interno del generico elemento finito sia omogenea, tale assunzione ovviamente introduce un errore che diminuisce all’aumentare del numero di elementi che compongono la mesh. Al fine di ottenere il modello (F.E.M) è necessario definire la mesh, ossia individuare gli elementi che comporranno la medesima. Tra tutte le scelte, quella che consente di ottenere un modello ad elementi finiti per il quale sia possibile specificare istante per istante la temperatura del liquido ingresso è la seguente: Elementi della Mesh relativa al liquido: ϑl1 := ϑl (l ) ϑl2 := ϑl (2 l ) ... ϑ N := ϑ ( N l ) l l Elementi della Mesh relativa al rame: ϑr1 := ϑr (l ) ϑr2 := ϑr (2 l ) ... ϑ N := ϑ ( N l ) r r (2.17) (2.18) 24 2.2 analisi f.e.m 25 A fronte della scelta fatta per la mesh, approssimando le derivate parziali rispetto alla variabile spaziale ξ con i relativi rapporti incrementali all’indietro1 , presenti nelle (2.6) si ottiene il seguente sistema (2.19) di 2 N (O.D.E) seguente, il quale rappresenta l’approssimazione (F.E.M) del sistema (2.6): dϑri (t) i i i =0 ( ϑ ( t ) − ϑ ( t ) − β ( ϑ ( t ) − ϑ ( t )) − τ v 2 r r l dt ϑli (t) − ϑli−1 (t) d(ϑli (t)) i i − τ + ( ϑ =0 ( t ) − ϑ ( t )) + v τ 1 1 r l dt l (2.19) i ϑri (0) = ϑr0 i ϑli (0) = ϑl0 ∀ i ∈ {1, . . . N } Tale sistema può essere riscritto in forma matriciale come segue: ˙ ϑl1 (t) .. . ˙N ϑl ( t ) ˙1 ϑr (t) .. . ˙ N ϑr (t) = −( vl + v l 0 β τ2 0 0 v l + 0 0 0 .. . 0 1 τ1 ) 0 0 .. . 0 .. v . −( l + 0 .. . 0 0 0 β τ2 1 τ1 1 τ1 ) 0 0 0 .. . 0 0 1 τ1 0 .. . 0 −( τβ2 + 0 0 1 τ2 ) 0 0 −( τβ2 0 + 1 ϑl ( t ) .. . ϑlN (t) + ϑ1 (t ) r .. . 1 N ϑr (t) τ2 ) 0 .. . 0 ϑli 1 ϑv τ2 .. . 1 τ2 (2.20) Per concludere, al fine di ottenere un modello che possieda come uscita la temm , è necessario introdurre la dinamica peratura misurata del liquido in uscita ϑlu della termocoppia ( 2.15 nella pagina precedente) caratterizzata dalla variabile di stato z, ottenendo il modello dinamico seguente: 1 Utilizzando un approssimazione basata su rapporto incrementale in avanti, scrivendo l’equazione differenziale che descrive l’evoluzione della temperatura del liquido nel primo tratto di condotta, si introdurrebbe un vincolo sulla derivata della temperatura del liquido in ingresso, contraddicendo l’ipotesi di arbitrarietà di tale variabile. 2.2 analisi f.e.m ˙1 ϑl ( t ) .. . ϑ˙N (t) l ˙1 ϑr (t) .. . ϑ˙rN (t) ż = −( vl + v l 0 β τ2 0 0 0 v l + 0 0 0 .. . 0 1 τ1 ) 0 .. . .. . 0 .. . 0 ... 1 τ1 0 0 −( vl + 0 1 τ1 ) 0 0 .. . 0 0 0 1 τ1 0 .. . 0 −( τβ2 + 0 β τ2 1 0 0 0 1 τ2 ) 0 ... 0 0 + −1 −( τβ2 1 ) τ2 26 ϑl1 (t) .. . ϑlN (t) + 1 ϑr (t) .. . ϑrN (t) z 0 .. . 0 ϑli 1 ϑv τ2 .. . 1 τ2 m ϑlu =z (2.21) Stabilita la struttura del sistema ad elementi finiti, al fine di poter utilizzare il modello per le relative simulazioni si deve stabilire quale sia il numero di elementi finiti opportuno nel caso in esame. Con lo scopo di ottenere un buon modello è conveniente fare fede ad un opportuno criterio di ottimo minimizzando una particolare cifra di merito che tenga conto dell’errore commesso nella discretizzazione del problema. Ebbene a fronte di numerose analisi statistiche è stato osservato che una buona cifra di merito è data dalla seguente espressione: J(N) = 1 2N N ∑ k =1 c (ϑlk − ϑl (k l ))2 + (1 − c)(ϑrk − ϑr (k l ))2 (2.22) In particolare scegliendo c = 0.8 si attribuisce un peso maggiore all’errore commesso sull’approssimazione del profilo ti temperatura del liquido. Ebbene osservando l’andamento dell’indice J, al variare del numero di elementi N, mostrato in (Figura 7) si nota che scegliendo N = 35 si ottiene un buon compromesso tra precisione del modello (almeno in condizioni statiche) e efficienza computazionale. A posteriori delle scelta fatta per N si ottiene il modello desiderato nella forma (2.21), nello specifico ponendo x := [ϑl1 , . . . , ϑl35 , . . . . . . , ϑr1 , . . . , ϑr35 , z] T (2.23) w := [ϑli , ϑv ] T 2.2 analisi f.e.m 100 90 80 70 J (N ) 60 50 40 30 20 10 0 0 10 20 30 40 50 60 70 80 90 100 N Figura 7.: Andamento dell’indice J ( N ) il modello ottenuto può essere riscritto nella forma: ẋ = A(v) x + B(v) w (2.24) y = Cx (2.25) avendo cura di indicare rispettivamente: −( vl + 1 τ1 ) v l A(v) = 0 β τ2 0 .. . 0 0 0 v l B(v) = 0 0 0 .. . 0 0 .. . .. . 0 ... 0 .. . 0 1 τ2 .. . 0 1 τ1 0 0 −( vl + 0 0 β τ2 1 1 τ1 ) 0 −( τβ2 + 0 0 0 0 .. . 0 1 0 τ1 1 ) 0 0 τ2 .. . 0 β 1 0 −( τ2 + τ2 ) ... −1 0 (2.26) 1 τ2 C = 0 ... 0 0 ... 1 Dunque al fine di valutare l’accuratezza del modello ad elementi finiti, si può determinare quale sia il profilo di temperatura che si stabilisce in regime stazionario ottenuto da tale modello. A tal scopo si può procedere come segue: 1. Fissate ϑli ,ϑv e ϑlu si determina il valore v̄ della velocità del fluido all’ equilibrio tramite la relazione (2.14) 27 2.2 analisi f.e.m 2. Si determina il valore dello stato di equilibrio x̄ risolvendo la seguente equazione lineare nell’incognita x̄ A(v̄) x̄ + B(v̄) w = 0 La procedura appena descritta può essere utilizzata nelle condizioni nominali, ossia: ϑli = 65◦ C ϑlu = 98◦ C (2.27) ◦ ϑv = 120 C ottenendo il profilo mostrato in (Figura 8). Quanto appena mostrato di fatto conclude la ricerca del modello ad elementi finiti necessario alla simulazione dello scambiatore di calore. 100 θ l (k l ) θ l (ξ ) 95 θ l (ξ ) ◦ C 90 85 80 75 70 65 0 0.5 1 1.5 2 2.5 ξ [m] (a) Profilo di temperatura del liquido 117 θ r (ξ ) θ r (k l ) 116 θ r (ξ ) ◦ C 115 114 113 112 111 0 0.5 1 1.5 2 2.5 ξ [m] (b) Profilo di temperatura del rame Figura 8.: Profili di temperatura del liquido e del rame in condizioni nominali confrontati con le relative approssimazione ottenute dal modello ad elementi finiti 28 3 N ANALISI NUMERICA DELLO S C A M B I ATO R E D I C A LO R E el capitolo precedente è stato ottenuto un modello ad elementi finiti in grado di descrivere in maniera accurata lo scambiatore di calore. In questo capitolo si effettueranno delle analisi sul modello matematico ad elementi finiti dello scambiatore di calore ottenuto in precedenza, atte a mettere in luce quali siano le caratteristiche dello scambiatore stesso dal punto di vista ingresso uscita. Tale caratterizzazione di fatto è necessaria ogni qual volta si voglia effettuare il controllo di un determinato processo. 3.1 caratterizzazione input-output dello scambiatore di calore Al fine di poter ottenere una caratterizzazione del processo in grado di fornire delle indicazioni sulla sintesi di una legge di controllo è necessario effettuare sul processo stesso delle analisi atte a mettere in evidenza quale sia la natura dei legami input-output in questione. Prima di ogni altra cosa è bene fare chiarezza su quale sia il ruolo delle diverse variabili presenti nel modello; pertanto considerando che in precedenza è stato indicato espressamente che l’unica variabile direttamente manipolabile è la velocità del fluido v, è evidente che tale variabile è la variabile di controllo del processo, a questo punto non resta allora che mettere in evidenza quale sia il ruolo delle restanti variabili di ingresso. Ebbene è conveniente attribuire il ruolo di variabili di disturbo alle grandezze ϑv , ϑli , infatti per esse non è possibile fissare arbitrariamente un valore. Avendo chiarito quale sia il ruolo delle diverse variabili si può passare ora ad analizzare le risposte del processo a fronte delle varie sollecitazioni rispetto ai diversi ingressi, il tutto sempre tenendo presente che in alcuni casi si analizzerà la risposta all’ingresso di controllo mentre in altri quello che si otterrà è la risposta a delle azioni disturbanti; tale circostanza è necessaria ad attribuire un peso diverso ai risultati che si otterranno nelle diverse analisi. Prima di procedere con le suddette analisi si tenga presente che ognuna di esse verrà effettuata a partire da diverse condizioni operative dell’impianto, condizioni che si ottengono partendo da una opportuno stato iniziale e mantenendo costanti di volta in volta alcune delle variabili di ingresso. m 3.1.1 Analisi del legame v − ϑlu Il primo legame che vuole studiare è quello tra la variabile di uscita e la variabile di controllo, certamente tale legame è molto importante in quanto esso influenza la procedura di sintesi della legge di controllo che verrà mostrata più avanti. 29 3.1 caratterizzazione input-output dello scambiatore di calore 3.1.1.1 Conduzione dell’impianto nello stato nominale La prima delle analisi che si vogliono effettuare consiste nel condurre l’impianto da una condizione operativa iniziale al punto di lavoro nominale. Per effettuare tale analisi si sollecita l’impianto con un ingresso di velocità del liquido v pari al valore nominale ossia 0.31 m s−1 a partire dalle seguenti condizioni iniziali: • Temperatura del rame: 90◦ C • Temperatura del liquido: 65◦ C • Temperatura del vapore: 120◦ C tali condizioni operative, significative dal punto di vista fisico, si traducono in dei valori opportuni dello stato iniziale del processo e dei valori temporali dei vari ingressi. In particolare si ha: • xi (0) = 90◦ C ∀ i ∈ {36, . . . , 70} elementi della mesh del rame • xi (0) = 65◦ C ∀ i ∈ {1, . . . , 35} elementi della mesh le liquido • x71 (0) = 65◦ C variabile di stato della termocoppia • ϑli (t) = 65◦ C ∀ t ≥ 0 • ϑv (t) = 120◦ C ∀ t ≥ 0 m ( t ) a fronte dell’ingresso di velocità del valore L’andamento della variabile ϑlu pari al valore nominale, è riportato in (Figura 9). Dall’andamento in (Figura 9) 98 94 90 θ lmu [ ◦ C ] 86 82 78 74 70 66 62 0 5 10 15 time [s] Figura 9.: Risposta ad un gradino di velocità pari a 0,31 m s−1 a partire dalla condizione indicata in 3.1.1.1 si può notare che la temperatura del fluido in uscita non raggiunge esattamente il valore di 98 ◦C in quanto è presente l’errore introdotto dalla approssimazione agli elementi finiti, si può inoltre notare che l’andamento non presenta né ritardi né oscillazioni. 30 3.1 caratterizzazione input-output dello scambiatore di calore Giunti a questo punto si vogliono effettuare ora delle analisi volte a mettere in luce quale sia il comportamento del sistema a fronte di piccole variazioni, del valore della velocità del liquido in ingresso, intorno al proprio valore nominale. Ebbene questo tipo di analisi è molto importante dal punto di vista controllistico in quanto, come ben noto, sollecitando un qualunque sistema dinamico non lineare con degli ingressi di ampiezza opportunamente piccola è possibile ottenere in uscita una risposta principalmente caratterizzata dalla componente lineare del sistema stesso, ottenendo in tal modo seppur in maniera qualitativa, delle indicazioni sulle difficoltà che si riscontreranno nella sintesi di una opportuna legge di controllo. Quanto appena detto di fatto può essere osservato conducendo due tipi di analisi. La prima analisi consiste nel valutare la distorsione armonica percentuale (T.H.D. %) della risposta del sistema, a fronte di una variazione sinusoidale di ampiezza opportuna della variabile v, se il sistema fosse lineare la distorsione armonica della risposta in uscita sarebbe chiaramente nulla, ovvero conducendo questa analisi si può ottenere una misura quantitativa dell’influenza dei termini non lineari, nella risposta del sistema. Per meglio comprendere quale sia il comportamento del sistema, l’analisi suddetta viene condotta per diversi valori di pulsazione dell’ingresso e per diversi valori di ampiezza del medesimo. Si riportano i risultati in (Figura 10) e (Figura 11). Osservando le figure di cui sopra si nota che all’aumentare dell’ampiezze dell’ingresso la distorsione armonica dell’uscita aumenta, tuttavia la distorsione risulta comunque contenuta. 31 3.1 caratterizzazione input-output dello scambiatore di calore ∆θ lmu(t ) [ ◦ C ] 1 0 −1 −2 10 15 20 25 30 35 40 30 35 40 time [s] v (t ) [ m s −1] 0.4 0.35 0.3 0.25 0.2 0.15 0.1 10 15 20 25 time [s] (a) Andamento della velocità del fluido e della relativa variazione di temperatura del liquido in uscita a partire da condizioni nominali T H D = 9.00% 1.6 1.4 1.2 Mag 1 0.8 0.6 0.4 0.2 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Harmonic order (b) Analisi armonica dell’andamento della variazione di temperatura del liquido in uscita Figura 10.: Sollecitazione dell’impianto con un ingresso di velocità v(t) = v̄ + 0.1 sin(0.6 t) [m s−1 ] In virtù della scarsa influenza dei termini non lineari nell’uscita dell’impianto, mostrata attraverso l’analisi appena condotta, risulta significativo ottenere una buona approssimazione lineare dell’impianto in grado di predire in maniera sufficientemente accurata l’evoluzione temporale della variazione rispetm ( t ) a fronte di piccole to al valore nominale, della temperatura in uscita ∆ϑlu variazioni della velocità del liquido in ingresso. Una prima metodologia che può essere utilizzata a tale scopo si basa sullo studio delle risposta indiciale dell’impianto ottenuta per variazioni a gradino della velocità del fluido di ampiezza modesta, a titolo d’esempio si osservi la risposta mostrata in (Figura 12). Più precisamente a partire dalla misura della risposta indiciale si identificano, attraverso l’utilizzo di metodi numerici iterativi basati sulla minimizzazione dello scarto quadratico medio, i parametri di opportuni modelli semplificati, in modo tale da approssimare la risposta indiciale del processo stesso; tale procedura di fatto rappresenta un metodo per descrivere il comportamento dinamico del processo, per piccole variazioni dell’ingresso, attraverso un numero sufficientemente ridotto di parametri. Ebbene un primo modello che può es- 32 3.1 caratterizzazione input-output dello scambiatore di calore ∆θ lmu(t) [ ◦ C ] 4 2 0 −2 −4 10 15 20 25 30 35 40 30 35 40 time [s] 0.5 v [m s −1] 0.4 0.3 0.2 0.1 0 −0.1 10 15 20 25 time [s] (a) Andamento della velocità del fluido e della relativa variazione di temperatura del liquido in uscita a partire da condizioni nominali T H D = 16.18% 3.5 3 2.5 Mag 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Harmonic order (b) Analisi armonica dell’andamento della variazione di temperatura del liquido in uscita Figura 11.: Sollecitazione dell’impianto con un ingresso di velocità v(t) = v̄ + 0.2 sin(0.6 t) [m s−1 ] 1.4 1.2 ∆θ lmu(t ) [ ◦ C ] 1 0.8 0.6 0.4 0.2 0 0 10 20 30 40 50 60 70 80 Ti me [ s] m ( t ) a fronte della variazioFigura 12.: Andamento della variazione di temperatura ∆ ϑlu ne della velocità del fluido rispetto al valore nominale: ∆ v(t) = −0.0203 δ−1 (t − 20) [m s−1 ] 33 3.1 caratterizzazione input-output dello scambiatore di calore sere considerato è conosciuto in letteratura come First Order Plus Dead Time (F.O.P.D.T.) ed è descritto dalla funzione di trasferimento seguente: P(s) = Kp e−Td s 1+τs (F.O.P.D.T.) (3.1) La struttura del modello (3.1) è motivata da considerazioni fisiche, in particolare il modello (F.O.P.D.T. cerca di approssimare la risposta indiciale del processo caratterizzando il comportamento statico del medesimo attraverso il guadagno statico K p e caratterizzando invece il comportamento dinamico, tenendo conto della “inerzia” del processo modellata attraverso la costante di tempo τ ed il ritardo Td . Pertanto avendo chiarito quale sia la giustificazione della scelta del modello (F.O.P.D.T.) ed avendo messo in luce quale sia la procedura di determinazione dei parametri dei medesimi, al fine di valutare la bontà del modello approssimante che è stato determinato, si riportano i parametri del modello di tipo (F.O.P.D.T.) identificati dalla risposta (Figura 12), ed i relativi andamenti delle risposte indiciali confrontati con quella dell’impianto (Figura 13). Modello (F.O.P.D.T.) • K p = −66.385 [◦ C m−1 s] • τ = 3.4801 [s] • Td = 2.0473 [s] P(s) = −66.385 −2.0473 s e 1 + 3.4801 s (3.2) Osservando la (Figura 13) si nota che il modello approssimato (3.2) commette un errore nella rappresentazione della risposta indiciale dell’impianto abbastanza evidente. Al fine di ridurre tale errore è possibile fare ricorso ad un ulteriore modello approssimato conosciuto in letteratura come Second Order Real Poles Plus Dead Time (S.O.R.P.P.D.T.) e descritto della seguente funzione di trasferimento: P(s) = Kp e−Td s (1 + τ1 s) (1 + τ2 s) (S.O.R.P.P.D.T.) (3.3) 34 3.1 caratterizzazione input-output dello scambiatore di calore v → θ lmu [ ◦ C ] 1.4 1.2 ∆θ lmu(t ) [ ◦ C ] 1 0.8 0.6 0.4 0.2 F. O. P. D. T. R eal P l ant 0 20 30 40 50 60 70 80 Ti me [ s] Figura 13.: Andamento reale ed approssimato con il modello (F.O.P.D.T.) della variam ( t ) rispetto al valore nominale a fronte della zione di temperatura ∆ ϑlu variazione della velocità del fluido rispetto al valore nominale ∆ v(t) = −0.0203 δ−1 (t − 20) [m s−1 ] Si noti che il modello (S.O.R.P.P.D.T.) di fatto prevede l’aggiunta di una ulteriore costante di tempo al modello (F.O.P.D.T.. In tal modo si cerca di ottenere una maggiore aderenza tra la risposta dell’impianto e quella del modello approssimante. Si riportano di seguito i parametri identificati del modello di tipo (S.OR.P.P.D.T.) identificati dalla risposta (Figura 12), ed i relativi andamenti delle risposte indiciali confrontati con quella dell’impianto (Figura 14). Modello (S.O.R.P.P.D.T.) I parametri del modello (S.O.R.P.P.D.T.) identificati a partire dalla risposta indiciale dell’impianto in (Figura 12) sono i seguenti: • K p = −67.38 [◦ C m−1 s] • τ1 = 2.2427 [s] • τ2 = 2.2427 [s] • Td = 0.90778 [s] tali parametri conducono alla funzione di trasferimento seguente per il modello in esame: P(s) = −67.38 e−0.90778 s (1 + 2.2427 s) (1 + 2.2427 s) (3.4) In conclusione, al fine di confrontare i risultati ottenuti nella loro totalità, si riportano in (Figura 15) le risposte indiciali dell’impianto reale, del modello (F.O.P.D.T. e del modello (S.O.R.P.P.D.T.) Osservando i risultati mostrati in (Figura 15) si nota che la migliore approssimazione dell’impianto è ottenuta adottando il modello (S.O.R.P.P.D.T.). 35 3.1 caratterizzazione input-output dello scambiatore di calore v → θ lmu [ ◦ C ] 1.4 1.2 ∆θ m l u(t ) [ ◦ C ] 1 0.8 0.6 0.4 0.2 S . O. R . P. P. D. T. R eal P l ant 0 20 30 40 50 60 70 80 Ti me [ s] Figura 14.: Andamento reale ed approssimato con il modello (S.O.R.P.P.D.T.) delm ( t ) rispetto al valore nominale a fronla variazione di temperatura ∆ ϑlu te della variazione della velocità del fluido rispetto al valore nominale ∆ v(t) = −0.0203 δ−1 (t − 20) [m s−1 ] v → θ lmu [ ◦ C ] 1.4 1.2 ∆θ m l u(t ) [ ◦ C ] 1 0.8 0.6 0.4 R eal P l ant F. O. P. D. T. 0.2 S . O. R . P. P. D. T. 0 20 30 40 50 60 70 80 Ti me [ s] Figura 15.: Andamento reale ed approssimato con i modelli (S.O.R.P.P.D.T.) e m ( t ) rispetto al valore no(F.O.P.D.T. della variazione di temperatura ∆ ϑlu minale a fronte della variazione della velocità del fluido rispetto al valore nominale ∆ v(t) = −0.0203 δ−1 (t − 20) [m s−1 ] m 3.1.2 Analisi del legame ϑli − ϑlu In questa sezione si vogliono effettuare delle analisi sul legame tra la temperatura misurata del fluido in uscita e la temperatura del fluido in ingresso, tale legame, fissata la velocità del fluido in ingresso v, risulta di tipo lineare come 36 3.1 caratterizzazione input-output dello scambiatore di calore m ◦ θ l i → θ lu [ C] 2 ∆θ lmu(t ) [ ◦ C ] 1.5 1 0.5 0 20 30 40 50 60 70 80 Ti me [ s] m ( t ) a fronte della variaFigura 16.: Andamento della variazione di temperatura ∆ ϑlu zione della temperatura del fluido in ingresso rispetto al valore nominale ∆ θli = 5 δ−1 (t − 20) [◦ C ] si può osservare dalla ( 2.24 a pagina 27), pertanto nel caso in esame non verrà effettuata nessuna analisi di distorsione armonica, come mostrato nel caso di cui sopra, viceversa la procedura di identificazione della risposta indiciale illustrata nel paragrafo precedente verrà condotta nuovamente in questa sede, considerando la risposta indiciale dell’impianto rispetto alla temperatura del liquido in ingresso θli , a partire dalle condizioni nominali, mostrata in (Figura 16). Adottando la medesima procedura di identificazione illustrata nella sezione precedente si determinano i parametri del modello in questione riportati di seguito: Modello (F.O.P.D.T.) • K p = 0.40474 • τ = 1.6055 [s] • Td = 7.1951 [s] I quali conducono alla seguente funzione di trasferimento: P(s) = 0.40474 e−7.1951 s (1 + 1.6055 s) (3.5) mentre in (Figura 17) si riporta il relativo andamento della risposta indiciale confrontato con quella dell’impianto. Osservando la (Figura 17) si può evincere che l’approssimazione fornita dal modello (F.O.P.D.T.) è accettabile. 37 3.1 caratterizzazione input-output dello scambiatore di calore m ◦ θ l i → θ lu [ C] 2 1.5 1 0.5 F. O. P. D. T. R eal P l ant 0 20 30 40 50 60 70 80 Ti me [ s] Figura 17.: Andamento reale ed approssimato con il modello (F.O.P.D.T.) della variam ( t ) rispetto al valore nominale a fronte della vazione di temperatura ∆ ϑlu riazione della temperatura del fluido in ingresso rispetto al valore nominale ∆ θli (t) = 5 δ−1 (t − 20) [◦ C ] m 3.1.3 Analisi del legame ϑv − ϑlu Per concludere la caratterizzazione del comportamento Input-Output del processo sotto analisi si effettua una analisi sul legame tra la temperatura del vapore e la temperatura misurata del fluido in uscita. Le considerazioni sopra esposte relativamente alla distorsione armonica restano valide anche in questo m . Ebbene sollecitando anche in caso in virtù della linearità del legame ϑv − ϑlu questo caso l’impianto a partire dalla condizione operativa nominale, con una variazione a gradino della temperatura del vapore, si ottiene per la risposta in uscita quanto mostrato in (Figura 18). Ebbene anche in questo caso, al fine di descrivere il legame in questione attraverso un numero sufficientemente ridotto di parametri, si ripete la procedura di identificazione illustrata in precedenza. In particolare considerando un modello di tipo (F.O.P.D.T.) si ottengono i seguenti valori per i parametri identificati: Modello (F.O.P.D.T.) • K p = 0.59525 • τ = 3.2133 [s] • Td = 1.68 [s] I quali conducono alla seguente funzione di trasferimento: P(s) = 0.59525 e−1.68 s (1 + 3.2133 s) (3.6) viceversa in (Figura 19) si riporta il relativo andamento della risposta indiciale confrontato con quella dell’impianto. 38 3.1 caratterizzazione input-output dello scambiatore di calore θ v → θ lmu [ ◦ C ] 3 ∆θ lmu(t ) [ ◦ C ] 2.5 2 1.5 1 0.5 0 20 30 40 50 60 70 80 Ti me [ s] m ( t ) a fronte della variaFigura 18.: Andamento della variazione di temperatura ∆ ϑlu zione della temperatura del fluido in ingresso rispetto al valore nominale ∆ θv = 5 δ−1 (t − 20) [◦ C ] θ v → θ lmu [ ◦ C ] 3 ∆ θ lmu(t ) [ ◦ C ] 2.5 2 1.5 1 0.5 F. O. P. D. T. R eal P l ant 0 20 30 40 50 60 70 80 Ti me [ s] Figura 19.: Andamento reale ed approssimato con il modello (F.O.P.D.T.) della variam ( t ) rispetto al valore nominale a fronte della vazione di temperatura ∆ ϑlu riazione della temperatura del fluido in ingresso rispetto al valore nominale ∆ θv (t) = 5 δ−1 (t − 20) [◦ C ] Anche in questo caso essendo l’approssimazione ottenuta tramite il modello (F.O.P.D.T.) poco soddisfacente si ricorre ad un modello di tipo (S.O.R.P.P.D.T.) per il quale la procedura di identificazione, basata sulla risposta indiciale in (Figura 18), fornisce i seguenti valori dei parametri in questione: 39 3.1 caratterizzazione input-output dello scambiatore di calore Modello (S.O.R.P.P.D.T.) • K p = 0.59525 • τ1 = 2.0403 [s] • τ2 = 2.0403 [s] • Td = 0.69411 [s] I quali conducono alla seguente funzione di trasferimento: P(s) = 0.59525 e−0.69 s (1 + 2.0403 s) (1 + 2.0403 s) (3.7) mentre la relativa risposta indiciale, confrontata con quella dell’impianto reale, è mostrata in (Figura 20) θ v → θ lmu [ ◦ C ] 3 ∆ θ lmu(t ) [ ◦ C ] 2.5 2 1.5 1 0.5 S . O. R . P. P. D. T. R eal P l ant 0 20 30 40 50 60 70 80 Ti me [ s] Figura 20.: Andamento reale ed approssimato con il modello (S.O.R.P.D.T.) della vam ( t ) rispetto al valore nominale a fronte delriazione di temperatura ∆ ϑlu la variazione della temperatura del vapore rispetto al valore nominale ∆ θv (t) = 5 δ−1 (t − 20) [◦ C ] 40 3.1 caratterizzazione input-output dello scambiatore di calore Anche in questo caso al fine di avere un quadro generale dei risultati ottenuti si mostrano in (Figura 21) le risposte indiciali dei modelli (F.O.P.D.T.), (S.O.R.P.D.T.) e dell’impianto reale. θ v → θ lmu [ ◦ C ] 3 ∆ θ lmu(t ) [ ◦ C ] 2.5 2 1.5 1 R eal P l ant F. O. P. D. T. 0.5 S . O. R . P. P. D. T. 0 20 30 40 50 60 70 80 Ti me [ s] Figura 21.: Andamento reale ed approssimato con il modello (F.O.P.D.T.) e con il mom ( t ) rispetto al dello (S.O.R.P.P.D.T.) della variazione di temperatura ∆ ϑlu valore nominale a fronte della variazione della temperatura del vapore rispetto al valore nominale ∆ θv (t) = 5 δ−1 (t − 20) [◦ C ] 41 3.2 3.2 analisi dell’impianto basata sul modello analisi dell’impianto basata sul modello Nella sezione precedente sono state condotte delle analisi, volte a mettere in risalto la componente lineare dell’impianto, basate sullo studio della risposta indiciale. Ebbene risulta evidente che tale tipo di analisi può essere condotta anche nei casi in cui non si disponga di un modello dell’impianto, tuttavia nel caso in esame disponendo di un modello del processo, ci si può giovare della conoscenza di tale modello al fine di condurre delle analisi del comportamento dinamico del medesimo in maniera più sistematica; ad esempio effettuando una linearizzazione dell’impianto attorno a diversi punti di equilibrio di interesse. A tal scopo si richiama la procedura generale di linearizzazione di un generico sistema non lineare forzato, attorno ad un generico equilibrio. Linearizzazione di un sistema non lineare Si consideri un generico sistema non lineare ẋ = f ( x, u) (3.8) y = h( x, u) e sia ( x̄ , ū) un qualunque punto di equilibrio per il sistema (3.8) ossia tale che: f ( x̄ , ū) = 0 si ponga allora: x̃ := x − x̄ ũ := u − ū ỹ := y − h( x̄, ū) pertanto assumendo f ( x , u) e h( x, u) sufficientemente regolari si ha: f ( x, u) = ∂ f ( x , u) ∂ f ( x , u) x̃ + ũ + h.o.t.( x, u) ∂x ∂u ( x ,u)=( x̄ ,ū) ( x ,u)=( x̄ ,ū) h( x, u) = h( x̄, ū) + ∂h( x , u) ∂h( x , u) x̃ + ũ + h.o.t.( x, u) ∂x ∂u ( x ,u)=( x̄ ,ū) ( x ,u)=( x̄ ,ū) Avendo indicato con h.o.t. i termini di ordine superiore al primo. Ebbene trascurando i termini di ordine superiore e osservando che x̃˙ = ẋ si ha: ∂ f ( x , u) ∂ f ( x , u) x̃˙ = x̃ + ũ ∂x ∂u ( x ,u)=( x̄ ,ū) ( x ,u)=( x̄ ,ū) ỹ = ∂h( x , u) ∂h( x , u) x̃ + ũ ∂x ∂u ( x ,u)=( x̄ ,ū) ( x ,u)=( x̄ ,ū) 42 3.2 analisi dell’impianto basata sul modello Il quale rappresenta l’approssimazione lineare del sistema di partenza, nell’intorno del punto ( x̄, ū). Per comodità di scrittura di può porre: ∂ f ( x , u) ∂x ( x ,u)=( x̄ ,ū) ∂ f ( x , u) B := ∂u ( x ,u)=( x̄ ,ū) ∂h( x , u) C := ∂x ( x ,u)=( x̄ ,ū) ∂h( x , u) D := ∂u ( x ,u)=( x̄ ,ū) A := (3.9) (3.10) (3.11) (3.12) (3.13) Potendo riscrivere in definitiva l’approssimazione lineare del sistema non lineare originario come segue: x̃˙ = A x̃ + B ũ ỹ = C x̃ + D ũ La procedura appena mostrata può essere condotta nel caso in esame prendendo in considerazione il modello dell’impianto (2.24), ossia nella forma: ẋ = A(v) x + B(v) w (3.14) y = Cx (3.15) Ebbene al fine ottenere una maggiore aderenza con la notazione indicata nella procedura di linearizzazione di cui sopra, si ponga: u := [v , w T ] T (3.16) in tal modo, indicando con: v̄, w̄, x̄ rispettivamente i valori della velocità del fluido, del vettore dei disturbi e dello stato, assunti nel generico punto di lavoro, si ottiene per l’impianto, un modello linearizzato, attorno al generico punto di lavoro, nella forma: x̃˙ = A x̃ + B ũ (3.17) ỹ = C x̃ (3.18) ove le matrici sopra indicate sono definite come segue: A := A(v̄) " ∂ A ∂ B x̄ + v̄ B := ∂v x̄, v̄, w̄ ∂v x̄, v̄, w̄ (3.19) # B(v̄) (3.20) Pertanto resta solo da chiarire quale o quali siano i punti di lavoro scelti per condurre l’analisi basata sulla linearizzazione dell’impianto. A tal scopo al 43 3.2 analisi dell’impianto basata sul modello fine di effettuare una analisi, il più possibile generale, del comportamento dell’impianto vengono presi in considerazione cinque punti di lavoro, compreso il nominale, scelti casualmente nella regione ammissibile in (Figura 5 a pagina 23). In (Tabella 1) sono riportati i vari punti operativi presi in considerazione. PerTabella 1.: Valore della temperatura del liquido in ingresso e del vapore per i punti di lavoro scelti. Il punto di lavoro nominale viene indicato con il numero 0. Punto di lavoro # 0 1 2 3 4 ϑli ◦ ( C) 65 75 50 75 50 ϑv ◦ ( C) 120 140 140 105 105 ciò, come già indicato in precedenza, al fine di determinare i diversi stati di equilibrio dell’impianto, dapprima si fissa la temperatura di uscita del liquido al valore nominale ossia, 98 ◦ C, successivamente si impone al vettore dei segnali esogeni w uno dei valori riportati in (Tabella 1), di seguito si determina la velocità del fluido all’equilibrio v̄, per i diversi punti di lavoro, tramite la ( 2.14 a pagina 22) ed infine si calcola lo stato di equilibrio x̄ risolvendo la seguente equazione algebrica nell’incognita x̄ A(v̄) x̄ + B(v̄) w = 0 (3.21) La soluzione della quale porge la seguente espressione generale per il valore dello stato di equilibrio x̄ x̄ = −A(v̄)−1 B(v̄) w̄ (3.22) A vantaggio di quanto appena esposto, si riportano nei paragrafi seguenti le analisi condotte, attraverso la procedura di linearizzazione effettuata attorno ai punti di lavoro indicati in (Tabella 1). 3.2.1 Analisi del comportamento dinamico del processo attorno al punto di lavoro nominale Lo studio che ci si appresta a mostrare in questo paragrafo è relativo al comportamento dinamico del processo in un intorno del punto di lavoro nominale. Ebbene tale punto di lavoro come già detto più volte è caratterizzato dai seguenti valori per le variabili di processo: • ϑli = 98◦ C • ϑv = 120◦ C • v = 0.31m.s−1 44 3.2 analisi dell’impianto basata sul modello a partire da tali informazioni e seguendo la procedura di linearizzazione sopra indicata si ottiene il modello linearizzato dell’impianto attorno al punto di lavoro nominale. Ottenuto il modello linearizzato del processo, si può passare ad effettuare un’analisi in frequenza del medesimo. In virtù di ciò è conveniente mettere in evidenza l’espressione della risposta forzata del sistema, nel dominio della variabile complessa s, ebbene indicando con U (s) e Y (s) rispettivamente, la trasformata di Laplace dell’ingresso vettoriale u(t) definito nella (3.16) e dell’uscita considerata, si ottiene, per la risposta forzata l’espressione seguente: Y ( s ) = C ( s I − A ) −1 B U ( s ) (3.23) In prima battuta, essendo il sistema a più ingressi, sembrerebbe opportuno effettuare un’analisi del massimo valore singolare della matrice di trasferimento tra gli ingressi e l’uscita considerata; tuttavia avendo considerato la temperatura del liquido di ingresso e la temperatura del vapore, come due disturbi è conveniente analizzare i singoli legami separatamente, considerando l’azione di un ingresso alla volta. Per questo è necessario riscrivere la risposta forzata (3.23), mettendo in evidenza l’azione distinta dei tre ingressi, in tal modo si ottiene l’espressione seguente: Y (s) = C (s I − A)−1 (b1 v(s) + b2 w1 (s) b3 w2 (s)) (3.24) Ove w1 (s) , w2 (s) sono le due componenti del vettore dei disturbi definito in ( 2.23 a pagina 26), mentre le tre matrici degli ingressi sopra indicate sono definite come segue: " # ∂ A ∂ B b1 := x̄ + v̄ ∂v x̄, v̄, w̄ ∂v x̄, v̄, w̄ (3.25) b2 := B1 (v̄) b3 := B2 (v̄) Avendo indicato con B1 e B2 rispettivamente la prima e la seconda colonna della matrice B(v̄). La relazione (3.24), al fine di snellire la notazione, può essere riscritta come segue: Y (s) = Fy v (s) v(s) + Fy w1 (s) w1 (s) + Fy w2 (s) w2 (s) dove si sono indicate rispettivamente: Fy v (s) = C (s I − A)−1 b1 Fy w1 (s) = C (s I − A)−1 b2 Fy w2 (s) = C (s I − A) −1 b3 (3.26) (3.27) (3.28) Avendo distinto le tre azioni è possibile ottenere i diagrammi di Bode dei tre legami, i quali vengono mostrati di seguito in (Figura 22), (Figura 23) e (Figura 24). Dalla analisi appena condotta si può stabilire, osservando il diagramma mostrato in (Figura 23), che il motivo per il quale è stato osservato nell’analisi nel tempo dell’impianto, un notevole ritardo nella risposta indiciale nel legame m , è giustificato dall’elevato grado relativo della risposta in questione, ϑli − ϑlu tale caratteristica del legame in esame di fatto discende dalla presenza di un 45 3.2 analisi dell’impianto basata sul modello Bode Diagram 60 Magnitude (dB) 40 20 0 −20 −40 −60 Phase (deg) 180 135 90 45 0 −3 10 −2 −1 10 0 1 10 10 Frequency (rad/sec) 10 2 10 Figura 22.: Diagramma di Bode della funzione di trasferimento Fy v (s) Bode Diagram Magnitude (dB) 0 −200 −400 −600 −800 Phase (deg) −1000 0 −360 −720 −1080 −1440 −1800 −2160 −2520 −2880 −3240 −3600 −2 10 −1 10 0 10 Frequency (rad/sec) 1 10 2 10 Figura 23.: Diagramma di Bode della funzione di trasferimento Fy w1 (s) ritardo di trasporto nell’impianto reale, il quale si manifesta nell’approssimazione ad elementi finiti attraverso un alto grado relativo. Una ulteriore analisi che può essere condotta, disponendo delle funzioni di trasferimento dei tre legami in questione, è la valutazione della bontà delle 46 3.2 analisi dell’impianto basata sul modello Bode Diagram Magnitude (dB) 0 −50 −100 Phase (deg) 0 −45 −90 −135 −180 −225 −270 −2 10 −1 0 10 1 10 Frequency (rad/sec) 10 Figura 24.: Diagramma di Bode della funzione di trasferimento Fy w2 (s) approssimazioni, dei legami dell’impianto, fornite dai modelli (F.O.P.D.T.) e (S.O.R.P.D.T.) e mostrate in (Capitolo 3.1 a pagina 29). Pertanto in (Figura 25), (Figura 26) e (Figura 27) vengono mostrati dei confronti tra le tre funzioni di risposta armonica e le relative migliori approssimazioni ottenute nel (Capitolo 3.1 a pagina 29). Bode Diagram Re al P l ant Magnitude (dB) 40 S. O. R. P. P. D . T. 20 0 −20 −40 Phase (deg) 180 135 90 45 0 −2 10 −1 10 0 10 Frequency (rad/sec) 1 10 Figura 25.: Diagramma di Bode della funzione di trasferimento Fy v (s) e della sua approssimante (3.3) 47 3.2 analisi dell’impianto basata sul modello Bode Diagram 10 F. O. P. D . T Magnitude (dB) 0 Re al P l ant −10 −20 −30 −40 −50 360 Phase (deg) 0 −360 −720 −1080 −1440 −1800 −2160 −2 10 −1 0 10 1 10 10 Frequency (rad/sec) Figura 26.: Diagramma di Bode della funzione di trasferimento Fy w1 (s) e della sua approssimante (3.5) Bode Diagram 20 Magnitude (dB) Re al P l ant 0 S. O. R. P. P. D . T. −20 −40 −60 90 Phase (deg) 0 −90 −180 −270 −360 −450 −540 −2 10 −1 0 10 10 1 10 Frequency (rad/sec) Figura 27.: Diagramma di Bode della funzione di trasferimento Fy w2 (s) e della sua approssimante (3.7) Osservazione 3.1: Osservando le (Figura 25), (Figura 26) e (Figura 27) si nota che le approssimazioni ottenute tramite identificazione dei parametri della risposta indiciale sono soddisfacenti solamente a bassa frequenza ciò è deriva in parte della natura semplificata dei modelli 48 3.2 analisi dell’impianto basata sul modello approssimanti considerati. 3.2.2 Analisi del comportamento dinamico del processo attorno ad alcuni punti di lavoro ammissibili In questo paragrafo si mostreranno solamente i diagrammi di Bode delle tre funzioni di trasferimento Fy v (s), Fy w1 (s), Fy w2 (s) ottenute linearizzando il modello dello scambiatore di calore attorno ai punti di lavoro indicati in (Tabella 1). Punto di lavoro 1 : θli = 75 ◦C, θv = 140 ◦C In (Figura 28),(Figura 29) e (Figura 30) si mostrano rispettivamente i diagrammi di Bode delle funzioni di trasferimento Fy v (s), Fy w1 (s), Fy w2 (s) relativamente al punto di lavoro di cui sopra. Bode Diagram 40 Magnitude (dB) 20 0 −20 −40 −60 −80 Phase (deg) 180 135 90 45 0 −45 −2 10 −1 10 0 10 Frequency (rad/sec) 1 10 2 10 Figura 28.: Diagramma di Bode della funzione di trasferimento Fy v (s) per il punto di lavoro θli = 75 ◦C θv = 140 ◦C 49 3.2 analisi dell’impianto basata sul modello Bode Diagram 200 Magnitude (dB) 0 −200 −400 −600 −800 −1000 Phase (deg) 0 −720 −1440 −2160 −2880 −3600 −2 10 −1 10 0 10 Frequency (rad/sec) 1 10 2 10 Figura 29.: Diagramma di Bode della funzione di trasferimento Fy w1 (s) per il punto di lavoro θli = 75 ◦C θv = 140 ◦C Bode Diagram Magnitude (dB) 0 −50 −100 −150 Phase (deg) 0 −45 −90 −135 −180 −225 −270 −2 10 −1 10 0 10 Frequency (rad/sec) 1 10 2 10 Figura 30.: Diagramma di Bode della funzione di trasferimento Fy w2 (s) per il punto di lavoro θli = 75 ◦C θv = 140 ◦C 50 3.2 analisi dell’impianto basata sul modello Punto di lavoro 2 : θli = 50 ◦C θv = 140 ◦C In (Figura 31), (Figura 32) e (Figura 33) si mostrano rispettivamente i diagrammi di Bode delle funzioni di trasferimento Fy v (s), Fy w1 (s), Fy w2 (s) relativamente al punto di lavoro di cui sopra. Bode Diagram 60 Magnitude (dB) 40 20 0 −20 −40 −60 Phase (deg) 180 135 90 45 0 −2 10 −1 10 0 10 Frequency (rad/sec) 1 10 2 10 Figura 31.: Diagramma di Bode della funzione di trasferimento Fy v (s) per il punto di lavoro θli = 50 ◦C θv = 140 ◦C 51 3.2 analisi dell’impianto basata sul modello Bode Diagram Magnitude (dB) 0 −200 −400 −600 −800 −1000 Phase (deg) 0 −720 −1440 −2160 −2880 −3600 −2 10 −1 10 0 10 Frequency (rad/sec) 1 10 2 10 Figura 32.: Diagramma di Bode della funzione di trasferimento Fy w1 (s) per il punto di lavoro θli = 50 ◦C θv = 140 ◦C Bode Diagram Magnitude (dB) 0 −50 −100 −150 45 Phase (deg) 0 −45 −90 −135 −180 −225 −270 −2 10 −1 10 0 10 Frequency (rad/sec) 1 10 2 10 Figura 33.: Diagramma di Bode della funzione di trasferimento Fy w2 (s) per il punto di lavoro θli = 50 ◦C θv = 140 ◦C 52 3.2 analisi dell’impianto basata sul modello Punto di lavoro 3 : θli = 75 ◦C θv = 105 ◦C In (Figura 34), (Figura 35) e (Figura 36) si mostrano rispettivamente i diagrammi di Bode delle funzioni di trasferimento Fy v (s), Fy w1 (s), Fy w2 (s) relativamente al punto di lavoro di cui sopra. Bode Diagram Magnitude (dB) 40 20 0 −20 −40 −60 −80 Phase (deg) 180 135 90 45 0 −2 10 −1 10 0 10 Frequency (rad/sec) 1 10 2 10 Figura 34.: Diagramma di Bode della funzione di trasferimento Fy v (s) per il punto di lavoro θli = 75 ◦C θv = 105 ◦C 53 3.2 analisi dell’impianto basata sul modello Bode Diagram 200 Magnitude (dB) 0 −200 −400 −600 −800 −1000 Phase (deg) 0 −720 −1440 −2160 −2880 −3600 −2 10 −1 10 0 10 Frequency (rad/sec) 1 10 2 10 Figura 35.: Diagramma di Bode della funzione di trasferimento Fy w1 (s) per il punto di lavoro θli = 75 ◦C θv = 105 ◦C Bode Diagram Magnitude (dB) 0 −50 −100 −150 45 Phase (deg) 0 −45 −90 −135 −180 −225 −270 −2 10 −1 10 0 10 Frequency (rad/sec) 1 10 2 10 Figura 36.: Diagramma di Bode della funzione di trasferimento Fy w2 (s) per il punto di lavoro θli = 75 ◦C θv = 105 ◦C 54 3.2 analisi dell’impianto basata sul modello Punto di lavoro 4 : θli = 50 ◦C θv = 105 ◦C In (Figura 37), (Figura 38) e (Figura 39) si mostrano rispettivamente i diagrammi di Bode delle funzioni di trasferimento Fy v (s), Fy w1 (s), Fy w2 (s) relativamente al punto di lavoro di cui sopra. Bode Diagram Magnitude (dB) 50 0 −50 −100 Phase (deg) 180 135 90 45 0 −2 10 −1 10 0 10 Frequency (rad/sec) 1 10 2 10 Figura 37.: Diagramma di Bode della funzione di trasferimento Fy v (s) per il punto di lavoro θli = 50 ◦C θv = 105 ◦C Infine di riportano in (Figura 40), (Figura 41) e (Figura 42) i diagrammi di Bode determinati linearizzando l’impianto intorno a tutti i punti di lavoro considerati. 55 3.2 analisi dell’impianto basata sul modello Bode Diagram 200 Magnitude (dB) 0 −200 −400 −600 −800 −1000 Phase (deg) 0 −720 −1440 −2160 −2880 −3600 −2 10 −1 10 0 10 Frequency (rad/sec) 1 10 2 10 Figura 38.: Diagramma di Bode della funzione di trasferimento Fy w1 (s) per il punto di lavoro θli = 50 ◦C θv = 105 ◦C Bode Diagram Magnitude (dB) 0 −50 −100 −150 45 Phase (deg) 0 −45 −90 −135 −180 −225 −270 −2 10 −1 10 0 10 Frequency (rad/sec) 1 10 2 10 Figura 39.: Diagramma di Bode della funzione di trasferimento Fy w2 (s) per il punto di lavoro θli = 50 ◦C θv = 105 ◦C 56 3.2 analisi dell’impianto basata sul modello Bode Diagram F0 Magnitude (dB) 50 F1 F2 0 F3 F4 −50 Phase (deg) 180 135 90 45 0 −45 −2 10 −1 10 0 1 10 10 Frequency (rad/sec) 2 10 3 10 Figura 40.: Diagrammi di Bode della funzione di trasferimento Fy v (s) per i diversi punti di lavoro. Le diverse funzioni di trasferimento nel diagramma sono indicate con Fi (s) con i ∈ {0, 1, 2, 3, 4} a seconda del punto di lavoro considerato, in maniera consistente alla numerazione riportata in (Tabella 1) Bode Diagram 200 F 0 (s) Magnitude (dB) 0 F 1 (s) −200 F 2 (s) −400 F 3 (s) −600 F 4 (s) −800 −1000 Phase (deg) 0 −720 −1440 −2160 −2880 −3600 −2 10 −1 10 0 1 10 10 Frequency (rad/sec) 2 10 3 10 Figura 41.: Diagrammi di Bode della funzione di trasferimento Fy w1 (s) per i diversi punti di lavoro. Le diverse funzioni di trasferimento nel diagramma sono indicate con Fi (s) con i ∈ {0, 1, 2, 3, 4} a seconda del punto di lavoro considerato, in maniera consistente alla numerazione riportata in (Tabella 1) 57 3.2 analisi dell’impianto basata sul modello Bode Diagram 50 F 0 (s) Magnitude (dB) 0 F 1 (s) F 2 (s) −50 F 3 (s) −100 F 4 (s) −150 45 Phase (deg) 0 −45 −90 −135 −180 −225 −270 −2 10 −1 10 0 1 10 10 Frequency (rad/sec) 2 10 3 10 Figura 42.: Diagrammi di Bode della funzione di trasferimento Fy w2 (s) per i diversi punti di lavoro. Le diverse funzioni di trasferimento nel diagramma sono indicate con Fi (s) con i ∈ {0, 1, 2, 3, 4} a seconda del punto di lavoro considerato, in maniera consistente alla numerazione riportata in (Tabella 1) 58 Parte II. Controllo dello scambiatore di calore 59 4 I 4.1 C O N T R O L LO D E L LO S C A M B I ATO R E DI CALORE descrizione del problema di controllo l problema che si vuole affrontare in questa parte è il problema della regolazione asintotica della temperatura misurata del fluido in uscita m , anche a fronte di variazioni della condizione operativa dallo scambiatore ϑlu del processo, ossia a fronte di variazioni ammissibili delle temperature del liquido in ingresso ϑli e del vapore ϑv . In particolare il problema appena esposto, con riferimento al modello dello scambiatore di calore indicato in ( 2.24 a pagina 27), può essere formalizzato come segue: Problema 4.1 (Problema di regolazione asintotica della temperatura): Sia assegnato il modello del processo P nella forma: ẋ (t) = A(v(t)) x (t) + B(v(t)) w(t) y(t) = C x (t) x (t ) = x 0 (4.1) 0 dove: x0 ∈ R71 t 0 ∈ R+ x (t) ∈ R71 C ∈ R1×71 u(·) : R+ → R w(·) : R+ → R2 Si supponga inoltre che x (t) e d(t) non siano misurabili. Trovare se esiste un controllore u(·) tale che: 1. lim ||y(t) − y∗ || = 0 ∀ y∗ ∈ R, ∀ w∗ ∈ Θy∗ , ∀ x0 ∈ R71 t→+∞ 2. x (·) , u(·) ∈ L∞ 3. Il controllo u(·) è tale che se: x0 : y ( t0 ) = y ∗ 60 (4.2) allora deve risultare: 4.2 problema di regolazione asintotica y ( t ) = y ∗ ∀ t ≥ t0 Osservando il problema appena esposto si può osservare che la specifica (1) riassume la proprietà di regolazione asintotica in maniera più formale, la specifica (2) impone che sia lo stato che il controllo siano limitati ed infine la (3) impone che se l’uscita all’istante iniziale è uguale alla temperatura di riferimento y∗ allora essa deve rimanere tale identicamente. E’ evidente che la struttura del problema esposto è del tutto generale e prescinde dalla realtà fisica in cui il problema centrale della presente trattazione di colloca. Pertanto, al fine di mostrare come sia possibile risolvere il problema (4.1), si farà riferimento ad un problema del tutto generale. Successivamente, una volta chiarito quali siano le diverse strategie per la soluzione del problema di regolazione, si mostreranno i risultati ottenuti applicando tali metodologie al caso in esame. 4.2 problema di regolazione asintotica Il problema di controllo che si vuole affrontare si colloca nell’ambito del problema di regolazione [2]. Tale problema rappresenta un caso particolare del problema di inseguimento di traiettoria. In particolare questo problema può essere, in generale, formalizzato come segue: Problema 4.2 (Problema di regolazione perfetta): Sia P un sistema dinamico assegnato del tipo: ẋ (t) = f ( x (t), u(t), d(t)) y(t) = h( x (t)) x (t ) = x 0 (4.3) 0 dove: x 0 ∈ Rn t 0 ∈ R+ x ( t ) ∈ Rn y(t) ∈ R u(·) : R+ → R d(·) : R+ → R p Si supponga inoltre che: f (·) ∈ C 1 (Rn × R × R p ), h(·) ∈ C 1 (Rn ) x (t) e d(t) sono non noti. Trovare se esiste un controllo u(·) tale che: (4.4) (4.5) 61 4.2 problema di regolazione asintotica 1. lim ||y(t) − y∗ || = 0 ∀ y∗ ∈ R ∀ d∗ ∈ R p , ∀ x0 ∈ Rn t→+∞ (4.6) 2. x (·) , u(·) ∈ L∞ 3. Il controllo u(·) garantisce che se: x0 : h ( x0 ) = y ∗ allora: h( x (t)) = y∗ ∀t ≥ t0 Per comodità di esposizione è conveniente introdurre la seguente grandezza ỹ(t) = y ∗ −y(t) detta errore di regolazione in tal modo si può riscrivere la (1) nella seguente forma equivalente: lim ||ỹ(t)|| = 0 ∀ y∗ ∈ R, ∀ d∗ ∈ R p , ∀ x0 ∈ Rn t→+∞ Il problema esposto di fatto consiste di determinare, se esiste, un controllo u(t) limitato in grado di annullare asintoticamente l’errore di regolazione a partire da qualunque condizione iniziale dell’impianto con l’aggiunta che la traiettoria nello stato, e quindi sotto le ipotesi di regolarità della funzione h( x ) indicate in (4.4) dell’uscita, sia limitata. Infine si osservi che la specifica (3) garantisce che se l’errore di regolazione ỹ è nullo nell’istante iniziale allora esso si mantiene nullo identicamente, l’aggiunta di tale specifica conferisce l’attributo perfetto al problema di regolazione (Problema 4.2). Al fine di proseguire la trattazione in maniera più agevole è conveniente indicare con Σ il sistema ottenuto controllando l’impianto P attraverso il generico controllo u(·), a tale sistema ci si riferirà nel corso della trattazione con l’accezione Sistema di controllo Σ. 4.2.1 Primi passi verso la soluzione del problema di regolazione e assunzioni sul problema Il problema di regolazione appena esposto equivale ad un problema di stabilizzazione di un particolare equilibrio, infatti si osservi che le traiettorie di un qualunque sistema dinamico, purché limitate, convergono sempre su insiemi invarianti, perciò affinché l’errore di regolazione possa annullarsi asintoticamente, risulta necessario che, fissati y∗ e d∗ , il sistema di controllo Σ possieda un punto di equilibrio per il quale si abbia y ≡ y∗ . Ebbene in virtù della ipo- 62 4.2 problema di regolazione asintotica tesi di stazionarietà dell’impianto P e di quanto appena asserito, si osservi che affinché il problema di regolazione ammetta soluzione è necessario che: ∃ ( x̄ , ū) ∈ Rn+1 : f ( x̄ , ū, d∗ ) = 0 , h( x̄ ) = y∗ . Tuttavia nella presente trattazione per semplicità si assumerà non solo che la coppia ( x̄, ū), dalla sopracitata proprietà esista, ma anche che essa sia unica. Ossia si farà la seguente assunzione: Assunzione 4.1: ∀ (y∗ , d∗ ) ∈ R p+1 ∃! ( x̄ , ū) ∈ Rn+1 : f ( x̄ , ū, d∗ ) = 0 , h( x̄ ) = y∗ Si noti che il fatto di aver assunto l’unicità di tale coppia non è fortemente restrittivo, ai fini primari della seguente trattazione, in quanto tale ipotesi può essere ragionevolmente ritenuta soddisfatta nell’ambito del controllo dei processi.1 Per comprendere quale sia la struttura della legge di controllo in grado di risolvere il problema di regolazione, si osservi che il valore del controllo ū è non noto e dipendente da y∗ e da d∗ , per cui sotto le ipotesi di non misurabilità di d∗ è immediato comprendere che la legge di controllo u(·) sia necessariamente in retroazione dall’uscita dell’impianto P , ma non solo infatti dovendo il sistema di controllo ammettere un equilibrio z̄ per il quale y ≡ y∗ sembra ragionevole utilizzare un controllore in retroazione dall’errore di regolazione ỹ, infatti in tal modo l’equilibrio z̄ è indipendente dai singoli valori di y∗ e y. In definitiva dunque il controllore K candidato per la soluzione del problema di regolazione è espresso nella forma seguente: ( ż = ϕ(z, ỹ) K: (4.7) u = ε(z, ỹ) E’ evidente che il problema di regolazione appena esposto, data la natura non lineare del sistema considerato non è di semplice soluzione. Tuttavia rilassando alcune specifiche è possibile determinare una legge di controllo lineare in grado di assicurare una regolazione di tipo locale. In altre parole attraverso l’uso di un controllore LTI in retroazione da ỹ è possibile, sotto (l’Assunzione 4.1), garantire che: lim ||ỹ(t)|| = 0 ∀ y∗ ∈ R, ∀ d∗ ∈ R p , ∀ x0 ∈ I( x̄ ) ⊂ Rn t→+∞ . Quando detto induce a formulare un ulteriore problema di regolazione semplificato. In particolare tale problema verrà denominato problema di regolazione locale e verrà formulato come segue: 1 Si noti che nel caso dello scambiatore di calore, questa ipotesi è certamente verifica in quanto m ,ϑ , ϑ esiste in corrispondenza un solo valore della velocità del fluido v fissato un valore per ϑlu v li e un solo stato di equilibrio x 63 4.2 problema di regolazione asintotica Problema 4.3 (Problema di regolazione perfetta locale): Sia P un sistema dinamico assegnato del tipo: ẋ (t) = f ( x (t), u(t), d(t)) x ( t0 ) = x0 y(t) = h( x (t)) (4.8) dove: x ( t ) ∈ Rn y(t) ∈ R u(·) : R+ → R d(·) : R+ → R p Si supponga inoltre che: • f (·) ∈ C 1 (Rn × R × R p ) h(·) ∈ C 1 (Rn ) • • x (t) e d(t) sono non noti. Trovare se esiste un controllo u(·) tale che: 1. lim ||ỹ(t)|| = 0 ∀ y∗ ∈ R, ∀d∗ ∈ R p , ∀ x0 ∈ I( x̄ (y∗ , d∗ )) ⊂ Rn (4.9) t→+∞ 2. x (·) , u(·) ∈ L∞ 3. Il controllo u(·) garantisce che se: x0 : h ( x0 ) = y ∗ allora: h( x (t)) = y∗ ∀t ≥ t0 Fatte le relative precisazioni si può passare a mostrare quale sia la soluzione del problema 4.3. In particolare la soluzione che si vuole adottare consiste nell’utilizzo di un controllore LTI in retroazione dall’errore di regolazione ỹ del tipo: ( ϕ̇ = Ψ ϕ + Λ ỹ K: (4.10) u = K1 ϕ + K2 ỹ 64 4.2 problema di regolazione asintotica con: ϕ ∈ R` K1 ∈ R K2 ∈ R (4.11) 1×` Ψ ∈ R (4.12) (4.13) `×` (4.14) Λ ∈ R` (4.15) L’adozione della legge di controllo (4.11) per il controllo dell’impianto P da luogo allo schema di controllo mostrato in (Figura 43). Ebbene al fine di stabili- Figura 43.: Schema di controllo adottato re sotto quali condizioni esista un controllore K che sia soluzione del problema di regolazione locale (4.3), è conveniente enunciare alcuni risultati. 65 4.2 problema di regolazione asintotica Lemma 4.1: In merito al sistema in (4.3), si fissino y∗ ∈ R, d∗ ∈ R p in modo qualunque ma a patto che valga la Assunzione 4.1 con ū 6= 0. Si osservi ora che il controllore (4.10) fornisce la seguente dinamica per il sistema di controllo: ẋ = f ( x, ϕ, d∗ , y∗ ) (4.16) ϕ̇ = Ψ ϕ + Λ (y∗ − h( x )) Ebbene se sotto l’azione di tale controllo esiste, per il sistema di controllo (A.11), uno stato di equilibrio ( x̂ , ϕ̂) ∈ R`+n tale che: h( x̂ ) = y∗ allora valgono le seguenti affermazioni: • ū = K1 ϕ̂. • la matrice Ψ ammette λ = 0 come autovalore ed inoltre ϕ̂ ∈ Ker (Ψ) • K1T ∈ / Im Ψ Dimostrazione. La prova del lemma è riportata nell’Appendice A Il lemma appena enunciato di fatto afferma che affinché il K controllore sia soluzione del problema di regolazione (4.3) esso deve necessariamente comprendere una azione integrale. Teorema 4.1: In merito al sistema in (4.3), si fissino y∗ ∈ R, d∗ ∈ R p arbitrariamente a patto che valga la Assunzione 4.1. Ebbene in si indichino ora rispettivamente: ∂ f ( x , u, d∗ ) A := ∂x ( x ,u)=( x̄ ,ū) ∂ f ( x , u , d∗ ) B := ∂u ( x ,u)=( x̄ ,ū) ∂h( x ) C := ∂x (x ,u)=( x̄ ,ū) (4.17) dunque se: • la terna ( A, B, C ) è stabilizzabile e rilevabile A B • rank = n+1 C 0 allora esiste un controllore K nella forma (4.10) che sia soluzione del problema (4.3) ed inoltre esso è espresso nella forma seguente: ϕ̇1 = λ1 ỹ K : ϕ̇2 = Ψ̂ ϕ2 + λ2 ỹ u = k11 ϕ1 + K12 ϕ2 + K2 ỹ (4.18) 66 4.2 problema di regolazione asintotica dove: ϕ1 ∈ R (4.19) λ1 ∈ R (4.21) k11 ∈ R (4.23) ϕ2 ∈ R`−1 λ2 ∈ R`−1 K12 ∈ R1×`−1 Dimostrazione. La prova del teorema è riportata nell’Appendice A (4.20) (4.22) (4.24) (4.25) Osservazione 4.1: La trattazione sviluppata fino ad ora di per se garantisce che fissati d∗ e y∗ è possibile progettare un regolatore in grado di annullare asintoticamente l’errore di regolazione a patto che lo stato iniziale dell’impianto sia sufficientemente prossimo allo stato di equilibrio che garantisce l’errore nullo di regolazione. Ciononostante nella realtà sia il riferimento y∗ sia il vettore dei disturbi d∗ non sono noti, ciò implica che di fatto il modello linearizzato dell’impianto non risulta univocamente determinato! Si noti che il problema che si presenta al variare delle condizioni operative è legato solamente alla stabilizzazione del processo allargato. Infatti la presenza dell’azione integrale conferisce al sistema di controllo di per se un elevata robustezza in quanto tale azione rende il punto di lavoro desiderato un equilibrio per il sistema di controllo, qualunque siano d∗ e y∗ , per cui se si riesce determinare una legge di controllo in grado di rendere il generico equilibrio asintoticamente stabile, il problema è risolto indipendentemente dal punto di lavoro considerato. E’ evidente però che utilizzando tecniche di stabilizzazione lineare la dipendenza dal punto di lavoro è un fattore critico in quanto tali tecniche sono basate sul modello linearizzato dell’impianto il quale dipende fortemente dal punto di lavoro considerato. Fortunatamente, sopratutto nell’ambito del controllo dei processi, è noto almeno un punto di lavoro nominale dell’impianto, perciò in virtù di tale conoscenza è possibile determinare un modello nominale sul quale basare la procedura di sintesi del regolatore. E’ evidente che il modello di cui si dispone in questo modo fornisce una conoscenza parziale dell’impianto per cui il progetto delle sintesi di controllo va effettuato garantendo che il controllore stabilizzante sia in grado di stabilizzare una famiglia di impianti e non il singolo impianto nominale. Pertanto attraverso l’introduzione dei due risultati di cui sopra, e partizionando il controllore K come segue: ( ϕ̇1 = λ1 ỹ I: u1 = k11 ϕ1 ( (4.26) K: ϕ̇2 = Ψ̂ ϕ2 + λ2 ỹ K : s u2 = K12 ϕ2 + K2 ỹ u = u +u 2 1 si perviene ad una legge di controllo che può essere rappresentata attraverso lo schema in (Figura 44) 67 4.2 problema di regolazione asintotica Figura 44.: Schema di controllo con azione integrale 4.2.2 Regolazione con azione feed-back+feed-forward Il problema di regolazione è stato risolto, almeno nella sua versione locale, con successo. Tuttavia si noti che la soluzione ottenuta, utilizzando il compensatore K in (4.10), potrebbe causare dei problemi in caso di variazioni a gradino del riferimento y∗, infatti in tali situazioni l’andamento della variabile di controllo u(t), a causa del legame diretto ingresso uscita presente nel controllore K , assumerebbe un andamento discontinuo provocando delle eccessive sollecitazioni dell’attuatore. Per risolvere tale problema è possibile effettuare una sistema di controllo comprensivo anche di un azione in avanti, ossia invece di utilizzare come riferimento il segnale y∗ si utilizza come riferimento l’uscita di un sistema LTI avente come ingresso il riferimento y∗ (azione in avanti), lo schema di controllo che deriva da quanto detto è mostrato in (Figura 45). Figura 45.: Schema di controllo con azione avanti Per comprendere se in tali condizioni il controllore K in (A) è ancora soluzione del problema di regolazione proposto, è conveniente enunciare il seguente risultato: 68 4.2 problema di regolazione asintotica Proposizione 4.1: Sia H il seguente sistema espresso nella forma: ε̇ = H ε + η y∗ m = Lε (4.27) e sia ΣH il sistema di controllo derivante dall’utilizzo del sistema H come azione in avanti per il sistema di controllo Σ, ossia: ẋ = f ( x, u , d∗ ) ε̇ = H ε + η y∗ ϕ̇ = Ψ ϕ + Λ ( L ε − h( x )) (4.28) u = K1 ϕ + K2 ( L ε − h( x )) e sia il controllore K (4.10) soluzione del problema di regolazione allora esso risolve anche il problema di regolazione con azione in avanti se e soltanto se: • H ≺ 0 • L H −1 η = − 1 Dimostrazione. Se il controllore (4.10) risolve il problema di regolazione allora ∃( x̄ , ϕ̄) ∈ Rn+` tale che: f ( x̄ , ϕ̄) = 0 (4.29) Λ ϕ̄ = 0 (4.30) y = y∗ (4.31) Tuttavia per l’ipotesi poste sul sistema H nell’enunciato si ha che ∃ ε̄ ∈ R tale che: ε̄ = y∗ Si osservi quindi che il punto di equilibrio stato ( x̄, ε̄, ϕ̄) è tale che: f ( x, z , d∗ ) = 0 Ψ ϕ̄ = 0 (4.32) Ossia l’azione in avanti non altera le posizione del punto di equilibrio necessario alla regolazione. Per concludere l’assunto occorre e basta osservare che linearizzando il sistema si controllo ΣH attorno al punto ( x̄, ε̄, ϕ̄) e utilizzando la notazione introdotta nel Teorema A.1 si ottiene il seguente sistema LTI: x̃˙ A − B K2 C B K1 B K2 L x̃ ϕ̃˙ = −Λ C (4.33) Ψ Λ L ϕ̃ ˙ε̃ 0 0 H ε̃ Ebbene è evidente che tale sistema è asintoticamente stabile in quanto gli autovalori della sua matrice dinamica sono l’unione degli autovalori del sistema di controllo Σ, i quali sono certamente a parte reale negativa dato che K è soluzione del problema di regolazione per ipotesi, e degli autovalori del sistema H i quali sono a parte reale negativa per ipotesi. Resta ora da mostrare la necessità delle condizioni sopra enunciate. A tal scopo si osservi che se H fosse singolare allora non esisterebbe nessun punto di equilibrio per il sistema di controllo 69 4.3 metodi basati sulla risposta armonica per la regolazione ΣH , viceversa se fosse H 0 allora il sistema di controllo ΣH ammetterebbe il punto di equilibrio desiderato ma esso in questo caso sarebbe instabile. Per quanto riguarda la necessità della seconda delle condizioni, si noti che se essa non fosse vera allora, come conseguenza del Lemma A.1, il sistema di controllo ammetterebbe come unico punto di equilibrio un punto per i quale si avrebbe: y = − L H −1 η y ∗ 6 = y ∗ Concludendo in tal modo la dimostrazione dell’assunto. 4.2.3 Controllo Observer-Based Per quanto è stato mostrato nelle sezione precedente, al fine di risolvere il problema di regolazione è sufficiente progettare un compensatore in grado di stabilizzare il processo allargato. Detto ciò osservando la struttura della legge di controllo stabilizzante, è immediato constatare che la strategia di controllo che può essere utilizzata allo scopo di stabilizzare il processo allargato può essere ad esempio di tipo observer-based. Infatti, come è noto, dal principio di separazione, si può progettare un compensatore stabilizzante in due passi. In particolare indicando con η lo stato del processo allargato si può sintetizzare una legge di controllo stabilizzante nel modo seguente: 1. si determina un controllo stabilizzante in retroazione dallo stato η 2. si sintetizza un osservatore asintotico dello stato del processo allargato e si utilizza la stima dello stato η̂ fornita dall’osservatore per effettuare la retroazione dallo stato determinata al passo precedente. La procedura observer-based indicata tuttavia richiede la conoscenza dei parametri dell’impianto da controllare. Come già detto una conoscenza parziale dell’impianto può essere ottenuta considerando un opportuno modello nominale. In questo modo è evidente che il progetto del compensatore deve essere effettuato utilizzando delle metodologie che garantiscono una certa robustezza. Ebbene una tecnica di controllo observer-based in grado di garantire una buona robustezza è il ben noto controllo LQG. Chiaramente l’utilizzo di una metodologia di tipo observer-based essendo basata sul modello del processo è intrinsecamente poco robusta. Per fare fronte a ciò nella sezione seguente si mostrerà un ulteriore metodologia di sintesi per il compensatore, attraverso la quale si può garantire un maggiore grado di robustezza. 4.3 metodi basati sulla risposta armonica per la regolazione In questo capitolo si vuole illustrare come sia possibile ricorrere a tecniche di sintesi basate sulla risposta armonica per sintetizzare una legge di controllo L.T.I. in grado di risolvere il problema delle regolazione locale di un sistema non lineare. Nello studio dei sistemi di controllo lineari, il criterio di Nyquist rappresenta uno strumento di notevole importanza soprattutto in fase di sintesi di un 70 4.3 metodi basati sulla risposta armonica per la regolazione controllore, in quanto attraverso l’utilizzo di tale criterio è possibile progettare un controllore che assicuri la stabilità asintotica del sistema di controllo facendo riferimento alla risposta armonica del sistema di controllo ad anello aperto, sfortunatamente,nel controllo dei sistemi non lineari non è possibile fare ricorso, di per se al criterio di Nyquist. Tuttavia nell’ambito della presente trattazione, essendo il progetto del regolatore fondato su tecniche lineari basate quest’ultime, sul teorema di approssimazione lineare di Lyapunov, sembra ragionevole pensare che si possa fare anche in questo, caso uso del criterio di Nyquist per garantire le proprietà di stabilità asintotica di un dato stato di equilibrio. Pertanto si osservi che sotto le ipotesi di utilizzo della legge di controllo (controllore lineare), come già detto, al fine di studiare le proprietà di stabilità esponenziale locale del punto di equilibrio ( x̄, ϕ¯1 ) per il sistema di controllo, è sufficiente studiare la stabilità asintotica dell’approssimazione lineare di esso (controllo) ottenuta linearizzando quest’ultimo attorno al punto ( x̄, ϕ¯1 ). Tale approssimazione lineare di fatto può essere pensata come il sistema LTI che si otterrebbe interconnettendo opportunamente il controllore (4.10), di per se già lineare, e l’approssimazione lineare del processo. Attraverso l’ausilio di questa osservazione è possibile utilizzare tutti gli strumenti noti della teoria lineare del controllo, basati sullo studio della risposta armonica considerando le funzioni di trasferimento del controllore, la quale risulta ben definita essendo quest’ultimo lineare, e la funzione di trasferimento dell’approssimazione lineare dell’impianto. Inoltre per quanto è stato detto in precedenza essendo a priori non noti y∗ e d∗ non è nota univocamente la funzione di trasferimento del processo allargato. 4.3.1 Controllo PID based Una delle leggi di controllo utilizzate soprattutto nell’ambito del controllo dei processi è certamente la legge di controllo P.I.D.. Ebbene per quanto esposto nei risultati presentati fin ora, sotto opportune condizioni sull’impianto da controllare, è sempre possibile determinare una legge di controllo L.T.I. in grado di risolvere il problema delle regolazione, nella forma: ϕ̇1 = λ1 ỹ ϕ̇2 = Ψ̂ ϕ2 + λ2 ỹ (4.34) u = k11 ϕ1 + K12 ϕ2 + K2 ỹ Pertanto per mostrare che un regolatore P.I.D. sia in grado di risolvere il problema di regolazione è sufficiente mostrare che tale controllore è un caso particolare di quella sopra indicata. A riguardo si riporta l’espressione della funzione di trasferimento di un regolatore P.I.D. causale: 1 s + Kd s 1+τs K p , Ki , K d ∈ R τ ∈ R+ G ( s ) = K p + Ki (4.35) Per mostrare l’equivalenza di tale legge di controllo con quella più generale, è conveniente determinare una rappresentazione con lo spazio di stato di tale fun- 71 4.3 metodi basati sulla risposta armonica per la regolazione zione di trasferimento. A tal fine, tramite opportune manipolazioni algebriche, si riscrive la (4.35) come segue: (Ki − τ −1 Kd0 ) s + Ki τ −1 s 2 + τ −1 s ove è stato posto per comodità G (s) = K p + Kd0 + (4.36) Kd0 = Kd τ −1 A questo punto, indicando rispettivamente con ξ,v, w, lo stato, l’ingresso e l’uscita del controllore P.I.D., si può determinare una realizzazione minima della funzione di trasferimento (4.36), ad esempio in forma canonica raggiungibile, ottenendo: 0 1 0 ξ+ v ξ̇ = − 1 0 −τ 1 w = Ki τ −1 Ki − τ −1 Kd0 ξ + (K p + Kd0 ) v (4.37) ξ ∈ R2 , v ∈ R, w ∈ R A questo punto di noti che la matrice dinamica del controllore (4.37) è banalmente diagonalizzabile, per cui per porre il controllore (4.37) nella forma desiderata è sufficiente esprimere le equazioni (4.37) nella base degli autovettori della matrice in questione. Tale trasformazione, come è noto, può essere effettuata ponendo: ϕ = Tξ (4.38) dove la matrice T è quadrata di dimensione 2, non singolare e tale che la sua inversa T −1 possieda come colonne una base di autovettori della matrice dinamica che si vuole rendere diagonale. In particolare nel caso in questione si ha: T −1 1 1 = − τ −1 0 (4.39) Perciò effettuando la trasformazione suddetta si ottiene: − τ −1 0 −τ ϕ+ v 0 0 τ h 0 i w = K2d Kτi ϕ + (K p + Kd0 ) v τ ϕ̇ = (4.40) ϕ ∈ R2 , v ∈ R, w ∈ R Infine esplicitando le equazioni differenziali che descrivono l’evoluzione delle due componenti del vettore ϕ si ottiene: ϕ̇1 = −τ −1 ϕ1 − τ v ϕ̇2 = τ v K0 Ki w = 2d ϕ1 + ϕ2 + (K p + Kd0 ) v τ τ (4.41) 72 4.4 valutazione delle prestazioni del sistema di controllo in piccolo Confrontando le espressioni appena ottenute con le (4.34) si hanno le seguenti corrispondenze tra le due diverse parametrizzazioni del controllore P.I.D. (4.35): k11 = Ki τ −1 K12 = Kd τ −1 λ1 = τ (4.42) λ2 = − τ K2 = K p + K d τ −1 E’ evidente che la trattazione appena svolta può essere ripetuta utilizzando un controllore PID affiancato ad un altro compensatore lineare, giungendo alle medesime conclusioni, pertanto un generico compensatore comprensivo di un azione integrale soddisfa in pieno i requisiti necessari per la soluzione del problema di regolazione locale. Osservazione 4.2: Quanto mostrato mette bene in evidenza quali siano i ruoli delle tre azioni di controllo del regolatore P.I.D. [1]. Infatti l’azione integrale rappresenta il modello interno, mentre l’azione proporzionale e l’azione derivativa rappresentano la parte stabilizzante del controllore. Si noti che il complesso dell’azione derivativa e dell’azione proporzionale di fatto può essere interpretata come l’azione di controllo ottenuta tramite l’utilizzo di un osservatore di ordine ridotto. In particolare se il modello del processo allargato fosse del secondo ordine allora l’azione proporzionale congiunta alla azione derivativa rappresenterebbe esattamente l’azione di un osservatore asintotico di ordine ridotto dello stato di tale processo. L’aggiunta di ulteriori azioni dinamiche nel controllore fa si che l’azione di stima dell’osservatore sia ricalcata il più possibile. 4.4 valutazione delle prestazioni del sistema di controllo in piccolo La trattazione affrontata fin ora, è stata condotta prevalentemente utilizzando la procedura di linearizzazione come metodo per garantire la stabilizzazione esponenziale locale dei punti di equilibrio di interesse. Tuttavia la procedura di linearizzazione può essere utilizzata anche come strumento per la valutazione delle prestazioni del sistema di controllo, in termini di risposta forzata, almeno ingressi di ampiezza modesta, ossia attraverso l’utilizzo della linearizzazione è possibile condurre una analisi sul comportamento del sistema di controllo per piccoli segnali. Per far ciò si indichi con r l’ingresso di riferimento e si consideri il sistema di controllo ΣH in (4.28), il quale può essere espresso nella forma seguente dopo aver sostituito nella dinamica del processo l’espressione del controllo u(t). ẋ = f ( x, ϕ , d) ε̇ = H ε + η r ϕ̇ = Ψ ϕ + Λ ( L ε − h( x )) (4.43) 73 4.4 valutazione delle prestazioni del sistema di controllo in piccolo ebbene, indicando con ( x̂, ϕ̂, ε̂ ) il punto di equilibrio per il sistema di controllo in corrispondenza ai valori costanti y∗ , d∗ per gli ingressi r, d e ponendo: x̃ = x − x̂ (4.44) ε̃ = ε − ε̂ (4.46) ϕ̃ = ϕ − ϕ̂ (4.45) ∗ r̃ = r − y d˜ = d − d∗ (4.47) (4.48) ỹ = y − y∗ (4.49) allora per piccoli scostamenti dalla condizione di equilibrio la dinamica del sistema di controllo ΣH è bene approssimata dal seguente sistema LTI2 : x̃˙ = ( A − B K2 C ) x̃ + B K1 ϕ̃ + B K2 Lε̃ + Bd d˜ ε̃˙ = H ε̃ + η r̃ ϕ̃˙ = Ψ ϕ̃ − Λ C x̃ + Λ L ε̃ (4.50) ỹ = −C x̃ dove le matrici A, B, C sono definite come di consueto mentre la matrice Bd è definita come segue: pertanto ponendo: Bd := ∂ f ∂d x̂, ϕ̂, dˆ (4.51) al fine di ottenere una rappresentazione più compatta è possibile descrivere il sistema di controllo ΣH linearizzato nel dominio di Laplace, ossia attraverso una rappresentazione con matrici di trasferimento. In particolare indicando con Ỹ (s), D̃ (s), R̃(s) rispettivamente le trasformate di Laplace delle funzioni ỹ(t), d˜(t), r̃ (t) e ponendo: P ( s ) = C ( s I − A ) −1 B G ( s ) = K1 ( s I − Ψ ) −1 Λ + K2 Pyd (s) = C (s I − A)−1 Bd H(s) = (s I − H )−1 η (4.52) (4.53) (4.54) (4.55) (4.56) si ottiene il seguente legame tra le trasformate di Laplace dell’uscita controllata e le trasformate di Laplace dei segnali esogeni ỹ(t) e r̃ (t): Ỹ (s) = P(s) G (s) 1 H(s) R̃(s) + P (s) D̃ (s) 1 + P(s) G (s) 1 + P(s) G (s) yd tali relazioni possono essere espresse in maniera più compatta introducendo due particolari funzioni di trasferimento, ossia la funzione di sensitività S(s) 2 Per giungere alla seguente rappresentazione si faccia riferimento al sistema (4.28) nel quale compare in maniera esplicita l’espressione della legge di controllo 74 4.4 valutazione delle prestazioni del sistema di controllo in piccolo e la funzione di sensitività complementare T (s), definite rispettivamente come segue: S(s) = 1 1 + P(s) G (s) T (s) = P(s) G (s) 1 + P(s) G (s) (4.57) ottenendo in definitiva il seguente legame tra le suddette grandezze: Ỹ (s) = T (s) H(s) R̃(s) + S(s) Pyd (s) D̃ (s) (4.58) in tal modo il comportamento del sistema di controllo, per piccoli segnali, può essere dedotto facendo riferimento al classico schema di controllo. Ottenuta questa descrizione del sistema di controllo, è possibile fare riscorso ad alcuni indicatori utilizzati nella teoria classica del controllo, per stabilire quali sia il comportamento in piccolo del sistema di controllo. In particolare, al fine di valutare la rapidità di risposta del sistema di controllo ad una variazione a gradino della temperatura di riferimento, si può fare riferimento alla banda passante delle funzione di risposta armonica H( j ω ) T ( j ω ), mentre si può fare ricorso allo studio del modulo della funzione di risposta armonica della funzione dei sensitività, ossia |S( j ω )| per valutare quale sia l’effetto delle variazioni della funzione di trasferimento P(s) sul comportamento del sistema di controllo. Questo ultimo fatto risulta molto importante nel caso in esame, in quanto a fronte di una variazione del punto di lavoro del processo, certamente corrisponderà una variazione della funzione di trasferimento P(s). Inoltre osservando la (4.58), si può evincere che lo studio della funzione di sensitività fornisce un indicazione, almeno in prima approssimazione, dell’effetto delle variazioni del disturbo d˜(t) sull’errore di regolazione ỹ. Si tenga bene presente infine che le indicazioni ottenute tramite lo studio dell’approssimazione lineare del sistema di controllo, rappresentano solamente una stima del comportamento reale del sistema di controllo, tuttavia sembra ragionevole cercare di ottenere delle buone prestazioni in piccolo al fine di ottenerle per il reale sistema di controllo. 75 5 N S I N T E S I D E I R E G O L ATO R I P E R LO S C A M B I ATO R E D I C A LO R E el capitolo precedente sono stati forniti alcuni strumenti di controllo necessari per la risoluzione del problema di regolazione asintotica locale di un sistema non lineare, principalmente tali strumenti sono incentrati attorno al Teorema A.1. Tuttavia questi strumenti possono essere utilizzati per risolvere il Problema 4.1, a patto di rilassare alcune specifiche e considerare i vincoli di ammissibilità imposti in tale problema. Infatti il Problema 4.1 differisce per alcuni aspetti dal Problema 4.3. Più precisamente le differenze tra i due problemi sono principalmente tre: 1. Il Problema 4.3 è di tipo locale, mentre il Problema 4.1 di per se è di tipo globale. 2. Per risolvere il Problema 4.3 è stata imposta l’Assunzione 4.1 che non è detto sia verificata nel Problema 4.1 3. Nel Problema 4.3 non ci sono vincoli sulla variabile di controllo mentre nel Problema 4.1 tali vincoli sono invece presenti A fronte delle constatazioni appena riportate è necessario fare delle opportune osservazioni. In particolare riguardo alla prima delle tre constatazioni, si noti che essa impone che la soluzione determinata a partire dagli strumenti di controllo presentati, sarà necessariamente di tipo locale. Ciononostante però questo fatto non è un grosso problema in quanto nella ( 3.1.1 a pagina 29) è stato mostrato che il comportamento dell’impianto non è fortemente non lineare, ovvero ciò suggerisce intuitivamente che: anche se la soluzione ottenuta per il problema di regolazione della temperatura è di tipo locale, la regione di attrazione del generico equilibrio può essere ritenuta sufficientemente ampia da ritenere questo tipo di soluzione ingegneristicamente accettabile. Anche il secondo dei punti indicati di fatto non rappresenta un problema in quanto nel caso in esame l’Assunzione 4.1 è banalmente soddisfatta dato che ∗ e una coppia ( ϑ , ϑ ), purchè tale fissata una temperatura di riferimento ϑlu v li scelta dia luogo ad un valore della velocità del liquido v che sia ammissibile, è sempre possibile trovare in corrispondenza una e una sola coppia ( x, v) che sia un equilibrio dell’impianto e che fornisca una temperatura del liquido in m uguale a quella di riferimento ϑ ∗ . Viceversa la terza delle consideuscita ϑlu lu razioni suddette, di per se rappresenta un problema non banale in quanto la generica legge di controllo, sintetizzata tramite gli strumenti presentati nel capitolo precedente, garantisce che il valore del controllo sia ammissibile solamente ∗ , ϑ , ϑ in modo tale che in in senso asintotico, ossia garantisce che fissati ϑlu v li corrispondenza vi sia un valore di velocità v̂ ammissibile allora: lim v(t) = v̂ t→ +∞ 76 5.1 verifica delle condizioni per la sintesi del controllore viceversa non è noto quello che accade nel generico istante di tempo finito. Questo fatto rappresenta un grosso problema in quanto l’effetto della saturazione introdotta dalla servovalvola, negli istanti in cui il controllo assume dei valori non ammissibili, non è stato in nessun modo preso in considerazione, ciò implica che possano instaurarsi dei fenomeni oscillatori a priori non previsti! Tuttavia se il punto di lavoro che si desidera raggiungere è ammissibile, è evidente che il fenomeno di saturazione può solamente presentarsi durante il transitorio del sistema di controllo, perciò se per avventura, non esiste nessun istante della risposta transitoria in cui il controllo assume dei valori non ammissibili, allora l’effetto della saturazione può essere del tutto ignorato garantendo che la strategia di controllo, purché lo stato iniziale del sistema di controllo sia sufficientemente prossimo al punto di lavoro desiderato, permetta di ottenere la regolazione asintotica. E’ evidente che analiticamente, con gli strumenti presentati, non si può garantire che il controllo sia ammissibile per ogni istante di tempo, d’altra parte è anche evidente che se lo stato iniziale del sistema di controllo è sufficientemente prossimo al punto di lavoro desiderato e se la legge di controllo è progettata in modo da non generare grosse sovraelongazioni della variabile di controllo, allora ingegneristicamente si può ritenere alquanto improbabile l’attivazione della saturazione. Quanto detto riguardo alla mancata conoscenza a priori del comportamento del sistema di controllo, mette in evidenza un fatto molto importante e di “stampo” prettamente ingengeristico, ossia la legge di controllo deve necessariamente essere validata attraverso numerose simulazioni, in modo da poter stabilire quale sia il comportamento del sistema di controllo nel maggior numero di situazioni. Avendo messo bene in chiaro quali siano le limitazioni derivanti dall’utilizzo delle metodologie di controllo presentate nel capitolo precedente, si può passare a l’utilizzo di tali metodologie per la soluzione del Problema 4.1. 5.1 verifica delle condizioni per la sintesi del controllore Si osservi prima di ogni altra cosa che la trattazione svolta fino ad ora è applicabile al processo in esame a patto che siano verificate sia l’assunzione 4.1 sia le condizioni dettate nel Teorema A.1. A tal scopo si noti che per il modello ( 2.5 a pagina 20), è stato mostrato che fissata una terna (ϑlu , ϑli , ϑv ) ammissibile esiste sempre un valore del velocità del fluido v tale da mantenere la temperatura del fluido in uscita ad un valore costante e pari a ϑlu , ed inoltre è stato mostrato che in corrispondenza a tale valore di velocità v esiste sempre un solo equilibrio per il sistema a dimensione finita ( 2.20 a pagina 25), pertanto assumendo che il modello ad elementi finiti approssimi in maniera sufficientemente accurata il modello ( 2.5 a pagina 20) si può ritenere ragionevolmente che il suddetto valore della velocità del fluido sia unico, ossia l’assunzione 4.1 può essere ritenuta ingegneristicamente verificata nel caso in esame. Per quanto riguarda le ipotesi dettate nel Teorema A.1 si tenga presente che esse dovrebbero essere verificate per infiniti punti di lavoro il che è chiaramente impossibile, tuttavia si osservi che le condizioni imposte dal teorema in questione, sono delle condizioni per la stabilizzazione, pertanto se esse sono verificate almeno dal punto di lavoro 77 5.2 controllo pid dello scambiatore di calore nominale è possibile sintetizzare una legge di controllo basandosi sul modello nominale, cercando di garantire il più possibile il mantenimento della proprietà della stabilità anche a fronte di variazioni del punto di lavoro, a posteriori delle sintesi del controllore sarà sufficiente effettuare delle analisi volte a valutare il margine di robustezza del sistema di controllo. Per quanto detto si consideri il ben noto modello dinamico dello scambiatore di calore: ẋ (t) = A(v(t)) x (t) + B(v(t)) w(t) x ( t0 ) = x0 y(t) = C x (t) x (t) ∈ R71 v(t) ∈ R (5.1) w(t) ∈ R2 A(v) ∈ R71×71 B(v) ∈ R71×2 con ciò linearizzando tale modello attorno al punto di lavoro nominale, ( x̄, v̄ , w̄), si ottiene il seguente modello LTI: x̃˙ = A x̃ + B ṽ (5.2) ỹ = C x̃ (5.3) ove sono state effettuare per comodità le seguenti posizioni: A := A(v̄) ∂ A B := ∂v x̄, v̄, w̄ (5.4) Pertanto le condizioni per verificare se le ipotesi del Teorema A.1 siano soddisfatte si deve verificare che: 1. La tripla ( A, B, C ) è stabilizzabile e rilevabile A B 2. rank = n+1 C D Ebbene verificando tali condizioni attraverso un programma di calcolo si evince che esse sono entrambi valide e con ciò è possibile sintetizzare un controllore K LTI del tipo indicato in (A.23) in grado di risolvere il problema di regolazione locale. Prima di procedere, per quanto detto in precedenza in merito all’analisi del processo per piccoli segnali, è conveniente anche linearizzare l’impianto attorno al punto di lavoro nominale considerando le relative variazioni delle temperature del liquido in ingresso e della temperatura del vapore rispetto ai loro valori nominali. 5.2 controllo pid dello scambiatore di calore Il primo dei controllori sintetizzati per la soluzione del problema di regolazione locale è un controllore PID, come già mostrato in precedenza, tale controllore possiede tutti i requisiti per essere un compensatore candidato per la 78 5.2 controllo pid dello scambiatore di calore soluzione del problema in questione. In particolare la sintesi di tale compensatore consiste nel determinare i quattro parametri della seguente funzione di trasferimento: Ki Kd C (s) = K p + + s 1+τs in modo tale da garantire che il modello nominale (5.2) sia asintoticamente stabile e in modo che si ottengano dei buoni margini di robustezza. In particolare la sintesi di tale controllore, come già detto in precedenza può essere fatta utilizzando le tecniche di sintesi in frequenza, basate sulla funzione di trasferimento del processo nominale (5.2), che verrà indicata con P(s). Dunque al fine di ottenere almeno una stima di primo tentativo dei parametri da determinare è possibile utilizzare per il modello nominale uno dei modelli approssimati determinati nel Capitolo 3, ad esempio il modello ( 3.3 a pagina 34). Effettuando vari tentativi e utilizzando gli strumenti messi a disposizione dai programmi di calcolo si perviene alle seguente funzione di trasferimento per il controllore PID: C (s) = −0.4339 s2 − 0.267 s − 0.03998 s2 + 15.15 s (5.5) Alla quale corrispondo i seguenti valori per i parametri sopra citati: K p = −0.017449 (5.6) Ki = −0.0026388 (5.7) τ = 0.066007 (5.9) Kd = −0.02749 (5.8) L’adozione di tale controllore nel controllo dell’impianto (5.2) consente di ottenere i seguenti valori per il margine di fase m ϕ e per il margine di guadagno mg : m ϕ = 99.2◦ (5.10) m g = +∞ (5.11) Per comprende quale sia l’azione espletata dal controllore PID sintetizzato sul processo è conveniente osservare il diagramma di Bode della risposta armonica del medesimo, mostrata in (46). Per valutare l’effetto del controllore PID sulla dinamica del sistema linearizzato, in (47) viene mostrato il diagramma di Bode della funzione di trasferimento di anello L(s) = C (s) P(s), mentre in (48) è mostrato il diagramma di Nyquist della funzione L(s). Per concludere la sintesi del sistema di controllo, è necessario sintetizzare il sistema H al fine di ridurre le sollecitazioni trasmesse all’attuatore a fronte di variazioni a gradino della temperatura di riferimento y∗ . Inoltre, come detto in precedenza, è possibile attraverso una opportuna scelta del sistema H ottenere anche dei benefici in termini di risposta transitoria dell’errore di regolazione, pertanto la sintesi di tale sistema verrà effettuata anche prendendo in considerazione tale obiettivo. Attraverso tecniche di tipo trial-&-error basate su ripetute simulazioni numeriche, si perviene alla seguente funzione di trasferimento per il sistema H. H(s) = s2 0.25 + 0.23 s + 0.25 (5.12) 79 5.2 controllo pid dello scambiatore di calore Bode Diagram Magnitude (dB) 0 −10 −20 −30 Phase (deg) −40 270 225 180 135 90 −2 10 −1 10 0 1 10 10 Frequency (rad/sec) 2 3 10 10 Figura 46.: Diagramma di Bode della risposta armonica del controllore PID Bode Diagram 100 L (s) Magnitude (dB) 50 P (s) 0 −50 −100 −150 Phase (deg) −200 180 0 −180 −2 10 −1 10 0 1 10 10 Frequency (rad/sec) 2 10 3 10 Figura 47.: Diagramma di Bode della funzione d’anello L(s) Al fine di valutare quale sia l’effetto dell’azione in avanti nel sistema di controllo si osservi in (Figura 49) che l’introduzione dell’azione in avanti produce un aumento della banda passante del legame tra r̃ e ỹ, in particolare con l’introduzione dell’azione in avanti, la banda passante passa da circa 0.17 rad s−1 a circa 0.58 rad s−1 . Per quanto detto prima al fine di valutare, in prima approssimazione l’effetto delle variazioni della temperatura del fluido in ingresso e della temperatura del vapore sull’errore di regolazione e l’effetto delle variazioni del 80 5.2 controllo pid dello scambiatore di calore Nyquist Diagram 1.5 1 Imaginary Axis 0.5 0 −0.5 −1 −1.5 −1.5 −1 −0.5 0 Real Axis 0.5 1 1.5 Figura 48.: Diagramma di Nyquist della funzione d’anello L(s) punto di lavoro sul comportamento del sistema di controllo, si può prendere in considerazione il modulo della funzione di risposta armonica associato alla funzione di sensitività S(s), il quale è mostrato in (50). Per valutare la robustezza del sistema di controllo sintetizzato si riportano in (Figura 51) i diagrammi di Bode delle funzioni di anello Li (s) ottenute considerando i modelli linearizzati attorno ai punti di equilibrio indicati nel corso della trattazione. Si riportano inoltre in (Tabella 2) i margini di guadagno, i margini di fase e le pulsazioni di attraversamento, ottenute per ogni punti di lavoro, considerato. 81 5.2 controllo pid dello scambiatore di calore Bode Diagram T (s) Magnitude (dB) 50 T (s) H (s) 0 −50 −100 −150 −200 Phase (deg) 0 −90 −180 −270 −360 −2 10 −1 10 0 10 Frequency (rad/sec) 1 2 10 10 Figura 49.: Diagrammi di Bode della funzione di sensitività complementare T (s) e della funzione H(s) T (s) Bode Diagram 10 0 Magnitude (dB) −10 −20 −30 −40 −50 −3 10 −2 10 −1 0 10 10 Frequency (rad/sec) 1 10 2 10 Figura 50.: Diagramma di Bode del modulo della funzione di sensitività S(s) 82 5.2 controllo pid dello scambiatore di calore Bode Diagram 60 L 1 (s) Magnitude (dB) 40 20 L 2 (s) 0 L 3 (s) −20 L 4 (s) −40 L 0 (s) −60 Phase (deg) −80 −45 −90 −135 −180 −2 10 −1 0 10 1 10 10 Frequency (rad/sec) 2 10 3 10 Figura 51.: Diagrammi di Bode delle funzioni anello ottenute per i diversi punti di lavoro. Le funzioni di trasferimento nel diagramma sono indicate con Li (s) con i ∈ {0, 1, 2, 3, 4} a seconda del punto di lavoro considerato, in maniera consistente alla numerazione riportata in (Tabella 1) Tabella 2.: Valore dei margini di guadagno e di fase ottenuti per i punti di lavoro scelti con il controllo PID. Il punto di lavoro nominale viene indicato con il numero 0. Punto di lavoro # mg db 0 1 2 3 4 ∞ ∞ ∞ ∞ ∞ mϕ (deg) 99 47 106 96 81 ωt (rad s−1 ) 0,221 0,208 0,078 0,414 0,143 83 5.3 5.3 controllo pid con rete compensatrice controllo pid con rete compensatrice I risultati ottenuti con il controllore PID (5.5) rappresentano un buon compromesso tra rapidità di risposta e robustezza alle variazioni parametriche indotte dalle variazioni del punto di lavoro. Tuttavia si nota in (Figura 51) che attraverso l’introduzione di una rete compensatrice si può aumentare lievemente la pulsazione di attraversamento della funzione L(s), ottenendo in tal modo un miglioramento delle prestazioni del sistema di controllo, anche in situazioni diverse dalla nominale. In particolare il miglioramento desiderato è in termini di attenuazione degli effetti indotti, sull’errore di regolazione, da parte dalle variazioni della temperatura del liquido in ingresso e della temperatura del vapore, tali benefici ovviamente devono essere ottenuti senza ridurre drasticamente i margini di stabilità del sistema di controllo ottenuti attraverso il controllore PID (5.5). In definitiva la funzione di compensazione sintetizzata è la seguente: G (s) = 0.7979 s + 0.1535 s + 0.1094 (5.13) in altre parole la funzione di compensazione è una funzione ritardatrice e viene progettata con lo scopo di ottenere un aumento del modulo della funzione d’anello L(s) solamente a bassa frequenza, in modo da lasciare praticamente inalterate le prestazioni nel caso nominale, in quanto esse già abbastanza soddisfacenti. Per capire quale sia l’effetto sulla funzione di anello L(s) del compensatore G (s) in (Figura 52) è mostrato in diagramma di Bode del medesimo. Mentre in (Figura 53) è mostrato il compensatore risultante dall’azione combinata dal controllore PID (5.5) e della rete ritardatrice (5.13). Per quanto riguarda l’azione di feed-forward H(s) si sceglie di utilizzare la medesima azione utilizzata nel controllo PID, ossia: H(s) = s2 0.25 + 0.23 s + 0.25 (5.14) Conclusa la sintesi del sistema di controllo, al fine di valutare i risultati ottenuti in (Figura 54) e in (Figura 55) si mostrano, rispettivamente i diagrammi di Bode e di Nyquist della funzione d’anello L(s) = G (s) C (s) P(s) nel caso nominale. Mentre al fine di valutare, in prima approssimazione l’effetto delle variazioni della temperatura del fluido in ingresso e della temperatura del vapore sull’errore di regolazione e l’effetto delle variazioni del punto di lavoro sul comportamento del sistema di controllo, si riporta in (Figura 56) il diagramma di Bode del modulo della funzione di sensitività nel caso nominale. Per quanto riguarda l’analisi della risposta nell’errore di regolazione del sistema di controllo a fronte di piccole variazioni della temperatura di riferimento si mostra in (Figura 57) il diagramma di Bode della funzione di trasferimento T (s) H(s) nel caso nominale. In particolare osservando la (Figura 57) si può constatare che la banda passante del sistema di controllo nominale, per quanto riguarda il legame tra la variazione delle temperatura di riferimento e l’errore di regolazione è pari a: ωb = 0.5550 rad s−1 Si riportano inoltre in (Tabella 3) i margini di guadagno, i margini di fase e le pulsazioni di attraversamento, ottenute per ogni punti di lavoro, considerato. 84 5.3 controllo pid con rete compensatrice Bode Diagram Magnitude (dB) 3 2 1 0 −1 −2 −3 5 Phase (deg) 0 −5 −10 −15 −20 −3 10 −2 10 −1 10 Frequency (rad/sec) 0 1 10 10 Figura 52.: Diagrammi di Bode della funzione di compensazione G (s) Magnitude (dB) Bode Diagram 20 G(s) C (s) 10 C (s) 0 −10 −20 −30 Phase (deg) −40 270 225 180 135 90 −2 10 −1 10 0 1 10 10 Frequency (rad/sec) 2 10 3 10 Figura 53.: Diagrammi di Bode del controllore C (s) G (s) Infine per valutare la robustezza del sistema di controllo sintetizzato si ripor- 85 5.3 controllo pid con rete compensatrice Bode Diagram 50 Magnitude (dB) L (s) 0 P (s) −50 −100 Phase (deg) −150 180 90 0 −90 −180 −2 10 −1 10 0 1 2 10 10 Frequency (rad/sec) 3 10 10 Figura 54.: Diagrammi di Bode della funzione di trasferimento P(s) a confronto con la funzione d’anello L(s) = P(s) C (s) G (s) Nyquist Diagram 1.5 Imaginary Axis 1 0.5 0 −0.5 −1 −1.5 −1.5 −1 −0.5 0 Real Axis 0.5 1 1.5 Figura 55.: Diagramma di Nyquist della funzione di anello L(s) nel caso nominale 86 5.3 controllo pid con rete compensatrice Bode Diagram 5 0 −5 −10 Magnitude (dB) −15 −20 −25 −30 −35 −40 −45 −50 −3 10 −2 −1 10 0 1 10 10 Frequency (rad/sec) 10 Figura 56.: Diagramma di Bode della funzione sensitività nel caso nominale Bode Diagram Magnitude (dB) 50 0 −50 −100 −150 −200 Phase (deg) 0 −90 −180 −270 −360 −2 10 −1 10 0 10 Frequency (rad/sec) 1 10 2 10 Figura 57.: Diagramma di Bode della funzione T (s) H(s) nel caso nominale tano in (Figura 58) i diagrammi di Bode delle funzioni di anello Li (s) ottenute 87 5.3 controllo pid con rete compensatrice Tabella 3.: Valore dei margini di guadagno e di fase ottenuti per i punti di lavoro scelti con il controllo PID+ Rete compensatrice. Il punto di lavoro nominale viene indicato con il numero 0. Punto di lavoro # mg db mϕ (deg) 0 1 2 3 4 ∞ ∞ ∞ ∞ ∞ 86,4 33,4 94,2 91,7 65,1 ωt (rad s−1 ) 0,208 0,204 0,095 0,326 0,15 considerando i modelli linearizzati attorno ai punti di equilibrio indicati nel corso della trattazione. Bode Diagram 100 Magnitude (dB) L 1 (s) L 2 (s) 50 L 3 (s) 0 L 4 (s) L 0 (s) −50 Phase (deg) −45 −90 −135 −180 −2 10 −1 10 0 1 10 10 Frequency (rad/sec) 2 10 3 10 Figura 58.: Diagrammi di Bode delle funzioni anello ottenute per i diversi punti di lavoro. Le funzioni di trasferimento nel diagramma sono indicate con Li (s) con i ∈ {0, 1, 2, 3, 4} a seconda del punto di lavoro considerato, in maniera consistente alla numerazione riportata in (Tabella 1) 88 5.4 controllo lqg 5.4 controllo lqg L’ultimo dei controllori che viene proposto è un controllore di tipo observer based, in particolare un controllore LQG, progettato sul modello del processo linearizzato attorno al punto di lavoro nominale. Tale metodologia di controllo, come illustrato precedentemente, consiste nel progettare un controllore, basato sulla metodologia ad assegnazione degli autovalori con osservatore asintotico dello stato per il processo allargato definito come segue: x̃˙ = A x̃ + B k11 ϕ̃1 + B v ϕ̃˙ 1 = −λ1 C x̃ (5.15) ỹ = −C x̃ Tale sistema avendo verificato in precedenza che le ipotesi del Teorema A.1 sono soddisfatte, è certamente stabilizzabile e rilevabile, perciò certamente esiste un controllore stabilizzante di tipo LQG. Ebbene, si osservi che come bene noto la sintesi di un controllore LQG avviene utilizzando degli algoritmi numerici che a partire da alcuni parametri di progetto e dal modello del processo sotto controllo, forniscono il compensatore cercato, pertanto dato che il processo allargato è di dimensione considerevole è necessario effettuare, se possibile, una riduzione dell’ordine di tale modello e dato che la dimensione elevata del processo allargato deriva naturalmente dalla dimensione elevata del processo soggetto al controllo, è evidente che la riduzione dell’ordine andrà effettuata relativamente al solo processo da controllare linearizzato attorno al punto di lavoro nominale. A tal scopo è possibile effettuare una analisi ai valori singolari del processo in esame, il quale esito è mostrato in (Figura 59). Osservando la (Figura 59) si nota che l’ordine del modello, data la presenza di numerosi modi naturali aventi scarsa influenza sul comportamento ingresso uscita del processo, può essere ridotto notevolmente attraverso una realizzazione in forma bilanciata. In particolare una scelta opportuna prevede di ridurre il modello del processo da dimensione 71 a dimensione 6, a seguito di tale scelta si può osservare in (Figura 60) come la risposta in frequenza del modello originario sia sovrapponibile a quella modello di ordine ridotto, la cui funzione di trasferimento viene indicata con Pr (s). Ottenuto il modello di ordine ridotto, è possibile sintetizzare attraverso una procedura trial & error e attraverso l’uso di opportuni algoritmi numerici il controllore LQG sopra citato, pertanto includendo in fine nel compensatore LQG sintetizzato la parte integrale del controllore come indicato in (Figura 44), ottenendo un controllore descritto dalla seguente funzione di trasferimento: C (s) = 10.4126(s + 1)(s + 0.4477)(s2 + 0.6461s + 0.3566)(s2 + 0.4299s + 1.169) s(s2 + 0.1777s + 0.6906)(s2 + 1.071s + 1.553)(s2 + 36.29s + 657.4) (5.16) Per comprendere l’azione del controllore progettato, si riporta in (Figura 61) il diagramma di Bode della funzione di trasferimento C (s) di cui sopra. Osservando il diagramma di Bode in (Figura 61) si nota che il controllore progettato è fortemente legato al modello nominale del processo, in quanto tale controllore tende a imporre una dinamica al sistema di controllo compensando alcuni fenomeni presenti nella dinamica del processo nominale, pertanto risulta necessario 89 5.4 controllo lqg Hankel Singular Values 45 40 35 30 abs 25 20 15 10 5 0 0 10 20 30 40 Order 50 60 70 80 Figura 59.: Valori singolari della funzione di trasferimento del processo linearizzato attorno al punto di lavoro nominale operare ad una riduzione dell’ordine di tale controllore in modo da evitare di ottenere un controllore che sia poco robusto. A tal scopo si può operare una riduzione dell’ordine basata sulle realizzazioni in forma bilanciata, come quanto è stato fatto in precedenza per il processo, in particolare in (Figura 62) si mostrano i valori singolari del controllore. osservando la (Figura 62) si nota che si può ad esempio ottenere un modello del controllore ridotto del terzo ordine, in tal modo si ottiene il controllore descritto dalla seguente funzione di trasferimento. Cr (s) = 10.42s2 + 12.47s + 1.802 s3 + 36.28s2 + 653.5s (5.17) Per comprendere l’effetto della riduzione dell’ordine in (Figura 63) vengono mostrati i diagrammi di Bode della funzione di trasferimento C (s) confrontata con la funzione di traferimento Cr (s). La sintesi del controllore LQG si può ritenere ultimata, tuttavia si nota attraverso numerose simulazioni che il controllore progettato, data l’assenza di un legame diretto ingresso uscita risulta poco performante rispetto ai controllori mostrati in precedenza. Per risolvere tale problema si può aggiungere un legame diretto ingresso-uscita al controllore, tuttavia si noti che tale modifica potrebbe comportare una perdita di stabilità del sistema di controllo in quanto come è stato mostrato la matrice dinamica del sistema di controllo linearizzato attorno ad una generico punto di lavoro, dipende anche dal legame diretto ingresso-uscita K2 del controllore. Tuttavia quanto detto non presenta un problema a livello progettuale, in quanto attraverso delle simulazioni è possibile 90 5.4 controllo lqg Bode Diagram 40 P r (s) Magnitude (dB) 20 P (s) 0 −20 −40 −60 Phase (deg) −80 180 135 90 45 0 −2 10 −1 10 0 10 Frequency (rad/sec) 1 10 2 10 Figura 60.: Diagrammi di Bode della funzione di trasferimento del processo nominale e della sua approssimazione di ordine 6 determinare un valore del suddetto legame diretto ingresso-uscita che non pregiudichi la stabilità del sistema di controllo ma che ne migliori le prestazioni. La procedura appena descritta conduce alla seguente scelta: K2 = 0.01 in tal modo la funzione di trasferimento del controllore LQG con legame diretto è la seguente: Cd (s) = 0.01s3 + 10.78s2 + 19.01s + 1.802 s3 + 36.28s2 + 653.5s (5.18) Per ultimare la sintesi del sistema di controllo non resta che progettare l’azione in avanti H(s), anche in questo caso tale azione viene progettata attraverso una procedura trial & error, giungendo alla seguente funzione di trasferimento: H(s) = s2 0.16 + 0.16s + 0.16 Completata la sintesi del sistema di controllo si mostrano, come di consueto i diagrammi di Bode del compensatore Cd (s) (Figura 64), della funzione d’anello L(s) (Figura 65),della funzione di sensitività S(s) (Figura 66), della funzione di trasferimento tra la variazione di riferimento e la variazione dell’errore di regolazione ossia H(s) T (s) (Figura 67) mentre in (Figura 68) si mostra in diagramma di Nyquist della funzione d’anello L(s) nel caso nominale. In particolare osservando la (Figura 67) si può constatare che la banda passante del 91 5.4 controllo lqg Bode Diagram Magnitude (dB) 0 −10 −20 −30 Phase (deg) −40 270 225 180 135 90 −2 10 −1 10 0 1 10 10 Frequency (rad/sec) 2 10 3 10 Figura 61.: Diagrammi di Bode della funzione di trasferimento del compensatore LQG con azione integrale sistema di controllo nominale, per quanto riguarda il legame tra la variazione delle temperatura di riferimento e l’errore di regolazione, è pari a: ωb = 0.5674 rad s−1 Ai fini di valutare la robustezza e in una qualche misura le perfomance del sistema di controllo, si riportano inoltre in (Tabella 4) i margini di guadagno, i margini di fase e le pulsazioni di attraversamento, ottenute per ogni punti di lavoro, considerato. Tabella 4.: Valore dei margini di guadagno e di fase ottenuti per i punti di lavoro scelti con il controllo LQG. Il punto di lavoro nominale viene indicato con il numero 0. Punto di lavoro # mg db 0 1 2 3 4 48,3 50,8 49,2 44,3 54,1 mϕ (deg) 65,2 33 123 52,3 69,6 ωt (rad s−1 ) 0,416 0,26 0,116 0,582 0,234 Infine per valutare maggiormente la robustezza del sistema di controllo sintetizzato si riportano in (Figura 69) i diagrammi di Bode delle funzioni di anello Li (s) ottenute considerando i modelli linearizzati attorno ai punti di equilibrio indicati nel corso della trattazione. 92 5.4 controllo lqg Hankel Singular Values (1 of jw−axis pole(s) = inf, not shown) 0.16 0.14 0.12 abs 0.1 0.08 0.06 0.04 0.02 0 1 2 3 4 Order 5 6 7 Figura 62.: Valori singolari associati ai modi naturali presenti nel legame ingressouscita del controllore LQG Bode Diagram 20 C r (s) Magnitude (dB) 10 C (s) 0 −10 −20 −30 Phase (deg) −40 270 225 180 135 90 −3 10 −2 10 −1 10 0 1 10 10 Frequency (rad/sec) 2 10 Figura 63.: Diagramma di Bode di Cr (s) e di C (s) a confronto 3 10 93 5.4 controllo lqg Bode Diagram 20 Magnitude (dB) 10 0 −10 −20 −30 −40 Phase (deg) −50 270 225 180 135 90 −2 10 0 2 10 10 Frequency (rad/sec) 4 10 Figura 64.: Diagramma di Bode della funzione di trasferimento Cd (s) Bode Diagram Magnitude (dB) 100 0 −100 −200 −300 −45 Phase (deg) −90 −135 −180 −225 −270 −2 10 0 2 10 10 Frequency (rad/sec) 4 10 Figura 65.: Diagramma di Bode della funzione di trasferimento L(s) 94 5.4 controllo lqg Bode Diagram 10 0 Magnitude (dB) −10 −20 −30 −40 −50 −3 10 −2 10 −1 0 1 10 10 Frequency (rad/sec) 2 10 10 Figura 66.: Diagramma di Bode del modulo della funzione di sensitività S(s) Bode Diagram 100 Magnitude (dB) 0 −100 −200 −300 −400 −500 0 Phase (deg) −90 −180 −270 −360 −450 −2 10 0 2 10 10 Frequency (rad/sec) 4 10 Figura 67.: Diagramma di Bode della funzione di trasferimento H(s) T (s) 95 5.4 controllo lqg Nyquist Diagram 1.5 Imaginary Axis 1 0.5 0 −0.5 −1 −1.5 −2 −1.5 −1 −0.5 0 Real Axis 0.5 1 1.5 2 Figura 68.: Diagramma di Nyquist della funzione d’anello L(s) Bode Diagram 100 L 1 (s) Magnitude (dB) 50 L 2 (s) 0 L 3 (s) −50 L 4 (s) −100 L 0 (s) −150 −200 −45 Phase (deg) −90 −135 −180 −225 −270 −4 10 −2 10 0 2 10 10 Frequency (rad/sec) 4 10 6 10 Figura 69.: Diagrammi di Bode delle funzioni anello ottenute per i diversi punti di lavoro. Le funzioni di trasferimento nel diagramma sono indicate con Li (s) con i ∈ {0, 1, 2, 3, 4} a seconda del punto di lavoro considerato, in maniera consistente alla numerazione riportata in (Tabella 1) 96 5.5 5.5 scelta del controllore più soddisfacente scelta del controllore più soddisfacente In questa sezione si vuole scegliere tra tutti i controllori progettati, quello che si ritiene possa essere il più soddisfacente. Per tale valutazione si prenderanno in considerazione, per i tre controllori, gli andamenti delle funzioni di sensitività S(s), delle funzioni di trasferimento tra la variazione del riferimento e l’errore di regolazione H(s) T (s), i margini di stabilità e la banda passante nel caso nominale. 97 5.5 scelta del controllore più soddisfacente Margini di stabilità e banda passante Tabella 5.: Valore dei margini di guadagno e di fase ottenuti per i punti di lavoro scelti con il controllo LQG. Il punto di lavoro nominale viene indicato con il numero 0. Controlore mg db mϕ (deg) PID LQG PID+Rete ∞ 48,3 ∞ 99,2 65,2 84,6 ωb (rad s−1 ) 0,5721 0,5674 0,5550 Funzione di sensitività Bode Diagram 20 LQG PID PID+Rete 10 0 Magnitude (dB) −10 −20 −30 −40 −50 −60 −4 10 −3 10 −2 10 −1 0 10 10 Frequency (rad/sec) 1 10 2 10 3 10 Figura 70.: Diagrammi di Bode del modulo della funzione di sensitività per i tre controllori nel caso nominale Funzione di trasferimento H(s) T (s) Osservando la (Figura 70), la (Figura 71) e la (Tabella 5) , si può evincere che il controllore PID con rete compensatrice, rappresenta tra tutte la scelta per la quale si ottiene il miglior compromesso tra robustezza e performance, in particolare tale controllore garantisce delle buone prestazioni, dei buoni margini di robustezza e una buona reiezione dei disturbi a bassa frequenza. Pertanto al fine di disporre di una realizzazione del controllore candidato che sia pronta all’uso, se ne deriva una versione a tempo discreto che può essere facilemente 98 5.5 scelta del controllore più soddisfacente Bode Diagram Magnitude (dB) 50 LQG PID PID+Rete 0 −50 −100 Phase (deg) −150 0 −90 −180 −270 −360 −3 10 −2 10 −1 10 Frequency (rad/sec) 0 10 1 10 Figura 71.: Diagrammi di Bode della funzione di trasferimento tra la variazione del riferimento di temperatura e l’errore di regolazione H(s) T (s) per i tre controllori nel caso nominale implementata con un dispositivo digitale. In particolare scegliendo un tempo di campionamento Ts = 0.01 s (un tempo certamente ragionevole considerando un dispositivo di elaborazione di costo medio-basso) e discretizzando il regolatore in questione con il metodo di Tustin, si ottiene il regolatore a tempo discreto descritto dalla seguente funzione di trasferimento. Cd (z) = 0.323z3 − 0.9663z2 + 0.9637z − 0.3204 z3 − 2.858z2 + 2.716z − 0.8582 Al fine di stabilire se le performance del sistema di controllo a ottenuto con il regolatore a tempo continuo siano sufficientemente recuperate in (72), si riportano i diagrammi di Bode delle funzioni di trasferimento Cd (z) e del regolatore originario C (s) 99 5.5 scelta del controllore più soddisfacente Bode Diagram 0 C (s) C d( z ) −5 −10 Magnitude (dB) −15 −20 −25 −30 −35 −40 −2 10 −1 10 0 1 10 10 Frequency (rad/sec) 2 10 3 10 Figura 72.: Diagrammi di Bode delle funzioni di trasferimento C (s) e Cd (z) a confronto 100 6 A L C U N I R I S U LT A T I N U M E R I C I In questo capitolo vengono mostrati alcuni risultati ottenuti utilizzando i controllori sintetizzati nel corso della presente trattazione. In particolare vengono valutate le risposte del sistema di controllo a variazioni a gradino della temperatura di riferimento e a variazioni a gradino delle temperatura del liquido in ingresso e della temperatura del vapore. Nello specifico vengono imposte le variazioni illustrate in (Figura 73) e in (Figura 74). 102 100 ◦ C 98 96 94 92 90 0 100 200 300 400 500 600 700 800 900 1000 ti me [ s] Figura 73.: Andamento del riferimento y∗ (t) 160 ϑ v (t ) ϑ l i(t ) 140 ◦ C 120 100 80 60 40 0 100 200 300 400 500 600 700 800 900 1000 ti me [ s] Figura 74.: Andamento delle variabili ϑli (t) e ϑv (t) A fronte degli andamenti imposti per le grandezze esogene, controllando il processo con i controllori ottenuti nel corso della trattazione, si ottengono i risultati mostrati in (Figura 75) e in (Figura 76). I risultati mostrati di cui sopra non evidenziano molto bene l’andamento della risposta transitoria del sistema di controllo, pertanto nei paragrafi seguenti 101 6.1 controllo pid con azione in avanti 115 PID PID+RETE LQG PID+RETE discreto 110 ◦ C 105 100 95 90 85 0 100 200 300 400 500 600 700 800 900 1000 ti me [ s] Figura 75.: Andamento dell’uscita controllata y(t) 0.8 PID PID+RETE LQG PID+RETE discreto 0.7 m s −1 0.6 0.5 0.4 0.3 0.2 0.1 0 100 200 300 400 500 600 700 800 900 1000 ti me [ s] Figura 76.: Andamento del controllo u(t) verrano mostrati singolarmente, gli andamenti delle risposte transitorie di cui sopra, ottenute con tutti i controllori sintetizzati. 6.1 controllo pid con azione in avanti Risultati ottenuti con il controllo PID. 102 6.1 controllo pid con azione in avanti 99 P l a nt Ou tp u t R eferen ce 98 [ ◦C ] 97 96 95 94 93 60 80 100 120 140 160 180 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.4 0.39 0.38 v (t )[ m s −1] 0.37 0.36 0.35 0.34 0.33 0.32 0.31 0.3 90 100 110 120 130 140 150 160 Ti me [ s] (b) Andamento della variabile di controllo v Figura 77.: Variazione delle temperatura di riferimento da 98 ◦ a 94 ◦ 103 6.1 controllo pid con azione in avanti 100 P l a nt Ou tp u t R eferen ce 99 98 [ ◦C ] 97 96 95 94 93 190 200 210 220 230 240 250 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.4 v (t )[ m s −1] 0.35 0.3 0.25 200 210 220 230 240 250 260 270 Ti me [ s] (b) Andamento della variabile di controllo v Figura 78.: Variazione delle temperatura di riferimento da 94 ◦ a 98 ◦ 104 6.1 controllo pid con azione in avanti 100 P l a nt Ou tp u t R eferen ce 99 98 [ ◦C ] 97 96 95 94 93 310 320 330 340 350 360 370 380 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.3 0.28 0.26 v (t )[ m s −1] 0.24 0.22 0.2 0.18 0.16 0.14 0.12 0.1 320 325 330 335 340 345 350 355 360 365 370 375 Ti me [ s] (b) Andamento della variabile di controllo v Figura 79.: Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (65◦ , 120◦ ) a (75◦ , 140◦ ) 105 6.1 controllo pid con azione in avanti 108 P l a nt Ou tp u t R eferen ce 106 [ ◦C ] 104 102 100 98 96 510 520 530 540 550 560 570 580 590 600 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.7 0.6 v (t )[ m s −1] 0.5 0.4 0.3 0.2 0.1 520 530 540 550 560 570 580 590 Ti me [ s] (b) Andamento della variabile di controllo v Figura 80.: Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (75◦ , 140◦ ) a (50◦ , 140◦ ) 106 6.1 controllo pid con azione in avanti P l a nt Ou tp u t 100 R eferen ce 98 [ ◦C ] 96 94 92 90 88 690 700 710 720 730 740 750 760 770 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.65 v (t )[ m s −1] 0.6 0.55 0.5 0.45 0.4 0.35 0.3 700 710 720 730 740 750 Ti me [ s] (b) Andamento della variabile di controllo v Figura 81.: Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (50◦ , 140◦ ) a (75◦ , 105◦ ) 107 6.1 controllo pid con azione in avanti 104 P l a nt Ou tp u t 102 R eferen ce 100 [ ◦C ] 98 96 94 92 90 88 835 840 845 850 855 860 865 870 875 880 885 890 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.6 0.55 0.5 v (t )[ m s −1] 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 835 840 845 850 855 860 865 870 875 880 885 890 Ti me [ s] (b) Andamento della variabile di controllo v Figura 82.: Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (75◦ , 105◦ ) a (50◦ , 105◦ ) 108 controllo pid con rete compensatrice con azione in avanti 99 P l a nt Ou tp u t R eferen ce 98 97 96 95 94 93 92 90 100 110 120 130 140 150 160 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.4 0.39 0.38 0.37 v (t )[ m s −1] 6.2 6.2 controllo pid con rete compensatrice con azione in avanti 0.36 0.35 0.34 0.33 0.32 0.31 0.3 90 100 110 120 130 140 150 160 170 180 190 Ti me [ s] (b) Andamento della variabile di controllo v Figura 83.: Variazione delle temperatura di riferimento da 98 ◦ a 94 ◦ 109 6.2 controllo pid con rete compensatrice con azione in avanti 100 P l a nt Ou tp u t 99 R eferen ce 98 97 96 95 94 93 92 190 200 210 220 230 240 250 260 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.4 0.38 0.36 v (t )[ m s −1] 0.34 0.32 0.3 0.28 0.26 0.24 0.22 0.2 190 200 210 220 230 240 250 260 270 Ti me [ s] (b) Andamento della variabile di controllo v Figura 84.: Variazione delle temperatura di riferimento da 94 ◦ a 98 ◦ 110 6.2 controllo pid con rete compensatrice con azione in avanti 101 P l a nt Ou tp u t 100 R eferen ce 99 98 97 96 95 94 93 320 330 340 350 360 370 380 390 400 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.3 v (t )[ m s −1] 0.25 0.2 0.15 0.1 310 320 330 340 350 360 370 380 390 400 Ti me [ s] (b) Andamento della variabile di controllo v Figura 85.: Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (65◦ , 120◦ ) a (75◦ , 140◦ ) 111 6.2 controllo pid con rete compensatrice con azione in avanti 110 P l a nt Ou tp u t 108 R eferen ce 106 104 102 100 98 96 94 520 530 540 550 560 570 580 590 600 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.7 0.6 v (t )[ m s −1] 0.5 0.4 0.3 0.2 0.1 510 520 530 540 550 560 570 Ti me [ s] (b) Andamento della variabile di controllo v Figura 86.: Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (75◦ , 140◦ ) a (50◦ , 140◦ ) 112 6.2 controllo pid con rete compensatrice con azione in avanti 102 P l a nt Ou tp u t R eferen ce 100 98 96 94 92 90 88 690 700 710 720 730 740 750 760 770 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.6 v (t )[ m s −1] 0.55 0.5 0.45 0.4 0.35 0.3 700 710 720 730 740 750 760 Ti me [ s] (b) Andamento della variabile di controllo v Figura 87.: Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (50◦ , 140◦ ) a (75◦ , 105◦ ) 113 6.2 controllo pid con rete compensatrice con azione in avanti 102 P l a nt Ou tp u t R eferen ce 100 98 96 94 92 90 840 850 860 870 880 890 900 910 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.4 0.35 v (t )[ m s −1] 0.3 0.25 0.2 0.15 0.1 840 850 860 870 880 890 Ti me [ s] (b) Andamento della variabile di controllo v Figura 88.: Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (75◦ , 105◦ ) a (50◦ , 105◦ ) 114 controllo lqg con azione in avanti 99 P l a nt Ou tp u t R eferen ce 98 97 96 95 94 93 90 100 110 120 130 140 150 160 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.4 0.38 0.36 0.34 v (t )[ m s −1] 6.3 6.3 controllo lqg con azione in avanti 0.32 0.3 0.28 0.26 0.24 0.22 0.2 90 100 110 120 130 140 150 160 170 180 190 Ti me [ s] (b) Andamento della variabile di controllo v Figura 89.: Variazione delle temperatura di riferimento da 98 ◦ a 94 ◦ 115 6.3 controllo lqg con azione in avanti 100 P l a nt Ou tp u t R eferen ce 99 98 97 96 95 94 93 200 210 220 230 240 250 260 270 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.4 v (t )[ m s −1] 0.35 0.3 0.25 200 210 220 230 240 250 260 270 Ti me [ s] (b) Andamento della variabile di controllo v Figura 90.: Variazione delle temperatura di riferimento da 94 ◦ a 98 ◦ 116 6.3 controllo lqg con azione in avanti 105 P l a nt Ou tp u t R eferen ce 100 95 90 320 340 360 380 400 420 440 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.4 0.35 v (t )[ m s −1] 0.3 0.25 0.2 0.15 0.1 320 340 360 380 400 420 440 Ti me [ s] (b) Andamento della variabile di controllo v Figura 91.: Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (65◦ , 120◦ ) a (75◦ , 140◦ ) 117 6.3 controllo lqg con azione in avanti P l a nt Ou tp u t R eferen ce 110 105 100 95 90 510 520 530 540 550 560 570 580 590 600 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.6 v (t )[ m s −1] 0.5 0.4 0.3 0.2 0.1 510 520 530 540 550 560 570 580 590 600 610 620 Ti me [ s] (b) Andamento della variabile di controllo v Figura 92.: Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (75◦ , 140◦ ) a (50◦ , 140◦ ) 118 6.3 controllo lqg con azione in avanti P l a nt Ou tp u t 104 R eferen ce 102 100 98 96 94 92 90 88 690 700 710 720 730 740 750 760 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.65 0.6 0.55 v (t )[ m s −1] 0.5 0.45 0.4 0.35 0.3 0.25 0.2 690 700 710 720 730 740 750 Ti me [ s] (b) Andamento della variabile di controllo v Figura 93.: Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (50◦ , 140◦ ) a (75◦ , 105◦ ) 119 6.3 controllo lqg con azione in avanti P l a nt Ou tp u t 104 R eferen ce 102 100 98 96 94 92 90 88 840 850 860 870 880 890 900 910 Ti me [ s] (a) Andamento della temperatura del fluido in uscita misurata ϑlu 0.5 0.45 v (t )[ m s −1] 0.4 0.35 0.3 0.25 0.2 0.15 0.1 840 850 860 870 880 890 Ti me [ s] (b) Andamento della variabile di controllo v Figura 94.: Variazione della temperatura del fluido in ingresso e della temperatura del vapore da (75◦ , 105◦ ) a (50◦ , 105◦ ) 120 7 L CONCLUSIONI ’esito del progetto si può ritenere positivo in quanto la regolazione del processo avviene con successo anche in condizioni operative molto diverse dalla nominale. Tuttavia l’approccio utilizzato non consente di certificare quale sia l’insieme delle condizioni operative per le quali l’obiettivo della regolazione venga conseguito. Un altro grosso problema è quello della saturazione dell’attuatore che non è stata presa in considerazione in fase di sintesi del sistema di controllo. Pertanto come sviluppi futuri si potrebbe pensare di studiare il problema della stabilizzazione robusta del processo con un approccio basato sui dei criteri di stabilità robusta in presenza di incertezze parametriche, mentre per quanto riguarda la saturazione si può sintetizzare uno schema di controllo con azione di anti-windup. 121 A STRUMENTI NECESSARI PER LA S I N T E S I D E I R E G O L ATO R I In questa appendice verranno illustrati gli strumenti concettuali utilizzati nella presente trattazione, per la sintesi dei regolatori mostrati. Lemma A.1: Sia K un generico sistema LTI nella forma: ż = A z + B u (A.1) z ∈ Rn (A.2) con: u∈R (A.3) dunque se A è singolare e se ( A, B) stabilizzabile, allora preso ( x̄, ū) tale che: A x̄ + B ū = 0 necessariamente deve risultate: ū = 0 Dimostrazione. Prima di ogni cosa si noti che essendo ( A, B) stabilizzabile si ha necessariamente che: rank[ A B] = n e quindi segue anche che: dim{Ker A} = 1 pertanto il seguente sistema equazioni ammette ∞1 soluzioni. Ossia: x̄ =0 A B ū (A.4) ebbene in virtù del fatto che la matrice A è singolare si ha che: ∗ x ∗ ∀ x ∈ Ker A : [ A B] =0 0 tuttavia dato che: dim{Ker A} = 1 si ha che tutte e sole le soluzioni del sistema (A.4) sono del tipo ( x ∗ , 0) con x ∗ ∈ Ker A, infatti se per assurdo esistesse un soluzione del sistema (A.4) del tipo: ( x̂ û) ∈ Rn con û 6= 0 allora vuol dire il sottospazio delle soluzioni per il sistema (A.4) ammetterebbe come base B il seguente insieme dei vettori: ( ) x∗ 0 B= , 0 û 122 strumenti necessari per la sintesi dei regolatori ma ciò contraddirebbe il fatto che il sottospazio delle soluzioni del sistema (A.4) sia di dimensione unitaria, contraddicendo in tal modo l’ipotesi di stabilizzabilità della coppia ( A, B). Lemma A.2: In merito al sistema in (4.3), si fissino y∗ ∈ R, d∗ ∈ R p in modo qualunque ma a patto che valga la Assunzione 4.1 con ū 6= 0. Si consideri ora la seguente legge di controllo L.T.I. in retroazione dall’errore di regolazione ỹ = y∗ − y: ϕ̇ = Ψ ϕ + Λ ỹ (A.5) u = K1 ϕ + K2 ỹ Con: ϕ ∈ R` K1 ∈ R (A.6) 1×` (A.7) K2 ∈ R (A.8) Λ ∈ R` (A.10) Ψ ∈ R`×` (A.9) La quale fornisce la seguente dinamica per il sistema di controllo: ẋ = f ( x, ϕ, d∗ y∗ ) (A.11) ϕ̇ = Ψ ϕ + Λ (y∗ − h( x )) Ebbene se sotto l’azione di tale controllo esiste, per il sistema di controllo (A.11), uno stato di equilibrio ( x̂ , ϕ̂) ∈ R`+n tale che: h( x̂ ) = y∗ allora valgono le seguenti affermazioni: • ū = K1 ϕ̂. • la matrice Ψ ammette λ = 0 come autovalore ed inoltre ϕ̂ ∈ Ker (Ψ) • K1T ∈ / Im Ψ Dimostrazione. Se lo stato ( x̂ , ϕ̂) è di equilibrio per il sistema (A.11) allora si ha: f ( x̂, ϕ̂, d∗ y∗ ) = 0 Ψ ϕ̂ + Λ (y∗ − h( x̂ )) = 0 (A.12) Tuttavia in vista dell’Assunzione 4.1 si ha necessariamente: x̄ = x̂ ū = K1 ϕ̂ (A.13) 123 strumenti necessari per la sintesi dei regolatori Provando in tal modo in primo e il secondo punto dell’enunciato. Allo scopo di completare la dimostrazione del lemma si osservi che le espressioni in (A.12) si modificano come segue: f ( x̂, ϕ̂, d∗ y∗ ) = 0 (A.14) Ψ ϕ̂ = 0 (A.15) ossia il vettore ϕ̂ deve soddisfare il seguente sistema di equazioni lineari: Ψ 0 ϕ̂ = `×1 K1 ū (A.16) Per cui affinché esista almeno una soluzione di tale sistema di equazioni, occorre e basta, per il teorema di Rouchè-Capelli che: Ψ rank K1 Ψ 0`×1 = rank K1 ū (A.17) ebbene si osservi ora che dato che per ipotesi ū allora: Ψ 0`×1 = rank Ψ + 1 rank K1 ū (A.18) pertanto dalla (A.17) si ha che le colonne della matrice Ψ sono linearmente indipendenti dal vettore K1T concludendo la dimostrazione del lemma. Osservazione A.1: Il lemma appena enunciato e provato, di fatto indica che un qualunque controllore L.T.I. che risolva il problema di regolazione contiene necessariamente una azione integrale. Tale proprietà è ben nota ed in particolare essa è denominata in alcuni testi con in nome di Principio del modello interno, facendo riferimento al fatto che il sistema di controllo contiene al suo interno un modello degli ingressi esogeni da compensare, siano essi riferimenti o disturbi, i quali nel caso del problema di regolazione sono di tipo costante, ossia descritti da un sistema dinamico di tipo integrale non forzato. Infine si noti che nelle ipotesi del lemma si assume che ū 6= 0, in caso contrario il lemma perde la sua validità di condizione necessaria, tuttavia essendo in generale ū qualunque, tale condizione rimane necessaria per gli scopi preposti. Si noti infine che la condizione sul rango enunciata nel lemma implica che se il sistema (A.5) è rilevabile allora può esistere uno e un solo valore di ϕ̂ per il quale vale la proprietà di regolazione. Prima di proseguire è conveniente enunciare un ulteriore risultato. Proposizione A.1: Il sistema (A.5) è rilevabile soltanto se: rank [Ψ − λ I] ≥ ` − 1 ∀ λ ∈ C+ Dimostrazione. La dimostrazione dell’assunto è immediata e discende direttamente dal P.B.H. test di rilevabilità. Infatti il sistema (A.5) è rilevabile se e solo se: Ψ−λI rank = ` ∀ λ ∈ C+ K1 124 strumenti necessari per la sintesi dei regolatori ma ciò è possibile solo se: rank Ψ − λ I ≥ ` − 1 ∀ λ ∈ C+ dimostrando l’assunto. Proposizione A.2: Sia data la legge di controllo (A.5), allora se sotto l’azione di tale controllo esiste, per il sistema di controllo (A.11), un solo stato di equilibrio ( x̂ , ϕ̂) ∈ R`+n tale che: h( x̂ ) = y∗ e se la coppia (Ψ, K1 ) è rilevabile allora esiste una trasformazione di coordinate invertibile ζ = T ϕ tramite la quale si può esprimere il sistema (A.5) nella seguente forma: 01×1 0 ζ + Λ J ỹ ζ̇ = 0 Ψ̂ u= K1J (A.19) ζ + K2 ỹ Dimostrazione. Si indichi con JΨ la forma di Jordan della matrice Ψ, ebbene è noto che certamente esiste una matrice quadrata e non singolare TJ tale che: TJ Ψ TJ−1 = JΨ Tuttavia dalle ipotesi per il Lemma 4.1 e per la Proposizione A.1, λ = 0 è un autovalore della matrice Ψ con molteplicità geometrica unitaria, quindi la forma di Jordan della matrice Ψ può essere scritta nel modo seguente: 0 1 0 ... 0 .. . . . . .. . . 1 . . .. O . 0 0 0 0 JΨ = .. . . . . . . . . . . 1 0 0 0 0 0 O ĴΨ dove il blocco di Jordan relativo all’autovalore nullo è di dimensione pari alla molteplicità algebrica di tale autovalore. A fronte di quanto mostrato si consideri il seguente cambio di coordinate: ζ̄ = TJ x il quale fornisce la seguente rappresentazione equivalente per il sistema (A.5): 0 1 .. . . . . 0 0 ζ̄˙ = .. . . . . 0 0 0 ... .. . 1 .. . 0 .. .. . . 0 0 O u = K̄1J ζ̄ + K2 ỹ 0 .. . 0 1 0 O ζ̄ + Λ̄ J ỹ ĴΨ (A.20) 125 strumenti necessari per la sintesi dei regolatori Infine attraverso l’ulteriore trasformazione: ζ = T̂ ζ̄ con: 0 T̂ = 0 0 0 1 .. . 0 .. . . . . 0 1 .. . 0 1 0 0 0 0 0 .. . 1 .. . .. 0 . .. . 0 0 1 ... 0 O O I Si ottiene per il sistema (A.5) la seguente rappresentazione: 0 0 0 0 0 . . .. . . 1 . . . ... . . 0 0 0 . 0 O .. .. .. ζ + Λ J ỹ . . . 0 ζ̇ = 1 0 0 0 0 0 0 1 0 . . . 0 O ĴΨ (A.21) u = K1J ζ + K2 ỹ In definitiva dato che la trasformazione: ζ = T̂ Tj x conduce il sistema (A.5) nella forma (A.21) la dimostrazione dell’enunciato è conclusa. Osservazione A.2: Si osservi che nella rappresentazione (A.21) ottenuta per il sistema (A.5) la matrice ĴΨ è in forma di Jordan, tuttavia ciò non da luogo a nessuna perdita di generalità in quanto è possibile, attraverso un ulteriore cambio di coordinate, portare la matrice ĴΨ in una forma del tutto generale. Teorema A.1: In merito al sistema in (4.3), si fissino y∗ ∈ R, d∗ ∈ R p arbitrariamente a patto che valga la Assunzione 4.1. Ebbene in si indichino ora rispettivamente: ∂ f ( x , u, d∗ ) A := ∂x ( x ,u)=( x̄ ,ū) ∗ ∂ f (x , u , d ) (A.22) B := ∂u ( x ,u)=( x̄ ,ū) ∂h( x ) C := ∂x (x ,u)=( x̄ ,ū) dunque se: 126 strumenti necessari per la sintesi dei regolatori • la terna ( A, B, C ) è stabilizzabile e rilevabile • rank A B = n+1 C 0 allora esiste un controllore K nella forma (4.10) che sia soluzione del problema (4.3) ed inoltre esso è espresso nella forma seguente: ϕ̇1 = λ1 ỹ (A.23) K : ϕ̇2 = Ψ̂ ϕ2 + λ2 ỹ u = k11 ϕ1 + K12 ϕ2 + K2 ỹ con: ϕ1 ∈ R ϕ2 ∈ R`−1 λ1 ∈ R (A.24) (A.25) (A.26) λ2 ∈ R`−1 (A.27) K12 ∈ R1×`−1 (A.29) k11 ∈ R (A.28) (A.30) Dimostrazione. La dimostrazione del risultato appena enunciato è di tipo costruttivo e conduce alla costruzione del controllore che risolve il problema di regolazione in questione. Si consideri l’espressione della dinamica ad anello chiuso, ottenuta controllando il sistema (4.3) con il controllore (A.23) : ẋ = f ( x, ϕ, d∗ y∗ ) ϕ̇ = Ψ ϕ + Λ (y∗ − h( x )) (A.31) Ebbene affinché l’errore di regolazione possa annullarsi asintoticamente come richiesto, è necessario che il sistema di controllo ammetta uno stato di equilibrio ( x̂ , ϕ̂) ∈ R`+n tale che: h( x̂ ) = y∗ Dunque a fronte del Lemma A.2 è necessario per la sintesi del controllore che: 0 ∈ σ (Ψ) Ebbene a questo punto essendo arbitraria la molteplicità algebrica dell’autovalore nullo della matrice Ψ si può scegliere essa unitaria e quindi imponendo evidentemente che il compensatore sia rilevabile, dal lemma A.2, partizionando i vettori ϕ e Λ come segue: ϕ := [ ϕ1 , ϕ2 T ] T (A.32) Λ := [λ1 , λ2 T ] T (A.33) K1 := [k11 , K12 ] (A.34) si può scrivere il controllore (A.23) in questo modo: ϕ̇1 = λ1 ỹ ϕ̇2 = Ψ̂ ϕ2 + λ2 ỹ u = k11 ϕ1 + K12 ϕ2 + K2 ỹ (A.35) 127 strumenti necessari per la sintesi dei regolatori dove: ϕ1 ∈ R ϕ2 ∈ R λ1 ∈ R λ2 ∈ R (A.36) `−1 (A.37) (A.38) `−1 k11 ∈ R K12 ∈ R (A.39) (A.40) 1×`−1 (A.41) (A.42) In tal modo la dinamica del sistema di controllo può essere riscritta in tal modo: Dinamica del sistema di controllo ẋ = f ( x, ϕ1 , ϕ2 , d∗ , y∗ ) ϕ̇1 = λ1 (y∗ − h( x )) (A.43) ∗ ϕ̇2 = Ψ̂ ϕ2 + λ2 (y − h( x )) Prima di procedere si noti che la scelta fatta per la struttura del controllore impone che se esso è stabilizzabile e rilevabile, allora k11 6= 0 e λ1 6= 0. Si osservi ora che a fronte dell’Assunzione 4.1 e del Lemma 4.1 ∃( x̄ , ϕ̄1 , ϕ̄2 ) tale che: f ( x̄, ϕ̄1 , ϕ̄2 , d∗ , y∗ ) = 0 Ψ̂ ϕ̄2 = 0 ū = k11 ϕ̄1 + K12 ϕ̄2 (A.44) h( x̄ ) = y∗ Ossia il sistema di controllo ammette come punto di equilibrio lo stato che garantisce un errore di regolazione identicamente nullo, pertanto per conseguire l’obiettivo della regolazione locale, occorre e basta provvedere alla stabilizzazione asintotica, almeno locale, di tale equilibrio. Pertanto come noto dalla teoria dei sistemi dinamici non lineari, al fine di conseguire l’obiettivo della stabilità esponenziale locale di un dato equilibrio, è sufficiente rendere asintoticamente stabile l’approssimazione lineare del sistema in esame attorno all’equilibrio considerato. A tal scopo si consideri la linearizzazione, attorno al punto ( x̄ , ϕ̄1 , ϕ̄2 ), dei termini a secondo membro della (A), la quale porge le seguenti espressioni: x̃˙ = ( A − B K2 C ) x̃ + B (k11 ϕ̃1 + K12 ϕ̃2 ) ϕ̃˙ 1 = −λ1 C x̃ ϕ̃˙ 2 = Ψ̂ ϕ̃2 − λ2 C x̃ (A.45) 128 strumenti necessari per la sintesi dei regolatori dove: ∂ f ( x , u d) ∂x ( x ,u)=( x̄ ,ū) ∂ f ( x , u) B := ∂u ( x ,u)=( x̄ ,ū) ∂h( x ) C := ∂x (x ,u)=(x̄ ,ū) A := (A.46) (A.47) (A.48) x̃ = x − x̄ (A.49) ϕ̃1 = ϕ1 − ϕ̄1 (A.51) ũ = u − ū (A.50) ϕ̃2 = ϕ2 − ϕ̄2 (A.52) (A.53) Al fine di stabilizzare il sistema ottenuto, ci si può giovare della arbitrarietà del sottosistema caratterizzato dallo stato ϕ2 . Infatti si può osservare che il sistema (A.45) può essere riscritto come l’interconnessione di due sistemi L.T.I., di cui uno dei due, il quale verrà denominato processo allargato, è composto dall’impianto e dalla componente della legge di controllo caratterizzata dallo stato ϕ1 ; mentre l’altro,che verrà indicato come controllore stabilizzante, è composto dalla componente della legge di controllo caratterizzata dallo stato ϕ2 e verrà sintetizzato arbitrariamente con il solo scopo di conseguire l’obiettivo della stabilità asintotica dell’intero sistema di controllo (A.45). Quanto detto conduce alla seguente scrittura per il sistema (A.45). Processo Allargato x̃˙ = A x̃ + B k11 ϕ̃1 + B v ϕ̃˙ 1 = −λ1 C x̃ (A.54) ỹ = −C x̃ Controllore stabilizzante ϕ̃˙ 2 = Ψ̂ ϕ̃2 + λ2 ỹ (A.55) Per concludere quindi la prova del teorema occorre e basta mostrare che il processo allargato (A.54) sia stabilizzabile e rilevabile. Per cui si consideri il processo allargato riscritto in forma compatta come segue: x̃˙ A B k11 x̃ B = + v ϕ̃1 0 ϕ̃˙ 1 − λ1 C 0 x̃ ỹ = −C 0 ϕ̃1 (A.56) a questo punto il processo allargato è stabilizzabile se e soltanto se: A − α I B k11 B rank = n + 1 ∀ α ∈ C+ − λ1 C − α 0 (A.57) 129 strumenti necessari per la sintesi dei regolatori come noto, la condizione sopra indicata va verificata solamente per i valori di α che sono autovalori dell’impianto allargato aventi parte reale non negativa, ossia in questo caso, per α = 0 e per tutti gli α ∈ C+ che sono autovalori di A. Perciò per α = 0 la condizione (A.57) diviene ha: A B k11 B rank = n+1 (A.58) − λ1 C 0 0 ma sotto l’ipotesi che: A B rank = n+1 − λ1 C 0 (A.59) ebbene si vede subito che tale condizione è certamente verificata in quanto: A B k11 B A B rank = rank = n+1 (A.60) − λ1 C 0 0 −C 0 Resta dunque da mostrare che la condizione di cui sopra sia verificata anche per tutti gli α 6= 0 ∈ C+ che sono autovalori di A. Ma ciò è immediato da constatare in quanto data l’ipotesi di stabilizzabilità della coppia ( A, B) tutte le righe della matrice: A − α I B B k11 (A.61) rank − λ1 C − α 0 sono linearmente indipendenti gli α 6= 0 ∈ C+ che sono autovalori di A. Infatti le prime n sono linearmente indipendenti, per l’ipotesi di stabilizzabilità della coppia ( A, B) e l’ultima riga di tale matrice è linearmente indipendente da quelle soprastanti dato che l’elemento di posto (n + 1 , n + 2) è nullo, dimostrando di fatto che il rango della matrice considerata sia proprio n + 1. Rimane dunque da dimostrare che l’impianto allargato sia rilevabile dall’errore di regolazione ỹ. Ossia bisogna verificare che: A − α I B k11 rank −λ1 C −α = n + 1 ∀ α ∈ C+ −C 0 (A.62) La verifica di tale proprietà va condotta per i valori di α indicati nel caso precedente. In particolare per α = 0 si ha che la condizione (A.62) è banalmente verificata sotto le ipotesi che: A B rank = n+1 (A.63) − λ1 C 0 in quanto: A B k11 A B rank −λ1 C = n+1 0 = rank −C 0 −C 0 (A.64) Mentre ragionando in maniera duale, rispetto al caso della proprietà di stabilizzabilità analizzata in precedenza, si giunge anche in questo caso, sotto l’ipotesi di rilevabilità della coppia ( A, C ) alla conclusione che il processo allargato sia 130 strumenti necessari per la sintesi dei regolatori rilevabile, pertanto è possibile sintetizzare un controllore stabilizzante in retroazione da ỹ, ottenendo in tal modo la stabilità esponenziale locale dell’equilibrio di regolazione. Per concludere la dimostrazione resta da provare che il controllo u(t) sia limitato, ma ciò è immediato in quanto tale controllo per struttura è combinazione lineare di variabili di stato che sono limitate per condizioni iniziali contenute nel bacino di attrazione del punto di equilibrio in questione, concludendo di fatto la dimostrazione dell’enunciato. Prima di concludere la dimostrazione sembra opportuno ribadire che le condizioni di rilevabilità e stabilizzabilità in questo caso sono solamente sufficienti per la soluzione del problema in quanto, nulla si può dire se la matrice dinamica del sistema di controllo contiene uno o più autovalori a parte reale nulla. 131 BIBLIOGRAFIA [1] Francesco Ferrante. “I regolatori standard PID: Il punto di vista delle teoria del controllo”. In: (2011). [2] Hassan Khalil. Nonlinear Systems, Third Edition. Prentice-Hall, 2002. [3] Riccardo Marino. Appunti del corso di Controllo dei Processi. 2011. 132