# Solver differences between sim!(::SimModel) and sim!(::PredictiveController) in ModelPredictiveControl.jl

**URL:** <https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819>\
**Category:** Modelling & Simulations\
**Created:** [October 2, 2024, 3:46pm UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819 "2024-10-02T15:46:40Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Bart\_van\_de\_Lint](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bart_van_de_lint/32/212161_2.png) [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)\
**Post date:** [October 2, 2024, 3:46pm UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/1 "2024-10-02T15:46:40Z")

</div>

I am trying to use NMPC in ModelPredictiveControl.jl, but while my model runs with `sim!(model)` it gets a boundserror with `sim!(nmpc)`.

To try to find out what is going on, I tried to run the two sim! functions on a simpler model. When running the simpler model, the results of `res` and `res_ry` are not exactly the same, even though I set the constraints of the NMPC so that `u = 0.5`. Any idea why these are different @franckgaga? Is there any difference between the solvers in the two `sim!` functions?

```julia
model = setname!(NonLinModel(f!, h!, Ts, nu, nx, ny, solver=RungeKutta(4; supersample=10)); u=vu, x=vx, y=vy)

println("sim")
u = [0.5]
N = 35
res = sim!(model, N, u, x_0 = [0, 0])
display(plot(res, plotu=false))

println("nmpc")
Hp, Hc, Mwt, Nwt = 20, 2, [0.5], [2.5]
nmpc = NonLinMPC(model; Hp, Hc, Mwt, Nwt, Cwt=Inf)
umin, umax = [0.5], [0.5]
nmpc = setconstraint!(nmpc; umin, umax)

res_ry = sim!(nmpc, N, [0.0], plant=model, x_0=[0, 0], x̂_0=[0, 0, 0])
display(plot(res_ry, plotu=false))

```

```julia
julia> println(res.X_data[1:10])
[0.0, 0.0, 0.8240082452793688, 0.04487962009129703, 1.2137241016911229, 0.15021540447023307, 1.23805728754557, 0.27538582798594163, 1.018826941414315, 0.3897021060356283]

julia> println(res_ry.X_data[1:10])
[0.0, 0.0, 0.0, 0.0, 0.8240082289015881, 0.0448796191992897, 1.2137240774825946, 0.1502154014804102, 1.2380572626538922, 0.27538582249114774]

```

---

<div class="post-metadata">

**Author:** ![franckgaga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/franckgaga/32/218241_2.png) [@franckgaga](https://discourse.julialang.org/u/franckgaga)\
**Post date:** [October 2, 2024, 4:06pm UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/2 "2024-10-02T16:06:19Z")

</div>

It’s presumably caused but the time limit in the solver. Do you see a warning like this ?

```julia
┌ Warning: MPC termination status not OPTIMAL or LOCALLY_SOLVED: keeping solution anyway
│ status = TIME_LIMIT::TerminationStatusCode = 12
└ @ ModelPredictiveControl ~/.julia/dev/ModelPredictiveControl/src/controller/execute.jl:459

```

I activated a time limit at `mpc.estim.model.Ts` by default using `JuMP.set_time_limit_sec(mpc.optim, mpc.estim.model.Ts)`. You can disable it with:

```julia
using JuMP; unset_time_limit_sec(nmpc.optim)

```

Tight constraint like above means a problem harder to solve for Ipopt (even if the solution seems obvious to us lol)

maybe this aspect should be more clear in the documentation

---

<div class="post-metadata">

**Author:** ![Bart\_van\_de\_Lint](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bart_van_de_lint/32/212161_2.png) [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)\
**Post date:** [October 2, 2024, 4:27pm UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/3 "2024-10-02T16:27:53Z")

</div>

> [@franckgaga](#):
>
> `using JuMP; unset_time_limit_sec(nmpc.optim)`

When running unset\_time\_limit and unset\_silent on the complex model, I get the following result. I still get NaNs in the sim!(nmpc), even though sim!(model) works fine.

```julia-auto
This is Ipopt version 3.14.16, running with linear solver MUMPS 5.7.3.

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

Number of Iterations....: 0

Number of objective function evaluations = 0
Number of objective gradient evaluations = 1
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 1
Number of equality constraint Jacobian evaluations = 0
Number of inequality constraint Jacobian evaluations = 1
Number of Lagrangian Hessian evaluations = 0
Total seconds in IPOPT = 0.373

EXIT: Invalid number in NLP function or derivative detected.
┌ Error: MPC terminated without solution: returning last solution shifted
│ status = INVALID_MODEL::TerminationStatusCode = 21
└ @ ModelPredictiveControl ~/.julia/packages/ModelPredictiveControl/OoBo3/src/controller/execute.jl:456

```

---

<div class="post-metadata">

**Author:** ![Bart\_van\_de\_Lint](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bart_van_de_lint/32/212161_2.png) [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)\
**Post date:** [October 2, 2024, 7:10pm UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/4 "2024-10-02T19:10:43Z")

</div>

In other words: in which situations does `sim!(model)` work, while `sim!(nmpc)` causes a `BoundsError`?

---

<div class="post-metadata">

**Author:** ![franckgaga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/franckgaga/32/218241_2.png) [@franckgaga](https://discourse.julialang.org/u/franckgaga)\
**Post date:** [October 2, 2024, 7:11pm UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/5 "2024-10-02T19:11:11Z")

</div>

This problem can be hard to debug.

From my experience, one possible cause is a `NonLinModel` in which the `f`/`h` function are not differentiable with ForwardDiff.jl at a specific point (but presumably differentiable at the first guess \mathbf{\Delta U = 0} and \epsilon=0, else Ipopt would have not been called at all). A sudden division by 0 is an exemple of code that is not differentiable with `ForwardDiff.jl`.

You can try put some:

```julia
if ~all(isfinite(ẋ))
    println(ẋ)
    # print other intermediary results that was used to compute ẋ
end

```

at the end of you `f` function, or a similar block of code at the end of you `h` function. It can help you to troubleshoot.

> [@Bart\_van\_de\_Lint](#):
>
> In other words: in which situations does `sim!(model)` work, while `sim!(nmpc)` causes a `BoundsError`?

`sim!(model)` implies no optimization at all, it’s just an open-loop simulation of the plant model at a constant \mathbf{u} value. `sim!(nmpc)` will call `JuMP.optimize!`, thus many automatic differentiation with ForwardDiff.jl and solving using the default Ipopt.jl solver.

edit:

can you copy-paste the stack trace of your `BoundsError` error here ?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [October 3, 2024, 12:13am UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/6 "2024-10-03T00:13:45Z")

</div>

You can try to solve the problem with a pure Julia solver in order to get a better stack trace. MadNLP.jl is a Julia-only version of Ipopt that I often use.

You can also try differentiating the observed function `h` with ForwardDiff manually, if that fails, I believe the problem is that JuliaSimCompiler hardcodes a Float64 argument into the out-of-place version of the built function, and in that case you can solve the issue by this kind of wrapper

```julia
h = build_explicit_observed_function(ssys, outputs; target = JuliaTarget())
h_oop = function (x, p, t)
    y = zeros(promote_type(eltype(x), eltype(p)), length(outputs))
    h(y, x, p, t)
    y
end

```

---

<div class="post-metadata">

**Author:** ![Bart\_van\_de\_Lint](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bart_van_de_lint/32/212161_2.png) [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)\
**Post date:** [October 3, 2024, 8:35am UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/7 "2024-10-03T08:35:43Z")

</div>

There are definitely a lot of NaN’s, long before the program crashes. So the problem probably lies in the forwarddiff.

The error message pasted below. I don’t call log at any point in my model.

```julia
┌ Error: MPC terminated without solution: returning last solution shifted
│ status = INVALID_MODEL::TerminationStatusCode = 21
└ @ ModelPredictiveControl ~/.julia/packages/ModelPredictiveControl/OoBo3/src/controller/execute.jl:456
ERROR: LoadError: DomainError with -1.0171433030139918e-5:
log was called with a negative real argument but will only return a complex result if called with a complex argument. Try log(Complex(x)).
DomainError detected in the user `f` function. This occurs when the domain of a function is violated.
For example, `log(-1.0)` is undefined because `log` of a real number is defined to only output real
numbers, but `log` of a negative number is complex valued and therefore Julia throws a DomainError
by default. Cases to be aware of include:

* `log(x)`, `sqrt(x)`, `cbrt(x)`, etc. where `x<0`
* `x^y` for `x<0` floating point `y` (example: `(-1.0)^(1/2) == im`)

Within the context of SciML, this error can occur within the solver process even if the domain constraint
would not be violated in the solution due to adaptivity. For example, an ODE solver or optimization
routine may check a step at `new_u` which violates the domain constraint, and if violated reject the
step and use a smaller `dt`. However, the throwing of this error will have halted the solving process.

Thus the recommended fix is to replace this function with the equivalent ones from NaNMath.jl
(https://github.com/JuliaMath/NaNMath.jl) which returns a NaN instead of an error. The solver will then
effectively use the NaN within the error control routines to reject the out of bounds step. Additionally,
one could perform a domain transformation on the variables so that such an issue does not occur in the
definition of `f`.

For more information, check out the following FAQ page:
https://docs.sciml.ai/Optimization/stable/API/FAQ/#The-Solver-Seems-to-Violate-Constraints-During-the-Optimization,-Causing-DomainErrors,-What-Can-I-Do-About-That?

Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float64)
    @ Base.Math ./math.jl:33
  [2] _log(x::Float64, base::Val{:ℯ}, func::Symbol)
    @ Base.Math ./special/log.jl:301
  [3] log
    @ ./special/log.jl:267 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/JuliaSimCompiler/y10jZ/src/backends/julia.jl:214 [inlined]
  [5] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:163 [inlined]
  [6] macro expansion
    @ ./none:0 [inlined]
  [7] generated_callfunc
    @ ./none:0 [inlined]
  [8] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::SubArray{…}, ::Vector{…}, ::Vector{…}, ::Float64)
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:150
  [9] f!
    @ ~/Code/KitePredictiveControl.jl/src/mtk_interface.jl:19 [inlined]
 [10] (::ModelPredictiveControl.var"#inner_solver_f!#19"{…})(xnext::SubArray{…}, x::SubArray{…}, u::Vector{…}, d::Vector{…}, p::Vector{…})
    @ ModelPredictiveControl ~/.julia/packages/ModelPredictiveControl/OoBo3/src/model/solver.jl:61
 [11] f!
    @ ~/.julia/packages/ModelPredictiveControl/OoBo3/src/model/nonlinmodel.jl:208 [inlined]
 [12] f̂!(x̂next0::SubArray{…}, û0::Vector{…}, estim::UnscentedKalmanFilter{…}, model::NonLinModel{…}, x̂0::Vector{…}, u0::Vector{…}, d0::Vector{…})
    @ ModelPredictiveControl ~/.julia/packages/ModelPredictiveControl/OoBo3/src/estimator/execute.jl:44
 [13] update_estimate!(estim::UnscentedKalmanFilter{Float64, NonLinModel{…}}, y0m::Vector{Float64}, d0::Vector{Float64}, u0::Vector{Float64})
    @ ModelPredictiveControl ~/.julia/packages/ModelPredictiveControl/OoBo3/src/estimator/kalman.jl:850
 [14] updatestate!(estim::UnscentedKalmanFilter{Float64, NonLinModel{…}}, u::Vector{Float64}, ym::Vector{Float64}, d::Vector{Float64})
    @ ModelPredictiveControl ~/.julia/packages/ModelPredictiveControl/OoBo3/src/estimator/execute.jl:277
 [15] updatestate!
    @ ~/.julia/packages/ModelPredictiveControl/OoBo3/src/controller/execute.jl:508 [inlined]
 [16] macro expansion
    @ ~/.julia/packages/ModelPredictiveControl/OoBo3/src/plot_sim.jl:303 [inlined]
 [17] macro expansion
    @ ~/.julia/packages/ProgressLogging/6KXlp/src/ProgressLogging.jl:470 [inlined]
 [18] sim_closedloop!(est_mpc::NonLinMPC{…}, estim::UnscentedKalmanFilter{…}, N::Int64, u_ry::Vector{…}, d::Vector{…}, ru::Vector{…}; plant::NonLinModel{…}, u_step::Vector{…}, u_noise::Vector{…}, y_step::Vector{…}, y_noise::Vector{…}, d_step::Vector{…}, d_noise::Vector{…}, x_noise::Vector{…}, x_0::Vector{…}, xhat_0::Nothing, lastu::Vector{…}, x̂_0::Nothing)
    @ ModelPredictiveControl ~/.julia/packages/ModelPredictiveControl/OoBo3/src/plot_sim.jl:285
 [19] #sim!#157
    @ ~/.julia/packages/ModelPredictiveControl/OoBo3/src/plot_sim.jl:244 [inlined]
 [20] sim!
    @ ~/.julia/packages/ModelPredictiveControl/OoBo3/src/plot_sim.jl:236 [inlined]
 [21] top-level scope
    @ ~/Code/KitePredictiveControl.jl/src/KitePredictiveControl.jl:50
 [22] include(fname::String)
    @ Base.MainInclude ./client.jl:489
 [23] top-level scope
    @ REPL[46]:1
in expression starting at /home/bart/Code/KitePredictiveControl.jl/src/KitePredictiveControl.jl:50
Some type information was truncated. Use `show(err)` to see complete types.

```

When using MadNLP I get the following message:

```julia
EXIT: Invalid number in NLP objective gradient detected.

```

So apparantly the gradients of my model are invalid / NaN. Is there any way to debug what parts of my model are causing this? I will try to run Forwarddiff manually.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [October 3, 2024, 9:43am UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/8 "2024-10-03T09:43:19Z")

</div>

Well, if the compression stiffness of your tether segments is different from the extension stiffness, in other words, if you are using **ifelse** statements in your equations the function is not longer differentiable.

This problem can be avoided by using differentiable versions of curtain functions, for example

```julia
# differentiable version of the sign function
function smooth_sign(x)
    EPSILON = 6
    x / sqrt(x * x + EPSILON * EPSILON)
end

```

---

<div class="post-metadata">

**Author:** ![Bart\_van\_de\_Lint](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bart_van_de_lint/32/212161_2.png) [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)\
**Post date:** [October 3, 2024, 10:24am UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/9 "2024-10-03T10:24:38Z")

</div>

Thanks, I am using ifelse statements. I will try using smooth\_sign.

---

<div class="post-metadata">

**Author:** ![Bart\_van\_de\_Lint](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bart_van_de_lint/32/212161_2.png) [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)\
**Post date:** [October 3, 2024, 3:57pm UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/10 "2024-10-03T15:57:49Z")

</div>

Calculating the jacobian works fine, without creating NaN, as long as x is set to x\_0:

```julia
f_grad(res, state) = f!(res, state, zeros(3), nothing, nothing)
@show ForwardDiff.jacobian(f_grad, zeros(length(x_0)), x_0)
70×70 Matrix{Float64}:
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 ⋮ ⋮ ⋮ ⋱ ⋮ ⋮                              
 0.0 0.0 34.4607 0.0 0.0 0.0 0.0 0.0 0.0 0.933865 1.0461 0.383042 0.0 -272.262 0.0 0.0 0.0 0.0 0.0 0.0 1185.76 34.6067 466.945
 0.0 0.0 23509.9 0.0 0.0 0.0 0.0 0.0 0.0 24.4376 0.383042 8.17345 0.0 -14309.5 0.0 0.0 0.0 0.0 0.0 0.0 29787.5 466.893 11746.8
 0.0 0.0 -359.463 0.0 0.0 0.0 0.0 0.0 0.0 -0.534668 -0.00838054 -0.218609 -108.622 263.334 0.0 0.0 0.0 0.0 0.0 0.0 -668.745 -10.4821 -263.334
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

```

When using zeros instead of x\_0, it creates a lot of NaN’s in the jacobian.

x\_0 is a known stable initial state:

```julia
julia> x_0
70-element Vector{Float64}:
 51.717090693272105
 57.697634043433354
 57.697634043463076
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
 55.33012009110322
  3.174296775796758e-11
  ⋮
  0.0
  0.0
  0.0
 10.57113449102748
 -5.606197393920564e-19
 -2.0908552687219248e-17
 -1.9994069460380887e-17
 -1.4220384836768373e-17
 -5.7396544513739685e-21
 -6.607803229225811e-19
 -2.832936859732727e-17
 -1.3222407170511217e-20
 -8.175880683050451e-18

```

---

<div class="post-metadata">

**Author:** ![Bart\_van\_de\_Lint](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bart_van_de_lint/32/212161_2.png) [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)\
**Post date:** [October 3, 2024, 4:03pm UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/11 "2024-10-03T16:03:49Z")

</div>

Removing the ifelse statements did not help. If I catch and discard the DomainError, I get the following error multiple times:

```julia
This is MadNLP version v0.8.4, running with umfpack

Number of nonzeros in constraint Jacobian............: 13
Number of nonzeros in Lagrangian Hessian.............: 0

Total number of variables............................: 4
                     variables with only lower bounds: 0
                variables with lower and upper bounds: 0
                     variables with only upper bounds: 0
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 13
        inequality constraints with only lower bounds: 0
   inequality constraints with lower and upper bounds: 0
        inequality constraints with only upper bounds: 13

EXIT: Invalid number in NLP objective gradient detected.
┌ Error: MPC terminated without solution: returning last solution shifted
│ status = INVALID_MODEL::TerminationStatusCode = 21
└ @ ModelPredictiveControl ~/.julia/packages/ModelPredictiveControl/OoBo3/src/controller/execute.jl:456

```

---

<div class="post-metadata">

**Author:** ![franckgaga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/franckgaga/32/218241_2.png) [@franckgaga](https://discourse.julialang.org/u/franckgaga)\
**Post date:** [October 3, 2024, 4:14pm UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/12 "2024-10-03T16:14:09Z")

</div>

You may have already figured that but `log` is called in `JuliaSimCompiler`. You can maybe temporarily add some print statements directly in `/.julia/packages/JuliaSimCompiler/y10jZ/src/backends/julia.jl` to guide you.

In my experience, calling an optimizer on a simulator is kinda the ultimate final test of its codebase. You may had done multiple “what if” simulations and experimental validations of the model with 0 bug in the past. There will always be some forgotten corner case left that the optimizer will manage to expose. The good news is that in many cases it only require the addition of statements like `x < 0 && #force assign another variable to NaN or 0`. Explicit `NaN` assignement can be useful: it tells the optimizer that the objective function is not defined in this region.

In the worst case it needs modifications in the plant model assumptions/equations, as mentioned by @ufechner7.

But yeah, hard to debug, as I said.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [October 3, 2024, 4:29pm UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/13 "2024-10-03T16:29:00Z")

</div>

To find issues like this, I usually copy-paste the code from the RuntimeGeneratedFunction returned by `generate_control_function` into a file and then `include` that file. That way I have a normal julia function I can debug as I want.

If you do that, you would get a stacktrace that points you to the offending line in the Julia code instead of the cryptic stack frames `[4]-[8]` above.

---

<div class="post-metadata">

**Author:** ![Bart\_van\_de\_Lint](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bart_van_de_lint/32/212161_2.png) [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)\
**Post date:** [October 3, 2024, 6:17pm UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/14 "2024-10-03T18:17:59Z")

</div>

How do you copy-paste the function? Using

```julia
    open("f.jl", "w") do file
        println(file, RuntimeGeneratedFunctions.get_expression(f_ip))
    end

```

causes syntax errors when trying to include “f.jl”.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [October 3, 2024, 6:23pm UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/15 "2024-10-03T18:23:16Z")

</div>

If you look at the written code, you might see some stuff in the beginning that needs to be pruned

---

<div class="post-metadata">

**Author:** ![Bart\_van\_de\_Lint](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bart_van_de_lint/32/212161_2.png) [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)\
**Post date:** [October 4, 2024, 6:55am UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/16 "2024-10-04T06:55:20Z")

</div>

Thanks for the help, I found the log function hidden in a getter function… I am not getting a domainerror anymore, but the` EXIT: Invalid number in NLP objective gradient detected.` are not gone yet. One thing I don’t understand: why aren’t there 3 variables with lower and upper bounds here:

```julia
This is MadNLP version v0.8.4, running with umfpack

Number of nonzeros in constraint Jacobian............: 13
Number of nonzeros in Lagrangian Hessian.............: 0

Total number of variables............................: 4
                     variables with only lower bounds: 0
                variables with lower and upper bounds: 0
                     variables with only upper bounds: 0
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 13
        inequality constraints with only lower bounds: 0
   inequality constraints with lower and upper bounds: 0
        inequality constraints with only upper bounds: 13

```

As I have 3 inputs that are bound to `-max, max` each.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [October 4, 2024, 7:01am UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/17 "2024-10-04T07:01:06Z")

</div>

Likely because ModelPredictiveControl.jl is using generic constraints instead of variable bounds. This, in turn, might be because there is an implicit integrator on the input, so that the optimization variables are really \Delta u instead of u. @franckgaga knows this better

---

<div class="post-metadata">

**Author:** ![franckgaga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/franckgaga/32/218241_2.png) [@franckgaga](https://discourse.julialang.org/u/franckgaga)\
**Post date:** [October 4, 2024, 1:03pm UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/18 "2024-10-04T13:03:24Z")

</div>

There is always a single additional \epsilon variable if constraint softening is activated. If you only have bound on the manipulated inputs, you can disable softening with the `Cwt=Inf` option.

I do not use pure lower and upper bounds on the decision variable for this reason (that’s also why you’re seeing `variables with lower and upper bounds: 0`). Supporting relaxation of \mathbf{u} and \mathbf{\Delta u} require linear inequality constraints \mathbf{A \Delta U \le b} (that’s the line `inequality constraints with only upper bounds: 13` above). Some nonlinear solver like Ipopt AFAIK are able to take advantage of linear inequality constraints. I’m not sure if MadNLP treats them as nonlinear…

edit: and yes the decision variables are the input increments over the prediction horizon \mathbf{\Delta U} (strictly speaking: \mathbf{\Delta \tilde{U}} that is augmented with \epsilon, if softening is enable). The constraint \mathbf{u\_{min}} and \mathbf{u\_{max}} are truly linear equality constraint in this context (because of the integrator for the conversion from input increments to inputs)

---

<div class="post-metadata">

**Author:** ![franckgaga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/franckgaga/32/218241_2.png) [@franckgaga](https://discourse.julialang.org/u/franckgaga)\
**Post date:** [October 4, 2024, 4:45pm UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/19 "2024-10-04T16:45:21Z")

</div>

You’re quick Fredrik! 🧠🙂 🐆

---

<div class="post-metadata">

**Author:** ![franckgaga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/franckgaga/32/218241_2.png) [@franckgaga](https://discourse.julialang.org/u/franckgaga)\
**Post date:** [June 24, 2026, 5:15pm UTC](https://discourse.julialang.org/t/solver-differences-between-sim-simmodel-and-sim-predictivecontroller-in-modelpredictivecontrol-jl/120819/20 "2026-06-24T17:15:12Z")

</div>

Just a quick follow-up for the posterity. As for [v2.5.0](https://discourse.julialang.org/t/ann-modelpredictivecontrol-jl/100147/75), the bounds on the slack variable \epsilon, the input increments \mathbf{\Delta U} and the terminal states \mathbf{\hat{x}}(k+H\_p) are now treated as box constraints. And pretty-printing an MPC object now shows the dimensions of the optimization problem, which is just clearer for the user.