系统辨识与自适应作业.docx

1、系统辨识与自适应作业 系统辨识与自适应控制学 院： 电气工程学院 专 业： 控制理论与控制工程 学 号： 2010020119 姓 名： 玉海超 系统辨识与自适应控制作业一、系统方程为如下方程，e(k)为零均值白噪声要进行系统辨识并讨论：定常系统a=0.8,b=0.5参数递推估计时变系统，取不同值时递推辨识结果并讨论解答：方法分析利用批处理法得到系统参数的最小二乘法估计式中Y=, , 采用遗忘因子递推最小二乘法，有所学理论知，其参数估计公式如下：其中： 程序代码如下：%批处理最小二乘参数估计（LS）clear all;a=1 0.8; b=0.5; d=1; %对象参数na=length(a)

2、-1; nb=length(b)-1; %na、nb为A、B阶次L=300; %数据长度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:L phi(k,:)=-yk;uk(d:d+nb); %此处phi(k,:)为行向量，便于组成phi矩阵 y(k)=phi(k,:)*theta+xi(k); %采集输出数据 IM=xor(S,x4);

3、 %产生逆M序列 if IM=0 u(k)=-1; else u(k)=1; end S=not(S); M=xor(x3,x4); %产生M序列 %更新数据 x4=x3; x3=x2; x2=x1; x1=M; %for i=d+nb:-1:1 %uk(i)=uk(i-1); %end uk(1)=u(k); %for i=na:-1:1 %yk(i)=yk(i-1); %end yk(1)=y(k);endthetae=inv(phi*phi)*phi*y %计算参数估计值thetae运行结果如下：程序代码如下：%遗忘因子递推最小二乘参数估计（FFRLS）clear all; close

4、all;a=1 0.8; b=0.5; d=1; %对象参数na=length(a)-1; nb=length(b)-1; %na、nb为A、B阶次L=600; %仿真长度uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i)yk=zeros(na,1); %输出初值u=randn(L,1); %输入采用白噪声序列xi=sqrt(0.1)*randn(L,1); %白噪声序列thetae_1=zeros(na+nb+1,1); %thetae初值P=106*eye(na+nb+1);lambda=0.95; %遗忘因子范围0.9 1for k=1:L if k=301 a=

5、1 0.6;b=0.3; %对象参数突变 end theta(:,k)=a(2:na+1);b; %对象参数真值 phi=-yk;uk(d:d+nb); y(k)=phi*theta(:,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:2 uk(i)=uk(i-1); end uk(1)=u

6、(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(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,:);hold on;plot(1:L,thetae(na+1:na+nb+1,:),k:);xlabel(k);ylabel(参数估计b);axis(0 L -0.5 2);仿真

7、结果：依次为=0.95、=0.98、=1的结果=0.95当遗忘因子=0.95时，从仿真图可以看出，辨识结果很明显。开始K=0，估计值a的值比较小，而b的值比较大。随着K的增大估计值a增大，b减小，并且变化率很大，当K=20左右变化比较缓慢了，一直到k=300，a围绕0.8上下波动，b围绕0.6上下波动，波动范围0.2左右；当k300的时候，由于a和b变化了，所以当变化稳定的时候a就围绕0.6变化，b围绕0.3变化。总的来说当=0.95时，系统达到辨识的要求，但是辨识效果不太理想。=0.98当遗忘因子=0.98时，可以从图中看出和=0.95的时候存在一定的差别。首先是当K=0的时候，a,b的值都

8、很大，随着K的增大，a,b开始减小，当K=300的时候，a,b接近真值，当K大于300的时候，两者围绕变化的数值变为0.6和0.3。可以看出当=0.98的时候辨识结果明显。=1当=1的时候我们称为此时为递推最小二乘法，此时就不带遗忘因子了。从图中可以从看出当K=0的时候出现的结果与=0.98时的一样，但是当K100后变化缓慢，最终a达到0.8左右，b达到0.5左右。但是当K300系统估计值a和b却没有变化到0.6和0.3。所以可以看出系统辨识结果也不太理想。二、已知系统方程为其中e(K)为白噪声，在输入信号为方波时，分析（白噪声为信号幅值的10%，最大为20%）系统开环响应在PID控制下系统闭

9、环响应在最小方差控制下，系统闭环响应（要有分析、比较说明）解：增量式PID控制算法，如下式所示u(k)=kp(e(k)-e(k-1)+kie(k)+kd(e(k)-2e(k-1)+e(k-2)u(k)=u(1)+ u(k);式中，kp，ki，kd为PID调节参数，且u(k)=u(k)-u(k-1) e(k)=yr(k)-y(k)最小方差控制算法：1) 输入初始数据2) 采样当前实际输出y(k)和期望输出Yr(k+d)(或Yr(k);3) 求解Diophantine方程，得到多项式E、F、G的系数4) 采用式计算并实施u(k);5) 返回第二步（kk+1）,继续循环程序代码如下：c%白噪声及有色

10、噪声序列的产生clear all; close all;L=400; %仿真长度d=1 1.2 0.35; j=1 0.4; c=1 1.1 0.28;a=2;%D、C多项式的系数(可用roots命令求其根)nd=length(d)-1; nc=length(c)-1; nj=length(j)-1; %nd、nc为D、C的阶次xik=zeros(nc,1); %白噪声初值，相当于(k-1).(k-nc)ek=zeros(nd,1); %输出Y初始值eiu=zeros(a+nj,1);%输入信号的初始值xi=randn(L,1); %randn产生均值为0，方差为1的高斯随机序列（白噪声序列）

11、err=zeros(L,1);u=10*ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4,1); %产生方波for k=1:L time(k)=k; e(k)=-d(2:nd+1)*ek+c*xi(k);xik+j*eiu(a:a+nj); %系统的实际输出 err=u(k)-e(k);%计算偏差 %数据更新 for i=nd:-1:2 ek(i)=ek(i-1); end ek(1)=e(k); for i=nc:-1:2 xik(i)=xik(i-1); end xik(1)=xi(k); for i=nj+a:-1:2 eiu(i)=eiu(i

12、-1); end eiu(1)=u(k);endsubplot(2,1,1);plot(time,u(1:L),r:,time,e);xlabel(t); ylabel(y_r(t)、y(t);legend(y_r(t),y(t);subplot(2,1,2);plot(time,err,r);%输出偏差图形xlabel(t); ylabel(e(t);仿真结果：程序代码：%PID控制clear all; close all;a=1 1.2 0.35; b=1 0.4; c=1 1.1 0.28; d=2; %对象参数na=length(a)-1; nb=length(b)-1; nc=len

13、gth(c)-1; %na、nb、nc为多项式A、B、C阶次kp=0.2; ki=0.5; kd=0.1; %PID控制器参数（试凑法L=400; %控制步数uk=zeros(d+nb,1); %控制初值：uk(i)表示u(k-i);yk=zeros(na,1); %输出初值xik=zeros(nc,1); %白噪声初值ek=zeros(2,1); %输出误差初值yr=10*ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4+d,1); %期望输出xi=sqrt(0.1)*randn(L,1); %白噪声序列for k=1:L time(k)=k; y

14、(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+c*xi(k);xik;%采集输出数据 e(k)=yr(k)-y(k); %增量式PID控制律 du=kp*(e(k)-ek(1)+ki*e(k)+kd*(e(k)-2*ek(1)+ek(2); u(k)=uk(1)+du; %更新数据 for i=d+nb:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k); for i=nc:-1:2 xik(i)=xik(i-1); end if nc0 xik(1)=xi(k); en

15、d for i=d+nb:-1:2 ek(i)=ek(i-1); end ek(1)=e(k);endsubplot(2,1,1);plot(time,yr(1:L),r:,time,y);xlabel(t); ylabel(y_r(t)、y(t);legend(y_r(t),y(t);subplot(2,1,2);plot(time,e);xlabel(t); ylabel(e(t);仿真结果：程序代码%最小方差控制（MVC）clear all; close all;a=1 1.2 0.35; b=1 0.4; c=1 1.1 0.28; d=2; %对象参数na=length(a)-1;

16、nb=length(b)-1; nc=length(c)-1; %na、nb、nc为多项式A、B、C阶次nf=nb+d-1; %nf为多项式F的阶次L=400; %控制步数uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i);yk=zeros(na,1); %输出初值yrk=zeros(nc,1); %期望输出初值xik=zeros(nc,1); %白噪声初值yr=10*ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4+d,1); %期望输出xi=sqrt(0.1)*randn(L,1); %白噪声序列ek=zeros(L,1)

17、;e,f,g=sindiophantine(a,b,c,d); %求解单步Diophantine方程for k=1:L time(k)=k; y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+c*xi(k);xik;%采集输出数据 ek(k)=yr(k)-y(k); u(k)=(-f(2:nf+1)*uk(1:nf)+c*yr(k+d:-1:k+d-min(d,nc);yrk(1:nc-d)-g*y(k);yk(1:na-1)/f(1);%求控制量 %更新数据 for i=d+nb:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2

18、 yk(i)=yk(i-1); end yk(1)=y(k); for i=nc:-1:2 yrk(i)=yrk(i-1); xik(i)=xik(i-1); end if nc0 yrk(1)=yr(k); xik(1)=xi(k); endendsubplot(2,1,1);plot(time,yr(1:L),r:,time,y);xlabel(k); ylabel(y_r(k)、y(k);legend(y_r(k),y(k);subplot(2,1,2);plot(time,ek);xlabel(k); ylabel(e(k);function e,f,g=sindiophantine(

19、a,b,c,d)%* %功能：单步Diophanine方程的求解 %调用格式：e,f,g=sindiophantine(a,b,c,d) %输入参数：多项式A、B、C系数（行向量）及纯滞后（共4个） %输出参数：Diophanine方程的解e,f,g（共3个）%*na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %A、B、C的阶次ne=d-1; ng=na-1; %E、G的阶次ad=a,zeros(1,ng+ne+1-na); cd=c,zeros(1,ng+d-nc); %令a(na+2)=a(na+3)=.=0e(1)=1;for i=2:ne

20、+1 e(i)=0; for j=2:i e(i)=e(i)+e(i+1-j)*ad(j); end e(i)=cd(i)-e(i); %计算eiendfor i=1:ng+1 g(i)=0; for j=1:ne+1 g(i)=g(i)+e(ne+2-j)*ad(i+j); end g(i)=cd(i+d)-g(i); %计算giendf=conv(b,e); %计算F仿真结果如下：最小方差控制放大图分析如下：从仿真图可以看出，系统的开环响应，实际的输出与期望的误差另很大，超调量也比较高，这主要是开环模型的一个效果。当系统引进闭环的PID控制时，系统的超调量比较小，上升时间也比较快，最后稳态

21、误差几乎为零，但在信号改变的那一刻，也即波峰与波谷转变的那瞬间，系统的误差也比较大。PID控制能够达到一定的控制效果，但PID参数的调整比较麻烦，需要反复的试奏才能达到一个满意的控制效果。从最小方差控制的仿真图来看，刚开始超调量有点大，误差也比较大。但很快的，随着控制的深入，系统的误差很小，几乎与输入信号的波形重合，在信号发生变化的瞬间，控制效果也很好，但将图形放大，还是有点小误差，但总体上，控制效果还是比较好的。三、进行基于波波夫稳定性理论的MARS设计及算法仿真算法的模型图如下 程序算法：1) 选择参考模型，即Gm(s);2) 选择输入信号yr(t)3) 采样当前参考模型输出ym(t)和系

22、统实际输出yp(t);4) 利用式;计算u(t)，其中r为输入信号，v=e*(d0+d1*s), 为的初始状态5) tt+h,返回第3步，继续循环考虑如下被控对象模型：kp未知（仿真时取kp=1）选择参考模型为：Km=1程序代码如下：%波波夫程序clear all; close all;h=0.01; L=100/h; %数值积分步长和仿真步数（减小h，可以提高积分精度）num=2 1; den=1 2 1; n=length(den)-1; %对象参数（严格正实）d0=1;d1=5;ek=zeros(n,1);kp=1; Ap,Bp,Cp,Dp=tf2ss(kp*num,den); %对象参

23、数（传递函数型转换为状态空间型）km=1; Am,Bm,Cm,Dm=tf2ss(km*num,den); %参考模型参数yr0=0; u0=0; e0=0;v0=0; %初值xp0=zeros(n,1); xm0=zeros(n,1); %状态向量初值kc0=0; %可调增益初值r=2; yr=r*ones(1,L/4) -ones(1,L/4) ones(1,L/4) -ones(1,L/4); %输入信号for k=1:L time(k)=k*h; xp(:,k)=xp0+h*(Ap*xp0+Bp*u0); yp(k)=Cp*xp(:,k); %计算yp xm(:,k)=xm0+h*(Am

24、*xm0+Bm*yr0); ym(k)=Cm*xm(:,k); %计算ym e(k)=ym(k)-yp(k); %e=ym-yp v(k)=d0*ek(n)+d1*(ek(n)-ek(n-1)/h); kc=kc0+1*h*v0*yr0+2*v0*yr0; %波波夫算法 u(k)=kc*yr(k); %控制量 %更新数据 yr0=yr(k); u0=u(k);v0=v(k); for i=n:-1:2 ek(n)=ek(n-1); end ek(1)=e(k); xp0=xp(:,k); xm0=xm(:,k); kc0=kc;endsubplot(2,1,1);plot(time,ym,r,time,yp,:);xlabel(t); ylabel(y_m(t)、y_p(t);legend(y_m(t),y_p(t);subplot(2,1,2);plot(time,u);xlabel(t); ylabel(u(t);仿真结果：