k c t t实验一一、 可增益Lyapunov-MRAC 算法1.1 步骤：已知：N （s ）D (s )第一步：选择参考模型，即Gm （s ）；第二步：选择输入信号 y r （t ）和自适应增益γ；第三步：采样当前参考模型输出 y m （t ）和系统实际输出 y p （t ）；第四步：利用公式 &( )= γe (t ) y r (t ) 和公式 u ( )= k c (t ) y r (t ) ；第五步：t t+h ，返回第三步，继续循环。

1.2 考虑如下被控方对象模型：G p (s )=选择参考模型为：k p (s + 1) s 2 + 5s + 1 ， k p 未知（仿真时取 k p =1）G m (s )= k m (s + 1) s 2 + 5s + 1， k m =1因为 G P (s )、 G m (s )均为严格正实函数。

取自适应增益γ=0.2，输入 y r 为 方波信号，幅值r=1，采用可调增益Lyapunov-MRAC 算法，仿真程序以及仿真结 果如下。

二、仿真程序%可调增益Lyapunov-MRACclear all ;close all ;h=0.1;L=100/h;%数值积分步长和仿真步数num=[1 1];den=[1 5 1];n=length(den)-1; kp=1;[Ap,Bp,Cp,Dp]=tf2ss(kp*num,den); km=1;[Am,Bm,Cm,Dm]=tf2ss(km*num,den); gamma=0.2;yr0=0;u0=0;e0=0;xp0=zeros(n,1);xm0=zeros(n,1); kc0=0;r=1;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;y m (t )、y p (t )u (t )xp(:,k)=xp0+h*(Ap*xp0+Bp*u0); yp(k)=Cp*xp(:,k);xm(:,k)=xm0+h*(Am*xm0+Bm*yr0); ym(k)=Cm*xm(:,k); e(k)=ym(k)-yp(k);kc=kc0+h*gamma*e0*yr0;%Lyapunov-MRAC 自适应律 u(k)=kc*yr(k); %更新数据yr0=yr(k);u0=u(k);e0=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)');三、仿真结果21 0-121-1-2y m (t) y p (t)0 10 20 30 40 50 60 70 80 90 100tt第一步：设置初值 和 ,输入初始数据；( )0θ( )0P∧ 参数 ,即 、 和 ；θGf C 取初值P(0)=I 、 =0.01，的下界为 =0.1，期望输出为610( )0θf min f( )ky r 第三步： 构造观测数据向量 ，利用递推算法在线实时估计控制器()dk -ϕ 实验二一、最小方差直接自校正控制算法1.1 步骤：已知：模型阶次 n a 、 n b 、 n c 及纯延时d 。

∧第二步：采用当前实际输出 y (k ) 和期望输出 y r (k + d ) ； ∧∧ ∧ ∧第四步：计算并实施 u (k )；第五步：k k+h ，返回第二步，继续循环。

1.2 考虑如下被控方对象模型：y(k)-1.7y(k-1)+0.7y(k-2)=u(k-4)+0.5u(k-5)+ξ(k)+0.2ξ(k-1)式中，ξ(k)为方差为0.1的白噪声。

∧ ∧幅值为10的方波信号，采用最小方差直接自校正控制算法，其控制效果如图 1、2.二、仿真程序%最小方差直接自校正控制a=[1 -1.7 0.7];b=[1 0.5];c=[1 0.2];d=4;na=length(a)-1;nb=length(b)-1;nc=length(c)-1; nf=nb+d-1;ng=na-1; L=400;uk=zeros(d+nf,1); yk=zeros(d+ng,1); yek=zeros(nc,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); %递推估计初值thetaek=zeros(na+nb+d+nc,d);P=10^6*eye(na+nb+d+nc);for k=1:Ltime(k)=k;y(k)=-a(2:na+1)*yk(1:na)+b*uk(d:d+nb)+c*[xi(k);xik];%递推增广最小二乘法phie=[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.9ce(2)=sign(ce(2))*0.9;endif fe(1)<0.1fe(1)=0.1;endu(k)=(-fe(2:nf+1)*uk(1:nf)+ce*[yr(k+d:-1:k+d-min(d,nc));yrk(1:nc-d)]-ge*[y(k);yk(1:na-1)])/fe(1);%控制量%更新数据for i=d:-1:2thetaek(:,i)=thetaek(:,i-1);endthetaek(:,1)=thetae(:,k);for i=d+nf:-1:2;uk(i)=uk(i-1);enduk(1)=u(k);for i=d+ng:-1:2yk(i)=yk(i-1);endyk(1)=y(k);for i=nc:-1:2;yek(i)=yek(i-1);u (k )yrk(i)=yrk(i-1); xik(i)=xik(i-1); endif nc>0yek(1)=ye; yrk(1)=yr(k); xik(1)=xi(k); end endfigure(1);subplot(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)');axis([0 L -20 20]); subplot(2,1,2); plot(time,u);xlabel('k');ylabel('u(k)');axis([0 L -40 40]); figure(2);subplot(2,1,1);plot([1:L],thetae(1:ng+1,:),[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(2,1,2);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]);三、仿真结果20 10-1-204020-2-40y (k)ry (k)050100150200250300350400k050100150200250300350400k图1直接算法控制结果、、、、g 、c 、、、、f42g 0 g 1 c 1-2 0 50100150200 250300350400k4321f 0 f 1 f 2 f 3 f 40 50 100 150 200 250 300 350 400k图2 控制器参数估计结果。