系统辨识大作业

系统辨识大作业系统辨识大作业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）。)答：选择如下形式的辨识模型)()4()3()2()1()(2121kkubkubkzakzakz设输入信号的取值是1．从1k到16k的M序列，采用4位移位寄存器产生，幅度为0.03；2.方差为1、长度为L的正态分布白噪声，用randn(L,1)命令产生；i)基本的迭代公式为1ˆˆˆ()(1)()[()()(1)]1()(1)()()(1)()()()[()()](1)TTTkkkzkkkkkkkkkkkIkkkKhKPhhPhPKhPii)程序流程框图Yz(k)工作间清零产生输出采样信号给被辨识参数和P赋初值计算P(k)计算被辨识参数的相对变化量迭代次数达到L2？停机计算K(k)计算(k)第四个移位寄存器的输出取反，并将幅值变为0.03得到辨识系统的输入信号样本值给M序列的长度L1和移位寄存器的输入赋初始值画出被辨识参数的各次递推估计值图形分离参数画出被辨识参数的相对误差的图形对z(0)、z(1)初始化图1最小二乘递推算法辨识的Malab程序流程图(输入信号为M序列)Yz(k)工作间清零产生输出采样信号给被辨识参数和P赋初值计算P(k)计算被辨识参数的相对变化量迭代次数达到L2？停机计算K(k)计算(k)用randn(L,1)命令产生方差为1的正态白噪声画出被辨识参数的各次递推估计值图形分离参数画出被辨识参数的相对误差的图形对z(0)、z(1)初始化图2最小二乘递推算法辨识的Malab程序流程图(输入信号为正态白噪声)iii)源代码与输入输出曲线、参数估计曲线、误差曲线%最小二乘递推算法辨识程序(输入信号采用M序列)clear%清理工作间变量L1=15;%M序列的周期L2=400;%迭代L2步后停止y1=1;y2=1;y3=1;y4=0;%四个移位寄存器的输出初始值%用4位移位寄存器生成M序列fori=1:L1;%开始循环，长度为L1x1=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);%第1个图形stem(u),gridon%显示出输入信号径线图并给图形加上网格z(2)=0;z(1)=0;%设z的前两个初始值为零fork=3:15;%循环变量从3到15z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);%输出采样信号endfigure(2);%第2个图形stem(z),gridon%显示出输出信号径线图并给图形加上网格%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);%相对误差的初始值及大小num=0;%迭代次数初始化fork=3:15;%开始求Kh1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';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;%给下次用num=num+1;%迭代次数加1ifnum==L2break;%如果参数收敛情况满足要求，终止计算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(3);%第3个图形i=1:15;%横坐标从1到15plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':')%画出a1，a2，b1，b2的各次辨识结果title('ParameterIdentificationwithRecursiveLeastSquaresMethod')%图形标题figure(4);%第4个图形i=1:15;%横坐标从1到15plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:')%画出a1，a2，b1，b2的各次辨识结果的收敛情况title('IdentificationPrecision')%图形标题输入输出曲线：参数估计曲线和误差曲线：%最小二乘递推算法辨识程序(输入信号采用正态白噪声)clear%清理工作间变量L2=400;%迭代L2步后停止%产生白噪声输入信号u=randn(L2,1);%产生长度为L2、方差为1的正态白噪声作为输入figure(1);%绘制输入信号图形stem(u),gridon%显示出输入信号径线图并给图形加上网格z(2)=0;z(1)=0;%设z的前两个初始值为零fork=3:15;%循环变量从3到15z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);%输出采样信号endfigure(2);%绘制输出信号图形stem(z),gridon%显示出输出信号径线图并给图形加上网格%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);%相对误差的初始值及大小num=0;%迭代次数初始化fork=3:15;%开始求Kh1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';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;%给下次用num=num+1;%迭代次数加1ifnum==L2break;%如果参数收敛情况满足要求，终止计算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(3);%第3个图形i=1:15;%横坐标从1到15plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':')%画出a1，a2，b1，b2的各次辨识结果title('ParameterIdentificationwithRecursiveLeastSquaresMethod')%图形标题figure(4);%第4个图形i=1:15;%横坐标从1到15plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:')%画出a1，a2，b1，b2的各次辨识结果的收敛情况title('IdentificationPrecision')%图形标题输入输出曲线：参数估计曲线：2.考虑如下确定性系统()1.5(1)0.7(2)(3)0.5(4)ykykykukuk采用梯度校正算法进行参数估计，通过相对误差准则停止计算。取初值ˆ(0)0，并取()()RkckI，其中212()()kiickhk。要求i)给出基本迭代公式；ii)画出程序流程框图；iii)画出输入输出数据曲线、参数估计曲线iv)自行设计适合的数据序列（M序列，白噪声等），并解释仿真结果。答：i)基本迭代公式：脉冲响应估计值)(k^θ的递推算式为][1(k)(k)y(k)(k)(k)(k))(kTθhhRθθ，初始值取0θ)(0。权矩阵)(kR取如下形式()()RkckI其中212()()kiickhk，所以可选取kiikhkc12)(1)(。ii)程序流程框图工作区间清零给出输入输出观测数据，并绘制其曲线给待辨识参数g、h赋初值计算权矩阵R(k)计算脉冲响应估计值g画出被辨识参数g的各次估计值停止iii)源代码和输入输出曲线、参数估计曲线%梯度校正参数辨识%给定的输入输出数据u=[-1,-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1,1,-1,-1,-1,-1,1];y=[0,-2,-6,-7,-7,-3,5,7,3,-1,5,3,-5,-3,1,-1,1,-5,-7,-7];%画出输入u和输出y的曲线figure(1),subplot(2,1,1),stem(u),subplot(2,1,2),stem(y),holdonk=1:20;plot(k,y)%给出初始值h1=[-1,0,0]';h2=[-1,-1,0]';g=[0,0,0]';I=[1,0,0;0,1/2,0;0,0,1/4];h=[h1,h2,zeros(3,16)];%计算样本数据h(k)fork=3:18h(:,k)=[u(k),u(k-1),u(k-2)]';end%计算权矩阵R(k)和g的估计值fork=1:18a=h(1,k)^2+(h(2,k)^2)+(h(3,k)^2);c=1/a;%计算得到选定的c(k)R=c*I;%计算权矩阵R(k)g(:,k+1)=g(:,k)+R*h(:,k)*(y(k+1)-h(:,k)'*g(:,k));%计算脉冲响应的估计值end%绘制待辨识参数g1、g2、g3的估计曲线g1=g(1,:);g2=g(2,:);g3=g(3,:);figure(2