Corso di Identificazione dei Modelli e Analisi dei Dati (AA 2009...
Transcript of Corso di Identificazione dei Modelli e Analisi dei Dati (AA 2009...
Corso di Identificazione dei Modelli e Analisi dei Dati (AA 2011-2012)
Prof. Sergio Bittanti
I esercitazione di laboratorio 21/11/2011
Marcello Farina
Sommario • Serie temporali e sistemi
– Generazione di sistemi e serie temporali – Analisi di sistemi – Analisi di serie temporali (“identificazione non parametrica”)
• Predittori
• Identificazione
– Identificazione parametrica – Selezione della struttura – Validazione dei modelli identificati
• Esercizi
Serie temporali e sistemi
• Generazione di sistemi a tempo discreto – LTI format (rappresentazione esterna)
>> system = tf(num,den,Ts)
– num: numeratore
– den: denominatore
– Ts: tempo di campionamento (-1 default)
>> system=tf([1 1],[1 -0.8],1)
Transfer function:
z + 1
-------
z - 0.8
Serie temporali e sistemi
• Generazione di segnali
– Risposta allo scalino >> T=[0:Ts:Tfinal];
>> [y,T]=step(sistema,Tfinal); % step response
>> step(sistema,Tfinal) % plot
– Risposta a un ingresso qualsiasi >> T=[0:Ts:Tfinal];
>> u=sin(0.1*T); % ingresso sinusoidale
>> [y,T]=lsim(sistema,u,T);
Serie temporali e sistemi – Rumore bianco
>> N=100; % numero di campioni
>> lambda=1; % deviazione standard
>> mu=1; % media
• I metodo >> u=mu+lambda*randn(N,1);
• II metodo >> u=wgn(M,N,P); % segnale a media nulla
M: numero di segnali
P: potenza del segnale (dBW) (10log10(lambda^2))
• Somma di segnali e WGN >> u=awgn(X,SNR,SIGPOWER);
X: segnale di partenza
SRN: signal-to-noise ratio (dB)
SIGPOWER: potenza del segnale X (dBW) (20log10(.))
- se SIGPOWER=‘measured’, AWGN misura la potenza el
segnale X.
Serie temporali e sistemi
– Risposta di un sistema (lineare) a segnali qualsiasi >> u=1+randn(3000,1);
>> [y,T]=lsim(sistema,u,[1:3000]);
– Il toolbox di identificazione si serve di istruzioni differenti per definire un modello (detto IDPOLY format) >> modello = idpoly(A,B,C,D,F,NoiseVariance,Ts)
A,B,C,D,F sono vettori (contengono coefficienti di
polinomi in ordine crescente di esponente di z-1)
Modello:
A(q) y(t) = [B(q)/F(q)] u(t-nk) + [C(q)/D(q)] e(t)
(Se ARMAX: F=[], D=[])
>> modello = idpoly(sistema)
Sistema è un sistema LTI (definito tramite ss, tf, etc.)
Serie temporali e sistemi
– Risposta di un sistema (lineare, in IDPOLY format) a segnali qualsiasi
>> y=sim(model,u);
– I dati possono essere forniti come IDDATA objects: >> y_id=sim(y,u,Ts);
u: input data
y: output data
Ts: sampling time
Serie temporali e sistemi
• Analisi del sistema (LTI e IDPOLY format):
– Zeri e poli: >> [Z,P,K]=zpkdata(sistema)
Z = [-1] % posizione degli zeri
P = [0.8000] % posizione dei poli
K = 1 % guadagno (non è pari a G(1)!)
– Diagramma di Bode: >> [A,phi] = bode(sistema,W)
W: valori della frequenza ‘w’ in cui calcolare ‘A(w)’ e
‘phi(w)’
A: amplificazione (A(w))
Phi: sfasamento (phi(w))
Serie temporali e sistemi
– Diagramma di Bode:
>> [A,phi] = bode(sistema,[0.1 1])
A(:,:,1) =9.1179
A(:,:,2) =1.9931
phi(:,:,1) =-24.2456
phi(:,:,2) =-78.5036
>> bode(sistema)
>> % disegna il diagramma
(modulo –dB- e fase)
Serie temporali e sistemi
• Analisi di serie temporali - m=mean(X) % Average or mean value.
– Y=detrend(X,’constant’) % Remove constant trend from data sets.
– Y=detrend(X,’linear’) % Remove linear trend from data sets.
– Funzione di covarianza (¿) >> gamma=covf(y,10) % tau=0,1,…,10-1
>> plot([0:9],gamma)
– Spettro ¡(!) – metodo del periodogramma >> periodogram(y)
>> [Pxx,w]=periodogram(y);
>> plot(w/pi,10*log10(Pxx))
° Attenzione: covf calcola la funzione di correlazione! BISOGNA PRIMA DEPOLARIZZARE IL SEGNALE
Predittori
• Dato il sistema
il predittore a k passi si calcola come
dove Ek(z) e Rk(z) sono quoziente e resto della lunga divisione, e
(si calcolano “a mano”)
y(t) =B(z)
A(z)u(t¡ ¿) +
C(z)
A(z)e(t), e(t) =WN(0; ¸2)
y(t=t¡ k) =B(z)Ek(z)
C(z)u(t¡ ¿) +
Rk(z)
C(z)y(t)
=B(z)Ek(z)
C(z)u(t¡ ¿) +
~Rk(z)
C(z)y(t¡ k)
R(z) = ~R(z)z¡k
Predittori
- Predittore
>> Yp = predict(model,data,k) % predittore a k passi.
model: IDPOLY format
data: IDDATA format
k: passi di predizione
– Errore di predizione
>> Epsilon = pe(model,data) % errore di predizione a k passi
– TEST DI ARDERSON PER L’ERRORE DI PREDIZIONE (regione di
confidenza del 99%)
>> resid(model,data)
Predittori
• Esempio:
>> modello=idpoly([1 0.5],[0 1 0.3],[1 0.6],[],[],1,1)
Discrete-time IDPOLY model: A(q)y(t) = B(q)u(t) + C(q)e(t)
A(q) = 1 + 0.5 q^-1
B(q) = q^-1 + 0.3 q^-2
C(q) = 1 + 0.6 q^-1
This model was not estimated from data.
Sampling interval: 1
>> dati=randn(100,2);
>> resid(modello,dati)
Identificazione Altre istruzioni Matlab (2 Identification toolbox)
• Parametric model estimation – ar - AR-models of signals using various approaches.
– armax - Prediction error estimate of an ARMAX model.
– arx - LS-estimate of ARX-models.
– bj - Prediction error estimate of a Box-Jenkins model.
– init - Initialize (randomize) the parameters of a model.
– ivar - IV-estimates for the AR-part of a scalar time series.
– iv4 - Approximately optimal IV-estimates for ARX-models.
– n4sid - State-space model estimation using a sub-space method.
– nlarx - Prediction error estimate of a nonlinear ARX model.
– nlhw - Prediction error estimate of a Hammerstein-Wiener model.
– oe - Prediction error estimate of an output-error model.
– pem - Prediction error estimate of a general model.
Identificazione • Model structure selection
– aic - Compute Akaike's information criterion.
– fpe - Compute final prediction criterion.
– arxstruc - Loss functions for families of ARX-models.
– selstruc - Select model structures according to
various criteria.
– idss/setstruc - Set the structure matrices for idss
objects.
– struc - Generate typical structure matrices for
ARXSTRUC.
• Per maggiori dettagli si veda:
>> help ident