Search for optimal ODE parameters

This example demonstrates a workflow for finding optimal parameters for an ordinary differential equation (ODE) describing a chain of chemical reactions. We will find the parameters and evaluate how closely they can reproduce the true dynamics of the system.

This problem is solved using the problem-based optimisation workflow.

In this example, we will use the functions of the JuMP.jl library to formulate the optimisation problem, the nonlinear optimisation library Ipopt.jl, the library for solving differential equations DifferentialEquations.jl and the Plots.jl library to visualise the results.

Problem description

Consider a model of a multi-step chain of chemical reactions involving several substances $C$. Substances react with each other with a certain rate $r$ and form new substances. The true rate of reactions for this problem is unknown. But we can make observations a certain number of times in a given time $t$. Based on these observations, we will find the optimal reaction rates $r$.

There are six substances in the model, from $C_1$ to $C_6$, that react as follows:

• one substance $C_1$ and one substance $C_2$ react to form one substance $C_3$ at a rate of $r_1$;
• one substance $C_3$ and one substance $C_4$ react to form one substance $C_5$ at a rate of $r_2$;
• one substance $C_3$ and one substance $C_4$ react to form one substance $C_6$ at a rate of $r_3$.

The reaction rate is proportional to the product of the quantities of the required substances. Thus, if $y_i$ represents the amount of substance $C_i$, then the reaction rate to produce $C_3$ is equal to $r_1y_1y_2$. Similarly, the reaction rate to produce $C_5$ is equal to $r_2y_3y_4$, and the reaction rate to produce $C_6$ is equal to $r_3y_3y_4$.

Thus, we can make a differential equation describing the process of change in our system:

$$ \frac{dy}{dt} = \begin{bmatrix} -r_1y_1y_2 \\ -r_1y_1y_2 \\ -r_2y_3y_4 + r_1y_1y_2 - r_3y_3y_4 \\ -r_2y_3y_4 - r_3y_3y_4 \\ r_2y_3y_4 \\ r_3y_3y_4 \\ \end{bmatrix} $$

Set initial values $y(0)=[1,1,0,1,0,0]$. Such initial values guarantee that all substances will completely react, as a result of which $C_1–C_4$ will approach zero as the time increases.

Setting the libraries

If your environment does not have the latest version of the JuMP package installed , uncomment and run the box below:

Pkg.add(["Ipopt", "DifferentialEquations", "JuMP"])

#Pkg.add("JuMP");

To launch a new version of the library after the installation is complete, click on the "My Account" button:

Then click on the "Stop" button:

Restart the session by pressing the "Start Engee" button:

Connecting libraries

Connect the library JuMP:

using JuMP;

Connect the nonlinear solver library Ipopt:

using Ipopt;

Connect the library for solving differential equations DifferentialEquations:

using DifferentialEquations;

Connect the graphing library Plots:

using Plots;

Configuration of visualisation parameters

Set customised chart size and background plotly:

plotly(size=(1500, 600), xlims=(0, :auto));

Describe a function that defines bounds for $min$ and $max$ axes to approximate the graphs:

function check_bounds(min::Number, max::Number)
    if min > max
        error("Значение min должно быть меньше чем значение max")
    elseif max < min
        error("Значение max должно быть больше чем значение min")
    end
end

check_bounds (generic function with 1 method)

Solving differential equations

Create a function dif_fun, which describes the ODE of our chemical reaction chain:

function dif_fun(du, u, p, t)
    r1, r2, r3 = p

    du[1] = -r1 * u[1] * u[2]
    du[2] = -r1 * u[1] * u[2]
    du[3] = -r2 * u[3] * u[4] + r1 * u[1] * u[2]- r3 * u[3] * u[4]
    du[4] = -r2 * u[3] * u[4] - r3 * u[3] * u[4]
    du[5] = r2 * u[3] * u[4]
    du[6] = r3 * u[3] * u[4]
end

dif_fun (generic function with 1 method)

Set the true values of $r$ and store them in the variable r_true:

r_true = [2.5, 1.2, 0.45];

The values $r_1$, $r_2$ and $r_3$ in the variable r_true are the benchmark that the optimisation algorithm should try to recreate.

Set the response simulation duration and initial values for $С_1 - С_6$:

t_span = (0.0, 5.0);
y_0 = [1.0, 1.0, 0.0, 1.0, 0.0, 0.0];

Create an ODE problem using the function ODEProblem:

prob = ODEProblem(dif_fun, y_0, t_span, r_true);

Solve the problem using the function solve. In brackets, give the name of the solver and the frequency of observations in the argument saveat:

sol_true = solve(prob, Tsit5(), saveat=0.1);

To solve this problem we use the Runge-Kutta Tsipouras solver of order 5 Tsit5(), which is recommended for most non-rigid problems.

Graph the obtained results of true solutions:

colors = [:red, :blue, :green, :orange, :purple, :brown]

plot(layout=(3,2))
for i in 1:6
    plot!(subplot=i, sol_true.t, sol_true[i,:], color=colors[i], title="y($i)", label="Истинное y($i)")
end
display(current())

Solution of the optimisation problem

Describe the objective function obj_fun, which calculates the difference between the true parameters $r_1$, $r_2$ and $r_3$ and the calculated values.

function obj_fun(r1, r2, r3)
    r = [r1, r2, r3]
    prob_new = ODEProblem(dif_fun, y_0, t_span, r)
    sol_new = solve(prob_new, Tsit5(), saveat=0.1)
    return sum(sum((sol_new[i,:] .- sol_true[i,:]).^2) for i in 1:6)
end

obj_fun (generic function with 1 method)

We will need to minimise this function to find optimal solutions.

Create an optimisation problem using the function Model and specify the name of the solver in brackets:

model = Model(Ipopt.Optimizer)

A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Ipopt

Set the model to "silent" mode to hide system messages:

set_silent(model)

Create a variable r, containing three elements - $r_1$, $r_2$ and $r_3$, and set the bounds from 0.1 to 10:

@variable(model, 0.1 <= r[1:3] <= 10)

3-element Vector{VariableRef}:
 r[1]
 r[2]
 r[3]

Since we are dealing with a complex user-defined nonlinear function obj_fun, it must be registered with the function register beforehand in order to be correctly integrated into the optimisation environment of the library JuMP. In brackets specify the number of incoming arguments ($r_1$, $r_2$ and $r_3$) and enable automatic differentiation.

register(model, :obj_fun, 3, obj_fun; autodiff = true)

Create a non-linear target function to find the minimum values of $r_1$, $r_2$ and $r_3$ and apply the previously created function obj_fun:

@NLobjective(model, Min, obj_fun(r[1], r[2], r[3]))

Solve the optimisation problem:

optimize!(model);

Store the optimised parameters in a variable:

optimized_params = value.(r)

3-element Vector{Float64}:
 2.5000119303476147
 1.200000509788957
 0.4499998880470981

Output the results of solving the problem. You can use the function round and the argument digits, to adjust the accuracy to the accuracy of the true parameter data.

optimized_params = round.(value.(r), digits=2)
println("Истинные параметры: ", r_true)
println("Оптимизированные параметры: ", optimized_params)

Истинные параметры: [2.5, 1.2, 0.45]
Оптимизированные параметры: [2.5, 1.2, 0.45]

Optimisation with limited observations

Let us simulate a situation where we cannot observe all components $y$, but only the final results $y(5)$ and $y(6)$, and calculate the rates of all reactions based on this limited information.

Describe the goal function obj_fun_lim, which is a copy of the function obj_fun from the previous section, but returns only the results for $y(5)$ and $y(6)$.

function obj_fun_lim(r1, r2, r3)
    r = [r1, r2, r3]
    prob_new = ODEProblem(dif_fun, y_0, t_span, r)
    sol_new = solve(prob_new, Tsit5(), saveat=0.1)
    return sum(sum((sol_new[i,:] .- sol_true[i,:]).^2) for i in 5:6)
end

obj_fun_lim (generic function with 1 method)

Create an optimisation problem using the function Model and give the name of the solver in brackets:

model_lim = Model(Ipopt.Optimizer)

A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Ipopt

Set the model to "silent" mode to hide system messages:

set_silent(model_lim)

Create a variable r, containing three elements - $r_1$, $r_2$ and $r_3$, and set the bounds from 0.1 to 10:

@variable(model_lim, 0.1 <= r_lim[1:3] <= 10)

3-element Vector{VariableRef}:
 r_lim[1]
 r_lim[2]
 r_lim[3]

Register the function obj_fun_lim with the function register:

register(model_lim, :obj_fun_lim, 3, obj_fun_lim; autodiff = true)

Create a non-linear target function to find the minimum values of $r_1$, $r_2$ and $r_3$ and apply the previously created function obj_fun_lim:

@NLobjective(model_lim, Min, obj_fun_lim(r[1], r[2], r[3]))

Solve the optimisation problem:

optimize!(model_lim)

Output the results of solving the problem. You can use the function round, to round off excessively accurate values if necessary.

println("Истинные параметры: ", r_true)
println("Оптимизированные параметры (ограниченные): ", value.(r_lim))

Истинные параметры: [2.5, 1.2, 0.45]
Оптимизированные параметры (ограниченные): [1.8085855441520966, 1.5505236127674895, 0.5814611991591475]

Create an ODE problem using the function ODEProblem, applying the optimised values r_lim:

prob_opt = ODEProblem(dif_fun, y_0, t_span, value.(r_lim));

Solve the ODE problem:

sol_opt = solve(prob_opt, Tsit5(), saveat=0.1);

Visualise the results for $y_5$ and $y_6$:

y_min = 0 # @param {type:"slider", min:0, max:1, step:0.01}
y_max = 0.6 # @param {type:"slider", min:0, max:1, step:0.01}
check_bounds(y_min, y_max)
x_min = 0 # @param {type:"slider", min:0, max:5, step:0.1}
x_max = 5 # @param {type:"slider", min:0, max:5, step:0.1}
check_bounds(x_min, x_max)

plot(sol_true.t, [sol_true[5,:] sol_true[6,:]], label=["Истинное y(5)" "Истинное y(6)"], ylimits=(y_min, y_max), xlimits=(x_min, x_max))
plot!(sol_true.t, [sol_opt[5,:] sol_opt[6,:]], label=["Оптимизированное y(5)" "Оптимизированное y(6)"], linestyle=:dash, ylimits=(y_min, y_max), xlimits=(x_min, x_max))
display(current())

Visualise and compare the results of the previously set true values and the obtained values $y$:

true_colors = [:darkred, :darkblue, :darkgreen, :darkorange, :darkviolet, :brown]
opt_colors = [:red, :blue, :green, :orange, :violet, :orange]

plot(layout=(3,2))
for i in 1:6
    plot!(subplot=i, sol_true.t, sol_true[i,:], color=true_colors[i], label="Истинное y($i)", title="y($i)")
    plot!(subplot=i, sol_opt.t, sol_opt[i,:], color=opt_colors[i], label="Оптимизированное y($i)", linestyle=:dash)
end
display(current())

We can find that at the new optimised parameters obtained with limited observations, substances $C_1$ and $C_2$ are consumed more slowly and substance $C_3$ accumulates to a lesser extent. However, substances $C_4$, $C_5$ and $C_6$ have exactly the same dynamics at both found and true parameters.

Conclusion

In this example, we solved the problem of finding the optimal parameters of an ODE describing a multistage chain of chemical reactions using the [problem-based optimisation] workflow (https://engee.com/community/catalogs/projects/algoritm-resheniia-zadach-optimizatsii). We modelled the situation under the condition of limited observations and compared the dynamics of the system.