share/calculus/optmiz.mac

   1 ttyoff:true $
   2 /* Functions and options for optimization using the algebraic
   3 manipulation language MACSYMA.  Programmed by STOUTE (David
   4 Stoutemyer), Electrical Engineering Department, University of Hawaii,
   5 4/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])
  23 matchdeclare([a1,a2,a3,a4], true) $
  24 tellsimp (stap(a1), stapoints(a1,[],[],[])) $
  25 tellsimp (stap(a1,a2), stapoints(a1,a2,[],[])) $
  26 tellsimp (stap(a1,a2,a3), stapoints(a1,a2,a3,[])) $
  27 tellsimp (stap(a1,a2,a3,a4), stapoints(a1,a2,a3,a4)) $
  28 */
  29 eval_when([translate,batch,demo,load,loadfile],rectform(sinh(1.2-.4*%i)),
  30 dv(a)::=buildq([a],define_variable(a,'a,any)))$
  31
  32 dv(lagrangian)$ dv(eigen)$ dv(modhessian)$
  33 dv(grad)$ dv(gradsub)$ dv(stapts)$ dv(mhess)$ dv(objsub)$
  34 dv(mhesssub)$ dv(cpolysub)$ dv(eigen)$ dv(decslkmults)$
  35 dv(ndec)$ dv(neqz)$ dv(nlez)$ dv(ntot)$ dv(eigens)$
  36
  37 gradient(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
  44 eval_when([translate,batch,demo,load,loadfile],
  45 modh(g,d)::=buildq([g,d],
  46 apply('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
  54 stapoints(objective, lezeros, eqzeros, decisionvars) := block( */
  55
  56 /* new definition */
  57 stapoints(objective,[list_args]):=
  58 block([lezeros, eqzeros, decisionvars],
  59 (if list_args=[] then  lezeros: eqzeros: decisionvars:[]
  60 else if length(list_args)=1 then (lezeros:first(list_args),
  61 eqzeros: decisionvars:[]) else if length(list_args)=2 then
  62 (lezeros:first(list_args),eqzeros:last(list_args),decisionvars:[])
  63 else if length(list_args)=3 then (lezeros:first(list_args),
  64 eqzeros:first(rest(list_args)), decisionvars:last(list_args)) else
  65 error("wrong number of args to stapoints"), block(
  66 /* This is the major function, which prints information about any
  67 stationary 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,
  73 if member('modhessian,arrays) then apply('remarray,['modhessian]),
  74
  75 modh(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
 127 ttyoff:false$