一、引言
在计算机参与的具有联系受控对象的控制系统中，有必要对联系控制系统设计数字控制器的必要，一般对于联系的控制对象设计数字控制器的方法有：第一种是应用联系系统理论得到的联系控制规律，再将控制规律离散化，用控制器实现，第二种是将联系的控制对象离散化，用离散控制理论设计控制器参数，数字再设计就是根据连续系统及相应的控制规律如何重新设计对应的离散系统与相应的离散控制规律。我们采用的是最优等价准则、双线性变换法、平均增益法进行数字再设计。 二、LQR控制器设计
(1) 二级倒立摆系统的状态空间模型 设线性定常系统为
x’=A*x（t）+B*u（t）,y=C*x（t）
其初始条件为x（t）=x0;
其中：A=[0,1,0,0;40,0,0,0;0,0,0,1;-6,0,0,0];B=[0;-2;0;0.8]; C=[1,0,0,0;0,0,1,0] (2) 系统的能控性判定
n=size(A); Tc=ctrb(A,B); nc=rank(Tc) n=6 6 nc=6
从运行结果可知，系统的阶次为6，能控性矩阵的秩也为6，因此系
统是能控的。 (3) 系统的能观性判定 To=obsv(A,C);no=rank(To) no=6
从运行结果可知，能观性矩阵的秩为6，与系统的阶次相等，因此系统是能观测的。 (4) LQR控制设计
基于一级倒立摆系统具有能控性和能观性，因此可采用LQR进行控制，经大量反复试验和仿真，选取R=0.2，
Q=[1 0 0 0 0 0;0 64 0 0 0 0;0 0 256 0 0 0;0 0 0 0 0 0;0 0 0 0 0 0;0 0 0 0 0 0];
F=lqr(A,B,Q,R)得到： F =
2.2361 106.6465 -155.4620 5.1719 4.9639 -24.5330 三、仿真曲线
采用LQR控制方式，设初始状态为x(0)=[1,-1,0,0]’，在相同采样周期T下应用数字再设计方法对一级倒立摆系统进行仿真，其中F(T)分别取为： 1. F(T)=F1(T)=F
2. F(T)=F2(T)=F[I+(A+BF)T/2] 3. F(T)=F3(T)=F[I-(A+BF)/2]-1
(1) T=0.013s， øc=e(A+BF)T时系统的极点、状态x1、x2、x3的离散仿
真曲线
A=[0,0,0,1,0,0;0,0,0,0,1,0;0,0,0,0,0,1;0,0,0,0,0,0;0,77.0642,-21.1927,0,0,0;0,-38.5321,37.8186,0,0,0]; B=[0;0;0;1;5.7012;-0.0728];
C=[1,0,0,0,0,0;0,1,0,0,0,0;0,0,1,0,0,0]; D=[0;0;0];
Q=[1 0 0 0 0 0;0 64 0 0 0 0;0 0 256 0 0 0;0 0 0 0 0 0;0 0 0 0 0 0;0 0 0 0 0 0]; R=0.2; F=lqr(A,B,Q,R) T=0.013;
[G,H]=c2d(A-B*F,B,T); %%离散一的函数 p0=eig(G), x0=[1 -1 0.5 0 0 0]'; [y,x t]=dinitial(G,B,C,D,x0); t=0:0.1:(t-1)/10;
subplot(3,1,1),x1=[1 0 0 0 0 0]*x'; %%响应曲线 plot(t,x1);grid;title('状态变量x1的响应曲线') subplot(3,1,2),x2=[0 1 0 0 0 0]*x';
plot(t,x2);grid;title('状态变量x2的响应曲线')
subplot(3,1,3),x3=[0 0 1 0 0 0]*x';
plot(t,x3);grid;title('状态变量x3的响应曲线') p0 =
0.8647 + 0.0473i 0.8647 - 0.0473i 0.9224 + 0.0618i 0.9224 - 0.0618i 0.9932 + 0.0066i 0.9932 - 0.0066i
图1 øc=e(A+BF)T
(2) T=0.013s， øc=ø +ΓF1(T)时系统的极点、状态x1、x2、x3的离散仿真曲线
A=[0,0,0,1,0,0;0,0,0,0,1,0;0,0,0,0,0,1;0,0,0,0,0,0;0,77.0642,-21.1927,0,0,
0;0,-38.5321,37.8186,0,0,0]; B=[0;0;0;1;5.7012;-0.0728];
C=[1,0,0,0,0,0;0,1,0,0,0,0;0,0,1,0,0,0]; D=[0;0;0];
Q=[1 0 0 0 0 0;0 64 0 0 0 0;0 0 256 0 0 0;0 0 0 0 0 0;0 0 0 0 0 0;0 0 0 0 0 0]; R=0.2; F=lqr(A,B,Q,R) T=0.013;
[Ad,B]=c2d(A,B,T); %%离散二的函数 Ad=Ad-B*F; p1=eig(Ad) x0=[1 -1 0.5 0 0 0]'; [y,x t]=dinitial(Ad,B,C,D,x0); t=0:0.1:(t-1)/10;
subplot(3,1,1),x1=[1 0 0 0 0 0]*x'; %%显示程序 plot(t,x1);grid;title('状态变量x1的响应曲线') subplot(3,1,2),x2=[0 1 0 0 0 0]*x';
plot(t,x2);grid;title('状态变量x2的响应曲线') subplot(3,1,3),x3=[0 0 1 0 0 0]*x';
plot(t,x3);grid;title('状态变量x3的响应曲线') p1 =
0.8349 + 0.0388i 0.8349 - 0.0388i 0.9247 + 0.0561i 0.9247 - 0.0561i 0.9932 + 0.0066i 0.9932 - 0.0066i
图2 øc=ø +ΓF1(T)
(3) T=0.013s，øc=ø+ΓF2(T)时系统的极点、F(T)值和状态x1、x2、x3的离散仿真曲线
A=[0,0,0,1,0,0;0,0,0,0,1,0;0,0,0,0,0,1;0,0,0,0,0,0;0,77.0642,-21.1927,0,0,0;0,-38.5321,37.8186,0,0,0];
B=[0;0;0;1;5.7012;-0.0728];
C=[1,0,0,0,0,0;0,1,0,0,0,0;0,0,1,0,0,0]; D=[0;0;0];
Q=[1 0 0 0 0 0;0 64 0 0 0 0;0 0 256 0 0 0;0 0 0 0 0 0;0 0 0 0 0 0;0 0 0 0 0 0]; R=0.2; F=lqr(A,B,Q,R) T=0.013;
P2=(A-B*F)*T/2; %%离散3的函数 F2=F*(eye(size(P2))+P2) [Add,B]=c2d(A,B,T); Ad=[Add-B*F2]; p2=eig(Ad) x0=[1 -1 0.5 0 0 0]'; [y,x,t]=dinitial(Ad,B,C,D,x0); t=0:0.1:(t-1)/10;
subplot(3,1,1),x1=[1 0 0 0 0 0]*x'; %%显示程序 plot(t,x1);grid;title('状态变量x1的响应曲线') subplot(3,1,2),x2=[0 1 0 0 0 0]*x';
plot(t,x2);grid;title('状态变量x2的响应曲线')
subplot(3,1,3),x3=[0 0 1 0 0 0]*x';
plot(t,x3);grid;title('状态变量x3的响应曲线') F2 =
1.7236 90.8365 -126.5481 4.0012 4.5195 -19.9211 p2 =
0.8676 + 0.0465i 0.8676 - 0.0465i 0.9224 + 0.0627i 0.9224 - 0.0627i 0.9932 + 0.0066i 0.9932 - 0.0066i
图3 øc=ø+ΓF2(T)
(4) T=0.013s，øc=ø+ΓF3(T)时系统的极点、F(T)值和状态x1、x2、x3的离散仿真曲线
A=[0,0,0,1,0,0;0,0,0,0,1,0;0,0,0,0,0,1;0,0,0,0,0,0;0,77.0642,-21.1927,0,0,0;0,-38.5321,37.8186,0,0,0]; B=[0;0;0;1;5.7012;-0.0728];
C=[1,0,0,0,0,0;0,1,0,0,0,0;0,0,1,0,0,0]; D=[0;0;0];
Q=[1 0 0 0 0 0;0 64 0 0 0 0;0 0 256 0 0 0;0 0 0 0 0 0;0 0 0 0 0 0;0 0 0 0 0 0]; R=0.2; F=lqr(A,B,Q,R) T=0.013;
P3=(A-B*F)*T/2; %%离散4的函数 F3=F*(eye(size(P3))-P3)^-1 [Add,B]=c2d(A,B,T); Ad=[Add-B*F3]; p3=eig(Ad),
[y,x,t]=dinitial(Ad,B,C,D,x0); t=0:0.1:(t-1)/10;
subplot(3,1,1),x1=[1 0 0 0 0 0]*x'; %%显示程序
plot(t,x1);grid;title('状态变量x1的响应曲线') subplot(3,1,2),x2=[0 1 0 0 0 0]*x';
plot(t,x2);grid;title('状态变量x2的响应曲线') subplot(3,1,3),x3=[0 0 1 0 0 0]*x';
plot(t,x3);grid;title('状态变量x3的响应曲线') F3 =
1.7779 92.1683 -129.2365 4.1238 p3 =
0.8655 + 0.0476i 0.8655 - 0.0476i 0.9222 + 0.0622i 0.9222 - 0.0622i 0.9932 + 0.0066i 0.9932 - 0.0066i
4.5459 -20.3464
图4 øc=ø+ΓF3(T)
由上面的1-4图我们可以知道：F(T)分别取F1(T)，F2(T)，F3(T)构成的闭环离散系统时仿真曲线基本一致，相应情况的闭环极点也基本相同，而取F(T)=F3(T)时，从系统的极点看，用øc=ø+ΓF3(T)代替øc=e(A+BF)T 构成闭环系统的精确度相当好。
(5) 当T=0.07s时，重新运行上述程序，运行结果及仿真曲线如下： F =
2.2361 106.6465 -155.4620 5.1719 4.9639 -24.5330 p0 =
0.4410 + 0.1336i 0.4410 - 0.1336i 0.6133 + 0.2309i
0.6133 - 0.2309i 0.9635 + 0.0344i 0.9635 - 0.0344i p1 =
-1.6200 0.4446 0.6537 + 0.1930i 0.6537 - 0.1930i 0.9634 + 0.0345i 0.9634 - 0.0345i F2 =
-0.5233 21.5158 0.3004 p2 =
1.5961 0.4379 0.6521 + 0.1966i 0.6521 - 0.1966i 0.9634 + 0.0345i 0.9634 - 0.0345i F3 =
0.6369 51.2148 -58.6229 1.4954 3.2064 -9.0394
0.2283 -1.1321 2.5709 p3 =
0.5960 + 0.3519i 0.5960 - 0.3519i 0.5336 + 0.1127i 0.5336 - 0.1127i 0.9636 + 0.0342i 0.9636 - 0.0342i
图5 T=0.07s，øc=e(A+BF)T仿真曲线
图6 T=0.07s，øc=ø+ΓF1(T)仿真曲线
图7 T=0.07s，øc=ø+ΓF2(T)仿真曲线
图8 T=0.07s，øc=ø+ΓF3(T)时仿真曲线
由以上图5-8可知，当F(T)取F1(T)或F2(T)构成闭环离散系统时，仿真曲线已经发散，系统变得不稳定了，而用F(T)= F3(T)构成闭环离散系统时，该仿真曲线与øc=e(A+BF)T时系统的仿真曲线基本一致，并且相应的闭环极点仍保持一定的精度。 五、结束语
从两次仿真结果可知，当采样时间很小时，人们通常用连续系统的状态反馈矩阵F构成闭环离散系统，这没有多大问题。但是随着采样时间的增大，仍采取这种方法，则闭环系统的状态响应变坏，甚至出现不稳定现象。这时应用闭环系统离散化的状态矩阵F(T)=F[I-(A+BF)/2]-1作为状态反馈矩阵构成离散化的闭环系统才是一种既简单而又具有较高精度的方法。 六、试验心得
通过此次试验对Matlab的一些应用有了一定的了解，在仿真的过程也碰到了不少问题，开始的有某些函数打错了一些字，在commend框的提示下解决一些问题，还有就是程序中的一些问题，开始参考的程序中漏了一句程序[Add,B]=c2d(A,B,T),使得运行出现错误，经过自己的调试和查找资料，最后终于顺利的完成了试验，于此同时对Matlab有了较好的兴趣，并且发现了其功能的强大，在往后自己会去深刻的学习这个软件，并且多和论坛中的高手交流交流。
因篇幅问题不能全部显示，请点此查看更多更全内容