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

**URL:** <https://discourse.julialang.org/t/parameter-estimation-of-ode-system-where-the-response-variables-are-not-the-state-variables/86171>\
**Category:** Optimization (Mathematical)\
**Tags:** question, diffeq\
**Created:** [August 23, 2022, 3:46am UTC](https://discourse.julialang.org/t/parameter-estimation-of-ode-system-where-the-response-variables-are-not-the-state-variables/86171 "2022-08-23T03:46:03Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![jonlym](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonlym/32/20007_2.png) [@jonlym](https://discourse.julialang.org/u/jonlym)\
**Post date:** [August 23, 2022, 3:46am UTC](https://discourse.julialang.org/t/parameter-estimation-of-ode-system-where-the-response-variables-are-not-the-state-variables/86171/1 "2022-08-23T03:46:03Z")

</div>

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](https://github.com/SciML/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](https://diffeqparamestim.sciml.ai/stable/tutorials/ODE_inference/), 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**

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/4/0/40f9d9e0b4de264b8615f131e9aa450a3580e772.png)

And the following CSV file:

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

```julia
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();

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 23, 2022, 9:47am UTC](https://discourse.julialang.org/t/parameter-estimation-of-ode-system-where-the-response-variables-are-not-the-state-variables/86171/2 "2022-08-23T09:47:22Z")

</div>

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

[https://sensitivity.sciml.ai/dev/](https://sensitivity.sciml.ai/dev/)

---

<div class="post-metadata">

**Author:** ![jonlym](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonlym/32/20007_2.png) [@jonlym](https://discourse.julialang.org/u/jonlym)\
**Post date:** [August 24, 2022, 3:02am UTC](https://discourse.julialang.org/t/parameter-estimation-of-ode-system-where-the-response-variables-are-not-the-state-variables/86171/3 "2022-08-24T03:02:41Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![jonlym](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonlym/32/20007_2.png) [@jonlym](https://discourse.julialang.org/u/jonlym)\
**Post date:** [August 24, 2022, 7:01pm UTC](https://discourse.julialang.org/t/parameter-estimation-of-ode-system-where-the-response-variables-are-not-the-state-variables/86171/4 "2022-08-24T19:01:00Z")

</div>

So I modified my code based on these examples:

- [Optimization of Ordinary Differential Equations · SciMLSensitivity.jl](https://sensitivity.sciml.ai/dev/ode_fitting/optimization_ode/)
- [Parameter Estimation on Highly Stiff Systems · SciMLSensitivity.jl](https://sensitivity.sciml.ai/dev/ode_fitting/stiff_ode_fit/)

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:

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

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

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

```

---

<div class="post-metadata">

**Author:** ![jonlym](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonlym/32/20007_2.png) [@jonlym](https://discourse.julialang.org/u/jonlym)\
**Post date:** [August 24, 2022, 8:52pm UTC](https://discourse.julialang.org/t/parameter-estimation-of-ode-system-where-the-response-variables-are-not-the-state-variables/86171/5 "2022-08-24T20:52:59Z")

</div>

Even though the error message said try/catch is problematic, I noticed the [Zygote documentation](https://fluxml.ai/Zygote.jl/latest/limitations/) 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**

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

```

---

<div class="post-metadata">

**Author:** ![jonlym](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonlym/32/20007_2.png) [@jonlym](https://discourse.julialang.org/u/jonlym)\
**Post date:** [August 26, 2022, 12:36pm UTC](https://discourse.julialang.org/t/parameter-estimation-of-ode-system-where-the-response-variables-are-not-the-state-variables/86171/6 "2022-08-26T12:36:06Z")

</div>

[See Chris’ post below for the actual reason the new code works.](https://discourse.julialang.org/t/parameter-estimation-of-ode-system-where-the-response-variables-are-not-the-state-variables/86171/7)  
~~I circumvented the issue by using [JuMP](https://jump.dev/JuMP.jl/stable/) 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**

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/f/f/ff81cb88e7cb0c022e13677b7ee8c1dca7c44165.png)

**Code**

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 26, 2022, 12:40pm UTC](https://discourse.julialang.org/t/parameter-estimation-of-ode-system-where-the-response-variables-are-not-the-state-variables/86171/7 "2022-08-26T12:40:59Z")

</div>

> [@jonlym](#):
>
> I circumvented the issue by using [JuMP](https://jump.dev/JuMP.jl/stable/) instead of SciML. The code converges to the true values (i.e. `logk1 = 1.0`, `logk2 = 2.0`).

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

> [@jonlym](#):
>
> ```julia
> ζ_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.

---

<div class="post-metadata">

**Author:** ![jonlym](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonlym/32/20007_2.png) [@jonlym](https://discourse.julialang.org/u/jonlym)\
**Post date:** [August 26, 2022, 12:45pm UTC](https://discourse.julialang.org/t/parameter-estimation-of-ode-system-where-the-response-variables-are-not-the-state-variables/86171/8 "2022-08-26T12:45:19Z")

</div>

I see! Thank you for the insight!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 26, 2022, 12:46pm UTC](https://discourse.julialang.org/t/parameter-estimation-of-ode-system-where-the-response-variables-are-not-the-state-variables/86171/9 "2022-08-26T12:46:50Z")

</div>

> [@jonlym](#):
>
> EDIT: Updated SciML to SciMLSensitivity

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)