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

一、系统辨识部分

1、SISO系统作为仿真对象

Z(k)?1.5Z(k?1)?0.7Z(k?2)?U(k?1)?0.5U(k?2)?e(k)二阶的离散的，

其中{e(k)}为服从N(0,1)分布的白噪声序列；

输入信号U(k)采用四阶逆重复m序列，其中幅值为1，数据信噪比β=14.3% 选择的辨识模型为：

Z(k)?a1Z(k?1)?a2Z(k?2)?b1U(k?1)?b20.5U(k?2)?ε(k)

用最小二乘估计的一次性完成算法和LS递推算法分别估计参数，选取数

?Q0?0.001据长度l?480，选取的初始值?（遗忘因子?=0.995）。 6p?10I2?2?0

解：

最小二乘估计的一次性完成算法程序代码：

clear clc

%-----产生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; end

figure(1);stem(u) grid on

title('输入信号')

%-----产生白噪声信号------- A=19;x0=12;M=500; for k=1:l x=A*x0;

x1=mod(x,M); v(k)=x1/512; x0=x1; end

figure(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); end

H=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); end

estimate=inv(H'*H)*H'*(zstar(3:479))' 辨识结果： estimate =

-1.5376 0.6938 -0.9780 -0.4565

最小二乘估计的递推算法的程序元代码： clear clc

%-----产生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; end

figure(1);stem(u) grid on

title('输入信号')

%-----产生白噪声信号------- A=19;x0=12;M=500; for k=1:l

x=A*x0;

x1=mod(x,M); v(k)=x1/512; x0=x1; end

figure(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); end

P=10^6*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)]; end

a1=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])

仿真结果：

输入信号10.80.60.40.20-0.2-0.4-0.6-0.8-1050100150200250300350400450500

白噪声信号10.90.80.70.60.50.40.30.20.10050100150200250300350400450500

最小二乘递推算法辨识曲线21.510.50-0.5-1-1.5-2050100150200250300350400450500

最小二乘递推算法辨识误差曲线109876543210050100150200250300350400450500

二、自适应控制部分

1、设有二阶系统，D(s)?a2s2?a1s?1，N(s)?1。 （1）按MIT规则设计的自适应控制系统

（2）取a1?a2?1，Km?1，??0.1。当yr(t)分别等于0.1、1.0、3.1、3.5时进行仿真并得出结论；

（3）采用修正的MIT，取a1?a2?1，Km?1，??0.1，??2，??0.01，当yr(t)分别等于0.1、1.0、3.1、3.5时进行仿真并得出结论； （4）改变?的值，并再对上边的模型进行仿真。 解：

（1）采用MIT规则设计二阶系统的稳定性分析 参考模型的微分算子方程为：

(a2p2?a1p?1)ym(t)?kmyr(t) (1)

并联可调增益系统的微分算子方程：

(a2p2?a1p?1)yp(t)?kckvyr(t) (2)

其数学模型：

?????a2e?a1e?e?(km?kckv)yr(t)?????a2ym?a1ym?ym?kmyr(t) (3) ???k?c??eym对(3)中的第一个式子对t求导：

a??????2e?a1e?e???eym(t)kvyr(t) (4) a??????2e?a1e?e??eym(t)kvyr(t)?0 (5)

ym(t)为参考模型的输出，参考模型稳定，因此当稳态时有：

ym?kvyr(t) (6)

根据上式可以得到：

a??????22e?a1e?e??ekmkvyr(t)?0 (7)

根据Huwits 判据，可得要使系统稳定的充要条件为：

2??a2?0,a1?0,?kmkvyr(t)?0 ?2??a1?1?a2?kmkvyr(t)(8)

满足式(8)，则系统稳定。

由(8)可知，若输入yr或自适应增益?较大时，系统可能出现不稳定现象。因此，必须对输入信号yr加以限制，同时自适应增益?也不能选择太大。

（2）二阶系统的MIT算法模型参考自适应控制系统仿真

图：MIT算法的matlab仿真结构图

0.120.10.080.060.040.020-0.02-0.04-0.06ypymyrye05101520253035404550

图：??0.1，yr(t)?0.1仿真曲线

1.2ypymyrye10.80.60.40.20-0.205101520253035404550

图：??0.1，yr(t)?1仿真曲线

6ypymyrye543210-105101520253035404550

图：??0.1，yr(t)?3.1仿真曲线

40ypymyrye3020100-10-20-3005101520253035404550

图：??0.1，yr(t)?3.5仿真曲线

仿真结论：

由以上几图可以看出，只有满足(8)式的情况下，系统才稳定。在这个仿真系统中，取a1?a2?1，Km?1，Kv?1??0.1。由（8）式可知，只有当

yr(t)<3.1时系统才是稳定的，而仿真结果也说明了这一点。当自适应增益?固定不变时，随着输入信号yr(t)的增大，系统的动态特性得到改善，当满足（8）式时，被控过程出现振荡，再增大，系统出现不稳定现象，输入信号对系统的影响较大，即系统对输入敏感；在选取适当自适应增益?的情况下，当输入信号yr(t)值较小时，系统响应时间快，随着yr(t)的增大，响应时间逐渐延长。 （3）二阶系统的修正的MIT算法模型参考自适应控制系统仿真

0.01a1km12s +s+1Transfer Fcneymym*ymAdd2e*ym0.015Band-LimitedWhite NoiseStepu(1)*u(2)Gain20.1FcnIntegratorSaturation1sb*xxDivide2bGain11AddGain12s +s+1Transfer Fcn1Add1ypym1Out1yr(t)

图：修正的MIT算法matlab仿真结构

仿真图形：

0.120.10.080.060.040.020-0.02-0.04-0.06ypymyrye05101520253035404550

图：??0.1，yr(t)?0.1仿真曲线

1.2ypymyrye10.80.60.40.20-0.205101520253035404550

图：??0.1，yr(t)?1仿真曲线

3ypymyrye2.521.510.50-0.505101520253035404550

图：??0.1，yr(t)?2.5仿真曲线

4.543.532.521.510.50-0.505101520253035404550ypymyrye

图：??0.1，yr(t)?3.5仿真曲线

仿真结论：

由上面的几个仿真图形可以看出，在自适应增益?固定时，该自适应调整律在yr(t)变化很大时，依然工作良好。由此可知修正的M IT算法对输入信号

yr(t)变的不敏感了。

（4）在不加修正的系统中，

当yr(t)=2.5时，分别取?值为0.01、0.06、0.1、0.16、0.2进行仿真

3ypymyrye2.521.510.50-0.505101520253035404550

图：??0.01，yr(t)?2.5仿真曲线

3.532.521.510.50-0.5ypymyrye05101520253035404550

图：??0.06，yr(t)?2.5仿真曲线

43.532.521.510.50-0.5ypymyrye05101520253035404550

图：??0.1，yr(t)?2.5仿真曲线

5ypymyrye43210-105101520253035404550

图：??0.16，yr(t)?2.5仿真曲线

302520151050-5-10-15-2005101520253035404550ypymyrye

图：??0.2，yr(t)?2.5仿真曲线

仿真结论：

在输入信号yr取固定的值时，当自适应增益?较小时，系统跟踪性能不好，随着?的增大，动态特性逐渐改善，当满足（8）式时，即??0.16时，被控过程出现等幅振荡，再增大，则系统出现不稳定现象。

在修正的系统中，

当yr(t)=0.1时，分别取?值为0.01、0.1、0.3、1.5、10进行仿真：

0.120.10.080.060.040.020-0.02-0.04-0.06ypymyrye05101520253035404550

图：??0.01，yr(t)?0.1仿真曲线

0.120.10.080.060.040.020-0.02-0.04-0.06ypymyrye05101520253035404550

图：??0.1，yr(t)?0.1仿真曲线

0.120.10.080.060.040.020-0.02-0.04-0.06ypymyrye05101520253035404550

图：??0.3，yr(t)?0.1仿真曲线

0.25ypymyrye0.20.150.10.050-0.05-0.105101520253035404550

图：??1.5，yr(t)?0.1仿真曲线

0.8ypymyrye0.60.40.20-0.2-0.4-0.605101520253035404550

图：??10，yr(t)?0.1仿真曲线

仿真结论：

有上面一组图可知，在输入信号yr(t)固定的情况下，存在一个适当的自适应增益?（不大不小的中间值）使得系统有较好的动态性能。

2、采用lyapunov稳定性理论设计图示系统中的自适应律即f0，f1和kc调节规律。

ym + e + r(t) Kms2?a1s?a0 yp _ Kp Kc s2?b1s?b0 + 解：

f1s+f0 自适应律 Km系统参考模型的传递函数为 s2?as?a10而可调系统的传递函数表示为：

2s??b1?Kpf1?s??b0?Kpf0?KpKc?Kp?s2?b1?s?b0?T显然，只要调整Kc，f1，f0，，总可以使被控对象与参考模型的传递函数相等。取系统参数误差向量???b0??a0b1??a1Kp??Km?， 广义误差向量为，

??T??? e??ee???

广义误差为e?ym?yp，满足微分方程

e+ a1 ?b1′ e+ a0?b0′ e= Km?Kp′ r

其等价的状态方程为：e?Ae?Br 其中

?

?01?A???????b0?a0b1?a1???0?B??????Km?Kp??12??diag???0?1???对角正定阵，则V是正定的，现对t求导得：

T??T??dV1??TT??ePe?ePe?θΛθ?θΛθ?dt2????1?21?2?TTV??ePe?????取Lyapunov函数为， ，其中P为正定对称阵，??

?T?TT?TTeAP?PAe?BrPe?ePBr+2θΛθ???????????T?TT?TeAP?PAe?2BrPe+2θΛθ??????为了使V?0，用以得出自适应律，在上式中可以取：

第一项eT?APT?PAe负定；

T?后两项之和2所以有：

?Br?Pe+2θΛθ为零。

T?? BrTPe??erP?erP?K?K????1222?mp?? ????TθΛθ?Kp?0f0b0??a0?Kp?1f1b1??a1?Kp?KcKp??Km

????????可调系统中Kc,f1,f0的自适应控制律为：

t???ep?epr?1222???d??K0c??Kc????0?Kpf0?t0???ep?epyp?1222???d??f00???0Kpf1?t0???ep?epyp?1222???d??f01???1Kp为讨论的方便，取 ???diag???0?1????diag??111，?,正定矩阵

P?? ?21???12?a1?3,b1?2,a0?Km?5,b0?Kp?3;f1?0??f0?0??0;Kc(0)?1并对自适应方案进行仿真。

利用MATLAB进行仿真，构建自适应控制系统仿真结构框图如下图所示：

仿真时，在命令窗口中输入如下的程序： plot(tout,yout)

title('阶跃输入下仿真效果图') legend('yr(t)','ym','yp','yp0')

当输入信号为阶跃信号时，参考模型对阶跃的响应、可调对象对阶跃的响应曲线如下图所示。从图上可看出：未采用自适应控制的，系统输出与参考模型输出有很大的差别，当采用模型参考自适应控制后，系统的输出与参考模型几乎一致。

阶跃输入下仿真效果图1.4yr(t)ymypyp01.210.80.60.40.20012345678910

若给系统输入幅值为0.1，频率为100的正弦信号，其仿真效果图如下图所示：

12x 10-4正弦输入下仿真效果图ymypyp01086420-2012345678910

从上图的仿真结果可以看出，随着时间的增加，系统误差逐渐趋近于0.

3.设对象用下面的CARMA模型描述：

y(k)?a1y(k?1)?a2y(k?2)?b0u(k?2)?bu1(k?3)?ew(k)?c1ew(k?1)?c2ew(k?2) 式中ew(k)是方差为?w的零均值白噪声，a1,a2,b0,b1,c1,c2未知并有缓慢时变，试设计一个最小方差控制器。 解：

①预测模型：

2m?2,nA?2,nB?1,nF?nA?1?1,nG?nB?m?1?2

A q?1 =1+a1 q?1+a2q?2 B q?1 =b0 +b1q?1

C q?1 =1+c1 q?1+c2q?2 E q?1 =1+(c1 ?a1)q?1

F q?1 =(c2?a2?a1c1+a12)?(a2c1 ?a1a2)q?1 G q?1 =b0+(b0c1 ?b0a1+b1)q?1+b1(c1?a1)q?2

y(k?2)?F(q?1)y(k)?G(q?1)u(k)??(k?2)

?f0y(k)?f1y(k?1)?g0u(k)?g1u(k?1)?g2u(k?2)??(k?2) ②LSE 模型

y(k)?f0y(k?2)?f1y(k?3)?g0u(k?2)?g1u(k?3)?g2u(k?4)??(k)??T???(k)

?T?[y(k?2),y(k?3),u(k?2),u(k?3),u(k?4)] ??[f0,f1,g0,g1,g2]T

③估计器

?(0)?0,p(0)?a2I

?P(k?1)?(k)T?(k)??(k?1)?[y(k)??(k)?(k?1)]???T(k)p(k?1)?(k) ???P(k?1)?(k)?T(k)P(k?1)P(k)?[P(k?1)?]/????T(k)P(k?1)?(k)

④控制器 u(k)?r(k)?F(q?1)y(k)G(q?1)????

????[r(k)?f0y(k)?f1y(k?1)?g1u(k?1)?g2u(k?2)]/g0

⑤程序流程图如下：

初始化数据 Y数据区[y(k),y(k?1),y(k?2),y(kU数据区[u(k),u(k?3),y(k?4)] ?1),u(k?2),u(k?3),u(k?4)] 右移一个时刻 采样，得y(k)，填入y数据 进行常规控制，将u(k)填入u数据 N 采样周期到否？ Y k=k+1 N k=0? Y ?置初值：?(0)?0,p(0)?a2I y,u数据区右移一个时刻 采样，得y(k)，填入y数据 按步骤（2）LSE模型公式，构造数据向量?(k) 按步骤（3）估计器公式，计算?(k) P(k) ?按步骤（4）控制器，计算u(k)并输出，u(k)填入u数据区 N 采样周期到否？ Y

根据程序流程图和所计算的模型编写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和p beta=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

%初始化theta

theta(:,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]; %计算控制率u

u(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)); end figure;

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 on

plot(theta(4,3:i-1),'-c') legend('a1','a2','b0','b1'); title('系数辨识');

各未知系数的仿真结果如下图所示：

系数辨识6050403020100-10-20-30a1a2b0b1020406080100120140160180200

系数辨识2a1a2b0b10-2-4-6-8-10-12020406080100120140160180200

系统中存在有噪声信号的影响，仿真的过程中噪声信号具有随机性，所以每次仿真的所得出的图形不尽相同，但只要最后未知参数达到一稳定值则表明实现了最小方差控制器的功能，使系统达到了稳定。

4.对于系统模型

y(k)?1.6y(k?1)?0.6y(k?2)?1.5u(k?1)?0.53u(k?2)?0.9u(k?3)?ew(k)?0.4ew(k?1) 求系统的一个广义最小方差控制律。

解 ： 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代入式

-1C(q-1)?E(q-1)A(q-1)?q-dD(q-1) P(q) 中 (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-1

代入式 u(k)?-??1(q?1)[?(q?1)y(k)??(q?1)w(k)]

?(q?1)?D(q?1)?1.2?0.6q?1 ?(q?1)?B(q?1)E(q?1)?C(q?1)Q(q?1)?3.5?1.33q?1?0.9q?2 ?(q?1)??C(q?1)R(q?1)??1?0.4q?1P q?1 =1,R q?1 =1,Q q?1 =2

式中w(k)为已知的参考输入量

P,R,和Q’分别为对实际输出、参考输入和控制输入的加权多项式 即可得广义最小方差控制律：

1.2?0.6q?1 y k + ?1+0.4q?1 w(k)

u k =— 3.5?1.33q?1+0.9q?2化简上式后广义最小方差控制律为：

u(k)=-0.34y(k)+0.17y(k-1)+0.38u(k-1)-0.26u(k-2)-0.28w(k)+0.12w(k-1); 将上式代入受控对象模型式

A(q-1)y(k)=B(q-1)u(k-d)+C(q-1)e(k) 可得到采用广义最小方差控制后的闭环系统模型： q?dB(q?1)R(q?1)y(k)?w(k)?1?1?1?1A(q)Q(q)?B(q)P(q)

B(q?1)E(q?1)?C(q?1)Q(q?1)?e(k) A(q?1)Q(q?1)?B(q?1)P(q?1)

本文来源：https://www.bwwdw.com/article/ao9x.html