Design of Experiments for Lotka-Volterra

Now, let us consider the well-known Lotka-Volterra system

\[\begin{aligned} \dot{x_1} &= p_1 x - p_2 x y\\  \dot{x_2} &= - p_3 y + p_4 x y\\ y(x) &= x \end{aligned} \]

where we are interested in estimating the parameters $p_2$ and $p_4$. We can measure the states directly with some fixed measurement rate.

We start by using ModelingToolkit.jl and DynamicOED.jl to model the system.

using DynamicOED
using ModelingToolkit
using Optimization, OptimizationMOI, Ipopt
using Plots

@variables t
@variables x(t)=0.5 [description = "Biomass Prey"] y(t)=0.7 [description ="Biomass Predator"]
@variables u(t) [description = "Control"]
@parameters p[1:2]=[1.0; 1.0] [description = "Fixed Parameters", tunable = false]
@parameters p_est[1:2]=[1.0; 1.0] [description = "Tunable Parameters", tunable = true]
D = Differential(t)
@variables obs(t)[1:2] [description = "Observed", measurement_rate = 96]
obs = collect(obs)

@named lotka_volterra = ODESystem(
    [
        D(x) ~   p[1]*x - p_est[1]*x*y;
        D(y) ~  -p[2]*y + p_est[2]*x*y
    ], tspan = (0.0, 12.0),
    observed = obs .~ [x; y]
)

\[ \begin{align} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& p_1 x\left( t \right) - p_{est_1} x\left( t \right) y\left( t \right) \\ \frac{\mathrm{d} y\left( t \right)}{\mathrm{d}t} =& - p_2 y\left( t \right) + p_{est_2} x\left( t \right) y\left( t \right) \end{align} \]

Like in the Design of Experiments for a simple system, we added important information:

• The observed variables are initialized with a measurement_rate. This time we use an integer measurement rate, resulting in $96$ subintervals of equal length.
• The parameters $p_2$ and $p_4$ are marked as tunable.

Now we can augment the system with the needed expressions for the sensitivities and the Fisher Information matrix by constructing an OEDSystem.

@named oed_system = OEDSystem(lotka_volterra)

\[ \begin{align} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& p_1 x\left( t \right) - p_{est_1} x\left( t \right) y\left( t \right) \\ \frac{\mathrm{d} y\left( t \right)}{\mathrm{d}t} =& - p_2 y\left( t \right) + p_{est_2} x\left( t \right) y\left( t \right) \\ 0 =& - \frac{\mathrm{d}}{\mathrm{d}t} G(t)_{1}ˏ_1 - x\left( t \right) y\left( t \right) + G(t)_{1}ˏ_1 \left( p_1 - p_{est_1} y\left( t \right) \right) - G(t)_{2}ˏ_1 p_{est_1} x\left( t \right) \\ 0 =& - \frac{\mathrm{d}}{\mathrm{d}t} G(t)_{2}ˏ_1 + G(t)_{1}ˏ_1 p_{est_2} y\left( t \right) + G(t)_{2}ˏ_1 \left( - p_2 + p_{est_2} x\left( t \right) \right) \\ 0 =& - \frac{\mathrm{d}}{\mathrm{d}t} G(t)_{1}ˏ_2 + G(t)_{1}ˏ_2 \left( p_1 - p_{est_1} y\left( t \right) \right) - G(t)_{2}ˏ_2 p_{est_1} x\left( t \right) \\ 0 =& - \frac{\mathrm{d}}{\mathrm{d}t} G(t)_{2}ˏ_2 + x\left( t \right) y\left( t \right) + G(t)_{1}ˏ_2 p_{est_2} y\left( t \right) + G(t)_{2}ˏ_2 \left( - p_2 + p_{est_2} x\left( t \right) \right) \\ \frac{\mathrm{d}}{\mathrm{d}t} F(t)_{1}ˏ_1 =& G(t)_{1}ˏ_1^{2} w_1 + G(t)_{2}ˏ_1^{2} w_2 \\ \frac{\mathrm{d}}{\mathrm{d}t} F(t)_{1}ˏ_2 =& G(t)_{1}ˏ_1 G(t)_{1}ˏ_2 w_1 + G(t)_{2}ˏ_1 G(t)_{2}ˏ_2 w_2 \\ \frac{\mathrm{d}}{\mathrm{d}t} F(t)_{2}ˏ_2 =& G(t)_{1}ˏ_2^{2} w_1 + G(t)_{2}ˏ_2^{2} w_2 \\ Q(t)_{1}ˏ_1 =& G(t)_{1}ˏ_1 \\ Q(t)_{2}ˏ_1 =& G(t)_{2}ˏ_1 \\ Q(t)_{1}ˏ_2 =& G(t)_{1}ˏ_2 \\ Q(t)_{2}ˏ_2 =& G(t)_{2}ˏ_2 \end{align} \]

With this augmented ODESystem we can set up the OEDProblem by specifying the criterion we want to optimize.

oed_problem = OEDProblem(structural_simplify(oed_system), DCriterion())

OEDProblem{ModelingToolkit.ODESystem, DCriterion, DynamicOED.Timegrid{Tuple{Symbol, Symbol}, Matrix{Int64}, Tuple{Vector{Tuple{Float64, Float64}}, Vector{Tuple{Float64, Float64}}}, Vector{Tuple{Float64, Float64}}}, OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, @NamedTuple{}}(ModelingToolkit.ODESystem(0x0000000000000004, Symbolics.Equation[Differential(t)(x(t)) ~ p[1]*x(t) - p_est[1]*x(t)*y(t), Differential(t)(y(t)) ~ -p[2]*y(t) + p_est[2]*x(t)*y(t), Differential(t)((G(t))[1, 1]) ~ -x(t)*y(t) + (G(t))[1, 1]*(p[1] - p_est[1]*y(t)) - (G(t))[2, 1]*p_est[1]*x(t), Differential(t)((G(t))[2, 1]) ~ (G(t))[1, 1]*p_est[2]*y(t) + (G(t))[2, 1]*(-p[2] + p_est[2]*x(t)), Differential(t)((G(t))[1, 2]) ~ (G(t))[1, 2]*(p[1] - p_est[1]*y(t)) - (G(t))[2, 2]*p_est[1]*x(t), Differential(t)((G(t))[2, 2]) ~ x(t)*y(t) + (G(t))[1, 2]*p_est[2]*y(t) + (G(t))[2, 2]*(-p[2] + p_est[2]*x(t)), Differential(t)((F(t))[1, 1]) ~ ((G(t))[1, 1]^2)*w₁ + ((G(t))[2, 1]^2)*w₂, Differential(t)((F(t))[1, 2]) ~ (G(t))[1, 1]*(G(t))[1, 2]*w₁ + (G(t))[2, 1]*(G(t))[2, 2]*w₂, Differential(t)((F(t))[2, 2]) ~ ((G(t))[1, 2]^2)*w₁ + ((G(t))[2, 2]^2)*w₂], t, Any[x(t), y(t), (G(t))[1, 1], (G(t))[2, 1], (G(t))[1, 2], (G(t))[2, 2], (F(t))[1, 1], (F(t))[1, 2], (F(t))[2, 2]], SymbolicUtils.BasicSymbolic[p_est[1], p[1], p_est[2], p[2], w₁, w₂], (0.0, 12.0), Dict{Any, Any}(:F => F(t), :p => p, :y => y(t), :G => G(t), :Q => Q(t), :w₂ => w₂, :obs => obs(t), :p_est => p_est, :w₁ => w₁, :x => x(t)…), Any[], Symbolics.Equation[(obs(t))[1] ~ x(t), (obs(t))[2] ~ y(t), (Q(t))[1, 1] ~ (G(t))[1, 1], (Q(t))[2, 1] ~ (G(t))[2, 1], (Q(t))[1, 2] ~ (G(t))[1, 2], (Q(t))[2, 2] ~ (G(t))[2, 2]], Base.RefValue{Vector{Symbolics.Num}}(Symbolics.Num[]), Base.RefValue{Any}(Matrix{Symbolics.Num}(undef, 0, 0)), Base.RefValue{Any}(Matrix{Symbolics.Num}(undef, 0, 0)), Base.RefValue{Matrix{Symbolics.Num}}(Matrix{Symbolics.Num}(undef, 0, 0)), Base.RefValue{Matrix{Symbolics.Num}}(Matrix{Symbolics.Num}(undef, 0, 0)), :oed_system, ModelingToolkit.ODESystem[], Dict{Any, Any}((G(t))[1, 2] => 0.0, (F(t))[1, 1] => 0.0, p_est[1] => 1.0, (F(t))[2, 2] => 0.0, (G(t))[1, 1] => 0.0, p[2] => 1.0, p_est[2] => 1.0, (G(t))[2, 2] => 0.0, x(t) => 0.5, w₂ => 0.0…), nothing, nothing, nothing, ModelingToolkit.SymbolicContinuousCallback[ModelingToolkit.SymbolicContinuousCallback(Symbolics.Equation[], Symbolics.Equation[])], ModelingToolkit.SymbolicDiscreteCallback[], nothing, nothing, TearingState of ModelingToolkit.ODESystem, ModelingToolkit.Substitutions(Symbolics.Equation[(Q(t))[1, 1] ~ (G(t))[1, 1], (Q(t))[2, 1] ~ (G(t))[2, 1], (Q(t))[1, 2] ~ (G(t))[1, 2], (Q(t))[2, 2] ~ (G(t))[2, 2]], [Int64[], [1], [1, 2], [1, 2, 3]], nothing), true, nothing, nothing, nothing, ModelingToolkit.ODESystem(0x0000000000000004, Symbolics.Equation[Differential(t)(x(t)) ~ p[1]*x(t) - p_est[1]*x(t)*y(t), Differential(t)(y(t)) ~ -p[2]*y(t) + p_est[2]*x(t)*y(t), 0.0 ~ -Differential(t)((G(t))[1, 1]) - x(t)*y(t) + (G(t))[1, 1]*(p[1] - p_est[1]*y(t)) - (G(t))[2, 1]*p_est[1]*x(t), 0.0 ~ -Differential(t)((G(t))[2, 1]) + (G(t))[1, 1]*p_est[2]*y(t) + (G(t))[2, 1]*(-p[2] + p_est[2]*x(t)), 0.0 ~ -Differential(t)((G(t))[1, 2]) + (G(t))[1, 2]*(p[1] - p_est[1]*y(t)) - (G(t))[2, 2]*p_est[1]*x(t), 0.0 ~ -Differential(t)((G(t))[2, 2]) + x(t)*y(t) + (G(t))[1, 2]*p_est[2]*y(t) + (G(t))[2, 2]*(-p[2] + p_est[2]*x(t)), Differential(t)((F(t))[1, 1]) ~ ((G(t))[1, 1]^2)*w₁ + ((G(t))[2, 1]^2)*w₂, Differential(t)((F(t))[1, 2]) ~ (G(t))[1, 1]*(G(t))[1, 2]*w₁ + (G(t))[2, 1]*(G(t))[2, 2]*w₂, Differential(t)((F(t))[2, 2]) ~ ((G(t))[1, 2]^2)*w₁ + ((G(t))[2, 2]^2)*w₂, (Q(t))[1, 1] ~ (G(t))[1, 1], (Q(t))[2, 1] ~ (G(t))[2, 1], (Q(t))[1, 2] ~ (G(t))[1, 2], (Q(t))[2, 2] ~ (G(t))[2, 2]], t, SymbolicUtils.BasicSymbolic{Real}[x(t), y(t), (G(t))[1, 1], (G(t))[2, 1], (G(t))[1, 2], (G(t))[2, 2], (F(t))[1, 1], (F(t))[1, 2], (F(t))[2, 2], (Q(t))[1, 1], (Q(t))[2, 1], (Q(t))[1, 2], (Q(t))[2, 2]], SymbolicUtils.BasicSymbolic[p_est[1], p[1], p_est[2], p[2], w₁, w₂], (0.0, 12.0), Dict{Any, Any}(:F => F(t), :p => p, :y => y(t), :G => G(t), :Q => Q(t), :w₂ => w₂, :obs => obs(t), :p_est => p_est, :w₁ => w₁, :x => x(t)…), Any[], Symbolics.Equation[(obs(t))[1] ~ x(t), (obs(t))[2] ~ y(t)], Base.RefValue{Vector{Symbolics.Num}}(Symbolics.Num[]), Base.RefValue{Any}(Matrix{Symbolics.Num}(undef, 0, 0)), Base.RefValue{Any}(Matrix{Symbolics.Num}(undef, 0, 0)), Base.RefValue{Matrix{Symbolics.Num}}(Matrix{Symbolics.Num}(undef, 0, 0)), Base.RefValue{Matrix{Symbolics.Num}}(Matrix{Symbolics.Num}(undef, 0, 0)), :oed_system, ModelingToolkit.ODESystem[], Dict{Any, Any}((G(t))[1, 2] => 0.0, (F(t))[1, 1] => 0.0, p[2] => 1.0, (G(t))[2, 2] => 0.0, (F(t))[1, 2] => 0.0, p[1] => 1.0, w₁ => 0.0, y(t) => 0.7, p_est[1] => 1.0, (F(t))[2, 2] => 0.0…), nothing, nothing, nothing, ModelingToolkit.SymbolicContinuousCallback[ModelingToolkit.SymbolicContinuousCallback(Symbolics.Equation[], Symbolics.Equation[])], ModelingToolkit.SymbolicDiscreteCallback[], nothing, nothing, nothing, nothing, true, nothing, nothing, nothing, nothing)), DCriterion(), DynamicOED.Timegrid{Tuple{Symbol, Symbol}, Matrix{Int64}, Tuple{Vector{Tuple{Float64, Float64}}, Vector{Tuple{Float64, Float64}}}, Vector{Tuple{Float64, Float64}}}((:w₁, :w₂), [1 2 … 95 96; 1 2 … 95 96], ([(0.0, 0.125), (0.125, 0.25), (0.25, 0.375), (0.375, 0.5), (0.5, 0.625), (0.625, 0.75), (0.75, 0.875), (0.875, 1.0), (1.0, 1.125), (1.125, 1.25)  …  (10.75, 10.875), (10.875, 11.0), (11.0, 11.125), (11.125, 11.25), (11.25, 11.375), (11.375, 11.5), (11.5, 11.625), (11.625, 11.75), (11.75, 11.875), (11.875, 12.0)], [(0.0, 0.125), (0.125, 0.25), (0.25, 0.375), (0.375, 0.5), (0.5, 0.625), (0.625, 0.75), (0.75, 0.875), (0.875, 1.0), (1.0, 1.125), (1.125, 1.25)  …  (10.75, 10.875), (10.875, 11.0), (11.0, 11.125), (11.125, 11.25), (11.25, 11.375), (11.375, 11.5), (11.5, 11.625), (11.625, 11.75), (11.75, 11.875), (11.875, 12.0)]), [(0.0, 0.125), (0.125, 0.25), (0.25, 0.375), (0.375, 0.5), (0.5, 0.625), (0.625, 0.75), (0.75, 0.875), (0.875, 1.0), (1.0, 1.125), (1.125, 1.25)  …  (10.75, 10.875), (10.875, 11.0), (11.0, 11.125), (11.125, 11.25), (11.25, 11.375), (11.375, 11.5), (11.5, 11.625), (11.625, 11.75), (11.75, 11.875), (11.875, 12.0)]), Tsit5(; stage_limiter! = trivial_limiter!, step_limiter! = trivial_limiter!, thread = static(false),), NamedTuple())

We choose the DCriterion, which minimizes the determinant of the inverse of the Fisher Information matrix. For constraining the time we can measure by defining a ConstraintSystem from ModelingToolkit on the optimization variables. We want to measure for at most $4$ units of time. Since we discretized the observed variables on $96$ subintervals on a time horizon of $12$ units of time, this translates to an upper limit on the measurements of $32$.

optimization_variables = states(oed_problem)

constraint_equations = [
    sum(optimization_variables.measurements.w₁) ≲ 32,
    sum(optimization_variables.measurements.w₂) ≲ 32,
]

@named constraint_system = ConstraintsSystem(
    constraint_equations, collect(optimization_variables), Num[]
)

Finally, we are now able to convert our OEDProblem into an OptimizationProblem and solve it.

optimization_problem = OptimizationProblem(
    oed_problem, AutoForwardDiff(), constraints = constraint_system,
    integer_constraints = false
)

optimal_design = solve(optimization_problem, Ipopt.Optimizer(); hessian_approximation="limited-memory")

u_opt = optimal_design.u + optimization_problem.u0

ComponentVector{Float64}(controls = Float64[], initial_conditions = Float64[], measurements = (w₁ = [8.150444691091382e-7, 8.163093694236809e-7, 8.189840660982695e-7, 8.233196800270488e-7, 8.296873643858106e-7, 8.386102058885147e-7, 8.508142555935555e-7, 8.673092593643241e-7, 8.895185223552471e-7, 9.194947024687153e-7  …  0.9999989892807954, 0.9999976003263097, 2.4075867388679947e-5, 2.6411920148780405e-6, 1.5899207116473755e-6, 1.23742983839834e-6, 1.0704250376404709e-6, 9.783547235911938e-7, 9.232085383209892e-7, 8.884586865318383e-7], w₂ = [6.082383847585814e-7, 6.086558131420219e-7, 6.094145554435602e-7, 6.104660817179761e-7, 6.11793544001915e-7, 6.134082806902657e-7, 6.153495754637033e-7, 6.176877087372267e-7, 6.20530885794599e-7, 6.240373305039844e-7  …  0.9999994678008942, 0.9999993693655769, 0.9999992146046348, 0.9999989479780069, 0.9999984062904621, 0.9999968096403158, 0.9999659040229282, 4.081571854515923e-6, 2.049937558183312e-6, 1.4270493773817238e-6]), regularization = 0.9999999900100008)

Now we want to visualize the found solution.

function plotoed(problem, res)

    predictor = DynamicOED.build_predictor(problem)
    x_opt, t_opt = predictor(res)
    timegrid = problem.timegrid

    state_plot = plot(t_opt, x_opt[1:2, :]', xlabel = "Time", ylabel = "States", label = ["x" "y"])

    measures_plot = plot()
    for i in 1:2
        t_measures = vcat(first.(timegrid.timegrids[i]), last.(timegrid.timegrids[i]))
        sort!(t_measures)
        unique!(t_measures)
        _measurements = getfield(res.measurements |> NamedTuple, timegrid.variables[i])
        plot!(t_measures,
            vcat(_measurements, last(_measurements)),
            line = :steppost,
            xlabel = "Time",
            ylabel = "Measurement",
            color = i == 2 ? :red : :blue,
            label = string(timegrid.variables[i]))
    end

    plot(state_plot, measures_plot, layout=(2,1))
end

plotoed(oed_problem, u_opt)