如何将方程式系统求解为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']))

此代码错误带有消息" Runtime Warning：在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

声明变量时，它正是您声明的变量类型。您给了它一个列表，该列表将其解释为列表，这不是支持差异化的数据类型，因此牛顿无能为力。您需要通过进行以下更改来使其成为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])))