DifferentialEquations.jl - part 2: decision from the model

In the last article, we explored different modeling options for a three-component systems which could represent the dynamics of a chemical reaction or a disease propagation in a population. Building on top of this model, we will formulate a desirable outcome and find a decision which maximizes this outcome.

In addition to the packages imported in the last post, we will also use BlackBoxOptim.jl:

import DifferentialEquations
const DiffEq = DifferentialEquations
import Plots
import BlackBoxOptim

The model

The same chemical system with three components, A, B and R will be used: $$A + B → 2B$$ $$B → R$$

The reactor where the reaction occurs must remain active for one minute. Let’s imagine that $B$ is our valuable component while $R$ is a waste. We want to maximize the quantity of $B$ present within the system after one minute, that’s the objective function. For that purpose, we can choose to add a certain quantity of new $A$ within the reactor at any point. $$t_{inject} ∈ [0,t_{final}]$$.

Implementing the injection

There is one major feature of DifferentialEquations.jl we haven’t explored yet: the event handling system. This allows for the system state to change at a particular point in time, depending on conditions on the time, state, etc…

# defining the problem
const α = 0.8
const β = 3.0
diffeq = function(t, u, du)
    du[1] = - α * u[1] * u[2]
    du[2] = α * u[1] * u[2] - β * u[2]
    du[3] = β * u[2]
u0 = [49.0;1.0;0.0]
tspan = (0.0, 1.0)
prob = DiffEq.ODEProblem(diffeq, u0, tspan)

const A_inj = 30
inject_new = function(t0)
    condition(t, u, integrator) = t0 - t
    affect! = function(integrator)
        integrator.u[1] = integrator.u[1] + A_inj
    callback = DiffEq.ContinuousCallback(condition, affect!)
    sol = DiffEq.solve(prob, callback=callback)

# trying it out with an injection at t=0.4
sol = inject_new(0.4)

Injection simulation

The ContinuousCallback construct is the central element here, it takes as information:

  • When to trigger the event, implemented as the condition function. It triggers when this function reaches 0, which is here the case when $t = t_0$.
  • What to do with the state at that moment. The state is encapsulated within the integrator variable. In our case, we add 30 units to the concentration in A.

As we can see on the plot, a discontinuity appears on the concentration in A at the injection time, the concentration in B restarts increasing.

Finding the optimal injection time: visual approach

From the previously built function, we can get the whole solution with a given injection time, and from that the final state of the system.

tinj_span = 0.05:0.005:0.95
final_b = [inject_new(tinj).u[end][2] for tinj in tinj_span]
Plots.plot(tinj_span, final_b)

Using a plain for comprehension, we fetch the solution of the simulation for the callback built with each $t_{inject}$.

Quantity of B

Injecting $A$ too soon lets too much time for the created $B$ to turn into $R$, but injecting it too late does not let enough time for $B$ to be produced from the injected $A$. The optimum seems to be around ≈ 0.82,

Finding the optimum using BlackBoxOptim.jl

The package requires an objective function which takes a vector as input. In our case, the decision is modeled as a single variable (the injection time), it’s crucial to make the objective use a vector nonetheless, otherwise calling the solver will just explode with cryptic errors.

compute_finalb = tinj -> -1 * inject_new(tinj[1]).u[end][2]
# trust the default algorithm
BlackBoxOptim.bboptimize(compute_finalb, SearchRange=(0.1,0.9), NumDimensions=1)
# use probabilistic descent
BlackBoxOptim.bboptimize(compute_finalb, SearchRange=(0.1,0.9), NumDimensions=1, Method=:probabilistic_descent)

The function inject_new we defined above returns the complete solution of the simulation, we get the state matrix u, from which we extract the final state u[end], and then the second component, the concentration in B: u[end][2]. The black box optimizer minimizes the objective, while we want to maximize the final concentration of B, hence the -1 multiplier used for

The bboptimize function can also be passed a Method argument specifying the optimization algorithm. In this case, the function is smooth, so we can suppose gradient estimation methods would work pretty well. We also let the default algorithm (differential evolution) be used. After some lines logging the progress on the search, we obtain the following for both methods:

Best candidate found: [0.835558]
Fitness: -24.039369448

More importantly, we can refer to the number of evaluations of the objective function as a measure of the algorithm performance, combined with the time taken.

For the probabilistic descent:

Optimization stopped after 10001 steps and 10.565439939498901 seconds
Steps per second = 946.5767689058794
Function evals per second = 1784.6866867802482
Improvements/step = 0.0
Total function evaluations = 18856

For the differential evolution:

Optimization stopped after 10001 steps and 5.897292137145996 seconds
Steps per second = 1695.863078751937
Function evals per second = 1712.8200138459472
Improvements/step = 0.1078
Total function evaluations = 10101

We found the best injection time (0.835558), and the corresponding final concentration (24.04).

Extending the model

The decision over one variable was pretty straightforward. We are going to extend it by changing how the $A$ component is added at $t_{inject}$. Instead of being completely dissolved, a part of the component will keep being poured in after $t_{inject}$. So the decision will be composed of two variables:

  • The time of the beginning of the injection
  • The part of $A$ to inject directly and the part to inject in a continuous fashion. We will note the fraction injected directly $\delta$.

Given a fixed available quantity $A₀$ and a fraction to inject directly $\delta$, the concentration in A is increased of $\delta \cdot A₀$ at time $t_{inject}$, after which the rate of change of the concentration in A is increased by a constant amount, until the total amount of A injected (directly and over time) is equal to the planned quantity.

We need a new variable in the state of the system, $u_4(t)$, which stands for the input flow of A being active or not.

  • $u(t) = 0$ if $t < t_{inject}$
  • $u(t) = 0$ if the total flow of A which has been injected is equal to the planned quantity
  • $u(t) = \dot{A}\ $ otherwise, with $\dot{A}\ $ the rate at which A is being poured.

New Julia equations

We already built the key components in the previous sections. This time we need two events:

  • A is directly injected at $t_{inject}$, and then starts being poured at constant rate
  • A stops being poured when the total quantity has been used
const inj_quantity = 30.0;
const inj_rate = 40.0;

diffeq_extended = function(t, u, du)
    du[1] = - α * u[1] * u[2] + u[4]
    du[2] = α * u[1] * u[2] - β * u[2]
    du[3] = β * u[2]
    du[4] = 0.0

u0 = [49.0;1.0;0.0;0.0]
tspan = (0.0, 1.0)
prob = DiffEq.ODEProblem(diffeq_extended, u0, tspan)

We wrap the solution building process into a function taking the starting time and the fraction being directly injected as parameters:

inject_progressive = function(t0, direct_frac)
    condition_start(t, u, integrator) = t0 - t
    affect_start! = function(integrator)
        integrator.u[1] = integrator.u[1] + inj_quantity * direct_frac
        integrator.u[4] = inj_rate
    callback_start = DiffEq.ContinuousCallback(
        condition_start, affect_start!, save_positions=(true, true)
    condition_end(t, u, integrator) = (t - t0) * inj_rate - inj_quantity * (1 - direct_frac)
    affect_end! = function(integrator)
        integrator.u[4] = 0.0
    callback_end = DiffEq.ContinuousCallback(condition_end, affect_end!, save_positions=(true, true))
    sol = DiffEq.solve(prob, callback=DiffEq.CallbackSet(callback_start, callback_end), dtmax=0.005)


We can notice callback_start being identical to the model we previously built, while condition_end corresponds to the time when the total injected quantity reaches inj_quantity. The first events activates $u_4$ and sets it to the nominal flow, while the second callback resets it to 0.

Constant rate

BlackBoxOptim.jl can be re-used to determine the optimal decision:

objective = function(x)
    sol = inject_progressive(x[1], x[2])
BlackBoxOptim.bboptimize(objective, SearchRange=[(0.1,0.9),(0.0,1.0)], NumDimensions=2)

The optimal solution corresponds to a complete direct injection ($\delta \approx 1$) with $t_{inject}^{opt}$ identical to the previous model. This means pouring the A component in a continuous fashion does not allow to produce more $B$ at the end of the minute.


We could still built on top of this model to keep refining it, taking more phenomena into account (what if the reactions produce heat and are sensitive to temperature?). The structures describing models built with DifferentialEquations.jl are transparent and easy to use for further manipulations.

One point on which I place expectations is some additional interoperability between DifferentialEquations.jl and JuMP, a Julia meta-package for optimization. Some great work was already performed to combine the two systems, one use case that has been described is the parameter identification problem (given the evolution of concentration in the system, identify the α and β parameters).

But given that the function I built from a parameter was a black box (without an explicit formula, not a gradient), I had to use BlackBoxOptim, which is amazingly straightforward, but feels a bit overkill for smooth functions as presented here. Maybe there is a different way to build the objective function, using parametrized functions for instance, which could make it transparent to optimization solvers.

If somebody has info on that last point or feedback, additional info you’d like to share regarding this post, hit me on Twitter. Thanks for reading!

Edits and improvements

Of course, BlackBoxOptim.jl was not the most appropriate algorithm as predicted. Patrick and Chris gave me some hints in this thread and I gave Optim.jl a try.

This package has a range of algorithms to choose from depending on the structure of the function and the knowledge of its gradient and Hessian. The goal is continuous optimization, (as opposed to BlackBoxOptim.jl which supports more exotic search spaces).

Finding the optimum $t_{inject}$ of the first problem is pretty simple:

import Optim
Optim.optimize(compute_finalb, 0.1, 0.9)

This yields the following information:

Results of Optimization Algorithm
 * Algorithm: Brent's Method
 * Search Interval: [0.100000, 0.900000]
 * Minimizer: 8.355891e-01
 * Minimum: -2.403824e+01
 * Iterations: 13
 * Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(1.5e-08*|x|+2.2e-16): true
 * Objective Function Calls: 14

14 calls to the objective function, pretty neat compared to the hundreds of BlackBoxOptim. We also confirm the optimum of 0.8355891.

[1] Cover image: Lorenz attractor on Wikimedia, again.