UNIVERSITÀ DEGLI STUDI DI PADOVA

FACOLTÀ DI INGEGNERIA —

LAUREA MAGISTRALE IN INGEGNERIA DELL’AUTOMAZIONE Progettazione di sistemi di controllo

IDENTIFICAZIONE DI MODELLI

TERMODINAMICI

Prof. Schenato Luca Prof. Chiuso Alessando

ALBERTON RICCARDO, AUSSERER MARKUS, BARAZZUOL ANDREA

ANNO ACCADEMICO 2010-2011

Indice

Introduzione  IV
0.1 Pem: Prediction Error Methods . . . . . . . . . . . . . . . . . . . 1
0.2 Definizione del problema . . . . . . . . . . . . . . . . . . . . . . . 4
0.3 Predizione in RKHS e approssimazione finito-dimensionale …. 7
0.4 Stabil Spline Group LAR (SSGLAR) algoritm . . . . . . . . . . . 8
0.4.1 Algoritmo SS-Glar . . . . . . . . . . . . . . . . . . . . . . 9
0.4.2 Stima dei hyperparametri . . . . . . . . . . . . . . . . . . 9
0.5 Stable Spline with Exponential Hyperparametri (SSEH) Algoritm 10
0.5.1 Hyperprobabilità a priori per i hyperparametri . . . . . . . 10
0.5.2 Il modello Bayesiano intero . . . . . . . . . . . . . . . . . 10
0.5.3 Stima degli hyperparametri . . . . . . . . . . . . . . . . . 11
0.5.4 Stima delle risposte impulsive del predittore per ζ noto . . 13
0.6 Validazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
Conclusioni  18
Bibliografia  19

Introduzione

IV

0.1 Pem: Prediction Error Methods

In questo capitolo viene spiegato brevemente il principio di funzionamento dei metodi a minimizzazione dell’errore di predizione. Dato un modello M(θ) appartenente a una certa classe parametrica M ≡ {M(θ); θ ∈ Θ} e dato una sequenza di N misure dei segnali d’ingresso e d’uscita

u^{N} := {u(t); t = 1, …, N} y^{N} := {y(t); t = 1, …, N}

si procede nel modo seguente:

1. Per un qualche valore di θ fissato si costruisce il miglior (secondo qualche criterio) predittore all’istante t − 1 dell’uscita y(t). Questo predittore è una funzione dei dati passati che viene denotata col simbolo M̂(θ):

M̂(θ) : (y^{t−1}, u^{t−1}) → ŷ_{θ}(t|t − 1).

La predizione ŷ_{θ}(t|t − 1) si può pensare come funzione dei dati passati e quindi come una funzione aleatoria che viene indicata usando il simbolo in grassetto ŷ_{θ}(t|t − 1).

2. Si formano gli errori di predizione

ε_{θ}(t) := y(t) − ŷ_{θ}(t) per t = 1, …, N.

Analogamente a quanto detto nel caso del predittore anche gli errori di predizione possono essere visti come quantità aleatorie e vengono quindi indicate col simbolo in grassetto ε_{θ}(t).

3. Si minimizza rispetto a θ una cifra di merito che descrive quanto bene il modello predice il dato successivo. Ad esempio si minimizza l’errore quadratico medio di predizione

V_{N}(θ) := \frac{1}{N} ∑_{t=1}^{N} ε_{θ}(t)^{2}

o più in generale una media degli errori quadratici di predizione pesati con una cifra di merito non negativa β(N, t) > 0,

V_{N}(θ) := \frac{1}{N} ∑_{t=1}^{N} β(N, t)ε_{θ}(t)^{2}

che per N piccoli dà peso minore agli errori di predizione compiuti nella fase iniziale dell’algoritmo quando l’influenza di condizioni iniziali stimate in modo approssimativo è più deleteria. Per N grande β(N, t) tende a essere inutile e quindi deve tendere a uno. In ogni caso si ricava la stima si θ dalla minimizzazione della cifra di merito

θ̂_{N} := Arg min_{θ} V_{N}(θ)

che è la stima PEM del parametro del modello. Lo stimatore θ̂_{N} viene chiamato stimatore PEM del parametro θ.

4. Infine si prende come stima della varianza dell’innovazione λ^{2} = var{e(t)} l’errore quadratico residuo

λ̂^{2} := V_{N}(θ̂_{N}).

Per quanto questa procedura sembra intuitivamente sensata, l’unica giustificazione per la sua adozione nei procedimenti di identificazione sta nelle sue proprietà statistiche. Per avere una trattazione completa delle proprietà statistiche dello stimatore PEM si consulti [1].

Finora si è visto come ricavare lo stimatore PEM nel caso si ha un parametro θ fissato. Poiché in generale però non si ha alcuna informazioni sulla struttura del sistema vero serve uno strumento per scegliere la classe parametrica M alla quale appartiene il modello identificato M(θ). Dato un set di misure e diverse classi di modelli con ordini p plausibili si calcolano degli indici che danno la possibilità di scegliere la classe che promette di ottenere un modello identificato con buone prestazioni. Questi indici sono sostanzialmente una generalizzazione del criterio FPE (Final Prediction Error, vedere [1] per la trattazione completa). Nelle simulazioni svolte per questo elaborato sono stati usati l’AICC (Akaike Information Criterion with a second-order Correction) e il BIC (Bayesian Information Criterion) definito dalle seguenti formule:

AICC = N · ln \frac{rss}{N} + N · \frac{N + p}{N − p − 2}
BIC = N · ln \frac{rss}{N} + p · ln(N)

dove rss (Residual Sum of Squares) è definito da rss = ∑_{t=1}^{N} ε_{θ}(t)^{2}, mentre N è il numero di misure e p è l’ordine del parametro θ del modello. Nel progetto in esame è stata implementata una routine nella quale viene considerato la classe dei modelli ARMAX. In un ciclo ad ogni passo viene aumentato l’ordine del parametro θ e vengono calcolati i valori dei due indici. Terminato il ciclo si controlla in relazione a quale ordine gli indici hanno raggiunto i valori minimi. Viene quindi rifatta l’identificazione PEM usando i rispettivi ordini ritenuti ottimi ottenendo due modelli che spesso si verificano aver ordine differente. Si noti però che per i due criteri considerati all’aumentare di N la probabilità che la scelta dell’ordine sia sbagliato tende a zero (e quindi si dice che la selezione dell’ordine è consistente).

0.2 Definizione del problema

Sia {z_{t}}_{t∈ℤ} z_{t} ∈ ℝ^{m} un processo stocastico, stazionario che descrive l’evoluzione temporale di qualche variabile di interesse. Con un piccolo abuso di notazione in questo contesto z_{t} indica sia la variabile aleatoria che la sua realizzazione. Come mostrato in figure 1 ogni componente del processo {z_{t}}_{t∈ℤ} può essere pensato come un nodo di una grafo. L’obiettivo di questi algoritmi di identificazione è di trovare modelli lineari dinamici che descrivono tutti i componenti di {z_{t}} come funzione degli altri componenti. A questo scopo si definisce y_{t} := z_{t}^{[i]} (la componente i-esima di z_{t}) come uscita e tutte le altre componenti u_{t} := z_{t}^{[−i]} ∈ ℝ^{m−1} come ingressi. Ovviamente il procedimento può essere ripetuto per ogni i = 1, …, m ottenendo in questo modo una descrizione completa di tutte le variabili di z_{t} come funzione degli altri. In questa sezione però per questioni di chiarezza si sceglie un valore specifica per i, per esempio i = 1. In questo caso z_{t} può essere rappresentato come

z_{t} = [\begin{matrix}y_{t} \\ u_{t} \end{matrix}]

Si definiscono gli insiemi delle variabili passate all’istante t

Y^{t} = [y_{t−1}, yt − 2, …] U^{t} = [u_{t−1}, ut − 2, …]

Dalla stazionarietà del processo {z_{t}}_{t∈ℤ} segue che anche i processi {y_{t}}_{t∈ℤ} e {u_{t}}_{t∈ℤ} sono processi stocastici, congiuntamente stazionari che possono quindi essere pensati rispettivamente come ingressi e uscita di un sistema dinamico tempo-invariante

y_{t} = ∑_{k=1}^{∞} f_{k}u_{t−k} + ∑_{k=0}^{∞} g_{k}e_{t−k}  (1)

dove f_{k} ∈ ℝ^{1×m} (avendo assunto assenza di retroazione diretta, cioè f_{0} = 0) e g_{k} ∈ ℝ sono coefficienti incogniti della risposta impulsiva e e_{t} rappresenta un processo di innovazione che è equivalente all’errore di predizione di un passo

e_{t} := y_{t} − ŷ_{t|t−1} = y_{t} − E[y_{t}|Y^{t}, U^{t}]
E[y_{t}|Y^{t}, U^{t}] := ∑_{j=1}^{m−1} [∑_{k=1}^{∞} h_{k}^{j}u_{t−k}^{[j]}] + ∑_{k=1}^{∞} h_{k}^{[m]}y_{t−k}.

Figura 1: Grafo che rappresenta l’interconnessione di m = 6 variabili. Gli archi solidi rappresentano le connessioni nel modello dinamico con il nodo y_{t} := z_{t}^{[1]} da tutti gli altri nodi. L’assenza degli archi da u − t^{[i]}, i = 1, 5 a y_{t} indica che h_{k}^{[i]} = 0, i = 1, 5, ∀k ∈ ℤ^{+}. Il nodo y_{t} ha l’arco entrante indicando la dipendenza dall’errore di predizione a un passo e_{t}. Gli archi punteggiati indicano dipendenze tra gli ingressi u_{t}.

Le sequenze h_{k} := [h_{k}^{[1]}, h_{k}^{[2]}, …, h_{k}^{[m]}] ∈ ℝ^{1×m}, k ∈ ℤ^{+} rappresentano i coefficienti della risposta impulsiva del predittore e servono a descrivere sistemi (BIBO) stabili e si ha quindi h^{[m]} ∈ l_{1}(ℤ^{+}).

Nell’ambito della PEM la predizione di sistemi dinamici può essere vista come predizione dei coefficienti della risposta impulsiva h_{k} del predittore. Si vuole prestare particolare attenzione a situazioni in cui non tutte le variabili sono necessari per predire y_{t}, il che in termini matematici equivale ad avere h_{k}^{[i]} = 0 per ogni k ∈ ℤ^{+}. In una rappresentazione grafica c’è una arco diretto u_{k}^{[i]} al nodo y_{k} se e solo se esiste k ∈ ℤ^{+} tale che h_{k}^{[i]} ≠ 0 per i = 1, …, m−1; inoltre non c’è un capo se e solo se esiste k ∈ ℤ^{+} tale che h_{k}^{[m]} ≠ 0. Per esempio per il grafo rappresentato in figura 1 si ha h_{k}^{[5]} = h_{k}^{[1]} = 0, ∀k ∈ ℤ+, mentre h_{k}^{[2]} = h_{k}^{[3]} = h_{k}^{[4]} e h_{k}^{[6]} non sono identicamente nulle. Questo significa che per predire y_{t} serve il passato di u^{[2]}, u^{[3]}, u^{[4]} e di y stesso.

In pratica non si sa a priori quali variabili saranno significativi per predire la uscita y_{t}. I metodi PEM e N4SID non effettuano alcuna selezione di ingressi e stimano un modello “completo” che tiene conto di tutti gli ingressi. Questo è un motivo per cui questi due algoritmi a volte producono stime scarse.

Qui di seguito si presentano degli algoritmi che favoriscono la sparsità e che sono quindi capaci di individuare una struttura del sistema dinamico simile a quella riportata in figura 1. Allo stesso momento riescono a stimare le risposte impulsive h_{k}^{[i]} nulle.

0.3. PREDIZIONE IN RKHS E APPROSSIMAZIONE FINITO-DIMENSIONALE 7

0.3 Predizione in RKHS e approssimazione finito-

dimensionale

Si assume che le risposte impulsive h^{[i]} siano (le versioni campionate) di funzioni in H_{K}, lo spazio di Hilbert riproduttivo associato al nucleo campionato K. Questo corrisponde a dire che le risposte impulsive h^{[i]} sono processi Gaussiani con funzione di covarianza K. Il problema di stimare le risposte impulsive h^{[i]} da dati misurati può essere formulato come il seguente problema di regolarizzazione di tipo Tikhonov:

{ĥ^{[i]}} = arg min_{h^{[i]}∈H_{K}} ∑_{t=t_{0}}^{N}(y_{t} − ŷ_{t|t−1})^{2} + γ_{2}^{2} (∑_{i=1}^{m} ∥h^{[i]}∥_{H_{K}}^{2})  (2)

soggetto a

ŷ_{t|t−1} = ∑_{i=1}^{m−1} [∑_{k=1}^{∞} h_{k}^{[i]}u_{t−k}^{[i]}] + ∑_{k=1}^{∞} h_{k}^{[m]}y_{t−k}.  (3)

Il parametro γ_{2} è il cosiddetto parametro di regolarizzazione che deve regolare il trade-off tra fit e regolarità di h^{[i]}.

Per portare questo problema in dimensione finita si modifica il nucleo K in modo che la realizzazione del nucleo modificato faccia soddisfare la risposta impulsiva alla condizione h_{k}^{[i]} = 0, ∀k > J. È facile vedere che il nucleo (la covarianza, indicata con K_{J}) che soddisfa la condizione è K_{J}(h, k) = 0, ∀(h, k) : k > J oppure h > J.

Raccogliendo i coefficienti della risposta impulsiva in un vettore colonna si ha

\underline{h}^{[i]} := [h_{1}^{[i]}, h_{2}^{[i]}, …, h_{t_{0}}^{[i]}]^{T} .

Il numero J non ha influenza sul trade-off tra bias e varianza, ma serve solo per ragioni di complessità computazionale. Con un piccolo abuso di notazione si indica con K_{J} ∈ ℝ^{j×j} la matrice simmetrica, definita positiva formato dai nuclei campionati fino al ritardo J. Sotto queste condizione la norma ∥h^{[i]}∥_{H_{K_{J}}}^{2} può essere scritta come

∥h^{[i]}∥_{H_{K_{J}}}^{2} = (\underline{h}^{[i]})^{T} K_{J}^{−1}\underline{h}^{[i]}.

Il predittore a un passo ŷ_{t|t−1} diventa

ŷ_{t|t−1} = ∑_{i=1}^{m−1} [∑_{k=1}^{J} h_{k}^{[i]}u_{t−k}^{[i]}] + ∑_{k=1}^{J} h_{k}^{[m]}y_{t−k}.  (4)

Risulta che sotto la restrizione h^{[i]} ∈ H_{K_{J}} la soluzione del problema 2 è

arg min_{\underline{h}^{[i]}∈ℝ^{t_{0}}} ∑_{t=t_{0}}^{N}(y_{t} − ŷ_{t|t−1})^{2} + γ_{2}^{2} ∑_{i=1}^{m} [(\underline{h}^{[i]})^{T} K_{J}^{−1}\underline{h}^{[i]}]  (5)

soggetto a 4. Si vede quindi che il problema può essere riformulato come quello di stimare un modello ARX con un peso di regolarizzazione incluso nel problema della stima.

0.4 Stabil Spline Group LAR (SSGLAR) algo-

ritm

Viene ora discusso come si può risolvere il problema 5 invocando sparsità nei gruppi h_{t_{0}}^{[i]} usando l’algoritmo GLAR. Per fare questo si assume che i parametri θ siano stati fissati; senza alcun’informazione a priori questo parametro viene posto uguale a zero. Il problema 5 è un problema di regressione lineare con una funzione peso in l_{2} che può essere riscritto nella forma

\overline{y}_{t_{0}}^{+} = ∑_{i=1}^{m} \overline{A}_{Ji}h_{t_{0}}^{[i]} + W  (6)

avendo definito

\overline{y}_{t_{0}}^{+} := [\begin{matrix}y_{t_{0}}^{+} \\ 0_{1×(t_{0}·m)} \end{matrix}]
\overline{A}_{Ji} := [A_{Ji} v_{i} ⊗ Λ] Λ := γ_{2}K_{J}^{−1/2}
v_{i} := [\underbrace{0 … 0}_{i−1} 1 \underbrace{0 … 0}_{m−i}]^{T}

La soluzione degli ingressi come descritto prima può essere eseguito dall’algoritmo GLAR applicato al problema di regressione. Si definisce SS-Glar il seguente algotimo:

0.4.1 Algoritmo SS-Glar

1. stabilire il parametro β del nucleo;

2. stabilire il parametro γ_{2} di 5 e formare \overline{A}_{Ji} in 6;

3. stima di h^{[i]} applicando l’algoritmo GLAR a 6;

0.4.2 Stima dei hyperparametri

Prima di poter eseguire l’algoritmo appena visto bisogna fissare i seguenti parametri:

1. la funzione peso γ_{2} in l_{2} di equazione 6 (regolarità di [h_{1}^{[i]}, h_{2}^{[i]}, …, h_{t_{0}}^{[i]}]^{T} nello spazio H_{K_{t_{0}}});

2. il parametro β del nucleo (velocità di decrescita del nucleo);

3. Il numero di blocchi diversi da zero stimati dall’algoritmo GLAR;

Questi parametri possono essere stimati usando un approccio di validazione come segue: Disponendo delle misure {y_{t}, u_{t}}_{t=1,...,N} si suddivide l’insieme in due parti. Il primo {y_{t}, u_{t}}_{t=1,...,2N/3} viene chiamato l’insieme di identificazione, i dati rimanenti l’insieme di validazione. Si esegue l’algoritmo per hyperparametri fissati usando l’insieme di identificazione e successivamente si valuta il modello ottenuto sull’insieme di validazione. Dall’insieme degli hyperparametri (β ∈ ℝ^{+}, γ ∈ ℝ^{+}) si sceglie un insieme con un numero finito (possibilmente piccolo) di alternative, su quale poi deve essere testato l’algoritmo. Inoltre vengono anche testati diversi livelli di sparsità. Alla fine viene scelto il modello identificato che ha dato le miglior prestazioni sul set di validazione in termini di RMS_{1} dell’errore di predizione a un passo:

RMS_{k} := \sqrt{\frac{1}{500} ∑_{i=1}^{500}(y_{t} − ŷ_{t|t−k})^{2}}.

Vengono quindi fissati gli hyperparametri e il livello di sparsità con i quali viene rifatto la stima del modello usando l’intero insieme di dati {y_{t}, u_{t}}_{t=1,...,N}.

0.5 Stable Spline with Exponential Hyperpara-

metri (SSEH) Algoritm

In questo sezione il predittore di equazione 3 viene formulato in un contesto Bayesiano assumendo che le risposte impulsive h^{[i]} siano processi Gaussiani con covarianza λ_{i}K(t, s). Serve specificare i hyperparametri β, θ e λ_{i} e la varianza di rumore σ^{2} che descrivono interamente le densità di probabilità delle misure di y e le risposte impulsive dei predittori h^{[i]}.

0.5.1 Hyperprobabilità a priori per i hyperparametri

La varianza di rumore σ^{2} viene stimata in un passo preliminare usando un modello ARX. Quindi questo parametro viene supposto noto nella stima Bayesiana. I hyperparametri β, θ e λ_{i} vengono invece modellati come variabili aleatorie indipendenti. β ha una densità di probabilità poco informativa su ℝ^{+}, mentre θ è distribuita uniformemente. Ogni λ_{i} è una variabile aleatoria esponenziale con l’inversa della media ξ ∈ ℝ^{+}, cioè

p(λ_{i}) = ξexp(−ξλ_{i})χ(λ_{i} ≥ 0), i = 1, …, m

con χ la funzione indicatrice. ξ può essere interpretato come una variabile aleatoria poco informativa su ℝ^{+}. Tutti i hyperparametri possono essere raccolti in vettore aleatorio indicato con ζ := [λ_{1}, …, λ_{m}, θ_{1}, θ_{2}, β, ξ]. Si noti che i parametri λ_{i}, che hanno lo stesso ruolo di σ^{2}/γ^{2} dove γ^{2} è definito in 2, possono ora essere tutti differenti, il che aumenta la flessibilità del modello.

0.5.2 Il modello Bayesiano intero

Sia A_{ti} ∈ ℝ^{N×∞} in modo che per le uscite future y^{+} := y_{1}^{+} risulta

y^{+} = [∑_{i=1}^{m} A_{1i}h^{[i]}] + e^{+},

0.5. STABLE SPLINE WITH EXPONENTIAL HYPERPARAMETRI (SSEH) ALGORITM11

dove e^{+} := e_{1}^{+}. In pratica y^{−} := [y_{0}, y_{−1}, y − 2, …] non è mai interamente noto e una soluzione consiste nel porre a zero i componenti non noti. Viene quindi imposto la seguente approssimazione:

p(y^{+}, {h^{[i]}}, y^{−}, u|ζ) ∝ [∏_{t=1}^{N} p(y_{t}|{h^{[i]}}, y_{t}^{−}, ζ, u_{t}^{−})] p(y^{−}, {h^{[i]}}, u^{−}|ζ)

≈ [∏_{t=1}^{N} p(y_{t}|{h^{[i]}}, y_{t}^{−}, ζ, u_{t}^{−})] p({h^{[i]}}|ζ)p(y^{−}, u^{−})

Il ∝ deriva dal assunzione che il predittore di u_{t} dato il passato di u_{t}^{−} e y_{t+1}^{−} non dipende da ζ, il ≈ invece dall’assunzione che il passato u_{t}^{−} e y_{t}^{−} porta informazione sulla risposta impulsiva del predittore e sui hyperparametri. Questo modello statistico è descritto dalla rete Bayesiana riportata in figura 2 (sinistra). Come visto sopra, per questioni di complessità computazionale bisogna inoltre troncare la risposta impulsiva h^{[i]} a una certa lunghezza t_{0}.

0.5.3 Stima degli hyperparametri

Il vettore degli hyperparametri ζ viene stimato minimizzando la probabilità marginale, cioè la densità di y^{+}, ζ e {h^{[i]}}. Può essere dimostrato (si guardi [2]) che la stima si ottiene da

ζ̂ = arg min_{ζ} J(y^{+}; ζ)  (7)

con ξ, β > 0, λ_{i} > 0, i = 1, …, m e J è dato da

J^{(}y^{+}; ζ) = \frac{1}{2}log(det[2πV [y^{+}]]) + \frac{1}{2}(y^{+})^{T} (V [y^{+}])y^{+} + ξ ∑_{i=1}^{m} λ_{i} − log(ξ) (8)
V [y^{+}] = σ^{2}I_{N} + ∑_{k=1}^{m} λ_{i}A_{1i}KA_{1i}^{T}.  (9)

Una questione molto importante per l’uso pratica dello schema matematico è l’utilizzo di un buon valore iniziale. Sotto viene descritto un schema con il quale si ottiene una soluzione sub ottimale risolvendo il problema di ottimizzazione in ℝ^{4} in relazione al modello Bayesiano in figura 2 (destra).

Figura 2: Rete Bayesiana che descrive il nuovo modello non parametrico per identificazione di sistemi dinamici sparsi, dove y^{−} := [y_{l−1}, y_{l−2}, …] e nel modello ridotto λ := λ_{1}, …, λ_{m}.

1. Si ottengono {λ̂_{i}}, θ̂ e β̂ risolvendo il problema 7 nella versione modificato

seguente

arg min_{ζ}[J(y^{+}; ζ) − ξ ∑_{i=1}^{m} λ_{i} + log(ξ)]

con β > 0, λ_{1} = … = λ_{m} ≥ 0

2. Si fissano ξ̂ = 1/λ̂_{1} e ζ̂ := [λ̂_{1}, …, λ̂_{m}, θ̂_{1}, θ̂_{2}, β̂, ξ̂]. Poi, per i = 1, …, m si fissa ζ̄ = ζ̂ tranne che la i-esima componente di ζ̄ che viene messa a zero e se J(y^{+}; ζ̄) ≤ J(y^{+}; ζ̂) si fissa ζ̂ = ζ̄.

La procedura descritta al passo 2 potrebbe non essere adatta nel caso in cui gli ingressi sono altamente correlati, perché potrebbe risultare sensitiva dal ordine con cui le componenti di ζ̂ convergono a zero. Per aggirare questo problema si propone di usare l’algoritmo SS-GLAR per effettuare questa selezione. In particolare si può sostituire il passo 2 con il passo seguente:

1. Si fissa γ_{2}^{2} = σ^{2}/λ̂_{1} e si aggiusta il nucleo visto in sezione 0.4 con K(:,:, θ̂, β̂). Si esegue poi l’algoritmo SS-GLAR selezionando il numero di componenti non nulli usando l’approccio di validazione. A questo punto i componenti {j_{1}, …, j_{k}} ⊆ {1, …, m} di ζ̂ vengono annullati, il che equivale a annullare ĥ^{[j]} = 0, j ∈ {j_{1}, …, j_{k}}.

Il valore ζ̂ finale può essere usato come punto iniziale in problema 9.

0.5.4 Stima delle risposte impulsive del predittore per ζ noto

Sotto certi condizioni (vedi [2]) si può dimostrare:

\{ĥ^{[i]}\}_{i=1}^{m} = arg min_{\{f^{i}∈H\}_{i=1}^{m}} ∥y^{+} − ∑_{i=1}^{m} A_{1if^{i}}∥^{2} + σ^{2} ∑_{i=1}^{m} \frac{∥f^{i}∥_{H_{K}}^{2}}{λ_{i}^{2}}

Inoltre per k = 1, …, m + 1 si ha

ĥ^{[i]} = λ_{i}^{2}KA_{1i}^{T}c, c = (σ^{2}I_{n} + ∑_{k=1}^{m+1} λ_{i}A_{1i}KA_{1i}^{T})^{−1} y^{+}

Dopo aver ottenute le stime di {h^{[i]}} semplice formule possono essere usate per ottenere le risposte impulsive f e g in 1 e quindi anche il predittore a k passi (per dettagli vedere [?]).

0.6 Validazione

Una parte molto importante dell’identificazione è la fase di validazione del modello ottenuto. Si cerca di capire la bontà del modello identificato nel descrivere il processo. L’analisi di prestazione deve essere fatta su dati diversi da quelli usati per l’identificazione. I dati disponibili vengono quindi suddivisi in due insiemi:

• dati per l’identificazione,

• dati per la validazione.

La validazione consiste nel confrontare l’uscita reale con l’uscita del modello identificato avendo come ingressi i dati reali e di calcolare successivamente diversi indici che danno informazione sulla bontà del modello. Di seguito vengono spiegati i diversi indici significativi alla fine di fare validazione.

Covarianza e correlazione: La covarianza (campionaria) e la correlazione (campionaria) si possono facilmente calcolare in Matlab cov(x_{1}, x_{2}) e corr(x_{1}, x_{2}) rispettivamente. La covarianza indica se esiste una dipendenza tra le due variabili. Valori elevati indicano una forte legame tra le due variabili, mentre per un legame debole la covarianza tende a zero. La covarianza dipende però della base scelta per rappresentare il modello (per esempio dipende dalla unità di misura utilizzata) e non è dunque un indice oggettivo per confrontare la bontà di diversi modelli Per migliori chiarimenti vedere [1]. La correlazione è invece un parametro adimensionale compreso tra -1 e 1 che indica la relazione lineare tra due variabili. Valori negativi della correlazione indicano che le due variabili sono inversamente correlati, cioè se una variabile cresce l’altra tende a decrescere. Due variabili indipendenti hanno correlazione nulla mentre non è vero il viceversa perché la correlazione può essere nulla anche se esiste una qualsiasi relazioni non lineari tra le due variabili.

Test di bianchezza: I due indici appena visti spiegano l’importanza del test di bianchezza. Serve per verificare se l’errore di predizione (residuo) tra i dati reali e l’uscita dal modello identificato è effettivamente un processo bianco. In base alla bianchezza di una realizzazione del processo vero si decide se il modello identificato è una buona descrizione del sistema dinamico reale. Per effettuare il test si bianchezza è utile il commando resid di Matlab che plotta un grafico per la auto-correlazione del errore di predizione e un grafico con la cross-correlazione tra l’errore di predizione e ogni ingresso reale che poi serve per il test di indipendenza. Inoltre nei grafici viene visualizzato l’intervallo di confidenza al 99%, che è un intervallo di valori plausibili ricavato dalla specifica distribuzione di probabilitá dei residui. Per avere bianchezza tutte le correlazioni (tranne la prima che in realtá è la varianza del residuo) dovrebbero trovarsi all’interno dell’intervallo di confidenza, o almeno non eccessivamente fuori. Nel caso contrario si ha che i campioni dell’errore di predizione sono correlati, il che si può avere per esempio se l’ordine del modello usato per l’identificazione non è adeguato e quindi il modello non è adeguato per descrivere il sistema dinamico vero. Assumendo che l’errore di predizione ε(t) sia bianco si ha che la covarianza r_{ε}(τ) = 0 per ogni τ ≠ 0 e quindi per la covarianza campionaria

r̂_{ε}(τ) = \frac{1}{N} ∑_{t=1}^{N−τ} ε(t + τ)ε(t)

si ha per N → ∞

r̂_{ε}(τ) → 0 ∀ τ ≠ 0,
r̂_{ε}(0) → Eε^{2}(t).

Per il test di bianchezza infine viene considerato la quantità normalizzata

x_{τ} = \frac{r̂_{ε}(τ)}{r̂_{ε}(0)}.

Test di indipendenza: Il test di indipendenza è relativo alla cross-correlazione tra l’errore di predizione e gli ingressi reali. Si verifica se i residui sono effettivamente indipendenti dalle osservazioni passate e quindi se è stato estratto tutta la informazione possibile dai dati. Come accennato prima per concludere che il modello identificato è adeguato si controlla nei grafici della cross-correlazione tra il residuo e gli ingressi reali ottenuto mediante il commando resid che non ci sono picchi significativamente fuori della regione di confidenza. Se per esempio si ha un picco al ritardo k-esimo che sta fuori della regione di confidenza questo significa che l’uscita y(t) non ha ancora estratto tutta l’informazione possibile dall’ingresso u(t − k) e un altro modello potrebbe spiegare meglio l’uscita all’istante t. Il commando resid assume che il residuo ε(t) sia indipendente dagli ingressi passati e sia a media nulla e che si ha quindi r_{εu}(τ) = Eε(t + τ)u(t) = 0. Per τ si considerano sia valori positivi che negativi. Mentre se il modello non è adeguato per descrivere il sistema per valori di τ positivi si può aspettare che ci sono picche di cross-correlazione che sono fuori dell’intervallo di confidenza questo non è detto per valori di τ negativi. Per esempio nel caso che u(t) sia un processo bianco questi picchi sono comunque nulle. Calcolando la correlazione campionaria come

r̂_{εu}(τ) = \frac{1}{N} ∑_{t=1−min(0,τ)}^{N−max(0,τ)} ε(t + τ)u(t)

per il test di indipendenza si considera la quantità normalizzata

x_{τ} = \frac{r̂_{εu}(τ)}{[r̂_{ε}(0)r̂_{ε}(0)]^{1/2}} .

Mappa zeri-poli: Tracciando la mappa zeri-poli di ogni funzione di trasferimento del modello identificato si può dare una prima indicazione sull’esattezza dell’ordine del modello. Se ci sono delle quasi-cancellazioni si può aspettare che è possibile ridurre l’ordine del modello identificato.

fit: Il fit da la percentuale della varianza spiegata ed è calcolato mediante la seguente formula:

fit = 100 ∗ [1 − \frac{norm(y_{h} − y)}{norm(y − mean(y))}]

dove y_{h} sono le uscite dal modello identificato, mentre y è l’uscita reale. Dalla formula si vede che più l’uscita del modello è vicina a quella reale, più vicino è il fit al 100% e piú sarrà quindi adeguato il modello identificato per spiegare il sistema reale.

Errore quadratico medio: Un indice che valuta la prestazione del modello identificato che viene spesso usato è l’MSE (mean square error). Esso è definito da

MSE = \frac{1}{N} ∑_{i=1}^{N}(y_{i} − ŷ_{i})^{2}
= V ar(ŷ) + mean^{2}(ŷ − y)

Questo indice da quindi informazioni sulla bontà dello stimatore in termini di varianza e dispersione della stima. Spesso si calcola anche la radice di questo indice, che viene chiamata con radice dell’errore di predizione (in inglese root mean square error, RMSE) che è quindi proporzionale alla deviazione standard della stima.

Sparsità La sparsità da la percentuale di funzioni di trasferimento nulle nel modello identificato. Avere una funzione di trasferimento nulla significa che la rispettiva misura (ingresso) non è stata utilizzata alla fine di calcolare l’uscita del modello identificato e quindi che il rispettivo sensore può essere eliminato nel calcolo della rispettiva uscita.

18 CONCLUSIONI

Conclusioni

Bibliografia

[1] Giorgio Picci. Metodi statistici per l’identificazione di sistemi lineari. 2011.

[2] Alessandro Chiuso and Gianluigi Pillonetto. Identification of large-scale sparse linear dynamic systems: a regularization based approach. 2010.

[3] Lennart Ljung. System identification - theory for the user. 1999.

[4] Alessandro Chiuso and Gianluigi Pillonetto. Nonparametric sparse estimators for identification of large scale linear systems.

[5] Alessandro Chiuso and Gianluigi Pillonetto. Learning sparse dynamic linear systems using stable spline kernels and exponential hyperpriors.

[6] Federico Bonollo, Daniele Cattai, and Marco Verga. Identificazione termodinamica di un edificio con l’utilizzo di pca e pls. 2007-2008.

[7] Brandly Efron, Trevor Hastie, Iain Johnstone, and Robert Tibshirani. Least angle regression. Journal of the Royal Statistical Society, 2004.

[8] Elena Toffoli, Giancarlo Baldan, and Guido Albertin. Identificazione termodinamica di un edificio. 2006-2007.

[9] Ming Yuan and Yi Lin. Model selection and estimation in regression with grouped variables. Journal of the Royal Statistical Society, 2006.

[10] Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society, 1996.

[11] Sensirion AG. Datasheet SHT1x - Humidity and temperature sensor, 2010.

20 BIBLIOGRAFIA [12] Moteiv Corporation. Tmote sky - humidity, light and temperature sensor

with USB, 2006.