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; 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]; # 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 == v_0 h == h_0 m == 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 (
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.
- OS Windows
- Julia Version 1.5.3
- JuMP v0.21.4
- Ipopt v0.6.5