1ttyoff:true $ 2/* Functions and options for optimization using the algebraic 3manipulation language MACSYMA. Programmed by STOUTE (David 4Stoutemyer), Electrical Engineering Department, University of Hawaii, 54/3/74. For a description of its usage, see text file OPTMIZ USAGE. */ 6 7/* First, we set some options that respectively cause automatic 8 printing of cpu time in milliseconds, force attempted equation 9 solution even when there are more variables than unknowns, 10 enables some (time consuming) techniques for solving equations that 11 contain logs and exponentials, and enables the solution 12 of consistent singular linear equations: */ 13 14/* Next, we set a switch to suppress the message that ordinarily 15 occurs whenever a floating-point number is replaced with a rational 16 number, and to prevent 1-by-1 matrices from being converted 17 to scalars: */ 18 19/* The following pattern-matching statements permit trailing 20[] arguments to be omitted: */ 21 22/* updated by function stapoints(objective,[list_args]) 23matchdeclare([a1,a2,a3,a4], true) $ 24tellsimp (stap(a1), stapoints(a1,[],[],[])) $ 25tellsimp (stap(a1,a2), stapoints(a1,a2,[],[])) $ 26tellsimp (stap(a1,a2,a3), stapoints(a1,a2,a3,[])) $ 27tellsimp (stap(a1,a2,a3,a4), stapoints(a1,a2,a3,a4)) $ 28*/ 29eval_when([translate,batch,demo,load,loadfile],rectform(sinh(1.2-.4*%i)), 30dv(a)::=buildq([a],define_variable(a,'a,any)))$ 31 32dv(lagrangian)$ dv(eigen)$ dv(modhessian)$ 33dv(grad)$ dv(gradsub)$ dv(stapts)$ dv(mhess)$ dv(objsub)$ 34dv(mhesssub)$ dv(cpolysub)$ dv(eigen)$ dv(decslkmults)$ 35dv(ndec)$ dv(neqz)$ dv(nlez)$ dv(ntot)$ dv(eigens)$ 36 37gradient(decslkmults) := /* This function recursively defines 38 the gradient of the Lagrangian, with respect to the decision 39 variables, rtslacks, and Lagrange multipliers. */ 40 if decslkmults = [] then [] 41 else cons(diff(lagrangian, first(decslkmults)), 42 gradient(rest(decslkmults))) $ 43 44eval_when([translate,batch,demo,load,loadfile], 45modh(g,d)::=buildq([g,d], 46apply('define,[arrayapply('modhessian,[i,j]), 47 /*internal array function for MODHESSIAN*/ 48 '('(if j>i then modhessian[j,i] /*(symmetric)*/ 49 else diff(g[i],d[j]) /*minus EIGEN from up-left diag*/ 50 - (if i=j and j<=ndec+nlez then eigen else 0)))])))$ 51 52 53/* old definition 54stapoints(objective, lezeros, eqzeros, decisionvars) := block( */ 55 56/* new definition */ 57stapoints(objective,[list_args]):= 58block([lezeros, eqzeros, decisionvars], 59(if list_args=[] then lezeros: eqzeros: decisionvars:[] 60else if length(list_args)=1 then (lezeros:first(list_args), 61eqzeros: decisionvars:[]) else if length(list_args)=2 then 62(lezeros:first(list_args),eqzeros:last(list_args),decisionvars:[]) 63else if length(list_args)=3 then (lezeros:first(list_args), 64eqzeros:first(rest(list_args)), decisionvars:last(list_args)) else 65error("wrong number of args to stapoints"), block( 66/* This is the major function, which prints information about any 67stationary points, then returns the value DONE. */ 68 69 [grindswitch, solveradcan, singsolve, ratprint, scalarmatrixp, 70 dispflag,eqmult,rtslack,lemult,i,j], /* declare local variables */ 71 grindswitch: solveradcan: singsolve: true, 72 ratprint: scalarmatrixp: false, 73if member('modhessian,arrays) then apply('remarray,['modhessian]), 74 75modh(grad,decslkmults) /*end MODHESSIAN*/, 76 77 if not listp(lezeros) then lezeros: [lezeros],/* ensure list args*/ 78 if not listp(eqzeros) then eqzeros: [eqzeros], 79 if decisionvars = [] /*default to all decision variables*/ 80 then decisionvars: listofvars([objective, lezeros, eqzeros]) 81 else if not listp(decisionvars) then decisionvars: [decisionvars], 82 83 ndec: length(decisionvars), /*determine number of decision vars. */ 84 nlez: length(lezeros),/*determine number of inequality constraints*/ 85 neqz: length(eqzeros), /*determine number of equality constraints*/ 86 lagrangian: objective + sum(eqzeros[i]*eqmult[i],i,1,neqz) 87 + sum((lezeros[i]+rtslack[i]**2)*lemult[i],i,1,nlez), 88 89 decslkmults: [], /*form list of dec.vars., rtslacks & multipliers*/ 90 for i:neqz step -1 thru 1 do decslkmults: cons(eqmult[i], 91 decslkmults), 92 for i:nlez step -1 thru 1 do 93 decslkmults: cons(lemult[i], decslkmults), 94 for i:nlez step -1 thru 1 do 95 decslkmults: cons(rtslack[i], decslkmults), 96 decslkmults: append(decisionvars, decslkmults), 97 grad: gradient(decslkmults), /* form gradient */ 98 dispflag: false, /* suppress automatic output from solve */ 99 stapts: solve(grad,decslkmults),/* solve GRAD=0*/ 100 101 if stapts = [] then apply('disp,["no stationary points found"]) 102 else( ntot: ndec + nlez + nlez + neqz, 103 mhess:'mhess, mhesssub:'mhesssub, /* unbind global matrices from 104 previous case to permit different sizes. */ 105 mhess: genmatrix(modhessian, ntot, ntot), /*form HESS*/ 106 apply('remarray,['modhessian]), 107 modhessian:'modhessian, /*destroy array to permit new 108 definition for next use of analyze. */ 109 dispflag: true, /* permit automatic output from SOLVE */ 110 for i thru length(stapts) do ( 111 objsub: apply('ev,[objective,stapts[i],'rectform]), 112 /*evaluate objective.*/ 113 gradsub:apply('ev,[grad,stapts[i],'rectform]), 114 /*eval. gradient at point*/ 115 apply('ldisplay,[arrayapply('stapts,[i]), 'objsub, 'gradsub]), 116 /* output */ 117 mhesssub:apply('ev,[mhess,stapts[i],'rectform]), 118 /*eval. modified Hessian */ 119 cpolysub:rectform(newdet(mhesssub)),/*form poly in EIGEN*/ 120 /* if CPOLYSUB is univariate use REALROOTS, otherwise 121 use the slower SOLVE function: */ 122 if listofvars(cpolysub) = [eigen] and freeof('%i,cpolysub) 123 then eigens: apply('ev,[realroots(cpolysub,rootsepsilon),'numer]) 124 else eigens: solve(cpolysub, eigen) ))))) 125 /* end of function stapoints. */ $ 126 127ttyoff:false$ 128