| blob | 89753d0bdad8af47d33348e9526b22cdecd0577b |
1 (C5) BATCH(OPTVAR, DEMO60, DSK,STOUTE);
3 (C6) /* This macsyma batch file illustrates how to use the
4 functions in OPTVAR > or OPTVAR LISP to help solve classical
5 variational or optimal control problems. These functions
6 use the calculus of variations and Pontryagin's
7 Maximum-Principle to derive the governing differential and
8 algebraic equations; and this demonstration also shows how
9 the governing equations may be solved using the built-in
10 SOLVE function and/or the functions in the SHARE files ODER
11 LISP, and DESOLN LISP. For more details, see the
12 corresponding SHARE text files OPTVAR USAGE, ODER USAGE, or
13 DESOLN USAGE.
15 The illustrative examples here are intentionally elementary.
16 For a more thorough discussion of the mathematical
17 principles demonstrated here, and for results with more
18 difficult examples, see the the ALOHA project technical
19 report by David Stoutemyer, "Computer algebraic manipulation
20 for the calculus of variations, the maximum principle, and
21 automatic control", University of Hawaii, 1974.
23 Many classical variational problems are analogous to special
24 cases of choosing a path which minimizes the transit time
25 through a region for which the speed is a function of
26 position. There are applications in, optics, acoustics,
27 hydrodynamics, and routing of aircraft or ships. In the
28 two-dimensional case, assuming the path may be represented
29 by a single-valued function, Y(X), the transit time between
30 Y(A) and Y(B) is given by the integral of
31 Q(X,Y)*SQRT(1+'D(y,x)**2), from X=A to X=B, where D is an
32 alias for DIFF and Q(X,Y) is the reciprocal of the speed,
33 and where we have used the fact that 'D(arclength,X) =
34 SQRT(1+'D(Y,X)**2). For simplicity, assume Q independent of
35 X. We may use the function EL, previously loaded by
36 LOADFILE(OPTVAR,LISP,DSK, SHARE), to derive the associated
37 Euler-Lagrange equation together with an associated energy
38 and/or momentum integral if they exist: */
40 EL(Q(Y)*SQRT(1+'D(Y,X)**2), Y, X);
42 DY 2
43 Q(Y) (--)
44 DY 2 DX
45 (E6) Q(Y) SQRT((--) + 1) - --------------- = K
46 DX DY 2 0
47 SQRT((--) + 1)
48 DX
50 DY
51 Q(Y) --
52 D DX DY 2 D
53 (E7) -- --------------- = SQRT((--) + 1) (-- Q(Y))
54 DX DY 2 DX DY
55 SQRT((--) + 1)
56 DX
57 TIME= 844 MSEC.
58 (D7) [E6, E7]
60 (C8) /* To expand the Euler-Lagrange equation: */
62 DEPENDENCIES(Y(X))$
63 TIME= 4 MSEC.
65 (C9) %TH(2)[2], EVAL, D;
66 TIME= 300 MSEC.
67 2 2
68 D Y DY 2 D Y
69 Q(Y) --- Q(Y) (--) ---
70 2 DX 2
71 DX DX
72 (D9) --------------- - --------------
73 DY 2 DY 2 3/2
74 SQRT((--) + 1) ((--) + 1)
75 DX DX
77 DY 2 D
78 = SQRT((--) + 1) (-- Q(Y))
79 DX DY
81 (C10) /* The routines in the SHARE file ODE LISP, previously
82 loaded, may solve an expanded first or second order
83 quasi-linear ordinary differential equations; so the
84 equation must be linear in its highest- order derivative.
85 The Euler-Lagrange equation is always of this form, but when
86 given a second-order equation, the ODE solver often returns
87 with a first-order equation which we must quasi-linearize
88 before proceeding; so it is usually most efficient to take
89 advantage of a first integral when one exists, even though
90 it requires a certain amount of manipulation. The SOLVE
91 function is currently somewhat weak with fractional powers;
92 so we must massage the above energy integral before solving
93 for 'D(Y,X): */
95 %TH(3)[1] * SQRT(1+'D(Y,X)**2), EXPAND, EVAL;
96 TIME= 211 MSEC.
97 DY 2
98 (D10) Q(Y) = K SQRT((--) + 1)
99 0 DX
101 (C11) SOLVE(%/K[0], 'D(Y,X));
102 SOLUTION
104 2 2
105 SQRT(Q (Y) - K )
106 DY 0
107 (E11) -- = - -----------------
108 DX K
109 0
111 2 2
112 SQRT(Q (Y) - K )
113 DY 0
114 (E12) -- = -----------------
115 DX K
116 0
117 TIME= 1543 MSEC.
118 (D12) [E11, E12]
120 (C13) ODE2(EV(%[2],EVAL), Y, X);
121 TIME= 9982 MSEC.
122 /
123 [ 1
124 X - K (I (-----------------) DY)
125 0 ] 2 2
126 / SQRT(Q (Y) - K )
127 0
128 (D13) - --------------------------------- = C
129 K
130 0
132 (C14) /* We must specify Q(Y) to proceed further. Q(Y)=1 is
133 associat- ed with the Euclidian shortest path (a straight
134 line), Q(Y)=SQRT(Y) is associated with the stationary
135 Jacobian-action path of a projectile (a parabola), Q(Y)=Y is
136 associated with the minimum-surface body of revolution (a
137 hyperbolic cosine), and Q(Y)=1/SQRT(Y) is associated with
138 the minimum-time path of a falling body, starting at Y=0,
139 measured down, (a cycloid). Of these, the cycloid presents
140 the most difficult integration. In fact, I have never seen
141 the the integration for this case performed directly except
142 by macsyma: */
144 SUBST(Q(Y)=1/SQRT(Y), %) $
145 TIME= 2961 MSEC.
147 (C15) %,INTEGRATE;
148 TIME= 2468 MSEC.
149 1 2
150 SQRT(- + K )
151 Y 0
152 (D15) - (X - K (-------------------
153 0 2 1 2 4
154 K (- + K ) - K
155 0 Y 0 0
157 1 2 1 2
158 LOG(SQRT(- + K ) - K ) LOG(SQRT(- + K ) + K )
159 Y 0 0 Y 0 0
160 + ----------------------- - -----------------------))/K
161 3 3 0
162 2 K 2 K
163 0 0
165 = C
167 (C16) /* This equation may be simplified by combining the
168 LOGs and clearing some fractions, but it is transcendental
169 in Y; so there is no hope for a closed-form representation
170 for Y(X). For completeness we should try the other
171 alternative for 'D(Y,X), but it turns out to lead to the
172 same result.
174 Most authors solve this problem by introducing the change of
175 variable 'D(Y,X)=TAN(T), useful also for other Q(Y), which
176 leads to the para- metric representation: Y=K[0]*COS(T)**2,
177 X=C+K[0]*(T+SIN(T)*COS(T)). Solving the former for T and
178 substituting into the latter gives the representation
179 equivalent to that found directly by macsyma.
181 Some straightforward investigation reveals that C is the
182 initial value of X, but the equation is transcendental in
183 K[0], so K[0] cannot be determined analytically. However,
184 it may be determined to arbitrary numerical accuracy by an
185 iterative algorithm such as that described in the SHARE file
186 ZEROIN USAGE.
188 To prove that the above solution is a strong global minimum,
189 we must prove that such a minimum exists and that no other
190 answer gives a smaller value to the functional. This may be
191 done by the direct approach of Caratheodory, or by checking
192 the classical Weierstrass, Weierstrass-Erdman, Legendre, and
193 and Jacobi necessary and sufficient conditions. Computer
194 algebraic manipulation can assist these investigations, but
195 they are beyond the scope of this demonstration. Instead,
196 to keep the demonstration brief, we will now consider other
197 kinds of problems.
199 Stationary values of a functional may be subject to
200 algebraic or differential constraints of the form
201 F(X,Y(X),'D(Y,X))=0 and/or isoperimetric constraints of the
202 form INTEGRATE(F(X,Y(X),'D(Y,X)),X,A,B) = J, where J is a
203 given constant. Sometimes it is possible to eliminate one or
204 more constraints, substituting them into the functional and
205 any remaining constraints; but if not, we may use Lagrange
206 multipliers. For example, suppose we wish to determine the
207 curve that a string of given length assumes to minimize its
208 potential energy. The potential energy is
209 INTEGRATE(Y*SQRT(1+'D(Y,X)**2), X, A, B), whereas the length
210 is INTEGRATE(SQRT(1+'D(Y,X)**2), X, A, B); so letting MULT
211 be the Lagrange multiplier, our augmented functional is of
212 the form that we have been studying, with Q(Y)=Y+MULT.
213 Consequently, we may simply substitute this in the general
214 integral that we found before: */
216 SUBST(Q(Y)=Y+MULT, %TH(3))$
217 TIME= 119 MSEC.
219 (C17) /* Type NONZERO; in response to the question
220 generated by the following statement: */
222 %,INTEGRATE;
223 IS K ZERO OR NONZERO?
224 0
226 NONZERO;
227 TIME= 1310 MSEC.
228 2
229 (Y + MULT) Y + MULT
230 X - K LOG(SQRT(----------- - 1) + --------)
231 0 2 K
232 K 0
233 0
234 (D17) - -------------------------------------------- = C
235 K
236 0
238 (C18) /* to get Y as a function of X: */
240 (%*K[0]+X)/K[0];
241 TIME= 92 MSEC.
242 2 X + K C
243 (Y + MULT) Y + MULT 0
244 (D18) LOG(SQRT(----------- - 1) + --------) = --------
245 2 K K
246 K 0 0
247 0
249 (C19) EXP(LHS(%))=EXP(RHS(%));
250 TIME= 29 MSEC.
251 X + K C
252 0
253 --------
254 2 K
255 (Y + MULT) Y + MULT 0
256 (D19) SQRT(----------- - 1) + -------- = %E
257 2 K
258 K 0
259 0
261 (C20) SOLVE((%*K[0]-Y-MULT)**2, Y);
262 SOLUTION
264 (E20) Y
266 X 2 X X
267 - -- - C --- + 2 C -- + C
268 K K K
269 0 0 0
270 %E (K %E - 2 %E MULT + K )
271 0 0
272 = ---------------------------------------------------
273 2
274 TIME= 8342 MSEC.
275 (D20) [E20]
277 (C21) %,EVAL,EXPAND,RATPRODEXPAND:TRUE;
278 TIME= 450 MSEC.
279 X X
280 - -- - C -- + C
281 K K
282 0 0
283 K %E K %E
284 0 0
285 (D21) [Y = -------------- + ----------- - MULT]
286 2 2
288 (C22) /* We may rewrite the above expression as: */
290 K[0]*COSH(X/K[0]+C) - MULT $
291 TIME= 65 MSEC.
293 (C23) /* We may substitute this into the integral constraint
294 to get 1 equation relating the constant K[0] and C to the
295 given length L: */
297 SQRT(1+D(%,X)**2);
298 TIME= 99 MSEC.
299 2 X
300 (D23) SQRT(SINH (-- + C) + 1)
301 K
302 0
304 (C24) /* This is simply COSH(X/K[0]+C), so: */
306 L = INTEGRATE(COSH(X/K[0]+C), X, A, B);
308 TIME= 1787 MSEC.
309 K C + B K C + A
310 0 0
311 (D24) L = K SINH(--------) - K SINH(--------)
312 0 K 0 K
313 0 0
315 (C25) /* Given A, B, Y(A), and Y(B), we may substitute into
316 the general solution to get two more relations among K[0], C
317 and MULT, but the three equations are transcendental; so a
318 numerical solution is generally necessary.
320 For completeness, we should also investigate the case of
321 K[0]=0, which yields the solution Y=0; but for brevity, we
322 will omit that analysis.
324 The functional may include derivatives of higher that first
325 order. For example, including the small-deflection energy
326 of bending, shear, and an elastic foundation, the elastic
327 energy of a nonuniform beam with loading W(X) is the
328 integral of the following function: */
330 A(X)*'D(Y,X,2)**2 + B(X)*'D(Y,X)**2 + C(X)*Y**2 + W(X)*Y$
331 TIME= 98 MSEC.
333 (C26) /* to get the Euler-Lagrange equation: */
335 EL(%, Y, X);
337 2 2
338 D DY D D Y
339 (E26) -- 2 B(X) -- - --- 2 A(X) --- = 2 C(X) Y + W(X)
340 DX DX 2 2
341 DX DX
342 TIME= 929 MSEC.
343 (D26) [E26]
345 (C27) /* If the shear & foundation energy are negligible, as
346 is often the case, we may solve this differential equation,
347 for specific A(X) and W(X), by two successive applications
348 of ODE2. For example: */
350 %[1], EVAL,A(X)=X, B(X)=0, C(X)=0, W(X)=1, 'D(Y,X,2)=DY2$
351 TIME= 341 MSEC.
353 (C28) DEPENDENCIES(DY2(X))$
354 TIME= 4 MSEC.
356 (C29) EV(%TH(2),D,EVAL);
357 TIME= 159 MSEC.
358 2
359 D DY2 DDY2
360 (D29) - 2 ----- X - 4 ---- = 1
361 2 DX
362 DX
364 (C30) ODE2(%, DY2, X);
365 TIME= 894 MSEC.
366 X K2 K1
367 (D30) DY2 = - - + -- - --
368 4 X 2
370 (C31) ODE2(SUBST([K1=K3,K2=K4,DY2='D(Y,X,2)],%), Y, X);
371 TIME= 1972 MSEC.
372 3 2
373 K4 X (24 - 24 LOG(X)) + X + 6 K3 X
374 (D31) Y = - ------------------------------------ + K2 X
375 24
377 + K1
379 (C32) /* Now let's impose the boundary conditions
380 Y='D(Y,X)=0 at X=1, Y=0, 'D(Y,X)=1 at X=2: */
382 IC(%, X=1, Y=0, 'D(Y,X)=0);
383 TIME= 4634 MSEC.
384 3 2
385 K4 X (24 - 24 LOG(X)) + X + 6 K3 X
386 (D33) [Y = - ------------------------------------
387 24
389 (4 K3 + 1) X 12 K4 - 3 K3 - 1
390 + ------------ + ----------------]
391 8 12
393 (C34) IC(SUBST([K3=K1,K4=K2],FIRST(%)), X=2, Y=0,
394 'D(Y,X)=1);
395 TIME= 3087 MSEC.
396 49 X (24 - 24 LOG(X)) 3
397 (D35) [Y = - ( - --------------------- + X
398 72 LOG(2) - 48
400 4 (110 LOG(2) - 57)
401 2 (1 - -------------------) X
402 6 (110 LOG(2) - 57) X 18 LOG(2) - 12
403 - ----------------------)/24 + ---------------------------
404 18 LOG(2) - 12 8
406 588 3 (110 LOG(2) - 57)
407 - -------------- + ------------------- - 1
408 72 LOG(2) - 48 18 LOG(2) - 12
409 + -------------------------------------------]
410 12
412 (C36) /* as another specific beam example: */
414 %TH(8),EVAL,A(X)=1,B(X)=2,C(X)=1,W(X)=1$
415 TIME= 271 MSEC.
417 (C37) %[1],D;
418 TIME= 180 MSEC.
419 2 4
420 D Y D Y
421 (D37) 4 --- - 2 --- = 2 Y + 1
422 2 4
423 DX DX
425 (C38) /* We cannot treat this as 2 successive second-order
426 equations, but the function DESOLVE in the SHARE file
427 DESOLN, already loaded, is applicable to some linear
428 arbitrary-order constant coefficient equations. First we
429 must convert the equation so that the dependency of Y is
430 explicitly indicated: */
432 CONVERT(%, Y, X);
433 TIME= 217 MSEC.
434 2 4
435 D D
436 (D38) 4 (--- Y(X)) - 2 (--- Y(X)) = 2 Y(X) + 1
437 2 4
438 DX DX
440 (C39) /* DESOLVE requires all boundary conditions to be at
441 X=0, and it works best if they are specified before calling
442 the function. However, we may overcome this restriction, as
443 illustrated by the following example, where the beam has
444 zero slope and deflection at both ends: */
446 (ATVALUE(Y(X),X=0,0), ATVALUE('D(Y(X),X),X=0,0),
447 ATVALUE('D(Y(X),X,2),X=0,K1), ATVALUE('D(Y(X),X,3),X=0,K2))$
448 TIME= 77 MSEC.
450 (C40) DESOLVE(%TH(2), Y(X));
451 TIME= 5740 MSEC.
452 - X - X
453 (2 K2 - 2 K1 + 1) X %E (K2 + 1) %E
454 (D40) Y(X) = ------------------------- + --------------
455 8 4
457 X X
458 (2 K2 + 2 K1 - 1) X %E (K2 - 1) %E 1
459 + ----------------------- - ------------ - -
460 8 4 2
462 (C41) IC(SUBST(Y(X)=Y,%), X=1, Y=0, 'D(Y,X)=0);
463 TIME= 4435 MSEC.
464 2 2
465 2 (%E - 2 %E - 1) 2 (%E - 2 %E + 1)
466 (D42) [Y = ((------------------ + ------------------ + 1) X
467 2 2
468 2 %E + 4 %E - 2 %E + 2 %E - 1
470 2
471 %E - 2 %E + 1 - X
472 (-------------- + 1) %E
473 2
474 - X %E + 2 %E - 1
475 %E )/8 + --------------------------
476 4
478 2 2
479 2 (%E - 2 %E - 1) 2 (%E - 2 %E + 1) X
480 + (( - ------------------ + ------------------ - 1) X %E )
481 2 2
482 2 %E + 4 %E - 2 %E + 2 %E - 1
484 2
485 %E - 2 %E + 1 X
486 (-------------- - 1) %E
487 2
488 %E + 2 %E - 1 1
489 /8 - ------------------------ - -]
490 4 2
492 (C43) /* We may also treat problems with more than one
493 dependent variable. For example, the charge, Q, and
494 displacement from equilibrium, Y, of a simple
495 electromagnetic loudspeaker, with M, K, L, S, E(T), and T
496 denoting mass, spring constant, inductance, electro-
497 magnetic work coefficient, voltage, and time, respectively,
498 are given by the stationary value of the integral of the
499 function in the first argument of EL below. (ref. S.H.
500 Crandall, D.C. Karnop, E.F. Kurtz Jr., & D.C. Pridmore-
501 Brown, "Dynamics of Mechanical and Electromechanical
502 Systems", McGraw-Hill). */
504 EL((M*'D(Y,T)**2 - K*Y**2 + L*'D(Q,T)**2)/2 + S*'D(Q,T)*Y +
505 E(T)*Q,
506 [Y, Q], T);
508 D DY DQ
509 (E43) -- M -- = -- S - K Y
510 DT DT DT
512 D DQ
513 (E44) -- (S Y + L --) = E(T)
514 DT DT
515 TIME= 1437 MSEC.
516 (D44) [E43, E44]
518 (C45) /* We may simplify these equations, replacing 'D(Q,T)
519 with the current, I; and we may introduce linear mechanical
520 resistance B and linear electrical resistance R as follows:
521 */
523 %, EVAL, 'D(Q,T)=I, S*I=S*I-B*'D(Y,T), E(T)=E(T)-R*I;
524 TIME= 477 MSEC.
525 D DY DY D
526 (D45) [-- M -- = - B -- - K Y + I S, -- (S Y + I L)
527 DT DT DT DT
529 = E(T) - I R]
531 (C46) /* DESOLVE also works for some systems of arbitrary-
532 order constant-coefficient equations; so let's try it with
533 the following parameter values, excitation E(T), and initial
534 conditions: */
536 SUBST([M=2,K=1,B=1,S=1,L=1,E(T)=SIN(T),R=1], %)$
537 TIME= 491 MSEC.
539 (C47) CONVERT(%,
540 [Y,I], T)$
541 TIME= 466 MSEC.
543 (C48) (ATVALUE(Y(T),T=0,0), ATVALUE(I(T),T=0,0),
544 ATVALUE('D(Y(T),T),T=0,0))$
545 TIME= 53 MSEC.
547 (C49) DESOLVE(%TH(2), [Y(T),I(T)]);
548 TIME= 10117 MSEC.
549 SQRT(3) T SQRT(3) T
550 SIN(---------) COS(---------)
551 - T/2 2 2
552 (D50) [Y(T) = %E (-------------- - --------------)
553 SQRT(3) 3
555 - T/2
556 2 SIN(T) COS(T) 8 %E
557 - -------- - ------ + ---------,
558 5 5 15
560 SQRT(3) T SQRT(3) T
561 SIN(---------) COS(---------)
562 - T/2 2 2
563 I(T) = %E ( - -------------- - --------------)
564 SQRT(3) 3
566 - T/2
567 3 SIN(T) COS(T) 8 %E
568 + -------- - ------ + ---------]
569 5 5 15
571 (C51) /* The functional may also have more than 1 independ-
572 ent variable. For example, the general 2-dimensional linear
573 elliptic partial differential equation is equivalent to the
574 variational problem:*/
576 EL(A*'D(U,X)**2 + B*'D(U,Y)**2 + C*U**2 + E*'D(U,X) +
577 F*'D(U,Y) + G*U, U, [X,Y]);
579 D DU D DU
580 (E51) -- (2 B -- + F) + -- (2 A -- + E) = 2 C U + G
581 DY DY DX DX
582 TIME= 1160 MSEC.
583 (D51) [E51]
585 (C52) /* Analytic solutions to even the simplest cases of
586 this equation are rare, but computer algebraic
587 manipulation has been used to construct series solutions to
588 such equations.
590 Many variational optimization problems are easier to treat
591 using Pontryagin's Maximum-Principle than by the calculus of
592 varia- tions. This is particularly true of optimal control
593 problems.
595 As an elementary example, suppose that we have a unit mass
596 at an arbitrary position X[0] with arbitrary velocity V[0]
597 at time T=0, and we wish to vary the force F, subject to
598 the constraint -1 <= F <= 1, such that the mass arrives at
599 position X=0 with velocity V=0, in minimum time. The force
600 equals the rate of change of momentum; so the motion is
601 governed by the pair of first-order differential equations:
602 */
604 ['D(V,T)=F, 'D(X,T)=V]$
605 TIME= 24 MSEC.
607 (C53) /* Using this list as the argument to the function HAM
608 results in output of the Hamiltonian followed by
609 differential equations for the so-called Auxiliary
610 variables, together with their solution when- ever the
611 differential equation is of the trivial form: 'D(aux[i],t) =
612 0, as is often the case: */
614 HAM(%);
616 (E53) AUX V + AUX F
617 2 1
619 D
620 (E54) -- AUX = - AUX
621 DT 1 2
623 D
624 (E55) -- AUX = 0
625 DT 2
627 (E56) AUX = C
628 2 2
629 TIME= 339 MSEC.
630 (D56) [E53, E54, E55, E56]
632 (C57) /* We may substitute the given solution to the last
633 differential equation into the one before it, then use
634 either DESOLVE or ODE2: */
636 %[4],EVAL$
637 TIME= 25 MSEC.
639 (C58) ODE2(EV(%TH(2)[2],%,EVAL), AUX[1], T);
640 TIME= 465 MSEC.
641 (D58) AUX = C - C T
642 1 2
644 (C59) /* Now we may substitute the values of the auxiliary
645 variables into the Hamiltonian: */
647 %TH(3)[1], %, EVAL$
648 TIME= 126 MSEC.
650 (C60) SUBST(%TH(3), %);
651 TIME= 35 MSEC.
652 (D60) C V + F (C - C T)
653 2 2
655 (C61) /* According to the maximum principle, we should vary
656 F to maximize this expression for all values of T in the
657 time-interval of interest. Considering the constraints on
658 F, clearly F should be SIGN(C-C[2]*T). (We could use the
659 function STAP in the SHARE file OPTMIZ > or OPTMIZ LISP
660 to determine this.)
662 We have found that every optimal control is F=1 and/or
663 F=-1, with at most one switch between them, when T=C/C[2].
664 Since F is piecewise constant, we may combine the two
665 first-order equations of motion and solve them as follows:
666 */
668 ODE2('D(X,T,2)=F, X, T);
669 TIME= 618 MSEC.
670 2
671 F T
672 (D61) X = ---- + K2 T + K1
673 2
675 (C62) /* Now we may choose the constants of integration
676 to sat- isfy any given boundary conditions. For example,
677 suppose that we start with V=0 and X=1 at T=0. Assume
678 we start with F=-1: */
680 IC(SUBST(F=-1,%), T=0, X=1, 'D(X,T)=0);
681 TIME= 635 MSEC.
682 2
683 T
684 (D63) [X = 1 - --]
685 2
687 (C64) /* Assuming that we terminate with F=+1: */
689 IC(SUBST(F=1,%TH(2)), T=TFINAL, X=0, 'D(X,T)=0);
690 TIME= 4104 MSEC.
691 2 2
692 TFINAL T
693 (D65) [X = ------- - T TFINAL + --]
694 2 2
696 (C66) /* To determine TFINAL and the time T at which F
697 switches sign, we may impose the condition that the two
698 solutions must agree in position and velocity at that time:
699 */
701 SUBST(%, FIRST(%TH(2)));
702 TIME= 36 MSEC.
703 2 2 2
704 TFINAL T T
705 (D66) ------- - T TFINAL + -- = 1 - --
706 2 2 2
708 (C67) SOLVE([%, D(%,T)]);
709 TIME= 1309 MSEC.
710 (D67) [[TFINAL = 2.0, T = 1.0], [TFINAL = - 2.0, T = - 1.0]]
712 (C68) /* For the first of these solutions, 0 < T < TFINAL;
713 so the assumption that F switches from -1 to 1 is
714 justified. F is now completely determined as a function of
715 time, but to represent it with our expression
716 SIGN(C-C[2]*T): */
718 SOLVE(SUBST(%[1], C-C[2]*T), C);
719 SOLUTION
721 (E68) C = C
722 2
723 TIME= 106 MSEC.
724 (D68) [E68]
726 (C69) /* As is always the case, one of the integration
727 constants for the auxiliary equations is redundant; so we
728 may set C[2] and C to the same arbitrary negative constant,
729 such as -1.
731 It is important to remember that so far, we have only
732 determined an EXTREMAL control. To prove that it is an
733 OPTIMAL control, we must prove that an optimal control
734 exists and that no more-optimal extremal control exists.
735 However, such considerations are beyond the scope of this
736 demonstration. */
739 TIME= 78843 MSEC.
741 (C70) CLOSEFILE(OPTVAR,OUT60)$