如何在 OpenMDAO 中将方程组求解为数组结构？

Question

我尝试使用 OpenMDAO 求解下面显示的隐式方程。等式如下所示，

x[i]*z[i] + z[i] - 4 = 0,  
y[i] = x[i] + 2*z[i]

The solution is (for i=2) z = [2.666667, 2], y = [5.833333, 5] for x = [0.5, 1.0].

对于这种情况，我使用了如下所示的代码，

from __future__ import print_function
from openmdao.api import Component, Group, Problem, Newton, ScipyGMRES, IndepVarComp

class SimpleEquationSystem(Component):
    """Solve the Equation 
            x[i]*z[i] + z[i] - 4 = 0
            y[i] = x[i] + 2*z[i]

       Solution: z = [2.666667, 2], y = [5.833333, 5] for x = [0.5, 1.0]
    """

    def __init__(self):
        super(SimpleEquationSystem, self).__init__()


        self.add_param('x', [0.5, 1.0])
        self.add_state('y', [0.0,0.0])
        self.add_state('z', [0.0,0.0])
        self.iter=0

    def solve_nonlinear(self, params, unknowns, resids):
        """This component does no calculation on its own. It mainly holds the
        initial value of the state. An OpenMDAO solver outside of this
        component varies it to drive the residual to zero."""
        pass

    def apply_nonlinear(self, params, unknowns, resids):
        """ Report the residual """
        self.iter+=1
        for i in range(2):
            x=params['x'][i]
            y = unknowns['y'][i]
            z = unknowns['z'][i]


            resids['y'][i] = x*z + z - 4
            resids['z'][i] = x + 2*z - y

        print('y_%d' % self.iter,'=%s' %resids['y'], 'z_%d' % self.iter, '=%s' %resids['z'])
        print('x' ,'=%s' %x, 'y', '=%s' %y, 'z', '=%s' %z)

top = Problem()
root = top.root = Group()
root.add('comp', SimpleEquationSystem())
root.add('p1', IndepVarComp('x', [0.5, 1.0]))
root.connect('p1.x', 'comp.x')
# Tell these components to finite difference
root.comp.deriv_options['type'] = 'fd'
root.comp.deriv_options['form'] = 'central'
root.comp.deriv_options['step_size'] = 1.0e-4
root.nl_solver = Newton()
root.nl_solver.options['maxiter']=int(200)
root.ln_solver = ScipyGMRES()

top.setup()
top.print_all_convergence(level=1, depth=2)
top.run()
print('Solution x=%s, y=%s, z=%s' % (top['comp.x'], top['comp.y'], top['comp.z']))

此代码出错并显示消息 "RuntimeWarning: invalid value encountered in double_scalars"。使用上述代码片段时，应该能够在 OpenMDAO 中得到此错误。

当残差方程被实现为标量而不是向量时，如下所示，它工作正常。

resids['z1']= params['x1'] + 2*unknowns['z1'] - unknowns['y1']
resids['z2']= params['x2'] + 2*unknowns['z2'] - unknowns['y2']
resids['y1']= params['x1']*unknowns['z1'] + unknowns['z1'] - 4
resids['y2']= params['x2']*unknowns['z2'] + unknowns['z2'] - 4

但我希望将残差作为向量，以便使用 for 循环更容易处理。你能帮我解决这个问题吗？

Answer 1

当你声明一个变量时，它就是你声明它的变量类型。你给了它一个列表，OpenMDAO 将其解释为一个列表，它不是支持微分的数据类型，所以 Newton 不能用它做任何事情。您需要通过进行以下更改使它们成为 numpy 数组：

import numpy as np

和

self.add_param('x', np.array([0.5, 1.0]))
self.add_state('y', np.array([0.0,0.0]))
self.add_state('z', np.array([0.0,0.0]))

还有：

root.add('p1', IndepVarComp('x', np.array([0.5, 1.0])))

如何在 OpenMDAO 中将方程组求解为数组结构？

How to solve an equation system as an array structure in OpenMDAO?

python

openmdao