使用 openMDAO 处理 SLSQP 优化中的约束

Question

使用 openMDAO，我正在使用 FD 导数并尝试使用 SLSQP 方法解决非线性约束优化问题。每当优化器到达违反其中一个约束的点时，它就会崩溃并显示以下消息：

优化失败。线搜索的正向导数

例如，如果我故意将初始点设置为不可行的设计点，优化器执行 1 次迭代并退出并出现上述错误（当我从可行点开始时也会发生同样的情况，但随后优化器到达不可行点经过几次迭代）。

根据In OpenMDAO, is there a way to ensure that the constraints are respected before proceeding with a computation? 中问题的答案，我假设在我的情况下引发AnalysisError 异常将不起作用，对吗？有没有其他方法可以防止优化器进入不可行的区域或至少回溯在线搜索并尝试不同的方向/距离？还是应该只在解析导数可用时才使用 SLSQP 方法？

可重现的测试用例：

import numpy as np
import openmdao.api as om


class d1(om.ExplicitComponent):

        def setup(self):

            # Global design variables
            self.add_input('r', val= [3,3,3])          
            self.add_input('T', val= 20)
            
            # Coupling output
            self.add_output('M', val=0)
            self.add_output('cost', val=0)
            
        def setup_partials(self):
            # Finite difference all partials.
            self.declare_partials('*', '*', method='fd')

        def compute(self, inputs, outputs):
            # define inputs
            r = inputs['r']
            T = inputs['T'][0]
            
            cost =  174.42 * T * (r[0]**2 + 2*r[1]**2 + r[2]**2 + r[0]*r[1] + r[1]*r[2])
            
            M =     456.19 * T * (r[0]**2 + 2*r[1]**2 + r[2]**2 + r[0]*r[1] + r[1]*r[2]) - 599718
                 
            outputs['M'] = M
            outputs['cost'] = cost
            

class MDA(om.Group):

    class ObjCmp(om.ExplicitComponent):

        def setup(self):
            # Global Design Variable
            self.add_input('cost', val=0)

            # Output
            self.add_output('obj', val=0.0)

        def setup_partials(self):
            # Finite difference all partials.
            self.declare_partials('*', '*', method='fd')

        def compute(self, inputs, outputs):
            
            outputs['obj'] = inputs['cost']
           
            
    class ConCmp(om.ExplicitComponent):

        def setup(self):
            # Global Design Variable
            self.add_input('M', val=0)
            
            # Output
            self.add_output('con', val=0.0)

        def setup_partials(self):
            # Finite difference all partials.
            self.declare_partials('*', '*', method='fd')

        def compute(self, inputs, outputs):
            
            # assemble outputs
            outputs['con'] = inputs['M']

    def setup(self):
        
        self.add_subsystem('d1', d1(), promotes_inputs=['r','T'],
                            promotes_outputs=['M','cost'])
                                     
        self.add_subsystem('con_cmp', self.ConCmp(), promotes_inputs=['M'],
                            promotes_outputs=['con'])
                          
        self.add_subsystem('obj_cmp', self.ObjCmp(), promotes_inputs=['cost'],
                           promotes_outputs=['obj'])
                           

# Build the model
prob = om.Problem(model=MDA())
model = prob.model

model.add_design_var('r', lower= [3,3,3], upper= [10,10,10])
model.add_design_var('T', lower= 20, upper= 220)
model.add_objective('obj', scaler=1)
model.add_constraint('con', lower=0)

# Setup the optimization
prob.driver = om.ScipyOptimizeDriver(optimizer='SLSQP', tol=1e-3, disp=True)

prob.setup()
prob.set_solver_print(level=0)
prob.run_driver()

# Printout
print('minimum found at')
print(prob.get_val('T')[0])
print(prob.get_val('r'))

print('constraint')
print(prob.get_val('con')[0])

print('minimum objective')
print(prob.get_val('obj')[0])

Answer 1

根据您提供的测试用例，这里的问题是您的目标和约束的规模确实很差（您还有一些非常奇怪的编码选择......我修改了）。

运行OpenMDAO scaling report 表明您的目标值和约束值都在 1e6 左右：

这是相当大的，是您的问题的根源。一个（非常粗略的）经验法则是您的目标和约束应该在 1 阶左右。这不是硬性规定，但通常是一个很好的起点。有时，如果您有非常非常大或非常小的导数，其他缩放会更好地工作……但是 SQP 方法的某些部分对目标和约束值的缩放直接敏感。所以尽量让它们大致保持在 1 的范围内是个好主意。

将ref=1e6 添加到目标和约束中，为数值方法提供了足够的分辨率来收敛问题：

            Current function value: [0.229372]
            Iterations: 8
            Function evaluations: 8
            Gradient evaluations: 8
Optimization Complete
-----------------------------------
minimum found at
20.00006826587515
[3.61138704 3.         3.61138704]
constraint
197.20821903413162
minimum objective
229371.99547899762

这是我修改的代码（包括删除组内似乎没有做任何事情的额外类定义）：

import numpy as np
import openmdao.api as om


class d1(om.ExplicitComponent):

        def setup(self):

            # Global design variables
            self.add_input('r', val= [3,3,3])          
            self.add_input('T', val= 20)
            
            # Coupling output
            self.add_output('M', val=0)
            self.add_output('cost', val=0)
            
        def setup_partials(self):
            # Finite difference all partials.
            self.declare_partials('*', '*', method='cs')

        def compute(self, inputs, outputs):
            # define inputs
            r = inputs['r']
            T = inputs['T'][0]
            
            cost =  174.42 * T * (r[0]**2 + 2*r[1]**2 + r[2]**2 + r[0]*r[1] + r[1]*r[2])
            
            M =     456.19 * T * (r[0]**2 + 2*r[1]**2 + r[2]**2 + r[0]*r[1] + r[1]*r[2]) - 599718
                 
            outputs['M'] = M
            outputs['cost'] = cost
            

class MDA(om.Group):

    def setup(self):
        
        self.add_subsystem('d1', d1(), promotes_inputs=['r','T'],
                            promotes_outputs=['M','cost'])
                                     
        # self.add_subsystem('con_cmp', self.ConCmp(), promotes_inputs=['M'],
        #                     promotes_outputs=['con'])
                          
        # self.add_subsystem('obj_cmp', self.ObjCmp(), promotes_inputs=['cost'],
        #                    promotes_outputs=['obj'])
                           

# Build the model
prob = om.Problem(model=MDA())
model = prob.model

model.add_design_var('r', lower= [3,3,3], upper= [10,10,10])
model.add_design_var('T', lower= 20, upper= 220)
model.add_objective('cost', ref=1e6)
model.add_constraint('M', lower=0, ref=1e6)

# Setup the optimization
prob.driver = om.ScipyOptimizeDriver(optimizer='SLSQP', tol=1e-3, disp=True)

prob.setup()
prob.set_solver_print(level=0)

prob.set_val('r', 7.65)
prob.run_driver()

# Printout
print('minimum found at')
print(prob.get_val('T')[0])
print(prob.get_val('r'))

print('constraint')
print(prob.get_val('M')[0])

print('minimum objective')
print(prob.get_val('cost')[0])

Answer 2

您使用的是哪种 SLSQP 方法？在 pyOptSparse 中有一种实现，在 ScipyOptimizer 中有一种实现。 pyoptsparse 中的那个较旧，并且不尊重边界约束。 Scipy中的那个更新并且确实如此。 （是的，它们具有相同的名称并共享一些血统，这非常令人困惑......但不再是同一个优化器）

当您超出范围时不应引发分析错误。如果您需要遵守严格的界限，我建议在 pyoptsparse 中使用 IPopt（如果可以编译的话）或切换到 ScipyOptimizerDriver 及其 SLSQP 实现。

根据您的问题，如果您在谈论边界约束或不平等/平等约束，我并不完全清楚。如果是后者，那么没有任何优化器可以保证您始终处于可行区域。像 IPopt 这样的内点方法会更好地留在区域内，但不是 100% 的时间。

一般来说，使用基于梯度的优化非常关键，即使问题在约束区域之外，您也可以使问题平滑和连续。如果有部分空间是你绝对不能进入的，那么你需要将这些变量变成设计变量并使用绑定约束。这有时需要稍微重新制定您的问题公式，可能通过添加一种兼容性约束，即“设计变量 = 计算值”。然后，您可以确保将设计变量传递给任何要求该值严格在界限内的东西，并且（希望）一个收敛的答案也将满足您的兼容性约束。

如果您提供某种测试用例或示例，我可以用更具体的建议修改我的答案。