Parameter Estimation of ODE system where the response variables are not the state variables

jonlym · August 23, 2022, 3:46am

I am new to Julia and I would like to write a program to estimate the kinetic parameters for a chemical system. I am aware of Catalyst.jl but my work involves systems with complex geometries so I often have to modify my equations accordingly.

The MWE below depicts a plug-flow reactor or PFR (which is essentially a cylinder) where reactants are flowed through one end and is converted to products along the axial direction, z. For a hypothetical reaction system:

A + B \rightarrow C ... 1
2A \rightarrow D ... 2

we can predict how much reaction has occurred by writing the ODEs:

\frac{d\zeta_1} {dz} = k_1 x_A x_B
\frac{d\zeta_2} {dz} = k_2 x_A^2
x_j = \frac{F_j} {\sum_k F_k}
F_j = F_{j0} + \sum_i\zeta_i \nu_{ij}

where:

\zeta_i is the extent of reaction i (the state variable),
k_i is the rate constant of reaction i (the parameter I want to estimate),
x_j is the mole fraction of species j,
F_j is the molar flow rate of species j (the response variables),
F_{j0} is the feed molar flow rate,
\nu is the stoichiometric matrix.

The crux of the problem is the state variable is \zeta_i but the response variables are F_j. From what I’ve seen on the tutorials for DiffEqParamEst.jl, the loss function is defined with respect to the state variables.

Would I have to write my own loss function? Or perhaps I could define F_j as state variables and calculate them from ζ_i using DAEs? Any help is greatly appreciated!

Code to solve ODEs

using CSV
using DataFrames
using DifferentialEquations
using Plots

gr();

"""ODE representing steady-state, isothermal, isobaric 1D PFR""" 
function pfr_ode!(dζ_dz, ζ, logk, z, v_mat, Fj0)
    # Calculate rate constants
    k = 10.0.^logk;
    # Calculate molar flow rate
    Fj = Fj0 .+ ζ * v_mat;
    # Calculate mole fractions
    xj = Fj / sum(Fj);
    # Calculate derivatives
    dζ_dz[1] = k[1] * xj[1] * xj[2];
    dζ_dz[2] = k[2] * (xj[1]^2.0);
    return nothing;
end

"""Converts extent DataFrame to species DataFrame"""
function get_species_df(ζ_df, Fj0, v_mat)
    ζ = Matrix(ζ_df[:, ["ζ1", "ζ2"]]);
    Fj = repeat(Fj0, size(ζ)[1]) .+ ζ * v_mat;
    Fj_df = DataFrame(z=ζ_df[:, "z"], F_A=Fj[:, 1], F_B=Fj[:, 2], F_C=Fj[:, 3],
                      F_D=Fj[:, 4]);
    return Fj_df;
end

function main()
     ζ0 = [0. 0.];      # Initial extent of reaction in mol/s. (The initial condition)
     zspan = (0.0,1.0); # Dimensionless axial position

     species_names = ["A" "B" "C" "D"];

     # Experimental settings
     Fj0 = [10. 20. 0. 0.]; # Feed molar flow rate (mol/s)

     # Stoichiometric matrix
     #         A   B   C   D
     v_mat = [-1. -1.  1.  0.;  # A + B -> C
              -2.  0.  0.  1];  # 2A -> D
     logk0 = [1., 2.];

     # Assigning variables that will be known when doing parameter estimation
     pfr_ode1!(dζ_dz, ζ, logk, z) = pfr_ode!(dζ_dz, ζ, logk, z, v_mat, Fj0);

     # Solve the ODE
     prob = ODEProblem(pfr_ode1!, ζ0, zspan, logk0);
     sol = solve(prob);
     println(sol.destats);

     # Converting solution to DataFrame
     ζ_df = DataFrame(sol);
     rename!(ζ_df, ["z", "ζ1", "ζ2"]);

     # Create DataFrame with 'observed' responses for parameter estimation
     Fj_df = get_species_df(ζ_df, Fj0, v_mat);

     # Write results
     CSV.write("pfr_results.csv", Fj_df);

     # Plotting
     fig = plot(Fj_df[:, "z"], Fj_df[:, "F_A"], label="A");
     plot!(Fj_df[:, "z"], Fj_df[:, "F_B"], label="B");
     plot!(Fj_df[:, "z"], Fj_df[:, "F_C"], label="C");
     plot!(Fj_df[:, "z"], Fj_df[:, "F_D"], label="D");
     xaxis!("Dimesionless axial position, z");
     yaxis!("Molar flow rate, F (mol/s)");
     savefig(fig, "pfr_results.png");
end

main();

This produces the following plot:
pfr_results

And the following CSV file:

z,F_A,F_B,F_C,F_D
0.0,10.0,20.0,0.0,0.0
9.999999999999999e-5,9.997556018444103,19.999777796295028,0.000222203704970252,0.001110888925463598
0.0010999999999999998,9.973167030717427,19.99755779461088,0.0024422053891180197,0.012195381946727549
0.011099999999999997,9.734270573617644,19.975559772073733,0.024440227926266128,0.12064459922804459
0.04290790248508085,9.030825578022126,19.907959532572843,0.09204046742715813,0.4385669772753575
0.09755188056087721,7.9951721383517205,19.79969740901408,0.20030259098592074,0.9022626353311793
0.1715344690632295,6.866946373020592,19.667254587787266,0.3327454122127348,1.4001541073833368
0.2710991757834323,5.708462873550774,19.510783102109617,0.4892168978903844,1.901160114279421
0.40051974176133837,4.617088219659976,19.337521890703986,0.6624781092960134,2.3602168355220052
0.5677240764121487,3.6402107192720203,19.15230381367386,0.8476961863261382,2.7560465472009206
0.7813853351060724,2.808067119889561,18.961835249559847,1.038164750440154,3.0768840648351423
1.0,2.234894364474867,18.80531841976786,1.1946815802321382,3.2852120276464976

What I have so far for my parameter estimation

using CSV
using DataFrames
using DifferentialEquations
using Optim
using DiffEqParamEstim
using NLopt

gr();

"""ODE representing steady-state, isothermal, isobaric 1D PFR""" 
function pfr_ode!(dζ_dz, ζ, logk, z, v_mat, Fj0)
    # Calculate rate constants
    k = 10.0.^logk;
    # Calculate molar flow rate
    Fj = Fj0 .+ ζ * v_mat;
    # Calculate mole fractions
    xj = Fj / sum(Fj);
    # Calculate derivatives
    dζ_dz[1] = k[1] * xj[1] * xj[2];
    dζ_dz[2] = k[2] * (xj[1]^2.0);
    return nothing;
end

function main()
    # Read CSV file for response variables
    Fj_obs_df = CSV.read("./pfr_results.csv", DataFrame);
    z = convert(Array, Fj_obs_df[:, "z"]);
    Fj_obs = Matrix(Fj_obs_df[:, r"F_"]);

    ζ0 = [0. 0.];      # Initial extent of reaction in mol/s. (The initial condition)
    zspan = (0.0,1.0); # Dimensionless axial position
    species_names = ["A" "B" "C" "D"];

    # Experimental settings
    Fj0 = [10. 20. 0. 0.]; # Feed molar flow rate (mol/s)

    # Stoichiometric matrix
    #         A   B   C   D
    v_mat = [-1. -1.  1.  0.;  # A + B -> C
             -2.  0.  0.  1];  # 2A -> D
    logk0 = [1.5, 1.5]; # Initial parameter guesses. True logk0 is [1., 2.]

    # Assigning variables that will be known when doing parameter estimation
    pfr_ode1!(dζ_dz, ζ, logk, z) = pfr_ode!(dζ_dz, ζ, logk, z, v_mat, Fj0);

    # Define the ODE
    prob = ODEProblem(pfr_ode1!, ζ0, zspan, logk0);

    # Build loss function wrt Fj (NOT ζ!)
    cost_func = build_loss_objective(ode_prob, Tsit5(),
                                     L2Loss(z, Fj_obs));

    # Optimize parameters
    result = Optim.optimize(cost_func, logk0)
end

main();

ChrisRackauckas · August 23, 2022, 9:47am

For this you’ll want to directly use SciMLSensitivity.jl and write the loss funciton.

https://sensitivity.sciml.ai/dev/

jonlym · August 24, 2022, 3:02am

Thank you for pointing me in the right direction, Chris! The documentation is very helpful!

I’ll post my code (hopefully by the end of the week) for an example for others.

jonlym · August 24, 2022, 7:01pm

So I modified my code based on these examples:

But when I run it, it throws an error (see code and complete traceback at the bottom of the post). My guess is this section is the most relevant:

warning: Linking two modules of different target triples: 'bcloader' is 'x86_64-w64-windows-gnu' whereas 'text' is 'x86_64-w64-mingw32'

warning: Linking two modules of different target triples: 'bcloader' is 'x86_64-w64-windows-gnu' whereas 'text' is 'x86_64-w64-mingw32'

warning: Linking two modules of different target triples: 'bcloader' is 'x86_64-w64-windows-gnu' whereas 'text' is 'x86_64-w64-mingw32'

┌ Warning: Using fallback BLAS replacements, performance may be degraded
└ @ Enzyme.Compiler C:\Users\Jon\.julia\packages\GPUCompiler\jVY4I\src\utils.jl:35
ERROR: LoadError: Compiling Tuple{Type{Dict}, Base.Generator{Base.Iterators.Enumerate{Vector{Symbol}}, SciMLBase.var"#516#517"}}: try/catch is not supported.
Refer to the Zygote documentation for fixes.
https://fluxml.ai/Zygote.jl/latest/limitations

My code doesn’t use any try/catch statements, so I think I need to resolve the two modules having different target triples? I’m not sure what that means. Do you have any advice?

Code

using CSV
using DataFrames
using DifferentialEquations
using DiffEqParamEstim
using Optimization
using OptimizationOptimJL
using LinearAlgebra
using OptimizationPolyalgorithms
using SciMLSensitivity

"""ODE representing steady-state, isothermal, isobaric 1D PFR""" 
function pfr_ode!(dζ_dz, ζ, logk, z, v_mat, Fj0)
    # Calculate rate constants
    k = 10.0.^logk;
    # Calculate molar flow rate
    Fj = Fj0 .+ ζ * v_mat;
    # Calculate mole fractions
    xj = Fj / sum(Fj);
    # Calculate derivatives
    dζ_dz[1] = k[1] * xj[1] * xj[2];
    dζ_dz[2] = k[2] * (xj[1]^2.0);
    return nothing;
end

"""Custom loss function for pfr_ode by calculating L2 norm of Fj"""
function pfr_loss_func(p, prob, z_vec, Fj_obs, v_mat, Fj0)
    sol = DifferentialEquations.solve(prob, Tsit5(), p=p, saveat=z_vec)

    tot_loss = 0.0;
    if (sol.retcode != :Success)
        tot_loss = Inf;
    else
        # Calculate molar flow rates
        ζ_df = DataFrame(sol);
        rename!(ζ_df, ["z", "ζ1", "ζ2"]);
        ζ = Matrix(ζ_df[:, ["ζ1", "ζ2"]]);
        Fj = repeat(Fj0, size(ζ)[1]) .+ ζ * v_mat;

        tot_loss = LinearAlgebra.norm(Fj .- Fj_obs, 2);
    end
    return tot_loss;
end

"""Callback function to observe training"""
function callback(p, l, pred)
    println("Loss: $l");
    println("Parameters: $p");
    println("-")
    return false;
end

function main()
    # Read CSV file for response variables
    Fj_obs_df = CSV.read("./pfr_results.csv", DataFrame);
    z_vec = convert(Array, Fj_obs_df[:, "z"]);
    Fj_obs = Matrix(Fj_obs_df[:, r"F_"]);

    ζ0 = [0. 0.];      # Initial extent of reaction in mol/s. (The initial condition)
    zspan = (0.0,1.0); # Dimensionless axial position
    species_names = ["A" "B" "C" "D"];

    # Experimental settings
    Fj0 = [10. 20. 0. 0.]; # Feed molar flow rate (mol/s)

    # Stoichiometric matrix
    #         A   B   C   D
    v_mat = [-1. -1.  1.  0.;  # A + B -> C
             -2.  0.  0.  1];  # 2A -> D
    logk0 = [1.5, 1.5]; # Initial parameter guesses. True logk0 is [1., 2.]

    # Assigning variables that will be known when doing parameter estimation

    pfr_ode1!(dζ_dz, ζ, logk, z) = pfr_ode!(dζ_dz, ζ, logk, z, v_mat, Fj0);

    # Define the ODE
    prob = DifferentialEquations.ODEProblem(pfr_ode1!, ζ0, zspan, logk0);

    # Assign values to loss function
    pfr_loss_func1(p) = pfr_loss_func(p, prob, z_vec, Fj_obs, v_mat, Fj0);

    # Optimize parameters
    adtype = Optimization.AutoZygote();
    optf = Optimization.OptimizationFunction((x,p)->pfr_loss_func1(x), adtype);
    optprob = Optimization.OptimizationProblem(optf, logk0)

    result_ode = Optimization.solve(optprob,
                                    OptimizationPolyalgorithms.PolyOpt(),
                                    callback=callback,
                                    maxiters=100);
    return nothing;
end

main();

Julia Version and Installed Packages

julia> VERSION
v"1.8.0"

julia> Pkg.status()
Status `C:\Users\Jon\.julia\environments\v1.8\Project.toml`
  [336ed68f] CSV v0.10.4
  [a93c6f00] DataFrames v1.3.4
  [1130ab10] DiffEqParamEstim v1.26.0
  [0c46a032] DifferentialEquations v7.2.0
  [7f7a1694] Optimization v3.8.2
  [36348300] OptimizationOptimJL v0.1.2
  [500b13db] OptimizationPolyalgorithms v0.1.0
  [91a5bcdd] Plots v1.31.7
  [1ed8b502] SciMLSensitivity v7.6.1
  [37e2e46d] LinearAlgebra

Full Traceback

PS C:\Users\Jon\Downloads\toy_julia_model\toy_julia_model> julia .\PfrParamEstMWE.jl 
warning: Linking two modules of different target triples: 'bcloader' is 'x86_64-w64-windows-gnu' whereas 'text' is 'x86_64-w64-mingw32'

warning: Linking two modules of different target triples: 'bcloader' is 'x86_64-w64-windows-gnu' whereas 'text' is 'x86_64-w64-mingw32'

warning: Linking two modules of different target triples: 'bcloader' is 'x86_64-w64-windows-gnu' whereas 'text' is 'x86_64-w64-mingw32'

┌ Warning: Using fallback BLAS replacements, performance may be degraded
└ @ Enzyme.Compiler C:\Users\Jon\.julia\packages\GPUCompiler\jVY4I\src\utils.jl:35
ERROR: LoadError: Compiling Tuple{Type{Dict}, Base.Generator{Base.Iterators.Enumerate{Vector{Symbol}}, SciMLBase.var"#516#517"}}: try/catch is not supported.
Refer to the Zygote documentation for fixes.
https://fluxml.ai/Zygote.jl/latest/limitations

Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:35
  [2] instrument(ir::IRTools.Inner.IR)
    @ Zygote C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\reverse.jl:121
  [3] #Primal#23
    @ C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\reverse.jl:205 [inlined]
  [4] Zygote.Adjoint(ir::IRTools.Inner.IR; varargs::Nothing, normalise::Bool)
    @ Zygote C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\reverse.jl:330
  [5] _generate_pullback_via_decomposition(T::Type)
    @ Zygote C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\emit.jl:101
  [6] #s2772#1066
    @ C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\interface2.jl:28 [inlined]
  [7] var"#s2772#1066"(::Any, ctx::Any, f::Any, args::Any)
    @ Zygote .\none:0
  [8] (::Core.GeneratedFunctionStub)(::Any, ::Vararg{Any})
    @ Core .\boot.jl:582
  [9] _pullback
    @ C:\Users\Jon\.julia\packages\SciMLBase\DPsGA\src\tabletraits.jl:13 [inlined]
 [10] _pullback(::Zygote.Context{false}, ::Type{SciMLBase.AbstractTimeseriesSolutionRows}, ::Vector{Symbol}, ::Vector{Type}, ::Vector{Float64}, ::Vector{Matrix{Float64}})
    @ Zygote C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\interface2.jl:0
 [11] _pullback
    @ C:\Users\Jon\.julia\packages\SciMLBase\DPsGA\src\tabletraits.jl:40 [inlined]
 [12] _pullback(ctx::Zygote.Context{false}, f::typeof(Tables.rows), args::ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, var"#pfr_ode1!#2"{Matrix{Float64}, Matrix{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, var"#pfr_ode1!#2"{Matrix{Float64}, Matrix{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats})
    @ Zygote C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\interface2.jl:0
 [13] _pullback
    @ C:\Users\Jon\.julia\packages\SciMLBase\DPsGA\src\tabletraits.jl:2 [inlined]
 [14] _pullback(ctx::Zygote.Context{false}, f::typeof(Tables.columns), args::ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, var"#pfr_ode1!#2"{Matrix{Float64}, Matrix{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, var"#pfr_ode1!#2"{Matrix{Float64}, Matrix{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats})
    @ Zygote C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\interface2.jl:0
 [15] _pullback
    @ C:\Users\Jon\.julia\packages\DataFrames\zqFGs\src\other\tables.jl:58 [inlined]
 [16] _pullback(::Zygote.Context{false}, ::DataFrames.var"##DataFrame#796", ::Nothing, ::Type{DataFrame}, ::ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, var"#pfr_ode1!#2"{Matrix{Float64}, Matrix{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, var"#pfr_ode1!#2"{Matrix{Float64}, Matrix{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats})
    @ Zygote C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\interface2.jl:0
 [17] _pullback
    @ C:\Users\Jon\.julia\packages\DataFrames\zqFGs\src\other\tables.jl:48 [inlined]
 [18] _pullback(ctx::Zygote.Context{false}, f::Type{DataFrame}, args::ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, var"#pfr_ode1!#2"{Matrix{Float64}, Matrix{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, var"#pfr_ode1!#2"{Matrix{Float64}, Matrix{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats})
    @ Zygote C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\interface2.jl:0
 [19] _pullback
    @ C:\Users\Jon\Downloads\toy_julia_model\toy_julia_model\PfrParamEstMWE.jl:34 [inlined]
 [20] _pullback(::Zygote.Context{false}, ::typeof(pfr_loss_func), ::Vector{Float64}, ::ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, var"#pfr_ode1!#2"{Matrix{Float64}, Matrix{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Vector{Float64}, ::Matrix{Float64}, ::Matrix{Float64}, ::Matrix{Float64})
    @ Zygote C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\interface2.jl:0
 [21] _pullback
    @ C:\Users\Jon\Downloads\toy_julia_model\toy_julia_model\PfrParamEstMWE.jl:79 [inlined]
 [22] _pullback
    @ C:\Users\Jon\Downloads\toy_julia_model\toy_julia_model\PfrParamEstMWE.jl:83 [inlined]
 [23] _pullback(::Zygote.Context{false}, ::var"#1#4"{var"#pfr_loss_func1#3"{ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, var"#pfr_ode1!#2"{Matrix{Float64}, Matrix{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Vector{Float64}}}, ::Vector{Float64}, ::SciMLBase.NullParameters)
    @ Zygote C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\interface2.jl:0
 [24] _apply(::Function, ::Vararg{Any})
    @ Core .\boot.jl:816
 [25] adjoint
    @ C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\lib\lib.jl:203 [inlined]
 [26] _pullback
    @ C:\Users\Jon\.julia\packages\ZygoteRules\AIbCs\src\adjoint.jl:65 [inlined]
 [27] _pullback
    @ C:\Users\Jon\.julia\packages\SciMLBase\DPsGA\src\scimlfunctions.jl:3024 [inlined]
 [28] _pullback(::Zygote.Context{false}, ::OptimizationFunction{true, Optimization.AutoZygote, var"#1#4"{var"#pfr_loss_func1#3"{ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, var"#pfr_ode1!#2"{Matrix{Float64}, Matrix{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Vector{Float64}}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, ::Vector{Float64}, ::SciMLBase.NullParameters)
    @ Zygote C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\interface2.jl:0
 [29] _apply(::Function, ::Vararg{Any})
    @ Core .\boot.jl:816
 [30] adjoint
    @ C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\lib\lib.jl:203 [inlined]
 [31] _pullback
    @ C:\Users\Jon\.julia\packages\ZygoteRules\AIbCs\src\adjoint.jl:65 [inlined]
 [32] _pullback
    @ C:\Users\Jon\.julia\packages\Optimization\6nIwk\src\function\zygote.jl:30 [inlined]
 [33] _pullback(ctx::Zygote.Context{false}, f::Optimization.var"#136#146"{OptimizationFunction{true, Optimization.AutoZygote, var"#1#4"{var"#pfr_loss_func1#3"{ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, var"#pfr_ode1!#2"{Matrix{Float64}, Matrix{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Vector{Float64}}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, SciMLBase.NullParameters}, args::Vector{Float64})
    @ Zygote C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\interface2.jl:0
 [34] _apply(::Function, ::Vararg{Any})
    @ Core .\boot.jl:816
 [35] adjoint
    @ C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\lib\lib.jl:203 [inlined]
 [36] _pullback
    @ C:\Users\Jon\.julia\packages\ZygoteRules\AIbCs\src\adjoint.jl:65 [inlined]
 [37] _pullback
    @ C:\Users\Jon\.julia\packages\Optimization\6nIwk\src\function\zygote.jl:32 [inlined]
 [38] _pullback(ctx::Zygote.Context{false}, f::Optimization.var"#139#149"{Tuple{}, Optimization.var"#136#146"{OptimizationFunction{true, Optimization.AutoZygote, var"#1#4"{var"#pfr_loss_func1#3"{ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, var"#pfr_ode1!#2"{Matrix{Float64}, Matrix{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Vector{Float64}}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, SciMLBase.NullParameters}}, args::Vector{Float64})
    @ Zygote C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\interface2.jl:0
 [39] pullback(f::Function, cx::Zygote.Context{false}, args::Vector{Float64})
    @ Zygote C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\interface.jl:44
 [40] pullback
    @ C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\interface.jl:42 [inlined]
 [41] gradient(f::Function, args::Vector{Float64})
    @ Zygote C:\Users\Jon\.julia\packages\Zygote\qGFGD\src\compiler\interface.jl:96
 [42] (::Optimization.var"#137#147"{Optimization.var"#136#146"{OptimizationFunction{true, Optimization.AutoZygote, var"#1#4"{var"#pfr_loss_func1#3"{ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, var"#pfr_ode1!#2"{Matrix{Float64}, Matrix{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Vector{Float64}}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, SciMLBase.NullParameters}})(::Vector{Float64}, ::Vector{Float64})
    @ Optimization C:\Users\Jon\.julia\packages\Optimization\6nIwk\src\function\zygote.jl:32
 [43] macro expansion
    @ C:\Users\Jon\.julia\packages\OptimizationOptimisers\XLPqT\src\OptimizationOptimisers.jl:35 [inlined]
 [44] macro expansion
    @ C:\Users\Jon\.julia\packages\Optimization\6nIwk\src\utils.jl:35 [inlined]
 [45] __solve(prob::OptimizationProblem{true, OptimizationFunction{true, Optimization.AutoZygote, var"#1#4"{var"#pfr_loss_func1#3"{ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, var"#pfr_ode1!#2"{Matrix{Float64}, Matrix{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Vector{Float64}}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}, SciMLBase.NullParameters, Nothing, Nothing, Nothing, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, opt::Optimisers.Adam{Float64}, data::Base.Iterators.Cycle{Tuple{Optimization.NullData}}; maxiters::Int64, callback::Function, progress::Bool, save_best::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OptimizationOptimisers C:\Users\Jon\.julia\packages\OptimizationOptimisers\XLPqT\src\OptimizationOptimisers.jl:33
 [46] #solve#509
    @ C:\Users\Jon\.julia\packages\SciMLBase\DPsGA\src\solve.jl:71 [inlined]
 [47] __solve(::OptimizationProblem{true, OptimizationFunction{true, Optimization.AutoZygote, var"#1#4"{var"#pfr_loss_func1#3"{ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, var"#pfr_ode1!#2"{Matrix{Float64}, Matrix{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Vector{Float64}}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}, SciMLBase.NullParameters, Nothing, Nothing, Nothing, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::PolyOpt; maxiters::Int64, kwargs::Base.Pairs{Symbol, typeof(callback), Tuple{Symbol}, NamedTuple{(:callback,), Tuple{typeof(callback)}}})
    @ OptimizationPolyalgorithms C:\Users\Jon\.julia\packages\OptimizationPolyalgorithms\5gDHf\src\OptimizationPolyalgorithms.jl:25
 [48] #solve#509
    @ C:\Users\Jon\.julia\packages\SciMLBase\DPsGA\src\solve.jl:71 [inlined]
 [49] main()
    @ Main C:\Users\Jon\Downloads\toy_julia_model\toy_julia_model\PfrParamEstMWE.jl:86
 [50] top-level scope
    @ C:\Users\Jon\Downloads\toy_julia_model\toy_julia_model\PfrParamEstMWE.jl:93
in expression starting at C:\Users\Jon\Downloads\toy_julia_model\toy_julia_model\PfrParamEstMWE.jl:93

jonlym · August 24, 2022, 8:52pm

Even though the error message said try/catch is problematic, I noticed the Zygote documentation said mutating functions can be problematic also. I converted my pfr_ode! function so it would return \frac{d\zeta} {dz} instead of mutating it (see code below) but I still receive the same error.

Converted pfr_ode function

"""ODE representing steady-state, isothermal, isobaric 1D PFR""" 
function pfr_ode(ζ, logk, z, v_mat, Fj0)
    # Calculate rate constants
    k = 10.0.^logk;
    # Calculate molar flow rate
    Fj = Fj0 .+ ζ * v_mat;
    # Calculate mole fractions
    xj = Fj / sum(Fj);
    # Calculate derivatives
    dζ_dz = zeros(size(ζ))
    dζ_dz[1] = k[1] * xj[1] * xj[2];
    dζ_dz[2] = k[2] * (xj[1]^2.0);
    return dζ_dz;
end

jonlym · August 26, 2022, 12:36pm

See Chris’ post below for the actual reason the new code works.
~~I circumvented the issue by using JuMP instead of SciMLSensitivity.~~ The code converges to the true values (i.e. logk1 = 1.0, logk2 = 2.0).

I’d still be interested to learn how to solve the triplets issue if anyone knows!

Solver Output

termination_status(jumodel) = MathOptInterface.LOCALLY_SOLVED
primal_status(jumodel) = MathOptInterface.FEASIBLE_POINT
dual_status(jumodel) = MathOptInterface.NO_SOLUTION
objective_value(jumodel) = 7.267740308589925e-6
value(logk1) = 0.9999995071557675
value(logk2) = 1.9999995819408962

Comparison of observed and predicted data
Param_est_result

Code

using CSV
using DataFrames
using DifferentialEquations
using DiffEqParamEstim
using LinearAlgebra
using NLopt
using JuMP
using OrdinaryDiffEq
using Plots

"""Converts extent DataFrame to species DataFrame"""
function get_species_df(ζ_df, Fj0, v_mat)
    ζ = Matrix(ζ_df[:, ["ζ1", "ζ2"]]);
    Fj = repeat(Fj0, size(ζ)[1]) .+ ζ * v_mat;
    Fj_df = DataFrame(z=ζ_df[:, "z"], F_A=Fj[:, 1], F_B=Fj[:, 2], F_C=Fj[:, 3],
                      F_D=Fj[:, 4]);
    return Fj_df;
end

"""ODE representing steady-state, isothermal, isobaric 1D PFR""" 
function pfr_ode!(dζ_dz, ζ, logk, z, v_mat, Fj0)
    # Calculate rate constants
    k = 10.0.^logk;
    # Calculate molar flow rate
    Fj = Fj0 .+ ζ * v_mat;
    # Calculate mole fractions
    xj = Fj / sum(Fj);
    # Calculate derivatives
    dζ_dz[1] = k[1] * xj[1] * xj[2];
    dζ_dz[2] = k[2] * (xj[1]^2.0);
    return nothing;
end

"""Custom loss function for pfr_ode by calculating L2 norm of Fj"""
function pfr_loss_func(prob, z_vec, Fj_obs, v_mat, Fj0, p)
    sol = DifferentialEquations.solve(prob, Tsit5(), p=p, saveat=z_vec)

    tot_loss = 0.0;
    if (sol.retcode != :Success)
        tot_loss = Inf;
    else
        # Calculate molar flow rates
        ζ_df = DataFrame(sol);
        rename!(ζ_df, ["z", "ζ1", "ζ2"]);
        ζ = Matrix(ζ_df[:, ["ζ1", "ζ2"]]);
        Fj = repeat(Fj0, size(ζ)[1]) .+ ζ * v_mat;

        tot_loss = LinearAlgebra.norm(Fj .- Fj_obs, 2);
    end
    return tot_loss;
end

function main()
    # Read CSV file for response variables
    Fj_obs_df = CSV.read("./pfr_results.csv", DataFrame);
    z_vec = convert(Array, Fj_obs_df[:, "z"]);
    Fj_obs = Matrix(Fj_obs_df[:, r"F_"]);

    ζ0 = [0. 0.]; # Initial extent of reaction in mol/s. (The initial condition)
    zspan = (0.0,1.0); # Dimensionless axial position
    species_names = ["A" "B" "C" "D"];

    # Experimental settings
    Fj0 = [10. 20. 0. 0.]; # Feed molar flow rate (mol/s)

    # Stoichiometric matrix
    #         A   B   C   D
    v_mat = [-1. -1.  1.  0.;  # A + B -> C
             -2.  0.  0.  1];  # 2A -> D
    logk0 = [1.5, 1.5]; # Initial parameter guesses. True logk0 is [1., 2.]

    # Assigning variables that will be known when doing parameter estimation
    pfr_ode1!(dζ_dz, ζ, logk, z) = pfr_ode!(dζ_dz, ζ, logk, z, v_mat, Fj0);

    # Define the ODE
    prob = DifferentialEquations.ODEProblem(pfr_ode1!, ζ0, zspan, logk0);

    # Assign values to loss function
    pfr_loss_func1(p) = pfr_loss_func(prob, z_vec, Fj_obs, v_mat, Fj0, p);

    # Define JuMP model
    juobj(args...) = pfr_loss_func1(args);
    jumodel = JuMP.Model(NLopt.Optimizer);
    set_optimizer_attribute(jumodel, "algorithm", :LD_MMA);

    JuMP.register(jumodel, :juobj, 2, juobj, autodiff=true);
    @variables jumodel begin
        logk1, (start=logk0[1])
        logk2, (start=logk0[2])
    end
    @NLobjective(jumodel, Min, juobj(logk1, logk2));
    
    # Optimize parameters
    JuMP.optimize!(jumodel);

    # Print parameter estimation final status
    @show termination_status(jumodel);
    @show primal_status(jumodel);
    @show dual_status(jumodel);
    @show objective_value(jumodel);
    @show value(logk1);
    @show value(logk2);

    # Create plots
    optim_sol = DifferentialEquations.solve(prob,
                                            Tsit5(),
                                            p=[value(logk1), value(logk2)],
                                            saveat=z_vec);
    ζ_df = DataFrame(optim_sol);
    rename!(ζ_df, ["z", "ζ1", "ζ2"]);
    Fj_optim_df = get_species_df(ζ_df, Fj0, v_mat);

    colors = ("#44B6FB","#EC9F84", "#7EC288", "#DBAAE5");
    plot(xlabel="Dimensionless axial position, z",
         ylabel="Molar flow rate, F (mol/s)");
    for (species_name, color_str) in zip(species_names, colors)
        color = parse(Colorant, color_str);
        Fj_str = "F_" * species_name;
        scatter!(z_vec, Fj_obs_df[:, Fj_str], color=color,
                 label=Fj_str * " obs");
        plot!(z_vec, Fj_optim_df[:, Fj_str], color=color,
              label=Fj_str * " pred");
    end
    savefig("param_est_result.png")

    return nothing;
end

main();

EDIT: Updated SciML to SciMLSensitivity
EDIT1: Highlighted Chris’ response.

ChrisRackauckas · August 26, 2022, 12:40pm

No, what actually happened here is you switch from Zygote to ForwardDiff. I think the issue really comes down to:

jonlym:

        ζ_df = DataFrame(sol);
        rename!(ζ_df, ["z", "ζ1", "ζ2"]);
        ζ = Matrix(ζ_df[:, ["ζ1", "ζ2"]]);
        Fj = repeat(Fj0, size(ζ)[1]) .+ ζ * v_mat;

        tot_loss = LinearAlgebra.norm(Fj .- Fj_obs, 2);

DataFrame operations use try/catch internally, which seems to break Zygote. You can probably build an MWE that is directly doing loss functions on DataFrames without an ODE solve. That would be nice to isolate and open an issue on Zygote.jl for.

But for this, ForwardDiff is going to be faster anyways, so adtype = Optimization.AutoForwardDiff(); is a better idea, and would fix your issue.

jonlym · August 26, 2022, 12:45pm

I see! Thank you for the insight!

ChrisRackauckas · August 26, 2022, 12:46pm

It’s still using SciMLSensitivity in some sense, it’s just ForwardDiff of the solver instead of adjoints (which is only good for small numbers of parameters, which is your case)