| blob | 412d768568cd1c78a4584ba4044f9948de5f4973 |
1 load('optmiz);
2 /*
3 Now let's try the famous post-office parcel problem: maximize the
4 volume of a rectangular parallelopiped parcel, subject to the
5 constraints that the length plus the girth can't exceed 72,
6 and that the length, width, and depth cannot be negative or
7 exceed 42. With three decision variables and 7 inequality constraints,
8 a direct approach using STAP would involve 17 simultaneous nonlinear
9 equations, which is beyond its present capabilities. However, since we
10 are only interested in stationary points that are maxima, it is clear
11 that none of the non-negativity constraints will be active, and the
12 length-plus-girth constraint will be active. Knowing this, we may
13 treat the length-plus-girth constraint as an equality, thus avoiding
14 one slack variable; and we may ignore the non-negativity constraints,
15 rejecting any negative solutions, avoiding three multipliers and three
16 more slacks. Moreover, neither the width nor depth can be 42, because
17 then the girth alone would exceed 72; so we may ignore these
18 constraints to avoid two more slacks and two more multipliers. These
19 simplifications result in a simple enough system of 6 simultaneous
20 nonlinear equations to be solved by SOLVE in a reasonable amount
21 of time: */
22 vol: l*w*d ;
23 eqcon: l + 2*w + 2*d - 72 ;
24 stapoints(vol, l-42, eqcon) ;
25 /* However, we may simplify the problem further, saving
26 additional time, by making the change of variable L=42-S**2, which
27 automatically satisfies the upper bound on L; and we may solve the
28 equality constraint for either W or D, then eliminate it from the
29 objective These changes reduce the problem to two simultaneous
30 nonlinear equations: */
31 solve(eqcon,w);
32 volsub: ev(vol,%);
33 ev(%, l=42-s**2);
34 stapoints(%) ;
35 /* It is almost always worth solving the largest possible subset
36 of the equality constraints for variables that enter them linearly,
37 then using this solution to eliminate these variables from the
38 remaining constraints and from the objective. This is also worth
39 doing for variables that enter nonlinearly, provided it introduces no
40 fractional powers in the objective or remaining constraints.
41 For example, to find the stationary points of X + 3*Y - 6*Z**2,
42 subject to X**2 + Y**2 + Z**2 = 1, it is well worth eliminating Z.
44 The change of variable for eliminating the inequality constraint
45 L <= 42, is equivalent to converting the inequality to an equality by
46 introducing a squared slack variable, solving for L, then eliminating
47 L from VOLSUB. From this viewpoint, the "change of variable"
48 technique is seen to be applicable to a great variety of inequality
49 constraints, not merely upper and lower bounds as is implied in most
50 textbooks. Together with applicable equality constraints, it is
51 generally worth including in the elimination the largest possible
52 subset of inequality constraints for which the eliminated variables
53 enter linearly, up to the number of decision variables minus the number
54 of equality constraints.
56 Another technique for treating an inequality constraint is to
57 solve the problem with the constraint assumed active, and also with
58 the constraint ignored, checking any stationary points for violation
59 of the constraint in the latter case: */
60 stapoints(ev(volsub,l=42)) ;
61 stapoints(volsub) ;
62 /* The second-order test is inadequate when
63 constraints have been artificially activated. For example, the one
64 eigenvalue for the case D = 15/2 above is negative, but this
65 stationary point is actually a saddle. To see this, we must check the
66 unactivated gradient at this point to see if it points into or out of
67 the feasible region: */
68 ev(grad, d=15/2, l=42);
69 decslkmults;
70 /* GRAD is a global variable bound by STAP to the
71 symbolic expression for the gradient, and DECSLKMULTS is bound to the
72 variables that the gradient is taken with respect to, in the order
73 of their components in GRAD. Thus, -405/4 points in, making the
74 point a saddle.
76 The report referenced at the beginning of this demonstration
77 explains how to generalize this test to more than one constraint.
79 When there is more than one inequality constraint, each feasible
80 combination must be activated, unless some additional convexity
81 requirements are satisfied. Nevertheless, this combinatorial
82 activation technique is probably capable of treating
83 larger problems than either the change-of-variable or the squared-
84 slack-variable-with-multiplier techniques of treating inequality
85 constraints.
87 We have already seen how STAP finds poles of the objective or gradient
88 only by accident. The following examples are included as reminders
89 about other limitations of using calculus to find extrema:*/
90 stapoints(x, [], x**3-y**2) ;
91 /* Lagrange multipliers require the constraints to have con-
92 tinuous tangents, and this is violated for the above example at X=0,
93 Y=0, where the objective is a minimum. A direct use of elimination
94 would fail too, because it results in an undefined derivative at the
95 minimum. In such cases, each piece between discontinuities must be
96 considered a separate constraint.
98 The next example has a minimum at X=0, Y=0, where the two
99 active constraints are parallel, making them linearly dependent.
100 Such cusps or "wiskers" on the feasible region violate the so-called
101 "constraint-qualification" requirement; so the extremum is not found:*/
102 stapoints(x, [-y, y-x**3]) ;
103 /* The following example doesn't have a strictly feasible point;
104 so it violates Slater's condition; and the extremum at X=0 is not
105 found: */
106 stapoints(x, x**2) ;
107 /* Finally, it is important to remember that unbounded regions
108 may have nonstationary or asymptotically stationary points
109 at infinity, which STAP will not find. Such situations are usually
110 obvious from qualitative considerations, provided the objective and
111 constraints are not too complicated and numerous, but the macsyma LIMIT
112 function can be of help. However, it is important to remember that a
113 multivariate limit depends upon the way it is taken. This is
114 illustrated by the following example, where you will have to type
115 NONZERO; in response to two interrupts: */
116 peano: (y - c - (a*x)**2)*(y - c - (b*x)**2) ;
117 limit(peano, y, inf);
118 limit(peano, x, inf);
119 radcan(ev(peano, y=c+(a**2+b**2)*x**2/2));
120 limit(%, x, inf);