(C13) BATCH (OPTMIZ, DEMO, DSK, SHARE);

(C14) /* This macsyma batch file illustrates how to use the function
STAP for determining the stationary points of a function of 1 or more
variables, either unconstrained or subject to equality &/or inequality
constraints.  As a prerequisite, see the text file  OPTMIZ USAGE.  The
examples here are chosen to be instructive, geometrically visualizable,
and easy enough to be solved by hand.  For a more thorough discussion
and results of some more demanding test cases, see the report 
"Automatic analytical optimization using computer algebraic
manipulation", by David R. Stoutemyer, (ALOHA project technical report,
University of Hawaii, Honolulu, June 1974).

Here is an example with a maximum: */

STAP(5*LOG(Y) - 3*X**2*Y - 4*Y) $

ALGSYS FASL DSK MACSYM BEING LOADED 
LOADING DONE

					5
(E14)                    STAPTS  = [Y = -, X = 0.0]
			       1        4

					 5
(E15)                     OBJSUB = 5 LOG(-) - 5.0
					 4

(E16)                       GRADSUB = [0.0, 0.0]

POLYRZ FASL DSK MACSYM BEING LOADED 
LOADING DONE

(E17)                          EIGEN = - 7.5

(E18)                          EIGEN = - 3.2
TIME= 4469 MSEC.

(C19) /* Here is an example with a saddle, which also illustrates 
that we may use subscripted variables: */

STAP(ATAN(X[1]) + ATANH(X[2]) + X[2]/X[1]) $

(E19)         STAPTS  = [X  = 0.409324944, X  = - 0.83245319]
		    1     2                 1

(E20)                      OBJSUB = - 0.75112805

(E21)              GRADSUB = [5.2154064E-8, 1.3411045E-7]

(E22)                       EIGEN = - 1.5897126

(E23)                        EIGEN = 1.93282488
TIME= 7350 MSEC.

(C24) /* Here is a famous example by Peano, for which the second-
order test reveals only that the stationary point is not a maximum.
It turns out that the stationary point is a saddle, unless A=B, in
which case it is a minimum, along with all other points on the curve
Y = C + (B*X)**2.  See "Theory of Maxima and Minima" by H. Hancock,
(Dover Press) for a discussion of how to analytically determine the
nature of a stationary point when the 2nd-order test is inconclusive. 
This example also illustrates how we may explicitly indicate the
decision variables. */

PEANO:  (Y - C - (A*X)**2)*(Y - C - (B*X)**2) $
TIME= 125 MSEC.

(C25) STAP(PEANO, [], [], [X,Y]) $

(E25)                    STAPTS  = [X = 0.0, Y = C]
			       1

(E26)                           OBJSUB = 0.0

(E27)                       GRADSUB = [0.0, 0.0]

(E28)                            EIGEN = 0

(E29)                            EIGEN = 2
TIME= 6177 MSEC.

(C30) /* Note how the answer is independent of A and B, but not C.

Now, note how when A=B, giving a non-isolated stationary point, an
extra parameter is automatically introduced to describe the stationary
curve: */

STAP(EV(PEANO,A=B), [], [], [X,Y]) $

				       2   2
				  - 2 B  R1  - 2 C
(E30)           STAPTS  = [Y = - -----------------, X = R1]
		      1                  2

				  2   2
			     - 2 B  R1  - 2 C    2   2     2
(E31)          OBJSUB = ( - ----------------- - B  R1  - C)
				    2

				     2   2
		      2         - 2 B  R1  - 2 C    2   2
(E32) GRADSUB = [- 4 B  R1 ( - ----------------- - B  R1  - C),
				       2

					     2   2
					- 2 B  R1  - 2 C    2   2
				 2 ( - ----------------- - B  R1  - C)]
					       2
SOLUTION

(E33)                            EIGEN = 0

				       4   2
(E34)                       EIGEN = 8 B  R1  + 2
TIME= 5647 MSEC.

(C35) /* Let's make certain that the GRADSUB expressions simplify to
zero: */

RADCAN(GRADSUB);
TIME= 323 MSEC.
(D35)                              [0, 0]

(C36) /* The limitations of STAP are primarily dependent upon the
limitations of the macsyma SOLVE command for solving nonlinear
equations.  SOLVE will attempt to find a closed-form solution to a
single equation that is irrational or involves elementary
transcendental functions; but as of June 1974, SOLVE is intended
to solve simultaneous equations only if they are polynomial
in the unknowns.  However, it generally converts rational equations
to this form by multiplying each equation by the least common
denominator of its terms; so it may solve simultaneous rational
equations too.  Thus, as in our first two examples, the objective may
contain any terms with rational gradients, such as ATAN(R), ATANH(R),
or LOG(R), where R is rational in the decision variables.  The clearing
of the denominator may result in a "solution" which is a pole for some
gradient components while a zero for the others.  Although not a true
solution to the equations, this is a bonus in our case, because we are 
interested in all extrema, not just stationary points.  In fact, we
are usually more interested in any poles of the proper sign than in
stationary points of the proper type with finite objective values.
Poles, if found, are usually revealed in an alarming way by one or
more large-magnitude components in GRADSUB or OBJSUB, or by an error
interrupt such as a zerodivide.

A surprising number of optimization problems can be converted to the
required form by a change of variable  For example, fractional powers
of a decision variable, X, may be eliminated by letting X = Z**Q,
where Q is the least common denominator of the fractional powers of x.
For example: */

FRAC: X**(2/3) - 2*X*Y + Y $
TIME= 55 MSEC.

(C37) EV(FRAC, X=Z**3);
TIME= 56 MSEC.
				     3    2
(D37)                         - 2 Y Z  + Z  + Y

(C38) STAP(%) $

(E38)            STAPTS  = [Z = 0.7937005, Y = 0.41997372]
		       1

(E39)                       OBJSUB = 0.62996052

(E40)             GRADSUB = [1.1920929E-7, - 8.9406967E-8]

(E41)                       EIGEN = - 4.9098093

(E42)                        EIGEN = 2.9098091
TIME= 6703 MSEC.

(C43) /* A similar technique may be used to eliminate fractional powers
of more complicated subexpressions, provided we include appropriate
supplementary equality constraints.  For example: */

SQRT(X+Y)*Y - X $
TIME= 31 MSEC.

(C44) STAP(EV(%,SQRT(X+Y)=Z), [], Z**2-X-Y) $

(E44)    STAPTS  = [X = 3.0, Y = - 2.0, Z = - 1.0, EQMULT  = - 1.0]
	       1                                         1

(E45)                          OBJSUB = - 1.0

(E46)                  GRADSUB = [0.0, 0.0, 0.0, 0.0]

(E47)                       EIGEN = - 1.4484026

(E48)                        EIGEN = 0.1150693
TIME= 9172 MSEC.

(C49) /* exponentials, hyperbolic functions, and trig functions may
often be converted to the required form by using the multiple-argument
formulas together with equality constraints such as the sum of the
squares of the sine and cosine equals 1.  For example: */

TRIGEXPAND: EXPTSUBST: TRUE $
TIME= 7 MSEC.

(C50) COS(2*X) + SIN(X)*EXP(2*Y) - EXP(Y);
TIME= 209 MSEC.
			    2 Y     Y      2         2
(D50)              SIN(X) %E    - %E  - SIN (X) + COS (X)

(C51) SUBST([SIN(X)=S, COS(X)=C, %E**Y=E], %);
TIME= 267 MSEC.
			       2    2          2
(D51)                       - S  + E  S - E + C

(C52) STAP(%, [], S**2+C**2-1) $

(E52) STAPTS  = [E = 1.25992104, C = - 0.91788341, S = 0.39685026,
	    1

						       EQMULT  = - 1.0]
							     1

(E53)                       OBJSUB = 0.055059291

(E54)        GRADSUB = [0, - 1.49011612E-8, 0.0, 8.19563865E-8]

(E55)                       EIGEN = - 4.4000479

(E56)                        EIGEN = 1.82370886

(E57) STAPTS  = [E = 1.25992104, C = 0.91788339, S = 0.39685026,
	    2

						       EQMULT  = - 1.0]
							     1

(E58)                       OBJSUB = 0.055059254

(E59)        GRADSUB = [0, - 1.49011612E-8, 0.0, 4.4703483E-8]

(E60)                       EIGEN = - 4.4000479

(E61)                        EIGEN = 1.82370886

				 9        1
(E62)       STAPTS  = [EQMULT  = -, E = - -, S = - 1.0, C = 0.0]
		  3          1   8        2

(E63)                         OBJSUB = - 0.75

(E64)                  GRADSUB = [0.0, 0.0, 0.0, 0.0]

(E65)                           EIGEN = - 2

(E66)                           EIGEN = 4.25

				   7      1
(E67)         STAPTS  = [EQMULT  = -, E = -, S = 1.0, C = 0.0]
		    4          1   8      2

(E68)                         OBJSUB = - 1.25

(E69)                  GRADSUB = [0.0, 0.0, 0.0, 0.0]

(E70)                            EIGEN = 2

(E71)                           EIGEN = 3.75
TIME= 30639 MSEC.

(C72) /* For any of the above answers, we may use LOG(E) or ?ATAN(S,C)
to compute the corresponding values of the original decision
variables, where E, S, or C are the right-hand-sides of the
appropriate components of the chosen component of STAPTS.

Now let's try the famous post-office parcel problem:  maximize the
volume of a rectangular parallelopiped parcel, subject to the
constraints that the length plus the girth can't exceed 72,
and that the length, width, and depth cannot be negative or 
exceed 42.  With three decision variables and 7 inequality constraints,
a direct approach using STAP would involve 17 simultaneous nonlinear
equations, which is beyond its present capabilities.  However, since we
are only interested in stationary points that are maxima, it is clear
that none of the non-negativity constraints will be active, and the
length-plus-girth constraint will be active.  Knowing this, we may
treat the length-plus-girth constraint as an equality, thus avoiding 
one slack variable; and we may ignore the non-negativity constraints,
rejecting any negative solutions, avoiding three multipliers and three
more slacks.  Moreover, neither the width nor depth can be 42, because
then the girth alone would exceed 72; so we may ignore these
constraints to avoid two more slacks and two more multipliers.  These
simplifications result in a simple enough system of 6 simultaneous
nonlinear equations to be solved by SOLVE in a reasonable amount
of time: */

VOL: L*W*D $
TIME= 14 MSEC.

(C73) EQCON: L + 2*W + 2*D - 72 $
TIME= 31 MSEC.

(C74) STAP(VOL, L-42, EQCON) $

(E74) STAPTS  = [RTSLACK  = 4.2426405, W = 12.0, D = 12.0, L = 24.0,
	    1           1

				      LEMULT  = 0.0, EQMULT  = - 144.0]
					    1              1

(E75)                         OBJSUB = 3456.0

(E76)       GRADSUB = [0.0, 0.0, 0.0, 0.0, - 1.66893005E-6, 0.0]

(E77)                           EIGEN = - 24

(E78)                        EIGEN = - 7.902439

				     15      15
(E79) STAPTS  = [RTSLACK  = 0.0, W = --, D = --, L = 42.0,
	    2           1            2       2

						  405              315
					LEMULT  = ---, EQMULT  = - ---]
					      1    4         1      2

(E80)                         OBJSUB = 2362.5

(E81)              GRADSUB = [0.0, 0, 0.0, 0.0, 0.0, 0.0]

(E82)                           EIGEN = - 42

(E83)                          EIGEN = 202.5
TIME= 47070 MSEC.

(C84) /* However, we may simplify the problem further, saving
additional time, by making the change of variable L=42-S**2, which
automatically satisfies  the upper bound on L; and we may solve the
equality constraint for either W or D, then eliminate it from the
objective These changes reduce the problem to two simultaneous
nonlinear equations: */

SOLVE(EQCON,W);
SOLUTION

				   L + 2 D - 72
(E84)                        W = - ------------
					2
TIME= 69 MSEC.
(D84)                              [E84]

(C85) VOLSUB: EV(VOL,%);
TIME= 66 MSEC.
			      D L (L + 2 D - 72)
(D85)                       - ------------------
				      2

(C86) EV(%, L=42-S**2);
TIME= 78 MSEC.
				 2       2
			D (42 - S ) ( - S  + 2 D - 30)
(D86)                 - ------------------------------
				      2

(C87) STAP(%) $

(E87)               STAPTS  = [S = 4.2426405, D = 12.0]
			  1

(E88)                         OBJSUB = 3456.0

(E89)            GRADSUB = [- 1.90734863E-5, 1.8310547E-4]

(E90)                       EIGEN = - 876.51387

(E91)                       EIGEN = - 35.486028

						15
(E92)                   STAPTS  = [S = 0.0, D = --]
			      2                 2

(E93)                         OBJSUB = 2362.5

(E94)                       GRADSUB = [0.0, 0.0]

(E95)                           EIGEN = - 84

(E96)                          EIGEN = 202.5
TIME= 18140 MSEC.

(C97) /* It is almost always worth solving the largest possible subset
of the equality constraints for variables that enter them linearly,
then using this solution to eliminate these variables from the
remaining constraints and from the objective.  This is also worth
doing for variables that enter nonlinearly, provided it introduces no
fractional powers in the objective or remaining constraints.
For example, to find the stationary points of  X + 3*Y - 6*Z**2,
subject to  X**2 + Y**2 + Z**2 = 1, it is well worth eliminating Z.

The change of variable for eliminating the inequality constraint
L <= 42, is equivalent to converting the inequality to an equality by
introducing a squared slack variable, solving for L, then eliminating
L from VOLSUB.  From this viewpoint, the "change of variable"
technique is seen to be applicable to a great variety of inequality
constraints, not merely upper and lower bounds as is implied in most
textbooks.  Together with applicable equality constraints, it is
generally worth including in the elimination the largest possible
subset of inequality constraints for which the eliminated variables
enter linearly, up to the number of decision variables minus the number
of equality constraints.

Another technique for treating an inequality constraint is to
solve the problem with the constraint assumed active, and also with
the constraint ignored, checking any stationary points for violation
of the constraint in the latter case: */

STAP(EV(VOLSUB,L=42)) $

				       15
(E98)                              D = --
				       2

					4725
(E99)                          OBJSUB = ----
					 2

(E100)                          GRADSUB = [0]

(E101)                          EIGEN = - 84
TIME= 1071 MSEC.

(C102) STAP(VOLSUB) $

(E102)                 STAPTS  = [L = 24.0, D = 12.0]
			     1

(E103)                         OBJSUB = 3456.0

(E104)                      GRADSUB = [0.0, 0.0]

(E105)                       EIGEN = - 51.633308

(E106)                      EIGEN = - 8.36669242

(E107)                  STAPTS  = [D = 36.0, L = 0.0]
			      2

(E108)                          OBJSUB = 0.0

(E109)                      GRADSUB = [0.0, 0.0]

(E110)                       EIGEN = - 58.249224

(E111)                       EIGEN = 22.2492234

(E112)                  STAPTS  = [L = 72.0, D = 0.0]
			      3

(E113)                          OBJSUB = 0.0

(E114)                      GRADSUB = [0.0, 0.0]

(E115)                      EIGEN = - 152.498447

(E116)                       EIGEN = 8.49844718

(E117)                  STAPTS  = [D = 0.0, L = 0.0]
			      4

(E118)                          OBJSUB = 0.0

(E119)                      GRADSUB = [0.0, 0.0]

(E120)                          EIGEN = - 36

(E121)                           EIGEN = 36
TIME= 11506 MSEC.

(C122) /* The second-order test is inadequate when
constraints have been artifically activated.  For example, the one
eigenvalue for the case  D = 15/2 above is negative, but this
stationary point is actually a saddle.  To see this, we must check the 
unactivated gradient at this point to see if it points into or out of
the feasible region: */

EV(GRAD, D=15/2, L=42);
TIME= 142 MSEC.
				       405
(D122)                           [0, - ---]
					4

(C123) DECSLKMULTS;
TIME= 1 MSEC.
(D123)                             [D, L]

(C124) /*  GRAD is a global varible bound by STAP to the
symbolic expression for the gradient, and DECSLKMULTS is bound to the
variables that the gradient is taken with respect to, in the order
of their components in GRAD.  Thus, -405/4 points in, making the
point a saddle.

The report referenced at the beginning of this demonstration
explains how to generalize this test to more than one constraint.

When there is more than one inequality constraint, each feasible
combination must be activated, unless some additional convexity
requirements are satisfied.  Nevertheless, this combinatorial
activation technique is probably capable of treating
larger problems than either the change-of-variable or the squared-
slack-variable-with-multiplier techniques of treating inequality 
constraints. 

We have already seen how STAP finds poles of the objective or gradient
only by accident.  The following examples are included as reminders
about other limitations of using calculus to find extrema:*/

STAP(X, [], X**3-Y**2) $

(E124)                   NO STATIONARY POINTS FOUND
TIME= 4727 MSEC.

(C125) /* Lagrange multipliers require the constraints to have con-
tinuous tangents, and this is violated for the above example at X=0,
Y=0, where the objective is a minimum.  A direct use of elimination
would fail too, because it results in an undefined derivative at the
minimum.  In such cases, each piece between discontinuities must be
considered a separate constraint.

The next example has a minimum at X=0, Y=0, where the two
active constraints are parallel, making them linearly dependent.
Such cusps or "wiskers" on the feasible region violate the so-called
"constraint-qualification" requirement; so the extremum is not found:*/

STAP(X, [-Y, Y-X**3]) $

(E125)                   NO STATIONARY POINTS FOUND
TIME= 17706 MSEC.

(C126) /* The following example doesn't have a strictly feasible point;
so it violates Slater's condition; and the extremum at X=0 is not
found: */

STAP(X, X**2) $

(E126)                   NO STATIONARY POINTS FOUND
TIME= 4668 MSEC.

(C127) /* Finally, it is important to remember that unbounded regions
may have nonstationary or asymptotically stationary points
at infinity, which STAP will not find.  Such situations are usually
obvious from qualitative considerations, provided the objective and
constraints are not too complicated and numerous, but the macsyma LIMIT
function can be of help.  However, it is important to remember that a
multivariate limit depends upon the way it is taken.  This is
illustrated by the following example, where you will have to type
NONZERO; in response to two interrupts: */

LIMIT(PEANO, Y, INF);

LIMIT FASL DSK MACSYM BEING LOADED 
LOADING DONE
TIME= 280 MSEC.
(D127)                               INF

(C128) LIMIT(PEANO, X, INF);
IS  A B  ZERO OR NONZERO?

NONZERO;
TIME= 402 MSEC.
(D128)                               INF

(C129) RADCAN(EV(PEANO, Y=C+(A**2+B**2)*X**2/2));
TIME= 818 MSEC.
			      4      2  2    4   4
			    (B  - 2 A  B  + A ) X
(D129)                    - ----------------------
				      4

(C130) LIMIT(%, X, INF);
IS  (B + A) (B - A)  ZERO OR NONZERO?

NONZERO;
TIME= 1085 MSEC.
(D130)                              MINF

TIME= 223257 MSEC.
(D131)                           BATCH DONE