如何使用 OpenMDAO 中的两个隐式且可能耦合的系统进行优化？
How to optimize with two implicit, and possibly coupled, systems in OpenMDAO?
我一直在参考 Sellar MDA 和 BalanceComp 教程，但我无法制定 OpenMDAO 架构来优化类型隐式方程的耦合系统：
Minimize the function F with design variables x, y such that the following nonlinear equations are simultaneously satisfied as constraints:
R(x, y, z) = f(x, y, c) - g(x, y, d) = 0
S(x, y, z) = f(x, y, c) - h(x, y, d) = e
with constants c, d, e.
所以我的想法是让F, f, g, h成为ExplicitComponent类，让R, S成为Group类，使用BalanceComp（每个里面一个！）来确定它们各自的解决方案，因为这些系统可能会在其他应用程序中独立使用。我相信，常数 c, d, e 可以通过 IndepVarComp 提供。
一起解决系统涉及创建一个 Group 组件（称之为 ImpMDA），其中 R, S 使用 Group 像 cycle 来自 Sellar 教程。
我的主要问题是如何在循环中跨两个系统全局“同步”x, y。是否应该将 R 的结果提供给 S，反之亦然？由于它们是同时求解的，所以我认为它们是耦合的。
请注意，x, y 可能来自其他一些独立的组件，因此很难将它们作为向量提供，除非创建了一些 ExplicitComponent，如 Lister，它接受输入并生成一个向量以馈送到 ImpMDA.
虽然在更简单的玩具箱（如果可能的话）上的具体实现会非常有帮助，但如果我能对在 OpenMDAO 中建模这个问题有一些一般性的了解，我将非常感激。
对于特定的用例，在尝试保持俯仰 trim（隐式系统 1) 和飞行中所有操作点的俯仰稳定性（隐式系统 2），从而确定水平尾翼面积和机翼位置。与重心匹配的压力中心（取决于机翼面积）是俯仰的要求trim，同样对于俯仰稳定的中性点，因此设计变量在系统中发挥作用。
编辑：
这是我正在寻找的模块化方法的 Python 示例：
import numpy as np
from scipy.optimize import fsolve, minimize
def f(xs, c): return -xs[0]**2 + xs[1] + c # ExplicitComponent 1
def g(xs, d): return -xs[0]**3 + xs[1]**2 - d # ExplicitComponent 2
def h(xs, d): return xs[0] + xs[1]**2 + d # ExplicitComponent 3
# This function is to quickly generate residual functions shown in the question.
# The OpenMDAO implementations for different residual functions will be done either
# via BalanceComp or ImplicitComponent, and may not share the same equation "setup" as here.
def residualer(f1, f2, a=6, b=4, const=0):
return lambda xs: f1(xs, a) - f2(xs, b) + const
def F(xs): # Objective function
return xs[0]**2 - xs[1]**2 - 5
x0 = [-2, 2] # Initial guess
# Problem 1
x1_opt = minimize(F, x0,
method='trust-constr',
constraints={ 'type': 'eq', 'fun': residualer(f, g, a=1, b=2) },
options = { 'disp': False } )
print(x1_opt, '\n\n')
# Problem 2
m, n = 5, 6
def coupled(xs): # Coupling the equations
return [residualer(f, g)(xs), residualer(f, h, a=4, b=3, const=m*n)(xs)]
x2_opt = minimize(F, x0,
method='trust-constr',
constraints={ 'type': 'eq', 'fun': coupled },
options = { 'disp': False } )
print(x2_opt)
哪个returns:
cg_niter: 22
cg_stop_cond: 2
constr: [array([-8.8817842e-16])]
constr_nfev: [102]
constr_nhev: [0]
constr_njev: [0]
constr_penalty: 3.2410415627552283
constr_violation: 8.881784197001252e-16
execution_time: 0.028232574462890625
fun: -10.302775637731996
grad: array([ 1.19209290e-07, -4.60555129e+00])
jac: [array([[-5.96046448e-08, -3.60555129e+00]])]
lagrangian_grad: array([ 1.95345288e-07, -8.88178420e-16])
message: '`xtol` termination condition is satisfied.'
method: 'equality_constrained_sqp'
nfev: 102
nhev: 0
nit: 23
niter: 23
njev: 0
optimality: 1.9534528835387254e-07
status: 2
success: True
tr_radius: 5.238777835236851e-09
v: [array([-1.2773501])]
x: array([1.17707733e-08, 2.30277564e+00])
cg_niter: 0
cg_stop_cond: 1
constr: [array([4.4408921e-14, 0.0000000e+00])]
constr_nfev: [36]
constr_nhev: [0]
constr_njev: [0]
constr_penalty: 2.1867764273118655
constr_violation: 4.440892098500626e-14
execution_time: 0.012656688690185547
fun: -24.59492687150548
grad: array([ 5.27636947, -10.30629829])
jac: [array([[15.60368701, -9.30629829],
[-6.27636948, -9.30629823]])]
lagrangian_grad: array([1.77635684e-15, 0.00000000e+00])
message: '`gtol` termination condition is satisfied.'
method: 'equality_constrained_sqp'
nfev: 36
nhev: 0
nit: 12
niter: 12
njev: 0
optimality: 1.7763568394002505e-15
status: 1
success: True
tr_radius: 4.9172480000000025
v: [array([-0.55882674, -0.54862737])]
x: array([2.63818474, 5.1531491 ])
我试着以模块化的方式来做这件事。由于脚本相当大，它在 github 要点 here.
这里的 Residualer 是一个为你的 residualer 函数服务的组。在它的初始化方法中，我定义了特定于每个实例化的选项。
每个实例化都为 f、g 和 h 创建自己的子系统。实际上，我可能会制作一个输出 f、g 和 h 的组件，因为它们是简单的方程式并且彼此独立。
problem1_optimizer() 和 problem2_optimizer() 使用 scipy 优化器对残差进行等式约束，就像您的示例一样。
problem2_optimizer 的 N2 是：
problem2_solver() 在 residualer 上使用 use_solver 选项来包含平衡组件和 NewtonSolver。这个平衡组件为x中的每个元素添加一个隐变量，它有一个对应的残差（R或S）。这两个标量值（x1 和 x2）被混合到 x 变量中，然后通过提升传递给适当的组件。按照编码，这要求我们对每个标量隐式变量都有一个残差，因此它只适用于你的第二个例子。请注意，以下未连接的输入是由于我们的平衡补偿仅采用一个变量并在等式的另一侧使用假设值 0.0。在第二种情况下没有优化，因为只有两个变量和两个未知数，因此没有更多的优化自由度。
我一直在参考 Sellar MDA 和 BalanceComp 教程，但我无法制定 OpenMDAO 架构来优化类型隐式方程的耦合系统：
Minimize the function
Fwith design variablesx, ysuch that the following nonlinear equations are simultaneously satisfied as constraints:
R(x, y, z) = f(x, y, c) - g(x, y, d) = 0
S(x, y, z) = f(x, y, c) - h(x, y, d) = ewith constants
c, d, e.
所以我的想法是让F, f, g, h成为ExplicitComponent类，让R, S成为Group类，使用BalanceComp（每个里面一个！）来确定它们各自的解决方案，因为这些系统可能会在其他应用程序中独立使用。我相信，常数 c, d, e 可以通过 IndepVarComp 提供。
一起解决系统涉及创建一个 Group 组件（称之为 ImpMDA），其中 R, S 使用 Group 像 cycle 来自 Sellar 教程。
我的主要问题是如何在循环中跨两个系统全局“同步”x, y。是否应该将 R 的结果提供给 S，反之亦然？由于它们是同时求解的，所以我认为它们是耦合的。
请注意，x, y 可能来自其他一些独立的组件，因此很难将它们作为向量提供，除非创建了一些 ExplicitComponent，如 Lister，它接受输入并生成一个向量以馈送到 ImpMDA.
虽然在更简单的玩具箱（如果可能的话）上的具体实现会非常有帮助，但如果我能对在 OpenMDAO 中建模这个问题有一些一般性的了解，我将非常感激。
对于特定的用例，在尝试保持俯仰 trim（隐式系统 1) 和飞行中所有操作点的俯仰稳定性（隐式系统 2），从而确定水平尾翼面积和机翼位置。与重心匹配的压力中心（取决于机翼面积）是俯仰的要求trim，同样对于俯仰稳定的中性点，因此设计变量在系统中发挥作用。
编辑：
这是我正在寻找的模块化方法的 Python 示例：
import numpy as np
from scipy.optimize import fsolve, minimize
def f(xs, c): return -xs[0]**2 + xs[1] + c # ExplicitComponent 1
def g(xs, d): return -xs[0]**3 + xs[1]**2 - d # ExplicitComponent 2
def h(xs, d): return xs[0] + xs[1]**2 + d # ExplicitComponent 3
# This function is to quickly generate residual functions shown in the question.
# The OpenMDAO implementations for different residual functions will be done either
# via BalanceComp or ImplicitComponent, and may not share the same equation "setup" as here.
def residualer(f1, f2, a=6, b=4, const=0):
return lambda xs: f1(xs, a) - f2(xs, b) + const
def F(xs): # Objective function
return xs[0]**2 - xs[1]**2 - 5
x0 = [-2, 2] # Initial guess
# Problem 1
x1_opt = minimize(F, x0,
method='trust-constr',
constraints={ 'type': 'eq', 'fun': residualer(f, g, a=1, b=2) },
options = { 'disp': False } )
print(x1_opt, '\n\n')
# Problem 2
m, n = 5, 6
def coupled(xs): # Coupling the equations
return [residualer(f, g)(xs), residualer(f, h, a=4, b=3, const=m*n)(xs)]
x2_opt = minimize(F, x0,
method='trust-constr',
constraints={ 'type': 'eq', 'fun': coupled },
options = { 'disp': False } )
print(x2_opt)
哪个returns:
cg_niter: 22
cg_stop_cond: 2
constr: [array([-8.8817842e-16])]
constr_nfev: [102]
constr_nhev: [0]
constr_njev: [0]
constr_penalty: 3.2410415627552283
constr_violation: 8.881784197001252e-16
execution_time: 0.028232574462890625
fun: -10.302775637731996
grad: array([ 1.19209290e-07, -4.60555129e+00])
jac: [array([[-5.96046448e-08, -3.60555129e+00]])]
lagrangian_grad: array([ 1.95345288e-07, -8.88178420e-16])
message: '`xtol` termination condition is satisfied.'
method: 'equality_constrained_sqp'
nfev: 102
nhev: 0
nit: 23
niter: 23
njev: 0
optimality: 1.9534528835387254e-07
status: 2
success: True
tr_radius: 5.238777835236851e-09
v: [array([-1.2773501])]
x: array([1.17707733e-08, 2.30277564e+00])
cg_niter: 0
cg_stop_cond: 1
constr: [array([4.4408921e-14, 0.0000000e+00])]
constr_nfev: [36]
constr_nhev: [0]
constr_njev: [0]
constr_penalty: 2.1867764273118655
constr_violation: 4.440892098500626e-14
execution_time: 0.012656688690185547
fun: -24.59492687150548
grad: array([ 5.27636947, -10.30629829])
jac: [array([[15.60368701, -9.30629829],
[-6.27636948, -9.30629823]])]
lagrangian_grad: array([1.77635684e-15, 0.00000000e+00])
message: '`gtol` termination condition is satisfied.'
method: 'equality_constrained_sqp'
nfev: 36
nhev: 0
nit: 12
niter: 12
njev: 0
optimality: 1.7763568394002505e-15
status: 1
success: True
tr_radius: 4.9172480000000025
v: [array([-0.55882674, -0.54862737])]
x: array([2.63818474, 5.1531491 ])
我试着以模块化的方式来做这件事。由于脚本相当大，它在 github 要点 here.
这里的 Residualer 是一个为你的 residualer 函数服务的组。在它的初始化方法中，我定义了特定于每个实例化的选项。
每个实例化都为 f、g 和 h 创建自己的子系统。实际上，我可能会制作一个输出 f、g 和 h 的组件，因为它们是简单的方程式并且彼此独立。
problem1_optimizer() 和 problem2_optimizer() 使用 scipy 优化器对残差进行等式约束，就像您的示例一样。
problem2_optimizer 的 N2 是：
problem2_solver() 在 residualer 上使用 use_solver 选项来包含平衡组件和 NewtonSolver。这个平衡组件为x中的每个元素添加一个隐变量，它有一个对应的残差（R或S）。这两个标量值（x1 和 x2）被混合到 x 变量中，然后通过提升传递给适当的组件。按照编码，这要求我们对每个标量隐式变量都有一个残差，因此它只适用于你的第二个例子。请注意，以下未连接的输入是由于我们的平衡补偿仅采用一个变量并在等式的另一侧使用假设值 0.0。在第二种情况下没有优化，因为只有两个变量和两个未知数，因此没有更多的优化自由度。