share/calculus/optmiz.dem

   1 load('optmiz);
   2 /* This macsyma batch file illustrates how to use the function
   3 STAP for determining the stationary points of a function of 1 or more
   4 variables, either unconstrained or subject to equality &/or inequality
   5 constraints.  As a prerequisite, see the text file  OPTMIZ USAGE.  The
   6 examples here are chosen to be instructive, geometrically visualizable,
   7 and easy enough to be solved by hand.  For a more thorough discussion
   8 and results of some more demanding test cases, see the report
   9 "Automatic analytical optimization using computer algebraic
  10 manipulation", by David R. Stoutemyer, (ALOHA project technical report,
  11 University of Hawaii, Honolulu, June 1974).
  12
  13 Here is an example with a maximum: */
  14 stapoints(5*log(y) - 3*x**2*y - 4*y) ;
  15 /* Here is an example with a saddle, which also illustrates
  16 that we may use subscripted variables: */
  17 /* this example is broken
  18 stapoints(atan(x[1]) + atanh(x[2]) + x[2]/x[1]) ; */
  19 /* Here is a famous example by Peano, for which the second-
  20 order test reveals only that the stationary point is not a maximum.
  21 It turns out that the stationary point is a saddle, unless A=B, in
  22 which case it is a minimum, along with all other points on the curve
  23 Y = C + (B*X)**2.  See "Theory of Maxima and Minima" by H. Hancock,
  24 (Dover Press) for a discussion of how to analytically determine the
  25 nature of a stationary point when the 2nd-order test is inconclusive.
  26 This example also illustrates how we may explicitly indicate the
  27 decision variables. */
  28 peano:  (y - c - (a*x)**2)*(y - c - (b*x)**2) ;
  29 stapoints(peano, [], [], [x,y]) ;
  30 /* Note how the answer is independent of A and B, but not C.
  31
  32 Now, note how when A=B, giving a non-isolated stationary point, an
  33 extra parameter is automatically introduced to describe the stationary
  34 curve: */
  35 stapoints(ev(peano,a=b), [], [], [x,y]) ;
  36 /* The limitations of STAP are primarily dependent upon the
  37 limitations of the macsyma SOLVE command for solving nonlinear
  38 equations.  SOLVE will attempt to find a closed-form solution to a
  39 single equation that is irrational or involves elementary
  40 transcendental functions; but as of June 1974, SOLVE is intended
  41 to solve simultaneous equations only if they are polynomial
  42 in the unknowns.  However, it generally converts rational equations
  43 to this form by multiplying each equation by the least common
  44 denominator of its terms; so it may solve simultaneous rational
  45 equations too.  Thus, as in our first two examples, the objective may
  46 contain any terms with rational gradients, such as ATAN(R), ATANH(R),
  47 or LOG(R), where R is rational in the decision variables.  The clearing
  48 of the denominator may result in a "solution" which is a pole for some
  49 gradient components while a zero for the others.  Although not a true
  50 solution to the equations, this is a bonus in our case, because we are
  51 interested in all extrema, not just stationary points.  In fact, we
  52 are usually more interested in any poles of the proper sign than in
  53 stationary points of the proper type with finite objective values.
  54 Poles, if found, are usually revealed in an alarming way by one or
  55 more large-magnitude components in GRADSUB or OBJSUB, or by an error
  56 interrupt such as a zerodivide.
  57
  58 A surprising number of optimization problems can be converted to the
  59 required form by a change of variable  For example, fractional powers
  60 of a decision variable, X, may be eliminated by letting X = Z**Q,
  61 where Q is the least common denominator of the fractional powers of x.
  62 For example: */
  63 frac: x**(2/3) - 2*x*y + y ;
  64 ev(frac, x=z**3);
  65 stapoints(%) ;
  66 /* A similar technique may be used to eliminate fractional powers
  67 of more complicated subexpressions, provided we include appropriate
  68 supplementary equality constraints.  For example: */
  69 sqrt(x+y)*y - x ;
  70 stapoints(ev(%,sqrt(x+y)=z), [], z**2-x-y) ;