CAPITOLO 1 . IL PROBLEMA TRATTATO
4
istante influenza le uscite successive) definito a tempo discreto che può essere rappresentato da
un’equazione alle differenze del tipo:
y((k + n)T) + a
n−1
y((k + n − 1)T) + · · · + a
1
y((k + 1)T) + a
0
y(kT) =
b
m
u((k + m)T) + b
m−1
u((k + m − 1)T) + · · · + b
1
u((k + 1)T) + b
0
u(kT)
1.2 Strumenti utilizzati
Il principale strumento utilizzato per l’analisi di sistemi dinamici ed anche per l’identificazione
dei modelli è il software Matlab®, un programma interattivo basato sui sistemi matriciali, per
applicazioni computazionali scientifiche di tipo numerico, gestibile principalmente da linea di
comando. Per quanto riguarda quest’ultimo campo è stata sviluppata una particolare
applicazione grafica detta GUI, Graphical User Interface, contenuta nel SIT, System
Identification Toolbox.
Digitando da linea di comando “ident” si accede comodamente alla GUI del SIT e da qui si
parte per importare ed analizzare i dati che si vuole modellizzare.
Per quanto riguarda invece l’analisi vera e propria dei sistemi dinamici e la progettazione di reti
correttrici/regolatori si utilizzano altre due principali applicazioni sviluppate sempre per
Matlab®: la prima, TFI, Transfer Function Interpreter sviluppato da A. Civolani e G. Marro è
un tool creato per gestire modelli generici, lineari e multivariabile; la seconda, LTIVIEW, un
altro tool creato per vedere il comportamento dinamico dei sistemi con vari ingressi (gradino o
impulso) e i diagrammi poli-zeri. Successivamente tratterò più in dettaglio il loro utilizzo.
1.3 Modelli utilizzati
Come già accennato, lo scopo dell’identificazione è trovare il modello che riproduce
ottimamente sequenze di campioni di uscita se soggetto a particolari ingressi.
Prima di tutto definiamo il concetto di modello: è una rappresentazione matematica del sistema
reale che pur non descrivendo esattamente il suo comportamento, permette di formulare in
maniera matematica il problema di controllo. La complessità del modello dipende dal sistema
che si deve descrivere e dall’uso che se ne intende fare.
CAPITOLO 1 . IL PROBLEMA TRATTATO
5
Esistono due tipi di modelli nell’identificazione: quelli parametrici e quelli non parametrici; i
primi consistono in un tentativo di stima legato ad specifico modello definito dall’utente,
mentre i secondi si riferiscono alla stima relativa a segnali generici come il segnale “gradino” o
quello “impulso”.
Solitamente l’analisi parte con una stima non parametrica per valutare il comportamento
generale del sistema. Poi si passa ai modelli parametrici iniziando dalla famiglia ARX,
acronimo di AutoRegressive with eXogeneous input, che significa un modello composto da una
parte autoregressiva ed una esogena, come vedremo nel prossimo paragrafo.
1.3.1 Modelli ARX
Questo modello è caratterizzato da un errore d'equazione e(t), rappresentato da un processo
stocastico bianco con valore medio nullo, e da un ordine n che ne indica la memoria:
y(t)= -α
n
y(t-1) + … -y(t-n) + β
n
u (t-1) + …… + β
1
u(t-n) + e(t)
dove i parametri α sono relativi alla memoria delle uscite y(...), mentre i β si riferiscono agli
ingressi u(...). Questi parametri sono rappresentabili in un singolo vettore colonna, chiamato
appunto vettore dei parametri, e definito come:
θ = [α
1
… α
n
| β
1
… β
n
]
Il predittore ARX y(t | t-1), detto predittore ad un passo in avanti, è ottimo se l'errore di
previsione
ε (t,θ) = y(t) - ŷ(t|t-1,θ)
eguaglia l'errore d'equazione
ε (t) = y(t) – y
m
(t).
Il nostro scopo è quello di minimizzare la funzione costo:
ƒ
N
t
t
N
nJ
1
2
),(
1
),( Τ Η Τ
costruita come la somma dei quadrati degli errori di previsione sui dati; il metodo usato
viene chiamato stima ai minimi quadrati (LS – Least Square) e si applica definendo un
vettore dei parametri stimati costruito imponendo:
θ
0
= H
T
y
0
= [( H
T
H )
-1
H
T
] y
0
CAPITOLO 1 . IL PROBLEMA TRATTATO
6
dove è indicata con H una particolare matrice, chiamata matrice di Hankel, che è costruita a
partire da due sottomatrici formate dai valori degli y e u .
La matrice θ
0
esiste ed è unica e minimizza la funzione costo se H ha rango massimo.
Indicando con θ* il valore vero dei dati generati ed N = m - n, in cui m è il numero di
campioni a disposizione, mentre n è l’ordine del modello, la stima LS risulta consistente se:
*0
lim
Τ Τ
φ οN
oppure sotto la condizione che l'errore d'equazione abbia proprietà di bianchezza ed
ergodicità, ed è una stima efficiente se l'errore d'equazione è un processo stocastico
gaussiano.
1.3.2 Modelli ARMAX
Il modello ARMAX è anch'esso appartenente alla famiglia ad equazione d'errore, ma in
questo caso e(t) è descritto tramite un processo MA (Mobile Average) a media mobile, in
questo modo:
e(t) = ω(t) + y
ny
ω(t-1) +… + y
1
ω(t-n)
dove ω (...) è modellato come un processo stocastico a valor medio nullo (µ =0) in cui
l’ordine n di γ
n
è permesso assumerlo pari a quello della parte autoregressiva.
Il modello ARMAX completo ha quindi la seguente struttura:
y(t) = -α
n
y(t-1) + …- α
1
y(t- n) + β
n
u(t-1) + … + β
1
u(t- n) + γ
n
ω(t-1) + … + γ
1
ω(t- n)
Il predittore ottimo per l'ARMAX è individuabile semplicemente osservando che l'unico
termine non predicibile all'istante t-1 è ω(t), da cui
y(t) - y(t|t-1)=ω(t).
Purtroppo non è più applicabile il metodo di stima LS, poiché l'errore d'equazione non è più
bianco e porterebbe ad una stima non polarizzata.
Quindi viene introdotto il metodo della variabile strumentale (IV), che sfrutta una matrice Z
delle stesse dimensioni di H, scritta come:
Z=[ H
n
(η) H
1
(u)]
ma con valori incorrelati da e(t).
CAPITOLO 1 . IL PROBLEMA TRATTATO
7
Le componenti di questa matrice ( strumenti) , possono essere ricavate in due modi:
ξ applicando una stima LS preliminare ottengo gli α
i
e β
i
che utilizzo per creare
η(t) = α
n
y(t-1) + … + α
1
y(t- n) + β
n
u(t-1) + … + β
1
u(t- n)
successivamente si trovano altri parametri α
i
e β
i
estrapolati da un θ
0
che è stato creato
utilizzando la matrice Z e si prosegue per passi successivi fino a giungere a convergenza;
ξ utilizzando come strumenti gli ingressi ai passi precedenti,posso scrivere η(t)=u(t-n), da
qui procedo riempiendo la sottomatrice H
n
( η ). Per determinare invece i valori di y
i
posso
utilizzare:
e(t) = y(t) - α
n
y(t-1) - … - α
1
y(t- n) - β
n
u(t-1) - … - β
1
u(t- n)
e stimando i parametri con il modello MA..
1.3.3 Rappresentazione polinomiale del modello Arx
Quando si opera con Matlab, i modelli vengono espressi non in forma estesa ma in una
particolare rappresentazione polinomiale più compatta; ad esempio il modello Arx descritto nel
seguente modo:
y(t)= α
n
y(t-1) + … + y(t-n) + β
n
u (t-1) + …… + β
1
u(t-n) + e(t)
può essere riscritto anche come:
y(t)- α
1
y(t-1) - … -α
na
y(t-n
a
) = β
1
u (t-n
k
) +… + β
nb
u(t-nk-nb+1) + e(t)
e più esplicitamente:
A(q) y(t) =B(q) u(t-n
k
) + e(t)
dove B ed A sono polinomi dell’operatore ritardo q
-1
, espressi così:
CAPITOLO 1 . IL PROBLEMA TRATTATO
8
A(q) = 1+ α
1
q
-1
+ … + α
na
q
-na
B(q) = β
1
+ β
2
q
-1
+ … + β
nb
q
-nb+1
I numeri n
a
e n
b
sono gli ordini dei rispettivi polinomi, quindi del modello. Il numero n
k
è il
ritardo fra ingresso e uscita.
I parametri da modificare per trovare il modello migliore sono perciò proprio
kba
nnn , ,: per
evitare di vagliare tutte le possibili combinazioni della terna
kba
nnn , , si ricorre alla
“valutazione dell’ordine a priori” che restituisce un valore indicativo di n nell’intorno del quale
limitare l’escursione dei vari
ba
nn , , nonché del ritardo
k
n fra ingresso e uscita. Nel prossimo
capitolo affronteremo meglio questo metodo, mentre nel prossimo paragrafo tratteremo con
precisione i criteri per la determinazione dell’ordine di un modello dal punto di vista teorico.
1.3.4 Struttura e rappresentazione polinomiale di altri modelli
Partendo dallo sviluppo in schema a blocchi del modello Arx:
sviluppiamo ora la struttura di un modello Armax :
A(q) y(t) = B(q) u(t-nk) + C(q) e(t)
dove A(q) e B(q) sono gli stessi definiti per il modello Arx, mentre:
C(q) = 1+ c
1
q
-1
+ …+ c
n c
q
-n c
y[ k ]
e[ k ]
B(q)
)(
1
qA
u[ k ]
][
)(
1
][
)(
)(
][ ke
qA
ku
qA
qB
ky
ARX
CAPITOLO 1 . IL PROBLEMA TRATTATO
9
quindi graficamente:
Un altro modello chiamato Box & Jenkins è invece ottenuto in questo modo:
con
F(q) = 1+ f
1
q
-1
+ …+ f
nf
q
-nf
D(q) = 1+ d
1
q
-1
+ …+ d
nd
q
-nd
che risulta essere la più complessa fra le strutture di modelli analizzati.
Invece la struttura più semplice realizzabile è quella del modello Output Error, data da:
y[ k ]
B(q)
)(
1
qA
u[ k ]
C(q)
e[ k ]
][
)(
)(
][
)(
)(
][ ke
qA
qC
ku
qA
qB
ky
ARMAX
e[ k ]
y[ k ]
)(
)(
qF
qB
)(
)(
qD
qC
u[ k ]
][
)(
)(
][
)(
)(
][ ke
qD
qC
ku
qF
qB
ky
BOX & JENKINS
][][
)(
)(
][ keku
qF
qB
ky
e[ k ]
y[ k ]
)(
)(
qF
qB
u[ k ]
OUPUT ERROR
CAPITOLO 1 . IL PROBLEMA TRATTATO
10
Generalizzando possiamo trovare la struttura di modelli generali parametrici:
][
)(
)(
][
)(
)(
][)( ke
qD
qC
ku
qF
qB
kyqA
L’ultima struttura che andremo ad analizzare sarà quella del modello Spazio degli Stati, che si
presenta così:
]1[]1[]1[
][][]1[
kDukCxky
kBukAxkx
Il System Identification Toolbox (SIT) di Matlab supporta due tipi di parametrizzazioni di
modelli nello Spazio degli Stati, ma a noi interessa il metodo black-box, nel quale
l’indice più importante della struttura è l’ordine del modello n, cioè la dimensione del
vettore degli stati x(t).
Solitamente si assume che per sistemi fisici si considerano modelli puramente dinamici.
Un’ulteriore stima può essere fatta attorno alla matrice di guadagno H : rappresenta il
guadagno dell’osservatore dinamico, in assenza di disturbi, o il filtro di Kalman in
ambiente stocastico.
Il SIT ci fornisce due metodi di base per la stima del modello nello Spazio degli Stati:
ξ Pem: è metodo standard basato sulla minimizzazione iterativa del criterio di
predizione d’errore. Le iterazioni partono da valori di parametri che sono calcolati
dal n4sid , una funzione del SIT. La parametrizzazione delle matrici A, B, C ed H
è effettuata in questo modo: considerati N campioni dei vettori ingresso-uscita
u(t) e y(t), si ricerca il minimo di
)(
1
1
2
te
N
J
N
t
ƒ
CAPITOLO 1 . IL PROBLEMA TRATTATO
11
ξ N4sid: è un metodo nel quale viene selezionato il modello migliore in termini di
predizione )(ty
)
delle uscite misurate y(t) (dato che e(t)=y(t) - )(ty
)
). Questo
metodo non usa la ricerca iterativa, ma la stima può risultare ugualmente di
qualità, a seconda di alcune opzioni scelte nella finestra di comando del SIT.
1.4 Determinazione dell’ordine dei modelli
Per la determinazione di un ordine appropriato del modello, è necessario implementare i
seguenti criteri per vari ordini k ed effettuare una scelta oculata tra i risultati dei diversi
criteri, basati su considerazioni fisiche o statistiche dei dati:
- PPCRE: è definito come un rapporto segnale/rumore (S/N) , precisamente come 100
volte il rapporto tra la deviazione standard dell'errore di previsione e la deviazione
standard dell'uscita. Il criterio viene calcolato per i vari potenziali ordini e viene
selezionato come esatto il minimo ordine per il quale per gli ordini immediatamente
successivi non avviene un decremento sensibile del PPCRE.
]det[
det
110100)(
2
2
k
T
kk
T
k
k
y
k
HHyy
S
kPPCRE
Ρ Ρ
Η
ς
ς
- FPE: consiste nella minimizzazione del valore atteso della varianza dei residui 2ε σ . Si
sceglie come esatto l'ordine associato al minimo valore del criterio.
2
1
0
)(
)(
)()( t
dNN
dN
J
dN
dN
kFPE
L
kt
N ƒ
Η Τ
- AIC: (AKAIKE) è un criterio che tiene conto anche della penalizzazione dei modelli di
ordine elevato; è asintoticamente equivalente all' FPE. Se N è elevato, questo criterio
tende a selezionare modelli di ordine più alto.
dJNkAIC
N
2)](log[)(
0
Τ
CAPITOLO 1 . IL PROBLEMA TRATTATO
12
- MDL: al contrario dei precedenti, che sono basati su considerazioni statistiche, questo
criterio è basato sulla minimizzazione delle informazioni necessarie a descrivere il
modello e il suo errore di previsione.
)](log[]log[
0
N
JNdNMDL Τ
Tutto ciò si tradurrà in questo studio nella stima attraverso SIT di determinati parametri, a
seconda dei modelli che si sta analizzando; più precisamente, per un Arx:
ξ n
a
relativo a A(q) = 1+ α
1
q
-1
+ … + α
na
q
-na
ξ n
b
relativo a B(q) = β
1
+ β
2
q
-1
+ … + β
nb
q
-nb+1
ξ n
k
, ritardo fra ingresso ed uscita.
Per un modello Armax:
ξ n
a
relativo a A(q) = 1+ α
1
q
-1
+ … + α
na
q
-na
ξ n
b
, relativo a B(q) = β
1
+ β
2
q
-1
+ … + β
nb
q
-nb+1
ξ n
c
, relativo a C(q) = 1+ c
1
q
-1
+ …+ c
n c
q
-n c
ξ n
k
, ritardo fra ingresso ed uscita.
Per un modello Box & Jenkins:
ξ n
b
, relativo a B(q) = β
1
+ β
2
q
-1
+ … + β
nb
q
-nb+1
ξ n
c
, relativo a C(q) = 1+ c
1
q
-1
+ …+ c
n c
q
-n c
ξ n
d
, relativo a F(q) = 1+ f
1
q
-1
+ …+ f
nf
q
-nf
ξ n
f
, relativo a D(q) = 1+ d
1
q
-1
+ …+ d
nd
q
-nd
ξ n
k
, ritardo fra ingresso ed uscita.
Per un modello Output Error, come per l’Arx:
ξ n
a
relativo a A(q) = 1+ α
1
q
-1
+ … + α
na
q
-na
ξ n
b
, relativo a B(q) = β
1
+ β
2
q
-1
+ … + β
nb
q
-nb+1
ξ n
k
, ritardo fra ingresso ed uscita.
Ed infine per un modello Stato degli Spazi:
ξ n , ordine generale del modello
ξ n
k
, ritardo fra ingresso ed uscita.
CAPITOLO 1 . IL PROBLEMA TRATTATO
13
1.4.1 Test di bianchezza dei residui
Sapendo che l'errore di previsione è
)()()(
0
Η Η ytyt
dove y° rappresenta il vettore delle uscite osservate del sistema, per un modello identificato
con ordine corretto ε(t) deve essere un processo bianco, ed è buona norma verificare questa
condizione effettuando un controllo sulla covarianza:
)()(
1
)(
1
Ω Η Η Ω
Η
ƒ
tt
N
r
N
t
N
con τ=1, … , M con M solitamente fissato a 8.
Se ε (t) è un rumore bianco allora la quantità:
2
1
2
,
)]([
)0(
1
Ω [
Η
Ω Η
N
M
MN
r
r
ƒ
converge asintoticamente ad una distribuzione di
)(
2
M Φ
Perciò la bianchezza dei residui è solitamente verificata nelle condizioni di:
)8(
2
99,08,
Φ [
N
In pratica si traduce nel fatto che le funzioni di autocorrelazione:
)()(
1
)(
ˆ
1
Ω Η Η Ω
Η
ƒ
tt
N
R
N
t
e cross-correlazione:
)()(
1
)(
ˆ
1
Ω Η Ω
Η
ƒ
ttu
N
R
N
t
u
rimangano entro i limiti del 95 % attorno allo 0, cosa facilmente visualizzabile con il SIT.
CAPITOLO 2. L’IDENTIFICAZIONE DEI MODELLI
14
Capitolo 2
L’identificazione dei modelli
Il particolare progetto che analizzeremo prevede la trattazione di un sistema SISO, singolo
ingresso e singola uscita, a dati campionati a distanza di un giorno; l’analisi si dividerà
principalmente in due aree tematiche principali:
ξ in simulazione ,
ξ in previsione ,
e, a differenza della teoria, sarà composta da solo quattro fasi principali (in quanto non c’è stato
bisogno di raccogliere i dati, avendoli già a disposizione):
a) elaborazione preliminare dei dati
b) scelta dei modelli da analizzare
c) stima del miglior modello
d) validazione e verifica.
Anticipo subito che analizzerò modello per modello poi confronterò i risultati ottenuti.
2.1 Elaborazione preliminare dei dati
Per sfruttare al meglio tutti le informazioni ricavabili dai dati in modo iterativo e semplice ho
costruito una funzione ( “data.m”, fruibile al capitolo 5) che in modo automatico fornisce una
serie di figure necessarie all’analisi preliminare.
La funzione innanzitutto carica i campioni di dati dal file “data_ac.txt”, poi li trasforma e li
rende utilizzabili dall’interfaccia “ident”, cioè in formato “iddata” .
Come consuetudine nel campo dell’identificazione ho diviso la serie temporale di dati in mio
possesso in due parti, una prima per la modellizzazione e le valutazioni, la seconda per la
validazione dei modelli trovati. Cioè le nostre coppie ingresso-uscita sono vettori di 2217
elementi; identificheremo i modelli a partire dai primi 1107 valori e utilizzeremo i restanti per
la validazione. Questo punto consisterà nel testare le capacità predittive del modello sulla
porzione dei dati cui non si è ricorsi in fase di identificazione.
Dopo aver sottoposto i dati alla sottrazione del valor medio gli ho visualizzati, ho operato
l’analisi spettrale dell’uscita riferita all’ingresso , dell’ingresso riferita all’uscita , il diagramma
di Bode ed infine la risposta all’impulso.
Ricevi informazioni sui nostri servizi, sulle offerte e non perdere news e consigli su università e lavoro.