Scipy.optimize.minimize约束优化失效：结果不满足约束且取变量下限

核心错误点 #

1. 约束函数未正确读取输入参数
   原约束函数中inputs = [a,b]完全错误，这行代码试图用未定义的a、b覆盖传入的优化参数，导致约束计算与优化变量无关，优化器无法根据约束调整变量，最终只能取边界值。必须改为从输入参数中解包变量：

   a, b = inputs
   

2. 未实现阶跃函数D(t)
   代码中多次调用的D(t)是单位阶跃函数（通常定义为x≥0时返回1，否则返回0），但未定义该函数，会直接导致运行错误。需要补充定义：

   def D(x):
       return 1.0 if x >= 0 else 0.0
   

   若使用numpy，也可以用np.heaviside(x, 1.0)替代。

3. 初始猜测值不合理
   初始猜测的b=0.005刚好是变量b的上限，优化器容易卡在边界无法迭代。建议将初始猜测调整到变量范围的中间值，比如：

   first_guess = np.array([0.05, 0.003])
   

4. 约束逻辑验证
   SLSQP的不等式约束要求fun(x) ≥ 0才视为满足约束。原约束逻辑是0.0031 - f_t10minus，对应f_t10minus ≤ 0.0031（挠度不超过最大值），这个逻辑是正确的，但需要确保f_t10minus的计算正确。

修正后的完整代码 #

import math
import numpy as np
from scipy import optimize

# 实现单位阶跃函数
def D(x):
    return 1.0 if x >= 0 else 0.0

# 修正后的约束函数
def deflection_constraint(inputs):
    a, b = inputs  # 正确读取优化参数
    
    mass_propeller_motor = 0.0075
    weight_propeller_motor = mass_propeller_motor * 9.81
    L = 0.035
    density_material = 1200
    density_air = 1.2
    mass_propeller_motor_battery = 0.150
    C_d = 0.7
    v_v = 3
    v_l = 7
    alpha = math.pi/6
    beta = math.pi/3
    
    # 惯性矩计算（矩形截面，b为竖直高度）
    I = (a * (b**3)) / 12
    
    # 静态弯矩计算
    M = weight_propeller_motor * L
    
    # 竖直工况下的弯矩
    V_arm = a * b * L
    V_box = 1.045 * 10**(-5)
    V_frame = 4 * V_arm + V_box
    M_frame = density_material * V_frame
    A_r = 0.0184
    F_w = (M_frame + mass_propeller_motor_battery) * 9.81
    F_d = 0.5 * density_air * A_r * C_d * v_v**2
    F_v = (F_w + F_d) / 4
    M_v = F_v * L
    
    # 侧向工况下的弯矩
    A_rl = A_r * math.cos(beta)
    F_dl = 0.5 * density_air * A_rl * C_d * v_l**2
    F_parallel = F_dl * math.sin(alpha) + F_w * math.cos(alpha)
    F_perpendicular = F_dl * math.cos(alpha) + F_w * math.sin(alpha) 
    F_l = math.sqrt(F_parallel**2 + F_perpendicular**2) * 0.25
    M_l = F_l * L
    
    # 时间参数转换为小时
    t1 = 13.33 / 3600
    t2 = 7186.67 / 3600
    t3 = 2
    t4 = 6
    t5 = 21613.33 / 3600
    t6 = 28786.67 / 3600
    t7 = 8
    t8 = 12
    t9 = 43213.33 / 3600
    t10 = 50386.67 / 3600
    
    # 挠度计算
    f_t10minus = (M_v * D(t10) + 
                  (M_l - M_v) * D(t10 - t1) - 
                  (M_l - M_v) * D(t10 - t2) - 
                  M_v * D(t10 - t3) + 
                  M_v * D(t10 - t4) + 
                  (M_l - M_v) * D(t10 - t5) - 
                  (M_l - M_v) * D(t10 - t6) - 
                  M_v * D(t10 - t7) + 
                  M_v * D(t10 - t8) + 
                  M_v * D(t10 - t9)) * L**2 / (2 * I)
    
    # 约束：挠度不超过0.0031，返回值≥0时满足约束
    total = 0.0031 - f_t10minus
    return total

# 目标函数：最小化横截面积a*b
def g(parameters):
    a, b = parameters
    return a * b

# 优化执行
first_guess = np.array([0.05, 0.003])  # 调整初始猜测到范围中间
my_constraints = ({'type': 'ineq', "fun": deflection_constraint})
result = optimize.minimize(g, 
                          first_guess, 
                          method='SLSQP', 
                          bounds=((0.005, 0.1), (0.001, 0.005)),
                          options={'disp': True},
                          tol=1e-5,
                          constraints=my_constraints)

print("优化结果：", result.x)
print("目标函数值：", result.fun)
print("约束函数值：", deflection_constraint(result.x))

额外注意事项 #

• 检查物理参数合理性：如果优化后仍然无法满足约束，可能是给定的变量范围（比如b的上限0.005）太小，导致即使取最大值也无法满足挠度要求，需要放宽变量边界。
• 数值稳定性：当a或b过小时，惯性矩I会很小，导致挠度计算值过大，约束函数返回负值。可以在约束函数中添加数值检查，避免除以零或溢出。
• 梯度计算：SLSQP默认使用数值梯度，若计算效率低或精度不够，可以手动实现约束函数和目标函数的梯度（雅可比矩阵），提升优化性能。

内容的提问来源于stack exchange，提问作者Giovanni Bugli