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
    # 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

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