Modellistica ed identificazione dei processi dinamici PDF

Preview

Summary

Relazione Modellistica ed identificazione dei processi dinamici...

Description

Modellistica ed identificazione dei processi dinamici RELAZIONE DI LABORATORIO Tommasi Mauro 1048462 Modellistica ed Identificazione dei Processi Dinamici Anna Maria Perdon – Andrea Bonci A.A. 2016/2017 – Ing. Informatica e dell’Automazione

Sommario Introduzione alla Modellistica Approccio BlackBox Fase ed Ampiezza Ipotesi comportamentali a semplici funzioni White Noise, Test di Anderson Test di laboratorio BlackBox Approccio GreyBox Equazione motore DC Labview (per acquisizione/calcolo della velocità) MATLAB (per la determinazione del miglior modello)

ii

Introduzione alla modellistica Definizioni e concetti base La modellistica viene definita come la tecnica, lo studio, l’attività di creare modelli. In ambito ingegneristico rappresenta la capacità di determinare una o più funzioni che approssimino un sistema a noi ignoto (o parzialmente conosciuto) analizzando gli ingressi e le correlative uscite del sistema stesso. Nasce dunque la necessità di conoscere quel che abbiamo definito sistema. Un sistema è un insieme di equazioni rappresentanti, generalmente, il funzionamento di un sistema fisico. Consideriamo ad esempio un semplice circuito:

Per le leggi di Kirchoff denotiamo che: 𝛿𝑥! 𝑡 + 𝑅! 𝑖! (𝑡) 𝛿𝑡 𝛿𝑥! 𝑡 + 𝑥! (𝑡) 𝑅! 𝑖! 𝑡 = 𝑅! 𝐶! 𝛿𝑡 𝛿𝑥! 𝑡 𝑖! 𝑡 = 𝑥! 𝑡 − 𝐶! 𝛿𝑡 𝑢 𝑡 =𝐿

Sostituendo l’ultima relazione nelle precedenti e ponendo 𝑥 𝑡 =

!! (!) !! (!)

, definendo A, B, C

e D come matrici di dimensioni opportune che premoltiplicano lo stato 𝑥 𝑡 e 𝑢 𝑡 avremo:

𝑅! 𝑅! 𝐿 𝑅! + 𝑅! 𝐴= 𝑅! 𝐶! 𝑅! + 𝑅! −

𝛿𝑥 𝑡 = 𝐴𝑥 𝑡 + 𝐵𝑢 𝑡 𝛿𝑡

𝑦 𝑡 = 𝐶𝑥 𝑡 + 𝐷𝑢 𝑡

𝑅! 𝐿 𝑅! + 𝑅! 1 − 𝐶! 𝑅! + 𝑅!

𝑅! 𝐿 𝑅! + 𝑅! 𝐶= 𝑅! 𝑅! + 𝑅!

−

1 𝐵= 𝐿 0

1 𝐿 𝑅! + 𝑅! 1 − 𝐶! 𝑅! + 𝑅!

𝐷=

0 0

Da cui si ottiene una funzione di trasferimento (nel caso continuo con Laplace) pari a: 𝐹𝐷𝑇 = 𝐶 𝑠𝐼 − 𝐴

!!

𝐵+𝐷

Essendo D nullo e le matrici A, B e C nel dominio di Laplace osservabili e raggiungibili (scomposizione di Kalman). Nel caso continuo la funzione di trasferimento si presenta sotto la forma di un rapporto tra poli e zeri aventi differenti molteplicità: 𝐹𝐷𝑇 = 𝐶 𝑠𝐼 − 𝐴

!!

𝐵=

𝑠 − 𝑧! 𝑠 − 𝑝!

! !

oppure nel caso discreto pari a: 𝐹𝐷𝑇 = 𝐶 𝑧𝐼 − 𝐴

!!

𝐵=

𝑏!!!! 𝑧 !!!! 𝑏! 𝑧 ! + ⋯ + 𝑏! 𝑧 !! = 𝑎!!!! 𝑧 !!!! 𝑎! 𝑧 ! + ⋯ + 𝑎! 𝑧 !!

Questa piccola introduzione ci è servita per introdurre il sistema di controllo, il quale oltre ad avere una forma del tipo ingresso-uscita, considereremo quello che definiamo disturbo, ossia un rumore bianco (vedremo poi le caratteristiche). Da qui poi, in base al tipo di configurazione scelta, utilizzeremo delle classi di modelli (AR, ARX, ARMA, ARMAX, ARXAR) che ci permetteranno dunque di “modellare” il nostro sistema.

Approccio BlackBox Analisi, modellazione, ipotesi, test di laboratorio L’approccio BlackBox (scatola nera) è il sistema più complicato per identificare un modello in quanto non conosciamo in alcun modo il funzionamento interno della scatola. Quindi ponendo un opportuno ingresso avremo altrettante uscite le quali verranno poi analizzate per determinare il miglior modello. Il primo passo è quello di capire innanzitutto la natura della blackbox (elettrica, meccanica, ecc) in maniera da studiarne poi i vari effetti in uscita. Nel nostro test analizzeremo una BlackBox costituita da un circuito elettrico a noi ignoto. La strumentazione necessaria per inviare degli ingressi e ricevere delle uscite è il seguente: -

una scheda myDAQ della National Instrument che permette la ricetrasmissione delle informazioni tramite i propri canali di I/O la BlackBox contenente un pin di input ed uno di output (oltre alla massa), collegato all’alimentatore un generatore di segnali (Funzione ARB del software Elvis per la myDAQ) un visualizzatore/registratore dei segnali di ingresso/uscita un sistema di analisi per la determinazione del miglior modello (o famiglia di modelli) tramite software MATLAB

Considerando la strumentazione di laboratorio necessaria per la lettura/scrittura e considerando la natura della BlackBox, possiamo iniziare l’analisi. Per prima cosa studio il comportamento del sistema per ingressi di varia natura: -

-

utilizzo la funzione gradino considerando l’ipotesi di un circuito a corrente continua. Bisogna denotare che, qualora utilizzassimo una funzione gradino, dovremmo determinare gli intervalli di massimo e di minimo affinché potremmo lavorare su tensioni consone alla circuito in questione. utilizzo la funzione cosinusoidale considerando l’ipotesi di un circuito a corrente alternata. Ovviamente in questa ipotesi potremmo lavorare: o sull’ampiezza della cosinusoide: 𝑣 𝑡 = 𝑨𝑐𝑜𝑠(𝜔𝑡 + 𝜑) o sulla frequenza della cosinusoide: 𝑣 𝑡 = 𝐴𝑐𝑜𝑠(𝝎𝑡 + 𝜑) o sulla fase della cosinusoide: 𝑣 𝑡 = 𝐴𝑐𝑜𝑠(𝜔𝑡 + 𝝋)

In un circuito elettrico l’uscita a regime non è mai immediata a causa di componenti derivativi ed integrativi (induttanza e capacità), si instaura dunque un periodo in cui le uscite saranno diverse dalla normale, motivo per il quale un certo intervallo non viene considerato. Considero la costante di tempo 𝜏 = 𝑅𝐶 (considerando un circuito RC) e considero le equazioni: 𝑡! 20% → 80% = 1,4 𝜏 𝑡! 80% → 90% = 2,2 𝜏 Rispettivamente il tempo di salita del condensatore. Altra informazione utile è la frequenza di taglio, che risulta essere pari a: 𝑓! =

1 2𝜋𝜏

La frequenza di taglio è un parametro di definizione delle proprietà dei filtri elettrici passa basso e passa alto. Spesso è anche definita come la frequenza alla quale il rapporto fra l'ampiezza del segnale di uscita e quello di ingresso (o attenuazione) vale ~ 0.707, ovvero 1/ 2. (wikipedia) Costanti tempo e frequenze di taglio considerate standard Per la pre-emphasis/de-emphasis dei filtri RC: Costante temporale Organizzazione Frequenza di taglio fc in Hz in µs RIAA 7958 20 RIAA, NAB 3183 50 1592 100 RIAA 318 500 200 796 140 1137 MC 120 1326 100 1592 MC 90 1768 RIAA 75 2122 FM 70 2274 NAB, PCM 50 3183 DIN 35 4547 25 6366 AES 17.5 9095 PCM 15 10610

Non conoscendo la quantità dei contributi capacitivi ed ipotizzando una struttura non eccessivamente complessa, avremo un tempo da considerare pari a 87,538ms (5 capacità RIAA al tempo di carica fino al 90%). Invece la costante di tempo di un induttore è: 𝜏=

! !" 𝐿 → 𝑣 𝑡 = 𝑉! 𝑒 !! = 𝑉! 𝑒 ! ! 𝑅

il quale scarica la corrente accumulata nel circuito nel tempo. Il valore della corrente e della tensione ai capi di un induttore è:

essendo la linea verde la tensione in ingresso, la linea rossa la corrente e la linea blu la tensione. In commercio esistono numerosi induttori aventi valori diversi, considerando che la media oscilla tra i 0,18𝜇𝐻 e i 100𝜇𝐻 avremo, prendendo come riferimento quest’ultimo valore, un ritardo nell’ordine di 0,5ms da dividere per una resistenza, che se la consideriamo unitaria, il valore massimo per una complessità di 5 induttanze totali è di 2,5ms. Questo valore è da aggiungere al ritardo avuto dai condensatori, avremo un totale di 90,038ms. Generalizzando il ritardo da considerare prima di effettuare un corretto campionamento è pari a: 𝑟 = 17,5076 ∗ 𝑛𝑢𝑚𝑒𝑟𝑜!"#$%#&'(")* + 0,5 ∗ 𝑛𝑢𝑚𝑒𝑟𝑜!"#$%%&'!

Considerando la stessa complessità 𝛾 per L e C (ossia lo stesso numero di componenti) avremo che: 𝑟 = 18,0056 ∗ 𝛾 (𝑚𝑠) Vale a dire che considereremo i dati dal valore 𝑟 e non da 0. Ovviamente un margine maggiore comporta una maggior sicurezza ed una maggior affidabilità (e veridicità) del modello in quanto avremo la risposta a regime dell’intero sistema senza eventuali ritardi di corrente o tensione (oppure è possibile eliminare i primi n elementi che vengono considerati non buoni, ad esempio il primo secondo). Fase ed ampiezza Le ipotesi e l’analisi della fase e dell’ampiezza hanno caratteristiche diverse rispetto a quelle in frequenza (le quali vedremo poi). Le componenti circuitali RLC hanno le seguenti caratteristiche di ampiezza/fase:

Ma ci sono altri componenti da considerare: -

Diodi ed applicazioni (rettificatori a semplice e a doppia semi-onda, ponte di diodi) Diodo Zener (direttamente ed inversamente polarizzati) Mosfet BJT Amplificatori operazionali Inverter

Ipotesi comportamentali a semplici funzioni Oltre alla componente cosinusoidale possiamo avere altre funzioni da poter inserire all’ingresso della nostra blackbox. Ipotizziamo di avere una funzione gradino, una tensione del tipo: 𝑣 𝑡 =𝑘 Le informazioni ricavabili da una funzione di questo tipo sono pressoché minime in quanto restituiscono in uscita un valore della tensione costante aumentata o diminuita di un determinato fattore (non viene considerata ne la fase che la frequenza). La consideriamo una sinusoide a frequenza nulla, fase pari a 90° ed ampiezza k. Risulta più utile invece la funzione rampa, in cui la tensione è funzione del tempo nella forma: 𝑣 𝑡 =𝑡 Questa funzione ci permette di stabilire l’intervallo di funzionamento del circuito, il quale necessita a volte delle tensioni di soglia minime al funzionamento dei vari componenti (ad esempio un diodo). White Noise, Test di Anderson Il rumore bianco è una funzione aperiodica ad ampiezza costante caratterizzata dalla distribuzione uniforme delle frequenze. E’ un ottimo strumento per la determinazione della banda di funzionamento del sistema. Viene definito un vettore casuale 𝑤 un rumore bianco se e solo se il suo valor medio è nullo e la sua matrice di autocorrelazione è multiplo della matrice identità per un valore pari alla varianza del segnale stesso: 𝜇! = 𝐸 𝑤 = 0 𝑅!! = 𝐸 𝑤𝑤 ! = 𝜎 ! 𝐼 Un processo continuo aleatorio 𝑤 𝑡 , 𝑡 ∈ ℝ è bianco se e solo se vengono rispettate le seguenti condizioni: 𝜇! = 𝐸 𝑤(𝑡) = 0 𝑅!! = 𝐸 𝑤 𝑡! 𝑤(𝑡! ) = 𝜎 ! 𝛿(𝑡! − 𝑡! )

Si ha un valor medio nullo ed una potenza infinita al tempo zero poiché la funzione autocorrelatrice è la funzione delta di Dirac. La trasformata di Fourier della funzione di autocorrelazione è costante pari a: 𝐹 𝑅!! = 𝑆!! = 𝜎 ! La figura mostra il comportamento del rumore bianco:

Una possibile configurazione elettronica del rumore bianco è:

limitato ovviamente alle frequenza fisiche delle componenti circuitali.

Un’altra configurazione è quella simulata tramite software. Per far fronte a ciò utilizzeremo MATLAB. % SEZIONE GENERAZIONE SEGNALE RANDOM % SEZIONE SETTAGGIO LIMITI DI AMPIEZZA DEL SEGNALE RANDOM n_campioni=1000; %definizione del numero di campioni del rumore che si desidera generare n=rand(n_campioni,1); %genera vettore di numeri random compresi tra 0 e 1%n=randn(10:1); ! %genera vettore di numeri random con distribuzione normale (gaussiana) non compresi tra 0 e 1 maxIN=1; %valore massimo del segnale generato da rand o randn(sempre 1) !minIN=0; %valore minimo del segnale generato da rand o randn(sempre 0)! maxOUT=+5; %valore massimo di ampiezza del segnale casuale desiderato (max. tensione di OUTPUT della scheda di acquisiz.)! minOUT=-5; %valore minimo di ampiezza del segnale casuale desiderato (min. tensione di OUTPUT della scheda di acquisiz.) syms t; %imposto la lettera t come variabile % GENERAZIONE SEGNALE RANDOM SENZA DEFINIRE TEMPO DI CAMPIONAMENTO (ASSE DEI TEMPI)% ! verrà" definito successivamente in fase di lettura dei dati dalla scheda di acquisizione n1=minOUT+(maxOUT-minOUT)*n; %trasla il valore dal valore compreso tra 0 e 1 al valore compreso tra minOUT e maxOUT! n2=detrend(n1); %elimina il valore medio del segnale Tc=0.01; % definizione tempo di campionamento in secondi = 100Hz !tmin=0; % definizione valore minimo della finestra temporale tmax=n_campioni*Tc-Tc; % definizione valore minimo della finestra temporale t=[tmin:Tc:tmax]'; !t_n2=[t,n2]; subplot(2,1,1) plot(t,n1) ylabel('n1') grid subplot(2,1,2) plot(t,n2) ylabel('n2')

Il rumore bianco può essere testato per vedere la banda di funzionamento di un circuito elettrico oltre che l’identificazione del modello. Per verificare la bianchezza del rumore bianco adotteremo il test di Anderson, il quale consiste nei seguenti passaggi: -

si fissa un valore 𝛼 compreso nello spazio probabilità (tra 0 e 1) calcolo un valore ±𝛽 tale per cui l’area sottesa alla Gaussiana sia pari ad 1 − 𝛼 : ! ! per 𝛽; +∞ → 𝐴𝑟𝑒𝑎 = per −∞; −𝛽) → 𝐴𝑟𝑒𝑎 = !

-

!

considero la funzione di covarianza campionaria (indico con n la finestra di campionamento): !!!

1 𝛾 𝜏 = 𝑛

𝜀 𝑡 𝜀 𝑡+𝜏 !!!

-

nel test di Anderson verifico la bianchezza tramite la covarianza campionaria normalizzata: 𝜌 𝜏 =

-

Se il rumore è bianco è approssimabile con una gaussiana a media nulla e varianza unitaria ridotto di un fattore dipendente dalla dimensione della finestra campionata 𝜌 𝜏 =

-

𝛾 𝜏 𝛾 0

𝐺 0,1 𝑛

Asintoticamente i valori di 𝜌(𝑖) e 𝜌(𝑗) per 𝑖 ≠ 𝑗 sono variabili scorrelate

-

Considerate 𝑚 valutazioni di 𝜌(𝜏) si valuta il numero di campioni di 𝜌(𝜏) contenuti fuori dall’intervallo normalizzato della Gaussiana, ossia il numero di campioni per cui: 𝜌 𝜏 ∉ −

-

𝛽 𝑛

;+

𝛽 𝑛

se il rumore è bianco allora: 𝑛 5){ //salta le prime 5 righe temp = readline(1); //legge riga per riga il file 1 tempArray = split(temp, “\t”); //divido la riga in un vettore number = split[tempArray[2], “E”]; // prendo il valore della terza colonna e divido a sua volta value = number[0]*(10 ^ number[1]); printfile(2,replace(value, “,”,”.”)); number = split[tempArray[5], “E”]; // prendo il valore della terza colonna e divido a sua volta value = number[0]*(10 ^ number[1]); printfile(3,replace(value, “,”,”.”)); } i++; }

In MATLAB inseriamo i valori nei due vettori (input e output) ed eseguiamo il System Identification Tools. Importiamo i dati indicando i vettori di ingresso e uscita ricordandoci di avere un sample time consono alla frequenza di campionamento altrimenti possiamo non soddisfare la condizione di Nyquist utili al corretto campionamento e alla successiva corretta ricostruzione (v. dispense TLC Ennio Gambi). Sono stati stimati i parametri di identificazione tramite dei modelli polinomiali ARX e ARMAX:

Andiamo ad analizzare i modelli ARX e ARMAX utilizzando le varie funzionalità offerte da MATLAB: -

Model output Model resids Transient resp Frequency resp Zero and poles Noise spectrum

MODEL RESIDS - MODELLI ARX

MODEL RESIDS - MODELLI ARMAX

MODEL OUTPUT / POLI E ZERI- MODELLI ARX

MODEL OUTPUT / POLI E ZERI - MODELLI ARMAX

TRANSIENT E FREQUECY STEP / POWER SPECTRUM– MODELLO ARX

TRANSIENT E FREQUECY STEP / POWER SPECTRUM– MODELLO ARMAX

Analizzando i modelli scelti si può notare che: -

-

i modelli ARX 111 e ARMAX 1111 sono i modelli da scartare a priori, hanno un’autocorrelazione ed una crosscorrelazione sicuramente peggiore rispetto a tutti gli altri modelli ed inoltre non rientrano nei limiti utili per la validazione del modello. I modelli ARX diversi da 111 e i modelli ARMAX diversi da 1111 hanno dei best fit pressoché simili, come simili sono le altre caratteristiche, si scartano tutti quei modelli che: o Hanno dei poli e zeri che si annullino a vicenda (il che implica che esistono modelli a complessità minore con ugual risultato) o Fuoriescono maggiormente rispetto ad altri modelli nei limiti utili per la validazione o Hanno una complessità troppo elevata considerando che esistono modelli di complessità minore che mi garantiscono prestazioni pressoché simili se non maggiori rispetto al modello scelto

Tra tutti i vari modelli analizzati abbiamo scelto il modello ARX 211 o il modello ARMAX 2111. Infatti considerato anche il grafico a pag.14 e con le ipotesi comportamentali di circuiti puramente R, L o C, deduciamo che la fase dell’uscita è uguale a quella d’entrata ma esistono componenti che mi creano un picco di tensione alto in uscita quando in entrata abbiamo una tensione negativa (caratteristiche proprie di una coppia condensatoreinduttanza in cui la fase si annulla). La nostra ipotesi comportava una complessità non molto elevata la quale è riscontrata nei modelli scelti. Tra il modello ARX e il modello ARMAX abbiamo scelto il secondo in quanto: -

-

lo sfasamento dovuto a componenti induttivi e capacitivi, anche se la loro differenza di fase tende a 0 mi genera un ritardo che devo tener in considerazione quindi il mio WN utilizzato per l’identificazione è considerato, oltre che nell’istante t, anche in un certo numero di istanti precedenti (risultando il numero di istanti precedenti pari a 1) L’ARX ha componente 𝑛𝑎 = 2, 𝑛𝑏 = 1, 𝑛𝑘 = 1 mentre nell’ARMAX 𝑛𝑎 = 2, 𝑛𝑏 = 1, 𝑛𝑐 = 1, 𝑛𝑘 = 1. Si può notare infatti che i coefficienti di na e nb hanno lo stesso grado, l’unica differenza è la parte media mobile dell’ARMAX

Da considerare anche il comportamento della fase e dell’ampiezza dei vari modelli, infatti se non consideriamo i modelli ARX111 o ARMAX1111 si evince che la posizione finale ed iniziale della fase tendono ad allinearsi in un unico punto:

il che significa che la somma delle fasi dei vari poli/zeri hanno un unico risultato: −𝜋. L’ampiezza invece figura risultati diversi in quanto il suo valore agli estremi è dato da: 𝐹 𝑗𝜔 𝐹 𝑗𝜔

!"#$%&

!"!#!$%&

= 20 log 𝐾 + 𝑟20log (𝜔)

= 20 log 𝐾 + 𝑟20 log 𝜔 +

𝑚!

𝑧! 20 log 𝜔 − 𝑧!

𝑛!

𝑝! 20 log 𝜔 𝑝!

Essendo K la costante di guadagno e 𝑧! , 𝑝! gli zeri e i poli della funzione di trasferimento del tipo: 𝐹 𝑠 = 𝐾𝑠 !

1 + 𝑧! 𝑠 1 + 𝑝! 𝑠

! !

semplificata a fattori esclusivamente binomiali e monomiali.

Ricordando che i valori della fase possono darmi un contributo iniziale/finale pari a: 𝜑!"!#!$%& =

𝜑!"#$%& =

𝜋 2

𝐾 𝑟𝜋 𝜋 + −1 + 2 2 𝐾

𝜋 2

𝐾 𝑟𝜋 𝜋 = −1 + 2 𝐾 2

𝑚!

𝑧! 𝜋 − 𝑧! 2

𝑛!

𝑝! 𝜋 = 𝑝! 2

𝐾 −1 +𝑟 𝐾 𝐾 −1 +𝑟+ 𝐾

𝑚!

𝑧! − 𝑧!

𝑛!

𝑝! 𝑝!

Essendo la funzione: 𝑥 𝑥 = −1 𝑠𝑒 𝑥 < 0 = 1 𝑠𝑒 𝑥 > 0, 𝑥 𝑥 esprimibile anche come funzione sign(x), letteralmente segno di x. Si nota che la fase iniziale dipende da K (il guadagno) e da r (la molteplicità del termine monomio), inoltre essa risulta costante per ogni combinazione di produttorie di termini binomici. Si evince che a parità di fase finale devono esistere termini che tendono ad annullarmi eventuali differenze di fase. Comunque sia ogni modello con questa fase a cui tutti tendono può essere ritenuto valido al fine della determinazione del modello. L’ampiezza ha fase finale diversa da caso a caso in quanto il suo valore finale è dato dalla somma di un certo numero di contributi. Altro fattore che teniamo in considerazione è l’ellisse di confidenza dei poli e zeri. Dalle dispense del docente Lucio Demeio (allegato in PDF) ricordiamo che l’ellisse ha equazione (considerando una bernoulliana): 𝑧!! 𝐹 𝑋! , 𝑝 = 1 + 𝑧!! 𝑝! − 2 𝑋! + !

!

2𝑛

!

𝑝 + 𝑋! = 0

Il quale mi restituisce due valori 𝑝! , 𝑝! (per un certo 𝑋! ).

Applicando lo stesso principio, considerando uno spostamento nel piano di Gauss, si evince che la funzione sia del tipo: 𝑧!! 𝐹 𝑥 − 𝑥! , 𝑦 − 𝑦! = 1 + 𝑧!! !

𝑦 − 𝑦!

!

−2

𝑥 − 𝑥! +

!

2𝑛

(𝑦 − 𝑦! ) + 𝑥 − 𝑥!

!

=0

Ossia un ellisse di confidenza centrato in 𝑃(𝑥! , 𝑦! ) essendo 𝑥! la parte reale e 𝑦! la parte immaginaria. La generazione di un contributo a forma di “ellisse” è dato da componenti aventi una parte reale ed una immaginaria, quindi saremo nella condizione di avere termini di secondo grado con soluzioni complesse e coniugate, ad es...