系统辨识作业

1、系统建模方法大作业1.考虑如下系统式中，为白噪声。取初值，。分别选择M序列和方差为1的正态分布白噪声作为输入信号，采用递推最小二乘算法进行参数估计，迭代L=400步停止计算。要求i) 给出基本迭代公式；ii) 画出程序流程框图；iii) 画出输入输出数据曲线、参数估计曲线、误差曲线；提示：产生长度为L方差为1的正态分布白噪声，相应的MATLAB命令为 randn（L,1）。（1） z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-3)+0.5*u(k-4) h1=-z(k-1),-z(k-2),u(k-3),u(k-4) c1=c0+k1*z(k)-h1'*c0（2）Yz(

2、k)工作间清零产生输出采样信号给被辨识参数和P赋初值计算P(k)计算被辨识参数的相对变化量满足迭代次数?停机计算K(k)计算(k)第四个移位寄存器的输出取反，并将幅值变为0.03得到辨识系统的输入信号样本值给M序列的长度L和移位寄存器的输入赋初始值 画出被辨识参数c的各次递推估计值图形分离参数画出被辨识参数c的相对误差的图形画出辨识的输入信号曲线图形图 最小二乘递推算法辨识的Malab程序流程图（3）1.当输入为M序列时L=15; % M序列的周期y1=1;y2=1;y3=1;y4=0; %四个移位寄存器的输出初始值for i=1:L; %开始循环，长度为L x1=xor(y3,y4); %第

3、一个移位寄存器的输入是第三个与第四个移位寄存器的输出的 x2=y1; %第二个移位寄存器的输入是第一个移位寄存器的输出 x3=y2; %第三个移位寄存器的输入是第二个移位寄存器的输出 x4=y3; %第四个移位寄存器的输入是第三个移位寄存器的输出 y(i)=y4; %取出第四个移位寄存器的幅值为"0"和"1"的输出信号，即M序列 if y(i)>0.5,u(i)=-0.03; %如果M序列的值为"1", 辨识的输入信号取“-0.03” else u(i)=0.03; %如果M序列的值为"0", 辨识的输入信号

4、取“0.03” end %小循环结束y1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号做准备end %大循环结束，产生输入信号u figure(1); %第一个图形title('输入M序列') %图形标题 stem(u),grid on %显示出输入信号径线图并给图形加上网格z(4)=0;z(3)=0;z(2)=0;z(1)=0; %设z的前四个初始值为零for k=5:15; %循环变量从3到15 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-3)+0.5*u(k-4); %输出采样信号 end%RLS递推最小二乘辨识c0=0.001 0.

5、001 0.001 0.001' %直接给出被辨识参数的初始值,即一个充分小的实向量p0=106*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005; %取相对误差E=0.000000005c=c0,zeros(4,14); %被辨识参数矩阵的初始值及大小e=zeros(4,15); for(n=1:400); %迭代次数for k=5:15; %开始求K h1=-z(k-1),-z(k-2),u(k-3),u(k-4)'x=h1'*p0*h1+1; x1=inv(x); %开始求K(k) k1=p0*h1*x1; %求出K

6、的值 d1=z(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 %小循环结束end %大循环结束c,e %显示被辨识参数及

7、其误差(收敛)情况%分离参数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:15; %横坐标从1到15plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':') %画出a1，a2，b1，b2的各次辨识结果title('递推最小二乘参数辨识') %图形标题figure(3); %第三个图形i=1:15; %横坐标从1到15p

8、lot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:') %画出a1，a2，b1，b2的各次辨识结果的收敛情况title('辨识精度') %图形标题 辨识曲线图：2当输入为随机正态分布时u=randn(L,1); %产生随机正态分布白噪声figure(1); %第一个图形title('输入正态分布白噪声') %图形标题 stem(u),grid on %显示出输入信号径线图并给图形加上网格z(4)=0;z(3)=0;z(2)=0;z(1)=0; %设z的前四个初始值

9、为零for k=5:15; %循环变量从3到15 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-3)+0.5*u(k-4); %输出采样信号 end%RLS递推最小二乘辨识c0=0.001 0.001 0.001 0.001' %直接给出被辨识参数的初始值,即一个充分小的实向量p0=106*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005; %取相对误差E=0.000000005c=c0,zeros(4,14); %被辨识参数矩阵的初始值及大小e=zeros(4,15); for(n=1:400); %迭代次数for k=

10、5:15; %开始求K h1=-z(k-1),-z(k-2),u(k-3),u(k-4)'x=h1'*p0*h1+1; x1=inv(x); %开始求K(k) k1=p0*h1*x1; %求出K的值 d1=z(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

11、-k1*k1'*h1'*p0*h1+1; %求出 p(k)的值 p0=p1; %给下次用 %如果参数收敛情况满足要求，终止计算 end %小循环结束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:15; %横坐标从1到15plot(i,a1,'r',i,a2,':',i,b1,'g',i,

12、b2,':') %画出a1，a2，b1，b2的各次辨识结果title('递推最小二乘参数辨识') %图形标题figure(3); %第三个图形i=1:15; %横坐标从1到15plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:') %画出a1，a2，b1，b2的各次辨识结果的收敛情况title('辨识精度') %图形标题 辨识曲线：2.考虑如下确定性系统采用梯度校正算法进行参数估计，通过相对误差准则停止计算。取初值，并取，其中。要求i) 给出基本

13、迭代公式；ii) 画出程序流程框图；iii) 画出输入输出数据曲线、参数估计曲线iv) 自行设计适合的数据序列（M序列，白噪声等），并解释仿真结果。(1)迭代公式： 准则函数：(2)程序流程框图：工作间清零给被辨识参数g和h赋初始值画出脉冲响应曲线计算模型输出值ym及系统输出与模型输出之间的误差Ey停机 计算出权矩阵R(k) 计算脉冲响应估计值g给出系统输入输出观测数据并画出图形画出被辨识参数g的各次估计值显示脉冲响应估计值、模型输出值、系统输出与模型输出之间的误差画出系统输出与模型输出之间的误差图二 确定性梯度校正辨识的Malab程序流程图（3）将题中的差分方程化成： 程序如下：1 当输入为

14、M序列时L=15; % M序列的周期y1=1;y2=1;y3=1;y4=0; %四个移位寄存器的输出初始值for i=1:L; %开始循环，长度为L x1=xor(y3,y4); %第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或” x2=y1; %第二个移位寄存器的输入是第一个移位寄存器的输出 x3=y2; %第三个移位寄存器的输入是第二个移位寄存器的输出 x4=y3; %第四个移位寄存器的输入是第三个移位寄存器的输出 m(i)=y4; %取出第四个移位寄存器的幅值为"0"和"1"的输出信号，即M序列 if m(i)>0.5,u(i)

15、=-1; %如果M序列的值为"1", 辨识的输入信号取“-0.03” else u(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的前四个初始值为零for k=5:200; %循环变量从5到200 y(k)=1.5*y(k-1)-0.7*y(k-2)+u(k-3)+0.5*u(k-4); %输出采样信号 endfigure(1), s

16、ubplot(2,1,1), stem(u),title('ijput')subplot(2,1,2), stem(y), hold onk=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)for k=4:198 h(:,k)=-y(k),-y(k-1),u(k-2

17、),u(k-3)'endE=0.0005;%计算权矩阵R(k)和g的估计值for k=1:198a=(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; if e2<=E,q=1;break; end %计算脉冲响应的估计值 end if q=1 for j=k:198 g(:,j)=g(:,k);end %绘

18、图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'), grid ontitle('参数辨识曲线') 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&

19、#39;,k,g3,'b',k,g4,'y'), grid on title('参数辨识曲线') end%计算模型输出ym及系统输出与模型输出之间的误差Eyfor k=1:198ym(k)=h(:,k)'*g(:,k); Ey(k)=y(k+1)-ym(k);end k=1:198subplot(1,2,2); plot(k,Ey), grid ontitle('误差曲线')g, ym, Ey %显示脉冲响应估计值、模型输出及系统输出与模型输出之间的误差 figure(3);x=0:1:4;y=0,g(1,198),g(2

20、,198),g(3,198),g(4,198);xi=linspace(0,4);yi=interp1(x,y,xi,'cubic');plot(x,y,'o',xi,yi,'m'), grid on title('三次插值曲线') %画出脉冲响应估计值及其三次插值曲线2 输入为白噪声时clear %清理工作间变量u=randn(200,1);y(4)=0;y(3)=0; y(2)=0;y(1)=0; %设y的前四个初始值为零for k=5:200; %循环变量从5到200 y(k)=1.5*y(k-1)-0.7*y(k-2)+u

21、(k-3)+0.5*u(k-4); %输出采样信号 endfigure(1), subplot(2,1,1), stem(u),title('方差为1的正态分布白噪声')subplot(2,1,2), stem(y), hold onk=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);%计算样本数据

22、h(k)for k=4:198 h(:,k)=-y(k),-y(k-1),u(k-2),u(k-3)'endE=0.000000000005;%计算权矩阵R(k)和g的估计值for k=1:198a=(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; if e2<=E,q=1;break; end %计算脉冲

23、响应的估计值 end if q=1 for j=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'), grid ontitle('参数辨识曲线') else g1=g(1,:); g2=g(2,:); g3=g(3,:);g4=g(4,:);figure(2); k=1:199; s

24、ubplot(1,2,1); plot(k,g1,'r',k,g2,'g',k,g3,'b',k,g4,'y'), grid on title('参数辨识曲线') end%计算模型输出ym及系统输出与模型输出之间的误差Eyfor k=1:198ym(k)=h(:,k)'*g(:,k); Ey(k)=y(k+1)-ym(k);endk=1:198subplot(1,2,2); plot(k,Ey), grid ontitle('误差曲线')g, ym, Ey %显示脉冲响应估计值、模型输出及系

25、统输出与模型输出之间的误差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'), grid on title('三次插值曲线') %画出脉冲响应估计值及其三次插值曲线辨识曲线： 3. 考虑如下系统式中，为白噪声。其中为方差为1的白噪声。取初值，。选择方差为1的正态分布白噪声作为输入信号，采用递推增广最小二乘算法进行参数估计，自行设计准则

26、停止计算。要求i) 画出程序流程框图；ii) 画出输入输出数据曲线、参数估计曲线提示：对于，其中。待估计参数，取数据向量为。由于中的不可测，采用估计值代替，即其中形成的递推公式为（1）程序框图:Yy(k)工作间清零产生输出采样信号给被辨识参数和P赋初值计算P(k)计算被辨识参数的相对变化量参数收敛满足要求？停机计算K(k)计算(k)第四个移位寄存器的输出取反，并将幅值变为1得到辨识系统的输入信号样本值给M序列的长度L和移位寄存器的输入赋初始值和白噪声画出被辨识参数的各次递推估计值图形分离参数画出被辨识参数的相对误差的图形画出辨识的输入信号曲线图形图 最小二乘递推算法辨识的Malab程序流程图N

27、（2）程序：L=20;u=randn(L,1); %u(k)为方差为1的白噪声y(4)=0;y(3)=0; %设z的前两个初始值为零nos=randn(L,1); %白噪声for k=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); %输出采样信号 endfigure(1), subplot(3,1,1), stem(u), title('输入u(k)为方差为1的白噪声')subplot(3,1,2), stem(nos),title('

28、误差为方差为1的白噪声')subplot(3,1,3), stem(y), hold onk=1:20;plot(k,y)title('输出') %RLS递推最小二乘辨识c0=0 0 0 0 0 0 0' %直接给出被辨识参数的初始值,即一个充分小的实向量p0=106*eye(7,7); %直接给出初始状态P0，即一个充分大的实数单位矩阵c=c0,zeros(7,19); %被辨识参数矩阵的初始值及大小e=zeros(7,20); %相对误差的初始值及大小for k=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