Merge pull request #332 from SciML/dds_mtk

Adding MTK interface for the discrete-time identifiability

docs/src/assets/Project.toml

@@ -7,5 +7,5 @@ StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"
 [compat]
 BenchmarkTools = "1.3"
 Documenter = "0.27, 1"
-ModelingToolkit = "8.74"
+ModelingToolkit = "9"
 StructuralIdentifiability = "0.5"

docs/src/tutorials/discrete_time.md

@@ -7,14 +7,7 @@ $\begin{cases}
 \mathbf{y}(t) = \mathbf{g}(\mathbf{x}(t), \mathbf{p}, \mathbf{u(t)}),
 \end{cases}$
 
-or in terms of **difference**
-
-$\begin{cases}
-\Delta\mathbf{x}(t) = \mathbf{f}(\mathbf{x}(t), \mathbf{p}, \mathbf{u}(t)),\\
-\mathbf{y}(t) = \mathbf{g}(\mathbf{x}(t), \mathbf{p}, \mathbf{u(t)}),
-\end{cases} \quad \text{where} \quad \Delta \mathbf{x}(t) := \mathbf{x}(t + 1) - \mathbf{x}(t).$
-
-In both cases,$\mathbf{x}(t), \mathbf{y}(t)$, and $\mathbf{u}(t)$ are time-dependent states, outputs, and inputs, respectively,
+where $\mathbf{x}(t), \mathbf{y}(t)$, and $\mathbf{u}(t)$ are time-dependent states, outputs, and inputs, respectively,
 and $\mathbf{p}$ are scalar parameters.
 As in the ODE case, we will call that a parameter or a states (or a function of them) is **identifiable** if its value can be recovered from
 time series for inputs and outputs (in the generic case, see Definition 3 in [^1] for details).
@@ -34,14 +27,6 @@ S(t + 1) = S(t) - \beta S(t) I(t),\\
 I(t + 1) = I(t) + \beta S(t) I(t) - \alpha I(t),\\
 R(t + 1) = R(t) + \alpha I(t),\\
 y(t) = I(t),
-\end{cases}
-\quad \text{or}
-\quad
-\begin{cases}
-\Delta S(t) = -\beta S(t) I(t),\\
-\Delta I(t) = \beta S(t) I(t) - \alpha I(t),\\
-\Delta R(t) = \alpha I(t),\\
-y(t) = I(t),
 \end{cases}$
 
 where the observable is `I`, the number of infected people.
@@ -87,6 +72,28 @@ assess_local_identifiability(dds; funcs_to_check = [β * S])
 As other main functions in the package, `assess_local_identifiability` accepts an optional parameter `loglevel` (default: `Logging.Info`)
 to adjust the verbosity of logging.
 
+Another way to defined a discrete-time system and its assess identifiability is to use the [`DiscreteSystem`](https://docs.sciml.ai/ModelingToolkit/dev/tutorials/discrete_system/) object from `ModelingToolkit`.
+The following code will perform the same computation as above (note that `ModelingToolkit` requires the shifts to be negative):
+
+```@example mtk
+using ModelingToolkit
+using StructuralIdentifiability
+
+@parameters α β
+@variables t S(t) I(t) R(t)
+k = ShiftIndex(t)
+
+eqs = [
+    S(k) ~ S(k - 1) - β * S(k - 1) * I(k - 1),
+    I(k) ~ I(k - 1) + β * S(k - 1) * I(k - 1) - α * I(k - 1),
+    R(k) ~ R(k - 1) + α * I(k - 1),
+]
+
+@mtkbuild sys = DiscreteSystem(eqs, t)
+
+assess_local_identifiability(sys, measured_quantities = [I])
+```
+
 The implementation is based on a version of the observability rank criterion and will be described in a forthcoming paper.
 
 [^1]: > S. Nõmm, C. Moog, [*Identifiability of discrete-time nonlinear systems*](https://doi.org/10.1016/S1474-6670(17)31245-4), IFAC Proceedings Volumes, 2004.

ext/ModelingToolkitSIExt.jl

@@ -159,7 +159,7 @@ end
 
 Input:
 - `de` - ModelingToolkit.AbstractTimeDependentSystem, a system for identifiability query
-- `measured_quantities` - array of input function in the form (name, expression)
+- `measured_quantities` - array of output function in the form (name, expression)
 
 Output:
 - `ODE` object containing required data for identifiability assessment
@@ -432,7 +432,7 @@ end
         prob_threshold::Float64=0.99)
 
 Input:
-- `dds` - the DiscreteSystem object from ModelingToolkit (with **difference** operator in the left-hand side)
+- `dds` - the DiscreteSystem object from ModelingToolkit
 - `measured_quantities` - the measurable outputs of the model
 - `funcs_to_check` - functions of parameters for which to check identifiability (all parameters and states if not specified)
 - `known_ic` - functions (of states and parameter) whose initial conditions are assumed to be known
@@ -443,7 +443,6 @@ Output:
 
 The result is correct with probability at least `prob_threshold`.
 """
-#=
 function StructuralIdentifiability.assess_local_identifiability(
     dds::ModelingToolkit.DiscreteSystem;
     measured_quantities = Array{ModelingToolkit.Equation}[],
@@ -490,13 +489,8 @@ function _assess_local_identifiability(
     # Converting the finite difference operator in the right-hand side to
     # the corresponding shift operator
     eqs = filter(eq -> !(ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(dds))
-    deltas = [Symbolics.operation(e.lhs).dt for e in eqs]
-    @assert length(Set(deltas)) == 1
-    eqs_shift = [e.lhs ~ e.rhs + first(Symbolics.arguments(e.lhs)) for e in eqs]
-    dds_shift = DiscreteSystem(eqs_shift, name = gensym())
-    @debug "System transformed from difference to shift: $dds_shift"
 
-    dds_aux_ode, conversion = mtk_to_si(dds_shift, measured_quantities)
+    dds_aux_ode, conversion = mtk_to_si(dds, measured_quantities)
     dds_aux = StructuralIdentifiability.DDS{QQMPolyRingElem}(dds_aux_ode)
     if length(funcs_to_check) == 0
         params = parameters(dds)
@@ -505,7 +499,7 @@ function _assess_local_identifiability(
         )
         funcs_to_check = vcat(
             [
-                x for x in states(dds) if
+                x for x in unknowns(dds) if
                 conversion[x] in StructuralIdentifiability.x_vars(dds_aux)
             ],
             union(params, params_from_measured_quantities),
@@ -529,7 +523,6 @@ function _assess_local_identifiability(
     end
     return out_dict
 end
-=#
 
 # ------------------------------------------------------------------------------
 

test/extensions/modelingtoolkit.jl

@@ -428,16 +428,19 @@ if GROUP == "All" || GROUP == "ModelingToolkitSIExt"
               correct
     end
 
-    #=
     @testset "Discrete local identifiability, ModelingToolkit interface" begin
         cases = []
 
         @parameters α β
         @variables t S(t) I(t) R(t) y(t)
-        D = Difference(t; dt = 1.0)
+        k = ShiftIndex(t)
 
-        eqs = [D(S) ~ -β * S * I, D(I) ~ β * S * I - α * I, D(R) ~ α * I]
-        @named sir = DiscreteSystem(eqs)
+        eqs = [
+            S(k) ~ S(k - 1) - β * S(k - 1) * I(k - 1),
+            I(k) ~ I(k - 1) + β * S(k - 1) * I(k - 1) - α * I(k - 1),
+            R(k) ~ R(k - 1) + α * I(k - 1),
+        ]
+        @mtkbuild sir = DiscreteSystem(eqs, t)
         push!(
             cases,
             Dict(
@@ -451,12 +454,11 @@ if GROUP == "All" || GROUP == "ModelingToolkitSIExt"
         )
 
         @parameters θ
-        @variables t x(t) y(t)
-        D = Difference(t; dt = 1.0)
+        @variables x(t) y(t)
 
-        eqs = [D(x) ~ θ * x^3 - x]
+        eqs = [x(k) ~ θ * x(k - 1)^3]
 
-        @named eqs = DiscreteSystem(eqs)
+        @mtkbuild eqs = DiscreteSystem(eqs, t)
         push!(
             cases,
             Dict(
@@ -470,12 +472,11 @@ if GROUP == "All" || GROUP == "ModelingToolkitSIExt"
         )
 
         @parameters θ β
-        @variables t x1(t) x2(t) y(t)
-        D = Difference(t; dt = 1.0)
+        @variables x1(t) x2(t) y(t)
 
-        eqs = [D(x1) ~ x1 + x2, D(x2) ~ θ + β]
+        eqs = [x1(k) ~ x1(k - 1) + x2(k - 1), x2(k) ~ x2(k - 1) + θ + β]
 
-        @named eqs = DiscreteSystem(eqs)
+        @mtkbuild eqs = DiscreteSystem(eqs, t)
         push!(
             cases,
             Dict(
@@ -489,12 +490,14 @@ if GROUP == "All" || GROUP == "ModelingToolkitSIExt"
         )
 
         @parameters a b c d
-        @variables t x1(t) x2(t) u(t) y2(t)
-        D = Difference(t; dt = 1.0)
+        @variables x1(t) x2(t) y2(t)
 
-        eqs = [D(x1) ~ a * x1 - b * x1 * x2 + u, D(x2) ~ -c * x2 + d * x1 * x2]
+        eqs = [
+            x1(k) ~ a * x1(k - 1) - b * x1(k - 1) * x2(k - 1),
+            x2(k) ~ -c * x2(k - 1) + d * x1(k - 1) * x2(k - 1),
+        ]
 
-        @named lv = DiscreteSystem(eqs)
+        @mtkbuild lv = DiscreteSystem(eqs, t)
         push!(
             cases,
             Dict(
@@ -565,13 +568,14 @@ if GROUP == "All" || GROUP == "ModelingToolkitSIExt"
         )
 
         # Example 1 from https://doi.org/10.1016/j.automatica.2008.03.019
+        #=
+        # Commented out because MTK does not seem to support inputs
         @parameters theta1 theta2
-        @variables t x1(t) x2(t) u(t) y(t)
-        D = Difference(t; dt = 1.0)
+        @variables x1(t) x2(t) u(t) y(t)
 
-        eqs = [D(x1) ~ theta1 * x1 + x2, D(x2) ~ (1 - theta2) * x1 + x2^2 + u - x2]
+        eqs = [x1(k) ~ theta1 * x1(k - 1) + x2(k - 1), x2(k) ~ (1 - theta2) * x1(k - 1) + x2(k - 1)^2 + u(k - 1)]
 
-        @named abmd1 = DiscreteSystem(eqs)
+        @mtkbuild abmd1 = DiscreteSystem(eqs, t)
         push!(
             cases,
             Dict(
@@ -583,15 +587,16 @@ if GROUP == "All" || GROUP == "ModelingToolkitSIExt"
                 :to_check => Array{}[],
             ),
         )
+        =#
 
         # Example 2 from https://doi.org/10.1016/j.automatica.2008.03.019
+        #=
         @parameters theta1 theta2 theta3
-        @variables t x1(t) x2(t) u(t) y(t) y2(t)
-        D = Difference(t; dt = 1.0)
+        @variables x1(t) x2(t) u(t) y(t) y2(t)
 
-        eqs = [D(x1) ~ theta1 * x1^2 + theta2 * x2 + u - x1, D(x2) ~ theta3 * x1 - x2]
+        eqs = [x1(k) ~ theta1 * x1(k - 1)^2 + theta2 * x2(k - 1) + u(k - 1), x2(k) ~ theta3 * x1(k - 1)]
 
-        @named abmd2 = DiscreteSystem(eqs)
+        @mtkbuild abmd2 = DiscreteSystem(eqs, t)
         push!(
             cases,
             Dict(
@@ -626,14 +631,13 @@ if GROUP == "All" || GROUP == "ModelingToolkitSIExt"
                 :to_check => Array{}[],
             ),
         )
-
+        =#
         @parameters a b
-        @variables t x1(t) y(t)
-        D = Difference(t; dt = 1.0)
+        @variables x1(t) y(t)
 
-        eqs = [D(x1) ~ a]
+        eqs = [x1(k) ~ x1(k - 1) + a]
 
-        @named kic = DiscreteSystem(eqs)
+        @mtkbuild kic = DiscreteSystem(eqs, t)
         push!(
             cases,
             Dict(
@@ -673,7 +677,7 @@ if GROUP == "All" || GROUP == "ModelingToolkitSIExt"
             ) == c[:res]
         end
     end
-    =#
+
     @testset "Exporting ModelingToolkit Model to SI Model" begin
 
         # Creates MTK model and assesses its identifiability.

test/local_identifiability_discrete.jl

@@ -23,8 +23,6 @@
 
     #---------------------
 
-    cases = []
-
     dds = @DDSmodel(x1(t + 1) = a * x1(t) * x2(t), x2(t + 1) = b * x1(t), y(t) = x2(t))
 
     push!(
@@ -48,6 +46,67 @@
 
     #---------------------
 
+    # Example 1 from https://doi.org/10.1016/j.automatica.2008.03.019
+    dds = @DDSmodel(
+        x1(t + 1) = theta1 * x1(t) + x2(t),
+        x2(t + 1) = (1 - theta2) * x1(t) + x2(t)^2 + u(t),
+        y(t) = x1(t)
+    )
+
+    push!(
+        cases,
+        Dict(
+            :dds => dds,
+            :res => OrderedDict(x1 => true, x2 => true, theta1 => true, theta2 => true),
+            :known => :none,
+        ),
+    )
+
+    # Example 2 from https://doi.org/10.1016/j.automatica.2008.03.019
+    dds = @DDSmodel(
+        x1(t + 1) = theta1 * x1(t)^2 + theta2 * x2(t) + u(t),
+        x2(t + 1) = theta3 * x1(t),
+        y(t) = x1(t)
+    )
+
+    push!(
+        cases,
+        Dict(
+            :dds => dds,
+            :res => OrderedDict(
+                x1 => true,
+                x2 => false,
+                theta1 => true,
+                theta2 => false,
+                theta3 => false,
+            ),
+            :known => :none,
+        ),
+    )
+
+    dds = @DDSmodel(
+        x1(t + 1) = theta1 * x1(t)^2 + theta2 * x2(t) + u(t),
+        x2(t + 1) = theta3 * x1(t),
+        y(t) = x1(t),
+        y2(t) = x2(t)
+    )
+    push!(
+        cases,
+        Dict(
+            :dds => dds,
+            :res => Dict(
+                x1 => true,
+                x2 => true,
+                theta1 => true,
+                theta2 => true,
+                theta3 => true,
+            ),
+            :known => :none,
+        ),
+    )
+
+    #---------------------
+
     for c in cases
         @test assess_local_identifiability(
             c[:dds];