OpenMDAO飞艇优化中Re数归零导致计算中断的问题求助

问题现象 #

在以两种重量计算方法差值最小化为目标的飞艇体积优化项目中，迭代过程中雷诺数公式 Re=(lb*Vs*pho_cruise_alt)/mu 突然降至0，导致计算摩擦系数 Cf=0.455/(log10(Re)**2.58) 时因对数定义域错误触发程序中断，排查确认变量lb同步归零。

可能原因 #

• 优化器探索行为：梯度类优化器（如SLSQP、IPOPT）在迭代时会试探设计变量的边界值，若lb未设置非零下限，优化器可能会尝试lb=0的极端点。
• 目标函数梯度特性：当lb趋近于0时，目标函数的差值可能出现突变，优化器误将lb=0判定为局部最优候选点。
• 物理约束缺失：未针对lb添加符合工程逻辑的非零约束，导致优化器可自由将其推向0值。

规避方法 #

1. 设置设计变量非零下限
   定义lb为设计变量时，根据飞艇最小物理尺寸设置大于0的下限，示例代码：

   prob.model.add_design_var('lb', lower=1e-3, upper=100.0)
   

2. 优化目标函数平滑性
   给目标函数添加正则项，惩罚lb过小的情况，避免优化器偏向极端值：

   def compute_objective(self, inputs, outputs):
       weight_diff = abs(inputs['weight1'] - inputs['weight2'])
       reg_term = 1e-6 / inputs['lb']  # 正则项，lb越小惩罚越大
       outputs['obj'] = weight_diff + reg_term
   

3. 添加物理约束
   根据工程需求添加不等式约束，强制lb满足最小尺寸要求：

   prob.model.add_constraint('lb', lower=0.1)  # 假设最小长度为0.1单位
   

强制变量非零的可行方案 #

• 硬边界限制：直接在设计变量定义时设置严格非零下限，这是最直接的方法，优化器会自动避开lb<=0的区域。
• 变量替换法：将lb替换为指数变换后的变量，确保其始终大于0：

  class LbTransform(om.ExplicitComponent):
      def setup(self):
          self.add_input('x', val=0.0)
          self.add_output('lb', val=1.0)
      
      def compute(self, inputs, outputs):
          outputs['lb'] = np.exp(inputs['x'])
      
      def compute_partials(self, inputs, partials):
          partials['lb', 'x'] = np.exp(inputs['x'])
  
  prob.model.add_subsystem('lb_transform', LbTransform())
  prob.model.add_design_var('lb_transform.x', lower=-5.0, upper=5.0)  # 对应lb范围约0.0067~148.4
  

• 兜底逻辑处理：在雷诺数计算代码中添加判断，当Re趋近于0时强制设置最小值，避免对数运算报错：

  Re = (lb * Vs * pho_cruise_alt) / mu
  Re = max(Re, 1e-3)  # 兜底，确保Re不小于1e-3
  Cf = 0.455 / (np.log10(Re) ** 2.58)
  

完整项目代码示例 #

import openmdao.api as om
import numpy as np

class WeightCalculator(om.ExplicitComponent):
    def setup(self):
        self.add_input('lb', val=10.0)
        self.add_input('Vs', val=5.0)
        self.add_input('pho_cruise_alt', val=1.225)
        self.add_input('mu', val=1.81e-5)
        self.add_output('weight1', val=0.0)
        self.add_output('weight2', val=0.0)
        self.add_output('Cf', val=0.0)
    
    def compute(self, inputs, outputs):
        lb = inputs['lb']
        Vs = inputs['Vs']
        pho = inputs['pho_cruise_alt']
        mu = inputs['mu']
        
        # 计算雷诺数并添加兜底
        Re = (lb * Vs * pho) / mu
        Re = max(Re, 1e-3)
        # 计算摩擦系数
        Cf = 0.455 / (np.log10(Re) ** 2.58)
        outputs['Cf'] = Cf
        
        # 两种重量计算示例逻辑
        outputs['weight1'] = lb * 100 + Cf * 50
        outputs['weight2'] = lb * 90 + Vs * 20

class Objective(om.ExplicitComponent):
    def setup(self):
        self.add_input('weight1', val=0.0)
        self.add_input('weight2', val=0.0)
        self.add_input('lb', val=10.0)
        self.add_output('obj', val=0.0)
    
    def compute(self, inputs, outputs):
        weight_diff = abs(inputs['weight1'] - inputs['weight2'])
        reg_term = 1e-6 / inputs['lb']  # 添加正则项避免lb过小
        outputs['obj'] = weight_diff + reg_term

# 构建优化问题
prob = om.Problem()
model = prob.model

# 添加子系统
model.add_subsystem('weight_calc', WeightCalculator())
model.add_subsystem('objective', Objective())

# 变量连接
model.connect('weight_calc.weight1', 'objective.weight1')
model.connect('weight_calc.weight2', 'objective.weight2')
model.connect('weight_calc.lb', 'objective.lb')

# 设置设计变量与约束
model.add_design_var('weight_calc.lb', lower=0.1, upper=50.0)
model.add_objective('objective.obj')

# 配置优化器
prob.driver = om.ScipyOptimizeDriver()
prob.driver.options['optimizer'] = 'SLSQP'
prob.driver.options['maxiter'] = 100

# 运行优化
prob.setup()
prob.run_driver()

# 输出结果
print(f"优化后lb值：{prob.get_val('weight_calc.lb')[0]:.4f}")
print(f"目标函数值：{prob.get_val('objective.obj')[0]:.6f}")

内容的提问来源于stack exchange，提问作者Simon Boudoux