Quick Ex#2: Moon Lander

Given:

A space-ship landing on the moon

Dynamic Constraints

Boundary Conditions

Control Limits

Find:

The track that minimizes time

Solution:

In this problem, we put the bounds directly into define(). Also, now we have constant limits on the control variables and those can be added as shown below This problem can be found here.

Packages that will be used

using NLOptControl, PrettyPlots

Define the Problem:

In this problem, we put the bounds directly into define(). Also, now we have constant limits on the control variables and those can be added as shown below

n = define(numStates=2,numControls=1,X0=[10.,-2],XF=[0.,0.],CL=[0.],CU=[3.])

State and Control Names

The state and control variables are by default, $:x1,:x2,..$ and $:u1,:u2,..$, but they can be changed with the following commands:

states!(n,[:h,:v];descriptions=["h(t)","v(t)"])
controls!(n,[:T];descriptions=["T(t)"])

Next, now that the problem is configured, all of the state and control variables are stored in JuMP Arrays, n.r.ocp.x[:,:] and n.r.u[:,:], respectively. For instance;

typeof(n.r.ocp.x)

Array{Array{Any,2},1}

Differential Equations

dx=[:(v[j]),:(T[j]-1.625)]
dynamics!(n,dx)

Configure the Problem:

configure!(n;(:finalTimeDV=>true));

Integral Terms in the Cost Function

integrate!() is used to make terms that can be added to the cost function that need to be integrated. When calling this function an expression must be passed:

In this example the first control variable Tneeds to be integrated, it must be passed in an expression :() with the index [j]. To do this, integrate!() can be used as:

obj = integrate!(n,:(T[j]))
# Now this term can be added as the objective function and the problem can be solved
@NLobjective(n.ocp.mdl, Min, obj);

Optimize

optimize!(n)

Post Process

allPlots(n)

Other Dynamic Constraint Methods

Currently there are three different methods to ensure that the dyanamic constraints are satisfied and they are set when configure!() is called using the :integrationScheme key. They are listed below:

The later two are time-marching methods and default number of points is 100, but that can be changed by setting N. So, the above problem can be solved using one of the time-marching schemes as:

n = define(numStates=2,numControls=1,X0=[10.,-2],XF=[0.,0.],CL=[0.],CU=[3.])
states!(n,[:h,:v];descriptions=["h(t)","v(t)"])
controls!(n,[:T];descriptions=["T(t)"])
dx=[:(v[j]),:(T[j]-1.625)]
dynamics!(n,dx)
configure!(n,N=200;(:integrationScheme=>:trapezoidal),(:finalTimeDV=>true))
obj = integrate!(n,:(T[j]))
@NLobjective(n.ocp.mdl, Min, obj)
optimize!(n)
allPlots(n)

Constraints

Often when building a model and using it to solve an optimal control problem, their are issues associated with infeasibility. NLOptControl has functionality to help deal with these issues. For instance, the dual infeasibility values can stored and quickly viewed. They are stored in a DataFrame which can be referenced with n.r.ocp.constraint.value as:

n.r.ocp.constraint.value

0-element Array{Any,1}

It is empty, because by default this data is not calculated and stored. This option can be turned on by modifying the settings for the problem:

n.s.ocp.evalConstraints

false

n.s.ocp.evalConstraints = true
optimize!(n)
n.r.ocp.constraint.value

3-element Array{Any,1}:
 2×2 DataFrames.DataFrame
│ Row │ step │ x0_con    │
├─────┼──────┼───────────┤
│ 1   │ 1    │ 0.396941  │
│ 2   │ 2    │ -0.492598 │
 2×2 DataFrames.DataFrame
│ Row │ step │ xf_con    │
├─────┼──────┼───────────┤
│ 1   │ 1    │ -0.396941 │
│ 2   │ 2    │ 2.18394   │
 200×3 DataFrames.DataFrame
│ Row │ step │ h        │ v         │
├─────┼──────┼──────────┼───────────┤
│ 1   │ 1    │ 0.396941 │ -0.496826 │
│ 2   │ 2    │ 0.396941 │ -0.505283 │
│ 3   │ 3    │ 0.396941 │ -0.513739 │
│ 4   │ 4    │ 0.396941 │ -0.522196 │
│ 5   │ 5    │ 0.396941 │ -0.530653 │
│ 6   │ 6    │ 0.396941 │ -0.53911  │
│ 7   │ 7    │ 0.396941 │ -0.547566 │
│ 8   │ 8    │ 0.396941 │ -0.556023 │
⋮
│ 192 │ 192  │ 0.396941 │ -2.11206  │
│ 193 │ 193  │ 0.396941 │ -2.12051  │
│ 194 │ 194  │ 0.396941 │ -2.12897  │
│ 195 │ 195  │ 0.396941 │ -2.13743  │
│ 196 │ 196  │ 0.396941 │ -2.14588  │
│ 197 │ 197  │ 0.396941 │ -2.15434  │
│ 198 │ 198  │ 0.396941 │ -2.1628   │
│ 199 │ 199  │ 0.396941 │ -2.17125  │
│ 200 │ 200  │ 0.396941 │ -2.17971  │

evalMaxDualInf(n)

2.1839379581754303

The last function called, searches through all of the dual infeasibilities to find the largest value. As, this problem is, it is feasible and optimal. But if there was an issue, often looking for high values in these DataFrame structures is the quickest way to figure out the constraints that are giving the solver trouble.

Tolerances

If there was an example where the dual infeasibility value for one or more of the variables was very high, but the actual constraint is only being violated slightly (by some reasonable amount) then the tolerances on the initial and terminal states can be adjusted. This will also improve the solve time, so it is good practice to set these to reasonable values. For instance, in the Moon Lander example, we can set them as:

n = define(numStates=2,numControls=1,X0=[10.,-2],XF=[0.,0.],CL=[0.],CU=[3.])
states!(n,[:h,:v];descriptions=["h(t)","v(t)"])
controls!(n,[:T];descriptions=["T(t)"])
dx = [:(v[j]),:(T[j]-1.625)]
dynamics!(n,dx)
XF_tol = [2.0,0.5]
X0_tol = [0.05,0.05]
defineTolerances!(n;X0_tol=X0_tol,XF_tol=XF_tol)
configure!(n,N=50;(:integrationScheme=>:bkwEuler),(:finalTimeDV=>true))
obj = integrate!(n,:(T[j]))
@NLobjective(n.ocp.mdl, Min, obj)
optimize!(n)
allPlots(n)