东北大学系统建模作业 12.docx
- 文档编号:3329777
- 上传时间:2022-11-21
- 格式:DOCX
- 页数:29
- 大小:378.70KB
东北大学系统建模作业 12.docx
《东北大学系统建模作业 12.docx》由会员分享,可在线阅读,更多相关《东北大学系统建模作业 12.docx(29页珍藏版)》请在冰豆网上搜索。
东北大学系统建模作业12
系统建模方法作业
姓名:
冯天宇
学号:
1001258
班级:
双控3班
1.考虑如下系统
式中,
为白噪声。
取初值
,
。
分别选择M序列和方差为1的正态分布白噪声作为输入信号
,采用递推最小二乘算法进行参数估计,迭代L=400步停止计算。
要求
i)给出基本迭代公式;
ii)画出程序流程框图;
iii)画出输入输出数据曲线、参数估计曲线、误差曲线;
提示:
产生长度为L方差为1的正态分布白噪声,相应的MATLAB命令为randn(L,1)。
i)迭代公式:
在计算过程中,增益矩阵K(k)可能产生误差,经过算法的迭代,造成误差传递和累积,最后将影响辨识算法的精度。
为此,可采用下面的计算式,以截断误差的传递,保证辨识精度。
ii)程序流程框图:
N
iii)程序、输入输出数据曲线、参数估计曲线、误差曲线、结果分析:
(1)当输入信号u(k)为M序列时,程序:
clear%清理工作间变量
L=400;%M序列的周期
m1=1;m2=1;m3=1;m4=0;%四个移位寄存器的输出初始值
fori=1:
L;%开始循环,长度为L
x1=xor(m3,m4);%第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“异或”
x2=m1;%第二个移位寄存器的输入是第一个移位寄存器的输出
x3=m2;%第三个移位寄存器的输入是第二个移位寄存器的输出
x4=m3;%第四个移位寄存器的输入是第三个移位寄存器的输出
m(i)=m4;%取出第四个移位寄存器的幅值为"0"和"1"的输出信号,即M序列
ifm(i)>0.5,u(i)=-1;%如果M序列的值为"1",辨识的输入信号取“-1”
elseu(i)=1;%如果M序列的值为"0",辨识的输入信号取“1”
end%小循环结束
m1=x1;m2=x2;m3=x3;m4=x4;%为下一次的输入信号做准备
end%大循环结束,产生输入信号u
y(4)=0;y(3)=0;y
(2)=0;y
(1)=0;%设y的前四个初始值为零
fork=5:
400%循环变量从5到400
y(k)=1.5*y(k-1)-0.7*y(k-2)+u(k-3)+0.5*u(k-4);%输出采样信号
end
figure
(1),subplot(2,1,1),stem(u),
title('input')
subplot(2,1,2),stem(y),holdon
k=1:
400;
plot(k,y)
title('output')
%RLS递推最小二乘辨识
c0=[0000]';%直接给出被辨识参数的初始值,即一个充分小的实向量
p0=10^6*eye(4,4);%直接给出初始状态P0,即一个充分大的实数单位矩阵
c=[c0,zeros(4,399)];%被辨识参数矩阵的初始值及大小
e=zeros(4,400);%相对误差的初始值及大小
fork=5:
400;%开始求K
h1=[-y(k-1),-y(k-2),u(k-3),u(k-4)]';
x=h1'*p0*h1+1;
x1=inv(x);%开始求K(k)
k1=p0*h1*x1;%求出K的值
d1=y(k)-h1'*c0;c1=c0+k1*d1;%求被辨识参数c
e1=c1-c0;%求参数当前值与上一次的值的差值
e2=e1./c0;%求参数的相对变化
e(:
k)=e2;%把当前相对变化的列向量加入误差矩阵的最后一列
c0=c1;%新获得的参数作为下一次递推的旧参数
c(:
k)=c1;%把辨识参数c列向量加入辨识参数矩阵的最后一列
p1=p0-k1*k1'*[h1'*p0*h1+1];%求出p(k)的值
p0=p1;%给下次用
end%循环结束
c,e%显示被辨识参数及其误差(收敛)情况
%分离参数
a1=c(1,:
);
a2=c(2,:
);
b1=c(3,:
);
b2=c(4,:
);
ea1=e(1,:
);
ea2=e(2,:
);
eb1=e(3,:
);
eb2=e(4,:
);
figure
(2);%第二个图形
i=1:
400;%横坐标从1到400
plot(i,a1,'r',i,a2,':
',i,b1,'g',i,b2,':
')%画出a1,a2,b1,b2的各次辨识结果
title('ParameterIdentificationwithRecursiveLeastSquaresMethod')%图形标题
figure(3);%第三个图形
i=1:
400;%横坐标从1到400
plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:
')%画出a1,a2,b1,b2的各次辨识结果的收敛情况
title('IdentificationPrecision')%图形标题
输入输出数据曲线:
参数估计曲线:
误差曲线:
结果分析:
表最小二乘递推算法的辨识结果
参数a1a2b1b2
真值-1.50.71.00.5
估计值-1.50.71.00.5
仿真结果表明,由于题意所选步长过大,不能明显看出辨识曲线走势,但从整体看,参数辨识在很短的步长变趋于稳定,由程序结果可得最终a1=-1.5,a2=0.7,b1=1.0,b2=0.5。
(2)当输入信号u(k)为方差为1的正态分布白噪声时,程序:
clear%清理工作间变量
u=randn(400,1);
y(4)=0;y(3)=0;y
(2)=0;y
(1)=0;%设y的前四个初始值为零
fork=5:
400;%循环变量从5到400
y(k)=1.5*y(k-1)-0.7*y(k-2)+u(k-3)+0.5*u(k-4);%输出采样信号
end
figure
(1),subplot(2,1,1),stem(u),
title('方差为1的正态分布白噪声')
subplot(2,1,2),stem(y),holdon
k=1:
400;
plot(k,y)
title('输出')
%RLS递推最小二乘辨识
c0=[0000]';%直接给出被辨识参数的初始值,即一个充分小的实向量
p0=10^6*eye(4,4);%直接给出初始状态P0,即一个充分大的实数单位矩阵
c=[c0,zeros(4,399)];%被辨识参数矩阵的初始值及大小
e=zeros(4,400);%相对误差的初始值及大小
fork=5:
400;%开始求K
h1=[-y(k-1),-y(k-2),u(k-3),u(k-4)]';
x=h1'*p0*h1+1;
x1=inv(x);%开始求K(k)
k1=p0*h1*x1;%求出K的值
d1=y(k)-h1'*c0;c1=c0+k1*d1;%求被辨识参数c
e1=c1-c0;%求参数当前值与上一次的值的差值
e2=e1./c0;%求参数的相对变化
e(:
k)=e2;%把当前相对变化的列向量加入误差矩阵的最后一列
c0=c1;%新获得的参数作为下一次递推的旧参数
c(:
k)=c1;%把辨识参数c列向量加入辨识参数矩阵的最后一列
p1=p0-k1*k1'*[h1'*p0*h1+1];%求出p(k)的值
p0=p1;%给下次用
end%循环结束
c,e%显示被辨识参数及其误差(收敛)情况
%分离参数
a1=c(1,:
);
a2=c(2,:
);
b1=c(3,:
);
b2=c(4,:
);
ea1=e(1,:
);
ea2=e(2,:
);
eb1=e(3,:
);
eb2=e(4,:
);
figure
(2);%第二个图形
i=1:
400;%横坐标从1到400
plot(i,a1,'r',i,a2,':
',i,b1,'g',i,b2,':
')%画出a1,a2,b1,b2的各次辨识结果
title('ParameterIdentificationwithRecursiveLeastSquaresMethod')%图形标题
figure(3);%第三个图形
i=1:
400;%横坐标从1到400
plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:
')%画出a1,a2,b1,b2的各次辨识结果的收敛情况
title('IdentificationPrecision')%图形标题
输入输出数据曲线:
参数估计曲线:
误差曲线:
结果分析:
表最小二乘递推算法的辨识结果
参数a1a2b1b2
真值-1.50.71.00.5
估计值-1.50.71.00.5
由于输入信号不一样,所以输入输出曲线略有不同,但参数估计曲线相差不多,很快便趋于稳定。
2.考虑如下确定性系统
采用梯度校正算法进行参数估计,通过相对误差准则停止计算。
取初值
,并取
,其中
。
要求
i)给出基本迭代公式;
ii)画出程序流程框图;
iii)画出输入输出数据曲线、参数估计曲线
iv)自行设计适合的数据序列(M序列,白噪声等),并解释仿真结果。
i)迭代公式:
ii)程序流程框图:
iii)程序、输入输出数据曲线、参数辨识曲线和误差曲线、三次插值曲线和结果分析:
(1)
当输入信号u(k)为M序列时,程序:
clear%清理工作间变量
L=200;%M序列的周期
y1=1;y2=1;y3=1;y4=0;%四个移位寄存器的输出初始值
fori=1:
L;%开始循环,长度为L
x1=xor(y3,y4);%第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”
x2=y1;%第二个移位寄存器的输入是第一个移位寄存器的输出
x3=y2;%第三个移位寄存器的输入是第二个移位寄存器的输出
x4=y3;%第四个移位寄存器的输入是第三个移位寄存器的输出
m(i)=y4;%取出第四个移位寄存器的幅值为"0"和"1"的输出信号,即M序列
ifm(i)>0.5,u(i)=-1;%如果M序列的值为"1",辨识的输入信号取“-0.03”
elseu(i)=1;%如果M序列的值为"0",辨识的输入信号取“0.03”
end%小循环结束
y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号做准备
end%大循环结束,产生输入信号u
y(4)=0;y(3)=0;y
(2)=0;y
(1)=0;%设y的前四个初始值为零
fork=5:
200;%循环变量从5到200
y(k)=1.5*y(k-1)-0.7*y(k-2)+u(k-3)+0.5*u(k-4);%输出采样信号
end
figure
(1),subplot(2,1,1),stem(u),
title('ijput')
subplot(2,1,2),stem(y),holdon
k=1:
200;
plot(k,y)
title('output')
%给出初始值
h1=[1,1,1,1]';h2=[1,1,1,1]';h3=[1,1,1,1]';g=[0,0,0,0]';I=[1/4,0,0,0;0,1/8,0,0;0,0,1/2,0;0,0,0,1/2];
h=[h1,h2,h3,zeros(4,197)];
%计算样本数据h(k)
fork=4:
198
h(:
k)=[-y(k),-y(k-1),u(k-2),u(k-3)]';
end
E=0.0005;
%计算权矩阵R(k)和g的估计值
fork=1:
198
a=(h(1,k)^2)/4+(h(2,k)^2)/8+(h(3,k)^2)/2+(h(4,k)^2)/2;a1=1/a;R=a1*I;
%计算权矩阵
g(:
k+1)=g(:
k)+R*h(:
k)*(y(k+1)-h(:
k)'*g(:
k));
e1=g(:
k+1)-g(:
k);
e2=e1./g(:
k);
q=0;
ife2<=E,q=1;break;
end
%计算脉冲响应的估计值
end
ifq==1
forj=k:
198
g(:
j)=g(:
k);
end
%绘图
g1=g(1,:
);g2=g(2,:
);g3=g(3,:
);g4=g(4,:
);
figure
(2);
k=1:
198;
subplot(1,2,1);
plot(k,g1,'r',k,g2,'g',k,g3,'b',k,g4,'y'),gridon
title('参数辨识曲线')
else
g1=g(1,:
);g2=g(2,:
);g3=g(3,:
);g4=g(4,:
);
figure
(2);
k=1:
199;
subplot(1,2,1);
plot(k,g1,'r',k,g2,'g',k,g3,'b',k,g4,'y'),gridon
title('参数辨识曲线')
end
%计算模型输出ym及系统输出与模型输出之间的误差Ey
fork=1:
198
ym(k)=h(:
k)'*g(:
k);
Ey(k)=y(k+1)-ym(k);
end
k=1:
198
subplot(1,2,2);
plot(k,Ey),gridon
title('误差曲线')
g,ym,Ey
%显示脉冲响应估计值、模型输出及系统输出与模型输出之间的误差
figure(3);
x=0:
1:
4;
y=[0,g(1,198),g(2,198),g(3,198),g(4,198)];
xi=linspace(0,4);
yi=interp1(x,y,xi,'cubic');
plot(x,y,'o',xi,yi,'m'),gridon
title('三次插值曲线')%画出脉冲响应估计值及其三次插值曲线
输入输出数据曲线:
参数辨识曲线和误差曲线:
三次插值曲线:
结果分析:
表梯度校正算法的辨识结果
参数g1g2g3g4
真值-1.50.71.00.5
估计值-1.49620.70301.05820.5665
(2)当输入信号u(k)为方差为1的正态分布白噪声时,程序:
clear%清理工作间变量
u=randn(200,1);
y(4)=0;y(3)=0;y
(2)=0;y
(1)=0;%设y的前四个初始值为零
fork=5:
200;%循环变量从5到200
y(k)=1.5*y(k-1)-0.7*y(k-2)+u(k-3)+0.5*u(k-4);%输出采样信号
end
figure
(1),subplot(2,1,1),stem(u),
title('方差为1的正态分布白噪声')
subplot(2,1,2),stem(y),holdon
k=1:
200;
plot(k,y)
title('输出')
%给出初始值
h1=[1,1,1,1]';h2=[1,1,1,1]';h3=[1,1,1,1]';g=[0,0,0,0]';I=[1/4,0,0,0;0,1/8,0,0;0,0,1/2,0;0,0,0,1/2];
h=[h1,h2,h3,zeros(4,197)];
%计算样本数据h(k)
fork=4:
198
h(:
k)=[-y(k),-y(k-1),u(k-2),u(k-3)]';
end
E=0.000000000005;
%计算权矩阵R(k)和g的估计值
fork=1:
198
a=(h(1,k)^2)/4+(h(2,k)^2)/8+(h(3,k)^2)/2+(h(4,k)^2)/2;a1=1/a;R=a1*I;
%计算权矩阵
g(:
k+1)=g(:
k)+R*h(:
k)*(y(k+1)-h(:
k)'*g(:
k));
e1=g(:
k+1)-g(:
k);
e2=e1./g(:
k);
q=0;
ife2<=E,q=1;break;
end
%计算脉冲响应的估计值
end
ifq==1
forj=k:
198
g(:
j)=g(:
k);
end
%绘图
g1=g(1,:
);g2=g(2,:
);g3=g(3,:
);g4=g(4,:
);
figure
(2);
k=1:
198;
subplot(1,2,1);
plot(k,g1,'r',k,g2,'g',k,g3,'b',k,g4,'y'),gridon
title('参数辨识曲线')
else
g1=g(1,:
);g2=g(2,:
);g3=g(3,:
);g4=g(4,:
);
figure
(2);
k=1:
199;
subplot(1,2,1);
plot(k,g1,'r',k,g2,'g',k,g3,'b',k,g4,'y'),gridon
title('参数辨识曲线')
end
%计算模型输出ym及系统输出与模型输出之间的误差Ey
fork=1:
198
ym(k)=h(:
k)'*g(:
k);
Ey(k)=y(k+1)-ym(k);
end
k=1:
198
subplot(1,2,2);
plot(k,Ey),gridon
title('误差曲线')
g,ym,Ey
%显示脉冲响应估计值、模型输出及系统输出与模型输出之间的误差
figure(3);
x=0:
1:
4;
y=[0,g(1,198),g(2,198),g(3,198),g(4,198)];
xi=linspace(0,4);
yi=interp1(x,y,xi,'cubic');
plot(x,y,'o',xi,yi,'m'),gridon
title('三次插值曲线')%画出脉冲响应估计值及其三次插值曲线
输入输出数据曲线:
参数辨识曲线和误差曲线:
三次插值曲线:
结果分析:
表梯度校正算法的辨识结果
参数g1g2g3g4
真值-1.50.71.00.5
估计值-1.50670.70621.00750.4904
由
(1)知当输入信号为M序列时,由仿真结果可知,程序设置的E=0.0005,所以程序检测到E=0.0005时,自动跳出。
估计参数的曲线则保持不变。
由于设置的E较大,所以误差曲线波动较明显,如果设置的E更精确,则误差曲线将变得平缓。
由
(2)知当输入信号为方差为1的白噪声时,由仿真结果可知,程序设置的E很小,但从参数估计曲线可以看出,参数估计在50步便趋于稳定,而此时误差也趋于0,误差曲线不再像第一次仿真时那样波动。
3.考虑如下系统
式中,
为白噪声。
其中
为方差为1的白噪声。
取初值
,
。
选择方差为1的正态分布白噪声作为输入信号
,采用递推增广最小二乘算法进行参数估计,自行设计准则停止计算。
要求
i)画出程序流程框图;
ii)画出输入输出数据曲线、参数估计曲线
提示:
对于
,其中
。
待估计参数
,
取数据向量为
。
由于
中的
不可测,采用估计值代替,即
其中
形成的递推公式为
i)程序框图:
N
ii)程序、输入输出曲线、参数估计曲线、误差曲线和结果分析:
程序:
clear%清理工作间变量
L=20;
u=randn(L,1);%u(k)为方差为1的白噪声
y(4)=0;y(3)=0;%设z的前两个初始值为零
nos=randn(L,1);%白噪声
fork=5:
20;%循环变量从5到20
y(k)=1.5*y(k-1)-0.7*y(k-2)+u(k-3)+0.5*u(k-4)+nos(k)-nos(k-1)+0.2*nos(k-2);
%输出采样信号
end
figure
(1),subplot(3,1,1),stem(u),
title('输入u(k)为方差为1的白噪声')
subplot(3,1,2),stem(nos),
title('nos为方差为1的白噪声')
subplot(3,1,3),stem(y),holdon
k=1:
20;
plot(k,y)
title('输出')
%RLS递推最小二乘辨识
c0=[0000000]';%直接给出被辨识参数的初始值,即一个充分小的实向量
p0=10^6*eye(7,7);%直接给出初始状态P0,即一个充分大的实数单位矩阵
c=[c0,zeros(7,19)];%被辨识参数矩阵的初始值及大小
e=zeros(7,20);%相对误差的初始值及大小
fork=5:
20;%开始求K
h1=[-y(k-1),-y(k-2),u(k-3),u(k-4),nos(k),nos(k-1),nos(k-2)]';
x=h1'*p0*h1+1;
x1=inv(x);%开始求K(k)
k1=p0*h1*x1;%求出K的值
d1=y(k)-h1'*c0;c1=c0+k1*d1;%求被辨识参数c
e1=c1-c0;%求参数当前值与上一次的值的差值
e2=e1./c0;%求参数的相对变化
e(:
k)=e2;%把当前相对变化的列向量加入误差矩阵的最后一列
c0=c1;%新获得的参数作为下一次递推的旧参数
c(:
k)=c1;%把辨识参数c列向量加入辨识参数矩阵的最后一列
p1=p0-k1*k1'*[h1'*p0*h1+1];%求出p(k)的值
p0=p1;%给下次用
end%循环结束
c,e%显示被辨识参数及其误差(收敛)情况
%分离参数
a1=c(1,:
);
a2=c(2,:
);
b1=c(3,:
);
b2=c(4,:
);
m1=c(5,:
);
m2=c(6,:
);
m3=c(7,:
);
ea1=e(1,:
);
ea2=e(2,:
);
eb1=e(3,:
);
eb2=e(4,:
);
em1=e(5,:
);
em2=e(6,:
);
em3=e(7,:
);
figure
(2);%第二个图形
i=1:
20;%横坐标从1到20
plot(i,a1,'r',i,a2,':
',i,b1,'g',i,b2,':
',i,m1,'b',i,m2,'y',i,m3,'b')%画出a1,a2,b1,b2的各次辨识结果
title('ParameterIdentificationwithRecursiveLeastSquaresMethod')%图形标题
figure(3);%第三个图形
i=1:
20;%横坐
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 东北大学系统建模作业 12 东北大学 系统 建模 作业
![提示](https://static.bdocx.com/images/bang_tan.gif)