# Nonlinear Objective Function Splatting

**URL:** <https://discourse.julialang.org/t/nonlinear-objective-function-splatting/51251>\
**Category:** Optimization (Mathematical)\
**Tags:** jump\
**Created:** [December 4, 2020, 1:25pm UTC](https://discourse.julialang.org/t/nonlinear-objective-function-splatting/51251 "2020-12-04T13:25:04Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Flesko](https://avatars.discourse-cdn.com/v4/letter/f/ba9def/32.png) [@Flesko](https://discourse.julialang.org/u/Flesko)\
**Post date:** [December 4, 2020, 1:25pm UTC](https://discourse.julialang.org/t/nonlinear-objective-function-splatting/51251/1 "2020-12-04T13:25:04Z")

</div>

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](https://jump.dev/JuMP.jl/stable/nlp/#). I noticed the difference at an own problem but could reproduce it with the [rocket control problem](https://nbviewer.jupyter.org/github/jump-dev/JuMPTutorials.jl/blob/master/notebook/modelling/rocket_control.ipynb) from JuMPTutorials.jl. Here is the slightly modified version to show the differences:

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

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [December 5, 2020, 3:00am UTC](https://discourse.julialang.org/t/nonlinear-objective-function-splatting/51251/2 "2020-12-05T03:00:30Z")

</div>

JuMP doesn’t support second-derivative information for multivariate user-defined functions (docs: [https://jump.dev/JuMP.jl/stable/nlp/#User-defined-Functions-1](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:

```julia
-------------------------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:

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

```

---

<div class="post-metadata">

**Author:** ![Flesko](https://avatars.discourse-cdn.com/v4/letter/f/ba9def/32.png) [@Flesko](https://discourse.julialang.org/u/Flesko)\
**Post date:** [December 7, 2020, 7:17am UTC](https://discourse.julialang.org/t/nonlinear-objective-function-splatting/51251/3 "2020-12-07T07:17:01Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [December 7, 2020, 8:30am UTC](https://discourse.julialang.org/t/nonlinear-objective-function-splatting/51251/4 "2020-12-07T08:30:01Z")

</div>

There’s an issue to track: [https://github.com/jump-dev/JuMP.jl/issues/1198](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.