使用 openMDAO 处理 SLSQP 优化中的约束
Treatment of constraints in SLSQP optimization with openMDAO
对于 openMDAO，我正在使用 FD 导数并尝试使用 SLSQP 方法解决非线性约束优化问题。每当优化器到达违反其中一个约束的点时，它就会崩溃并显示消息：
Optimization FAILED. Positive directional derivative for linesearch
例如，如果我故意将初始点设置为不可行的设计点，优化器将执行 1 次迭代并退出并出现上述错误（当我从可行点开始时也会发生同样的情况，但随后优化器到达了几次迭代后的不可行点）。
根据 In OpenMDAO, is there a way to ensure that the constraints are respected before proceeding with a computation? 中问题的答案，我假设引发 AnalysisError 异常对我的情况不起作用，对吗？有没有其他方法可以防止优化器进入不可行的区域或至少在 linesearch 上回溯并尝试不同的 direction/distance？还是仅在解析导数可用时才使用 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])
您使用的是哪种 SLSQP 方法？ pyOptSparse 中有一种实现，ScipyOptimizer 中有一种实现。 pyoptsparse 中的那个比较旧，不遵守边界约束。 Scipy 中的那个较新并且确实如此。 （是的，这很令人困惑，因为它们具有相同的名称并共享一些血统......但不再是同一个优化器）
你不应该在超出范围时引发分析错误。如果您需要严格遵守边界，我建议在 pyoptsparse 中使用 IPopt（如果您可以编译）或切换到 ScipyOptimizerDriver 及其 SLSQP 实现。
根据你的问题，我不太清楚你是在谈论边界约束还是 inequality/equality 约束。如果是后者，那么就没有任何优化器可以保证您始终保持在可行区域中。像 IPopt 这样的内部点方法会更好地留在区域内，但不是 100% 的时间。
一般来说，使用基于梯度的优化，即使在约束区域之外，也能使问题平滑和连续，这一点非常重要。如果space中有部分是你绝对不能进去的，那你就需要把那些变量做成设计变量，使用bound constraints。这有时需要稍微重新制定您的问题公式，可能通过添加一种兼容性约束来表示“设计变量 = computed_value”。然后您可以确保将设计变量传递给任何需要值严格在界限内的东西，并且（希望）收敛的答案也将满足您的兼容性约束。
如果您提供某种测试用例或示例，我可以用更具体的建议修改我的答案。
根据您提供的测试用例，这里的问题是您的缩放 objective 和约束真的很差（您还有一些非常奇怪的编码选择......我修改了）。
运行 OpenMDAO scaling report 显示您的 objective 和约束值的大小都在 1e6 左右：
这非常大，是您问题的根源。一个（非常粗略的）经验法则是你的 objectives 和约束应该在顺序 1 附近。这不是硬性规定，但通常是一个很好的起点。有时其他缩放会更好，如果你有非常大或小的导数......但是 SQP 方法的某些部分对 objective 的缩放和约束值直接敏感。所以尽量让它们大致保持在 1 的范围内是个好主意。
向 objective 和约束添加 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
这是我修改的代码（包括删除您组中似乎没有做任何事情的额外 class 定义）：
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])
对于 openMDAO，我正在使用 FD 导数并尝试使用 SLSQP 方法解决非线性约束优化问题。每当优化器到达违反其中一个约束的点时，它就会崩溃并显示消息：
Optimization FAILED. Positive directional derivative for linesearch
例如，如果我故意将初始点设置为不可行的设计点，优化器将执行 1 次迭代并退出并出现上述错误（当我从可行点开始时也会发生同样的情况，但随后优化器到达了几次迭代后的不可行点）。
根据 In OpenMDAO, is there a way to ensure that the constraints are respected before proceeding with a computation? 中问题的答案，我假设引发 AnalysisError 异常对我的情况不起作用，对吗？有没有其他方法可以防止优化器进入不可行的区域或至少在 linesearch 上回溯并尝试不同的 direction/distance？还是仅在解析导数可用时才使用 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])
您使用的是哪种 SLSQP 方法？ pyOptSparse 中有一种实现，ScipyOptimizer 中有一种实现。 pyoptsparse 中的那个比较旧，不遵守边界约束。 Scipy 中的那个较新并且确实如此。 （是的，这很令人困惑，因为它们具有相同的名称并共享一些血统......但不再是同一个优化器）
你不应该在超出范围时引发分析错误。如果您需要严格遵守边界，我建议在 pyoptsparse 中使用 IPopt（如果您可以编译）或切换到 ScipyOptimizerDriver 及其 SLSQP 实现。
根据你的问题，我不太清楚你是在谈论边界约束还是 inequality/equality 约束。如果是后者，那么就没有任何优化器可以保证您始终保持在可行区域中。像 IPopt 这样的内部点方法会更好地留在区域内，但不是 100% 的时间。
一般来说，使用基于梯度的优化，即使在约束区域之外，也能使问题平滑和连续，这一点非常重要。如果space中有部分是你绝对不能进去的，那你就需要把那些变量做成设计变量，使用bound constraints。这有时需要稍微重新制定您的问题公式，可能通过添加一种兼容性约束来表示“设计变量 = computed_value”。然后您可以确保将设计变量传递给任何需要值严格在界限内的东西，并且（希望）收敛的答案也将满足您的兼容性约束。
如果您提供某种测试用例或示例，我可以用更具体的建议修改我的答案。
根据您提供的测试用例，这里的问题是您的缩放 objective 和约束真的很差（您还有一些非常奇怪的编码选择......我修改了）。
运行 OpenMDAO scaling report 显示您的 objective 和约束值的大小都在 1e6 左右：
这非常大，是您问题的根源。一个（非常粗略的）经验法则是你的 objectives 和约束应该在顺序 1 附近。这不是硬性规定，但通常是一个很好的起点。有时其他缩放会更好，如果你有非常大或小的导数......但是 SQP 方法的某些部分对 objective 的缩放和约束值直接敏感。所以尽量让它们大致保持在 1 的范围内是个好主意。
向 objective 和约束添加 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
这是我修改的代码（包括删除您组中似乎没有做任何事情的额外 class 定义）：
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])