# Dowling21.wxmx:   Optimal Control Theory

load(draw)$  set_draw_defaults(line_width=2, grid = [2,2], point_type = filled_circle,
          head_type = 'nofilled, head_angle = 20, head_length = 0.5,
             background_color = light_gray,  draw_realpart=false)$           

load ("Econ2.mac");

## Preface

Dowling21.wxmx is one of a number of wxMaxima files available in the section
        Economic Analysis with Maxima
on my CSULB webpage. 

In Dowling21.wxmx, we use Maxima to discuss the maximization of an integral over a finite period
of time containing both a time dependent state variable x(t) subject to constraints and a time 
dependent control variable y(t) whose forms are sought using the calculus of variations, following 
Dowling's Chapter 21: "Optimal Control Theory". The problems considered are restricted to the
context of continuous time, a finite time horizon, and fixed endpoint constraints on the state
variable.

We have changed some of the symbols used in particular problems. We use p(t) to represent the
co-state variable instead of Dowling's λ(t). An approximate pdf translation (using Microsoft print 
to pdf) is available as Dowling21fit.pdf. That pdf file can be searched using Ctrl-F.

Ted Woollett
https://home.csulb.edu/~woollett/
woollett@charter.net
April 7, 2022

## References

Introduction to Mathematical Economics, 3rd ed., Edward T. Dowling, 2012, Schaum's 
Outline Series, McGraw-Hill.

Fundamental Methods of Mathematical Economics, Alpha C. Chiang and Kevin Wainwright, 
4th ed., 2005, McGraw-Hill

Elements of Dynamic Optimization, Alpha C. Chiang, 2012, Waveland Press

## Terminology and Statement of the Problem

Roughly quoting Dowling:
"In optimal control theory, the aim is to find the optimal time path for a 'control variable', which
we shall denote as y. The variable for which we previously sought an optimal time path in the
calculus of variations [Ch. 20], known as a 'state variable', we shall continue to designate as x.
The goal of optimal control theory is to select a stream of values over time for the control variable
that will optimize a functional (J) subject to the constraints set on the state variable."

"Optimal control theory problems involving 1) continuous time,  2) a finite time horizon, 
and  3) fixed endpoints are generally written:

    Maximize the functional:    J = integrate (f[x(t), y(t), t], t, 0, T)

                  subject to:   dx/dt = g[x(t), y(t), t],   x(0) = x0,   x(T) = xT,                                  (21.1)

where J = the value of the functional to be optimized; y(t) = the control variable, so called because
its value is selected or 'controlled' to optimize J; x(t) = the state variable, which changes over time
according to the differential equation implied by the constraint dx/dt = g[x(t), y(t), t] , and whose 
value is 'controlled' by the presence of the control variable y(t) in that constraint; and 
t = time. The solution to an optimal control problem demarcates the optimal dynamic time path for 
the control variable y(t)."

## Necessary Conditions for Max. with Fixed End Pts.

"Dynamic optimization of a functional subject to a constraint on the state variable in optimal 
control involves a 'Hamiltonian function H' similar to the Lagrangian function in concave 
programming.

In terms of (21.1) the Hamiltonian is defined as"

      H[x(t), y(t), p(t), t] = f[x(t), y(t),t ] + p(t)*g[x(t), y(t), t]                                                     (21.2)

where p(t) is called the 'co-state variable' [also the 'generalized momentum', the 'marginal 
valuation' of x(t), the 'shadow price', and 'how much a unit increment in x at time t contributes 
to the optimal  objective functional J' ".
 
Dowlingl uses λ(t) instead of p(t) for the co-state variable. In a problem in which the state 
variable is  a price, we can use symbols p(t) as the state variable and λ(t) as the co-state 
variable as a more  natural notation.

"Similar to the Lagrangian multiplier, the co-state variable p(t) estimates the marginal value or 
shadow price of the associated state variable x(t)."

"Assuming the Hamiltonian is differentiable in y and strictly concave, so there is an interior 
solution and not an endpoint solution, the necessary conditions for maximization are"

1.                                                ∂H/∂y = 0

2.         a)              dp/dt =  - ∂H/∂x                    b)       dx/dt =  ∂H/∂p

3.         a)               x(0) = x0                             b)        x(T) = xT

"The first two conditions are known as 'the [Pontryagin] maximum principle' (PMP) and the 
third is called the 'boundary condition.' The two equations of motion in the second condition 
are  generally referred to as the 'Hamiltonian system' or  the 'canonical system'."

"For minimization (instead of maximization) the objective functional can simply be multiplied 
by -1, as in concave programming. If the solution involves an end point, ∂H/∂y  need not equal 
zero in the first condition, but H must still be maximized with respect to y. See Ch. 13, Ex. 9 
and Fig. 13-1(b)-(c) for clarification."

## Sufficient Conditions for Maximization

Quoting Dowling:
"Assuming the maximum principle representing the necessary conditions for maximization in
optimal control are satisfied, the sufficiency conditions will be fulfilled if

   1. Both the objective functional [integrand] f(x,y,t) and the constraint g(x,y,t) are differentiable
and jointly concave in x and y, and

   2. p(t) >= 0, if the constraint g(x,y,t) is nonlinear in x or y. If the constraint is linear, then p(t) 
may assume any sign."

"Linear functions are always both concave and convex, but neither strictly concave nor strictly
convex. For nonlinear functions, an easy test for joint concavity is the simple discriminant test. 

Given the discriminant matrix D of the second order derivatives of a function [say, f(x,y)]:
       D = matrix ([fxx, fxy], [fyx, fyy]),
a function f(x,y) will be 'strictly concave' if the discriminant is 'negative definite':
       D1 = D[1,1] = fxx < 0  and  D2 = determinant(D) > 0
and 'simply concave' if the discriminant is 'negative semidefinite':
       D1 = fxx <= 0    and   D2 >= 0.

A negative definite discriminant indicates a global maximum and is, therefore, always sufficient
for a maximum. A negative semidefinite discriminant is indicative of a relative maximum and is
sufficient for a maximum if the test is conducted for every possible ordering of the variables (x,y)
with similar results."

## ConcaveTest(f, x, y)

The Maxima function ConcaveTest (f, x, y), defined in Econ2.mac, determines if f(x,y) is jointly
concave in x and y using discriminant tests.
 
The matrix d1 which is displayed is equivalent to
matrix ([fxx, fxy], [fyx, fyy]).  d11 is fxx and dd1 is determinant(d1).

 The matrix d2 which is displayed is equivalent to matrix ([fyy, fyx], [fxy, fxx]). d21 is fyy and
dd2 is determinant (d2).

ConcaveTest (4*x - 5*y^2, x, y);

ConcaveTest (8*y, x, y)$

## Problems with Fixed End Points

### Example 1

Use the necessary conditions for maximization to solve the problem:

     Maximize the functional:   J = integrate ( 4*x - 5*y^2, t, 0, 3)

        subject to:     dx/dt = 8*y,    x(0) = 2,  and  x(3) = 117.2.

### Interior Solution equations of motion

A.   Set up the Hamiltonian H + p*g:      

H : 4*x - 5*y^2 + p*8*y;

B. Assuming an interior solution, apply the maximum principle (1 & 2):
    1.   ∂H/∂y = 0: 

dy : diff (H, y);

sy : solve(dy,y);

Let sy stand for y(t) in terms of p(t):

sy : rhs (sy[1]);

so y(t) = 4*p(t)/5.

2.  use dx/dt = diff(H,p) and dp/dt = -diff(H,x).

Let xd stand for dx/dt =  ∂H/∂p.

xd : diff(H,p);

Make use of our result that sy stands for  4*p/5.

xd : at (xd, y = sy);

So dx/dt = 32*p/5, which means to find x(t) we need to know p(t).
Let pd stand for dp/dt.

pd : - diff (H, x);

We can immediately integrate this first order ode:
      dp/dt = -4, 
to get p(t) in terms of an arbitrary constant k1.

Let sp stand for p(t) in terms of the constant k1.

sp : integrate(pd, t) + k1;

Now replace p(t) in expression dx/dt (xd).

xd : at(xd, p = sp), expand;

We can then integrate dx/dt = F(t) in terms of another arbitrary integration constant k2.
Let sx stand for x(t) in terms of k1 and k2.

sx : integrate (xd, t) + k2;

Apply the given boundary condtions at t = 0 and t = 3 to determine the values of k1 and k2.

ksoln : solve ([ at(sx, t = 0) = 2, at (sx, t = 3) = 117.2] );

ksoln : ksoln[1];

Update sx, which is x(t) to reflect the consequences of the endpoint boundary conditions.

sx : at (sx, ksoln );

Obviously, x(0) = 2.

at (sx, t = 0);

at (sx, t = 3), numer;

Update sp which is p(t), the co-state variable, to reflect the consequences of the boundary
conditions:

sp : at(sp, ksoln );

Obviously, p(0) = 12, p(3) = 0.

at(sp, t = 3);

sy is y(t), the control variable, which is proportional to p. Replace p with sp.

sy : at (sy, p = sp), expand;

The optimal path of the control variable y(t) is linear, starting at (0, 9.6) and ending at (3,0), with
a slope of -16/5 = - 3.2.

Value of y(0) in terms of a floating point number:

at(sy, t = 0), numer;

Value of y(3) in terms of a floating point number. Maxima's numerical calculations carry along 
16 digits, so the following answer is approx. 0.

at(sy, t = 3), numer;

Doing the evaluation symbolically, we get exactly 0.

at(sy, t = 3);

So y(0) = 9.6, y(3) = 0.

### Plot of solutions of Example 1

wxdraw2d (xlabel = "t", key = " state var x(t)", yrange = [0, 140], key_pos = top_left,
    explicit (sx, t, 0, 3), color = red, key = "control var y(t)", explicit (sy, t, 0, 3),
    color = brown, key = " co-state var p(t)", explicit (sp, t, 0, 3))$

If we relace x,y,and p with our solutions x(t), y(t) and p(t) we can evaluate H(t):

Ht : at (H, [x = sx, y = sy, p = sp]), expand;

We evaluate the extremal value of J = integrate(f, t, 0, 3), all we need are the x(t) and y(t) 
solutions. We first define f as an expression in terms of x and y, and then replace x and y 
by sx and sy,  respectively [redefining f as an expression depending only on t].

f : 4*x - 5*y^2;
f : at (f, [x = sx, y = sy]), expand;

Integrating this expression over the time interval [0, 3] then gives us the extremum value of J.

Jextr : integrate (f, t, 0, 3);

float (%);

### Sufficient Conditions for Example 1

Check the sufficient conditions for a maximum in J. 

1. The constraint function g = 8*y is linear in y and is thus always both concave and convex.
That implies we don't need to worry about the sign of p(t), the co-state variable.  

2. We need the integrand f(x,y) of the functional J to be jointly concave in x and y.

ConcaveTest (f, x, y), discussed above, finds that f(x,y) is jointly (simply) concave in both x and y, 
and the maximum is therefore a relative maximum.

ConcaveTest (4*x - 5*y^2, x, y)$

### Prob. 21.1

Maximize the functional:    J = integrate (6*x - 4*y^2, t, 0, 2)
  subject to:     dx/dt =  16*y,  x(0) = 24,  x(2) = 408.

A.   Set up the Hamiltonian, H = f + p*g.       

H : 6*x - 4*y^2 + p*16*y;

B. Assuming an interior solution, apply the maximum principle (1 & 2):
    1.   ∂H/∂y = 0: 

dy : diff (H, y);

sy : solve(dy,y);

sy : rhs (sy[1]);

so y(t) = 2*p(t)

2.  use dx/dt = diff(H,p) and dp/dt = -diff(H,x).

xd : diff(H,p);

xd : at (xd, y = sy);

dx/dt depends on p(t).

pd : - diff (H, x);

We can immediately integrate this first order ode:
      dp/dt = - 6, 
to get p(t) in terms of an arbitrary constant k1.

sp is p(t) in terms of constant k1.

sp : integrate(pd, t) + k1;

Replace p(t) in the expression xd (which stands for dx/dt)

xd : at(xd, p = sp), expand;

We can then integrate dx/dt = F(t) in terms of another arbitrary integration constant k2.
Let sx stand for x(t) in terms of k1 and k2.

sx : integrate (xd, t) + k2;

Apply the boundary condtions at t = 0 and t = 3 to determine the values of k1 and k2.

ksoln : solve ([ at(sx, t = 0) = 24, at (sx, t = 2) = 408] );

ksoln : ksoln[1];

Update sx, which is x(t).

sx : at (sx, ksoln );

Check  x(0) = 24, x(2) = 408.

[ at (sx, t = 0), at (sx, t = 2)];

Update sp which is p(t), the co-state variable.

sp : at(sp, ksoln );

We now find that p(0) = 12, p(2) = 0.

[at (sp, t = 0),  at(sp, t = 2) ];

sy is y(t), the control variable, which is proportional to p, so update sy:

sy : at (sy, p = sp), expand;

The optimal path of the control variable y(t) is linear, starting at (0, 24) and ending at (2,0), with
a slope of -12.

Check initial and final values of the control variable y(t):

[ at(sy, t = 0), at(sy, t = 2)];

Value of y(2).

wxdraw2d (xlabel = "t", key = " state var x(t)", yrange = [0, 450], key_pos = top_left,
    explicit (sx, t, 0, 2), color = red, key = "control var y(t)", explicit (sy, t, 0, 2),
    color = brown, key = " co-state var p(t)", explicit (sp, t, 0, 2))$

Check the sufficient conditions for a maximum in J for Prob. 21.1.

1. The constraint function g = 16*y is linear in y and is always both concave and convex, and we 
hence don't need to worry about the sign of p(t), the co-state variable.  

2. We need the integrand f(x,y) of the functional J to be jointly concave in x and y.
ConcaveTest (f, x, y) finds f(x,y) is jointly (simply) concave in both x and y, and the maximum is
therefore a relative maximum. (Note: f above ended up as an expression depending only of t.)

ConcaveTest (6*x - 4*y^2, x, y)$

### Prob. 21.2

Maximize the functional:    J = integrate (5*x + 3*y - 2*y^2, t, 0, 1)
subject to:     dx/dt =  6*y,  x(0) = 7,  x(1) = 70.

A.   Set up the Hamiltonian, H = f + p*g.       

H : 5*x + 3*y - 2*y^2 + p*6*y;

B. Assuming an interior solution, apply the maximum principle (1 & 2):
    1.   ∂H/∂y = 0: 

dy : diff (H, y);

sy : solve(dy,y), expand;

sy is y(t) solution in terms of p(t).

sy : rhs (sy[1]);

2. Require: dx/dt = diff(H,p) and dp/dt = -diff(H,x).

xd : diff(H,p);

xd : at (xd, y = sy), expand;

dx/dt  depends on p(t).

pd : - diff (H, x);

We can immediately integrate this first order ode:
      dp/dt = -5, 
to get p(t) in terms of an arbitrary constant k1.

sp is p(t) in terms of the constant k1.

sp : integrate(pd, t) + k1;

replace p(t) in expression xd, which stands for dx/dt.

xd : at(xd, p = sp), expand;

We can then integrate dx/dt = F(t) in terms of another arbitrary integration constant k2.
sx is x(t) in terms of k1 and k2.

sx : integrate (xd, t) + k2;

Apply the boundary condtions at t = 0 and t = 1 to determine the values of k1 and k2.

ksoln : solve ([ at(sx, t = 0) = 7, at (sx, t = 1) = 70] );

ksoln : ksoln[1];

Update sx, which is x(t).

sx : at (sx, ksoln );

float (sx);

[at(sx, t = 0), at (sx, t = 1) ];

Update sp which is p(t), the co-state variable.

sp : at(sp, ksoln );

[at (sp, t = 0), at(sp, t = 1) ];

sy stands for y(t), the control variable, which is proportional to p, so update sy:

sy : at (sy, p = sp), expand;

float (sy);

[at (sy, t = 0),  at(sy, t = 1)], numer;

The optimal path of the control variable y(t) is linear, starting at (0, 14.25) and ending at 
(1, 6.75), with a slope of - 7.5.

wxdraw2d (xlabel = "t", key = " state var x(t)", yrange = [0, 100], key_pos = top_left,
    explicit (sx, t, 0, 1), color = red, key = "control var y(t)", explicit (sy, t, 0, 1),
    color = brown, key = " co-state var p(t)", explicit (sp, t, 0, 1))$

Check the sufficient conditions for a maximum in J for Prob. 21.2.

1. The constraint function g = 6*y is linear in y and is always both concave and convex, and we 
hence don't need to worry about the sign of p(t), the co-state variable.  

2. We need the integrand f(x,y) of the functional J to be jointly concave in x and y.
ConcaveTest (f, x, y) finds f(x,y) is jointly (simply) concave in both x and y, and the maximum is
therefore a relative maximum.

ConcaveTest (5*x + 3*y - 2*y^2, x, y)$

## Neccessary Conditions for Max. with a Free End Point

Quoting Dowling:
"The general format for an optimal control problem involving continuous time with a finite time
horizon and a free endpoint is:

    Maximize the functional  J = integrate (f[x(t), y(t), t], t, 0, T)
              subject to:  dx/dt = g[x(t), y(t), t],
                                 x(0) = x0,   x(T) =  free and unrestricted.

"Assuming an interior solution, the first two conditions for maximization, comprising the
maximum condition, remain the same, but the third or boundary condition changes."

The Hamiltonian function is:   H[x(t), y(t), p(t), t] = f[x(t), y(t),t ] + p(t)*g[x(t), y(t), t]  with
p(t) the co-state variable, x(t) the state variable, y(t) the control variable.

"Assuming the Hamiltonian is differentiable in y and strictly concave so there is an interior 
solution and not an endpoint solution, the necessary conditions for maximization are

1.                                                ∂H/∂y = 0

2.         a)              dp/dt =  - ∂H/∂x                    b)       ∂x/∂t =  ∂H/∂p

3.         a)               x(0) = x0                             b)        p(T) = 0

"The very last condition, p(T) = 0 is known as the 'transversality condition' for a free endpoint.
The rationale for the transversality condition follows straightforward from what we learned in
concave programming. If the value of x at T is free to vary, the constraint must be nonbinding 
and the shadow price p evaluated at T must equal 0, i.e., p(T) = 0. See Prob. 21.4 to 21.6."

"For minimization (instead of maximization), the objective functional can simply be multiplied by 
-1, as in concave programming. If the solution involves an end point, ∂H/∂y  need not equal zero 
in the first condition, but H must still be maximized with respect to y."

### Example 3

  Maximize the functional J = integrate ( 3*x - 2*y^2,  t, 0, 2)
      subject to:   dx/dt = 8*y,  x(0) = 5,  x(2) = free and unrestricted.

A.   Set up the Hamiltonian, H = f + p*g.       

H : 3*x - 2*y^2 + p*8*y;

B. Assuming an interior solution, apply the maximum principle (1 & 2):
    1.   ∂H/∂y = 0: 

dy : diff (H, y);

sy : solve(dy,y);

sy is y(t) solution in terms of p(t).

sy : rhs (sy[1]);

2.  use dx/dt = diff(H,p) and dp/dt = -diff(H,x).

xd : diff(H,p);

xd : at (xd, y = sy);

dx/dt  depends on p(t).

pd : - diff (H, x);

We can immediately integrate this first order ode:
      dp/dt = -3, 
to get p(t) in terms of an arbitrary constant k1.

sp is p(t) in terms of constant k1.

sp : integrate(pd, t) + k1;

replace p(t) in expression dx/dt (xd).

xd : at(xd, p = sp), expand;

We can then integrate dx/dt = F(t) in terms of another arbitrary integration constant k2.
sx is x(t) in terms of k1 and k2.

sx : integrate (xd, t) + k2;

Apply the boundary condtions at t = 0 and t = 2 to determine the values of k1 and k2.
Instead of x(2) value we use p(2) value (transversality condition).

ksoln : solve ([ at(sx, t = 0) = 5, at (sp, t = 2) = 0] );

ksoln : ksoln[1];

Update sx, which is x(t).

sx : at (sx, ksoln );

[at (sx, t = 0),  at (sx, t = 2)];

Update sp which is p(t), the co-state variable.

sp : at(sp, ksoln );

p(0) = 6, p(2) = 0.

[at (sp, t = 0),  at(sp, t = 2) ];

sy is y(t), the control variable, which is proportional to p, so update sy:

sy : at (sy, p = sp), expand;

The optimal path of the control variable y(t) is linear, starting at (0, 12) and ending at 
(2, 0), with a slope of - 6.

wxdraw2d (xlabel = "t", key = " state var x(t)", yrange = [0, 120], key_pos = top_left,
    explicit (sx, t, 0, 2), color = red, key = "control var y(t)", explicit (sy, t, 0, 2),
    color = brown, key = " co-state var p(t)", explicit (sp, t, 0, 2))$

Check the sufficient conditions for a maximum in J for Example 3.

1. The constraint function g = 8*y is linear in y and is always both concave and convex, and we 
hence don't need to worry about the sign of p(t), the co-state variable.  

2. We need the integrand f(x,y) of the functional J to be jointly concave in x and y.
ConcaveTest (f, x, y) finds f(x,y) is jointly (simply) concave in both x and y, and the maximum is
therefore a relative maximum.

ConcaveTest (3*x - 2*y^2, x, y)$

### Prob. 21.4

Maximize the functional  J = integrate (8*x - 10*y^2, t, 0, 4)]
   subject to:  dx/dt = 24*y,  x(0)  = 7,  x(4) = free and unrestricted.

A.   Set up the Hamiltonian, H = f + p*g.       

H : 8*x - 10*y^2 + p*24*y;

B. Assuming an interior solution, apply the maximum principle (1 & 2):
    1.   ∂H/∂y = 0: 

dy : diff (H, y);

sy : solve(dy,y);

sy is y(t) solution in terms of p(t).

sy : rhs (sy[1]);

2.  use dx/dt = diff(H,p) and dp/dt = -diff(H,x).

xd : diff(H,p);

xd : at (xd, y = sy);

dx/dt  depends on p(t).

pd : - diff (H, x);

We can immediately integrate this first order ode:
      dp/dt = -8, 
to get p(t) in terms of an arbitrary constant k1.

sp is p(t) in terms of constant k1.

sp : integrate(pd, t) + k1;

replace p(t) in expression dx/dt (xd).

xd : at(xd, p = sp), expand;

We can then integrate dx/dt = F(t) in terms of another arbitrary integration constant k2.
sx is x(t) in terms of k1 and k2.

sx : integrate (xd, t) + k2;

Apply the boundary condtions at t = 0 and t = 4 to determine the values of k1 and k2.
Instead of x(4) value we use p(4) = 0 (transversality condition).

ksoln : solve ([ at(sx, t = 0) = 7, at (sp, t = 4) = 0] );

ksoln : ksoln[1];

Update sx, which is x(t).

sx : at (sx, ksoln );

float(sx);

[ at (sx, t = 0), at (sx, t = 4) ];

float(%);

Update sp which is p(t), the co-state variable.

sp : at(sp, ksoln );

[at (sp, t = 0), at (sp, t = 4)];

sy is y(t), the control variable, which is proportional to p, so update sy:

sy : at (sy, p = sp), expand;

float(sy);

at (sy, t = 4);

The optimal path of the control variable y(t) is linear, starting at (0, 38.4) and ending at 
(4, 0), with a slope of - 9.6.

wxdraw2d (xlabel = "t", key = " state var x(t)", yrange = [0, 2e3], key_pos = top_left,
    explicit (sx, t, 0, 4), color = red, key = "control var y(t)", explicit (sy, t, 0, 4),
    color = brown, key = " co-state var p(t)", explicit (sp, t, 0, 4))$

Check the sufficient conditions for a maximum in J for Prob. 21.4.

1. The constraint function g = 24*y is linear in y and is always both concave and convex, and we 
hence don't need to worry about the sign of p(t), the co-state variable.  

2. We need the integrand f(x,y) of the functional J to be jointly concave in x and y.
ConcaveTest (f, x, y) finds f(x,y) is jointly (simply) concave in both x and y, and the maximum is
therefore a relative maximum.

ConcaveTest (8*x - 10*y^2, x, y)$

### Prob. 21.6  Use of Zindef_float (A, B)

Maximize  the functional J  = integrate (4*y - y^2 - x - 2*x^2, t, 0, 1)
   subject to:  dx/dt = x + y,  x(0) = 6.15,  x(1) = free and unrestricted.

A.   Set up the Hamiltonian, H = f + p*g.       

H : 4*y - y^2 - x - 2*x^2 + p*(x + y);

B. Assuming an interior solution, apply the maximum principle (1 & 2):
    1.   ∂H/∂y = 0: 

dy : diff (H, y);

sy : solve(dy,y), expand;

sy is y(t) solution in terms of p(t).

sy : rhs (sy[1]);

2.  use dx/dt = diff(H,p) and dp/dt = -diff(H,x).

xd : diff(H,p);

xd : at (xd, y = sy);

dx/dt  depends on both x(t) and p(t).

pd : - diff (H, x);

dp/dt depends on both x(t) and p(t). 
So we have a pair of coupled first order odes to solve.
Let Z = cvec ([x, p]) and write the set of first order odes in matrix form:

       dZ/dt = A . Z + B.

The Maxima function Zindef_float (A, B), defined in Econ2.mac, returns the indefinite solution 
as a matrix column vector.

A : matrix ( [1, 1/2], [4, -1] );
B : cvec ( [2, 1]);

eigenvalues(A);

float(%);

Zs : Zindef_float (A, B);

xs : Zs[1,1];
ps : Zs[2,1];

To determine the constants %k1 and %k2 we use x(0) = 6.15 and p(1) = 0, the latter 
corresponding to x(1) = free and unrestricted.

ksoln : solve ([at(xs, t = 0) = 6.15, at (ps, t = 1) = 0]), numer;

sx : at (xs, ksoln[1]);
sp : at (ps, ksoln[1]);

To compare with Dowling's solutions we need to replace some pieces of this with floating point
values, which we can do with the Maxima function substpart (zz, expr, n1, n2,...).

part (sp, 3);

sp : substpart (float (-7/3), sp, 3);

[at (sp, t = 0), at (sp, t = 1)];

sx : substpart (float (- 5/6), sx, 3);

at (sx, t = 1);

Finally, we can update sy, which is y(t), and which depends on p(t):

sy : at (sy, p = sp), expand;

[at (sy, t = 0), at (sy, t = 1)];

wxdraw2d (xlabel = "t", key = " state var x(t)", yrange = [-40, 20], key_pos = top_left,
    explicit (sx, t, 0, 1), color = red, key = "control var y(t)", explicit (sy, t, 0, 1),
    color = brown, key = " co-state var p(t)", explicit (sp, t, 0, 1), color = black,
    key = "", line_width = 1, explicit (0, t, 0, 1))$

Check the sufficient conditions for a maximum in J for Prob. 21.6.

1. The constraint function g = x + y is linear in x and y and is always both concave and convex, 
and hence we don't need to worry about the sign of p(t), the co-state variable.  

2. We need the integrand f(x,y) of the functional J to be jointly concave in x and y.
ConcaveTest (f, x, y) finds f(x,y) is strictly concave, and the maximum is therefore a 
global maximum.

ConcaveTest (4*y - y^2 - x - 2*x^2, x, y)$

## Case of Endpoint Inequality Constraints

Quoting Dowling:
"If the terminal value of the state variable is subject to an inequality constraint, x(T) >= xmin,
the optimal value x_opt(T) may be chosen freely as long as it does not violate the value set by the
constraint xmin. If x_opt(T) > xmin, the constraint is nonbinding and the problem reduces to a
free endpoint problem. So
     p(T) = 0   when x_opt(T) > xmin.

If x_opt(T) < xmin, the constraint is binding and the optimal solution will involve setting x(T) = xmin,
which is equivalent to a fixed-end problem with
   p(T) >= 0   when  x_opt(T) = xmin.

For conciseness, the endpoint conditions are sometimes reduced to a single statement analogous
to the Kuhn-Tucker condition:
     p(T) >= 0,    x(T) >= xmin,    and  [x(T) - xmin]*p(T) = 0."

"In practice, solving problems with inequality constraints on the endpoints is straightforward.

First solve the problem as if it were a free endpoint problem. If the optimal value of the state
variable is greater than the minimum required by the endpoint condition, i.e., if x_opt(T) >= xmin,
the correct solution has been found. If x_opt(T) < xmin, set the terminal endpoint equal to the
value of the constraint, x(T) = xmin, and solve as a fixed endpoint problem."

### Example 5

Maximize J = integrate (3*x - 2*y^2, t, 0, 2)
  subject to:  dx/dt = 8*y,  x(0) = 5,   x(2) >= 95.

To solve an optimal control problem involving an inequality constraint, solve it first as an 
unconstrained problem with a free end point. This we did in Example 3 where we found
x_opt(t) = - 24*t^2 + 96*t + 5,

at (-24*t^2 + 96*t + 5, t = 2);

Since x_opt(2) = 101  > 95, the free endpoint solution satisfies the inequality constraint required,
and the constraint is nonbinding, so we have found the correct solution for the state variable x.
The solution for the control variable y(t) found in Ex. 3 also is still valid:
  y(t) = 12 - 6*t.

### Example 6

Redo the same problem in Ex. 5 above with the new required boundary conditions:
  x(0) = 5,  x(2) >= 133. 

Step 1 is to solve Ex. 5 as an unconstrained problem with a free endpoint, which we did in Ex. 3, 
finding x_opt(2) = 101, which is less than the required value of 133. Hence we need to solve 
as a fixed end point problem with x(2) = 133.

From our work above on Ex. 3, we found y = 2*p, dx/dt = 16*p, dp/dt = -3, p = k1 - 3*t,
dx/dt = 16*k1 - 48*t, sx = -24*t^2 + 16*k1*t + k2.

sx : -24*t^2 + 16*k1*t + k2;
ksoln : solve ([ at(sx, t = 0) = 5, at (sx, t = 2) = 133])[1];

sx : at (sx, ksoln);

[at (sx, t = 0),  at(sx, t = 2)];

sp : at(k1 - 3*t, ksoln);

[at (sp, t = 0),  at (sp, t = 2)];

sy : expand (2*sp);

[ at (sy, t = 0),  at (sy, t = 2)];

### Example 7

Quoting Dowling:
"With an inequality constraint as a terminal endpoint...we first optimize the Hamiltonian subject to
a free endpoint. With a free endpoint we set p(T) = 0, allowing the marginal value of the state 
variable to be taken down to zero. This, in effect, means that as long as the minimum value set by
the constraint is met, the state variable is no longer of any value to us. Our interest in the state
variable does not extend beyond time T."

"Most variables have value, however, and our interest generally extends beyond some narrowly
limited time horizon. In such cases we will not treat the state variable as a free good by permitting
its marginal value to be reduced to zero. We will rather require some minimum value of the state
variable to be preserved for use beyond time T. This means maximizing the Hamiltonian subject 
to a fixed endpoint determined by the minimum value of the constraint. In such cases, p(T) > 0, 
the constraint is binding, and we will not use as much of the state variable as we would if it were 
a free good."

### Prob. 21.7

Maximize J = integrate (8*x - 10*y^2, t, 0, 4)
  subject to:  dx/dt = 24*y, x(0) = 7,  x(4) >= 2000.

This problem was previously solve in Prob. 21.4 as a free endpoint, where we found
  y = 6*p/5,  dx/dt = 144*p/5, dp/dt = -8, p = k1 - 8*t, dx/dt = 144*k1/5 - 1152*t/5,
  x = - 576*t^2/5 + 144*k1*t/5 + k2, and with p(4) = 0, we got sx = -115.2*t^2 + 921.6*t + 7,
at(sx, t = 4) = 1850.2 < 2000. So the free endpoint solution does not provide a large enough
value for x(4). We instead use x(0) = 7 and x(4) = 2000 to determine the constants k1 and k2.

sx : - 576*t^2/5 + 144*k1*t/5 + k2;
ksoln : solve ([ at (sx, t = 0) = 7,  at (sx, t = 4) = 2000] )[1], numer;

sx : at (sx, ksoln), numer;

sp : at (k1 - 8*t, ksoln);

at (sp, t = 4);

sy : float (6*sp/5), expand;

at (sy, t = 4);

### Prob. 21.10  Use of Zindef_float (A, B)

Maximize  J = integrate (4*y - y^2 - x - 2*x^2, t, 0, 1)
  subject to  dx/dt  = x + y,  x(0) = 6.15,  x(1) >= 8.

We have already optimized this function subject to free endpoint conditions in Prob. 21.6 above.
There we found for the value of x(1):

float( 0.98291*exp(sqrt(3))  + 6.0004*exp(-sqrt(3))  - 0.83333);

which is less than 8, so we must evaluate the constants k1 and k2 using fixed endpoint conditions
x(0) = 6.15, x(1) = 8.

Since the path to a solution involved the pair of coupled first order ode's dx/dt = u(x,p) and 
dp/dt = v(x,p) and we used matrix methods, we repeat our solution from scratch.

H : 4*y - y^2 - x - 2*x^2 + p*(x + y);
dy : diff (H, y);

solve (dy, y);

sy : rhs(%[1]);

xd : diff (H, p);

xd : at (xd, y = sy), expand;

pd : - diff (H, x);

dp/dt depends on both x(t) and p(t), as does dx/dt.
So we have a pair of coupled first order odes to solve.
Let Z = cvec ([x, p]) and write the set of first order odes in matrix form:
       dZ/dt = A . Z + B.
The Maxima function Zindef (A, B) returns the indefinite solution as a matrix column vector.

A : matrix ( [1, 1/2], [4, -1] );
B : cvec ( [2, 1]);

eigenvalues (A);

float (%);

Zs : Zindef_float (A, B);

xs : Zs[1,1];
ps : Zs[2,1];

To determine the constants %k1 and %k2 we use x(0) = 6.15 and x(1) = 8.

ksoln : solve ([at(xs, t = 0) = 6.15, at (xs, t = 1) = 8]);

ksoln : float (ksoln);

sx : at (xs, ksoln[1]);
sp : at (ps, ksoln[1]);

To compare with Dowling's solutions we need to replace some pieces of this with floating point
values, which we can do with the Maxima function substpart (zz, expr, n1, n2,...).

part (sp, 3);

sp : substpart (float (-7/3), sp, 3);

[at (sp, t = 0), at (sp, t = 1)];

sx : substpart (float (- 5/6), sx, 3);

Check x(1) value:

at (sx, t = 1);

Finally, we can update sy, which is y(t), and which depends on p(t):

sy : at (sy, p = sp), expand;

[at (sy, t = 0), at (sy, t = 1)];

wxdraw2d (xlabel = "t", key = " state var x(t)", yrange = [-40, 20], key_pos = top_left,
    explicit (sx, t, 0, 1), color = red, key = "control var y(t)", explicit (sy, t, 0, 1),
    color = brown, key = " co-state var p(t)", explicit (sp, t, 0, 1), color = black,
    key = "", line_width = 1, explicit (0, t, 0, 1))$

Check the sufficient conditions for a maximum in J for Prob. 21.10.
1. The constraint function g = x + y is linear in x and y and is always both concave and convex, 
and hence we don't need to worry about the sign of p(t), the co-state variable.  
2. We need the integrand f(x,y) of the functional J to be jointly concave in x and y.
ConcaveTest (f, x, y) finds f(x,y) is strictly concave, and the maximum is therefore a 
global maximum.

ConcaveTest (4*y - y^2 - x - 2*x^2, x, y)$

## The Current-Valued Hamiltonian Hc

### Modified Maximum Principle with Hc

We seek to maximize J = integrate (exp(-r*t)*f(x(t),y(t),t), t, 0, T),  with r = constant,
     subject to: dx/dt = g(x(t), y(t), t), x(0) = x0, x(T) = free.
Using p(t) as the co-state variable again, the Hamiltonian is
   H = exp (-r*t)*f(x,y,t) + p(t)*g(x,y,t).

Let μ = p*exp (r*t), then p = μ*exp(-r*t) and ∂μ/∂p = exp(r*t).

Now define the current valued Hamiltonian Hc as:
Hc = H*exp(r*t) = f + μ*g, and then, of course, H = Hc*exp(-r*t).

Then dp/dt = -∂H/∂x = - ∂/∂x ( Hc*exp(-r*t) ) = - exp(-r*t)*∂Hc/∂x.

But using p = μ*exp(-r*t), we also have
   dp/dt = -r*μ*exp(-r*t) + exp(-r*t)*dμ/dt = - exp(-r*t)*∂Hc/∂x.

Cancelling the common factor exp (-r*t), we get:

            dμ/dt = r*μ - ∂Hc/∂x.

Moreover, the endpoint condition p(T) = 0 implies μ(T)*exp(-r*T) = 0.

Assuming an interior solution, the requirement
   ∂H/∂y = 0 implies  ∂( Hc*exp(-r*t) )/∂y = 0, or, because exp(-r*t) > 0,
     
           ∂Hc/∂y = 0

  We next work on dx/dt.
   dx/dt = ∂H/∂p  =  ∂ (Hc*exp(-r*t)) /∂p = exp(-r*t)*∂Hc/∂p = exp(-r*t)*∂μ/∂p*∂Hc/∂μ,
    using the chain rule of differentiation. But ∂μ/∂p = exp(r*t), so we finally have

         dx/dt = ∂Hc/∂μ.
      

Summarizing the current-valued Hamiltonian approach to the first order equations of motion:

  1. For an interior solution, require ∂Hc/∂y = 0.
  2. Require the pair of equations:
           dx/dt = ∂Hc/∂μ,   dμ/dt = r*μ - ∂Hc/∂x.
  3. For x(T) = free boundary condition, use   μ(T)*exp(-r*T) = 0.

### Example 8  Use of Zindef_float (A, B)

Maximize J = integrate (exp (-0.02*t)*(x - 3*x^2 - 2*y^2), t, 0, 2)
  subject to:  dx/dt = y - 0.5*x,  x(0) = 93.91, x(2) = free.

A. Set up the current valued Hamiltonian Hc

Hc : x - 3*x^2 - 2*y^2  + μ*(y - 0.5*x);

B. Assuming an interior solution, apply the modified Maximum principle.
1. ∂Hc/∂y = 0.

dy : diff (Hc, y);

solve (dy, y);

sy : rhs (%[1]);

2. dx/dt = ∂Hc/∂μ and dμ/dt = r*μ - ∂Hc/∂x:

xd : diff (Hc, μ);

xd : at (xd, y = sy), numer;

μd : 0.02*μ - diff (Hc, x);

dx/dt and dμ/dt each depend on x and μ, so we have a pair of first order ode's which we solve 
using matrix methods, as in Prob. 21.6 and Prob. 21.10 above. Our solution matrix column vector
Zs has the form:

cvec([ μ, x ]);

A : matrix ([ 0.52, 6], [0.25, -0.5]);
B : cvec ([ -1, 0]);

eigenvalues (A);

float(%);

cvec ([ dμ/dt, dx/dt]) = A . cvec ([μ, x]) + B;

Zs : Zindef_float (A, B);

μs : Zs[1,1];
xs : Zs[2,1];

With T = 2, we require exp(-r*T)*μ(T) = 0 as part of the next step.

ksoln : float (solve ([at (xs, t = 0) = 93.91, exp(-0.02*2)*at (μs, t = 2) = 0]));

ksoln : ksoln[1];

sμ : at(μs, ksoln);
sx : at (xs, ksoln);

[at (sx, t = 0), at (sx, t = 2)];

[at (sμ, t = 0), at (sμ, t = 2)];

sy : at (sy, μ = sμ), expand;

[at( sy, t = 0), at (sy, t = 2)];

wxdraw2d (xlabel = "t", key = " state var x(t)", yrange = [-350, 150], key_pos = bottom_right,
    explicit (sx, t, 0, 2), color = red, key = "control var y(t)", explicit (sy, t, 0, 2),
    color = brown, key = " co-state var mu(t)", explicit (sμ, t, 0, 2), color = black,
    key = "", line_width = 1, explicit (0, t, 0, 2))$

Applying the sufficiency test we look at the matrix

matrix([Fxx, Fxy], [Fyx, Fyy]);

with F = exp(-r*t)*f(x,y) and with r a constant, so Fxy = exp(-r*t)*fxy, etc, the sign of Fxx will
be the same as the sign of fxx, and the determinant of the Hessian matrix of F will be equal to
exp(-2*r*t) times the determinant of the Hessian matrix of f, with both determinants having the
same sign.

Hence we can use ConcaveTest(f(x,y), x, y) and ignore the exponential factor.

ConcaveTest (x - 3*x^2 - 2*y^2, x, y)$

With g = x + y linear in both x and y, we therefore have found a global maximum.

### Prob. 21.11  Use of Zindef_float (A,B)

Maximize  J = integrate (exp(- 0.05*t)*( x*y - x^2 - y^2), t, 0, 3)
  subject to:    dx/dt = x + y,  x(0) = 134.35,  x(3) = free.

A. Set up the current valued Hamiltonian Hc

Hc : x*y - x^2 - y^2  + μ*(x + y);

B. Assuming an interior solution, apply the modified Maximum principle.
1. ∂Hc/∂y = 0.

dy : diff (Hc, y);

solve (dy, y);

sy : rhs (%[1]);

sy : float (expand (sy));

2. dx/dt = ∂Hc/∂μ and dμ/dt = r*μ - ∂Hc/∂x:

xd : diff (Hc, μ);

xd : at (xd, y = sy);

μd : 0.05*μ - diff (Hc, x);

μd : at (μd, y = sy);

A : matrix ([ - 1.45, 1.5], [0.5, 1.5]);
B : cvec ([ 0, 0]);

eigenvalues (A);

float(%);

cvec ([ dμ/dt, dx/dt]) = A . cvec ([μ, x]) + B;

Zs : Zindef_float (A, B);

μs : Zs[1,1];
xs : Zs[2,1];

ksoln : float (solve ([at (xs, t = 0) = 134.35, exp(-0.05*3)*at (μs, t = 3) = 0]));

ksoln : ksoln[1];

sμ : at(μs, ksoln);
sx : at (xs, ksoln);

[at (sx, t = 0), at (sx, t = 3)];

[at (sμ, t = 0), at (sμ, t = 3)];

sy : at (sy, [μ = sμ, x = sx]), expand;

[at( sy, t = 0), at (sy, t = 3)];

wxdraw2d (xlabel = "t", key = " state var x(t)", yrange = [-900, 150], key_pos = bottom_right,
    explicit (sx, t, 0, 3), color = red, key = "control var y(t)", explicit (sy, t, 0, 3),
    color = brown, key = " co-state var mu(t)", explicit (sμ, t, 0, 3), color = black,
    key = "", line_width = 1, explicit (0, t, 0, 3))$

ConcaveTest (x*y - x^2 - y^2, x, y)$

With g = x + y linear in both x and y, we therefore have found a global maximum.