MODELE REGRESYJNE
Model AR
Model AR (Auto Regresive) zdefiniowany jest jako:
n
y(t)= - y(t - i)+ e(t) (1)
"ai
i=1
Gdzie:
" y(t) jest wyjściem układu,
" ai - parametrami modelu.
" n  jest rzędem modelu
" y(t-1) jest wyjściem modelu z chwili poprzedniej
" e(t)  jest białym szumem.
Inaczej równanie 1, można zapisać
A(q)y(t)= e(t), A(q)= 1+ a1q-1 +K+ anq-n
Gdzie oraz to operator przesunięcia
wstecz.
Błąd estymacji parametrów modelu zdefiniowany jest jako
µ(t)= y(t)- w
(t)
Gdzie y(t) jest sygnałem wejściowym do procesu estymacji a ( ) jest sygnałem wyestymowanymi.
Model ARMA
Model ARMA (AutoRegresiveMovingAverage) układu jest rozszerzeniem modelu AR o tzw. Średnia
ruchomą i może być zapisany jako:
y(t)+ a1y(t -1)+K+ an y(t - n)= e(t)+ c1e(t -1)+K+ cn y(t - nc)
c
lub
A(q)y(t)= C(q)e(t), A(q)= 1+ a1q-1 +K+ anq-n
c
C(q)= 1+ c1q-1 +K+ cncq-n
Gdzie A(q) i C(q) SA operatorami przesunięcia wstecz
Modele ARX
Model ARX (AutoRegressive with eXogenous input) jest zdefiniowany jako model autoregresyjny z
zewnętrznym wejściem:
Gdzie
" y(t) jest wyjściem układu,
" a1,& ,an b1, & bn- parametrami modelu.
" na  jest liczbą biegunów modelu
" nb  jest liczba zer modelu +1
" nk  jest opóz
nieniem ()
" y(t-1), & y(t-na) jest wyjściem modelu z chwili poprzedniej
" u(t- nk), & u(t- nk  nb +1) jest wejściem modelu z chwili poprzedniej
" e(t)  jest białym szumem.
Inaczej model ten można zapisać:
Gdzie:
a q-1 jest operatorem przesunięcia wstecz.
Analogicznie zdefiniowany jest model typu ARMAX:
Lub
Gdzie:
" y(t) jest wyjściem układu,
" a1,& ,an b1, & bn- parametrami modelu.
" na  jest liczbą biegunów modelu
" nb  jest liczba zer modelu +1
" nk  jest opóz
nieniem ()
" y(t-1), & y(t-na) jest wyjściem modelu z chwili poprzedniej
" u(t- nk), & u(t- nk  nb +1) jest wejściem modelu z chwili poprzedniej
" e(t- nk), & e(t- nk  nb +1) jest białym szumem.
a A, B i C są wyrażone jako funkcja q
Przykład
Model AR
Załóżmy sygnał odpowiedzi układu o jednym stopniu swobody. Korzystając z modelu AR dokonajmy
estymacji modelu układu, korzystając z pakietu MATLAB. Przykładowy kod programu poniżej.
clear
close all
t = 0:0.01:10;
y = exp(-0.02*2*pi*2*t).*sin(2*pi*2*t);
[mod] = ar(y,2);
Y_es = predict(mod,y');
figure
plot(t,y,'LineWidth',2)
hold on
plot(t,Y_es{1},'--r','LineWidth',2)
xlabel('Czas [s]','FontSize',16)
ylabel('Amplituda ','FontSize',16)
grid on
legend('Analizowany sygnał','Sygnał z wyestymowanego modelu')
efekt działania algorytmu przedstawia rysunek poniżej
Dodatkowo, wprowadz
my zakłócenie do sygnału
y_nois = exp(-0.02*2*pi*2*t).*sin(2*pi*2*t)+rand(1,length(t))*0.25;
Przeprowadz
my ponownie proces estymacji modelu i sprawdz
my wynik (rysunek poniżej)
Model ARMA/ARMAX
Załóżmy układ o jednym stopniu swobody zdefiniowany wg wskazówek z poprzednich zajęć.
Przykładowy kod poniżej:
m=1;
k=3000;
c=2;
A = [0 1;-k/m -c/m];
B = [-1/m 0]';
C = [0 1];
D = 0;
sys = ss(A,B,C,D);
aby wyznaczyć parametry modalne modelu użyjemy funkcji damp
[Wn,Z] = damp(sys);
Gdzie Wn to częstości własne układu, a Z  współczynniki tłumienia.
Przeprowadz
my symulacje, podając na wejście układu wymuszenie w postaci białego szumu.
t = 0:0.01:50;
u = 0.5*rand(1,length(t));
Y_out = lsim(sys,u,t);
Wynik syulacji przedstawia rysunek poniżej
Dokonajmy estymacji parametrów modelu I porównajmy parametry wyznaczone z
wyestymowanymi:
model = armax([Y_out],[2 2]);
model_z = zpk(model)
[Wn_es, Z_es P] = damp(model_z);
[Wn(2)/2/pi Wn_es(1)/2/pi*100;Z(1) Z_es(1)]
ans =
Częstości wyznaczone [Hz 8.7173 7.6344 Częstości wyestymowane [Hz]
Wsp. tłumienia wyestymowane 0.0183 0.0912 Wsp. tłumienia wyestymowane
Jak można zauważyć model nie został poprawnie wyznaczony. Zwiększmy wiec rząd modelu
model = armax([Y_out],[4 4]);
damp(model_z);
Wybieramy rozwiÄ…zanie 3 i 4 (odpowiedz dlaczego)
[Wn_es, Z_es P] = damp(model_z);
[Wn(1)/2/pi Wn_es(3)/2/pi*100;Z(1) Z_es(3)]
Częstości wyznaczone [Hz 8.7173 8.7478 Częstości wyestymowane [Hz]
Wsp. tłumienia wyestymowane 0.0183 0.0199 Wsp. tłumienia wyestymowane
Sformułuj wnioski na temat powyższego zadania.
Rekursywne metody identyfikacji.
Dla układów niestacjonarnych bardzo często wykorzystuje się techniki oparte o modele regresyjne
działające w sposób rekursywny. Przykładem takiego podejścia może być model RARMAX (Recursive
ARMAX). Schemat identyfikacji układu z wykorzystaniem takiego algorytmu można sformuować
następująco:
Krok 1: Pobranie bieżącej odpowiedzi systemu y(i) z przetwornika A/C.
Krok 2: Obliczenie estymatora bÅ‚Ä™du predykcji µ (i) :
T
Ć
µ (t) = y(i) -¸ (i -1)Õ(i)
Krok 3: Utworzenie pomocniczych wektorów f i g
T
f = U (i -1)Õ(i)
g = D(i -1) f
oraz zapamiętanie współczynnika zapominania 0 =  .
Dla j = 1,& , n¸ , kroki 4 ÷ 5zostajÄ… powtórzone.
Krok 6: Obliczenie poszczególnych elementów macierzy D(i) :
 =  + f g ,
j j-1 j j
D (i) = ( )-1 D(i -1) ,
j j j-1
v = g .
j j
gdzie: j jest skalarem, f , g , Dj i v są elementami poszczególnych wektorów.
j j j
Krok 4: Obliczenie poszczególnych elementów macierzy U (i) dla k = 1, & , j  1:
jeśli j > 1
vk f
j
Ukj (i) = Ukj (i -1) -

j-1
vk = vk + Ukj (i -1)v
j
gdzie: Ukj (i) oznacza odpowiedni element macierzyU (i) .
Krok 5: Uaktualnienie wektora wzmocnień:
v
L(i) =
d
Krok 6: Aktualizacja wektora estymowanych parametrów modelu ¸Ć(i) :
Krok 7: Oczekiwanie na następną próbkę i powrót do Kroku 1.
przy czym:
¸ (i) - wektor parametrów modelu o wymiarze n¸ ×1,
Õ(i) - wektor obserwacji wyjÅ›ciowo-wejÅ›ciowej (ang. regressor vector) o wymiarze n¸ ×1,
Ważnym elementem jest w tym przypadku współczynnik zapominania decydujący o szybkości reakcji
algorytmu na zmienność parametrów oraz jakości wyestymowanego modelu. Współczynnik ten
dobiera się w zależności od dynamiki procesu. Praktyczny przedział stosowalności współczynnika
zapominania to 0.9 do 1. Dla wartości 1 algorytm przyjmuje postać modelu ARMA (ARMAX).
Przykład
Załóżmy sygnał harmoniczny o zmiennej częstotliwości (przykładowy kod poniżej).
t = 0:0.01:40;
frq = [2*ones(1,1001) 3*ones(1,1000) 4*ones(1,1000) 3*ones(1,1000)];
for i = 1:1:length(t)
y(i) = sin(2*pi*t(i)*frq(i));
end
Dodatkowo dodajmy zakłócenie.
y = y+rand(1,length(y))*0.25;
Przebieg sygnału oraz przebieg zmian częstotliwości w czasie pokazano poniżej.
Przeprowadz
my proces identyfikacji
na = 2;
nb = 2;
theta = rarmax(y',[na nb],'ff',0.97);
i wyznaczmy parametry modalne (częstotliwość)
for i=1:length(theta)
aa = [1 theta(i,1:na)];
[wnout_v,z,r]=damp(aa, 1/100);
OMEGA(i,:) = [wnout_v(1:2:na)'/(2*pi)];
end
Przebieg czasowy wyestymowanego sygnału pokazano poniżej.
Zbadajmy wpływ wartości współczynnika zapominania na proces identyfikacji.
theta = rarmax(y',[na nb],'ff',0.99);
wyniki poniżej
Przeprowadz
proces identyfikacji dla różnych współczynników tłumienia. Co zauważasz?
Przykład
Poniżej przeprowadzimy proces identyfikacji układu o jednym stopniu swobody o zmiennych
parametrach. Zmiana parametrów występować może na wskutek np. uszkodzenia układu.
Załóżmy układ o jednym stopniu swobody, o zmiennych parametrach. Uproszczoną symulacje takiego
układu przedstawia kod poniżej:
t = 0:0.01:10;
m=1;
k= [3000 14000 5000 4000];
c=[2 10 20 5];
y = [];U = [];
for i = 1:1:4
A = [0 1;-k(i)/m -c(i)/m];
B = [-1/m 0]';
C = [0 1];
D = 0;
sys = ss(A,B,C,D);
u = 0.5*rand(1,length(t));
y_tmp = lsim(sys,u,t);
y = [y y_tmp'];
U = [U u];
[Wn,Z] = damp(sys);
Om(i) = Wn(1)/2/pi;
Da(i) = Z(1);
end
Proces identyfikacji:
na = 2;
nb = 2;
theta = rarmax([y' U'],[na nb 2 1],'ff',0.95);
for i=1:length(theta)
aa = [1 theta(i,1:na)];
[wnout_v,z,r]=damp(aa, 1/100);
DAMPING(i,:) = [z(1:2:na)'];
OMEGA(i,:) = [wnout_v(1:2:na)'/(2*pi)];
end
t1 = 0:0.01:40.03
figure
subplot(2,1,1)
plot(t1,OMEGA,'LineWidth',2)
hold on
% plot(t,frq,'r','LineWidth',2)
ylim([0 20])
OO = [Om(1)*ones(1,1001) Om(2)*ones(1,1001) Om(3)*ones(1,1001)
Om(4)*ones(1,1001)];
plot(t1,OO,'--r','LineWidth',2)
xlabel('Czas [s]','FontSize',16)
ylabel('Częstotliwość [Hz] ','FontSize',16)
grid on
legend('Parametry wyestymowanego modelu','Zało\
ona zmiana czestotliwość')
subplot(2,1,2)
plot(t1,DAMPING*100,'LineWidth',2)
hold on
DD = [Da(1)*ones(1,1001) Da(2)*ones(1,1001) Da(3)*ones(1,1001)
Da(4)*ones(1,1001)];
plot(t1,DD*100,'--r','LineWidth',2)
grid on
xlabel('Czas [s]','FontSize',16)
ylabel('Wsp. TÅ‚umienia [-] ','FontSize',16)
legend('Parametry wyestymowanego modelu','Zało\
ona zmiana tłumienia')
Wynik
Przeprowadz
proces identyfikacji dla różnych wartości współczynnika zapominania. Co zauważasz?
Zadania do samodzielnego wykonania:
1. Na podstawie znanych metod modelowania przeprowadz
proces identyfikacji układu o
dwóch stopniach swobody (dwumasowy model sprężysto - tłumiący) o stałych parametrach.
Porównaj wartości parametrów założonych oraz wyestymowanych. Dokonaj identyfikacji dla
różnego poziomu zakłóceń sygnału odpowiedzi.
2. Dokonaj identyfikacji parametrów modalnych układu z punktu 1 z uwzględnieniem zmian
parametrów modelu w trakcji symulacji (podobnie jak w ostatnim przykładzie). Analizy
wykonaj dla różnych wartości współczynnika zapominania.