APPROACH TO CONSTRAINT SET: SIMPLE OPTIMAL CONTROL PROBLEM
				
by

Paul Scott

Introduction

This notebook performs the numerical calculations proposed in a paper with same title.  References below refere to this paper.  This is a shooting method in which one varies the inital values of the costate until the constraint set is reached.

Activate Part One initially.  Then part two can be repeated with new initial costate values until the constraint set is reached.

Part One

Declaring Variables

Number of y state variables (called n).

n  =  2

Declaration of y state variables.

yy  =  Table[ (Subscript[y, i]) ,  { i ,  1 ,  n } ]

Declaration of y- costate variables (called λ).

λλ  =  
      Table[  (Subscript[λ, i])  ,  { i , 1 , n }  ]

Number of z state variables (called m).

m  =  2

Declaration of z-state variables.

zz  =  Table[  (Subscript[z, i])  ,  { i , 1 , m }  ]

Declaration of z-costate variables (called ξ).

ξξ  =  Table[  (Subscript[ξ, i])  ,  { i , 1 , m}  ]

Declaration of control variables (called u).  See paper (1.0.3) and (1.0.4).

uu  =  Table[  (Subscript[u, i])  ,  { i , 1 , m }  ]

Defining the differential equation and constraint

Definition of rhs of ODE   F.  See paper (1.0.1)

Each component of F is defined.

(Subscript[f, 1])   =   -(Subscript[y, 1])   +   (Subscript[z, 1])

(Subscript[f, 2])   =   -(Subscript[y, 2])   +   (Subscript[z, 2])

The components are listed in vector F.

F  =  Table[  (Subscript[f, i])  ,  { i , 1 , n }  ]

Definition of constraint function G.  See paper (1.0.2)

Each component of g is defined.

(Subscript[g, 1])  =   
      (Subscript[y, 1])^2   +   (Subscript[y, 2])^2   -  ((Subscript[z, 1])^2   +   (Subscript[z, 2])^2 )

(Subscript[g, 2])  =    (Subscript[y, 1])   -   (Subscript[z, 1])

The components are listed in a vector G.

G  =  Table[  (Subscript[g, i])  ,  { i , 1 , m }  ]

Definition of Hamiltonian.  See paper (3.0.7)

H  =   
      G . G   +   λλ . F   +   
        ξξ . uu

The following are the partial derivatives that occur in the costate differential equations.  This is to show them.

(∂_((Subscript[y, 1]))) H

(∂_((Subscript[y, 2]))) H

(∂_((Subscript[z, 1]))) H

(∂_((Subscript[z, 2]))) H

The feedback control law.  See paper (3.0.10)

Do[ uu[[i]]   =    -Sign[  ξξ[[i]]  ] , 
       {i, 1, m} ]

uu

Define variables for numerical calculations

These variables are functions of t for substitution into the differential equations.

YY  =  Table[ (Subscript[Y, i])[t] ,  { i ,  1 ,  n } ]

ΓΓ  =  
      Table[  (Subscript[Γ, i])[t]  ,  { i , 1 , n }  ]

ZZ  =  Table[ (Subscript[Z, i])[t]  ,  { i ,  1 , m }  ]

ΔΔ  =  
      Table[  (Subscript[Δ, i])[t]  ,  { i , 1 , m }  ]

UU  =  Table[  (Subscript[U, i])[t]  ,  { i , 1 , m }  ]

The following list called "sset" is the list of unknown functions to be found by solving the differential equations set up below.

sset   =   
      Union[ YY ,  ZZ ,  ΓΓ ,  
        ΔΔ ]

Substitution lists to insert numerical variables

suby  =  
      Table[ yy[[i]]  ->  YY[[i]]  ,  {i, 1, n} ]

subλ  =  
      Table[ λλ[[i]]  ->  
          ΓΓ[[i]]  ,  {i, 1, n} ]

subz  =  
      Table[ zz[[i]]  ->  ZZ[[i]]  ,  {i, 1, m} ]

subξ  =  
      Table[ ξξ[[i]]  ->  
          ΔΔ[[i]]  ,  {i, 1, m} ]

The following list called "sub" is the union of all the above substitution lists.

sub  =  
      Union[  suby ,  subz ,  subλ ,  subξ   ]

For graphic analysis latter

The variable "GG" is the constraint computed along the trajectories, and should go to zero as the trasjectory approaches the constraint set.

GG   =   G  /.  sub

Substitute numerical variables into the feedback law

This substitution replaces the variable ξξ in the feedback by the corresponding function ΔΔ[t] for numerical solution.

Do[ 
      UU[[i]]   =    
        uu[[i]]  /.  
          ξξ[[i]]  ->  
            ΔΔ[[i]]  , 
                  {i, 1, m}
                                                                                                        ]

UU

Set initial values for y-state and z-state.

initialy  =  {   2 ,  0 }

initialz  =  { 0 ,  1 }

Set up state equations for y and z

Below ystate is the set of differential equations for y, see paper (1.0.1)  and initial conditions.  The command makes a table with each entry containing tow parts.  The first part is the DE in which the Derivative of y is set equal to F and the numerical variables are substituted in.  The second part sets the initial value of y.

\(ystate\  = \ 
      Table[\ {\ 
          D[\ YY[\([i]\)]\ \ , \ t\ ]\  == \ F[\([i]\)]\ \  /. \ sub\ \ , \ 

          				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          \((\ \ YY[\([i]\)]\  /. \ t\  -> \ 0\ \ )\)\  == \ 
            initialy[\([i]\)]\ }\ , 

        				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\ i\ , \ 
          1\ , \ n\ }
          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\)

Below zstate is the set of differential equations for z, see paper (1.0.4)  and initial conditions.

\(zstate\  = \ 
      Table[\ {\ 
          D[\ ZZ[\([i]\)]\ \ , \ t\ ]\  == \ 
            UU[\([i]\)]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ , \ 

          				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          \((\ \ ZZ[\([i]\)]\  /. \ t\  -> \ 0\ \ )\)\  == \ 
            initialz[\([i]\)]\ }\ , 

        				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\ i\ , \ 1
          \ , \ m\ }
          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\)

Part Two

Set initial costate values

Adjust initial costate values until constraint set is approached.

initialλ  =  { -200 ,  -300 }

initialξ  =  { -400 ,  -200 }

"tend" is the end time for solving the trajectories.  You might want to change this sometimes.

tend   =   1.15

tinterval   =   { t, 0, tend }

Set the costate equations

The following sets the λ-costate equations and initial conditions.  See paper (3.0.8)

\(λstate\  = \ 
      Table[\ {\ 
          D[\ ΓΓ[\([i]\)]\ \ , \ t\ ]\  == \ 
              \(-\ ∂\_\(y\_i\)\ H\)\ \  /. \ sub\ \ , \ 

          				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          \((\ \ ΓΓ[\([i]\)]\  /. \ 
                t\  -> \ 0\ \ )\)\  == \ initialλ[\([i]\)]\ }\ , 

        				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\ i\ , \ 
          1\ , \ n\ }
          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\)

The following sets the ξ-costate equations and initial conditions.  See paper (3.0.8)

\(ξstate\  = \ 
      Table[\ {\ 
          D[\ ΔΔ[\([i]\)]\ \ , \ t\ ]\  == \ 
              \(-\ ∂\_\(z\_i\)\ H\)\ \  /. \ sub\ \ , \ 

          				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          \((\ \ ΔΔ[\([i]\)]\  /. \ 
                t\  -> \ 0\ \ )\)\  == \ initialξ[\([i]\)]\ }\ , 

        				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\ i\ , \ 
          1\ , \ m\ }
          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\)

Putting all the DE's and initial conditions in one list called "state".

state  =  
      Flatten[Union[ ystate ,  zstate ,  λstate ,  ξstate
           ]]

Solving equations

"NDSolve" numerically solves the differential equations and initial conditions in list "state" for all the unknown functions in list "sset", over the time interval "tinterval".  The solution is kept in a list called "sol".

sol   =   NDSolve[  state ,  sset ,  tinterval ]

Plots of the components of the constraint GG

Below the plot of the first component of GG.

Plot[ Evaluate[ GG[[1]]  /.  sol ] ,  { t, 0, tend } ,  
      AxesOrigin  ->  {0, 0} ]

Below the plot of the second component of GG.

Plot[ Evaluate[ GG[[2]]  /.  sol ] ,  { t, 0, tend } ,  
      AxesOrigin  ->  {0, 0} ]

Click the hypertext "il" to go back to inital costate values and try again.

il

Numerical values at a point near the constraint set

Set the following variable to the value of t when the graphs are both close to zero

tset   =   1.07

The following command shows the list of unknown functions

sset

The following command gives the values of all the functions in list "sset" at time "tset".

( sset  /.  sol )  /.  t  ->  tset