Un sentito ringraziamento per l’assistenza prestata ... Un sentito ringraziamento per...

1

Un sentito ringraziamento per l’assistenza

prestata durante lo svolgimento del seguente

lavoro a:

• Mrs Jacqueline Wilkie,

ICC, Strathclyde University, Glasgow.

• Proff. Mario Innocenti, Aldo Balestrino, Alberto Landi,

Ing. Francesco Nasuti,

DSEA, Ingegneria, Università di Pisa.

2

1. INTRODUZIONE..................................................................................................................................... 3

2. MODELLO NON LINEARE................................................................................................................... 5

2.1 MODELLO ............................................................................................................................................ 72.2 ATTUATORI.......................................................................................................................................... 82.3 RIDUZIONE A 5 INGRESSI ..................................................................................................................... 92.4 DISACCOPPIAMENTO.......................................................................................................................... 11

2.4.1 Piano orizzontale. ..................................................................................................................... 132.4.2 Piano verticale .......................................................................................................................... 162.4.3 Matrice di disaccoppiamento.................................................................................................... 19

3. MODELLO LINEARE .......................................................................................................................... 22

3.1 STATI................................................................................................................................................. 223.2 CONTROLLABILITÀ ED OSSERVABILITÀ ............................................................................................. 233.3 AUTOSTRUTTURA .............................................................................................................................. 263.4 VALORI SINGOLARI ........................................................................................................................... 343.5 RIDUZIONE DEL MODELLO ................................................................................................................ 38

4. CONTROLLO H∞∞∞∞ ................................................................................................................................. 45

4.1 FORMULAZIONE STANDARD DEL PROBLEMA..................................................................................... 454.2 SELEZIONE DELLE FUNZIONI DI PESO ................................................................................................. 494.3 OTTIMIZZAZIONE DELLE FUNZIONI DI PESO........................................................................................ 53

4.3.1 Funzione “hidx”. ...................................................................................................................... 564.3.2 Funzione “hopt”. ...................................................................................................................... 71

5. RISULTATI ............................................................................................................................................ 76

5.1 SENSITIVITÀ E CONTROLLO DELLA SENSITIVITÀ. ............................................................................... 765.2 IMPIANTO A CICLO CHIUSO LINEARIZZATO......................................................................................... 775.3 IMPIANTO NON LINEARE CONTROLLATO ............................................................................................ 80

6. CONTROLLO LQG............................................................................................................................... 90

6.1 TEORIA .............................................................................................................................................. 906.2 ANALOGIE CON IL PROBLEMA H∞ ...................................................................................................... 926.3 DIFFERENZE....................................................................................................................................... 946.4 RISULTATI ......................................................................................................................................... 94

6.4.1 Sensitività e sensitività di controllo .......................................................................................... 956.4.2 Impianto a ciclo chiuso linearizzato ......................................................................................... 96

6.5 IMPIANTO NON LINEARE CONTROLLATO ............................................................................................ 98

7. CONTROLLO ADATTIVO ................................................................................................................ 102

7.1 IDENTIFICAZIONE :........................................................................................................................... 1037.2 ALGORITMI DI IDENTIFICAZIONE MIMO : ....................................................................................... 104

7.2.1 Funzione arx3 ......................................................................................................................... 1107.2.2 Funzione rarx3........................................................................................................................ 112

7.3 CONFRONTO .................................................................................................................................... 1167.4 FUNZIONI PER L’ADATTAMENTO...................................................................................................... 1207.5 RISULTATI ....................................................................................................................................... 125

8. CONCLUSIONI.................................................................................................................................... 127

3

1. Introduzione

I veicoli subacquei a rimorchio (T.U.V.), correntemente utilizzati sono per lo più passivi.

A causa di errori nella rotta della nave madre, correnti, e vari fattori di disturbo non

prevedibili, desta interesse la possibilità di munirli di un controllo attivo.

Un lavoro per lo sviluppo di un modello matematico non lineare per la simulazione

dinamica di un tale veicolo, è stato svolto in [i].

Lo scopo della presente tesi è sviluppare un sistema di controllo robusto ed affidabile che

consente a un T.U.V. stabilità dell’assetto di navigazione ed un accurato inseguimento delle

traiettorie di riferimento.

Il problema non è semplice per vari motivi:

1) Il comportamento non lineare sia del veicolo che del cavo di traino.

2) La natura multivariabile del sistema, che si manifesta con un alto livello di interazione

tra i canali di comando.

3) L’alto livello di incertezza dovuto ad approssimazioni nel calcolo delle forze

idrodinamiche , a variazioni dei parametri, e a fattori di disturbo non conosciuti

introdotti nella linearizzazione del modello non lineare.

L’obiettivo prefissato è stato comunque raggiunto utilizzando il metodo H∞,che insieme

alla µ-sinthesys è diventato il metodo standard per affrontare problemi di controllo nei

sistemi multivariabili.

Per trovare un controllore H∞, per il sistema, è stato studiato il problema generale di

selezione delle funzioni di peso, ed è stato ricavato un algoritmo in codice Matlab per la

loro ottimizzazione.

Tale algoritmo minimizza una funzione costo che include, per ciascun canale, sia termini

appartenenti al dominio della frequenza (sensitività, frequenza di taglio, ecc), sia termini

appartenenti al dominio del tempo (settling time, overshoot, ecc).

4

Un controllore H∞ è stato quindi trovato con l’aiuto dell’algoritmo sopra detto, ed è stata

verificata la sua validità sull’impianto non-lineare, definendo il range del campo operativo

in cui, essendo completamente affidabile, può essere utilizzato.

Allo scopo di confrontare il controllo H∞ con un controllo ottimo standard, sono state

messe in luce le analogie e le differenze tra la metodologia H∞ e quella LQG, quindi un

controllo di tipo LQG è stato realizzato e le sue prestazioni sull’ impianto sono state

analizzate e confrontate con le prestazioni del controllore H∞ precedentemente ottenuto.

Infine, al fine di realizzare un controllo H∞ di tipo adattivo, il problema dell’ identificazione

ricorsiva di impianti multivariabili è stato preso in considerazione, realizzando due diversi

algoritmi di natura parallela per l’identificazione multivariabile.

Tali algoritmi, pur essendo confrontabili con quelli presenti nel System Identification

Toolbox, non hanno tuttavia dato risultati soddisfacenti ai fini del controllo adattivo, percui

sono stati indicati altri possibili metodi per realizzare tale tipo di controllo.

5

2. Modello Non lineare

Il modello non-lineare utilizzato è il risultato del lavoro fatto in [i] ed è riferito a un veicolo

con due gruppi di pinne trasversali, ciascuna delle quali può muoversi indipendentemente

dalle altre.( vedi Fig. 2-1 ).

Fig. 2-1

6

• VEICOLO

Il veicolo viene considerato come corpo rigido tridimensionale.Il verso di rotazione di

ogni pinna si assume positivo quando è in senso antiorario lungo il proprio asse,ed il

relativo angolo viene indicato con δ; con δw1, per esempio, si indica l’angolo

dell’pinna relativo all’asse w1.

Indicheremo con “B” la terna di riferimento solidale al corpo.

• CAVO

Il cavo è modellato con un numero variabile di segmenti inestensibili collegati da

giunti sferici; la massa di ogni segmento si assume essere alla fine del segmento stesso,

nell’ estremità più vicina al veicolo.

Tutte le forze relative al segmento si conviene siano applicate a questa massa.

È stato considerato un cavo di 100mt diviso in 5 segmenti.

• NAVE

Il punto di attacco della linea di rimorchio alla nave si assume essere nella posizione

[0,0,0] della terna inerziale “E”, fissata a terra; si assume inoltre costante la velocità

della nave: “v-ship”.

7

2.1 Modello

Il modello può essere scritto in questa forma (con riferimento alle coordinate E):

& ( ) ( , , , , )u M q F q u v vc= δ Eq 2-1

dove:

• u è il vettore contenente la velocità angolare del veicolo (3 numeri) e la velocità di

ogni elemento del cavo, riferita quest’ultima al relativo segmento. Da notare che

vista l’inestensibilità dei segmenti, saranno sufficienti soltanto due numeri per

indicare la velocità di ogni elemento ( vedi [i] ).

• q è il vettore contenente l’assetto del veicolo (3 numeri) e, l’ orientamento di ogni

elemento del cavo rispetto alla terna E. Anche qui, sono sufficienti soltanto due

numeri per indicare la posizione di ogni elemento ( vedi [i] ).

• δδδδ è il vettore contenente la posizione di ogni pinna.

• v è la velocità della nave.

• vc è la velocità della corrente.

• M(q) è la matrice contenente i termini inerziali del sistema.

• F è il vettore delle forze che operano sul sistema.

Questo modello non-lineare viene scritto come Matlab s-function.

• ingressi :

a) δδδδ (vettore 8 per 1)

b) velocità e accelerazione della nave (vettore 6 per 1).

• stato :

Contiene lo stato del sistema, formato da entrambi i vettori q ed u descritti

sopra, è quindi, un vettore di 2*(3+2*N) elementi dove N è il numero dei

segmenti del cavo.

Poichè il cavo è diviso in 5 segmenti, il modello ha 26 stati.

Lo stato iniziale viene automaticamente calcolato a partire dal vettore q nel

punto di equilibrio desiderato.

8

Il punto di equilibrio si trova memorizzato in un file.Sono disponibili le funzioni

per il calcolo del punto di equilibrio e per l’avvio della s-function.

• uscite:

a) posizione del veicolo rispetto alla terna E.

b) assetto del veicolo rispetto alla terna E.

c) velocità e accelerazione del veicolo rispetto alla terna B.

d) velocità angolare ed accelerazione angolare del veicolo rispetto alla terna B.

e) posizione di ogni elemento del cavo.

f) ampiezza, angolo di attacco e angolo di slittamento della forza trainante

agente sul veicolo, rispetto alla terna B.

La posizione del veicolo rispetto alla terna E è stata considerata misurabile nelle

sue 3 componenti.

2.2 Attuatori

Ogni attuatore è modellato come in Fig. 2-2:

Fcn = 0.5 − ( ((u[2] <= (−0.4363))*(u[1] < 0)) + ((u[2] >= 0.4363)*(u[1] > 0)) )

1

out_1

1

in_1

f(u)

Fcn

Mux

Mux 0

Constant SwitchSaturation

86.0285

s+20Transfer Fcn

1/s

Integrator

−+

Sum

Fig. 2-2

Il primo blocco di saturazione limita la velocità angolare (±10 rad/sec).

Il secondo gruppo, composto dalla funzione e dall’interruttore, limita l’angolo della relativa

pinna (±0.43 rad).

Se angolo e velocità si trovano entro questi limiti, il guadagno che si ottiene è molto vicino

ad uno.

9

2.3 Riduzione a 5 ingressi

Nel modello base ogni pinna può muoversi indipendentemente dalle altre e pertanto si hanno

8 superfici manovrabili.

Tuttavia, per ottenere un parziale disaccoppiamento del sistema, e per ridurre il numero di

attuatori è stato deciso di considerare accoppiate tutte le pinne contrapposte, eccetto la

coppia orizzontale anteriore, come descritto in seguito.

Il blocco Simulink chiamato “5_2_8” fornisce in uscita l’angolo di ogni pinna (δw1 ... δt4 )

prendendo in ingresso solamente cinque numeri (δ1..δ5):

1. δw1 = -δ1 -δ4

2. δw2 = -δ1 +δ4

3. δw3 = δ2

4. δw4 = -δ2

5. δt1 = -δ5

6. δt2 = δ5

7. δt3 = δ3

8. δt4 = -δ3

Da notare che in questo modo diventa possibile utilizzare per le ultime 6 pinne soltanto 3

attuatori invece di 6.

Per esempio, lo stesso attuatore può fornire gli angoli δ3 e -δ3 per la coppia verticale di

pinne δt3 e δt4.

Quindi il numero totale di attuatori può essere ridotto da 8 a 5.

10

Il significato dei cinque comandi è quindi questo:

1. δ1 : comando della rotazione attorno all’asse x.

2. δ2 : comando delle pinne verticali anteriori.

3. δ3 : comando delle pinne verticali posteriori.

4. δ4 : comando delle pinne orizzontali anteriori.

5. δ5 : comando delle pinne orizzontali posteriori.

Da notare che i comandi δ2 e δ3 causano un movimento nel piano orizzontale (deriva lungo

y ed imbardata), ed i comandi δ4 e δ5 causano un movimento nel piano verticale (deriva

lungo z e beccheggio).

Il blocco Simulink chiamato “attuatori” fornisce quindi 5 attuatori necessari per le 8 pinne,

come mostrato in Fig. 2-3

−1

Gain5

Single Actuator angle & speed saturation

Act_4

−1

Gain6


Act_5

−1

Gain4


Act_3

1

in_1

Demux

Demux

1

out_1

Mux

Mux2


Act_2


Act_1

Fig. 2-3

11

2.4 Disaccoppiamento

Per le variabili di posizione ed assetto saranno utilizzate le seguenti notazioni:

• “roll”(rollio): la rotazione ϕ lungo l’asse x nella terna B.

• “sway”(slittamento): lo spostamento Y lungo l’asse y nella terna E.

• “yaw”(beccheggio): la rotazione ψ lungo l’asse z nella terna B.

• “heave”(sollevamento): lo spostamento Z lungo l’asse z nella terna E.

• “pitch”(imbardata): la rotazione θ lungo l’asse y nella

terna B.

È evidente, che avendo un numero n di ingressi, non è possibilie controllare

indipendentemente più di n di uscite, dovremo quindi selezionare tra tutte le uscite

disponibili, le cinque che vogliamo controllare.

Per applicazioni quali il controllo delle tubazioni o dei cavi sul fondo del mare, è molto

importante fornire una buona manovrabilità nell’ assetto (roll, pitch e yaw) e un buon

controllo della posizione, almeno lungo gli assi x ed y della terna E.

A questo punto è necessario verificare se sia veramente possibile, utilizzando i cinque

ingressi descritti sopra, controllare indipendentemente roll, pitch, yaw, sway ed heave.

Questo può essere appurato linearizzando l’impianto con questo insieme di ingressi e uscite

intorno ad un punto di lavoro, e controllando il rango (o il numero di condizione) della

funzione di trasferimento nell’intervallo di frequenza in cui vogliamo controllare il veicolo.

Se il rango non è pieno (o se il numero di condizione, che è il rapporto tra il più alto ed il

più basso dei valori singolari, è molto alto) significa che vicino al punto operativo scelto, il

problema è mal condizionato, vale a dire che non saremo in grado di controllare

indipendentemente tutte le uscite con gli ingressi selezionati, e quindi dovremo considerare

la possibilità di aumentare il numero di ingressi oppure diminuire il numero di uscite da

controllare.

12

Il punto di lavoro prescelto è quello in cui sia la nave che il veicolo si muovono ad una

velocità di 3 m/s nella direzione dell’asse x, della terna E.

Come potremo vedere nel capitolo 3, la funzione di trasferimento in questo punto operativo

ha rango pieno (ed un numero di condizione di 44 dB), questo significa che vicino a questo

punto sarà possibile controllare l’assetto e le due variabili di posizione.

Ad ogni modo, per fornire una visione più ampia del comportamento del veicolo, faremo un

analisi dell’assetto e della posizione, per un certo numero di punti di lavoro, ognuno dei

quali è un punto di equilibrio risultato da una differente combinazione di valori di ingresso

δ2..δ5.

13

2.4.1 Piano orizzontale.

In questo esempio, ogni punto di equilibrio si ottiene con un valore diverso degli ingressi

δ2 e δ3 mentre δ1, δ4 e δ5 sono fissati a zero.

Per tutti questi punti operativi, roll, pitch ed heave rimangono fissati a zero, quindi gli

ingressi δ2 e δ3 incidono soltanto sulla posizione nel piano orizzontale.

Il valore assoluto dello sway è rappresentato in Fig. 2-4

−20

−10

0

10

20

−20−10

010

200

10

20

30

40

delta 2delta 3

abs(

Y)

sway

Fig. 2-4

Si nota facilmente che c’è una linea, (chiamata linea dello“zero-sway”), in cui lo sway è

zero, quindi se noi muovessimo le pinne in modo da spostare il punto di equilibrio su

questa linea, il veicolo sarebbe soggetto ad una rotazione (yaw) senza nessun altro

movimento.

Per lo stesso punto di equilibrio, la Fig. 2-5 mostra il valore assoluto dello yaw.

14

−20−10

010

20

−20

−10

0

10

200

5

10

15

20

delta 2delta 3

abs(

Psi

)

yaw

Fig. 2-5

Anche qui, c’è una linea (linea “zero yaw”), in cui lo yaw è zero, quindi, se muovessimo le

pinne in modo da far rimanere su questa linea il punto di equilibrio, il veicolo sarebbe

soggetto ad uno slittamento (sway) senza nessun altro movimento.

La possibilità di avere lo yaw senza sway o lo sway senza yaw non si sarebbe ottenuta se il

veicolo avesse avuto soltanto una coppia di pinne per il movimento sul piano orizzontale,

(per esempio soltanto la coppia di pinne di coda o soltanto la coppia di pinne anteriori),

invece di due coppie.

Questo fatto suggerisce due cose:

a) se vogliamo controllare sway e yaw indipendentemente, non possiamo unire gli ingressi

δ2 e δ3 per avere un solo ingresso indipendente (per esempio non possiamo utilizzare

solamente δ2 e porre δ3 = K*δ2).

15

b) Possiamo provare a disaccoppiare yaw e sway utilizzando due comandi diversi

direttamente collegati con δ2 e δ3:

δs (sway command) → causa un movimento lungo la linea“zero yaw”,

cioè uno sway senza yaw.

δy (yaw command) → causa un movimento lungo la linea “zero sway”

cioè, uno yaw senza sway.

Questo si può ottenere utilizzando una matrice che trasforma questi nuovi comandi nella

corretta combinazione di δ2 e δ3:

[ ] [ ][ ]δδ

δδ

2

3 1 1=

a b s

yEq 2-2

δ2 = a * δ3 è l’equazione della linea zero yaw

δ2 = b * δ3 è l’equazione della linea zero sway

Per il punto di equilibrio di riferimento il valore dei due coefficenti è:

a = -2.56

b = 2.63

16

2.4.2 Piano verticale

Si può procedere con la stessa analisi per il piano verticale.

In questo esempio, ogni punto di equilibrio viene ottenuto con diversi valori per gli ingressi

δ4 e δ5 mentre δ2, δ3 e δ1 sono fissati a zero.

Per tutti questi punti operativi, roll sway e yaw rimangono fissati a zero, e quindi gli

ingressi δ2 e δ3 incidono soltanto sulla posizione nel piano verticale.

Il valore assoluto dell’ heave è mostrato in Fig. 2-6:

−20−10

010

20

−20

−10

0

10

200

10

20

30

40

50

delta 4delta 5

abs(

Z)

heave

Fig. 2-6

C’è una linea “zero-heave” dove l’heave è zero, e quindi se muoviamo le pinne in modo

tale da spostare il punto di equilibrio su questa linea, il veicolo sarà soggetto ad una

rotazione attorno all’asse y (pitch) senza nessun altro movimento.

17

Per lo stesso punto di equilibrio, la Fig. 2-7 mostra il valore assoluto del pitch.

−20−10

010

20

−20

−10

0

10

200

5

10

15

20

delta 4delta 5

abs(

The

ta)

pitch

Fig. 2-7

Si può vedere una linea “zero-pitch”, dove il pitch è zero, quindi, se muoviamo le pinne in

modo da far rimanere il punto di equilibrio su questa linea, il veicolo sarà soggetto ad uno

spostamento lungo l’asse z (heave) senza alcun altro movimento.

Come già visto nel piano orizzontale, abbiamo che:

a) Se vogliamo controllare indipendentemente gli angoli di heave e di pitch, non possiamo

collegare gli ingressi δ4 e δ5 sotto un unico ingresso indipendente (per esempio

utilizzando soltanto δ4 e ponendo δ5 = K*δ4).

b) Possiamo provare a disaccoppiare pitch ed heave utilizzando due comandi differenti

collegati direttamente con δ4 e δ5:

18

δh (heave command) → causa un movimento lungo la linea “zero

pitch”, cioè, uno heave senza pitch.

δp (pitch command) → causa un movimento lungo la linea “zero

heave” , esattamente, un pitch senza heave.

Questo si può ottenere utilizzando una matrice che trasformi questi nuovi comandi in una

corretta combinazione di δ4 e δ5:

δδ

δδ

4

5 1 1

=

c d h

p

Eq 2-3

δ4 = c * δ5 è l’equazione della linea “zero pitch”

δ4 = d * δ5 è l’equazione della linea “zero heave”

Per il punto di equilibrio di riferimento il valore dei due coefficenti è:

c = -2.56

d = 1.69

19

2.4.3 Matrice di disaccoppiamento

Se mettiamo insieme le due matrici per il disaccoppiamento nei piani orizzontale e

verticale, otteniamo una matrice 5 per 5 che potrà essere utilizzata così:

δδδδδ

δδδδδ

1

2

3

4

5

1 0 0 0 0

0 0 0

0 1 1 0 0

0 0 0

0 0 0 1 1

=a b

c d

r

s

y

h

p

Eq 2-4

In questo modo, il nuovo vettore di ingressi δr.. δp, viene trasfomato nel desiderato vettore

di comandi δ1..δ5.

Si può dimostrare che la sensitività dei coefficienti (a, b, c, e d) rispetto alle condizioni di

equilibrio è molto ridotta (la variazione è minore del 10% quando il punto di equilibrio

viene scelto entro un raggio di 40 m2 intorno al punto di equilibrio di riferimento).

Questo significa che l’orientamento di tutte le linee prima descritte, rimane quasi lo stesso

per un raggio di punti di lavoro ragionevolmente ampio, e quindi, probabilmente saremo in

grado di utlizzare controllori molto simili (o forse lo stesso) per punti operativi differenti.

20

La Fig. 2-8 mostra la posizione del punto di equilibrio ottenuto lavorando sul nuovo

insieme di ingressi.

−100

10 −10

0

10−50

0

50

delta sdelta y

Y

sway

−100

10 −10

0

10−50

0

50

delta h delta p

Z

heave

−100

10 −10

0

10−20

0

20

delta h delta p

The

tapitch

−10

0

10 −10

0

10

−20

0

20

delta s delta y

Psi

yaw

Fig. 2-8

Come si può notare, ogni comando muove il punto di equilibrio su un solo asse, cioè il

veicolo si muove soltanto in una direzione senza alcuna rotazione (o ruota intorno ad un

solo asse senza alcun altro movimento).

Si può mostrare anche che l’ingresso δ1 incide soltanto sul roll.

La matrice di disaccoppiamento è inclusa nel blocco simulink DM.

L’impianto preceduto da questa matrice si presta ad un controllo fatto attraverso metodi più

“classici”: per esempio, si potrebbe ricercare la possibilità di 5 diversi controllori SISO

(single-input single-output), uno per ogni canale.

21

Lo scopo principale di questa sezione, ad ogni modo, è quello di stabilire gli insiemi di

ingresso ed uscita che saranno utilizzati, piuttosto che provvedere ad un vero

disaccoppiamento del sistema; infatti se il problema non è mal condizionato, il regolatore è

in grado di disaccoppiare il sistema indipendentemente dalla presenza della matrice di

disaccoppiamento.

Effettivamente, sono stati trovati controllori per entrambi gli impianti, con e senza matrice

di disaccoppiamento: le prestazioni sono risultate uguali eccetto alcune differenze molto

piccole imputabili a questioni di tipo numerico.

L’impianto che, quindi, andremo a considerare sarà quello descritto in Fig. 2-9.

2

d23

d34

d4

1

d1

5

d5

Mux

Mux

Demux

Demux1

R2D

rad2deg

Demux

Demux2

2

y4

z1

fi

3

psi

5

th uequ(9:14)

Ship vel & acc

Mux

yequ

Eq. values

CABLE − VEHICLESYSTEM

−

+

Sum

Actuators 5_2_8

Demixing

D2R

deg2rad

Fig. 2-9

I blocchi “deg2rad” e “rad2deg” provvedono alle trasformazioni da gradi in radianti e

viceversa, per cui gli ingressi del sistema vengono misurati in gradi, e le uscite in gradi e

metri.

Il blocco “costante” yequ e uequ(9:14) fornisce le condizioni iniziali per il sistema.

Le variabili uequ e yequ vengono caricate nello spazio di lavoro dalla funzione “nli.m”,

che deve essere lanciata prima di utilizzare lo schema (vedere l’ help della funzione per

ulteriori dettagli).

22

3. Modello Lineare

Il modello lineare viene calcolato utilizzando la funzione Matlab “linmod”, per l’estrazione

di un modello lineare da uno schema simulink o da una s-function.

Come risulato avremo un modello lineare instabile con 5 ingressi, 5 uscite, 32 stati, le cui

caratteristiche saranno descritte in dettaglio nella sezione che segue.

3.1 Stati

In base al modo in cui il simulink ordina gli stati delle s-functions riferite ai vari blocchi in

uno schema, il significato fisico degli stati è :

stati 1 .. 5 : stati degli attuatori: derivano dal blocco integratore di ogni attuatore,

nell’ordine in cui appaiono nella Fig. 2-3.

stati 6 .. 8 : roll, pitch e yaw del veicolo.

stati 9 .. 11 : componenti della velocità angolare del veicolo, lungo gli assi x,y,x, di B.

stati 12 .. 13 : θ15 e θ25, orientamento dell’ultimo punto del quinto segmento del cavo (il

punto di aggancio tra il veicolo e il cavo), rispetto alla terna E (vedere [i]

per dettagli).

stati 14 .. 15 : u15 e u25, velocità del punto suddetto rispetto alla terna di riferimento del

segmento (vedere [i] per dettagli).

stati 16 .. 17 : θ14 e θ24, orientamento dell’ultimo punto del quarto segmento del cavo

rispetto alla terna E.


segmento.

23

stati 20 .. 21 : θ13 e θ23, orientamento dell’ultimo punto del terzo segmento del cavo



segmento.

stati 24 .. 25 : θ12 e θ22, orientamento dell’ultimo punto del secondo segmento del cavo


stati 26 .. 27 : u12 e u22, velocità del punto suddetto rispetto al riferimento del segmento.

stati 28 .. 29 : θ11 e θ21, orientamento dell’ultimo punto del primo segmento del cavo


stati 30 .. 31 : u11 e u21, velocità del punto suddetto rispetto alla terna E.

stati 32 .. 36 : stati degli attuatori, più precisamente dei blocchi “transfer-fcn” di ogni

attuatore, nell’ordine in cui appaiono nella Fig. 2-3.

3.2 Controllabilità ed osservabilità

I grahmiani di controllabilità ed osservabilità sono le matrici Gc e Go tali che:

0

0

=++

=++

CCAGGA

BBAGAGT

OOT

TTCC

Eq. 3-1

Indichiamo con Sc e So le matrici diagonali che contengono lungo la diagonale i valori

singolari dei grahmiani di controllabilità e osservabilità, e con Uc e Uo le matrici le cui

colonne sono i vettori singolari di questi grahmiani.

24

Da notare che essendo i grahmiani simmetrici, i vettori singolari di ingresso sono uguali a

quelli di uscita, ( e sono uguali agli autovettori ) percui si puo’ parlare di vettore singolare i-

esimo senza generare confusione.

Il grahmiano di controllabilità descrive la regione dello spazio degli stati raggiungibile con

un ingresso di una data energia, e con condizioni iniziali nulle.

Il grahmiano di osservabilità descrive la regione di punti nello spazio degli stati che, se

presi come condizioni iniziali, provocano un uscita di una data energia, con un ingresso

nullo.

Queste regioni sono iper-ellissoidi con semiasse i-esimo dato da √(S(i))*U(i), per la

controllabilità, e (1/√(S(i)))*U(i) per l’osservabilità, dove S(i) e’ l’i-esimo valore singolare

ed U(i) il rispettivo vettore singolare.

Volendo una misura dell’osservabilità, faremo riferimento al suo ellissoide inverso, con

semiassi dati da √(S(i))*U(i).

Le equazioni di questi ellissoidi, negli spazi di coordinate xc= Uc-1 x ed xo= Uo

-1 x sono:

1

11

1

=

=−

−

ooT

o

ccT

c

xSx

xSxEq. 3-2

Per ogni stato, chiameremo controllabilità dello stato stesso il modulo dell’ intersezione tra

la direzione dello stato, e l’ellisse di controllabilità, analogamente chiameremo osservabilità

dello stato il modulo dell’ intersezione tra la direzione dello stato e l’ellissoide inverso di

osservabilità.

Per ogni polo, chiameremo controllabilità del polo stesso il modulo dell’ intersezione tra la

direzione del suo autovettore e l’ellissoide di controllabilità, analogamente chiameremo

osservabilità del polo il modulo dell’ intersezione tra la direzione del suo autovettore , e

l’ellissoide inverso di osservabilità.

25

La Fig. 3-1 mostra la controllabilità e osservabilità degli stati del sistema.

Fig. 3-1

Alcuni stati dell’attuatore (4,5,35,36 dai due attuatori di coda) hanno una controllabilità

molto limitata, ed è interessante vedere che alcuni stati del cavo diventano meno

controllabili ed osservabili quando il segmento relativo si allontana dal veicolo.

0 5 10 15 20 25 30 35 4010

−8

10−7

10−6

10−5

10−4

10−3

10−2

10−1

100

101

102

state

ctrb

(re

d), o

bsv

(blu

e)ctrb & obsv gram. sv

26

3.3 Autostruttura

Il sistema ha 32 poli e 30 zeri come mostrato nella Fig. 3-2.

−14 −12 −10 −8 −6 −4 −2 0 2−2.5

−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

2.5

Real Axis

Imag

Axi

s

system poles & zeros

Fig. 3-2

Come possiamo vedere, ci sono molti poli quasi sull’ asse immaginario, si vedrà in seguito

che questi derivano dal cavo, ed alcuni di essi sono scarsamente controllabili, o osservabili,

come si può dedurre dal fatto che sono coperti da zeri.

27

In dettaglio, i poli del sistema sono:

-10.78267997416636

-12.20943288835049

-12.81857043722292

-12.48110950043442

-8.88568277132924

-0.00000598044715 + 2.43005154563980i

-0.00000598044715 - 2.43005154563980i

-1.35719650928428 + 1.68805104150768i

-1.35719650928428 - 1.68805104150768i

-0.00001694702355 + 1.95889578230110i

-0.00001694702355 - 1.95889578230110i

-1.93370877010762

0.00006274829235 + 1.41675327418150i

0.00006274829235 - 1.41675327418150i

-1.21489866154187

0.00060264088379 + 0.74644629911264i

0.00060264088379 - 0.74644629911264i

-0.70866401731913

28

-0.62541170112530

-0.49805540421703 + 0.13924393897637i

-0.49805540421703 - 0.13924393897637i

-0.31816541278695

-0.12121713292313

-0.04119107837663

-0.03183408600548 + 0.05569078955429i

-0.03183408600548 - 0.05569078955429i

-6.26215302613939

-6.26215302613923

-6.26215302613926

-13.73784697386058

-13.73784697386077

-13.73784697386074

-13.73784697386074

-13.73784697386075

-6.26215302613925

-6.26215302613923

29

mentre gli zeri di trasmissione sono:

1.348921368961114e+16 - 4.020845700183954e+14i

-7.619196471470863e+15 + 5.873919545104043e+13i

2.004570579182056e+14 + 1.345073802265565e+12i

-1.884457552285398e+14 + 2.250560231041280e+12i

-1.664662554998474e+15 + 4.168019134461320e+13i

-4.504071436680867e+14 - 3.597053425848050e+13i

1.849889636566270e+05 - 2.486281833308713e+04i

-1.846195790399146e+05 + 2.479175889954268e+04i

-7.443550783881528e+04 - 7.286231184565877e+04i

7.486074650008259e+04 + 7.284723447664162e+04i

7.306837910591389e+04 - 7.465468171851218e+04i

-7.264314122837997e+04 + 7.463961612719821e+04i

-2.464227403396731e+04 - 1.848377944059751e+05i

2.501163846042297e+04 + 1.847667370775895e+05i

-12.81835142914654 - 0.00000000000000i

-12.48074869583046 - 0.00000000000000i

-12.20887247085085 + 0.00000000000004i

-10.75071561219304 - 0.00000000000002i

-0.00000508190568 - 2.43004624109670i

-0.00000508190568 + 2.43004624109669i

30

-0.00000508190568 - 1.95875225000447i

-0.00000508190569 + 1.95875225000446i

-0.00000508190567 - 1.41604232476610i

-0.00000508190569 + 1.41604232476610i

-0.00000508190569 - 0.74355371345502i

-0.00000508190569 + 0.74355371345502i

-0.69953819552865 + 0.00000000000000i

-0.48317205044712 + 0.00000000000000i

-0.08852299104955 - 0.00000000000000i

-0.26240694205817 + 0.00000000000000i

I primi 14 zeri derivano probabilmente da problemi numerici.

Il sistema è instabile in quanto ci sono due coppie di poli complessi-coniugati appena oltre

l’asse immaginario.

Gli ultimi dieci poli sono quelli introdotti dai cinque attuatori, infatti, come abbiamo visto

nel secondo capitolo linearizzando il modello dell’attuatore, ogni attuatore ha una dinamica

del secondo ordine con i poli in -13.73 and -6.26.

Avendo a disposizione alcune informazioni sul significato fisico degli stati, è interessante

utilizzarle per vedere in che modo i poli del nostro sistema sono correlati agli stati.

Questo potrebbe dirci di più sul significato fisico dei modi del sistema.

31

I grafici da Fig. 3-3 sino alla Fig. 3-14 mostrano, per ogni polo, un diagramma a barre

rappresentante il valore assoluto delle componenti dell’ autovettore relativo al polo, lungo

le direzioni degli stati, o in altre parole, il modo in cui ogni polo viene ripartito tra gli stati.

Si mostrano anche la controllabilità e osservabilità del polo.

Per ciascuna coppia di autovalori complessi coniugati si mostra soltanto un grafico.

Le Fig. 3-3, 3-4 mostrano le due coppie di autovalori instabili.

Fig. 3-3 Fig. 3-4

Fig. 3-5 Fig. 3-6

È evidente che entrambi modi sono correlati agli stati del cavo.

I due modi rappresentati nelle Fig 3-5, 3-6 sono ugualmente correlati agli stati del cavo.

0 5 10 15 20 25 30 35 400

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

states

abs(

eige

nvec

tor)

Pole : 0.0006026+0.7464i ctrb : 4.596 obsv : 0.0001904

0 5 10 15 20 25 30 35 400

0.1

0.2

0.3

0.4

0.5

0.6

states

abs(

eige

nvec

tor)

Pole : 0.00006275+1.417i ctrb : 15.68 obsv : 0.0006203

0 5 10 15 20 25 30 35 400

0.1

0.2

0.3

0.4

0.5

0.6

0.7

states

abs(

eige

nvec

tor)

Pole : −5.98e−006+2.43i ctrb : 7.211 obsv : 0.01292

0 5 10 15 20 25 30 35 400

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

states

abs(

eige

nvec

tor)

Pole : −0.00001695+1.959i ctrb : 20.06 obsv : 0.002055

32

I quattro modi rappresentati. da Fig 3-7 a Fig 3-10 invece sono chiaramente correlati al

veicolo, e sono molto più osservabili che controllabili.

Fig. 3-7 Fig. 3-8

Fig. 3-9 Fig. 3-10

Da notare nella Fig. 3-9 la forte influenza degli stati 6 e 9, cioè del roll e della componente

lungo l’asse x della velocità angolare del veicolo.

0 5 10 15 20 25 30 35 400

0.1

0.2

0.3

0.4

0.5

0.6

states

abs(

eige

nvec

tor)

Pole : −1.357+1.688i ctrb : 0.05416 obsv : 4.05e−006

0 5 10 15 20 25 30 35 400

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

states

abs(

eige

nvec

tor)

Pole : −0.4981+0.1392i ctrb : 0.004519 obsv : 3.344e−007

0 5 10 15 20 25 30 35 400

0.1

0.2

0.3

0.4

0.5

0.6

0.7

states

abs(

eige

nvec

tor)

Pole : −1.934 ctrb : 0.05976 obsv : 1.072e−006

0 5 10 15 20 25 30 35 400

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

states

abs(

eige

nvec

tor)

Pole : −1.215 ctrb : 0.02857 obsv : 5.668e−006

33

Le Fig. 3-11 e 3-12 mostrano due modi correlati agli attuatori, il secondo sembra avere più

influenza sul cavo e sul veicolo rispetto al primo.

Fig. 3-11 Fig. 3-12

Le Fig 3-13, 3-14 mostrano due modi correlati al veicolo e alla parte del cavo vicino al

veicolo.

Fig. 3-13 Fig. 3-14

Tutti i poli degli attuatori (quelli in -6.26 e in -13.73), hanno dei grafici abbastanza simili.

Si potrebbe anche mostrare che i gruppi di poli in -10 e -12 rappresentano modi

scarsamente osservabili, correlati al cavo.

Tutti gli altri modi, sono molto più osservabili che controllabili, e, paiono coinvolgere sia il

cavo che il veicolo, ma è difficile ottenere un’interpretazione più dettagliata.

0 5 10 15 20 25 30 35 400

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

states

abs(

eige

nvec

tor)

Pole : −6.262 ctrb : 0.004594 obsv : 4.722e−007

0 5 10 15 20 25 30 35 400

0.1

0.2

0.3

0.4

0.5

0.6

0.7

states

abs(

eige

nvec

tor)

Pole : −13.74 ctrb : 0.00002181 obsv : 4.321e−007

0 5 10 15 20 25 30 35 400

0.1

0.2

0.3

0.4

0.5

0.6

0.7

states

abs(

eige

nvec

tor)

Pole : −0.04119 ctrb : 0.04383 obsv : 1.401e−006

0 5 10 15 20 25 30 35 400

0.1

0.2

0.3

0.4

0.5

0.6

0.7

states

abs(

eige

nvec

tor)

Pole : −0.03183+0.05569i ctrb : 2.394 obsv : 0.0006888

34

3.4 Valori Singolari

Ad una data frequenza il massimo valore singolare è il guadagno maggiore possibile del

sistema, ed il minimo valore singolare è il minore guadagno possibile [ii].

In dettaglio, i valori singolari della funzione di trasferimento rappresentano il guadagno del

sistema lungo le proprie direzioni principali e cioè, se il segnale di ingresso è nella

direzione di un vettore singolare di ingresso, il segnale di uscita sarà nella direzione del

relativo vettore singolare di uscita e il guadagno del sistema sarà il relativo valore

singolare.

La Fig. 3-15 mostra il comportamento in frequenza dei valori singolari del sistema.

10−3

10−2

10−1

100

101

102

−150

−100

−50

0

50

Frequency (rad/sec)

Sin

gula

r V

alue

s dB

system singular values

Fig. 3-15

35

La banda passante del sistema è di circa 3 rad/sec, e il guadagno (in ogni direzione)

diminuisce velocemente dopo i 10 rad/sec.

Esistono due direzioni con guadagno alto, ed altre tre con guadagno di circa 0 dB.

Vedremo adesso con un’ analisi dei vettori singolari, quali ingressi e uscite sono coinvolti

in questi guadagni principali.

Dalla Fig 3-16 alla 3-27 sono mostrati, per ciascun valore singolare, due grafici a barre

rappresentanti il valore assoluto delle componenti dei vettori singolari corrispondenti al

valore singolare.

La prima analisi viene fatta per una frequenza di 0.01 rad/sec, e ci mostra una esatta

corrispondenza tra direzioni principali e “canali” del sistema.

0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

output

abs(

s. v

ec.)

Singular Value : 178

0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

input

abs(

s. v

ec.)

0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

output

abs(

s. v

ec.)

Singular Value : 44.4

0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

input

abs(

s. v

ec.)

Fig. 3-16 Fig. 3-17

A questa frequenza, gli ingressi δ2 e δ4 (coppia di pinne verticali e orizzontali anteriori),

sono fortemente correlati con sway ed heave (misurati in metri), ed in questi canali ci sono i

due guadagni più alti del sistema (vedi Fig 3-16, 3-17).

36

Le Fig da 3-18 a 3-21 ci mostrano, come già noto,che l’ingresso 1 è fortemente correlato al

roll e che gli ingressi 3 e 5 (coppia di pinne verticali e orizzontali posteriori) sono correlati

alle angolazioni di imbardata e di beccheggio (misurate in gradi).

In questi canali il guadagno è quasi 1.

0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

output

abs(

s. v

ec.)


0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

input

abs(

s. v

ec.)

0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

output

abs(

s. v

ec.)


0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

input

abs(

s. v

ec.)

Fig. 3-18 Fig. 3-19

0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

output

abs(

s. v

ec.)


0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

input

abs(

s. v

ec.)

1 1.5 2 2.5 3 3.5 4 4.5 5−5

−4.5

−4

−3.5

−3

−2.5

−2

−1.5

−1

input

outp

ut

max.sv.(G(jw,i,j))

Fig. 3-20 Fig. 3-21

Dal momento che le direzioni sono fortemente correlate ai canali del sistema, possiamo

dedurre che a questa frequenza la struttura della funzione di trasferimento dovrebbe essere

piuttosto diagonale, con una “dominanza” nella matrice dei termini (2,2) e (4,4).

La Fig 3-20 ci mostra che questo è vero; infatti,è stata calcolata l’ampiezza di tutti i 25

elementi della funzione di trasferimento e le linee in figura rappresentano le curve di livello

della funzione da N2 a R+ così ottenuta.

37

Per frequenze maggiori, la situazione è comunque differente e la stessa analisi per una

frequenza di 0.1 rad/sec mostra che il sistema non è più così diagonale.

0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

output

abs(

s. v

ec.)


0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

input

abs(

s. v

ec.)

0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

output

abs(

s. v

ec.)


0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

inputab

s(s.

vec

.)

Fig. 3-22 Fig. 3-23

0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

output

abs(

s. v

ec.)


0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

input

abs(

s. v

ec.)

0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

output

abs(

s. v

ec.)


0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

input

abs(

s. v

ec.)

Fig. 3-24 Fig. 3-25

38

0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

output

abs(

s. v

ec.)


0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

input

abs(

s. v

ec.)

1 1.5 2 2.5 3 3.5 4 4.5 5−5

−4.5

−4

−3.5

−3

−2.5

−2

−1.5

−1

input

outp

ut

max.sv.(G(jw,i,j))

Fig. 3-26 Fig. 3-27

Come possiamo vedere nei grafici da fig 3-22 a fig 3-27, c’è un forte accoppiamento tra i

comandi per il movimento sullo stesso piano.

3.5 Riduzione del Modello

Il fatto che ci siano molti stati che mostrano una controllabilità o osservabilità molto scarsa,

suggerisce che probabilmente è opportuno provare una riduzione del modello del sistema

per semplificarlo, senza però perdere nessuna delle sue caratteristiche.

Sappiamo infatti che generalmente, tutti gli algoritmi numerici lavorano meglio e più

velocemente al diminuire delle dimensioni delle matrici coinvolte.

Il modo più comune di ridurre un sistema è quello di ottenerne una realizzazione bilanciata

e dopo, procedere ad un troncamento in modo da liberarsi degli stati meno utili. Esistono

diverse funzioni nell’ambiente Matlab che permettono di fare una cosa del genere.

La funzione “balreal” del “control toolbox” serve per ottenere una realizzazione bilanciata,

ma non funziona se la realizzazione non è minima o se ci sono degli stati che sono molto

poco controllabili o osservabili come nel nostro caso.

La funzione “minreal” del “control toolbox”, serve per ottenere una realizzazione minima,

e, applicata al nostro sistema con una tolleranza di default rimuove soltanto uno stato.

39

Il risultato non è ancora tale la permettere a “balreal” di funzionare.

Se la tolleranza in “minreal” va oltre 0.3 allora questa funzione rimuove altri stati lasciando

lavorare balreal, ma il risultato non è ancora soddisfacente.

I risultati migliori si ottengono con le funzioni del µ-tools toolbox.

La funzione “sysbal”, elimina automaticamente gli stati non-controllabili o inosservabili.

Non funziona se il sistema è instabile come nel nostro caso; comunque questa non è

una limitazione, poichè non sarebbe prudente eliminare qualche polo instabile anche se la

sua controllabilità ed osservabilità è scarsa.

Si è quindi deciso di procedere in questo modo:

1. Dividere il sistema in due parti, una instabile e l’altra stabile.

2. Bilanciare e ridurre la parte stabile.

3. Aggiungere la parte instabile al sistema ridotto e bilanciato.

4. Confrontare il sistema risultante con quello originale e quindi ripetere tutto il

procedimento fino ad ottenere un risultato soddisfacente.

Come risultato, è stato ottenuto un sistema di 19 stati che non conviene ulteriormente

ridurre per non perdere alcune caratteristiche importanti come ad esempio il numero di

condizionamento e la forma dei vettori singolari.

40

La Fig. 3-28 mostra la controllabilità e l’osservabilità degli stati del sistema ridotto. Da

notare che ogni informazione sul significato fisico degli stati viene persa perchè vengono

cambiate le coordinate nello spazio degli stati.

Fig. 3-28

0 2 4 6 8 10 12 14 16 18 2010

−1

100

101

102

state

ctrb

(re

d), o

bsv

(blu

e)

ctrb & obsv gram. sv

41

La Fig. 3-29 mostra i poli e gli zeri del sistema.

−5 0 5 10 15 20−2.5

−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

2.5

Real Axis

Imag

Axi

s


Fig. 3-29

Più in dettaglio i poli sono:

0.00006274829235 + 1.41675327418149i

0.00006274829235 - 1.41675327418149i

0.00060264088379 + 0.74644629911264i

0.00060264088379 - 0.74644629911264i

-0.00000573233871 + 2.43005191447068i

42

-0.00000573233871 - 2.43005191447068i

-0.00001700486128 + 1.95889589084189i

-0.00001700486128 - 1.95889589084189i

-0.03175503176227 + 0.05568076358300i

-0.03175503176227 - 0.05568076358300i

-0.04561044561980

1.64169096029987 + 2.17650615654612i

-1.64169096029987 - 2.17650615654612i

-1.05403974832363 + 1.50335622090630i

-1.05403974832363 - 1.50335622090630i

-0.71813672411532 + 0.40721327794409i

-0.71813672411532 - 0.40721327794409i

-0.28119580369960 + 0.14405349626301i

-0.28119580369960 - 0.14405349626301i

Insieme ai due modi instabili, anche i quattro successivi modi appaiono nel sistema

originale, (vedere Fig da 3-5 a 3-14) gli ultimi quattro modi invece sono propri del modello

ridotto e non trovano riscontro nel modello originale.

Gli zeri del sistema sono:

5.611781738910854e+17 - 8.358571600612663e+18i

-1.347886426525800e+14+ 3.591267196515551e+14i

43

1.025405654024770e+16+ 1.036423806824973e+14i

-3.734676705916108e+15 - 5.555069415363998e+15i

-3.282825130228596e+15 + 5.355217793995528e+15i

17.53398985391323 - 0.00000000000004i

8.97968262968278 + 0.00000000000001i

5.54758531982268 - 0.93135381773286i

5.54758531982278 + 0.93135381773283i

5.59680561739371 - 0.00000000000000i

0.00006271452900 + 1.41619016832985i

0.00006271452900 - 1.41619016832985i

0.00003173065577 + 1.95879210243769i

0.00003173065578 - 1.95879210243769i

0.00003203827808 + 0.74369685772168i

0.00003203827808 - 0.74369685772168i

-0.00000312190411 + 2.43004852576467i

-0.00000312190411 - 2.43004852576467i

-0.17889854127239 - 0.00000000000000i

I primi tre zeri sono probabilmente dovuti a problemi numerici.

44

La Fig. 3-30 mostra i valori singolari del sistema.

10−3

10−2

10−1

100

101

102

−80

−60

−40

−20

0

20

40

60

Frequency (rad/sec)

Sin

gula

r V

alue

s dB


Fig. 3-30

Non esiste alcuna differenza apprezzabile rispetto agli andamenti dei valori singolari del

sistema originale non ridotto.

Infine, per una frequenza di 0.1 rad/sec, la “struttura” del sistema, cioè la grandezza di tutti

i 25 elementi della matrice, è sostanzialmente identica a quella dell’ impianto originale (Fig

3-21, 3-27).

45

4. Controllo H∞∞∞∞

Nei capitoli seguenti adotteremo la nomenclatura standard descritta in [iii].(Vedere [iv] per

una dettagliata introduzione all’ argomento, e [v] per la soluzione completa del problema

nello spazio degli stati,).

4.1 Formulazione Standard del Problema

Il diagramma a blocchi utilizzato per la formulazione standard del problema H∞ è riportato

in Fig. 4-1:

P(s)

K(s)

w z

u e

Fig. 4-1

• Il segnale w contiene tutti gli ingressi esterni, inclusi i

disturbi; il segnale z contiene tutte le uscite verso

l’esterno, ed è normalmente utilizzato come un segnale

di errore; e contiene le variabili misurate usate come

ingressi del controllore; u contiene gli ingressi di

controllo, ossia le uscite del controllore.

• P è l’ impianto generalizzato, cioè quello che

generalmente viene chiamato impianto in un problema

di controllo, e comprende anche tutte le funzioni di

peso. Ci riferiremo ad esso con una delle matrici delle

due equazioni seguenti, rispettivamente nel dominio

della frequenza e dello spazio degli stati.

46

z

y

P P

P P

w

u

=

11 12

21 22

&x

z

y

A B B

C D D

C D D

x

w

u

=

1 2

1 11 12

2 21 22

Eq 4-1 Eq 4-2

• K è il controllore da trovare.

Ci si riferisce all impianto di Fig. 4-1 anche come ad una trasformazione lineare frazionaria

(LFT) su K con coefficiente P.

LFT P K T s P s P s K s P s K s P szw( , ) ( ) ( ) ( ) ( )( ( ) ( )) ( )= = + + −11 12 22

1211 Eq. 4-3

La funzione di trasferimento risultante a circuito chiuso da w a z viene chiamata Tzw.

Per una data matrice hamiltoniana M, denoteremo con Ric(M) la soluzione (se esistente)

della relativa equazione Riccati, come in Eq. 4-4, vedi [iv] per maggiori dettagli.

Ric M X

XM M X XM X MT

( ) =

+ + − =

X t.c.

11 11 12 21 0Eq. 4-4

Si chiama “norma-infinito” di una funzione di trasferimento l’estremo superiore, al variare

della frequenza, del suo massimo valore singolare:

T s T ss j

( ) sup( ( ( ))∞=

=ω

σ Eq. 4-5

Tale norma rappresenta in pratica il massimo guadagno possibile del sistema al variare

della direzione e della frequenza del segnale in ingresso.

Un sistema si dice “stabilizzabile” se tutti i suoi poli a parte reale positiva sono

controllabili, si dice “individuabile” se tutti i suoi poli a parte reale positiva sono

osservabili.

47

Infine, chiameremo ρ(N) il raggio spettrale della matrice N, cioè il valore assoluto del suo

autovalore maggiore in modulo.

48

Il problema del controllo H∞ è quello di trovare, tra tutti i controllori che portano ad un

sistema a ciclo chiuso internamente stabile, un controllore K che minimizza la norma

infinito di Tzw

Il risultato matematico basilare per questo problema è espresso sia in [iv] che in [v] e segue

questo schema:

Si suppone che P sia tale che:

1. (A,B2) sia stabilizzabile, e (A,C2) sia individuabile.

2. Rango di D12 uguale al suo numero di colonne,

rango di D21 uguale al suo numero di righe.

Sia γ un numero reale positivo.

Siano H∞ e J∞ due matrici hamiltoniane, funzioni di P e γ-2, costruite come spiegato in [iv]

e in [v].

Esiste un controllore Tzw che è internamente stabile, e la cui norma infinito è minore di γ,

se, e soltanto se, si verificano le tre condizioni che seguono:

1. Esistenza di una soluzione semidefinita positiva di Ric(H∞), chiamata X∞.

2. Esistenza di una soluzione semidefinita positiva di Ric(J∞), chiamata Y∞.

3. ρ(X∞Y∞) < γ2.

Inoltre, quando queste condizioni sono soddisfatte, viene data una formula nello spazio

degli stati per calcolare un controllore che risolve il problema.

In pratica, questo teorema viene utilizzato nella maniera che segue:

Per prima cosa viene selezionato un intervallo di valori di γ in cui vogliamo cercare il

valore minimo della norma infinito di Tzw, dopodichè viene adottato il metodo della

bisezione per iterare sul valore di γ e trovare così il controllore H∞ “ottimo” cioè quello che

minimizza la norma infinito di Tzw.

49

Più dettagliatamente, per ogni γ l’impianto viene scalato per ottenere D11=0, D22=0,

D21D12T=I, D21D12

T=I il che semplifica molto le due matrici hamiltoniane, (senza perdita di

generalità).

Le due matrici hamiltoniane suddette vengono quindi calcolate, e le tre condizioni vengono

testate per stabilire se il controllore esiste o no, dopodichè, il nuovo intervallo di ricerca per

γ sarà delimitato a destra dall’ ultimo valore per cui il controllore esiste ed a sinistra dall’

ultimo valore in cui il controllore non esiste.

Si ripete il procedimento per il γ che divide in due il nuovo intervallo di ricerca

concludendo quando la differenza tra il valore γ ed il precedente è minore della tolleranza

stabilita.

Di solito, quando γ raggiunge il proprio valore ottimale, la soluzione dell’equazione di

Riccati diventa numericamente sensibile per via del fatto che gli autovalori delle matrici

Hamiltoniane corrispondenti si muovono verso l’asse jωdove la soluzione non dà più

controllori stabilizzanti.

L’algoritmo completo sopra descritto viene utilizzato nelle funzioni “hinfsyn” e “hinfsyne”

del “µ tools toolbox” per l’utilizzazione con Matlab [iii].

4.2 Selezione delle funzioni di peso

La Fig. 4-2 sotto mostra il convenzionale ciclo chiuso “ad un grado di libertà” che

adotteremo nel progetto del nostro controllo.

y

++

Sum

r

G(s)

K(s)eu

Fig. 4-2

50

In particolare chiameremo:

• Sensitività: la matrice S(s) tale che e(s)=S(s)r(s), cioè

(I+G(s)K(s))-1.

• Sensitività di controllo:la matrice M(s) tale che u(s)=M(s)r(s), cioè S(s)K(s).

• Sensitività complementare: la matrice T(s) tale che y(s)=T(s)r(s), cioè G(s)M(s).

L’aspetto più cruciale nel modello di controllo H∞ ,è la definizione delle specifiche di

progetto in termini di funzioni di peso.

Normalmente il sistema a circuito chiuso richiesto si deve comportare in maniera molto

diversa a basse o ad alte frequenze[ii].

Alle basse frequenze, dove si assume una conoscenza soddisfacente dell’ impianto, le

necessità più importanti sono l’inseguimento del segnale di riferimento ed una buona

reiezione dei segnali di disturbo.

Queste necessità possono entrambe essere soddisfatte mantenendo la sensitività bassa, visto

che questo significa una buona reiezione dei segnali di disturbo, ed una sensitività

complementare vicina all’ identità, il che implica che il segnale uscita sia vicino a quello di

riferimento.

Al contrario, alle alte frequenze, le necessità più importanti sono il rigetto del rumore dei

sensori ed una stabilità robusta a fronte delle incertezze dovute a dinamiche non modellate,

non linearità, troncamenti del sistema, che di solito sono fattori che diventano rilevanti

all’aumentare della frequenza di lavoro.

Inoltre, generalmente vogliamo che l’energia spesa per il controllo venga ridotta il più

possibile a queste frequenze, visto che non dobbiamo seguire nessun segnale di riferimento.

Si può mostrare facilmente (vedi [ii]) che tutte queste necessità possono essere ottenute

tenendo la sensitività di controllo, e di conseguenza la sensitività complementare più bassa

possibile.

La sensitività di controllo infatti è la funzione trasferimento vista da una incertezza additiva

sull’ impianto, e quindi, per evitare un circuito chiuso instabile tra l’ impianto controllato e

l’incertezza, deve essere ridotta.

51

Inoltre tenere la sensitività di controllo bassa significa tenere basso il segnale di controllo e

la sensitività complementare, avendo così una buona reiezione del rumore dei sensori.

In termini di valori singolari, queste specifiche di progetto significano minimizzare

σmax(S(jω)) alle basse frequenze e σmax(M(jω)) alle alte frequenze.

Se prendiamo

W s M s

W s S s

m

s

( ) ( )

( ) ( )

max

max

=

=

−

−

1

1Eq. 4-6

dove Smax(jω) e Mmax(jω) esprimono in qualche modo la forma desiderata per S ed M,

allora possiamo esprimere le nostre specifiche di progetto in questi termini:

σ ω ωσ ω ω

( ( ) ( ))

( ( ) ( ))

W j S j

W j M js

m

<<1

1Eq. 4-7

Si può mostrare [ii] che se:

W s S s

W s M s

s

m

( ) ( )

( ) ( )

∞< 1 Eq. 4-8

allora le Eq. 4-7 sono soddisfatte su tutte le frequenze.

Quindi, se il controllore è tale da soddisfare la condizione 4-8, le nostre specifiche saranno

rispettate.

Inoltre, tanto più la norma infinito nella 4-8 sarà minore di 1,tanto più saranno ampi i

margini con cui le specifiche verranno rispettate, e quindi migliore sarà il controllore.

52

Ora, se costruiamo l’impianto P in modo tale che:

T sW s S s

W s M swz

s

m( )

( ) ( )

( ) ( )=

Eq. 4-9

Il problema di trovare un controllore tale che le 4-6 siano soddisfatte diventa quello di

trovare un controllore tale che:

[ ]T swz ( )∞

< 1 Eq. 4-10

ossia diventa un problema H∞ standard.

La procedura H∞ infatti riesce a trovare il controllore che minimizza la norma infinito di

Twz(s), se questa norma è anche minore di 1 allora tale controllore è tale da soddisfare la 4-

11 e quindi le 4-6.

Rimane da specificare come costruire l’impianto P in modo da soddisfare la 4-10.

Ciò si può fare facilmente “portando fuori” dall’ impianto i punti tra i quali possiamo

misurare S(s) ed M(s), e moltiplicare le uscite per le funzioni di peso Ws(s) e Wm(s) scelte

come detto sopra.

53

La costruzione dell’ impianto è mostrata in fig. 4-3.

y

++

Sum

reu

Wm(s) Ws(s)

z2 z1

G(s)

Fig. 4-3

4.3 Ottimizzazione delle funzioni di peso.

Anche se controllori H∞ sono facili da ottenere utilizzando pacchetti software esistenti, la

maggior parte dei “controllisti” ha non poche difficoltà nella scelta delle funzioni di peso.

Questo è essenzialmente dovuto al fatto che esistono aspetti familiari dei problemi di

controllo che non vengono completamente catturati dall’approccio nel dominio della

frequenza descritto sopra.

Inoltre, qualche volta, a causa di problemi numerici, la sensitività del controllore finale a

fronte di piccoli cambiamenti delle funzioni di peso, è veramente alta.

Quindi, sarà necessario spendere molto tempo provando varie funzioni di peso per ottenere

un controllore soddisfacente.

Questo è vero specialmente in caso di complessi impianti multivariabili dove le relazioni

tra specifiche di progetto, matrici di sensitività, e funzioni di peso, spesso diventano difficili

da gestire.

L’algoritmo descritto sotto, potrebbe essere uno strumento utile ai progettisti nel lavoro di

selezione delle funzioni di peso.

54

Le principali idee da cui deriva l’ algoritmo sono:

• Se i canali dell’ impianto non sono molto differenti tra

di loro, e richiediamo più o meno lo stesso tipo di

prestazioni da ogni canale possiamo selezionare una

matrice di peso diagonale con la stessa funzione per

ogni canale ( altrimenti possiamo sempre scalare l’

impianto per ottenere la condizione di somiglianza tra i

canali ).In questo modo il problema si riduce alla

selezione di due funzioni scalari di peso, una per S(s) e

una per M(s) o T(s).

• Se utilizziamo una forma standard per queste funzioni

scalari di peso (per esempio usando funzioni di primo

ordine), allora sarà possibile descriverle con pochi

numeri, per esempio la frequenza di taglio e il

guadagno. In questo modo le due matrici delle funzioni

di peso vengono descritte completamente da un vettore

in uno spazio quadrimensionale.

• È possibile costruire una funzione che prenda in

ingresso questo vettore, calcoli le relative matrici delle

funzioni di peso,il controllore H∞, compia varie

simulazioni per testare il controllore ( a questo punto è

possibile anche una simulazione su un impianto non-

lineare), e alla fine come risultato delle simulazioni,

fornisca un indice di qualità del controllore basato sui

termini di dominio temporale e sulle quantità, come ad

esempio il valore finale γ o la posizione dei poli del

controllore.

• È possibile applicare a questa funzione qualche

algoritmo standard di ottimizzazione vincolata per

minimizzare l’indice di qualità e quindi trovare il

55

migliore controllore possibile (semprechè l’indice di

qualità rifletta veramente le specifiche di progetto e le

intenzioni del progettista).

L’algoritmo è scritto completamente in codice Matlab, e nella sezione che segue ne daremo

una spiegazione dettagliata.

56

4.3.1 Funzione “hidx”.

Questa funzione è il cuore del programma di ottimizzazione.La sua sintassi di chiamata è:

[f,g]=hidx(x) dove

• x è il vettore che descrive le matrici delle funzionidi

peso.

• f è l’indice di qualità che deve essere minimizzato nella

regione ammissibile.

• g definisce la regione ammissibile come quella regione

nello spazio delle x in cui ogni elemento di g è minore

di zero.

Come già detto, adotteremo matrici diagonali delle funzioni di peso con la stessa funzione

della frequenza per ogni canale.

Sia:

σ ω ωσ ω ω

( ( )) ( )

( ( )) ( )

M j lm j

S j ls j

<<

Eq.4-11

Per il modo in cui definiremo le funzioni di peso (vedere sezione 4.2), dovremo prima

definire ls(s) e lm(s) come funzioni invertibili di primo ordine, e quindi per ognuna di esse,

avremo bisogno di tre numeri:

1. frequenza del polo.

2. frequenza dello zero.

3. guadagno per s=∞ per s=0.

Per ls la frequenza del polo è fissata a 10^(-6) rad/sec e dovremo fornire la frequenza dello

zero (maggiore di quella del polo) e per il guadagno in s=∞.

I primi due elementi di x si riferiscono ad ls ed hanno questo significato:

x(1) = guadagno per ls(s=inf) in decibel.

x(2) tale che 10^x(2) = frequenza dello zero di ls in rad/sac.

57

Per lm la frequenza dello zero è fissata a 10^(6) rad/sec e dovremo fornire la frequenza del

polo (minore di quella dello zero) e per il guadagno in s=0.

Gli ultimi due elementi di x hanno questo significato:

x(3) = guadagno per lm(s=0) (decibel).

x(4) tale che 10^x(2) = frequenza del polo di ls(rad/sac).

Prima di chiamare questa funzione, le variabili che seguiranno devono essere definite nello

spazio di lavoro come globali; normalmente questa funzione è chiamata da una funzione

più esterna, la quale provvede a ciò, e quindi non sarà necessario preoccuparsi di questo, a

meno che, non vogliamo utilizzare la funzione indipendentemente dall’algoritmo di

ottimizzazione.

Queste variabili sono:

plt, se è diversa da zero abilita l’esposizione automatica di alcuni risultati e deve

essere tenuta a zero durante il procedimento di ottimizzazione per evitare di

rallentarlo eccessivamente.

nl, se è diversa da zero abilita una simulazione sull’impianto non lineare e viene

generalmente mantenuta a zero se tale simulazione necessita di tempi

piuttosto lunghi.

sz, numero di righe e colonne dell’ impianto in A B C D (quadrato).

A B C D, modello lineare del sistema G in

Ak Bk Ck Dk, sono le variabili destinate ad essere riempite dalle matrici del controllore.

a1 b1 c1 d1, sono le variabili destinate ad essere riempite dal sistema che rappresenta la

matrice delle funzioni di peso per la sensitività, cioè Ws(s).

58

a2 b2 c2 d2, sono le variabili destinate ad essere riempite dal sistema che rappresenta la

matrice delle funzioni di peso per la sensitività di controllo, cioè Wm(s).

Le variabili uequ yequ e N, vengono riempite rispettivamente con u e y, nel punto di

equilibrio, e con il numero dei segmenti nel cavo; quindi queste sono variabili che

riguardano il sistema particolare e che non dovrebbero più essere nello spazio di lavoro per

un altro sistema.

Le variabili che seguono sono necessarie soltanto per la procedura di ottimizzazione:

X, È una matrice contenente in ogni colonna il vettore x corrispondente ai

quattro migliori controllori ottenuti (la prima colonna individua il miglior

controllore ottenuto).

F, È un numero che rappresenta il migliore indice

ottenuto { F=max(0,1e5*g(X(:,1)))+f(X(:,1)) }.

CT, È un contatore che indica quante volte questa funzione è stata chiamata

Vediamo il codice con maggiori dettagli:

• Per prima cosa c’è l’inizializzazione e l’incremento del

contatore globale:

global plt nl sz A B C D Ak Bk Ck Dk a1 b1 c1 d1 a2 b2 c2 d2 uequ yequ N X F CT

CT=CT+1,

59

• L’ istruzione:

w=logspace(-6,6,100);

definisce un vettore di 100 punti spaziato logaritmicamente tra 10^-6 e 10^6, che

utilizzeremo come frequenza.

• Le due linee:

[a1,b1,c1,d1]= hfot(-6,x(2),x(1),1,sz);

[a2,b2,c2,d2]= hfot (6,x(4),x(3),1,sz);

servono per la creazione delle matrici delle funzioni di peso a partire dal vettore x

(vedere l’ help nella funzione hfot per maggiori dettagli).

• Le linee:

[aw1,bw1,cw1,dw1]=hfot(-6,x(2),x(1),0,1);

[aw2,bw2,cw2,dw2]=hfot(6,x(4),x(3),0,1);

sw1=max([sigma(aw1,bw1,cw1,dw1,w);-inf*ones(size(w))]);

sw2=max([sigma(aw2,bw2,cw2,dw2,w);-inf*ones(size(w))]);

costruiscono ls(s) e lm(s) e li immagazzinano nei vettori sw1 e sw2.

A cosa serve la riga -inf*ones(size(w)) ?.

Supponiamo che essa non ci sia e che sigma(aw2,bw2,cw2,dw2,w) esprima l’andamento

in frequenza di un solo valore singolare.In tal caso l’istruzione “max”estrarrebbe il

massimo valore in frequenza di questa risposta.Quello che noi vogliamo,invece,è la

risposta in frequenza tutta intera,e per questo motivo è necessaria la presenza della riga -

inf.*ones(size(w)).

60

Adesso che abbiamo sistemi e funzioni di peso, dobbiamo costruire l’ impianto nella forma

generale come mostrato in Fig. 4-3 per la procedura H∞.

Esistono vari metodi in Matlab per l’interconnessione dei sistemi, il µ-tools toolbox

fornisce anche un m-file chiamato “sysic” piuttosto semplice da utilizzare.

Dopo aver fatto alcuni confronti, è stato deciso di interconnettere i blocchi utilizzando

l’interfaccia grafica Simulink, e dopo utilizzare “linmod” per estrarre il modello lineare di

tutto il sistema ottenuto.

Questo viene fatto tramite le istruzioni “[ap,bp,cp,dp]=linmod('hloop');”

dove ap,bp,cp,dp è l’impianto da controllare e hloop è lo schema in Fig. 4-4.

emu

Demux2

6

e67

e78

e89

e910

e10

5

e5

4

e4

3

e3

2

e2

1

e1

emu

Demux1

++

Sum

10

u10

8

u8

Mux

Mux2

9

u9

7

u7

6

u6

15

e15

14

e14

13

e13

12

e12

11

e11

emu

Demux3

x’ = Ax+Bu y = Cx+Du

System

1

u12

u2

5

u5

4

u4

3

u3

Mux

Mux1x’ = Ax+Bu y = Cx+Du

W2(s)


W1(s)

Fig. 4-4

61

Le ragioni di questa scelta tra i vari metodi di interconnessione possibili sono tre:

1) È il metodo più adatto per l’utente ed è facile da gestire.

2) È il metodo più veloce, e questo è importante visto che utilizzeremo questa funzione

diverse volte durante l’ottimizzazione.

3) È stato provato essere il metodo numericamente più affidabile, specialmente per sistemi

complessi di ordine elevato.

Svantaggi possibili di questo metodo sono:

1) Necessità di utilizzare variabili globali.

2) Necessità di cambiare il numero delle porte di input ed output nello schema Simulink, se

vogliamo utilizzare l’algoritmo di ottimizzazione per un altro impianto con un numero

diverso di ingressi o uscite

• L’istruzione:

plant=pck(ap,bp,cp,dp);

Impacchetta l’impianto nella forma ” packed state space” utilizzata in µ-tools.

• Con le istruzioni

gmax=1;

[K,P,gf]=hinfsyne(plant,sz,sz,0,gmax,.01,inf,1);

viene fatto un primo tentativo di trovare la soluzione per γ nell’intervallo 0-1.

62

• Se il tentativo fallisce,l’intervallo di ricerca per γ viene

successivamente esteso fino a che non si trova un

controllore:

while K==[]:

gmax=10*gmax;

[K,P,gf]=hinfsyne(plant,sz,sz,0,gmax,.01,inf,1);

end

• Con l’istruzione:

[Ak,Bk,Ck,Dk]=unpck(K);

il controllore viene poi memorizzato nelle variabili globali Ak,Bk,Ck,Dk.

La sezione che segue fornisce la sensitività dell’ impianto a circuito chiuso ottenuto.


[ap1,bp1,cp1,dp1]=linmod('hloop1');

dove hloop1 è lo schema della Fig. 4-5, serve per calcolare le variabili ap1,bp1,cp1,dp1,

che contengono la funzione di trasferimento della sensitività.

1

u12

u2

5

u5

4

u4

3

u3

Mux

Mux1

++

Sum


System


Control

5

e5

4

e4

3

e3

2

e2

1

e1

emu

Demux1

Fig. 4-5

63


sv1=max([sigma(ap1,bp1,cp1,dp1,w);-inf*ones(size(w))]);

calcola la risposta in frequenza della sensitività.

• L’istruzione successiva è:

ow1=20*log10(min(sw1./sv1));

dove la variabile ow1 è una misura in decibel della distanza tra il valore singolare

massimo della sensitività e il limite superiore ls(s), questa variabile potrebbe essere

utilizzata nel calcolo dell’indice di qualità.

• Il blocco di istruzioni:

if plt,

f1=fig.;

subplot(2,1,1);

semilogx(w,20*log10(sw1),'w');

hold on

semilogx(w,20*log10(sv1),'r');

title('sensitivity');

end

serve per disegnare l’andamento della sensitività e di ls(s) nel caso che la variabile plt

sia diversa da zero.

64

La sezione che segue calcola la sensitività di controllo dell’impianto a circuito chiuso.


[ap2,bp2,cp2,dp2]=linmod('hloop2');

dove hloop2 è lo schema della Fig. 4-5, serve per calcolare le variabili ap2,bp2,cp2,dp2,

che contengono la funzione di trasferimento della sensitività di controllo.

emu

Demux1

1

e12

e23

e34

e45

e5

1

u12

u2

5

u5

4

u4

3

u3

Mux

Mux1


System


Control

++

Sum

Fig. 4-6

• L’istruzione:

sv2=max([sigma(ap2,bp2,cp2,dp2,w);-inf*ones(size(w))]);

calcola la risposta in frequenza della sensitività di controllo.

65

• L’istruzione successiva è:

ow2=20*log10(min(sw2./sv2));

dove la variabile ow2 è una misura in decibel della distanza tra il valore singolare

massimo della sensitività di controllo e il limite superiore lm(s), questa variabile

potrebbe essere utilizzata nel calcolo dell’indice di qualità.

La sezione successiva calcola le risposte nel tempo (una per ogni canale), del segnale di

controllo con un gradino come comando di riferimento nell’ingresso relativo:

• Le relative istruzioni sono:

mxu=zeros(1,sz);

for n=1:sz;

[y,xsys,t]=step(ap2,bp2(:,n),cp2(n,:),dp2(n,n),1);

mxu(n)=max(y);

end

• L’istruzione:

mxu=max(abs(mxu));

estrae il valore massimo del segnale di controllo, anche questo valore potrà essere

utilizzato nel calcolo dell’indice di qualità.

66


if plt,

subplot(2,1,2);

semilogx(w,20*log10(sw2),'w');

hold on

semilogx(w,20*log10(sv2),'b');

title('control sensitivity');

end

serve per disegnare l’andamento della sensitività di controllo e di lm(s),se plt è diversa

da zero.

La sezione che segue analizzerà completamente l’ impianto a ciclo chiuso.

Le matrici nello spazio degli stati dell’ impianto a ciclo chiuso potrebbero essere calcolate

interconnettendo il sistema A,B,C,D ed il controllore Ak Bk Ck Dk, ad ogni modo è stato

deciso di ottenere queste matrici applicando il controllore direttamente all’ impianto non

lineare e linearizzando poi tutto l’impianto, poichè questa è la maniera più affidabile, nel

senso che l’impianto lineare a ciclo chiuso così ottenuto è più vicino a quello vero di quanto

lo sia quello calcolato direttamente interconnettendo i sistemi lineari tra loro.

• L’istruzione:

s4=lmc2('nl3c',[0;0;0;0;0],[3;0;0],[0;0;0]);

Restituisce il modello lineare impacchettato nella variabile s4.

La funzione “lmc2” linearizza l’ impianto nel punto operativo di riferimento (vedere l’

help della funzione per maggiori dettagli).

67

• Le istruzioni:

mp4=max(real(spoles(s4)));

mf4=abs(max(spoles(s4)));

calcolano le variabili mp4 ed mf4 che rappresentano rispettivamente la massima parte

reale di un polo, e il massimo valore assoluto di un polo, tra tutti i poli dell’ impianto.

La sezione seguente calcola le risposte del tempo dell’ impianto a ciclo chiuso per ogni

canale con un gradino come comando di riferimento nel relativo ingresso.

• Le istruzioni:

[A4,B4,C4,D4]=unpck(s4);

g4=D4-C4*inv(A4)*B4;

“spacchettano” l’impianto a ciclo chiuso, e calcolano la matrice di trasferimento per s=0

ogni colonna della quale, come è noto, rappresenta anche il valore delle uscite per t=∞

quando all’ ingresso corrispondente alla colonna viene applicato un gradino.

• Istruzione di inizio ciclo:

for n=1:sz;

68

• Le istruzioni seguenti:

[y,xsys,t]=step(A4,B4(:,n),C4(n,:),D4(n,n),1);

g0=g4(n,n);

Calcolano la risposta al gradino per l’ n-esimo canale ed il suo valore in t=∞.

• L’istruzione:

ts4(n)=max(t'.*(abs(y-g0)>.05*abs(g0)));

calcola il settling time della risposta y(t) (notare che g0=y(∞)):

• La linea che calcola l’ overshoot delle stesse risposte è:

os4(n)=abs(max(y-g0)/g0);

• Istruzione di fine ciclo:

end

• Le istruzioni:

ts4=max(ts4);

os4=max(os4);

Calcolano il valore massimo dell’ overshoot e del setting time tra i diversi canali:

69

C’è poi una parte della funzione che, se la variabile nl è diversa da zero, permette di fare

una simulazione sull’ impianto non-lineare.

Questa parte richiede tempi di calcolo molto lunghi e per questo non è stata utilizzata

nell’ottimizzazione, ma soltanto per analizzare i risultati a procedura di ottimizzazione

conclusa.

Essendo poi abbastanza complessa ed irrilevante per la comprensione dell’algoritmo, non è

stata qui presa in considerazione.

L’indice di qualità è ottenuto come prodotto di un vettore q contenente tutte le variabili

che possono essere utilizzate, per un vettore di coefficenti scelti dal progettista, e

dipendenti dalle specifiche di impianto.

• Le istruzioni

q =[ x(1); x(2); x(3); x(4); ow1; ow2; gf; mxu; mp4; mf4; os4; ts4]

f = [0, 0, 0, 0, 0, 0, 1e-1, 0, 0, 0, 0, 1 ]*q

servono appunto per definire il vettore delle variabili che possono essere usate, e

calcolare l’ indice di qualità f come prodotto di questo vettore per un vettore di

coefficienti scelti dall’ utente.

I vincoli sono espressi dalla relazione g <= 0.

• L’istruzione che definisce il vettore g è:

g=[mxu-2;1e7*mp4]

70

La funzione costo e i vincoli mostrati sopra, sono proprio quelli utilizzati nel nostro caso.

Il vincolo è fatto in modo che il sistema a ciclo chiuso sia stabile (1e7*mp4 < 0), che il

segnale di controllo sia sempre minore di 2 con un gradino in ingresso (mxu-2 < 0).

La quantità da minimizzare è una combinazione tra il settling time (1*ts4), ed il valore

finale di γ (1e-1*gf).

• L’istruzione:

if min([ finite(f) finite(g') ])==0, f=inf;g=inf; end

pone f e g ad infinito se qualcuno dei loro valori non è un numero finito.

Cio’ viene fatto per evitare possibili propagazioni di NaN (Not A Number) attraverso la

procedura di ottimizzazione.

• L’istruzione:

fc=max(max(g),0)*1e5+f;

calcola l’indice utilizzato per confrontare i vari controllori ( che è uguale ad f se g <= 0

).


if fc<F,

F=fc;

X=[reshape(x,4,1) X(:,1:3)]

end

sposta a destra le prime tre colonne di X, perdendo l’ultima e salva il vettore x relativo al

miglior controllore, come prima colonna di X.

71

Il mantenimento in memoria dei migliori controllori trovati rende possibile fermare la

routine di ottimizzazione e vedere quale è il controllore migliore trovato fino a quel

momento.

• Istruzione di fine funzione:

end

4.3.2 Funzione “hopt”.

Questo m-file è l’interfaccia esterna dell’algoritmo di ottimizzazione.

La variabile s0 è destinata a contenere il sistema che deve essere controllato, ed è l’unica

variabile che deve essere messa nello spazio di lavoro prima di far funzionare hopt.

• L’istruzione:

global plt nl sz A B C D Ak Bk Ck Dk a1 b1 c1 d1 a2 b2 c2 d2 X F N uequ yequ CT

inizializza come “globali” le variabili in essa contenute.

• L’istruzione:

plt=0;nl=0;

pone a zero le variabili plt ed nl , in modo da saltare la routine di visualizzazione e la

routine per calcolare la risposta dell’ impianto non lineare in hidx.

72

• Le istruzioni:

[A,B,C,D]=unpck(s0);

sz=size(B,2);

format long

CT=0;

F=inf;

X=zeros(4,4);

servono per inizializzare le variabili globali.

La variabile options contiene opzioni per la funzione “constr” utilizzata per

l’ottimizzazione ( vedere help FOPTIONS per ulteriori dettagli ).

• Le istruzioni:

options=[];

options(14)=40;

servono per settare a 40 il massimo numero delle iterazioni, e lasciare le altre opzioni ai

valori di default.

Dopo la scelta della funzione da minimizzare, la scelta del punto di partenza è la parte

più critica della procedura di ottimizzazione.

Il punto di partenza dovrebbe essere scelto secondo alcuni metodi come quello presentato

in sezione 4.2; è sicuramente utile, se non necessario, lanciare alcune volte la funzione

hidx, per conoscere le quantità con cui abbiamo a che fare e per migliorare l’idoneità degli

indici, dei punti di partenza e dei limiti.

73

I limiti definiscono la regione in cui stiamo cercando la soluzione, che non dovrebbe essere

troppo grande se scegliamo un buon punto di partenza.

Per il nostro problema, punto di partenza e limiti, sono stati trovati come risultato di varie

ricerche per il punto minimo e sono così definiti:

• L’istruzione:

x=[6.63071609934561 -0.00000001000000 29.24136343802332 0.67538629717924]';

definisce il punto di partenza.

• Le istruzioni:

vlb=x-[5 .5 10 1]';

vub=x+[5 .5 10 1]';

definiscono i limiti inferiore e superiore.

La parte che segue serve per la ricerca di un punto di partenza migliore nella regione

definita dai limiti, e può essere saltata se siamo sicuri che non esistono punti di partenza

migliori nelle vicinanze di quello che abbiamo scelto.

• L’istruzione:

eval('[f,g]=hidx(x);','');

esegue hidx sul punto di partenza.

74

La funzione hidx viene chiamata indirettamente con la sintassi eval('[f,g]=hidx(x);','');

perchè così, se questa si arresta per qualche motivo, il programma non si ferma e considera

un altro punto.

In ognuna delle quattro direzioni dello spazio di x , avremo un segmento (indichiamo con

lseg la sua lunghezza), che và dal limite inferiore a quello superiore, ed avrà il punto di

partenza fornito all’ inizio come suo punto medio.

In questo segmento vengono considerati due punti vicini alle estremità, essi sono:

limite inferiore più 1/5*lseg, e limite superiore meno 1/5*lseg.

Avremo quindi un totale di 16 punti che individuano un ipercubo nello spazio delle x.


stp1=(vub(1)-vlb(1))/5;

for i1=[vlb(1)+stp1 vub(1)-stp1],

stp2=(vub(2)-vlb(2))/5;


stp3=(vub(3)-vlb(3))/5;


stp4=(vub(4)-vlb(4))/5;


x=[i1 i2 i3 i4];

eval('[f,g]=hidx(x);','');

end

end

end

end

effettua le ricerca del miglior punto tra i 16 suddetti.

75

• Il seguente blocco di istruzioni:

error=0;

x=X(:,1);

while (error < 3)&(CT < options(14)),

eval('constr(''hidx'',x,options,vlb,vub);',...

'x=X(:,1)+rand(size(X(:,1)));error=error+1');

end

seleziona come punto di partenza il punto migliore trovato ed esegue la funzione

“constr”, se qualcosa va male nel corso dell’ottimizzazione, essa viene rieseguita (per

non più di due volte) partendo dal migliore tra punti disponibili in quel momento.

• Le istruzioni:

nl=1;plt=1;

[f,g]=hidx(X(:,1));

save sc3 Ak Bk Ck Dk X s0

visualizzano i risultati e salvano i dati relativi al miglior controllore.

76

5. RisultatiLa procedura di ottimizzazione descritta nel capitolo precedente, applicata al modello

lineare ridotto, con il punto di partenza e la funzione costo precedentemente mostrati, ha

dato come risultato un controllore H∞ (5 ingressi, 5 uscite, 29 stati, stabile), per il punto di

equilibrio scelto.

5.1 Sensitività e controllo della sensitività.

La Fig. 5-1 mostra il valore singolare massimo della sensitività, e della sensitività di

controllo, ed i relativi limiti superiori ls(s) e lm(s).

10−6

10−4

10−2

100

102

104

106

−150

−100

−50

0

50sensitivity

10−6

10−4

10−2

100

102

104

106

−200

−100

0

100control sensitivity

Fig. 5-1

Come possiamo vedere, alle basse frequenze la sensitività, che si trova leggermente al di

sotto del suo limite superiore, è abbastanza piccola, e discende con 20 dB/dec,il che

77

dovrebbe garantire buone prestazioni in termini di reiezione dei disturbi ed inseguimento

del segnale di riferimento.

Alle alte frequenze la sensitività di controllo è molto al di sotto dei propri limiti, e discende

con 40 dB/dec.,il che dovrebbe garantire un buona reiezione del rumore dei sensori, e

buona robustezza a fronte di incertezze additive.

5.2 Impianto a ciclo chiuso linearizzato.

Chiameremo impianto “a ciclo chiuso linearizzato” l’ impianto ottenuto dalla

linearizzazione del circuito chiuso tra il sistema non lineare di Fig. 2-9 e controllore.

Questo è un impianto stabile di 5 ingressi, 5 uscite, 65 stati, le cui caratteristiche saranno

descritte in questa sezione.

Per via del modo in cui il Simulink ordina gli stati, i primi 31 stati sono esattamente quelli

descritti in sezione 3.1, ci sono poi i 29 stati del controllore, e gli ultimi 5 stati

appartengono alla funzione di trasferimento del blocco fcn degli attuatori.

La Fig. 5-2 mostra i poli e gli zeri del sistema, i valori singolari, controllabilità e

osservabilità degli stati, e guadagno di ogni elemento della funzione di trasferimento ad una

frequenza di 0.1 rad/sec.

78

−2 −1 0 1

x 106

−20

−10

0

10

20

Real Axis

Imag

Axi

s


10−4

10−2

100

102

−200

−100

0

100

Frequency (rad/sec)

Sin

gula

r V

alue

s dB


0 20 40 60 8010

−5

100

105

1010

1015

state

ctrb

(re

d), o

bsv

(blu

e)


1 2 3 4 5−5

−4

−3

−2

−1

input

outp

ut

max.sv.(G(jw,i,j))

Fig. 5-2

A causa dei tre zeri attorno a -1*10^(6) non possiamo vedere bene le posizioni degli altri

zeri e dei poli.

Per ciò che riguarda controllabilità e osservabilità, è interessante notare che tutti gli stati del

controllore sono molto più controllabili che osservabili, e ad ogni modo la controllabilità e

l’osservabilità hanno circa lo stesso valore medio (il che non si verificava nell’ impianto a

ciclo aperto).

Dal diagramma dei valori singolari possiamo vedere che fino a 1 rad/sec il guadagno del

sistema è 1 in ogni direzione, questo dovrebbe garantire un buon tracking del segnale di

riferimento a queste frequenze.

Possiamo anche vedere che il guadagno del sistema diminuisce velocemente dopo 10

rad/sec.

79

Infine, a 0.1 rad/sec la struttura del sistema è molto diagonale, come si potrebbe mostrare

anche con una analisi dei vettori singolari: questo dovrebbe garantire un buon

disaccoppiamento (almeno intorno al punto di lavoro ).

La Fig. 5-3, mostra le risposte al gradino relative ad ogni elemento della funzione di

trasferimento.

0 5 100

1

2

0 1 2

x 106

−5

0

5x 10

−4

0 1 2

x 106

−2

0

2x 10

−4

0 1 2

x 106

−1

0

1x 10

−5

0 1 2

x 106

−2

0

2x 10

−4

0 5 10

x 105

−2

0

2x 10

−5

0 100 2000

1

2

0 5 10

x 105

−2

0

2x 10

−5

0 5 10

x 105

−2

0

2x 10

−6

0 5 10

x 105

−2

0

2x 10

−5

0 1 2

x 106

−5

0

5x 10

−4

0 1 2

x 106

−5

0

5x 10

−4

0 5 100

1

2

0 50 100−0.05

0

0.05

0 5 10

x 105

−5

0

5x 10

−4

0 5 10

x 105

−1

0

1x 10

−5

0 1000 2000−0.02

0

0.02

0 5 10

x 105

0

x 10−4

0 20 400

1

2

0 5 10

x 105

−5

0

5x 10

−6

0 1 2

x 106

0

x 10−4

0 1 2

x 106

−2

0

2x 10

−4

0 5 10

x 105

−2

0

2x 10

−4

0 5 10

x 105

−2

0

2x 10

−5

0 10 200

1

2

Fig. 5-3

Si può vedere che il settling time va dai 2 ai 4 secondi, e l’ accoppiamento tra i canali è

molto limitato.

80

5.3 Impianto non lineare controllato

In questa sezione analizzeremo tutte le risposte dell’ impianto non lineare controllato.

Nelle figure dalla Fig. 5-4 alla 5-17, sono riportate le risposte del sistema quando si insegue

una rampa di pendenza 0,1 e ampiezza 15, applicata in un unico ingresso.

Questo è un modo per capire quanto è ampia la regione circostante il punto di lavoro dove

il controllore lavora bene, e come la capacità del sistema di seguire comandi piuttosto

semplici, diminuisce quando ci muoviamo lontano dal punto di lavoro.

• Primo canale (regolazione del roll):

0 50 100 150 200 250−2

0

2

4

6

8

10

12

14

16

time

chan

nel1

nonlinear plant response

Fig. 5-4

l’ inseguimento del segnale di comando è perfetto.

81

• Secondo canale (sway):

0 50 100 150 200 2500

2

4

6

8

10

12

14

16

time

chan

nel2


Fig. 5-5

l’inseguimento è molto buono fino a 14 metri.

• Terzo canale (yaw):

0 50 100 150 200 250−2

0

2

4

6

8

10

12

14

16

18

time

chan

nel3


Fig. 5-6

82

Il terzo canale sembra essere il più critico: dopo 13 gradi, la non-linearità dell’ impianto

inizia ad agire, anche se l’attuatore è ancora in grado di seguire il riferimento ma con un

errore che può raggiungere i 2 gradi.

• Canali quarto (heave), e quinto (pitch):

0 50 100 150 200 2500

2

4

6

8

10

12

14

16

time

chan

nel4


Fig. 5-7

0 50 100 150 200 2500

2

4

6

8

10

12

14

16

time

chan

nel5


L’inseguimento del comando è molto buono.

83

Si potrebbe mostrare che in tutti questi casi l’interazione fra i canali è scarsa.

Nelle figure dalla Fig. 5-8 alla Fig. 5-12 è riportata la risposta del sistema ad un comando

più complesso.

Viene dato come comando una rampa di ampiezza 5 metri con una pendenza di 0.2 ,prima

per lo sway , dopo 150 secondi per l’heave e dopo intervalli di 150 secondi lo stesso

comando (5 gradi di ampiezza, pendenza 0.2) viene dato per pitch, yaw e roll.

Tutto ciò allo scopo di vedere come si comporta il sistema quando è soggetto

contemporeanamente a comandi su canali differenti.

• Comando dopo 1 secondo: 5 metri di sway:

0 100 200 300 400 500 600 700 800−1

0

1

2

3

4

5

6

Fig. 5-8

84

• Comando dopo 150 secondi: 5 metri di heave.

0 100 200 300 400 500 600 700 800−1

0

1

2

3

4

5

6

Fig. 5-9

• Comando dopo 300 secondi: 5 gradi di pitch.

0 100 200 300 400 500 600 700 800−1

0

1

2

3

4

5

6

Fig. 5-10

85

• Comando dopo 450 secondi: 5 gradi di yaw.

0 100 200 300 400 500 600 700 800−1

0

1

2

3

4

5

6

Fig. 5-11

• Comando dopo 600 secondi: 5 gradi di roll.

0 100 200 300 400 500 600 700 800−1

0

1

2

3

4

5

6

Fig. 5-12

86

Come possiamo notare, il comportamento del sistema è molto buono, può seguire tutti i 5

comandi indipendentemente con un’interazione molto limitata.

Nei grafici daFig. 5-13 a Fig. 5-17 è riportata la risposta del sistema relativa ad un comando, simile al

precedente, avente però ampiezza 10 anzichè 5.

• Comando dopo 1 secondo: 10 metri di sway.

0 100 200 300 400 500 600 700 800−2

0

2

4

6

8

10

12

Fig. 5-13

87

• Comando dopo 150 secondi: 10 metri di heave:

0 100 200 300 400 500 600 700 800−2

0

2

4

6

8

10

12

Fig. 5-14

• Comando dopo 300 secondi: 10 gradi di pitch:

0 100 200 300 400 500 600 700 800−2

0

2

4

6

8

10

12

Fig. 5-15

88


0 100 200 300 400 500 600 700 800−2

0

2

4

6

8

10

12

Fig. 5-16


0 100 200 300 400 500 600 700 800−0.5

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Fig. 5-17

89

Come prevedibile il comportamento del sistema non è stato buono come nel caso

precedente; infatti l’interazione tra canali è maggiore, ed in particolare l’interazione tra

pitch ed heave.

Il sistema segue difficilmente comandi più severi e contemporanei per heave, pitch, e roll.

Da notare però che, in pratica, di un comando così articolato ne avremo bisogno raramente,

e che comunque i comandi isolati di heave e pitch funzionano molto bene.

In conclusione, il funzionamento di questi controllori è molto buono nell’arco dei 5 metri-5

gradi , e potrebbe essere accettabile nell’arco di 10 metri-10 gradi, a meno di non voler

inseguire contemporeanamente 5 comandi abbastanza ampi..

Fuori dall’arco dei 15 metri-15gradi le prestazioni del sistema diminuiscono rapidamente,

quindi, se vogliamo controllare il veicolo fuori da questi limiti, dovremo affrontare la non-

linearità dell’ impianto e progettare un controllore gain-scheduling o self-tuning.

90

6. Controllo LQG

Allo scopo di confrontare il controllo H∞ con un controllo di tipo standard è stato realizzato

un controllo di tipo lineare quadratico gaussiano (LQG).

6.1 TeoriaIl problema del controllo LQG è il seguente:

Supponiamo di avere un sistema nello spazio degli stati:

[ ] [ ][ ] [ ]&( )

( )

( )

( )

( )

( )x t

y t

A B

C D

x t

u tW

V

d t

d t= +

1

2

1

2

0

0

1

2Eq 6-1

dove il vettore d(t) è un rumore bianco con densità spettrale di potenza:

Φdd s I( ) = Eq 6-2

percui i rumori w(t) e v(t) entranti rispettivamente sugli stati e sulle uscite hanno densità

spettrali di potenza:

Φ

ΦΦ

wwT

vvT

wv

s W W

s V V

s

( )

( )

( )

= = ≥

= = >=

0

0

0Eq 6-3

Il problema è allora trovare un controllo tale da minimizzare l’indice di costo:

J E z t z t dtt

Tt

= ∫

→∞

lim ( ) ( )0

Eq 6-4

91

dove:

z tz t

z tQ

R

x t

u t( )

( )

( )

( )

( )=

=

1

2

1

2

1

2

0

0Eq 6-5

per cui, all’ interno dell’ integrale le matrici che pesano stati e ingressi sono:

Q Q

R R

T

T

= ≥

= >

0

0Eq 6-6

Il controllore che risolve il problema ha la seguente funzione di trasferimento:

K s K sI A BK K C K DK Kc c f f c f( ) ( )= − − + + + −1 Eq 6-7

con

K R B Ric H

K Ric J C V

cT

fT

=

=

−

−

12

21

( )

( )Eq 6-8

dove H2 ed J2 sono le seguenti matrici hamiltoniane:

HA BR B

Q A

JA C V C

W A

T

T

T T

T

2

1

2

1

=−

− −

=−

− −

−

− Eq 6-9

Da notare che K(s) contiene tanti stati quanti ne contiene l’impianto.

92

6.2 Analogie con il problema H∞∞∞∞

Data una funzione di trasferimento T(s) si chiama “norma-due” della funzione la quantità:

( )T s tr T s T s dsT

j

j

( ) ( ) ( )2

1

2= ∫

− ∞

∞

π Eq 6-10

dove con tr(M) si indica la traccia della matrice M, ossia la somma degli elementi sulla

diagonale principale.

La norma-due è la radice quadrata della potenza media in uscita dal sistema quando in

ingresso è presente un rumore bianco, essa quindi dà precise informazioni sul guadagno di

potenza medio del sistema ( la norma-infinito riguardava invece il guadagno massimo ).

L’indice dell’ Eq. 6-4 ha un analogo corrispondente nel dominio della frequenza[ii]:

( )J tr s dszz

j

j

= ∫− ∞

∞1

2πΦ ( ) Eq 6-11

Applicando la legge di controllo:

u s K s y s( ) ( ) ( )= Eq 6-12

dove K(s) è quello descritto nell Eq. 6-7, l’unico segnale entrante del sistema rimane d(s).

Tra d(s) e z(s) ci sarà una funzione di trasferimento Tzd(s) (che ovviamente dipenderà anche

dal controllo K(s)), percui si ha:

Φ Φzz zdT

dd zd zdT

zds T s s T s T s T s( ) ( ) ( ) ( ) ( ) ( )= = Eq 6-13

e quindi l’indice che viene minimizzato dal controllo LQG è:

( )J T szd= ( ) 2

2Eq 6-14

93

Da notare che Tzd(s)=LFT(K,P) dove P è l’impianto descritto dall’ Eq. 6-15:

&( )

( )

( )

( )

( )

( )

( )

( )

x t

z t

z t

y t

A W B

Q

R

C V D

x t

d t

d t

u t

1

2

1

2

1

2

1

2

1

2

1

2

0

0 0 0

0 0 0

0

=

Eq 6-15

graficamente:

+

+

+

Sum_y

+

+

Sum_x

1

d12

d23

u

3

y

B

W^(1/2)

V^(1/2)

R^(1/2)

Q^(1/2)

2

z21

z1C

D

sI−A)^(−1)

Fig. 6-1

Si può dimostrare [4] che il metodo H∞ per valori di γ molto elevati di fatto minimizza la

norma-due della LFT, cioè:

lim ( ( ), ) ( , )γ

γ→∞ ∞ =LFT K P LFT K P 2 Eq 6-16

A questo punto, possiamo concludere che il problema LQG è un caso particolare del

problema H∞ nel senso che, se aggiungiamo al sistema da controllare gli ingressi fittizi d(t),

e le uscite fittizie z(t) “pesandoli” con le matrici W1/2 V1/2 Q1/2 R1/2,come in Fig. 6-1, ed

applichiamo la procedura illustrata nel capitolo 4 con un γ → ∞, il controllore risultante

sarà un controllore LQG, che minimizzerà l’indice dell’ Eq. 6-4.

94

6.3 DifferenzeOltre al fatto che la procedura LQG è indubbiamente molto più semplice, ci sono tre

differenze che vale la pena di notare, in particolare le ultime due sono fondamentali e

rappresentano la vera ragione del successo del metodo H∞.

1) Il controllo LQG minimizza la norma-due di una funzione di trasferimento, percui

minimizza un guadagno medio di potenza. La procedura H∞ invece minimizza la

norma-infinito e quindi un guadagno massimo.

2) La f.d.t minimizzata da LQG non è la f.d.t. vista da una qualche incertezza applicata

all’ impianto, ecco perchè LQG non è un controllo robusto. La f.d.t minimizzata da H∞

può essere scelta a piacere e normalmente comprende la f.d.t. vista da una incertezza

additiva applicata all’ impianto (sensitività di controllo).

3) Nel metodo H∞ normalmente le matrici di peso sono funzioni della frequenza e questo

vuol dire che:

A) Il controllore H∞ avrà degli stati in più rispetto al controllore LQG e precisamente

tutti gli stati delle matrici di peso.

B) Con un controllo LQG non riusciamo a modellare l’andamento in frequenza della

funzione di trasferimento che viene minimizzata, cosa che invece riesce

facilmente con il controllo H∞.

6.4 Risultati

Dopo diverse prove sono state scelte queste matrici di peso:

R V I

Q I

W I

= =

=

= −

5

191 37

192 39

10

10

*

*

.

.

si è calcolato il controllore descritto dalla 6-7 e lo si è applicato all’ impianto, con i seguenti

risultati:

95

6.4.1 Sensitività e sensitività di controllo

La Fig. 6-2 mostra il massimo valore singolare della sensitività e della sensitività di

controllo:

10−6

10−4

10−2

100

102

104

106

0

2

4

6sensitivity

10−6

10−4

10−2

100

102

104

106

−150

−100

−50

0

50control sensitivity

Fig. 6-2

La sensitività di controllo potrebbe essere accettabile, la sensitività però ha un picco attorno

ad 1 rad/sec e anche a basse frequenze si mantiene costante sopra 0dB anzichè scendere.

Questo ci fa supporre che le prestazioni del controllore in termini di inseguimento dei

comandi non saranno buone.

96

6.4.2 Impianto a ciclo chiuso linearizzato

In questa sezione chiameremo impianto “a ciclo chiuso linearizzato” l’ impianto ottenuto

dalla linearizzazione del circuito chiuso tra il sistema non lineare di Fig. 2-9 e controllore

LQG. Questo è un impianto stabile di 5 ingressi, 5 uscite, 55 stati, le cui caratteristiche

saranno descritte.

Per via del modo in cui il Simulink ordina gli stati, i primi 31 stati sono esattamente quelli

descritti in sezione 3.1, ci sono poi i 19 stati del controllore, e gli ultimi 5 stati

appartengono alla funzione di trasferimento del blocco fcn degli attuatori.

La Fig. 6-3 mostra i poli e gli zeri del sistema, i valori singolari, controllabilità e

osservabilità degli stati, e guadagno di ogni elemento della funzione di trasferimento ad una

frequenza di 0.1 rad/sec; ne derivano le seguenti considerazioni:

−40 −20 0 20−4

−2

0

2

4

Real Axis

Imag

Axi

s


10−4

10−2

100

102

−300

−200

−100

0

100

Frequency (rad/sec)

Sin

gula

r V

alue

s dB


0 20 40 6010

−10

10−5

100

105

state

ctrb

(re

d), o

bsv

(blu

e)


1 2 3 4 5−5

−4

−3

−2

−1

input

outp

ut

max.sv.(G(jw,i,j))

Fig. 6-3

97

Non ci sono zeri a grande distanza dall’ origine.

Per ciò che riguarda controllabilità e osservabilità, non si riscontra una sostanziale

differenza tra l’osservabilità e la controllabilità degli stati del controllore.

L’ andamento dei valori singolari è molto meno “liscio” di quello che si aveva con il

controllo H∞, inoltre anche a basse frequenze il guadagno del sistema non è lo stesso in

ogni direzione, e questo è una conferma del fatto che probabilmente il tracking del segnale

di riferimento non sarà molto buono.

Anche in questo caso il guadagno del sistema diminuisce velocemente dopo 10 rad/sec.

Infine, come si poteva immaginare, a 0.1 rad/sec la struttura del sistema non è così

diagonale come nel caso precedente.

La Fig. 6-4, mostra le risposte al gradino relative ad ogni elemento della funzione di

trasferimento.

0 50000

0.05

0.1

0 100 2000

0.1

0.2

0 2000 4000−0.2

0

0.2

0 100 200−0.02

0

0.02

0 2000 4000−0.2

−0.1

0

0 50 1000

0.5

0 50 1000

0.5

1

0 50 1000

0.5

1

0 50 1000

0.1

0.2

0 50 100−0.4

−0.2

0

0 200 400−0.4

−0.2

0

0 50 100−0.5

0

0 2000 4000−0.5

0

0.5

0 50 100−0.1

0

0.1

0 5000

0.1

0.2

0 1 2

x 104

−5

0

5x 10

−3

0 50 100−0.01

−0.005

0

0 5000 10000−0.01

0

0.01

0 50 1000

0.2

0.4

0 50 1000

0.1

0.2

0 500 1000−0.2

0

0.2

0 50 100−0.1

0

0.1

0 500 1000−0.5

0

0.5

0 500

0.2

0.4

0 200 4000

0.2

0.4

Fig. 6-4

98

Si può vedere che il settling time va dai 20 ai 40 secondi, ( un ordine di grandezza in più

rispetto al caso precedente ) e l’ accoppiamento tra i canali è considerevole.

6.5 Impianto non lineare controllato

In questa sezione analizzeremo la risposta dell’ impianto non lineare controllato con il

controllore LQG, allo stesso comando dato come ultimo caso nel capitolo 5 cioè la rampa

di ampiezza 5 metri con una pendenza di 0.2 , prima per lo sway , dopo 150 secondi per

l’heave e dopo intervalli di 150 secondi lo stesso comando (5 gradi di ampiezza, pendenza

0.2) viene dato per pitch, yaw e roll.

• Comando dopo 1 secondo: 15 metri di sway:

0 100 200 300 400 500 600 700 800−25

−20

−15

−10

−5

0

5

10

15

Fig. 6-5

99

• Comando dopo 150 secondi: 15 metri di heave.

0 100 200 300 400 500 600 700 800−0.5

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Fig. 6-6

• Comando dopo 300 secondi: 15 gradi di pitch.

0 100 200 300 400 500 600 700 800−3

−2

−1

0

1

2

3

4

Fig. 6-7

100


0 100 200 300 400 500 600 700 800−50

−40

−30

−20

−10

0

10

Fig. 6-8


0 100 200 300 400 500 600 700 800−1

0

1

2

3

4

5

6

Fig. 6-9

101

La Fig. 6-10 mostra , per un confronto più facile, le risposte che si avevano nel caso H∞

(tutte insieme):

0 100 200 300 400 500 600 700 800−2

0

2

4

6

8

10

12

Fig. 6-10

Al di la del fatto che il comando è molto impegnativo e che probabilmente è utile solo per

confrontare i tipi di controllo, si vede subito che il confronto con il controllo H∞ è

improponibile: il controllore LQG non riesce a disaccoppiare i canali, e dopo 450 secondi,

in corrirpondenza del comando sullo yaw, il sistema esce si esce definitivamente dalla zona

di linearità.

102

7. Controllo Adattivo

Allo scopo di gestire meglio le nonlinearità dell’ impianto si è provato a realizzare un

controllo di tipo adattivo, ossia un meccanismo che identifica alcuni parametri dell’

impianto durante il funzionamento e si adatta agli eventuali cambiamenti di tali parametri

ridisegnando il controllore.

Lo schema risultante è mostrato in Fig. 7-1:

Demux

Demux2

5T

4Z

3P

2Y

1F

Yaw

Pitch

Roll

Heave

Sway

1 f 2 y 3 p 5 t

Mux Mux2

4 z

Mux

yequ

Eq. values

CABLE − VEHICLESYSTEM

Demux

Demux1

R2D

rad2deguequ(9:14)

Ship vel & acc

Actuators −

+

Sum1

Mux Mux

1

Gainreference

roll

sway

yaw

heave

pitch

1/s

Lint11/s

Lint21/s

Lint31/s

Lint41/s

Lint5

Mux

Mux3

Clock zoh4

tTws4

ada1

Adaptation


Control

+++−

Sum2

blwnSign

1

Gain1

Mux Mux4

5_2_8

Demixing

D2R

deg2rad

y Tws2

zoh2zoh3

u Tws3

Load Data and model initialization

Fig. 7-1

In definitiva si cerca di realizzare un controllore che possa adattarsi ai cambiamenti dell’

impianto dovuti a fattori di tipo esterno (disturbi) o interni (spostamento dello stato dal

punto operativo).

Il blocco ‘ada1’ riceve in entrata ingressi ed uscite dell’

impianto, li campiona, e li usa per costruirsi un nuovo modello che a sua volta serve per

calcolare un nuovo controllore.

Sono necessarie quindi funzioni per l’identificazione (magari ricorsiva) di sistemi

multivariabili.

103

7.1 Identificazione :

Il problema dell’ identificazione per sistemi lineari tempoinvarianti tempodiscreti [vi] è:

[ ]h g V e h g

y t g t u t h t e t

h g

n(. ), (. ) arg min ( (., (. ), (. ))

( ) ( ) ( ) ( ) ( )

( ) , ( )

== ⊗ + ⊗

= =0 1 0 0

Eq 7-1

Nel caso di identificazione “parametrica” le funzioni h(.) e g(.) dipendono da un vettore di

parametri θ , e la notazione che si usa é questa :

y t G q u t H q e t( ) ( , ) ( ) ( , ) ( )= +ϑ ϑ Eq 7-2

G ed H sono le z-trasformate delle risposte impulsive g ed h:

G q g k q

H q h k q

k

k

k

k

( , ) ( , )

( , ) ( , )

ϑ ϑ

ϑ ϑ

= ∑

= +∑

=

∞ −

=

∞ −

1

11

Eq 7-3

Le quantità X(q)y(t) rappresentano x(t)⊗ y(t) nel dominio del tempo, o X(q)Y(q) nel

dominio della frequenza.

La quantità da minimizzare è in genere la norma del segnale e(t), il quale rappresenta la

cosiddetta “innovazione” ossia la parte di y(t) che non puo’ essere predetta dai valori

passati di y(t) ed u(t).

Vn en

e t e tt

nT( (. )) ( ) ( )= ∑

=

1

2 1 Eq 7-4

104

7.2 Algoritmi di identificazione MIMO :

Il system identification toolbox del Matlab [vii], supporta le seguenti funzioni per

l’identificazione multivariabile:

Consideriamo la funzione arx:

Si suppone che il modello abbia la seguente struttura:

G q A q B q

H q A q

( ) = ( ) ( )

( ) = ( )

-1

-1 Eq 7-5

cioè obbedisca all’ equazione:

A q y t B q u t e t( ) ( ) = ( ) ( ) + ( ) Eq 7-6

dove A(q) e B(q) sono matrici polinomiali i cui elementi sono:

a a q a q

b b q b q

ij ij ij ij

naij naij

ij ij

nkij

ij

nbij nbij

= + + ... +

= + ... +

δ 1 1

1

− −

− − Eq 7-7

la struttura è quindi individuata dalle 3 matrici na nb ed nk,

la prima è quadrata di dimensione ny, le altre due rettangolari di dimensione ny*nk.

la sintassi è la seguente:

th=arx([[[[y u]]]],[[[[na nb nk]]]]);

dove ogni riga di y (u) è il vettore delle uscite (ingressi) ad un certo istante, th è una matrice

che contiene la stima dei parametri, la loro varianza,ed altre utili informazioni.

Il modello arx è a regressione lineare, cioè le uscite sono una funzione lineare del vettore di

parametri θ il quale puo’ essere calcolato quindi con una semplice pseudoinversione di una

matrice contenente i dati, l’algoritmo è pertanto molto veloce.

105

Un altra funzione adatta all’identificazione multivariabile è iv4

Il modello si suppone avere la stessa struttura vista in arx.

La sintassi è esattamente la stessa, si usa pero’ un altro criterio di minimizzazione, quello

delle variabili strumentali: in sostanza si tenta di imporre che alcune proiezioni della norma

di e(t) vadano a zero.

L’algoritmo è più lento di arx e da risultati leggermente inferiori rispetto a quest’ultimo.

Il problema con questo tipo di algoritmi è che i risultati sono fortemente dipendenti dal tipo

di struttura che si vuole identificare, e la struttura del sistema (che è comunque soggetta a

variazioni durante il funzionamento), non si conosce a priori.

Per sistemi SISO questo problema può essere risolto mediante le funzioni arxstruc e

selstruc:

Formando una matrice nn in modo che ogni sua riga sia del tipo [na nb nk], ed individui

quindi una struttura, la funzione arxstruc calcola per ogni struttura un indice v (loss-

function), rappresentante il valore della funzione minimizzata.

La sintassi è:

V=arxstruc([ye ue],[yv uv],nn);

[ye ue] sono i dati usati per calcolare il modello,

[yv uv] sono i dati usati per calcolare la funzione errore, usando il modello calcolato

precedentemente (cross validation).

V è il vettore contenente i valori di v per ogni struttura.

La funzione selstruc seleziona tra tutte le strutture quella con il minor valore di v:

nn=selstruc(V,0);

se il secondo argomento è ‘AIC’ o ‘FPÈ anzichè 0, si usa un criterio per compensare la

diminuzione di v all’ aumentare del numero di parametri della struttura nel caso in cui si

usino per la validazione gli stessi dati usati per la stima.

106

La funzione spa (spectral analisys) è quella che, nel campo dell’identificazione non

parametrica, da i risultati migliori, ed è inoltre l’unica a poter essere usata per i sistemi

MIMO.

La risposta in frequenza del sistema viene calcolata direttamente nel seguente modo:

G e juu yu( ) ( ) ( )ω ω ω= −Φ Φ1 Eq 7-8

dove Φuu e Φyu sono le stime delle densità spettrali di potenza in ingresso,e mutua uscita

ingresso.

Si puo’ dimostrare che facendo cio’ viene implicitamente minimizzata una certa funzione

della norma di e(t).

I risultati ottenuti sono comunque inferiori di molto a quelli ottenuti con l’identificazione

parametrica.

La funzione forse più adatta all’ identificazione multivariabile (off line) è la funzione pem

(prediction error method).

Il modello è rappresentato nello spazio degli stati:

G q C qI A B

H q C qI A K I

( ) = ( )

( ) = ( ) +

-1

-1

−

−Eq 7-9

e si aggiunge ad y(t) il contributo di uno stato iniziale x(0) :

y t G q u t H q e t C qI A x t( ) = ( ) ( ) + ( ) ( ) + ( ) ( )-1− ( )0 δ Eq 7-10

Il tutto si puo esprimere più sinteticamente così:

[ ] [ ]x t

y t

A B K x

C I

x t

u t

e t

t

( )

( )

( )

( )

( )

( )

( )

=

0

0 0

δ

Eq 7-11

107

La struttura delle matrici è questa:

• Ogni elemento di B, K, x(0) è un parametro da

identificare.

• A è una matrice di 0 ed 1 con ny righe piene di

parametri da identificare, la posizione della prima riga e

le ny-1 distanze tra una di queste righe e la precedente si

dicono indici di pseudo-osservabilità e rappresentano in

sostanza il max numero di ritardi usati per la rispettiva

uscita; la loro somma è uguale al numero di righe di A,

cioè all’ ordine del sistema da identificare.

• C è una matrice di zeri, con degli 1 posizionati in

funzione degli indici di pseudo-osservabilità, e quindi

completamente nota.

Questo tipo di struttura risponde all’esigenza di avere una realizzazione minima di un certo

ordine, con il minimo numero di parametri da stimare, e di essere identificabile.

Per maggiori dettagli vedi Ljung 1987,pag.115-125.

La sintassi del comando pem è la seguente:

th=pem([y u],canstart([y u],[[[[v1 ... vny]]]],nu));

th y u hanno lo stesso significato visto in precedenza,

v1 ... vny, sono gli indici di pseudo osservabilità, nu è il numero di ingressi.

L’algoritmo minimizza la norma dell’errore usando un metodo di tipo Gauss-Newton.

I risultati sono ottimi con indici di pseudo osservabilità tra 8 e 12, tuttavia ciò da luogo ad

un numero di parametri dell’ordine di 20*nu*ny, e quindi ad una generale perdita di

affidabilità all’ aumentare della grandezza del sistema, e comunque ad una calcolo

108

estremamente lento (1h per un sistema 2*3, 486dx33), percui l’algoritmo non è applicabile

on-line, almeno nel nostro caso(sistema 5*5).

Esiste una versione ricorsiva degli algoritmi di identificazione in cui si compie una

iterazione nello stesso tempo in cui una nuova osservazione si aggiunge alle precedenti:

ϑ ϑ( ) ( ) ( ) ( ) ( )t t Q t t e t= −1 + ϕ Eq 7-12

e(t) è l’errore di predizione ϕ(t) è il suo gradiente rispetto a θ.

La scelta di Q(t) condiziona la direzione ed il passo della ricerca, ci sono 3 diverse filosofie

di scelta di Q(t):

1) Metodi del gradiente, in cui Q(t)=γI o Q(t)=(γ/|ϕ(t)|)I.

2) Q è calcolata dalle equazioni del filtro di kalman.

3) Forgetting factor, nella formula si aggiunge un valore λ che porta alla minimizzazione

della seguente funzione:

( )V e e k e ktt k

k

tT( (.)) ( ) ( )= ∑ −

=λ

1Eq 7-13

In pratica le osservazioni vengono pesate in modo da far contare quelle vecchie sempre di

meno, moltiplicandole per λ ad ogni passo, il valore 1/(1-λ) può essere pensato come

l’orizzonte di osservazione.

Questo approccio sembra particolarmente adatto per un uso on-line dove occorre

dimenticare le uscite troppo vecchie.

Esiste una versione ricorsiva di arx e pem; queste funzioni in pratica compiono un singolo

passo della ricerca ogni volta che vengono chiamate, la loro sintassi è la seguente:

[th,yh,P,phi]=rarx([y u], [na nb nk],adm,adg,th0,P0,phi0);

109

[y u] e [na nb nk] hanno il significato visto.

adm,adg specificano il metodo di scelta della Q(t),per esempio, adm=‘ff’, adg=0.99,

specificano un algoritmo di tipo forgetting factor con fattore , λ=0.99.

th0,P0,phi0, contengono rispettivamente i parametri , la loro matrice di covarianza,e i dati

al tempo zero; questi ultimi sono ottenuti con una prima chiamata della funzione:

[th0,yh0,P0,phi0]=rarx([y u], [na nb nk],adm,adg);

th,P,phi sono i valori aggiornati di un passo rispetto a th0,P0 ...

yh è l’uscita predetta al tempo t.

L’identificazione ricorsiva è particolarmente adatta ad essere applicata on-line, e quindi

sfruttata in un controllo adattivo; sfortunatamente però nessuno di tali algoritmi è in grado

di maneggiare sistemi MIMO.

Verranno ora presentate due funzioni arx3, e la sua versione ricorsiva rarx3, che effettuano

una stima sfruttando un modello di tipo arx e selezionando automaticamente la migliore

struttura anche nel caso di sistemi multivariabili.

L’idea di fondo che queste funzioni sfruttano è semplice:

tutte le funzioni di trasferimento di cui la matrice di trasferimento è composta vengono

trattate separatamente,(in parallelo), per ognuna di esse si seleziona la struttura migliore e

poi si estrae un modello; alla fine tutti i modelli vengono messi insieme ottenendo un unico

modello nello spazio degli stati.

Nonostante l’apparente complessità le funzioni, basandosi su algoritmi lineari sono molto

veloci.

110

7.2.1 Funzione arx3

Segue il listato, in linguaggio Matlab, della funzione arx3:

function th=arx3(y,u,mxs)

% th=arx3(y,u,mxs) MIMO model obtained joining all the MISO arx models% between inputs and each output.% Each arx model is chosen with the structure that achieve the smallest AIC,% according to a simple search algorithm, in which the first half of data% is used for estimation, the second for cross validation.% Each column of y,(z) contains the story of the related output,(input).% the stable part of the whole system is balanced and truncated in order% to reduce the states to mxs if possible (mxs=inf skips reduction).

% Giampiero Campa 12-12-95

% initializationno=size(y,2);ni=size(u,2);

f=size(y,1);l=floor(f/2);

% outer cycle (constructs matrices row by row)A=[];B=[];C=[];D=[];K=[];X0=[];for i=1:no,

% computes optimal structuren0=struc(1:3,1:3,0:3);n0=[n0(:,1) n0(:,2)*ones(1,ni) n0(:,3)*ones(1,ni)];n1=sls3(arxstruc([y(1:l,i) u(1:l,:)],[y(l+1:f,i) u(l+1:f,:)],n0),'AIC');vtx=(n1==3);c=0;

while any(vtx)&(c<9),n0=struc(n1(1)-vtx(1):n1(1)+vtx(1),n1(2)-vtx(2):n1(2)+vtx(2),n1(2+ni)-vtx(2+ni):n1(2+ni)+vtx(2+ni));n0=[n0(:,1) n0(:,2)*ones(1,ni) n0(:,3)*ones(1,ni)];n2=sls3(arxstruc([y(1:l,i) u(1:l,:)],[y(l+1:f,i) u(l+1:f,:)],n0),'AIC');vtx=(n2==(n1+1));n1=n2;c=c+1;end

th=arx([y(:,i) u],n1);[a,b,c,d,k,x0]=th2ss(th);

% state space matrices construction

111

A=[A zeros(size(A,1),size(a,2)); zeros(size(a,1),size(A,2)) a];B=[B;b];C=[C zeros(size(C,1),size(c,2)); zeros(size(c,1),size(C,2)) c];D=[D;d];K=[K zeros(size(K,1),size(k,2)); zeros(size(k,1),size(K,2)) k];X0=[X0; x0];

% outer cycleend

% balance and truncation if requiredif ~isinf(mxs),

% system in continuous time domains0=bilinz2s(pck(A,[B K X0(1:size(A,1))],C,[D eye(size(K,2)) zeros(size(D,1),1)]),1);

% extract unstable part[ss0,su0] = sdecomp(s0,-1e-12);

% stable part balanced realization[s6,h6]=sysbal(ss0);[ty6,no6,ni6,ns6]=minfo(s6);if ty6=='cons', s6=[s6,zeros(size(s6,1));zeros(size(s6,2)),-inf]; end

% stable part truncation and reduced system constructions2=madd(strunc(s6,max(0,mxs-su0(1,size(su0,2)))),su0);

% modal form with sorted eigenvalues[ty2,no2,ni2,ns2]=minfo(s2);if ty2=='cons', s2=[s2,zeros(size(s2,1));zeros(size(s2,2)),-inf]; end[T,E]=eig(s2(1:ns2,1:ns2));[l,i]=sort(real(diag(E)));

% final discrete time system% s3=bilins2z(statecc(s2,T(:,i)),1);s3=bilins2z(s2,1);

A=s3(1:ns2,1:ns2);B=s3(1:ns2,ns2+[1:ni]);C=s3(ns2+[1:no],1:ns2);D=s3(ns2+[1:no],ns2+[1:ni]);K=s3(1:ns2,ns2+ni+[1:no]);X0=s3(1:ns2,ns2+ni+no+1);

% end ifend

112

% final theta formatth=ss2th(A,B,C,D,K,X0);

%bodeplot([th2ff(th) spa([y u],30*ones(1,no))]);end

7.2.2 Funzione rarx3

La funzione rarx3 è la versione ricorsiva di arx3, segue il suo listato:

function [th,datan]=rarx3(y,u,mxs,adm,adg,data0)

% [th,datan]=rarx3(y,u,mxs,adm,adg,data0) MIMO model obtained joining% all the MISO recursive arx models between inputs and each output.% Each rarx model is chosen on the first time, with the structure that achieve% the smallest AIC, according to a simple search algorithm, in which the first% half of data is used for estimation, the second for cross validation.% Each column of y,(z) contains the story of the related output,(input).% The stable part of the whole system is balanced and truncated in order% to reduce the states to mxs if possible.% mxs=inf simply skips balance and truncation, it is advisable to skip it% during recursive use, since it requires many computations.% data0 should be zero in the first execution of this function, then,% once the structure is selected, use data0 as the output datan of the last% execution, it contains parameters, covariance matrix, regression vectors, for% each of the transfer functions.% adm: Adaptation mechanism. adg: Adaptation gain% adm='ff', adg=lam: Forgetting factor algorithm, with forg factor lam% adm='kf', adg=R1: The Kalman filter algorithm with R1 as covariance% matrix of the parameter changes per time step% adm='ng', adg=gam: A normalized gradient algorithm, with gain gam% adm='ug', adg=gam: An Unnormalized gradient algorithm with gain gam

% Giampiero Campa 12-12-95

% initializationno=size(y,2);ni=size(u,2);

dt1=[];dt2=[];

fex=size(data0,1);f=size(y,1);l=floor(f/2);cnt=0;

113

% outer cycle (constructs matrices row by row)A=[];B=[];C=[];D=[];K=[];X0=[];for i=1:no,

cnt=cnt+1;

if fex<2, % first time: computes optimal structure

n0=struc(1:3,1:3,1:3);n0=[n0(:,1) n0(:,2)*ones(1,ni) n0(:,3)*ones(1,ni)];n1=sls3(arxstruc([y(1:l,i) u(1:l,:)],[y(l+1:f,i) u(l+1:f,:)],n0),'AIC');vtx=(n1==3);c=0;

while any(vtx)&(c<9),n0=struc(n1(1)-vtx(1):n1(1)+vtx(1),n1(2)-vtx(2):n1(2)+vtx(2),n1(2+ni)-vtx(2+ni):n1(2+ni)+vtx(2+ni));n0=[n0(:,1) n0(:,2)*ones(1,ni) n0(:,3)*ones(1,ni)];n2=sls3(arxstruc([y(1:l,i) u(1:l,:)],[y(l+1:f,i) u(l+1:f,:)],n0),'AIC');vtx=(n2==(n1+1));n1=n2;c=c+1;end

% execution of algorithm for first time[tmn,yh,Pn,Phn]=rarx([y(:,i) u],n1,adm,adg);

else% NOT FIRST TIME, it is all in data0

n1=data0(1,(8+2*ni)*(cnt-1)+[1:1+2*ni]);

la1=data0(1,(8+2*ni)*(cnt-1)+2+2*ni);la2=data0(1,(8+2*ni)*(cnt-1)+3+2*ni);lb1=data0(1,(8+2*ni)*(cnt-1)+4+2*ni);lb2=data0(1,(8+2*ni)*(cnt-1)+5+2*ni);lc1=data0(1,(8+2*ni)*(cnt-1)+6+2*ni);lc2=data0(1,(8+2*ni)*(cnt-1)+7+2*ni);ldt=data0(1,(8+2*ni)*(cnt-1)+8+2*ni);

% construction of matrices from second row of data0tm0=reshape(data0(2,ldt+[1:la1*la2]),la1,la2);P0=reshape(data0(2,ldt+la1*la2+[1:lb1*lb2]),lb1,lb2);Ph0=reshape(data0(2,ldt+la1*la2+lb1*lb2+[1:lc1*lc2]),lc1,lc2);

% execution of recursive algorithm[tmn,yh,Pn,Phn]=rarx([y(:,i) u],n1,adm,adg,tm0(la1,:),P0,Ph0);

114

% end if first timeend

% takes only last row of tmn to avoid data explosion% anyway could be skipped without many consequencestmn=tmn(size(tmn,1),:);

% new data constructionla1=size(tmn,1);la2=size(tmn,2);lb1=size(Pn,1);lb2=size(Pn,2);lc1=size(Phn,1);lc2=size(Phn,2);ldt=size(dt2,2);

dt1=[dt1 n1 la1 la2 lb1 lb2 lc1 lc2 ldt];dt2=[dt2 reshape(tmn,1,la1*la2) reshape(Pn,1,lb1*lb2) reshape(Phn,1,lc1*lc2)];

% puts theta in state space formBB=[];for j=1:niBB=[BB zeros(1,n1(1+ni+j)) tmn(la1,n1(1)+[1:n1(1+j)])];endAA=[1 tmn(la1,1:n1(1))];th=arx2th(AA,BB,1,ni);[a,b,c,d,k,x0]=th2ss(th);

% state space matrices constructionA=[A zeros(size(A,1),size(a,2)); zeros(size(a,1),size(A,2)) a];B=[B;b];C=[C zeros(size(C,1),size(c,2)); zeros(size(c,1),size(C,2)) c];D=[D;d];K=[K zeros(size(K,1),size(k,2)); zeros(size(k,1),size(K,2)) k];X0=[X0; x0];

% outer cycleend

% datan from dt1 and dt2l1=size(dt1,2);l2=size(dt2,2);if l2>l1,datan=[dt1 zeros(1,l2-l1);dt2];elsedatan=[dt1; dt2 zeros(1,l1-l2)];end

115

% balance and truncation if requiredif ~isinf(mxs),

% system in continuous time domains0=bilinz2s(pck(A,[B K X0(1:size(A,1))],C,[D eye(size(K,2)) zeros(size(D,1),1)]),1);

% extract unstable part[ss0,su0] = sdecomp(s0,-1e-12);

% stable part balanced realization[s6,h6]=sysbal(ss0);[ty6,no6,ni6,ns6]=minfo(s6);if ty6=='cons', s6=[s6,zeros(size(s6,1));zeros(size(s6,2)),-inf]; end

% stable part truncation and reduced system constructions2=madd(strunc(s6,max(0,mxs-su0(1,size(su0,2)))),su0);

% modal form with sorted eigenvalues[ty2,no2,ni2,ns2]=minfo(s2);if ty2=='cons', s2=[s2,zeros(size(s2,1));zeros(size(s2,2)),-inf]; end[T,E]=eig(s2(1:ns2,1:ns2));[l,i]=sort(real(diag(E)));

% final discrete time system% s3=bilins2z(statecc(s2,T(:,i)),1);s3=bilins2z(s2,1);

A=s3(1:ns2,1:ns2);B=s3(1:ns2,ns2+[1:ni]);C=s3(ns2+[1:no],1:ns2);D=s3(ns2+[1:no],ns2+[1:ni]);K=s3(1:ns2,ns2+ni+[1:no]);X0=s3(1:ns2,ns2+ni+no+1);

% end ifend

% final theta formatth=ss2th(A,B,C,D,K,X0);

% Comparation with spectral analysis if is the first time or a reduction is required%if (~isinf(mxs) | fex<2),%bodeplot([th2ff(th) spa([y u],30*ones(1,no))]);%end

end

116

7.3 Confronto

Segue un confronto tra tutte le funzioni viste, il modello utilizzato è un sistema con 2 uscite

e 3 ingressi, proposto come esempio nel manuale del System Identification toolbox, pag 1-

32.

Nelle figure dalla 7-2 alla 7-7 vengono mostrati i diagrammi in ampiezza e fase della f.d.t.

tra il primo ingresso e la prima uscita.

Su ciascun diagramma vengono riportati gli andamenti sia della f.d.t. vera che della f.d.t.

identificata.

1) spa:

10−2 10−1 100 10110−4

10−2

100

102

104

frequency (rad/sec)

AMPLITUDE PLOT, input # 1 output # 1

10−2 10−1 100 101−1500

−1000

−500

0

500PHASE PLOT, input # 1 output # 1

frequency (rad/sec)

phas

e

Fig. 7-2

117

2) arx (struttura nota).

10−2 10−1 100 10110−4

10−2

100

102

frequency (rad/sec)


10−2 10−1 100 101−1000

−500

0


frequency (rad/sec)

phas

e

Figure 7-3

3) iv4 (struttura nota).

10−2

10−1

100

101

10−4

10−2

100

102

frequency (rad/sec)


10−2

10−1

100

101

−400

−200

0


frequency (rad/sec)

phas

e

118

4) pem (ps-idx=3):

10−2

10−1

100

101

10−4

10−2

100

102

frequency (rad/sec)


10−2

10−1

100

101

−400

−200

0


frequency (rad/sec)

phas

e

Fig. 7-4

5) pem (ps-idx=10):

10−2

10−1

100

101

10−4

10−2

100

102

frequency (rad/sec)


10−2

10−1

100

101

−600

−400

−200

0


frequency (rad/sec)

phas

e

Fig. 7-5

119

5) arx3:

Fig. 7-6

6) rarx3:

Fig. 7-7

10−2

10−1

100

101

10−4

10−2

100

102

frequency (rad/sec)


10−2

10−1

100

101

−600

−400

−200

0


frequency (rad/sec)

phas

e

10−2

10−1

100

101

10−4

10−2

100

102

frequency (rad/sec)


10−2

10−1

100

101

−400

−200

0


frequency (rad/sec)

phas

e

120

7.4 Funzioni per l’adattamentoVediamo infine la struttura delle due funzioni che realizzano il controllo adattivo.

La seconda, ada2 , è la versione ricorsiva della prima, ed è un po più complessa.

Si da il loro listato in Matlab:

function [ret,x0] = ada1(t,x,u,flag,name,ts,tn,nst,buf)

% [ret,x0] = ada1(t,x,u,flag,name,ts,tn,nst,buf) adaptive control function,% computes a new controller every tn seconds in the scheme 'name'% the first 5 entries of u are supposed to be inputs of system,% the last 5 its outputs, both are sampled every ts seconds.% This function finds a new model based on the story of inputs% and outputs stored in a buffer of size buf ( buf*ts=horizon ).% it execute the function arx3 every tn seconds.% The model can be balanced and truncated% in order to have a whole linear system of nst states,% (if nst==NaN the balance and truncation is disabled).% An H infinity controller with weighting matrices% (depending on vector x) stored in a1 .. d2 is then found% and stored in the block "control".% This controller is in any case "fixed" to have the same% number of states of the initial controller.

global N uequ yequ A B C D Ak Bk Ck Dk Df yd ud

if flag==0,[ni,no]=size(Dk);ret=[0 1 0 ni+no 0 0];x0 = 0;ud=[];yd=[];

elseif abs(flag)==2,ret = x;

elseif abs(flag)==4,ret = t+ts; % next step time.

elseret = [];

end

if (abs(flag)==2) & (abs(round(t/ts) - t/ts) < 1e-8),% step time

% size of the (square) plant

121

[ni,no]=size(Dk);

if x==0,% not enough data yetud=[ud; u(1:ni)'];yd=[yd; u(ni+[1:no])'];if size(ud,1)>buf,x=1;end

else% enough dataud=[ud(2:size(ud,1),:); u(1:ni)'];yd=[yd(2:size(yd,1),:); u(ni+[1:no])'];

% end if x==0 or notend

if (abs(round(t/tn) - t/tn) < 1e-8) & (t>=0.9*tn),% time to find a new controller

% new plantth=arx3(yd,ud,nst);[A,B,C,D,K,Xe]=th2ss(thd2thc(sett(th,ts),'del','NoP'));

% creation of weighting functionsxc=[6.63071609934561 -0.00000001000000 29.24136343802332 0.67538629717924]';[a1,b1,c1,d1]=hfot(-6,xc(2),xc(1),1,no);[a2,b2,c2,d2]=hfot(6,xc(4),xc(3),1,ni);

G=pck(A,B,C,D);s_out=mmult([-eye(no) eye(no)],daug(G,eye(no)));s_kos=abv(eye(no),pck(a1,b1,c1,d1));s_kks=abv(pck(a2,b2,c2,d2),eye(ni));s_kin=mmult(daug(eye(ni),s_out),daug(s_kks,eye(no)));s_fin=mmult(daug(eye(ni),s_kos),s_kin);

plant=mmult([zeros(no,ni+no) eye(no); eye(ni+no) zeros(ni+no,no)],...s_fin,...[zeros(ni,no) eye(ni); eye(no) zeros(no,ni)]);

% new controller computationgmax=1;[K,P,gf]=hinfsyne(plant,no,ni,0,gmax,.01,inf,1);while K==[],gmax=10*gmax;[K,P,gf]=hinfsyne(plant,no,ni,0,gmax,.01,inf,1);endK=nsfix(K,size(Ak,1),1,1);

122

[Ak,Bk,Ck,Dk]=unpck(K);

% [ak,bk,ck,dk]=reg(A,B,C,D,...% lqr(A,B,10^(1.37464278889085)*eye(size(A,1)),eye(size(B,2))),...% lqe(A,eye(size(A,1)),C,10^(-2.39996276978935)*eye(size(A,1)),eye(size(C,1))));% [Ak,Bk,Ck,Dk]=unpck(nsfix(pck(ak,bk,ck,dk),size(Ak,1),1,1));

s_l3=starp(mmult([eye(no);eye(no)],G,K),[-eye(no) eye(no)],no,no);[A3,B3,C3,D3]=unpck(s_l3);mp3=max(real(eig(A3)));mf3=max(abs(eig(A3)));g3=D3-C3*inv(A3)*B3;

% feedforward matrixDf=inv(g3);B3=B3*Df;D3=D3*Df;

% set parameters in the blockset_param([name,'/','Control'],'A',Ak,'B',Bk,'C',Ck,'D',Dk);set_param([name,'/','F-fwd'],'D',Df);

% end controller or notend

ret=x;% end step timeend

123

function [ret,x0] = ada2(t,x,u,flag,name,ts,tn,nst,adm,adg)

% [ret,x0] = ada2(t,x,u,flag,name,ts,tn,nst,adm,adg) adaptive control function,% computes a new controller every tn seconds in the scheme 'name'.% the first 5 entries of u are supposed to be inputs of system,% the last 5 its outputs, both are sampled every ts seconds.% This function finds a new model based on the story of these inputs% and outputs, executing a step of the function rarx3 every ts seconds.% the stable part of this model can be balanced and truncated% in order to have a whole linear system of nst states,% (if nst==NaN the balance and truncation is disabled).% An H infinity controller with weighting matrices% (depending on vector x) stored in a1 .. d2 is then found% and stored in the global variables Ak .. Dk% This controller is in any case "fixed" to have the same% number of states of the initial controller.

global N uequ yequ A B C D Ak Bk Ck Dk Df data yd ud

if flag==0,[ni,no]=size(Dk);ret=[0 1 0 ni+no 0 0];x0 = 0;ud=[];yd=[];

elseif abs(flag)==2,ret = x;

elseif abs(flag)==4,ret = t+ts; % next step time.

elseret = [];

end

if (abs(flag)==2) & (abs(round(t/ts) - t/ts) < 1e-8),% step time

% size of the plant[ni,no]=size(Dk);

if (abs(round(t/tn) - t/tn) < 1e-8) & (t>=0.9*tn),% time to find a new controller

% new plant[th,data]=rarx3(u(ni+[1:no])',u(1:ni)',nst,adm,adg,data);

124

[A,B,C,D,K,Xe]=th2ss(thd2thc(sett(th,ts),'del','NoP'));

% creation of weighting functionsxc=[6.63071609934561 -0.00000001000000 29.24136343802332 0.67538629717924]';[a1,b1,c1,d1]=hfot(-6,xc(2),xc(1),1,no);[a2,b2,c2,d2]=hfot(6,xc(4),xc(3),1,ni);

G=pck(A,B,C,D);s_out=mmult([-eye(no) eye(no)],daug(G,eye(no)));s_kos=abv(eye(no),pck(a1,b1,c1,d1));s_kks=abv(pck(a2,b2,c2,d2),eye(ni));s_kin=mmult(daug(eye(ni),s_out),daug(s_kks,eye(no)));s_fin=mmult(daug(eye(ni),s_kos),s_kin);

plant=mmult([zeros(no,ni+no) eye(no); eye(ni+no) zeros(ni+no,no)],...s_fin,...[zeros(ni,no) eye(ni); eye(no) zeros(no,ni)]);

% new controller computationgmax=1;[K,P,gf]=hinfsyne(plant,no,ni,0,gmax,.01,inf,1);while K==[],gmax=10*gmax;[K,P,gf]=hinfsyne(plant,no,ni,0,gmax,.01,inf,1);end

K=nsfix(K,size(Ak,1),1,1);[Ak,Bk,Ck,Dk]=unpck(K);

% [ak,bk,ck,dk]=reg(A,B,C,D,...% lqr(A,B,10^(1.37464278889085)*eye(size(A,1)),eye(size(B,2))),...% lqe(A,eye(size(A,1)),C,10^(-2.39996276978935)*eye(size(A,1)),eye(size(C,1))));% [Ak,Bk,Ck,Dk]=unpck(nsfix(pck(ak,bk,ck,dk),size(Ak,1),1,1));

s_l3=starp(mmult([eye(no);eye(no)],G,K),[-eye(no) eye(no)],no,no);[A3,B3,C3,D3]=unpck(s_l3);mp3=max(real(eig(A3)));mf3=max(abs(eig(A3)));g3=D3-C3*inv(A3)*B3;

% feedforward matrixDf=inv(g3);B3=B3*Df;D3=D3*Df;

% set parameters in the blockset_param([name,'/','Control'],'A',Ak,'B',Bk,'C',Ck,'D',Dk);set_param([name,'/','F-fwd'],'D',Df);

125

else% not yet time to find a new controller

if x==0,% there is no structure yet

% update ud & ydud=[ud; u(1:ni)'];yd=[yd; u(ni+[1:no])'];

if t>0.9*tn,% first time of execution[th,data]=rarx3(yd,ud,inf,adm,adg,0);ud=[];yd=[];x=1;

end

else% there is the structure% just a step needed[th,data]=rarx3(u(ni+[1:no])',u(1:ni)',inf,adm,adg,data);

% end if x==0 or notend

% end controller or notend

ret=x;% end step timeend

7.5 Risultati

Purtroppo il controllo adatttivo non da i risultati sperati in quanto l’identificazione in

closed-loop non è altrettanto buona, con nessuno dei metodi visti sopra, per cui succede che

a volte, specialmente in presenza di segnali di riferimento costanti, il sistema identificato è

diverso da quello vero, e ed il controllore non riesce a stabilizzare l’impianto.

L’algoritmo di identificazione che funziona meglio, pem, per un sistema 5*5 occupa

rapidamente tutta la memoria del sistema, bloccando la simulazione.

126

La combinazione tra identificazione pem e controllo adattivo H∞ potrebbe quindi

funzionare su un sistema con maggiori risorse di memoria RAM tuttavia essa sarebbe di

gran lunga troppo complessa per essere usata in un controllo on-line.

Risulta chiaro quindi che l’unico modo di affrontare la non-linearità dell’ impianto è quello

di sfruttare il più possibile la conoscenza fisica dell’ impianto stesso.

Ciò può essere fatto sia con un controllo gain-scheduling, o ancora con un controllo

adattivo.

Tuttavia, volendo seguire la strada del controllo adattivo, bisognerà provvedere ad un

meccanismo di identificazione che, identificando alcuni parametri fisici del sistema sia in

grado, sfruttando la conoscenza fisica del veicolo, di costruirne un modello affidabile, in

modo da permettere la realizzazione di un controllore idoneo.

127

8. Conclusioni

L’ obiettivo di dotare il veicolo sottomarino trainato (TUV) di un sistema di controllo atto a

consentire stabilità dell’ assetto di navigazione ed un accurato inseguimento delle traiettorie

di riferimento, è stato raggiunto tramite un controllore sviluppato con il metodo H∞.

Il range del campo operativo in cui il controllore può essere adoperato risulta esteso fino a

spostamenti di 15 metri e rotazioni di 15 gradi, ed è quindi abbastanza ampio da

comprendere la maggior parte delle situazioni che normalmente si verificano durante la

navigazione.

Oltre tale range le prestazioni del sistema di controllo decadono, od occorre provvedere ad

un controllo di tipo gain scheduling o self-tuning.

Il risultato suddetto è stato raggiunto tramite i seguenti passi:

• Studio del modello non lineare di un veicolo subacqueo

trainato.

• Scelta degli insiemi di ingresso e uscita appropriati.

• Linearizzazione del modello ed analisi del modello

linearizzato, con particolare attenzione per l’

autostruttura ed i vettori singolari.

• Riduzione del modello con ulteriori analisi del modello

ridotto.

• Messa a punto di un algoritmo per la selezione

automatica e l’ottimizzazione delle funzioni di peso per

un controllore H∞.

• Uso dell’ algoritmo suddetto per l’ottenimento di un

controllore.

128

• Verifica della validità del suddetto controllore sull’

impianto non lineare, e definizione del suo range di

validità.

129

Inoltre, è stata verificata la migliore performance del controllo H∞ rispetto a quello LQG,

mettendo in rilievo le analogie e le differenze tra i due metodi.

Un ulteriore miglioramento del sistema di controllo potrebbe ottenersi rendendo adattivo il

sistema di controllo trovato.

Tale possibilità è stata studiata, e sono stati scritti due di algoritmi di natura parallela, per

l’identificazione multivariabile.

Tali algoritmi, pur essendo confrontabili con quelli presenti nel System Identification

Toolbox, non hanno tuttavia dato risultati soddisfacenti ai fini del controllo adattivo, percui

sono stati indicati altri possibili metodi per realizzare tale tipo di controllo, migliorando

quello già trovato.

130

Bibliografia:

[i] F.Nasuti, “Sviluppo di un simulatore per la dinamica di veicoli sottomarini trainati”, Electrical Engineering Thesis,Apr-1994. University of Pisa.

[ii] J.M.Maciejowski, “Multivariable Feedback Design”, Addison-Wesley 1989.[iii] J.C.Doyle, K.Glover, A.Packard, R.Smith, G.Balas, “µ-analysis and synthesis toolbox for use with MATLAB”, Jan

1994 The Math Works Inc.

[iv] J.C.Doyle, K.Glover, P.Khargonekar, B.Francis “State space solutions for standard H2 and H∞ control problems”IEEE transactions on automatic control, Aug 1989.

[v] J.C.Doyle, K.Glover, “State space formulae for all stabilizing controllers that satisfy an H∞ norm bound and relactionto risk sensitivity” system and control letters, vol 11, 1988.

[vi] L.Ljung, “System Identification: Theory for the User” Prentice Hall 1987.

[vii] L.Ljung “System Identification Toolbox for use with Matlab” July 1991 The Math Works Inc.

Un sentito ringraziamento per l’assistenza prestata ... Un sentito ringraziamento per...

Documents

Transcript of Un sentito ringraziamento per l’assistenza prestata ... Un sentito ringraziamento per...