2. 线性化与状态空间模型提取

任务定义

给定某型飞行器在平衡飞行状态（配平点）附近的纵向运动方程组，通过小扰动线性化方法提取状态空间模型，得到形如 ẋ = Ax + Bu 的线性时不变系统，并验证线性化模型在配平点附近的近似精度。

问题描述与输入数据

飞行器纵向运动的非线性方程组（力与力矩方程）如下：

dV/dt   = (T·cosα - D)/m - g·sinγ
dγ/dt   = (T·sinα + L)/(m·V) - g·cosγ/V
dα/dt   = q - dγ/dt
dq/dt   = M/I_yy
dh/dt   = V·sinγ

其中状态变量为 [V, γ, α, q, h]，控制输入为 [δ_e, δ_T]（升降舵偏角与推力）。在配平点 V₀=250 m/s，γ₀=0，α₀=2°，q₀=0，h₀=5000 m 处，气动力与力矩系数为：

L = ½ρV²·S·C_L,  C_L = C_L0 + C_Lα·α + C_Lδe·δ_e
D = ½ρV²·S·C_D,  C_D = C_D0 + C_Dα·α + C_Dδe·δ_e
M = ½ρV²·S·c̄·C_M,  C_M = C_M0 + C_Mα·α + C_Mq·(q·c̄/(2V)) + C_Mδe·δ_e
T = T₀ + k_T·δ_T

配平条件要求 dγ/dt=0 且 dq/dt=0，由此确定配平时的 δ_e0 与 δ_T0。物理参数：m=12000 kg，S=50 m²，c̄=4 m，I_yy=180000 kg·m²，ρ=0.736 kg/m³（5000m高度），g=9.8 m/s²。气动导数：C_L0=0.2，C_Lα=5.2（1/rad），C_Lδe=0.4（1/rad），C_D0=0.03，C_Dα=0.15（1/rad），C_Dδe=0，C_M0=0.02，C_Mα=-0.8（1/rad），C_Mq=-10（1/rad），C_Mδe=-1.2（1/rad），T₀=50000 N，k_T=2000 N/°。

线性化方法

采用小扰动线性化：设 x = x₀ + Δx，u = u₀ + Δu，将非线性函数 f(x,u) 在配平点处做一阶泰勒展开，忽略高阶项，得到：

Δẋ = A·Δx + B·Δu

其中雅可比矩阵 A = ∂f/∂x|₀，B = ∂f/∂u|₀。由于状态方程中 h 仅通过 V·sinγ 影响，且对纵向短周期与长周期模态影响较弱，本案例提取包含 [ΔV, Δγ, Δα, Δq] 的四阶状态空间模型（忽略高度通道），以聚焦核心纵向动力学。

Python实现

import numpy as np

# 物理参数
m = 12000.0
S = 50.0
cbar = 4.0
Iyy = 180000.0
rho = 0.736
g = 9.8

# 气动导数 (SI单位, 角度转弧度)
C_L0 = 0.2; C_Lalpha = 5.2; C_Lde = 0.4
C_D0 = 0.03; C_Dalpha = 0.15; C_Dde = 0.0
C_M0 = 0.02; C_Malpha = -0.8; C_Mq = -10.0; C_Mde = -1.2
T0 = 50000.0; kT = 2000.0  # N/deg

# 配平点状态
V0 = 250.0
gamma0 = 0.0
alpha0 = np.deg2rad(2.0)
q0 = 0.0
h0 = 5000.0

# 配平计算: 由 dgamma/dt=0 和 dq/dt=0 求 delta_e0, delta_T0
# 推力 T = T0 + kT*delta_T (delta_T 单位度)
# 升力 L = 0.5*rho*V^2*S*(C_L0 + C_Lalpha*alpha + C_Lde*delta_e)
# 阻力 D = 0.5*rho*V^2*S*(C_D0 + C_Dalpha*alpha)
# 力矩 M = 0.5*rho*V^2*S*cbar*(C_M0 + C_Malpha*alpha + C_Mq*(q*cbar/(2V)) + C_Mde*delta_e)

qbar = 0.5 * rho * V0**2  # 动压

# 由 dq/dt=0: M=0
# C_M0 + C_Malpha*alpha0 + C_Mq*(q0*cbar/(2V0)) + C_Mde*delta_e0 = 0
delta_e0 = -(C_M0 + C_Malpha*alpha0 + C_Mq*(q0*cbar/(2*V0))) / C_Mde

# 由 dgamma/dt=0: (T*sin(alpha0) + L)/(m*V0) - g*cos(gamma0)/V0 = 0
# 且 gamma0=0, cos(gamma0)=1
L0 = qbar * S * (C_L0 + C_Lalpha*alpha0 + C_Lde*delta_e0)
D0 = qbar * S * (C_D0 + C_Dalpha*alpha0)
# T*sin(alpha0) + L0 = m*g
T_trim = (m*g - L0) / np.sin(alpha0) if abs(np.sin(alpha0)) > 1e-6 else 0
delta_T0 = (T_trim - T0) / kT

print(f"配平舵偏 delta_e0 = {np.rad2deg(delta_e0):.4f} deg")
print(f"配平推力 delta_T0 = {delta_T0:.4f} deg")

# 定义非线性函数 f(x,u) = [dV/dt, dgamma/dt, dalpha/dt, dq/dt]
def f_nonlin(x, u):
    V, gamma, alpha, q = x
    de, dT = u
    T = T0 + kT * dT
    rho_local = rho  # 假设高度不变
    qbar_loc = 0.5 * rho_local * V**2
    L = qbar_loc * S * (C_L0 + C_Lalpha*alpha + C_Lde*de)
    D = qbar_loc * S * (C_D0 + C_Dalpha*alpha + C_Dde*de)
    M = qbar_loc * S * cbar * (C_M0 + C_Malpha*alpha + C_Mq*(q*cbar/(2*V)) + C_Mde*de)
    dV = (T*np.cos(alpha) - D)/m - g*np.sin(gamma)
    dgamma = (T*np.sin(alpha) + L)/(m*V) - g*np.cos(gamma)/V
    dalpha = q - dgamma
    dq = M / Iyy
    return np.array([dV, dgamma, dalpha, dq])

# 配平点
x0 = np.array([V0, gamma0, alpha0, q0])
u0 = np.array([delta_e0, delta_T0])

# 数值雅可比 (中心差分)
def jacobian(f, x, u, eps=1e-6):
    n = len(x); m = len(u)
    A = np.zeros((n, n))
    B = np.zeros((n, m))
    for i in range(n):
        xp = x.copy(); xm = x.copy()
        xp[i] += eps; xm[i] -= eps
        A[:, i] = (f(xp, u) - f(xm, u)) / (2*eps)
    for j in range(m):
        up = u.copy(); um = u.copy()
        up[j] += eps; um[j] -= eps
        B[:, j] = (f(x, up) - f(x, um)) / (2*eps)
    return A, B

A, B = jacobian(f_nonlin, x0, u0)

print("\n状态矩阵 A (状态: [ΔV, Δγ, Δα, Δq]):")
print(np.array2string(A, precision=4, suppress_small=True))
print("\n输入矩阵 B (输入: [Δδ_e, Δδ_T]):")
print(np.array2string(B, precision=4, suppress_small=True))

# 验证: 在配平点附近给一个小扰动, 比较非线性与线性模型
delta_x = np.array([2.0, 0.01, 0.02, 0.01])  # ΔV=2m/s, Δγ=0.01rad, Δα=0.02rad, Δq=0.01rad/s
delta_u = np.array([0.5, 100.0])  # Δδ_e=0.5deg, Δδ_T=100N(对应0.05deg)

# 非线性响应
x_nl = x0 + delta_x
u_nl = u0 + delta_u
f_nl = f_nonlin(x_nl, u_nl)
f0 = f_nonlin(x0, u0)
nonlinear_delta = f_nl - f0

# 线性响应
linear_delta = A @ delta_x + B @ delta_u

print("\n非线性增量: ", np.array2string(nonlinear_delta, precision=6))
print("线性增量:   ", np.array2string(linear_delta, precision=6))
print("误差:       ", np.array2string(nonlinear_delta - linear_delta, precision=6))

# 判定: 相对误差
rel_err = np.linalg.norm(nonlinear_delta - linear_delta) / (np.linalg.norm(nonlinear_delta) + 1e-12)
print(f"\n相对误差: {rel_err:.4%}")
print("若相对误差 < 5%, 线性化模型在扰动范围内有效")

运行与验证

直接运行上述Python脚本（需安装NumPy），输出配平点、A矩阵、B矩阵以及线性化误差验证结果。判定标准：

配平计算得到的 δ_e0 与 δ_T0 应使非线性方程在配平点处的导数为零（数值上接近0）；
在给定小扰动下，线性模型输出与非线性模型输出的相对误差小于5%，表明线性化模型在配平点附近有效；
A矩阵特征值实部为负（稳定）或按需分析模态特性。

预期输出示例

配平舵偏 delta_e0 = 1.2345 deg
配平推力 delta_T0 = 3.4567 deg

状态矩阵 A (状态: [ΔV, Δγ, Δα, Δq]):
[[-0.0123  9.8100 -0.5000  0.0000]
 [-0.0001 -0.0010  0.8000  0.0000]
 [ 0.0001  0.0010 -0.8000  1.0000]
 [ 0.0000  0.0000 -2.5000 -0.6000]]

输入矩阵 B (输入: [Δδ_e, Δδ_T]):
[[ 0.0000  0.8000]
 [ 0.0000  0.0000]
 [ 0.0000 -0.0001]
 [-1.2000  0.0000]]

非线性增量: [ 0.123456 -0.001234  0.056789 -0.045678]
线性增量:   [ 0.123400 -0.001230  0.056700 -0.045600]
误差:       [ 0.000056 -0.000004  0.000089 -0.000078]

相对误差: 0.1234%
若相对误差 < 5%, 线性化模型在扰动范围内有效

以上数值为示例格式，实际运行结果取决于具体参数与扰动选取。判定标准为相对误差小于5%即认为线性化模型在给定扰动范围内有效。

📚 课程目录