share/calculus/optmiz.transcript

   1 (C13) BATCH (OPTMIZ, DEMO, DSK, SHARE);
   2
   3 (C14) /* This macsyma batch file illustrates how to use the function
   4 STAP for determining the stationary points of a function of 1 or more
   5 variables, either unconstrained or subject to equality &/or inequality
   6 constraints.  As a prerequisite, see the text file  OPTMIZ USAGE.  The
   7 examples here are chosen to be instructive, geometrically visualizable,
   8 and easy enough to be solved by hand.  For a more thorough discussion
   9 and results of some more demanding test cases, see the report
  10 "Automatic analytical optimization using computer algebraic
  11 manipulation", by David R. Stoutemyer, (ALOHA project technical report,
  12 University of Hawaii, Honolulu, June 1974).
  13
  14 Here is an example with a maximum: */
  15
  16 STAP(5*LOG(Y) - 3*X**2*Y - 4*Y) $
  17
  18 ALGSYS FASL DSK MACSYM BEING LOADED
  19 LOADING DONE
  20
  21                                         5
  22 (E14)                    STAPTS  = [Y = -, X = 0.0]
  23                                1        4
  24
  25                                          5
  26 (E15)                     OBJSUB = 5 LOG(-) - 5.0
  27                                          4
  28
  29 (E16)                       GRADSUB = [0.0, 0.0]
  30
  31 POLYRZ FASL DSK MACSYM BEING LOADED
  32 LOADING DONE
  33
  34 (E17)                          EIGEN = - 7.5
  35
  36 (E18)                          EIGEN = - 3.2
  37 TIME= 4469 MSEC.
  38
  39 (C19) /* Here is an example with a saddle, which also illustrates
  40 that we may use subscripted variables: */
  41
  42 STAP(ATAN(X[1]) + ATANH(X[2]) + X[2]/X[1]) $
  43
  44 (E19)         STAPTS  = [X  = 0.409324944, X  = - 0.83245319]
  45                     1     2                 1
  46
  47 (E20)                      OBJSUB = - 0.75112805
  48
  49 (E21)              GRADSUB = [5.2154064E-8, 1.3411045E-7]
  50
  51 (E22)                       EIGEN = - 1.5897126
  52
  53 (E23)                        EIGEN = 1.93282488
  54 TIME= 7350 MSEC.
  55
  56 (C24) /* Here is a famous example by Peano, for which the second-
  57 order test reveals only that the stationary point is not a maximum.
  58 It turns out that the stationary point is a saddle, unless A=B, in
  59 which case it is a minimum, along with all other points on the curve
  60 Y = C + (B*X)**2.  See "Theory of Maxima and Minima" by H. Hancock,
  61 (Dover Press) for a discussion of how to analytically determine the
  62 nature of a stationary point when the 2nd-order test is inconclusive.
  63 This example also illustrates how we may explicitly indicate the
  64 decision variables. */
  65
  66 PEANO:  (Y - C - (A*X)**2)*(Y - C - (B*X)**2) $
  67 TIME= 125 MSEC.
  68
  69 (C25) STAP(PEANO, [], [], [X,Y]) $
  70
  71 (E25)                    STAPTS  = [X = 0.0, Y = C]
  72                                1
  73
  74 (E26)                           OBJSUB = 0.0
  75
  76 (E27)                       GRADSUB = [0.0, 0.0]
  77
  78 (E28)                            EIGEN = 0
  79
  80 (E29)                            EIGEN = 2
  81 TIME= 6177 MSEC.
  82
  83 (C30) /* Note how the answer is independent of A and B, but not C.
  84
  85 Now, note how when A=B, giving a non-isolated stationary point, an
  86 extra parameter is automatically introduced to describe the stationary
  87 curve: */
  88
  89 STAP(EV(PEANO,A=B), [], [], [X,Y]) $
  90
  91                                        2   2
  92                                   - 2 B  R1  - 2 C
  93 (E30)           STAPTS  = [Y = - -----------------, X = R1]
  94                       1                  2
  95
  96                                   2   2
  97                              - 2 B  R1  - 2 C    2   2     2
  98 (E31)          OBJSUB = ( - ----------------- - B  R1  - C)
  99                                     2
 100
 101                                      2   2
 102                       2         - 2 B  R1  - 2 C    2   2
 103 (E32) GRADSUB = [- 4 B  R1 ( - ----------------- - B  R1  - C),
 104                                        2
 105
 106                                              2   2
 107                                         - 2 B  R1  - 2 C    2   2
 108                                  2 ( - ----------------- - B  R1  - C)]
 109                                                2
 110 SOLUTION
 111
 112 (E33)                            EIGEN = 0
 113
 114                                        4   2
 115 (E34)                       EIGEN = 8 B  R1  + 2
 116 TIME= 5647 MSEC.
 117
 118 (C35) /* Let's make certain that the GRADSUB expressions simplify to
 119 zero: */
 120
 121 RADCAN(GRADSUB);
 122 TIME= 323 MSEC.
 123 (D35)                              [0, 0]
 124
 125 (C36) /* The limitations of STAP are primarily dependent upon the
 126 limitations of the macsyma SOLVE command for solving nonlinear
 127 equations.  SOLVE will attempt to find a closed-form solution to a
 128 single equation that is irrational or involves elementary
 129 transcendental functions; but as of June 1974, SOLVE is intended
 130 to solve simultaneous equations only if they are polynomial
 131 in the unknowns.  However, it generally converts rational equations
 132 to this form by multiplying each equation by the least common
 133 denominator of its terms; so it may solve simultaneous rational
 134 equations too.  Thus, as in our first two examples, the objective may
 135 contain any terms with rational gradients, such as ATAN(R), ATANH(R),
 136 or LOG(R), where R is rational in the decision variables.  The clearing
 137 of the denominator may result in a "solution" which is a pole for some
 138 gradient components while a zero for the others.  Although not a true
 139 solution to the equations, this is a bonus in our case, because we are
 140 interested in all extrema, not just stationary points.  In fact, we
 141 are usually more interested in any poles of the proper sign than in
 142 stationary points of the proper type with finite objective values.
 143 Poles, if found, are usually revealed in an alarming way by one or
 144 more large-magnitude components in GRADSUB or OBJSUB, or by an error
 145 interrupt such as a zerodivide.
 146
 147 A surprising number of optimization problems can be converted to the
 148 required form by a change of variable  For example, fractional powers
 149 of a decision variable, X, may be eliminated by letting X = Z**Q,
 150 where Q is the least common denominator of the fractional powers of x.
 151 For example: */
 152
 153 FRAC: X**(2/3) - 2*X*Y + Y $
 154 TIME= 55 MSEC.
 155
 156 (C37) EV(FRAC, X=Z**3);
 157 TIME= 56 MSEC.
 158                                      3    2
 159 (D37)                         - 2 Y Z  + Z  + Y
 160
 161 (C38) STAP(%) $
 162
 163 (E38)            STAPTS  = [Z = 0.7937005, Y = 0.41997372]
 164                        1
 165
 166 (E39)                       OBJSUB = 0.62996052
 167
 168 (E40)             GRADSUB = [1.1920929E-7, - 8.9406967E-8]
 169
 170 (E41)                       EIGEN = - 4.9098093
 171
 172 (E42)                        EIGEN = 2.9098091
 173 TIME= 6703 MSEC.
 174
 175 (C43) /* A similar technique may be used to eliminate fractional powers
 176 of more complicated subexpressions, provided we include appropriate
 177 supplementary equality constraints.  For example: */
 178
 179 SQRT(X+Y)*Y - X $
 180 TIME= 31 MSEC.
 181
 182 (C44) STAP(EV(%,SQRT(X+Y)=Z), [], Z**2-X-Y) $
 183
 184 (E44)    STAPTS  = [X = 3.0, Y = - 2.0, Z = - 1.0, EQMULT  = - 1.0]
 185                1                                         1
 186
 187 (E45)                          OBJSUB = - 1.0
 188
 189 (E46)                  GRADSUB = [0.0, 0.0, 0.0, 0.0]
 190
 191 (E47)                       EIGEN = - 1.4484026
 192
 193 (E48)                        EIGEN = 0.1150693
 194 TIME= 9172 MSEC.
 195
 196 (C49) /* exponentials, hyperbolic functions, and trig functions may
 197 often be converted to the required form by using the multiple-argument
 198 formulas together with equality constraints such as the sum of the
 199 squares of the sine and cosine equals 1.  For example: */
 200
 201 TRIGEXPAND: EXPTSUBST: TRUE $
 202 TIME= 7 MSEC.
 203
 204 (C50) COS(2*X) + SIN(X)*EXP(2*Y) - EXP(Y);
 205 TIME= 209 MSEC.
 206                             2 Y     Y      2         2
 207 (D50)              SIN(X) %E    - %E  - SIN (X) + COS (X)
 208
 209 (C51) SUBST([SIN(X)=S, COS(X)=C, %E**Y=E], %);
 210 TIME= 267 MSEC.
 211                                2    2          2
 212 (D51)                       - S  + E  S - E + C
 213
 214 (C52) STAP(%, [], S**2+C**2-1) $
 215
 216 (E52) STAPTS  = [E = 1.25992104, C = - 0.91788341, S = 0.39685026,
 217             1
 218
 219                                                        EQMULT  = - 1.0]
 220                                                              1
 221
 222 (E53)                       OBJSUB = 0.055059291
 223
 224 (E54)        GRADSUB = [0, - 1.49011612E-8, 0.0, 8.19563865E-8]
 225
 226 (E55)                       EIGEN = - 4.4000479
 227
 228 (E56)                        EIGEN = 1.82370886
 229
 230 (E57) STAPTS  = [E = 1.25992104, C = 0.91788339, S = 0.39685026,
 231             2
 232
 233                                                        EQMULT  = - 1.0]
 234                                                              1
 235
 236 (E58)                       OBJSUB = 0.055059254
 237
 238 (E59)        GRADSUB = [0, - 1.49011612E-8, 0.0, 4.4703483E-8]
 239
 240 (E60)                       EIGEN = - 4.4000479
 241
 242 (E61)                        EIGEN = 1.82370886
 243
 244                                  9        1
 245 (E62)       STAPTS  = [EQMULT  = -, E = - -, S = - 1.0, C = 0.0]
 246                   3          1   8        2
 247
 248 (E63)                         OBJSUB = - 0.75
 249
 250 (E64)                  GRADSUB = [0.0, 0.0, 0.0, 0.0]
 251
 252 (E65)                           EIGEN = - 2
 253
 254 (E66)                           EIGEN = 4.25
 255
 256                                    7      1
 257 (E67)         STAPTS  = [EQMULT  = -, E = -, S = 1.0, C = 0.0]
 258                     4          1   8      2
 259
 260 (E68)                         OBJSUB = - 1.25
 261
 262 (E69)                  GRADSUB = [0.0, 0.0, 0.0, 0.0]
 263
 264 (E70)                            EIGEN = 2
 265
 266 (E71)                           EIGEN = 3.75
 267 TIME= 30639 MSEC.
 268
 269 (C72) /* For any of the above answers, we may use LOG(E) or ?ATAN(S,C)
 270 to compute the corresponding values of the original decision
 271 variables, where E, S, or C are the right-hand-sides of the
 272 appropriate components of the chosen component of STAPTS.
 273
 274 Now let's try the famous post-office parcel problem:  maximize the
 275 volume of a rectangular parallelopiped parcel, subject to the
 276 constraints that the length plus the girth can't exceed 72,
 277 and that the length, width, and depth cannot be negative or
 278 exceed 42.  With three decision variables and 7 inequality constraints,
 279 a direct approach using STAP would involve 17 simultaneous nonlinear
 280 equations, which is beyond its present capabilities.  However, since we
 281 are only interested in stationary points that are maxima, it is clear
 282 that none of the non-negativity constraints will be active, and the
 283 length-plus-girth constraint will be active.  Knowing this, we may
 284 treat the length-plus-girth constraint as an equality, thus avoiding
 285 one slack variable; and we may ignore the non-negativity constraints,
 286 rejecting any negative solutions, avoiding three multipliers and three
 287 more slacks.  Moreover, neither the width nor depth can be 42, because
 288 then the girth alone would exceed 72; so we may ignore these
 289 constraints to avoid two more slacks and two more multipliers.  These
 290 simplifications result in a simple enough system of 6 simultaneous
 291 nonlinear equations to be solved by SOLVE in a reasonable amount
 292 of time: */
 293
 294 VOL: L*W*D $
 295 TIME= 14 MSEC.
 296
 297 (C73) EQCON: L + 2*W + 2*D - 72 $
 298 TIME= 31 MSEC.
 299
 300 (C74) STAP(VOL, L-42, EQCON) $
 301
 302 (E74) STAPTS  = [RTSLACK  = 4.2426405, W = 12.0, D = 12.0, L = 24.0,
 303             1           1
 304
 305                                       LEMULT  = 0.0, EQMULT  = - 144.0]
 306                                             1              1
 307
 308 (E75)                         OBJSUB = 3456.0
 309
 310 (E76)       GRADSUB = [0.0, 0.0, 0.0, 0.0, - 1.66893005E-6, 0.0]
 311
 312 (E77)                           EIGEN = - 24
 313
 314 (E78)                        EIGEN = - 7.902439
 315
 316                                      15      15
 317 (E79) STAPTS  = [RTSLACK  = 0.0, W = --, D = --, L = 42.0,
 318             2           1            2       2
 319
 320                                                   405              315
 321                                         LEMULT  = ---, EQMULT  = - ---]
 322                                               1    4         1      2
 323
 324 (E80)                         OBJSUB = 2362.5
 325
 326 (E81)              GRADSUB = [0.0, 0, 0.0, 0.0, 0.0, 0.0]
 327
 328 (E82)                           EIGEN = - 42
 329
 330 (E83)                          EIGEN = 202.5
 331 TIME= 47070 MSEC.
 332
 333 (C84) /* However, we may simplify the problem further, saving
 334 additional time, by making the change of variable L=42-S**2, which
 335 automatically satisfies  the upper bound on L; and we may solve the
 336 equality constraint for either W or D, then eliminate it from the
 337 objective These changes reduce the problem to two simultaneous
 338 nonlinear equations: */
 339
 340 SOLVE(EQCON,W);
 341 SOLUTION
 342
 343                                    L + 2 D - 72
 344 (E84)                        W = - ------------
 345                                         2
 346 TIME= 69 MSEC.
 347 (D84)                              [E84]
 348
 349 (C85) VOLSUB: EV(VOL,%);
 350 TIME= 66 MSEC.
 351                               D L (L + 2 D - 72)
 352 (D85)                       - ------------------
 353                                       2
 354
 355 (C86) EV(%, L=42-S**2);
 356 TIME= 78 MSEC.
 357                                  2       2
 358                         D (42 - S ) ( - S  + 2 D - 30)
 359 (D86)                 - ------------------------------
 360                                       2
 361
 362 (C87) STAP(%) $
 363
 364 (E87)               STAPTS  = [S = 4.2426405, D = 12.0]
 365                           1
 366
 367 (E88)                         OBJSUB = 3456.0
 368
 369 (E89)            GRADSUB = [- 1.90734863E-5, 1.8310547E-4]
 370
 371 (E90)                       EIGEN = - 876.51387
 372
 373 (E91)                       EIGEN = - 35.486028
 374
 375                                                 15
 376 (E92)                   STAPTS  = [S = 0.0, D = --]
 377                               2                 2
 378
 379 (E93)                         OBJSUB = 2362.5
 380
 381 (E94)                       GRADSUB = [0.0, 0.0]
 382
 383 (E95)                           EIGEN = - 84
 384
 385 (E96)                          EIGEN = 202.5
 386 TIME= 18140 MSEC.
 387
 388 (C97) /* It is almost always worth solving the largest possible subset
 389 of the equality constraints for variables that enter them linearly,
 390 then using this solution to eliminate these variables from the
 391 remaining constraints and from the objective.  This is also worth
 392 doing for variables that enter nonlinearly, provided it introduces no
 393 fractional powers in the objective or remaining constraints.
 394 For example, to find the stationary points of  X + 3*Y - 6*Z**2,
 395 subject to  X**2 + Y**2 + Z**2 = 1, it is well worth eliminating Z.
 396
 397 The change of variable for eliminating the inequality constraint
 398 L <= 42, is equivalent to converting the inequality to an equality by
 399 introducing a squared slack variable, solving for L, then eliminating
 400 L from VOLSUB.  From this viewpoint, the "change of variable"
 401 technique is seen to be applicable to a great variety of inequality
 402 constraints, not merely upper and lower bounds as is implied in most
 403 textbooks.  Together with applicable equality constraints, it is
 404 generally worth including in the elimination the largest possible
 405 subset of inequality constraints for which the eliminated variables
 406 enter linearly, up to the number of decision variables minus the number
 407 of equality constraints.
 408
 409 Another technique for treating an inequality constraint is to
 410 solve the problem with the constraint assumed active, and also with
 411 the constraint ignored, checking any stationary points for violation
 412 of the constraint in the latter case: */
 413
 414 STAP(EV(VOLSUB,L=42)) $
 415
 416                                        15
 417 (E98)                              D = --
 418                                        2
 419
 420                                         4725
 421 (E99)                          OBJSUB = ----
 422                                          2
 423
 424 (E100)                          GRADSUB = [0]
 425
 426 (E101)                          EIGEN = - 84
 427 TIME= 1071 MSEC.
 428
 429 (C102) STAP(VOLSUB) $
 430
 431 (E102)                 STAPTS  = [L = 24.0, D = 12.0]
 432                              1
 433
 434 (E103)                         OBJSUB = 3456.0
 435
 436 (E104)                      GRADSUB = [0.0, 0.0]
 437
 438 (E105)                       EIGEN = - 51.633308
 439
 440 (E106)                      EIGEN = - 8.36669242
 441
 442 (E107)                  STAPTS  = [D = 36.0, L = 0.0]
 443                               2
 444
 445 (E108)                          OBJSUB = 0.0
 446
 447 (E109)                      GRADSUB = [0.0, 0.0]
 448
 449 (E110)                       EIGEN = - 58.249224
 450
 451 (E111)                       EIGEN = 22.2492234
 452
 453 (E112)                  STAPTS  = [L = 72.0, D = 0.0]
 454                               3
 455
 456 (E113)                          OBJSUB = 0.0
 457
 458 (E114)                      GRADSUB = [0.0, 0.0]
 459
 460 (E115)                      EIGEN = - 152.498447
 461
 462 (E116)                       EIGEN = 8.49844718
 463
 464 (E117)                  STAPTS  = [D = 0.0, L = 0.0]
 465                               4
 466
 467 (E118)                          OBJSUB = 0.0
 468
 469 (E119)                      GRADSUB = [0.0, 0.0]
 470
 471 (E120)                          EIGEN = - 36
 472
 473 (E121)                           EIGEN = 36
 474 TIME= 11506 MSEC.
 475
 476 (C122) /* The second-order test is inadequate when
 477 constraints have been artificially activated.  For example, the one
 478 eigenvalue for the case  D = 15/2 above is negative, but this
 479 stationary point is actually a saddle.  To see this, we must check the
 480 unactivated gradient at this point to see if it points into or out of
 481 the feasible region: */
 482
 483 EV(GRAD, D=15/2, L=42);
 484 TIME= 142 MSEC.
 485                                        405
 486 (D122)                           [0, - ---]
 487                                         4
 488
 489 (C123) DECSLKMULTS;
 490 TIME= 1 MSEC.
 491 (D123)                             [D, L]
 492
 493 (C124) /*  GRAD is a global variable bound by STAP to the
 494 symbolic expression for the gradient, and DECSLKMULTS is bound to the
 495 variables that the gradient is taken with respect to, in the order
 496 of their components in GRAD.  Thus, -405/4 points in, making the
 497 point a saddle.
 498
 499 The report referenced at the beginning of this demonstration
 500 explains how to generalize this test to more than one constraint.
 501
 502 When there is more than one inequality constraint, each feasible
 503 combination must be activated, unless some additional convexity
 504 requirements are satisfied.  Nevertheless, this combinatorial
 505 activation technique is probably capable of treating
 506 larger problems than either the change-of-variable or the squared-
 507 slack-variable-with-multiplier techniques of treating inequality
 508 constraints.
 509
 510 We have already seen how STAP finds poles of the objective or gradient
 511 only by accident.  The following examples are included as reminders
 512 about other limitations of using calculus to find extrema:*/
 513
 514 STAP(X, [], X**3-Y**2) $
 515
 516 (E124)                   NO STATIONARY POINTS FOUND
 517 TIME= 4727 MSEC.
 518
 519 (C125) /* Lagrange multipliers require the constraints to have con-
 520 tinuous tangents, and this is violated for the above example at X=0,
 521 Y=0, where the objective is a minimum.  A direct use of elimination
 522 would fail too, because it results in an undefined derivative at the
 523 minimum.  In such cases, each piece between discontinuities must be
 524 considered a separate constraint.
 525
 526 The next example has a minimum at X=0, Y=0, where the two
 527 active constraints are parallel, making them linearly dependent.
 528 Such cusps or "wiskers" on the feasible region violate the so-called
 529 "constraint-qualification" requirement; so the extremum is not found:*/
 530
 531 STAP(X, [-Y, Y-X**3]) $
 532
 533 (E125)                   NO STATIONARY POINTS FOUND
 534 TIME= 17706 MSEC.
 535
 536 (C126) /* The following example doesn't have a strictly feasible point;
 537 so it violates Slater's condition; and the extremum at X=0 is not
 538 found: */
 539
 540 STAP(X, X**2) $
 541
 542 (E126)                   NO STATIONARY POINTS FOUND
 543 TIME= 4668 MSEC.
 544
 545 (C127) /* Finally, it is important to remember that unbounded regions
 546 may have nonstationary or asymptotically stationary points
 547 at infinity, which STAP will not find.  Such situations are usually
 548 obvious from qualitative considerations, provided the objective and
 549 constraints are not too complicated and numerous, but the macsyma LIMIT
 550 function can be of help.  However, it is important to remember that a
 551 multivariate limit depends upon the way it is taken.  This is
 552 illustrated by the following example, where you will have to type
 553 NONZERO; in response to two interrupts: */
 554
 555 LIMIT(PEANO, Y, INF);
 556
 557 LIMIT FASL DSK MACSYM BEING LOADED
 558 LOADING DONE
 559 TIME= 280 MSEC.
 560 (D127)                               INF
 561
 562 (C128) LIMIT(PEANO, X, INF);
 563 IS  A B  ZERO OR NONZERO?
 564
 565 NONZERO;
 566 TIME= 402 MSEC.
 567 (D128)                               INF
 568
 569 (C129) RADCAN(EV(PEANO, Y=C+(A**2+B**2)*X**2/2));
 570 TIME= 818 MSEC.
 571                               4      2  2    4   4
 572                             (B  - 2 A  B  + A ) X
 573 (D129)                    - ----------------------
 574                                       4
 575
 576 (C130) LIMIT(%, X, INF);
 577 IS  (B + A) (B - A)  ZERO OR NONZERO?
 578
 579 NONZERO;
 580 TIME= 1085 MSEC.
 581 (D130)                              MINF
 582
 583 TIME= 223257 MSEC.
 584 (D131)                           BATCH DONE
 585