系统辨识作业

系统建模方法大作业1.考虑如下系统y()1.5(1)0.7(2)(3)0.5(4)()kykykukukk式中，()k为白噪声。取初值6(0)10PI，ˆ(0)0。分别选择M序列和方差为1的正态分布白噪声作为输入信号()uk，采用递推最小二乘算法进行参数估计，迭代L=400步停止计算。要求i)给出基本迭代公式；ii)画出程序流程框图；iii)画出输入输出数据曲线、参数估计曲线、误差曲线；提示：产生长度为L方差为1的正态分布白噪声，相应的MATLAB命令为randn（L,1）。110**[*0*1]kphhph（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(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;%四个移位寄存器的输出初始值fori=1:L;%开始循环，长度为Lx1=xor(y3,y4);%第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的x2=y1;%第二个移位寄存器的输入是第一个移位寄存器的输出x3=y2;%第三个移位寄存器的输入是第二个移位寄存器的输出x4=y3;%第四个移位寄存器的输入是第三个移位寄存器的输出y(i)=y4;%取出第四个移位寄存器的幅值为0和1的输出信号，即M序列ify(i)0.5,u(i)=-0.03;%如果M序列的值为1,辨识的输入信号取“-0.03”elseu(i)=0.03;%如果M序列的值为0,辨识的输入信号取“0.03”end%小循环结束y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号做准备end%大循环结束，产生输入信号ufigure(1);%第一个图形title('输入M序列')%图形标题stem(u),gridon%显示出输入信号径线图并给图形加上网格z(4)=0;z(3)=0;z(2)=0;z(1)=0;%设z的前四个初始值为零fork=5:15;%循环变量从3到15z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-3)+0.5*u(k-4);%输出采样信号end%RLS递推最小二乘辨识c0=[0.0010.0010.0010.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(4,4);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005;%取相对误差E=0.000000005c=[c0,zeros(4,14)];%被辨识参数矩阵的初始值及大小e=zeros(4,15);for(n=1:400);%迭代次数fork=5:15;%开始求Kh1=[-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;%求被辨识参数ce1=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%显示被辨识参数及其误差(收敛)情况%分离参数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到15plot(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),gridon%显示出输入信号径线图并给图形加上网格z(4)=0;z(3)=0;z(2)=0;z(1)=0;%设z的前四个初始值为零fork=5:15;%循环变量从3到15z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-3)+0.5*u(k-4);%输出采样信号end%RLS递推最小二乘辨识c0=[0.0010.0010.0010.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(4,4);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005;%取相对误差E=0.000000005c=[c0,zeros(4,14)];%被辨识参数矩阵的初始值及大小e=zeros(4,15);for(n=1:400);%迭代次数fork=5:15;%开始求Kh1=[-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;%求被辨识参数ce1=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%显示被辨识参数及其误差(收敛)情况%分离参数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到15plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:')%画出a1，a2，b1，b2的各次辨识结果的收敛情况title('辨识精度')%图形标题辨识曲线：2.考虑如下确定性系统()1.5(1)0.7(2)(3)0.5(4)ykykykukuk采用梯度校正算法进行参数估计，通过相对误差准则停止计算。取初值ˆ(0)0，并取()()RkckI，其中212()()kiickhk。要求i)给出基本迭代公式；ii)画出程序流程框图；iii)画出输入输出数据曲线、参数估计曲线iv)自行设计适合的数据序列（M序列，白噪声等），并解释仿真结果。(1)迭代公式：ˆˆˆ1()()[()()()]TkkkkykkkggRhhg)](,),(),([)()(1)(2112kkkdiagkhkkNNiiiR准则函数：(2)程序流程框图：工作间清零给被辨识参数g和h赋初始值画出脉冲响应曲线计算模型输出值ym及系统输出与模型输出之间的误差Ey停机计算出权矩阵R(k)计算脉冲响应估计值g给出系统输入输出观测数据并画出图形画出被辨识参数g的各次估计值显示脉冲响应估计值、模型输出值、系统输出与模型输出之间的误差画出系统输出与模型输出之间的误差图二确定性梯度校正辨识的Malab程序流程图min|),(21|)()(2)(kkkJ)()(),(khkyk1,2hkykyk1212[,,,]aabb（3）将题中的差分方程化成：程序如下：1当输入为M序列时L=15;%M序列的周期y1=1;y2=1;y3=1;y4=0;%四个移位寄存器的输出初始值fori=1:L;%开始循环，长度为Lx1=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%大循环结束，产生输入信号uy(4)=0;y(3)=0;y(2)=0;y(1)=0;%设y的前四个初始值为零fork=5:200;%循环变量从5到200y(k)=1.5*y(k-1)-0.7*y(k-2)+u(k-3)+0.5*u(k-4);%输出采样信号endfigure(1),subplot(2,1,1),stem(u),title('ijput')subplot(2,1,2),stem(y),holdonk=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:198h(:,k)=[-y(k),-y(k-1),u(k-2),u(k-3)]';endE=0.0005;%计算权矩阵R(k)和g的估计值fork=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;ife2=E,q=1;break;end%计算脉冲响应的估计值endifq==1forj=k:198g(:,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'),gri