LOAD('OPTMIZ); /* 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 ; eqcon: l + 2*w + 2*d - 72 ; stapoints(vol, l-42, eqcon) ; /* 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); volsub: ev(vol,%); ev(%, l=42-s**2); stapoints(%) ; /* 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: */ stapoints(ev(volsub,l=42)) ; stapoints(volsub) ; /* 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); decslkmults; /* 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:*/ stapoints(x, [], x**3-y**2) ; /* 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:*/ stapoints(x, [-y, y-x**3]) ; /* 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: */ stapoints(x, x**2) ; /* 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: */ peano: (y - c - (a*x)**2)*(y - c - (b*x)**2) ; limit(peano, y, inf); limit(peano, x, inf); radcan(ev(peano, y=c+(a**2+b**2)*x**2/2)); limit(%, x, inf);