Globally Identifiable Functions (a.k.a. Observation Field)

In addition to assessing identifiability of given functions of parameters or states, StructuralIdentifiability.jl provides the function find_identifiable_functions which finds a set of identifiable functions such that any other identifiable function can be expressed using them. This allows to find out what actually is identifiable and what non-identifiability in the model at hand looks like.

For example, consider the following model^[1].

LLW1987 = @ODEmodel(
    x1'(t) = -p1 * x1(t) + p2 * u(t),
    x2'(t) = -p3 * x2(t) + p4 * u(t),
    x3'(t) = -(p1 + p3) * x3(t) + (p4 * x1(t) + p2 * x2(t)) * u(t),
    y1(t) = x3(t)
)

x1'(t) = -x1(t)*p1 + u(t)*p2
x2'(t) = -x2(t)*p3 + u(t)*p4
x3'(t) = x1(t)*u(t)*p4 + x2(t)*u(t)*p2 - x3(t)*p1 - x3(t)*p3
y1(t) = x3(t)

Several decades ago, this model was introduced to demonstrate the presence of nontrivial unidentifiability in nonlinear systems of ODEs. Nowadays, we can automatically find the identifiable combinations of parameters:

find_identifiable_functions(LLW1987)

3-element Vector{AbstractAlgebra.Generic.FracFieldElem{Nemo.QQMPolyRingElem}}:
 p1 + p3
 p2*p4
 p1*p3

From these expressions, we see that the values of p1 and p3 are not identifiable but an unordered pair of numbers {p1, p3} is uniquely determined since p1 + p3 and p1 * p3 are known. Furthermore, we see that, for fixed input and output, p2 and p4 can take infinitely many values but knowing one of them, we would also be able to determine the other.

observation_field(LLW1987)

Moreover, we can find generators of all identifiable functions in parameters and states:

find_identifiable_functions(LLW1987, with_states = true)

7-element Vector{AbstractAlgebra.Generic.FracFieldElem{Nemo.QQMPolyRingElem}}:
 p1 + p3
 p2*p4
 p1*p3
 x3(t)
 x1(t)*x2(t)
 x1(t)*p4 + x2(t)*p2
 x1(t)*p1*p4 + x2(t)*p2*p3

By default, find_identifiable_functions tries to simplify the output functions as much as possible, and it has simplify keyword responsible for the degree of simplification. The default value is :standard but one could use :strong to try to simplify further (at the expense of heavier computation) or use :weak to simplify less (but compute faster).

When several generating sets describe the same field, the cmp keyword can be used to control which rational functions are considered simpler and end up in the final result. It is possible to provide a custom function cmp(f, g) that returns true when f should be preferred to g. For example, the following comparator prefers functions with shorter printed representations:

find_identifiable_functions(
    LLW1987,
    with_states = true,
    cmp = (f1, f2) -> length(string(f1)) < length(string(f2)),
)

7-element Vector{AbstractAlgebra.Generic.FracFieldElem{Nemo.QQMPolyRingElem}}:
 p1*p3
 p2*p4
 x3(t)
 p1 + p3
 x1(t)*x2(t)
 x1(t)*p4 + x2(t)*p2
 x1(t)*p1*p4 + x2(t)*p2*p3

As assess_identifiability and assess_local_identifiability, find_identifiable_functions accepts an optional parameter loglevel (default: Logging.Info) to adjust the verbosity of logging.

Finally, as for assess_identifiability, an experimental feature allows to provide an additional keyword argument known_ic to indicate functions of states and parameters for which the initial conditions are assumed to be known (while the initial conditions of the system are still assumed to be generic). In this case, the function will find identifiable functions of parameters and initial conditions rather than of parameters and states. Let us add an assumption that the initial condition x1(0) is known:

find_identifiable_functions(LLW1987, known_ic = [x1])

8-element Vector{AbstractAlgebra.Generic.FracFieldElem{Nemo.QQMPolyRingElem}}:
 p1 + p3
 p2*p4
 p1*p3
 x3(0)
 x2(0)
 x1(0)
 p2*x2(0) + p4*x1(0)
 p1*p4*x1(0) + p2*p3*x2(0)

We see that x2(0) becomes an identifiable function as well, which is natural since x1(t) * x2(t) was an identifiable function before.