自适应控制结课作业

1、自适应控制结课作业2016年1月目录1 遗忘因子递推最小二乘法11.1最小二乘理论11.2带遗忘因子的递推最小二乘法11.2.1白噪声与白噪声序列11.2.2遗忘因子递推最小二乘法21.3仿真实例4程序1 带遗忘因子的递推最小二乘法：62 广义最小方差自校正控制82.1广义最小方差自校正控制82.1.1 逆M序列82.1.2 广义最小方差自校正控制82.2仿真实例10程序2 广义最小方差自校正控制：133 参考模型自适应控制163.1参考模型自适应控制163.2仿真实例183.2.1数值积分183.2.2仿真结果18程序3 参考模型自适应控制：2321 遗忘因子递推最小二乘法1.1最小二乘理论

2、最小二乘最早的想法是高斯在1795年预测行星和彗星运动轨道时提出来的，“未知量的最大可能的值是这样一个数值，它使各次实际观测和计算值之间的差值的平方乘以度量其精确度的数值以后的和为最小”1。这一估计方法原理简单，不需要随机变量的任何统计特性，目前已经成为动态系统辨识的主要手段。最小二乘辨识方法使其能得到一个在最小方差意义上与实验数据最好拟合的数学模型。由最小二乘法获得的估计在一定条件下有最佳的统计特性，即统计结果是无偏的、一致的和有效的。1.2带遗忘因子的递推最小二乘法1.2.1白噪声与白噪声序列系统辨识中所用到的数据通常含有噪声。从工程实际出发，这种噪声往往可以视为具有理想谱密度的平稳随机过

3、程。白噪声是一种最简单的随机过程，是由一系列不相关的随机变量组成的理想化随机过程。白噪声的数学描述如下：如果随机过程均值为0，自相关函数为，即式中，为单位脉冲函数(亦称为Dirac函数)，即，且则称该随机过程为白噪声，其离散形式是白噪声序列。如果随机序列均值为零，且两两互不相关，即对应的相关函数为：则这种随机序列称为白噪声序列。其谱密度函数为常数。白噪声序列的功率在到的全频段内均匀分布2。建立系统的数学模型时，如果模型结构正确，则模型参数辨识的精度将直接依赖于输入信号，因此合理选用辨识输入信号是保证能否获得理想的辨识结果的关键之一。理论分析表明，白噪声作为被辨识系统的输入时，可以激发系统的所有

4、模态，可对系统充分激励，可防止数据病态，保证辨识精度，可以保证获得较好的辨识效果。图1-1 白噪声序列1.2.2遗忘因子递推最小二乘法假设被辨识的系统为一单入单出的离散时间系统，且已知为CAR模型，如图2所示：(k)u(k)z-dB(z-1)A(z-1)y(k)图1-2 辨识系统模型即式中，为输入变量，为输出变量，为白噪声，且则上式可转化为如下最小二乘格式：式中，为数据向量，待估参数向量，且取算法的性能指标为式中，为遗忘因子（）。带遗忘因子的递推最小二乘估计的算法公式为：公式表明，新的参数估计是用新的实际测量值与基于老模型进行预测得到的量之偏差，对前面的参数估计加以修正得到的，修正系数阵为。的

5、物理意义是参数估计误差的方差，作为参数估计精度的一种度量。遗忘因子的作用是削弱过去数据的作用，通常选择0.95到0.998之间的数。带遗忘因子的递推最小二乘估计算法属于在线辨识所用方法的一种，它既能克服离线辨识的缺点，也能克服递推最小二乘估计中的“数据饱和”现象。遗忘因子最小二乘法的算法：已知式阶次、式阶次以及延迟。步骤1：设置初值和及遗忘因子，输入初始数据；步骤2：采样当前输出和输入；步骤3：利用递推公式，计算、和；步骤4：，返回步骤2，继续循环。1.3仿真实例系统模型如下：其中，输入为方差为1的白噪声，为方差为0.1的白噪声。由于和的选择可按如下方法：式中，为充分大的正实数，为零向量或充分

6、小的正的实向量。因此，取初值、。仿真结果如下：图 1-3 参数估计结果图 1-4实际输出与辨识输出对比图 1-5 实际输出与辨识输出对比（局部放大）附录程序1 带遗忘因子的递推最小二乘法：clear all; close all;a=1 1.2 0.8 0.1' b=1 0.5 2 1.5' d=3;na=length(a)-1; nb=length(b)-1;L=1000;uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i)yk=zeros(na,1);%输出初值ymk=zeros(na,1); %辨识输出的初值u=randn(L,1); %输入采用白噪

7、声序列xi=sqrt(0.1)*randn(L,1); %白噪声序列thetaek=zeros(na+nb+1,1); %thetae初值P=106*eye(na+nb+1);lambda=0.998;for k=1:L theta(:,k)=a(2:na+1);b;%对象参数真实值 phi=-yk;uk(d:d+nb); phim=-ymk;uk(d:d+nb); y(k)=phi'*theta(:,k)+xi(k); ym(k)=phim'*thetaek; %递推公式： K=P*phi/(lambda+phi'*P*phi); thetae(:,k)=thetae

8、k+K*(y(k)-phi'*thetaek); P=(eye(na+nb+1)-K*phi')*P/lambda; thetaek=thetae(:,k); %产生新的phi和phim 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); ymk(i)=ymk(i-1); end ymk(1)=ym(k); yk(1)=y(k);endfigure(1)plot(1:L,thetae(1:na+nb+1,:); hold on; plot(1:L,theta(1:na+nb+1,

9、:),'k:');legend('a_1','a_2','a_3','b_0','b_1','b_2','b_3'); axis(0 L -0.5 2.5);%x轴范围，y轴范围figure(2)plot(1:L,y);hold onplot(1:L,ym,'-r');legend('y','ym');axis(0 L -30 30);2 带有辅助变量的最小方差控制器2.1带有辅助变量的最小方差控制器2.1.1 逆M序列

10、设M(k)是周期为Np bit、元素取值为0或1的M序列，S(k)是周期为2bit、元素依次取值为0或1的方波序列，将这两个序列按位进行异或运算，得到的复合序列就是周期为2Np bit、元素取值为0或1的逆M序列，记作IM(k)，即有IM(k)=M(k)S(k)将上述逆M序列的逻辑值“0”或“1”分别变换为1或1，此时逆M序列均值为0.虽然逆M序列只是M序列与方波序列简单复合的结果，但其性质却优于M序列，使其在辨识领域有着更为广泛的应用。2.1.2 广义最小方差自校正控制器当考虑干扰对系统的作用时，控制器的设计就是要最大限度的减小干扰对系统的影响。鉴于一般被控对象或过程都存在不同程度的纯迟延，

11、控制对系统的作用要到时刻才有响应。在这段纯迟延的时间内，干扰仍会作用于系统，所以在时刻预测时刻的输出，并按照预测误差的方差最小的原则，设计现时控制，并加以实施。当过程参数未知，或者时变时，用递推最小二乘法估计，或者直接估计控制器参数，然后算出控制量来，这就是最小方差自校正控制的基本思想。 而广义最小方差自校正控制就是在性能指标函数中引入加权多项式。被控对象为：其中：为延迟因子，为输入变量，为输出变量，为白噪声。为Hurwitz多项式。选择性能指标函数为：式中，为期望输出；为第拍的输出；为第拍的控制；、和分别为输出、期望输出和控制的加权多项式，它们分别具有改善闭环系统性能，软化输入和约束控制量的

12、作用。并且上述多项式的阶次及参数根据实际需要确定。由此，据文献1知，广义最小方差控制律为：uk=Cz-1Rz-1yrk+d-G(z-1)P(z-1)y(k)q0b0Cz-1Qz-1+F(z-1)P(z-1) (2-1)在进行控制系统设计时，一般可以取加权多项式、和，而大小的选取需要在快速性和稳定性方面进行权衡。当被控对象的参数未知时，也可以利用递推算法直接估计广义最小方差控制的参数，即广义最小方差自校正控制直接算法。由式21知，在选定加权多项式P(z-1)、Rz-1、Qz-1后，只需再确定多项式G(z-1)、F(z-1)、Cz-1，即可确定广义最小方差控制律u(k)。控制器参数的估计模型及递推

13、公式如下：k=k-1+Kkyk-Tk-dk-1Kk=PK-1k-d+Tk-dPk-1k-dPk=1I-K(k)Tk-dP(k-1) (2-2)式中k-dyk-d,yk-d-ng,uk-d,uk-d-nf,-y*k-1,-y*k-ncTRng+nf+nc+2×1=g0,gng,f0,fnf,c1,cncTRng+nf+nc+2×1y*k=Tk-dk-d或y*k=yrk（2-3）由式21得广义最小方差控制律为uk=1q02f0+f0q0f0q0-CQ+f0-FPuk+CRyrk+d-GPy(k) (2-4)由上式可知，在控制算法实施过程中，如果f0趋于零，则会出现零除现象。为此

14、，应对f0的最小值加以约束，这就意味着应该事先知道f0的符号和下界，或事先确定f0的值。2.1.3 带有辅助变量的最小方差控制器广义最小方差尚需解决偏差问题，因为采用广义最小方差性能指标中含有u2(k)，稳态时必将使y2(k+d)产生某一数值，来补偿u2(k)的值，而在稳态时y(k+d)0产生了稳态偏差，而且不便于处理设定值问题。定义一种辅助系统的广义输出为：xk=yk-yrk-d+uk-d (2-5)求使性能指标函数：J=Ez2(k+d)为最小的控制u(k)为：uk=-1g0Hz-1yrk-Mz-1yk+i=1nggiuk-i (2-6)采用带遗k=k-1+Kkzk-Tk-dk-1Kk=PK

15、-1k-d+Tk-dPk-1k-dPk=1I-K(k)Tk-dP(k-1) (2-7)式中：k-dyk-d,yk-d-nm,uk-d,uk-d-ng,-yrk-d,-yrk-d-nhTRnm+ng+nh+2×1=m0,mnm,g0,gng,h0,hnhTRnm+ng+nh+2×1（2-8）带有辅助变量的最小方差控制器的算法：已知：模型阶次、以及延迟。步骤1：设置初值和，输入初始数据，并设置加权多项式、；步骤2：采样当前实际输出和期望输出yr(k)；步骤3：根据式28构造观测数据向量k-d，利用遗忘因子递推最小二乘法递推公式27在线实时估计控制器参数，即m、G和h;步骤4：利

16、用式26计算并实施u(k);步骤5：返回步骤2（），继续循环。2.2仿真实例设系统模型如下：yk-1.8yk-1+0.7yk-2+0.1yk-3=uk-4+0.2uk-5+k+0.1(k-1)其中，为取值为0.1的逆M序列，采用带有辅助变量的最小方差控制器。取初值、；设置加权多项式、Rz-1=1、2。期望输出采用幅值为10的方波型号，其控制结果如下：图21期望输出与实际输出对比图22控制量u(k)图23 参数辨识附录程序2 带有辅助变量的最小方差控制器：clear all;close all;clca=1 -1.8 0.7 0.1;b=1 1.8;c=1 0.2;d=4;na=length(a

17、)-1; nb=length(b)-1; nc=length(c)-1;nf=nb+d-1; ng=na-1; Pw=1;R=1;Q=2;np=length(Pw)-1; nh=length(R)-1+nc; nq=length(Q)-1;nm=np+ng;nG=max(np+nf,nc+nq);L=400;uk=zeros(d+nG,1);yk=zeros(d+nm,1);yrk=zeros(d+nh,1);xmk=zeros(nc,1);yr=10*ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4+d,1);x1=1;x2=1;x3=1;x4=0

18、;S=1;for k=1:L IM=xor(S,x4); if IM=0 xm(k)=-0.1; else xm(k)=0.1; end S=not(S); M(k)=xor(x3,x4); x4=x3;x3=x2;x2=x1;x1=M(k);endthetae_1=zeros(nm+nG+nh+3,d);P=106*eye(nm+nG+nh+3);lambda=0.998;for k=1:L time(k)=k; y(k)=-a(2:na+1)*yk(1:na)+b*uk(d:d+nb)+c*xm(k);xmk; phi=yk(d:d+nm);uk(d:d+nG);-yrk(d:nh+d);

19、 z(k)=Pw*y(k)-R*yrk(d)+Q/b(1)*Q*uk(d); K=P*phi/(lambda+phi'*P*phi); thetae(:,k)=thetae_1(:,1)+K*(z(k)-phi'*thetae_1(:,1); P=(eye(nm+nG+nh+3)-K*phi')*P/lambda; me=thetae(1:nm+1,k)'ge=thetae(nm+2:nm+nG+2,k)'he=thetae(nm+nG+3:nm+nG+3+nh,k)' if ge(1)<0.1 ge(1)=0.1; end u(k)=(h

20、e(1)*yr(k)+he(2:nh+1)*yrk(1:nh)-me(1)*y(k)-me(2:nm+1)*yk(1:nm)-ge(2:nG+1)*uk(1:nG)/ge(1); thetae; for i=d:-1:2 thetae_1(i)=thetae_1(i-1); end thetae_1(:,1)=thetae(:,k); for i=nG+d:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=nm+d:-1:2 yk(i)=yk(i-1); end yk(1)=y(k); for i=d+nh:-1:2 yrk(i)=yrk(i-1); end

21、yrk(1)=yr(k); for i=nc:-1:2 xmk(i)=xmk(i-1); end xmk(1)=xm(k);endfigure(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); figure(2)plot(time,u);xlabel('k');ylabel('u(k)');ax

22、is(0 L -10 10); figure(3)plot(1:L,thetae(1:nm+1,:),1:L,thetae(nm+2:nm+nG+2,:),1:L,thetae(nm+nG+3:nm+nG+3+nh,:);xlabel('k');ylabel('参数估计m,G,h');legend('m_0','m_1','m_2','G_0','G_1','G_2','G_3','G_4','h_0','h_1

23、')axis(0 L -3 10);3 参考模型自适应控制3.1参考模型自适应控制模型参考自适应控制器（MRAC，model reference adaptive control），即为利用可调系统（包含被控对象）的各种信息，度量或测出某种性能指标，把它与参考模型期望的性能指标相比较；用性能指标偏差（广义误差）通过非线性反馈的自适应律来调节可调系统，以削弱可调系统因“不确定性”所造成的性能指标的偏差，最后达到使被控的可调系统获得较好的性能指标的目的。MARC中的可调系统一般包括被控对象和调节器，它们形成一常规的反馈控制系统。 这个系统相对于MARC 系统来说是一个子系统或称“内回路”。

24、另外， MARC系统还有一个自适应反馈回路，称为外回路，它用来调节可调系统。有内外回路组成双回路系统是MARC 系统的结构特点5。模型参考自适应控制可以处理缓慢变化的不确定性对象的控制问题。 它由于可以不必经过系统辨识而度量性能指标，因而有可能或得快速跟踪控制。由于被控对象的全部状态要准确得到很困难， 按被控对象输入和输出直接设计自适应控制系统更有价值，一般有直接法和间接法。所谓直接和间接，指的是对未知的被控对象进行直接控制和间接控制。 间接控制的基本思想是用未知的被控对象的输入输出数据来估计被控对象的参数， 并用这些参数估计值产生一个反馈函数去调整调节器参数。直接控制和间接控制不同，在产生反

25、馈控制信号之前没有明显的被控对象的辨识。所以二者之间主要不同在于：在直接控制中要有一个显式的理想特性的参考模型， 而间接控制则需要被控对象模型进行在线辨识并用隐式方法去产生自适应律。此外，间接控制是用辨识误差，而直接控制则用控制误差去修改调节器参数6。现用一种直接法，即 K.S.Narendra 提出的稳定性自适应控制方案。原理框图如图3-1所示：图3-1 模型参考自适应控制原理框图由于系统的对象的分母阶次比分子大两阶，自适应控制器可设计如下。(1)选择，并构造辅助信号状态方程(2)可调参数自适应律为：式中，为正定矩阵。(3)自适应控制律为：参考模型自适应控制的算法：已知：被控对象的阶数、。步

26、骤1：选择参考模型为严格正实、稳定最小相位系统，与阶数及相对阶相同，并具有理想的动态性能；并利用构造辅助信号发生器状态矩阵；步骤2：设置初值，选择自适应增益矩阵和输入信号，并初始化数据；步骤3：采样当前参考模型输出和系统实际输出，并计算；步骤4：辅助信号状态方程计算和；步骤5：利用自适应律计算；步骤6：组件，并有自适应控制律计算；步骤7：，返回步骤3，继续循环。3.2仿真实例3.2.1数值积分对于如下连续系统的一阶微分方程：在离散系统中可采用欧拉法求解，其递推公式如下式中，称为计算步长或者步距。该方法简单、计算量小，由前一点即可推出后一点的值，属于单步法。适当减小计算步长h有助于提高计算精度。

27、3.2.2仿真结果采用的对象模型为模型静态增益为1，故采用的参考模型为改模型是严格正实的最小相位系统，并且具有良好的动态性能。取自适应增益矩阵，输入信号为方波信号，使用模型参考自适应控制后的仿真结果如下：图 3-1 实际输出与参考输出对比(左图无噪声，右图有噪声) 图 3-2 偏差值e(t) 图 3-3 控制量u(t)图 3-4 可调参数kc、c1、c2、c3曲线 图 3-5 可调参数d0、d1、d2、d3曲线附录程序3 参考模型自适应控制：clear all; close all;h=0.1; L=800/h; F=1 2;nump=1 4 4; denp=1 8 22 24 9; Ap,B

28、p,Cp,Dp=tf2ss(nump,denp); n=length(denp)-1; %对象参数 xi=sqrt(0.01)*randn(L,1); %白噪声 numm=3.375 10.125 6.75; denm=1 8 21.75 22.5 6.75; Am,Bm,Cm,Dm=tf2ss(numm,denm); Df=conv(F,numm); %辅助信号发生器传递函数分母Af=zeros(n-2,1),eye(n-2);-Df(n:-1:2); %辅助信号发生器状态矩阵Bf=zeros(n-2,1);1; %辅助信号发生器输入矩阵 yr0=0; yp0=0; u0=0; e0=0;

29、%初值v10=zeros(n-1,1); v20=zeros(n-1,1); %辅助信号发生器状态初值xp0=zeros(n,1); xm0=zeros(n,1); %状态向量初值theta0=zeros(2*n,1); %可调参数向量初值zeta0=zeros(2*n,1); yr=ones(1,L/80) -ones(1,L/80);for i=1:39 yr=yr ones(1,L/80) -ones(1,L/80);end Gamma=diag(0.8 1 0.2 0.2 0.3 0.1 0.2 10); %自适应增益矩阵 for k=1:L time(k)=k*h; xp(:,k)=xp0+h*(Ap*xp0+Bp*u0); yp(k)=Cp*xp(:,k)+Dp*u0; yp(k)=yp(k)+xi(k); %加入白噪声 xm(:,k)=xm0+h*(Am*xm0+Bm*yr0); ym(k)=Cm*xm(:,k)+