share/calculus/optvar_2.dem

   1 load('optvar)$
   2 /* We may also treat problems with more than one
   3 dependent variable.  For example, the charge, Q, and
   4 displacement from equilibrium, Y, of a simple
   5 electromagnetic loudspeaker, with M, K, L, S, E(T), and T
   6 denoting mass, spring constant, inductance, electro-
   7 magnetic work coefficient, voltage, and time, respectively,
   8 are given by the stationary value of the integral of the
   9 function in the first argument of EL below.  (ref. S.H.
  10 Crandall, D.C. Karnop, E.F. Kurtz Jr., & D.C. Pridmore-
  11 Brown, "Dynamics of Mechanical and Electromechanical
  12 Systems", McGraw-Hill). */
  13 el((m*'d(y,t)**2 - k*y**2 + l*'d(q,t)**2)/2 + s*'d(q,t)*y +
  14 e(t)*q,
  15    [y, q], t);
  16 /* We may simplify these equations, replacing 'D(Q,T)
  17 with the current, I; and we may introduce linear mechanical
  18 resistance B and linear electrical resistance R as follows:
  19 */
  20 %, eval, 'd(q,t)=i, s*i=s*i-b*'d(y,t), e(t)=e(t)-r*i;
  21 /* DESOLVE  also works for some systems of arbitrary-
  22 order constant-coefficient equations; so let's try it with
  23 the following parameter values, excitation E(T), and initial
  24 conditions: */
  25 subst([m=2,k=1,b=1,s=1,l=1,e(t)=sin(t),r=1], %)$
  26 convert(%,
  27 [y,i], t)$
  28 (atvalue(y(t),t=0,0), atvalue(i(t),t=0,0),
  29 atvalue('d(y(t),t),t=0,0))$
  30 desolve(%th(2), [y(t),i(t)]);
  31 /* The functional may also have more than 1 independ-
  32 ent variable.  For example, the general 2-dimensional linear
  33 elliptic partial differential equation is equivalent to the
  34 variational problem:*/
  35 el(a*'d(u,x)**2 + b*'d(u,y)**2 + c*u**2 + e*'d(u,x) +
  36 f*'d(u,y) + g*u,  u, [x,y]);
  37 /* Analytic solutions to even the simplest cases of
  38 this equation are rare, but computer algebraic
  39 manipulation has been used to construct series solutions to
  40 such equations.
  41
  42 Many variational optimization problems are easier to treat
  43 using Pontryagin's Maximum-Principle than by the calculus of
  44 varia- tions.  This is particularly true of optimal control
  45 problems.
  46
  47 As an elementary example, suppose that we have a unit mass
  48 at an arbitrary position X[0] with arbitrary velocity V[0]
  49 at time T=0, and we wish to vary the force  F, subject to
  50 the constraint -1 <= F <= 1, such that the mass arrives at
  51 position X=0 with velocity V=0, in minimum time.  The force
  52 equals the rate of change of momentum; so the motion is
  53 governed by the pair of first-order differential equations:
  54 */
  55 ['d(v,t)=f, 'd(x,t)=v]$
  56 /* Using this list as the argument to the function HAM
  57 results in output of the Hamiltonian followed by
  58 differential equations for the so-called Auxiliary
  59 variables, together with their solution when- ever the
  60 differential equation is of the trivial form: 'D(aux[i],t) =
  61 0, as is often the case: */
  62 cc:ham(%);
  63 /* We may substitute the given solution to the last
  64 differential equation into the one before it, then use
  65 either DESOLVE or ODE2: */
  66 part(%,4),eval$
  67 ode2(ev(part(cc,2),%,eval), aux[1], t);
  68 /* Now we may substitute the values of the auxiliary
  69 variables into the Hamiltonian: */
  70 first(cc), %, eval$
  71 subst(%th(3), %);
  72 /* According to the maximum principle, we should vary
  73 F to maximize this expression for all values of T in the
  74 time-interval of interest.  Considering the constraints on
  75 F, clearly F should be SIGN(C-C[2]*T).  (We could use the
  76 function  STAP in the SHARE file OPTMIZ >  or OPTMIZ LISP
  77 to determine this.)
  78
  79 We have found that every optimal control is  F=1  and/or
  80 F=-1, with at most one switch between them, when T=C/C[2].
  81 Since  F  is piecewise constant, we may combine the two
  82 first-order equations of motion and solve them as follows:
  83 */
  84 ode2('d(x,t,2)=f, x, t);
  85 /* Now we may choose the constants of integration
  86 to sat- isfy any given boundary conditions.  For example,
  87 suppose that we start with  V=0  and  X=1  at  T=0.  Assume
  88 we start with  F=-1: */
  89 ic2(subst(f=-1,%), t=0, x=1, 'd(x,t)=0);
  90 /* Assuming that we terminate with  F=+1: */
  91 ic2(subst(f=1,%th(2)), t=tfinal, x=0, 'd(x,t)=0);
  92 /* To determine  TFINAL  and the time  T  at which  F
  93 switches sign, we may impose the condition that the two
  94 solutions must agree in position and velocity at that time:
  95 */
  96 %,%th(2);
  97 solve([%, d(%,t)]);
  98 /* For the first of these solutions,  0 < T < TFINAL;
  99 so the assumption that  F  switches from  -1  to  1  is
 100 justified.  F is now completely determined as a function of
 101 time, but to represent it with our expression
 102 SIGN(C-C[2]*T): */
 103 solve(subst(first(%), c-c[2]*t), c);
 104 /* As is always the case, one of the integration
 105 constants for the auxiliary equations is redundant; so we
 106 may set C[2] and C to the same arbitrary negative constant,
 107 such as -1.
 108
 109 It is important to remember that so far, we have only
 110 determined an EXTREMAL control.  To prove that it is an
 111 OPTIMAL control, we must prove that an optimal control
 112 exists and that no more-optimal extremal control exists.
 113 However, such considerations are beyond the scope of this
 114 demonstration.  */
 115 "end of demo"$