parasol.Optimize

1 # Aerojet FraMES 2 # (FraMES) Framework for Modeling Engineering Systems 3 # 4 # Written by Charlie Taylor <charles.taylor@aerojet.com> (916) 355-2773 5 # Oct,21 2005 6 from scipy.optimize.cobyla import fmin_cobyla 7 import os 8 import cast 9 10 optFSys = None # set this global variable to current object being optimized 11 _optLoopCount = 0 12 _findmin = 0 13 _constraintLoopCount = 0 14

15 -def objFunction(x):

16 global _optLoopCount 17 #global _fResult 18 19 for i in range(len(x)): 20 dvStr = optFSys.optDesVars[i] 21 dv = optFSys.desVarDict[ dvStr ] 22 #print dvStr,dv.val, # debug print 23 24 optFSys.setDesignVar( dvStr, x[i]*dv.scaleFactor) 25 26 #print "calling control Routine from optimizer" 27 optFSys.evaluate() 28 29 # f[0] is the object's figure of merit attribute value 30 # after the control routine implements the new value(s) of design vars. 31 32 try: 33 fObj = optFSys.getResultVar(optFSys.figureOfMerit) 34 #print 'fObj=',fObj # debug print 35 except: 36 fObj = 0.0 37 print "ERROR... calling optimizer" 38 39 if not _findmin: 40 fObj *= -1.0 41 42 #print "in fvals", x, _fResult, 43 print ".", 44 _optLoopCount += 1 45 return fObj

46

47 -def constraintFunc(x, iconstraint):

48 global _constraintLoopCount 49 _constraintLoopCount += 1 50 51 for i in range(len(x)): 52 dvStr = optFSys.optDesVars[i] 53 dv = optFSys.desVarDict[ dvStr ] 54 #print dvStr,dv.val, # debug print 55 56 optFSys.setDesignVar( dvStr, x[i]*dv.scaleFactor) 57 58 #print "calling control Routine from optimizer" 59 optFSys.evaluate() 60 61 row = optFSys.constraintList[iconstraint] 62 constrCond, key = row 63 64 rv = optFSys.resultVarDict[ key ] 65 66 #print '... Calling Constraint',iconstraint,key,constrCond, 67 # constraint held to greater than zero 68 if constrCond == "<": 69 #print rv.hiLimit 70 return (rv.hiLimit - rv.val) # / rv.scaleFactor 71 elif constrCond == ">": 72 #print rv.loLimit 73 return (rv.val - rv.loLimit) # / rv.scaleFactor 74 else: 75 #print 'on error' 76 return 1.0 # mark as satisfied if not ">" or "<"

77 78 79

80 -def getOptVarSummary(PS, optVarList=None):

81 summary = PS.getSummary() + \ 82 '\n====================== OPTIMIZATION DESIGN VARIABLES ==' +\ 83 '=====================\n' +\ 84 ' name value minimum maximum\n' 85 86 for key in optVarList: 87 dv = PS.desVarDict[key] 88 desc = dv.description 89 if len(dv.units) > 0: 90 desc += ' (%s)'%dv.units 91 summary += '%10s %12g %12g %12g %s\n'%\ 92 (key,dv.val,dv.minVal,dv.maxVal, desc ) 93 94 if _findmin: 95 mLabel = 'Minimize' 96 else: 97 mLabel = 'Maximize' 98 fmerit = PS.resultVarDict[PS.figureOfMerit] 99 label = fmerit.description + ' (' + fmerit.name + ')' 100 summary += '\n Figure of Merit: %s = %g %s <== %s'%(label, 101 PS.getResultVar(PS.figureOfMerit), fmerit.units, mLabel) 102 103 showCon = 0 104 for key,rv in PS.resultVarDict.items(): 105 if rv.isConstrained(): 106 if not showCon: 107 showCon = 1 108 summary += '\n\n========================== OPTIMIZATION CONSTRAINTS ===='\ 109 + '=====================\n'\ 110 + ' name value minimum maximum' 111 112 minVal = '------------' 113 if rv.hasLowConstraint() : 114 minVal = '%12g'%rv.loLimit 115 maxVal = '------------' 116 if rv.hasHighConstraint(): 117 maxVal = '%12g'%rv.hiLimit 118 119 if len(rv.units) > 0: 120 desc += ' (%s)'%rv.units 121 summary += '\n%10s %12g %12s %12s %s'%\ 122 (key,rv.val,minVal,maxVal, rv.description ) 123 124 125 return summary + '\n======================================'\ 126 + '======================================\n'

127

128 -def saveOptVarSummary(PS, optVarList=None, title=''):

129 print 'saving Optimization Variable Summary to',os.path.split(PS.summFileName)[-1] 130 131 summStr = getOptVarSummary(PS, optVarList) 132 PS.summFile.write( title + summStr + '\n' ) 133 print title, summStr 134 135 PS.optimizeHistoryL.append( [title, summStr] ) 136 137 PS.htmlFile.write('<center><table class="mytable">') 138 if title: 139 PS.htmlFile.write('<th bgcolor="#CCCCCC"> %s </th>'%title) 140 PS.htmlFile.write('<tr><td nowrap>'+ '<pre>' + summStr + '</pre>' ) 141 PS.htmlFile.write('</td></tr></table></center>') 142 143 if PS.userOptions.ppt: 144 try: 145 PS.pptDoc.addTextSlide( text=summStr.replace('\n','\r'), title=title.replace('\n','\r'), 146 textFont='Courier New', textFont='Courier New', textFontSize=14, 147 noBullets=1) 148 except: 149 print "ERROR... FAILED to put plot in PowerPoint file" 150 print traceback.print_exc() 151 152 if PS.userOptions.word: 153 tableStr = [(title,),(' ',)] 154 wordTable1 = PS.wordDoc.addTable( tableStr, Range=PS.wordDoc.selectCharacter(-2) ) 155 wordTable1.Style = PS.tblstyl 156 PS.wordDoc.setCellStyle(wordTable1,1,1, just='c',bold=True, 157 fontName='Courier New', fontSize=14, bgcolor='15') 158 PS.wordDoc.setCellStyle( wordTable1, 2, 1, 159 text=summStr ) 160 PS.wordDoc.selectCharacter(-1) 161 PS.wordDoc.addText(' ') 162 163 if PS.userOptions.excel: # save all excel output later from optimizeHistoryL 164 pass

165 166

167 -def minimize(PS, figureOfMerit="mass_lbm", desVars=None, MaxLoop=500, printLevel=1):

168 optimize(PS, figureOfMerit=figureOfMerit, desVars=desVars, 169 findmin=1, MaxLoop=MaxLoop, printLevel=printLevel)

170

171 -def maximize(PS, figureOfMerit="mass_lbm", desVars=None, MaxLoop=500, printLevel=1):

172 optimize(PS, figureOfMerit=figureOfMerit, desVars=desVars, 173 findmin=0, MaxLoop=MaxLoop, printLevel=printLevel)

174

175 -def optimize(PS, figureOfMerit="mass_lbm", desVars=None, 176 findmin=1, MaxLoop=500, printLevel=1):

177 '''to find max, set findmin = 0''' 178 179 PS.firstEvaluateIfReqd() 180 181 global optFSys, _optLoopCount, _findmin, _constraintLoopCount 182 optFSys = PS 183 _optLoopCount = 0 184 _constraintLoopCount = 0 185 _findmin = findmin # use global flag in objective function 186 187 PS.figureOfMerit = figureOfMerit 188 PS.optDesVars = [] 189 190 for dvStr in desVars: 191 dv = PS.desVarDict[dvStr] 192 if PS.hasFeasibleControlVar( dv.name ): 193 print 'WARNING... %s is a feasible design variable\n It can NOT be used as an optimization variable.'%dvStr 194 pass 195 else: 196 PS.optDesVars.append( dvStr ) 197 198 #PS.summFile.write('\nPRIOR TO OPTIMIZATION\n') 199 if _findmin: 200 mLabel = 'MINIMIZE' 201 else: 202 mLabel = 'MAXIMIZE' 203 if printLevel: 204 saveOptVarSummary(PS, PS.optDesVars,'\nPRIOR TO %s OPTIMIZATION\n'%mLabel) 205 206 207 ndv = len( PS.optDesVars ) 208 ncon = 0 209 210 211 if ndv <= 0: 212 print "ERROR... there are no legal design variables for optimization" 213 return 214 215 xlow=[] 216 xhigh=[] 217 xinit = [] 218 219 # get limits of design variables 220 i = 0 221 for dvStr in PS.optDesVars: 222 dv = PS.desVarDict[dvStr] 223 xlow.append( dv.minVal / dv.scaleFactor ) 224 xhigh.append( dv.maxVal / dv.scaleFactor ) 225 xinit.append( PS.getDesignVar( dvStr) / dv.scaleFactor ) 226 i += 1 227 228 # count any result variables that might be constraints 229 PS.constraintList[:] = [] 230 for key,rv in PS.resultVarDict.items(): 231 # only handle inequality constraints, equality constraints handled by feasibility variables 232 if rv.hasLowConstraint() : 233 PS.constraintList.append( ['>',key] ) 234 ncon += 1 235 if rv.hasHighConstraint(): 236 PS.constraintList.append( ['<',key] ) 237 ncon += 1 238 239 #xOut = mdot.findminormax(findmin,ndv,ncon,MaxLoop, xlow,xhigh,xinit,objAndConstraints) 240 241 # set up constraint calling functions 242 cons = [] 243 244 # first do bounds on x array 245 def makeLo(n): 246 return lambda x: x[n] - xlow[n]

247 def makeHi(n): 248 return lambda x: xhigh[n] - x[n] 249 for n in range( len(xlow) ): 250 conLo = makeLo(n) 251 conHi = makeHi(n) 252 cons.append( conLo ) 253 cons.append( conHi ) 254 255 # then do any constraints on result variables 256 def makeCon(n): 257 return lambda x: constraintFunc(x, n) 258 for n in range( len(PS.constraintList) ): 259 con = makeCon(n) 260 cons.append( con ) 261 262 263 # rhobeg can be 1.0 because x values are scaled. 264 265 xOut = fmin_cobyla(objFunction, xinit, cons, rhobeg=1.0, rhoend=1e-6, 266 iprint=cast.intDammit(printLevel), maxfun=MaxLoop) 267 268 if printLevel: 269 print 270 print " ==> OPTIMIZATION (Loops = %i)"%_optLoopCount, "(Max = %g)"%MaxLoop, 271 print " (%i constraint loops)\n"%_constraintLoopCount 272 273 for i in range( len(PS.optDesVars) ): 274 dvStr = PS.optDesVars[i] 275 dv = PS.desVarDict[ dvStr ] 276 277 dvVal = xOut[i] * dv.scaleFactor 278 if printLevel: 279 print 'optimum',dvStr,'=',dvVal 280 #PS.summFile.write( '\noptimum %10s = %g'%(dvStr, dvVal) ) 281 if printLevel: 282 print '\nscaled desVars',PS.optDesVars,'=',xOut[:ndv] 283 284 print 'gives',figureOfMerit,'=',PS.getResultVar(figureOfMerit),'\n' 285 286 #PS.summFile.write('\n\ngives '+ str(figureOfMerit)+\ 287 # ' = '+ str(PS.getResultVar(figureOfMerit))+'\n') 288 289 290 for i in range(ndv): 291 dvStr = PS.optDesVars[i] 292 dv = PS.desVarDict[ dvStr ] 293 dvVal = xOut[i] * dv.scaleFactor 294 PS.setDesignVar( dvStr, dvVal) 295 PS.evaluate() 296 297 #PS.summFile.write('\nAFTER OPTIMIZATION\n') 298 if printLevel: 299 saveOptVarSummary(PS, PS.optDesVars,'\nAFTER %s OPTIMIZATION\n'%mLabel) 300 301 302 if __name__ == "__main__": #self test 303 304 from parasol_main import ParametricSoln 305 306 # create system object (make sure author is correct... it's used for report) 307 PS = ParametricSoln(subtaskName="Cardboard Box", 308 author="Charlie Taylor", taskName="Sample Optimization", constraintTolerance=0.001) 309 310 density = 0.025 #: density of cardboard (lbm/cuin) 311 312 # add design variables to the system (these variables may be used to 313 # optimize the system or to create plots) 314 # design vars have: 315 # name, value, minVal, maxVal, NSteps, units, description 316 PS.addDesVars( 317 ['L',12,1,48,20,'in','Length of Cardboard'], 318 ['W',12,1,48,20,'in','Width of Cardboard'], 319 ['h',4,1,6,20,'in','Height of Box'], 320 ['thk',0.1,0.02,0.4,20,'in','Thickness of Cardboard'], 321 ) 322 323 324 # now add any Result Variables That might be plotted 325 # result variables have: 326 # name, units, description 327 PS.addResultVars( 328 ['Volume','cuin','Box Volume'], 329 ['sysMass','lbm','Box Mass'], 330 ) 331 332 # the following control routine ties together the system components 333 # with the system design variables

334 - def myControlRoutine(PS):

335 # get current values of design variables 336 L,W,h,thk = PS.getDesVars("L","W","h","thk") 337 338 Volume = max(0.,(W-2*h)) * max(0.,(L-2*h)) * h 339 340 sysMass = max(0.,(W*L - 4*h**2)) * thk * density 341 342 # "sysMass" is required 343 PS.setResultVar("sysMass", sysMass) 344 PS.setResultVar("Volume", Volume)

345 346 # need to tell system the name of the control routine 347 PS.setControlRoutine(myControlRoutine) 348 349 PS.evaluate() 350 351 352 # now optimize the system... 353 optimize(PS, figureOfMerit="Volume", desVars=['h'], findmin=0) 354 355 print getOptVarSummary(PS, PS.optDesVars ) 356 print ' ===> Answer Should Be h = 2, Volume = 128' 357 358 # now save summary of system 359 PS.saveFullSummary() 360

Source Code for Module parasol.Optimize