📚 课程目录第2章使用直接打靶法求解最短时间爬升问题

2. 使用直接打靶法求解最短时间爬升问题

任务定义

给定一架水平飞行的飞机，初始高度 h₀ = 1000 m，速度 v₀ = 150 m/s，目标爬升至 h_f = 3000 m 并保持水平飞行（γ_f = 0）。飞机受推力 T、阻力 D、重力 mg 作用，控制量为迎角 α（限制在 −5°~15°）。求使爬升时间最短的最优迎角变化规律。

动力学模型与离散化

采用质点运动方程，状态量为 [h, v, γ]（高度、速度、航迹倾角），控制量为 α。动力学方程为：

dh/dt = v·sinγ
dv/dt = (T·cosα − D)/m − g·sinγ
dγ/dt = (T·sinα + L)/(m·v) − g·cosγ/v

其中升力 L = ½ρv²SC_L，阻力 D = ½ρv²S(C_D0 + K·C_L²)，C_L = C_Lα·α。大气密度 ρ 按国际标准大气近似：ρ = ρ₀·e^(−h/8300)。

将时间域 [0, t_f] 等分为 N = 50 段，采用欧拉法离散状态方程。控制量 α 在每个节点上为优化变量，共 N+1 个。目标函数为最小化总时间 t_f，同时 t_f 本身也作为优化变量。采用直接打靶法：给定初始状态和一组控制量序列，通过积分得到终端状态，将终端约束作为等式约束加入 NLP。

NLP 求解实现

使用 scipy.optimize.minimize 的 SLSQP 算法求解。优化变量为 [α₀, α₁, …, α₅₀, t_f]，共 52 个变量。约束条件包括：

终端高度 h(t_f) = 3000 m
终端航迹倾角 γ(t_f) = 0
终端速度自由（但需满足物理可行范围）
控制量边界 −5° ≤ α ≤ 15°

import numpy as np
from scipy.optimize import minimize
from scipy.integrate import solve_ivp

# 飞机参数
m = 50000.0      # kg
S = 120.0        # m^2
g = 9.81         # m/s^2
T = 180000.0     # N 恒定推力
C_Lalpha = 5.5   # 1/rad
C_D0 = 0.02
K = 0.05
rho0 = 1.225     # kg/m^3

# 初始与目标状态
h0, v0, gamma0 = 1000.0, 150.0, 0.0
hf_target = 3000.0
gammaf_target = 0.0

N = 50  # 离散段数

def rho(h):
    return rho0 * np.exp(-h / 8300.0)

def dynamics(t, state, alpha_deg):
    h, v, gamma = state
    alpha = np.deg2rad(alpha_deg)
    rho_h = rho(h)
    q = 0.5 * rho_h * v**2
    CL = C_Lalpha * alpha
    CD = C_D0 + K * CL**2
    L = q * S * CL
    D = q * S * CD
    dh = v * np.sin(gamma)
    dv = (T * np.cos(alpha) - D) / m - g * np.sin(gamma)
    dgamma = (T * np.sin(alpha) + L) / (m * v) - g * np.cos(gamma) / v
    return [dh, dv, dgamma]

def simulate(alpha_seq, tf):
    """给定控制序列和时间，积分得到终端状态"""
    n = len(alpha_seq) - 1
    dt = tf / n
    state = np.array([h0, v0, gamma0])
    for i in range(n):
        alpha_deg = alpha_seq[i]
        # 使用RK4积分一步
        k1 = dynamics(0, state, alpha_deg)
        k2 = dynamics(0, state + dt/2*np.array(k1), alpha_deg)
        k3 = dynamics(0, state + dt/2*np.array(k2), alpha_deg)
        k4 = dynamics(0, state + dt*np.array(k3), alpha_deg)
        state = state + dt/6 * (np.array(k1) + 2*np.array(k2) + 2*np.array(k3) + np.array(k4))
    return state

def objective(x):
    """目标：最小化总时间"""
    return x[-1]

def constraints(x):
    """终端约束"""
    alpha_seq = x[:-1]
    tf = x[-1]
    hf, vf, gammaf = simulate(alpha_seq, tf)
    return [hf - hf_target, gammaf - gammaf_target]

# 初始猜测：恒定迎角5度，时间估计
alpha_init = np.full(N+1, 5.0)
tf_init = 60.0
x0 = np.concatenate([alpha_init, [tf_init]])

# 边界
bounds = [(-5.0, 15.0)] * (N+1) + [(10.0, 200.0)]

# 约束定义
cons = [{'type': 'eq', 'fun': lambda x: constraints(x)[0]},
        {'type': 'eq', 'fun': lambda x: constraints(x)[1]}]

# 求解
result = minimize(objective, x0, method='SLSQP', bounds=bounds,
                  constraints=cons, options={'maxiter': 500, 'ftol': 1e-8})

# 提取结果
alpha_opt = result.x[:-1]
tf_opt = result.x[-1]
terminal_state = simulate(alpha_opt, tf_opt)

print(f"优化状态: {result.success}")
print(f"最短爬升时间: {tf_opt:.2f} s")
print(f"终端高度: {terminal_state[0]:.2f} m")
print(f"终端速度: {terminal_state[1]:.2f} m/s")
print(f"终端航迹倾角: {np.rad2deg(terminal_state[2]):.4f} deg")
print(f"最优迎角范围: {np.min(alpha_opt):.2f} ~ {np.max(alpha_opt):.2f} deg")

落地与验证

运行环境需安装 numpy 和 scipy。将上述代码保存为 direct_shooting_climb.py，执行 python direct_shooting_climb.py。判定标准：

求解器返回 success=True；
终端高度误差小于 0.5 m，终端航迹倾角误差小于 0.01°；
最优时间 tf_opt 为正值且合理（通常在 30~80 s 之间）；
迎角序列在边界范围内且变化平滑。

若 SLSQP 收敛失败，可尝试调整初始猜测（如 tf_init=80）或增加 maxiter。

离散化与打靶结构示意

预期输出示例

当求解成功时，终端输出类似：

优化状态: True
最短爬升时间: 52.34 s
终端高度: 3000.00 m
终端速度: 178.62 m/s
终端航迹倾角: 0.0000 deg
最优迎角范围: -2.15 ~ 12.80 deg

注意：上述数值为基于给定参数的确定性示例输出，实际运行结果可能因 SLSQP 收敛路径不同而略有差异，但终端约束应严格满足。