Assertion error with `FastDifferentiation` and `ModelPredictiveControl`

[franckgaga]

I'm raising an issue here, related this discussion : JuliaControl/ModelPredictiveControl.jl#174

I think it's a bug/limitation in FastDifferentiation. Here's a not-so-minimal exemple:

using ModelPredictiveControl
using DifferentiationInterface
import FastDifferentiation
function f!(ẋ, x, u, _ , p)
    g, L, K, m = p       # [m/s²], [m], [kg/s], [kg]
    θ, ω = x[1], x[2]    # [rad], [rad/s]
    τ = u[1]             # [Nm]
    ẋ[1] = ω
    ẋ[2] = -g/L*sin(θ) - K/m*ω + τ/m/L^2
end
h!(y, x, _ , _ ) = (y[1] = 180/π*x[1])   # [°]
p = [9.8, 0.4, 1.2, 0.3]
nu = 1; nx = 2; ny = 1; Ts = 0.1
model = NonLinModel(f!, h!, Ts, nu, nx, ny; p)

σQ = [0.1, 1.0]; σR=[5.0]; nint_u=[1]; σQint_u=[0.1]
estim = UnscentedKalmanFilter(model; σQ, σR, nint_u, σQint_u)

p_plant = copy(p); p_plant[3] = 1.25*p[3]
plant = NonLinModel(f!, h!, Ts, nu, nx, ny; p=p_plant)
N = 35; u = [0.5]; 
res = sim!(estim, N, u; plant, y_noise=[0.5]) # no error by commenting out this line ! 😮😮😮

Hp, Hc, Mwt, Nwt, Cwt = 20, 2, [0.5], [2.5], Inf
backend = AutoFastDifferentiation() 
transcription = MultipleShooting()
nmpc = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, transcription, jacobian=backend)

giving this error:

ERROR: AssertionError: Should only be one path from root 5 to variable 1. Instead have 2 children from node 122 on the path
Stacktrace:
  [1] follow_path(a::FastDifferentiation.DerivativeGraph{Int64}, root_index::Int64, var_index::Int64)
    @ FastDifferentiation ~/.julia/packages/FastDifferentiation/kJ1qL/src/Factoring.jl:525
  [2] evaluate_path
    @ ~/.julia/packages/FastDifferentiation/kJ1qL/src/Factoring.jl:555 [inlined]
  [3] _symbolic_jacobian!(graph::FastDifferentiation.DerivativeGraph{…}, partial_variables::Vector{…})
    @ FastDifferentiation ~/.julia/packages/FastDifferentiation/kJ1qL/src/Jacobian.jl:47
  [4] _symbolic_jacobian
    @ ~/.julia/packages/FastDifferentiation/kJ1qL/src/Jacobian.jl:61 [inlined]
  [5] jacobian
    @ ~/.julia/packages/FastDifferentiation/kJ1qL/src/Jacobian.jl:90 [inlined]
  [6] prepare_jacobian_nokwarg(::Val{…}, ::Function, ::Vector{…}, ::AutoFastDifferentiation, ::Vector{…}, ::Cache{…}, ::Cache{…}, ::Cache{…}, ::Cache{…}, ::Cache{…}, ::Cache{…}, ::Cache{…}, ::Cache{…}, ::Cache{…}, ::Cache{…}, ::Cache{…})
    @ DifferentiationInterfaceFastDifferentiationExt ~/.julia/packages/DifferentiationInterface/7eD1K/ext/DifferentiationInterfaceFastDifferentiationExt/twoarg.jl:320
  [7] #prepare_jacobian#48
    @ ~/.julia/packages/DifferentiationInterface/7eD1K/src/first_order/jacobian.jl:23 [inlined]
  [8] prepare_jacobian
    @ ~/.julia/packages/DifferentiationInterface/7eD1K/src/first_order/jacobian.jl:15 [inlined]
  [9] get_optim_functions(mpc::NonLinMPC{…}, ::JuMP.Model)
    @ ModelPredictiveControl ~/.julia/dev/ModelPredictiveControl/src/controller/nonlinmpc.jl:706
 [10] init_optimization!(mpc::NonLinMPC{…}, model::NonLinModel{…}, optim::JuMP.Model)
    @ ModelPredictiveControl ~/.julia/dev/ModelPredictiveControl/src/controller/nonlinmpc.jl:543
 [11] (NonLinMPC{…})(estim::UnscentedKalmanFilter{…}, Hp::Int64, Hc::Int64, M_Hp::Matrix{…}, N_Hc::Matrix{…}, L_Hp::Matrix{…}, Cwt::Float64, Ewt::Float64, JE::ModelPredictiveControl.var"#159#163", gc!::ModelPredictiveControl.var"#160#164", nc::Int64, p::Vector{…}, transcription::MultipleShooting, optim::JuMP.Model, gradient::AutoForwardDiff{…}, jacobian::AutoFastDifferentiation)
    @ ModelPredictiveControl ~/.julia/dev/ModelPredictiveControl/src/controller/nonlinmpc.jl:126
 [12] NonLinMPC(estim::UnscentedKalmanFilter{…}; Hp::Int64, Hc::Int64, Mwt::Vector{…}, Nwt::Vector{…}, Lwt::Vector{…}, M_Hp::Matrix{…}, N_Hc::Matrix{…}, L_Hp::Matrix{…}, Cwt::Float64, Ewt::Float64, JE::ModelPredictiveControl.var"#159#163", gc!::ModelPredictiveControl.var"#160#164", gc::ModelPredictiveControl.var"#160#164", nc::Int64, p::Vector{…}, transcription::MultipleShooting, optim::JuMP.Model, gradient::AutoForwardDiff{…}, jacobian::AutoFastDifferentiation)
    @ ModelPredictiveControl ~/.julia/dev/ModelPredictiveControl/src/controller/nonlinmpc.jl:393
 [13] top-level scope
    @ ~/Dropbox/Programmation/Julia/TestMPC/src/debug_FD.jl:28
Some type information was truncated. Use `show(err)` to see complete types.

Now comes the very weird part: if I comment out the line where I simulate the state estimator estim, there is no error and the mpc objects work as expected. I really don't know what's happening here.

[gdalle]

@brianguenter any idea what the error means?

[brianguenter]

Unfortunately, yes I do. It's a bug in the graph factorization code, or maybe several bugs. This is the most complex part of FastDifferentiation and debugging it is a pain so I've been putting this off. Too long. Can you send me a MWE? I'll see If I can fix the problem.

[franckgaga]

Above it's the most minimal I could get to reproduce. The error seems to come and go randomly, as mentioned above: just commenting out one line makes the error disappear. I can open an issue on your repo with this exemple.

[brianguenter]

What is going on inside NonLinMPC? Presumably it's computing a Jacobian over multiple time steps. Do you know if the symbolic Jacobian is computed just once at the beginning of the iteration and is then evaluated many times for different inputs? Or is the structure of the Jacobian changing at each time step?

[franckgaga]

In the code above there is no (explicit) Jacobian computation at all, since NonLinMPC object is constructed but never executed. I call prepare_jacobian of DifferentiationInterface.jl during construction to pre-allocate everything, and that's in prepare_jacobian that FastDifferentiation crashes.

[brianguenter]

just looked at the docs for NonLinMPC. The jacobian is for the non-linear constraints. Is this jacobian computed just once on entry to NonLinMPC and then evaluated multiple times? If this is the case it will be easier to extract the derivative graph that is causing problems.

Here's the equation from NonLinMPC docs:

is the Jₑ the function that encapsulates the non-linear constraints? In which case I could take the Jacobian of Jₑ directly rather than having to run NonLinMPC and debug through that function call.

[brianguenter]

Ok. I'll look at prepare_jacobian docs then. This is called inside NonLinMPC I assume?

[franckgaga]

This is called inside NonLinMPC constructor yes, in the last steps of construction. Here:

https://github.com/JuliaControl/ModelPredictiveControl.jl/blob/37058219cd0de86df95ec6e912ebf70514e1e479/src/controller/nonlinmpc.jl#L706

[brianguenter]

Reproduced the error on my machine.

The derivative graph for your Jacobian has 1339 nodes. Debugging an example that large is difficult and time consuming. If you need this code to work right now you should probably switch to another AD algorithm because it could take days, maybe weeks to track this down.

[franckgaga]

No, not needed right now, I was just comparing the performance of various AD backends in my package.

I'll raise an issue on your github if I'm able to produce a simpler MRE. Thanks for your time!