系统辨识与自适应控制学院：专业：学号：姓名：系统辨识与自适应控制作业一、 对时变系统进行参数估计。

系统方程为：y(k)+a(k)y(k-1)=b(k)u(k-1)+e(k) 其中：e(k)为零均值噪声，a(k)= b(k)=要求：1对定常系统（a=0.8,b=0.5）进行结构（阶数）确定和参数估计；2对时变系统，λ取不同值（0.9——0.99）时对系统辨识结果和过程进行比较、讨论3对辨识结果必须进行残差检验 解：一（1）：分析：采用最小二乘法（LS ）：最小二乘的思想就是寻找一个θ的估计值θˆ，使得各次测量的),1(m i Z i =与由估计θˆ确定的量测估计θˆˆi i H Z =之差的平方和最小，由于此方法兼顾了所有方程的近似程度，使整体误差达到最小，因而对抑制误差是有利的。

在此，我应用批处理最小二乘法，收敛较快，易于理解，在系统参数估计应用中十分广泛。

作业程序：clear all;a=[1 0.8]'; b=[ 0.5]'; d=3; %对象参数na=length(a)-1; nb=length(b)-1; %na 、nb 为A 、B 阶次 L=500; %数据长度uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i) yk=zeros(na,1); %输出初值x1=1; x2=1; x3=1; x4=0; S=1; %移位寄存器初值、方波初值 xi=randn(L,1); %白噪声序列theta=[a(2:na+1);b]; %对象参数真值 for k=1:Lphi(k,:)=[-yk;uk(d:d+nb)]'; %此处phi(k,:)为行向量，便于组成phi 矩阵y(k)=phi(k,:)*theta+xi(k); %采集输出数据IM=xor(S,x4); %产生逆M 序列 if IM==0u(k)=-1; elseu(k)=1; endS=not(S); M=xor(x3,x4); %产生M 序列%更新数据x4=x3; x3=x2; x2=x1; x1=M;for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endthetae=inv(phi'*phi)*phi'*y' %计算参数估计值thetae结果：thetae =0.7787 ，0.5103真值=0.8，0.5解：一（2）：采用遗忘因子递推最小二乘参数估计；其仿真算法如下：Step1：设置初值、，及遗忘因子，输入初始数据；Step2：采样当前输入和输出数据；Step3：利用含有遗忘因子的递推公式计算、和；Step4：k=k+1，返回Step2继续循环。

其中：对时变系统，λ取不同值（0.9——0.99）时对系统辨识结果和过程进行比较、讨论作业程序：clear all; close all;a=[1 0.8]’;b=[0.5]’;d=2; %对象参数na=length(a)-1;nb=length(b)-1; %na、nb为A、B阶数L=500; %仿真长度uk=zeros(d+nb,1); %输入初值yk=zeros(na,1); %输出初值u=randn(L,1); %输入采用白噪声序列xi=sqrt(0.1)*randn(L,1); %白噪声序列thetae_1=zeros(na+nb+1,1);p=10^6*eye(na+nb+1);lambda=0.95; %遗忘因子for k=1:L;if k>300;a=[1 0.6]';b=[0.3]';endthetae(:,k)=[a(2:na+1);b];phi=[-yk;uk(d:d+nb)];y(k)=phi'*thetae(:,k)+xi(k); %采集输出数据%递推最小二乘法K=p*phi/(lambda+phi'*p*phi);thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1);p=(eye(na+nb+1)-K*phi')*p/lambda;%更新数据thetae_1=thetae(:,k);for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endsubplot(1,2,1)plot([1:L],thetae(1:na,:));hold on; plot([1:L], thetae(1:na,:),'k:'); xlabel('k');ylabel('参数估计a');axis([0 L -0.5 2]);subplot(1,2,2)plot([1:L],thetae(na+1:na+nb+1,:));holdon;plot([1:L],thetae(na+1:na+nb+1,:),'k:');xlabel('k');ylabel('参数估计b');axis([0 L -0.5 2]);仿真结果图1-1遗忘因子递推最小二乘法的参数估计结果（=0.95）图1-2遗忘因子递推最小二乘法的参数估计结果（λ=0.99）综上所述，可知从仿真图可以看出，当0<k<300时，a的值大约在0.8左右波动，b的值大约在0.5左右波动。

当k>300时，a、b的值发生变化，此时a的值大约在0.6左右波动，b的值大约在0.3左右波动。

所以在此种情况下，可以达到辨识目的，但是由于波动较大，辨识效果不理想。

同时，当λ=0.99时比前面所述λ=0.96时的辨识结果更为明显。

这种情况下当估计值变化稳定时，估计值的变化率比较小，辨识效果较为理想。

一（3）：根据残差平方和估计模型的阶次SISO 过程的差分方程模型的输出残差为)(~k z ，数据长度L ，n H ˆ为n ˆ阶时的数据矩阵，nˆˆθ为n ˆ阶时的参数的估计量，n ˆ为模型阶次估计值，0n 为真实阶次，则残差平方和函数J ： )(~1)ˆ()ˆ(1~~1)ˆ(12ˆˆˆˆˆˆ00k z L H z H z L z z L n J L n n k n n n T n n n n T n ∑++==--==θθ 残差平方和有这样的性质：当L 足够大时，随着n ˆ增加)ˆ(n J 先是显著地下降，当nˆ>0n 时，)ˆ(nJ 值显著下降的现象就终止。

这就是损失函数法来定阶的原理。

程序：N=15;%4位移位寄存器产生的M 序列的周期 y1=1;y2=1;y3=1;y4=0; for i=1:N;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; endy1=x1;y2=x2;y3=x3;y4=x4; end%白噪声的产生A=19;x0=12;M=1024; for k=1:N x=A*x0; x1=mod(x,M); v(k)=x1/512; x0=x1; end%辨识主程序z=zeros(7,N);zs=zeros(7,N);zm=zeros(7,N);zmd=zeros(7,N); z(1)=0;z(2)=0; zs(1)=0; zs(2)=0; zm(1)=0; zm(2)=0; zmd(1)=0;zmd(2)=0;theta=zeros(7,1);p0=10^6*eye(7,7);the=[theta,zeros(7,14)];e=zeros(7,15);for k=3:N;h=[-z(k-1),-z(k-2),u(k-1),u(k-2),v(k),v(k-1),v(k-2)]';x=inv(h'*p0*h+1);q=p0*h*x;d1=z(k)-h'*theta;theta1=theta+q*d1;zs(k)=-1.5*z(k-1)+0.7*z(k-2)+u(k-1)+0.5*u(k-2);zm(k)=[-z(k-1),-z(k-2),u(k-1),u(k-2)]*[theta1(1);theta1(2);theta1(3); theta1(4)];zmd(k)=h'*theta1;e(:,k)=theta1-theta;theta=theta1;the(:,k)=theta1;p1=p0-q*q'*[h'*p0*h+1];p0=p1;endfigure(1);i=1:N;plot(i,the(1,:),'r',i,the(2,:),'r:',i,the(3,:),'b',i,the(4,:),'b:',i, the(5,:),'g',i,the(6,:),'g:',i,the(7,:),'g+')title('辨识曲线')figure(2);i=1:N;plot(i,e(1,:),'r',i,e(2,:),'r:',i,e(3,:),'b',i,e(4,:),'b:',i,e(5,:),' g',i,e(6,:),'g:',i,e(7,:),'r+')title('辨识误差曲线')图1-3-1辨识曲线图1-3-2辨识误差曲线由以上对比曲线可以看出，经过最小二乘法估计得到的数据与实际数据之间虽然存在区别，但是基本符合要求, 说明辨识效果较好。

二、已知系统方程：y(k)+1.2 y(k-1)+0.35 y(k-2)=u(k-2)+0.4 u(k-3)+e(k)+1.1e(k-1)+0.28e(k-2)其中e(k)为白噪声序列。

要求：1、对系统进行参数递推估计；2、进行最小方差控制器设计；3、进行基于最小方差控制的自适应控制系统仿真分析（设输入信号为方波信号）。

二（1）：程序：clear all; close all;a=[1 1.2 0.35]'; b=[1 0.4]'; c=[1 1.1 0.28]'; d=3; %对象参数na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %na、nb、nc为A、B、C阶次L=1000; %仿真长度uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i)yk=zeros(na,1); %输出初值xik=zeros(nc,1); %噪声初值xiek=zeros(nc,1); %噪声估计初值u=randn(L,1); %输入采用白噪声序列xi=sqrt(0.1)*randn(L,1); %白噪声序列theta=[a(2:na+1);b;c(2:nc+1)]; %对象参数thetae_1=zeros(na+nb+1+nc,1); %na+nb+1+nc为辨识参数个数P=10^6*eye(na+nb+1+nc);for k=1:Lphi=[-yk;uk(d:d+nb);xik];y(k)=phi'*theta+xi(k); %采集输出数据phie=[-yk;uk(d:d+nb);xiek]; %组建phie%递推增广最小二乘法K=P*phie/(1+phie'*P*phie);thetae(:,k)=thetae_1+K*(y(k)-phie'*thetae_1);P=(eye(na+nb+1+nc)-K*phie')*P;xie=y(k)-phie'*thetae(:,k); %白噪声的估计值%更新数据thetae_1=thetae(:,k);for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);for i=nc:-1:2xik(i)=xik(i-1);xiek(i)=xiek(i-1);endxik(1)=xi(k);xiek(1)=xie;endfigure(1)plot([1:L],thetae(1:na,:));xlabel('k'); ylabel('参数估计a');legend('a_1','a_2'); axis([0 L -2 2]);figure(2)plot([1:L],thetae(na+1:na+nb+1,:));xlabel('k'); ylabel('参数估计b');legend('b_0','b_1'); axis([0 L -2 2]);figure(3)plot([1:L],thetae(na+nb+2:na+nb+nc+1,:));xlabel('k'); ylabel('参数估计c');legend('c_1','c_2'); axis([0 L -2 2]);对系统进行最广最小二乘的递推估计（RELS），仿真结果如下图2.1所示。