I am trying to solve the problems in MDO Test suite using OpenMDAO as a means of validation. The issue I am facing is selecting the optimizer/algorithms for a given problem.
Following is the code for the Propane Combustion Problem from the suite formulated through IDF method. I tried solving it using SLSQP & NSGA2 available in pyOptSparseDriver. But both optimizer are giving sub-optimum solution.
from __future__ import print_function, division
import numpy as np
import time
from openmdao.api import *
class PropaneDis1(Component):
def __init__(self):
super(PropaneDis1, self).__init__()
self.add_param('z1', 1.0)
self.add_param('y3', np.ones(3))
self.add_output('y1', val=np.ones(2))
def solve_nonlinear(self, params, unknowns, resids):
z1 = params['z1']
y31 = params['y3'][0]
y12 = 3.0 - z1
y11 = z1*y31/y12
unknowns['y1'] = np.array([y11, y12])
class PropaneDis2(Component):
def __init__(self):
super(PropaneDis2, self).__init__()
self.add_param('z1', val= 0.0)
self.add_param('y1', val= np.ones(2))
self.add_param('y3', val= np.ones(3))
self.add_output('y2', val=np.ones(2))
def solve_nonlinear(self, params, unknowns, resids):
z1 = params['z1']
y12 = params['y1'][1]
y33 = params['y3'][2]
y21 = 0.1*z1*y33/(40.0*y12)
y22 = 0.1*(z1**2)*y33/(40.0*(y12**2))
unknowns['y2'] = np.array([y21, y22])
class PropaneDis3(Component):
def __init__(self):
super(PropaneDis3, self).__init__()
self.add_param('z1', val=0.0)
self.add_param('y1', val=np.ones(2))
self.add_param('y2', val=np.ones(2))
self.add_param('x3', val=np.ones(3))
self.add_output('y3', val=np.ones(3))
def solve_nonlinear(self, params, unknowns, resids):
z1 = params['z1']
y11 = params['y1'][0]
x31 = params['x3'][0]
y12 = params['y1'][1]
x32 = params['x3'][1]
x33 = params['x3'][2]
y21 = params['y2'][0]
y22 = params['y2'][1]
y31 = (8.0 - 2*y11 - x32 - x33)/2.0
y32 = 4*10.0 - 2*x31
y33 = z1 + y11 + x31 + y12 + x32 + x33 + y21 + y22 + y31 + y32
unknowns['y3'] = np.array([y31, y32, y33])
class PropaneIDF(Group):
def __init__(self):
super(PropaneIDF, self).__init__()
self.add('pz1', IndepVarComp('z1', val=1.0), promotes=['*'])
self.add('px3', IndepVarComp('x3', val=np.array([20.0,1.5,1.0])), promotes=['*'])
self.add('py1t', IndepVarComp('y1t', val=np.array([1.0,1.0])), promotes=['*'])
self.add('py2t', IndepVarComp('y2t', val=np.array([1.0,1.0])), promotes=['*'])
self.add('py3t', IndepVarComp('y3t', val=np.array([1.0,1.0,1.0])), promotes=['*'])
self.add('d1', PropaneDis1(), promotes=['z1', 'y1'])
self.add('d2', PropaneDis2(), promotes=['z1', 'y2'])
self.add('d3', PropaneDis3(), promotes=['z1', 'x3', 'y3'])
self.connect('y1t', 'd2.y1')
self.connect('y1t', 'd3.y1')
self.connect('y2t', 'd3.y2')
self.connect('y3t', 'd1.y3')
self.connect('y3t', 'd2.y3')
self.add('con_cmp1', ExecComp('f2=2*z1 + y1[0] + y1[1] + x3[2] + y2[0] + y3[1] + 2*y2[1] - 10.0'
, z1=0.0, y1=np.zeros(2), y2=np.zeros(2), y3=np.zeros(3), x3=np.zeros(3))
, promotes=['*'])
self.add('con_cmp2', ExecComp('f6=abs(y1[0]*y1[1])**0.5 - x3[1]*abs(40.0*z1/y3[2])**0.5'
, z1=0.0, y1=np.zeros(2), y3=np.ones(3), x3=np.zeros(3))
, promotes=['*'])
self.add('con_cmp3', ExecComp('f7=abs(z1*y1[0])**0.5 - x3[2]*abs(40.0*y1[1]/y3[2])**0.5'
, z1=0.0, y1=np.zeros(2), y3=np.ones(3), x3=np.zeros(3))
, promotes=['*'])
self.add('con_cmp4', ExecComp('f9=z1*(x3[0]**0.5) - y1[1]*y3[1]*(40.0/y3[2])**0.5'
, z1=0.0, y1=np.zeros(2), y3=np.zeros(3), x3=np.zeros(3))
, promotes=['*'])
self.add('con_cmp51', ExecComp('con51=y1t[0] - y1[0]', y1t=np.zeros(2), y1=np.zeros(2)), promotes=['*'])
self.add('con_cmp52', ExecComp('con52=y1t[1] - y1[1]', y1t=np.zeros(2), y1=np.zeros(2)), promotes=['*'])
self.add('con_cmp61', ExecComp('con61=y2t[0] - y2[0]', y2t=np.zeros(2), y2=np.zeros(2)), promotes=['*'])
self.add('con_cmp62', ExecComp('con62=y2t[1] - y2[1]', y2t=np.zeros(2), y2=np.zeros(2)), promotes=['*'])
self.add('con_cmp71', ExecComp('con71=y3t[0] - y3[0]', y3t=np.zeros(3), y3=np.zeros(3)), promotes=['*'])
self.add('con_cmp72', ExecComp('con72=y3t[1] - y3[1]', y3t=np.zeros(3), y3=np.zeros(3)), promotes=['*'])
self.add('con_cmp73', ExecComp('con73=y3t[2] - y3[2]', y3t=np.zeros(3), y3=np.zeros(3)), promotes=['*'])
self.add('obj_cmp', ExecComp('obj=2*z1 + y1[0] + y1[1] + x3[2] + y2[0] + y3[1] + 2*y2[1] - 10.0 +\
abs(y1[0]*y1[1])**0.5 - x3[1]*abs(40.0*z1/y3[2])**0.5+\
abs(z1*y1[0])**0.5 - x3[2]*abs(40.0*y1[1]/y3[2])**0.5+\
z1*(x3[0]**0.5) - y1[1]*y3[1]*(40.0/y3[2])**0.5 '
, z1=0.0, y1=np.zeros(2), y2=np.zeros(2), y3=np.zeros(3), x3=np.zeros(3))
, promotes=['*'])
self.fd_options['force_fd']=True
self.fd_options['step_size']=1.0e-12
top = Problem()
top.root = PropaneIDF()
top.driver = pyOptSparseDriver()
top.driver.options['optimizer']='SLSQP'
top.driver.add_desvar('z1', lower=0.0, upper=20.0)
top.driver.add_desvar('x3', lower=np.array([0.0, 0.0, 0.0]), upper=np.array([20.0, 20.0, 20.0]))
top.driver.add_desvar('y1t', lower=np.array([0.0, 0.0]), upper=np.array([20.0, 20.0]))
top.driver.add_desvar('y2t', lower=np.array([0.0, 0.0]), upper=np.array([20.0, 20.0]))
top.driver.add_desvar('y3t', lower=np.array([0.0, 0.0, 0.0]), upper=np.array([20.0, 20.0, 20.0]))
top.driver.add_objective('obj')
top.driver.add_constraint('con51', equals=0.0)
top.driver.add_constraint('con52', equals=0.0)
top.driver.add_constraint('con61', equals=0.0)
top.driver.add_constraint('con62', equals=0.0)
top.driver.add_constraint('con71', equals=0.0)
top.driver.add_constraint('con72', equals=0.0)
top.driver.add_constraint('con73', equals=0.0)
top.driver.add_constraint('f2', lower=0.0)
top.driver.add_constraint('f6', lower=0.0)
top.driver.add_constraint('f7', lower=0.0)
top.driver.add_constraint('f9', lower=0.0)
# top.driver.add_constraint('f2', equals=0.0)
# top.driver.add_constraint('f6', equals=0.0)
# top.driver.add_constraint('f7', equals=0.0)
# top.driver.add_constraint('f9', equals=0.0)
top.setup()
tt = time.time()
top.run()
Solution for SLSQP -
Optimization using pyOpt_sparse
================================================================================
Objective Function: _objfunc
Solution:
--------------------------------------------------------------------------------
Total Time: 0.1535
User Objective Time : 0.0183
User Sensitivity Time : 0.0806
Interface Time : 0.0534
Opt Solver Time: 0.0013
Calls to Objective Function : 10
Calls to Sens Function : 7
Objectives:
Name Value Optimum
obj 14.3997 0
Variables (c - continuous, i - integer, d - discrete):
Name Type Value Lower Bound Upper Bound
z1_0 c 2.130827 0.00e+00 2.00e+01
x3_0 c 20.000000 0.00e+00 2.00e+01
x3_1 c 0.000000 0.00e+00 2.00e+01
x3_2 c -0.000000 0.00e+00 2.00e+01
y1t_0 c 3.645641 0.00e+00 2.00e+01
y1t_1 c 0.869179 0.00e+00 2.00e+01
y2t_0 c 0.180723 0.00e+00 2.00e+01
y2t_1 c 0.679352 0.00e+00 2.00e+01
y3t_0 c 1.691004 0.00e+00 2.00e+01
y3t_1 c -0.000000 0.00e+00 2.00e+01
y3t_2 c 20.000000 0.00e+00 2.00e+01
Constraints (i - inequality, e - equality):
Name Type Bounds
con51 i 0.00e+00 <= -0.499951 <= 0.00e+00
con52 i 0.00e+00 <= 0.000006 <= 0.00e+00
con61 i 0.00e+00 <= 0.058146 <= 0.00e+00
con62 i 0.00e+00 <= 0.378850 <= 0.00e+00
con71 i 0.00e+00 <= 1.336646 <= 0.00e+00
con72 i 0.00e+00 <= -0.000000 <= 0.00e+00
con73 i 0.00e+00 <= -7.860081 <= 0.00e+00
f2 i 0.00e+00 <= -0.000000 <= 1.00e+20
f6 i 0.00e+00 <= 1.898220 <= 1.00e+20
f7 i 0.00e+00 <= 2.972127 <= 1.00e+20
f9 i 0.00e+00 <= 9.529347 <= 1.00e+20
--------------------------------------------------------------------------------
Solution for NSGA2 -
Optimization using pyOpt_sparse
================================================================================
Objective Function: _objfunc
Solution:
--------------------------------------------------------------------------------
Total Time: 107.5302
User Objective Time : 62.5869
User Sensitivity Time : 0.0000
Interface Time : 41.6987
Opt Solver Time: 3.2446
Calls to Objective Function : 99981
Calls to Sens Function : 0
Objectives:
Name Value Optimum
obj 7.67208 0
Variables (c - continuous, i - integer, d - discrete):
Name Type Value Lower Bound Upper Bound
z1_0 c 2.184051 0.00e+00 2.00e+01
x3_0 c 20.000000 0.00e+00 2.00e+01
x3_1 c 0.000000 0.00e+00 2.00e+01
x3_2 c 0.000000 0.00e+00 2.00e+01
y1t_0 c 3.757993 0.00e+00 2.00e+01
y1t_1 c 0.815949 0.00e+00 2.00e+01
y2t_0 c 0.193249 0.00e+00 2.00e+01
y2t_1 c 0.746901 0.00e+00 2.00e+01
y3t_0 c 1.592243 0.00e+00 2.00e+01
y3t_1 c -0.000000 0.00e+00 2.00e+01
y3t_2 c 20.000000 0.00e+00 2.00e+01
Constraints (i - inequality, e - equality):
Name Type Bounds
con51 i 0.00e+00 <= 0.173682 <= 0.00e+00
con52 i 0.00e+00 <= -0.520903 <= 0.00e+00
con61 i 0.00e+00 <= -0.122582 <= 0.00e+00
con62 i 0.00e+00 <= -0.292454 <= 0.00e+00
con71 i 0.00e+00 <= -0.120713 <= 0.00e+00
con72 i 0.00e+00 <= 0.000372 <= 0.00e+00
con73 i 0.00e+00 <= -7.729249 <= 0.00e+00
f2 i 0.00e+00 <= 0.143730 <= 1.00e+20
f6 i 0.00e+00 <= 1.641686 <= 1.00e+20
f7 i 0.00e+00 <= 1.957046 <= 1.00e+20
f9 i 0.00e+00 <= 3.929613 <= 1.00e+20
----------------------------------------------------------------------
I know the problem isn't with the definition as I tried running the problem with the constraints f2, f6, f7 & f9 as equality constraints (=0) instead of lower bound, I am getting the optima.
So, what other optimizers should I be looking into? Or is there any other issue that I am confusing for the optimizer's problem.
Thanks in advance.
