0

I am trying to implement some distributed MDO architectures following the tips that I found here and here

When I try to implement the CO for the Sellar problem it simply does not converge, it ends up violating the constraints, according to the answer below:

Minimum target found at (5.000000, 1.964134, 0.821694)
Coupling vars target: 3.514927, 1.876374
('Minimum objective: ', 6.3073869839326475)
('constraints: ', 1.4712589277025103, 5.8788877332794067)

I think it's not an optimizer problem, since it's a simple math challenge. But I can not identify the error.

Depending on the initial values the answer is NAN when using 'FD' instead of 'CS' for gradient calculations.

My code is ( OpenMDAO 1.7):

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=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=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
        J['y2', 'z'] = np.array([[1.0, 1.0]])

        return J



class SubOptimization1(Component):
    ''' minimize differences between target and local variables of the first disipline of the sellar problem '''
    def __init__(self):
        super(SubOptimization1, self).__init__()

        # Inputs to this subprob
        self.add_param('z', val=np.array([5.0, 2.0]))
        self.add_param('x', val=1.0)
        self.add_param('y2', val=1.0)
        self.add_param('zt', val=np.array([5.0, 2.0]))
        self.add_param('xt', val=1.0)
        self.add_param('y2t', val=1.0)
        self.add_param('y1t', val=3.16)

        # Unknowns for this sub prob
        self.add_output('y1', val=1.0)

        self.problem = s1prob = Problem()
        s1prob.root = Group()
        s1prob.root.add('p1x', IndepVarComp('x', 1.0), promotes=['x'])
        s1prob.root.add('p1z', IndepVarComp('z', np.array([5.0, 2.0])), promotes=['z'])
        s1prob.root.add('p1y2', IndepVarComp('y2', 1.0), promotes=['y2'])
        s1prob.root.add('d1', SellarDis1(), promotes=['x','z','y2','y1'])

        s1prob.root.add('obj_cmp1', ExecComp('obj1 = (xt-x)**2 + (zt[0]-z[0])**2+(zt[1]-z[1])**2+(y1t-y1)**2+(y2t-y2)**2',
                                z=np.array([5.0, 2.0])), promotes= ['obj1','x','z','y1','y2','xt','zt','y1t','y2t'])

    self.deriv_options['type'] = 'cs'
        #self.fd_options['force_fd'] = True



    s1prob.driver =pyOptSparseDriver()# ScipyOptimizer() #
        s1prob.driver.options['optimizer'] = 'SLSQP'

        s1prob.driver.add_desvar('x', lower=0., upper=10.0)
        s1prob.driver.add_desvar('z', lower=np.array([-10.0, 0.0]),upper=np.array([10.0, 10.0]))
        s1prob.driver.add_desvar('y2', lower=-10.00, upper=10.00)
        s1prob.driver.add_objective('obj1')
        s1prob.driver.add_constraint('y1', lower=3.16)

        s1prob.setup()


    def solve_nonlinear(self, params, unknowns, resids):

        s1prob = self.problem

        # Pass values into our problem
        s1prob['x'] = params['x']
        s1prob['z'] = params['z']
        s1prob['y2'] = params['y2']
        s1prob['xt'] = params['xt']
        s1prob['zt'] = params['zt']
        s1prob['y1t'] = params['y1t']
        s1prob['y2t'] = params['y2t']

        # Run problem
        s1prob.run()

        # Pull values from problem
        unknowns['y1'] = s1prob['y1']



class SubOptimization2(Component):
    ''' minimize differences between target and local variables of the second disipline of the sellar problem '''
    def __init__(self):
        super(SubOptimization2, self).__init__()

        # Inputs to this subprob
        self.add_param('z', val=np.array([5.0, 2.0]))
        self.add_param('y1', val=3.16)
        self.add_param('zt', val= np.array([5.0,2.0]))
        self.add_param('y2t', val=1.0)
        self.add_param('y1t', val=3.26)

        # Unknowns for this sub prob
        self.add_output('y2', val=1.0)

        self.problem = s2prob = Problem()
        s2prob.root = Group()
        s2prob.root.add('p2z', IndepVarComp('z', np.array([5.0, 2.0])), promotes=['z'])
        s2prob.root.add('p2y1', IndepVarComp('y1', 0.0), promotes=['y1'])
        s2prob.root.add('d2', SellarDis2(), promotes=['z','y1','y2'])

        s2prob.root.add('obj_cmp2', ExecComp('obj2 =(zt[0]-z[0])**2+(zt[1]-z[1])**2+(y1t-y1)**2+(y2t-y2)**2',
                                z=np.array([5.0, 2.0])), promotes= ['obj2','z','y1','y2','zt','y1t','y2t'])
    self.deriv_options['type'] = 'cs'
        #self.fd_options['force_fd'] = True



    s2prob.driver =pyOptSparseDriver()# ScipyOptimizer() #
        s2prob.driver.options['optimizer'] = 'SLSQP'

        s2prob.driver.add_desvar('z', lower=np.array([-10.0, 0.0]),upper=np.array([10.0, 10.0]))
        s2prob.driver.add_desvar('y1', lower=-10.00, upper=10.00)
        s2prob.driver.add_objective('obj2')
        s2prob.driver.add_constraint('y2', upper=24.00)

        s2prob.setup()


    def solve_nonlinear(self, params, unknowns, resids):

        s2prob = self.problem

        # Pass values into our problem
        s2prob['z'] = params['z']
        s2prob['y1'] = params['y1']
        s2prob['zt'] = params['zt']
        s2prob['y1t'] = params['y1t']
        s2prob['y2t'] = params['y2t']

        # Run problem
        s2prob.run()

        # Pull values from problem
        unknowns['y2'] = s2prob['y2']



class SellarMDO(Group):
    ''' optimize top objective function of the sellar problem with the target variables '''
    def __init__(self):
        super(SellarMDO, self).__init__()

    #### target variables initialization ####        
    self.add('pxt', IndepVarComp('xt', 1.0), promotes= ['xt'])
        self.add('pzt', IndepVarComp('zt', np.array([5.0, 2.0])),promotes= ['zt'])
        self.add('py2t', IndepVarComp('y2t',1.0),promotes= ['y2t'])
    self.add('py1t', IndepVarComp('y1t', 3.16),promotes= ['y1t'])

        self.add('d3', SubOptimization1(), promotes= ['xt','zt','y1t','y2t'])
        self.add('d4', SubOptimization2(), promotes= ['zt','y1t','y2t'])


    #### sellar problem objective function ####
        self.add('obj_cmp', ExecComp('obj = xt**2 + zt[1] + y1t + exp(-y2t)',
                                     zt=np.array([5.0, 2.0])), promotes= ['obj','xt','zt','y1t','y2t'])


    #### First discipline constraint ####
        self.add('con1_cmp', ExecComp('con1 = (xt-x)**2 + (zt[0]-z[0])**2+(zt[1]-z[1])**2+(y1t-y1)**2+(y2t-y2)**2',
                                     z=np.array([5.0, 2.0]), x=1.0, y2=1.0,zt=np.array([5.0, 2.0])), promotes= ['con1','xt','zt','y1t','y2t'])
        self.connect("d3.x", "con1_cmp.x")
        self.connect("d3.z", "con1_cmp.z")
        self.connect("d3.y1", "con1_cmp.y1")
        self.connect("d3.y2", "con1_cmp.y2")


    #### Second discipline constraint ####
        self.add('con2_cmp', ExecComp('con2 = (zt[0]-z[0])**2+(zt[1]-z[1])**2+(y1t-y1)**2+(y2t-y2)**2',
                                     z=np.array([5.0, 2.0]), y1=3.16 ,zt=np.array([5.0, 2.0])), promotes= ['con2','zt','y1t','y2t'])
        self.connect("d4.z", "con2_cmp.z")
        self.connect("d4.y1", "con2_cmp.y1")
        self.connect("d4.y2", "con2_cmp.y2")



    self.deriv_options['type'] = 'cs'
        #self.fd_options['force_fd'] = True


if __name__ == '__main__':

    from openmdao.api import Problem, ScipyOptimizer, SqliteRecorder

    top = Problem()
    top.root = SellarMDO()

    top.driver = pyOptSparseDriver()#ScipyOptimizer() #
    top.driver.options['optimizer'] = 'SLSQP'#'NSGA2'




    top.driver.add_desvar('zt', lower=np.array([-10.0, 0.0]),upper=np.array([10.0, 10.0]))
    top.driver.add_desvar('xt', lower=0.0, upper=10.0)
    top.driver.add_desvar('y1t', lower=-10.0, upper=10.0)
    top.driver.add_desvar('y2t', lower=-10.0, upper=10.0)



    top.driver.add_objective('obj')
    top.driver.add_constraint('con1', upper=0.005)
    top.driver.add_constraint('con2', upper=0.005)



    top.setup()
    top.run()
    print("\n")
    print( "Minimum target found at (%f, %f, %f)" % (top['zt'][0],
                                             top['zt'][1],
                                             top['xt']))

    print("Coupling vars target: %f, %f" % (top['y1t'], top['y2t']))
    print("Minimum objective: ", top['obj'])
    print("constraints: ", top['con1'] , top['con2'])

Thanks very much!

aL_eX
  • 1,453
  • 2
  • 15
  • 30

1 Answers1

0

I got a CO implementation to work for the sellar problem, but I did it in OpenMDAO 2.0.

The apis are a little different, but you should be able to back convert to 1.0 if you want to. Or make the jump to 2.0!

OpenMDAO - CO (Collaborative Optimization) on Sellar test case

Justin Gray
  • 5,605
  • 1
  • 11
  • 16
  • thank you so much! I'm a bit resiliant in jumping to OpenMdao2.0 because I'm using openaerostruct too, which was built for version 1.7. Could you please tell me why it was necessary to use apply_nonlinear with an embedded optimizer? Thank you very much, – Luiz Tibério Jan 23 '18 at 14:10
  • my original implementation used ImplicitComponent for the suboptimizations because the sub-optimziers behave a bit like nonlinear solvers in this case. If you ever wanted to try post-optimality-sensitivities here, you would need to use this approach. There is a version of OpenAeroStruct available for OpenMDAO 2.0 (https://github.com/mdolab/OpenAeroStruct/tree/v2). The code in that other post should be pretty easy to back port to OpenMDAO1.7 though. – Justin Gray Jan 23 '18 at 15:12
  • Thank you! I just translated the code for version 1.7 and everything went well. I also made a brief change to use subproblem instead of component, as you suggested in the topic in the link. As for post-optimality-sensitivities, to be honest I thought I could solve the problem with: fd_options ['force_fd'] = True Now that I've seen that this option is deprecated, I'll need to think of another strategy for that ... – Luiz Tibério Jan 25 '18 at 12:50
  • my implementation is also solving things with FD. Its very possible, I was just saying that in theory you could use post optimality sensitivities, but to do so would mean posing things as an ImplicitComponent. You'd also need access to the Lagrange multipliers though, which most optimizers don't provide. – Justin Gray Jan 26 '18 at 02:19