
\documentclass[journal]{IEEEtran}
% correct bad hyphenation here
%\hyphenation{op-tical net-works semi-conduc-tor}
\hyphenation{op-tical net-works semi-conduc-tor IEEEtran}
\usepackage[latin1]{inputenc}
\usepackage{graphicx}  
\usepackage{subfigure} 
\usepackage{url}       
\usepackage{array}
\usepackage[colorlinks=false]{hyperref}
\usepackage{amsfonts}
\usepackage{amssymb}

\begin{document}
% paper title
\title{Identificazione termodinamica di un edificio con l'utilizzo di PCA e PLS}
%
\author{Bonollo~Federico,~Cattai~Daniele,~Verga~Marco}
%\author{Bonollo~Federico}

% The paper headers
\markboth{Progettazione di sistemi di controllo a.a.~2007-2008}{Shell \MakeLowercase{\textit{et al.}}: Bare Demo of IEEEtran.cls for Journals}
% The only time the second header will appear is for the odd numbered pages
% after the title page when using the twoside option.
% 
% make the title area
\maketitle


\begin{abstract}
Una delle applicazioni più' comuni per reti di sensori wireless e' quella del monitoraggio di quantità di interesse quali temperature, livelli di sostanze biochimiche, etc. 
Spesso questi sensori sono disposti in modo casuale e in soprannumero in termini di qualità di stima del campo di interesse. Nel caso di un edificio e al  monitoraggio della temperatura interna, in linea di principio, sarebbe possibile ottenere delle stime di temperatura con una accuratezza elevata tramite l'utilizzo di un modello dinamico di stato in cui solamente poche misure sono necessarie per ricostruirlo. In questa trattazione si studiano due particolari tecniche: PCR e PLS, rispettivamente \emph{Principal Component Regression} e \emph{Partial Least Squares}. Queste due metodologie  vengono inserite in un contesto di identificazione a sottospazi che permette di ricavare un modello dell'edificio in questione. Alla loro applicazione su dati rilevati sperimentalmente, vengono affiancate delle simulazioni effettuate utilizzando il software \emph{DesignBuilder} con lo scopo di fornire diversi set di esperimenti su qui creare e testare il modello, oltre che avere un qualitativo confronto con le tecniche sopracitate
\end{abstract}
\begin{keywords}
identificazione a sottospazi, PCA, PCR, PLS, DesignBuilder.
\end{keywords}
% Note that keywords are not normally used for peerreview papers.

% For peerreview papers, inserts a page break and creates the second title.
% Will be ignored for other modes.
\IEEEpeerreviewmaketitle



\section{Introduzione}

\PARstart{N}{egli} ultimi anni l'interesse verso sistemi abitativi e industriali ecocompatibili \'e andato via via crescendo in relazione alla stessa crescita della sensibilit\'a ambientale. Il motore principale di questo crescente interesse \'e il massiccio utilizzo di risorse energetiche non rinnovabili e del conseguente impatto ambientale che a livello globale si \'e creato. \'E nata dunque l'esigenza di coordinare le ricerche in più ambiti per ottenere strategie di utilizzo energetico valide e nel rispetto della salvaguardia dell'ambiente. Purtroppo a livello domestico si ha una scarsa efficienza nell'utilizzo delle risorse energetiche in particolare, il riscaldamento e il  condizionamento, che nelle abitazioni risentono spesso dell'efficienza strutturale dell'edificio. \\
In tale contesto un'efficace identificazione del modello termodi\-namico si rivela utile per un mirato controllo dei parametri fisici. Ciò permette inoltre di  testare e classificare il rendimento termodinamico dell'edificio.\\ Questa trattazione si focalizza nel determinare attraverso  misure reali e simulative, un modello  per l'evoluzione della temperatura nello spazio e nel tempo dell'edificio. Conoscendo tale modello si vuole verificare e fornire un'idea qualitativa sul posizionamento ottimo dei sensori per un futuro controllo o classificazione dell'edificio, applicando i risultati ottenuti anche ad altri contesti data la grande versatilità delle reti di sensori wireless. Per la tematica della classificazione,  a livello europeo esiste gi\'a una normativa in materia la 2002/91/EC per l'introduzione della classe energetica anche al settore edilizio. A tale scopo gi\'a diverse aziende e imprese di costruzioni si sono rivolte alla ricerca di metodologie e software per la determinazione di tale classe, e una di queste viene descritta e sfruttata in questo progetto, ossia $DesignBuilder$. Il progetto cercherà inoltre di migliorare e comparare i risultati ottenuti nel lavoro \cite{baldan}.  In tale lavoro la selezione di ingressi/uscite significativi, \'e stata effettuata utilizzando algoritmi di tipo Greedy Search, e Genetic Search, senza focalizzarsi sulla fase identificativa. Con l'uso delle tecniche di selezione $PCR$ e $PLS$, \'e possibile constatare un effettivo miglioramento nella bont\'a e consistenza del model\-lo identificato, permettendo ampie conclusioni. Per quanto riguarda l'identificazione vera e propria sono state utilizzate tecniche a sottospazi del tipo: $N4SID$, $PBSID$ e una versione ottimizzata e gentilmente fornita dal prof. A.Chiuso, denominata $PBSIDopt$, mentre attraverso l'algoritmo PLS è stato identificato un modello MIMO-ARX. La teoria a sottospazi \'e stata ampiamente sviluppata in letteratura e vi sono algoritmi, tra cui quelli sopracitati, che forniscono risultati robusti. Si \'e studiata l'identificazione a sottospazi in maniera approfondita quanto bastava per poter intervenire sul problema, descrivendone alcuni dettagli e basi teoriche.\\ L'ottimalità della soluzione, ovvero il miglior sistema identificato possibile con un sottoinsieme di dati, verrebbe raggiunta comparando tutti i sistemi indentificati scegliendo di volta in volta tutte le possibili combinazioni di ingressi possibili per il sistema stesso. Naturalmente, a causa dell'elevato numero di sensori presenti, non è possibile procedere con questa metodologia: la complessità di un'algoritmo che risolve questo problema sarebbe difatti esponenziale. Nel precedente lavoro, sono state affrontate le tecniche di selezione ingressi/uscite con metodologie greedy o con algoritmi genetici. Queste tipologie di tecniche non forniscono però una soluzione ottima, bensì partono con una soluzione sub-ottima cercando di migliorarla nell'intorno. L'obbiettivo è riuscire a compiere un'identificazione a sottospazi migliore, selezionando i dati in ingresso mediante algoritmi più efficienti come PCR e PLS, che sfruttano strumenti quali l'analisi delle componenti principali.
Per la verifica della bontà di selezione e del relativo modello identificato si sono sviluppati diversi criteri di validazione.
  
% You must have at least 2 lines in the paragraph with the drop letter
% (should never be an issue)
%\hfill Marzo 28, 2008

\section{I Simulatori}

\subsection{Introduzione}

A seguito dell'introduzione della certificazione energetica degli edifici da parte della Comunit\'a Europea allo scopo di ridurne i consumi e l'inquinamento, molti software sono stati sviluppati per calcolarne l'\emph{EPB}, ossia \emph{"l'energy performance of buildings"}. Essi considerano nel calcolo svariati parametri in condizioni statiche o dinamiche, il primo dei quali \'e la posizione geografica. Nel Nord d'Europa il clima \'e freddo durante gran parte dell'anno, per cui l'EPB dipende essenzialmente dall'energia del calore disperso. Di conseguenza le progettazione degli edifici richiede un cospicuo isolamento termico e soluzioni adatte alla massima cattura dell'energia solare per il relativo sfruttamento. Viceversa nei paesi molto pi\'u a sud, dove il clima \'e caldo e umido, le soluzioni tecniche richieste al fine di massimizzare la dispersione del calore sono:
\begin{itemize}
	\item ventilazione passiva,
	\item progetti di climatizzazione e sfruttamento dell'inerzia termica.
\end{itemize}
     Al centro di queste due filosofie di progettazione, si inserisce l'area mediterranea caratterizzata da andamenti climatici intermedi, tali da rendere necessaria una maggiore flessibilit\'a delle soluzioni tecniche, tradotta nella realt\'a in una diversificazione dell'architettura storica di questi paesaggi. 

Lo sviluppo delle tecnologie costruttive agli inizi del 20°secolo non prendevano in considerazione il carattere ener\-getico degli edifici. Solo a partire dagl'anni '50 sono stati introdotti i primi accorgimenti nell'ambito, per esempio:
\begin{itemize}
	\item  riduzione dell'inerzia termica coibentando gli edifici,
	\item riduzione dello spessore dei muri.
\end{itemize}
Queste prime soluzioni, per\'o, si sono dimostrate subito insufficienti con l'avvento della crisi economica del 1973 e con l'adozione di nuove politiche climatiche. Durante gli anni '80 l'Europa a tale proposito ha prodotto una serie di standard costruttivi fino alla recente EN 832. Nello stesso periodo la societ\'a ha subito un forte sviluppo economico, e una seguente modifica della percezione del comfort, sia abitativo che lavorativo.\\ 
Purtroppo questo sviluppo nel benessere abitativo unito ad altri fattori, ha causato l'incremento dei consumi di energia nel settore edilizio, specialmente nell'area mediterranea dove si utilizza molta energia sia per riscaldare durante l'inverno, sia per il raffreddamento d'estate. La Comunit\'a Europea stima che pi\'u del 40\% del consumo di energia sia attribuibile agli edifici. Per una riduzione e una maggiore razionalizzazione delle risorse, sono state emanate leggi e normative all'interno del direttivo 91/2002 EPBD \emph{"Energy Performance of Buildings Directive"}. Le normative tecniche collegate al EPBD sono state sviluppate dal C.E.N.(\emph{"European Committee for Standardization"}) e vanno sotto il nome di CEN-Umbrella. L'EPBD introduce l'obbligo della certificazione energetica degli edifici, il cui approccio \'e molto simile alle etichette di classe energetica delle apparecchiature elettriche. La CEN-Umbrella include una revisione della EN 832, della EN 13790, introducendo nuovi paramteri di valutazione: 

\begin{itemize}
\item l'energia richiesta durante l'inverno e l'estate per:  riscaldamento,  raffreddamento e la ventilazione.
\item la produzione di acqua calda ad uso domestico (DHW- "`Domestic Hot Water"').
\item consumi energia elettrica.
\end{itemize}

In tale contesto i software e i modelli per il calcolo dell'EPB ricoprono un valore importante ai fini di una corretta identificazione della classe, potendo contribuire anche alla ricerca della migliore soluzione per incrementare l'efficienza energetica. Teoricamente si dovrebbe garantire una forte corrispondenza tra i dati simulati e quelli reali, per permettere un altrettanto miglioramento nella fase progettuale e migliorandone il controllo della stessa efficienza anche da parte dell'utente. Evidentemente questo si scontra con la richiesta di software con capacit\'a globali, ossia corrette per diverse aree climatiche e caratteristiche costruttive, e robusti all'aumentare della complessit\'a dei fattori e delle variabili in gioco che non sono comparabili e standardizzabili.

\subsection{Modelli di calcolo e riferimenti standard del simulatore}

Il metodo di calcolo utilizzato dal software DesignBuilder, e adottato nelle simulazioni in questo progetto, segue i riferi\-menti CEN con gli standard: prEN 13790, prEN 15603 e relative normative sugli impianti. Si tenga presente che DesignB. non corrisponde perfettamente allo standard CEN-Umbrella per la mancanza di aggiornamenti su vai coefficienti e nell'utilizzo fonti di energia rinnovabili.

\subsubsection{Descrizione calcoli}

Il software per l'energia termica usa la formula dettata dagli standard europei.

\begin{eqnarray}
	Q_{H,nd} & = & Q_{L,H} - \eta_{g,H} Q_{G,H}
	\label{heat}
\end{eqnarray}

Dove $Q_{H,nd}$ \'e l'energia richiesta per il riscaldamento dell'edificio (in MJ), $Q_{L,H}$ \'e la quantit\'a di calore totale trasferito e $Q_{G,H}$ \'e il calore totale fornito dalle risorse energe\-tiche, infine $\eta_{g,H}$ \'e il fattore di guadagno di utilizzazione adimensionale. 
Per il raffreddamento la formula diviene:

\begin{eqnarray}
	Q_{C,nd} & = & Q_{G,C} - \eta_{l,C} Q_{L,C}
	\label{cool}
\end{eqnarray}

Dove $Q_{C,nd}$ \'e l'energia richiesta per il raffrescamento dell'edificio (in MJ), $Q_{G,C}$ \'e la quantit\'a di raffreddamento totale trasferito e $Q_{L,C}$ \'e il calore totale fornito dalle risorse energe\-tiche, infine $\eta_{g,H}$ \'e il fattore di guadagno di utilizzazione adimensionale. 
Nelle equazioni \ref{heat}, \ref{cool} le lettere H e C fanno riferimento rispettivamente alle modalit\'a di riscaldamento e climatizzazione, mentre i rispettivi addendi per entrambe le equazioni si scompongono in:

\begin{eqnarray}
	Q_{L} & = & Q_{tr} + Q_{ve}
	\label{L}
\end{eqnarray}

\begin{eqnarray}
	Q_{G} & = & Q_{int} + Q_{sol}
	\label{G}
\end{eqnarray}

Dove il calore totale trasferito $Q_{L}$ \'e composto da $Q_{tr}$ calore trasferito per trasmissione, dipendente dal valore della trasmittanza per ogni componente e $Q_{ve}$ per la ventilazione.
In analogia le totali fonti di calore $Q_{G}$ per ogni funzionamento, sono espresse come la somma di tutte le fonti di calore interne sul periodo $Q_{int}$ e la somma delle fonti solari sempre sul periodo $Q_{sol}$, quest'ultimo dipendente dai dati meteorologici. La $Q_{int}$ \'e ottenuta con metodi stocastici, non comparabile con i reali consumi. Questo valore include tutta l'energia derivante da cucine, luce, consumi elettrici e generali dell'unit\'a.

\begin{eqnarray}
	Q_{int} & = & \sum_{k}Q_{int,k} + \sum_{l}(1-b_l)\phi_{int,u,l}
	\label{int}
\end{eqnarray}

Con $Q_{int,k}$ il rate orario di flusso di calore dalla k-esima fonte di calore interna (computer, luci artificiali, apparecchi domestici,.. sia nel caso reale che da tabelle normalizzate), $\phi_{int,u,l}$ rate orario di flusso di calore dalla l-esima fonte di calore interna nello spazio incondizionato W, infine con $b_l(W)$ il relativo fattore di riduzione (ISO 13789). Il calore disperso e la fornitura dello stesso, sono calcolate in relazione con le variazioni della temperatura esterna e con l'irraggiamento medio durante tutto l'anno, basato sui dati meteorologici. 
La valutazione per il calcolo dei consumi DHW sfruttano la formula:

\begin{eqnarray}
	Q_{h,w} & = & V_w \rho C n (t_w - t_0)
	\label{hw}
\end{eqnarray}

dove V \'e il volume di DHW consumato, $\rho$ la densit\'a dell'acqua, C il calore specifico dell'acqua (4186 J/Kg K), $t_W$ la temperatura dell'acqua in uscita e $t_0$ quella in ingresso. Si tenga presente che l'uso di tutte queste standardizzazioni (DHW, luce, energia per cucinare, consumo elettrico e ventilazione naturale) rappresentano un punto molto critico per la modelizzazione e il calcolo. Questo perch\'e sono innanzitutto ideali, estremamente variabili, e non sono completamente valide per la normativa italiana (Appendice \ref{calcolo}).
DHW (\emph{Domestic hot water}) consumata dipende da molteplici fattori: destinazione d'uso, numero di utenti, numero di utenze finali. Le altre udm si riassumono in guadagni di energia interna ottenuti dai computer, apparecchi domestici, usi del gas per la cucina standardizzati con il parametro $W/m^2$. Una nota va' ai valori dei fattori di ventilazione naturale che differiscono anch'essi dalla normativa italiana.

\subsection{DesignBuilder}

Il mondo dell'edilizia stà rapidamente evolvendo verso un'industrializzazione sia della fase costruttiva dei componenti che della costruzione degli edifici nei cantieri. La disponibilit\'a di nuovi strumenti software per la gestione ottimale di questa spinta evolutiva ha in effetti migliorato la qualit\'a complessiva delle costruzioni e ridotto le tempistiche nei cantieri. Non tutto per\'o \'e ancora stato ottimizzato anche perch\'e, nel mondo delle costruzioni, operano diverse discipline non sempre fra loro omogenee. La differenza maggiore di impostazione riguarda il mondo degli impianti, la cui importanza, alla luce dei costi energetici, stà crescendo in modo esponenziale. Tanto che la qualit\'a energetica di un edificio sar\'a l'elemento discriminante, a breve, nelle compravendite di immobili. In questa integrazione si inserisce il software DesignBuilder suggerito per le simulazioni dell'edificio di questo progetto.
DesignBuilder nasce da un progetto europeo e si propone come uno strumento di grande qualit\'a per analisi energetiche serie e attendibili di edifici semplici o estremamente complessi.
\'E un software che unisce alla rapidit\'a e alla facilit\'a di rappresentazione tridimensionale dell'edificio, l'analisi energetica del sistema Edificio/Impianto che si basa sul motore di simulazione energetica "`Energyplus"' \footnote{EnergyPlus \'e il motore di calcolo energetico sviluppato dal ministero degli Stati Uniti(DOE) a seguito della prima crisi energetica del 1974. Nato alla Berkely University continua lo sviluppo anche in organizzazioni scientifiche e tecniche come la ASHARE che ne valida le procedure di calcolo tramite un apposito protocollo internazionale}. Questo software fornisce una serie di dati di rendimento ambientali come: consumo annuo di energia, temperature massime estive e potenze dei componenti \emph{HVAC (Heating Ventilation and Air Conditioning)}. Il software pu\'o essere utilizzato efficacemente in qualunque fase del processo di disegno dalla fase di concetto dove si pu\'o voler provare svariate tipologie di strutture in ogni singola stanza, oppure attraverso tutta la fase di realizzazione. L'utilizzo di "`Template"' e dati di "`Default"' permettono di  guidare l'utente attraverso il processo di inserimento dei dati, delle analisi e delle varianti in corso d'opera. I modelli di dati (Template) permettono all'utente di caricare geometrie, strutture e impianti di edifici tipici, ottenendo immediatamente calcoli attendibili per valutazioni di massima. Successivamente si pu\'o passare a fasi di dettaglio, modificando opportunamente i singoli dati, per definire al meglio i valori di potenza e di consumi dello specifico progetto. \'E inoltre possibile aggiungere i propri modelli se si lavora spesso su tipologie di edifici simili. L'interazione con la vista tridimensionale renderizzata permette all'utente di capire immediatamente la differenza tra architettura e altri dati pratici come il surriscaldamento di porzioni non climatizzate dell'edificio. \'E dunque adatto all'utilizzo di architetti, ingegneri edili, termotecnici e studenti.

%%%%%%%%%%%%%%%%%%%%%%%%%%% fine simulatori %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{L'impostazione dell'esperimento}

Tenendo in considerazione il lavoro svolto lo scorso anno \cite{baldan}, e con l'obiettivo di comparare i risultati ottenuti si sono utilizzati come dati sperimentali quelli raccolti in \cite{baldan}. Si ritiene utile tuttavia fornire indicazioni su come impostare al meglio la fase di test e raccolta dati sul campo, in modo da permettere una successiva fase di identificazione che sia la pi\'u efficace possibile. \'E possibile definire in maniera adatta le modalit\'a e il set di dati da raccogliere nell'esperimento principale \cite{manual}.\\
Le fasi preliminari si possono distinguere in:
\begin{itemize}
\item \textbf{Tempo invarianza} Molti metodi di identificazione assumono modelli tempo-invarianti e gli esperimenti/misure riescono a determinare tale propriet\'a. Per\'o \'e possibile che vi sia una valutazione errata. Infatti delle dinamiche lente appaiono come dinamiche tempo varianti se consideriamo una scala dei tempi troppo corta. \'E necessario qualora si voglia identificare delle dinamiche veloci, eliminare quelle lente assimilandole a trend.
\item \textbf{Linearit\'a} In Matlab, il toolbox di identificazione contiene routine per identificare solo sistemi dinamici. La linearit\'a pu\'o essere verificata attraverso l'analisi delle risposte del sistema a differenti ampiezze degli ingressi, con la simmetria della risposta ad un'onda quadra o derivando propriet\'a spettrali.\\ Essendo molti sistemi non lineari, \'e comunque possibi\-le, sotto certe ipotesi, ricavarsi la linearizzazione del model\-lo. 
\item \textbf{Transitori} L'analisi dei transitori attraverso risposte all'impulso o al gradino, fornisce valutazioni su constanti di tempo, ritardi, guadagni stazionari. 
\item \textbf{Risposta in frequenza} Sono metodologie relegate alla stima non parametrica, in comune con l'analisi delle cor\-relazioni. Esse forniscono la stima della funzione di trasferi\-mento del sistema. In tale contesto il diagramma di Bode si rivela utile per le informazioni sulla locazione di poli e zeri, guadagni e ritardi. 
\end{itemize}

Per quanto riguarda la parte principale dell'esperimento, si sa' che l'identificazione fornisce accurati modelli alle frequenze dove il segnale di ingresso contiene maggiore energia. Questo si traduce in una sufficiente eccitazione dell'ingresso a tali frequenze. Allora \'e naturale cercare di concentrare le frequenze del segnale di ingresso in modo tale da rendere piccolo l'errore di stima dove ci\'o si rende pi\'u necessario. Una scelta comune \'e l'ingresso PRBS \emph{pseudo-random binary sequence} che ha come caratteristica di avere molta energia su uno spettro ampio. Se questa scelta non fosse possibile, come in questo caso di ingresso esterno, \'e auspicabile avere un segnale con un buon rapporto segnale/rumore e una condizione di persistente eccitazione sufficiente ad una buona identificazione.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  L'identificazione   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{L'identificazione}

Per quanto riguarda l'identificazione di modelli a variabili di stato, il procedimento ideale che si basa sulla fattorizzazione della matrice di Hankel, raramente porta a risultati soddisfacenti nelle applicazioni con dati reali. Questo tipo di identificazione comporta generalmente grossolani errori di stima relazionati a:
\begin{itemize}
	\item non è possibile disporre di dati esatti, conseguenza di disturbi nelle misure (errori di misura),
	\item inevitabili arrotondamenti nelle elaborazioni numeriche (approssimazioni numeriche),
	\item il modello lineare a dimensioni finite ipotizzato per spiegare i dati, si dimostra inadatto a descrivere la realtà (semplificazioni di modello).
\end{itemize}
 Questa serie di cause comporta un aumento di rango della matrice di Hankel, e con esso una diretta sovrastima del modello a fronte di un complessità reale ben minore (\cite{Bittanti} p.132]). Questo cruciale problema trova un adeguato rimedio nella scomposizione matriciale di Beltrami, nota come scomposizione ai valori singolari (SVD). Si scompone la matrice di Hankel $H$ in tre matrici reali: $U$ ($s$ x $s$ con $s$ le righe di $H$), $V$ ($r$x$r$ con $r$ colonne di $H$) e $S$ matrice diagonale con medesime dimensioni di $H$. Si tenga presente che la nozione di valore singolare non và confusa con quella di autovalore di una matrice, poichè i primi sono tutti reali non negativi, possono essere associati a matrici con dimensione qualsiasi e sono di facile e robusto computo.
Un'altra caratteristica significativa della scomposizione ai valori singolari \'e che la determinazione del rango coincide con il numero di autovalori singolari non nulli della matrice $S$. \'E possibile capire cos\'i l'ordine del sistema o l'ordine della realizzazione minima. Tenendo poi presente che i campioni della risposta impulsiva sono affetti da rumore, gli stessi valori singolari si presentano in soprannumero ai fini identificativi, per cui in tale contesto vengono classificati in due famiglie: valori singolari dominanti e secondari. La stima dell'ordine $n$ viene associata ai valori dominanti, una scelta che in alcuni casi si pu\'o presentare soggettiva, a causa di una non evidente suddivisione nei due insiemi. Una volta scelto l'ordine si procede al calcolo delle matrici $A,B,C,D$ del sistema, per esempio con il metodo a sottospazi di Kung (\cite{Bittanti}). Vi sono poi diverse altre tecniche descritte in questa relazione di identificazione a sottospazi, non basate sull'identificazione della risposta impulsiva, bens\'i a partire da dati ingresso-uscita generici, e supportati sempre da tecniche numeriche. Sono il PBSID, PBSID\_opt fornito dal prof.Chiuso, N4SID \footnote{Una tecnica di identificazione a sottospazi generica \'e descritta in \cite{Bittanti}}.

\subsection{Sottospazi e catena chiusa}

I metodi subspace identification trovano il loro utilizzo  nell'identificazione di sistemi MIMO dinamici lineari a tempo discreto e dimensione finita di grandi dimensioni, dove i metodi a minimizzazione dell'errore di predizione sarebbero proibitivi per motivi tecnici. Da una prima analisi dei dati, valutati con l'uso degli strumenti forniti da Matlab, si evince quanto gli ingressi e le uscite in questo progetto siano fortemente correlati. A tale proposito \'e stato dunque necessario introdurre una estensione degli algoritmi di identificazione a sottospazi per operare in condizioni di anello chiuso, il cui studio ha subito uno sviluppo solo recentemente. Questa nuova classe di algoritmi consistenti per la catena chiusa, si basa sull'identificazione del predittore, in maniera molto simile al metodo ad errore di predizione ($PEM$), avvicinandosi in alcuni casi al suo comportamento. Le condizioni richieste per questo tipo di identificazione hanno richiesto di colmare un vuoto nelle procedure di realizzazione stocastica in presenza di feedback, cercando di eliminarne i termini indesiderati. 
Le procedure per l'identificazione con i metodi a sottospazi sono composte essenzialmente da due step:

\begin{itemize}
\item[1)] costruzione di una base per lo spazio di stato attraverso una proiezione adatta delle sequenze dei dati.
\item[2)] data questa base, risolvere al tempo $t$ e $t+1$ nel senso dei minimi quadrati ("least squares") il sistema:

          \begin{equation}
						\left\{
						\begin{array}{rcl}
							\hat{X}_{t+1} & \simeq & A\hat{X}_t + B\hat{U}_t + KE_t\\
							 Y_{t}        & \simeq & C\hat{X}_t + E_t
						\end{array}\right.
						\end{equation}
\end{itemize}

In parole povere, i metodi subspace permettono di ricavare gli stati del sistema a partire dai dati I/O attraverso opportune proiezioni ortogonali o oblique degli sottospazi generati, defi\-nite in questi casi da matrici di Hankel a blocchi. Dopo di che, a partire dalla proiezione attraverso una decomposizione ai valori singolari (SVD) si ottiene una stima dell'ordine del sistema, la matrice di osservabilità estesa e la sequenza degli stati. Da questo si ricavano le matrici A,B,C,D come soluzione esatta o ai minimi quadrati di un opportuno sistema,  e le matrici di covarianza Q,S,R dai residui. Si tenga presente che gli algoritmi a sottospazi differiscono essenzialmente sull'implementazione del primo punto, mentre il secondo \'e virtualmente uguale per tutti. Nelle prossime sezioni si descrivono qualitativamente gli algoritmi $PBSID$ e la sua versione ottimizzata $PBSID_{opt}$ data l'utilit\'a che ha avuto in questo progetto (vedi \cite{2007b} per una trattazione esaustiva). 

\subsection{$PBSID$}

Il PBSID ("`Prediction based subspace identification"') a differenza degli altri metodi a sottospazi tipo N4SID, utilizza differenti matrici di Hankel a blocchi, per operare in presenza di feedback. Si cerca cos\'i di risolvere il problema, introdotto dagli ingressi colorati, della non eliminazione del rumore additivo stocastico nella proiezione obliqua. Si utilizzano poi queste proiezioni oblique per ottenere una base per lo stato e la matrice di osservabilit\'a estesa; segue una decomposizione SVD per separarle e la risoluzione ai minimi quadrati per ricavare le matrici del sistema. 

Dato il sistema lineare stocastico nella forma di in\-novazione:

\begin{equation}
 \left\{
	\begin{array}{rcl}
			  \textbf{x}(t+1) & = & A\textbf{x}(t) + B\textbf{u}(t) + K\textbf{e}(t)\\
				\textbf{y}(t)   & = & C\textbf{x}(t) + D\textbf{u}(t) + \textbf{e}(t)
	\end{array}\right.
	\label{stoc}
\end{equation}
\\
Per la trattazione si assume che:

\begin{itemize}
\item presenza di feedback tra I/O. 
\item la rappresentazione \'e minima.
\item supposto D=0.
\item si definisce $\overline{A} := A - KC$.
\item il processo ergodico stazionario del secondo ordine $\textbf{z}(t) := [\textbf{y}^\top \textbf{u}^\top]^T$ abbia densit\'a spettrale razionale e zeri tutti all'interno del cerchio 					unitario 
\item Il processo di rumore bianco $\textbf{e}$, di innovazione di $\textbf{y}$ dato il passato di $\textbf{z}$, \'e definito come l'errore di predizione di un passo di 											$\textbf{y}(t)$  dato il passato stretto di \textbf{z}(t) al tempo $t$.
\item la differenza tra l'istante presente $t$ e quello iniziale $t_0$ tende all'infinito rendendo trascurabile quest'ultimo termine. \'E un assunto necessario per assicurare la 						consistenza dello stimatore in catena chiusa,  e la stima di infinite covarianze  per ottenere l'asintotica efficienza degli stimatori, anche se questo mediamente pu\'o non 						essere automatico.
\item per semplificare l'analisi nell'articolo \cite{2007b}, si adotta un ordine $\nu$ di orizzonte futuro fissato e finito, sufficientemente grande.
\end{itemize}

Durante la trattazione, si useranno diversi strumenti la cui trattazione completa esula dallo scopo prefisso. La proiezione obliqua \'e la scomposizione in due sottospazi tra loro non ortogonali. Si denota come $\hat{E}_{||B}[\cdot|C_t]$, ossia la proiezione lungo lo spazio generato delle righe di $B$ sullo spazio generato dalle righe di $C_t$. La formula estesa \'e:

\begin{equation}
 \hat{E}_{||B_t}[A_t|C_t] = \hat{\Sigma}_{ac|b}\hat{\Sigma}_{cc|b}^{-1}C_t
\label{obl}
\end{equation}

Dove sotto le ipotesi di ergodicit\'a le matrici di covarianza campionaria (condizionata e non) sono:

\begin{equation}
\begin{array}{lll}
  \hat{\Sigma}_{ab} &:=& \frac{1}{N}\sum_{i=0}^{N-1}a_{t+i}b_{t+i}^{\top}\\\\
  \hat{\Sigma}_{ab|c} &:=& \hat{\Sigma}_{ab} - \hat{\Sigma}_{ac}\hat{\Sigma}_{cc}^{-1}\hat{\Sigma}_{cb}
   \label{sigma}
\end{array}
\end{equation}

Infine si definisce la matrice di osservabilit\'a estesa:

\begin{equation}
 \overline{\Gamma}_{\nu}^{\top} := [C^{\top} \quad \overline{A}^{\top}C^{\top} \quad \ldots \quad (\overline{A}^{\top})^{\nu-1}C^{\top}]
\label{obv}
\end{equation}


\subsubsection{L'algoritmo}

L'algoritmo PBSID per la ricerca di una base per lo spazio di stato usa diverse proiezioni oblique. La proiezione di ogni riga dell'uscita $Y_{t+h} \quad h=0,...,\nu$ pu\'o essere vista con un modello VARX del tipo:

\begin{eqnarray} \label{varx}
	\begin{array}{rcl}
	\hat{Y}_{t+h} &:=& \hat{E}[Y_{t+h}|Z_{[t_0,t+h)}]\\\\
	              &=& \hat{\Psi}_{1,h}Z_{t+h-1} + \ldots + \hat{\Psi}_{t+h-t_0,h}Z_{t_0}
	\end{array}
\end{eqnarray}\\

dalla quale attraverso l'algoritmo "`predictor-based"'(la P in apice) le proiezioni oblique asintoticamente in $N$, forniscono:

\begin{eqnarray} \label{varx1}
	\begin{array}{rcllll}
	\hat{Y}_{t+h}^P &:=& \hat{E}_{||Z_{[t,t+h)}}[Y_{t+h}|Z_{[t_0,t)}]\\\\
	                 &=& \sum_{i=h+1}^{t-t_0+h} \hat{\Psi}_{i,h}Z_{t+h-i} & \simeq & C\overline{A}^{h-1}X_t
	\end{array}
\end{eqnarray}

Tutti i predittori fino a $h=T$ con tale approssimazione si possono raggruppare come:

\begin{equation} \label{varx2}
  \begin{array}{cccccc}
     \hat{Y}_{[t,T)}^P &:=& \left[%
     	\begin{array}{c}
     		\hat{Y}_t^P \\\\
     		\hat{Y}_{t+1}^P \\\\
     		\vdots\\\\
     		\hat{Y}_{T-1}^P \\
     	\end{array}
     \right]
     & \simeq & \overline{\Gamma}_{\nu}X_t\\  
	\end{array}
\end{equation}

Dalla decomposizione ai valori singolari:

\begin{eqnarray}
			W_p^{-1}\hat{Y}_{[t,T)}^P &:=& PDQ^{\top}\\
			&=& \left[P_n \quad \tilde{P}_n \right]
     \left[
     	\begin{array}{cc}
			D_n & 0\\
			0 & \tilde{D}_n
     	\end{array}
     \right]
      \left[ Q_n^{\top} \quad \tilde{Q}_n^{\top} \right]
      \label{varx3}
\end{eqnarray}\\
  
con una opportuna matrice peso $W_p$, \'e possibile stimare una matrice di osservabilit\'a, ottenuta ignorando i valori singolari che descrivono meno il sistema (ossia ponendo nell'equazione \ref{varx3} la $\tilde{D}_n \simeq 0$) \footnote{questa analisi \'e ripresa nella descrizione del metodo $PCA$}. Allora:

\begin{equation}
\hat{\overline{\Gamma}}_{\nu} = W_PP_nD_n^{1/2} 
\end{equation}

quindi dalla \ref{varx2} e sfruttando le equazioni precedenti si trova una base per lo spazio di stato

\begin{equation}
\hat{X}_t^{PBSID} := \hat{\overline{\Gamma}}_{\nu}^{-L}\hat{Y}_{[t,T)}^P 
\end{equation}

\begin{equation}
\hat{X}_{t+1}^{PBSID} := \hat{\overline{\Gamma}}_{\nu}^{-L}\hat{Y}_{[t+1,T)}^P 
\end{equation}

Dove $\hat{\overline{\Gamma}}_{\nu}^{-L}$ \'e l'inversa sinistra definita come

\begin{equation}
\hat{\overline{\Gamma}}_{\nu}^{-L} := (\hat{\overline{\Gamma}}_{\nu}^{-L}W_p^{-\top}W_p^{-1}\hat{\overline{\Gamma}}_{\nu})^{-1}\hat{\overline{\Gamma}}_{\nu}^{-L}W_p^{-\top}W_p^{-1}
\end{equation}\\

\subsection{$PBSID_{opt}$}

Il metodo $PBSID_{opt}$ \'e numericamente equivalente alla stima di un modello VARX, utilizzante i medesimi passi per l'identificazione a sottospazi e con il pregio di avere una minore complessit\'a computazionale. Unisce cos\'i i pregi dello stimatore consistente CCA per la catena aperta e migliora il comportamento in retroazione rispetto a modelli SSARX/PBSID, nel senso della asintotica varianza. Infatti il PBSID \'e asintoticamente equivalente al CCA ("Canonical Correlation Analysis") per l'identificazione di serie di dati e in presenza di ingressi bianchi. Mentre la sua versione ottimizzata fornisce una minore o uguale asintotica varianza per ogni scelta del segnale di ingresso, anche colorati, (dimostrato in \cite{2007}~). Comunque non \'e ancora del tutto chiaro il comportamento dell'efficienza in generale. Questa versione ottimizzata per i sistemi MIMO, presenta migliori prestazioni di identificazione in sistemi con contenuto alle basse frequenze (come le temperature in un edificio), la sostituzione del calcolo esplicito delle proiezioni oblique con passaggi di minimizzazione iterativa, la velocit\'a di calcolo e la produzione di modelli stabili e non affetti da bias. La selezione automatica dell'algoritmo dell'ordine del modello segue il criterio di rapporto tra valori singolari adiacenti e tra il valore massimo e minimo. Per quanto riguarda l'algoritmo, esso differisce dal PBSID originale, per il fatto che la stima del predittore \ref{varx} \'e formulato come un problema ai minimi quadrati pesati, dove la proiezione equivale, dopo opportune manipolazioni, a risolvere riga per riga il seguente sistema vettoriale:

\begin{eqnarray*} 
     Y &:=& \left[%
     	\begin{array}{c}
     		vec(Y_t)\\\\
     		vec(Y_{t+1}) \\\\
     		\vdots\\\\
     		vec(Y_T)
     	\end{array}
     \right]
      =  S^P\Omega^P +  
     \left[%
     	\begin{array}{c}
     		vec(E_t)\\\\
     		vec(E_{t+1}) \\\\
     		\vdots\\\\
     		vec(E_T)
     	\end{array}
     \right]\\\\
     &=& S^P\Omega^P+E
\end{eqnarray*}

Dove si sono definite le nuove matrici:

\begin{equation}
\begin{array}{lll}
  \Xi_i &=& C\overline{A}^i\Delta \quad \forall i=1,\ldots,\nu \\
  \Delta &:=& [\overline{A}^{(t-t_0-1)}[K B]\overline{A}^{(t-t_0-2)}[K B]\ldots[K B]] \\
  \Psi_{ij} &=& C\overline{A}^{j-1}[K B] \quad \forall i,j=1,\ldots,\nu 
 \label{matrici}
\end{array}
\end{equation}

A noi basti sapere che attraverso uno stimatore di Markov si trova iterativamente una soluzione ottima $\hat{\Omega}^{Popt}$ come stima di $\Omega^P$ che comporta la pi\'u piccola varianza asintotica basata sul modello. 

\begin{equation}
 Y := S^P\Omega^P+E
 \label{sistema}
\end{equation}

Con la soluzione $\hat{\Omega}^{Popt}$ di \ref{sistema} \'e possibile rivedere i passaggi dell'algoritmo PBSID precedente, pervenendo al nuovo stato:

\begin{equation} \label{varx5}
  \begin{array}{cccccc}
     \hat{X}_{t+1}^{Popt} &:=& (\hat{\overline{\Gamma}}_{\nu}^{Popt})^{-L} \left[%
     	\begin{array}{cc}
     		\hat{\Xi}_1^{Popt} & \hat{\Psi}_{11}^{Popt}\\\\
     		\hat{\Xi}_2^{Popt} & \hat{\Psi}_{22}^{Popt}\\\\
     		\vdots & \vdots\\\\
     		\hat{\Xi}_{\nu}^{Popt} & \hat{\Psi}_{\nu\nu}^{Popt}\\\\
     	\end{array}
     \right]
     Z_{[t_0,t+1)}\\  
	\end{array}
\end{equation}

Dove $\hat{\overline{\Gamma}}_{\nu}^{Popt}$ \'e la matrice estesa di osservabilit\'a stimata (sostituendo $\hat{Y}_{t+k}^{Popt}$ a $\hat{Y}_{t+k}^{P}$ nella decomposizione \ref{varx3}).
Similmente \'e possibile stimare la sequenza $E_t$, come:

\begin{equation}
 \hat{E}_t^{Popt} := Y_t - \hat{E}[\hat{Y}_t^{Popt}|\hat{X}_t^{Popt}]
 \label{E}
\end{equation}

\section{Il problema del modello nella regressione lineare\\}

\PARstart{N}{}ei problemi di controllo di sistemi dinamici, uno dei passi fondamentali è quello di poter ottenere un modello matematico che descriva il comportamento e le relazioni tra le variabili di ingresso o predditive e le variabili di uscita. Nel caso ideale il modello matematico è disponibile come set di equazioni ai principi primi esplicitato mediante equazioni di bilancio di materia, e di energia in cui i parametri sono coefficienti reali del processo, aventi significato fisico.\\
Nella stragrande maggioranza dei casi invece non si dispone del modello deterministico, questo significa che non sono noti o ricavabili alcuni parametri del modello. In questi casi si dispone invece di una collezione di misure delle variabili di interesse o di alcune grandezze ad esse collegate e si cerca di approssimare il reale comportamento del processo mediante un'equazione empirica di \emph{regressione}.\\
Le tecniche di analisi statistica multivariata basate sul concetto di variabile latente forniscono gli strumenti necessari per estrarre da una moltitudine di dati le informazioni necessarie separandole dal rumore e dai disturbi presenti nel processo.\\
L'applicazione dei metodi di regressione ed analisi statistica multivariata permette di identificare o meglio stimare la relazione esistente tra i due set di variabili, di ingresso e di uscita; in questo modo si ottiene un modello matematico empirico del processo che ci permette di progettare un sistema di controllo. Tali tipi di modelli vengono denominati in letteratura come modelli black-box.\\
L'identificazione del modello è quindi una parte importante e ingegniristicamente onerosa nel controllo di qualunque sistema, come indicato nell'abstract noi cercheremo di identificare un modello termodinamico che descriva l'evolversi della temperatura all'interno di un edificio sotto l'effetto della temperatura esterna e della radiazione luminosa che interessa l'edificio su cui sono state svolte le misure sperimentali.
Cercheremo di fornire un quadro di insieme della teoria della regressione lineare multivariata per capire bene dove l'applicheremo nella soluzione del nostro problema evidenziandone i punti di forza che presenta.
Presentiamo due metodi di regressione lineare multipla:
\begin{itemize}
	\item la regressione ai minimi quadrati classici \emph{Ordinary Least Squares Regression};
	\item la regressione sulle componenti principali \emph{Principal Component Regression};
	\item la proiezione su Strutture Latenti \emph{Projectin on Latent Structures}.
\end{itemize}

\subsection{Modelli Lineari Multivariati}
Supponiamo di disporre di $n$ misure sperimentali di $k$ variabili indipendenti $x_{1},\ldots,x_{k}$ e di altrettante misure della variabile dipendente $y$, con \mbox{$n > k$}. Un  modello multivariato coinvolge una o più variabili dipendenti, rappresentate come funzione lineare di più variabili indipendenti o predittive.\\
Il modello di equazione lineare multivariata è:
   
\begin{equation}
\label{eq1}
\begin{array}{ccc}
 y_{i} &=& \alpha_{0} + \alpha_{1}x_{i1} + \ldots + \alpha_{k}x_{ik} + \epsilon_{i} \\
				&=&  \sum^{k}_{j=1}\alpha_{j}x_{ij} + \epsilon_{i} \quad \forall i = 1,\ldots,n
\end{array}
\end{equation}

dove $y_{i}$ è l'i-esimo valore della variabile dipendente, $x_{1},\ldots,x_{k}$ sono le $k$ variabili indipendenti e $\alpha_{0},\ldots,\alpha_{k}$ sono $k+1$ parametri matematici o fisici il cui valore viene stimato con metodi di regressione statistica in modo da minimizzare $\epsilon_{i}$ che è un errore stocastico sistematico che tiene conto del rumore di misura e di altre cause erratiche nella costruzione del modello.\\
 L'applicazione del modello di regressione lineare ai dati di processo è resa complicata da numerose cause: sia la variabile dipendente che le variabili preditive sono affette da errore di misura, le variabili sono molto correlate tra loro, i dati sono imprecisi o in alcuni casi mancanti, a causa della rottura o malfunzionamento di un sensore. Inolre il rapporto segnale/rumore delle variabili può essere molto basso, introducendo rumore sistematico all'interno dell'equazione di regressione: questo fatto è molto importante, perchè induce nell'equazione di regressione, attraverso i coefficienti di regressione, a modellare  anche l'errore misurato dalle variabili predittive e dipendenti per minimizzare l'errore $\epsilon_{i}.$\\
Ciò significa che anche un'equazione di regressione statisticamente significativa che ottenga un buon fitting dei dati sperimentali e minimizzi l'errore in (\ref{eq1}) può essere non utilizzabile se prima non viene separata l'informazione utile dal rumore.

\subsection{Regressione lineare ai minimi quadrati}
\label{sec:OLS}
La regressione ai minimi quadrati (\emph{Ordinary Least Squares}) risolve un sistema di $n$, equazioni lineari in cui sono incogniti $k+1$ coefficienti di regressione $\alpha_{0}\ldots\alpha_{k}$ corrispondenti ad ogni variabile predittiva più l'intercetta, cercando di minimizzare la varianza dell'errore  $\epsilon$.\\
La regressione ai minimi quadrati è possibile solo se il numero di campioni disponibili è superiore al numero di variabili predittive;~ infatti se $n<k$ il sistema non è risolvibile in modo univoco, ma si hanno infinite soluzioni. Nel caso invece si abbia $n=k$ la soluzione esatta nel sistema lineare è possibile ed è unica: si ottiene per sostituzione o mediante algoritmi numerici come la decomposizione di Gauss.
L'equazione di regressione multivariata  ai minimi quadrati è data da:
	\[ \hat{y_{i}}= b_{0} + b_{1}x_{i1} + \ldots + b_{k}x_{ik} = b_{0} + \sum^{k}_{j=1}b_{j}x_{ij}\quad \forall i = 1,\ldots,n \] 
Dove $b_{0}\ldots b_{k}$ sono le stime dei coefficienti di regressione parziale ottenute risolvendo il sistema di equazioni normali ricavate differenziando l'errore totale rispetto ad ogni incognita $\alpha_{0}\ldots\alpha_{k}$.
Il residuo $e_{i}$ tiene conto dell'imprecisione che si introduce utilizzando valori stimati non esatti per i coefficienti di regressione.\\
I residui, nel caso il modello sia esatto, hanno queste caratteristiche:
\begin{itemize}
	\item si distribuiscono normalmente e non sono tra loro autocorrelati; 
	\item hanno media nulla;
	\item sono non correlati con la variabile calcolata: $cov(e_{i},\hat{y_{i}})$
\end{itemize}
Il metodo di regressione ai minimi quadrati pùo essere applicato anche al caso di più variabili ottenendo un'espressione matriciale.Si ha dunque:
\begin{equation}
\label{sistlineare}
	\mathbf{Y} = \mathbf{X}\mathbf{B}+\mathbf{E}
\end{equation}
Il sistema di equazioni normali del metodo di regressione ai minimi quadrati si scrive: 
	\[ \mathbf{X}^{T}\mathbf{X}\mathbf{B} = \mathbf{X}^{T}\mathbf{Y}\] 
la cui soluzione è 
\begin{equation}
 \label{PEM}
	\mathbf{\hat{B}} = (\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{Y} 
\end{equation}
 
Abbiamo cosi l'equazione di regressione lineare multivariata ai minimi quadrati in forma matricile:
	\[ \hat{\mathbf{Y}} = \mathbf{X}(\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{Y}\] 
Tale soluzione minimizza la funzione MSE, dove le stime $b_{j}$ dei coefficienti di regressione parziali $\alpha_{j}$	si caratterizzano perchè:
\begin{itemize}
	\item minimizzano MSE indipendentemente dalla distribuzione degli errori $e_{i}$;
	\item più la deviazione standard è piccola più  la stima dei coefficienti è precisa, se il modello soffre di multicollinearità l'errore standard è molto alto ;
	\item sono stime attendibili di $\alpha_{j}$ se le variabili predittive sono tra loro non correlate o se i residui $e_{i}$ sono statisticamente indipendenti.
\end{itemize}
La soluzione del problema ai minimi quadrati è possibile solo se la matrice $(\mathbf{X}^{T}\mathbf{X})^{-1}$ esiste, cioè se la matrice $X$ ha rango massimo. Nel caso le variabili siano tra loro correlate e non sia possibile effettuare l'inversione della matrice $\mathbf{X}^{T}\mathbf{X}$, si deve rinunciare ad ottenere una stima ai minimi quadrati dei coefficienti di regressione, utilizzzando metodi alternativi per la risoluzione di sistemi lineari di equazioni come la decomposizione QR.L'ideale sarebbe dunque disporre di dati linearmente indipendenti:in questo modo, non solo si avrebbe una soluzione, ma la soluzione trovata sarebbe molto più consistente. Basti pensare che in (\ref{PEM}), se i dati fossero linearmente dipendenti, l'espressione $(\mathbf{XX}^{T})^{-1}$ fornirebbe dei valori molto elevati, quindi i parametri identificati avrebbero poca consistenza, dato che nel sitema piccole variazioni di $\mathbf{X}$ in (\ref{sistlineare}) darebbero forti variazioni delle uscite $\mathbf{Y}$, invalidando l'identificazione effettuata.
Va rilevato infine che la regressione ai minimi quadrati classici si basa su un'assunzione che viene raramente rispettata dai dati di processo, ovvero l'esattezza delle misure sperimemtali delle variabili predittive: l'essenza di questo requisito fa cadere la validità dell'equazione di regressione ai minimi quadrati, perchè il modello non è in grado di verificare statisticamente l'esattezza delle stime dei coefficienti di regressione e l'attendibilità delle stime con dati di convalida.

\section{Metodi di regressione multivariati}
I metodi di regressione multivariata sono strumenti di regressione nati dall'esigenza di superare i limiti della regressione ai minimi quadrati nel trattare variabili predittive fortemente correlate tra loro e affette da rumore. Le tecniche statistiche multivariate più conosciute sono la \emph{Regressione delle Componenti Principali} (PCR) e la \emph{Proiezione su Strutture Latenti} (PLS), che sfruttano entrambe la Analisi delle Componenti Principali (PCA).
I metodi di regressione multivariata vengono definiti anche come \emph{Latent variable Multivariate Regression} in quanto presuppongono esista una struttura sottostante i dati descritta dalle variabili latenti; i metodi PLS e PCA sono strumenti per catturare questa struttura (Mggiori dettagli in\cite{Burnham}).  
compito dei metodi PCA e PLS é:
\begin{enumerate}
	\item  determinare la dimensione effettiva del sistema e cogliere la struttura delle correlazioni tra variabili, questo avviene mediante la definizione delle variabili latenti, tali variabili sono in  numero inferiore alle variabili originali e tra loro ortogonali, esse costituiscono le direzioni del nuovo iperpiano, di dimensione ridotte rispetto alle dimensioni del processo, ma comunque in grado di descrivere i dati.
	\item l'estrazione delle informazioni utili dai dati disponibili e la proiezione di tali informazioni lungo le direzioni prestabilite dalle variabili latenti.
\end{enumerate}
Rispetto ad un metodo OLS classico, i metodi PCA e PLS presentano molti vantaggi,(vedi \cite{kourti},\cite{Kresta},\cite{Tonizzo})
\begin{enumerate}
	\item Sono capaci di gestire l'inclusione nel modello di variabili altamente correlate tra loro. In questo modo si evita di dover escludere una variabile importante dal modello solo perchè fortemente correlata con un'altra variabile già compresa nel modello, perdendo parte dell'informazione essenziale di un processo. Questo è un vantaggio rispetto al metodo OLS, che invece è incapace di gestire variabili tra loro correlate a causa della non invertibilità dell'inversa di una matrice singolare.
	\item riducono le dimensioni del problema
	\item permettono la regressione anche con dati di variabili predittive affette da errore, cosa non possibili con la regressione ai minimi quadrati classichi,
	\item minimizzano il rumore presente nelle variabili latenti mediando il rapporto segnale/rumore proveniente dalle variabili preditive.
	\item forniscono stime più robuste dei coefficienti di regressione parziale.
\end{enumerate}


\section{Analisi delle Componenti Principali}
\label{sec:PCA}
Il problema principale che si può verificare in una procedura di identificazione che sfrutta i minimi quadrati è che i dati in ingresso formino una famiglia linearmente dipendente. In queste circonstanze, l'inversa di $\mathbf{X}\mathbf{X}^{T}$ in \ref{PEM} non esiste; in realtá i problemi sono gli stessi anche nel caso in cui la matrice $\mathbf{X}\mathbf{X}^{T}$ è invertibile ma si è in presenza di forti dipendenze. In questi casi, i parametri identificati che risultano non sono rappresentativi del sistema, e portano a un sistema che non supera la validazione o addirittura instabile in presenza di set di dati diversi da quelli utilizzati in fase identificativa.
Una tecnica che ci permette di evitare tali problematiche è l'analisi delle componenti principali (\emph{PCA}).
Questa tecnica si basa sull'analisi degli autovettori e degli autovalori della matrice di convarianza dei dati, oppure equivalentemente sulla decomposizione a valori singolari (SVD).
Tale metodologia sará utilizzata in seguito come strumento di analisi su dati multivariati. PCA puó esser considerato infatti uno strumento per analizzare la struttura di un insieme di dati multivariati, in particolare con il proposito di ridurne la dimensione: infatti, tramite PCA é possibile prendere i dati e proiettarli/ruotarli in uno spazio di dimen\-sione minore, selezionando le direzioni nello spazio dei dati che forniscono la massima variazione, o equivalentemente, la maggior informazione\footnote{Questo tipo di approccio sará utilizzato in seguito per evitare vari problemi indotti dalla \textit{collinearitá} dei dati.}.
In breve, si ha un vettore
\[
  \mathbf{X}=\left\{ \mathbf{x_{1}};\mathbf{x_{2}} ... \mathbf{x_{m}} \right\}
\]
in cui ogni componente $\mathbf{x_{i}}$ ha dimensione $n$, e la matrice $\mathbf{X}$ non é a rango pieno.
Si vuole ridurre la dimensione $m$ di $\mathbf{X}$, ovvero trovare una trasformazione lineare:
\[
  \mathbf{KLT}: \mathbb{R}^{m x n} \rightarrow \mathbb{R}^{l x n}
\]

\[
  \mathbf{X}\longmapsto \bar{\mathbf{X}}
\]
in modo che la matrice $\bar{\mathbf{X}}$ contenga le stesse informazioni di $\mathbf{X}$ ma sia di rango pieno, mettendo in evidenza le differenti caratteristiche tra i vettori delle misure.\\
Si richiamano in seguito due importanti strumenti teorici, necessari a comprendere e giustificare la metodologia PCA, qui brevemente introdotta.

\subsection{SVD}
\label{sec:SVD}
Si consideri una matrice nxk $\mathbf{X}$ e sia $t = min(n,k)$. La decomposizione a valori singolari SVD di $\mathbf{X}$ é la fattorizzazione:
\begin{equation}
\label{SVDeq}
	\mathbf{X}=\mathbf{U}\mathbf{D}\mathbf{P}^{T}
\end{equation}
    
dove
%elenco puntato

\begin{enumerate}
\item $\mathbf{U}$ é una matrice $n$ x $t$ ortogonale \footnote{ricordiamo che una matrice U é ortogonale se vale: $U \cdot U^{T}=I$, con I matrice identitá di dimensione $n$ x $n$. Questo porta che, per le matrici ortogonali, $U^{-1}=U^{T}$};
\item $\mathbf{P}$ é una matrice kxt ortogonale;
\item $\mathbf{D}$=diag $\left\{d_{1},...,d_{t}\right\}$é una matrice diagonale con elementi $d_{1}\geq...\geq d_{t}\geq 0$;
\item Gli $r\leq t$ valori non nulli della diagonale della matrice $\mathbf{D}$ son i valori singolari di $\mathbf{X}$;
\item Gli ultimi $t-r$ valori della matrice $\mathbf{D}$ possono esser scartati come le ultime $t-r$ colonne delle matrici $\mathbf{U}$ e $\mathbf{P}$.
\end{enumerate}
Si ha infatti:
\[
\mathbf{X}=\left[ \begin{array}{cc}
\mathbf{U1} &\mathbf{U2}	\end{array} \right] \cdot \left[ \begin{array}{cc}
\mathbf{D1}& 0\\0& 0	
\end{array} \right] \cdot \left[
\begin{array}{c}
\mathbf{P1}^{T} \\ \mathbf{P2}^{T}	
\end{array}\right] \simeq \mathbf{U1}\cdot \mathbf{D1}\cdot \mathbf{P1}
\]
La matrice $\mathbf{P}$ é usualmente chiamata \textbf{matrice dei loadings}, e le sue colonne \textbf{$p_{1},...,p_{t}$} sono i loadings.
La matrice $\mathbf{T}=\mathbf{UD}$ viene chiamata invece \textbf{matrice degli scores}, e le sue rispettive colonne \textbf{$t_{1},...,t_{t}$} scores.


\subsection{Decomposizione autovalori-autovettori}
\label{sec:EIG}
Si consideri una matrice quadrata $\mathbf{A}$ di dimensioni $n$ x $n$, e si consideri la sua decomposizione autovalori-autovettori:
\begin{equation}
  \label{eig}
  \mathbf{A}=\mathbf{P} \mathbf{\Lambda} \mathbf{P}^{T}= \Sigma_{i=1}^{n} \lambda_{i} p_{i} p_{i}^{T}	
\end{equation}

Si ha che $\mathbf{P}$ é matrice ortogonale di colonne \textbf{$p_{1},...,p_{n}$}, e si consideri la matrice degli autovettori $\mathbf{\Lambda}$ in modo che gli autovettori siano ordinati in ordine crescente:
\[
\lambda_{1} \geq ...\geq \lambda_{n}
\]

Si consideri una fattorizzazione di $\mathbf{A}$ come $A=X^{T} X$ e si consideri la decomposizione SVD di $\mathbf{X}$ trattata in precedenza.Risulta:
   
\begin{equation}
  \label{eigDeceq}
	 \mathbf{A} = \mathbf{X}^{T} \mathbf{X} = \mathbf{PDU}^{T}\mathbf{UDP}=\mathbf{PD}^{2}\mathbf{P}
\end{equation}

dove si é sfruttata la proprietá di ortogonalitá della matrice $\mathbf{U}$.
Risulta che i valori sulla diagonale principale della matrice $\mathbf{D}^{2}$, rispettivamente $d_{1}^{2},...,d_{r}^{2}$ sono gli autovalori non nulli della matrice $\mathbf{X}^{T} \mathbf{X}$ di rango $r$. Le colonne di $\mathbf{P}$ sono autovettori di $\mathbf{X}^{T} X$, corrispondenti agli autovalori non nulli.
\footnote{\textbf{N.B.}: se $\mathbf{X}$ rappresenta un set di dati centrati, la matrice $\mathbf{X}^{T} \mathbf{X}/(n-1)$ é la matrice di covarianza di $\mathbf{X}$.}


L'Analisi delle Componenti Principali (PCA) individua, all'interno di un set di variabili, la struttura correlativa esistente tra le variabili stesse ed esplicita la dimensione effettiva dello spazio rappresentato da queste variabili mediante la costruzione delle variabili latenti.
Data la matrice $\mathbf{X}$ dei dati delle variabili di processo, di dimensione n x k, dove $n$ sono i campioni o misure disponibili e $k$ le variabili predittive, il metodo PCA calcola dei vettori, detti componenti principali(PC) o variabili latenti (LV), che esprimono le direzioni di massima variabilità complessiva dei dati disponibili, cioè le direzioni lungo le quali si dispongono la maggior parte delle traiettorie delle variabili. Le variabili sono, per costruzione ortogonali tra loro; la prima variabile latente indica la direzione di massima variabilità delle variabili, la seconda variabile latente indica la seconda direzione di massima variabilità ortogonale alla prima. Per ottenere le componenti principali, il metodo PCA crea due nuovi vettori per ogni PC:
\begin{itemize}
	\item un vettore \emph{loading,} indicato con $p$~;
	\item un vettore \emph{score,}indicato con il simbolo $t.$
\end{itemize}
Per ottenere i vettori \emph{loadings}e gli \emph{scores} si scompone la matrice $X$ dei dati in una serie di $k$ matrici di rango 1, ognuna delle quali rappresenta una componente principale e al cui somma è la matrice $X$ stessa: ognuna di queste matrici è il prodotto vettoriale dei due vettori di \emph{loading} e \emph{score}di ogni PC.
Intuitivamente, si posson riprendere alcune nozioni trattate in precedenza, partendo dalla decomposizione della matrice $\mathbf{X^{T}X}$ come $\mathbf{X^{T}X} = \mathbf{P} \mathbf{\Lambda }\mathbf{P}^{T}$ trattata in \ref{eigDeceq} e la decomposizione SVD citata in \ref{SVDeq}.Ricordando che in \ref{sec:SVD} $\mathbf{P}$ é stata definita \textit{matrice dei loadings} e $\mathbf{T}=\mathbf{UD}=\mathbf{XP}$\footnote{sfruttando la proprietá di ortonormalitá di $\mathbf{P}$  e sostituendo $\mathbf{UD}$ in luogo $\mathbf{T}$ nella decomposizione SVD di $\mathbf{X}$;} \textit{matrice degli scores}, si ha:
	\[ \mathbf{X} = t_{1}p^{T}_{1} + \ldots t_{a}p^{T}_{a} +\ldots+t_{k}p^{T}_{k} = \mathbf{TP}^{T}.\]
Come detto, l'informazione ottenuta dalle variabili di processo è corrotta da rumore, di conseguenza non tutte le $k$ direzioni principali saranno necessarie per rappresentare l'informazione contenuta nei dati. Se vengono mantenute le prime $A$ direzioni principali o variabili latenti, la matrice $\mathbf{X}$ può essere rappresentata come:
	
\begin{equation}
\label{eq2}
	\mathbf{X} = t_{1}p^{T}_{1}+\ldots+t_{a}p^{T}_{a}+t_{r}p^{T}_{r}+E_{A} = \mathbf{TP}^{T}+\mathbf{E}_{A}
\end{equation}
con $A < k$ ed $E_{A}$	matrice dell'informazione residua.
Si presume che l'informazione utile sia contenuta nelle prime $A$ variabili latenti, e il rumore nella matrice residua $E_{A}$.\footnote{In realtà questo raramente accade, la determinazione del numero opportuno di variabili latenti da includere in un modello PLS rimane una questione ancora irrisolta. In assenza di certezze, il metodo più comunemente impiegato rimane quello della \emph{ cross-validation} (maggiori dettagli in \cite{Baibing}) \\}
L'equazione (\ref{eq2}) può essere così interpretata: il processo monitorato attraverso le variabili contenute nella matrice $X$ può essere descritto attraverso $A$ variabili fittizie artificiali, dette variabili latenti, lequali assicurano di rappresentare la maggior parte della variabilità utile del processo e di non essere correlate tra loro grazie all'ortogonalità dei \emph{loadings}.
Le variabili latenti costituiscono le nuove variabili del sistema e sono combinazione lineare delle variabili predittive, colonne di $X$. Le variabili latenti hanno un'importanza diversa nel modello, proporzionale alla quantità di variabilità della matrice $X$ che riescono a spiegare.
Al vettore dei \emph{loadings} e degli \emph{scores} è possibile associare un'interpretazione geometrica:
\begin{itemize}
	\item La matrice dei \emph{loadings} \[ \mathbf{P} = \left[_{1},\ldots,p_{a},\ldots p_{A}\right]\] ha dimensioni A x p ed ogni elemento $p_{aj}$ rappresenta il contributo della j-esima variabile predittiva alla a-esima variabile latente; geometricamente, $p_{aj}$ è il coseno, compreso tra 0 e 1 che quantifica l'angolo  di rotazione della a-esima variabile latente rispetto alla j-esima variabile predittiva nello spazio originario. In pratica ogni ogni vettore $p_{a}$,di dimensioni 1 x p, definisce una nuova direzione nello spazio cartesiano in funzione delle variabili originali.\\
	\item la matrice degli \emph{scores} \[ \mathbf{T} = \left[t_{1},\ldots,t_{a},\ldots,t_{A}\right]\] ha dimensione n x A; gli \emph{scores} sono il risultato della proiezione della matrice $\mathbf{X}$ lungo le direzioni principali individuate dai vettori \emph{loadings}, secondo la relazione:

\begin{equation}
\label{TXP}
	\mathbf{T} = \mathbf{XP}
\end{equation}.
\end{itemize}
Una problematica molto importante che resta da analizzare è la scelta del numero di componenti principali. In letteratura sono presenti molteplici metodologie che permettono di definire il numero di componenti in base a diversi criteri: la scelta è stata attuata confrontando i vari metodi in base alle conoscenze teoriche di cui si disponeva, privilegiando la semplicità sia di comprensione che, guardando più avanti, di implementazione. Il criterio che verrà presentato in seguito rispecchia tali scelte, e pertanto è uno tra i più semplici sia dal punto di vista logico che implementativo.


\subsection{ Scelta del numero di componenti principali}
\label{sec:SelezPC}
Come presentato fin qui, PCA permette di selezionare le componenti principali di una matrice $\mathbf{X}$ e di decomporla come:
\[
  \mathbf{X}= \mathbf{T}_g \mathbf{P}_{g}^{T} + \mathbf{X}_{g+1}
\]
 dove $\mathbf{T}_g \mathbf{P}_{g}^{T}$ rappresenta la parte rilevante della matrice, mentre $\mathbf{X}_{g+1}$ è la parte trascurabile dei residui dati dagli eventuali $r-g$ termini rimanenti\footnote{trascurabili , si intende, rispettivamente al problema che si sta affrontando. Si ricorda che $r$ è pari al numero totale delle componenti principali di $\mathbf{X}$}.
La problematica che si pone a questo punto, alquanto rilevante per i successivi trattamenti della matrice, è trovare una metodologia di scelta di $g$, ovvero quante componenti principali si devono selezionare in $\mathbf{X}$ per il problema che si sta trattando. Esiste un ampia scala metodologie di scelta di $g$, che variano a seconda della difficoltà di implementazione e del grado più o meno buono dell'ottimalità della soluzione che forniscono (vedi PRESS in \cite{Baibing}).
Un semplice metodo di scelta si basa sugli autovalori di $\mathbf{X}$ \footnote{Per le seguenti affermazioni, si ipotizza che gli autovalori siano stati ordinati in maniera decrescente.In caso di SVD, tale ordinamento viene automaticamente eseguito, in caso di utilizzo della scomposiczione autovalori-autovettori, è necessario invece ordinarli.}. Gli autovalori  vengono normalizzati in valore e ad ognuno viene assegnato un peso percentuale della somma di tutti gli autovalori:
\[
  peso(\lambda_{i})= \frac{\lambda_{i}}{\lambda_{1}+...+\lambda_{r}} \times 100
\]
e tale percentuale viene interpretata come la \textit{variazione percentuale spiegata} dalla corrispondente componente principale. Tali pesi vengono quindi usati in modo cumulativo, così che il peso dato dalle prime $g$ componenti è dato da:
\[
  \frac{\lambda_{1}+\lambda_{2}+...+\lambda_{g}}{\lambda_{1}+...+\lambda_{r}} \times 100
\]
Viene quindi fissata una regola di scelta di $g$ come l'insieme di componenti che \emph{spiegano} il \emph{tot\%} della matrice $\mathbf{X}$:in questo modo la scelta del numero di componenti da selezionare può esser automatizzata al variare della percentuale data. Un buon valore di percentuale da scegliere si assesta tra l'$80$ e il $90\%$.\\
Una giustificazione a questa scelta può esser facilmente spiegata: se si considera la varianza del j-esimo score $\mathbf{t}_{j}$ come
\[
  s^{2}(\mathbf{t}_{j})= \frac{1}{n-1} \left\| \mathbf{t}_{j} \right\|^{2}\\
  = \frac{1}{n-1} \mathbf{t}_{j}^{T} \cdot \mathbf{t}_{j}\\
  = \frac{1}{n-1} \lambda_j
\]
si nota come $\lambda_j$ sia proporzionale alla varianza del corrispondente score. In particolare, tutte le componenti che hanno il rispettivo $\lambda_j$ nullo non danno alcun peso e posson esser quindi scartate. Considerando quindi la varianza matriciale di $\mathbf{X}$
\[
  \mathbf{V}_{x} = \frac{1}{n-1} \mathbf{X}^{T} \mathbf{X}\\
  = \frac{1}{n-1} \sum^{j=1}_{k} \lambda_{j} \mathbf{p}_{j}^{T} \mathbf{p}_{j}
\]
si può considerare $\lambda_{j}$ come contributo (pesato) di $\mathbf{t}_{j}$ alla covarianza totale di  $\mathbf{X}$ \footnote{Si sottolinea come la somma degli autovalori sia pari alla varianza scalare di $\mathbf{X}$:
\[
tr\left\{\mathbf{X}^{T} \mathbf{X}\right\}= (n-1) \left\{ s^{2}(\mathbf{x}_{1})+...+s^{2}(\mathbf{x}_{k}) \right\}
\] 
}.
\\


\section{Regressione delle Componenti Principali}
\label{sec:PCR}
\textbf{PCR} é l'acronimo di \textit{Principal Component Regression}, ovvero regressione delle componenti principali. Con tale metodo si intende prelevare le componenti principali da una matrice di dati e regredire tali componenti mediante metodi lineari o minimi quadrati. A questo scopo si considerano le componenti principali di una matrice $\mathbf{X^{T}X}$, dove $\mathbf{X}$ é una matrice $n$ x $k$ di dati centrati. Ci sono fondamentalmente due metodologie per estrarre le componenti principali dalla matrice $\mathbf{X^{T}X}$, perfettamente equivalenti : SVD e decomposizione autovalori-autovettori (\textit{eigenvalue decomposition}). Applicando il primo dei due metodi, si ottiene la forma ridotta:
\[
X= U_{r} D_{r} P_{r}^{T}
\]
dove con $r$ si voglion indicare le matrici risultanti da SVD tolte le colonne/righe coorrispondenti agli ultimi $t-r$ autovalori nulli. Si definiscono gli \textit{scores} come:
\[
T_{r} = U_{r} D_{r}
\]
una matrice con le colonne ortogonali, ma non necessariamente ortonormali
\footnote{Infatti, \[ 
T^{T} T =D U^{T} U D =D^{2}=\Lambda
 \]}.
Da $X=T P^{T}$ si ha che:
\[
X^{T} X = P_{r} T_{r}^{T} T_{r} P_{r}^{T} = P_{r} \Lambda_r P_{r}^{T}
\]
che é la decomposizione autovalori-autovettori di $X^{T} X$, in cui son state eliminate le colonne di $P$ relative agli autovalori nulli di $\Lambda$.

Si può perciò pensare che la matrice $\mathbf{X}$ sia scomposta come:\\
\[  
\mathbf{X}= \mathbf{T}_g \mathbf{P}_{g}^{T} 
\]
e le caratteristiche importanti di $\mathbf{X}$ sian racchiuse in $\mathbf{T}_g$, matrice degli scores delle prime $g$ componenti principali. Si ipotizza che $\mathbf{X}$ sia la matrice contenente $m$ ingressi campionati in $n$ istanti differenti, e $\mathbf{Y}$ sia la relativa matrice delle $p$ uscite di un sistema lineare del tipo:
\[
\mathbf{Y}=\mathbf{X} \cdot \mathbf{C} + \mathbf{F}
\]
con $\mathbf{F}$ rumore del sistema.

Si puó a questo punto compiere PCR, ovvero \textit{Multiple Linear Regression} utilizzando peró non l'intera matrice $\mathbf{X}$, ma $\mathbf{T}_g$ in luogo di questa. Si ottiene dunque:
 
\begin{equation}
\mathbf{Y}=\mathbf{T}_g \cdot \mathbf{C} + \mathbf{F_g}	
\end{equation}
 
e la relativa soluzione ai minimi quadrati é data da:
\begin{equation}
	\mathbf{\hat{C}} =\left( \mathbf{T}_{g}^{T} \mathbf{T}_{g} \right)^{-1} \mathbf{T}_{g}^{T} \mathbf{Y}
\end{equation}


Nel caso che si sta trattando, tale approccio va esteso, perché non adatto al problema: si vuole infatti trovare un sistema \textbf{dimamico} e non un legame \textbf{statico} tra ingressi e uscite del sistema. L'idea é dunque quella di compiere \textit{PCA} sui dati d'ingresso e passare quindi ad una routine di identificazione a sottospazi l'insieme delle uscite e gli scores degli ingressi in luogo di questi ultimi. In questo modo, infatti, le problematiche trattate in \ref{sec:OLS}, comuni anche per tali tecniche ad identificazione, vengono a mancare permettendo così di ottenere un sistema più robusto e flessibile, sia dal punto di vista predittivo che in catena aperta. In realtà, si desidera un sistema in grado di accettare $m$ ingressi, pari all'insieme totale di ingressi di cui si dispone inizialmente, e non l'insieme $g$ ridotto.\\ Il sistema identificato si presenta nella forma:
\begin{equation}
\label{ss}
	\left\{
	\begin{array}{rcl}
	\mathbf{x}(t+1)= A \mathbf{x}(t)+ B \mathbf{u}(t)+Ke(t);\\
	\mathbf{y}(t)=C \mathbf{x}(t) + D \mathbf{u}(t)+e(t);
	\end{array}\right.
\end{equation}
con $A\in R^{n*n}$,$B\in \mathbb{R}^{n*g}$,$C\in \mathbb{R}^{p*n}$ e $D\in \mathbb{R}^{p*g}$, sistema che accetta $g$ ingressi. Dati i loading $\mathbf{P}_g$ ricavati tramite il metodo \textit{PCA} descritto in \ref{sec:PCA} si compie quindi di volta in volta un cambio di base sugli ingressi $u(t) \in \mathbb{R}^{m}$:
\[
  \hat{u}(t)=\mathbf{P}^{T}_{g} u(t)
\]
e $\hat{u}(t) \in \mathbb{R}^{g}$.



\section{Algoritmo NIPALS per PCA} 
NIPALS é l'acronimo di \textit{Nonlinear Iterative PArtial Least Square}, e si tratta di un algoritmo che implementa l'analisi delle componenti principali. Si vogliono estrarre le prime $g$ componenti principali di $X^{T} X$, iniziando con l'autovalore piú grande ($\lambda_1$) e analizzando via via gli autovalori in maniera decrescente\footnote{La metodologia con la quale si definisce il numero $g$ di componenti principali da selezionare é riportata in (\ref{sec:SelezPC}). Naturalmente dev'essere $g\leq r$, con $r$ pari al numero massimo delle componenti principali.}.La comprensione di tale algoritmo é necessaria perché verrá utilizzato ampiamente nelle seguenti sezioni:NIPALS permette di ottenere i vettori degli \emph{scores} e dei \emph{loadings} così come indicato in \cite{Wold},\cite{Tonizzo}.L'algoritmo viene inizializzato con $\mathbf{E}_{j}=\mathbf{X}$ e $j=0$.Successivamente vengono compiuti i seguenti passi:
\begin{enumerate}
	\item si prende un vettore colonna, cioè una variabile $\mathbf{e}_{j}$ di $\mathbf{E}_j$ e la si denomina $\mathbf{t}_{j}$: $\mathbf{e}_{j}=\mathbf{t}_{j}$ 
	\item si calcola il vettore \emph{loading} corrispondente: \[ \mathbf{p}_{j} =\mathbf{E}_{j} \mathbf{t}_{j}/\mathbf{t}^{T}_{j}\mathbf{t}_{j} \]
questo equivale a risolvere una regressione ai minimi quadrati delle colonne di $\mathbf{X}$ su $\mathbf{t}_{j}$ senza problemi di correlazione, perchè $\mathbf{t}^{T}_{j}\mathbf{t}_{j}$ è uno scalare; in questo modo si cattura anche la direzione di massima variabilità della matrice $\mathbf{X}$.
	\item si normalizza $\mathbf{p}_{j}$ a norma unitaria: $ \mathbf{p}_{j}= \mathbf{p}_{j}/\left\|\mathbf{p}_{j}\right\|$ ;
	\item si ricalcola il vettore \emph{score}: \[ \mathbf{t}_{j} = \mathbf{E}_{j}\mathbf{p}_{j}/\mathbf{p}^{T}_{j}\mathbf{p}_{j} \] anche in questo caso si ha una regressione ai minimi quadrati senza problemi di inversione, perchè $p^{T}_{j}p_{j}$; è uno scalare 
	\item si confronta $\mathbf{t}_{j}$ nel punto 2 con $\mathbf{t}_{j}$ al punto 4: se la differenza relativa è nulla (con uno scarto di $10^{-8}-10^{-6}$),si passa al punto 6, altrimenti si ritorna a 2;
	\item si calcola il residuo (della matrice $X$): $\mathbf{E}_{j+1}=\mathbf{E}_{j}-\mathbf{t}_{j}\mathbf{p}^{T}_{j}$ 
\end{enumerate}

Se $g=r$ allora l'algoritmo trova tutte le componenti principali. A questo punto le matrici $\mathbf{T}$ e $\mathbf{P}$ sono conformi alla relazione:
\begin{equation}
	\mathbf{X} = \mathbf{T} \mathbf{P}^{T}
\end{equation}

\subsection{ \textit{Proprietà}}
Per comprendere meglio come agisce l'algoritmo, vengono analizzati i passi in dettaglio, utilizzando le proprietà sopracitate per le matrici di \textit{loadings} e \textit{scores}.
Si ponga $\lambda_{j}=\left\|\mathbf{E}_j \mathbf{t}_j \right\|$ e si scriva il passo 2 dell'algoritmo come:
\[    \mathbf{E}_j \mathbf{t}_j =\lambda_{j} \mathbf{p}_j  \]
Dalla relazione del passo 4 segue:
\[  \mathbf{E}_j^{T} \mathbf{E}_j \mathbf{p}_j =\lambda_{j} \mathbf{p}_j \]
equazione che viene soddisfatta iterando dal passo 2 al passo 5 fino a che questi convergano ad una soluzione in norma d'errore minima. Da qui si nota facilmente come $\lambda_{j}$ e $\mathbf{p}_{j}$ siano rispettivamente gli autovalori e gli autovettori della matrice $\mathbf{E}_j^{T} \mathbf{E}_j$.\\
Si dimostra ora come, nel passo 6, si formi una matrice $\mathbf{E}_{j+1}=\mathbf{E}_{j} - \mathbf{t}_{j} \mathbf{p}_{j}^{T}$ che è ortogonale a $\mathbf{t}_{j}$, j-esimo score. Si ha infatti:
\[     
\begin{array}{ccc}
	\left(\mathbf{E}_{j} - \mathbf{t}_{j} \mathbf{p}_{j}^{T}\right)^{T} \cdot \mathbf{t}_{j} &=&\mathbf{E}_{j}^{T} \mathbf{t}_{j} - \mathbf{p}_{j} \mathbf{t}_{j}^{T} \mathbf{t}_{j} \\
	&=& \mathbf{E}_{j}^{T} \mathbf{E}_{j} \mathbf{p}_{j} - \mathbf{p}_{j} \lambda_j \\
	  &=&  \lambda_{j} \mathbf{p}_{j} - \mathbf{p}_{j} \lambda_j =\mathbf{0}
\end{array}
\]
Si può dunque facilmente comprendere come, prendendo $\mathbf{t}_{j+1}$, una colonna a caso di $\mathbf{E}_{j+1}$ nel passo 1 dell'algoritmo, questa sia ortogonale, qualunque essa sia, ai precedenti scores $\mathbf{t}_{j}$ selezionati.
Dopo $g$ cicli si avrà dunque:
\[ \mathbf{E}= \mathbf{t}_{1} \mathbf{p}_{1} + \mathbf{t}_{2} \mathbf{p}_{2} + ...+ \mathbf{t}_{g} \mathbf{p}_{g} +     
   \mathbf{E}_{g+1}
\]
con $\mathbf{E}_{g+1}=\mathbf{0}$ se son state selezionate tutte le componenti, ovvero $g=r$.
\\
\\

\section{Proiezioni su strutture latenti: PLS}
La proiezione su strutture latenti (PLS) è un metodo di regressione che costruisce un'equazione di regressione lineare tra gli \emph{scores} delle variabili predittive e gli \emph{scores} delle variabili dipendenti; analogamente al metodo PCA, il metodo PLS decompone sia la matrice delle variabili predittive $\mathbf{X}$ che la matrice delle variabili dipendenti $\mathbf{Y}$, creando direzioni principali che descrivano la massima variabilità di $\mathbf{X}$ cogliendo allo stesso tempo la massima covarianza tra gli \emph{scores} di $\mathbf{X}$ e di $\mathbf{Y}$.
Le due matrici $\mathbf{X}$ e $\mathbf{Y}$ possono essere decomposte con due PCA secondo le relazioni esterne:
	\[ \mathbf{X} = \mathbf{T}\mathbf{P}^{T} + \mathbf{E},\]
	\[ \mathbf{Y} = \mathbf{U}\mathbf{Q}^{T} + \mathbf{F}.\]
Se si utilizzano gli \emph{scores} di $\mathbf{X}$ come variabili predittive e gli \emph{scores} di $\mathbf{Y}$ come nuove variabili dipendenti, è possibile ottenere una equazione di regressione interna tra gli \emph{scores} cioè si pone 
	\[ \mathbf{U} = \mathbf{T}\mathbf{b_{pls}} + \mathbf{F}^{*}\]
dove il vettore $\mathbf{b_{pls}} = \left[b_{1},\ldots,b_{a},\ldots,b_{A} \right]^{T}$ è il vettore dei coefficienti di regressione tra le variabili latenti, $\mathbf{F}^{*}$ è la matrice degli scarti e $\mathbf{b_{pls}}$ è stato ottenuto con la regressione ai minimi quadrati tra gli \emph{scores} di $\mathbf{X}$ e di $\mathbf{Y}$: 
\[\mathbf{b}_{a,pls} = \mathbf{u}^{T}_{a}\mathbf{t}_{a}/\mathbf{t}^{T}_{a}\mathbf{t}_{a}.\] e 
\[ \hat{\mathbf{u}} = \mathbf{b}_{a,PLS}\mathbf{t}_{a}. \]
Il modello così calcolato non è però il migliore possibile, perchè non sempre gli \emph{scores} di $\mathbf{X}$ e $\mathbf{Y}$ sono correlati tra loro; si può quindi ottenere un'equazione di regressione che rappresenta solo una minima parte della variabilità di $\mathbf{Y}$ o che non minimizza $\mathbf{F}^{*}$.
Si è quindi cercato di decomporre simultaneamente le matrici $\mathbf{X}$ e $\mathbf{Y}$, ruotando parzialmente gli \emph{scores} $\mathbf{t}_{a}$ in modo da renderli più correlati agli \emph{scores} $\mathbf{u}_{a}$. Per fare questo si sono introdotti anche i vettori \emph{weights}, simili nel significato ai vettori \emph{loadings}.
I \emph{loadings} e gli \emph{scores} sono calcolati con l'algoritmo NIPALS, adattato alla regressione PLS:
\begin{enumerate}
	\item prendere un vettore colonna o meglio una variabile $\mathbf{y}_{j}$ di $\mathbf{Y}$ e denominarlo $\mathbf{u}_{a}:\mathbf{y}_{j} = \mathbf{u}_{a}$~;
	\item calcolare il vettore \emph{weight} $\mathbf{X}$ corrispondente: 
	\[\mathbf{w}_{a}=\mathbf{X}\mathbf{u}_{a}/\mathbf{u}^{T}_{a}\mathbf{u}_{a}\]
	 regredendo $\mathbf{X}$ su $\mathbf{u}_{a}$ trovo la relazione lineare che permette di spiegare una variabile dipendente, $\mathbf{u}_{a}$ utilizzando l'informazione contenuta nella matrice delle variabili predittive. L'utilizzo del modello ai minimi quadrati, inoltre, garantisce che ciò è stato fatto prendendo la direzione di massima variabilità dei dati di $\mathbf{X}$ e minimizzando la matrice residua $\mathbf{F}^{*}$~.
	\item normalizzare $\mathbf{w}_{a} = \mathbf{w}_{a}/\left\|\mathbf{w}_{a}\right\|$;
	\item calcolare il vettore \emph{score} di $\mathbf{X}: \mathbf{t}_{a} = \mathbf{X}\mathbf{w}_{a}/\mathbf{w}^{T}_{a}\mathbf{w}_{a}$~;
	\item calcolare il vettore \emph{loading} di $\mathbf{Y}: \mathbf{q}_{a} = \mathbf{Y}\mathbf{t}_{a}/\mathbf{t}^{T}_{a}\mathbf{t}_{a}$~;
  \item normalizzare $\mathbf{q}_{a}$ a norma unitaria: $\mathbf{q}_{a} = \mathbf{q}_{a}/\left\|\mathbf{q}_{a}\right\|$~;
  \item ricalcolare il vettore \emph{score} di $\mathbf{Y}: \mathbf{u}_{a} = \mathbf{Y}\mathbf{q}_{a}/\mathbf{q}^{T}_{a}\mathbf{q}_{a}$~;
  \item confrontare il valore di $\mathbf{t}_{a}$ al punto 4 tra un'iterazione e l'altra: se la differenza relativa è nulla (con uno scarto di $10^{-8} \div 10^{-6}$) passare al punto 9, altrimenti tornare al punto 2;
  \item calcolare il vettore \emph{loading} di 
	\[ \mathbf{X}: \mathbf{p}_{a} = \mathbf{X}\mathbf{t}_{a}/\mathbf{t}^{T}_{a}\mathbf{t}_{a}; \]
  \item ortogonalizzare il \emph{weight}: 
  \[\mathbf{w}_{a} = \mathbf{w}_{a}\mathbf{p}_{a}/\mathbf{p}^{T}_{a}\mathbf{p}_{a};\]
  \item calcolare il coefficiente di regressione: 
	\[ \mathbf{b}_{a,pls} = \mathbf{u}^{T}_{a}\mathbf{t}_{a}/\mathbf{t}^{T}_{a}\mathbf{t}_{a};\]
  la relazione interna è data da:
  \[ \hat{\mathbf{u}_{a}} = \mathbf{b}_{a,pls}\mathbf{t}_{a};\]
 ma la forma degli \emph{scores} e dei \emph{loadings} è cambiata;
  \item aggiornare le variabili: 
  \[ \mathbf{E}_{a} = \mathbf{E}_{a-1} - \mathbf{t}_{a}\mathbf{p}^{T}_{a};\]
  \[ \mathbf{F}_{a} = \mathbf{F}_{a-1} - \mathbf{u}_{a}\mathbf{q}^{T}_{a};\]
  \[ \mathbf{F}^{*}_{a} = \mathbf{F}^{*}_{a-1} - \hat{\mathbf{u}_{a}}q^{T}_{a};\]
\end{enumerate}


In pratica si combinano gli algoritmi di decomposizione PCA per le due matrici $\mathbf{X}$ e $\mathbf{Y}$ tra loro, scambiando informazioni tra gli \emph{scores} e i \emph{loadings} delle variabili predittive e dipendenti, di modo che sia $\mathbf{X}$ che $\mathbf{Y}$ siano modellate dalle stesse variabili latenti.
La matrice $\mathbf{F}^{*}$~minimizza grazie alla regressione ai minimi quadrati, la differenza tra errori calcolati e sperimentali: $\mathbf{F}^{*} = \mathbf{Y} - \mathbf{T}\mathbf{Q}^{T}$~.\\
A differenza dell'algoritmo PCA gli \emph{scores} $\mathbf{t}$ sono combinazioni lineari delle variabili predittive di $\mathbf{X}$ pesate secondo i \emph{weights} $\mathbf{w}^{*}$, e non secondo i \emph{loadings} $\mathbf{p}$. 
\[ \mathbf{T} = \mathbf{X}\mathbf{W}^{*}; \]
La differenza tra \emph{weights} e \emph{loadings} della matrice $\mathbf{X}$ è che: i primi forniscono i pesi di ogni variabile predittiva a formare la variabile latente, i secondi trovano la direzione della variabile latente e spiegano la varianza delle singole matrici; a differenza del PCA, le due caretteristiche, contributo alla variabile latente e formazione della direzione principale, non coincidono.\\
Si nota inoltre che esiste una differenza tra i \emph{weights} $\mathbf{w}$ e $\mathbf{w}^{*}$ legati dalla relazione: $\mathbf{W}^{*} = \mathbf{W}(\mathbf{P}^{T}\mathbf{W})^{-1}$.\\
Un modello PLS opportunamente costruito contiene variabili latenti che descrivono la varianza di $\mathbf{X}$ e $\mathbf{Y}$ con valori decrescenti e cumulativamente simili; se ad esempio, un modello spiega la maggior parte della varianza di $\mathbf{Y}$ ma poca varianza di $\mathbf{X}$~, significa che la matrice $\mathbf{X}$ contiene variabili poco correlate ad $\mathbf{Y}$, che eventualmente danneggiano la modellizzazione.\\
La regressione condotta con metodo PLS è rigorosamente lineare: per questo motivo, spesso vengono suggeriti metodi per eliminare il carattere dinamico delle traiettorie, come il centramento, la normalizzazione, la trasformazione di variabili.

\subsection{Interpretazione geometrica di PLS}
PLS è un metodo di proiezione che geometricamente si può interpretare come la proiezione della matrice $\mathbf{X}$ (Un insieme di punti in uno spazio a K dimensioni) su un iperpiano A-dimen\-sionale; in questo modo le coordinate dei proiettori (\emph{scores}) $t_{a}: \quad a = 1,\ldots ,A$ sono dei buoni predittori per $\mathbf{Y}$, figura \ref{fig:iperpiano}. 
La direzione del piano è espressa come inclinazione, $\mathbf{p_{ak}}$, di ogni componente rispetto alle coordinate, $\mathbf{x_{k}}.$ La pendenza è il coseno dell'angolo tra la direzione del piano di proiezione e l'asse delle ascisse.\\
PLS costruisce iterativamente un iperpiano nello spazio generato dai dati di $\mathbf{X}$; in questo modo tale iperpiano è la migliore approssimazione di $\mathbf{X}$, e allo stesso tempo, la posizione di dati proiettati sul piano, descritta dal valore degli \emph{scores} $\mathbf{t_{a}}$ sono utilizzati per la predizione della risposta del sistema, $\mathbf{Y}.$
 
\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.4in,height=3in]{./figurePLS/iperpiano.eps}
	\caption{Iperpiano generato dall'algoritmo PLS nello spazio dei punti di $\mathbf{X}$, si evidenzia l'ortogonalità delle componenti}
	\label{fig:iperpiano}
\end{figure}

\subsubsection{Confronto geometrico tra OLS e PLS}
\label{sec:ConfrontoGeometricoTraOLSEPLS}
La robustezza dell'algoritmo di regressione è riferita all'insensibilità del modello stimato al condizionamento errato dei dati e al rumore. Il grande vantaggio delle tecniche di proiezione su strutture latenti rispetto alla regressione lineare ordinaria sono illustrati geometricamente in figura \ref{fig:geometriaOLS} e in figura \ref{fig:geometriaPLS}.\\
Tali figure evidenziano il caso di due ingressi e una uscita in presenza di collinearità dei dati e di rumore sovrapposto. Tutti i dati sono esattamente sulla stessa retta tranne un punto $x,$ il quale è deviato per la presenza di rumore. Questi dati generano uno sottospazio a due dimensioni chiamiamolo, $\mathbf{X}$.\\
La regressione lineare ordinaria, proietta l'uscita $\mathbf{Y}$ ortogonalmente a $\mathbf{X}$. Purtroppo dato che il dato $x$ è corrotto da rumore l'orientazione del piano $\mathbf{X}$ è pesantemente influenzata dalla collocazione del punto  $x.$ Il risultato è che la proiezione  $\hat{Y}_{OLS}$  è altamente sensibile alla collocazione di $x$ e quindi al rumore casuale.\\
Nel caso invece di algoritmo PLS si nota come sia richiesto un  fattore, cioè una  proiezione artogonale sul sottospazio uno-dimensionale $\mathbf{t_{1}}$ di $\mathbf{X}$. In  questo caso la proiezione $\hat{Y}_{PLS}$ non è influenzata dalla collocazione di $x,$ cioè è robusta al rumore. (maggiori dettagli in \cite{Qin1}).

\begin{figure}[htb!]
	\centering
		 \includegraphics[width=3in,height=1.8in]{./figurePLS/geometriaOLS.eps}
		\caption{Interpretazione geometrica di OLS. Un dato corrotto da rumore condiziona negativamente la regressione lineare}
	  \label{fig:geometriaOLS}
\end{figure}


\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.1in,height=1.8in]{./figurePLS/geometriaPLS.eps}
	\caption{Interpretazione geometrica di PLS. Il metodo PLS è robusto a collinearità dei dati e agli effetti del rumore sovrapposto}
	\label{fig:geometriaPLS}
\end{figure}

\section{Selezione dei sensori con PLS: caso statico}
Ciò che il nostro lavoro si propone, oltre al fatto di identificare un modello termodinamico è anche quello di ridurre il numero di sensori utilizzati mantenendo buona la qualità della stima della temperatura interna. 
La regressione PLS è negativamente influenzata dalla presenza di variabili predittive non correlate con la variabile dipendente. La presenza di variabili "inutili", rende necessario un maggior numero di variabili latenti per modellare i dati. Anche se questo metodo di regressione non è influenzato dalla presenza di variabili tra loro correlate, diventa necessario operare una scelta sulle variabili indipendenti, in modo da selezionare quelle maggiormente esplicative della variabile dipendente e aventi un elevato rapporto segnale/rumore, in modo da rendere le stime di qualità del prodotto più robuste e il modello meno influenzato dal rischio di avarie dei sensori.
La scelta delle variabili è contrastata tra due esigenze:
\begin{enumerate}
	\item la mecessità di rappresentare nel modo migliore possibile i dati sperimentali, minimizzando il MSE, porterebbe a introdurre nel modello il maggior numero possibilie di variabili predittive. Con questo criterio vengono utilizzate variabili tra loro correlate nel modello e può condurre a \emph{overfitting}, perchè introduce eccessive variabili nel modello che spiegano solo la parte della varianza dovuta al rumore e impedisce l'estrapolazione dei risultati dell'equazione di regressione al di fuori dell'intervallo in cui sono stati calibrati i dati;
	\item la necassità di ottenere un modello il più possibile semplice e che minimizzi la varianza media del dato calcolato $var(\hat{y})$ richiederebbe di inserire ill minor numero possibile di variabili predittive.
\end{enumerate}
Oltre a queste esigenze, di tipo statistico, si desidera che le variabili predittive prescelte siano ottenibili con apparecchiature poco costose e soggette a piccoli errori di misura e rumore.
Esistono inn letteratura diverse soluzioni per risolvere problemi diversi:
\begin{itemize}
	\item Quale criterio si deve adottare per determinare la superiorità di un'equazione di regressione rispetto ad un'altra e il sottoinsieme "migliore" di variabili predittive che determina il modello. Il criterio migliore rimane la minimizzazione dell'errore di predizione MSEP. Questo criterio, però pone il problema della distinzione tra modellizzazione dell'informazione e del rumore: un'equazione può riprodurre ottimamente i dati reali, ma essere inutile perchè dii fatto riproduce ile oscillazioni casuali introdotte dalla strumentazione di misura nei dati raccolti;
	\item Quale algoritmo è capace di selezionare il set più appropriato di variabili in modo da risulatare superiore agli altri rispetto al criterio di giudizio scelto, sapendo che la selezione dell'insieme ottimo dipende anche dallo stimatore che si intende costruire.
\end{itemize}
\subsection{Selezione di variabili per PLS: VIP}
In letteratura sono stati svilippati diversi metodi di eliminazione di variabili implementati sull'algoritmo NIPALS, tra cui: OSC(\emph{Orthogonal Scatter Correction}) sviluppato in \cite{Ohman}; O-PLS (\emph{orthogonal-PLS}) sviluppato in \cite{Trygg}. \\
Tra quelli proposti il più utilizzato è il VIP, acronimio di \emph{Variable Importance in Projection,} trattato in \cite{Chong} il quale calcola un parametro di importanza per ogni variabile predittiva in base al valore pesato di ogni \emph{loadings} nelle variabili latenti secondo:
	\[ VIP_{j}=\sqrt{\frac{p\sum^{A}_{a=1}\left(b^{2}_{a}t^{T}_{a}t_{a}\left(w_{aj}/\left\|w_{aj}\right\|\right)^{2}\right)}{\sum^{A}_{a=1}\left(b^{2}_{a}t^{T}_{a}t_{a}\right)}} \]
Tutte le variabili aventi il criterio $VIP_{j} > 1$ sono ammesse al set ridotto.
VIP effettua sempre le stesse scelte indipendentemente dal rapporto segnale/rumore delle variabili, mentre varia se cambiano il numero di variabili latenti utilizzate dal modello.
Il limite intrinseco di questi metodi è quello di costiture una relazione \emph{statica} di ricerca delle variabili "ottime" per costruire un modello ridotto, noi abbiamo già lavorato in questa direzione ottenendo buoni risultati con l'analisi delle componenti principali e la relativa regressione, cercheremo invece di costruire un metodo dinamico che sfrutti le potenzialità di PLS per identificare un modello e scegliere dinamicamente un set opportuno di sensori pre creare un modello ridotto adatto al controllo e che dia buoni risultati inn fase di predizione.
 


	
\section{PLS: caso dinamico, modello ARX}
Si esamina in questa sessione l'utilizzo e l'implementazione di un algoritmo di PLS dinamico. I medoti di selezione delle variabili fin qui analizzati come PCA  permettono di ottenre delle informazioni sul numero di ingressi e uscite utilizzati in fase di identificazione attraverso un  procedimento di carattere statico. 
Si cercherà qui invece di costruire una relazione dinamica che permetta di ridurre la dimensione del problema mantenendo tuttavia la qualità del modello identificato. Il modello ottenuto con questo tipo di tecnica sarà un modello ARX diversamente da quanto fatto con metodi PCA.\\
Per testarne l'acuratezza si utilizzeranno metodi di validazione comuni. Descriviamo l'idea proposta da Qin J. in \cite{Qin} per costruire un modello ARX utilizzando dinamicamente l'algoritmo PLS e le sue potenzialità.
Nei processi di identificazione caratterizzati da una grande mole di dati, l'algoritmo PLS viene utilizzato per creare un modello dei che spieghi bene i dati. Questo tipo di approccio ha delle forti limitazioni:
\begin{itemize}
	\item  il modello che ne risulta è complicato da aggiornare e difficilmente riesce ad inseguire i cambiamenti dinamici del modello,
	\item  quando si devono manipolare grosse quantità di dati, con sistemi MIMO, come nel nostro caso la disponibilità di memoria fisica di un  calcolatore diviene un fattore limitante.
\end{itemize}
Si cerca quindi di modificare l'approccio statico dell'algoritmo PLS a vantaggio di un algoritmo ricorsivo per un'identificazione dinamica del modello. L'idea è di fondere assieme due concetti, quello della regressione su strutture latenti di dati altamente correlati e affetti da rumore, e la possibiltà di operare con un modello ARX che tiene conto di $n$ ritardi temporali antecedente all'istante $k-esimo$, in modo da avere una struttura dinamica. 
Diversi benefici vengono raggiunti utilizzando un algoritmo ricorsivo di PLS:
\begin{enumerate}
	\item Il PLS ricorsivo aggiorna il modello quando un nuovo campione è disponibile,
	\item Il PLS ricorsivo ha l'abilità di adattare la struttura correlativa del modello a seconda della dinamica del processo,
	\item Il PLS ricorsivo tende a comprimere l'informazione e contenuta nei dati e risolvere il problema della disponibilità di memoria.
\end{enumerate}
Un semplice modello stocastico auto-regressivo con ingressi esogeni è dato da:
\begin{equation}
	\mathbf{y(k)} = \sum^{n_{y}}_{i=1}\mathbf{A_{i}}\mathbf{y}(k-i) + \sum^{n_{u}}_{j=1}\mathbf{B_{j}}\mathbf{u}(k-j) + \mathbf{v}(k);
\end{equation}
dove $\mathbf{y}(k)$, $\mathbf{u}(k)$, e $\mathbf{v}(k)$ sono il processo di uscita, di ingresso e il vettore di rumore con appropriate dimensioni per un sistema multi-input, multi-output. \\
$A_{i}$ e $B_{j}$ sono matrici dei soefficienti del modello che devono essere identificati. $n_{y}$, e $n_{u}$, sono i ritardi, o meglio quanti campioni precedenti all'istante $k-esimo$ vengono considerati per l'uscita e l'ingresso.\\
Per poter elaborare i dati con l'algoritmo PLS è necesasario definire il vettore:

$\mathbf{x}^{T}(k)=\left[\mathbf{y}^{T}(k-1),\mathbf{y}^{T}(k-2),\ldots,\mathbf{y}^{T}(k-n_{y}),\right.$
\begin{equation}
	 \left.\mathbf{u}^{T}(k-1), \mathbf{u}^{T}(k-2),\ldots,\mathbf{u}^{T}(k-n_{u})\right]
\end{equation}

dove la dimensione è $m,$ cioè la somma di tutte le uscite e gli ingressi del sistema per il rispettivo ordine di ritardo $n_{y}$ e $n_{u}.$
Noi nella fase sperimentale cercheremo di trovare l'ordine ottimo del modello, ponendoci in questa prima fase nel caso in cui $n_{y} = n_{u} = n_{uy}.$\\
Una volta definita la matrice dei dati in questo modo si si costruiscono le matrici per n campioni:
\[ \mathbf{X} = \left[ \mathbf{x}(1),\mathbf{x}(2),\ldots, \mathbf{x}(n) \right]^{T} \quad \in \mathbf{R}^{nxm}\]
\[ \mathbf{Y} = \left[ \mathbf{y}(1), \mathbf{y}(2),\ldots, \mathbf{y}(n)\right]^{T} \quad \in \mathbf{R}^{nxp}\]
dove $p$ è la dimensione del vettore delle uscite  $\mathbf{y}(k)$.\\
Definiamo la matrice dei valori da identificare del modello ARX come, 
\[ \mathbf{C} = \left[ \mathbf{A}_{1},\mathbf{A}_{2},\ldots,\mathbf{A}_{n_{y}},\mathbf{B}_{1} ,\mathbf{B}_{2},\ldots,\mathbf{B}_{n_{u}}\right]^{T}\quad \in \mathbf{R}^{mxp},\] 
le due matrici dei dati $\mathbf{Y}$, $\mathbf{X}$ sono in relazione secondo:
\[ \mathbf{Y} = \mathbf{X}\mathbf{C} + \mathbf{V}.\]
Ciò che l'algoritmo PLS restituisce e che ci serviranno per ricostruire la matrice $\mathbf{C}$ dei coefficienti del modello ARX, sono le matrici,
\[ \mathbf{T}_{n} = \left[ \mathbf{t}_{1},\mathbf{t}_{2},\ldots, \mathbf{t}_{m} \right], \]
\[ \mathbf{P}_{n} = \left[ \mathbf{p}_{1}, \mathbf{p}_{2},\ldots, \mathbf{p}_{m} \right], \]
\[ \mathbf{B}_{n} = diag{(\mathbf{b}_{1},\mathbf{b}_{2},\ldots,\mathbf{b}_{m})}, \]
\[ \mathbf{Q}_{n} = \left[ \mathbf{q}_{1},\mathbf{q}_{2},\ldots,\mathbf{q}_{m}\right].\]
\footnote{Va notato che il numero dei fattori qui considerato è pari a alla dimensione dell'ingresso cioè $m$ e corrisponde a condiderare tutte le variabili latenti calcolate dall'algoritmo PLS. Si cercherà nella fase sperimentale di limitare il numero di variabili latenti utilizzate attraverso i calcolo del RMSECV e RMSEP per diversi set di ingressi e uscite.}
Le matrici di ingresso e di uscita possono essere riscritte in relazione alle matrici ottenute con PLS inn modo che:
\[ \mathbf{X}_{n} = \mathbf{T}_{n}\mathbf{P}^{T}_{n}\]
\[ \mathbf{Y}_{n} = \mathbf{T}_{n}\mathbf{B}_{n}\mathbf{Q}^{T}_{n} + \mathbf{F}_{n}\]
Dato che il numero di fattori considerato in questo caso è pari alla dimensione dell'ingresso, il residuo per la matrice $\mathbf{X}_{n}$ è scomparso, mentre in genarale la matrice residua $\mathbf{F}_{n}$ non è nulla per la presenza di rumore. Come già discusso se la matrice di covarianza $\mathbf{X}^{T}_{n}\mathbf{X}^{T}_{n}$ è singolare, e la soluzione ai minimi quadrati è mal condizionata.
Per minimizzare il residuo nell'equazione di uscita utilizziamo il criterio ai minimi quadrati generalizzandolo nella forma:
\begin{equation}
(\mathbf{X}^{T}_{n}\mathbf{X}^{T}_{n})\mathbf{C}_{n} = \mathbf{X}^{T}_{n}\mathbf{Y}_{n}
\label{xnxn}
\end{equation}
 
dato che le colonne di $\mathbf{T}_{n}$ sono mutuamente ortogonali, si può scrivere:
\[ (\mathbf{X}^{T}_{n}\mathbf{X}_{n}) = \mathbf{P}_{n}\mathbf{T^{T}_{n}}\mathbf{T}_{n}\mathbf{P}^{T}_{n} = \mathbf{P}_{n}\mathbf{P}^{T}_{n}\]

\begin{equation}
	\mathbf{X}^{T}_{n}\mathbf{Y}^{T}_{n} = \mathbf{X}^{T}_{n}\mathbf{T}_{n}\mathbf{B}_{n}\mathbf{Q}^{T}_{n} + \mathbf{X}^{T}_{n}\mathbf{F}_{n} = \mathbf{P}_{n}\mathbf{B}_{n}\mathbf{Q}^{T}_{n}
\label{PBQ} 
\end{equation}

Il fatto che la matrice residua $\mathbf{F}_{n}$ è ortogonale a $\mathbf{X}_{n}$ nella regressione ai minimi quadrati è usato per derivare l'equazione (\ref{PBQ}). Sostituendo le precedenti relazioni nell'equazione (\ref{xnxn}), otteniamo:
\[ (\mathbf{P}_{n}\mathbf{P}^{T}_{n})\mathbf{C}_{n} = \mathbf{P}_{n}\mathbf{B}_{n}\mathbf{Q}^{T}_{n}\]
se la matrice di covarianza $\mathbf{X}^{T}\mathbf{X}$ è singolare la matrice $\mathbf{P}\mathbf{P}^{T}_{n}$ è pure singolare. Per cui la matrice dei coefficienti del modello ARX ottenuta con l'algoritmo PLS è data generalizzando l'inversa di una matrice, attraverso la pseudo-inversa:
\[ \mathbf{C}^{PLS}_{n} = (\mathbf{P}_{n}\mathbf{B}^{T}_{n})^{+}\mathbf{P}_{n}\mathbf{B}_{n}\mathbf{Q}^{T}_{n}.\]
Si è così arrivati ad identificare il modello attraverso una procedura di tipo dinamico, che sfrutta le potenzialità dell'azione dell'algoritmo di PLS nel ricercare gli \emph{scores} di ingresso che forniscono la massima correlazione con il set di uscite.

%%%%%%%%%%%%%%%%%%%%%daniele%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



\section{La validazione}

\label{capitolovalidazione}

L'ultima, ma non meno importante, fase dell'identificazione \'e la validazione del modello ricavato. Ci si chiede se il modello identificato spieghi i dati in maniera soddisfacente o se addirittura sia il caso di riconsiderare l'intero processo di identificazione. L'analisi delle prestazioni viene effettuata in condizioni sperimentali diverse da quelle di identificazione. Una metodologia \'e la suddivisione dei dati in due insiemi:
\begin{itemize}
	\item l'uno per l'identificazione,
	\item il secondo per la validazione
\end{itemize}
Nel secondo set di dati si confrontano le uscite reali con quelle del modello con gli ingressi reali. Se tale fase comporta risultati insoddisfacenti \'e necessario rivedere il processo di identificazione, limitando il ricorso ad una nuova fase di raccolta dati che potrebbe risultare notevolmente onerosa, o non possibile.\\
 In questo caso \'e comprensibile quale rilevanza assume il design dell'esperimento di raccolta dati, in modo tale da fornire dati il pi\'u possibile ricchi di informazione per l'identificazione.

Lo scopo del simulatore in questo progetto \'e in parte diverso da quello precedente. \'E utile dal punto di vista di generazione di set di esperimenti (estate, inverno, \ldots) per validare il model\-lo, ma \'e anche vero che partendo dal presupposto che il set di rilevazioni di temperatura sia molto limitato, ossia unico per zona\footnote{zona \'e intesa come una stanza, una entrata o un'area pi\'u vasta non frazionata in tutti i lati da divisori}, e per il fatto che gli ingressi esterni siano un andamento medio e non ripartibile nelle quattro coordinate geografiche, $DesignBuilder$ non offre possibilit\'a di selezione, nel confronto con i modelli identificati con i dati reali. 
A correlazione di quello appena detto, si tenga presente che la divisione in due metodologie della fase di validazione \'e giustificata dall'uso futuro del modello. Ossia se il modello verr\'a sfruttato come si presenta in catena aperta o come predittore. Nel nostro caso, se si vorr\'a conoscere gli andamenti all'interno dell'edificio avendo a disposizione anche un solo termostato all'interno si utilizza il predittore. Se invece, si avranno a disposizione solo gli andamenti della temperatura esterna, si usa il modello identificato nella maniera migliore possibile, ossia con tutti i dati disponibili.

Nel caso di un predittore si definisce ottimo se rispetta le seguenti condizioni:

\begin{itemize}
\item predittore e l'errore di predizione devono essere totalmente scorrelati.
\item la varianza dell'errore di predizione sia minima (fra tutti i predittori che soddisfano il primo punto).
\end{itemize}
 
\subsection{metriche di valutazione}

Come precedentemente accennato, la validazione di un modello si divide in due tronconi: la cross-validazione o la validazione del modello con metodi fpe, aic. La prima \'e valida qualora la quantit\'a di dati sia elevata e il limitare l'intervallo di identificazione non comporta errori apprezzabili, la seconda invece viene sfruttata in quei casi in cui il set di dati \'e limitato, ed \'e necessario, quindi usarli tutti per l'identificazione. La verifica della bont\'a del modello avviene con i criteri fpe o akaike su un set di modelli identificati dai stessi dati a cui vengono individualmente cambiati i parametri e ordine.

\subsubsection{covarianza e correlazione}

La covarianza $cov(x_1,x_2)$ in questo caso campionaria, \'e un indice della relazione lineare esistente tra due variabili, cio\'e della loro variabilit\'a relativa. La covarianza pu\'o assumere valori diversi a seconda dell'unit\'a di misura con cui vengono riportate le variabili, e di conseguenza non \'e un buon indicatore. In generale, comunque, un valore assoluto elevato indica l'esistenza di una forte correlazione tra due variabili, viceversa nel caso il valore della covarianza tenda a zero. Mentre la correlazione $corr(x_1,x_2)$ di due variabili aleatorie \'e un indice adimensionale, compreso tra $-1$ e $1$, della relazione lineare esistente tra due variabili normalizzate. Un valore di correlazione positivo indica che le due variabili tendono a crescere o a decrescere simultaneamente, mentre un valore negativo indica che tendenzialmente le due variabili hanno carattere opposto. Infine, due variabili linearmente indipendenti l'una dall'altra hanno correlazione nulla. Ci\'o non significa che le due variabili sono indipendenti tra loro, cioè che sono ortogonali, ma solo che non esiste nessuna relazione lineare in grado di esprimere una variabile in funzione dell'altra; può però esistere una funzione di qualsiasi altro tipo in grado di mettere in relazione le due variabili. Queste definizioni permettono di capire ulteriormente il test di bianchezza e una sua contro-verifica. Quest'ultima viene attuata sempre per via simulativa con lo scopo di verificare se effettivamente l'errore di predizione fosse bianco, o almeno il suo andamento per ogni sensore tenda a rimanere all'interno dell'intervallo di confidenza.

\subsubsection{test bianchezza}

Il test di bianchezza \'e un test per decidere se un processo $X()$ è bianco, sulla base di una sua realizzazione $x(1),x(2),\ldots,x(N)$. Nel contesto dell'identificazione
il test di bianchezza \'e utile per verificare la bont\'a di un modello identificato mediante l'analisi del residuo, cio\'e dell'errore di predizione associato al modello identificato.
Se il residuo \'e un processo bianco (a media nulla), allora il modello identificato \'e una buona descrizione della serie temporale/sistema dinamico in esame. Risulta utile anche nel valutare sulla base di una realizzazione se, una volta rimosse le eventuali componenti deterministiche, la componente stocastica di una serie temporale è un processo bianco.

La funzione in Matlab \emph{resid} calcola delle stime della funzione di correlazione del residuo, e plotta un intervallo di confidenza \footnote{L'intervallo di confidenza corrisponde al range di valori dei residui che data la loro specifica probabilit\'a associata, per il sistema si possono ritenere statisticamente insignificanti.} all'interno del quale tutte le correlazioni eccetto la prima (la varianza) si dovrebbero trovare. Se cos'\'i risulta, possiamo ritenerci soddisfatti del modello identificato. Nel caso invece, alcune residui siano significativamente al di sopra della soglia, il modello relativo non \'e corretto, i residui non sono incorrelati ed \'e auspicabile una modifica. \'E il caso in cui l'ordine del modello non sia corretto, o la stessa struttura o la lunghezza della sequenza dei dati non sia appropriata.  

Il test dei residui deve essere usato con cautela ponendo particolare attenzione a un dettaglio pratico. Se si usano metodi a sottospazi (come N4sid  appartente alla classe IV \emph{instrumental variable} che si focalizzano sulla parte ingresso-uscita $F(z)$ nella relazione \ref{Gi}):

\begin{equation}
 \textbf{e}(t) = G(z)^{-1}[\textbf{y}(t)-F(z)\textbf{u}(t)]
 \label{Gi}
\end{equation}

viene consigliato in questa situazione di verificare l'indipendenza dell'errore $e$ e di $u$, dando meno importanza al test di bianchezza per l'errore. Si sposta cos\'i l'attenzione sulla dinamica di $G(z)$ e non sulle propriet\'a del segnale attraverso $F(z)$ (maggiori informazioni su \cite{dispensePicci}). Come conseguenza di ci\'o \'e lecito aspettarsi che i residui non soddisfino la propriet\'a menzionata prima di modello bianco. Non \'e comunque il nostro caso perch\'e il sistema che usiamo \ref{stoc} per identificare \'e di tipo \emph{forward innovation} dove il disturbo \'e esplicitamente modellato come rumore bianco $\textbf{e}(t)$, per cui il test di bianchezza sar\'a uno strumento di analisi ampiamente sfruttato in validazione. 

Partendo dall'assunto che $\epsilon(t)$ sia bianco e la sua funzione di covarianza $r_{\varepsilon}(\tau)=0 \quad \tau\neq0$, una volta che si costruisce la stima della funzione di covarianza:

\begin{equation}
 \hat{r}_{\varepsilon}(\tau) = \frac{1}{N}\sum_{t=1}^{N-\tau}\epsilon(t+\tau)\epsilon(t)
 \label{cov}
\end{equation}

Essendo validi gli assunti per $N \rightarrow \infty$:
\begin{eqnarray}
  \hat{r}_{\varepsilon}(\tau) \rightarrow 0 \qquad \tau\neq0\\
 \hat{r}_{\varepsilon}(0) \rightarrow \lambda^2 = E\varepsilon^2(t) 
  \label{epsilon}
\end{eqnarray}

Si considera cos\'i la quantit\'a normalizza per il test il rapporto:

\begin{equation}
 x_{\tau}=\frac{\hat{r}_{\varepsilon}(\tau)}{\hat{r}_{\varepsilon}(0)}
 \label{cov1}
\end{equation}


\subsubsection{test indipendenza}

In relazione con il test di bianchezza, il test di indipendenza fornisce un esito positivo riguardo la bont\'a del modello, nel caso i residui sia incorrelati con gli ingressi passati. Eventuali picchi fuori dall'intervallo di confidenza indicano che il modello non \'e descritto come parte di quell'uscita relativa al corrispondente ingresso analizzato. Per esempio, un picco al ritardo (lag) $k$-esimo significa che l'uscita $y(t)$ originata dall'ingresso $u(t-k)$ non \'e propriamente descritta dal modello. Questo test si basa sull'assunto $\varepsilon$ sia indipendente dagli ingressi passati e di media nulla. Se i residui e gli ingressi sono indipendenti, allora

\begin{equation}
 r_{\varepsilon u}(\tau)=E\varepsilon(t+\tau)u(t) = 0
 \label{cov1}
\end{equation}

Si considerano valori sia positivi che negativi di $\tau$ per il test. Come affermato precedentemente se il modello non \'e accurato, \'e ragionevole che l'andamento di $ r_{\varepsilon u}(\tau)$ per $\tau\geq0$ sia lontano da zero, mentre per $\tau<0$ pu\'o non essere nullo dipendentemente dall'autocorrelazione di $u(t)$. Per esempio, se $u(t)$ \'e bianco allora $r_{\varepsilon u}(\tau)=0$ per $\tau<0$ anche se il modello \'e inaccurato, ma pu\'o anche accadere che, nel caso il modello sia l'esatta descrizione del sistema, ossia  
$r_{\varepsilon u}(\tau)=0$ per $\tau\geq0$ si manifesti una $r_{\varepsilon u}(\tau)\neq0$ per $\tau<0$: questa \'e la manifestazione della presenza di feedback.

La quantit\'a normalizzata usata dal test \'e nella forma:

\begin{equation}
 x_{\tau}=\frac{\hat{r}_{\varepsilon u}(\tau)}{[\hat{r}_{\varepsilon}(0)\hat{r}_u(0)]^{1/2}}
 \label{cov1}
\end{equation}

Dove:

\begin{equation}
 \hat{r}_{\varepsilon u}(\tau) = \frac{1}{N}\sum_{t=1-min(0,\tau)}^{N-max(\tau,0)}\varepsilon(t+\tau)u(t)
 \label{cov}
\end{equation}

\subsubsection{errore quadratico medio}
\label{sec:MSE}
Un indice della qualit\'a del modello stimato diffusamente usato \'e l'MSE, ossia "`mean square error"'. Dato che per scarto o residuo viene definito come $\epsilon(t) = y_i - \hat{y}_i$ Esso \'e descritto dalla formula:

\begin{eqnarray}
 MSE &=& \frac{1}{N}\!\sum_{i=1}^N (y_i - \hat{y}_i)^2 \\
     &=& Var\{\hat{y}\} + mean^2(\hat{y}-y)
 \label{err}
\end{eqnarray}

MSE \'e un buon criterio di giudizio della bont\'a di un'equazione, perch\'e \'e un indice della dispersione degli errori rispetto alla retta o iperpiano di regressione. Si ponga attenzione al fatto che si \'e fatta l'assunzione $\frac{1}{N-k-1}\simeq\frac{1}{N}$, per cui la somma dei quadrati dei residui \'e mediata sul numero elevato dei campioni senza sottrarre ad essi i gradi di libert\'a del sistema come in realt\'a la vera definizione di MSE richiama. 
L'equazione \ref{err} come traduzione dell'MSE in media e varianza, esprime come queste propriet\'a della media e della varianza siano utilizzate nell'inferenza statistica per fornire informazioni sui parametri stimati di un'equazione di regressione; ad esempio, se si ha una stima degli errori, dalla varianza degli errori si pu\'o ricavare un valore che indichi la probabilità di ogni punto di essere uguale al valore atteso e un intervallo, detto intervallo di confidenza, all'interno del quale avere il 95\% di probabilità di trovare il valore atteso. Infine si tenga presente che questo criterio assume nomi diversi, ma con medesima formulazione, nel caso siano usati con set di dati diversi, ossia vengono utilizzati i criteri MSEP (Mean Square Error of Prediction) e MSECV (Mean Square Error of Cross-Validation) che sono computazionalmente identici, ma calcolati su dati di convalida o attraverso la cross-validation dei dati.

\subsubsection{mappa zero-polo}

Partendo dal presupposto che la stima del modello possa essere intrapreso come una allocazione zero-polo della fdt dei dati, la verifica della mappa zero-polo di ogni uscita o del modello identificato, fornisce una prima indicazione sull'esattezza dell'ordine del modello usato. Infatti in caso di zeri e poli molto vicini, suggerisce la possibilit\'a di una ulteriore riduzione del modello. 

\subsubsection{fit}

Un'altro indice di valutazione sfruttato nella fase di validazione \'e la percentuale di varianza spiegata, ossia il fit in percentuale fornito in maniera automatica dalla funzione $compare.m$ di Matlab:

\begin{equation}
 fit = 100*\left[1-\frac{norm(y_h-y)}{norm(y-\overline{y})}\right]
 \label{fit}
\end{equation}

Dove $y_h$ rappresenta l'uscita del modello con gli ingressi di validazione, mentre $y$ sono le uscite di validazione e la norma di solito \'e la norma due. Da \ref{fit} si evince che se l'uscita $h$-esima di validazione \'e uguale a quella in uscita del modello, la differenza tra le due \'e nulla e l'attinenza \' del 100\%. Essendo il sistema MIMO questo calcolo viene effettuato per ogni uscita in maniera automatica sempre dalla funzione $compare.m$, e verr\'a tenuto in considerazione nella valutazioni sperimentali, anche nell'analisi predittiva. 

\subsubsection{Fpe e Aic} 

I criteri per la scelta dell'ordine del modello in cui si cerchi di valutare la bont\'a del modello rispetto alle diverse possibili realizzazioni dei dati, assumono diverse connotazioni in relazioni al modello e alla tecnica di stima usata. In via generale, tali criteri vengono ottenuti penalizzando con una decrescita di una funzione perdita $V_N(\hat{\Theta}_N)$ in funzione dell'ordine o del set di modelli. La scelta del modello \'e associata al valore minimo dell'indice relativo. Una forma generale \cite{stoika} di questo criterio \'e:

\begin{equation}
 W_N = V_N(\hat{\Theta}_N)[1+\beta(N,p)]
 \label{fpe}
\end{equation}

Dove $\beta(N,p)$ \'e funzione di N e dei parametri p, penalizzando modelli troppo complessi. Nel caso PEM, una forma realizzabile \'e di tale criterio \'e:

\begin{equation}
 W_N = V_N(\hat{\Theta}_N)\frac{1+p/N}{1-p/N}
 \label{fpe2}
\end{equation}

Nell'ipotesi di N sufficientemente grande le approssimazioni portano ad un corrispondente criterio, denominato \emph{Akaike information criterion} (AIC).

\begin{equation}
 W_N = NlogV_N(\hat{\Theta}_N) + 2p
 \label{fpe3}
\end{equation}

Si tratta di un criterio di scelta della complessit\'a del modello di tipo statistico. Viene determinato minimizzando la distanza di Kullback tra la densit\'a di probabilit\'a dei dati osservati e quella che produrrebbe il modello in esame. 
Nel progetto in esame, esiste una routine che sfrutta tali criteri per determinare l'ordine del modello iterando pi\'u modelli ARX con differenti ordini. Con tale ordine si setta l'orizzonte passato/futuro da fornire al funzione PBSIDopt. Si richiama il fatto che la regola per la selezione dell'ordine viene detta ed \'e consistente se la probabilit\'a di selezionare l'ordine sbagliato tende a zero quanto pi\'u i dati tendono all'infinito.
 

\section{Prove simulative}
\label{sec:sim}

\subsection{L'identificazione e la validazione con DesignBuilder}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% parte che non compilo %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%fine Dany simulazione %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Per identificare il modello grazie al simulatore si hanno a disposizione svariati set di esperimenti. Diverse sono le prove che si apprestano ad eseguire, con lo scopo di fornire risposte alle seguenti domande:

\begin{enumerate}
	\item In quale periodo dell'anno l'identificazione permette di ottenere un buon modello?
	\item Qual'\'e il periodo minimo per una buona identificazione? Con quale campionamento?
	\item Quanto il modello identificato con gli esperimenti della stagione, per esempio estiva, \'e utile e corretto per predire le condizioni in altro periodo?
\end{enumerate}

All'atto dell'identificazione, nel caso le variabili siano tra loro linearmente correlate l'inversione di alcune matrici all'interno dell'algoritmo, pu\'o portare a un cattivo condizionamento numerico creando problemi. \'E possibile risolvere ove questo accada, utilizzando metodi alternativi per la risoluzione di sistemi lineari di equazioni come la decomposizione QR.

\subsection{modello estivo}

Le simulazioni per il periodo estivo sono state effettuate nelle stesse due settimane di raccolta dati reali. Ovviamente lo stesso simulatore comporta delle approssimazioni, la prima della quale \'e riconducibile al fatto che i dati luogo che il modello sfrutta sono relativi ad una zona nel trevigiano. Si \'e cercato invano di immettere i dati relativi delle stazioni meteorologici patavini, ma anche se fosse stato possibile non \'e comunque necessario ai nostri fini. Non si confrontano infatti gli andamenti, bens\'i gli ordini del modello identificati e i relativi MSE in fase di validazione. 

La prima analisi effettuata riguarda i segnali di ingresso. Per esempio l'ingresso temperatura esterna si dimostra un processo stazionario, come le altre grandezze. \'E persistentemente eccitante di ordine 50 (dato fornito dalla funzione matlab \emph{advice()} \cite{IDtool}), per cui non dovrebbero sorgere problemi di identificabilit\'a, dato l'ordine del modello limitato. Infatti per identificare i parametri di un modello di ordine $n$, di solito l'ingresso deve essere persistentemente eccitante (PE) di ordine $2n$, ma questo valore indicativo varia da metodo a metodo. Sfruttando la trasformata di Fourier del segnale si mette in luce la presenza di componenti stagionali come somma di sinusoidi. Infatti si visualizzano dei picchi presenti nella trasformata che indicano in questo modo le pulsazioni delle sinusoidi che contribuiscono alla componente deterministica. Il grafico \ref{fft_u2} rappresenta la temperatura esterna come ingresso e la sua densit\'a spettrale di potenza. Si verifica che il segnale \'e PE di ordine $n$ se la sua densit\'a spettrale \'e diversa da zero per almeno $n$ punti. Se poi il segnale ha densit\'a spettrale per tutte le frequenze strettamente positiva, garantisce persistente eccitazione per ordini ben maggiori. Il medesimo discorso  

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/due_ingressi/fft_u2.eps}
	\caption{andamento delle temperatura esterna estiva e relativa densità spettrale}
	\label{fft_u2}
\end{figure}

In figura \ref{identif_yu} vengono visualizzati la porzione dei dati utilizzati per la validazione: gli ingressi di temperatura e di umidit\'a esterna e tutte le uscite , entrambi scalati per della media. 

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/due_ingressi/identif_yu.eps}
	\caption{ingressi e uscite per l'identificazione}
	\label{identif_yu}
\end{figure}

\subsubsection{differenze tra $PBSID$ e $PBSID_opt$}

Durante il processo di identificazione, si \'e notata una traccia del diverso comportamento dei metodi nella determinazione delle matrici. Infatti per la matrice $A$ il modulo degli autovalori presenta diversit\'a per valori bassi del modulo, vedi \ref{conf_autov}. \'E possibile che il discostamento sia il risultato della diverso modo di procedere del metodo $PBSID\_opt$, che perviene ad una soluzione in maniera iterativa del sistema \ref{sistema}.

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/due_ingressi/autov_confronto.eps}
	\caption{confronto modulo degli autovalori di $A$ identificati con i metodi PBSID E PBSID\_opt}
	\label{conf_autov}
\end{figure}

Si tenga presente che per definire la bont\'a di uno stimatore dalla teoria siamo a conoscenza che la sola propriet\'a di asintotica correttezza non \' sufficiente a classificarne il comportamento, \'e necessaria anche la propriet\'a di consistenza. La propriet\'a di consistenza permette di affermare che la qualit\'a della stima migliora all'aumentare dei campioni rilevati ($N\rightarrow\infty$). Si \'e verificato che, nonostante alcune naturali discordanze come la sopracitata, entrambi i modelli assicurino tale caratteristica.

%%%%%%%%%%%%%%%%%%%%%%%% validazione in pratica %%%%%%%%%%%%%%%%%%%
%\graphicspath{{../}} % cerca le immagini nella cartella principale
%\graphicspath{{./figs/}} % cerca le immagini nella sotto-cartella "figs"

\subsubsection{validazione modello estivo}

Una volta effettuata l'identificazione del modello \'e necessario validarlo, ossia verificarne la bont\'a come ampiamente descritto nel capitolo \ref{capitolovalidazione}. 
In prima analisi nella cross-validazione si confrontano le uscite del modello con le uscite di validazione. Le prime sono le risposte del modello identificato con gli ingressi di validazione, mentre le seconde siano la parte delle uscite fornite in simulazione troncate e non utilizzate per l'identificazione. 

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/due_ingressi/uscite_validazione.eps}
	\caption{Confronto uscite di validazione con quelle di modello}
	\label{uscite_validazione}
\end{figure}

La figura \ref{uscite_validazione} fornisce solo un'idea indicativa della qualit\'a della stima. Si tenga presente che le uscite del modello sono in catena aperta, per cui \'e realistico aspettarsi simili andamenti. Infatti il modello soffre nel transitorio iniziale, dove le condizioni iniziali sono state stimate, e non riesce a descrivere le fluttuazioni intermedie. In tabella \ref{uscite} si relazionano le uscite con le stanze relative.

\begin{table}
	\centering
	\caption{relazione uscite con le stanze: mse e fit del modello in catena aperta}
		\begin{tabular}{|c|c|c|c|}
		\hline
			\textbf{uscita} n° & \textbf{stanza} & fit \tiny{[\%]} & MSE \tiny{[$10^{-3}$]}\\
			\hline
			1 & centrale termica 						& 83,46 & 39,5\\
			2 & ingresso e atrio 						& 72,98 & 71,4\\
			3 & laboratorio ricerca 2°piano & 82,57 & 68,0\\
			4 & laboratorio ricerca 1°piano & 77,13 & 34,7\\
			5 & pre-wc 											& 88,66 & 16,7\\
			6 & scala 											& 82,8  & 66,7\\
			7 & ufficio											& 77,59 & 66,5\\
			8 & bagno 2°piano 							& 69,79 & 166,2\\
			9 & bagno 1°piano 							& 77,96 & 54,1\\
		\hline
		\end{tabular}
	\label{uscite}
\end{table}

Per avere maggiori dettagli si usa l'indice MSE, come in figura \ref{MSE_catena_aperta}:

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/due_ingressi/MSE_catena_aperta.eps}
	\caption{Indice Mse per il modello estivo in catena aperta}
	\label{MSE_catena_aperta}
\end{figure}

L'errore quadratico medio per ogni uscita in condizioni in catena aperta fornisce errori molto limitati. Per esempio il massimo errore medio si localizza sull'uscita 8 relativa al bagno situato al secondo piano. Per quanto riguarda la percentuale di varianza spiegata, a titolo esemplificativo si raffigura in \ref{compare_catena_aperta} il risultato della funzione \emph{compare.m}. Da essa \'e possibile confrontare le uscite di validazione e di modelli singolarmente e trarre il relativo fit.  

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/due_ingressi/compare_catena_aperta.eps}
	\caption{Visualizzazione del comando \emph{compare} per la catena aperta}
	\label{compare_catena_aperta}
\end{figure}

%%%%%% test bianchezza catena aperta su dati validazione %%%%%%%%%%%%%%
\paragraph{test bianchezza e test indipendenza modello estivo}

Con il modello estivo identificato, le stime della funzione di correlazione ottenuta a partire dai dati e considerando i processi stazionari, non verificano il test di bianchezza. Si tenga presente che l'intervallo di confidenza \'e settato al 99\%. In figura \ref{testb_estivo}  rimangono al di fuori dell'intervallo di confidenza ,invece di tendere a zero (o al limite permanere sotto la soglia). Anche questa analisi \'e sinonimo della presenza di una componente deterministica, nonostante sia stata preventivamente tolta la media dei segnali.

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/due_ingressi/test_bianchezza.eps}
	\caption{Particolare del test di bianchezza modello estivo}
	\label{testb_estivo}
\end{figure}

Con l'ordine del modello a 21, come il toolbox di identificazione ritorna di default, il test di bianchezza fornisce quasi esito positivo. Comunque da un'analisi della mappa zero polo incrociata della fdt tra ogni uscita e ogni ingresso come in figura \ref{pzmap}, si evince che comunque ci sono cancellazioni zero-polo, indice della possibilit\'a di una riduzione dell'ordine. Il modello presenta poi i poli vicino al cerchio di raggio unitario, ossia alla pulsazione che contribuisce alla componente deterministica, che in statistica si traduce in presenza di un trend lineare. Resta comunque il fatto che la vicinanza dei poli/zeri nell'intorno di uno \'e associabile anche al tempo di campionamento, essendo che i modi del sistemi campionato, vanno come $e^{\Delta t}$.

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/due_ingressi/pzmap.eps}
	\caption{Mappa zero polo tra ingresso $u_1$ e uscita $y_1$}
	\label{pzmap}
\end{figure}

In figura \ref{testi_estivo} viene fornito una frazione del test di indipendenza. In generale l'esito \'e positivo. Le eventuali correlazioni tra i residui e gli ingressi per ritardi negativi non \'e necessariamente un indicazione di un modello non accurato. Mentre le correlazioni che dipendano dagli ingressi futuri, si traducono nella presenza di feedback nel sistema.

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/due_ingressi/test_indipendenza.eps}
	\caption{Particolare del test di indipendenza modello estivo}
	\label{testi_estivo}
\end{figure}

\paragraph{predittore modello estivo}

Nel problema dell'identificazione viene usato molto spesso, anche un approccio "predittivo"; ossia il modello \'e "buono" se \'e capace di "predire" il processo in un'arco di passi o di tempo richiesti. Dalla teoria il predittore \'e ottimo se rispetta le seguenti condizioni:

\begin{enumerate}
	\item predittore e errore di predizione devono essere totalmente scorrelati $E[\hat{y}(t|t-k)\varepsilon(t)]=0$ dove $\varepsilon=y(t)-\hat{y}(t|t-k)$.
	\item var[$\varepsilon(t)$] minima fra tutti i predittori che soddisfino il punto 1.
\end{enumerate}

Osserviamo che questa \'e una condizione necessaria ma non sufficiente di ottimalit\'a. Infatti due predittori diversi possono avverare entrambi le ipotesi di ottimalit\'a, nel caso $\varepsilon(t)$ e $\hat{y}(t|t-k)$ contengono termini che non si sovrappongono nel tempo, ovvero:

\begin{itemize}
	\item in $\varepsilon(t)$ ci sono solo termini del "futuro".
	\item in $\hat{y}(t|t-k)$ ci sono solo termini del "presente" e del "passato".
\end{itemize}

Per questa tematica \'e scelto di sviluppare il predittore a $k$ passi utilizzando sostanzialmente le funzioni di Matlab che lo implementano. Infatti la funzione $compare.m$ sfrutta la routine $predict.m$. Essa fornisce un modello in spazio di stato, dove le equazioni di stato per il predittore $k$ passi avanti hanno come valore iniziale $x(t-k)=\hat{x}(t-k)$, con $\hat{x}(t-k)$ \' stimato dal filtro di Kalman. 
Il predittore a un passo \'e utile per analizzare l'andamento dell'errore di innovazione e verificarne ulteriormente il test di bianchezza  attraverso la stima della funzione di covarianza degli stessi errori. Dalla figura \ref{e_inn estivo} si constata il fatto che questo particolare modello non avveri il test di bianchezza, proprio perch\'e alcune stime delle covarianza, non rimangono all'interno dell'intervallo, simile a quello di confidenza.

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/due_ingressi/covarianza_errore_innovazione.eps}
	\caption{Stima della funzione di covarianza dell'errore di innovazione}
	\label{e_inn estivo}
\end{figure}

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/due_ingressi/andamenti_err_innovazione.eps}
	\caption{Andamenti dell'errore di innovazione per le $n$ uscite}
	\label{e_inn_k estivo}
\end{figure}


Per i predittori $k$ passi, il loro confronto viene effettuato con l'utilizzo dell'indice MSE e della percentuale di varianza spiegata (ossia il fit). Data la notevole bont\'a del modello si dimostra inutile la visualizzazione delle uscite che per $k$ piccoli risultano pressoché sovrapposte come in figura \ref{validVspredk}. In tabella \ref{tabfitk} vengono riassunti i risultati per il predittore a diversi passi. I dati sono campionati con il minimo intervallo fornito dal simulatore, ossia 10 minuti. Per cui, per esempio, una predizione di 24 passi, equivale a una stima di quattro ore in avanti. Dal grafico \ref{fit_k_volte} \'e sensato giustificare un peggioramento della stima all'aumentare di $k$, piuttosto incuriosiscono i diversi andamenti. Le uscite $y_2$, $y_4$, $y_8$ risentono pesantemente della predizione, al contrario del set $y_5$, $y_6$, $y_7$, $y_9$. Le prime fanno riferimento ad andamenti delle temperature con diverse fluttuazioni, che nel modello simulativo della casetta sono rappresentate dalle stanze che si affacciano a est, a sud e a sud-ovest . Mentre le seconde sono gli andamenti pi\'u stazionari, riferiti alle stanze situate a nord e nord-ovest le quali non ricevendo un forte irraggiamento diretto come le prime, risentono di un'oscillazione termica pi\'u limitata. Per i fini predittivi, si vorrebbe che il modello riuscisse a spiegare meglio le dinamiche pi\'u veloci con minori errori, ma \'e anche vero che dipende da quale orizzonte di predizione venga richiesto e con quale accuratezza. Se \'e sufficiente una predizione di mezz'ora o un'ora, il fit peggiore \'e sull'ordine dell'93\%, per un periodo pi\'u lungo, per esempio fino a tre ore, si \'e costretti ad accettare un fit minimo sull'ordine del 65\%. La figure \ref{validaz_Vs_pred_15_passi} e \ref{k_pred_1y} forniscono i dettagli.

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/due_ingressi/validaz_Vs_pred_2_passi.eps}
	\caption{Confronto dati di validazione e predittore a 2 passi}
	\label{validVspredk}
\end{figure}

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/due_ingressi/fit_k_volte.eps}
	\caption{Fit di ogni uscita all'aumentare dell'intervallo di predizione}
	\label{fit_k_volte}
\end{figure}

\begin{table}
	\centering
	\caption{fit e mse di ogni uscita per diversi predittori}
		\begin{tabular}{|c|c|c|c|c|c|c|}
		\hline
			 & \multicolumn{2}{c|}{\textbf{k=1}} & \multicolumn{2}{c|}{\textbf{k=9}} & \multicolumn{2}{c|}{\textbf{k=15}}\\
			 \hline
			\textbf{n°} & fit \tiny{[\%]} & \tiny{MSE [$10^{-3}$]} & fit \tiny{[\%]} & \tiny{MSE [$10^{-3}$]} & fit \tiny{[\%]} & \tiny{MSE  [$10^{-3}$]}\\
			\hline
			%% inseriti a mano guardando le variabili MSEerr_pred e fit_predict !!
			1 &  99,62& 0,0086 & 90,87& 4,908& 85,23& 12,86\\
			2 &  99,52& 0,0130 & 90,02& 5,537& 78,65& 25,31\\
			3 &  99,58& 0,0160 & 94,11& 3,162& 87,73& 13,71\\
			4 &  99,48& 0,0085 & 91,64& 2,213& 82,57& 9,624\\
			5 &  99,78& 0,0022 & 98,44& 0,114& 96,47& 0,577\\
			6 &  99,64& 0,0124 & 97,26& 0,700& 95,27& 2,089\\
			7 &  99,48& 0,0250 & 96,02& 1,476& 93,06& 4,473\\
			8 &  99,37& 0,0412 & 87,45& 16,39& 73,30& 74,19\\
			9 &  99,57& 0,0137 & 96,54& 0,888& 94,06& 2,616\\
			\hline
		\end{tabular}
	\label{tabfitk}
\end{table}

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/due_ingressi/validaz_Vs_pred_15_passi.eps}
	\caption{Prima uscita di validazione e predittore a 15 passi}
	\label{validaz_Vs_pred_15_passi}
\end{figure}

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/due_ingressi/k_pred_1y.eps}
	\caption{Confronto validazione e predittori per l'ottava uscita}
	\label{k_pred_1y}
\end{figure}

\paragraph{variazione ingressi}

Le precedenti simulazioni sono state effettuate con due ingressi: le temperatura esterna e l'umidit\'a. In questa sezione si vuole indagare sulla qualit\'a del modello e della stima in funzione dell'utilizzo di uno, o pi\'u ingressi. Di tutti i risultati si riportano quelli pi\'u interessanti. L'ordine del modello identificato fornito di default dal tool varia da venti per l'ingresso singolo a ventiquattro per tutti e quattro ingressi a disposizione. Presupponendo la scelta dell'ordine, con il criterio di rapporto tra valori singolari adiacenti e tra il valore massimo e minimo, rimanga costante, si associa la crescita al solo aumento della quantit\'a di ingressi utilizzati.

Per un unico ingresso, ossia la temperatura esterna, l'MSE in catena aperta aumenta in generale, tranne per le uscite $y_1$, $y_4$, $y_5$. Il fit diminuisce sensibilente sull'ordine di uno o due punti percentuali. Il test di bianchezza non fornisce esito positivo, con il permanere di una situazione anomala sull'uscita $y_5$ non indagata. La mappa zero-polo delle ftd visualizza un insieme pi\'u fitto di questi attorno all'ascissa unitaria, mentre per quanto riguarda la predizione a un passo, ovviamente la stima \'e agli stessi livelli qualitativi. Comunque la stima per $k$ passi in avanti fornisce un andamento peggiore per i fit piuttosto che con quella a due ingressi.

Se si sfrutta oltre ai due ingressi precedenti, l'informazione della radiazione solare normale alla superficie (espressa in $[kW/m^2]$), la qualit\'a del modello in catena aperta peggiora, mentre ne migliora la predizione poich\'e la retroazione dello stato apporta maggiore informazione di controllo. Nel caso si usino tutti e quattro gli ingressi a disposizione, aggiungendo l'ingresso di radiazione solare diffusa nell'ambiente (espressa in $[kW/m^2]$) i risultati sostanzialmente migliorano (vedi figura \ref{yu_4u}). La combinazione lineare delle due entit\'a di radiazione solare essendo fortemente correlate, ovviamente non avverano il test di bianchezza. Mentre l'indice MSE rivela una limitazione dei valori massimi a patto che si accetti una sua redistribuzione sulle uscite con MSE piccolo (figura \ref{MSE4}). Decisamente efficace risulta la predizione, la quale mantiene un buon fit dei dati all'aumentare dei passi di predizione, come si illustrato in figura \ref{prek4}. Infine sono state effettuate simulazioni con la temperatura esterna, mantenuta sempre come riferimento minimo e le radiazioni solari. Il modello identificato peggiora ulteriormente, mentre la predizione si mantiene su buoni livelli. Con queste osservazioni in questo particolare modello estivo, \'e possibile concludere che: per il modello in catena aperta sono sufficienti la temperatura esterna e l'umidit\'a, mentre per la predizione la scelta del numero di ingressi efficace comprende sempre la temperatura esterna e necessariamente la radiazione normale con/senza la radiazione diffusa.

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/quattro_ingressi/identif_yu.eps}
	\caption{Uscite di validazione e di modello con $m = 4$ ingressi }
	\label{yu_4u}
\end{figure}

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/quattro_ingressi/MSE_catena_aperta.eps}
	\caption{MSE catena aperta con $m = 4$ ingressi }
	\label{MSE4}
\end{figure}

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_estivo_2weeks/quattro_ingressi/fit_k_volte.eps}
	\caption{Andamento fit delle uscite per k predittori nel modello con $m = 4$ ingressi }
	\label{prek4}
\end{figure}

%% quanto migliora la predizione con più dati?
\paragraph{campionatura e troncatura dati}

Per quanto riguarda il periodo minimo di campionamento per una fase di identificazione che non comporti la stima di un modello instabile, da varie prove si evince che dieci minuti \'e un periodo sufficiente per questo processo. Nel caso di periodi pi\'u lunghi, il modello identificato incorre in instabilit\'a. Per quanto riguarda il numero di campioni necessari, sempre da varie simulazioni, si aggira sempre attorno ai 500 campioni. Un minor numero di campioni produce delle forti oscillazioni delle uscite negli istanti iniziali, mentre gli indici MSE e fit si mantengono entro limiti accettabili. Sempre in tale tematica sono state effettuate delle simulazioni per l'intero arco estivo, dal 1° Luglio al 30 settembre. Non \'e utile per il caso predittivo, dove la bont\'a della stima varia con i passi di predizione richiesti, bens\'i per la catena aperta. Con una mole consistente di dati (4416 in identificazione) l'MSE permane al di sotto della soglia 0,8, mentre i fit si aggirano nell'intorno del 70\%. La sovrabbondanza dei campioni si constata a livello qualitativi nel test di indipendenza. Si conclude che non sono perci\'o necessari periodi di simulazione lunghi, piuttosto periodi sull'ordine di due settimane in stagioni diverse. Sar\'a interessante in seguito verificaree la bont\'a del modello, per esempio estivo, con i dati invernali, autunnali o primaverili.

\subsection{modello invernale}

Ora si indaga sulle differenze che possano nascere nell'identificare l'edificio con i dati delle simulazioni nel periodo invernale (dal 1 al 12 Gennaio). Il toolbox \footnote{d'ora in poi useremo la versione \emph{PBSID}} sceglie ancora un modello di ordine 21.

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_inverno/MSE_catena_aperta.eps}
	\caption{Indice Mse per il modello invernale in catena aperta}
	\label{MSE_ca_inverno}
\end{figure}

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_inverno/uscite_validazione.eps}
	\caption{Uscite di validazione e di modello invernale con ingressi $T_e$ e umidit\'a}
	\label{mod_val}
\end{figure}

\begin{figure}
	\centering
		\includegraphics[width=3.4in,height=2.7in]{./fig_modello_inverno/covarianza_errore_innovazione_uscita_1.eps}
	\caption{Andamento covarianza errore di innovazione 1° uscita modello inverno}
	\label{cov_val}
\end{figure}

Una volta identificato con la relativa routine, si convalida il modello focalizzando l'attenzione su alcune differenze sostanziali dal modello estivo. Nel caso si utilizzino gli ingressi $T_e$ e umidit\'a, il modello permane stabile, per\'o in fase di validazione fornisce risultati insoddisfacenti gi\'a da un confronto preliminare delle uscite di validazione e di modello come in figura \ref{mod_val}. L' MSE fornisce errori decisamente alti (fig.\ref{MSE_ca_inverno}) rispetto gl'altri modelli e il fit cala drasticamente all'aumentare dei passi di predizione. Tale comportamento \'e probabilmente riconducibile a forti correlazioni e al non verificarsi del test di bianchezza. A tal proposito si include una controverifica del test sempre utilizzata anche se non descritta esplicitamente nelle simulazioni; l'andamento delle covarianze dell'errore di predizione. In questo caso l'andamento dell'errore di innnovazione e della sua covarianza (fig.\ref{cov_val}), oltre che non permanente sotto la soglia del 99\%, contiene delle oscillazioni. Quest'ultime sono il possibile sinonimo della presenza di modi all'interno del modello dell'errore, cosa che idealmente non dovrebbe accadere, poich\'e per verificare il test il rumore di predizione dovrebbe essere bianco! Come nota implementativa si verificano problemi di condizioni iniziali con il metodo $PBSID\_opt$ ossia non è stato possibile stimarle e di conseguenza la validazione non funziona correttamente. Diversamente da quello che ci si poteva aspettare nel caso si utilizzino gli ingressi $T_e$ e i due relativi all'irraggiamento, la situazione si capovolge. Anche se non avvera affatto sia il test di bianchezza che quello di indipendenza, il modello ottiene buoni indici MSE e fit, permanendo per il primo sotto la soglia di 0,08 e ad di sopra dell'83\% per il secondo. Tale fatto non ha una spiegazione immediata, si pensa possa essere il risultato di una diversa e pi\'u sensibile identificazione nell'aggiungere andamenti periodici come l'irraggiamento, piuttosto che oscillatori aperiodici come l'umidit\'a.

\subsection{modello estivo nelle altre stagioni}

Come ultime simulazioni con l'ausilio del software DesignBuilder, si cerca di rispondere alla domanda antecedentemente posta, ossia quanto \'e capace di predire, per esempio il modello estivo, gli andamenti in diversi stagioni. Sfruttando la possibilit\'a di avere a disposizioni diversi set di esperimenti dal simulatore, ottengo dati delle grandezze per le stagioni invernali, autunnali e primaverili. Si usano gli ingressi di temperatura esterna e umidit\'a che forniscono maggiori prestazioni in catena aperta e leggermente inferiori per la predizione. Senza alcun indugio possiamo concludere che le procedure di identificazione e di validazione fin qui proposte hanno validit\'a solo "locale". Infatti il modello estivo, non riesce a descrivere in maniera soddisfacente gli andamenti molto diversi dal proprio set di dati usati per l'identificazione. Per il modello inverno l'MSE in catena aperta comporta errori sui 0,5 e fit attorno al 40\%. Tale situazione migliora con l'utilizzo degli ingressi di irraggiamento, l'MSE rimane al di sotto della soglia 0,2 e aumenta il fit fino al 60\% di media, cos\'i anche la predizione si mantiene su buoni livelli come si vede in figura \ref{predizione_inv} per l'ottava uscita, che in media \'e quella che introduce maggiori errori. Purtroppo si innescano indesiderati fenomeni di oscillazione per lungo periodo (in figura \ref{confronto_inv}). Medesime conclusioni si giunge per i dati autunnali e primaverili, il modello estivo non riesce a descrivere gli andamenti. Dove si raggiungono indici di qualit\'a soddisfacenti si instaurano oscillazioni indesiderate e viceversa.

\begin{figure}
	\centering
		\includegraphics[width=3in]{./fig_modello_inverno/uscite_validazione_inverno.eps}
	\caption{Confronto uscite della stagione invernale e le uscite del modello estivo con ingressi invernali}
	\label{confronto_inv}
\end{figure}


\begin{figure}
	\centering
		\includegraphics[width=3in]{./fig_modello_inverno/k_pred_1y_inverno.eps}
	\caption{Predizione delle uscite invernali con modello estivo}
	\label{predizione_inv}
\end{figure}

\subsection{DesignBuilder in queste simulazioni}

Al fine di identificare il modello con cui DesignBuilder determina le simulazioni, limitatamente all'andamento delle temperature, le scelte tipologiche e topologiche degli I/O sono la conseguenza dei seguenti ragionamenti. Come uscite il simulatore fornisce tre andamenti per le temperature: 
\begin{itemize}
\item[-]temperatura dell'aria $T_a$: calcolo della media delle temperature dell'aria. 
\item[-]temperatura radiante $T_r$: la media della temperatura radiante (MRT) della zona, calcolata assumendo la persona al centro della stanza, senza dare peso a particolari superfici.
\item[-]temperature operativa $T_{op}$: media tra la temperatura dell'aria e quella radiante.
\end{itemize}
Essendo del tutto simili a meno della media, si sceglie l'andamento $T_{op}$ della temperatura operativa, che ne riassume la media delle prime due, come espresso dalla formula $T_{op} = \frac{T_a+T_r}{2}$. Le temperature operative vengono fornite per ogni stanza, ma anche quelle medie per piano e per l'intero edificio cosa non necessaria e da non considerare all'atto dell'esportazione dei dati. Quest'ultime vengono prese come uscite del modello, mentre per quanto riguarda gli ingressi il primo utile \'e la temperatura esterna a bulbo secco, poi nel tentativo di aderenza alla metodologia con cui sono stati rilevati i dati nella realt\'a, si \'e estrapolano dal simulatore i dati di umidit\'a e irraggiamento. Quest'ultimo viene fornito scomposto in due entit\'a, irraggiamento diretto e diffuso, entrambe utilizzate in fase di identificazione. Nelle figure \ref{edita}, \ref{visualization}, \ref{temperature}, \ref{datiluogo}, vengono forniti alcuni frazioni della fase simulativa.

\begin{figure}[ht]
	\begin{center}
	\includegraphics[width=1.8in, angle=90]{./fig_designbuilder/planimetrie}
	\caption{planimetria complesso in esame.}
	\label{planimetrie}
	\end{center}
\end{figure}

\begin{figure}[ht]
	\begin{center}
	\includegraphics[width=2in, angle=90]{./fig_designbuilder/edita}
	\caption{modello edificio in DesignBuilder.}
	\label{edita}
	\end{center}
\end{figure}

\begin{figure}[ht]
	\begin{center}
	\includegraphics[width=2.3in, angle=90]{./fig_designbuilder/visualization}
	\caption{visualizzazione edificio in simulazione.}
	\label{visualization}
	\end{center}
\end{figure}

\begin{figure}[ht]
	\begin{center}
	\includegraphics[width=2.3in, angle=90]{./fig_designbuilder/temperature}
	\caption{esempio di grafici simulativi in DesignBuilder.}
	\label{temperature}
	\end{center}
\end{figure}

\begin{figure}[ht]
	\begin{center}
	\includegraphics[width=2.5in, angle=90]{./fig_designbuilder/datiluogo}
	\caption{andamento delle grandezze del luogo di riferimento.}
	\label{datiluogo}
	\end{center}
\end{figure}


\subsection{Identificazione termodinamica dell'edificio}
Nel caso che ci troviamo ad affrontare, un utilizzo delle tecniche fin'ora esposte é quanto mai necessario. Si é in presenza infatti di dati altamente correlati sia in ingresso che in uscita. Basta pensare infatti che
\begin{itemize}
\item per i \emph{sensori interni} (\emph{USCITE} del sistema): all'interno di una stanza, se come in questo caso le porte vengono tenute chiuse, le temperature misurate sono molto simili tra loro;  \\
\item per i \emph{sensori esterni} (\emph{INGRESSI} del sistema): all'esterno, le temperature dell'aria si posson ipotizzare molto simili tra loro; probabilmente l'irradiazione solare é la tipologia di ingressi selezionati che hanno maggior variabilitá, data la variazione di esposizione delle pareti della casa lungo l'arco della giornata.
\end{itemize}

É possibile affermare che, mentre la correlazione tra le uscite non crea nessun problema per quanto riguarda l'identificazione, la correlazione tra gli ingressi ne crea, e di molto influenti. Le problematiche sono quelle esposte in precedenza: se i dati sono dipendenti, la matrice degli ingressi non é invertibile, o lo é ma é anche molto vicina all'esser singolare, e questo porta a identificare un sistema molto sensibile e quindi inutilizzabile da ogni punto di vista.

Per tali simulazioni si sono utilizzate le misurazioni eseguite nell'anno precedente, prendendo in considerazione ingressi e uscite del sistema come trattati nel precedente progetto. Non è stato infatti possibile eseguire nuove misurazioni perchè la stuttura non è più agibile:tale procedura poteva esser di grande utilità perchè rilevazioni di misurazioni in un periodo dell'anno totalmente diverso dal precedente. 
In \cite{baldan} viene riportato che le misurazioni son state eseguite tenendo le porte chiuse per ogni stanza:tale procedura è stata resa necessaria perchè la parte non-lineare di un sistema di interscambio di calore è dovuta al circolo di masse d'aria, ovvero la parte relativa alla convezione. La parte di procedura che non è stata compresa è perché si sono mantenuti chiusi gli scuri e tenute aperte le finestre:per le considerazioni fatte in precedenza ciò apporta spostamento di masse d'aria anche considerevoli, e induce perciò non linearità. Si sottolinea inoltre come nella maggior parte dei sistemi di simulazione venga considerato lo scambio di calore attraverso i vetri delle finestre e l'effetto che l'irradiazione solare provoca su queste per le variazioni di temperatura interne:era stato consigliato quindi di rieseguire le misurazioni tenendo finestre chiuse e scuri aperti, e tenendo tirate le tende per permettere ai raggi solari di penetrare nelle stanze.
Inoltre, per il tempo che si disponeva, sarebbe stato interessante selezionare un periodo di grande variabilità climatica, allo scopo di trovare ingressi persistentemente eccitanti di ordine superiore a quello utilizzato: la settimana in cui sono state eseguite le rilevazioni, infatti, viene descritta per lo più con un tempo persistentemente sereno e con cambi di temperature per lo piú dolci, e non repentine. Si consiglia quindi, in caso del proseguimento di questo lavoro, di rieseguire le misurazioni in periodi di cambi stagionali (es. metà Febbraio-metà Marzo) in cui, almeno negli ultimi anni, si verificano molteplici variazioni di condizioni climatiche, anche all'interno di brevi periodi.

Si hanno quindi a disposizione:
\begin{itemize}
\item $53$ misure di temperatura interne, disposti uniformemente all'interno delle stanze (un sensore in ogni angolo + uno/due in centro per le stanze piú grandi);
\item $14$ misure di temperatura e luminosita' esterne, suddivise in $5$ di radiazione solare (uno per ogni facciata, $2$ per il lato est della casa) e $4$ sensori di temperatura. 
\end{itemize}
Già nel lavoro precedente erano state rilevati degli errori dovuti per lo più a problemi hardware dei \emph{T-mote}, che avevano costretto alla rimorzione di $2$ sensori tra i $65$ utilizzati. Da un esame più attento dei dati in uscita sono stati trovati altri sensori che soffrono di problematiche simili, nella fattispecie $5$ sensori, ma sono stati comunque mantenuti all'interno del set di misure. 
La differenza sostanziale che si vuole sottolineare da subito è che, mentre nel precedente progetto si è cercato di rimuovere "fisicamente" dei sensori, siano essi di ingresso o di uscita, tramite tecniche greedy, nel presente si sono utilizzate tecniche che consentono di rimuovere il "ridondante" dai set di misure, ma non per questo di rimuovere totalmente l'utilizzo di una specifica misura. Entrambe le tecniche che esporremo, infatti, riducono il set di ingresso o i set ingresso/uscita dei dati mediante trasformazioni lineari: i set di dati risultanti saranno frutto quindi di combinazioni lineari dell'insieme dei dati in considerazione, ma si continuerà a considerare comunque il set di dati nella sua interezza.

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.4in,height=2.6in]{./PCR/mappaSens.eps}
	\caption{Mappa tridimensionale della posizione dei sensori all'interno dell'edificio.In giallo e in verde, rispettivamente i sensori esterni di luminosità e di temperatura;in blu e in magenta i sensori disposti nello studio e nel bagno;in nero i rimanenti sensori.}
	\label{mappaSens}
\end{figure}


In figura (\ref{mappaSens}) viene riportata la disposizione fisica tridimensionale dei sensori, come specificato nella struttura dati fornita con il precedente progetto.

\subsection{PCR: l'implementazione}
\label{sec:PCRsim}
Si é implementata una funzione denominata \emph{KLT()} che compie \textit{PCA} su un set di dati qualsiasi, e tale funzione é stata utilizzata principalmente per la riduzione delle variabili in ingresso.

\begin{figure}
	\centering
	\includegraphics[width=3.4in,height=2.7in]{./PCR/ingIdent.eps}
	\caption{Ingressi del sistema, corrispondenti ai sensori esterni all'edificio.}
	\label{ingIdent}
\end{figure}


I passi che l'algoritmo esegue sono\footnote{per l'implementazione della routine si é utilizzata la \textit{scomposizione autovalori autovettori} della matrice di covarianza di$\mathbf{X}$ invece della \textit{scomposizione SVD} della matrice dei dati $\mathbf{X}$, come presentato nella sezione \ref{sec:EIG} : si é notato infatti che \emph{eig()} si comporta computazionalmente molto meglio di \emph{svd()} in presenza di una quantitá consistente di dati da trattare. In particolare, con la mole di dati a disposizione non è stato possibile compiere \textit{SVD} perché necessita di troppo dispendio di memoria fisica (almeno l'implementazione fornita da \texttt{Matlab}).}:
\begin{itemize}

\item calcola le medie campionarie e le leva ai singoli vettori di misura:
\[
 \left\{ \begin{array}{ccc}
	\mathbf{x}_c&=&\mathbf{x}-\bar{\mathbf{x}}\\
	\bar{\mathbf{x}}&=&\frac{1}{N} \sum^{N}_{i=1}x_i
  \end{array}\right.
\]

\item calcola la matrice di covarianza di $\mathbf{X}$ come:
\[
  Var{\mathbf{X}_{c}}=\frac{1} {n-1}  \mathbf{X}_{c} \mathbf{X}_{c}^{T};
\]

\item trova autovalori e autovettori di $Var{\mathbf{X}}$:
\[
  Var{\mathbf{X}_{c}} \cdot \mathbf{P} = \mathbf{P} \cdot \mathbf{D}^2;
\]
dove $\mathbf{P}$\footnote{che corrisponde poi precisamente alla matrice dei loadings (rif.\ref{sec:PCA})} é la matrice contenente gli autovettori disposti sulle colonne, e $\mathbf{D}^2$ é la matrice diagonale contenente gli autovalori della matrice $Var{\mathbf{X}}$ come spiegato in (\ref{sec:EIG});
\item ordina in maniera decrescente gli autovalori, mantenendo l'ordinamento con i relativi autovettori;
\item sceglie il numero $g$ di componenti principali secondo una determinata metrica , ad esempio quella che pesa gli autovalori trattata in (\ref{sec:SelezPC})
\item seleziona le prime $g$ colonne di $\mathbf{P}$, che rappresentano le direzioni lungo le quali i dati $\mathbf{X}$ hanno maggior rilevanza
\item proietta i dati lungo le $g$ direzioni date dai primi $g$ autovettori:
\[
  \mathbf{T}_{g}=\mathbf{X}\cdot \mathbf{P}_{g}
\]
come indicato in (\ref{TXP}).
\end{itemize}
La matrice risultante rappresenta gli \emph{scores} di $\mathbf{X}$\footnote{da qui in poi, quando si parlerà di $\mathbf{X}$ si terrà conto comunque degli ingressi centrati, anche se non esplicitato.} e fungerá da ingresso per le routine di identificazione a sottospazi al posto dell'intera matrice $\mathbf{X}$.
In questa maniera si ottiene un sistema in forma di stato come in \ref{ss} che ha peró $g$ ingressi e $p$ uscite:in caso il sistema venga utilizzato poi a scopo di controllo bisognerà comunque tener conto del cambio di base indotto dalla matrice degli scores in \ref{TXP}.

\begin{figure}
	\centering
	\includegraphics[width=3.7in,height=2.7in]{./PCR/ingPC.eps}
	\caption{componenti principali degli ingressi del sistema. Viene riportata la percentuale spiegata reale che la funzione ritorna come parametro.}
	\label{ingPC}
\end{figure}

In figura \ref{ingPC} vengono riportate le $g$ componenti principali selezionate:confrontate con la figura \ref{ingIdent}, che presenta gli ingressi del sistema rispetto al tempo, si puó notare come queste riportino le caratteristiche essenziali comuni tra l'insieme degli ingressi, e permettano quindi di fornire la stessa informazione con una mole di dati minori. Il nuovo set viene passato cosí, assieme a tutto il set di uscite, alla routine \texttt{PBSID\_opt()}. Si sono precedentemente eseguite numerose simulazioni per scegliere i parametri che indicano quanti passi considerare di passato e futuro: la scelta ottimale si sarebbe ottenuta eseguendo un ciclo e simulando vari sistemi arx, variando di volta in volta gli ordini relativi e scegliendo infine quello che riporta il minor \emph{fpe}.Purtroppo tale approccio è troppo oneroso dal punto di vista computazionale, e si è preferito un approccio greedy:i risultati migliori sono stati ottenuti quindi con $7$ passi per il passato e $4$ per il futuro.


\begin{figure}
	\centering
	\includegraphics[width=3.4in,height=2.6in]{./PCR/uscite_validazione.eps} 
	\caption{Uscite in catena aperta VS Uscite di validazione.}
	\label{outVal}
\end{figure}


\begin{figure}
	\centering
	\includegraphics[width=3.4in,height=2.6in]{./PCR/MSE_catena_aperta.eps} 
	\caption{Mean Square Error a catena aperta per il sistema identificato.}
	\label{MSEca}
\end{figure}


Viene dunque scelto un ordine del sistema leggermente inferiore di quello proposto da default, in modo da evitare l'overfitting\footnote{infatti di default la routine imposterebbe un valore pari al numero di sensori in uscita, e questo porterebbe porbabilmente ad \emph{overfitting} il sistema e ad errori molto alti in fase di validazione}. Dato il sistema identificato, si é passati alla fase di validazione del sistema \textit{in catena aperta}: innanzitutto, si é preso il set di validazione e si é simulato il sistema identificato con gli ingressi di validazione, relazionandolo successivamente con le uscite di validazione del sistema, il risultato é in figura (\ref{outVal}). Da questi si é calcolato il \textit{MSE} per ogni uscita (\ref{sec:MSE}) e si puó apprezzarne il risultato in figura \ref{MSEca}.

\begin{figure}
	\centering
	\includegraphics[width=3.4in,height=2.6in]{./PCR/covarianza_errore_innovazione.eps} 
	\caption{}
	\label{CovErrInn}
\end{figure}

I risutati che variano di più al variare degli indici di passato-futuro sono quelli riguardanti il test di bianchezza,e la correlazione errori d'uscita-ingressi. Succede difatti che il primo non viene mai superato se gli indici passato e futuro sono differenti da quelli scelti in ultima analisi, mentre viene superato, anche se solo da qualche sensore, nel caso i valori siano impostati a $7$ per il passato e $4$ per il futuro ;gli ingressi e i residui in uscita risultano comunque incorrelati e ciò indica che non si è in presenza di retroazione. Si é poi simulato l'andamento della covarianza dell'errore di innovazione:gli andamenti riportati in figura (\ref{CovErrInn}), in cui gli indici passato e futuro sono settati entrambi a $7$ passi , sottolineano la presenza di modi sinusoidali che non sono stati identificati e sono per questo presenti nel modello di innovazione\footnote{tali modi sono molto meno visibili nel caso di corretto settaggio degli indici "passato" e "futuro"}. 

%PLS FEDEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE simulazioneeeeeeeeeeeeeeeee
\subsection{PLS DINAMICO: scelta della complessità del modello ARX}
Nell'algoritmo di identificazione analizzato precedentemente vi sono due parametri che giocano un ruolo chiave nella complessità di calcolo e nella bontà di stima del modello identificato, sono:
\begin{itemize}
	\item numero di variabili latenti da considerare,
	\item numero di campioni precedenti l'istante k-esimo che vengono considerati nel modello ARX, nel nostro caso $n_{u} = n_{y} = n_{uy}$.
\end{itemize}

\subsubsection{Numero di variabili latenti}
Nell'analisi del numero di variabili latenti da considerare nel modello vengono utilizzati i criteri MSEP (Mean Square Error of Prediction) e MSECV (Mean Square Error of Cross-validation) che sono computazionalmente identici, ma calcolati sui dati di convalida o attraverso la cross-validation dei dati di calibrazione.  (maggiori dettagli in \cite{Baibing} e in \cite{Hubert})
un altro criterio di scelta è il PRESS (PRediction Error Sum of Squares) che si ottiene come somma dei quadrati delle differenze tra ogni osservazione $y_{i}$ e il valore calcolato $\hat{y}_{(i)}$ ottenuto dall'equazione di regressione costruita su tutti i punti tranne $y_{i}$. Formalmente
	\[ PRESS = \sum^{n}_{i=1}(y_{i} - \hat{y}_{(i)})^{2}\] 
Il criterio di scelta indica come modello migliore quello che presenta il valore più piccolo di PRESS. Il criterio PRESS è molto simile a MSE e MSEP ma viene utilizzato solitamente nella cross-validation, infatti si ha 
	\[ MSECV = PRESS/(n - k - 1);\] versioni diverse dello stesso criterio si hanno eliminando blocchi di dati invece che un singolo valore alla volta per la fase di calibrazione (dettagli in \cite{Tonizzo}).\\
E' da rilevare come la scelta del numero delle variabili latenti utilizzate in un modello PLS sia una questione ancora aperta in letteratura, lo dimostra la grande mole di articoli reperibili, ma tutti caratterizzati da una diversa impostazione delle ipotesi di lavoro e di verifica, in relazione allo scopo per il quale si cerca di ridurre il numero di variabili latenti.
Noi abbiamo sviluppato alcune funzioni in  \emph{Matlab} che analizzano il comportamento del sistema, in dettaglio sono: \emph{ipls} che sfrutta delle sottofunzioni per calcolare il  RMSE di ogni modello, in base agli intervalli di dati che si desidera considerare e al numero di variabili predittive del modello; e \emph{myplscv} che calcola e grafica il RMSEP. A tal proposito va sottolineato che la dimensione del problema e quindi del numero di variabili latenti massimo che possiamo considerare nel modello è per noi variabile, in quanto nel caso dinamico è funzione del:
\begin{itemize}
	\item numero di sensori di ingresso al sistema,
	\item numero di sensori di uscita considerato,
	\item ordine del modello ARX considerato $n_{uy}$.
\end{itemize}
Si riportano di seguito alcuni grafici ottenuti in fase di simulazione:

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.6in,height=2.6in]{./figurePLS/RMSECV.eps}
	\caption{Grafico di Root Mean Square Error of Cross -  Validation con 14 ingressi e 54 uscite}
	\label{fig:RMSECV}
\end{figure}
%%
In figura (\ref{fig:RMSECV}) è graficato il RMSECV nel caso di 14 ingressi e 54 uscite ottenuto applicando l'algoritmo in modo statico al set di dati disponibili. L'obiettivo è quello di un'analisi a priori della complessità del modello prima che i dati vengano organizzati in modo da formare la matrice $ \mathbf{X} \in \mathbf{R}^{nxm}$ e di un successivo confronto.

La figura (\ref{fig:RMSEP}) è ottenuta considerando i dati come in figura (\ref{fig:RMSECV}) ma in questo caso viene calcolato il RMSEP che è un indice di predizione del modello PLS. 
\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.6in,height=2.6in]{./figurePLS/RMSEP.eps}
	\caption{Grafico di Root Mean Square Error of Prediction con 14 ingressi e 54 uscite}
	\label{fig:RMSEP}
\end{figure}

Le figure (\ref{fig:RMSECV36}) e  (\ref{fig:RMSEP36}) graficano il RMSECV e il RMSEP nel caso descritto dall'algoritmo dinamico, sistemando cioè la matrice di ingresso del PLS come: 
\[ \mathbf{X} = \left[ \mathbf{x}(1),\mathbf{x}(2),\ldots, \mathbf{x}(n) \right]^{T} \quad \in \mathbf{R}^{nxm} \]

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.6in,height=2.6in]{./figurePLS/RMSECV36.eps}
	\caption{RMSECV nel caso PLS dinamico troncato a circa 36 variabili latenti}
	\label{fig:RMSECV36}
\end{figure}

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.6in,height=2.6in]{./figurePLS/RMSEP36.eps}
	\caption{RMSEP nel csao dinamico troncato a circa 36 variabili latenti}
	\label{fig:RMSEP36}
\end{figure}

Confrontando il caso statico con il caso dinamico si nota come nel caso dinamico le prime $5 \div 6$ variabili latenti spieghino "meglio" il sistema di quanto non accada nel caso statico, infatti nel caso statico il RMSECV varia da circa 1,3 a 1 mentre nel caso dinamico la variazione va da circa 1,3 a 0.05,  a dimostrazione del fatto che identificare il sistema con un procedimento dinamico porta a dei risultati più accurati.\\
Considerazioni analoghe possono essere fatte per il RMSEP. Nel seguito identificheremo il sistema considerando le prime 5 o 10 variabili latenti inn modo da percepire la variazione nell'errore di identificazione.

\subsubsection{Ordine del modello ARX, $n_{uy}$}
La scelta dell'ordine da utlizzare nel modello ARX è uno dei passi dell'algoritmo PLS di non facile analisi in quanto:
\begin{itemize}
	\item utilizzare la funzione \emph{FPE} di \emph{Matlab} richiede di disporre di un set di modelli ARX di ordine diverso costruiti con la funzione \emph{ARX}, la quale crea una struttura per ogni modello contenente l'indice FPE del modello,
	\item considerando un modello MIMO di 14 ingressi reali e 54 uscite la memoria di sistema rappresenta il collo di bottiglia in fase di simulazione, di fatto rendendola impossibile.
\end{itemize}
Si è quindi cercato, in un primo momento un approccio diverso che permettesse di superare tali difficoltà. \\
Si sono simulati 20 modelli ARX di ordine crescente da $n_{u} = n_{y} = n_{uy} = 1$ a $n_{u} = n_{y} = n_{uy} = 20$ con 14 ingressi effettivi e 54 uscite. Si sono scelti 5 sensori campione per ognuno dei quali si è calcolato il MSE al variare dell'ordine del modello e si sono quindi graficati i risultati.

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.6in,height=2.6in]{./figurePLS/MSE20ordiniVERO.eps}
	\caption{MSE di un sensore per ordine di modello }
	\label{fig:MSE20ordiniVERO}
\end{figure}

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.6in,height=2.6in]{./figurePLS/MSE20ordiniVERO2.eps}
	\caption{MSE di un sensore per ordine di modello }
	\label{fig:MSE20ordiniVERO2}
\end{figure}

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.6in,height=2.6in]{./figurePLS/MSE20ordiniVERO3.eps}
	\caption{MSE di un sensore per ordine di modello }
	\label{fig:MSE20ordiniVERO3}
\end{figure}

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.6in,height=2.6in]{./figurePLS/MSE20ordiniVERO4.eps}
	\caption{MSE di un sensore per ordine di modello }
	\label{fig:MSE20ordiniVERO4}
\end{figure}

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.6in,height=2.6in]{./figurePLS/MSE20ordiniVERO5.eps}
	\caption{MSE di un sensore per ordine di modello }
	\label{fig:MSE20ordiniVERO5}
\end{figure}

Le figure(\ref{fig:MSE20ordiniVERO}, \ref{fig:MSE20ordiniVERO2}, \ref{fig:MSE20ordiniVERO3}, \ref{fig:MSE20ordiniVERO4}, \ref{fig:MSE20ordiniVERO5} ) evidenziano come sia minima la variabilità del MSE tra i diversi modelli e in particolare come il minimo ordine vari in relazione al sensore scelto. Reputando scarso utilizzo solo di tale metodo si è cercato un modo per calcolare l'FPE dei diversi modelli identificati con l'algoritmo PLS dinamico.
Decidendo di utilizzare le funzioni fornite nel IDtoolbox di \emph{Matlab}, si è operato come segue:
\begin{enumerate}
	\item identificato il modello attraverso il PLS dinamico,
	\item estratte le matrici $A$ e $B$ dei coefficienti del modello ARX
	\item risimulato le uscite del sistema identificato utilizzando come ingressi gli stessi ingressi di identificazione del PLS,
	\item a questo punto si sono identificati nuovamente i diversi modelli con la funzione \emph{arx} utilizzando tutti gli ingressi di identificazione e come uscite quelle ottenute dalla simulazione del punto precedente.
\end{enumerate}
Tale procedimento è stato eseguito su un set ridotto di 10 uscite in luogo di 54, questo per risolvere i problemi di disponibilità di memoria in fase di simulazione.
L'analisi della figura (\ref{fig:FPE14ing10uscZOOM}) in relazioene ai precedenti grafici del MSE ci permette di dire che un buon ordine per il modello ARX è di $n_{u} = n_{y} = 9.$ \\ Va tuttavia considerato che a priori ci si aspettava una maggiore variabilità nella bontà del modello in relazione all'ordine, osservazione parzialmente smentita sia dal calcolo di MSE che di FPE.

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.7in,height=2.7in]{./figurePLS/FPE14ing10uscZOOM.eps}
	\caption{Grafico dell'FPE per ordine di modello}
	\label{fig:FPE14ing10uscZOOM}
\end{figure}

Si presentano i risultati ottenuti identificando il modello con l'algoritmo PLS dinamico simulando i risultati su dati di validazione e variando il numero di variabili latenti utilizzate. La funzione che calcola il modello ARX è \emph{ARX\_PLS} che sfrutta al suo interno la funzione \emph{PLSF2}.\\
Va notato che nella funzione \emph{ARX\_PLS} la validazione del modello è stata ottenuta con due metodologie diverse:
\begin{enumerate}
	\item Una utilizzando i dati di validazione e organizzandoli in modo da avere la matrice \[ \mathbf{X_{VALIDAZIONE}} = \left[       \mathbf{x}(1),\mathbf{x}(2),\ldots, \mathbf{x}(n) \right]^{T} \quad \in \mathbf{R}^{nxm},\] tale che l'uscita del sistema sia data come moltiplicazione per la matrice dei coefficienti
\item L'altra ricavando dalla matrice $C\_pls$ le due matrici dei coefficienti $A$ e $B$ che permettono di organizzare il modello ARX in forma canonica e simularne l'uscita utilizzando solo gli ingressi si validazione, $Uval$. Questa seconda metodologia rispecchia il reale utilizzo del modello anche per un futuro sviluppo di controllo.
\end{enumerate}
Pur essendo la prima metodologia poco significativa si è scelto di inserirla ugualmente, in modo da confrontare il grado di acuratezza del modello in presenza di conoscenza a priori del valore dell'uscita negli istanti precedenti al k-esimo, e nel caso invece più realistico di disponibilità dei soli ingressi di validazione.\\
Di seguito si presentano i risultati ottenuti in fase di validazione settando:
\begin{itemize}
	\item ordine del modello  $n_{u}=n_{y}=9,$
	\item utilizzando \textbf{5} variabili latenti
\end{itemize}
In figura(\ref{fig:54sensor9vl5met1MSE}) e figura(\ref{fig:54sensor9vl5met2MSE}) sono riportati i valori di MSE per tutti i 54 sensori di uscita, ottenuti attraverso i due metodi di validazione. 
\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.7in,height=2.7in]{./figurePLS/54sensor9vl5met1MSE.eps}
	\caption{MSE dei 54 sensori di uscita, \textbf{primo} metodo di validazione e \textbf{5~V.L.}}
	\label{fig:54sensor9vl5met1MSE}
\end{figure}

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.7in,height=2.7in]{./figurePLS/54sensor9vl5met2MSE.eps}
	\caption{MSE dei 54 sensori di uscita, \textbf{secondo} metodo di validazione e \textbf{5~V.L.} }
	\label{fig:54sensor9vl5met2MSE}
\end{figure}


In figura(\ref{fig:3sensor9vl5met1}) e in figura(\ref{fig:MSE3senor9vl5met1}) sono riportati i grafici dell'andamento di 3 sensori scelti casualmente con il relativo MSE ottenuti sistemando i dati di validazione secondo il primo metodo.
\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.7in,height=2.7in]{./figurePLS/3sensor9vl5met1.eps}
	\caption{Grafico di 3 sensori con andamento reale e  in uscita dal modello ARX ottenuto con il \textbf{primo} metodo di validazione e \textbf{5~V.L.}}
	\label{fig:3sensor9vl5met1}
\end{figure}

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.7in,height=2.7in]{./figurePLS/MSE3senor9vl5met1.eps}
	\caption{MSE dei tre sensori di figura(\ref{fig:3sensor9vl5met1}) relativo al metodo di \textbf{validazione~1} e \textbf{5~V.L.}}
	\label{fig:MSE3senor9vl5met1}
\end{figure}

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.7in,height=2.7in]{./figurePLS/1sensor9vl5met1.eps}
	\caption{Grafico di 1 sensore con andamento reale e  in uscita dal modello ARX ottenuto con il \textbf{primo} metodo di validazione e \textbf{5~V.L.}}
	\label{fig:1sensor9vl5met1}
\end{figure}


In figura(\ref{fig:3sensor9vl5met2}) e in figura(\ref{fig:MSE3senor9vl5met2}) sono riportati i grafici dell'andamento di 3 sensori scelti casualmente con il relativo MSE ottenuti con il secondo metodo di validazione.

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.7in,height=2.7in]{./figurePLS/3sensor9vl5met2.eps}
	\caption{Grafico di 3 sensori con andamento reale e  in uscita dal modello ARX ottenuto con il \textbf{secondo} metodo di validazione e \textbf{5~ V.L.}}
	\label{fig:3sensor9vl5met2}
\end{figure}

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.7in,height=2.7in]{./figurePLS/MSE3senor9vl5met2.eps}
	\caption{MSE dei tre sensori di figura(\ref{fig:3sensor9vl5met2}) relativo al metodo di \textbf{validazione~2} e \textbf{5~V.L.}}
	\label{fig:MSE3senor9vl5met2}
\end{figure}
 
\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.7in,height=2.7in]{./figurePLS/1sensor9vl5met2.eps}
	\caption{Grafico di 1 sensore con andamento reale e  in uscita dal modello ARX ottenuto con il \textbf{secondo} metodo di validazione e \textbf{5~ V.L.}}
	\label{fig:1sensor9vl5met2}
\end{figure}

Di seguito si presentano i risultati ottenuti in fase di validazione settando:
\begin{itemize}
	\item ordine del modello  $n_{u}=n_{y}=9,$
	\item utilizzando \textbf{10} variabili latenti
\end{itemize}
In figura(\ref{fig:3sensor9vl10met2}) e figura(\ref{fig:MSE3senor9vl10met2}) è riportato l'andamento di 3 sensori scelti casualmente con il relativo MSE ottenuti con il secondo metodo di validazione.

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.7in,height=2.7in]{./figurePLS/3sensor9vl10met2.eps}
	\caption{Grafico di 3 sensori con andamento reale e  in uscita dal modello ARX ottenuto con il \textbf{secondo} metodo di validazione e \textbf{10~V.L.}}
	\label{fig:3sensor9vl10met2}
\end{figure}

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.7in,height=2.7in]{./figurePLS/MSE3senor9vl10met2.eps}
	\caption{MSE dei tre sensori di figura(\ref{fig:3sensor9vl10met2}) relativo al metodo di \textbf{validazione~2} e \textbf{10~V.L.}}
	\label{fig:MSE3senor9vl10met2}
\end{figure}

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.7in,height=2.7in]{./figurePLS/1sensor9vl10met2.eps}
	\caption{Grafico di 1 sensori con andamento reale e  in uscita dal modello ARX ottenuto con il \textbf{secondo} metodo di validazione e \textbf{10~V.L.}}
	\label{fig:1sensor9vl10met2}
\end{figure}



%SELEZIONE USCITE
\subsection{Selezione dei sensori d'uscita}
Una delle problematiche fin qui non ancora affrontate riguarda la selezione dei sensori in uscita. Nel precedente lavoro sono stati utilizzati algoritmi euristici che sfruttano tecniche di ricerca greedy, locale e genetica e che sfruttano le metriche date dalla varianza dell'errore di predizione o il fit del modello risultante.


\begin{figure}
	\centering
	\includegraphics[width=3.7in,height=2.7in]{./PCR/uscite.eps} 
	\caption{Temperature interne del sistema.}
	\label{uscite}
\end{figure}


Si è più volte sottolineato come i metodi precedentemente esposti in questo lavoro permettano la riduzione del set di dati, in particolare quelli d'ingresso, ed evitino così numerosi problemi in fase d'identificazione:tali metodi non permettono però di ridurre \textsl{fisicamente} i sensori presenti nel sistema, perchè le trasformazioni indotte da questi metodi coinvolgono algebricamente tutti i processi. Richiamando la definizione di loading e score data in (\ref{sec:PCA}) e l'equazione \ref{eq2} si può infatti capire come la riduzione di dimensione sia algebrica e non fisica.


\begin{figure}
	\centering
	\includegraphics[width=3.7in,height=2.7in]{./PCR/uscitePC.eps} 
	\caption{Componenti principali delle uscite che spiegano la percentuale indicata dei dati.Da notare come queste rappresentino gli andamenti principali che compongono le uscite stesse.}
	\label{uscitePC}
\end{figure}


\begin{figure}
	\centering
	\includegraphics[width=3.7in,height=2.7in]{./PCR/uscLoading.eps} 
	\caption{Loadings definiti per le uscite selezionando 5 componenti principali.}
	\label{uscLoading}
\end{figure}


L'idea è stata a questo punto di analizzare la matrice $\mathbf{P}$ (fig. \ref{uscLoading}), che definisce il cambio di base, con riferimento ai singoli valori di ogni loadings: il k-esimo valore del loadings, difatti, pesa l'importanza che il k-esimo sensore ha per la q-esima componente principale (score) definita. La tecnica greedy adottata è stata dunque quella di selezionare, per ogni componente principale, il sensore che ha nella relativa PC il peso più alto, quindi il valore più alto all'interno del loading. \emph{Così facendo, però, viene tralasciata la restante informazione all'interno del loading, e può succedere che esistano più sensori rilevanti che definiscono il q-esimo score.Inoltre in questo modo il numero di sensori selezionati è sempre pari al numero di componenti principali definite.Ciò che fa variare il numero di sensori selezionati di volta in volta è dunque la percentuale fissata per la selezione delle componenti principali in (\ref{sec:SelezPC}).}

\begin{figure}[htb!]
	\centering
		\includegraphics[width=3.6in,height=2.6in]{./PCR/mappaSensSelez.eps}
	\caption{Mappa dei sensori.In giallo e in verde, rispettivamente i sensori di luminosità e di temperatura esterni, in rosso i sensori selezionati e in nero i rimanenti.}
	\label{mappaSensSelez}
\end{figure}

Successivamente alla selezione, si è simulato un modello predittivo tramite un filtro di Kalman\footnote{\'E stata riutilizzata la funzione implementata lo scorso anno, fk(), rispecchiando le modalità di lavoro.La scelta di lavorare in predizione è giustificata dal fatto che la selezione viene fatta a scopo di ridurre i costi di un sistema di controllo.}, allo scopo di analizzare come la scelta eseguita permetta di ricostruire il sistema originario. 


\begin{figure}
	\centering
	\includegraphics[width=3.7in,height=2.7in]{./PCR/uscKF.eps} 
	\caption{Uscite di validazione del sistema Vs Uscite simulate con filtro di Kalman sulle 'n' uscite selezionate.}
	\label{uscKF}
\end{figure}



In pratica, dopo l'identificazione eseguita in (\ref{sec:PCR}) si ha il modello dell'equazione \ref{ss};dopo la selezione si ha:

\begin{equation}
\label{sistID}
\left\{\begin{array}{ccc}
	 x(t+1)&=&Ax(t)+Bu(t)+Ke(t)\\
	 y_{l}(t)&=&C_{l}x(t)+Du(t)+e_{l}(t).
\end{array}\right.
\end{equation}
in cui si è indicato con $C_{l}$ la matrice $C$ di sistema in cui le righe relative ai sensori non selezionati son nulle.
Il relativo sistema con gli errori di ingresso e di misura scorrelati è:
\begin{equation}
\label{sistIDscorr}
\left\{\begin{array}{ccc}
	 x(t+1)&=&(A-SR^{-1}C)x(t)+Bu(t)+\tilde{w}(t)\\
	 y_{l}(t)&=&C_{l}x(t)+Du(t)+e_{l}(t).
\end{array}\right.
\end{equation}
Il filtro diventa dunque:\\\\
\textit{predizione}:\\
	\[
   \begin{array}{ccc}
	     \hat{x} (t+1|t)&=& F\hat{x} (t|t) +SR^{-1}y(t)+Bu(t);\\
	     P(t+1|t)&=& FP(t|t)F^{T}+\tilde{Q};
   \end{array}
	\]\\
\\
\textit{aggiornamento}: \\
	\[	  
	      \hat{x} (t+1|t+1)= \hat{x} (t+1|t) +L(t+1)\left[ y(t+1)-C_{l}\hat{x}(t+1|t)\right];\\
	\]      
	\[      P(t+1|t+1)= P(t+1|t)-L(t+1)C_{l}P(t+1|t);\\ \]
	\[     L(t)=P(t|t-1)C^{T}_{l}(C_{l}P(t|t-1)C^{T}_{l}+R)^{-1};
	\]\\
\textit{condizioni iniziali}:\\
	\[
   \begin{array}{ccc}
	 \hat{x} (0|-1)&=& \mu_{0};\\
	 P(0|-1)&=& P_{0}.
   \end{array}
	\]\\
Infine, le stime delle temperature non misurate vengono calcolate tramite:
\[
\hat{y}_{l}(t)=C_{l}\hat{x}(t).
\]
Per quanto riguarda le condizioni iniziali, $\hat{x} (0|-1)$ è stato settato tramite la creazione di un modello predittivo e l'utilizzo della funzione \texttt{pe()} di \texttt{Matlab}, come fatto in validazione; mentre per la $P_{0}$ si sono mantenuti i valori del precedente lavoro. Si sottolinea che si è utilizzato il sistema di stato identificato nella precedente analisi tramite PCR, e quindi un sistema che considera come ingresso gli scores degli ingressi.


\begin{figure}
	\centering
	\includegraphics[width=3.7in,height=2.7in]{./PCR/MSEselez.eps} 
	\caption{MSE per le uscite simulate con Kalman filter riportate in figura \ref{uscKF}.}
	\label{MSEselez}
\end{figure}

\newpage

\section{Conclusioni}
L'identificazione del modello avvenuta attraverso l'utilizzo di un algoritmo dinamico e sfruttando le potenzialità del PLS, ha il grande vantaggio di essere computazionalmente molto efficiente rispetto all'utilizzo della funzione di identificazione \emph{arx} del IDToolbox di \emph{Matlab} tanto più che tale funzione non porta ad alcun risultanto per scarsità di memoria di sistema in fase simulativa.\footnote{si dispone di un PC con 2 GByte di memoria RAM.} \\
Pur avendo ridotto la complessità del problema utilizzando di fatto al massimo 10 variabili latenti, il passo successivo è quello di costruire una relazione dinamica che permetta di limitare il numero di sensori utilizzati in fase di identificazione di un ordine di grandezza. L'obiettivo resta quello che la riduzione non deve peggiorare la qualità della stima della grandezza di interesse, nel nostro caso la temperatura interna.\\
Gli sviluppi futuri quindi di tale modello è quello di porlo in forma di predizione come fatto per il caso statico con PCA e valutarne la bontà predittiva, costruendo e ricercando un modo dinamico di selezionare un set ridotto di sensori sia di ingresso che di uscita, catturando la massima variabilità possibile delle grandezze che ci interessano.\\
In questa fase la nostra idea è quella di combinare le due metodologie che abbiamo analizzato nel seguente modo:
\begin{enumerate}
	\item Si esegue un'analisi delle componenti principali,
	\item si riduce attraverso PCR il numero di sensori di ingresso,
	\item si applica al set ridotto di sensori algoritmo di PLS dinamico,
	\item si testa il modello facendolo lavorare in forma di predizione.
\end{enumerate}
In letteratura si trovano soluzioni che risolvono il problema di ridurre il numero di sensori, utilizzando tecniche di posizionamento ottimo a priori per il successivo impiego di tecniche PLS, come in \cite{Muradore1}.
Questo tuttavia si allontana dalla nostra idea e dal nostro problema, noi cerchiamo un modo di ridurre i sensori utilizzati in fase di identificazione ma la cui collocazione è casuale, noi cerchiamo il set di sensori che porta il più alto grado di informazione per la fase di controllo e cerchiamo di farlo con un metodo dinamico.\\
Il lavoro svolto ha l'idea di fornire un quadro di insieme dello stato dell'arte nella soluzione del problema di identificazione termodinamica. Il lavoro si è concentrato nella ricerca di metodologie alternative rispetto al lavoro svolto lo scorso anno allo scopo di migliorare i risultati ottenuti sia in fase di identificazione che in fase predittiva. L'obiettivo che ci eravamo prefissi era quello di comparare parte dei risultati per capire dove le metodologie da noi applicate avevano reali vantaggi e dove invece il lavoro dello scorso anno permetteva di ottenere risultati più consistenti.
Purtroppo un'analisi approfondita del lavoro dello scorso anno non è stata possibile, ma d'altra parte, anche solo tramite riscontro visivo dell'andamento degli errori in catena aperta, è possibile apprezzare il netto miglioramento del sistema, sia con l'utilizzo di PCR che con PLS. Rapportando i due metodi in termini di MSE (fig.(\ref{fig:54sensor9vl5met2MSE}) e (\ref{MSEca})), si può constatare come il metodo PLS presenti MSE minori, come ci si poteva aspettare: il metodo PCR, infatti, compie elaborazioni solo sul set d'ingresso, mentre PLS estrae gli \emph{scores} degli ingressi che sono più rilevanti per le uscite.\\
Segnaliamo come sviluppo futuro l'idea avuta da R. Muradore e F. Bezzo in \cite{Muradore} di fondere il concetto e le potenzialità del PLS con il metodo di identificazione a sottospazi, cioè la regressione statica del PLS viene sostituita da il metodo di identificazione a sottospazi creando un algoritmo dinamico di identificazione di un modello in spazio di stato con ingressi e uscite altamente correlati e affetti da rumore di processo. 
Va considerato che la portata teorica e applicativa di questa idea è molto pesante e richiede una conoscenza dell'identificazione a sottospazi che al momento non possediamo. Altra questione è quella di poterlo implementare e testare nel nostro problema, passo anche questo di non facile attuazione.
Altri lavori futuri potrebbero essere lo studio e la valutazione di una versione non lineare dell'analisi delle componenti principali, come KPCA (Kernel PCA) trattato in \cite{KPCA}, oppure una versione di PCR che tiene conto della dinamica dell'errore di modello (MLPRC), riportata in \cite{Canada}. Numerose sono comunque le implementazioni trovate per PCA: nel lavoro svolto sono state trattate le piú semplici sia dal punto di vista teorico che implementativo, evitando di imbattersi in ulteriori teorie complicate.\\
Vogliamo infine sottolineare la difficoltà incontrata nell'affrontare gli argomenti trattati nel lavoro svolto, tanto più che la letteratura specialistica si basa su articoli di ricerca nell'ambito di applicazioni chimiche. Questo ci ha portato a considerare una grande mole di materiale per poter disporre di un quadro esaustivo della teoria di nostro interesse. Inoltre abbiamo dovuto ridefinire numerose volte il problema in relazione alle conoscenze che andavamo aquisendo: ritenevamo inizialmente che i metodi PCR e PLS permettessero una riduzione "fisica" del numero di sensori, mentre tali tecniche permettono, allo stato attuale, di ridurre il numero di variabili latenti o componenti principali che sono una combinazione lineare comprendente l'intero set di sensori.\\
Ci riteniamo soddisfatti del lavoro svolto, in quanto condotto in modo ordinato e descrivendo accuratamente i risultati ottenuti. Questo permetterà di riprendere agilmente il lavoro, valutandone i risultati ottenuti e i miglioramenti applicabili. 
Per quanto riguarda le possibili conclusioni inerenti alle simulazioni con l'ausilio di DesignBuilder \'e possibile affermare che per lo scopo prefisso l'utilizzo di tale simulatore si \'e dimostrato marginale. Si \'e cercato di immettere i dati meteorologici strettamente relativi alla settimana di raccolta dati, ma \'e risultato lungi dall'essere fattibile. Poi essendo il quantitativo di misurazioni limitato a una misura per stanza, non \'e possibile confrontarne qualitativamente l'andamento con quello risultante dalle eventuali selezioni dei sensori. Del resto questo simulatore risulta particolarmente utile qualora si voglia determinare un modello dell'edificio simulato, e con esso predire in maniera soddisfacente gli andamenti delle grandezze. In tal senso con il set di grandezze simulate, i risultati si sono dimostrati soddisfacenti sia in catena aperta, che per predizioni di molte ore. Si ricordi che comunque, il modello determinato ha validit\'a "locale", ossia nell'intorno del set di dati usato per identificarlo. Un possibile ampliamento sarebbe nel simulare il modello per esempio estivo, e validarne il comportamento con i dati reali dello stesso periodo. Infine si auspica l'analisi di altri software, tali da permettere, oltre che a maggiori dettagli di costruzione, l'estrapolazione di maggiori informazioni (come l'andamento dell'irraggiamento solare per almeno le quattro facciate in modo da confrontarle con le letture dei sensori)

\appendices

\section{Caratteristiche tecniche principali di DesignBuilder:}

\begin{itemize}
\item Interfaccia OPEN GL per grafica 3D fotorealistica.
\item Motore di calcolo EnergyPlus v 1.4.
\item Calcolo dinamico carichi termici e frigoriferi (ASHRAE Test certificate).
\item Simulazione dei consumi di combustibile ed energia elettrica per periodi temporali a scelta.
\item Collegamento a banche dati climatici mondiali.
\item	Calcolo della produzione della CO2 e verifica dei limiti di legge (Norma Inglese).
\item	Calcolo delle temperature reali nelle zone dell'edificio non climatizzate.
\item	Analisi corretto posizionamento pannelli solari e FV in funzione della posizione dinamica del sole nei vari periodi dell'anno.
\item	Archivi di strati elementari, strutture opache verticali e orizzontali, strutture vetrate verticali e orizzontali.
\item	Archivi di "Template" riferiti alle tipologie di attività sia di edifici,di impianti HVAC e di Illuminazione artificiale.
\item	Archivi riferiti a sistemi di schermatura interni ed esterni, ecct.
\item	Produzione di filmati AVI per simulazione ombre sull'edificio e rotazione del medesimo.
\item	Stampe in formato leggibile EXCEL e JPEG, Bitmap, DXF, ecct..
\item	Importazione/esportazione planimetrie DXF per agevolare costruzione edifici.
\item	Esportazione file per manipolazione esterna su EnergyPlus (casi complessi).
\item	Creazione zone termiche differenziate (attività,temperature,carichi endogeni,ecct).
\item	Check finale per Certificazione rendimento energetico edificio/impianto secondo normativa europea prEN 15217.
\end{itemize}

\subsection{Tipici usi del programma sono:}

\begin{itemize}
\item Valutare una serie di tipologie di facciate per quantificare gli effetti di surriscaldamento, utilizzo dell'energia e impatto visivo dell'edificio.
\item Valutare l'uso ottimale della luce naturale. Modellare sistemi di controllo di illuminazione e calcolare i risparmi in illuminazione elettrica.
\item Visualizzazione delle ombre solari in qualsiasi periodo dell'anno.
\item Simulazione termica di edifici ventilati naturalmente.
\item Calcolare le potenze di impianti di riscaldamento e condizionamento in tempi brevissimi.
\item Aiuto nella comunicazione durante conferenze o presentazioni.
\end{itemize}

% you can choose not to have a title for an appendix
% if you want by leaving the argument blank
\section{Cenni su modello di calcolo nello standard italiano}

\label{calcolo} Il modello di calcolo per la legislazione italiana fa' riferimento alla UNI EN 832, alla prEN 13790 e alla "`CTI Recommendation R30/3"'. CTI \'e il comitato termo-tecnico italiano,  l'organismo chiamato a legiferare normative sull'energia e gli impianti termo-tecnici, mantenendo sempre un riferimento al gruppo CEN. La R03/3 rappresenta una semplificazione nella valutazione dei rendimenti energetici degli edifici, tale da permettere un adeguato intervallo di confidenza. Perci\'o non permette simulazioni dinamiche, mentre semplifica le tipologie di ingresso-uscita dei simulatori. Un esempio di applicazione di questa metodologia \'e il software "BestClass" sviluppato dal politecnico di Milano. Esso stima il consumo di energia per il riscaldamento durante la stagione invernale e il DHW (consumo acqua calda). Purtroppo queste uniche caratteristiche limitano la completezza tecnica del pacchetto,  non considerando, per esempio, i consumi di energia elettrica. Per una descrizione dei metodi di calcolo del sofware \cite{performance}.

% use section* for acknowledgement
\section{Ringraziamenti}
% optional entry into table of contents (if used)
%\addcontentsline{toc}{section}{Acknowledgment}
 Ringraziamo per la disponibilità e la competenza dimostrata il prof. A. Chiuso a cui abbiamo fatto riferimento numerose volte.
 Un ringraziamento all'Ing. F.Bezzo e al prof. M. Barolo per materiale fornito e l'interesse dimostrato anche in prospettiva futura.
 
\section{Frasi celebri}
\textbf{Questo procedimento può essere criticabile: ma di certo conduce a risultati interessanti}\\
\emph{A.N. Kolmogorov}\\
\textbf{Se chiudete la porta agli errori, anche la verità ne resterà fuori.}\\
\emph{R. Tagore}\\
\textbf{\ldots esiste una mole imponente di materiale scritto da personaggi dalle dubbie credenziali scientifiche \ldots Ci\'o nonostante \'e invalso in questo settore l'uso di un linguaggio di tipo mistico-biologico che poco o niente ha a che fare col soggetto e apparentemente serve unicamente a fare audience.}\\
\emph{G. Picci}

\begin{thebibliography}{1}
\bibitem{performance} L. Tronchin, K. Fabbri, \emph{Energy performance building evaluation in Mediterranean countries: Comparison between software simulations and operating rating simulation}, Energy \& Buildings, 2007
\bibitem{2007b} A. Chiuso  \emph{The role of Vector AutoRegressive modeling in Predictor Based Subspace Identification} Automatica, 43(6), June 2007,  pp. 1034-1048.
\bibitem{2007} A. Chiuso \emph{On the relation between  CCA and predictor-based subspace identification} IEEE Transactions on Automatic Control, 52(10), 2007 
\bibitem{dispensePicci} G. Picci \emph{dispense corso di identificazione dei modelli} A.A. 2007/2008 
\bibitem{Bittanti} S.Bittanti \emph{Identificazione dei modelli e sistemi adattativi} Pitagora Editrice Bologna 2004 
\bibitem{manual} Lennart Andersson, Ulf Jönsson, Karl Henrik Johansson, and Johan Bengtsson \emph{A manual for system identification} Gennaio 2006
\bibitem{stoika}T.S\"oderstr\"om P.Soica \emph{System Identification} Prentice Hall
\bibitem{Burnham} Burnham,~A.~J.,~J.~F. MacGregor e R.~Viveros (1999). \emph{Latent Variable Multivariate Regression Modeling}, Chem.~Intll.~Lab. ~Syst.~\textbf{48},~167-180 
\bibitem{kourti}Kourti,T.e J.F. MacGregor. (1995)\emph{Process Analysis, Monitoring and Diagnosis, Using Multivariate Projection Methods,}  Chem.~Intll.~Lab. ~Syst.~\textbf{28},~3-21. 
\bibitem{Kresta} Kresta, J. V., J.F. MacGregor e T.E. Marlin (1991), \emph{Multivariate Statistical Monitoring of Process operating Performance,}Can.J.~Chem.~Engng.~\textbf{69},~35-47.
\bibitem{Tonizzo} Tonizzo~M.(2006-2007), \emph{Metodologie per la selezione delle variabili in ingresso nello sviluppo di sensori software basati su tecniche statistiche multivariate,} Barolo~M., Faccotesi~P., Tesi di laurea in Ingegneria Chimica anno accademico (2006-2007).
\bibitem{Qin} Qin S.Joe,(1993)\emph{Partial Least Squares Regression for Recursive System Identification} Fisher-Rosemount System, Inc. 1712 Centre Creek Drive, Austin
\bibitem{baldan} \emph{Identificazione termodinamica di un edificio} Toffoli, Baldan e Albertin  A.A. 2006/2007.
\bibitem{Wold} S. Wold, Erikson L. Trygg J. Kettaneh N. (2004)\emph{The PLS method - partial least squares projections to latent structures and its application in industrial RDP}
\bibitem{Baibing} Baibing Li, Julian M., Elaine B.M.(2002) \emph{Model selection for partial least squares regression} Chem.~Intll.~Lab. ~Syst.~\textbf{64},~79-89. 
\bibitem{Ohman} Wold, S., Annti, F. Lindgren and J. Ohman (1998). \emph{OSC of Near-infrared Spectra}, Chem.~Intll.~Lab. ~Syst.~\textbf{44},~ 175-185.
\bibitem{Trygg} Trygg, J. e S.Wold (2002). \emph{Orthogonalizated Projection to latent Structure} J. Chemom.~\textbf{16},~ 4389-4401.
\bibitem{Chong} Chong, I. G. e C.H. Jun(2005). \emph{Performance} of Some Variable selection Methods when Multicollinearity is Present, Chem.~Intll.~Lab. ~Syst.~\textbf{78},~ 103-112.
\bibitem{Qin1} S.J. Qin (1997) \emph{Recursive PLS algorithms for adaptive data modeling}, Computers and Chemical Engineering \textbf{22},~503-514.
\bibitem{Muradore} R. Muradore, F. Bezzo~(2006), \emph{A framework for PLS-SIM integration}, $14^{th}$ ~IFAC~Symposium~on~System~Identification.
\bibitem{Hubert} M.Hubert,~K. V. Branden~(2003),~\emph{Robust Methods for Partial Least Squares Regression}.
\bibitem{Muradore1} R.Muradore, F.Bezzo, M.Barolo (2006) \emph{Optimal sensor location for distributed sensor systems using multivariate regression}, Computers and Chemical Engineering ~\textbf{30}, %521-534.
\bibitem{Picci} G. Picci, \emph{Filtraggio statistico (Wiener, Levinson,Kalman) e applicazioni}
\bibitem{Stone} M. Stone, R.J. Brooks, ~\emph{Continuum regression:Cross-validated Sequentially constructed prediction embracing Ordinary Least square, Partial Least Square and Principal Component Regression}, J.R. Statist. Soc. B(1990) \textbf{52}, No. 2, pp. 237-269
\bibitem{Kemsley} E.K. Kemsley, \emph{Discriminant analysis of high-dimensional data:a comparison of principal component analysis and partial least squares data reduction methods}, Chemometrics and ILS 33(1996) 47-61
\bibitem{Roweis} S. Roweis, \emph{EM algorithms for PCA and SPCA}, Computation \& Neural systems, California Institute of Tech.
\bibitem{Hwang} J.T. Gene Hwang, Dan Nettleton, \emph{Principal Component Regression with Data-chosen components and Related methods}.
\bibitem{KPCA}B. Scholkopf, A. J. Smola, and K. R. MÄuller.\emph{Nonlinear Component Analysis as a Kernel Eigen-
value Problem}, Neural Computation, 10:1299-1319, (1998).
\bibitem{Sufian}Abu Jafar Mohammad Sufian, \emph{Analyzing Collinear Data by Principal Component Regression Approach: An Example from Developing Countries}, Journal of Data Science \textbf{3}(2005), 221-232
\bibitem{PCRfinl} P. A. Karjalainen, M. P. Tarvainen, T. Laitinen, \emph{Principal Component Regression Approach for QT
Variability Estimation}, Department of Applied Physics University of Kuopio, Finland.
\bibitem{Canada} Suzanne K. Schreyer, Michael Bidinosti,Peter D. Wentzell \emph{Application of Maximum Likelihood Principal Components Regression to Fluorescence Emission Spectra},Trace Analysis Research Centre, Department of Chemistry, Dalhousie University, Halifax, Nova Scotia, Canada B3H 4J3.
\bibitem{IDtool} \emph{Identification toolbox} di Matlab.

\end{thebibliography}

\end{document}

