extcode_paraboloid.py

if __name__ == '__main__':
    import sys
    import numpy as np
    input_filename = sys.argv[1]
    output_filename = sys.argv[2]
  
    file_contents = np.loadtxt(input_filename)
    x = file_contents [0] 
    y = file_contents [1] 
   
    f_xy = (x-3.0)**2 + x*y + (y+4.0)**2 - 3.0
    g_xy = (x-3.0)**2 + x*y + (y+4.0)**2 - 3.0 
    g_xy = g_xy *3

    np.savetxt(output_filename,[f_xy,g_xy])    






    

ObjComp.py

from __future__ import division
from openmdao.api import ExplicitComponent
import numpy as np

class ObjComp(ExplicitComponent):

    def setup(self):
        
        self.add_input('f_xy')
        self.add_input('g_xy')
     
        self.add_output('obj')        

        self.declare_partials(of='obj', wrt='f_xy')
        self.declare_partials(of='obj', wrt='g_xy')


    def compute(self, inputs, outputs):
        outputs['obj'] = 3*126*inputs['f_xy'] + 3*500*inputs['g_xy']
        print('OBJECTIVE',outputs['obj'][0])
        print('----------------------------------------------------------------------')        
        
    def compute_partials(self, inputs, partials):
        partials['obj', 'f_xy'] = 3*126
        partials['obj', 'g_xy'] = 3*500

        
        

Optimizer.py

from openmdao.api import Problem, Group, IndepVarComp
from openmdao.api import  ScipyOptimizeDriver
from ParaboloidExternalCode import ParaboloidExternalCode
from ObjComp import ObjComp

top = Problem()
top.model = model = Group()

# create and connect inputs
model.add_subsystem('p1', IndepVarComp('x', 3.0))
model.add_subsystem('p2', IndepVarComp('y', -4.0))
model.add_subsystem('p', ParaboloidExternalCode())
model.add_subsystem('ObjComp', ObjComp())

model.connect('p1.x', 'p.x')
model.connect('p2.y', 'p.y')
model.connect('p.f_xy', 'ObjComp.f_xy')
model.connect('p.g_xy', 'ObjComp.g_xy')

top.driver = ScipyOptimizeDriver()
top.driver.options['optimizer'] = 'SLSQP'

top.model.add_design_var('p1.x', lower=-50, upper=50)
top.model.add_design_var('p2.y', lower=-50, upper=50)
top.model.add_objective('ObjComp.obj')
#top.model.add_constraint('p.f_xy', lower=1)

top.driver.options['tol'] = 1e-5
top.driver.options['disp'] = True
top.setup()
#top.setup(check=True, mode='fwd')

top.run_driver()
# minimum value
# location of the minimum
print(top['p1.x'])
print(top['p2.y'])
print(top['ObjComp.obj'])

ParaboloidExternalCode.py

from openmdao.api import ExternalCode
import numpy as np
class ParaboloidExternalCode(ExternalCode):
    def setup(self):
        self.add_input('x', val=0.0)
        self.add_input('y', val=0.0)

        self.add_output('f_xy', val=0.0)
        self.add_output('g_xy', val=0.0)

        self.input_file = 'paraboloid_input.dat'
        self.output_file = 'paraboloid_output.dat'

        # providing these is optional; the component will verify that any input
        # files exist before execution and that the output files exist after.
        self.options['external_input_files'] = [self.input_file,]
        self.options['external_output_files'] = [self.output_file,]

        self.options['command'] = [
            'python', 'extcode_paraboloid.py', self.input_file, self.output_file
        ]
        self.declare_partials(of='*', wrt='*', method='fd', form='central', step=1e-6)
        

    def compute(self, inputs, outputs):
        x = inputs['x']
        y = inputs['y']
        
        # generate the input file for the paraboloid external code
        with open(self.input_file, 'w') as input_file:
            input_file.write('%f\n%f\n' % (x,y))

        # the parent compute function actually runs the external code
        super(ParaboloidExternalCode, self).compute(inputs, outputs)

        file_contents = np.loadtxt(self.output_file)

  
        outputs['f_xy'] = file_contents [0] 
        outputs['g_xy'] = file_contents [1] 

        print('ExtInX = {}, ExtInY = {},, ExtOutF_XY = {}, ExtOutG_XY = {} '.format(x,y,file_contents [0],file_contents[1]) )