OpenMDAO добавляет ограничение a‹=X‹=b для констант a,b и массива X

Я хочу ограничить (не привязать) переменную дизайна z, чтобы каждая запись была меньше или равна 5 и больше или равна -5: -5<=Z<=5. Я пытаюсь сделать это, определяя ограничение так, чтобы оно просто возвращало значения z (self.add('con_cmp1', ExecComp('con1 = z'), promotes=['z', 'con1'])), а затем определяя ограничение так, чтобы оно имело верхний предел 5.0 и нижний предел -5.0 (top.driver.add_constraint('con1', lower=np.array([-5., -5.]), upper=np.array([5.,5.]))).

Когда я это делаю, я получаю ошибку Type <type 'numpy.ndarray'> of source 'pz.z' (z) must be the same as type <type 'float'> of target 'con_cmp1.z' (z). Что означает эта ошибка? Как я могу правильно установить это ограничение?

from __future__ import print_function
from openmdao.api import ExecComp, IndepVarComp, Group, NLGaussSeidel, \
                         ScipyGMRES, Problem, ScipyOptimizer
import numpy as np

from openmdao.api import Component


class SellarDis1(Component):
    """Component containing Discipline 1."""

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

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

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

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

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

    def solve_nonlinear(self, params, unknowns, resids):
        """Evaluates the equation
        y1 = z1**2 + z2 + x1 - 0.2*y2"""

        z1 = params['z'][0]
        z2 = params['z'][1]
        x1 = params['x']
        y2 = params['y2']

        unknowns['y1'] = z1**2 + z2 + x1 - 0.2*y2

    def linearize(self, params, unknowns, resids):
        """ Jacobian for Sellar discipline 1."""
        J = {}

        J['y1','y2'] = -0.2
        J['y1','z'] = np.array([[2*params['z'][0], 1.0]])
        J['y1','x'] = 1.0

        return J


class SellarDis2(Component):
    """Component containing Discipline 2."""

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

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

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

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

    def solve_nonlinear(self, params, unknowns, resids):
        """Evaluates the equation
        y2 = y1**(.5) + z1 + z2"""

        z1 = params['z'][0]
        z2 = params['z'][1]
        y1 = params['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
        y1 = abs(y1)

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

    def linearize(self, params, unknowns, resids):
        """ Jacobian for Sellar discipline 2."""
        J = {}

        J['y2', 'y1'] = .5*params['y1']**-.5

        #Extra set of brackets below ensure we have a 2D array instead of a 1D array
        # for the Jacobian;  Note that Jacobian is 2D (num outputs x num inputs).
        J['y2', 'z'] = np.array([[1.0, 1.0]])

        return J

class SellarDerivatives(Group):
    """ Group containing the Sellar MDA. This version uses the disciplines
    with derivatives."""

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

        self.add('px', IndepVarComp('x', 1.0), promotes=['x'])
        self.add('pz', IndepVarComp('z', np.array([5.0, 2.0])), promotes=['z'])

        self.add('d1', SellarDis1(), promotes=['z', 'x', 'y1', 'y2'])
        self.add('d2', SellarDis2(), promotes=['z', 'y1', 'y2'])

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

        self.add('con_cmp1', ExecComp('con1 = z'), promotes=['z', 'con1'])

        self.nl_solver = NLGaussSeidel()
        self.nl_solver.options['atol'] = 1.0e-12

        self.ln_solver = ScipyGMRES()

top = Problem()
top.root = SellarDerivatives()

top.driver = ScipyOptimizer()
top.driver.options['optimizer'] = 'SLSQP'
top.driver.options['tol'] = 1.0e-8

top.driver.add_desvar('z', lower=np.array([-10.0, 0.0]),
                     upper=np.array([10.0, 10.0]))
top.driver.add_desvar('x', lower=0.0, upper=10.0)

top.driver.add_objective('obj')
top.driver.add_constraint('con1', lower=np.array([-5., -5.]), upper=np.array([5.,5.]))

top.setup()

# Setting initial values for design variables
top['x'] = 1.0
top['z'] = np.array([5.0, 2.0])

top.run()

print("\n")
print( "Minimum found at (%f, %f, %f)" % (top['z'][0], \
                                         top['z'][1], \
                                         top['x']))
print("Coupling vars: %f, %f" % (top['y1'], top['y2']))
print("Minimum objective: ", top['obj'])

person kilojoules    schedule 13.02.2017    source источник


Ответы (1)


Фундаментальная проблема заключается в том, что ExecComp предполагает, что все является скалярным, если вы не скажете иначе. Таким образом, вы можете заставить код, который у вас есть, работать, если вы дадите ему некоторую информацию о размере.

self.add('con_cmp1', ExecComp('con1 = z', inits={'z':np.zeros(2), 'con1': np.zeros(2)}), promotes=['z', 'con1'])

Однако на самом деле вам вовсе не нужен ExecComp. Вы можете напрямую ограничить переменную проекта следующим образом:

top.driver.add_constraint('z', lower=np.array([-5., -5.]), upper=np.array([5.,5.]))

person Justin Gray    schedule 13.02.2017