Hey everyone,
I was playing around with nonlinear functions for optimal control problems utilizing JuMP with IPOPT as the solver. I realized convergence problems as soon as I use the splatting operator for passing vectors to the objective function as suggested in the JuMP documentation for Nonlinear Optimization. I noticed the difference at an own problem but could reproduce it with the rocket control problem from JuMPTutorials.jl. Here is the slightly modified version to show the differences:
using JuMP, Ipopt
function optControlProb(objFcnSelector)
# Create JuMP model, using Ipopt as the solver
rocket = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 5, "mu_init" => 1/100))
# Constants
# Note that all parameters in the model have been normalized
# to be dimensionless. See the COPS3 paper for more info.
h_0 = 1 # Initial height
v_0 = 0 # Initial velocity
m_0 = 1 # Initial mass
g_0 = 1 # Gravity at the surface
T_c = 3.5 # Used for thrust
h_c = 500 # Used for drag
v_c = 620 # Used for drag
m_c = 0.6 # Fraction of initial mass left at end
c = 0.5 * sqrt(g_0 * h_0) # Thrust-to-fuel mass
m_f = m_c * m_0 # Final mass
D_c = 0.5 * v_c * m_0 / g_0 # Drag scaling
T_max = T_c * g_0 * m_0 # Maximum thrust
n = 799 # Time steps
@variables(rocket, begin
Δt ≥ 0, (start = 1/(n+1)) # Time step
# State variables
v[0:n] ≥ 0 # Velocity
h[0:n] ≥ h_0 # Height
m_f ≤ m[0:n] ≤ m_0 # Mass
# Control
0 ≤ T[0:n] ≤ T_max # Thrust
end)
if objFcnSelector == 1
# use splatting to pass a vector to the objective function for the user defined function
function objFcnRocket(x...)
result = x[1];
return result;
end
register(rocket, :objFcnRocket, 2, objFcnRocket; autodiff = true)
# define a vector with two entries, where only the first one is used within the objective function
temp = [h[n];h[0]];
# Objective: maximize altitude at end of time of flight
@NLobjective(rocket, Max, objFcnRocket(temp...))
elseif objFcnSelector == 2
# use the objective from the tutorial
@objective(rocket, Max, h[n])
elseif objFcnSelector == 3
# use the objective from the tutorial, however define it as a nonlinear one
@NLobjective(rocket, Max, h[n])
elseif objFcnSelector == 4
# define a user defined function for the objective which does not use splatting
function objFcnRocketScalar(x)
result = x;
return result;
end
register(rocket, :objFcnRocketScalar, 1, objFcnRocketScalar; autodiff = true)
# use the scalar user defined function for the optimization
@NLobjective(rocket, Max, objFcnRocketScalar(h[n]) )
end
# Initial conditions
@constraints(rocket, begin
v[0] == v_0
h[0] == h_0
m[0] == m_0
m[n] == m_f
end)
# Forces
# Drag(h,v) = Dc v^2 exp( -hc * (h - h0) / h0 )
@NLexpression(rocket, drag[j = 0:n], D_c * (v[j]^2) * exp(-h_c * (h[j] - h_0) / h_0))
# Grav(h) = go * (h0 / h)^2
@NLexpression(rocket, grav[j = 0:n], g_0 * (h_0 / h[j])^2)
# Time of flight
@NLexpression(rocket, t_f, Δt * (n+1))
# Dynamics
for j in 1:n
# h' = v
# Rectangular integration
# @NLconstraint(rocket, h[j] == h[j - 1] + Δt * v[j - 1])
# Trapezoidal integration
@NLconstraint(rocket,
h[j] == h[j - 1] + 0.5 * Δt * (v[j] + v[j - 1]))
# v' = (T-D(h,v))/m - g(h)
# Rectangular integration
# @NLconstraint(rocket, v[j] == v[j - 1] + Δt *(
# (T[j - 1] - drag[j - 1]) / m[j - 1] - grav[j - 1]))
# Trapezoidal integration
@NLconstraint(rocket,
v[j] == v[j-1] + 0.5 * Δt * (
(T[j] - drag[j] - m[j] * grav[j]) / m[j] +
(T[j - 1] - drag[j - 1] - m[j - 1] * grav[j - 1]) / m[j - 1]))
# m' = -T/c
# Rectangular integration
# @NLconstraint(rocket, m[j] == m[j - 1] - Δt * T[j - 1] / c)
# Trapezoidal integration
@NLconstraint(rocket,
m[j] == m[j - 1] - 0.5 * Δt * (T[j] + T[j-1]) / c)
end
# Solve for the control and state
println("-------------------------", objFcnSelector, "-------------------------")
println("Solving...")
status = optimize!(rocket)
# Display results
# println("Solver status: ", status)
println("Max height: ", objective_value(rocket))
println("Steptime: ", value.(Δt))
end
optControlProb(1);
optControlProb(2);
optControlProb(3);
optControlProb(4);
The first call of to solve the optimal control problem (optControlProb(1);) uses the splatting operator in a user defined function for the NLobjective, but only the first value from the vector is taken. A restoration error occurs. The second call (optControlProb(2);) just uses the linear objective function from the example itself. The third call (optControlProb(3);) uses NLobjective instead of objective (even though it is linear, just for checking the results). The fourth call (optControlProb(4);) also uses a user defined function for the NLobjective, however it does not applies the splatting operator.
I am curious about the different results if the splatting operator is utilized. The version with the splatting operator has lg(mu) = 0 for the barrier parameter, whereas the other ones have lg(mu) = -1 (which is supposed to be the default one) for iter 0 (taken from the IPOPT prints). Changing mu_init does not seem to affect the version with the splatting operator but does for the other ones.
To be honest, I lack the knowledge to find the cause for this behavior and I hope that someone could help me out here.
Kind regards
Additional information:
- OS Windows
- Julia Version 1.5.3
- JuMP v0.21.4
- Ipopt v0.6.5