1 Raccolta dati

1.1 Impostazione

Per eseguire le campagne di racca

1.2 I dispositivi utilizzati

Per eseguire le misurazioni delle grandezze di interesse sono stati utilizzati i dispositivi MoteIV Tmote Sky equipaggiati con sensori per il rilevamento di temperatura, umidità e luminosità. Sebbene questi sensori siano predisposti per una comunicazione dei dati mediante rete wireless si è scelto di effettuare l’acquisizione off-line sfruttando la memoria flash per immagazzinare i dati raccolti più delle informazioni aggiuntive che sono stati letti alla fine mediante usb. La scelta di utilizzare questa modalità di lettura è dipesa principalmente dall’esigenza di un attento utilizzo del consumo energetico, infatti sebbene questi dispositivi siano progettati per garantire un basso consumo ed abbiano un processore in grado di funzionare con tensioni che arrivano fino a 1.9 V i componenti necessari all’esperimento (i sensori e la memoria flash) necessitano di una tensione costante di circa 3 V fornita da due batterie di tipo AA alcaline. Se la tensione in ingresso diminuisse eccessivamente potrebbero accadere dei malfunzionamenti al convertitore ADC con una perdita di bit significativi (da 14 bit a 12 bit upper) e alla memoria flash che per una scrittura sicura richiede una almeno 2.6 V. Si fa notare che la scelta di una alimentazione mediante porta usb per questa esperienza non è possibile poichè in questo caso l’alimentazione a 5 V provoca un surriscaldamento del processore posto in prossimità del sensore di temperatura, sfalsando così la misura di temperatura. In questo modo si è riusciti ad avere una alimentazione corretta durante tutta la raccolta dei dati, in Figura ??? è riportato l’andamento della tensione di ogni dispositivo raccolta proprio per monitorare il buon funzionamento di ogni dispositivo. Viene riportata qui di seguito una breve descrizione dei sensori utilizzati dai Tmote e successivamente del software necessario:

1. Umidità e temperatura: Il sensore di temperatura e umidità utilizzato sul mote è il sensore Sensiron SHT11. Esso consta di un sensore di temperatura il cui principio di funzionamento si basa sulla dipendenza dalla temperatura della banda di non conduzione di un semiconduttore ed ha range di funzionamento: -40/+123.8 C. La sua risoluzione è a 14 bit (può essere ridotta a 12 bit per applicazioni che richiedono velocità elevate o un utilizzo con potenza in ingresso molto bassa). La calibrazione è gestita internamente in base a delle tabelle, la formula di conversione dalla parola a 14 bit T_{raw} presente nella memoria flash a gradi centigradi è

T = b_{1} + b_{2}T_{raw} − b_{3}(T_{raw} − b_{4})^{2} [^{∘}C]  (1.1)

con costanti b_{1} = −39.6, b_{2} = 0.01, b_{3} = −2 · 10^{−8} e b_{4} = 7000.

Il sensore per rilevare l’umidità è basato sulla variazione di capacità di un polimero in funzione dell’umidità presente nell’ambiente circostante e ha range di funzionamento 0/100 % RH. La misura viene registrata in una parola a 12 bit H_{raw} che può scendere ad 8 bit se necessario. La formula per la conversione è:

RH_{true} = (T − 25)(c_{1} + c_{2}H_{raw}) + c_{3} + c_{4}H_{raw} + c_{5}H_{raw}^{2}[%]  (1.2)

con constanti c_{1} = 0.01, c_{2} = 8 · 10^{−6}, c_{3} = −4,c_{4} = 0.0405 e c_{5} = −2.8 · 10^{−6}. Per maggiori informazioni si rimanda al datasheet ???.

2. Luminosità I sensori presenti nel mote per misurare la quantità di luce sono due fotodiodi prodotti dalla Hamamatsu. In particolare è presente il modello S1087 sensibile alla sola luce visibile (320 - 730nm) e il modello S1087-01 che rivela anche lo spettro infrarosso (320 - 1100nm). Le misure sono memorizzate in parole da 16 bit unsigned e la conversione in lux é data da:

AR = 625AR_{raw} [lux]  (1.3)

con AR_{raw} il valore in bit fornito dal sensore S1087. Analogamente a partire dal dato in bit T R_{raw} del sensore S1087-01 si implementa

T R = 76.9T R_{raw} [lux] .  (1.4)

Per la fase di identificazione la misura che viene utilizzata è quella fornita dal secondo sensore, ovvero quella comprensiva della radiazione visibile e dell’infrarosso poichè entrambe influiscono sul comportamento dinamico dell’edificio.

3. Memoria flash Il mote è equipaggiato con una unità di memoria flash esterna di tipo seriale di capienza 1024kb. Essendo tuttavia configurata in lettura e scrittura solo nei primi quindici settori lo spazio effettivamente disponibile è di 983040 byte. La struttura dati presente consta di 6 campi da 16 bit ciascuno, il che significa una capacità totale di 81920 campioni. Questo valore è in realtà ridotto a 65535 campioni poichè parte dello spazio è utilizzato per controllare se dei dati non vengono scritti in memoria (di fatto l’ultimo dei sei campi è occupato da un numero che aumenta di un passo ad ogni istante di campionamento). Contando che il tempo di campionamento adatto per le dinamiche termiche di un edificio è dell’ordine dei 5-10 minuti si capisce che la memoria disponibile nel mote è più che sufficiente per lo svolgimento dell’esperienza.

Si fornisce qui di seguito una breve descrizione del codice che viene utilizzato per la gestione dei mote. La sessione di acquisizione viene suddivisa in tre fasi dal firmaware caricato su ogni mote. La prima parte è dedicata al setup che viene effettuato comunicando con il mote tramite usb e pc. In questa fase vengono gestite le operazioni principali di lettura della memoria flash e di cancellazione della stessa. La seconda fase è dedicata alla sincronizzazione e il passaggio avviene allo scadere di un timer opportunamente impostato. In questa fase le comunicazioni usb sono interrotte e gli unici eventi che il mote è può rilevare sono ricezione del pacchetto radio di sincronizzazione o la pressione del tasto user (che causa la partenza di un nuovo timer e l’invio di un nuovo pacchetto di sincronizzazione). Con la terza fase si inizia l’acquisizione vera e propria dei dati. Il modulo radio viene disattivato ed un timer posto uguale al tempo di campionamento comanda la lettura dei dati rilevati dai sensori e la scrittura su flash. Il principio di funzionamento è riportato in Figura 1.1. L’interfaccia con il PC è sviluppata in Java su piattaforma Linux ed utilizza le librerie di comunicazione seriale e le classi per la gestione dei pacchetti fornite con TinyOS 2.x. La lettura del mote fornisce in outuput un file di testo contenente i valori registati dai sensori che vengono convertiti in unità ingegneristiche da uno script in Matlab secondo (1.2)-(1.4). La programmazione dei mote è avvenuta eseguendo lo script EggMote.h in cui si precedentemente impostato il tempo di campionamento.

Figure 1.1: Struttura degli stati e delle funzionalità del mote

2 Metodi di identificazione ai sottospazi

I metodi di identificazione a sottospazi sono particolarmente adatti per l’identificazione di sistemi MISO (o MIMO) di grandi dimensioni per la loro efficienza computazionale. Vengono implementati utilizzando strumenti derivati dall’algebra lineare (per effettuare opportune proiezioni ortogonali e oblique) come la decomposizione QR e SVD di cui esistono implementazioni molto efficienti. Si noti che in questo modo si evita l’utilizzo di minimizzazioni ricorsive e quindi si evitano i problemi di convergenza che esse comportano (e che possono capitare con i metodi PEM). Il problema è dunque, date le osservazioni {y_{1}, …, y_{T} , u_{1}, …, u_{T} }, y_{t} ∈ ℝ^{p}, u_{t} ∈ ℝ^{m} costruire il modello:

\{ \begin{matrix}x(t + 1) & = & Ax(t) + Bu(t) + Ke(t) \\ y(t) & = & Cx(t) + e(t) \end{matrix}  (2.1)

che descriva i dati (con e(t) rumore bianco). A partire dal processo stostico y(t) si scrivono i dati nella forma {y_{t}, …, y_{t+N−1}} = Y_{t} per creare una corrispondenza biunivoca (al crescere del numero N di dati) tra lo spazio astratto di variabili aleatorie span{y(t), …, y(s)} e lo spazio generato dalle righe Y_{t} della matrice di Henkel diagonale a blocchi:

Y_{[t,s]} = [\begin{matrix}y_{t} & y_{t+1} & … & y_{t+N−1} \\ y_{t+1} & y_{t+2} & … & y_{t+N} \\ ⋮ & ⋮ & ... & ⋮ \\ y_{s} & y_{s+1} & … & y_{s+N−1} \end{matrix}]  (2.2)

di dimensione Y_{[t,s]} ∈ ℝ^{p(s−t)xN} e da questa si definiscono le matrici delle uscite passate e future

Y_{[0,2t−1]} = [\begin{matrix}Y_{p} \\ Y_{f} \end{matrix}] = [\begin{matrix}Y_{[0,t−1]} \\ Y_{[t,2t−1]} \end{matrix}]  (2.3)

e si svolgono delle operazioni analoghe con i dati d’ingresso u(t). A questo punto si può costruire una base per lo spazio di stato X_{t} attraverso una proiezione adatta delle sequenze dei dati e successivamente da questa base si risolve nel senso dei minimi quadrati il sistema:

\{ \begin{matrix}X̂_{t+1} & ≃ & AX̂_{t} + BÛ_{t} + KE_{t} \\ Y_{t} & ≃ & CX̂_{t} + E_{t} \end{matrix}  (2.4)

dove Ê_{t} := Y_{t} − ĈX_{t}. Si definisce poi la matrice di osservabilità del sistema come:

Γ := [ C ; CA ; CA^{2} ; ⋮ ; CA^{f}]  (2.5)

e matrici triangolari inferiori a blocchi di Toeplitz con i parametri di Markov del sistema:

H_{d} = [\begin{matrix}0 & 0 & … & 0 \\ CB & 0 & … & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ CA^{f−1}B & CA^{f−2}B & … & 0 \end{matrix}] ,  (2.6)

H_{s} = [\begin{matrix}I & 0 & … & 0 \\ CK & I & … & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ CA^{f−1} & CA^{f−2}K & … & I \end{matrix}] .  (2.7)

Utilizzando ora la relazione

Y_{[p,p+f]} = ΓX_{p} + H_{d}U_{[p,p+f]} + H_{s}E_{[p,p+f]}  (2.8)

si può effettuare la proiezione obliqua ^{1} (mediante decomposizione QR)

ΓX_{p} ≃ Ŷ_{[p,p+f]} := E_{||U_{[p,p+f]}} [Y_{[p,p+f]}|Z_{[0,p)}] .  (2.10)

Figure 2.1: Rappresentazione geometrica della proiezione obliqua effettutata in (2.10)

1 La proiezione obliqua é la scomposizione in due sottospazi tra loro non ortogonali. Si denota come Ê_{||B}[·|C_{t}], ossia la proiezione lungo lo spazio generato delle righe di B sullo spazio generato dalle righe di C_{t}. La formula estesa é:

Ê_{∥B_{t}}[A_{t}|C_{t}] = Σ̂_{ac|b}Σ̂_{cc|b}^{−1}C_{t}  (2.9)

Dove sotto le ipotesi di ergodicitá le matrici di covarianza campionaria (condizionata e non) sono: Σ̂_{ab} := \frac{1}{N} ∑_{i=0}^{N−1} a_{t+i}b_{t+i}^{⊤} e Σ̂_{ab|c} := Σ̂_{ab} − Σ̂_{ac}Σ̂_{cc}^{−1}Σ̂_{cb}

Per stimare Γ ( e se necessario X_{p} ) si utilizza una decomposizione "‘SVD pesata"’

U_{n}S_{n}V_{n}^{T} ≃ W_{r}Ŷ_{[p,p+f]}W_{c}  (2.11)

e si ottiene

Γ̂ := W_{r}^{−1}U_{n}S_{n}^{1/2}  X̂_{p} := Γ̂_{r}^{−L}Ŷ_{[p,p+f]}  (2.12)

3 Stato dell’arte

L’identificazione black-box è ampiamente utilizzata per ricostruire la dinamica di modelli da un insieme finito di dati ingresso-uscita. In questo scenario una situazione particolare che si verifica in svariati ambiti (ad esempio in ingegneria chimica, in econometria, nella visione computazionale, e così via) è quella in cui si richiede l’identificazione di sistemi di grandi dimensioni in cui entrano in gioco un gran numero di variabili. In un contesto simile un punto cardine delle procedure di identificazione dovrebbe essere quello di favorire la sparsità dei risultati, ovvero di estrarre dal gran numero di sottosistemi presenti nel sistema principale solo il sottoinsieme che influenza l’uscita del sistema in maniera più significativa. Il fatto che ogni variabile possa dipendere solo da un sottoinsieme di altre variabili si esprime nel caso Gaussiano statico con delle condizioni di indipendenza condizionata fra sottoinsiemi di variabili. Nel caso dinamico, cioè quando i dati osservati sono traiettorie di processi stocastici (possibilmente stazionari) le condizioni di indipendenza condizionale si esprimono nel fatto che la predizione dell’andamento di una variabile (che verrà chiamata uscita) richieda solo il passato di qualche altra variabile (che verranno chiamate ingressi) più eventualmente il proprio. Si può rappresentare ciò con un grafo in cui i nodi sono le variabili, gli archi diretti sono le funzioni di trasferimento non nulle e i cappi (archi che si chiudono sullo stesso nodo da cui partono) codificano la dipendenza della variabile dal proprio passato, come si può vedere in Figura ??. In generale sia il sistema dinamico sia la struttura che rappresenta il legame fra le variabili dovranno essere estratti dai dati. In questo contesto di identificazione gli usuali metodi di minimizzazione dell’errore di predizione (PEM) possono diventare computazionalmente problematici poichè la teoria sottostante ai metodi di stima dell’ordine del modello (ad esempio AICC e BIC) inizia a vacillare. Questo problema insieme al problema legato alla presenza di non-convessità è stato risolto con l’introduzione degli algoritmi ai sottospazi che invece di fornire la stima utilizzando metodi iterativi di ottimizzazione non lineare (come in PEM), si basano su degli strumenti efficienti e robusti derivati dall’algebra lineare come la decomposizione QR e SVD. Tuttavia anche per questi metodi la presenza di un gran numero di parametri da stimare rispetto al numero di dati disponibili resta un problema. In questa situazione dunque è possibile che sia gli algoritmi parametrici tradizionali che quelli ai sottospazi possano fornire risultati non soddisfacenti. Oggigiorno sono stati introdotti dei nuovi approcci parametrici (come LAR e LASSO) che inducono sparsità nel risultato sfruttando la norma l_{1}. Un siffatto termine peso porta con sè la cosiddetta proprietà di bi-separazione, cioè la preferenza di soluzioni con molti termini nulli a spesa di pochi termini ma di modulo maggiore. Sono inoltre stati ampliati con delle versioni "a gruppi" (GLars e GLasso) in cui interi gruppi di termini sono posti a zeri ottenendo così una sparsità in termini di funzioni di trasferimento piuttosto che di singoli termini. Tuttavia la maggior parte del lavoro in questo ambito è stato svolto in uno scenario statico a discapito del caso dinamico. Gli algoritmi che verranno testati maggiormente in questo lavoro adottano un punto di vista Bayesiano e prendono come punto di partenza un nuovo paradigma sviluppato in [??] che si basa sulla stima non parametrica delle risposte impulsive. In questo modo invece di dover postulare una struttura finito dimensionale per le funzioni di trasferimento del sistema (con ordine quindi da stimare) si cerca la risposta impulsiva all’interno di uno spazio infinito dimensionale. Il mal condizionamento intrinseco a un problema così formulato viene risolto utilizzando dei metodi di regolarizzazione Bayesiana. Più nel dettaglio si può dire che le risposte impulsive sono modellate come dei processi Gaussiani la cui autocovarianza viene chiamata stable spline kernel ed include il vincolo di BIBO stabilità del risultato.

4 Preliminari: nuclei per l’identificazione

4.1 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 4.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}} in funzione degli altri componenti. A questo scopo si definisce y_{t} := z_{t}^{[i]} (la componente i-esima di z_{t}) come l’uscita e tutte le altre componenti u_{t} := z_{t}^{[−i]} ∈ ℝ^{m−1} come gli 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} in funzione delle rimanenti. In questa sezione però per questioni di chiarezza si sceglie un valore specifico 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}, y_{t−2}, …], U^{t} = [u_{t−1}, u_{t−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. Possono essere quindi descritti dalla seguente relazione:

y_{t} = ∑_{k=1}^{∞} f_{k}u_{t−k} + ∑_{k=0}^{∞} g_{k}e_{t−k}  (4.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. Esso è equivalente all’errore di predizione di un passo

e_{t} := y_{t} − ŷ_{t|t−1} = y_{t} − E[y_{t}|Y^{t}, U^{t}]  (4.2)
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}.  (4.3)

Figure 4.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 necessarie 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 il grafo rappresentato in Figura 4.1 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 significative per predire l’uscita y_{t}. I metodi PEM e N4SID non effettuano alcuna selezione di ingressi e stimano un modello “completo” che tiene conto di tutte le misure. 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 4.1 ed allo stesso momento riescono a stimare le risposte impulsive h_{k}^{[i]} nulle.

4.2 La regolarizzazione

Un approccio ampiamente utilizzato per ricostruire una funzione da un insieme di misure {y_{t}} consiste nel minimizzare un funzionale di regolarizzazione appartenente ad un RKHS (Reproducing Kernel Hilbert Space) H associato ad un nucleo K simmetrico e definito positivo. Date N misure, la regolarizzazione ai minimi quadrati in H stima la funzione sconosciuta come

ĥ = arg min_{h} ∑_{t=1}^{N}(y_{t} − Γ_{t}[h])^{2} + η∥h∥_{H}^{2}  (4.4)

dove {Γ_{t}} sono dei funzionali lineari e limitati in H legati al modello di misura mentre il parametro η è uno scalare positivo che permette di settare il trade-off fra l’errore di predizione e la regolarità della soluzione. La funzione che minimizza (4.4) può essere interpretata da un punto di vista Bayesiano come la stima a minima varianza di h quando h è visto come un processo Gaussiano a media nulla con autocovarianza K e {y_{t} − Γ_{t}[h]} un rumore bianco indipendente da h. Spesso la conoscenza a priori è limitata al fatto che i segnali in gioco siano continui e a energia limitata. In questo caso la h viene modellata con un integrale pdimensionale di un rumore bianco di intensità unitaria. L’autocorrelazione W_{p} di h si può così scrivere come:

W_{p}(s, t) = ∫_{0}^{1} G_{p}(s, u)G_{p}(t, u)du  (4.5)
G_{p}(r, u) = \frac{(r − u)_{+}^{p−1}}{(p − 1)!} , (u)_{+} = \{ \begin{matrix}u if u ≥ 0 \\ 0 if u < 0 \end{matrix}  (4.6)

dove G_{p} è l’autocovarianza associata all’interpretazione Bayesiana di una spline di regolarizzazione di ordine p. Se p = 2 si ottiene un nucleo a spline cubica. Il problema principale nell’utilizzo di un nucleo come (4.6) è che non tiene conto della stabilità della risposta impulsiva. Infatti se l’autocovarianza di h è proporzionale a W_{p} (come nel caso cubico) la varianza di h(t) è nulla nell’istante iniziale e poi diverge all’infinito. Per questo motivo è stato definito un nucleo apposta per l’identificazione di sistemi lineari in cui h rappresenti la risposta impulsiva di un sistema stabile (con covarianza non nulla nell’instante iniziale che converge esponenzialmente a zero). Questo viene garantito prendendo la covarianza proporzionale alla classe di nuclei

K_{p}(s, t) = W_{p}(e^{−βs}, e^{−βt}) s, t ∈ ℝ^{+}  (4.7)

dove β è uno scalare positivo che regola il decadimento della varianza; essendo in pratica un parametro sconosciuto sarà trattato come un iperparametro da stimare a partire dai dati. Tenendo presente le formule precedenti si vede che con p = 2 l’autocovarianza diventa ciò che in letteratura viene chiamato Stable Spline kernel:

K_{2}(t, τ) = \frac{e^{−β(t+τ)}e^{−β max(t,τ)}}{2} − \frac{e^{−3β max(t,τ)}}{6}  (4.8)

E’ stato dimostrato (vedi[??]) che se h è un processo Gaussiano a media nulla con autocovarianza K_{2} allora con probabilità uno le realizzazioni di h sono risposte impulsive di sistemi BIBO stabili. E’ chiaro che lavorando in pratica con segnali a tempo discreto si parlerà di versioni campionate h_{k} di h(t) così che {h_{k}}_{k∈ℤ^{+}} verranno viste come realizzazioni campionate del nucleo K(i, j).

Si considera ora il problema di stimare il parametro θ ∈ ℝ^{m} dal modello lineare

Y = Xθ + W  (4.9)

dove si ha Y ∈ ℝ^{N} è il vettore dei dati di uscita, X ∈ ℝ^{N×m} è il vettore di regressione e W ∈ ℝ^{N} è un rumore bianco vettoriale a media nulla con 𝔼[W W^{T}] = σ^{2}I. Ciò che verrà fatto è considerare ogni componente di θ appartenente ad uno spazio infinito dimensionale RKHS con nucleo K(s, t) come in (4.7). In altre parole θ_{i} viene visto come un processo gaussiano con covarianza K(s, t). La matrice di regressione X sarà un operatore lineare contente la storia passata di u e y.

Il problema di selezionare le variabili più significative può essere affrontato come un problema di ricerca sparsità nella soluzione con alcune importanti differenze. In particolare si tratta ora di un problema a gruppi, poichè per escludere dalla procedura di identificazione delle variabili è necessario che vengano poste a zero delle funzioni di trasferimento. Ogni gruppo di parametri descrive dunque una risposta impulsiva h^{[i]}. Per evitare il problema relativo alla selezione del modello e della parametrizzazione non lineare si preferisce non utilizzare algoritmi a gruppi parametrici scegliendo di lavorare in un contesto non parametrico. I parametri vivono dunque in un RKHS che come già detto possono essere trattati mediante probabilità a priori e regolarizzazione. I due approcci che verranno utilizzati sono:

1. SS-GLAR: Versione a gruppi del problema posto dal Lasso ^{2} portato in un contesto nonparametrico dove i gruppi h^{[i]} appartengono ad uno spazio infinito dimensionale. Per includere il vincolo del problema è necessario risolvere un problema di regolarizzazione misto l_{1} − l_{2}.

2. SSEH: un modello gerarchico dove h^{[i]} è un processo Gaussiano con covarianza λ_{i}K(s, t) e il fattore di scala λ_{i} ha una distribuzione esponenziale per favorire la sparsità sullo spazio dei fattori di scala.

A questo punto si rende necessario introdurre della notazione aggiuntiva per poter spiegare gli algoritmi che verranno utilizzati. Si definisce dunque

y_{t}^{+} := [ y_{t} ; ⋮ ; y_{t+N−1}] e_{t}^{+} := [ e_{t} ; ⋮ ; e_{t+N−1}] h := [ h^{[1]} ; ⋮ ; h^{[m]}] ,  (4.11)

dove h^{[i]} ∈ l_{1}(ℤ^{+}), i = 1, …, m. Il predittore (4.3) può essere così riscritto ^{3} come:

y_{t}^{+} = \underbrace{[A_{t1} … A_{tm}]}_{:=A_{t}} h + e_{t}^{+}  (4.12)

con A_{ti} ∈ ℝ^{N×∞}, i = 1, …, m definiti da:

\begin{aligned}A_{ti}^{[jk]} & := & u_{t−j−k}^{[i]} , i = 1, …, m − 1 , j, k ∈ ℤ^{+} \\ A_{tm}^{[jk]} & := & y_{t−j−k}, j , k ∈ ℤ^{+} \end{aligned}  (4.13)
θ̂:= arg min_{θ} ∥Y − Xθ∥_{2}^{2} + γ_{1}∥θ∥_{1}  (4.10)

2Nel Lasso la scelta del regressore avviene risolvendo il problema nella forma:

con γ_{1} parametro di regolarizzazione compreso fra 0 e 1.

3Il prodotto tra matrici semi-infinite va inteso come il limite della sequenza finita, che sotto l’ipotesi h^{[i]} ∈ l_{1}ℤ^{+}, è ben definito.

e il problema di identificazione può essere dunque visto come la stima di h in (4.11),(4.12) soggetto al vincolo che h^{[i]} ∈ l_{1}(ℤ^{+}), i = 1, …, m. Come già detto lo stimatore dovrà automaticamente selezionare tra {u^{[1]}, …, u^{[m]}, y} le variabili necessarie alla predizione di y e quali no e questo farà in modo che alcune risposte impulsive ĥ^{[i]} siano poste esattamente a zero.

5 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 al fatto che le risposte impulsive h^{[i]} sono processi Gaussiani con funzione di covarianza K. Il problema di stimare le risposte impulsive h^{[i]} dai dati misurati può essere formulato come il seguente problema di tipo Tikhonov di regolarizzazione:

{ĥ^{[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})  (5.1)

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}.  (5.2)

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 soddisfi 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} ,

dove t_{0} = J. 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}.  (5.3)

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

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.4)

soggetto a 5.3. 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.

6 Stabil Spline Group LAR (SSGLAR) algoritm

Viene ora discusso come si può risolvere il problema 5.4 invocando sparsità nei gruppi h_{t_{0}}^{[i]} usando l’algoritmo GLAR. Per fare questo si assume che i parametri θ siano stati fissati; senza alcuna informazione a priori questo parametro viene posto uguale a zero. Il problema 5.4 è 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.1)

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 selezione degli ingressi come descritto prima può essere eseguito dall’algoritmo GLAR applicato al problema di regressione. Nella prossima sezione si definisce nel dettaglio questo algoritmo che viene chiamato SS-GLAR.

6.1 Algoritmo SS-GLAR

I passi che compie l’algoritmo SS-GLAR sono:

1. stabilire il parametro β del nucleo;

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

3. stima di h^{[i]} applicando l’algoritmo GLAR al problema 6.1;

6.2 Stima degli iperparametri

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

1. la funzione peso γ_{2} in l_{2} nell’equazione 6.1 (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 il seguente approccio di validazione: 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, il secondo con i dati rimanenti viene chiamato l’insieme di validazione. Si esegue l’algoritmo per iperparametri fissati usando l’insieme di identificazione e successivamente si valuta il modello ottenuto sull’insieme di validazione. Dall’insieme degli iperparametri (β ∈ ℝ^{+}, γ ∈ ℝ^{+}) si sceglie un insieme con un numero finito (possibilmente piccolo) di alternative, sul quale poi si deve testare l’algoritmo. Vengono inoltre 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} (Root Mean Square) dell’errore di predizione a un passo, dove il generico RMS_{k} (k numero di passi) è definito da:

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

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

6.3 Stable Spline with Exponential Hyperpprior (SSEH) Algoritm

In questa sezione il predittore di equazione 5.2 viene formulato in un contesto Bayesiano assumendo che le risposte impulsive h^{[i]} siano processi Gaussiani con covarianza λ_{i}K(t, s). Serve specificare gli iperparametri β, θ 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]}.

6.3.1 Iperprobabilità a priori per gli iperparametri

La varianza di rumore σ^{2} viene stimata in un passo preliminare usando un modello ARX per poi supporla nota nella stima Bayesiana. Gli iperparametri β, θ e λ_{i} vengono invece modellati come variabili aleatorie indipendenti. β ha una densità di probabilità non informativa (cioè inizialmente non si hanno informazioni a priori sul parametro) su ℝ^{+}, mentre θ è distribuita uniformemente. Ogni λ_{i} è una variabile aleatoria esponenziale con media 1/ξ ∈ ℝ^{+}, quindi vale

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

con χ la funzione indicatrice. ξ può essere interpretata come una variabile aleatoria non informativa su ℝ^{+}. Tutti gli iperparametri possono essere raccolti in un vettore aleatorio indicato con ζ := [λ_{1}, …, λ_{m}, θ_{1}, θ_{2}, β, ξ]. Si noti che i parametri λ_{i}, che hanno lo stesso ruolo di σ^{2}/γ^{2} con γ^{2} definito in 5.1, possono ora essere tutti differenti aumentando così la flessibilità del modello.

6.3.2 Il modello Bayesiano intero

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

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

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 le componenti non note. Viene quindi imposta 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^{−})

La proporzionalità deriva dall’assunzione che il predittore di u_{t} dato il passato u_{t}^{−} e y_{t+1}^{−} non dipenda da ζ mentre l’approssimazione deriva dall’assunzione che il passato u_{t}^{−} e y_{t}^{−} porti informazione sulla risposta impulsiva del predittore e sugli iperparametri. Questo modello statistico è descritto dalla rete Bayesiana riportata in Figura 6.1 (sinistra). Come visto sopra, per questioni di complessità computazionale bisogna troncare la risposta impulsiva h^{[i]} a una certa lunghezza t_{0}.

6.3.3 Stima degli iperparametri

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

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

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

J(y^{+}; ζ) = \frac{1}{2}log(det[2πV [y^{+}]]) + \frac{1}{2}(y^{+})^{T} (V [y^{+}])^{−1}y^{+} + ξ ∑_{i=1}^{m} λ_{i} − log(ξ)  (6.3)

V [y^{+}] = σ^{2}I_{N} + ∑_{k=1}^{m} λ_{i}A_{1i}KA_{1i}^{T} .  (6.4)

Una questione molto importante per l’utilizzo pratico di questo procedimento matematico è l’utilizzo di un buon valore iniziale. Nel seguito viene descritto un schema con il quale è possibile ottenere una soluzione sub-ottimale risolvendo il problema di ottimizzazione in ℝ^{4} in relazione al modello Bayesiano di Figura 6.1 (destra).

(i) Si ottengono {λ̂_{i}}, θ̂ e β̂ risolvendo il problema 6.2 modificato nella seguente versione:

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

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

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

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

La procedura descritta al passo (ii) potrebbe non essere adatta nel caso in cui gli ingressi sono altamente correlati, perché potrebbe risultare sensibile all’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 (ii) con il passo seguente:

(ii’) Si fissa γ_{2}^{2} = σ^{2}/λ̂_{1} e si aggiusta il nucleo visto in sezione 6 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 6.4.

6.3.4 Stima delle risposte impulsive del predittore per ζ noto

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

\{ĥ^{[i]}\}_{i=1}^{m} = arg min_{\{f^{i}∈H\}_{i=1}^{m}} ∥y^{+} − ∑_{i=1}^{m} A_{1i}f^{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} + ∑_{i=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 4.1 e quindi anche il predittore a k passi (per dettagli vedere [?]).