<figure><img src="image_1.png" /></figure>

HAL
open science

# Identifiability of nonlinear ODE Models with Time-Varying Parameters: the General Analytical Solution and Applications in Viral Dynamics

Agostino Martinelli

► To cite this version:

Agostino Martinelli. Identifiability of nonlinear ODE Models with Time-Varying Parameters: the General Analytical Solution and Applications in Viral Dynamics. 2024. hal-04420973

HAL Id: hal-04420973

https://inria.hal.science/hal-04420973v1

Preprint submitted on 27 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Identifiability of nonlinear ODE Models with Time-Varying
Parameters: the General Analytical Solution and
Applications in Viral Dynamics

Agostino Martinelli<fnref n="1" />

December 4, 2023

<footnote n="1">¹A. Martinelli is with INRIA Rhone Alpes, Montbonnot, France e-mail: agostino.martinelli@inria.fr</footnote>
# Contents

<table><tbody><tr><td>1 Introduction</td><td>3</td></tr><tr><td>2 Basic definitions and operations</td><td>7</td></tr><tr><td>2.1 System characterization</td><td>7</td></tr><tr><td>2.2 Problem statement</td><td>8</td></tr><tr><td>2.2.1 Observability analysis</td><td>8</td></tr><tr><td>2.2.2 Identifiability analysis</td><td>8</td></tr><tr><td>2.3 Basic algebraic operations and notions</td><td>8</td></tr><tr><td>2.4 Minimal Invariant codistributions and distributions</td><td>10</td></tr><tr><td>2.5 Unknown input degree of reconstructability</td><td>11</td></tr><tr><td>2.6 Definition of the case study</td><td>11</td></tr><tr><td>3 Systematic procedure to perform the observability analysis</td><td>13</td></tr><tr><td>3.1 Basic operations in Algorithm 1</td><td>13</td></tr><tr><td>3.1.1 deg<sub>w</sub>(Ω) operation</td><td>13</td></tr><tr><td>3.1.2 S(Σ, Ω) operation</td><td>13</td></tr><tr><td>3.1.3 R(Σ, λ<sub>1</sub>,..., λ<sub>m</sub>) operation</td><td>15</td></tr><tr><td>3.1.4 A(Σ, m) operation</td><td>15</td></tr><tr><td>3.2 Algorithm 1 for systems canonic with respect to the unknown inputs</td><td>16</td></tr><tr><td>3.2.1 μ and ν</td><td>16</td></tr><tr><td>3.2.2 ĝ<sup>α</sup></td><td>16</td></tr><tr><td>3.2.3 The codistribution Õ</td><td>17</td></tr><tr><td>3.2.4 Final step</td><td>17</td></tr><tr><td>3.3 Algorithm 1 for systems that are not in canonical form</td><td>17</td></tr><tr><td>3.3.1 m<sup>α</sup>μ and m<sup>α</sup>ν</td><td>17</td></tr><tr><td>3.3.2 m<sup>α</sup>ĝ</td><td>18</td></tr><tr><td>3.3.3 The nested loop</td><td>19</td></tr><tr><td>3.4 Illustration by the case study</td><td>19</td></tr><tr><td>4 Systematic procedure to perform the identifiability analysis</td><td>22</td></tr><tr><td>4.1 Identifiability when the state is observable</td><td>22</td></tr><tr><td>4.2 Unobservability of the State, State Symmetries and Indistinguishable States</td><td>24</td></tr><tr><td>4.3 Identifiability for any system</td><td>24</td></tr><tr><td>4.3.1 Unknown input symmetries due to the system non canonicity</td><td>25</td></tr><tr><td>4.3.2 Possible unknown input symmetries due to the state unobservability</td><td>25</td></tr><tr><td>4.3.3 The systematic procedure</td><td>26</td></tr><tr><td>4.4 Illustration by the case study</td><td>27</td></tr><tr><td>4.4.1 Symmetries and Indistinguishable states</td><td>27</td></tr><tr><td>4.4.2 Characterization by an observable state</td><td>29</td></tr><tr><td>4.4.3 Identifiability</td><td>29</td></tr><tr><td>5 Indistinguishable states and parameters and minimal missing information to achieve observability/identifiability</td><td>31</td></tr><tr><td>5.1 Indistinguishable states and parameters</td><td>31</td></tr><tr><td>5.2 Minimal missing information to achieve observability/identifiability</td><td>33</td></tr><tr><td>5.3 Illustration by the case study</td><td>34</td></tr></tbody></table>
<table>
  <tbody>
    <tr>
      <td><b>6 HIV infection</b></td>
      <td><b>36</b></td>
    </tr>
    <tr>
      <td>6.1 Observability analysis</td>
      <td>37</td>
    </tr>
    <tr>
      <td>6.2 Symmetries and indistinguishable states</td>
      <td>38</td>
    </tr>
    <tr>
      <td>6.3 Characterization by an observable state</td>
      <td>39</td>
    </tr>
    <tr>
      <td>6.4 Identifiability analysis</td>
      <td>39</td>
    </tr>
    <tr>
      <td>6.5 Indistinguishable states and unknown inputs</td>
      <td>40</td>
    </tr>
    <tr>
      <td>6.6 Comparison with the state of the art results</td>
      <td>41</td>
    </tr>
    <tr>
      <td>6.7 Minimal external information</td>
      <td>44</td>
    </tr>
    <tr>
      <td><b>7 Compartmental SEIAR model of the epidemic dynamics of Covid-19</b></td>
      <td><b>46</b></td>
    </tr>
    <tr>
      <td>7.1 Observability analysis</td>
      <td>47</td>
    </tr>
    <tr>
      <td>7.2 Symmetries and indistinguishable states</td>
      <td>48</td>
    </tr>
    <tr>
      <td>7.2.1 First symmetry</td>
      <td>48</td>
    </tr>
    <tr>
      <td>7.2.2 Second symmetry</td>
      <td>49</td>
    </tr>
    <tr>
      <td>7.3 Characterization by an observable state</td>
      <td>50</td>
    </tr>
    <tr>
      <td>7.4 Identifiability analysis</td>
      <td>50</td>
    </tr>
    <tr>
      <td>7.4.1 First symmetry</td>
      <td>51</td>
    </tr>
    <tr>
      <td>7.4.2 Second symmetry</td>
      <td>51</td>
    </tr>
    <tr>
      <td>7.5 Indistinguishable states and unknown inputs</td>
      <td>52</td>
    </tr>
    <tr>
      <td>7.5.1 First set of solutions</td>
      <td>52</td>
    </tr>
    <tr>
      <td>7.5.2 Second set of solutions</td>
      <td>53</td>
    </tr>
    <tr>
      <td>7.6 Numerical results</td>
      <td>54</td>
    </tr>
    <tr>
      <td>7.6.1 β = 1</td>
      <td>55</td>
    </tr>
    <tr>
      <td>7.6.2 β(t) = cos(π/400t)</td>
      <td>58</td>
    </tr>
    <tr>
      <td><b>8 Genetic toggle switch</b></td>
      <td><b>61</b></td>
    </tr>
    <tr>
      <td>8.1 Implementation of Algorithm 1</td>
      <td>61</td>
    </tr>
    <tr>
      <td>8.2 Identifiability analysis and indistinguishable states and parameters</td>
      <td>62</td>
    </tr>
    <tr>
      <td>8.3 Direct proof</td>
      <td>64</td>
    </tr>
    <tr>
      <td>8.3.1 First set</td>
      <td>64</td>
    </tr>
    <tr>
      <td>8.3.2 Second set</td>
      <td>64</td>
    </tr>
    <tr>
      <td>8.3.3 Third set</td>
      <td>64</td>
    </tr>
    <tr>
      <td>8.3.4 Fourth set</td>
      <td>65</td>
    </tr>
    <tr>
      <td>8.3.5 Fifth set</td>
      <td>65</td>
    </tr>
    <tr>
      <td>8.3.6 Sixth set</td>
      <td>65</td>
    </tr>
    <tr>
      <td><b>9 Further simple examples</b></td>
      <td><b>66</b></td>
    </tr>
    <tr>
      <td>9.1 First Example</td>
      <td>66</td>
    </tr>
    <tr>
      <td>9.1.1 Observability analysis</td>
      <td>66</td>
    </tr>
    <tr>
      <td>9.1.2 Identifiability analysis</td>
      <td>67</td>
    </tr>
    <tr>
      <td>9.2 Second example</td>
      <td>67</td>
    </tr>
    <tr>
      <td>9.2.1 Observability analysis</td>
      <td>67</td>
    </tr>
    <tr>
      <td>9.2.2 Identifiability analysis</td>
      <td>68</td>
    </tr>
    <tr>
      <td>9.3 Third example</td>
      <td>69</td>
    </tr>
    <tr>
      <td>9.3.1 Observability analysis</td>
      <td>69</td>
    </tr>
    <tr>
      <td>9.3.2 Identifiability analysis</td>
      <td>70</td>
    </tr>
    <tr>
      <td><b>10 Conclusion</b></td>
      <td><b>71</b></td>
    </tr>
    <tr>
      <td>10.1 First contribution</td>
      <td>71</td>
    </tr>
    <tr>
      <td>10.2 Second contribution</td>
      <td>71</td>
    </tr>
    <tr>
      <td>10.3 Third contribution</td>
      <td>72</td>
    </tr>
    <tr>
      <td>10.4 Fourth contribution</td>
      <td>72</td>
    </tr>
    <tr>
      <td><b>A Proof of Theorem 4.2</b></td>
      <td><b>75</b></td>
    </tr>
    <tr>
      <td><b>B Proof of Theorem 4.3</b></td>
      <td><b>77</b></td>
    </tr>
  </tbody>
</table>
Abstract

Identifiability is a structural property of any ODE model characterized by a set of unknown parameters. It describes the possibility of determining the values of these parameters from fusing the observations of the system inputs and outputs. This paper finds the general analytical solution of this fundamental problem and, based on this, provides a general and automated analytical method to determine the identifiability of the unknown parameters. In particular, the method can handle any model, regardless of its complexity and type of non-linearity, and provides the identifiability of the parameters even when they are time-varying. In addition, it is automatic as it simply needs to follow the steps of a systematic procedure that only requires to perform the calculation of derivatives and matrix ranks. Time-varying parameters are treated as unknown inputs and their identification is based on the very recent analytical solution of the unknown input observability problem [1, 2]. The method is used to determine the identifiability of the unknown time-varying parameters that characterize two non-linear models in the field of viral dynamics (HIV and Covid-19) and a non-linear model that characterizes the genetic toggle switch. New fundamental properties that characterize these models are determined and discussed in detail through a comparison with the state-of-the-art results. In particular, regarding the very popular HIV ODE model and the genetic toggle switch model, the method automatically finds new important results that are in contrast with the results in the current literature.

**Keywords:** Structural Identifiability; System Identification; Nonlinear observability; Unknown Input Observability; Unknown Input Reconstruction; Viral dynamics; HIV dynamics; Covid-19 dynamics; Genetic toggle switch.
# Highlights

1. Algorithm that solves the most general Unknown Input Observability problem usable by a non specialist user.

2. Algorithm that determines the time-varying parameter identifiability of any ODE model

3. Determination of the continuous transformations (one parameter Lie groups) which allow us to construct the set of indistinguishable states and indistinguishable unknown inputs in the presence of unobservability and unidentifiability.

4. Observability and identifiability of a very popular HIV model, a Covid-19 model, and a model that characterizes the genetic toggle switch.

5. Detection of errors in the state of the art about the HIV ODE model identifiability and the genetic toggle switch identifiability. In particular, detection of a serious error in Section 6.2 of [3]<fnref n="1" />.

<footnote n="1">¹The interested reader can find a detailed explanation of the serious error made by the authors of [3] in [5]. Note that Section 6.2 of [3] contains two distinct errors. The former, which is a simple typo, was recently acknowledged by the authors (see [4]). However, in [4], the authors still concluded that the system is observable/identifiable. The very serious error committed by the authors is highlighted in Section 3.2 of [5] and lies in an incorrect use of the Inverse Function Theorem. In particular, when using this theorem, the authors did not take into account the constraints due to the dynamics of the system, as shown in Section 3.2 of [5].</footnote>
# Chapter 1

# Introduction

Ordinary Differential Equations (ODEs) are used to model a plethora of phenomena in many scientific domains, ranging from the natural and applied sciences up to the social sciences.

An ODE model is characterized by a state that consists of several scalar quantities. Its time evolution is precisely described by a set of ordinary differential equations (one per each state component). In addition to the state components, an ODE model can also include a set of parameters which may be known or unknown. Finally, an ODE model is characterized by a set of outputs and, very often, by a set of inputs. The inputs are functions of time that act on the state dynamics and that can be assigned (with some restrictions that depend on the specific case). The outputs are available quantities and can be expressed in terms of the state and, in some cases, also in terms of the above inputs and the above model parameters.

*State observability* and *parameter identifiability* are two structural properties of an ODE model. The former characterizes the possibility of inferring the values that the state components take, starting from the knowledge of the system inputs and outputs (e.g., [6, 7, 8, 9, 10]). The latter characterizes the possibility of inferring the values of the model parameters, again, from the knowledge of the system inputs and outputs (e.g., [11, 12, 13]).

Performing an observability and an identifiability analysis is fundamental in order to set up an ODE model, i.e., in order to correctly define the state and the parameters that the information contained in the sensor measurements allows us to estimate. Characterizing a given system with a state that is unobservable and/or with parameters that are unidentifiable can have dramatic consequences during the estimation and/or the identification process. Note that, for many complex sensor fusion problems studied in the past, this task (observability and identifiability analysis) was accomplished and played a crucial role to set up a suitable ODE model. For example, this has been the case for the following sensor fusion problems: visual and inertial (e.g., see the works [14, 15, 16, 17, 18]), visual and pressure (e.g., [19]), GPS and odometry (e.g., [20]), wheel encoders and ultra-wideband (e.g., [21]), radio beacons and inertial (e.g., [22]). Additionally, this has also been the case in the framework of distributed estimation over a low-cost sensor network (e.g., see [23] and references therein). In most previous works (and in many others), the observability analysis of the considered sensor fusion problem was essential to set up a proper ODE model (note that in the above works, the concept of identifiability was not distinguished from observability and, with observability analysis, it was actually intended observability and identifiability analysis). More in general, the key role of observability in information fusion has recently been investigated in [24].

An ODE model may also include the presence of *unknown* inputs, i.e., inputs that act on the system dynamics precisely as the aforementioned inputs. However, they differ from them because they are unknown. From a mathematical point of view, an unknown input plays exactly the same role of an unknown time-varying parameter. In the literature, the possibility of inferring their values has been called in different manners: *input observability*, *input reconstruction*, and sometimes *system invertibility* (e.g., [25, 26, 27, 28, 29, 30, 31, 32]).

Throughout this paper, we adopt the following terminology. *Observability* refers to the components of the state. *Unknown Input Observability* (abbreviated UIO) still refers only to the state components, but when the ODE model also includes the presence of unknown inputs (or time-varying parameters). *Identifiability* refers to all the unknown parameters. When a parameter is
time-varying we also use the term *unknown input reconstruction* to mean its identifiability.

We emphasize that the identifiability of constant parameters can be addressed using methods to study observability (e.g., through the use of the observability rank condition, [8, 9], when the system does not contain unknown time-varying parameters). Indeed, the time evolution of a constant parameter is trivially known (its time derivative vanishes, by definition). As a result, by including it in the state, we get a new ODE model in which the presence of this parameter is eliminated as it has become a component of the state (see the last paragraph of Section 2.2). In the presence of time-varying unknown parameters, the same approach does not provide the same immediate result. If a time-varying parameter is included in the state, the resulting ODE model will be characterized by a new unknown parameter, which is its time derivative (in general time-varying).

ODE identifiability analysis is present in a large variety of scientific domains and several interesting methods have been proposed. An exhaustive review of the methods proposed up to 2011 can be found in [3]. Some of these methods are very interesting and, in many cases, they can successfully be used to detect all the identifiability properties of a given ODE model. However, they have the following fundamental limitations:

* They are not general. They are based on several properties that sometimes cannot be exploited in the presence of a given type of system nonlinearity (i.e., a given function *f* and/or *h* in Equation (2.1)).

* They cannot be executed automatically. In most of cases, they must be adapted to the specific case under investigation and this adjustment requires inventiveness by the user. In particular, this process cannot be carried out by simply following the steps of a systematic procedure (e.g., by running a code without human intervention).

The second limitation becomes particularly relevant in the case of time varying parameters where, in some cases, an erroneous application of these methods provided wrong results. Specifically, in our study about the identifiability of the ODE HIV model given in Section 6.4, we obtain a fundamental new result that contradicts the result obtained in Section 6.2 of [3], which was obtained by using one of these methods (a detailed explanation of the serious error made by the authors of [3] is available in [5]).

In the last decade, thanks to the progress in UIO, several new interesting approaches have been proposed. Specifically, the nonlinear UIO problem was approached by introducing an extended state that includes the original state together with the unknown inputs and their time derivatives up to a given order [33, 34, 35, 36, 37]. In [35, 36, 37], this was used to study the identifiability of the unknown time-varying parameters. In particular, in [35, 36, 37], several automatic iterative algorithms were introduced. These algorithms work automatically and are able to check the identifiability of time-varying parameters. On the other hand, they suffer from the following limitations:

* They do not converge, automatically. In particular, at each iterative step, the state is extended by including new time derivatives of the unknown inputs. Consequently, if the extended state is observable at a given step, convergence is achieved. However, if this were not the case, we can never exclude that, at a later step, the extended state becomes observable. Therefore, in the presence of unobservability, all these algorithms remain in-conclusive.

* Due to the previous state augmentation, the computational burden can easily become prohibitive after a few steps.

Very recently, we introduced the analytical solution of the nonlinear UIO [2], starting from the results/derivations presented in [1]<fnref n="1" />. In particular, the algorithm introduced in [2] does not encounter the aforementioned limitations<fnref n="2" /> and provides a full answer to the problem of the state

<footnote n="1">¹This analytical solution is also based on the results obtained in [38] and [39], which deal with special systems, namely characterized by a single unknown input, and dynamics linear with respect to this unknown input.</footnote>

<footnote n="2">²Note that, the solution introduced in [2] could need the inclusion of some of the UIs in the state. However, this inclusion is targeted and terminates in a finite number of steps, as soon as the extended system achieves its *highest unknown input degree of reconstructability* (see Section 5 in [2]).</footnote>
observability for any nonlinear ODE model, in the presence of unknown inputs (or time-varying parameters). On the other hand, the algorithm introduced in [2] deals with the observability of the state and not with the identifiability of the time-varying parameters. In [2], we exploited this analytical solution to provide a preliminary result also about the identifiability of the time-varying parameters (unknown input reconstruction). One of the goals of this paper is the extension of this preliminary and partial result, and the introduction of a complete systematic procedure to study the identifiability of the time-varying parameters. Then, we adopt this systematic procedure to study the identifiability of a very popular HIV model (e.g., [3, 35, 40, 41, 42]) and a Covid-19 model [43, 44].

The contributions of this paper are the following four:

1. Compact presentation of the analytical solution introduced in [2] (Chapter 3) together with its basic mathematical foundation introduced in [1] (Chapter 2) in order to make this solution usable by a non specialist user.

2. Introduction of the systematic procedure to determine the identifiability of all the time-varying parameters of any ODE model (Chapter 4).

3. Determination of the continuous transformations (one parameter Lie groups) that allow us to build the set of indistinguishable states and indistinguishable unknown inputs in the presence of unobservability and unidentifiability (Chapter 5 and, in particular, the two systems of differential equations in (5.5) and (5.7)).

4. Determination of the observability and the identifiability properties of a very popular HIV model, a Covid-19 model, and a genetic toggle switch model (Chapters 6, 7, and 8, respectively), obtained by using the three above contributions.

Note that the first two contributions consist of two systematic procedures that fully allow us to perform the observability analysis and the identifiability analysis, respectively. Specifically, they are Algorithm 1 in Chapter 3 and Algorithm 4 in Chapter 4. In other words, given any nonlinear system, we can obtain the observability of the state and the identifiability of the parameters by simply following the steps of these systematic procedures.

In contrast, the third contribution does not consist of a systematic procedure. In particular, all that can be done automatically is the determination of the set of differential equations in (5.5) and (5.7) and the determination of the analytical solution of (5.5), which is trivial. However, the method to solve (5.7) depends on the specific case and it may also be not possible to obtain an analytical solution. Note that this third contribution is unnecessary for the observability and the identifiability analyses (which are fully achieved by running Algorithm 1 and Algorithm 4, respectively). On the other hand, the determination of the analytical solution of (5.7), when possible, provides very useful insights and additional very interesting properties of the ODE model under investigation.

This third contribution adopts precisely the same definition of Lie symmetry introduced in [45, 46]. However, in [45, 46] this definition was then adopted by using the aforementioned extended state, which includes the original state and the unknown inputs with their time derivatives up to a given order (i.e., the augmented state introduced in [33, 34, 35, 36, 37] to investigate the state observability in the presence of unknown inputs).

Regarding the fourth contribution, our results improve previous results in the state of the art. In particular, in all the three cases, Algorithm 4 automatically provides the unidentifiability of the time-varying parameters (for the HIV model and the genetic toggle switch model this result is in contrast with the state of the art). In other words, in all the three cases, the time-varying parameters cannot be uniquely identified. In addition, for these ODE models, it is possible to obtain an analytical solution of the differential equation in (5.7). This allows us to determine indistinguishable unknown inputs that agree with the same outputs of the model (these models do not contain known inputs and, consequently, the only source of information about the state and the model parameters comes from the outputs). Regarding the HIV ODE model, in Section 6.6 we use a data set available in the literature and, by using the solution of (5.7), we generate infinite different unknown inputs that produce the same outputs, starting from the true unknown input and the true states. This unequivocally shows the unidentifiability of the system, in contrast with the results available in the state of the art (e.g., see Section 6.2 of [3] and a detailed explanation
of the serious error made by its authors, available in [5]). Finally, we determine the minimal
external information (external to the knowledge of the outputs) requested to uniquely determine
the system parameters. In Section 7.6, we perform a similar analysis for the Covid-19 ODE
model here investigated.

The paper is organized as follows. Chapter 2 provides a basic mathematical characterization
of the problem together with basic mathematical concepts and definitions. It also introduces an
elementary and illustrative example from robotics (from now on, the *case study*). We refer to
the case study throughout all the paper to better illustrate all the theoretical concepts and the
analytical procedures here introduced. In particular, each section ends with a subsection where
the theoretical concepts provided by the section are applied to the case study. Chapter 3 provides
a summary of our recent solution of the nonlinear UIO introduced in [2]. Chapter 4 introduces the
systematic procedure to obtain the identifiability of any ODE model (i.e., Algorithm 4). Chapter
5 introduces the set of differential equations to generate indistinguishable unknown inputs and
initial states, simultaneously (i.e., Equations (5.5) and (5.7)). Chapters 6, 7, and 8 investigate
the HIV model discussed in [3, 35, 40, 41, 42], the Covid-19 model introduced in [43, 44], and
the genetic toggle switch model discussed in [35]. Chapter 9 provides the study of further very
simple examples. Finally, Chapter 10 provides our conclusion.

A short version of this paper is also available in [47].
# Chapter 2

# Basic definitions and operations

## 2.1 System characterization

We start from the following very general ODE model:

$$ \begin{cases} \dot{X} = f(X(t), t, U(t), Q, W(t)) \\ y = h(X(t), t, U(t), Q, W(t)), \end{cases} \qquad (2.1) $$

where: $X(t) \in \mathbb{R}^m$ is the state<fnref n="1" />, $y(t) \in \mathbb{R}^p$ the output vector, $U(t) \in \mathbb{R}^{m_u}$ the known system input vector, $Q \in \mathbb{R}^q$ the set of the unknown constant parameters, and $W(t) \in \mathbb{R}^{m_w}$ the set of the unknown parameters that depend on time. Note that $W(t)$ can be regarded as a system input vector that, precisely as $U(t)$, acts on the system dynamics. However, it differs from $U(t)$ in two fundamental respects: (i) it cannot be assigned, and (ii) it is unknown.

The problem that we solve in this paper is the introduction of the analytic condition that fully characterizes the identifiability of the system parameters, both constant and time-varying. This is obtained by exploiting recent results on the unknown input observability problem. To exploit these results, instead of (2.1) we adopt the following system characterization:

$$ \begin{cases} \dot{x} = g^0(x, t) + \sum_{k=1}^{m_u} f^k(x, t)u_k(t) + \sum_{j=1}^{m_w} g^j(x, t)w_j(t) \\ y = [h_1(x, t), \dots, h_p(x, t)], \end{cases} \qquad (2.2) $$

where:

* $x \in \mathbb{R}^n$ is the state.

* $y \in \mathbb{R}^p$ is the output vector.

* $u_1(t), \dots, u_{m_u}(t)$ are the known inputs.

* $w_1(t), \dots, w_{m_w}(t)$ are the unknown inputs or the unknown time-varying parameters.

* $f^1, \dots, f^{m_u}, g^0, g^1, \dots, g^{m_w}$ are $m_u+m_w+1$ vector fields, which are assumed to be smooth functions of $x$ and $t$.

From now on, we denote the system in (2.2) by $\Sigma$. Note that the above characterization can easily account for the presence of constant parameters ($Q$) by including all of them in the state $x$ and by suitably setting $f^1, \dots, f^{m_u}, g^0, g^1, \dots, g^{m_w}$, namely, by setting to zero all their components that yield $\dot{Q}$ in (2.2). Note that many nonlinear systems have the structure in (2.2). When this is not directly the case, it is possible to easily convert (2.1) to (2.2). This is obtained by setting $u(t) = \dot{U}, w(t) = \dot{W}$ (with $u = [u_1, \dots, u_{m_u}]$ and $w = [w_1, \dots, w_{m_w}]$) and by including $U(t), W(t)$ and $Q$ in the state (i.e., $x = [X^T, U^T, W^T, Q^T]^T$). In the rest of this paper, we directly refer to the characterization given in (2.2).

<footnote n="1">¹In general, instead of $\mathbb{R}^m$, the state belongs to a differential manifold of dimension $m$. This also holds for the other quantities that appear in (2.1) (and also in (2.2)), with the appropriate dimension. In this paper, for the sake of simplicity, we avoid the concept of differential manifold (the reader familiar with this concept can simply replace the Euclidean space of a given dimension with the differential manifold of the same dimension).</footnote>
## 2.2 Problem statement

Given the system characterized by (2.2) (which is equivalent to (2.1)), the problem we solve in this paper is the introduction of the systematic procedures to automatically perform the observability analysis and the identifiability analysis.

### 2.2.1 Observability analysis

As we mentioned in the introduction, the paper summarizes (Chapter 3) the solution introduced in our previous paper [2]. The solution provides the analytical condition to easily check whether the state that characterizes the system in (2.2) is observable or not. It actually tells us more than this. It builds the entire observation space, i.e., the set of all the observable functions. An observable function is a scalar function of the state $x$ whose value, at a given time $t_0$, can be expressed in terms of the values that the known inputs and the system outputs take on a given time interval that includes $t_0^2$. In practice, even when the state is unobservable, the solution here summarized provides all the functions of the state that are observable. For instance, let us suppose to have a vehicle that moves on a plane and let us suppose that, with respect to a reference frame, we cannot obtain from the inputs and the outputs its Cartesian coordinates $(x_1, x_2)$ but we can obtain its distance from the origin. If the state $x$ that characterizes this system contains $x_1$ and $x_2$, the state is unobservable. However, the function of the state $\theta(x) := \sqrt{x_1^2 + x_2^2}$ is observable. Note that this is precisely the case of our illustrative example introduced in Section 2.6 (the case study), where, however, we adopt polar coordinates $(\rho$ and $\phi$ instead of $x_1$ and $x_2$).

### 2.2.2 Identifiability analysis

Chapter 4 provides the systematic procedure to check if a given time-varying parameter (or unknown input) can be identified. If this is not the case, we analytically compute all the set of continuous transformations that transform the value taken by a given unidentifiable parameter in other indistinguishable values for that parameter. This is carried out in Chapter 5 and the set of continuous transformations are the differential equations in (5.5) and (5.7).

Note that, the identifiability of constant parameters is actually obtained by performing the above observability analysis. Indeed, the constant parameters are included in the state. To check if a given constant parameter $q$ is identifiable, we need to trivially check if the function $\theta(x) := q$ belongs to the observation space (i.e., we need to check if its gradient $(\nabla\theta)$ belongs to the observability codistribution $(\mathcal{O})$, which is automatically computed by Algorithm 1).

## 2.3 Basic algebraic operations and notions

In accordance with the control theory literature, we use the term *distribution* to denote the span of a set of $d \le n$ vector fields$^3$, $\Delta =\operatorname{span}\{\tau^1, \tau^2, \dots, \tau^d\}$. It is basically a vector space that depends on $x \in \mathbb{R}^n$. Similarly, we also use the term *codistribution* to denote the span of a set of $s \le n$ covector fields$^4$, $\Omega =\operatorname{span}\{\omega_1, \omega_2, \dots, \omega_s\}$. Again, it is a vector space that depends on $x \in \mathbb{R}^n$. In this paper, we always consider non singular distributions and codistributions. In other words, we always refer to an open set of $\mathbb{R}^n$ where the dimension of the considered distribution (or codistribution) takes the same value.

<footnote>$^2$The definition of observation space is provided in [38], where Definition 2 defines this concept starting from the concept of indistinguishability. An observable function is precisely an element of this function space.</footnote>

<footnote>$^3$Simply speaking, a vector field is a column vector of dimension $n$ whose components are functions of $x$. Actually, the correct definition of vector field must take into account its tensor nature, i.e., how its components change under a generic change of coordinates. In particular, a vector is a tensor of type (0, 1) (see Section 2.2 of [1]).</footnote>

<footnote>$^4$A covector field is the dual of a vector field. Simply speaking, it is a row vector of dimension $n$ whose components are functions of $x$. Actually, the correct definition of covector field must take into account its tensor nature, i.e., how its components change under a generic change of coordinates. In particular, a covector is a tensor of type (1, 0) (see Section 2.2 of [1]).</footnote>
We remind the reader of the definitions of the Lie derivative along a vector field $f$ of a scalar field $\lambda$, of a vector field $\tau$, and of a covector field $\omega$. We have, respectively:

$$ \mathcal{L}_f \lambda = \frac{\partial \lambda}{\partial x} \cdot f \qquad (2.3) $$

$$ \mathcal{L}_f \tau = [f, \tau] = \frac{\partial \tau}{\partial x} \cdot f - \frac{\partial f}{\partial x} \cdot \tau \qquad (2.4) $$

$$ \mathcal{L}_f \omega = f^T \cdot \left( \frac{\partial \omega^T}{\partial x} \right)^T + \omega \cdot \left( \frac{\partial f}{\partial x} \right) \qquad (2.5) $$

The square brackets in (2.4) are called *Lie brackets*. In the special case when the covector field $\omega$ is the gradient of a scalar field (i.e., $\omega = \nabla \lambda$), Equations (2.5) and (2.3) yield:

$$ \mathcal{L}_f \nabla \lambda = \nabla \mathcal{L}_f \lambda. $$

We use the following definitions/notation:

* Given a distribution $\Delta$ and a vector field $f$, we set $\mathcal{L}_f \Delta$ the distribution that is the span of all the vectors $\mathcal{L}_f \tau = [f, \tau]$, for any $\tau \in \Delta$. Similarly, given a codistribution $\Omega$ and a vector field $f$, we set $\mathcal{L}_f \Omega$ the codistribution that is the span of all the covectors $\mathcal{L}_f \omega$ for any $\omega \in \Omega$.

* We use the term *integrable codistribution* to denote the special codistribution that can be generated by the gradients of a set of scalar fields (i.e., $\Omega = \text{span}\{\nabla \lambda_1, \nabla \lambda_2, \dots, \nabla \lambda_s\}$). In this case, we often use the term *generators* to denote the scalar functions $\lambda_1, \dots, \lambda_s$ (instead of their gradients).

* Given two vector spaces $V_1$ and $V_2$, $V_1+V_2$ is their sum, i.e., the span of all the generators of $V_1$ and $V_2$. When we have $k(> 2)$ vector spaces $V_1, \dots, V_k$, we denote their sum in the compact notation $\sum_{i=1}^k V_i$.

* Greek indices take non negative values (e.g., $\alpha = 0, 1, 2, \dots$). Latin indices take positive values (e.g., $i = 1, 2, \dots$).

In addition, we remind the reader of a further algebraic operation introduced in [2] and called the *autobracket*. Given a set of $L + 1$ vector fields, $\tau^0, \tau^1, \dots, \tau^L$, and a non singular two-index tensor $\sigma^\alpha_\beta, \alpha, \beta = 0, 1, \dots, L$, the autobracket of a vector field $f$ along $\gamma$ ($\gamma = 0, 1, \dots, L$), denoted by $\sigma_\tau[f]^\gamma$, is defined as follows:

$$ \sigma_{\tau}[f]^{\gamma} = \sum_{\beta=0}^{L} \sigma_{\beta}^{\gamma} [\tau^{\beta}, f] + \delta_{0}^{\gamma} \frac{\partial f}{\partial t}, \quad \gamma = 0, 1, \dots, L, \qquad (2.6) $$

where the square brackets on the right hand side are the Lie brackets and $\delta_0^\gamma$ (the Kronecker delta) is 1 when $\gamma = 0$ and is 0 otherwise. In [2], and in this paper, we use this operation for two settings (note that in [2] the operation was directly defined for these settings):

1. $L = m_w$, $\tau^0 = g^0$, $\tau^1 = g^1$, $\dots$, $\tau^{m_w} = g^{m_w}$, which are the $m_w + 1$ vector fields that appear in (2.2)), and $\sigma = \nu$, which is the inverse of $\mu$ defined in Section 3.2.1 (see Equation (3.3)). For simplicity, we denote this operation by $[\cdot]^\gamma$ (instead of $g^\gamma[\cdot]^\gamma$)

2. $L = m$, $\tau^0 = g^0$, $\tau^1 = g^1$, $\dots$, $\tau^m = g^m$, which are the first $m + 1$ vector fields that appear in (2.2) (once the unknown inputs are re-ordered, in accordance with the $\mathcal{R}$ operation, defined in Section 3.1.3), and $\sigma = m\nu$, which is the inverse of $m\mu$ defined in Section 3.3.1 (see Equation (3.6)). For simplicity, we denote this operation by $[\cdot]^\gamma$ (instead of $g^\gamma[\cdot]^\gamma$).

We have $L+1$ possible first order autobrackets ($\gamma = 0, 1, \dots, L$). By applying the autobracket $k$ consecutive times, we obtain $(L+1)^k$ vector fields that we denote by $\sigma_\tau[f]^{(\gamma_1, \dots, \gamma_k)}$ and we have $\sigma_\tau[f]^{(\gamma_1, \dots, \gamma_{k-1}, \gamma_k)} = \sigma_\tau[\sigma_\tau[f]^{(\gamma_1, \dots, \gamma_{k-1})}]^{\gamma_k}$.

Finally, given a distribution $\Delta$, we set $\sigma_\tau[\Delta]^\gamma$ the span of all the vectors $\sigma_\tau[f]^\gamma$, for any $f \in \Delta$.
## 2.4 Minimal Invariant codistributions and distributions

Given a distribution $\Delta$ and a set of vector fields $\tau^1, \dots, \tau^d$ defined on $\mathbb{R}^n$, we denote by

$$ \langle \tau^1, \dots, \tau^d | \Delta \rangle $$

the smallest distribution that contains $\Delta$ and such that, for any $f \in \Delta$, we have $\mathcal{L}_{\tau^i} f = [\tau^i, f] \in \langle \tau^1, \dots, \tau^d | \Delta \rangle$, for any $i = 1, \dots, d$. This minimal distribution can be easily computed by a simple recursive algorithm, which is [10]:

$$ \left\{ \begin{array}{l} \Delta_0 = \Delta \\ \Delta_{k+1} = \Delta_k + \sum_{i=1}^{d} \mathcal{L}_{\tau^i} \Delta_k \end{array} \right. , \qquad (2.7) $$

where we provide the initialization and the recursive step. This algorithm converges at the smallest integer $j$ for which $\Delta_j = \Delta_{j-1}$ and $j \le n - \dim\{\Delta\} + 1$. In addition, $\langle \tau^1, \dots, \tau^d | \Delta \rangle = \Delta_{j-1}^5$.

Similarly, given a codistribution $\Omega$ and a set of vector fields $\tau^1, \dots, \tau^d$ defined on $\mathbb{R}^n$, we denote by

$$ \langle \tau^1, \dots, \tau^d | \Omega \rangle $$

the smallest codistribution that contains $\Omega$ and such that, for any $\omega \in \Omega$, we have $\mathcal{L}_{\tau^i} \omega \in \langle \tau^1, \dots, \tau^d | \Omega \rangle$, for any $i = 1, \dots, d$. This minimal codistribution can be easily computed by a simple recursive algorithm, which is [10]:

$$ \left\{ \begin{array}{l} \Omega_0 = \Omega \\ \Omega_{k+1} = \Omega_k + \sum_{i=1}^{d} \mathcal{L}_{\tau^i} \Omega_k \end{array} \right. , \qquad (2.8) $$

This algorithm converges at the smallest integer $j$ for which $\Omega_j = \Omega_{j-1}$ and $j \le n - \dim\{\Omega\} + 1$. In addition, $\langle \tau^1, \dots, \tau^d | \Omega \rangle = \Omega_{j-1}$.

Note that the execution of the above algorithms only requires to compute the Lie derivatives along $\tau^1, \dots, \tau^d$ of the generators of the previous distribution / codistribution.

In the following, we also need two modified versions of the two above algorithms. In particular, these new algorithms compute the following two quantities:

1. Given a distribution $\Delta$, a set of vector fields $\tau^0, \tau^1, \dots, \tau^L$ defined on $\mathbb{R}^n$, and a non singular two-index tensor $\sigma$, we denote by

$$ \langle \sigma_{\tau} [\cdot] | \Delta \rangle \qquad (2.9) $$

the smallest distribution that contains $\Delta$ and such that, for any $f \in \Delta$, we have $\sigma_{\tau} [f]^{\gamma} \in \langle \sigma_{\tau} [\cdot] | \Delta \rangle$, for any $\gamma = 0, 1, \dots, L$ (where the operation $\sigma_{\tau} [\cdot]$ is defined at the end of Section 2.3, by Equation (2.6)). This minimal distribution can be easily computed by a simple recursive algorithm, which is:

$$ \left\{ \begin{array}{l} \Delta_0 = \Delta \\ \Delta_{k+1} = \Delta_k + \sum_{\gamma=0}^{L} \sigma_{\tau} [\Delta_k]^{\gamma} \end{array} \right. , \qquad (2.10) $$

where we provide the initialization and the recursive step. This algorithm converges at the smallest integer $j$ for which $\Delta_j = \Delta_{j-1}$ and $j \le n - \dim\{\Delta\} + 1$. In addition, $\langle \sigma_{\tau} [\cdot] | \Delta \rangle = \Delta_{j-1}$. Note that also the execution of the above algorithm only requires to compute the autobracket along $\gamma = 0, 1, \dots, L$ of the generators of the previous distribution.

2. Given the codistribution $\Omega$ and the three sets of vector fields: $\tau^1, \dots, \tau^{d_1}, \xi^1, \dots, \xi^{d_2}$, and $\zeta^1, \dots, \zeta^{d_3}$, we denote by

$$ \langle \tau^1, \dots, \tau^{d_1}, \left\{ \begin{array}{c} \xi^1, \dots, \xi^{d_2} \\ \zeta^1, \dots, \zeta^{d_3} \end{array} \right\} | \Omega \rangle \qquad (2.11) $$

the codistribution computed by the algorithm in (2.8), where the sum $\sum_{i=1}^d$ at the recursive step includes all the $d_1$ vector fields $\tau^1, \dots, \tau^{d_1}$. In addition, by denoting with $\omega$ a generator

<footnote n="5">⁵Throughout this paper, when we provide the convergence properties of iterative algorithms that build a family of distributions or codistributions (one distribution or codistribution constructed at each iterative step), we always refer to an open set where these distributions or codistributions are non singular.</footnote>
of $\Omega_k$, the codistribution $\Omega_{k+1}$ also includes, among its generators, all the $d_2$ Lie derivatives of $\omega$ along the vector fields $\xi^1, \dots, \xi^{d_2}$ if and only if $\mathcal{L}_{\zeta^l} \omega$ vanishes for all $l = 1, \dots, d_3$. The algorithm is then interrupted at the smallest integer $j$ such that $\Omega_j = \Omega_{j-1}$.

## 2.5 Unknown input degree of reconstructability

The derivations in [2] and in this paper need the definition of the unknown input degree of reconstructability (Definition 2.2 below), which is based on the definition of the unknown input reconstructability matrix (Definition 2.1 below). These definitions were introduced in [2]. For the sake of completeness, we report them below.

**Definition 2.1** (Unknown input reconstructability matrix). Given the system $\Sigma$ in (2.2) and $k$ scalar functions of the state, $\lambda_1(x), \dots, \lambda_k(x)$, the unknown input reconstructability matrix of $\Sigma$ from $\lambda_1, \dots, \lambda_k$ is defined as follows:

$$ \mathcal{RM}(\lambda_1, \dots, \lambda_k) := \begin{bmatrix} \mathcal{L}_{g^1} \lambda_1 & \mathcal{L}_{g^2} \lambda_1 & \dots & \mathcal{L}_{g^{m_w}} \lambda_1 \\ \mathcal{L}_{g^1} \lambda_2 & \mathcal{L}_{g^2} \lambda_2 & \dots & \mathcal{L}_{g^{m_w}} \lambda_2 \\ \dots & \dots & \dots & \dots \\ \mathcal{L}_{g^1} \lambda_k & \mathcal{L}_{g^2} \lambda_k & \dots & \mathcal{L}_{g^{m_w}} \lambda_k \end{bmatrix} \quad (2.12) $$

**Definition 2.2** (Unknown input degree of reconstructability). Given the system $\Sigma$ in (2.2), and the functions $\lambda_1, \dots, \lambda_k$, the unknown input degree of reconstructability of $\Sigma$ from $\lambda_1, \dots, \lambda_k$ is the rank of $\mathcal{RM}(\lambda_1, \dots, \lambda_k)$.

By construction, given the system characterized by (2.2), its unknown input degree of reconstructability from any set of scalar functions cannot exceed $m_w$. In [2], we introduced the definitions of *Canonic system with respect to its unknown inputs* and *system in Canonical Form with respect to its unknown inputs*. The former is a system such that its unknown input degree of reconstructability from all the observable functions is $m_w$. The latter is a system such that its unknown input degree of reconstructability from the output functions $(h_1, \dots, h_p)$ is $m_w$. As the output functions are observable, a system in canonical form is certainly canonic. However, the viceversa does not hold, in general. Given a canonic system that is not in canonical form, we say that it has been set in canonical form, as soon as we are able to somehow determine a set of observable functions such that the unknown input degree of reconstructability from them is $m_w$.

## 2.6 Definition of the case study

<figure><img src="image_3.png" /><figcaption>Figure 2.1: The vehicle polar coordinates $(\rho, \phi)$, the orientation $(\theta)$, and the angle measured by the camera $(\beta)$.</figcaption></figure>
We introduce an elementary example of an ODE model that will be used in the rest of the paper to illustrate all the concepts and methods that will be introduced.

We consider a vehicle that moves on a 2D-environment. The configuration of the vehicle in a global reference frame can be characterized through the vector [ρ, φ, θ]ᵀ, where ρ and φ are the polar vehicle coordinates, and θ is the vehicle orientation (see Figure 2.1 for an illustration). We assume that the dynamics of this vector satisfy the unicycle model. In polar coordinates we have:

$$ \begin{bmatrix} \dot{\rho} \\ \dot{\phi} \\ \dot{\theta} \end{bmatrix} = \begin{matrix} v \cos(\theta - \phi) \\ \frac{v}{\rho} \sin(\theta - \phi) \\ \omega \end{matrix} \qquad (2.13) $$

where v and ω are the linear and the rotational vehicle speed, respectively, and they are the system inputs. We assume that the vehicle is equipped with a monocular camera that provides the bearing angle of the origin in the vehicle local frame (at the origin there is a landmark perceived by the camera). In other words, the camera provides the angle β in Figure 2.1. We can express β in terms of the vehicle configuration. We have: β = π - (θ - φ). Throughout this paper we will refer to the two scenarios defined by the following settings:

1. v is known and ω is unknown.

2. ω is known and v is unknown.

In both scenarios, the system is directly a special case of (2.2). In particular, x = [ρ, φ, θ]ᵀ, g° = [0, 0, 0]ᵀ, m_u = m_w = 1, p = 1, h₁(x) = φ - θ (we ignore the constant term π in β, which does not affect the observability of the state), and, in the first scenario f¹ = [cos(θ - φ), sin(θ - φ)/ρ, 0]ᵀ, and g¹ = [0, 0, 1]ᵀ, while, in the second scenario, f¹ = [0, 0, 1]ᵀ, and g¹ = [cos(θ - φ), sin(θ - φ)/ρ, 0]ᵀ.
# Chapter 3

# Systematic procedure to perform
the observability analysis

This section provides a very concise summary of the solution of the nonlinear unknown input
observability problem. In particular, it is written to make a non specialist able to implement it.
The interested reader can find all the theoretical foundation and proofs in [2].

The solution consists of a recursive algorithm (Algorithm 1') that, given a system characterized by (2.2), automatically builds the entire observation space in a finite number of iterations. In particular, it builds the so-called *observability codistribution*, which is, by definition, the span of the gradients of all the observable functions. Algorithm 1 builds the generators of the observability codistribution. To this regard, note that, in the algorithm, statements like "$\Omega = \operatorname{span}\{\nabla h_1, \dots, \nabla h_p\}$" (e.g., Line 2 of Algorithm 1) require no action. The algorithm builds the codistribution by building a set of its generators. When it recursively updates the codistribution it simply computes its new generators (and this is obtained by applying certain operations directly on the generators of the previous codistribution or by computing new covectors starting from the quantities that define the original system $\Sigma$ in (2.2)). In the following subsections, we provide all the ingredients necessary to implement Algorithm 1.

## 3.1 Basic operations in Algorithm 1

### 3.1.1 deg<sub>w</sub>(Ω) operation

In Section 2.5, we provided the definition of unknown input degree of reconstructability from a set of scalar functions. Given an integrable codistribution $\Omega = \operatorname{span}\{\nabla\lambda_1, \dots, \nabla\lambda_k\}$, we define the unknown input degree of reconstructability of $\Sigma$ from $\Omega$ the unknown input degree of reconstructability from a set of generators of $\Omega$ (e.g., from $\lambda_1, \dots, \lambda_k$). We denote it by $\operatorname{deg}_w(\Omega)$.

The main loop of Algorithm 1 starts by checking if $\operatorname{deg}_w(\Omega) == m_w$. If this is the case, the main loop is interrupted. Then, Lines 45-47 of Algorithm 1 are executed and Algorithm 1 ends. For the clarity sake, we devote Section 3.2 to describe the behaviour of Algorithm 1 in this case. Before, we provide the definition of further operations adopted by Algorithm 1.

### 3.1.2 $S(\Sigma, \Omega)$ operation

Let us consider a system $\Sigma$ that is characterized by (2.2) and an integrable codistribution $\Omega$. Let us set $m = \operatorname{deg}_w(\Omega)$. The $S(\Sigma, \Omega)$ operation provides a set of $m$ scalar functions, denoted by $\tilde{h}_1, \dots, \tilde{h}_m$, that make the unknown input reconstructability matrix full rank. The operation is denoted by $S$ because is the *Selection* of the aforementioned $\tilde{h}_1, \dots, \tilde{h}_m$ from the generators of $\Omega$. The execution of this operation is immediate because the generators of $\Omega$ are always available when executing Algorithm 1. In this paper, we use the notation $[\tilde{h}_1, \dots, \tilde{h}_m] = S(\Sigma, \Omega)$ to denote the outputs of this operation.

<footnote n="1">¹Algorithm 1 uses the automatic procedures provided by Algorithm 2 and Algorithm 3. Often, when we refer to Algorithm 1, we actually mean all the three algorithms (1, 2, and 3).</footnote>
**Algorithm 1:** The algorithm that builds the observable codistribution ($\mathcal{O}$) for any non-linear system characterized by (2.2). Line 19 runs Algorithm 2 and 3. All the functions/operations that appear in the algorithm are defined in Sections 3.1-3.3. Note that, for TV systems, the Lie derivative operator along $\hat{g}^0$ and $^m\hat{g}^0$, i.e., $\mathcal{L}_{\hat{g}^0}$ and $\mathcal{L}_{^m\hat{g}^0}$, must be replaced by $\mathcal{L}_{\hat{g}^0} := \mathcal{L}_{\hat{g}^0} + \frac{\partial}{\partial t}$ and $\mathcal{L}_{^m\hat{g}^0} := \mathcal{L}_{^m\hat{g}^0} + \frac{\partial}{\partial t}$, respectively.

1. Set $\Sigma$ the system in (2.2).
2. $\Omega = \operatorname{span}\{\nabla h_1, \dots, \nabla h_p\}$.

3. boolean Continue=true.

4. **while** 1 { **This is the MAIN loop**} **do**

5. &nbsp;&nbsp;&nbsp;&nbsp;**if** $\deg_w(\Omega) == m_w$ **then**

6. &nbsp;&nbsp;&nbsp;&nbsp;**break** main loop.

7. &nbsp;&nbsp;&nbsp;&nbsp;**end if**

8. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$m = \deg_w(\Omega) - 1$.

9. **while** $\deg_w(\Omega) > m$ & $\deg_w(\Omega) \neq m_w$ **do**

10. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$m = \deg_w(\Omega)$.

11. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[$\tilde{h}_1, \dots, \tilde{h}_m$] = $\mathcal{S}(\Sigma, \Omega)$, $\Sigma = \mathcal{R}(\Sigma, \tilde{h}_1, \dots, \tilde{h}_m)$.

12. Compute $^m\mu$, $^m\nu$, $^m\hat{g}^\alpha$ (Eqs (3.6-3.7)).

13. $$ \Omega = \langle f^1, \dots, f^{m_u}, \begin{Bmatrix} {}^m\hat{g}^0, & \dots, & {}^m\hat{g}^m \\ g^{m+1}, & \dots, & g^{m_w} \end{Bmatrix} \mid \Omega \rangle $$

14. **end while**

15. **if** $\deg_w(\Omega) == m_w$ **then**

16. **break** main loop.

17. **end if**

18. **if** $m_u \neq 0$ **then**

19. Run Algorithms 2 and 3 to compute $\hat{k}^m = s_x^m + r^m$.

20. **else**

21. $$ \hat{k}^m = 0. $$

22. **end if**

23. Set $k = 0$ and $\Omega_0 = \Omega$

24. **while** 1 { **This is the NESTED loop**} **do**

25. &nbsp;&nbsp;&nbsp;&nbsp;$k = k + 1$.

26. Reset $[\Sigma, \Omega_{k-1}] = \mathcal{A}(\Sigma, m, \Omega_{k-1})$.

27. Compute $^m\mu$, $^m\nu$, $^m\hat{g}^\alpha$ (Eqs (3.6-3.7)).

28. $$ \Omega^* = \Omega_{k-1}. $$

29. **if** ($k \le \hat{k}^m$)&($k \ge 2$) **then**

30. $$ \Omega^* = \Omega^* + \sum_{q=1}^m \sum_{i=1}^{m_u} \sum_{\alpha_1=0}^m \dots \sum_{\alpha_{k-1}=0}^m \operatorname{span}\{\mathcal{L}_{m[f^i]^{(\alpha_1, \dots, \alpha_{k-1})}} \nabla \tilde{h}_q\} $$

31. **end if**

32. $$ \Omega_k = \langle f^1, \dots, f^{m_u}, \begin{Bmatrix} {}^m\hat{g}^0, & \dots, & {}^m\hat{g}^m \\ g^{m+1}, & \dots, & g^{m_w} \end{Bmatrix} \mid \Omega^* \rangle. $$

33. **if** ($k - 1 \ge \hat{k}^m$)&($\Omega_k == \Omega_{k-1}$){ The highest unknown input degree of reconstructability is $m < m_w$ (system not canonic), and the observability codistribution is $\Omega_k$} **then**

34. Set $\mathcal{O} = \Omega_k$.

35. **Continue**=**false**.

36. **break** main and nested loop

37. **end if**

38. **if** $\deg_w(\Omega_k) > m$ **then**

39. Set $\Omega = \Omega_k$.

40. **break** nested loop

41. **end if**

42. **end while**

43. **end while**

44. **if** Continue **then**

45. &nbsp;&nbsp;&nbsp;&nbsp;[$\tilde{h}_1, \dots, \tilde{h}_{m_w}$] = $\mathcal{S}(\Sigma, \Omega)$.

46. Compute $\mu, \nu, \hat{g}^\alpha$, and $\tilde{\mathcal{O}}$ (Eqs (3.3-3.5)).

47. $$ \mathcal{O} = \langle f^1, \dots, f^{m_u}, \hat{g}^0, \hat{g}^1, \dots, \hat{g}^{m_w} \mid \Omega + \tilde{\mathcal{O}} \rangle $$

48. **end if**
**Algorithm 2:** The algorithm that returns the integer $s_x^m$ at Line 19 of Algorithm 1.

1: Set $\Sigma' = \Sigma$, $m, \tilde{h}_1, \dots, \tilde{h}_m$, as at Line 11 of Algorithm 1.

2: Set $k = 0$, $\Omega_0^m = \text{span}\{\nabla\tilde{h}_1, \dots, \nabla\tilde{h}_m\}$, and $^x\Omega_0^m = \Omega_0^m$.

3: **loop**

4: Set $k = k + 1$.

5: Reset $[\Sigma', \Omega_{k-1}^m] = \mathcal{A}(\Sigma', m, \Omega_{k-1}^m)$

6: $\Omega_k^m = \Omega_{k-1}^m + \mathcal{L}_{g^0} \Omega_{k-1}^m + \sum_{j=1}^m \mathcal{L}_{g_j} \Omega_{k-1}^m$.

7: $^x\Omega_k^m = \mathcal{D}_x(\Omega_k^m)$.

8: **if** $^x\Omega_k^m == ^x\Omega_{k-1}^m$ **then**

9: Set $s_x^m = k$ **then** exit.

10: **end if**

11: **end loop**

**Algorithm 3:** The algorithm that returns the integer $r^m$ at Line 19 of Algorithm 1.

1: Set $\Sigma' = \Sigma$ and $m$ as at Line 11 of Algorithm 1.

2: Set $k = 0$, $^m\Delta_0 = \text{span}\{f^1, \dots, f^{m_u}\}$.

3: **loop**

4: Set $k = k + 1$.

5: Reset $[\Sigma', ^m\Delta_{k-1}] = \mathcal{A}(\Sigma', m, ^m\Delta_{k-1})$

6: $^m\Delta_k = ^m\Delta_{k-1} + \sum_{\beta=0}^m {}^m\Delta_{k-1}^\beta$

7: **if** $^m\Delta_k == ^m\Delta_{k-1}$ **then**

8: Set $r^m = k - 1$ **then** exit.

9: **end if**

10: **end loop**

### 3.1.3 $\mathcal{R}(\Sigma, \lambda_1, \dots, \lambda_m)$ operation

Let us consider a system $\Sigma$ that is characterized by (2.2) and a set of $m < m_w$ scalar functions, $\lambda_1, \dots, \lambda_m$, such that rank($\mathcal{RM}(\lambda_1, \dots, \lambda_m)$) = $m$. The $\mathcal{R}(\Sigma, \lambda_1, \dots, \lambda_m)$ operation provides a new system that is obtained from $\Sigma$ by reordering the unknown inputs. In particular, they are reordered in such a way that the square submatrix that consists of the first $m$ columns of the unknown input reconstructability matrix from $\lambda_1, \dots, \lambda_m$, is non singular. The operation is denoted by $\mathcal{R}$ because it is a *Reordering* of the unknown inputs, as explained above. Its execution is immediate. It suffices to extract from $\mathcal{RM}(\lambda_1, \dots, \lambda_m)$ a set of $m$ independent columns. In this paper, we use the notation $\Sigma' = \mathcal{R}(\Sigma, \lambda_1, \dots, \lambda_m)$ to denote the new ordered system.

### 3.1.4 $\mathcal{A}(\Sigma, m)$ operation

Let us consider a system $\Sigma$ that is characterized by (2.2) and an integer $m < m_w$. This operation provides a new extended system, defined as follows. It is obtained by introducing a new extended state that includes the last $d = m_w - m$ unknown inputs. We have:

$$x \to [x^T, w_{m+1}, \dots, w_{m_w}]^T \qquad (3.1)$$

Starting from (2.2) we obtain the dynamics of the above extended state. We obtain a new system that still satisfies (2.2). It is still characterized by $m_u$ known inputs and $m_w$ unknown inputs. All the $m_u$ known inputs and the first $m$ unknown inputs coincide with the original ones. The last $d$ unknown inputs are $\dot{w}_{m+j}, j = 1, \dots, d$. Regarding the new vector fields that describe the dynamics we obtain:

$$f^i \to \begin{bmatrix} f^i \\ 0_d \end{bmatrix}, \quad g^0 \to \begin{bmatrix} g^0 + \sum_{l=m+1}^{m_w} g^l w_l \\ 0_d \end{bmatrix}, \qquad (3.2)$$

$$g^k \to \begin{bmatrix} g^k \\ 0_d \end{bmatrix}, \quad g^{m+j} \to \begin{bmatrix} 0_n \\ e^j \end{bmatrix},$$
with $i = 1, \dots, m_u$, $k = 1, \dots, m$, $j = 1, \dots, d$. In addition, $0_d$ and $0_n$ denote the zero $d$-column vector and the zero $n$-column vector, respectively, and $e^j$ denotes the $d$-column vector with the $j^{th}$ entry equal to 1 and the remaining $d-1$ entries equal to 0.

The operation is denoted by $\mathcal{A}$ because it consists of a state Augmentation and the consequent re-definition of all the key vector fields that characterize the new extended system, as specified above. In this paper, we use the notation $\Sigma' = \mathcal{A}(\Sigma, m)$ to denote the new extended system.

In Algorithms 1, 2, and 3, we also use this operation by including further outputs and further inputs. Specifically, we consider the following two cases:

1. $[\Sigma', \Omega'] = \mathcal{A}(\Sigma, m, \Omega)$, where the second output $(\Omega')$ is trivially the augmented codistribution obtained by extending all the covectors of $\Omega$ with $d$ zero entries.

2. $[\Sigma', \Delta'] = \mathcal{A}(\Sigma, m, \Delta)$, where the second output $(\Delta')$ is trivially the augmented distribution obtained by extending all the vectors of $\Delta$ with $d$ zero entries.

We can apply this operation multiple consecutive times (and this is the case of Algorithms 1, 2, and 3, as the operation appears in a loop). The resulting systems still satisfy (2.2) and they are all characterized by $m_u$ known inputs and $m_w$ unknown inputs. In [2], we called these extended systems Finite Unknown Inputs Extensions. In addition, in [2], we introduced the following definition:

**Definition 3.1** (Highest UI Degree of Reconstructability). Given the system in (2.2), the highest unknown input degree of reconstructability is the largest unknown input degree of reconstructability of all its finite unknown input extensions.

Finally, in [2] we called a system *Canonizable with respect to its unknown inputs* if its highest unknown input degree of reconstructability is equal to $m_w$.

## 3.2 Algorithm 1 for systems canonic with respect to the unknown inputs

When a system is canonic with respect to its unknown inputs, Algorithm 1 ends with the execution of Lines 45-47 (the boolean variable *Continue* remains set to "true"). The system could be directly in canonic form with respect to its unknown inputs or $\Omega$ was obtained in other parts of the algorithm (e.g., during the execution of the nested loop) and the system has been set in canonical form after the determination of one or more observable functions that are not among the outputs. The execution of Line 45 of Algorithm 1 provides the $m_w$ scalar functions $\tilde{h}_1, \dots, \tilde{h}_{m_w}$ ($[\tilde{h}_1, \dots, \tilde{h}_{m_w}] = S(\Sigma, \Omega)$). Then, the algorithm computes the following quantities: $\mu, \nu, \hat{g}^{\alpha}$, and the codistribution $\mathcal{O}$.

### 3.2.1 $\mu$ and $\nu$

The two-index tensor $\mu$ is defined as follows:

$$ \mu_j^i = \mathcal{L}_{g^i} \tilde{h}_j, \quad i, j = 1, \dots, m_w \qquad (3.3) $$

$$ \mu_0^0 = 1, \mu_0^i = 0, \mu_i^0 = \frac{\partial \tilde{h}_i}{\partial t} + \mathcal{L}_{g^0} \tilde{h}_i, \quad i = 1, \dots, m_w. $$

Note that the entries with $i, j = 1, \dots, m_w$ are the same entries of the unknown input reconstructability matrix from $\tilde{h}_1, \dots, \tilde{h}_{m_w}$, which is full rank. As a result, the tensor $\mu$ is also non singular and we denote by $\nu$ its inverse.

### 3.2.2 $\hat{g}^{\alpha}$

Starting from $\nu$ the algorithm builds the new vector fields $\hat{g}^0, \dots, \hat{g}^{m_w}$, defined as follows:

$$ \hat{g}^{\alpha} = \sum_{\beta=0}^{m_w} \nu_{\beta}^{\alpha} g^{\beta}, \quad \alpha = 0, 1, \dots, m_w. \qquad (3.4) $$
### 3.2.3 The codistribution $\tilde{\mathcal{O}}$

The codistribution $\tilde{\mathcal{O}}$ is defined as follows:

$$ \tilde{\mathcal{O}} := \sum_{q=1}^{m_w} \sum_{j=0}^{s+r} \sum_{\alpha_1=0}^{m_w} \dots \sum_{\alpha_j=0}^{m_w} \sum_{i=1}^{m_w} \mathrm{span} \left\{ \nabla \mathcal{L}_{[f^i]^{(\alpha_1, \dots, \alpha_j)}} \tilde{h}_q \right\}, \quad (3.5) $$

where the second sum is up to $s+r$, and the integers $s$ and $r$ are defined as follows:

* $s$ is the smallest integer such that $\Omega_s = \Omega_{s-1}$, where $\Omega_k$ is the codistribution at the $k^{th}$ step of the algorithm in (2.8), when computing:

$$ \langle g^0, g^1, \dots, g^{m_w} \mid \mathrm{span} \{ \nabla \tilde{h}_1, \dots, \nabla \tilde{h}_{m_w} \} \rangle $$

Note that, for TV systems, the Lie derivative operator along $g^0$, i.e., $\mathcal{L}_{g^0}$, must be replaced by $\dot{\mathcal{L}}_{g^0} := \mathcal{L}_{g^0} + \frac{\partial}{\partial t}$. Note that $s \le n - m_w + 1$.

* $r$ is the smallest integer such that $\Delta_{r+1} = \Delta_r$, where $\Delta_k$ is the distribution at the $k^{th}$ step of the algorithm in (2.10), when computing:

$$ \langle [\cdot] \mid \mathrm{span} \{ f^1, \dots, f^{m_u} \} \rangle $$

Note that $r \le n - 1$ (it is even $r \le n - \dim \{\mathrm{span} \{f^1, \dots, f^{m_u}\}\}$).

### 3.2.4 Final step

The last operation of Algorithm 1, when the system is canonical with respect to its unknown inputs, is the computation of the observability codistribution:

$$ \mathcal{O} = \langle f^1, \dots, f^{m_u}, \hat{g}^0, \hat{g}^1, \dots, \hat{g}^{m_w} \mid \Omega + \tilde{\mathcal{O}} \rangle $$

This is obtained by running the algorithm in (2.8) for the specific case (note that, for TV systems, the Lie derivative operator along $\hat{g}^0$, i.e., $\mathcal{L}_{\hat{g}^0}$, must be replaced by $\dot{\mathcal{L}}_{\hat{g}^0} := \mathcal{L}_{\hat{g}^0} + \frac{\partial}{\partial t}$). The convergence of the algorithm in (2.8) is attained at the smallest integer $j$ such that $\Omega_j = \Omega_{j-1}$. Note that $j$ cannot exceed $n - \dim(\Omega + \tilde{\mathcal{O}}) + 1$.

## 3.3 Algorithm 1 for systems that are not in canonical form

Let us back to the first part of Algorithm 1, at the beginning of the main loop, and let us consider now the case when $\mathrm{deg}_w(\Omega) < m_w$. The system is not in canonical form with respect to its unknown inputs. Algorithm 1 proceeds in this case with the loop at lines 9-14. The body of this loop is executed at least once. The set of scalar functions $\tilde{h}_1, \dots, \tilde{h}_m$ are obtained by executing the operation $\mathcal{S}(\Sigma, \Omega)$ defined in Section 3.1.2, and the unknown inputs are reordered accordingly ($\mathcal{R}(\Sigma, \tilde{h}_1, \dots, \tilde{h}_m)$, Section 3.1.3). Then, the algorithm computes the following quantities: $^m\mu$, $^m\nu$, and $^m\hat{g}^\alpha$.

### 3.3.1 $^m\mu$ and $^m\nu$

The two-index tensor $^m\mu$ is defined as follows:

$$ {}^m\mu_j^i = \mathcal{L}_{g^i} \tilde{h}_j, \quad i, j = 1, \dots, m \quad (3.6) $$

$$ {}^m\mu_0^0 = 1, \quad {}^m\mu_0^i = 0, \quad {}^m\mu_i^0 = \frac{\partial \tilde{h}_i}{\partial t} + \mathcal{L}_{g^0} \tilde{h}_i, \quad i = 1, \dots, m. $$

Note that the entries with $i, j = 1, \dots, m$ are the same entries of the matrix that consists of the first $m$ columns of the unknown input reconstructability matrix from $\tilde{h}_1, \dots, \tilde{h}_m$. By construction, this submatrix is full rank. As a result, the tensor $^m\mu$ is also non singular and we denote by $^m\nu$ its inverse.
### 3.3.2 $m\hat{g}^\alpha$

Starting from $^m\nu$ we build the new vector fields $^m\hat{g}^0, \dots, ^m\hat{g}^m$ defined as follows:

$$ ^m\hat{g}^\alpha = \sum_{\beta=0}^{m} {}^m\nu_\beta^\alpha g^\beta, \quad \alpha = 0, 1, \dots, m. \qquad (3.7) $$

Once all the quantities $^m\mu, ^m\nu$, and $^m\hat{g}^\alpha$ are computed, the loop at lines 9-14 of Algorithm 1 ends by computing

$$ \Omega = \langle f^1, \dots, f^{m_u}, \left\{ \begin{array}{ccc} {}^m\hat{g}^0, & \dots, & {}^m\hat{g}^m \\ g^{m+1}, & \dots, & g^{m_w} \end{array} \right\} \middle| \Omega \rangle. $$

(note that, for TV systems, the Lie derivative operator along $^m\hat{g}^0$, i.e., $\mathcal{L}_m\hat{g}^0$, must be replaced by $\mathcal{L}_m\hat{g}^0 := \mathcal{L}_m\hat{g}^0 + \frac{\partial}{\partial t}$). The loop is executed when the conditions given at Line 9 are satisfied. Note that, as $\deg_w(\Omega)$ cannot exceed $m_w$, the loop cannot be repeated more than $m_w - 1$ times.

After the execution of this loop, Algorithm 1 computes the integer $\hat{k}^m = s_x^m + r^m$, where the integers $r^m$ and $s_x^m$ are defined similarly to the integers $r$ and $s$ introduced in Section 3.2.3:

* $s_x^m$ is the smallest integer such that $^x\Omega_{s_x^m}^m = {}^x\Omega_{s_x^m-1}^m$, where $^x\Omega_k^m$ is the codistribution at the $k^{th}$ step computed by Algorithm 2. Note that $s_x^m \le n - m + 1$.

* $r^m$ is the smallest integer such that $^m\Delta_{r^m+1} = {}^m\Delta_{r^m}$, where $^m\Delta_k$ is the distribution at the $k^{th}$ step computed by Algorithm 3. Note that $r^m \le n - 1$ (it is even $r^m \le n - \dim\{\operatorname{span}\{f^1, \dots, f^{m_u}\}\}$).

We have $\hat{k}^m = s_x^m + r^m \le 2n - m + 1 - \dim\{\operatorname{span}\{f^1, \dots, f^{m_u}\}\} \le 2n - m$. Below, we provide a description of Algorithms 2 and 3.

#### Algorithm 2

It computes the codistribution $^x\Omega^m$. This algorithm uses a new operation denoted by $\mathcal{D}_x(\Omega)$. The input of this operation is a codistribution that is defined in a given augmented space (in the algorithm, before this operation, the $\mathcal{A}(\Sigma, m)$ operation is applied). The operation $\mathcal{D}_x(\Omega)$, first detects a basis of $\Omega$. As $\Omega$ is an integrable codistribution, it detects a set of scalar functions: $\theta_1, \dots, \theta_D$, where $D$ is the dimension of $\Omega$. In other words, $\Omega = \operatorname{span}\{\nabla\theta_1, \dots, \nabla\theta_D\}$, where $\nabla$ is the gradient with respect to the new extended state (i.e., the state defined at the last $\mathcal{A}(\Sigma, m)$ operation). The output of the operation $\mathcal{D}_x(\Omega)$, is $^x\Omega = \operatorname{span}\{\partial_x\theta_1, \dots, \partial_x\theta_D\}$, where $\partial_x$ is the gradient with respect to the original state, i.e., the state of $\Sigma$ defined at the first line of Algorithm 2.

The initialization step sets the codistribution equal to the span of the gradients of the selected functions $\tilde{h}_1, \dots, \tilde{h}_m$. As a result, the dimension of $\Omega_0^m$ is $m$. Then, each iteration of the loop executes the following operations:

1. System augmentation, as explained in Section 3.1.4.

2. Computation of $\Omega_k^m$ by adding to $\Omega_{k-1}^m$ the term $\mathcal{L}_{g^0}^m \Omega_{k-1}^m + \sum_{j=1}^m \mathcal{L}_{g^j} \Omega_{k-1}^m$.

3. Computation of $^x\Omega_k^m = \mathcal{D}_x(\Omega_k^m)$.

The convergence of Algorithm 2 occurs at the smallest integer $j$ such that $^x\Omega_j^m = {}^x\Omega_{j-1}^m$ and $j \le n - m + 1$. We set $s_x^m = j$. Note that, because of the presence of the state augmentation in the loop (Line 5 of Algorithm 2), proving the validity of the above convergence property is very demanding (see Appendix E of [2] and in particular the proof of Proposition 11 in that appendix).
### Algorithm 3

It computes the distribution $^m\Delta$. The initialization step sets the distribution equal to the span of the vector fields $f^1, \dots, f^{m_u}$. Then, the recursive step adds to $^m\Delta_{k-1}$ the term $\sum_{\beta=0}^{m_w} [^m\Delta_{k-1}]^\beta$. The convergence is attained at the smallest integer $j$ such that $^m\Delta_j = ^m\Delta_{j-1}$ and $j \le n - \dim\{\operatorname{span}\{f^1, \dots, f^{m_u}\}\} + 1$. We set $r^m + 1 = j$. Note that, because of the presence of the state augmentation in the loop (Line 5 of Algorithm 3), proving the validity of the above convergence property is non trivial (see Section 6.2.1 of [2] and in particular the proof of Proposition 1 given in Appendix D of [2]).

Note that, the system augmentation performed by the $\mathcal{A}$ operation executed by Algorithms 2 and 3, does not reset the system in Algorithm 1. The system defined in Algorithm 1 can only be reset (augmented) by the operation $\mathcal{A}$ that is executed at Line 26 of Algorithm 1. The execution of Algorithm 2 and Algorithm 3 only returns the integer $\hat{k}^m = s_x^m + r^m$. All remaining quantities remain unchanged.

#### 3.3.3 The nested loop

Algorithm 1 continues by computing the codistributions $\Omega_k$, $k = 0, 1, \dots$. For positive $k$, this is obtained by executing the nested loop (Lines 24-42). Each iteration, first augments the system (Line 26 of Algorithm 1). Note that, this operation updates the set of vector fields $g^0, g^1, \dots, g^m$, in accordance with (3.2). As a result, it is necessary to recompute $^m\hat{g}^0, ^m\hat{g}^1, \dots, ^m\hat{g}^m$ by using (3.7). This is carried out at the next line (Line 27). Then, Algorithm 1 computes $\Omega_k$ (Lines 29-32). The last part of the loop contains two if statements (Line 33, and Line 38). Each of them, when the condition is met, interrupts the loop. In [2], we proved that at least one of these two conditions is met in at most

$$ \hat{k}^m + n - \dim(\Omega_0) \le \hat{k}^m + n - m \le 3n - 2m $$

steps, where $n$ is the dimension of the original state, i.e., the state of the system before entering the nested loop. Note that, in [2] (Section 7), we actually proved that one of the two if conditions is met in at most $n - m$. This because the first $\hat{k}^m$ steps were included in the initialization step (this is less efficient because the condition of the second if can be met during the execution of these first $\hat{k}^m$ steps).

When the condition of the first if statement (Line 33) is met, not only the loop is interrupted but also the execution of Algorithm 1 terminates. In this case, the system is not canonic with respect to its unknown inputs. Algorithm 1 provides the entire observable codistribution. In addition, the augmented system is the one with the highest unknown input degree of reconstructability (see Definition 3.1). When the loop is interrupted by the second if statement (Line 38), the execution of the nested loop is interrupted. Again, as $m \le m_w$, the condition at Line 38 cannot be met more than $m_w - 1$ times.

## 3.4 Illustration by the case study

We execute Algorithm 1 to obtain the observability codistribution in the two scenarios of our case study defined in Section 2.6. Note that, in both scenarios the system is in canonical form. The systems investigated in Chapters 6 and 7 are not in canonical form and, to see the implementation of Algorithm 1 on a system that is not in canonical form, the reader is addressed to Section 6.1 or to Section 7.1.

### First scenario

We have:

$$ g^0 = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}, \quad f^1 = \begin{bmatrix} \cos(\theta - \phi) \\ \frac{\sin(\theta - \phi)}{\rho} \\ 0 \end{bmatrix}, \quad g^1 = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}, $$
and $h_1 = \phi - \theta$. We run Algorithm 1. Line 2 provides $\Omega = \text{span}\{\nabla h_1\} = \text{span}\{[0, 1, -1]\}$, and $\text{deg}_w(\Omega) = 1$, as the unknown input observability matrix from $h_1$ is the $1 \times 1$ matrix equal to $\mathcal{L}_{g^1}h_1 = -1 \neq 0$ (its rank is 1). The system is in canonical form with respect to its unknown input. Line 45 trivially selects $\tilde{h}_1 = h_1$. Then, Line 46 provides:

$$ \mu = \nu = \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix}, $$

obtained from (3.3). In addition, from (3.4) we obtain: $\hat{g}^0 = \nu_0^0 g^0 + \nu_1^0 g^1 = g^0 = [0, 0, 0]^T$, and $\hat{g}^1 = \nu_0^1 g^0 + \nu_1^1 g^1 = -g^1 = [0, 0, -1]^T$. Finally, we need to compute $\tilde{\mathcal{O}}$ from (3.5). This requires, first of all, to compute the two integers $s$ and $r$:

* $s$ is the smallest integer such that $\Omega_s = \Omega_{s-1}$, where $\Omega_k$ is the codistribution at the $k^{th}$ step of the algorithm in (2.8), when computing:

$$ \langle g^1 | \text{span}\{\nabla \tilde{h}_1\} \rangle = \left\langle \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} \right\rangle \left| \text{span}\{[0, 1, -1]\} \right\rangle $$

We immediately obtain $s = 1$.

* $r$ is the smallest integer such that $\Delta_{r+1} = \Delta_r$, where $\Delta_k$ is the distribution at the $k^{th}$ step of the algorithm in (2.10), when computing:

$$ \langle [\cdot] | \text{span}\{f^1\} \rangle = \left\langle [\cdot] \left| \text{span} \left\{ \begin{bmatrix} \frac{\cos(\theta - \phi)}{\rho} \\ 0 \end{bmatrix} \right\} \right\rangle \right\rangle. $$

We obtain $r = 1$ (see the calculation details below).

## Detail of computation of $r$

The algorithm in (2.10), for the specific case, provides: $\Delta_0 = \text{span}\{f^1\}$. We need to compute $[\Delta_0]^\alpha$, for $\alpha = 0, 1$. On the other hand, $[f^1]^0 = \nu_0^0[g^0, f^1] + \nu_1^0[g^1, f^1]$ vanishes (the first term vanishes because $g^0$ is null, the second term vanishes because $\nu_1^0 = 0$). $[f^1]^1 = \nu_0^1[g^0, f^1] + \nu_1^1[g^1, f^1] = -[g^1, f^1] = \chi^1$, with $\chi^1 := [\sin(\theta - \phi), -\cos(\theta - \phi)/\rho, 0]^T$. Hence, $\Delta_1 = \text{span}\{f^1, \chi^1\}$. Finally, $[\chi^1]^0$ vanishes and $[\chi^1]^1 = -f^1$. Therefore, $\Delta_2 = \Delta_1$ and $r+1 = 2$.

Let us compute $\tilde{\mathcal{O}}$. From (3.5), and from the above computation, we obtain that the generators of $\tilde{\mathcal{O}}$ are $\mathcal{L}_{f^1}h_1 = \frac{\sin(\theta-\phi)}{\rho}$ and $\mathcal{L}_{\chi^1}h_1 = -\frac{\cos(\theta-\phi)}{\rho}$. Hence:

$$ \tilde{\mathcal{O}} = \text{span} \left\{ \begin{bmatrix} -\frac{\sin(\theta - \phi)}{\rho^2}, & -\frac{\cos(\theta - \phi)}{\rho}, & \frac{\cos(\theta - \phi)}{\rho} \\ \frac{\sin(\theta - \phi)}{\rho} & 0 & 0 \end{bmatrix}, \begin{bmatrix} \frac{\cos(\theta - \phi)}{\rho^2}, & -\frac{\sin(\theta - \phi)}{\rho}, & \frac{\sin(\theta - \phi)}{\rho} \end{bmatrix} \right\} $$

The final step of Algorithm 1 is the execution of Line 47, namely:

$$ \mathcal{O} = \left\langle \begin{bmatrix} \frac{\cos(\theta - \phi)}{\rho} \\ \frac{\sin(\theta - \phi)}{\rho} \\ 0 \end{bmatrix}, \begin{bmatrix} 0 \\ 0 \\ -1 \end{bmatrix} \right\rangle \left| \text{span}\{[0, 1, -1]\} + \tilde{\mathcal{O}} \right\rangle $$

as $\hat{g}^0$ vanishes and $\hat{g}^1 = -g^1$. The algorithm in (2.8) converges at the first step ($\Omega_1 = \Omega_0$), and $\mathcal{O} =$

$$ \text{span} \left\{ [0, 1, -1], \begin{bmatrix} -\frac{\sin(\theta - \phi)}{\rho^2}, & -\frac{\cos(\theta - \phi)}{\rho}, & \frac{\cos(\theta - \phi)}{\rho} \end{bmatrix} \right\}. $$

The above codistribution contains all the observability properties of our system. We postpone a discussion about them to Section 4.4.1.
## Second scenario

With respect to the previous scenario, we now have:

$$g^1 = \begin{bmatrix} \cos(\theta - \phi) \\ \frac{\sin(\theta - \phi)}{\rho} \\ 0 \end{bmatrix}, \quad f^1 = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}$$

We run Algorithm 1. Line 2 provides $\Omega = \text{span}\{\nabla h_1\} = \text{span}\{[0, 1, -1]\}$, and $\text{deg}_w(\Omega) = 1$, as the unknown input observability matrix from $h_1$ is the $1 \times 1$ matrix equal to $\mathcal{L}_{g^1} h_1 = \frac{\sin(\theta - \phi)}{\rho} \neq 0$ (its rank is 1). The system is in canonical form with respect to its unknown input. Line 45 trivially selects $\tilde{h}_1 = h_1$. Then, Line 46 provides:

$$\mu = \begin{bmatrix} 1 & 0 \\ 0 & \sin(\theta - \phi)/\rho \end{bmatrix}, \quad \nu = \begin{bmatrix} 1 & 0 \\ 0 & \rho/\sin(\theta - \phi) \end{bmatrix}, \qquad (3.8)$$

obtained from (3.3). In addition, from (3.4) we obtain: $\tilde{g}^0 = \nu_0^0 g^0 + \nu_1^0 g^1 = g^0 = [0, 0, 0]^T$, and $\tilde{g}^1 = \nu_0^1 g^0 + \nu_1^1 g^1 = [\rho \frac{\cos(\theta - \phi)}{\sin(\theta - \phi)}, 1, 0]^T$. Finally, we need to compute $\tilde{\mathcal{O}}$ from (3.5). This requires, first of all, to compute the two integers $s$ and $r$:

* $s$ is the smallest integer such that $\Omega_s = \Omega_{s-1}$, where $\Omega_k$ is the codistribution at the $k^{th}$ step of the algorithm in (2.8), when computing:

$$\langle g^1 | \text{span}\{\nabla \tilde{h}_1\} \rangle = \left\langle \begin{bmatrix} \cos(\theta - \phi) \\ \frac{\sin(\theta - \phi)}{\rho} \\ 0 \end{bmatrix} \middle| \text{span}\{[0, 1, -1]\} \right\rangle$$

By a direct computation, we obtain $\Omega_2 = \Omega_1 = \text{span}\{\nabla h_1, \nabla \frac{\sin(\theta - \phi)}{\rho}\}$ and, consequently, $s = 2$.

* $r$ is the smallest integer such that $\Delta_{r+1} = \Delta_r$, where $\Delta_k$ is the distribution at the $k^{th}$ step of the algorithm in (2.10), when computing:

$$\langle [\cdot]^r | \text{span}\{f^1\} \rangle = \langle [\cdot]^r | \text{span}\left\{ \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} \right\} \rangle.$$

We obtain $r = 2$. We do not provide the detail of the computation. We only mention that, in this case, the generators of $\Delta$ are: $f^1$, $[\rho, -\frac{\cos(\theta - \phi)}{\sin(\theta - \phi)}, 0]^T$, and $[\rho \frac{\cos(\theta - \phi)}{\sin(\theta - \phi)}, -\frac{\cos^2(\theta - \phi)}{\sin^2(\theta - \phi)}, 0]^T$

Finally, Line 46 computes $\tilde{\mathcal{O}}$. From (3.5), by a direct computation we obtain:

$$\tilde{\mathcal{O}} = \text{span}\left\{\nabla \frac{\cos(\theta - \phi)}{\sin(\theta - \phi)}\right\}.$$

The final step of Algorithm 1 is the execution of Line 47, namely:

$$\mathcal{O} = \left\langle \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}, \begin{bmatrix} \rho \frac{\cos(\theta - \phi)}{\sin(\theta - \phi)} \\ 1 \\ 0 \end{bmatrix} \middle| \text{span}\{[0, 1, -1]\} \right\rangle,$$

as $\tilde{\mathcal{O}} \subseteq \text{span}\{[0, 1, -1]\}$, and $\tilde{g}^0$ is null. The algorithm in (2.8) converges at the first step ($\Omega_1 = \Omega_0$), and:

$$\mathcal{O} = \text{span}\{[0, 1, -1]\}.$$

The above codistribution contains all the observability properties of our system. We postpone a discussion about them to Section 4.4.1.
# Chapter 4

# Systematic procedure to perform
the identifiability analysis

This section is devoted to the problem of the identifiability of the time-varying parameters, or, equivalently, to the reconstructability of the system unknown inputs. As we will see, this problem is strongly related with the problem of state observability in the presence of unknown inputs (whose solution was provided in Chapter 3). We provide the general solution to obtain the reconstructability of the unknown inputs.

This section consists of four subsections. In 4.1, we extend a preliminary result provided in [2] and we obtain the answer to our problem when the state that characterizes our system is observable. In 4.2, we remind the reader of the concept of continuous symmetry introduced in [1]. In 4.3, we extend this concept to the unknown input vector and we provide the complete answer to our problem. In practice, we obtain the general systematic procedure that solves the problem of the identifiability of the time-varying parameters. The procedure is Algorithm 4. Finally, in 4.4, we illustrate all these concepts by using our case study.

## 4.1 Identifiability when the state is observable

In [2] we proved that when the state that characterizes the system is observable, if the system is canonic with respect to its unknown inputs, then all the unknown inputs can be reconstructed (Theorem 4 in [2]). It also holds the viceversa and, consequently, when the state is observable, we obtain a complete answer to our problem which is given by the following theorem:

**Theorem 4.1.** If the state of the system characterized by (2.2) is observable, the unknown inputs can be reconstructed if and only if the system is canonic with respect to its unknown inputs.

Note that, checking if a system is canonic with respect to its unknown input is immediate. It suffices to compute the rank of the reconstructability matrix from any set of generators of the observability codistribution (and Algorithm 1 provides a set of generators of this codistribution). Therefore, Theorem 4.1 provides a full answer to the problem of unknown input reconstruction when the state is observable.

*Proof*. Based on the result stated by Theorem 4 in [2], it remains to prove that if the unknown inputs can be reconstructed then the system is canonic with respect to its unknown inputs.

We proceed by contradiction. We assume that the system is not canonic with respect to its unknown inputs.

As the state is observable, all its components belong to the observation space and we have: $O = \text{span}\{\nabla x_1, \dots, \nabla x_n\}$. As we assumed that the system is not canonic, $\text{deg}_w(O) < m_w$. Therefore, the matrix $RM(x_1, \dots, x_n)$ has rank smaller than $m_w$. On the other hand, by an explicit computation we obtain
$$ \mathcal{RM}(x_1, \dots, x_n) = \begin{bmatrix} [g^1]_1 & [g^2]_1 & \dots & [g^{m_w}]_1 \\ [g^1]_2 & [g^2]_2 & \dots & [g^{m_w}]_2 \\ \vdots & \vdots & \ddots & \vdots \\ [g^1]_n & [g^2]_n & \dots & [g^{m_w}]_n \end{bmatrix} = \begin{bmatrix} g^1 & g^2 & \dots & g^{m_w} \end{bmatrix} $$

where $[g^i]_j$ is the $j^{th}$ component of the vector field $g^i$ ($i = 1, \dots, m_w$ and $j = 1, \dots, n$). Therefore, as rank($\mathcal{RM}(x_1, \dots, x_n)$) < $m_w$, the vectors $g^1, g^2, \dots, g^{m_w}$ are linearly dependent. We denote by $\hat{n}$ a non trivial vector that belongs to the null space of the above matrix. Let us consider the dynamics:

$$ \dot{x} = g^0 + \sum_{k=1}^{m_w} f^k u_k + \sum_{j=1}^{m_w} g^j w_j. $$

For any $w = [w_1, w_2, \dots, w_{m_w}]$, we obtain the same dynamics by transforming the UI as follows:

$$ w \rightarrow w' = w + \hat{n} $$

As result, $w$ cannot be distinguished from $w'$ and, consequently, cannot be reconstructed. ◀

As we said, when the state is observable, the above theorem provides a full answer to our problem. It remains to find the solution when the state is not observable. Unfortunately, we need a complete different approach. In other words, in the presence of the state unobservability, we cannot exploit the result stated by Theorem 4.1. We provide an explanation of this fact.

Let us suppose that the state is not observable. In order to exploit the result stated by Theorem 4.1, we need to characterize our system by an observable state. Algorithm 1 provides the entire observable space. We could introduce a new state whose components are the generators of this space, i.e., the generators of $O^1$. This new state is certainly observable and its components are expressed in terms of the old state. To use Theorem 4.1, we need to describe our system by using this new state starting from the description that is characterized by the old state. This means to obtain a set of equations of the same format of (2.2) where the state $x$ is the new state. In other words, the dynamics of $x$, precisely as in (2.2), must be expressed only in terms of the components of the new $x$ (and not in terms of the old state). On the other hand, when we try to obtain this description, we often encounter two fundamental difficulties that make it critical, and even useless:

* **Critical issues:** This description cannot be determined automatically. In some cases, its determination can be a very demanding and laborious task, based on enormous inventiveness (e.g., see Section 6.3 where, although the state is small, the determination of this description requires ability and inventiveness).

* **Uselessness:** In many cases, obtaining this description requires a re-definition of the unknown inputs. In other words, this task cannot be performed by maintaining the same original unknown inputs (see the discussion of the second scenario of our case study in Section 4.4.2, and our applications in Sections 6.3 and 7.3). When the final goal is precisely the reconstruction of the original unknown inputs, performing this task (i.e., describing the system by an observable state) can be useless.

The above discussion showed that, in general, we cannot exploit the result stated by Theorem 4.1 to deal with systems characterized by an unobservable state. In this case, we need a different approach. This method will be introduced in Section 4.3. It is based on the concept of continuous symmetry introduced in [48] and here summarized in Section 4.2.

<footnote n="1">¹Clearly, with generators, here we mean the scalar functions whose gradients generate $O$.</footnote>
## 4.2 Unobservability of the State, State Symmetries and Indistinguishable States

In [48], and more exhaustively in [1], we introduced the concept of system symmetry in the presence of unobservability. When we are in the presence of unobservability, the dimension of the observability codistribution ($\mathcal{O}$) is smaller than the state dimension, $n$. Let us denote it by $n_o$. The null space of $\mathcal{O}$ is not empty and it is a distribution (the orthogonal distribution, from now on denoted by $\mathcal{O}^{\perp}$), with dimension $n - n_o$. In [1, 48], we defined symmetry any vector field that belongs to $\mathcal{O}^{\perp}$. In Section 4.6 of [1] we showed that the generators of $\mathcal{O}^{\perp}$ are the generators of a Lie algebra associated to a Lie group. Given a state $x$, all its indistinguishable states can be obtained by the action of this Lie group on $x$. Specifically, in [1], we showed that, starting from a system symmetry $\xi \in \mathcal{O}^{\perp}$ and a given state $\bar{x}$ we can generate a set of states which are indistinguishable from $\bar{x}$ by solving the following differential equation:

$$ \begin{cases} \frac{dx}{d\tau} = \xi(x(\tau)) \\ x(0) = \bar{x} \end{cases} \qquad (4.1) $$

In particular, if this equation admits solution on the interval $\tau \in [0, T]$, then all the states $x(\tau)$, obtained by integrating the above equation up to $\tau$, are indistinguishable from each other. The above procedure defines a one parameter Lie group. The parameter is $\tau$ and its action on the state $\bar{x}$ returns the state $x(\tau)$:

$$ \bar{x} \to x(\tau) \qquad (4.2) $$

As usual, when dealing with Lie groups, it is very useful to refer to the action of the group when the parameter takes infinitesimal values. The above transformation becomes:

$$ \bar{x} \to \bar{x} + \epsilon\xi \qquad (4.3) $$

where $\epsilon$ is an infinitesimal<fnref n="2" />. Given a symmetry $\xi$ of the system in (2.2), the two states $\bar{x}$ and $\bar{x} + \epsilon\xi$ are indistinguishable. Note that we have a symmetry if and only if the distribution $\mathcal{O}^{\perp}$ is non trivial, i.e., if and only if the dimension of $\mathcal{O}$ is strictly smaller than $n$. Hence, we can say that the state is observable if and only we do not have symmetries. The symmetries can be computed automatically (it suffices to compute the nullspace of $\mathcal{O}$). Starting from them, we can also try to obtain the indistinguishable states, by solving the differential equation in (4.1). On the other hand, this last step cannot be performed automatically. When it can be performed analytically, it provides further useful insights on the ODE model. Finally, note that when the dimension of $\mathcal{O}^{\perp}$ is strictly larger the 1, the states indistinguishable from a given state $\bar{x}$ can be obtained by the action of a multi parameters Lie group (instead of a simple one-parameter Lie group). We do not provide details on this Lie group. We only mention that a set of generators of its associated Lie algebra is precisely the set of generators of $\mathcal{O}^{\perp}$ (note that, as $\mathcal{O}$ is an integrable codistribution, because of the Frobenius Theorem, $\mathcal{O}^{\perp}$ is involutive, i.e., closed with respect to the Lie brackets). We also emphasize that, in this case of multiple independent symmetries, we can still consider each symmetry independently and obtaining 1-dimensional indistinguishable sets by the action of the corresponding one-parameter Lie group (i.e., by solving the differential equation in (4.1)). We get a 1-dimensional indistinguishable set for each symmetry.

In Section 4.3, we extend the concept of symmetry to the unknown input vector and we obtain a full answer to the identifiability problem.

## 4.3 Identifiability for any system

We start by running Algorithm 1. By doing this, we automatically obtain the observability codistribution, $\mathcal{O}$. In addition, the algorithm provides a new system, which is an unknown input extension of the original system with the highest unknown input degree of recostructability. From now on, we refer to this system. Its structure still satisfies (2.2). We denote the dimension of its state by $n$ and we denote its unknown input degree of recostructability (which is the highest) by

<footnote n="2">²The reader is addressed to [49, 50, 51] for the definition of infinitesimal (or infinitesimal number) and its use in calculus.</footnote>
m. Note that, a set of m observable functions, $\tilde{h}_1, \dots, \tilde{h}_m$, is also available, after the execution of Algorithm 1. The new system could coincide with the original system, meaning that the UI degree of reconstructability of the original system is already the highest. In this case, Algorithm 1 returns the original system where only the order of the unknown inputs could have changed. Indeed, when $m < m_w$, the system returned by Algorithm 1 is ordered in accordance with the $\mathcal{R}$ operation, defined in Section 3.1.3. In other words, the rank of the first m columns of $\mathcal{RM}(\tilde{h}_1, \dots, \tilde{h}_m)$ is m.

In this section we provide the expression of all the symmetries that characterize the unknown input vector of the system returned by Algorithm 1. In other words, we provide the expression of all the vectors $\chi = [\chi_1, \dots, \chi_{m_w}]^T$ such that:

$$w' = w + \epsilon \chi \quad (4.4)$$

is indistinguishable from w. In Section 4.3.1, we provide the expression of a first set of such symmetries that exist when the system is not canonic with respect to its unknown inputs (i.e., $m < m_w$). They will be denoted by $^c\chi$, where c stands for system non-Canonicity. In Section 4.3.2 we provide the expression of a second set of such symmetries that may exist when the state is unobservable. They will be denoted by $^u\chi$, where u stands for state Unobservability.

Note that, the characterization of unknown input symmetry provided by (4.4) was also very recently introduced in [45, 46]. In these works, this definition was used by working in the augmented space, i.e., by characterizing the system with a state that includes the unknown inputs together with their time derivatives up to a given order. As we mentioned in Chapter 1, using this extended state has two limitations: (i) the computational burden can easily become prohibitive, and (ii) some of the symmetries could be actually spurious because we did not include enough time derivatives of the unknown inputs (and no criterion to establish which is the maximum order to be included was provided). In addition, working with the augmented state does not allow us to distinguish between the two aforementioned categories of unknown input symmetries (i.e., to distinguish between $^c\chi$ and $^u\chi$) and can also make more difficult the physical interpretation of a given symmetry.

### 4.3.1 Unknown input symmetries due to the system non canonicity

We provide the expression of a set of symmetries of the unknown input vector that certainly exist when the system is not canonic with respect to its unknown inputs. In particular, this set consists of $m_w - m$ independent symmetries.

**Theorem 4.2.** Let us consider the system returned by Algorithm 1 and let us suppose that $m < m_w$ (system not canonic with respect to the unknown inputs). For this system, the following $m_w - m$ column vectors:

$$^c\chi^i = \begin{bmatrix} 0_m \\ e_i \end{bmatrix}, \quad i = 1, \dots, m_w - m, \quad (4.5)$$

with $0_m$ the zero m-column vector and $e_i$ the unit vector of dimension $m_w - m$ with zero everywhere except that the $i^{th}$ component which is 1, are $m_w - m$ independent symmetries of the unknown input vector.

*Proof.* The proof is given in Appendix A.

### 4.3.2 Possible unknown input symmetries due to the state unobservability

When the state is unobservable, we have one or more nonnull symmetries of the state. We have the following fundamental result.

**Theorem 4.3.** Let us consider the system returned by Algorithm 1, and let us denote by $\xi \in \mathcal{O}^\perp$ a given symmetry of the state. The $m_w$- column vector $^u\chi$, with components:
$$ \left\{ \begin{array}{ll} {}^{u}\chi_k = -\sum_{i=1}^{m} {}^{m}\nu_k^{i} (\xi_i^0 + \sum_{j=1}^{m} \xi_j^j w_j) & k = 1, \dots, m \\ {}^{u}\chi_k = 0 & k = m + 1, \dots, m_w \end{array} \right. \qquad (4.6) $$

where, for any $i = 1, \dots, m$:

$$ \xi_i^\alpha = \begin{cases} \left( \frac{\partial}{\partial x} \left( \mathcal{L}_{g^0} \tilde{h}_i + \frac{\partial \tilde{h}_i}{\partial t} \right) \right) \cdot \xi & \alpha = 0 \\ \left( \frac{\partial}{\partial x} \mathcal{L}_{g^\alpha} \tilde{h}_i \right) \cdot \xi & \alpha = 1, \dots, m \end{cases} \qquad (4.7) $$

is a symmetry of w.

When the system is canonic with respect to its unknown inputs ($m = m_w$), the unknown input symmetry becomes:

$$ {}^{u}\chi_k = -\sum_{i=1}^{m_w} \nu_k^i \left( \xi_i^0 + \sum_{j=1}^{m_w} \xi_j^j w_j \right), \quad k = 1, \dots, m_w. $$

*Proof.* The proof is given in Appendix B.

We provide the following remarks.

**Remark 4.1.** Given a symmetry $\xi$ of the state, the corresponding symmetry ${}^u\chi$ in (4.6) can be trivial, i.e., can be a null vector. This is certainly the case when all the quantities $\xi_i^\alpha$ in (4.7) vanish. If for all the symmetries of the state the corresponding ${}^u\chi$ vanish, the first $m$ component of $w$ can be reconstructed. This means that there are cases where the first $m$ unknown inputs can be reconstructed even in the presence of state unobservability. Hence, it is better to call the vector ${}^u\chi$ in (4.6) a potential symmetry of the unknown input vector.

**Remark 4.2.** We have a potential symmetry of the unknown input vector for any symmetry of the state.

**Remark 4.3.** If the state is observable, there are not state symmetries and the first $m$ unknown inputs can be reconstructed. This is consistent with the result stated by Theorem 4.1. Indeed, when the state is observable we only have the unknown input symmetries that characterize the last $m_w - m$ unknown inputs. These symmetries disappear if and only if $m = m_w$ (i.e., the system is canonic with respect to the UIs).

### 4.3.3 The systematic procedure

The systematic procedure to automatically obtain the identifiability of all the time-varying parameters (or the reconstructability of all the unknown inputs) is Algorithm 4. It is based on the results stated by Theorems 4.2 and 4.3. It starts by running Algorithm 1 on the system under investigation (Line 1). This provides the observability codistribution $\mathcal{O}$, the highest UI degree of reconstructability with respect to the unknown inputs, $m$, and a new system $\Sigma'$ which is an unknown input extension of the original system with the highest UI degree of reconstructability. Note that, the identifiability analysis regards the new system $\Sigma'$. Hence, we investigate the reconstructability of the unknown inputs of this system. These unknown inputs may differ from the original ones for the following two reasons:

1. Their order may have been changed.

2. Some of them could be the time derivatives (of a given order) of the original ones.

Note that, in any case, we also obtain the reconstructability properties of the original unknown inputs. Indeed, if a given UI of $\Sigma'$ is the time derivative of a given UI of the original system, it certainly means that the original UI was included in the state. To check its identifiability it suffices to verify if its differential belongs to $\mathcal{O}$. Algorithm 4 will tell us if the corresponding UI of $\Sigma'$ is identifiable (i.e., if the time derivative of a given order of the original UI is identifiable).

Line 3 of Algorithm 4 provides a set of $m_w - m$ symmetries of the unknown input, which are the $m_w - m$ unit vectors defined in (4.5). This set is empty if and only if the system is canonic with respect to its unknown inputs ($m = m_w$).
Lines 5 and 6 compute the orthogonal distribution $O^\perp$ and a basis of it. In particular, its generators are $\xi^1, \dots, \xi^S$. Then, the loop at Lines 7-9 computes the corresponding symmetries of the unknown input vector, i.e., $^u\chi^1, \dots, ^u\chi^S$, on the basis of Theorem 4.3.

At this point, we have obtained all the symmetries that characterize the unknown input vector. They are the following $S + m_w - m$ vectors:

$$ ^c\chi^1, \dots, ^c\chi^{m_w-m}, ^u\chi^1, \dots, ^u\chi^S. $$

Finally, the loop at Lines 10-16 checks the reconstructability of all the unknown inputs $w_1, \dots, w_{m_w}$. Specifically, the $j^{th}$ unknown input ($w_j$) can be reconstructed if and only if:

$$ ^c\chi_j^1 = \dots = ^c\chi_j^{m_w-m} = ^u\chi_j^1 = \dots = ^u\chi_j^S = 0, $$

i.e., the $j^{th}$ components of all the unknown input symmetries, vanish, simultaneously.

Note that, when $S = 0$ (and this is the case when the state that characterizes $\Sigma'$ is observable) the first $m$ unknown inputs, $w_1, \dots, w_m$, can be reconstructed. If in addition the system is canonic with respect to its unknown inputs, all the unknown inputs can be reconstructed. Note that the observability of the state is a sufficient (but not a necessary) condition to have the reconstructability of the first $m$ unknown inputs. Indeed, there are cases where the state is unobservable (i.e., $S > 0$) but all the symmetries $^u\chi^1, \dots, ^u\chi^S$ vanish, simultaneously (e.g., this occurs in the first scenario of our case study, as discussed in Section 4.4.3).

**Algorithm 4:** The algorithm that builds the symmetries of the unknown input vector and check the reconstructability of all its components ($w_j, j = 1, \dots, m_w$).

1: Run Algorithm 1 and compute $O, m, \Sigma'$.

2: **if** $m < m_w$ **then**

3: Set $^c\chi^1, \dots, ^c\chi^{m_w-m}$, as in (4.5).

4: **end if**

5: Compute $O^\perp$.

6: Set $O^\perp = \text{span}\{\xi^1, \dots, \xi^S\}$.

7: **for** $i_s = 1 : S$ **do**

8: Set $\xi = \xi^{i_s}$, and compute the corresponding $^u\chi$ (from (4.6) and (4.7)), and denote it by $^u\chi^{i_s}$.

9: **end for**

10: **for** $j = 1 : m_w$ **do**

11: **if** $^c\chi_j^1 = \dots = ^c\chi_j^{m_w-m} = ^u\chi_j^1 = \dots = ^u\chi_j^S = 0$ **then**

12: $w_j$ can be reconstructed

13: **else**

14: $w_j$ cannot be reconstructed

15: **end if**

16: **end for**

## 4.4 Illustration by the case study

We illustrate all the concepts introduced by this section by applying them to our case study. Note that, for educational purposes, we provide the indistinguishable states and we also obtain a description of our system by an observable state. On the other hand, the implementation of Algorithm 4 requires neither the determination of the indistinguishable states nor a description of the system by an observable state. It only requires the computation of all the independent symmetries of the state, i.e., a basis of the orthogonal distribution $O^\perp$.

### 4.4.1 Symmetries and Indistinguishable states

#### First scenario

From the expression of $O$, provided in Section 3.4, we obtain:
$$ \mathcal{O}^{\perp} = \operatorname{span} \left\{ \begin{bmatrix} 0 \\ 1 \\ 1 \end{bmatrix} \right\}, $$

which is the generator of the rotations around the vertical axis. This agrees with our intuition as all the measurements (the linear velocity, $v$, and the bearing angle, $\beta$) are invariant under any rotation of the system around the vertical axis.

Let us consider a given state, $[\bar{\rho}, \bar{\phi}, \bar{\theta}]^T$ and let us compute its indistinguishable states. This computation cannot be performed automatically, as it requires to solve a differential equation. On the other hand, it is only provided for educational purposes. From Equation (4.1), we need to solve the following differential equation:

$$ \left\{ \begin{aligned} \frac{d\rho}{d\tau} &= 0 \\ \frac{d\phi}{d\tau} &= 1 \\ \frac{d\bar{\theta}}{d\tau} &= 1 \\ \rho(0) &= \bar{\rho}, \quad \phi(0) = \bar{\phi}, \quad \theta(0) = \bar{\theta} \end{aligned} \right. $$

We obtain that all the states:

$$ \begin{bmatrix} \bar{\rho} \\ \bar{\phi} + \tau \\ \bar{\theta} + \tau \end{bmatrix} $$

are indistinguishable. For a given $\tau$, the above state is obtained by rotating the state $[\bar{\rho}, \bar{\phi}, \bar{\theta}]^T$ by an angle of magnitude $\tau$ around the vertical axis.

## Second scenario

From the expression of $\mathcal{O}$, provided in Section 3.4, we obtain:

$$ \mathcal{O}^{\perp} = \operatorname{span} \left\{ \begin{bmatrix} 0 \\ 1 \\ 1 \end{bmatrix}, \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} \right\}. $$

The first generator is again the generator of the rotations around the vertical axis. It agrees with our intuition as also the measurements of the angular speed $\omega$ are invariant under any rotation of the system around the vertical axis. The new symmetry is the generator of the absolute scale transformations. Its presence is not surprising. When the linear speed is unknown, all the available measurements are angular measurements ($\omega$ and $\beta$) and there is no metric information. In other words, all the measurements are also invariant to the scale.

Let us consider a given state, $[\bar{\rho}, \bar{\phi}, \bar{\theta}]^T$ and let us compute the indistinguishable states obtained by using the second generator. Equation (4.1) becomes:

$$ \left\{ \begin{aligned} \frac{d\rho}{d\tau} &= 1 \\ \frac{d\phi}{d\tau} &= 0 \\ \frac{d\bar{\theta}}{d\tau} &= 0 \\ \rho(0) &= \bar{\rho}, \quad \phi(0) = \bar{\phi}, \quad \theta(0) = \bar{\theta} \end{aligned} \right. $$

We obtain that all the states:

$$ \begin{bmatrix} \bar{\rho} + \tau \\ \bar{\phi} \\ \bar{\theta} \end{bmatrix} $$

are indistinguishable. As $\rho$ is the only metric component of the state (the other two, $\theta$ and $\phi$, are angles), the above state can be regarded as a scale transformation of $[\bar{\rho}, \bar{\phi}, \bar{\theta}]^T$ (in particular, it is rescaled by the factor $\frac{\bar{\rho}+\tau}{\bar{\rho}}$). Note that, a symmetry is defined up to a multiplying factor (for any vector in $\mathcal{O}^{\perp}$ also the vector obtained by multiplying it by a given scalar is in $\mathcal{O}^{\perp}$). Hence, the transformation on the state defined by $\xi = [1, 0, 0]^T$ can also be defined by the same symmetry multiplied by $\rho$. In other words, we could consider the symmetry: $[\rho, 0, 0]^T$. By
using this symmetry, Equation (4.1) provides, $\frac{d\rho}{d\tau} = \rho$, whose solution is $\rho(\tau) = \bar{\rho}e^{\tau}$. This second option makes more readable that it is a scale transformation (in particular, the scale factor is $e^{\tau}$).

### 4.4.2 Characterization by an observable state

This step cannot be performed automatically but, as the system is very simple, its execution is trivial. In addition, as we mentioned, this step is not required for the implementation of Algorithm 4.

#### First scenario

The dimension of the observability codistribution is 2 and two independent generators are: $\rho$, and $\psi := \theta - \phi$. Hence, a possible observable state is precisely $x = [\rho, \psi]^T$. From (2.13), we easily obtain:

$$ \begin{bmatrix} \dot{\rho} = v \cos \psi \\ \dot{\psi} = \omega - \frac{v}{\rho} \sin \psi \end{bmatrix} \qquad (4.8) $$

where $\omega(=w)$ is the unknown input and $v(=u)$ the known input.

#### Second scenario

The dimension of the observability codistribution is 1 and one generator is $\psi$. Hence, a possible observable state is precisely $x = [\psi]$. From (2.13), we easily obtain: $\dot{\psi} = \omega - \frac{v}{\rho} \sin \psi$, where $v(=w)$ is the unknown input and $\omega(=u)$ the known input. Note that, by using $w$, we cannot obtain a description of the system in terms of the solely observable state, because the expression of the dynamics needs the quantity $\rho$. We must redefine the unknown input. A possible choice is $\tilde{w} := \frac{w}{\rho} = \frac{v}{\rho}$. The system dynamics become:

$$ \left[ \begin{matrix} \dot{\psi} = \omega - \tilde{w} \sin \psi \end{matrix} \right. \qquad (4.9) $$

### 4.4.3 Identifiability

We use Algorithm 4 to obtain the identifiability of the time-varying parameters for our case study. In particular, for the first scenario, the unknown time-varying parameter is the angular speed $\omega(t)$ while, in the second scenario, is the linear speed $v(t)$. In both scenarios, the system is canonic with respect to its unknown input. Hence, we can only have the unknown input symmetries due to the state unobservability ($^u\chi$).

#### First scenario

We have a single symmetry of the state that is $\xi = [0, 1, 1]^T$. We compute the corresponding symmetry of the unknown input, $^u\chi$, by using (4.6). We need, first of all, to compute the terms $\xi_i^\alpha$ by using (4.7), for $\alpha = 0, 1$, and $i = 1$. We have $\tilde{h}_1 = \tilde{h} = \phi - \theta$. As $g^0$ is null, $\mathcal{L}_{g^0} \tilde{h} = 0$. Hence:

$$ \xi_1^0 = 0 $$

Let us compute $\xi_1^1$. We have $g^1 = [0, 0, 1]^T$. We obtain: $\mathcal{L}_{g^1} \tilde{h} = -1$. As a result,

$$ \xi_1^1 = 0 $$

From (4.6), we obtain that the symmetry of the unknown input vanishes.

Therefore, we do not have non trivial symmetries of the unknown input. This means that the unknown input can be reconstructed (or, in other words, the time-varying parameter $\omega(t)$ is identifiable).

Note that the same result is obtained by using Theorem 4.1 on the system in (4.8). Indeed, the state that characterizes this system is observable and the system is canonic with respect to its unknown input. Hence, the unknown input can be reconstructed.
## Second scenario

We have two independent symmetries: $\xi^1 = [0, 1, 1]^T$, and $\xi^2 = [\rho, 0, 0]^T$. Note that, for the second symmetry, on the basis of the remark at the end of Section 4.4.1, we use $[\rho, 0, 0]^T$ instead of $[1, 0, 0]^T$. To compute the corresponding $^u\chi$, we need the terms $\xi_1^0$ and $\xi_1^1$. For both symmetries, we need to compute $\mathcal{L}_{g^0}\tilde{h}_1$ and $\mathcal{L}_{g^1}\tilde{h}_1$. On the other hand, $g^0$ is null, $g^1 = [\cos(\theta - \phi), \frac{\sin(\theta - \phi)}{\rho}, 0]^T$, and $\tilde{h}_1 = \tilde{h} = \phi - \theta$. Hence, we obtain:

$$ \mathcal{L}_{g^0}\tilde{h} = 0, \quad \mathcal{L}_{g^1}\tilde{h} = \frac{\sin(\theta - \phi)}{\rho}. $$

Let us consider the first symmetry: $\xi^1 = [0, 1, 1]^T$. From (4.7), and the above expressions of $\mathcal{L}_{g^0}\tilde{h}$ and $\mathcal{L}_{g^1}\tilde{h}$, we obtain:

$$ \xi_1^0 = \xi_1^1 = 0. $$

From (4.6), we obtain that the corresponding symmetry of the unknown input vanishes, i.e., $^u\chi^1 = 0$.

Let us consider the second symmetry $\xi^2 = [\rho, 0, 0]^T$. From (4.7), and the above expressions of $\mathcal{L}_{g^0}\tilde{h}$ and $\mathcal{L}_{g^1}\tilde{h}$, we obtain:

$$ \xi_1^0 = 0, \quad \xi_1^1 = \frac{\partial \mathcal{L}_{g^1}\tilde{h}}{\partial x} \xi^2 = -\frac{\sin(\theta - \phi)}{\rho}. $$

We use (4.6) to compute $^u\chi$. Note that, in this case $m = m_w = 1$, and, from (3.8), $^m\nu_1^1 = \nu_1^1 = \frac{\rho}{\sin(\theta - \phi)}$. We obtain:

$$ ^u\chi^2 = v. $$

This means that, although we have two symmetries of the state, we have a single non trivial symmetry of the unknown input. The unknown input cannot be reconstructed (or, in other words, the time-varying parameter $v(t)$ is not identifiable).

In this case, by applying Theorem 4.1 on (4.9), we only conclude that the unknown input $\tilde{w} = \frac{v(t)}{\rho(t)}$ can be reconstructed. However, the Theorem does not allow us to take a conclusion about the identifiability of the original time-varying parameter $w = v(t)$.
# Chapter 5

## Indistinguishable states and parameters and minimal missing information to achieve observability/identifiability

In this chapter we focus our attention to the systems that are unobservable and/or unidentifiable. We start by providing a general procedure that generate indistinguishable states and parameters (Section 5.1). Then, in Section 5.2 we exploit the determination of these indistinguishable states and parameters to characterize the minimal missing information to achieve observability and/or identifiability. In particular, we provide a preliminary result that allows us to obtain the above characterization in the presence of a single symmetry.

### 5.1 Indistinguishable states and parameters

Algorithm 4 can be executed automatically and provides a full answer to the problem of the identifiability of the time-varying parameters (or unknown input reconstruction). It uses the results stated by Theorems 4.2 and 4.3, which are based on the powerful concept of infinitesimal transformation. In this section, we want to discuss the case of finite transformations. In other words, we want to extend the same analysis provided in Section 4.2 to the unknown inputs. In that case, the infinitesimal transformation was given by Equation (4.3) and the finite transformation was given by Equation (4.2), where $x(\tau)$ is the solution of the differential equation in (4.1). Note that, the solution of the differential equation in (4.1) cannot be obtained automatically, and, in many cases, it can only be determined numerically. On the other hand, its determination is unnecessary to perform both the observability and the identifiability analysis. The same remarks hold for the differential equations given in this section (i.e., Equations (5.5) and (5.7)), which provide the finite transformations to build indistinguishable unknown inputs and initial states, simultaneously. On the other hand, their solution can provide useful insights on a given ODE model (in Sections 6.4 we provide the analytical solution for a model in the framework of HIV dynamics and, in Section 7.4, we provide the analytical solution for a model in the framework of Covid-19 dynamics).

Let us consider the system $\Sigma'$ returned by Algorithm 1. It still satisfies Equation (2.2). Let us assume that the initial time is $t = 0$ and the initial state is $x(0) = x_0$. The state at time $t$ is obtained by solving the following differential equation (where all the quantities are the ones that characterize the system $\Sigma'$ returned by Algorithm 1):

$$ \left\{ \begin{aligned} \dot{x}(t) &= g^0(x, t) + \sum_{k=1}^{m_u} f^k(x, t)u_k(t) + \sum_{j=1}^{m_w} g^j(x, t)w_j(t) \\ x(0) &= x_0 \end{aligned} \right. \qquad (5.1) $$
We consider any choice of the known inputs, $u_1(t), \dots, u_{m_u}(t)$, such that the above differential equation admits a solution in the time interval $[0, T]$.

Now, we introduce the state $x'(t)$ obtained by solving the following differential equation:

$$ \begin{cases} \dot{x}'(t) = g^0(x', t) + \sum_{k=1}^{m_u} f^k(x', t)u_k(t) \\ x'(0) = x'_0 \end{cases} + \sum_{j=1}^{m_w} g^j(x', t)w'_j(t) \quad (5.2) $$

We emphasize that the state $x'(t)$ is obtained by solving the same differential equation that provides $x(t)$ but with the following two differences:

* The initial state is $x'_0$ instead of $x_0$.

* The unknown input is $w'(t) = [w'_1(t), \dots, w'_{m_w}(t)]^T$ instead of $w(t) = [w_1(t), \dots, w_{m_w}(t)]^T$ (and this second difference involves the unknown input at any time in the considered time interval and not simply the initial time).

We are interested in detecting all the initial states $x'_0$ and all the unknown input functions $w'(t)$ such that, independently of the choice of the known input functions, the states $x'(t)$ produce exactly the same outputs of the states $x(t)$, i.e.:

$$ h_i(x(t), t) = h_i(x'(t), t), \quad i = 1, \dots, p, $$

for any $t \in [0, T]$.

On the basis of Theorems 4.2 and 4.3, we have two sets of solutions, which are provided by the two types of unknown input symmetries, $^c\chi$ and $^u\chi$. Specifically, on the basis of the result stated by Theorem 4.2, we obtain the following first set of solutions:

$$ x'_{0} = x_{0}, \quad w'(t) = w(t) + \epsilon \, {}^{c}\chi^{i}, \quad i = 1, \dots, m_w - m \quad (5.3) $$

On the basis of the result stated by Theorem 4.3, we obtain the following second set of solutions:

$$ x'_{0} = x_{0} + \epsilon \xi^{i}, \quad w'(t) = w(t) + \epsilon \, {}^{u}\chi^{i}, \quad i = 1, \dots, S \quad (5.4) $$

where $\xi^i \in \mathcal{O}^\perp$ is the symmetry of the state that provides $^u\chi^i$, in accordance with (4.6) (i.e., $^u\chi^i$ is obtained by setting $\xi = \xi^i$ in (4.7)).

Starting from the infinitesimal transformations given in (5.3) and in (5.4), we can compute the corresponding finite transformations. For the $m_w - m$ infinitesimal transformations given in (5.3), the corresponding finite transformations are the following. Regarding $x'_0$, it trivially remains the same for all of them, i.e.:

$$ x'_{0} = x_{0} $$

Regarding the unknown inputs, we have $m_w - m$ independent $w' = w'(t, \tau)$, which are the solutions of the following differential equations:

$$ \begin{cases} \frac{dw'}{d\tau} = {}^c\chi^i \\ w'(t, 0) = w(t), \end{cases} \quad (5.5) $$

$i = 1, \dots, m_w - m$. In the above equation, the time $t$ is a fixed parameter. In other words, to obtain $w'(t, \tau)$, we first set the time to the desired value $t$. Due to the trivial expression of $^c\chi$ in (4.5) (and in particular its independence of the specific case), this differential equation can be solved once forever and provides:

$$ \begin{cases} w'_{j}(t, \tau) = w_{j}(t), & j \neq m + i \\ w'_{j}(t, \tau) = w_{j}(t) + \tau, & j = m + i \end{cases} $$

For the $S$ infinitesimal transformations given in (5.4), the corresponding finite transformations now involve simultaneously the initial state and the unknown inputs. To analytically compute the corresponding finite transformations, we need the following condition to be satisfied:
$$ \left[ \xi, g^0 + \sum_{i=1}^{m_u} f^i u_i + \sum_{j=1}^{m_w} g^j w_j \right] + \sum_{j=1}^{m} g^j {}^u \chi_j = 0, \qquad (5.6) $$

where the square brackets are the Lie brackets. This condition expresses the commutativity between the temporal shift and the shift generated by the considered symmetry.

If the condition in (5.6) is honoured, the result of many consecutive applications of the S infinitesimal transformations in (5.4) can be analytically computed. Specifically, for the $i^{th}$ transformation in (5.4) we obtain the new $x'_0 = x'_0(\tau)$, and $w' = w'(t, \tau)$ by solving the following differential equations:

$$ \left\{ \begin{aligned} \frac{dx'}{d\tau'} &= \xi^i(x'(t, \tau)) \\ \frac{dw'}{d\tau'} &= {}^u\chi^i(x'(t, \tau), w'(t, \tau)) \\ x'(t, 0) &= x(t), \quad w'(t, 0) = w(t), \end{aligned} \right. \qquad (5.7) $$

i = 1, ..., S. Again, in the above equation, the time t is a fixed parameter. In other words, to obtain $x'(t, \tau)$, and $w'(t, \tau)$, we first set the time to desired value t. If, for any $t \in [0, T]$, the above equation admits a solution in the interval $[0, T]$, then, for any $\tau \in [0, T]$, the unknown input vector $w'(t, \tau)$ is indistinguishable from the unknown input vector $w(t)$. Note that, in the above equation, the value of $x(t)$, which appears at the initial condition $x'(t, 0) = x(t)$, is the solution of (5.1) at the time t. From $x'(t, \tau)$ we obtain the finite transformation of the initial state, which is the value of this function at the initial time, i.e., $x'_0 = x'_0(\tau) = x'(0, \tau)$.

The solution of (5.7) is obtained by first determining $x'(t, \tau)$ (i.e., by first solving $\frac{dx'}{d\tau} = \xi^i(x'(t, \tau))$, with initial condition $x'(t, 0) = x(t)$), and then by using it in $\frac{dw'}{d\tau} = {}^u\chi^i(x'(t, \tau), w'(t, \tau))$, which becomes a differential equation of only $w'(t, \tau)$. Note that the solution of $\frac{dx'}{d\tau} = \xi^i(x'(t, \tau))$, with initial condition $x'(t, 0) = x(t)$, is precisely the solution of (4.1) with $\bar{x} = x(t)$ and for the same state symmetry (i.e., $\xi = \xi^i$).

## 5.2 Minimal missing information to achieve observability/identifiability

We discuss a fundamental issue that arises every time we are in the presence of unobservability and/or unidentifiability. We wish to answer to the following practical question: What is the minimal external information (external to the available measurements of the system outputs) required to make observable the state and identifiable all the model parameters?

In particular, we would like an automatic answer, which holds for any system, as in the case of the observability and identifiability analysis that are obtained by simply running Algorithm 1 and 4.

Unfortunately, we do not have a general answer yet (this will be the matter of a future investigation). For the moment, we have a preliminary result that provides an answer to our question when the state is characterized by a single symmetry (i.e., the orthogonal distribution $O^\perp$ has a single generator). The key is that, in the presence of a single symmetry, we are able to determine all the indistinguishable states by solving the differential equations in (4.1) (or even to determine all the indistinguishable states and unknown inputs by solving (5.7)).

Let us denote the components of the true state at time t by $x_1(t), x_2(t), \dots, x_n(t)$. As the system is characterized by a single symmetry, all the indistinguishable states are:

$$ \left[ \begin{array}{l} x'_1(t, \tau) = s_1(x_1(t), x_2(t), \dots, x_n(t), \tau) \\ x'_2(t, \tau) = s_2(x_1(t), x_2(t), \dots, x_n(t), \tau) \\ \vdots \\ x'_n(t, \tau) = s_n(x_1(t), x_2(t), \dots, x_n(t), \tau) \end{array} \right. \qquad (5.8) $$

where $s_1, \dots, s_n$ are specific functions determined by solving the differential equations in (4.1) (or (5.7), in which case the system in (5.8) includes the components of the state and the unknown inputs).

We know that, by using all the available measurements, we cannot determine $x_1(t), x_2(t), \dots, x_n(t)$, but we can determine $x'_1(t, \tau), x'_2(t, \tau), \dots, x'_n(t, \tau)$. The latter are related to the former by the equation system in (5.8). In particular, this system can be regarded as a system of n
equations in the $n+1$ unknowns, which are the $n$ components of the true state and $\tau$. Note that we have such a system of equations for any $t$ in our time interval. On the other hand, while the components of the state are in general time dependent, $\tau$ is independent of $t$. Let us suppose that we add a measurement, at a single time $t^*$, on one of the unobservable components of the true state. Without loss of generality, let us assume that this is the last component of the state, (i.e., we are measuring $x_n(t^*)$). The fact that $x_n$ is unobservable, means that the function $s_n$ depends on $\tau$. Now we consider our equation system in (5.8) at $t = t^*$. This system consists of $n$ equations in $n$ unknowns, which are the first $n-1$ components of the true state at time $t^*$, and $\tau$. In general, this allows us to obtain the $n$ unknowns, and in particular, to obtain $\tau$. But once $\tau$ has been determined, our equation system in (5.8) for any $t \neq t^*$ consists of $n$ equations in $n$ unknowns, which are the $n$ components of the true state at time $t$. In general, this allows us to obtain the $n$ unknowns, i.e., to obtain the components of the state at any time $t$.

In Section 5.3 we provide a trivial application of this procedure for the case study. In Section 6.7 we will see a more complex application of this procedure for a model that characterizes the HIV dynamics.

## 5.3 Illustration by the case study

We determine the set of indistinguishable unknown inputs for our case study. Note that, when the system parameters are identifiable, this set is trivial and consists of a single (the true) unknown input. Hence, we only consider the second scenario for which we obtained a single non trivial symmetry of the unknown input vector. Specifically, in Section 4.4.3 we obtained $\xi = [\rho, 0, 0]^T$ and $^u\chi = v$. By an explicit computation, it is possible to check that the condition in (5.6) is honoured. The differential equation in (5.7) becomes:

$$ 
\left\{ 
\begin{aligned}
\frac{d\rho'}{d\tau} &= \rho'(t, \tau) \\
\frac{d\phi'}{d\tau} &= 0 \\
\frac{d\theta'}{d\tau} &= 0 \\
\frac{dv'}{d\tau} &= v'(t, \tau) \\
\rho'(t, 0) &= \rho(t), \quad \phi'(t, 0) = \phi(t), \quad \theta'(t, 0) = \theta(t) \\
v'(t, 0) &= v(t)
\end{aligned}
\right.
\qquad (5.9) $$

and its solution is:

$$ 
\left[
\begin{array}{l}
\rho'(t, \tau) = \rho(t) e^{\tau}, \\
\phi'(t, \tau) = \phi(t), \\
\theta'(t, \tau) = \theta(t), \\
v'(t, \tau) = v(t) e^{\tau},
\end{array}
\right.
\qquad (5.10) $$

which is a scale transformation with factor $e^\tau$. As expected, the unknown input $v(t)$ can only be reconstructed up to a scale. By computing the above $\rho'(t, \tau)$ at $t=0$ we obtain that the initial state is only observable up to the same scale, $\rho'_0 = \rho'(0, \tau) = \rho_0 e^\tau$.

Note that, by using $\xi = [1, 0, 0]^T$ instead of $\xi = [\rho, 0, 0]^T$, we would obtain the same result. The corresponding $^u\chi$ is $\frac{v}{\rho}$ and the differential equation in (5.7), becomes:

$$ 
\left\{ 
\begin{aligned}
\frac{d\rho'}{d\sigma} &= 1 \\
\frac{d\phi'}{d\sigma} &= 0 \\
\frac{d\theta'}{d\sigma} &= 0 \\
\frac{dv'}{d\sigma} &= \frac{v'(t, \sigma)}{\rho'(t, \sigma)} \\
\rho'(t, 0) &= \rho(t), \quad \phi'(t, 0) = \phi(t), \quad \theta'(t, 0) = \theta(t) \\
v'(t, 0) &= v(t)
\end{aligned}
\right.
$$

where we used $\sigma$ instead of $\tau$. The solution is:

$$ 
\left[
\begin{array}{l}
\rho'(t, \sigma) = \rho(t) + \sigma, \\
\phi'(t, \sigma) = \phi(t), \\
\theta'(t, \sigma) = \theta(t), \\
v'(t, \sigma) = \frac{v(t)}{\rho(t)}(\rho(t) + \sigma),
\end{array}
\right.
$$
which coincides with the previous solution by setting $\sigma = \rho(t)(e^{\tau} - 1)$.

We conclude by remarking that, in accordance with the above results, the ratio $\frac{v}{\rho}$ is unaffected by the above scale transformation. Indeed, we have:

$$ \frac{v}{\rho} \rightarrow \frac{ve^{\tau}}{\rho e^{\tau}} = \frac{v}{\rho}. $$

This means that this ratio is observable. This result is consistent with the result that we obtain by applying Theorem 4.1 on the system in (4.9), which is characterized by an observable state. As this system is also canonic with respect to its unknown input, it means that the unknown input $\tilde{w} = \frac{v}{\rho}$ can be reconstructed.

We conclude this section by characterizing the minimal missing information to achieve observability and identifiability. Clearly, due to the simplicity of the system, we obtain an answer by proceeding intuitively. Specifically, we know that to have observability and identifiability we need to know the absolute scale. This is provided by a single measurement of the distance $\rho$ or by a single measurement of the speed $v$. We show that we obtain the same result by applying the procedure provided in Section 5.2. Note that we are allowed to use this procedure because the system is characterized by a single symmetry.

Equation (5.8) becomes, in this case, (5.10). In order to determine the true state and the true parameter, we must determine the value of $\tau$. This can be determined by knowing $\rho$ or $v$ at a single time ($t^*$). In the first case, $\tau$ is obtained from the first equation in (5.10), which provides $\rho(t^*)e^{\tau}$. The knowledge of $\rho$ at $t^*$, i.e. $\rho(t^*)$, makes immediate the determination of $e^{\tau}$, and then $\tau$. Once $\tau$ has been determined, we determine the remaining quantities at any $t \in \mathcal{I}$.
# Chapter 6

## HIV infection

We investigate the observability and the identifiability properties of a simple ODE model widely used to describe HIV dynamics in HIV-infected patients with antiretroviral treatment [3, 35, 40, 41, 42]. The model is characterized by the following three equations and the following two outputs:

$$ 
\left\{
\begin{aligned}
\dot{T}_U &= \lambda - \rho T_U - \eta(t) T_U V \\
\dot{T}_I &= \eta(t) T_U V - \delta T_I \\
\dot{V} &= N \delta T_I - cV \\
y &= [V, T_U + T_I],
\end{aligned}
\right.
\qquad (6.1) $$

where, $T_U$ is the concentration of uninfected cells, $T_I$ the concentration of infected cells, and $V$ the viral load. The model is also characterized by the time-varying parameter $\eta(t)$ and the five constant parameters $\lambda, \rho, \delta, N, c$. They are defined as follows:

* $\eta(t)$ is the infection rate, which is a function of the antiviral treatment efficacy.

* $\lambda$ is the source rate of uninfected cells.

* $\rho$ is the death rate of uninfected cells.

* $\delta$ is the death rate of infected cells.

* $N$ is the average number of virions produced by a single infected cell during its lifetime.

* $c$ is the clearance rate of free virions.

All these parameters are assumed to be unknown.

To proceed, we need, first of all, to introduce a state that includes both the time-varying quantities (i.e., $T_U, T_I, V$) and the constant parameters (i.e., $\lambda, \rho, \delta, N, c$). We set:

$$ x = [T_U, T_I, V, \lambda, \rho, \delta, N, c]^T \qquad (6.2) $$

The system is directly a special case of (2.2). In particular, $m_u = 0, m_w = 1, p = 2, h_1(x) = V, h_2(x) = T_U + T_I,$

$$ g^0 = \begin{bmatrix} \lambda - \rho T_U \\ -\delta T_I \\ N \delta T_I - cV \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}, \quad g^1 = \begin{bmatrix} -T_U V \\ T_U V \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} \qquad (6.3) $$
## 6.1 Observability analysis

We run Algorithm 1 to obtain the observability codistribution. Line 2 provides:

$$ \begin{aligned} \Omega &= \operatorname{span}\{\nabla h_1, \nabla h_2\} \\ &= \operatorname{span}\{[0, 0, 1, 0, 0, 0, 0], [1, 1, 0, 0, 0, 0, 0]\}. \end{aligned} $$

$\deg_w(\Omega) = 0$ ($< m_w$), as $\mathcal{L}_{g^1}h_1 = \mathcal{L}_{g^1}h_2 = 0$ (and the rank of the unknown input reconstructability matrix from $h_1, h_2$ vanishes). The system is not in canonical form with respect to its unknown input. We enter the loop at lines 9-14.

Line 11 provides no functions ($m = 0$). Line 12 provides: $^m\mu_0^0 = ^m\nu_0^0 = 1$ obtained from (3.6) (these tensors have a single entry as the indices only take the value 0). In addition, from (3.7) we obtain: $^m\hat{g}^0 = ^m\nu_0^0g^0 = g^0$. Line 13 computes

$$ \Omega = \langle \left\{ \begin{array}{c} g^0 \\ g^1 \end{array} \right\} \mid \operatorname{span}\{\nabla h_1, \nabla h_2\} \rangle, $$

as $^m\hat{g}^0 = g^0$. By a direct computation, we obtain:

$$ \Omega = \operatorname{span}\{\nabla h_1, \nabla h_2, \nabla\mathcal{L}_{g^0}h_1, \nabla\mathcal{L}_{g^0}h_2\} = \operatorname{span}\{\nabla h_1, \nabla h_2, \nabla h_3, \nabla h_4\} $$

where we adopted the following notation:

$$ h_3 := \mathcal{L}_{g^0}h_1, \quad h_4 := \mathcal{L}_{g^0}h_2 $$

Now, $\deg_w(\Omega) = 1 = m_w$. Hence, the second condition at Line 9 is not met and the loop is interrupted. In addition, as the condition at Line 15 is met, also the main loop is interrupted. The execution of Algorithm 1 continues with Line 45. The operation at this line ($[\tilde{h}_1, \dots, \tilde{h}_{m_w}] = \mathcal{S}(\Sigma, \Omega)$) sets

$$ \tilde{h}_1 = h_4 = \mathcal{L}_{g^0}h_2 = \lambda - \rho T_U - \delta T_I $$

(we would obtain the same final result by setting $\tilde{h}_1 = h_3 = \mathcal{L}_{g^0}h_1$). Then, Line 46 provides:

$$ \begin{aligned} \mu &= \begin{bmatrix} 1 & 0 \\ T_I \delta^2 - \rho(\lambda - T_U \rho) & -T_U V(\delta - \rho) \end{bmatrix}, \\ \nu &= \begin{bmatrix} 1 & 0 \\ \frac{T_I \delta^2 + T_U \rho^2 - \lambda \rho}{T_U V(\delta - \rho)} & -\frac{1}{T_U V(\delta - \rho)} \end{bmatrix}, \end{aligned} \quad (6.4) $$

obtained from (3.3). In addition, from (3.4) we obtain:

$$ \hat{g}^0 = \begin{bmatrix} -\frac{\delta(T_I \delta - \lambda + T_U \rho)}{\rho(T_I \delta - \lambda + T_U \rho)} \\ \frac{\delta - \rho}{\delta - \rho} \\ NT_I \delta - Vc \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}, \quad \hat{g}^1 = \begin{bmatrix} \frac{1}{\delta - \rho} \\ -\frac{1}{\delta - \rho} \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$

Finally, from (3.5) we obtain that, when $m_u = 0$, the codistribution $\tilde{\mathcal{O}}$ vanishes. The final step of Algorithm 1 is the execution of Line 47, namely:

$$ \mathcal{O} = \langle \hat{g}^0, \hat{g}^1 | \operatorname{span}\{\nabla h_1, \nabla h_2, \nabla h_3, \nabla h_4\} \rangle. $$

The algorithm in (2.8) converges at the third step ($\Omega_3 = \Omega_2$), and $\mathcal{O} =$

$$ \operatorname{span}\{\nabla h_1, \nabla h_2, \nabla h_3, \nabla h_4, \nabla\mathcal{L}_{\hat{g}^0}h_3, \nabla\mathcal{L}_{\hat{g}^1}h_3, \nabla\mathcal{L}_{\hat{g}^2}h_3\}. $$
Its dimension is 7, which is smaller than the dimension of the state in (6.2). Hence, the state is not observable. In particular, regarding the constant parameters, it is immediate to check that two of them ($\delta$ and $N$) are unidentifiable (it suffices to show that their gradient, i.e., the two covectors [0, 0, 0, 0, 0, 1, 0, 0] and [0, 0, 0, 0, 0, 0, 1, 0], does not belong to $\mathcal{O}$). This result contradicts the result available in the state of the art (e.g., see Section 6.2 of [3] or Section 3.3 of [35]). In Section 6.6, we discuss the above result by referring to the same data set adopted in [35].

## 6.2 Symmetries and indistinguishable states

Once we have $\mathcal{O}$, the symmetries are immediately obtained by computing the orthogonal distribution. We have:

$$ \mathcal{O}^{\perp} = \operatorname{span} \left\{ \begin{bmatrix} T_I \delta \\ -T_I \delta \\ 0 \\ 0 \\ 0 \\ \delta(\delta - \rho) \\ N\rho \\ 0 \end{bmatrix} \right\}. \qquad (6.5) $$

As for the case study, starting from the generators of $\mathcal{O}^{\perp}$ we compute the indistinguishable states (although the determination of the indistinguishable states is not requested by Algorithm 4). We need to solve the differential equation in (4.1). Given a state $[\bar{T}_U, \bar{T}_I, \bar{V}, \bar{\lambda}, \bar{\rho}, \bar{\delta}, \bar{N}, \bar{c}]^T$, the solution of this differential equation provides states which are indistinguishable. For the specific case, by using the generator in (6.5) for $\xi$, Equation (4.1) becomes:

$$ \left\{ \begin{aligned} \frac{dT_U}{d\tau} &= T_I \delta \\ \frac{dT_I}{d\tau} &= -T_I \delta \\ \frac{d\bar{V}}{d\tau} &= 0 \\ \frac{d\bar{\lambda}}{d\tau} &= 0 \\ \frac{d\bar{\rho}}{d\tau} &= 0 \\ \frac{d\bar{\delta}}{d\tau} &= \delta(\delta - \rho) \\ \frac{d\bar{N}}{d\tau} &= N\rho \\ \frac{d\bar{c}}{d\tau} &= 0 \\ T_U(0) &= \bar{T}_U, \quad T_I(0) = \bar{T}_I, \quad V(0) = \bar{V}, \quad \lambda(0) = \bar{\lambda}, \\ \rho(0) &= \bar{\rho}, \quad \delta(0) = \bar{\delta}, \quad N(0) = \bar{N}, \quad c(0) = \bar{c} \end{aligned} \right. \qquad (6.6) $$

This equation can be solved analytically. The solution is:

$$ x(\tau) = \begin{bmatrix} \bar{T}_U + \bar{T}_I - \frac{\bar{\lambda}}{\bar{\rho}} (\bar{\delta}e^{-\bar{\rho}\tau} + \bar{\rho} - \bar{\delta}) \\ \frac{\bar{T}_I}{\bar{\rho}} (\bar{\delta}e^{-\bar{\rho}\tau} + \bar{\rho} - \bar{\delta}) \\ \bar{V} \\ \bar{\lambda} \\ \bar{\rho} \\ \frac{\bar{\delta}\bar{\rho}}{\bar{N}e^{\bar{\rho}\tau}} \\ \bar{c} \end{bmatrix} \qquad (6.7) $$

All the above states ($x(\tau)$) are indistinguishable for any $\tau$ (in particular, they are indistinguishable from $[\bar{T}_U, \bar{T}_I, \bar{V}, \bar{\lambda}, \bar{\rho}, \bar{\delta}, \bar{N}, \bar{c}]^T$). In Section 6.6 we explicitly show that they produce the same outputs.
## 6.3 Characterization by an observable state

Note that this task is not requested by Algorithm 4, to obtain the identifiability of the parameters. In addition, this task cannot be performed automatically. On the other hand, we were able to obtain this description that can be useful for further investigations.

The original state has dimension 8 and the observability codistribution 7. This is the reason why we have a single symmetry ($8 = 7 + 1$). A possible observable state consists precisely of a set of generators of the observability codistribution. Algorithm 1 automatically returns a set of generators. We obtained:

$$h_1, h_2, \mathcal{L}_{g^0}h_1, \mathcal{L}_{g^0}h_2, \mathcal{L}_{\tilde{g}^0}\mathcal{L}_{g^0}h_1, \mathcal{L}_{\tilde{g}^1}\mathcal{L}_{g^0}h_1, \mathcal{L}_{\tilde{g}^0}^2\mathcal{L}_{g^0}h_1.$$

On the other hand, by setting the state entries equal to these functions, obtaining the expression of its dynamics in terms of the same state components is definitely prohibitive (although possible). To simplify this task, first of all, we try to determine 7 independent and simple functions, which are observable and that generate the observability codistribution. By using the expression of the above generators and the expression of the above symmetry, by using our inventiveness, we were able to select the following observable functions:

$$\begin{aligned}
- \Psi_1 &:= V \\
- \Psi_2 &:= T_U + T_I \\
- \Psi_3 &:= N\delta T_I \\
- \mu &:= \frac{\delta - \rho}{N\delta} \\
- \lambda & \\
- c & \\
- \rho &
\end{aligned}
\qquad (6.8)$$

(by a direct calculation it is possible to verify that their gradients generate $\mathcal{O}$). As their expression is much easier than the expression of the generators automatically selected by Algorithm 1, we have the possibility to describe our system by using them. After some trials, by using (6.1), we were able to obtain the following dynamics:

$$\left\{
\begin{aligned}
\dot{\Psi}_1 &= \Psi_3 - c\Psi_1 \\
\dot{\Psi}_2 &= \lambda - \rho\Psi_2 - \Psi_3 \mu \\
\dot{\Psi}_3 &= \tilde{\eta} \\
\dot{\mu} &= 0 \\
\dot{\lambda} &= 0 \\
\dot{c} &= 0 \\
\dot{\rho} &= 0 \\
y &= [\Psi_1, \Psi_2]
\end{aligned}
\right.
\qquad (6.9)$$

The correctness of the above expression can be easily verified by using in (6.9) the expressions of $\Psi_1$, $\Psi_2$, $\Psi_3$, and $\mu$ given in (6.8) and by using Equation (6.1). Note that, in order to achieve this description, we had to redefine the unknown input ($\tilde{\eta} := N\delta T_U \Psi_1 \eta - \delta \Psi_3$). By using the original unknown input $\eta$, we would obtain $\dot{\Psi}_3 = N\delta T_U \Psi_1 \eta - \delta \Psi_3$. This equation cannot be expressed only in terms of the components of the new observable state. It needs the product $N\delta T_U$ and $\delta$, which cannot be expressed in terms of the selected observable functions because both $N\delta T_U$ and $\delta$ are unobservable.

## 6.4 Identifiability analysis

We execute Algorithm 4. First, as the system is canonic with respect to its unknown inputs, we do not have symmetries of the unknown input due to the system non canonicity ($c\chi$). On the other hand, we have a single symmetry of the state ($\xi$), which is the single generator of $\mathcal{O}^{\perp}$ in (6.5), and, as a result, we may have a non vanishing symmetry $^u\chi$ of the unknown input. Let us compute it. In accordance with Equation (4.6), we need to compute $\xi_1^0$ and $\xi_1^1$. We have: $\tilde{h} = \tilde{h}_1 = \lambda - \rho T_U - \delta T_I$. From (6.3) we obtain:
$$ \mathcal{L}_{g^0} \tilde{h} = T_I \delta^2 - \rho(\lambda - T_U \rho), \quad \mathcal{L}_{g^1} \tilde{h} = -T_U V(\delta - \rho). $$

From Equation (4.7), by using the generator of $O^{\perp}$ in (6.5) for $\xi$, we obtain:

$$ \xi_1^0 = T_I \delta(\delta - \rho)^2, \quad \xi_1^1 = -V \delta(T_I + T_U)(\delta - \rho). $$

In addition, from Equation (6.4) and by knowing that $m = m_w = 1$ (canonic system), we have $m\nu_1^1 = \nu_1^1 = -\frac{1}{T_U V(\delta - \rho)}$.

By using (4.6) and the above expressions, we obtain:

$$ \begin{aligned} {}^u\chi = & -\nu_1^1 (\xi_1^0 + \xi_1^1 w_1) = \frac{1}{T_U V(\delta - \rho)} T_I \delta(\delta - \rho)^2 \\ & - \frac{1}{T_U V(\delta - \rho)} V \delta(T_I + T_U)(\delta - \rho) \eta. \end{aligned} $$

Hence:

$$ {}^u\chi = \frac{T_I \delta(\delta - \rho)}{T_U V} - \frac{\delta(T_I + T_U)}{T_U} \eta. \quad (6.10) $$

We have obtained a non trivial symmetry ($^u\chi$) of the unknown input. This means that the unknown input cannot be reconstructed (or, in other words, the time-varying parameter $\eta(t)$ is not identifiable). This result contradicts the result available in the state of the art (e.g., see Section 6.2 of [3] or Section 3.3 of [35]). In Section 6.6, we prove the validity of the above result by quantitatively providing more unknown inputs that are consistent with the same outputs.

## 6.5 Indistinguishable states and unknown inputs

In order to obtain the indistinguishable states and unknown inputs we need to solve the differential equation in (5.7), for the specific case. We first need to check if the condition in (5.6) is honoured. By using the expression of $\xi$ (which is the generator of $O^{\perp}$ in (6.5)) and the expression of $^u\chi$ in (6.10), it is possible to check that this condition is satisfied and the differential equation in (5.7) becomes:

$$ \left\{ \begin{aligned} \frac{dT'_I}{d\tau'} &= T'_I \delta' \\ \frac{dT'_I}{d\tau'} &= -T'_I \delta' \\ \frac{dV'}{d\tau} &= 0 \\ \frac{d\lambda}{d\tau} &= 0 \\ \frac{d\rho'}{d\tau} &= 0 \\ \frac{d\delta'}{d\tau} &= \delta'(\delta' - \rho') \\ \frac{dN'}{d\tau} &= N' \rho' \\ \frac{dc'}{d\tau} &= 0 \\ \frac{d\eta'}{d\tau} &= \frac{T'_I \delta'(\delta' - \rho')}{T'_U V'} - \frac{\delta'(T'_I + T'_U)}{T'_U} \eta' \\ T'_I(t, 0) &= T_U(t), \quad T'_I(t, 0) = T_I(t), \\ V'(t, 0) &= V(t), \\ \lambda'(0) = \lambda, \quad \rho'(0) &= \rho, \quad \delta'(0) = \delta, \quad N'(0) = N, \quad c'(0) = c \\ \eta'(t, 0) &= \eta(t) \end{aligned} \right. \quad (6.11) $$

where we only provide the explicit dependence on $t$ and $\tau$ at the initial conditions and we omit it for the differential equations, for the simplicity sake. These equations can be solved analytically. The solution for the first eight components (which are the components of the state) is the same solution of (6.6) given in (6.7), with the appropriate values for the initial conditions. In other words, $x'(t, \tau)$ is:
$$
\left\{
\begin{aligned}
T'_{U}(t, \tau) &= T_{U}(t) + T_{I}(t) - \frac{T_{I}(t)}{\rho} (\delta e^{-\rho\tau} + \rho - \delta) \\
T'_{I}(t, \tau) &= \frac{T_{I}(t)}{\rho} (\delta e^{-\rho\tau} + \rho - \delta) \\
V'(t, \tau) &= V(t) \\
\lambda'(\tau) &= \lambda \\
\rho'(\tau) &= \rho \\
\delta'(\tau) &= \delta\rho / ((\rho - \delta)e^{\rho\tau} + \delta) \\
N'(\tau) &= N e^{\rho\tau} \\
c'(\tau) &= c
\end{aligned}
\right.
\qquad (6.12)
$$

Regarding the unknown input, we obtain:

$$
\eta'(t, \tau) = \frac{\eta(t)T_U(t)V(t)\rho e^{\rho\tau} + [T_I(t)\delta^2 - T_I(t)\delta\rho - \eta(t)T_U(t)V(t)\delta][e^{\rho\tau} - 1]}{V(t)[T_I(t)\delta + T_U(t)\rho]e^{\rho\tau} - V(t)T_I(t)\delta} \quad (6.13)
$$

We have obtained the following fundamental result. Let us consider the ODE model in (6.1) and let us consider a given time interval $\mathcal{I} := [t_0, T]$. Let us suppose that the initial conditions and the values of the model parameters are:

• $T_U(t_0) = T_{U_0}$, $T_I(t_0) = T_{I_0}$, and $V(t_0) = V_0$.

• The values of the constant parameters are $\lambda, \rho, \delta, N$, and $c$.

• The time-varying parameter on $\mathcal{I}$ is $\eta(t)$.

Then, the following initial conditions and the values of the model parameters:

$$
\bullet \quad
\left[
\begin{aligned}
T'_{U}(t_0, \tau) &= T_{U_0} + T_{I_0} - \frac{T_{I_0}}{\rho} (\delta e^{-\rho\tau} + \rho - \delta) \\
T'_{I}(t_0, \tau) &= \frac{T_{I_0}}{\rho} (\delta e^{-\rho\tau} + \rho - \delta) \\
V'(t_0, \tau) &= V_0 \\
\lambda'(\tau) &= \lambda \\
\rho'(\tau) &= \rho \\
\delta'(\tau) &= \delta\rho / ((\rho - \delta)e^{\rho\tau} + \delta) \\
N'(\tau) &= N e^{\rho\tau} \\
c'(\tau) &= c
\end{aligned}
\right.
\qquad (6.14)
$$

• the time-varying parameter given in (6.13), with $T_U(t)$, $T_I(t)$, and $V(t)$ the values of $T_U$, $T_I$, and $V$ obtained by integrating (6.1) with initial conditions $T_{U_0}$, $T_{I_0}$, and $V_0$, and unknown input $\eta(t)$,

produce exactly the same outputs $y_1(t)$, $y_2(t)$, on the entire time interval $\mathcal{I}$. This holds for any value of the parameter $\tau$ that belongs to an interval $\mathcal{U} \subseteq \mathbb{R}$, which includes $\tau = 0$. In addition, the interval $\mathcal{U}$ is such that, for any $\tau \in \mathcal{U}$, the solution in (6.14) must exist. On the other hand, even if the above requirements on $\mathcal{U}$ ensure the validity of our result (i.e., that the solution in (6.13) and (6.14) produces the same outputs at any time) we have a further requirement, due to the physical meaning of our states and parameters. For any $\tau \in \mathcal{U}$ the quantities $\delta'(\tau)$, $N'(\tau)$, $\eta'(t, \tau)$, $T'_U(t, \tau)$, and $T'_I(t, \tau)$ must be positive.

In Section 6.6, we test the validity of the above fundamental result by using the same data set adopted in [35].

## 6.6 Comparison with the state of the art results

As we aforementioned, our results contradict the result available in the state of the art for exactly the same ODE model characterized by (6.1) (e.g., see Section 6.2 of [3]<fnref n="1" /> or Section 3.3 of [35]). Specifically, in contrast with the state of the art, we obtained that both the state and

<footnote n="1">¹The reader can find a detailed explanation of the serious error made by the authors of [3] in [5].</footnote>
the unknown input (or the time-varying parameter) $\eta(t)$, cannot be recovered from the outputs. In addition, we also obtained a one dimensional set of indistinguishable states, $x'_0(\tau)$, and a one parameter set of functions $\eta'(t, \tau)$ such that, by only using the outputs, it is not possible to distinguish among the elements of these sets. In other words, the outputs obtained by using the initial state $x'_0(\tau)$ and by using the unknown input function $\eta'(t, \tau)$, are independent of the parameter $\tau$, and this holds for any $\tau \in \mathcal{U}$.

In this section, we explicitly provide these sets by referring to the same data available in the literature (e.g., see the electronic supplementary material of [35]) and we show that any choice in these sets produces the same outputs.

The data set is characterized as follows:

* The time interval is $\mathcal{I} = [0, 201]$ (in days).

* $T_U(0) = 600$, $T_I(0) = 0$, $V(0) = 10^5$, $\lambda = 36$, $\rho = 0.108$, $\delta = 0.5$, $N = 10^3$, and $c = 3$.

* The time-varying parameter is:

$$ \eta(t) = k \left( 1 - 0.9 \cos \left( \frac{\pi}{1000} t \right) \right), \qquad (6.15) $$

with $k = 0.00009$.

<figure><img src="image_4.png" /><figcaption>Figure 6.1: Several indistinguishable profiles for the time-varying parameter $\eta'(t, \tau)$ for several values of the parameter $\tau$ ranging from -3 up to $\tau_* \approx 2.2532$. The dashed purple line is the profile given by (6.15), i.e., $\eta'(t, \tau)$ for $\tau = 0$.</figcaption></figure>

We provide indistinguishable initial states and indistinguishable unknown inputs for several values of the parameter $\tau \in \mathcal{U}$. First of all, let us characterize $\mathcal{U}$.

From the expression of $\delta'(\tau)$ in (6.14) and the values of our data set, we obtain that

$$ \lim_{\tau \to \tau_*^-} \delta'(\tau) = +\infty, \quad \lim_{\tau \to \tau_*^+} \delta'(\tau) = -\infty $$

with $\tau_* := \frac{\log(\frac{\delta}{\delta-\rho})}{\rho} \approx 2.2532$.
<figure><img src="image_5.png" /><figcaption>Figure 6.2: The profiles $T_U'(t, \tau)$ for several values of the parameter $\tau$ ranging from -3 up to $\tau_* \approx 2.2532$. The dashed purple line is the profile for $\tau = 0$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 6.3: The profiles $T_I'(t, \tau)$ for several values of the parameter $\tau$ ranging from -3 up to $\tau_* \approx 2.2532$. The dashed purple line is the profile for $\tau = 0$.</figcaption></figure>
In addition, we found that, for $\tau = -4$ the function $\eta'(t, \tau)$ in (6.13) takes negative values for $t \in [2, 6]$s. Finally, we verified that in the interval

$$ \mathcal{U}_0 := [-3, \tau_*) $$

all the quantities $\delta'(\tau)$, $N'(\tau)$, $\eta'(t, \tau)$, $T'_U(t, \tau)$, and $T'_I(t, \tau)$ take positive values (the last three on the entire time interval $\mathcal{I}$). The condition $\mathcal{U}_0 \subset \mathcal{U}$ is ensured and we are allowed to use any $\tau \in \mathcal{U}_0$. Hence, we consider the indistinguishable initial states and unknown inputs obtained by varying $\tau$ in this interval $\mathcal{U}_0$.

Figure 6.1 displays the profiles of $\eta'(t, \tau)$. The profile for $\tau = 0$ (purple dashed line) is precisely the one set in (6.15). By varying $\tau$ we obtain a significant change of the profile meaning that this parameter is "strongly" unidentifiable.

Figures 6.2 and 6.3 display the profiles of $T'_U(t, \tau)$ and $T'_I(t, \tau)$, respectively. The initial conditions of $T_U$ and $T_I$ are obtained by these profiles at $t = 0$. By chance, in this case, these initial conditions are independent of $\tau$. This because the part that depends on $\tau$ is linear in $T_{I_0}$ (see the first two equations in (6.14)) and we set $T_{I_0} = 0$. Note that, by changing the values of $T_{I_0}$, the initial conditions significantly depend on $\tau$ meaning that also the initial state is "strongly" unobservable. In addition, by varying $\tau$, the two constant parameters that are not identifiable (i.e., $\delta$ and $N$) significantly change. For instance, in accordance with (6.14), for $\tau = -3$ and $\tau = \tau_*$ we obtain the following changes:

$$ \begin{aligned} \delta &= \delta'(0) = 0.5 \quad \rightarrow \delta'(-3) \approx 0.25 \quad \rightarrow \delta'(\tau \rightarrow \tau_*) = +\infty \\ N &= N'(0) = 10^3 \quad \rightarrow N'(-3) \approx 723 \quad \rightarrow N'(\tau_*) \approx 1276. \end{aligned} $$

The variation of $\delta'$ is even divergent.

We verified that the two outputs, $y_1(t) = V$ and $y_2(t) = T_U + T_I$, are independent of $\tau$, at any $t$. This is obtained by proceeding as follows.

We considered several values of $\tau$ in $\mathcal{U}_0$. For each $\tau$, we integrated the differential equation in (6.1) with the initial conditions given in (6.14) for that $\tau$ and by setting the unknown time-varying parameter as in (6.13) for the same value of $\tau$ (that is one of the profiles plotted in Figure 6.1). Both the initial conditions in (6.14) and the unknown input significantly depend on $\tau$. However, the two outputs (which are displayed in Figure 6.4) are independent, in accordance with our results.

We also check the validity of the same results by considering other settings (e.g., by setting $T_{I_0} \neq 0$).

## 6.7 Minimal external information

Our analysis showed that only three parameters are identifiable. They are $\lambda$, $\rho$, and $c$. The remaining constant parameters, $\delta$, and $N$, and the time varying parameter, $\eta(t)$, cannot be identified. It is possible to prove that, by also including the output $y_3 = T_U$ (or $y_3 = T_I$), the state becomes observable and all the parameters identifiable. This can be proved by repeating the analysis for the new system characterized by this further output. On the other hand, from our analysis, it is immediate to obtain a more interesting result.

The system is characterized by a single symmetry. As a result, we can apply the general procedure described in Section 5.2 and we can conclude that it suffices to have $T_U$ (or $T_I$) at a single time $t^* \in \mathcal{I}$ (and not necessarily on the entire time interval) to make possible the identification of all the parameters (and the observability of the state). Equation (5.8) becomes, in this case, (6.12) together with (6.13). In order to determine the true state and the true parameters, we must determine the value of $\tau$. This can be determined by knowing $T_U$ (or equivalently $T_I = y_2 - T_U$) at a single time. First of all, Equation (6.12) tells us that we have the true $V$, $\lambda$, $\rho$, and $c$. Now let us suppose that we know $T_I$ at a single time. From the second equation in (6.12) we obtain the quantity $\delta e^{-\rho\tau} - \delta = \frac{T'_I(t^*, \tau) \rho}{T_I(t^*)} - \rho$, which is an equation in the two unknowns $\delta$ and $\tau$. On the other hand, the sixth equation in (6.12) is a further equation in the same two unknowns. From them, we obtain both $\delta$ and $\tau$. Finally, once $\tau$ has been determined, we determine the remaining quantities at any $t \in \mathcal{I}$.
<figure><img src="image_7.png" /><figcaption>Figure 6.4: The two outputs $y_1(t) = V$ (solid blue line) and $y_2(t) = T_U + T_I$ (dashed red line). The profiles are independent of the parameter $\tau$, as expected.</figcaption></figure>
# Chapter 7

## Compartmental SEIAR model of the epidemic dynamics of Covid-19

We consider a simple model that generalizes the SEIR model commonly used for virus disease. This more general ODE model includes a further compartment, A(t), that is the asymptomatic infected, [44]. The other compartments are the same of a SEIR model, i.e.: S, Susceptible; E, Exposed; I, Infected (symptomatic infected); R, Removed (i.e., healed or dead).

The differential equations of this model describe the movement of individuals of the population between the five aforementioned classes, i.e.: S(t), E(t), I(t), A(t) and R(t).

The dynamics and the outputs are given by the following set of equations (e.g., see [44]):

$$ 
\begin{cases}
\dot{S} = -\beta S(I + A) \\
\dot{E} = \beta S(I + A) - \gamma E \\
\dot{I} = \gamma p E - \mu_1 I \\
\dot{A} = \gamma (1 - p) E - \mu_2 A \\
\dot{R} = \mu_1 I + \mu_2 A \\
y = [I, A, S + E + R],
\end{cases}
\quad (7.1) $$

This model includes the following constant parameters, which are assumed to be unknown:

* $\mu_1$ and $\mu_2$ that are the remove rates from the infected symptomatic individuals ($I(t)$) and the infected asymptomatic individuals ($A(t)$), respectively.

* $\gamma$ that is the rate at which the exposed individuals ($E(t)$) are infected.

* $p$ that is the probability that infected individuals are symptomatic (and $1-p$ is the probability that infected individuals are asymptomatic).

The model also includes the following unknown time-varying parameter:

* $\beta = \beta(t)$ that is the probability of disease transmission in a single contact times the average number of contacts per person, due to contacts between an individual of the class S and an individual that belongs to one of the two classes I and A. This parameter is assumed to be time-varying as it can be modified by government measures (e.g., using masks, closing schools, remote working).

We assume that we can measure the infected symptomatic individuals ($I(t)$), and the infected asymptomatic individuals ($A(t)$). Hence, our model is characterized by the first two outputs in (7.1). Finally, as we know the total population, $S(t)+E(t)+I(t)+A(t)+R(t)$, by removing from this the above two measurements, we obtain the third output in (7.1), i.e., $S(t)+E(t)+R(t)$.

To proceed, we need, first of all, to introduce a state that includes both the time-varying quantities (i.e., S, E, I, A, R) and the constant parameters (i.e., $\mu_1$, $\mu_2$, $\gamma$, $p$). We set:
$$x = [S, E, I, A, R, \mu_1, \mu_2, \gamma, p]^T \qquad (7.2)$$

The system is directly a special case of (2.2). In particular, $m_u = 0, m_w = 1, p = 3$, $h_1(x) = I, h_2(x) = A, h_3(x) = S + E + R$,

$$g^0 = \begin{bmatrix} 0 \\ -\gamma E \\ \gamma p E - \mu_1 I \\ \gamma(1-p)E - \mu_2 A \\ \mu_1 I + \mu_2 A \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}, \quad g^1 = \begin{bmatrix} -S(I+A) \\ S(I+A) \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} \qquad (7.3)$$

## 7.1 Observability analysis

We run Algorithm 1 to obtain the observability codistribution. Line 2 provides $\Omega = \text{span}\{\nabla h_1, \nabla h_2, \nabla h_3\}$, and $\text{deg}_w(\Omega) = 0$ ($< m_w$), as $\mathcal{L}_{g^1}h_1 = \mathcal{L}_{g^1}h_2 = \mathcal{L}_{g^1}h_3 = 0$ (and the rank of the unknown input reconstructability matrix from $h_1, h_2, h_3$ vanishes). The system is not in canonical form with respect to its unknown input. We enter the loop at lines 9-14. Line 11 provides no functions ($m = 0$). Line 12 provides: $^m\mu_0^0 = ^m\nu_0^0 = 1$ obtained from (3.6) (these tensors have a single entry as the indices only take the value 0). In addition, from (3.7) we obtain: $^m\hat{g}^0 = ^m\nu_0^0 g^0 = g^0$.

Line 13 computes

$$\Omega = \langle \left\{ \begin{matrix} g^0 \\ g^1 \end{matrix} \right\} | \text{span}\{\nabla h_1, \nabla h_2, \nabla h_3\} \rangle,$$

as $^m\hat{g}^0 = g^0$. By a direct computation, we obtain:

$$\Omega = \text{span}\{\nabla h_1, \nabla h_2, \nabla h_3, \nabla \mathcal{L}_{g^0}h_1, \nabla \mathcal{L}_{g^0}h_2\} = \\ \text{span}\{\nabla h_1, \nabla h_2, \nabla h_3, \nabla h_4, \nabla h_5\}$$

where we adopted the following notation:

$$h_4 := \mathcal{L}_{g^0}h_1, \quad h_5 := \mathcal{L}_{g^0}h_2$$

Now, $\text{deg}_w(\Omega) = 1 = m_w$. Hence, the second condition at Line 9 is not met and the loop is interrupted. In addition, as the condition at Line 15 is met, also the main loop is interrupted. The execution of Algorithm 1 continues with Line 45. The operation at this line ($[\tilde{h}_1, \dots, \tilde{h}_{m_w}] = S(\Sigma, \Omega)$) sets

$$\tilde{h}_1 = h_4 = \mathcal{L}_{g^0}h_1 = \gamma p E - \mu_1 I$$

(we would obtain the same final result by setting $\tilde{h}_1 = h_5 = \mathcal{L}_{g^0}h_2$). Then, Line 46 provides:

$$\mu = \begin{bmatrix} 1 & 0 \\ \mu_1(I\mu_1 - E\gamma p) - E\gamma^2 p & S\gamma p(A+I) \end{bmatrix}, \\ \nu = \begin{bmatrix} 1 & 0 \\ \frac{Ep\gamma^2 + Ep\gamma\mu_1 - I\mu_1^2}{S\gamma p(A+I)} & \frac{1}{S\gamma p(A+I)} \end{bmatrix},$$

obtained from (3.3). In addition, from (3.4) we obtain:
$$ \hat{g}^0 = \begin{bmatrix} -\frac{Ep\gamma^2+Ep\gamma\mu_1-I\mu_1^2}{\gamma p} \\ -\frac{\mu_1(I\mu_1-E\gamma p)}{\gamma p} \\ E\gamma p - I\mu_1 \\ -A\mu_2 - E\gamma(p-1) \\ A\mu_2 + I\mu_1 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}, \quad \hat{g}^1 = \begin{bmatrix} -\frac{1}{\gamma p} \\ \frac{1}{\gamma p} \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$

Finally, from (3.5) we obtain that, when $m_u = 0$, the codistribution $\tilde{\mathcal{O}}$ vanishes. The final step of Algorithm 1 is the execution of Line 47, namely:

$$ \mathcal{O} = \langle \hat{g}^0, \hat{g}^1 | \operatorname{span}\{\nabla h_1, \nabla h_2, \nabla h_3, \nabla h_4, \nabla h_5\} \rangle. $$

The algorithm in (2.8) converges at the second step ($\Omega_2 = \Omega_1$), and

$$ \mathcal{O} = \operatorname{span}\{\nabla h_1, \nabla h_2, \nabla h_3, \nabla h_4, \nabla h_5, \nabla \mathcal{L}_{\hat{g}^0} h_5, \nabla \mathcal{L}_{\hat{g}^1} h_5\}. $$

Its dimension is 7, which is smaller than the dimension of the state in (7.2). The state is unobservable.

## 7.2 Symmetries and indistinguishable states

Once we have $\mathcal{O}$, the symmetries are immediately obtained by computing the orthogonal distribution. We have:

$$ \mathcal{O}^{\perp} = \operatorname{span} \left\{ \begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \\ -1 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}, \begin{bmatrix} E \\ -E \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \gamma \\ 0 \end{bmatrix} \right\}. \qquad (7.4) $$

As for the case study, starting from the generators of $\mathcal{O}^{\perp}$ we compute the indistinguishable states (although the determination of the indistinguishable states is not requested by Algorithm 4). We need to solve the differential equation in (4.1). Given a state $[\bar{S}, \bar{E}, \bar{I}, \bar{A}, \bar{R}, \bar{\mu}_1, \bar{\mu}_2, \bar{\gamma}, \bar{p}]^T$, the solution of this differential equation provides states which are indistinguishable. For the specific case, we have two distinct independent symmetries that are the two generators of $\mathcal{O}^{\perp}$ in (7.4). Let us discuss the two cases separately.

### 7.2.1 First symmetry

Equation (4.1) becomes:
$$
\left\{
\begin{array}{l}
\displaystyle \frac{dS}{d\tau} = 1 \\
\displaystyle \frac{dE}{d\tau} = 0 \\
\displaystyle \frac{dI}{d\tau} = 0 \\
\displaystyle \frac{dA}{d\tau} = 0 \\
\displaystyle \frac{dR}{d\tau} = 0 \\
\displaystyle \frac{d\mu_1}{d\tau} = 0 \\
\displaystyle \frac{d\mu_2}{d\tau} = 0 \\
\displaystyle \frac{d\gamma}{d\tau} = 0 \\
\displaystyle \frac{dp}{d\tau} = 0 \\
S(0) = \bar{S}, E(0) = \bar{E}, I(0) = \bar{I}, A(0) = \bar{A}, R(0) = \bar{R}, \\
\mu_1(0) = \bar{\mu}_1, \mu_2(0) = \bar{\mu}_2, \gamma(0) = \bar{\gamma}, p(0) = \bar{p}
\end{array}
\right.
\qquad (7.5)
$$

This equation can be trivially solved analytically. The solution is:

$$
x(\tau) = \begin{bmatrix} \frac{\bar{S} + \tau}{\bar{E}} \\ \frac{\bar{I}}{\bar{A}} \\ \bar{R} - \tau \\ \frac{\bar{\mu}_1}{\bar{\mu}_2} \\ \frac{\bar{\gamma}}{\bar{p}} \end{bmatrix} \qquad (7.6)
$$

All the above states ($x(\tau)$) are indistinguishable for any $\tau$ (in particular, they are indistinguishable from $[\bar{S}, \bar{E}, \bar{I}, \bar{A}, \bar{R}, \bar{\mu}_1, \bar{\mu}_2, \bar{\gamma}, \bar{p}]^T$).

7.2.2 Second symmetry

Equation (4.1) becomes:

$$
\left\{
\begin{array}{l}
\displaystyle \frac{dS}{d\tau} = E \\
\displaystyle \frac{dE}{d\tau} = -E \\
\displaystyle \frac{dI}{d\tau} = 0 \\
\displaystyle \frac{dA}{d\tau} = 0 \\
\displaystyle \frac{dR}{d\tau} = 0 \\
\displaystyle \frac{d\mu_1}{d\tau} = 0 \\
\displaystyle \frac{d\mu_2}{d\tau} = 0 \\
\displaystyle \frac{d\gamma}{d\tau} = \gamma \\
\displaystyle \frac{dp}{d\tau} = 0 \\
S(0) = \bar{S}, E(0) = \bar{E}, I(0) = \bar{I}, A(0) = \bar{A}, R(0) = \bar{R}, \\
\mu_1(0) = \bar{\mu}_1, \mu_2(0) = \bar{\mu}_2, \gamma(0) = \bar{\gamma}, p(0) = \bar{p}
\end{array}
\right.
\qquad (7.7)
$$

This equation can be solved analytically. The solution is:

$$
x(\tau) = \begin{bmatrix} \frac{-\bar{E} e^{-\tau} + \bar{S} + \bar{E}}{\bar{E} e^{-\tau}} \\ \frac{\bar{I}}{\bar{A}} \\ \bar{R} \\ \frac{\bar{\mu}_1}{\bar{\mu}_2} \\ \frac{\bar{\gamma} e^{\tau}}{\bar{p}} \end{bmatrix} \qquad (7.8)
$$

All the above states ($x(\tau)$) are indistinguishable for any $\tau$ (in particular, they are indistinguishable from $[\bar{S}, \bar{E}, \bar{I}, \bar{A}, \bar{R}, \bar{\mu}_1, \bar{\mu}_2, \bar{\gamma}, \bar{p}]^T$).
## 7.3 Characterization by an observable state

Note that this task is not requested by Algorithm 4, to obtain the identifiability of the parameters. In addition, this task cannot be performed automatically. On the other hand, we were able to obtain this description that can be useful for further investigations.

The original state has dimension 9 and the observability codistribution 7. This is the reason why we have two independent symmetries ($9 = 7 + 2$). A possible observable state consists precisely of a set of generators of the observability codistribution. Algorithm 1 automatically returns a set of generators. We obtained: $h_1, h_2, h_3, h_4, h_5, \mathcal{L}_{\hat{g}^0}h_5, \mathcal{L}_{\hat{g}^1}h_5$. On the other hand, by setting the state entries equal to these functions, obtaining the expression of its dynamics in terms of the same state components is definitely prohibitive (although possible). To simplify this task, first of all, we try to determine 7 independent and simple functions, which are observable and that generate the observability codistribution. By using the expression of the generators of $\mathcal{O}$ and the expression of the generators of $\mathcal{O}^{\perp}$, by using our inventiveness we were able to select the following observable functions:

$$ 
\begin{array}{l}
- \Psi_1 := S + E + R \\
- \Psi_2 := \gamma E \\
- I \\
- A \\
- \mu_1 \\
- \mu_2 \\
- p
\end{array} 
\qquad (7.9) $$

(by a direct calculation it is possible to verify that their gradients generate $\mathcal{O}$). As their expression is much easier than the expression of the generators automatically selected by Algorithm 1, we have the possibility to describe our system by using them. After some trials, by using (7.1), we were able to obtain the following dynamics:

$$ 
\left\{ 
\begin{array}{l}
\dot{\Psi}_1 = -\Psi_2 + \mu_1 I + \mu_2 A \\
\dot{\Psi}_2 = \tilde{\beta} \\
\dot{I} = p\Psi_2 - \mu_1 I \\
\dot{A} = (1-p)\Psi_2 - \mu_2 A \\
y = [I, A, \Psi_1],
\end{array} 
\right. 
\qquad (7.10) $$

Note that, in order to achieve this description, we had to redefine the unknown input. By using the original unknown input $\beta$, we would obtain $\tilde{\Psi}_2 = \gamma \tilde{E} = \beta S \gamma (I + A) - \gamma \Psi_2$. This equation cannot be expressed only in terms of the components of the new observable state. It needs the product $S\gamma$ and $\gamma$, which cannot be expressed in terms of the selected observable functions because they are unobservable.

## 7.4 Identifiability analysis

We execute Algorithm 4. First, as the system is canonic with respect to its unknown inputs, we do not have symmetries of the unknown input due to the system non canonicity ($c\chi$). On the other hand, we have two independent symmetries of the state, which are the generators of $\mathcal{O}^{\perp}$ in (7.4), and, as a result, we may have non vanishing symmetries $^u\chi$ of the unknown input. Let us compute them. We use Equation (4.6) twice, by setting $\xi$ equal to the two generators of $\mathcal{O}^{\perp}$ in (7.4), respectively. In both cases, we have:

$$ \tilde{h} := \tilde{h}_1 = \gamma p E - \mu_1 I. $$

and from (7.3), we obtain:

$$ \mathcal{L}_{g^0} \tilde{h} = \mu_1(I\mu_1 - E\gamma p) - E\gamma^2 p, \quad \mathcal{L}_{g^1} \tilde{h} = S\gamma p(A + I). \qquad (7.11) $$
### 7.4.1 First symmetry

We set

$$ \xi = \xi^1 = \begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \\ -1 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}. $$

In accordance with Equation (4.6), we need to compute $\xi_1^0$ and $\xi_1^1$. From Equation (4.7), by using the above $\xi$ and (7.11), we obtain:

$$ \xi_1^0 = 0, \quad \xi_1^1 = \gamma p(A + I). $$

and, by substituting in (4.6), we obtain:

$$ ^u\chi = - \stackrel{m}{\nu}_1^1(\xi_1^0 + \xi_1^1 w_1) = -\frac{\xi_1^1}{S\gamma p(A+I)}w_1 = -\frac{1}{S}\beta $$

because $\stackrel{m}{\nu}_1^1 = \nu_1^1 = \frac{1}{S\gamma p(A+I)}$. Hence, we have the following first symmetry of the unknown input:

$$ ^u\chi^1 = -\frac{1}{S}\beta \qquad (7.12) $$

### 7.4.2 Second symmetry

We set

$$ \xi = \xi^2 = \begin{bmatrix} E \\ -E \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \gamma \\ 0 \end{bmatrix}. $$

In accordance with Equation (4.6), we need to compute $\xi_1^0$ and $\xi_1^1$. From Equation (4.7), by using the above $\xi$ and (7.11), we obtain:

$$ \xi_1^0 = -\gamma^2 p E, \quad \xi_1^1 = \gamma p(A + I)(E + S). $$

and, by substituting in (4.6), we obtain:

$$ ^u\chi = - \stackrel{m}{\nu}_1^1(\xi_1^0 + \xi_1^1 w_1) = -\frac{\xi_1^0}{S\gamma p(A+I)} - \frac{\xi_1^1}{S\gamma p(A+I)}w_1 = \frac{\gamma E}{S(A+I)} - \frac{E+S}{S}\beta $$

because $\stackrel{m}{\nu}_1^1 = \nu_1^1 = \frac{1}{S\gamma p(A+I)}$. Hence, we have the following second symmetry of the unknown input:

$$ ^u\chi^2 = \frac{\gamma E}{S(A+I)} - \frac{E+S}{S}\beta \qquad (7.13) $$

We conclude by remarking that the unknown input is characterized by two symmetries. Hence, $\beta(t)$ is unidentifiable.
## 7.5 Indistinguishable states and unknown inputs

We have two generators of $O^{\perp}$ and we discuss the two cases separately.

### 7.5.1 First set of solutions

By using the first generator of $O^{\perp}$ in (7.4)) and the expression of $^u\chi$ in (7.12), it is possible to check that the condition in (5.6) is satisfied and the differential equation in (5.7) becomes:

$$ 
\left\{ 
\begin{aligned}
\frac{dS'}{d\tau'} &= 1 \\
\frac{dE'}{d\tau'} &= 0 \\
\frac{dI'}{d\tau'} &= 0 \\
\frac{dA'}{d\tau'} &= 0 \\
\frac{dR'}{d\tau'} &= -1 \\
\frac{d\mu_1'}{d\tau'} &= 0 \\
\frac{d\mu_2'}{d\tau'} &= 0 \\
\frac{d\gamma'}{d\tau'} &= 0 \\
\frac{dp'}{d\tau'} &= 0 \\
\frac{d\beta'}{d\tau'} &= -\frac{\beta'}{S'} \\
S'(t, 0) &= S(t), \quad E'(t, 0) = E(t), \quad I'(t, 0) = I(t), \\
A'(t, 0) &= A(t), \quad R'(t, 0) = R(t), \quad \mu_1'(0) = \mu_1, \\
\mu_2'(0) &= \mu_2, \quad \gamma'(0) = \gamma, \quad p'(0) = p \\
\beta'(t, 0) &= \beta(t)
\end{aligned}
\right. 
\qquad (7.14) 
$$

where we only provide the explicit dependence on $t$ and $\tau$ at the initial conditions and we omit it for the differential equations, for the simplicity sake. These equations can be solved analytically. The solution for the first nine components (which are the components of the state) is the same solution of (7.5) given in (7.6), with the appropriate values for the initial conditions. In other words, $x'(t, \tau)$ is:

$$ 
\left\{ 
\begin{aligned}
S'(t, \tau) &= S(t) + \tau \\
E'(t, \tau) &= E(t) \\
I'(t, \tau) &= I(t) \\
A'(t, \tau) &= A(t) \\
R'(t, \tau) &= R(t) - \tau \\
\mu_1'(\tau) &= \mu_1 \\
\mu_2'(\tau) &= \mu_2 \\
\gamma'(\tau) &= \gamma \\
p'(\tau) &= p
\end{aligned}
\right. 
\qquad (7.15) 
$$

Regarding the unknown input, we obtain:

$$ \beta'(t, \tau) = \frac{\beta(t)S(t)}{S(t) + \tau}. \qquad (7.16) $$

We have obtained the following fundamental result. Let us consider the ODE model in (7.1) and let us consider a given time interval $[t_0, T]$. Let us suppose that the initial conditions and the values of the model parameters are:

* $S(t_0) = S_0$, $E(t_0) = E_0$, $I(t_0) = I_0$, $A(t_0) = A_0$, and $R(t_0) = R_0$.

* The values of the constant parameters are $\mu_1$, $\mu_2$, $\gamma$, and $p$.

* The time-varying parameter on the considered time interval $[t_0, T]$ is $\beta(t)$.

Then, the following initial conditions and the values of the model parameters:
<figure><img src="image_8.png" /></figure>

$$
\left\{
\begin{array}{l}
S'(t_0, \tau) = S_0 + \tau \\
E'(t_0, \tau) = E_0 \\
I'(t_0, \tau) = I_0 \\
A'(t_0, \tau) = A_0 \\
R'(t_0, \tau) = R_0 - \tau \\
\mu_1'(\tau) = \mu_1 \\
\mu_2'(\tau) = \mu_2 \\
\gamma'(\tau) = \gamma \\
p'(\tau) = p
\end{array}
\right.
\qquad (7.17)
$$

• the time-varying parameter given in (7.16), with $S(t)$, $E(t)$, $I(t)$, $A(t)$, and $R(t)$ the values of $S$, $E$, $I$, $A$, and $R$ obtained by integrating (7.1) with initial conditions $S_0$, $E_0$, $I_0$, $A_0$, and unknown input $\beta(t)$,

produce exactly the same outputs $y_1(t)$, $y_2(t)$, $y_3(t)$, on the entire time interval $[t_0, T]$. This
holds for any value of the parameter $\tau$ that belongs to an interval $\mathcal{U} \subseteq \mathbb{R}$, which includes $\tau = 0$,
and all the values of $\tau$ for which $\beta'(t, \tau)$, $S'(t, \tau)$, and $R'(t, \tau)$ are positive.

7.5.2 Second set of solutions

By using the second generator of $O^{\perp}$ in (7.4)) and the expression of $^u\chi$ in (7.13), it is possible
to check that the condition in (5.6) is satisfied and the differential equation in (5.7) becomes:

$$
\left\{
\begin{array}{l}
\displaystyle \frac{dS'}{d\tau'} = E' \\[1.5ex]
\displaystyle \frac{dE'}{d\tau'} = -E' \\[1.5ex]
\displaystyle \frac{dI'}{d\tau'} = 0 \\[1.5ex]
\displaystyle \frac{dA'}{d\tau'} = 0 \\[1.5ex]
\displaystyle \frac{dR'}{d\tau'} = 0 \\[1.5ex]
\displaystyle \frac{d\mu_1'}{d\tau'} = 0 \\[1.5ex]
\displaystyle \frac{d\mu_2'}{d\tau'} = 0 \\[1.5ex]
\displaystyle \frac{d\gamma'}{d\tau'} = \gamma' \\[1.5ex]
\displaystyle \frac{dp'}{d\tau'} = 0 \\[1.5ex]
\displaystyle \frac{d\beta'}{d\tau'} = \frac{\gamma' E'}{S'(A'+I')} - \frac{(E'+S')\beta'}{S'} \\[1.5ex]
S'(t, 0) = S(t), E'(t, 0) = E(t), I'(t, 0) = I(t), \\
A'(t, 0) = A(t), R'(t, 0) = R(t), \mu_1'(0) = \mu_1, \\
\mu_2'(0) = \mu_2, \gamma'(0) = \gamma, p'(0) = p \\[1.5ex]
\beta'(t, 0) = \beta(t)
\end{array}
\right.
\qquad (7.18)
$$

where we only provide the explicit dependence on *t* and *τ* at the initial conditions and we
omit it for the differential equations, for the simplicity sake. These equations can be solved
analytically. The solution for the first nine components (which are the components of the state)
is the same solution of (7.7) given in (7.8), with the appropriate values for the initial conditions.
In other words, *x*'(*t*, *τ*) is:

$$
\left\{
\begin{array}{l}
S'(t, \tau) = S(t) + E(t) (1 - e^{-\tau}) \\
E'(t, \tau) = E(t) e^{-\tau} \\
I'(t, \tau) = I(t) \\
A'(t, \tau) = A(t) \\
R'(t, \tau) = R(t) \\
\mu_1'(\tau) = \mu_1 \\
\mu_2'(\tau) = \mu_2 \\
\gamma'(\tau) = \gamma e^{\tau} \\
p'(\tau) = p
\end{array}
\right.
\qquad (7.19)
$$

Regarding the unknown input, we obtain:

$$
\beta'(t, \tau) = \frac{\mu(t)(1 - e^{\tau}) - S(t)\beta(t)}{E(t) - (E(t) + S(t))e^{\tau}}, \quad \mu := \frac{\gamma E(t)}{A(t) + I(t)}. \qquad (7.20)
$$
We have obtained the following fundamental result. Let us consider the ODE model in (7.1) and let us consider a given time interval $[t_0, T]$. Let us suppose that the initial conditions and the values of the model parameters are:

* $S(t_0) = S_0, E(t_0) = E_0, I(t_0) = I_0, A(t_0) = A_0$, and $R(t_0) = R_0$.

* The values of the constant parameters are $\mu_1, \mu_2, \gamma$, and $p$.

* The time-varying parameter on the considered time interval $[t_0, T]$ is $\beta(t)$.

Then, the following initial conditions and the values of the model parameters:

$$ \bullet \quad \left[ \begin{array}{lcl} S'(t_0, \tau) & = & S_0 + E_0 (1 - e^{-\tau}) \\ E'(t_0, \tau) & = & E_0 e^{-\tau} \\ I'(t_0, \tau) & = & I_0 \\ A'(t_0, \tau) & = & A_0 \\ R'(t_0, \tau) & = & R_0 \\ \mu_1'(\tau) & = & \mu_1 \\ \mu_2'(\tau) & = & \mu_2 \\ \gamma'(\tau) & = & \gamma e^{\tau} \\ p'(\tau) & = & p \end{array} \right. \qquad (7.21) $$

* the time-varying parameter given in (7.20), with $S(t), E(t), I(t), A(t)$, and $R(t)$ the values of $S, E, I, A$, and $R$ obtained by integrating (7.1) with initial conditions $S_0, E_0, I_0, A_0$, and unknown input $\beta(t)$,

produce exactly the same outputs $y_1(t), y_2(t), y_3(t)$, on the entire time interval $[t_0, T]$. This holds for any value of the parameter $\tau$ that belongs to an interval $\mathcal{U} \subseteq \mathbb{R}$, which includes $\tau = 0$, and all the values of $\tau$ for which $\gamma'(\tau), \beta'(t, \tau), S'(t, \tau)$, and $E'(t, \tau)$ are positive.

## 7.6 Numerical results

In this section, we provide some numerical results. We consider a given data set for our ODE model in (7.1) and we compute the outputs. Then, by using our results in Section 7.5, we produce other initial states and model parameters that produce exactly the same outputs. This unequivocally prove the state unobservability and the parameters unidentifiability.

The considered data set is characterized as follows. For the ground truth of the model parameters, we use the values available in [44]. We set:

* The time interval is $\mathcal{I} = [0, 200]$ (in days).

* $S(0) = 1 - 10^{-10}, E(0) = 0, I(0) = 10^{-10}, A(0) = 0, R(0) = 0, \mu_1 = \frac{1}{3}, \mu_2 = \frac{1}{10}, \gamma = \frac{1}{4}$, and $p = 0.14$.

* We consider separately two distinct profiles of the time-varying parameter, which are:

1. $\beta(t) = 1$, i.e., independent of time.

2.

$$ \beta(t) = \cos\left(\frac{\pi}{400}t\right). \qquad (7.22) $$

We provide indistinguishable initial states and indistinguishable unknown inputs. We have two distinct sets of solutions, which are the ones provided in Sections 7.5.1, and 7.5.2, respectively. However, we only provide the results obtained by using the second set. This because, by using the first set, independently of the choice of $\tau \in \mathbb{R}$, the solution in (7.15) provides a negative first component ($S(t) + \tau$) and/or a negative fifth component ($R(t) - \tau$). Indeed, with our data set, by integrating (7.1) in the time interval $\mathcal{I}$, the solutions $S(t)$ and $R(t)$ take all the values in the interval $[0, 1]$. As a result, there exists $t \in \mathcal{I}$ such that, for any $\tau \in \mathbb{R}$, at least one of the two components $S(t) + \tau$ and $R(t) - \tau$ is negative. In few words, for the first set of solutions given in
Section 7.5.1, the interval $\mathcal{U}$ that characterizes the parameter $\tau$ is empty. Note that this holds for the specific data set.

Let us consider the second set of solutions given in Section 7.5.2. We will discuss the two cases of constant $\beta$ ($\beta = 1$) and $\beta(t)$ given in (7.22), separately.

### 7.6.1 $\beta = 1$

First of all, let us characterize $\mathcal{U}$. By performing several runs, we obtain that all the states and parameters in (7.19) and $\beta'(t, \tau)$ in (7.20) take positive values, for any $t \in \mathcal{I}$, when $\tau \in [-0.05, \infty)$. Hence, we set:

$$ \mathcal{U}_0 := [-0.05, 10]. $$

The condition $\mathcal{U}_0 \subset \mathcal{U}$ is ensured and we are allowed to use any $\tau \in \mathcal{U}_0$. We also find that the plots obtained with $\tau > 10$ cannot be distinguished one each other, showing that we achieved a convergence. We consider the indistinguishable initial states and unknown inputs obtained by varying $\tau$ in $\mathcal{U}_0$.

<figure><img src="image_9.png" /><figcaption>Figure 7.1: Several indistinguishable profiles for the time-varying parameter $\beta'(t, \tau)$ for several values of the parameter $\tau$ ranging from -0.05 up to 10. The dashed purple line is the profile given by $\beta = 1$, i.e., $\beta'(t, \tau)$ for $\tau = 0$.</figcaption></figure>

Figure 7.1 displays the profiles of $\beta'(t, \tau)$. The profile for $\tau = 0$ (purple dashed line) is precisely the constant value $\beta(t) = 1$. By varying $\tau$ we obtain a significant change of the profile meaning that this parameter is "strongly" unidentifiable.

Figures 7.2 and 7.3 display the profiles of $S'(t, \tau)$ and $E'(t, \tau)$, respectively. The initial conditions of $S$ and $E$ are obtained by these profiles at $t = 0$. In addition, by varying $\tau$, the constant parameter that is not identifiable (i.e., $\gamma$) significantly changes. For instance, in accordance with (7.21), for $\tau = -0.05$ and $\tau = 10$ we obtain the following changes:

$$ \gamma = \gamma'(0) = 0.25 \rightarrow \gamma'(-0.05) \cong 0.2378 \rightarrow \gamma'(10) \cong 5506 $$
<figure><img src="image_10.png" /><figcaption>Figure 7.2: The profiles $S'(t, \tau)$ for several values of the parameter $\tau$ ranging from −0.05 up to 10. The dashed purple line is the profile for $\tau = 0$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 7.3: The profiles $E'(t, \tau)$ for several values of the parameter $\tau$ ranging from −0.05 up to 10. The dashed purple line is the profile for $\tau = 0$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 7.4: The three outputs $y_1(t) = I$ (solid blue line), $y_2(t) = A$ (dashed red line), and $y_3(t) = S + E + R$ (dotted yellow line). The three profiles are independent of the parameter $\tau$, as expected.</figcaption></figure>
We verified that the three outputs, $y_1(t) = I$, $y_2(t) = A$, and $y_3(t) = S + E + R$ are independent of $\tau$, at any $t$. This is obtained by proceeding as follows.

We considered several values of $\tau$ in $U_0$. For each $\tau$, we integrated the differential equation in (7.1) with the initial conditions given in (7.21) for that $\tau$ and by setting the unknown time-varying parameter as in (7.20) for the same value of $\tau$ (that is one of the profiles plotted in Figure 7.1). Both the initial conditions in (7.21) and the unknown input significantly depend on $\tau$. However, the three outputs (which are displayed in Figure 7.4) are independent, in accordance with our results.

$$7.6.2 \quad \beta(t) = \cos\left(\frac{\pi}{400}t\right)$$

We proceed exactly as in the previous case. In this case we are allowed to also include smaller $\tau$ and we set:

$$U_0 := [-0.10, 10].$$

<figure><img src="image_13.png" /><figcaption>Figure 7.5: Several indistinguishable profiles for the time-varying parameter $\beta'(t, \tau)$ for several values of the parameter $\tau$ ranging from -0.1 up to 10. The dashed purple line is the profile given by $\beta(t) = \cos(\frac{\pi}{400}t)$, i.e., $\beta'(t, \tau)$ for $\tau = 0$.</figcaption></figure>

Figure 7.5 displays the profiles of $\beta'(t, \tau)$. The profile for $\tau = 0$ (purple dashed line) is precisely the one set in (7.22). By varying $\tau$, we obtain a significant change of the profile meaning that this parameter is "strongly" unidentifiable.

Figures 7.6 and 7.7 display the profiles of $S'(t, \tau)$ and $E'(t, \tau)$, respectively.

Also in this case we verified that the three outputs, $y_1(t) = I$, $y_2(t) = A$, and $y_3(t) = S+E+R$ are independent of $\tau$, at any $t$. They are displayed in Figure 7.8.
<figure><img src="image_14.png" /><figcaption>Figure 7.6: The profiles $S'(t, \tau)$ for several values of the parameter $\tau$ ranging from -0.1 up to 10. The dashed purple line is the profile for $\tau = 0$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 7.7: The profiles $E'(t, \tau)$ for several values of the parameter $\tau$ ranging from -0.1 up to 10. The dashed purple line is the profile for $\tau = 0$.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 7.8: The three outputs $y_1(t) = I$ (solid blue line), $y_2(t) = A$ (dashed red line), and $y_3(t) = S + E + R$ (dotted yellow line). The three profiles are independent of the parameter $\tau$, as expected.</figcaption></figure>
# Chapter 8

## Genetic toggle switch

This example is provided to illustrate the use of our results when the dynamics of the system are not affine in the inputs. We study the observability and the identifiability properties of the ODE model characterized by the following dynamics and outputs:

$$ egin{cases} 
\dot{x}_1 = k_{01} + \frac{k_1}{1 + \left(\frac{x_2}{W_1}\right)^{n_1}} - x_1 \\ 
\dot{x}_2 = k_{02} + \frac{k_2}{1 + \left(\frac{x_1}{W_2}\right)^{n_2}} - x_2 \\ 
y = [x_1, x_2], 
\end{cases} 
\quad (8.1) $$

This model includes two TV-parameters, $W_1$ and $W_2$, and six constant parameters, $k_{01}$, $k_{02}$, $k_1$, $k_2$, $n_1$, $n_2$. All of them are assumed to be unknown. This model is commonly adopted to characterize the genetic toggle switch. The interested reader is addressed to Section 3.2 of [35] to find a biological description of this system, which includes the physical meaning of all the above states and parameters. Note that, (8.1) coincides with Eq. (3.5) in [35] with $n_1 = \eta_{TetR}$, $n_2 = \eta_{LacI}$, $W_1 = 1 + w_1$, and $W_2 = 1 + w_2$.

The model in (8.1) is a special case of (2.1) but not directly a special case of (2.2). In accordance with what we mentioned at the end of Section 2.1, we need to introduce a state that, in addition to the six constant parameters, also includes $W_1$ and $W_2$. In other words:

$$ x = [x_1, x_2, W_1, W_2, k_{01}, k_1, n_1, k_{02}, k_2, n_2]^T \quad (8.2) $$

and the dynamics are:

$$ egin{cases} 
\dot{x}_1 = k_{01} + \frac{k_1}{1 + \left(\frac{x_2}{W_1}\right)^{n_1}} - x_1 \\ 
\dot{x}_2 = k_{02} + \frac{k_2}{1 + \left(\frac{x_1}{W_2}\right)^{n_2}} - x_2 \\ 
\dot{W}_1 = w_1 \\ 
\dot{W}_2 = w_2 \\ 
\dot{k}_{01} = \dot{k}_1 = \dot{n}_1 = \dot{k}_{02} = \dot{k}_2 = \dot{n}_2 = 0. \\ 
y = [x_1, x_2], 
\end{cases} 
\quad (8.3) $$

The system is now a special case of (2.2). In particular, $m_u = 0$, $m_w = 2$, $p = 2$, $h_1(x) = x_1$, $h_2(x) = x_2$,

$$ g^0 = \begin{bmatrix} k_{01} + \frac{k_1}{1 + \left(\frac{x_2}{W_1}\right)^{n_1}} - x_1 \\ k_{02} + \frac{k_2}{1 + \left(\frac{x_1}{W_2}\right)^{n_2}} - x_2 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}, \quad g^1 = \begin{bmatrix} 0 \\ 0 \\ 1 \\ 0 \\ 0 \\ 0 \end{bmatrix}, \quad g^2 = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 1 \\ 0 \\ 0 \end{bmatrix} \quad (8.4) $$

### 8.1 Implementation of Algorithm 1

We run Algorithm 1 to obtain the observability codistribution. Line 2 provides:
$$ \Omega = \operatorname{span}\{\nabla h_1, \nabla h_2\} = \operatorname{span}\{[1, 0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0]\}. $$

$\deg_w(\Omega) = 0$ ($< m_w$), as $\mathcal{L}_{g^1}h_1 = \mathcal{L}_{g^1}h_2 = \mathcal{L}_{g^2}h_1 = \mathcal{L}_{g^2}h_2 = 0$ (and the rank of the unknown input reconstructability matrix from $h_1, h_2$ vanishes). The system is not in canonical form with respect to its unknown input. We enter the loop at lines 9-14.

Line 11 provides no functions ($m = 0$). Line 12 provides: $^m\mu_0^0 = ^m\nu_0^0 = 1$ obtained from (3.6) (these tensors have a single entry as the indices only take the value 0). In addition, from (3.7) we obtain: $^m\hat{g}^0 = ^m\nu_0^0g^0 = g^0$. Line 13 computes

$$ \Omega = \left\langle \left\{ \begin{array}{cc} g^0 & \\ g^1 & g^2 \end{array} \right\} \middle| \operatorname{span}\{\nabla h_1, \nabla h_2\} \right\rangle, $$

as $^m\hat{g}^0 = g^0$. By a direct computation, we obtain:

$$ \Omega = \operatorname{span}\{\nabla h_1, \nabla h_2, \nabla \mathcal{L}_{g^0}h_1, \nabla \mathcal{L}_{g^0}h_2\} = \\ \operatorname{span}\{\nabla h_1, \nabla h_2, \nabla h_3, \nabla h_4\} $$

where we adopted the following notation:

$$ h_3 := \mathcal{L}_{g^0}h_1, \quad h_4 := \mathcal{L}_{g^0}h_2 $$

Now, $\deg_w(\Omega) = 2 = m_w$. Hence, the second condition at Line 9 is not met and the loop is interrupted. In addition, as the condition at Line 15 is met, also the main loop is interrupted.

The execution of Algorithm 1 continues with Line 45. The operation at this line ($[\tilde{h}_1, \dots, \tilde{h}_{m_w}] = S(\Sigma, \Omega)$) sets

$$ \tilde{h}_1 = h_3 = k_{01} + \frac{k_1}{1 + \left(\frac{x_2}{W_1}\right)^{n_1}} - x_1, \quad \tilde{h}_2 = h_4 = k_{02} + \frac{k_2}{1 + \left(\frac{x_1}{W_2}\right)^{n_2}} - x_2 \quad (8.5) $$

Then, Line 46 provides:

$$ \mu = \begin{bmatrix} 1 & 0 & 0 \\ \mathcal{L}_{g^0}\tilde{h}_1 & \mathcal{L}_{g^1}\tilde{h}_1 & \mathcal{L}_{g^2}\tilde{h}_1 \\ \mathcal{L}_{g^0}\tilde{h}_2 & \mathcal{L}_{g^1}\tilde{h}_2 & \mathcal{L}_{g^2}\tilde{h}_2 \end{bmatrix}, $$

whose expression is obtained by using (8.4) and (8.5). From $\mu$ we obtain $\nu$ (by computing its inverse). Finally, from $\nu$, we obtain $\hat{g}^0, \hat{g}^1$, and $\hat{g}^2$, by using (3.4).

From (3.5) we obtain that, when $m_u = 0$, the codistribution $\tilde{\mathcal{O}}$ vanishes. The final step of Algorithm 1 is the execution of Line 47, namely:

$$ \mathcal{O} = \langle \hat{g}^0, \hat{g}^1, \hat{g}^2, \middle| \operatorname{span}\{\nabla h_1, \nabla h_2, \nabla h_3, \nabla h_4\} \rangle. $$

The algorithm in (2.8) converges at the first step ($\Omega_1 = \Omega_0$), and

$$ \mathcal{O} = \operatorname{span}\{\nabla h_1, \nabla h_2, \nabla h_3, \nabla h_4\}. $$

The rank of $\mathcal{O}$ is 4, which is smaller than the dimension of the state in (8.2). Hence, this state is not observable. Note that $x_1$ and $x_2$ are observable and we know this even without the implementation of Algorithm 1, as they are the outputs.

## 8.2 Identifiability analysis and indistinguishable states and parameters

Regarding the constant parameters, it is immediate to check that none is identifiable (it suffices to show that their gradient does not belong to $\mathcal{O}$). This result contradicts the result available in the state of the art (e.g., see Section 3.2 of [35]). In Section 8.3, we explicitly prove the validity of our result.

Let us study the identifiability of the time varying parameters.
First of all we remark that the two original TV-parameters, $W_1$ and $W_2$, were included in the state (defined in Eq. (8.2)). As result, their identifiability is obtained by simply checking if their differential belongs to $O$. It is immediate to check that this is not the case. This result contradicts the result available in Section 3.2 of [35]. In Section 8.3, we explicitly prove the validity of our result.

Let us check the identifiability of their time derivatives, i.e., $w_1$ and $w_2$.

We execute Algorithm 4. First, as the system is canonic with respect to its unknown inputs, we do not have symmetries of the unknown input due to the system non canonicity ($c^c\chi$).

Let us compute the distribution orthogonal to the observable codistribution. We have:

$$O^{\perp} = \operatorname{span}\{\xi_1, \xi_2, \xi_3, \xi_4, \xi_5, \xi_6\}. \qquad (8.6)$$

with:

$$\xi_1 = \begin{bmatrix} 0_2 \\ -\frac{W_1((x_2/W_1)^{n_1+1})^2}{n_1 k_1 (x_2/W_1)^{n_1}} \\ 0 \\ 1 \\ 0_5 \end{bmatrix}, \quad \xi_2 = \begin{bmatrix} 0_2 \\ -\frac{W_1((x_2/W_1)^{n_1+1})}{n_1 k_1 (x_2/W_1)^{n_1}} \\ 0_2 \\ 1 \\ 0_4 \end{bmatrix}, \quad \xi_3 = \begin{bmatrix} 0_2 \\ W_1 \log\left(\frac{x_2}{W_1}\right) \\ 0_3 \\ n_1 \\ 0_3 \end{bmatrix},$$

$$\xi_4 = \begin{bmatrix} 0_3 \\ -\frac{W_2((x_1/W_2)^{n_2+1})^2}{n_2 k_2 (x_1/W_2)^{n_2}} \\ 0_3 \\ 1 \\ 0_2 \end{bmatrix}, \quad \xi_5 = \begin{bmatrix} 0_3 \\ -\frac{W_2((x_1/W_2)^{n_2+1})}{n_2 k_2 (x_1/W_2)^{n_2}} \\ 0_4 \\ 1 \\ 0 \end{bmatrix}, \quad \xi_6 = \begin{bmatrix} 0_3 \\ W_2 \log\left(\frac{x_1}{W_2}\right) \\ 0_4 \\ 0 \\ n_2 \end{bmatrix}$$

Hence, we have six symmetries of the state. From them, it is immediate to compute the symmetries of the unknown inputs (i.e., by using Eq. (4.6)). It is easy to check that they do not vanish. We conclude that also $w_1$ and $w_2$ are unidentifiable, even locally.

As for the case study, starting from the generators of $O^{\perp}$ we compute the indistinguishable states (although the determination of the indistinguishable states is not requested by Algorithm 4). We need to solve the differential equation in (4.1).

As we have six generators, we obtain six distinct sets of solutions. We directly provide the results.

Let us denote the true state and parameters as follows: $x_1(t)$, $x_2(t)$, $W_1(t)$, $W_2(t)$, $w_1(t)$, $w_2(t)$, $k_{01}$, $k_1$, $n_1$, $k_{02}$, $k_2$, and $n_2$. The new parameters will be denoted by

$$W'_1(t), W'_2(t), w'_1(t), w'_2(t), k'_{01}, k'_1, n'_1, k'_{02}, k'_2, n'_2.$$

They will depend on a given parameter, denoted by $\tau$, that can take all the values in a given interval $(-\epsilon, \epsilon)$ (for a given $\epsilon \in \mathbb{R}^+$).

## First set

$$W_1(t)' = x_2 \left( \frac{1 - n_1 \left(1 + \left(\frac{x_2}{W_1(t)}\right)^{n_1}\right) \tau}{\left(\frac{x_2}{W_1(t)}\right)^{n_1} + n_1 \left(\frac{x_2}{W_1(t)}\right)^{n_1} \tau + n_1 \tau} \right)^{\frac{1}{n_1}}, \qquad (8.7)$$
$$k'_{01} = k_{01} + n_1 k_1 \tau,$$

and the remaining parameters do not change, i.e.: $W'_2(t) = W_2(t)$, $k'_1 = k_1$, $n'_1 = n_1$, $k'_{02} = k_{02}$, $k'_2 = k_2$, $n'_2 = n_2$.

## Second set

$$W'_1(t) = x_2 \left( \frac{1}{\left(1 + \left(\frac{x_2}{W_1(t)}\right)^{n_1}\right) e^{n_1 \tau} - 1} \right)^{\frac{1}{n_1}} \qquad (8.8)$$
$$k'_1 = k_1 e^{n_1 \tau}$$

and the remaining parameters do not change, i.e.: $W'_2(t) = W_2(t)$, $k'_{01} = k_{01}$, $n'_1 = n_1$, $k'_{02} = k_{02}$, $k'_2 = k_2$, $n'_2 = n_2$.
### Third set

$$W'_{1}(t) = x_{2} e^{\left(e^{-\tau} \log \frac{W_{1}(t)}{x_{2}}\right)} \qquad (8.9)$$

$$n'_{1} = n_{1}e^{\tau}$$

and the remaining parameters do not change, i.e.: $W'_{2}(t) = W_{2}(t)$, $k'_{01} = k_{01}$, $k'_{1} = k_{1}$, $k'_{02} = k_{02}$, $k'_{2} = k_{2}$, $n'_{2} = n_{2}$.

### Fourth set

As the first set, with the change $1 \leftrightarrow 2$.

### Fifth set

As the second set, with the change $1 \leftrightarrow 2$.

### Sixth set

As the third set, with the change $1 \leftrightarrow 2$.

## 8.3 Direct proof

We provide an unequivocal proof of the validity of our results. In particular, we directly prove that the values of the parameters in the six sets are indistinguishable from the true values.

The indistinguishability will be proved by showing that, for any $\tau$, the new parameters provide exactly the same dynamics of $x_1$ and $x_2$, obtained with the true parameters, i.e, from (8.1), by proving the following two equalities:

$$k'_{01} + \frac{k'_{1}}{1 + \left(\frac{x_{2}(t)}{W'_{1}(t)}\right)^{n'_{1}}} = k_{01} + \frac{k_{1}}{1 + \left(\frac{x_{2}(t)}{W_{1}(t)}\right)^{n_{1}}} \qquad (8.10)$$

$$k'_{02} + \frac{k'_{2}}{1 + \left(\frac{x_{1}(t)}{W'_{2}(t)}\right)^{n'_{2}}} = k_{02} + \frac{k_{2}}{1 + \left(\frac{x_{1}(t)}{W_{2}(t)}\right)^{n_{2}}} \qquad (8.11)$$

### 8.3.1 First set

By a direct substitution of (8.7) in (8.10), we obtain the validity of (8.10). The validity of (8.11) is trivial. The existence of this set of indistinguishable parameters proves the unidentifiability of $W_1(t)$ and $k_{01}$. In addition, as $w_1(t) = \dot{W}_1(t) \neq \dot{W}'_1(t) = w'_1(t)$, it also proves the unidentifiability of $w_1(t)$. Finally, as $\tau$ can take infinitesimal values, and $\lim_{\tau \to 0} W_1(t)' = W_1(t)$, and $\lim_{\tau \to 0} k'_{01} = k_{01}$, the unidentifiability is even local.

### 8.3.2 Second set

By a direct substitution of (8.8) in (8.10), we obtain the validity of the equality in (8.10). The validity of (8.11) is trivial. The existence of this set of indistinguishable parameters proves the unidentifiability of $W_1(t)$, $k_1$, and $w_1(t)$. Finally, as $\tau$ can take infinitesimal values, and $\lim_{\tau \to 0} W_1(t)' = W_1(t)$, and $\lim_{\tau \to 0} k'_1 = k_1$, the unidentifiability is even local.

### 8.3.3 Third set

By a direct substitution of (8.9) in (8.10), we obtain the validity of the equality in (8.10). The validity of (8.11) is trivial. The existence of this set of indistinguishable parameters proves the unidentifiability of $W_1(t)$, $n_1$, and $w_1(t)$. Finally, as $\tau$ can take infinitesimal values, and $\lim_{\tau \to 0} W_1(t)' = W_1(t)$, and $\lim_{\tau \to 0} n'_1 = n_1$, the unidentifiability is even local.
### 8.3.4 Fourth set

In this case the validity of (8.10) is trivial while to prove the validity of (8.11) we need to use (8.7) with the change $1 \leftrightarrow 2$. The existence of this set of indistinguishable parameters proves the unidentifiability of $W_2(t)$, $k_{02}$, and $w_2(t)$. Finally, as $\tau$ can take infinitesimal values, and $\lim_{\tau \to 0} W_2(t)' = W_2(t)$, and $\lim_{\tau \to 0} k_{02}' = k_{02}$, the unidentifiability is even local.

### 8.3.5 Fifth set

In this case the validity of (8.10) is trivial while to prove the validity of (8.11) we need to use (8.8) with the change $1 \leftrightarrow 2$. The existence of this set of indistinguishable parameters proves the unidentifiability of $W_2(t)$, $k_2$, and $w_2(t)$. Finally, as $\tau$ can take infinitesimal values, and $\lim_{\tau \to 0} W_1(t)' = W_2(t)$, and $\lim_{\tau \to 0} k_2' = k_2$, the unidentifiability is even local.

### 8.3.6 Sixth set

In this case the validity of (8.10) is trivial while to prove the validity of (8.11) we need to use (8.9) with the change $1 \leftrightarrow 2$. The existence of this set of indistinguishable parameters proves the unidentifiability of $W_2(t)$, $n_2$, and $w_2(t)$. Finally, as $\tau$ can take infinitesimal values, and $\lim_{\tau \to 0} W_2(t)' = W_2(t)$, and $\lim_{\tau \to 0} n_2' = n_2$, the unidentifiability is even local.
# Chapter 9

## Further simple examples

This chapter provides three simple examples. In particular, the former is a time varying system. We only study the observability of the state and the identifiability of the parameters, constant and time varying. In other words, we simply run Algorithms 1 and 4.

### 9.1 First Example

We consider the system characterized by the following equations:

$$ \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \\ y \end{bmatrix} = \begin{bmatrix} 1 \\ w_1 \\ x_1 + x_2 t \end{bmatrix} \qquad (9.1) $$

The system is directly a special case of (2.2). In particular, $m_u = 0$, $m_w = 1$, $p = 1$, $h_1(x) = x_1 + x_2t$,

$$ g^0 = \begin{bmatrix} 1 \\ 0 \end{bmatrix}, \quad g^1 = \begin{bmatrix} 0 \\ 1 \end{bmatrix}. \qquad (9.2) $$

#### 9.1.1 Observability analysis

We run Algorithm 1 to obtain the observability codistribution. Line 2 provides:

$$ \Omega = \operatorname{span}\{\nabla h_1\} = \operatorname{span}\{[1, t]\}. $$

$\deg_w(\Omega) = 1$ ($= m_w$), as $\mathcal{L}_{g^1}h_1 = t$ (and the unknown input reconstructability matrix from $h_1$ is full rank). The system is in canonical form with respect to its unknown input and the main loop is immediately interrupted. The execution of Algorithm 1 continues with Line 45. The operation at this line ($[\tilde{h}_1, \dots, \tilde{h}_{m_w}] = \mathcal{S}(\Sigma, \Omega)$) sets

$$ \tilde{h}_1 = h_1 = x_1 + x_2 t \qquad (9.3) $$

Then, Line 46 provides:

$$ \mu = \begin{bmatrix} \dot{\mathcal{L}}_{g^0} \tilde{h}_1 & 0 \\ 0 & \mathcal{L}_{g^1} \tilde{h}_1 \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 1 + x_2 & t \end{bmatrix}, \quad \nu = \begin{bmatrix} 1 & 0 \\ -\frac{1+x_2}{t} & \frac{1}{t} \end{bmatrix} $$

Finally, from $\nu$, we obtain $\hat{g}^0$, and $\hat{g}^1$, by using (3.4):

$$ \hat{g}^0 = g^0 - \frac{1+x_2}{t} g^1 = \begin{bmatrix} 1 \\ -\frac{1+x_2}{t} \end{bmatrix}, \quad \hat{g}^1 = \frac{1}{t} g^1 = \begin{bmatrix} 0 \\ \frac{1}{t} \end{bmatrix} $$

From (3.5) we obtain that, when $m_u = 0$, the codistribution $\tilde{\mathcal{O}}$ vanishes. The final step of Algorithm 1 is the execution of Line 47, namely:

$$ \mathcal{O} = \langle \hat{g}^0, \hat{g}^1 | \operatorname{span}\{\nabla h_1\} \rangle. $$
Note that, $\dot{\mathcal{L}}_{\hat{g}^0}h_1 = \mathcal{L}_{\hat{g}^0}h_1 + \frac{\partial h_1}{\partial t} = 0$, and $\mathcal{L}_{\hat{g}^1}h_1 = 1$. The algorithm in (2.8) converges at the first step ($\Omega_1 = \Omega_0$), and

$$ \mathcal{O} = \operatorname{span}\{\nabla h_1\} = \operatorname{span}\{[1, t]\}. $$

The rank of $\mathcal{O}$ is 1, which is smaller than the dimension of the state. Hence, the state is not observable.

## 9.1.2 Identifiability analysis

Let us check the identifiability of $w_1$.

We execute Algorithm 4. First, as the system is canonic with respect to its unknown inputs, we do not have symmetries of the unknown input due to the system non canonicity ($^c\chi$).

Let us compute the distribution orthogonal to the observable codistribution. We have:

$$ \mathcal{O}^{\perp} = \operatorname{span}\left\{\begin{bmatrix} -t \\ 1 \end{bmatrix}\right\}, \qquad (9.4) $$

We have a single symmetry of the state. It is immediate to compute the symmetry of the unknown input.

In accordance with Equation (4.6), we need to compute $\xi_1^0, \xi_1^1$. We have:

$$ \dot{\mathcal{L}}_{g^0} \tilde{h}_1 = \mathcal{L}_{g^0} \tilde{h}_1 + \frac{\partial \tilde{h}_1}{\partial t} = 1 + x_2, \quad \mathcal{L}_{g^1} \tilde{h}_1 = t. $$

From Equation (4.7), by using the generator of $\mathcal{O}^{\perp}$ in (9.4) for $\xi$, we obtain:

$$ \xi_1^0 = 1, \quad \xi_1^1 = 0, $$

and, from (4.6),

$$ u_{\chi_1} = -\frac{1}{t} \neq 0 $$

Hence, $w_1$ is not identifiable, even locally.

## 9.2 Second example

We study the system characterized by the following equations:

$$ \left\{ \begin{aligned} \dot{x}_1 &= x_3 \\ \dot{x}_2 &= x_3 w_1 \\ \dot{x}_3 &= w_2 + x_1 w_3 \\ y &= [x_1, x_2]. \end{aligned} \right. \qquad (9.5) $$

It is a special case of (2.2) with: $m_u = 0, m_w = 3, p = 2, h_1(x) = x_1, h_2(x) = x_2$,

$$ g^0 = \begin{bmatrix} x_3 \\ 0 \\ 0 \end{bmatrix}, \quad g^1 = \begin{bmatrix} 0 \\ x_3 \\ 0 \end{bmatrix}, \quad g^2 = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}, \quad g^3 = \begin{bmatrix} 0 \\ 0 \\ x_1 \end{bmatrix}. \qquad (9.6) $$

### 9.2.1 Observability analysis

We run Algorithm 1 to obtain the observability codistribution. Line 2 provides:

$$ \Omega = \operatorname{span}\{\nabla h_1, \nabla h_2\} = \operatorname{span}\{[1, 0, 0], [0, 1, 0]\}. $$

$\deg_w(\Omega) = 1$ ($< m_w$), as $\mathcal{L}_{g^1}h_1 = \mathcal{L}_{g^2}h_1 = \mathcal{L}_{g^3}h_1 = \mathcal{L}_{g^2}h_2 = \mathcal{L}_{g^3}h_2 = 0$, and $\mathcal{L}_{g^1}h_2 = x_3$ (and the rank of the unknown input reconstructability matrix from $h_1, h_2$ is 1). The system is not in canonical form with respect to its unknown input. We enter the loop at lines 9-14.

Line 11 provides: