Design of Experiments for a simple system

As a pedagogical example, lets start with a simple system.

\[\begin{aligned} \dot{x} &= p_1 x \\  y &= x \end{aligned} \]

where the state $x$ is dependent on a single, uncertain parameter $p_1$. We assume we can observe the state 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

@variables t x(t) = 1.0 [description = "State"]
@variables y(t) [description="Observations", measurement_rate = 0.1]
@parameters p = -2.0 [description= "Uncertain parameter", tunable = true]

D = Differential(t)

@named simple_system = ODESystem(
    [
        D(x) ~ p*x
    ], tspan = (0., 1.),
    observed = [
        y  ~ x
    ]
)

\[ \begin{align} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& p x\left( t \right) \end{align} \]

Note that in order to prepare the system for optimal experimental design, we added two important information:

• The parameter $p_1$ is marked as tunable.
• The observed variable $y$ has a measurement_rate.

This way, we can limit the set of tunable parameters of the system and derive time grids for the observations independently. To derive the associated system for experimental design, we simply construct an OEDSystem.

@named oed_system = OEDSystem(simple_system)

\[ \begin{align} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& p x\left( t \right) \\ 0 =& x\left( t \right) - \frac{\mathrm{d}}{\mathrm{d}t} G(t)_{1}ˏ_1 + G(t)_{1}ˏ_1 p \\ \frac{\mathrm{d}}{\mathrm{d}t} F(t)_{1}ˏ_1 =& G(t)_{1}ˏ_1^{2} w_1 \\ Q(t)_{1}ˏ_1 =& G(t)_{1}ˏ_1 \end{align} \]

This step augments the initial system with all necessary information, e.g. the sensitivity equations and the dynamics of the fisher information matrix. The output is simply another ODESystem, hence we can use available transformations and simplification routines from ModelingToolkit.

Next, we associate an optimization criterion to the system and instantiate it over a specific time grid.

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

OEDProblem{ModelingToolkit.ODESystem, DCriterion, DynamicOED.Timegrid{Tuple{Symbol}, Matrix{Int64}, Tuple{Vector{Tuple{Float64, Float64}}}, Vector{Tuple{Float64, Float64}}}, OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, @NamedTuple{}}(ModelingToolkit.ODESystem(0x0000000000000001, Symbolics.Equation[Differential(t)(x(t)) ~ p*x(t), Differential(t)((G(t))[1, 1]) ~ x(t) + (G(t))[1, 1]*p, Differential(t)((F(t))[1, 1]) ~ ((G(t))[1, 1]^2)*w₁], t, Any[x(t), (G(t))[1, 1], (F(t))[1, 1]], SymbolicUtils.BasicSymbolic[p, w₁], (0.0, 1.0), Dict{Any, Any}(:F => F(t), :p => p, :G => G(t), :y => y(t), :Q => Q(t), :w₁ => w₁, :x => x(t)), Any[], Symbolics.Equation[y(t) ~ x(t), (Q(t))[1, 1] ~ (G(t))[1, 1]], 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}(p => -2.0, x(t) => 1.0, (F(t))[1, 1] => 0.0, w₁ => 0.0, (G(t))[1, 1] => 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]], [Int64[]], nothing), true, nothing, nothing, nothing, ModelingToolkit.ODESystem(0x0000000000000001, Symbolics.Equation[Differential(t)(x(t)) ~ p*x(t), 0.0 ~ x(t) - Differential(t)((G(t))[1, 1]) + (G(t))[1, 1]*p, Differential(t)((F(t))[1, 1]) ~ ((G(t))[1, 1]^2)*w₁, (Q(t))[1, 1] ~ (G(t))[1, 1]], t, SymbolicUtils.BasicSymbolic{Real}[x(t), (G(t))[1, 1], (F(t))[1, 1], (Q(t))[1, 1]], SymbolicUtils.BasicSymbolic[p, w₁], (0.0, 1.0), Dict{Any, Any}(:F => F(t), :p => p, :G => G(t), :y => y(t), :Q => Q(t), :w₁ => w₁, :x => x(t)), Any[], Symbolics.Equation[y(t) ~ x(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}(p => -2.0, x(t) => 1.0, (F(t))[1, 1] => 0.0, w₁ => 0.0, (G(t))[1, 1] => 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}, Matrix{Int64}, Tuple{Vector{Tuple{Float64, Float64}}}, Vector{Tuple{Float64, Float64}}}((:w₁,), [1 2 … 9 10], ([(0.0, 0.1), (0.1, 0.2), (0.2, 0.3), (0.3, 0.4), (0.4, 0.5), (0.5, 0.6), (0.6, 0.7), (0.7, 0.8), (0.8, 0.9), (0.9, 1.0)],), [(0.0, 0.1), (0.1, 0.2), (0.2, 0.3), (0.3, 0.4), (0.4, 0.5), (0.5, 0.6), (0.6, 0.7), (0.7, 0.8), (0.8, 0.9), (0.9, 1.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. Given that the optimal solution would be to measure all occurrences, we also need to add constraints to the measurement function.

optimization_variables = states(oed_problem)

\[ \begin{equation} \left[ \begin{array}{c} w_{1 1} \\ w_{1 2} \\ w_{1 3} \\ w_{1 4} \\ w_{1 5} \\ w_{1 6} \\ w_{1 7} \\ w_{1 8} \\ w_{1 9} \\ w_{1 1 0} \\ \tau \\ \end{array} \right] \end{equation} \]

Now we have access to all optimization variables as a ComponentArray. We are interested in choosing at most 3 measurements, so we add a ConstraintsSystem from ModelingToolkit.

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

@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: 11-element Vector{Float64}:
 -6.442615705818314e-9
 -4.519395856235989e-9
  2.9261803704474444e-9
  1.5958188117575518e-7
  0.9999999732471323
  0.9999999743442556
  0.9999998731278948
  2.9335593476398022e-8
  5.497538495725719e-9
 -3.6755641831694274e-10
  0.9901904800443652

We see that we have indeed recovered a bang-bang solution to our problem, even though we only solve a relaxed problem. The interested reader can find a more in depth explanation here.