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{(), Tuple{}}}(ModelingToolkit.ODESystem(0x0000000000000005, 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(0x0000000000000005, 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, optimization_variables, []
+)
The optimization_states contain several groups of variables, namely measurements, controls, initial_conditions, and regularization. measurements represent the decision to observe at a specific time point at the grid. We currently work with the naming convention w_i for the i-th observed equation.
Finally, we are now able to convert our OEDProblem into an OptimizationProblem and solve it.
Currently we only support AutoForwardDiff() as an AD backend.
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.150444691085876e-7, 8.163093694231291e-7, 8.189840660977134e-7, 8.233196800264876e-7, 8.296873643852418e-7, 8.386102058879344e-7, 8.508142555929592e-7, 8.673092593637045e-7, 8.895185223545978e-7, 9.194947024680233e-7 … 0.9999989892807954, 0.9999976003263097, 2.4075867388525625e-5, 2.641192014873092e-6, 1.5899207116454061e-6, 1.2374298383970999e-6, 1.0704250376395544e-6, 9.783547235904154e-7, 9.232085383202934e-7, 8.884586865311912e-7], w₂ = [6.082383847588166e-7, 6.08655813142258e-7, 6.094145554437966e-7, 6.104660817182134e-7, 6.117935440021537e-7, 6.134082806905058e-7, 6.153495754639441e-7, 6.1768770873747e-7, 6.205308857948444e-7, 6.24037330504233e-7 … 0.9999994678008942, 0.9999993693655769, 0.9999992146046348, 0.9999989479780069, 0.9999984062904621, 0.9999968096403159, 0.999965904022932, 4.081571854526682e-6, 2.049937558186124e-6, 1.4270493773830689e-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)
