Introduzione a System Identification Toolbox di MATLAB

1. POLITECNICO DI MILANO Facoltà di Ingegneria Dipartimento di Elettronica e Informazione Introduzione aSystem Identification Toolbox di MATLAB Corso: Identificazione dei Modelli e Analisi dei Dati Prof. Sergio BITTANTI

2. ORGANIZZAZIONE DELLA PRESENTAZIONE • Esempio guida • Rappresentazione e trattamento dei dati sperimentali • Calcolo del predittore • Stima non parametrica di un modello: • Stima di correlazione • Stima spettrale • Stima parametrica di un modello • Stima dei modelli col metodo dei Minimi Quadrati • Stima dei modelli col metodo della Massima Verosimiglianza • Scelta del modello ottimo e sua validazione • Analisi della complessità • Analisi dei residui: bianchezza dell’errore • Confronto di diversi modelli • Algoritmi di identificazione ricorsivi e con oblio • L’interfaccia grafica ident - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

3. Introduzione • Tipicamente un problema di identificazione inizia dopo aver progettato un set significativo di prove sperimentali ed aver raccolto una serie di dati dal sistema in esame. • In questa presentazione supporremo sempre un approccio “Black Box” al problema dell’identificazione; in altre parole, ci muoveremo come se non sapessimo nulla della fisica del sistema che genera i dati. • Se i dati sono contenuti in un file .mat (estensione dei file di dati in matlab), il primo comando da digitare nel prompt di Matlab è: “» load nomefile” • Ciò renderà i dati sperimentali disponibili nel Workspace della sessione di Matlab attiva. - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

4. Esempio Guida • Nel seguito, tutti gli esempi presentati faranno riferimento ad uno specifico sistema, che ora introduciamo. • I dati sono stati raccolti simulando un sistema AR(3), asintoticamente stabile: con . • I dati raccolti dalla simulazione sono stati memorizzati in due vettori chiamati ingresso (per il segnale e(t)) ed uscita (per il segnale y(t)), ognuno contenente 10000 campioni. • Sulla base di tali dati si calcolerà il predittore e, nell’ipotesi che il vero meccanismo di generazione dei dati sia ignoto, si vedrà come impostare la stima parametrica e non parametrica di un modello. Infine si passerà alla validazione del modello scelto. - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

5. Rappresentazione e trattamento dei dati sperimentali • Per iniziare, viene creato un oggetto iddata dove vengono salvati i dati. Questo oggetto consentirà di porre in un’unica struttura le diverse informazioni, cui si potrà accedere (come in un’usuale struttura C) con la notazione puntata nomemodello.proprietà • Il comando Matlab è: >> mymod=iddata(uscita,ingresso,1); • Con il comando >> get(mymod) si visualizzano a video tutte le proprietà dei dati. • Analogamente, il comando >> set(mymod) consente all’utente di modificare le proprietà dei dati. - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

6. Rappresentazione e trattamento dei dati sperimentali • La prima operazione da farsi sui dati prima di intraprendere la fase di identificazione vera e propria è la rimozione di eventuali trend presenti nei dati stessi. Nel caso dell’esempio guida il trend si riduce al valor medio non nullo della serie temporale. I comandi che seguono rimuovono il trend e visualizzano il risultato di tale operazione: >> figure% creo una figura in cui plottare i dati >> subplot(2,1,1)%divido la figura in 2 parti in cui plotterò uscita e ingresso originali e senza trend >> plot(ingresso(1:100));% grafico sulla finestra di grafico creata i primi 100 campioni del segnale in ingresso >> hold on;% mantengo aperta la figura corrente per plottare l’ingresso dopo il detrend nella stessa finestra >> y=detrend(ingresso);%assegno alla variabile y i valori dell’ingresso da cui viene rimosso il trend >> plot(y(1:100),'r')%grafico nella stessa finestra i campioni dell’ingresso dopo il detrend in rosso >> hold off;% rilascio la figura corrente >> subplot(2,1,2)%creo una seconda finestra nella figura >> plot(uscita(1:100));% grafico sulla finestra di grafico creata i primi 100 campioni del segnale in uscita >> hold on;% mantengo aperta la figura corrente per plottare l’usicta dopo il detrend nella stessa finestra >> u=detrend(uscita);%assegno alla variabile u i valori dell’uscita da cui viene rimosso il trend >> plot(u(1:100),'r')%gafico nella stessa finestra i campioni dell’ingresso dopo il detrend in rosso >> hold off;% rilascio la figura corrente - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

7. Rappresentazione e trattamento dei dati sperimentali- la funzione detrend - • In generale, la funzione detrend (utilizzata nell’esempio precedente per rimuovere il valor medio) offre la possibilità di rimuovere diversi tipi di trend • La sintassi completa della funzione è: >> new_mod=detrend(mod,o,breakpoint) dove: mod è un oggetto iddata; o è l’ordine del trend da rimuovere. Di default o=0: in questo caso si rimuove il trend “di ordine 0”, ossia la media del segnale. Con o=1 si possono eliminare i trend lineari, con o=2 i trend quadratici e così via; aumentando il valore del parametro o vengono rimossi trend polinomiali di ordine crescente. E’ possibile decidere anche su quale parte dei dati la funzione deve agire (è possibile infatti che la presenza del trend interessi solo una porzione dei dati). Specificando nel parametro breakpoint gli estremi della porzione di dati a cui applicare il detrend si può applicare la funzione solo alla parte dei dati di interesse. L’output della funzione new_mod è un oggetto iddata che conterrà tutti i segnali contenuti in quello di partenza mod privati del trend. - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

8. Rappresentazione e trattamento dei dati sperimentali Nella figura vediamo il risultato dell’operazione di detrend sull’esempio guida. Nella finestra in alto sono mostrati i campioni del segnale in ingresso originali (blu) e dopo il detrend (rosso). Nella finestra in basso gli analoghi andamenti dell’uscita. - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

9. Calcolo del predittore • Si vuole ora ottenere, a partire dai dati raccolti, una predizione dei valori futuri dell’uscita del processo mediante un modello (MA,AR(X),ARMA(X)) che l’utente può scegliere per descrivere i dati raccolti. • La sintassi della funzione Matlab che fa questo tipo di analisi è: >> yhat = predict(modello,dati,k,init) dove: modelloè il tipo di modello che si è scelto per descrivere il meccanismo di generazione dei dati; datiè un oggetto di tipo iddata; kè l’orizzonte predittivo; initè la modalità di inizializzazione; se init ='estimate‘ lo stato iniziale è fissato al valore che minimizza l’errore di predizione associato al modello ed ai dati; init = 'zero' fissa lo stato iniziale a zero; init = 'model' usa lo stato inziale del modello che viene automaticamente salvato alla sua creazione. L’output è: yhat, un ogetto iddata i cui OutputData sono i valori predetti. - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

10. Calcolo del predittore • A partire dai primi 5000 dati raccolti nel vettore uscita, generatiper simulazione tramite l’esempio guida, si identifica un modello AR(3). Per l’identificazione del Matlab utilizza, in questo caso, l’algoritmo dei Minimi Quadrati. Se si fosse desiderato identificare, invece, un modello ARMA, Matlab avrebbe ricorso all’algoritmo della Massima Verosimiglianza. • >> modello = ar(uscita(1:5000),[3]); • In questo modo si è identificato il modello AR(3). Questo modello può essere utilizzato per formare la predizione. In una prima fase, la predizione viene effettuata sui primi 5000 dati, gli stessi usati per l’identificazione: • >> yhat = predict(modello,[ingresso(1:5000);uscita(1:5000)],2,’estimate’); • In una seconda fase, la predizione viene effettuata sui restanti dati,così da avere un criterio oggettivo di valutazione della bontà predittiva del modello: >> compare(yhat,modello,2,5001: end); - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

11. Calcolo del predittore Nella figura si vede il confronto tra la seconda metà dei dati e valori predetti dell’uscita in base al modello AR(3) stimato. - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

12. Stima non parametrica • In questa sezione ci occuperemo di stimare alcune caratteristiche del meccanismo di generazione dati (a noi noto solo attraverso una sua realizzazione) in modo non parametrico, senza passare attraverso l’identificazione preventiva di un modello. • Il primo passo per avere informazioni sulle caratteristiche del sistema è fare una stima della funzione di covarianza campionaria. Ciò si ottiene col comando: >>Gammahat = covf(mymod,1); dove: mymod è un iddata che contiene i dati in ingresso e in uscita; 1 indica che . Il valore di questo parametro (nell’esempio 1) decrementato di 1 sarà dunque il valore di t. Con 1, dunque, stimeremo la varianza campionaria. In questo caso, con i dati dell’esempio guida, si ottiene: - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

13. Stima non parametrica • Il passo successivo della stima non parametrica è l’analisi spettrale diretta. • Un primo strumento per stimare lo spettro in questo modo è il periodogramma. A partire dai dati in esame, coi comandi >> [Pyy,w] = periodogram(uscita(1:1024)); >> plot(w,Pyy) viene calcolato e graficato il periodogramma a partire da 1024 campioni dell’uscita del sistema. Dato che, numericamente, il calcolo del periodogramma si basa sulla Fast Fourier Transform (FFT), per ottimizzare l’algoritmo è consigliabile utilizzare un numero di campioni che sia una potenza di 2. - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

14. Stima non parametrica In questa figura si vede il grafico del periodogramma calcolato su 1024 campioni dell’uscita del processo. La densità spettrale di potenza stimata viene graficata come funzione delle frequenze w, la cui unità di misura, in Matlab,è [rad/sample]. - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

15. Stima non parametrica • Com’è noto, la stima dello spettro via periodogramma soffre di diversi problemi, da cui deriva l’estrema irregolarità del grafico. • Uno dei rimedi più noti per ottenere stime più regolari dello spettro è utilizzare il metodo di Bartlett, ovvero suddividere i dati in diverse sottosequenze, calcolare il periodogramma di ognuna di esse, e ricavare la stima globale dello spettro come media di questi nuovi periodogrammi. Ciò si ottiene, in Matlab, con il comando: >> [Pyy,w]= pwelch(uscita(1:1024),bartlett(256),0); dove bartlett(256) indica che il tipo di finestra scelta è quella di bartlett (molte altre sono disponibili, ognuna ha delle caratteristiche peculiari che non discuteremo in questa sede) e che i campioni in ogni sottosequenza sono 256; 0 è il parametro con cui scegliere il numero di campioni di overlap. Il valore 0 indica che non si utilizzano finestre parzialmente sovrapposte. • Matlab offre molte altre funzioni per la stima non parametrica della densità spettrale di potenza (per una panoramica esaustiva si rimanda alla System Identification Toolbox User’s Guide). - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

16. Stima non parametrica In questa figura si vede il periodogramma stimato col metodo di Bartlett, con 4 sottosequenze da 256 dati l’una. Si può apprezzare la riduzione delle fluttuazioni. Aumentando il numero delle sottosequenze il diagramma è via via più regolare. - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

17. Stima non parametrica della risposta in frequenza • Matlab consente anche la stima della risposta in frequenza del sistema a partire dai dati. Con la funzione spa : >>g = spa (data) dove data è un vettore che contiene campioni di una realizzazione del processo in esame, si ottiene in output g, che contiene tale stima: Le frequenze cui viene calcolata la stima sono fissate di default ai valori w = [1:128]/128*pi/Ts, dove Ts è il tempo di campionamento (1 nel contesto in esame), ma possono essere anche scelte dall’utente (nel qual caso dovranno essere passate come ulteriore parametro alla funzione spa). • Con il comando >> bode (g) si ottiene il diagramma di Bode di modulo e fase della risposta in frequenza stimata. - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

18. Stima parametrica • Passiamo ora ad analizzare le tecniche di identificazione parametrica, ovvero i metodi di identificazione che mirano a descrivere il meccanismo di generazione dati ignoto mediante un modello di una data famiglia (MA, AR(X), ARMA(X), ecc.). • Cominciamo ad analizzare la procedura per la stima di un modello AR(n): >> modello = ar(data ,n, metodo) dove: data è un oggetto iddata che contiene i campioni di ingresso e uscita; n è l’ordine del modello; metodo permette di scegliere l’algoritmo da utilizzare. Le opzioni più comuni sono ‘ls’, ovvero il metodo dei Minimi Quadrati , e ‘yw’ che, per trovare la stima, risolve le equazioni di Yule Walker a partire dal calcolo delle funzioni di covarianza campionarie. - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

19. Stima parametrica • Per la stima di modelli MA e ARMA(X), si utilizza la funzione: >> modello = armax(dati,'na',na,'nb',nb,'nc',nc,'nk',nk) in cui i parametri sono quelli del modello teorico: Con un’opportuna scelta degli ordini delle diverse parti (AR, MA, X), si possono descrivere modelli AR,MA,ARMAX degli ordini desiderati. La procedura di stima si basa sull’algoritmo della Massima Verosimiglianza (MLE). - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

20. Scelta del modello ottimo e sua validazione • In base alle funzioni viste nella precedente sezione, procederemo ora a stimare diversi modelli per descrivere i dati dell’esempio. Si tratterà poi di scegliere il modello ottimo, sia dal punto di vista della complessità sia dal punto di vista della famiglia di modelli. Si sappia che i dati utilizzati nell’esempio sono stati generati da un AR(3). • Stimiamo 4 diversi modelli: un AR(1), un AR(3), un MA(3) ed un ARMA(1,2): >>m1=ar(uscita, 1, ‘ls’); >>m2=ar(uscita, 3, ‘ls’); >>m3=armax(uscita, ‘na’,0,’nc’,3 ); >>m4=armax(uscita, ‘na’,1,’nc’,2); - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

21. Scelta del modello ottimo e sua validazione • Prima di tutto occorre comprendere quale sia la famiglia di modelli cui è più probabile appartenga il meccanismo di generazione dati. Per far ciò, analizziamo l’errore di predizione, sapendo che il caso ottimo è quello in cui l’errore di predizione è un rumore bianco. A tal fine, utilizziamo una particolare funzione Matlab: >> resid(m1,uscita) >> resid(m2,uscita) >> resid(m3,uscita) >> resid(m4,uscita) Tale funzione calcola l’errore di predizione commesso dal predittore associato al modello e ne valuta la funzione di autocorrelazione. Il grafico mostra i valori di tale funzione rispetto ad un intervallo di confidenza del 99%. Se l’errore è bianco, pertanto, ci aspettiamo che tutti i valori della funzione di autocorrelazione, ad eccezione ovviamente del primo, siano nulli, o comunque sufficientemente piccoli da restare all’interno dell’intervallo di confidenza. - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

22. Scelta del modello ottimo e sua validazione La figura mostra l’analisi dei residui prima descritta nel caso del modello AR(1). Come si vede l’errore di predizione non è bianco, coerentemente col fatto che il meccanismo di generazione dati reale è un AR(3). - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

23. Scelta del modello ottimo e sua validazione La figura mostra l’analisi dei residui prima descritta nel caso del modello MA(3). Come si vede l’errore di predizione non è bianco, coerentemente col fatto che il meccanismo di generazione dati reale è un AR(3). - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

24. Scelta del modello ottimo e sua validazione La figura mostra l’analisi dei residui prima descritta nel caso del modello ARMA(1,2). Come si vede l’errore di predizione non è bianco, coerentemente col fatto che il meccanismo di generazione dati reale è un AR(3). - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

25. Scelta del modello ottimo e sua validazione La figura mostra l’analisi dei residui prima descritta nel caso del modello AR(3). Come si vede l’errore di predizione in questo caso è bianco (si vede che l’unico campione della funzione di autocorrelazione diverso da zero è il primo), coerentemente col fatto che il meccanismo di generazione dati reale è proprio un AR(3). - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

26. Scelta del modello ottimo e sua validazione • Una volta individuata la famiglia di modelli, nel nostro caso AR(n), occorre determinare la complessità ottima. Esaminiamo le funzioni Matlab che consentono questa analisi, applicandole ai modelli della famiglia AR con ordine n = 1 ,...,10. • Per fare ciò utilizziamo i criteri noti, ovvero il Final Prediction Error (FPE) e l’Akaike Information Criterion (AIC). • I comandi Matlab >> aic(modello) >>fpe (modello) consentono di calcolare i valori di tali indici per il modello identificato. - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

27. Scelta del modello ottimo e sua validazione La figura mostra l’andamento di FPE al variare dell’ordine. Coerentemente con quanto ci si aspetta, l’ordine che ottimizza questo criterio nel caso dell’esempio guida è n =3. - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

28. Algoritmi di identificazione ricorsivi e con oblio • Spesso, per ottimizzare i tempi di calcolo, è opportuno utilizzare algoritmi ricorsivi per l’identificazione, ad esempio il Recursive Least Square Algorithm (RLS). Allo stesso modo, spesso si vuole poter agire sul peso relativo dei dati, attribuendo importanza minore ai dati lontani nel tempo. Ciò si realizza introducendo il coefficiente d’oblio. • Matlab offre la possibilità di scegliere questa strada. Per esempio: >> [est,yhat] = rarx (uscita(1:length(uscita)/2),[3],’ff’,0.98 ); fornisce una stima dei parametri di un modello AR(3), costruito a partire dalla prima metà dei dati utilizzando l’algoritmo RLS. In questo caso si è usata l’opzione ‘ff’ (leggi: forgetting factor), che utilizza un coefficiente di oblio, fissato, in questo caso, a 0.98. Gli output di questa funzione sono est, che è una matrice la cui k-esima riga contiene i parametri stimati al tempo k (ovvero stimati in base ai dati fino al k-esimo incluso); yhat è il vettore che contiene la predizione dell’uscita secondo il modello stimato. Si può valutarne la bontà analizzando l’errore di tale predizione sulla seconda metà dei dati (non utilizzati nella stima del modello). - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

29. L’interfaccia grafica • Per concludere questa breve panoramica sul System Identification Toolbox, introduciamo la GUI che questo tool offre. • Per aprirla occorre scrivere al prompt >> ident Questa è l’interfaccia che si presenta all’utente - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

30. L’interfaccia grafica • Per analizzare i dati occorre importarli: il menù a tendina Data consente di importare i dati dal Workspace di Matlab. • E’ offerta all’utente la possibilità di pretrattare i dati importati, secondo le modalità indicate nel menù posto sotto alla scritta Operations (si possono rimuovere i trend, selezionare solo alcune variabili da analizzare, effettuare dei prefiltraggi, ecc.). • Dopo aver scelto i dati da analizzare si possono costruire delle stime sia parametriche sia non parametriche, come specificato nelle opzioni del menù a tendina Estimate • E’ possibile, inoltre, importare dei modelli creati precedentemente in Matlab (mediante l’opzione import del menù Models ). - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

31. L’interfaccia grafica • Scegliendo l’opzione di identificazione parametrica, per esempio, si apre una finestra di dialogo in cui scegliere il tipo di modello e il suo ordine: - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-

32. L’interfaccia grafica • Una volta creati i modelli, essi vengono visualizzati nei riquadri sotto la scritta Models. Selezionando una dellecaselle sotto la scritta Model Views si possono valutare e confrontare le caratteristiche dei diversi modelli (analizzare gli errori di predizione, valutare la risposta all’impulso e allo scalino dei vari modelli, le risposte in frequenza dei modelli stimati, ecc.) • Dunque questa interfaccia è uno strumento utile ed intuitivo per la fase iniziale di analisi in una procedura reale di indentificazione. E’ ovvio, peraltro, che ogni singolo problema avrà delle peculiarità che, quasi sempre, richiederanno un’analisi studiata ad hoc, per la quale sarà necessario creare nuove funzioni Matlab che sfruttino le potenzialità del Toolbox ma siano orientate alla specifica applicazione. - IMAD - Politecnico di Milano: System Identification Toolbox di Matlab-