pulsipher
commented
3 years ago Thank you for reaching out and trying InfiniteOpt.jl! The behavior described above can be explained via how the time supports are managed in building the model, and how this behavior differs between using trapezoid rule and Gauss quadrature. In the above example this leads to certain time points not being penalized in the objective which allows large thrust values to occur at those instances leading to a result one would likely not expect. Let me illustrate this with the cases below, notice the added comments on how supports are added and used.
Default Model (using trapezoid rule):
using InfiniteOpt, JuMP, Ipopt
# Waypoints
xw = [1 4 6 1; 1 3 0 1] # positions
tw = [0, 25, 50, 60] # times
# Initialize the model
m = InfiniteModel(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
# Set the parameters and variables
@infinite_parameter(m, t in [0, 60], num_supports = 61) # adds 61 supports at [0, 1, ..., 60]
@infinite_variable(m, x[1:2](t), start = 1)
@infinite_variable(m, v[1:2](t), start = 0)
@infinite_variable(m, u[1:2](t), start = 0)
# Specify the objective
@objective(m, Min, integral(u[1]^2 + u[2]^2, t)) # uses trapezoid rule which will incorporate all supports added before `optimize!`
# Set the constraints (this will incorporate all supports added before `optimize!` is invoked)
@BDconstraint(m, initial_velocity[i = 1:2](t == 0), v[i] == 0) # adds a support 0 to t if not already there (which it is in this case)
@constraint(m, [i = 1:2], deriv(x[i], t) == v[i])
@constraint(m, [i = 1:2], deriv(v[i], t) == u[i])
@BDconstraint(m, [i = 1:2, j = eachindex(tw)](t == tw[j]), x[i] == xw[i, j])
# Optimize the model (will make fresh transcription model formatted as a JuMP model and then solved)
optimize!(m)
# Plot results (not shown here for conciseness)
Thus, invoking this and plotting obtains the trajectory we expect.
Adding Supports and Using Quadrature (leads to many ignored terms in the objective):
using InfiniteOpt, JuMP, Ipopt
# Waypoints
xw = [1 4 6 1; 1 3 0 1] # positions
tw = [0, 25, 50, 60] # times
# Initialize the model
m = InfiniteModel(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
# Set the parameters and variables
@infinite_parameter(m, t in [0, 60], num_supports = 61) # adds 61 supports at [0, 1, ..., 60]
@infinite_variable(m, x[1:2](t), start = 1)
@infinite_variable(m, v[1:2](t), start = 0)
@infinite_variable(m, u[1:2](t), start = 0)
# Specify the objective using Gauss quadrature
# Note this will only incorporate the support points it generates and will ignore the 61 supports added above
@objective(m, Min, integral(u[1]^2 + u[2]^2, t, eval_method = Quadrature)) # adds 10 Gauss quadrature nodes to the supports
println(num_supports(t)) # prints 71 (now our model has the 61 supports we added + the 10 Gauss quadrature ones)
# Set the constraints (this will incorporate all 71 supports)
@BDconstraint(m, initial_velocity[i = 1:2](t == 0), v[i] == 0) # would add a support 0 to t if not already there, but this has already been added
@constraint(m, [i = 1:2], deriv(x[i], t) == v[i])
@constraint(m, [i = 1:2], deriv(v[i], t) == u[i])
@BDconstraint(m, [i = 1:2, j = eachindex(tw)](t == tw[j]), x[i] == xw[i, j])
# Optimize the model (will make fresh transcription model formatted as a JuMP model and then solved)
optimize!(m)
# Plot results (not shown here for conciseness)
Now we obtain a very different result, however it does satisfy all the constraints (including the derivatives). This v shape occurs because we are only penalized to use thrust on the 10 Gauss quadrature points [0.782804144485, 4.04809899933, ..., 59.2171958555] thus the thrust given at these instances is 0 and any amount is given at the other time points without penalty allowing the hovercraft to hit all the needed waypoints. Visually this is shown in the plot below. Notice that all the Gauss time points have 0 thrust and large amounts of thrust are used liberally at the other points.
Only Use Quadrature Supports:
using InfiniteOpt, JuMP, Ipopt
# Define function to find the closest time point
function findnearest(a, x)
if isempty(a)
return
end
a = sort(a)
idx = findfirst(v -> v >= x, a)
if idx isa Nothing
return last(a)
elseif idx == 1
return a[idx]
elseif abs(a[idx] - x) <= abs(a[idx - 1] - x)
return a[idx]
else
return a[idx - 1]
end
end
# Waypoints
xw = [1 4 6 1; 1 3 0 1] # positions
tw = [0, 25, 50, 60] # times (we will find the closest quadrature points further down)
# Initialize the model
m = InfiniteModel(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
# Set the parameters and variables
@infinite_parameter(m, t in [0, 60]) # don't add any supports yet
@infinite_variable(m, x[1:2](t), start = 1)
@infinite_variable(m, v[1:2](t), start = 0)
@infinite_variable(m, u[1:2](t), start = 0)
# Specify the objective using Gauss quadrature and specify we want 61 supports
# Note this will only incorporate the support points it generates and will ignore any others we add elsewhere
@objective(m, Min, integral(u[1]^2 + u[2]^2, t, eval_method = Quadrature, num_supports = 61)) # adds 61 Gauss quadrature nodes to the supports
# Set the constraints (this will incorporate all supports)
@BDconstraint(m, initial_velocity[i = 1:2](t == 0), v[i] == 0) # this adds a support of 0 to t (now we have 62 supports)
@constraint(m, [i = 1:2], deriv(x[i], t) == v[i])
@constraint(m, [i = 1:2], deriv(v[i], t) == u[i])
twc = [findnearest(supports(t), w) for w in tw] # match the waypoint times to the closest support
@BDconstraint(m, [i = 1:2, j = eachindex(twc)](t == twc[j]), x[i] == xw[i, j])
# Optimize the model (will make fresh transcription model formatted as a JuMP model and then solved)
optimize!(m)
# Plot results (not shown here for conciseness)
This works as we would intuitively expect since all the supports (except the first) are included in the objective integral. Excluding the first is ok since the initial condition does not correspond to a degree of freedom in the transcripted problem.
Pulling this all together, the takeaway is that quadrature can be used to evaluate integrals, but this should be done with caution knowing that any supports added to the model elsewhere will not be included in the integral approximation since quadrature exacts that the included supports be at particular points that are orthogonal to the polynomial basis. For this reason, we employ trapezoid rule by default since it can and will include all the supports that are added to the model.
I hope this clarifies your issues you were running into. Suggestions and/or contributions to make this behavior more clear in our documentation would be appreciated.
mzagorowska
In brief, if I want to use any other quadrature than trapezoid in the hovercraft example, the constraints with derivatives are not considered in the model. The minimal working example contains 4 cases (to comment/uncomment):
The code with all the cases:
Desktop:
I am attaching a document with all the plots I got ( InfiniteOpt_quadratures.pdf). I am not getting any errors when running the cases, but let me know if I should provide more information.