如何解决solve_bvp中奇异雅可比矩阵导致的收敛失败问题？

问题概述 #

模拟火箭接近火星减速进入最终轨道时，使用scipy.integrate.solve_bvp求解位置与速度变化，但求解器无法收敛，提示**"A singular Jacobian encountered when solving the collocation system."**。仅修改边界条件为Sa[0]（实验用）时可求解，但不符合实际需求，需针对边界条件及物理模型进行修复。

原始代码 #

import numpy as np 
from scipy.integrate import solve_ivp, solve_bvp 
import matplotlib.pyplot as plt

# 定义必要常数
Cd = .2 
A = 11.4  # m^2 
G = 6.673 * 10**-11 # Nm**2/kg**2 
m1 = 5.97219 * 10**24  # kg 
re = 6.371 * 10**6  # m 
Isp = 300  # sec 
Isp2 = 450 # sec 
g0 = 9.81  # m/s^2 
p0 = 101325  # Pa 
M = .0289652  # kg/mol 
R = 8.31446  # J/(mol*K) 
T0 = 288.15  # K 
L = .0065  # K/m 
dry_mass = 21000 # kg 
slv_mass = 2300 # kg 
payload_mass = 2180 # kg 
stage1_fuel_mass = 200000 # kg 
stage2_fuel_mass = 50000 # kg 
deceleration_fuel_mass = stage2_fuel_mass 
total_mass = dry_mass + stage1_fuel_mass + stage2_fuel_mass + payload_mass + slv_mass + deceleration_fuel_mass 
total_mass_2 = stage2_fuel_mass + payload_mass + slv_mass + deceleration_fuel_mass 
m2dot = 2100 # kg/s 
m2dot2 = 400 # kg/s

total_mass_3 = payload_mass + slv_mass + deceleration_fuel_mass

required_velocity = 3000
ta,tb = 0,900
va,vb = 40000, required_velocity

def fun3(t,S): 
    h,v = S 
    dhdt = v 
    dvdt = (m2dot2*(-Isp2*g0+v))/((total_mass_3)-m2dot2*t) 
    return np.vstack([dhdt,dvdt])

def bc(Sa,Sb): 
    bc1 = Sa[1] - va 
    bc2 = Sb[1] - vb 
    return np.array([bc1,bc2])

t_bvp = np.linspace(ta,tb,1000) 
S = np.zeros((2,t_bvp.size)) 
sol3 = solve_bvp(fun3, bc, t_bvp, S, max_nodes = 10000)

print(sol3)

问题分析 #

1. 边界条件欠定：状态变量包含高度h和速度v两个维度，但当前仅约束了速度的初始值与终值，高度无任何边界约束，导致系统存在无穷多解，雅可比矩阵奇异。
2. 物理模型错误：
   • 使用了地球参数而非火星参数，引力计算完全不符合场景；
   • 微分方程未考虑火星引力，仅包含推力项，模型脱离物理实际；
   • 推力公式符号错误，未正确匹配火箭减速的受力方向。
3. 初始猜测不合理：初始猜测全为0，与真实解差距过大，增加了求解器收敛难度。

修复方案与代码 #

关键修复点 #

• 替换为火星物理参数；
• 补充高度边界条件，使系统适定；
• 修正微分方程，加入火星引力项，调整推力公式符号；
• 提供更贴合真实解的初始猜测。

修复后代码 #

import numpy as np 
from scipy.integrate import solve_ivp, solve_bvp 
import matplotlib.pyplot as plt

# 火星物理参数
G = 6.673 * 10**-11       # 万有引力常数 Nm²/kg²
m_mars = 6.39e23          # 火星质量 kg
r_mars = 3389.5e3         # 火星半径 m

# 火箭参数
Isp2 = 450                # 比冲 sec
g0 = 9.81                 # 地球表面重力加速度 m/s²
dry_mass = 21000          # 干质量 kg
slv_mass = 2300           # 上面级质量 kg
payload_mass = 2180       # 有效载荷质量 kg
deceleration_fuel_mass = 50000  # 减速段燃料质量 kg
total_mass_3 = payload_mass + slv_mass + deceleration_fuel_mass
m2dot2 = 400              # 燃料消耗率 kg/s

# 任务参数
required_velocity = 3000  # 目标轨道速度 m/s
ta, tb = 0, 900           # 减速时间区间 s
va = -40000               # 初始接近速度（负号表示朝向火星）m/s
ha = r_mars + 100000      # 初始高度（火星表面上方100km）m

def fun3(t, S): 
    h, v = S 
    # 火星引力加速度（指向火星中心，与远离方向相反）
    g_mars = G * m_mars / (h ** 2)
    # 当前火箭质量
    m = total_mass_3 - m2dot2 * t
    # 排气速度
    exhaust_vel = Isp2 * g0
    # 推力加速度（与接近方向相反，用于减速）
    thrust_acc = exhaust_vel * m2dot2 / m
    # 速度微分：推力减速 + 引力加速（接近火星时引力助力）
    dvdt = thrust_acc - g_mars
    dhdt = v
    return np.vstack([dhdt, dvdt])

def bc(Sa, Sb): 
    # 边界条件：初始高度固定，终态速度达到目标值
    bc1 = Sa[0] - ha
    bc2 = Sb[1] - required_velocity
    return np.array([bc1, bc2])

# 生成初始猜测：高度维持初始值，速度从初始值线性过渡到目标值
t_bvp = np.linspace(ta, tb, 100) 
S_guess = np.zeros((2, t_bvp.size))
S_guess[0, :] = ha
S_guess[1, :] = np.linspace(va, required_velocity, t_bvp.size)

# 求解BVP
sol3 = solve_bvp(fun3, bc, t_bvp, S_guess, max_nodes=10000)

# 输出结果与绘图
print(sol3)
if sol3.success:
    plt.figure(figsize=(10, 6))
    plt.subplot(211)
    plt.plot(sol3.t, sol3.y[0] - r_mars, label='火星表面上方高度')
    plt.ylabel('高度(m)')
    plt.legend()
    plt.subplot(212)
    plt.plot(sol3.t, sol3.y[1], label='速度')
    plt.xlabel('时间(s)')
    plt.ylabel('速度(m/s)')
    plt.legend()
    plt.show()
else:
    print("求解失败：", sol3.message)

修复效果说明 #

• 边界条件补充后，系统变为适定问题，雅可比矩阵奇异问题解决；
• 修正后的物理模型符合火星减速场景，求解器可收敛到合理解；
• 初始猜测贴近真实变化趋势，进一步提升了求解成功率。

内容的提问来源于stack exchange，提问作者Owen Wieland