python - Optimizer selection for Multi-Disciplinary Problems -
i trying solve problems in mdo test suite using openmdao means of validation. issue facing selecting optimizer/algorithms given problem.
following code propane combustion problem suite formulated through idf method. tried solving using slsqp & nsga2 available in pyoptsparsedriver. both optimizer giving sub-optimum solution.
__future__ import print_function, division import numpy np import time 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 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 objective function : 10 calls sens function : 7 objectives: name value optimum obj 14.3997 0 variables (c - continuous, - 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 0.00e+00 <= -0.499951 <= 0.00e+00 con52 0.00e+00 <= 0.000006 <= 0.00e+00 con61 0.00e+00 <= 0.058146 <= 0.00e+00 con62 0.00e+00 <= 0.378850 <= 0.00e+00 con71 0.00e+00 <= 1.336646 <= 0.00e+00 con72 0.00e+00 <= -0.000000 <= 0.00e+00 con73 0.00e+00 <= -7.860081 <= 0.00e+00 f2 0.00e+00 <= -0.000000 <= 1.00e+20 f6 0.00e+00 <= 1.898220 <= 1.00e+20 f7 0.00e+00 <= 2.972127 <= 1.00e+20 f9 0.00e+00 <= 9.529347 <= 1.00e+20 -------------------------------------------------------------------------------- solution 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 objective function : 99981 calls sens function : 0 objectives: name value optimum obj 7.67208 0 variables (c - continuous, - 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 0.00e+00 <= 0.173682 <= 0.00e+00 con52 0.00e+00 <= -0.520903 <= 0.00e+00 con61 0.00e+00 <= -0.122582 <= 0.00e+00 con62 0.00e+00 <= -0.292454 <= 0.00e+00 con71 0.00e+00 <= -0.120713 <= 0.00e+00 con72 0.00e+00 <= 0.000372 <= 0.00e+00 con73 0.00e+00 <= -7.729249 <= 0.00e+00 f2 0.00e+00 <= 0.143730 <= 1.00e+20 f6 0.00e+00 <= 1.641686 <= 1.00e+20 f7 0.00e+00 <= 1.957046 <= 1.00e+20 f9 0.00e+00 <= 3.929613 <= 1.00e+20 ---------------------------------------------------------------------- i know problem isn't definition tried running problem constraints f2, f6, f7 & f9 equality constraints (=0) instead of lower bound, getting optima.
so, other optimizers should looking into? or there other issue confusing optimizer's problem.
thanks in advance.
firstly, looking @ slsqp output, optimizer unable deliver feasible solution @ all. specifically, constraints con51 through con73 have numerous violations. these compatibility equations replace cyclic dependencies between discipline components, seems optimizer unable converge multidisciplinary coupling.
i have taken code , tried few things it, including converting on mdf configuration , switching nonlinear solver newton see how handled converging single mda cycle (so took out optimizer). assuming other design variables stayed @ initial conditions, got converge nice 1e-10 tolerance:
y1 [ 1. 2.] y2 [ 0.05886036 0.02943018] y3 [ 2. 38. 47.08829054] notice converged value y3 way outside upper bound specified when using y3t design variable under idf:
top.driver.add_desvar('y3t', lower=np.array([0.0, 0.0, 0.0]), upper=np.array([10.0, 10.0, 10.0])) i not sure got upper limits, didn't see them specified in propane problem. more importantly, haven't been able match see in code specified in mdo test suite link. suspect there may errors, wasn't able verify anything.
Comments
Post a Comment