系统辨识与自适应控制作业.doc

一、系统辨识部分1、SISO系统作为仿真对象二阶的离散的，其中e(k)为服从N(0,1)分布的白噪声序列；输入信号采用四阶逆重复m序列，其中幅值为1，数据信噪比=14.3%选择的辨识模型为：用最小二乘估计的一次性完成算法和LS递推算法分别估计参数，选取数据长度，选取的初始值（遗忘因子=0.995）。解：最小二乘估计的一次性完成算法程序代码：clearclc%-产生M序列输入信号-l=480;y1=1;y2=1;y3=1;y4=0;for i=1:l; x1=xor(y3,y4); x2=y1;x3=y2;x4=y3;y(i)=y4; if y(i)0.143,u(i)=-1; else u(i)=1; end y1=x1;y2=x2;y3=x3;y4=x4;endfigure(1);stem(u)grid ontitle(输入信号)%-产生白噪声信号-A=19;x0=12;M=500;for k=1:l x=A*x0; x1=mod(x,M); v(k)=x1/512; x0=x1;endfigure(2);stem(v)title(白噪声信号)z=zeros(479,1);z(2)=0;z(1)=0;w=0.995;l=477;for k=3:479; z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+v(k); zstar(k)=z(k)*w(l-k+2);endH=zeros(477,4);for k=1:477 H(k,1)=-z(k+1)*w(l-k); H(k,2)=-z(k)*w(l-k); H(k,3)=-u(k+1)*w(l-k); H(k,4)=-u(k)*w(l-k);endestimate=inv(H*H)*H*(zstar(3:479)辨识结果： estimate = -1.5376 0.6938 -0.9780 -0.4565最小二乘估计的递推算法的程序元代码：clearclc%-产生M序列输入信号-l=480;y1=1;y2=1;y3=1;y4=0;for i=1:l; x1=xor(y3,y4); x2=y1;x3=y2;x4=y3;y(i)=y4; if y(i)0.143,u(i)=-1; else u(i)=1; end y1=x1;y2=x2;y3=x3;y4=x4;endfigure(1);stem(u)grid ontitle(输入信号)%-产生白噪声信号-A=19;x0=12;M=500;for k=1:l x=A*x0; x1=mod(x,M); v(k)=x1/512; x0=x1;endfigure(2);stem(v)title(白噪声信号)z=zeros(479,1);z(2)=0;z(1)=0;for k=3:479; z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+v(k);endP=106*eye(4,4);e=zeros(4,478);e(:,1)=P(1,1),P(2,2),P(3,3),P(4,4);c=zeros(4,478);c(:,1)=0.001 0.001 0.001 0.001;K=10;10;10;10;w=0.995;for k=3:479; h=-z(k-1),-z(k-2),u(k-1),u(k-2); K=P*h*inv(h*P*h+w); c(:,k-1)=c(:,k-2)+K*(z(k)-h*c(:,k-2); P=(eye(4)-K*h)*P/w; e(:,k-1)=P(1,1),P(2,2),P(3,3),P(4,4);enda1=c(1,:);a2=c(2,:);b1=c(3,:); b2=c(4,:); ea1=e(1,:);ea2=e(2,:);eb1=e(3,:); eb2=e(4,:); figure(3);i=1:478;plot(i,a1,r,i,a2,y:,i,b1,g,i,b2,:)title(最小二乘递推算法辨识曲线)axis(0,500,-2,2)figure(4);i=1:478;plot(i,ea1,r,i,ea2,:,i,eb1,g,i,eb2,:)title(最小二乘递推算法辨识误差曲线)axis(0,500,0,10)仿真结果：二、自适应控制部分1、设有二阶系统，。（1）按MIT规则设计的自适应控制系统 （2）取，。当分别等于0.1、1.0、3.1、3.5时进行仿真并得出结论； （3）采用修正的MIT，取，当分别等于0.1、1.0、3.1、3.5时进行仿真并得出结论； （4）改变的值，并再对上边的模型进行仿真。解：（1）采用MIT规则设计二阶系统的稳定性分析参考模型的微分算子方程为： (1)并联可调增益系统的微分算子方程： (2) 其数学模型： (3)对(3)中的第一个式子对t求导： (4) (5)为参考模型的输出，参考模型稳定，因此当稳态时有： (6)根据上式可以得到： (7)根据Huwits 判据，可得要使系统稳定的充要条件为： (8)满足式(8)，则系统稳定。 由(8)可知，若输入或自适应增益较大时，系统可能出现不稳定现象。因此，必须对输入信号加以限制，同时自适应增益也不能选择太大。（2）二阶系统的MIT算法模型参考自适应控制系统仿真 图：MIT算法的matlab仿真结构图图：，仿真曲线图：，仿真曲线图：，仿真曲线图：，仿真曲线仿真结论： 由以上几图可以看出，只有满足(8)式的情况下，系统才稳定。在这个仿真系统中，取，。由（8）式可知，只有当yr(t)3.1时系统才是稳定的，而仿真结果也说明了这一点。当自适应增益固定不变时，随着输入信号yr(t)的增大，系统的动态特性得到改善，当满足（8）式时，被控过程出现振荡，再增大，系统出现不稳定现象，输入信号对系统的影响较大，即系统对输入敏感；在选取适当自适应增益的情况下，当输入信号yr(t)值较小时，系统响应时间快，随着yr(t)的增大，响应时间逐渐延长。（3）二阶系统的修正的MIT算法模型参考自适应控制系统仿真 图：修正的MIT算法matlab仿真结构 仿真图形：图：，仿真曲线图：，仿真曲线图：，仿真曲线 图：，仿真曲线仿真结论：由上面的几个仿真图形可以看出，在自适应增益固定时，该自适应调整律在yr(t)变化很大时，依然工作良好。由此可知修正的M IT算法对输入信号变的不敏感了。（4）在不加修正的系统中，当yr(t)=2.5时，分别取值为0.01、0.06、0.1、0.16、0.2进行仿真图：，仿真曲线 图：，仿真曲线 图：，仿真曲线图：，仿真曲线图：，仿真曲线仿真结论：在输入信号取固定的值时，当自适应增益较小时，系统跟踪性能不好，随着的增大，动态特性逐渐改善，当满足（8）式时，即时，被控过程出现等幅振荡，再增大，则系统出现不稳定现象。在修正的系统中，当yr(t)=0.1时，分别取值为0.01、0.1、0.3、1.5、10进行仿真：图：，仿真曲线图：，仿真曲线图：，仿真曲线图：，仿真曲线图：，仿真曲线仿真结论：有上面一组图可知，在输入信号固定的情况下，存在一个适当的自适应增益（不大不小的中间值）使得系统有较好的动态性能。2、采用lyapunov稳定性理论设计图示系统中的自适应律即，和调节规律。yme+ypr(t)_Kc +_f1s+f0自适应律解：系统参考模型的传递函数为 而可调系统的传递函数表示为：显然，只要调整，总可以使被控对象与参考模型的传递函数相等。取系统参数误差向量， 广义误差向量为， 广义误差为，满足微分方程e+a1 -b1e+a0-b0e=Km-Kpr其等价的状态方程为： 其中 取Lyapunov函数为， ，其中P为正定对称阵，对角正定阵，则V是正定的，现对t求导得： 为了使，用以得出自适应律，在上式中可以取：第一项负定；后两项之和为零。所以有：可调系统中的自适应控制律为：为讨论的方便，取 ，,正定矩阵 并对自适应方案进行仿真。利用MATLAB进行仿真，构建自适应控制系统仿真结构框图如下图所示：仿真时，在命令窗口中输入如下的程序：plot(tout,yout)title(阶跃输入下仿真效果图)legend(yr(t),ym,yp,yp0)当输入信号为阶跃信号时，参考模型对阶跃的响应、可调对象对阶跃的响应曲线如下图所示。从图上可看出：未采用自适应控制的，系统输出与参考模型输出有很大的差别，当采用模型参考自适应控制后，系统的输出与参考模型几乎一致。 若给系统输入幅值为0.1，频率为100的正弦信号，其仿真效果图如下图所示：从上图的仿真结果可以看出，随着时间的增加，系统误差逐渐趋近于0.3.设对象用下面的CARMA模型描述：式中是方差为的零均值白噪声，a1,a2,b0,b1,c1,c2未知并有缓慢时变，试设计一个最小方差控制器。解：预测模型： Aq-1=1+a1 q-1+a2q-2 Bq-1=b0 +b1q-1 Cq-1=1+c1 q-1+c2q-2 Eq-1=1+(c1 -a1)q-1 Fq-1=(c2-a2-a1c1+a12)-(a2c1 -a1a2)q-1 Gq-1=b0+(b0c1 -b0a1+b1)q-1+b1(c1-a1)q-2 LSE 模型估计器控制器程序流程图如下：初始化数据Y数据区U数据区 右移一个时刻采样，得，填入数据进行常规控制，将填入数据k=k+1y,u数据区右移一个时刻置初值：按步骤（2）LSE模型公式，构造数据向量采样，得，填入数据按步骤（3）估计器公式，计算 按步骤（4）控制器，计算u(k)并输出，u(k)填入u数据区采样周期到否？YN采样周期到否？YNYNk=0?根据程序流程图和所计算的模型编写MATLABE程序。编程的过程中，设定参考输入的方波信号利用MATLABE命令square直接产生。程序代码如下：clc clear%*设定参考输入周期方波信号*i=200;r=zeros(i,1);t = 0:.1:20;r =square(2*pi*0.03*t); %*设定噪声信号w*w=idinput(i,rgs,0 1,-1 1);%*初始化参数*%预先定维y=zeros(i,1);u=zeros(i,1);x=zeros(1,4);%初始化theta，使其第一列为设定的初始值theta=zeros(4,i);theta(:,1)=2;-10;10;5;%初始化遗忘因子beta和pbeta=0.99;p=2.50*eye(4);T=1;%采样周期K=1;%比比例环节%初始化输入和控制信号y(3)=1;y(2)=2;u(1)=1;u(2)=-1;for m=4:i%构造预测函数 x(1,1)=y(m-1); x(1,2)=y(m-2); x(1,3)=u(m-2); x(1,4)=u(m-3); y(m)=x*theta(:,m-3)+w(m); %利用常规控制手段-PI控制获取输入输出数据 u(m)=K*(r(m)-y(m)+(1/T)*(r(m-1)-y(m-1);end%初始化thetatheta(:,1)=0;0;0;0;for m=4:i %输入LSE模型 u(m-1)=u(m); y(m-1)=y(m); x(1,1)=y(m-1); x(1,2)=y(m-2); x(1,3)=u(m-2); x(1,4)=u(m-3); y(m)=x*theta(:,m-3)+w(m); %利用递推最小二乘法对未知参数进行辨识，并输出系数a1,a2,b0,b1结果 k=p*x*inv(beta+x*p*x); theta(:,(m-1)=theta(:,(m-2)+k*(y(m)-x*theta(:,(m-2); p=(1/beta)*p-k*x*p; %计算控制率uu(m)=(r(m)-theta(1,m-1)*y(m)-theta(2,m-1)*y(m-1)*inv(theta(3,m-1)*u(m-2)+theta(4,m-1)*u(m-3);endfigure;plot(-theta(1,3:i-1),-);hold on;plot(-theta(2,3:i-1),-.k);hold on;plot(theta(3,3:i-1),:r);hold onplot(theta(4,3:i-1),-c)legend(a1,a2,b0,b1);title(系数辨识);各未知系数的仿真结果如下图所示：系统中存在有噪声信号的影响，仿真的过程中噪声信号具有随机性，所以每次仿真的所得出的图形不尽相同，但只要最后未知参数达到一稳定值则表明实现了最小方差控制器的功能，使系统达到了稳定。4.对于系统模型求系统的一个广义最小方差控制律。解 ： A(q-1)= (1-1.6q-1+0.6q-2)，C(q-1)=(1-0.4q-1)，B(q-1)=1.5-0.53 q-1+0.9 q-2 , d=1. 将A,C,d代入式 中 (1-0.4q-1)= (1-1.6q-1+0.6q-2)E(q-1)+q-1D(q-1)可解得 E(q-1)=1， D(q-1)=1.2-0.6q-