华工自动化线性系统第一次大作业

专业：控制理论与控制工程     学号与姓名：

一、飞行器原理及结构和空间坐标系

为了进行控制系统设计的目的，飞机动力学经常称为飞行姿态的一些操作状态进行线性化，它假设飞机的速度（马赫数）和姿态是不变的。控制面（The control surfaces）和发动机推力装置设置或修改，以达到这些状态，我们设计控制系统就是为了维护这些条件，例如，强制将到这些状态的扰动（偏差）变为零。

    通常地，我们将横向运动和纵向运动分开，因为在许多情况下，横向运动和纵向运动的耦合很弱。因此，控制系统的设计可以将各个通道进行考虑。

    飞机机体的轴由X、Y和Z表示，原点选在某些参考点（通常是飞机的重心）。侧倾、俯仰和偏航运动组成围绕x，y和z轴的旋转。θ表示在侧视图向上旋转俯仰角的增加；在飞机的顶视图上，ψ表示偏航角，以逆时针方向增加；在飞机的前视图上，滚动角φ（以逆时针方向增加）。我们让分别表示偏航率、俯仰速率和滚转率。速度矢量V是以U、V和W为投影分别投射到X、Y和Z轴上的。迎角α是速度向量与X轴俯仰方向的夹角，侧滑角β是它与X轴偏航方向的夹角。注意，当角很小的时候，

     飞机的俯仰运动通常由叫做升降舵的控制面控制，翻滚由副翼控制，偏航由方向舵控制。

我们只对其中一些变量感兴趣。作为一个具体的例子，考虑实际飞机的数字数据，一架AFTI-16（F16战斗机的一个修改版本）正在进行着陆（速度V=16英里/小时）。描述飞机纵向运动的状态空间方程给出如下：

控制输入是升降舵角度和向量的状态变量分别是速度的变化，迎角，俯仰速率和俯仰度。

该飞机的纵向模式称为短周期和长周期。在长周期特征值，这也是一种复杂的共轭特征值接近虚轴，造成长周期运动，在水平面缓慢地震荡。

二、状态转移矩阵的重要性与意义

线性时不变系统的响应=零输入响应+零状态响应。即：，其中，为初始条件，为输入，所以在输入矩阵、输入信号、初始状态已知的情况下，只要求出来状态转移阵，就可以求出系统的输出。例如：线性系统输入为一个脉冲信号时，系统的响应为脉冲响应。脉冲响应可以描述系统的输入输出特性。当线性系统初始条件为零时，脉冲响应：。它本质上也是求的系统的状态转移阵。

所以对一个系统而言，的求解是很重要的，只要求出的值就可以计算出系统的输出，也就可以了解到系统的特性。显然，求解状态方程的解的实质是求解状态转移矩阵的问题。

为了得到系统的零输入响应，需要求出状态转移阵，同时可以得到每个状态在任意时刻的输出，由于输出是各个状态的线性组合，所以就可以知道了系统在不同的C阵和D阵下的输出。因此，求解对系统很重要，是求解系统的核心，所以先讨论的求解方法。

三、求状态转移矩阵的几种方法

求的方法有时域的求解方法和频域的求解方法。

方法1：根据或者的定义直接计算：

=I+++…++…=

从公式可以看出，右边是一个无穷项的和，要精确计算出结果是很困难的，所以无论是手工计算还是利用电脑计算，都不可能取无穷项计算，通常是取有限项,得到一个近似的值，以满足不同的精度要求即可.对于不同的精度要求，n的值会不同。在工程上，只要取它的前几项就可以满足要求，本方法易于理解，适合计算机编程。

方法2：利用拉氏反变换法求：

     ==

因为，两边取拉氏变换得：，即

，所以。对上式两边取拉氏反变换，可得：   。

方法3：利用凯莱—哈密顿定理求：

（1）由凯莱—哈密顿定理，方阵A满足其自身的特征方程，即

所以有：

它是的线性组合。

同理，可得，以此类推，都可以用线性表示。

（2）在的定义中，用上面的方法可以消去的n及n以上的幂次项，即

=

   =

当A的特征根互异时，一般计算方法为：

A的特征根均相同时，为时，则

方法4：变换为约旦标准型：

（1）的特征根互异，此时，其中是使变换为对角矩阵的变换阵，。设A的特征根为，i=1,2,3……，n。是A关于的特征向量，满足。可得到基于n个特征向量构造变换阵T，，则状态方程可通过奇异变换而化为Jordan标准型。

，

，。

当A的特征根存在重根时，若A的特征根为，代数重数为4，几何重数为1，则只有一个无关解，则需要3个线性无关向量构成的一组 基。这3个 向量满足下面的性质。。如果A的个特征向量v满足，。v称为generalized eigenvector of grade n。

（2）的特征值有重根，此时，，其中。

,.

四、仿真结果及分析

下面通过飞机的例子说明对系统的状态转移矩阵的求解方法，比较各种方法的不同之处，并对例子进行一定的分析。

例如: 描述AFTI-16战斗机纵向运动的状态空间方程为：

1、利用定义法求及仿真：

的求解及求系统单位脉冲响应的Matlab程序如下：

clear

syms M t X0 B XT

X0=[0;0;0;0];%初始状态。

disp('矩阵A为:')

A=[-0.0507 -3.861 0 -32.2;-0.00117 -0.51 1 0;-0.000129 1.4168 -0.4932 0;0 0 1 0]

%求解2!3!4！5！.

n2=factorial(2);

n3=factorial(3);

n4=factorial(4);

n5=factorial(5);

I=eye(4);

disp('取最高次为5时M的值为：')

M=I+A*t+(1/n2)*A^2*t^2+(1/n3)*A^3*t^3+(1/n4)*A^4*t^4+(1/n5)*A^5*t^5

%求解本系统的状态转移矩阵。

B=[0;-0.0717;-1.5;0];

XT=M*(X0+B) %求解系统的状态响应。

%画图

subplot(2,2,1)

ezplot(XT(1,1),[0,2])%画出x(t)d第一个分量X1(t),并把它显示在左上角。

subplot(2,2,2)

ezplot(XT(2,1),[0,2]) %画出x(t)d第二个分量X2(t),并把它显示在右上角。

subplot(2,2,3)

ezplot(XT(3,1),[0,2]) %画出x(t)d第三个分量X3(t),并把它显示在左下角。

subplot(2,2,4)

ezplot(XT(4,1),[0,2]) %画出x(t)d第四个分量X1(t),并把它显示在右下角。

运行结果为：

图1为系统的单位脉冲响应x（t）的四个状态的仿真波形图。

图1 定义法求单位脉冲响应仿真波形

2、利用拉氏反变换法求及仿真

的求解及求系统单位脉冲响应的Matlab程序如下：

syms M s d1 t XT X0;

A=[-0.0507 -3.861 0 -32.2;-0.00117 -0.51 1 0;-0.000129 1.4168 -0.4932 0;0 0 1 0];

disp('矩阵A的行列式如下:');

d1=det(A);

I=eye(4);

disp('[sI-A]^(-1)为:');

B=(s*I-A);

C=inv(B);

digits(4)

C=vpa(C)

disp('状态转移阵为');

D=ilaplace(C);

digits(4);

M=vpa(D)

X0=[0;0;0;0];

B=[0;-0.0717;-1.5;0];

XT=M*(X0+B) %求解系统的状态响应。

%画图

subplot(2,2,1) %画出x(t)d第一个分量X1(t),并把它显示在左上角。

ezplot(XT(1,1),[0,2])

subplot(2,2,2) %画出x(t)d第二个分量X2(t),并把它显示在右上角。

ezplot(XT(2,1),[0,2])

subplot(2,2,3) %画出x(t)d第三个分量X3(t),并把它显示在左下角。

ezplot(XT(3,1),[0,2])

subplot(2,2,4) %画出x(t)d第一个分量X4(t),并把它显示在右下角。

ezplot(XT(4,1),[0,2])

运行结果为：

图2为系统的单位脉冲响应x（t）的四个状态的仿真波形图。

图2 拉氏反变换法求单位脉冲响应仿真波形

3、利用凯莱—哈密顿定理求及仿真

的求解及求系统单位脉冲响应的Matlab程序如下：

clear

clc

syms s t a M a0 a1 a2 a3 b0 b1 b2 b3 X0 B XT

I=eye(4);

A=[-0.0507 -3.861 0 -32.2;-0.00117 -0.51 1 0;-0.000129 1.4168 -0.4932 0;0 0 1 0];

[V,D]=eig(A);

b1=D(1,1);b2=D(2,2);b3=D(3,3);b4=D(4,4);

c=[1 b1 (b1)^2 (b1)^3;1 b2 (b2)^2 (b2)^3;1 b3 (b3)^2 (b3)^3;1 b4 (b4)^2 (b4)^3];

d=inv(c);

disp('系数矩阵如下:');

f=[exp((b1)*t);exp((b2)*t);exp((b3)*t);exp((b4)*t)];

a=d*f;

a0=a(1,1); %求解公式中的第一个系数

a1=a(2,1); %求解公式中的第二个系数

a2=a(3,1); %求解公式中的第三个系数

a3=a(4,1); %求解公式中的第四个系数

disp('状态转移阵位');

M=a0*I+a1*A+a2*A^2+a3*A^3; %求系统的状态方程

X0=[0;0;0;0];

B=[0;-0.0717;-1.5;0];

XT=M*(X0+B) 

subplot(2,2,1) %画出x(t)d第一个分量X1(t),并把它显示在左上角。

ezplot(XT(1,1),[0,2])

subplot(2,2,2) %画出x(t)d第二个分量X2(t),并把它显示在右上角。

ezplot(XT(2,1),[0,2])

subplot(2,2,3) %画出x(t)d第三个分量X3(t),并把它显示在左下角。

ezplot(XT(3,1),[0,2])

subplot(2,2,4) %画出x(t)d第一个分量X4(t),并把它显示在右下角。

ezplot(XT(4,1),[0,2])

运行结果为：

图3为系统的单位脉冲响应x（t）的四个状态的仿真波形图。

图3 凯莱—哈密顿定理求单位脉冲响应仿真波形

4、利用约旦标准型法求及仿真

的求解及求系统单位脉冲响应的Matlab程序如下：

clear

clc

syms M B N t a0 a1 a2 a3 X0 XT

A=[-0.0507 -3.861 0 -32.2;-0.00117 -0.51 1 0;-0.000129 1.4168 -0.4932 0;0 0 1 0];

[T,D]=eig(A);%求特征值和特征向量.

H=inv(T)*A*T;

a0=D(1,1);

a1=D(2,2);

a2=D(3,3);

a3=D(4,4);

N=[exp(a0*t) 0 0 0;0 exp(a1*t) 0 0;0 0 exp(a2*t) 0;0 0 0 exp(a3*t)];

M=T*N*inv(T);

X0=[0;0;0;0];

B=[0;-0.0717;-1.5;0];

XT=M*(X0+B) 

%其中，M即为所求的结果，B为A矩阵的约旦标准型，T为变换矩阵。

subplot(2,2,1) %画出x(t)d第一个分量X1(t),并把它显示在左上角。

ezplot(XT(1,1),[0,2])

subplot(2,2,2) %画出x(t)d第二个分量X2(t),并把它显示在右上角。

ezplot(XT(2,1),[0,2])

subplot(2,2,3) %画出x(t)d第三个分量X3(t),并把它显示在左下角。

ezplot(XT(3,1),[0,2])

subplot(2,2,4) %画出x(t)d第一个分量X4(t),并把它显示在右下角。

ezplot(XT(4,1),[0,2])

运行结果为：

图4为系统的单位脉冲响应x（t）的四个状态的仿真波形图。

图4 约旦标准型法求单位脉冲响应仿真波形

由上述的四种方法的输出的波形可以看出，四种方法的结果是基本一致的，对于第一个，由于不能取完无穷项，所以和其他的三种有一点的差异，但是大体生仿真波形是一致的，只是存在很小的误差。只要取足够多的项，就能把误差控制在允许的范围之内。可以说明四种方法的结果是相同的，即四种方法都可以求出系统的状态转移矩阵。

由此可知，飞机的纵向运动和速度的变化，迎角，俯仰速率和俯仰度都有关系，飞行器的速度的变化，或者是迎角，俯仰速率和俯仰度的变化都能影响它的状态。

又如：一架时速190英里的飞行在海平面附近的波音747客机的数据。纵向运动的的状态空间描述如下：

利用拉氏反变换法对其进行计算及仿真：

Matlab程序如下：

syms M s d1 t XT X0;

A=[-0.0188 11.5959 0 -32.2;-0.0007 -0.5357 1 0;0.000048 -0.4944 -0.4935 0;0 0 1 0];

disp('矩阵A的行列式如下：');

d1=det(A);

I=eye(4);

disp('[sI-A]^(-1)为：');

B=(s*I-A);

C=inv(B);

digits(4)

C=vpa(C)

disp('状态转移矩阵为');

D=ilaplace(C);

digits(4);

M=vpa(D)

X0=[0;0;0;0];

B=[0;0;-0.5632;0];

XT=M*(X0+B) 

subplot(2,2,1)

ezplot(XT(1,1),[0,2])

subplot(2,2,2)

ezplot(XT(2,1),[0,2])

subplot(2,2,3)

ezplot(XT(3,1),[0,2])

subplot(2,2,4)

ezplot(XT(4,1),[0,2])

运行结果为：

图5为系统的单位脉冲响应x（t）的四个状态的仿真波形图。

图5 波音飞机利用拉氏反变换法求单位脉冲响应仿真波形

同理，用另外三种方法计算状态转移矩阵及仿真结果都是一样的，这里就不再贽述。

对波音747客机和AFTI-16战斗机的脉冲响应的四个状态进行对比，可知，两个飞机的纵向运动和速度的变化，迎角，俯仰速率和俯仰度都有关系，飞行器的速度的变化，或者是迎角，俯仰速率和俯仰度的变化都能影响它的状态。

同时，两个又有不同之处，两者的不稳定的分量的个数不相同，表明对于不同的飞机，受速度的变化，迎角，俯仰速率和俯仰度变化的影响不相同。对于第一个分量，波音747客机和AFTI-16战斗机的变化方向是相同的，但是波音747客机的变化更平稳一些。对于第二个和第四个分量，波音747客机的也是跟平稳一些，这更加适合于乘客的舒适性，符合客机主要强调稳定性的要求；而对于AFTI-16战斗机则是变化更快一点，这对于快速反应的需要比较好，符合战斗机强调快速性的要求。对于第三个分量，AFTI-16战斗机和波音747客机的变化方向不同，后者是趋向于稳定的，而前者则是不稳定的。两者的不同之处反应了飞机设计的时候对于不同性能的要求，一个是追求快速性，一个则是追求舒适性。

对于AFTI-16战斗机的状态空间方程，利用拉氏反变换的方法计算状态转移矩阵，及其单位阶跃响应计算的Matlab程序如下：

clear

clc

syms M s d1 t XT X0 W;

A=[-0.0507 -3.861 0 -32.2;-0.00117 -0.51 1 0;-0.000129 1.4168 -0.4932 0;0 0 1 0];

disp('矩阵A的行列式为：');

d1=det(A);

W=inv(A);

I=eye(4);

disp('[sI-A]^(-1)为：');

B=(s*I-A);

C=inv(B);

digits(4)

C=vpa(C)

disp('状态转移阵为：');

D=ilaplace(C);

digits(4);

M=vpa(D)

X0=[0;0;0;0];

B= [0;-0.0717;-1.5;0];

XT=W*(M-I)*B+M*X0;

subplot(2,2,1)

ezplot(XT(1,1),[0,2])

subplot(2,2,2)

ezplot(XT(2,1),[0,2])

subplot(2,2,3)

ezplot(XT(3,1),[0,2])

subplot(2,2,4)

ezplot(XT(4,1),[0,2])

运行结果为：

图6为系统的单位阶跃响应x（t）的四个状态的仿真波形图。

图6 AFTI-16战斗机利用拉氏反变换法求单位阶跃响应仿真波形

利用定义的的方法计算状态转移矩阵，AFTI-16战斗机单位阶跃响应计算的Matlab程序如下：

clear

clc

syms M t X0 B XT

X0=[0;0;0;0];

disp('矩阵A为:')

A=[-0.0507 -3.861 0 -32.2;-0.00117 -0.51 1 0;-0.000129 1.4168 -0.4932 0;0 0 1 0]

n2=factorial(2);

n3=factorial(3);

n4=factorial(4);

n5=factorial(5);

I=eye(4);

W=inv(A);

disp('取最高次为5：')

M=I+A*t+(1/n2)*A^2*t^2+(1/n3)*A^3*t^3+(1/n4)*A^4*t^4+(1/n5)*A^5*t^5

B=[0;-0.0717;-1.5;0];

XT=W*(M-I)*B+M*X0;

subplot(2,2,1)

ezplot(XT(1,1),[0,2])

subplot(2,2,2)

ezplot(XT(2,1),[0,2])

subplot(2,2,3)

ezplot(XT(3,1),[0,2])

subplot(2,2,4)

ezplot(XT(4,1),[0,2])

运行结果为：

图7为系统的单位阶跃响应x（t）的四个状态的仿真波形图。

图7 AFTI-16战斗机利用定义法求单位阶跃响应仿真波形

通过对比AFTI-16战斗机的单位脉冲响应（图4）和单位阶跃响应（图7）可知，波形有一些不同之处，但是其基本的稳定性是相同的。这是因为其稳定性是系统的基本属性，与所加入的信号的类型无关。它只与系统矩阵A有关系。

对于AFTI-16战斗机，利用Simulink对其单位阶跃响应进行仿真，并与上面的图7进行比较，其单位阶跃响应的Simulink仿真模型和仿真波形如下图所示：

输入是一个单位阶跃信号。C分别取[1 0 0 0]和[0 1 0 0]和[0 0 1 0],[0 0 0 1]，则y的输出就是对应于单位阶跃响应的四个状态x1(t)，x2(t)，x3(t)和x4(t)的四个状态波形。

图8 AFTI-16战斗机单位阶跃响应Simulink仿真模型

图9 AFTI-16战斗机单位阶跃响应Simulink仿真波形

对于AFTI-16战斗机，利用Simulink对其单位脉冲响应进行仿真，并与上面的图4进行比较，其单位脉冲响应的Simulink仿真模型和仿真波形如下图所示：

输入是一个单位脉冲信号，它是由两个阶跃信号叠加而成的。C分别取[1 0 0 0]和[0 1 0 0]和[0 0 1 0],[0 0 0 1]，则y的输出就是对应于单位脉冲响应的四个状态x1(t)，x2(t)，x3(t)和x4(t)的四个状态波形。

图10 AFTI-16战斗机单位脉冲响应Simulink仿真模型

图11 AFTI-16战斗机单位脉冲响应Simulink仿真波形

将Simulink仿真的图形与上面的用Matlab语句的输出的波形比较，两者是相同的。

说明了利用上述四种方法求解状态转移矩阵的方法是正确的。

当输入为方波信号时，AFTI-16战斗机系统的Simulink仿真如下：

图12 AFTI-16战斗机方波信号响应Simulink仿真模型

如上图所示，输入是一个占空比为50%的单位方波信号。

C分别取[1 0 0 0]和[0 1 0 0]和[0 0 1 0],[0 0 0 1]，则y的输出就是对应于方波的响应的x1(t)，x2(t)，x3(t)和x4(t)四个状态波形。如下图所示。

图13 AFTI-16战斗机方波信号响应Simulink仿真波形

因为系统的输出就是4个状态的线性组合，所以输出的特性可以用四个状态所表示。由图可知，系统在阶跃信号和方波信号输入下，均有不稳定的状态，所以当系统的输出含有不稳定的状态时，系统是不稳定的，当不含有那个不稳定的状态时，系统可以有稳定的输出，但系统的本质上是不稳定的，这是由系统的内部结构所决定的。

五、  结论

系统的输出就是4个状态的线性组合，所以输出的特性可以用四个状态所表示。对一个系统而言，的求解是很重要的，只要求出的值就可以计算出系统的输出，也就可以了解到系统的特性。显然，求解状态方程的解的实质是求解状态转移矩阵的问题。

四种求解状态转移矩阵的方法各自有不同的特点，但是结果是相同的。

利用matlab语句生成的输出的各个分量的波形和利用Simulink仿真的波形是相同的，这证明了四种方法都可以得到正确的结果。

对于同一个系统，当输入是不同的时候，波形不同，但是其稳定性只与系统矩阵A有关。由于波音747客机和AFTI-16战斗机用途不同，其对平稳性和快速性的要求也有所差别，波音747客机主要强调平稳性和舒适性，而AFTI-16战斗机则主要是强调快速性，所以它们对同一输入的响应就要分别体现出它们各自的要求。