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
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 T
needs 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);
Change Solver Options
Solver options can be easily changed, for example, consider changing the print level in Ipopt as
defineSolver!(n,((:name=>:Ipopt),(:print_level=>8)))
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:
:integrateScheme | Description |
---|---|
:lgrExplicit | default scheme; implementation derivative constraints in hp-pseudospecral method |
:lgrImplicit | implementation of integral constraints in hp-pseudospecral method |
:bkwEuler | approximate using backward euler method |
:trapezoidal | approximate using trapezoidal method |
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)