专业：控制理论与控制工程 学号与姓名：
一、飞行器原理及结构和空间坐标系
为了进行控制系统设计的目的，飞机动力学经常称为飞行姿态的一些操作状态进行线性化，它假设飞机的速度（马赫数）和姿态是不变的。控制面（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战斗机则主要是强调快速性，所以它们对同一输入的响应就要分别体现出它们各自的要求。下载本文
