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