Imad

Tesina di

Identificazione dei Modelli E Analisi dei Dati [96-009] Data from a flexible robot arm “Data from a flexible robot arm. The arm is installed on an electrical motor. We have modeled the transfer function from the measured reaction torque of the structure on the ground to the acceleration of the flexible arm. The applied input is a periodic sine sweep.”

Di Elisa Benetti Anno Accademico 2007/2008

1

1. Verifica delle condizioni di persistente eccitazione: Per prima cosa carichiamo i campioni sperimentali sul workspace e creiamo i vettori di ingressi, che chiamiamo u, ed uscite reali, chiamate y, calcolandone inoltre lunghezza e varianza: >>load robot_arm.dat –ascii; >>u=robot_arm(:,1); >>y=robot_arm(:,2); >> length(u) ans = 1024 >> length(y) ans = 1024 >> var(u) ans = 0.0271 >> var(y) ans = 0.0758

Verifichiamo quindi che l’ingresso sia eccitante persistentemente (PE). Per far ciò occorre creiamo la matrice di Hankel del modello. Per la creazione della matrice sfruttiamo la funzione: newhank(,,) che compie in automatico l’algoritmo. Come ordine della matrice verifichiamo i casi di ordine 10,20,30 e 40. >>PhiU10=newhank(u,1024,10); >>PhiU20=newhank(u,1024,20); >>PhiU30=newhank(u,1024,30); >>PhiU30=newhank(u,1024,30);

Creiamo quindi le matrici Cn per verificare le condizioni di non singolarità >>Cu10=PhiU10’*PhiU10/1024; >>Cu20=PhiU20’*PhiU20/1024; >>Cu30=PhiU30’*PhiU30/1024; >>Cu40=PhiU40’*PhiU40/1024;

Dalla teoria sappiamo che il segnale in ingresso u è PE di ordine n se per tale ordine la matrice Cn risulta definita e positiva. Verifichiamolo visualizzando il rango delle matrici Cn appena calcolate: >> rank(Cu10) ans = 10 >> rank(Cu20) ans = 20 >> rank(Cu30) ans = 30 >> rank(Cu40) 2

ans = 40

Risulta verificata la condizione nei casi esaminati. Calcoliamo ora gli autovalori di tali matrici >> eig(Cu10) ans = 0.0000 0.0000 0.0000 0.0001 0.0020 0.0186 0.0473 0.0539 0.0701 0.0764

>> eig(Cu20) ans = 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0003 0.0029 0.0165 0.0408 0.0436 0.0520 0.0585 0.0757 0.0768 0.0824 0.0828

>> eig(Cu30) ans = 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0004 0.0032 0.0155 0.0286 0.0343 0.0466 0.0509 0.0632 0.0673 0.0719 0.0725 0.0788 0.0793 0.0883 0.0887

>> eig(Cu40) ans = 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0005 0.0032 0.0140 0.0244 0.0262 0.0418 0.0448 0.0536 0.0558 0.0664 0.0687 0.0705 0.0706 0.0753 0.0770 0.0834 0.0840 0.0898 0.0899

Possiamo notare come le matrici siano non singolari, ma gli autovalori siano in ordine crescente, e questo è indice che sono prossime alla singolarità. Si potranno quindi riscontrare errori o warning durante la fase di stima del modello, soprattutto cercando di stimare un modello con ordine elevato.

3

2. Grafico dei dati ingresso-uscita: Apriamo l’interfaccia grafica all’identification toolbox di Matlab con il comando >> ident Opening ident ....... done.

Importiamo poi i dati:

Inseriamo come ingresso il vettore u e come uscita il vettore y, rinominiamo i dati in “robotarm” e fissiamo il tempo iniziale a 0. L’intervallo di campionamento scegliamo di lasciarlo ad 1 secondo, non essendo specificato diversamente tra i dati del Systa.

Così carichiamo il dato in Data Views:

Spuntando la casella Time Plot possiamo visualizzare il grafico dei dati Ingresso-uscita: 4

3. Grafico degli spettri: Spuntando invece Data Spectra possiamo averne un’ analisi spettrale e un Periodogramma:

4. Generazione delle sotto-sequenze di identificazione e validazione: 5

Suddividiamo ora i dati in due metà, la prima da usare per stimare il modello e la seconda per le operazione di validazione e confronto tra il modello e il sistema reale. Per fare questo selezioniamo dal menu a tendina Preprocessing, l’opzione Select Range, e selezioniamo nella finestra che si aprirà i 2 range, tra 1 e 512 il primo che chiamiamo robotarme, tra 513 e 1024 il secondo, robotarmv:

Con Time plot sui verifichiamo la divisione tra evalutation and validation data, che sommati ripresentano il segnale di partenza:

Trasciniamo quindi robotarme nella casella del Working Data e robotarmv in Validation Data e cominciamo la stima dei modelli.

6

5. Prefiltraggio dei dati: Selezionando Preprocess, effettuiamo il Remove Means e il Remove Trends:

Valutiamone quindi la differenza nel grafico degli spettri:

7

6. Stima dei modelli parametrici: 6.1. ARX: Selezionando Estimate e Parametic models si apre la seguente finestra in cui consideriamo l’opzione ARX nel menu a tendina:

Occorre inserire i valori dell’ordine na nb nk, corrispondenti all’ordine della matrice di Hankel dipendente da y, da x e dall’ordine di e. Order selection dà la possibilità che vengano analizzati i n modo automatico molti modelli differenti. Sfruttiamo questa opportunità selezionando come range di ordini 1:40, ovvero quelli valutati precedentemente.

Cambiamo il range del’asse y del grafico in quanto il modello di ordine 1 porta un Undexplained della varianza di uscita percentuale elevatissimo, da evitare. Il toolbox ident individua automaticamente, nell’intervallo da noi inserito, i modelli migliori secondo i criteri MDL, AIC e quello a Best-fit più elevato (con MDL e AIC coincidenti). Selezionando i due modelli, li inseriamo nel models views con il bottone “Insert”. Inseriamone poi un terzo, corrispondente al valore per cui l’istogramma si stabilizza

I modelli quindi che andremo a valutare sono: MDL/AIC [21 18 1], Best-Fit [32 34 1], stabilizzazione istogramma [20 17 1]. 8

Visualizziamo il loro Best-fit simulato dove il migliore risulta il modello di ordine [32 34 1], seguito dal [20 17 1], con lieve differenza di best fit, ma notevole differenza nell’ordine:

Visualizziamo poi il Best fit previsto. In questo caso i risultati sono molto simili, prenderemo allora in considerazione i modelli a ordine piu basso [20 17 1] e [18 18 1]

9

Analizziamone ora la loro bianchezza dei residui. Quelli che danno una bianchezza ‘migliore’ pur non rispettandola sono i modelli [32 34 1] e [21 18 1]. Il primo non rispetta però la cross correlazione, il secondo si. Teniamo quindi presente questo, al posto del modello [20 17 1] che avevamo tenuto in considerazione prima, data comunque la loro lieve differenza nei risultati di Best-fit :

Gli ordini valutati sono molto elevati, testando ordini più bassi troviamo che per un modello [9 9 1] risulta il migliore Best-fit simulato, al 97.34%. La bianchezza anche in questo caso non è verificata, ma la cross correlazione si. Terremo quindi conto anche di questo modello nell’analisi finale.

10

11

6.2.

ARMAX

Seleziono alla voce “Structure” nella form “Parametric models” ARMAX per stimare un modello di tipo ARMAX. Non posso stavolta utilizzare la order selection, fisso quindi gli ordini na nb nc nk.

Generalmente, si può preferire un ARMAX ad un ARX se questo stima in modo migliore il modello, con ordine inferiore a quello che faceva l’ARX. A parità di risultato e di ordine si preferisce l’ARX quindi, in quando più semplice dal punto di vista computazionale. In ogni caso decidiamo di stimare diversi ARMAX, ovvero [6 6 6 1], [8 8 8 1], [9 9 9 1], [10 10 10 1], [15 15 15 1], [21 21 21 1] e per sicurezza proviamo anche per ordini molto elevati come [30 30 30 1] fino a [40 40 40 1] che lo prendiamo come ‘limite di valutazione’. Nel caso del Best-fit Simulato abbiamo come modello migliore quello di ordine 40, ma subito dopo il [9 9 9 1]:

Nel caso del Best fit previsto abbiamo ai primi posti l’ordine 15, il 21 e di nuovo il [9 9 9 1] 12

Valutiamo infine la bianchezza. L’unico modello che presenta bianchezza è il [21 21 21] che non rispetta però la cross correlazione, il [9 9 9 1] invece rispetta la cross correlazione ed è molto vicino alla bianchezza:

13

6.3.

OE:

Nel caso dell’ Output Error proviamo vari ordini a partire dal piu basso [8 8 1], provandone di intermedi di ordine 15, 20, 30 fino al massimo [40 40 1]. Il Best fit simulato risulta il seguente, dove vediamo che al primo posto abbiamo il modello di ordine 30, e per secondo, anche se con netto calo di best-fit, il [8 8 1]:

Per quanto riguarda la bianchezza questa non è mai rispettata, ma è rispettata la cross correlazione. Prenderemo quindi i modelli migliori come Best-fit:

14

6.4.

BJ:

Ber quanto riguarda il Box-jenkins valutiamo i risultati per due ordini bassi, 2 e 8, e due ordini piu elevati, 20 e 30. Per il Best-fit simulato otteniamo un risultato accettabile solo dal modello [8 8 8 1], seguito dall’ ordine 2 ma con notevole differenza:

Per il Best-fit previsto abbiamo per primo sempre il modello [8 8 8 1] mentre quello che piu si discosta come risultato è proprio la nostra seconda scelta al punto precedente, il [2 2 2 1]:

15

Infine la bianchezza, mai verificata. Nemmeno la cross correlazione è mai verificata, ma i due modelli che possiamo definire ‘meno disastrosi’ sono il [8 8 8 1] e [30 30 30 1]:

16

7. Stima di modelli non parametrici: 7.1.

Analisi spettrale

Andiamo in Estimate, Spectral Model, ci appare la seguente finestra con possibilità di scegliere tra metodo Blackman-Tukey e Smoothed Fourier Transform:

7.2.

Correlazione

In questo caso basta semplicemente selezionare Estimate e Correlation Model per ottenere il grafico. Concludendo, i vari modelli scelti risultano quindi essere i seguenti:

17

8. Validazione Del Modello: 8.1. Confronto delle risposte al transitorio: La risposta al gradino e quella all’impulso risultano di difficile interpretazione o comunque non sono utili per fornire quelle informazioni consistenti. Infatti il robot arm è una struttura flessibile: quando viene sollecitata la struttura incomincia da oscillare, e ciò spiega il tipo di risposta visualizzata. Pertanto non ricaviamo nessuna conclusione da tale risposta, per qualsiasi modello di qualsiasi ordine.

E’ stata eliminata la risposta del modello OE [30 30 1] in quanto la sua oscillazione troppo amplificata rispetto alle altre, ne impediva un confronto anche solo visivo:

18

8.2.

Confronto delle risposte frequenziali:

19

8.3.

Confronto diagrammi poli-zeri:

OE ARX

BJ ARMAX

Nella tabella finale si trovano nello specifico la stabilità e le possibili sovra parametrizzazione per ognuno di questi modelli valutati

20

9. Calcolo incertezza parametrica dei modelli migliori: Importiamo i modelli scelti nel workspace di Matlab, trascinandoli nella casella della GUI denominata “To Workspace”. Possiamo ora estrapolare le matrici parametriche del modello attraverso il comando predefinito in Matlab “polydata”, e calcolare poi l’ incertezza parametrica relativa percentuale con max(100*abs(dX./X). Riportiamo come esempio gli output ottenuti con il primo modello, per gli altri risultati, riferirsi alla tabella finale: >> [A,B,C,D,F,dA,dB,dC,dD,dF]=POLYDATA(arx20181) A= Columns 1 through 12 1.0000 -5.7014 17.8517 -39.7568 70.8240 -107.4803 144.0422 -173.3674 188.6903 186.1854 166.9097 -136.0082 Columns 13 through 21 100.5362 -66.9765 39.7923 -20.7547 9.2959 -3.4710 1.0349 -0.2316 0.0298 B= Columns 1 through 12 0 -0.5366 3.3342 -11.1951 26.2058 -47.3174 69.7236 -86.8094 93.3896 -87.9007 72.6355 -52.4457 Columns 13 through 18 32.6482 -17.1340 7.3035 -2.3766 0.5224 -0.0553 C= 1 D= 1 F= 1 dA = Columns 1 through 12 0 0.0459 0.2635 0.8441 1.9438 3.5954 5.6659 7.8685 9.8003 11.0364 11.2731 10.4573 Columns 13 through 21 8.8059 6.7068 4.5762 2.7478 1.4145 0.6038 0.2056 0.0512 0.0071 dB = Columns 1 through 12 0 0.0040 0.0389 0.1962 0.6273 1.4669 2.7095 4.1451 5.4112 6.1392 6.1054 5.3190 Columns 13 through 18 4.0197 2.5828 1.3650 0.5617 0.1628 0.0258 dC = 0 dD = 0 dF = 0 >> max(100*abs(dA./A)) ans = 23.9532 >> max(100*abs(dB./B)) Warning: Divide by zero. ans = 46.7497 21

10.

Tabella di riepilogo: Modello Parametrico

Model

Struct

Best fit Val.

Best fit Prev.

Loss function

 

V1 ˆN

Altri criteri (FPE)

Bianche zza residuo (cross corr.)

Modello non parametrico C Spectr A analys

100 

Oss.

max |  | 

Solo a bassissime frequenze Solo a bassissime frequenze Solo a bassissime frequenze

INSTABILE SOVRAPARAM

A 23.95 B 46,75

STABILE SOVRAPARAM

A 35.57 B 260.70

Migliore best-fit

INSTABILE SOVRAPARAM

Si bianchezza

-

Solo a bassissime frequenze

STABILE SOVRAPARAM

NO,SI

-

INSTABILE SOVRAPARAM

1.81E-5

NO,SI

-

INSTABILE SOVRAPARAM

B 8.87 F 0.89

2.80E-7

4.52E-7

NO,NO

-

Solo a bassissime frequenze Solo a bassissime frequenze Solo a bassissime frequenze

A 29.06 B 49.93 C 527.41 A 35.57 B 260.69 C 70.47 B 0.06 F 0.16

STABILE NO SOVRAP

3.75E-7

4.53E-7

NO,NO

-

Solo a bassissime frequenze

INSTABILE SOVRAPARAM

B 301.81 C 18274 D 32.75 F 1047 B 68.24 C 673.20 D 46.38 F 24.08

ARX

20 17 1

97.27

99.98

2.01E-9

2.60E-9

NO,SI

-

ARX

991

97.34

99.92

4.00e-8

4.30E-8

NO,SI

-

ARMAX

21 21 21 1

96.94

99.98

2.24E-9

3.12E-9

SI,NO

-

ARMAX

9991

97.21

99.98

2.40E-9

2.77E-9

NO,SI

OE

30 30 1

97.8

97.81

8.96E-6

1.28E-5

OE

881

97.13

97.13

1.65E-5

BJ

20 20 20 20 1

-56.29

99.81

BJ

88881

95.58

99.94

 

Diag. poli-zeri

Bianchezza marginale Errore minore

No sovrapara metrizzazio ne



22



11.

Osservazioni Conclusive:



Non è possibile in nessun caso, con i dati acquisiti, ottenere una bianchezza del residuo, ma al più si più raggiungere solamente una bianchezza marginale. Pertanto il predittore non riesce ad essere ottimale.



Dobbiamo infine fare ipotesi sull’applicazione finale:

_ Se il modello serve per il progetto di un controllore del sistema, allora è necessario un buon adattamento dei dati, quindi un buon Best-fit, nel nostro caso quindi un ARX [20 17 1], o , prediligendo una complessità minore a scapito di un lieve peggioramento di best-fit, l’ ARX [9 9 9 1], considerando inoltre che a ordini maggiori corrisponde instabilità del modello. _ Se il modello è stato creato per il filtraggio dei dati, allora occorre un predittore ottimo e si darà priorità alla bianchezza dei residui, ottenibile marginalmente con un ARMAX [9 9 9 1] Data la natura del problema, l’applicazione ssarà probabilmente la prima, ovvero il progetto di un controllore che controlla l’input, ovvero la coppia fornita al braccio del robot , allo scopo effettuare al braccio un profilo di moto prestabilito. Pertanto, risulta sia sufficiente un ARX [9 9 1].

23