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