Identification/RELS_method.m at main · AliAlaei76/Identification · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
%% RELS method
%% Estimate of Theta0 and P0

clear;
clc;
close all;
s=tf('s');
G=-6.2447/((s+0.2423)*(s^2+0.35*s+77.37));
t=0:0.1:20;
N=numel(t);
u1=wgn(N,1,1);
[y,t]=lsim(G,u1,t);
y1=y+0.005*rand(N,1); % Measured output
n=5;
p=16;
u1e=[0;0;0;0;0;u1];
y1e=[0;0;0;0;0;y1];
U=zeros(20,11);
for i=1:20
    U(i,:)=[-y1e(i+4) -y1e(i+3) -y1e(i+2) -y1e(i+1) -y1e(i) u1e(i+5) u1e(i+4) u1e(i+3) u1e(i+2) u1e(i+1) u1e(i)];
end
theta0=(U'*U)^(-1)*U'*y1(1:20);
P0=100*eye(p);
v0=y1(1:20)-U*theta0;
theta0=[theta0;zeros(n,1)];
%% Implemention of RELS method

shift=20;
theta=[theta0 zeros(p,N)];
P=zeros(p,p,N);
P(:,:,1)=P0;
v=[v0;zeros(N-20,1)];
ut=zeros(N,p);
for i=21:N
    ut(i-shift,:)=[-y1(i-1) -y1(i-2) -y1(i-3) -y1(i-4) -y1(i-5) u1(i) u1(i-1) u1(i-2) u1(i-3) u1(i-4) u1(i-5) v(i-1) v(i-2) v(i-3) v(i-4) v(i-5)];
    P(:,:,i-shift+1)=P(:,:,i-shift)-(1/(1+ut(i-shift,:)*P(:,:,i-shift)*ut(i-shift,:)'))*(P(:,:,i-shift)*ut(i-shift,:)'*ut(i-shift,:)*P(:,:,i-shift));
    theta(:,i-shift+1)=theta(:,i-shift)+P(:,:,i-shift+1)*ut(i-shift,:)'*(y1(i)-ut(i-shift,:)*theta(:,i-shift));
    v(i)=y1(i)-ut(i-shift,:)*theta(:,i-shift+1);
    if norm(theta(:,i-shift+1)-theta(:,i-shift))<10^(-6)
        break
    end
end
numberofiteration=i-shift
%% Plots and results

m=0:1:numberofiteration;
plot(m,theta(1,1:numberofiteration+1),m,theta(2,1:numberofiteration+1),m,theta(3,1:numberofiteration+1),m,theta(4,1:numberofiteration+1),m,theta(5,1:numberofiteration+1),'linewidth',1.5)
grid on
xlabel('Iteration')
ylabel('$\hat{\theta}$','Interpreter',"latex",'FontSize',13)
legend('a_1','a_2','a_3','a_4','a_5')

plot(m,theta(6,1:numberofiteration+1),m,theta(7,1:numberofiteration+1),m,theta(8,1:numberofiteration+1),m,theta(9,1:numberofiteration+1),m,theta(10,1:numberofiteration+1),m,theta(11,1:numberofiteration+1),'linewidth',1.5)
grid on
xlabel('Iteration')
ylabel('$\hat{\theta}$','Interpreter',"latex",'FontSize',13)
legend('b_0','b_1','b_2','b_3','b_4','b_5')

plot(m,theta(12,1:numberofiteration+1),m,theta(13,1:numberofiteration+1),m,theta(14,1:numberofiteration+1),m,theta(15,1:numberofiteration+1),m,theta(16,1:numberofiteration+1),'linewidth',1.5)
grid on
xlabel('Iteration')
ylabel('$\hat{\theta}$','Interpreter',"latex",'FontSize',13)
legend('c_1','c_2','c_3','c_4','c_5')

yhat=zeros(i,1);
for j=21:i
    yhat(j)=ut(j-shift,:)*theta(:,j-shift+1);
end
e=y1(21:i)-yhat(21:i);
s=e'*e
sigma=s/(N-p)

plot(t(21:N),y1(21:N),t(21:N),yhat(21:N))
grid on
xlabel('Time(second)')
ylabel('y(t)')
legend('y(t)','yhat(t)')

plot(t(21:N),e)
grid on
xlabel('Time(second)','FontSize',15);
ylabel('$\hat{e}$','Interpreter',"latex",'FontSize',15)