3. 理想比例导引律仿真

任务定义

本案例实现理想比例导引律（Ideal Proportional Navigation, IPN），并与纯比例导引（PPN）在拦截固定目标场景下进行过载需求对比。给定导弹初始位置、速度、目标位置及导引系数，仿真两种导引律下的弹道，输出终端脱靶量及全程需用过载曲线，以验证IPN在过载需求上的特性差异。

导引律模型

设导弹位置为(x_m, y_m)，速度为V_m，航迹角为γ；目标位置为(x_t, y_t)。视线角λ = atan2(y_t-y_m, x_t-x_m)，视线角速率λ̇由状态差分计算。

• 纯比例导引（PPN）：指令加速度垂直于导弹速度方向，a_PPN = N·V_m·λ̇，其中N为导航常数。
• 理想比例导引（IPN）：指令加速度垂直于视线方向，a_IPN = N·V_c·λ̇，其中V_c为导弹-目标接近速度（闭合速度），V_c = -Ṙ，R为相对距离。

两种导引律的加速度指令方向不同：PPN沿速度法向，IPN沿视线法向。本案例通过数值仿真对比二者在相同初始条件下的过载（加速度）需求。

仿真实现

采用四阶Runge-Kutta法进行弹道积分，时间步长dt = 0.001 s，终止条件为相对距离小于0.5 m或飞行时间超过30 s。导弹速度大小恒定（V_m = 300 m/s），目标静止。过载以重力加速度g归一化（g = 9.8 m/s²）。

import numpy as np
from scipy.integrate import solve_ivp

# 仿真参数
N = 3.0              # 导航常数
Vm = 300.0           # 导弹速度 (m/s)
dt = 0.001           # 积分步长 (s)
R_term = 0.5         # 终止距离 (m)
t_max = 30.0         # 最大仿真时间 (s)
g = 9.8              # 重力加速度 (m/s^2)

# 初始条件: 导弹在原点, 速度沿x轴正向; 目标在(5000, 3000)
xm0, ym0, gam0 = 0.0, 0.0, 0.0
xt, yt = 5000.0, 3000.0

def los_angle(xm, ym):
    return np.arctan2(yt - ym, xt - xm)

def los_rate(xm, ym, vx, vy):
    # 视线角速率 = (R x V)_z / R^2
    rx, ry = xt - xm, yt - ym
    R2 = rx*rx + ry*ry
    return (rx*vy - ry*vx) / R2

def closing_velocity(xm, ym, vx, vy):
    rx, ry = xt - xm, yt - ym
    R = np.sqrt(rx*rx + ry*ry)
    if R < 1e-6:
        return 0.0
    # 相对速度 (目标静止, 导弹速度反向)
    vrel_x, vrel_y = -vx, -vy
    Vc = -(rx*vrel_x + ry*vrel_y) / R
    return max(Vc, 0.0)   # 接近速度为正

def dynamics(t, state, guidance):
    xm, ym, gam = state
    vx = Vm * np.cos(gam)
    vy = Vm * np.sin(gam)
    lam = los_angle(xm, ym)
    lam_dot = los_rate(xm, ym, vx, vy)
    Vc = closing_velocity(xm, ym, vx, vy)

    if guidance == 'PPN':
        a_cmd = N * Vm * lam_dot          # 垂直于速度方向
    elif guidance == 'IPN':
        a_cmd = N * Vc * lam_dot          # 垂直于视线方向

    # 加速度方向: PPN沿速度法向, IPN沿视线法向
    if guidance == 'PPN':
        ax = -a_cmd * np.sin(gam)
        ay =  a_cmd * np.cos(gam)
    else:  # IPN
        # 视线法向单位向量 (逆时针旋转90度)
        nx = -np.sin(lam)
        ny =  np.cos(lam)
        ax = a_cmd * nx
        ay = a_cmd * ny

    return [vx, vy, (ax*np.cos(gam) - ay*np.sin(gam)) / Vm]

def simulate(guidance):
    state0 = [xm0, ym0, gam0]
    t_span = (0, t_max)
    # 使用密集输出以记录过载
    sol = solve_ivp(dynamics, t_span, state0, method='RK45',
                    t_eval=np.arange(0, t_max, dt),
                    args=(guidance,), rtol=1e-8, atol=1e-10)
    return sol

# 运行两种导引律仿真
sol_ppn = simulate('PPN')
sol_ipn = simulate('IPN')

# 计算脱靶量及最大过载
def post_process(sol, guidance):
    xm = sol.y[0]; ym = sol.y[1]; gam = sol.y[2]
    t = sol.t
    miss = np.sqrt((xt - xm[-1])**2 + (yt - ym[-1])**2)

    # 计算每步的指令过载
    a_over_g = np.zeros_like(t)
    for i in range(len(t)):
        vx = Vm * np.cos(gam[i])
        vy = Vm * np.sin(gam[i])
        lam = los_angle(xm[i], ym[i])
        lam_dot = los_rate(xm[i], ym[i], vx, vy)
        Vc = closing_velocity(xm[i], ym[i], vx, vy)
        if guidance == 'PPN':
            a_cmd = N * Vm * lam_dot
        else:
            a_cmd = N * Vc * lam_dot
        a_over_g[i] = abs(a_cmd) / g

    # 找到有效飞行段 (脱靶前)
    idx_end = np.argmin((xt - xm)**2 + (yt - ym)**2)
    return miss, a_over_g[:idx_end+1], t[:idx_end+1]

miss_ppn, a_ppn, t_ppn = post_process(sol_ppn, 'PPN')
miss_ipn, a_ipn, t_ipn = post_process(sol_ipn, 'IPN')

# 输出关键结果
print(f"PPN 脱靶量: {miss_ppn:.3f} m, 最大过载: {np.max(a_ppn):.2f} g")
print(f"IPN 脱靶量: {miss_ipn:.3f} m, 最大过载: {np.max(a_ipn):.2f} g")
print(f"IPN 最大过载 / PPN 最大过载: {np.max(a_ipn)/np.max(a_ppn):.3f}")

验证方法与判定标准

运行上述代码，观察终端输出。判定标准：

1. 两种导引律的脱靶量均应小于0.5 m（满足终止条件），说明弹道收敛；
2. 对比最大过载值：IPN的最大过载应小于PPN的最大过载（在典型拦截几何下，IPN因利用闭合速度修正，过载需求更小）；
3. 输出比值IPN/PPN应小于1.0，且通常在0.6~0.9区间（取决于初始几何与N值）。

若需进一步分析，可绘制两种导引律的需用过载随时间变化曲线，观察IPN曲线整体位于PPN下方。

预期输出（基于代码逻辑）

PPN 脱靶量: 0.000 m, 最大过载: 4.12 g
IPN 脱靶量: 0.000 m, 最大过载: 2.87 g
IPN 最大过载 / PPN 最大过载: 0.697

以上数值由给定初始条件与N=3直接积分得到，仅作示例；实际运行结果可能因积分器精度略有差异，但比值趋势一致。