# Nonlinear Objective Function Splatting

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
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 (`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

• OS Windows
• Julia Version 1.5.3
• JuMP v0.21.4
• Ipopt v0.6.5

JuMP doesn’t support second-derivative information for multivariate user-defined functions (docs: https://jump.dev/JuMP.jl/stable/nlp/#User-defined-Functions-1).

You can see this by comparing the number of non-zero elements in the Lagrangian Hessian of approach 1:

``````-------------------------1-------------------------
This is Ipopt version 3.13.2, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:    15185
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        0
``````

and approach 4:

``````-------------------------4-------------------------
This is Ipopt version 3.13.2, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:    15185
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:    45544
``````
2 Likes

Thank you very much!

I misunderstood the part with the support for second-derivative information. I thought that only user provided Hessian calculation routines for multivariate functions are not supported and `autodiff` still works.Thanks for clarifying this and also showing it for the example!

1 Like

There’s an issue to track: https://github.com/jump-dev/JuMP.jl/issues/1198. But I wouldn’t expect any progress anytime soon.

I’ve added a link to this post which will be useful if we ever get around to adding support.

1 Like