从飞机的状态矩阵A阵求运动模态的特征根、求传递函数矩阵
1 从状态矩阵求运动模态的特征根、求传递函数矩阵
求特征根
系统的状态空间方程为：
x
˙
=
A
x
+
B
u
(1)
\dot{x}=Ax + Bu \tag{1}
x˙=Ax+Bu(1)
y
=
C
x
+
D
u
(2)
y=Cx+Du \tag{2}
y=Cx+Du(2)
对式（1）两端取Laplace变换，写为运动方程：
(
s
I
−
A
)
x
=
B
u
(3)
(sI-A)x=Bu \tag{3}
(sI−A)x=Bu(3)方程的特征多项式为
ρ
(
s
)
=
d
e
t
(
s
I
−
A
)
(4)
\rho(s)=det(sI-A) \tag{4}
ρ(s)=det(sI−A)(4)即可求A的特征根
求传递函数矩阵
根据 G ( s ) = C ( s I − A ) − 1 B + D (5) G(s) = C(sI-A)^{-1}B+D\tag{5} G(s)=C(sI−A)−1B+D(5)推导过程见参考文献[2]
2 编写MATLAB程序
使用参考文献[3]中的算例：
x
=
(
Δ
V
Δ
α
Δ
θ
Δ
q
)
x=\begin{pmatrix} \Delta V \\ \Delta \alpha \\ \Delta \theta \\ \Delta q\\ \end{pmatrix}
x=
ΔVΔαΔθΔq
u
=
(
δ
e
)
u=\begin{pmatrix} \delta_e \end{pmatrix}
u=(δe)
%% Ex_6_2
% Liang
% 20250401
%% Clear
clear;
%% 定义系统矩阵
A = [-0.014 , 3.5432, -32.2 , 0 ;
-0.0001, -0.806 , 0 , 1 ;
0 , 0 , 0 , 1 ;
0.0007, -8.8077, 0 , -1.3442];
B = [0;
-0.042;
0;
-4.5724];
C = [1, 0, 0, 0;
0, 1, 0, 0;
0, 0, 1, 0;
0, 0, 0, 1];
D = zeros(4, 1);
%% 求特征值
% 第一种方法
egi_A1 = eig(A);
% 第二种方法
I4 = eye(4);
syms s;
Ts = s*I4 - A;
R = solve(det(Ts), s);
egi_A2 = vpa(R);
%% 求传递函数矩阵
Gs = C*(inv(Ts))*B;
%% 增加初始扰动，观察纵向参数响应
% x0 = [0; 5/180*pi; 0; 0]; % 初始扰动
x0 = [0; 0; 0; 0]; % 初始扰动
% u = [0];
% u = @(t) 1/180*pi;
u = @(t) (t>=1) * 1/180*pi;
%% 定义微分方程
odefun = @(t, x) A * x + B * u(t);
%% 求解时间范围
tspan = [0, 300];
%% 使用 ode45 求解
[t, x] = ode45(odefun, tspan, x0);
%% 绘制响应曲线
figure(40101);
xlabel('Time (s)');
ylabel('State Variables');
title('State Response via ODE45');
subplot(4,1,1);
plot(t, x(:,1), 'b'),
legend('\Delta V');
grid on;
subplot(4,1,2);
plot(t, x(:,2)*180/pi, 'r');
legend('\Delta \alpha');
grid on;
subplot(4,1,3)
plot(t, x(:,3)*180/pi, 'c');
legend('\Delta \theta');
grid on;
subplot(4,1,4)
plot(t, x(:,4)*180/pi, 'm');
legend('q');
grid on;
% 设置全屏显示
set(gcf,'outerposition',get(0,'screensize'));
[1] 方振平，陈万春，张曙光. 航空飞行器飞行动力学[M]. 北京：北京航空航天大学出版社，2005.
[2] 吴麒，王诗宓. 自动控制原理（第二版）（上册）[M]. 北京：清华大学出版社，2006.
[3] 刘世前. 现代飞机飞行动力学与控制（第2版）[M]. 上海：上海交通大学出版社，2018.