如何对含有u(k)的表达式进行系统辨识。

希望大神们可以教教我。。。
书上的例子都是y(k)=a1*y(k-1)+a2*y(k-2)+b1*u(k-1)+b2*u(k-2)+v(k)这种类似形式的。
如果是y(k)=a1*y(k-1)+a2*y(k-2)+b1*u(k)+b2*u(k-1)+b3*u(k-2)+v(k)这种形式的话。利用matlab进行递推最小二乘辨识和递推极大似然辨识的话程序该怎么编写呢?小白一个,希望大佬教教我图片说明例如上图的差分方程,我该怎么弄呢?
下面是我的代码。。。存在错误,但不懂怎么改正。很懵逼
clear all
close all
%%%%%%%%%产生仿真数据%%%%%%%%%
n = 2;
total = 1000;
sigma = 1;%噪声变量的均方根

%M序列作为输入
y1 = 1; y2 = 1; y3 = 1; y4 = 0;
for i = 1:total;
x1 = xor(y3,y4);
x2 = y1;
x3 = y2;
x4 = y3;
y(i) = y4;
if y(i)>0.5,u(i) = -1;
else u(i) = 1;
end
y1 = x1; y2= x2; y3 = x3; y4 = x4;
end
figure(1);
stem(u),grid on
%%%%%%%%%%系统输出%%%%%%%%%%
y(1)=0;y(2)=0;
v=sigma*randn(total,1);%噪声
y(1)=1;y(2)=0.01;
for k=3:total
y(k) = (-0.1538)*y(k-1)-0.3846*y(k-2)+0.3846*u(k)+0.7692*u(k-1)+0.3846*u(k-2)+1*v(k)+0*v(k-1)+0*v(k-2);
end

%%%%%%%%%初始化%%%%%%%%%
theta0=0.001*ones(7,1);
e1(1)=(-0.1538)-theta0(1);
e2(1)=(-0.3846)-theta0(2);%误差初始化
e3(1)=0.3846-theta0(3); e4(1)=0.7692-theta0(4);
e5(1)=0.3846-theta0(5);
e6(1)=0-theta0(6); e7(1)=0-theta0(7);
a_hat(1)=theta0(1);
a_hat(2)=theta0(2);
b_hat(1)=theta0(3);
b_hat(2)=theta0(4);
b_hat(3)=theta0(5);
c_hat(1)=theta0(6);
c_hat(2)=theta0(7);

P0=eye(7,7);%矩阵P初始化
for i=1:n
yf(i)=0.1;uf(i)=0.1;vf(i)=0.1;
fai0(i,1)=-yf(i);
fai0(n+i,1)=uf(i);
fai0(2*n+i,1)=vf(i);
end
e(1)=1.0;
e(2)=1.0;

%%%%%%递推算法%%%%%%
for i=n+1:total
pusai=[-y(i-1);-y(i-2);u(i);u(i-1);u(i-2);e(i-1);e(i-2)];

C=zeros(n*3,n*3);
Q=zeros(3*n,1);
Q(1)=-y(i-1);
Q(n+1)=u(i-1);
Q(2*n+1)=e(i-1);
for j=1:n
    C(1,j)=-c_hat(j);
    c(n+1,n+j)=-c_hat(j);
    c(2*n+1,2*n+j)=-c_hat(j);
    if j>1
        C(j,j-1)=1.0;
        C(n+j,n+j-1)=1.0;
        C(2*n+j,2*n+j-1)=1.0;
    end
end
fai=C*fai0+Q;
K=P0*fai*inv(fai'*P0*fai+1);
P=[eye(7,7)-K*fai']*P0;

e(i)=y(i)-pusai'*theta0;
theta=theta0+K*e(i);

P0=P;
theta0=theta;
fai0=fai;

a_hat(1)=theta(1);
a_hat(2)=theta(2);
b_hat(1)=theta(3);
b_hat(2)=theta(4);
b_hat(3)=theta(5);
c_hat(1)=theta(6);
c_hat(2)=theta(7);
e1(i)=(-0.1538)-a_hat(1); e2(i)=(-0.3846)-a_hat(2);
e3(i)=0.3846-b_hat(1); e4(i)=0.7692-b_hat(2);e5(i)=b_hat(3);
e6(i)=0-c_hat(1);e7(i)=0-c_hat(2);
end

figure(2)
plot(e1);
hold on
plot(e2);
hold on
plot(e3);
hold on
plot(e4);
hold on
plot(e5);
hold on
plot(e6);
hold on
plot(e7);
title('Parameter Estimation Error');
xlabel('times');
ylabel('error');
hold off
figure(3)
plot(e);
title('Output Error');
xlabel('times');
ylabel('error');

Csdn user default icon
上传中...
上传图片
插入图片
抄袭、复制答案,以达到刷声望分或其他目的的行为,在CSDN问答是严格禁止的,一经发现立刻封号。是时候展现真正的技术了!
立即提问