兰州理工大学线性系统理论期末MATLAB大作业

1、线性系统理论的Matlab实践1、在造纸流程中，投料箱应该把纸浆流变成2cm的射流，并均匀喷洒在网状传送带上。为此，要精确控制喷射速度和传送速度之间的比例关系。投料箱内的压力是需要控制的主要变量，它决定了纸浆的喷射速度。投料箱内的总压力是纸浆液压和另外灌注的气压之和。由压力控制的投料箱是个耦合系统，因此，我们很难用手工方法保证纸张的质量。在特定的工作点上，将投料箱线性化，可以得到下面的状态空间模型： x = -0.8+0.02-0.020 x+0.0510.0010 u y =x1 , x2 其中，系统的状态变量x1=液面高度，x2=压力，系统的控制变量u1=纸浆流量u2=气压阀门的开启量。在

2、上述条件下，试设计合适的状态变量反馈控制器，使系统具有实特征根，且有一个根大于5解：下面是对此设计的MATLAB程序实现：>> A=-0.8 0.02;-0.02 0;>> B=0.05 1;0.001,0;>> r=rank(ctrb(A,B)r =2>> C=1 1;>> P=1 6;>> K=place(A,B,P)K = 1.0e+003 * -0.0200 -6.0000 -0.0008 0.30002、描述恒速制导导弹的运动方程为： x = 01000-0.1-0.50000.500000 010000.510

3、00 x + 01000 u y= 0 0 0 1 0 x (a) 运用ctrb函数计算系统的能控型矩阵，并验证系统是不可控的；(b) 计算从u到Y的传递函数，并消去传递函数中的分子和分母公因式，由此可以得到能控的状态空间模型。在消去了公因子之后，请用tf2ss函数确定新的状态变量模型；(c) 证明(b)中得到的状态变量模型是能控的；(d) 说明恒速制导导弹是否稳定？(e) 讨论状态变量模型的能控性和复杂性的关系（假设用状态变量的数目来度量复杂性）解程序如下:clearA=input('请输入系统矩阵:');B=input('请输入输入矩阵:');C=input

4、('请输入输出矩阵:');Qc1=ctrb(A,B)N1=size(A);n1=N1(1) %判断状态方程维数rc1=rank(Qc1)if rc1=n1 disp('系统可控')elseif rc1<n1 disp('系统不可控')endsyms s I=eye(n1);Q=inv(s*I-A);sys=collect(C*Q*B) %求解原状态方程的频域传递函数并化简num=500 250 50;den=1 0 0;A1 B1 C1 D1=tf2ss(num,den)Qc2=ctrb(A1,B1)N2=size(A1);n2=N2(1)

5、 %判断状态方程维数rc2=rank(Qc2)if rc2=n2 disp('系统可控')elseif rc2<n2 disp('系统不可控')endd1=eig(A)'d2=eig(A1)'flag1=0;flag2=0;for i=1:n1 if real(d1(i)>0 flag1=1; endendif flag1=1 disp('原系统不稳定')else disp('原系统稳定')endfor j=1:n2 if real(d2(j)>0 flag2=1; endendif flag2=

6、1 disp('新系统不稳定')else disp('新系统稳定')end运行结果：请输入系统矩阵:0 1 0 0 0;-0.1 -0.5 0 0 0;0.5 0 0 0 0;0 0 10 0 0;0.5 1 0 0 0请输入输入矩阵:0;1;0;0;0请输入输出矩阵:0 0 0 1 0Qc1 = 0 1.0000 -0.5000 0.1500 -0.0250 1.0000 -0.5000 0.1500 -0.0250 -0.0025 0 0 0.5000 -0.2500 0.0750 0 0 0 5.0000 -2.5000 0 1.0000 0 -0.100

7、0 0.0500n1 = 5rc1 = 4系统不可控 sys = 50/s2/(10*s2+5*s+1) A1 = 0 0 1 0B1 = 1 0C1 = 250 50D1 = 500Qc2 = 1 0 0 1n2 = 2rc2 = 2系统可控d1 = 0 0 0 -0.2500 - 0.1936i -0.2500 + 0.1936id2 = 0 0原系统稳定新系统稳定分析：由上述分析结果可知原系统和新系统均稳定，而实际上由系统的极点可知，原系统是稳定的，新系统实际上处于临界稳定状态也可认为是不稳定的；若以状态变量的数目来度量复杂性，可知系统的完全可控性与复杂性存在类似反比的关系，及复杂性越高

8、系统完全可控的难度越大，复杂性越低系统完全可控的难度越低。3、垂直起降的飞机的线性化模型为：x=Ax+B1u1+B2u2其中A = -0.03660.02710.0188-0.04550.0482-0.01000.0214-4.02080.10020.3681-0.70701.42000000B1 = 0.44223.5466-5.52000 ， B2 = 0.1761-7.59224.49000系统的状态变量为水平速度（节）、垂直速度（节）、倾斜率（度/秒）和倾斜角（度）；系统的控制输入为和，其中用于控制垂直运动，用于控制水平运动。(a) 计算系统矩阵的特征值，并由此判断系统是否稳定；(b)

9、 利用poly函数确定的特征多项式，计算特征根，并与(a)中得到的特征根相比较；(c) 当只有发挥作用时，系统能控吗？当只有发挥作用时，结果又如何？请比较解释你的结论。解：矩阵A的特征值由下列方式实现：>> A=-0.0366 0.0271 0.0188 -0.4555;0.0482 -1.0100 0.0024 -4.0208;0.1002 0.3681 -0.7070 1.4200;0 0 1 0;>> Lambda=eig(A)Lambda = 0.2758 + 0.2576i 0.2758 - 0.2576i -0.2325 -2.0727 由上可知，系统有两个

10、负实根-0.2325和-2.0727，两个实部为正的共轭复根0.2758 + 0.2576i和0.2758 - 0.2576i，而要使系统渐进稳定，所有的特征根必须都具有负实部，所以此系统为不稳定系统。由poly函数确定系统特征多项式的实现如下：>> A=-0.0366 0.0271 0.0188 -0.4555;0.0482 -1.0100 0.0024 -4.0208;0.1002 0.3681 -0.7070 1.4200;0 0 1 0;>> b=poly(A)b = 1.0000 1.7536 -0.6472 0.0625 0.0686上述b的值为系统特征多项

11、式的系数，则系统特征多项式为a(s)= s4+1.7536s3-0.6472s2+0.0625s+0.0686，计算此特征多项式的根有如下实现过程：>> b=1 1.7536 -0.6472 0.0625 0.0686;>> roots(b)ans = -2.0727 0.2758 + 0.2576i 0.2758 - 0.2576i -0.2324 由特征多项式所得特征根为两个负实根 -2.0727和 -0.2324，两个具有正实部的共轭复根0.2758 + 0.2576i和0.2758 - 0.2576i，将其与(a)中所得特征根比较如下：(a): -2.0727

12、-0.2325 0.2758 + 0.2576i 0.2758 - 0.2576i (b): -2.0727 -0.2324 0.2758 + 0.2576i 0.2758 - 0.2576i 可以看出，(a)和(b)所得的系统特征值只有一个负实根不相同之外其他的特征根都相同，而不相同的两个负实根-0.2325和-0.2324只相差0.0001，相差甚微，仅仅是这么小的差距对分析系统性能并没有很大的影响，完全可以忽略。(c)当只有u1作用于系统时，即输入矩阵只有B1矩阵，由系统矩阵A和输入矩阵B1确定的能控性判别矩阵可知系统的能控性，对此有如下实现过程：>> A=-0.0366 0

13、.0271 0.0188 -0.4555;0.0482 -1.0100 0.0024 -4.0208;0.1002 0.3681 -0.7070 1.4200;0 0 1 0;>> B1=0.4422;3.5466;-5.5200;0;>> Co1=ctrb(A,B1)Co1 = 0.4422 -0.0238 2.5171 -2.0270 3.5466 -3.5740 25.8160 -47.1028 -5.5200 5.2525 -12.8699 26.3126 0 -5.5200 5.2525 -12.8699>> r=rank(Co1)r =4系统状态

14、维数n=4，又有r=n=4，所以只有u1单独作用时系统完全能控。当只有u2作用于系统时，相应的输入矩阵只有B2，可以由系统矩阵A和输入矩阵B2确定的能控性判别矩阵判断系统的能控性，具体的实现如下过程：>> B2=0.1761;-7.5922;4.4900;0;>> Co2=ctrb(A,B2)Co2 = 0.1761 -0.1278 -1.9441 2.3338 -7.5922 7.6874 -25.8381 49.9646 4.4900 -5.9515 13.4004 -27.6310 0 4.4900 -5.9515 13.4004>> r=rank(C

15、o2)r =4同样有r=n=4，所以只有u2单独作用于系统时系统也是完全可控的。4、为了探究月球背面（远离地球的一面）的奥秘，人们付出了不懈的努力。例如，在地球-太阳-月球系统中，人们希望通信卫星能定点在不受月球遮挡的轨道上，并为此开展了广泛的论证研究工作。图中给出了预期卫星轨道的示意图，从地球上看上去，卫星轨道的光影恰似环绕月球的外层光晕，因此这种轨道又称为光晕轨道。轨道控制的目的是，使通信卫星在地球可见的光晕轨道上运行，从而保证通信链路的畅通，所需的通信链路包括从地球到卫星和从卫星到月球背面共两段线路。卫星绕定点位置运动时，经过标准化和线性化的漂移运动方程为：其中，状态变量是卫星在三个方向

16、上的位置和速度漂移，输入分别是轨控发动机在、和方向上产生的加速度。(a) 卫星的定点位置是否稳定？(b) 如果只有发挥作用，卫星是否能控？(c) 如果只有发挥作用，卫星是否能控？(d) 如果只有发挥作用，卫星是否能控？(e) 如果能够测得方向的位置漂移，请确定由到该位置漂移量的传递函数。（提示：可以令观测输出为）(f) 用tf2ss函数，计算(e)中得到的传递函数的状态变量模型，并验证该轨迹子系统是能控系统；(g) 采用状态反馈，设计合适的反馈控制器，使(f)中得到的系统的闭环极点为和。解：(a)由系统矩阵A的特征值可以判断系统的稳定性，如下可知系统含有正实数特征根，故系统是不稳定的>&

17、gt; A=0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;7.3809 0 0 0 2 0;0 -2.1904 0 -2 0 0;0 0 -3.1904 0 0 0>> Lambda=eig(A)Lambda = -2.1587 2.1587 0.0000 + 1.8626i 0.0000 - 1.8626i 0 + 1.7862i 0 - 1.7862i(b) 系统只有u1发挥作用时，利用系统矩阵A和输入矩阵B1得到的能控性判别矩阵Co1判断系统的能控性，有如下实现过程：>> A=0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0

18、 0 1;7.3809 0 0 0 2 0;0 -2.1904 0 -2 0 0;0 0 -3.1904 0 0 0;>> B1=0;0;0;1;0;0;>> Co1=ctrb(A,B1)Co1 = 0 1.0000 0 3.3809 0 20.1921 0 0 -2.0000 0 -2.3810 0 0 0 0 0 0 0 1.0000 0 3.3809 0 20.1921 0 0 -2.0000 0 -2.3810 0 -35.1688 0 0 0 0 0 0>> r1=rank(Co1)r1 =4系统状态变量的维数n=5，由于r1=4<n=5，所

19、以只有u1发挥作用时，卫星定点系统不完全能控。(c)当系统只有u2发挥作用时，由系统矩阵A和输入矩阵B2构成能控性判别矩阵Co2，判断系统能控性的实现过程如下：>> A=0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;7.3809 0 0 0 2 0;0 -2.1904 0 -2 0 0;0 0 -3.1904 0 0 0;>> B2=0;0;0;0;1;0;>> Co2=ctrb(A,B2)Co2 = 0 0 2.0000 0 2.3810 0 0 1.0000 0 -6.1904 0 8.7975 0 0 0 0 0 0 0 2.

20、0000 0 2.3810 0 35.1688 1.0000 0 -6.1904 0 8.7975 0 0 0 0 0 0 0>> r2=rank(Co2) r2=4同样有r2=4<n=6，故只有u2作用时，卫星定点系统也是不完全能控的。(d)当系统只有u3发挥作用时，由系统矩阵A和输入矩阵B3构成能控性判别矩阵Co3，判断系统能控性的实现过程如下：>> A=0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;7.3809 0 0 0 2 0;0 -2.1904 0 -2 0 0;0 0 -3.1904 0 0 0;>> B3=0;

21、0;0;0;0;1;>> Co3=ctrb(A,B3)Co3 = 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0000 0 -3.1904 0 10.1787 0 0 0 0 0 0 0 0 0 0 0 0 1.0000 0 -3.1904 0 10.1787 0>> r3=rank(Co3)r3 =2可知r3=2<n=6,所以当只有u3作用于系统时，卫星定点系统也是不可控的。(e) 以U2为输入到该位置漂移量求得系统的传递函数MATLAB实现如下：>> syms s>> A=0 0 0 1 0 0;0 0 0 0 1 0;0

22、0 0 0 0 1;7.3809 0 0 0 2 0;0 -2.1904 0 -2 0 0;0 0 -3.1904 0 0 0；>> B=0;0;0;0;0;1;>> C=0 1 0 0 0 0;>> I=eye(6)；>> F=inv(s*I-A)；>> G=simple(simple(C*F*B)；最后得到系统传递函数为：G =-(6250000*s2 - 46130625)/(7440625*s2 6250000*s4 + 101044521)(f)由上述所得传递函数获取状态变量模型有如下实现：>> num=-625

23、0000 0 46130625;>> den=-6250000 0 7440625 0 101044521;>> G=tf(num,den)>> sys=ss(G)a = x1 x2 x3 x4x1 0 0.5952 0 2.021x2 2 0 0 0x3 0 2 0 0x4 0 0 2 0b = u1 x1 1 x2 0 x3 0 x4 0c = x1 x2 x3 x4 y1 0 0.5 0 -0.9226d = u1 y1 0 上述系统的能控性判别矩阵及其秩如下：QC = 1.0000 0 1.1904 0 0 2.0000 0 2.3808 0 0 4

24、.0000 0 0 0 0 8.0000>> r=rank(QC)r = 4r=4=n故系统是完全能控的(g)由于是单输入系统，故在极点配置的过程中对于重极点的情形会在实现时程序出错，故将给定的两个重极点-10中的一个取为与其很接近的-10.0001，这样的取法对系统的性能影响甚微，故具备可行性，对此的极点配置的状态反馈矩阵有如下实现：>> A=0 0.5952 0 2.021;2 0 0 0;0 2 0 0;0 0 2 0;>> B=1;0;0;0;>> C=0 0.5 0 -0.9266;>> P=-1+j -1-j -10 -1

25、0.0001;>> K=place(A,B,P)K = 22.0001 71.5958 60.0005 27.0212则有u2=Kx所得控制输入为：u2 = 22.0001 71.5958 60.0005 27.0212x5、在8.2风力机的一阶模型中，采用浆距角控制风力机的转速，风速的变化视为扰动，设计风力机转速的闭环PI控制，使转速恒定。解：要设计风力机转速的闭环PI控制，首先绘制出带有PI控制器的风力机控制系统，如图所示：带有PI控制器的风力机控制器系统框图图中PI控制器的模型为K1+K2/s，从而可以求出该系统闭环系统的特征方程为S2+0.3397K1*S+0.3397K2

26、=0要使得该风力机稳定运行，则需让特征方程的根具有负实部。设计程序如下：syms s k1 k2;>> s=solve(s2+0.339*k1*s+0.339*k2) s = - (339*k1)/2000 - (114921*k12)/1000000 - (339*k2)/250)(1/2)/2 (114921*k12)/1000000 - (339*k2)/250)(1/2)/2 - (339*k1)/2000使其为负数，从而可得出K1<=180，K2<=86,得到的PI控制器就能满足要求。6、在8.2风力机的三阶模型中，采用浆距角控制风力机的转速，风速的变化视为扰

27、动，电磁转矩视为常数，采用状态反馈和极点配置算法，设计风力机转速的闭环控制系统。解：即确定状态反馈矩阵K，从而使得系统矩阵的特征根为期望特征根，不妨设期望的系统矩阵特征根为-2，-3，-4；设计程序如下：A=-0.8468 -8.1598e-008 0.5071;867637000 0 -867637000;1.1636e+004 0.0019 -1.1636e+004;>> B=-0.0262;0;0;>> Qc=ctrb(A,B)Qc = 1.0e+011 * -0.0000 0.0000 -0.0000 0 -0.0002 2.6453 0 -0.0000 0.0

28、000>> Rc=rank(Qc)Rc = 3>> P=-2 -3 -4;>> K=place(A,B,P)K = 1.0e+005 *4.4381 0.0000 -4.4380按照所求K所构成的状态反馈矩阵从而形成的状态反馈能够满足要求。7、问题同上，设计风力机转速的LQR控制器。解：以浆距角为控制输入的三阶模型如下：=x+v要对系统进行状态反馈极点配置需先判断系统能控性，能控性判别矩阵为：Qc = 1.0e+011 * -0.0000 0.0000 -0.0000 0 -0.0002 2.6453 0 -0.0000 0.0000可以看出rankQc=3

29、，系统完全能控，则可以进行任意极点的配置。如下计算系统的特征值：>> A=-0.8468 -8.1598e-008 0.5071;867637000 0 -867637000;1.1636e+004 0.0019 -1.1636e+004;>> Lambda=eig(A)Lambda = 1.0e+004 * -1.1493 -0.0000 -0.0143为了满足控制要求，为系统取定几个位于左半S平面的几个闭环极点，首先取离虚轴较近的一对共轭主导极点，其余的极点取到离虚轴的距离为主导极点距离虚轴3-6倍，则系统性能主要取决于主导极点，而受其他非主导极点的影响不大，基此取

30、定如下三个期望的系统闭环极点：S1,2=-2±4j，S3=-10对期望极点进行配置算法如下：A=-0.8468 -8.1598e-008 0.5071;867637000 0 -867637000;1.1636e+004 0.0019 -1.1636e+004;>> B=-0.0262;0;0;>> P=-2+4j -2-4j -10;>> K=place(A,B,P)K = 1.0e+005 * 4.4362 0.0000 -4.4361状态反馈之后闭环系统的方框图如下所示：X=AX+B1+B2V U -K8、风力机三阶模型中表示系统轴上的扭转弹

31、力第二个状态是不容易测量的状态变量，设计关于该状态的观测器，并构成状态反馈系统。解：本题目要设计风力机三阶模型的状态反馈及其状态观测器，并且为降维观测器，由于输出矩阵C的秩为1，故而降维观测器的维数为3-1=2；不妨设状态反馈的期望极点为-3，-4，-5，状态观测器的期望极点为-7，-8则相应的设计程序如下：设计程序：>> A=-0.8468 -8.1598e-008 0.5071;867637000 0 -867637000;1.1636e+004 0.0019 -1.1636e+004;>> B=-0.0262;0;0;>> C=0 0 1;>> Qc=ctrb(A,B)Qc = 1.0e+011 * -0.0000 0.0000 -0.0000 0 -0.0002 2.6453 0 -0.0000 0.0000>> Rc=rank(Qc)Rc = 3>> Qo=obs