自适应大作业(共21页).docx

《自适应大作业(共21页).docx》由会员分享，可在线阅读，更多相关《自适应大作业(共21页).docx（21页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上自适应控制大作业给定如下双容水箱对象，其中阀门开度为0100%；液位高度单位为cm（下水箱液位最高为80cm）。其数学模型为：，采样周期为1s.其中对应连续传递函数：，描述测量噪声对系统的影响：。要求完成如下任务：（1） 图中 有均值为0,方差为0.005的白噪声, 施加M序列信号,采集过程输入输出数据；（2） 使用递推增广最小二乘法辨识过程模型参数，并分析其辨识效果；（3） 自行设置某个过程参数随时间发生变化，编程实现遗忘因子增广递推最小二乘法，分析其与递推增广最小二乘法的效果差异；（4） 以下水箱液位 为被控对象，阀门开度u为控制量，施加PID控制,分析控制效果；

2、（5） 针对上述过程设计最小方差控制控制器 分析其原理并编程实现，并与PID控制效果进行对比分析；（6） 针对上述过程设计最小方差自适应控制控制器 分析其原理并编程实现，最后并与PID控制进行对比分析。（7） 不考虑噪声影响，设计一种模型参考自适应控制器，并编程实现；解：（1）由双容水箱对象数学模型，最终确定对象的递进函数关系如下：利用matlab进行编程，得到输入数据为：图一 输入m序列输出数据为：图二 输出随机数据（2）利用递推最小二乘法进行matlab仿真曲线为：图三 递推最小二乘法最终所得辨识参数为：a1=-1.7655，a2=0.7749，b6=0.3062，b7=0.2794。由此

3、可以看出，所辨识得到参数与原模型参数相差较小，辨识效果较好。（3）令原模型参数b6在500s时增大，则通过带遗忘因子的最小二乘方法得到结果如下：图四 遗忘因子最小二乘法所辨识得到参数为a1=-1.7639，a2=0.7732，b6=0.5048，b7=0.2763。将相同变化施加于递推最小二乘法辨识，所得曲线为：图五 递推最小二乘法所辨识得到参数为a1=-1.7704，a2= 0.7792，b6=0.3910，b7=0.2663。综上可以比较出，对于过程参数随时间发生变化的情况，带遗忘因子的递推最小二乘法能够较好的适应改变的过程条件，辨识得到参数与改变后的参数相差较小。但是递推最小二乘法不能适

4、应改变的过程条件，无法辨识出改变后的过程参数，所以，带遗忘因子的递推最小二乘法辨识效果相较于递推最小二乘法要好。（4）利用PID控制器控制该二容水箱下水箱液位曲线如下：图六 PID控制曲线由上图可以看出，PID控制器能够使下水箱液位最终稳定于稳态点，且超调量、调节时间和上升时间都较小，控制效果较好。（5）最小方差控制器原理简介：被控对象数学模型：其中：这是一个 d 步延时系统，t 时刻控制 u(t) 只能影响 t+d 时刻及以后的输出。要想得到t时刻系统的最优输出量u(t)，就必须对t+d时刻的输出y( t+d)进行预测（预测它与u(t)的关系）。通过“输出方差最小”这个性能指标的约束，求得u

5、(t)的表达是即为最优输出。在t时刻，根据模型计算t+d时刻的输出即 。由与 可得由上式及可得又因为，所以最终其中，且。所以。预测误差为：优化目标函数为预测误差方差最小化：所以最优预测模型为：预测方差为：。利用最小方差控制器设计该系统，被控对象数学模型为得到最小方差控制效果为：图7 最小方差控制仿真曲线由上图可以看出，最小方差控制器可以跟踪方波的输入变化而变化，但有较大的毛刺产生，控制效果不是太好。将其与PID进行比较：图8 最小方差与PID控制效果比较由上图可以看出，由于最小方差虽然能跟踪方波的变化而变化，但由于产生较多毛刺，使得控制效果较差，与PID效果差的较大。下面将设计最小方差直接自校

6、正控制器。（6）在最小方差控制器的基础上设计最小方差直接自校正控制器，即将参数估计与最小方差控制器结合起来，通过在每一次循环时利用递推最小二乘估计模型参数以达到自适应的效果。编程步骤： 输入对象的结构和参数，为做对象仿真做准备； 为RELS算法做准备，设定 生成白噪声和期望输出 进入循环，计算对象的最新时刻的输出，并得到期望输出； 利用RELS算法估计模型参数A,B,C 求解Diophantine方程，得到E,F,G； 计算最优控制律 t=t+1，返回4。利用最小方差直接自校正控制器控制仿真效果为：图9 输入输出仿真曲线图10 参数辨识仿真曲线与PID的控制效果进行对比曲线为：图11 自校正与

7、PID控制效果对比曲线由上图可以看出，最小方差直接自校正控制效果相比与最小方差控制都能够跟踪设定值方波的变化而变化，但最小方差直接自校正控制产生的毛刺比最小方差控制产生的毛刺要稀疏一点且毛刺的幅值要更小，可能是因为最小方差直接自校正控制能够利用递推最小二乘法在线辨识模型参数，使得控制性能较原来更好。但就总体控制效果而言，都不如PID的控制效果好。（7）设计一种模型参考自适应控制器仿真结果如下所示：1）当期望输出为方波时：图12 模型参考自适应仿真图图13 参数估计1图14 参数估计22）当期望输出为随机序列时：图15模型参考自适应仿真图图16 参数估计1图17 参数估计2程序附录：（1）cle

8、ar allclcclose allL=1000;n = 10;p = 2n -1; ms = idinput(p, prbs,0,1,0,1);u=ms(1:L,1);v=gauss_rand(0,0.005,L);z=zeros(1,L); z(1,1:7)=0;for i=8:Lz(1,i)=1.765*z(1,i-1)-0.7745*z(1,i-2)+0.3064*u(1,i-6)+0.2814*u(1,i-7)+v(1,i)+0.5*v(1,i-1)+0.3*v(1,i-2);endfigurestairs(z)title(输出数据 )（2）clear allclcclose all

9、L=1000;lambda=1*ones(1,L);lambda_L=diag(lambda);n = 10; p = 2n -1; ms = idinput(p, prbs,0,1,0,1);u=ms(1:L,1);v=gauss_rand(0,0.005,L);z=zeros(1,L); z(1,1:7)=0;for i=8:L z(1,i)=1.765*z(1,i-1)-0.7745*z(1,i-2)+0.3064*u(1,i-6)+0.2814*u(1,i-7)+v(1,i)+0.5*v(1,i-1)+0.3*v(1,i-2);endZ_L=z(1,1:L);V_L=v(1,1:L);

10、phi_L=zeros(L,4);for j=8:L phi_L(j,:)=-1*z(1,j-1) -1*z(1,j-2) 0 0+0 0 u(1,j-6) u(1,j-7);endtheta_WLS=inv(phi_L*lambda_L*phi_L)*phi_L*lambda_L*Z_Lz0=z(1,1);phi0=phi_L(1,:);theta0=zeros(4,1);theta=theta0;a=1000;p0=(a2)*eye(4,4);p=p0;phi=phi0;zk=z0;for k=1:999 l=p*phi/(phi*p*phi+1/lambda(1,k); theta=th

11、eta+l*(zk-phi*theta); p=(eye(4,4)-l*phi)*p; phi=phi_L(k+1,:); zk=z(1,k+1); theta1(:,k)=theta;endthetasubplot(2,2,1)plot(theta1(1,:)title(a1)subplot(2,2,2)plot(theta1(2,:)title(a2)subplot(2,2,3)plot(theta1(3,:)title(b6)subplot(2,2,4)plot(theta1(4,:)title(b7)（3）clear allclcclose allL=1000;lambda=1*one

12、s(1,L);lambda_L=diag(lambda);n = 10; p = 2n -1; ms = idinput(p, prbs,0,1,0,1);u=ms(1:L,1);v=gauss_rand(0,0.005,L);z=zeros(1,L); z(1,1:7)=0;for i=8:Lz(1,i)=1.765*z(1,i-1)-0.7745*z(1,i-2)+0.3064*u(1,i-6)+0.2814*u(1,i-7)+v(1,i)+0.5*v(1,i-1)+0.3*v(1,i-2); if i500 z(1,i)=1.765*z(1,i-1)-0.7745*z(1,i-2)+0.

13、5064*u(1,i-6)+0.2814*u(1,i-7)+v(1,i)+0.5*v(1,i-1)+0.3*v(1,i-2); endendZ_L=z(1,1:L);V_L=v(1,1:L);phi_L=zeros(L,4);for j=8:L phi_L(j,:)=-1*z(1,j-1) -1*z(1,j-2) 0 0+0 0 u(1,j-6) u(1,j-7);endtheta_WLS=inv(phi_L*lambda_L*phi_L)*phi_L*lambda_L*Z_Lz0=z(1,1);phi0=phi_L(1,:);theta0=zeros(4,1);theta=theta0;a=

14、1000;p0=(a2)*eye(4,4);p=p0;phi=phi0;zk=z0;alpha=0.99;for k=1:999 l=p*phi/(phi*p*phi+alpha); theta=theta+l*(zk-phi*theta); p=(eye(4,4)-l*phi)*p/alpha; phi=phi_L(k+1,:); zk=z(1,k+1); theta1(:,k)=theta;endthetasubplot(2,2,1)plot(theta1(1,:)title(a1)subplot(2,2,2)plot(theta1(2,:)title(a2)subplot(2,2,3)p

15、lot(theta1(3,:)title(b6)subplot(2,2,4)plot(theta1(4,:)title(b7)（4）clcclear allclose allTs=1;L=1000;u=0.6 0.6 0.6 0.6 0.6 0.6 0.6;v(1:L)=0;z=zeros(1,L); z(1:7)=0;uu(1:7)=u(1:7);H=z;%PID%Kp=0.04;Ti=60;Td=0;e_1=0;e_2=0;for i=8:L if i=350&i=800&i1001 H2_sp(i)=25; ende(i)=H2_sp(i)-H(i-1);du(i)=Kp*(e(i)-e

16、_1)+Ts/Ti*e(i)+Td/Ts*(e(i)-2*e_1+e_2);uu(i)=uu(i-1)+du(i);if uu(i)1 uu(i)=1;endu(i)=uu(i);e_2=e_1;e_1=e(i);z(i)=1.765*z(i-1)-0.7745*z(i-2)+0.3064*u(i-6)+0.2814*u(i-7)+v(i)+0.5*v(i-1)+0.3*v(i-2);H(i)=z(i); endt=1:Ts:L;subplot(211)plot(t,H,t,H2_sp)xlabel(t/s)title(下水箱H2)subplot(212),plot(t,uu)title(控制

17、作用u)（5）clear all; close all;a=1 -1.765 0.7745; b=0.3064 0.2814; c=1 0.5 0.3;d=6; %对象参数na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %nanbnc为多项式A,B,C的阶次nf=nb+d-1; %nf多项式F的阶次L=1000; %控制步数uk=zeros(d+nb,1); %输入初值yk=zeros(na,1); %输出初值yrk=zeros(nc,1); %期望输出初值xik=zeros(nc,1); %白噪声初值yr=40*ones(L/4,1);20*o

18、nes(L/4,1);40*ones(L/4,1);20*ones(L/4+d,1); %期望输出xi=sqrt(0.1)*randn(L,1); %白噪声序列e,f,g=sindiophantine(a,b,c,d); %求解Diophantine方程for k=1:L t(k)=k; y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+c*xi(k);xik;%采集输出数据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

19、=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 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(t,yr(1:L),r:,t,y);xlabel(k); ylabel(h2_r(k)h(k);legend(h2_r,h2);axis(0 L 10 60);subplot(2,1,2);plot(

20、t,u);xlabel(k); ylabel(u(k);（6）clear all; close all;a=1 -1.765 0.7745; b=0.3064 0.2814; c=1 0.5 0.3; d=6; %对象参数na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %nanbnc为多项式A,B,C的阶次nf=nb+d-1; ng=na-1; %nfng多项式F,G的阶次 L=1000; %控制步数uk=zeros(d+nf,1); %输入初值yk=zeros(d+ng,1); %输出初值yek=zeros(nc,1); %最优预测输出估计初值

21、yrk=zeros(nc,1); %期望输出初值xik=zeros(nc,1); %白噪声初值yr=40*ones(L/4,1);20*ones(L/4,1);40*ones(L/4,1);20*ones(L/4+d,1); %期望输出xi=sqrt(0.1)*randn(L,1); %白噪声序列%递推估计初值thetaek=zeros(na+nb+d+nc,d);P=106*eye(na+nb+d+nc);for k=1:L time(k)=k; y(k)=-a(2:na+1)*yk(1:na)+b*uk(d:d+nb)+c*xi(k);xik; %采集输出数据 %递推增广最小二乘法 phi

22、e=yk(d:d+ng);uk(d:d+nf);-yek(1:nc); K=P*phie/(1+phie*P*phie); thetae(:,k)=thetaek(:,1)+K*(y(k)-phie*thetaek(:,1); P=(eye(na+nb+d+nc)-K*phie)*P; ye=phie*thetaek(:,d); %预测输出估计值 %提取辨识参数 ge=thetae(1:ng+1,k); fe=thetae(ng+2:ng+nf+2,k); ce=1 thetae(ng+nf+3:ng+nf+2+nc,k); if abs(ce(2)0.9 ce(2)=sign(ce(2)*0

23、.9; end if fe(1)0 yek(1)=ye; yrk(1)=yr(k); xik(1)=xi(k); endendfigure(1);subplot(2,1,1);plot(time,yr(1:L),r:,time,y);xlabel(k); ylabel(h2_r(k)h2(k);legend(h2_r(k),h2(k); axis(0 L 10 60);subplot(2,1,2);plot(time,u);xlabel(k); ylabel(u(k); axis(0 L -40 40);figure(2)subplot(211)plot(1:L,thetae(1:ng+1,:

24、),1:L,thetae(ng+nf+3:ng+2+nf+nc,:);xlabel(k); ylabel(参数估计g、c);legend(g_0,g_1,c_1); axis(0 L -3 4);subplot(212)plot(1:L,thetae(ng+2:ng+2+nf,:);xlabel(k); ylabel(参数估计f);legend(f_0,f_1,f_2,f_3,f_4); axis(0 L 0 4);（7）clear all; close all;am=1.765 -0.7745; bm=0.3064 0.2814;d=6; %对像参数thetam=am;bm; %参考模型参数

25、向量na=length(am); nb=length(bm); %可调参数个数L=700; %仿真长度ypk=zeros(na,1); ymk=zeros(na,1); yrk=zeros(d+nb-1,1); epsilonk=zeros(na,1); %初值，ypk(i)表示yp(k-i)yr=0.6*ones(L/4,1);0.3*ones(L/4,1);0.6*ones(L/4,1);0.3*ones(L/4+d,1); %期望输出G=1*eye(na+nb); G1=0.1*G; D=-am+0.1; thetapr_1=zeros(na+nb,1); %可调参数初值for k=1:

26、L time(k)=k; xm_1=ymk;yrk(d:d+nb-1); %构造参考模型数据向量 ym(k)=thetam*xm_1; %采集参考模型输出 xp_1=ypk;yrk(d:d+nb-1); %构造对象数据向量 v0(k)=ym(k)-thetapr_1*xp_1+D*epsilonk; thetapp(:,k)=G1*xp_1*v0(k)/(1+xp_1*(G+G1)*xp_1); thetapr(:,k)=thetapr_1+G*xp_1*v0(k)/(1+xp_1*(G+G1)*xp_1); thetap(:,k)=thetapr(:,k)+thetapp(:,k); yp(

27、k)=thetap(:,k)*xp_1; %利用最新可调参数计算ypepsilon=ym(k)-yp(k);%更新数据 thetapr_1=thetapr(:,k);for i=d+nb-1:-1:2 yrk(i)=yrk(i-1); end if nb1 yrk(1)=yr(k); end for i=na:-1:2 ypk(i)=ypk(i-1); ymk(i)=ymk(i-1); epsilonk(i)=epsilonk(i-1); end ypk(1)=yp(k); ymk(1)=ym(k); epsilonk(1)=epsilon;endfigure(1)plot(time,thetap(1:na,:);xlabel(k); ylabel(可调参数ap);legend(a_p_1,a_p_2,a_p_3); axis(0 L -2 2.5);figure(2)plot(time,thetap(na+1:na+nb,:);xlabel(k); ylabel(可调参数bp);legend(b_p_0,b_p_1); axis(0 L 0 1);figure(3)plot(time,ym,r,time,yp,:);xlabel(k); ylabel(y_m(k)y_p(k);legend(y_m(k),y_p(k);专心-专注-专业