OpenMDAO使用说明-多目标优化

OpenMDAO多目标优化

Case:Sellar - A Two-Discipline Problem with a Nonlinear Solver

link.
现在，我们将解决一个稍微复杂的问题，涉及两个学科，因此涉及两个主要的component。 您将学习如何将组件组合到一个更大的模型中，以及如何使用NonlinearBlockGaussSeidel非线性求解器将具有耦合组件的组收敛。
Sellar问题是一个两学科的问题，每个学科都用一个方程式描述。 这些方程式没有任何物理意义，它只是一个示例，可用来理解简单的耦合模型。 每个组件的输出都传入另一个组件的输入，这会创建一个耦合模型，该模型需要收敛才能使输出有效。 您可以在下面描述问题结构的XDSM图中通过y1和y2变量看到两个学科之间的耦合：

1、 构建各学科的component

在以下组件定义中，setup方法中有一个对declare_partials的调用，如下所示：

self.declare_partials('*', '*', method='fd')

构建两个学科的component

class SellarDis1(om.ExplicitComponent):
    """
    Component containing Discipline 1 -- no derivatives version.
    """

    def setup(self):

        # Global Design Variable
        self.add_input('z', val=np.zeros(2))

        # Local Design Variable
        self.add_input('x', val=0.)

        # Coupling parameter
        self.add_input('y2', val=1.0)

        # Coupling output
        self.add_output('y1', val=1.0)

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

    def compute(self, inputs, outputs):
        """
        Evaluates the equation
        y1 = z1**2 + z2 + x1 - 0.2*y2
        """
        z1 = inputs['z'][0]
        z2 = inputs['z'][1]
        x1 = inputs['x']
        y2 = inputs['y2']

        outputs['y1'] = z1**2 + z2 + x1 - 0.2*y2
class SellarDis2(om.ExplicitComponent):
    """
    Component containing Discipline 2 -- no derivatives version.
    """

    def setup(self):
        # Global Design Variable
        self.add_input('z', val=np.zeros(2))

        # Coupling parameter
        self.add_input('y1', val=1.0)

        # Coupling output
        self.add_output('y2', val=1.0)

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

    def compute(self, inputs, outputs):
        """
        Evaluates the equation
        y2 = y1**(.5) + z1 + z2
        """

        z1 = inputs['z'][0]
        z2 = inputs['z'][1]
        y1 = inputs['y1']

        # Note: this may cause some issues. However, y1 is constrained to be
        # above 3.16, so lets just let it converge, and the optimizer will
        # throw it out
        if y1.real < 0.0:
            y1 *= -1

        outputs['y2'] = y1**.5 + z1 + z2

2、 Grouping and Connecting Components

现在，我们要在OpenMDAO中构建由上面的XDSM图表示的模型。 我们已经定义了两个学科组成部分，但是仍然需要计算模型的三个输出。 另外，由于我们将计算分为多个部分，因此需要将它们全部分组，并告诉OpenMDAO如何在它们之间传递数据。

class SellarMDA(om.Group):
    """
    Group containing the Sellar MDA.
    """

    def setup(self):
        indeps = self.add_subsystem('indeps', om.IndepVarComp(), promotes=['*'])
        indeps.add_output('x', 1.0)
        indeps.add_output('z', np.array([5.0, 2.0]))

        cycle = self.add_subsystem('cycle', om.Group(), promotes=['*'])
        cycle.add_subsystem('d1', SellarDis1(), promotes_inputs=['x', 'z', 'y2'],
                            promotes_outputs=['y1'])
        cycle.add_subsystem('d2', SellarDis2(), promotes_inputs=['z', 'y1'],
                            promotes_outputs=['y2'])

        # Nonlinear Block Gauss Seidel is a gradient free solver
        cycle.nonlinear_solver = om.NonlinearBlockGS()

        self.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
                                                  z=np.array([0.0, 0.0]), x=0.0),
                           promotes=['x', 'z', 'y1', 'y2', 'obj'])

        self.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])
        self.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2'])

我们可以创建一个新的类：Group。 这是使您可以构建复杂的模型层次结构的容器。 Group可以包含其他group，component或group和component的组合。

您可以直接创建要使用的Group实例，也可以从其子类定义您自己的自定义组。 我们正在上面做这两个事情。 首先，我们定义了自己的自定义Group子类，称为SellarMDA。 在我们的运行脚本中，我们创建了一个SellarMDA实例来实际运行它。 然后，在SellarMDA的设置方法中，我们还通过添加子系统循环来直接处理Group实例：

cycle = self.add_subsystem('cycle', om.Group(), promotes=['*'])
cycle.add_subsystem('d1', SellarDis1(), promotes_inputs=['x', 'z', 'y2'], promotes_outputs=['y1'])
cycle.add_subsystem('d2', SellarDis2(), promotes_inputs=['z', 'y1'], promotes_outputs=['y2'])

# Nonlinear Block Gauss-Seidel is a gradient-free solver
cycle.nonlinear_solver = NonlinearBlockGS()

实例化后，我们的SellarMDA组将具有一个三级层次结构，其自身为最高层次，如下图所示：

3、Why do we create the cycle subgroup?

d1和d2之间存在循环数据相关性，需要使用非线性求解器对其进行收敛才能获得有效答案。 将这两个组件放到各自的子组中会更有效率，这样我们就可以自行迭代地对其进行收敛，然后再继续进行模型中的其他计算。 其中包含循环的模型通常称为“多学科分析”或MDA。 您可以选择要使用哪种求解器来收敛MDA。 最常见的选择是：
NonlinearBlockGaussSeidel
NewtonSolver
NonlinearBlockGaussSeidel求解器，有时也称为“定点迭代求解器”，是一种无梯度方法，在许多情况下都可以很好地工作。 耦合更紧密的问题，或者ImplicitComponent实例的问题未实现自己的solve_nonlinear方法，则需要使用牛顿求解器。
名为cycle的subgroup在这里很有用，因为它包含塞拉尔问题的多学科耦合。 这使我们可以将非线性求解器分配给循环以仅收敛这两个分量，然后再继续进行obj_cmp，con_cmp1和con_cmp2的最终计算以计算问题的实际输出。

4、Promoting variables with the same name connects them

此模型中的数据连接是通过promotion建立的。 OpenMDAO将查看层次结构的每个级别，并连接所有具有相同名称的输出-输入对。

5、ExecComp is a helper component for quickly defining components for simple equations

在模型中，很多时候，您需要定义一个新变量作为其他变量的简单函数。 OpenMDAO提供了一个辅助组件，称为ExecComp，可简化此工作。 它相当灵活，可让您处理标量或数组，处理单位以及调用基本数学函数（例如sin或exp）。 我们在此模型中使用了ExecComp来计算目标和约束。

Method:Linking Variables with Promotion vs. Connection

在上一教程中，我们使用两个学科的component和一些ExecComps建立了Sellar问题的模型。 为了使OpenMDAO在所有组件之间传递数据，我们使用提升的变量将所有内容链接在一起，以使数据从输出传递到具有相同提升名称的输入。

提取变量通常是建立从输出到输入的数据传递链接的便捷方法。 但是，您也可以使用对connect方法的调用，以将输出链接到输入，而无需进行任何升级。 使用以下方法定义相同的Sellar模型：
Variable promotion
Connect statements
Both variable promotion and connect statements
这三个方法都将给出完全相同的答案，但是每个变量的寻址方式都将略有不同。

Variable Promotion

添加子系统时，可以提取输入和输出变量：

import numpy as np

import openmdao.api as om
from openmdao.test_suite.components.sellar import SellarDis1, SellarDis2

class SellarMDA(om.Group):
    """
    Group containing the Sellar MDA.
    """

    def setup(self):
        indeps = self.add_subsystem('indeps', om.IndepVarComp(), promotes=['*'])
        indeps.add_output('x', 1.0)
        indeps.add_output('z', np.array([5.0, 2.0]))

        cycle = self.add_subsystem('cycle', om.Group(), promotes=['*'])
        cycle.add_subsystem('d1', SellarDis1(),
                            promotes_inputs=['x', 'z', 'y2'],
                            promotes_outputs=['y1'])
        cycle.add_subsystem('d2', SellarDis2(),
                            promotes_inputs=['z', 'y1'],
                            promotes_outputs=['y2'])

        # Nonlinear Block Gauss Seidel is a gradient free solver
        cycle.nonlinear_solver = om. NonlinearBlockGS()

        self.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
                                                  z=np.array([0.0, 0.0]), x=0.0),
                           promotes=['x', 'z', 'y1', 'y2', 'obj'])

        self.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'),
                           promotes=['con1', 'y1'])
        self.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'),
                           promotes=['con2', 'y2'])


prob = om.Problem()
prob.model = SellarMDA()

prob.setup()

prob['x'] = 2.
prob['z'] = [-1., -1.]

prob.run_model()
=====
cycle
=====
NL: NLBGS Converged in 9 iterations

print((prob['y1'][0], prob['y2'][0], prob['obj'][0], prob['con1'][0], prob['con2'][0]))

(2.10951650608596, -0.5475825303701556, 6.838584498046503, 1.0504834939140402, -24.547582530370157)

另外，配置group时可以提取变量：

import numpy as np

import openmdao.api as om
from openmdao.test_suite.components.sellar import SellarDis1, SellarDis2

class SellarMDA(om.Group):
    """
    Group containing the Sellar MDA.
    """

    def setup(self):
        # set up model hierarchy
        indeps = self.add_subsystem('indeps', om.IndepVarComp(), promotes=['*'])
        indeps.add_output('x', 1.0)
        indeps.add_output('z', np.array([5.0, 2.0]))

        cycle = self.add_subsystem('cycle', om.Group())
        cycle.add_subsystem('d1', SellarDis1())
        cycle.add_subsystem('d2', SellarDis2())

        cycle.nonlinear_solver = om. NonlinearBlockGS()

        self.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
                                                  z=np.array([0.0, 0.0]), x=0.0))

        self.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'))
        self.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'))

    def configure(self):
        # connect everything via promotes
        self.cycle.promotes('d1', inputs=['x', 'z', 'y2'], outputs=['y1'])
        self.cycle.promotes('d2', inputs=['z', 'y1'], outputs=['y2'])
        self.promotes('cycle', any=['*'])

        self.promotes('obj_cmp', any=['x', 'z', 'y1', 'y2', 'obj'])
        self.promotes('con_cmp1', any=['con1', 'y1'])
        self.promotes('con_cmp2', any=['con2', 'y2'])


prob = om.Problem()
prob.model = SellarMDA()

prob.setup()

prob['x'] = 2.
prob['z'] = [-1., -1.]

prob.run_model()
=====
cycle
=====
NL: NLBGS Converged in 9 iterations

print((prob['y1'][0], prob['y2'][0], prob['obj'][0], prob['con1'][0], prob['con2'][0]))

(2.10951650608596, -0.5475825303701556, 6.838584498046503, 1.0504834939140402, -24.547582530370157)

区别就是前一个是在添加子系统的时候promote了变量，这个是先添加子系统，再用self.promote设置变量。

有一些重要的细节要注意：

不能有两个具有相同名称的输出；可以将多个输入升级为相同的名称，但是要建立连接，还必须有一个具有相同名称的输出。否则，将无法建立连接。

您可以使用全局模式来提升许多变量，而无需全部指定它们，但请少使用promotes = [’’]。尽管这似乎是一种方便的处理方法，但是它可能会使正在阅读您的代码的其他人难以理解哪些变量相互关联。如果在不会引起混淆的情况下（例如与周期混淆），可以使用promots = [’’]，这种情况仅存在于非线性求解器可以收敛两个分量的情况下。何时可以安全使用promots = [’*’]的另一个示例是：有ExecComps来确保可以清楚地知道该组件的I / O是什么。

Connect Statements

使用connect语句而不是Promotions可以获得完全相同的模型结果。 但是，请注意在这些connect和print语句中如何处理变量。

import numpy as np

import openmdao.api as om
from openmdao.test_suite.components.sellar import SellarDis1, SellarDis2

class SellarMDAConnect(om.Group):
    """
    Group containing the Sellar MDA. This version uses the disciplines without derivatives.
    """

    def setup(self):
        indeps = self.add_subsystem('indeps', om.IndepVarComp())
        indeps.add_output('x', 1.0)
        indeps.add_output('z', np.array([5.0, 2.0]))

        cycle = self.add_subsystem('cycle', om.Group())
        cycle.add_subsystem('d1', SellarDis1())
        cycle.add_subsystem('d2', SellarDis2())
        cycle.connect('d1.y1', 'd2.y1')
        cycle.connect('d2.y2', 'd1.y2')

        # Nonlinear Block Gauss Seidel is a gradient free solver
        cycle.nonlinear_solver = om.NonlinearBlockGS()

        self.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
                                                  z=np.array([0.0, 0.0]), x=0.0))

        self.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'))
        self.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'))

        self.connect('indeps.x', ['cycle.d1.x', 'obj_cmp.x'])
        self.connect('indeps.z', ['cycle.d1.z', 'cycle.d2.z', 'obj_cmp.z'])
        self.connect('cycle.d1.y1', ['obj_cmp.y1', 'con_cmp1.y1'])
        self.connect('cycle.d2.y2', ['obj_cmp.y2', 'con_cmp2.y2'])

prob = om.Problem()
prob.model = SellarMDAConnect()

prob.setup()

prob['indeps.x'] = 2.
prob['indeps.z'] = [-1., -1.]

prob.run_model()
=====
cycle
=====
NL: NLBGS Converged in 9 iterations

print((prob['cycle.d1.y1'][0], prob['cycle.d2.y2'][0], prob['obj_cmp.obj'][0], prob['con_cmp1.con1'][0], prob['con_cmp2.con2'][0]))

(2.10951650608596, -0.5475825303701556, 6.838584498046503, 1.0504834939140402, -24.547582530370157)

Variable Promotion and Connect Statements

也可以在单个模型中组合promotion和connect。 在这里，请注意，我们不必在任何内容前添加“循环”，因为我们从这个group中提取了所有变量。

import numpy as np

import openmdao.api as om
from openmdao.test_suite.components.sellar import SellarDis1, SellarDis2

class SellarMDAPromoteConnect(om.Group):
    """
    Group containing the Sellar MDA. This version uses the disciplines without derivatives.
    """

    def setup(self):
        indeps = self.add_subsystem('indeps', om.IndepVarComp(), promotes=['*'])
        indeps.add_output('x', 1.0)
        indeps.add_output('z', np.array([5.0, 2.0]))

        cycle = self.add_subsystem('cycle', om.Group(), promotes=['*'])
        cycle.add_subsystem('d1', SellarDis1())
        cycle.add_subsystem('d2', SellarDis2())
        cycle.connect('d1.y1', 'd2.y1')
        cycle.connect('d2.y2', 'd1.y2')

        # Nonlinear Block Gauss Seidel is a gradient free solver
        cycle.nonlinear_solver = om.NonlinearBlockGS()

        self.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
                                                  z=np.array([0.0, 0.0]), x=0.0))

        self.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'))
        self.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'))

        self.connect('x', ['d1.x', 'obj_cmp.x'])
        self.connect('z', ['d1.z', 'd2.z', 'obj_cmp.z'])
        self.connect('d1.y1', ['con_cmp1.y1', 'obj_cmp.y1'])
        self.connect('d2.y2', ['con_cmp2.y2', 'obj_cmp.y2'])


prob = om.Problem()
prob.model = SellarMDAPromoteConnect()

prob.setup()

prob['x'] = 2.
prob['z'] = [-1., -1.]

prob.run_model()
=====
cycle
=====
NL: NLBGS Converged in 9 iterations

print((prob['d1.y1'][0], prob['d2.y2'][0], prob['obj_cmp.obj'][0], prob['con_cmp1.con1'][0], prob['con_cmp2.con2'][0]))

(2.10951650608596, -0.5475825303701556, 6.838584498046503, 1.0504834939140402, -24.547582530370157)

Optimizing the Sellar Problem

在先前的教程中，我们向您展示了如何定义Sellar模型并直接运行它。 现在，让我们看看如何优化Sellar问题以最小化目标函数。 这是塞拉优化问题的数学问题公式：

前面我们按以下方式构建了Sellar模型：

class SellarMDA(om.Group):
    """
    Group containing the Sellar MDA.
    """

    def setup(self):
        indeps = self.add_subsystem('indeps', om.IndepVarComp(), promotes=['*'])
        indeps.add_output('x', 1.0)
        indeps.add_output('z', np.array([5.0, 2.0]))

        cycle = self.add_subsystem('cycle', om.Group(), promotes=['*'])
        cycle.add_subsystem('d1', SellarDis1(), promotes_inputs=['x', 'z', 'y2'],
                            promotes_outputs=['y1'])
        cycle.add_subsystem('d2', SellarDis2(), promotes_inputs=['z', 'y1'],
                            promotes_outputs=['y2'])

        # Nonlinear Block Gauss Seidel is a gradient free solver
        cycle.nonlinear_solver = om.NonlinearBlockGS()

        self.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
                                                  z=np.array([0.0, 0.0]), x=0.0),
                           promotes=['x', 'z', 'y1', 'y2', 'obj'])

        self.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])
        self.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2'])

优化所需的变量均已存在。 因此，现在我们只需要运行脚本来执行优化。

import openmdao.api as om
from openmdao.test_suite.components.sellar_feature import SellarMDA

prob = om.Problem()
prob.model = SellarMDA()

prob.driver = om.ScipyOptimizeDriver()
prob.driver.options['optimizer'] = 'SLSQP'
# prob.driver.options['maxiter'] = 100
prob.driver.options['tol'] = 1e-8

prob.model.add_design_var('x', lower=0, upper=10)
prob.model.add_design_var('z', lower=0, upper=10)
prob.model.add_objective('obj')
prob.model.add_constraint('con1', upper=0)
prob.model.add_constraint('con2', upper=0)

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

# Ask OpenMDAO to finite-difference across the model to compute the gradients for the optimizer
prob.model.approx_totals()

prob.run_driver()

print('minimum found at')
print(prob['x'][0])
print(prob['z'])

print('minumum objective')
print(prob['obj'][0])

Optimization terminated successfully.    (Exit mode 0)
            Current function value: 3.183393951729169
            Iterations: 6
            Function evaluations: 6
            Gradient evaluations: 6
Optimization Complete
-----------------------------------
minimum found at
0.0
[1.97763888e+00 8.83056605e-15]
minumum objective
3.183393951729169

Controlling the Solver Print Output

请注意对prob.set_solver_print（）的调用，该调用将求解器的输出设置为0级，在该设置中，只有在求解器收敛失败时才会通知您。 有很多方法可以在模型中配置求解器打印输出以适合您的需求。

Approximate the total derivatives with finite difference

在这个案例中，我们使用的是SLSQP算法，这是一种基于梯度的优化方法。 到目前为止，我们的组件都没有提供任何解析导数，因此我们将对整个模型进行有限差分以近似导数。 这是通过以下代码行完成的：

prob.model.approx_totals()

为了简单起见，我们在这里使用有限差分，但是对于较大的模型，有限差分会导致较高的计算成本，并且精度可能会受到限制。 在模型中使用解析导数会更好。 您可以在《高级用户指南》中了解有关此内容的更多信息。