# Algorithmic study of the algebraic parameter estimation problem for a class of perturbations

Maya Chartouny, Thomas Cluzeau, Alban Quadrat

► To cite this version:

Maya Chartouny, Thomas Cluzeau, Alban Quadrat. Algorithmic study of the algebraic parameter estimation problem for a class of perturbations. Maple Transactions, 2023, 3, 10.5206/mt.v2i2.14467. hal-04203089

HAL Id: hal-04203089

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

Submitted on 11 Sep 2023

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Algorithmic study of the algebraic parameter estimation problem for a class of perturbations

CHARTOUNY MAYA, Thales DIS - University of Paris-Saclay Versailles, France

CLUZEAU THOMAS, Univ. Limoges, CNRS, XLIM, UMR 7252, F-87000 Limoges, France

QUADRAT ALBAN, IMJ - PRG, Sorbonne University, Inria Paris-Ouragan project-team, France

We consider the algebraic parameter estimation problem for a class of standard perturbations. We assume that the measurement $z(t)$ of a solution $x(t)$ of a linear ordinary differential equation, whose coefficients depend on a set $\theta := \{\theta_1, \dots, \theta_r\}$ of unknown constant parameters, is affected by a perturbation $\gamma(t)$ whose structure is supposed to be known (e.g., an unknown bias, an unknown ramp), i.e., $z(t) = x(t, \theta) + \gamma(t)$. We investigate the problem of obtaining closed-form expressions for the parameters $\theta_i$'s in terms of repeated indefinite integrals or convolutions of $z$. We illustrate the different results with explicit examples computed using the NonA package, developed in Maple, in which we have implemented our main contributions.

Additional Key Words and Phrases: Parameter estimation problem, inverse Cauchy problem, algebraic systems, elimination, annihilators, rings of ordinary differential operators.

**Recommended Reference Format:**

Chartouny Maya, Cluzeau Thomas, and Quadrat Alban. 2022. Algorithmic study of the algebraic parameter estimation problem for a class of perturbations. *Maple Trans.* 2, 2, Article 14467 (November 2022), 26 pages.
https://doi.org/10.5206/mt.v2i2.14467

## 1 Introduction

In many real applications studied in engineering sciences, applied mathematics, mathematical biology, etc., the estimation of constant parameters of a dynamical system (e.g., a mass, a spring constant, a damping coefficient, a resistance) is a fundamental issue. Hence, parameter estimation has widely been studied in different scientific fields. In particular, this problem has been actively studied in both control theory and signal processing. See, e.g., [Kailath et al. (2000), Fliess et al. (2003), Poor (1994), Van Trees (2004)].

In the present paper, we study the *algebraic parameter estimation problem* following the approach initiated in [Fliess et al. (2003)] and continued in [Belkoura et al. (2009), Mboup (2009), Quadrat (2017), Ushirobira et al. (2016)]. For more references, see the references therein. This problem aims at estimating the constant parameters of the solutions of linear ordinary differential equations (ODEs) with polynomial coefficients. More precisely, if we assume that $x(t)$ is a solution of a linear ODE depending on a set $\theta := \{\theta_1, \dots, \theta_r\}$ of unknown constant parameters, then we look for explicit expressions for the parameters $\theta_i$'s in terms of a measured function $z(t) = x(t, \theta) + \gamma(t) + \omega(t)$, where $\gamma$ is a perturbation (whose global structure can sometimes be supposed to be known) and $\omega$ is a noise (e.g., a zero-mean Gaussian noise). Moreover, we are interested in closed-form expressions of the parameters $\theta_i$'s involving repeated indefinite integrals or convolutions of $z$ which allow us to filter a part of the influence of noise $\omega$.

Authors' addresses: Chartouny Maya, Thales DIS - University of Paris-Saclay Versailles, France, mayachartouny@hotmail.com; Cluzeau Thomas, Univ. Limoges, CNRS, XLIM, UMR 7252, F-87000 Limoges, France, thomas.cluzeau@unilim.fr; Quadrat Alban, IMJ - PRG, Sorbonne University, Inria Paris-Ouragan project-team, Paris, France, alban.quadrat@inria.fr.

Permission to make digital or hard copies of all or part of this work for any use is granted without fee, provided that copies bear this notice and the full citation on the first page. Copyrights for third-party components of this work must be honored. © 2022 Copyright held by the owner/author(s). Publication rights licensed to Maple Transactions, under Creative Commons CC-BY 4.0 License.

https://doi.org/10.5206/mt.v2i2.14467
Note that the standard estimation approaches yield an estimation $\hat{\theta}$ of the parameters which asymptotically converges to $\theta$. In this sense, these methods are “asymptotic methods”. They mainly use the resolution of an optimization problem. The first main interest of the approach studied in this paper is that it gives an explicit exact formula for $\theta$. This formula can then be used to develop real-time numerical schemes for the estimation of $\theta$ [Fliess et al. (2003), Mboup (2009)]. In this sense, this approach is a “non-asymptotic method”. The noise sensitivity analysis of the real-time numerical schemes can be studied afterward. A second interest of this approach is that its formulation is exact and algebraic, and thus, the use of algebraic methods is natural to investigate the possibility to obtain closed-form solutions for $\theta$. Hence, the algebraic parameter estimation problem belongs to the realm of computer algebra. The first goal of this paper is to popularize the algebraic parameter estimation problem among the computer algebra community.

The first case to be considered in the study of the algebraic parameter estimation problem is that of a perfect measurement, i.e., $\gamma(t) = \omega(t) = 0$. As explained in [Chartouny et al. (2021)], this can be seen as an *inverse Cauchy problem*. Indeed, the *Cauchy problem* characterizes the solutions of an ODE that satisfies fixed initial conditions. Conversely, given a function that is known to satisfy a linear ODE of fixed order with polynomial coefficients, the inverse Cauchy problem studies when the constant coefficients of these polynomials as well as the initial conditions can be expressed with repeated indefinite integrals or convolutions of the solution. To achieve our goal, we have to combine *operational calculus* (Laplace transform, convolution product - see Appendix A.1) with *algebraic methods* (elimination methods, linear algebra). This leads to Theorems 2.3 and 2.5, first obtained in [Chartouny et al. (2021)], where it is shown that if $t = 0$ is an ordinary point of the ODE, then we can obtain explicit expressions for the coefficients of the ODE and the initial conditions in terms of repeated indefinite integrals of the measured solution.

The second goal of the paper is to consider the case where the signal $x$ is corrupted by a perturbation $\gamma$ so that we measure $z(t) := x(t) + \gamma(t)$. Here, the function $\gamma$ is a *structured perturbation* of the form $\gamma(t) = c t^r H(t)$, where $c$ is an unknown constant, $r \in \mathbb{Z}_{\ge 0}$ and $H$ is the Heaviside function. These structured perturbations are standard, e.g., in the *disturbance rejection problem* (see, e.g., [Quadrat et al. (2014)] and the references therein). For instance, if $r = 0$, then $\gamma(t) = c H(t)$ is a *step function* which corresponds to a unknown *bias*. If $r = 1$, then $\gamma(t) = c t H(t)$ is a *ramp function* and if $r = 2$, then $\gamma(t) = c t^2 H(t)$ is a *parabolic function*. We tackle the problem following the approach of [Fliess et al. (2003)] (see also [Mboup (2009), Quadrat (2017), Ushirobira et al. (2016)], which corresponds to the particular case $r = 0$). Applying the Laplace transform to the ODE, we obtain a new equation in the frequency domain which involves the initial conditions, the coefficients of the ODE, and the unknown constant parameter $c$. Using Gröbner basis techniques (see Appendix A.3) to compute the annihilator of a polynomial in the first Weyl algebra (see Appendix A.2), we manage to get rid of the term involving the unknown constant parameter $c$. Then, we use elimination techniques to obtain explicit expressions of the initial conditions in terms of the Laplace transform of $z(t)$, its derivatives, and the coefficients of the ODE. Applying the inverse Laplace transform, we can obtain explicit expressions of the initial conditions in terms of $z(t)$, its repeated indefinite integrals (convolutions), and the coefficients of the ODE. This is the first main result of the paper, stated in Theorem 3.8. Finally, as in the case of a perfect measurement, the second main result of the paper, stated in Theorem 3.10, shows that we can use elimination techniques to obtain the desired explicit expressions of the coefficients of the ODE in terms of $z(t)$ and its repeated indefinite integrals (convolutions).
We illustrate our results by explicit examples where the computations were done using the NonA package, developed in Maple, in which we have implemented the results of this paper. Details on the package as well as a complete worked example are given in Section 4.3.

The paper is organized as follows. In Section 2, we recall the main results of [Chartouny et al. (2021)] concerning the case of a perfect measurement. Section 3 contains the main contribution of the present paper for a class of standard structured perturbations. Section 4 contains the details concerning the implementation in Maple of the results developed in the paper and a worked example that illustrates the main commands. Finally, all the mathematical prerequisites (Laplace transform, convolution product, Weyl algebra, Gröbner basis, annihilator) are recalled in Section A.

## 2 The algebraic parameter estimation problem without perturbation

In this section, we first recall the approach developed in [Chartouny et al. (2021)] for the exact algebraic parameter estimation problem, i.e., for the case without perturbation, and we then state again the main results obtained in [Chartouny et al. (2021)], namely, Theorems 2.3, 2.4 and 2.5. These results will be generalized in Section 3 to include an important class of perturbations.

As explained in the introduction, the general parameter estimation problem consists in estimating constant parameters $\theta := (\theta_1, \dots, \theta_r)$, defining a solution $x$ of a Cauchy problem for an ODE depending on $\theta$, from the measurement of a signal of the form

$$z(t) = x(\theta, t) + \gamma(t) + \omega(t), \quad (1)$$

where $\gamma$ is a perturbation (whose global structure is sometimes supposed to be known) and $\omega$ is a noise (e.g., a zero-mean Gaussian noise). It is natural to first ask whether or not the system parameters $\theta$ can always be estimated in the exact case, i.e., when $\gamma = 0$ and $\omega = 0$. Indeed, if it is not possible, then the algebraic parameter estimation problem cannot be solved. Hence, in this section, we ask whether or not the coefficients of a linear ODE with polynomial coefficients and the initial conditions of the Cauchy problem can be exactly recovered from the knowledge of a “generic solution” and its repeated indefinite integrals. In other words, we aim to explicitly study the inverse Cauchy problem for linear ODEs with polynomial coefficients.

### 2.1 Estimation of the initial conditions

Let $\mathbb{K}$ be a field of characteristic 0. For instance, $\mathbb{K} = \mathbb{Q}$ or $\mathbb{K} = \mathbb{Q}(\theta_1, \dots, \theta_r)$ is the field of rational functions in the indeterminates $\theta_1, \dots, \theta_r$. Let us suppose that $x$ satisfies the following ODE

$$\sum_{i=0}^{n} a_i(t) x^{(i)}(t) = 0, \quad n \ge 1, \quad \forall i = 0, \dots, n, \quad a_i(t) = \sum_{j=0}^{d_i} a_{ij} t^j, \quad a_n \ne 0, \quad (2)$$

where $x^{(i)}$ stands for the $i^{th}$ derivative of $x$ with respect to $t$. Let us set $m := \max_{0 \le i \le n} d_i$. We first suppose that the coefficients $a_{ij}$'s are known, i.e., $a_{ij} \in \mathbb{K}$, and we focus on the possibility to recover the initial conditions of the Cauchy problem for (2), i.e., the values $x^{(i)}(0)$ for $i = 0, \dots, n-1$. In Section 2.2, we shall show how the coefficients $a_{ij}$'s can be explicitly determined from $x$ and its repeated indefinite integrals.

We first apply the Laplace transform to (2) and using the notation $\hat{x} = \mathcal{L}(x)$, we obtain

$$\sum_{i=0}^{n} a_i(-\partial_s) \left( s^{i} \hat{x}(s) - \sum_{j=0}^{i-1} s^{i-j-1} x^{(j)}(0) \right) = 0,$$
where $\partial_s = \frac{d}{ds}$ is the differential operator with respect to $s$ and $s$ is the Laplace variable (see Appendix A). This can be rewritten as

$$ \sum_{i=0}^{n} a_i(-\partial_s) s^i \hat{x}(s) - \sum_{i=0}^{n} \sum_{j=0}^{i-1} a_i(-\partial_s) s^{i-j-1} x^{(j)}(0) = 0. \quad (3) $$

Expanding the second term of the above identity, we get

$$ \sum_{i=0}^{n} \sum_{j=0}^{i-1} a_i(-\partial_s) s^{i-j-1} x^{(j)}(0) = \sum_{j=0}^{n-1} \sum_{i=j+1}^{n} a_i(-\partial_s) s^{i-j-1} x^{(j)}(0) = \sum_{k=0}^{n-1} \sum_{i=k+1}^{n} a_i(-\partial_s) s^{i-k-1} x^{(k)}(0). $$

Then, (3) becomes

$$ \sum_{i=0}^{n} a_i(-\partial_s) s^i \hat{x}(s) - \sum_{k=0}^{n-1} \sum_{i=k+1}^{n} a_i(-\partial_s) s^{i-k-1} x^{(k)}(0) = 0. $$

In what follows, we shall note

$$ \left\{ \begin{array}{l} P(s, \partial_s) := \sum_{i=0}^{n} a_i(-\partial_s) s^i, \\ S_k(s) := - \sum_{i=k+1}^{n} a_i(-\partial_s) s^{i-k-1}, \quad \partial_k := x^{(k)}(0), \quad k = 0, \dots, n-1, \\ Q(s) := \sum_{k=0}^{n-1} S_k \partial_k. \end{array} \right. \quad (4) $$

Notice first that $P(s, \partial_s)$ is an ordinary differential (OD) operator in $\partial_s$ with polynomial coefficients in $s$, i.e., is an element of the first Weyl algebra $A_1(\mathbb{K}) := \mathbb{K}[s] \langle \partial_s | \partial_s s = s \partial_s + 1 \rangle$ (see, e.g., [Coutinho (1995)] and Appendix A), which applies to the function $\hat{x}(s)$. Moreover, in the term $S_k$, $a_i(-\partial_s) \in A_1(\mathbb{K})$ applies to $s^{i-k-1}$, which shows that $S_k \in \mathbb{K}[s]$ and $Q \in \mathbb{K}[\vartheta_0, \dots, \vartheta_{n-1}][s]$. Equation (3) then becomes

$$ P(s, \partial_s) \hat{x}(s) + Q(s) = 0. \quad (5) $$

Let us now study when the $\partial_k$'s can be explicitly characterized. If we note $\Theta := (\vartheta_0 \dots \vartheta_{n-1})^T$, then (5) can be rewritten as

$$ (S_0 \dots S_{n-1}) \Theta = -P(s, \partial_s) \hat{x}(s). \quad (6) $$

**Definition 2.1.** The valuation of $a_i \in \mathbb{K}[t] \setminus \{0\}$ at $t = 0$, denoted by $v_0(a_i)$, is the degree of the maximum power of $t$ that divides $a_i$. Moreover, we set $v(0) = -\infty$.

In general, for all $k = 0, \dots, n-1$, we have $\deg_S(S_k) = \max_{i=k+1, \dots, n} \{i-k-1-v_0(a_i)\}$. However, in what follows, we shall suppose that $v_0(a_n) = 0$, i.e., $a_{n0} \neq 0$, which then implies $\deg_S(S_k) = n-k-1$ for $k = 0, \dots, n-1$. The latter assumption $a_{n0} \neq 0$ means that $t = 0$ is a regular point of (2), and thus, it makes sense to consider the Cauchy problem for (2) at $t = 0$.

**Example 2.2.** If $t = 0$ is a singular point for (2), i.e., $a_{n0} = 0$, then we have $a_n(t) = p(t)t$ for a certain polynomial $p$, and thus, $S_{n-1} = -a_n(-\partial_s)s^0 = p(-\partial_s)\partial_s 1 = 0$ and $Q = \sum_{k=0}^{n-2} S_k \partial_k$. Hence, $\vartheta_{n-1} = x^{(n-1)}(0)$ cannot be estimated as it does not appear in $Q$. Similarly, if $v_0(a_n) \ge 2$ and $v_0(a_{n-1}) \ge 1$ (e.g., $t^2 x^{(2)}(t) + t x^{(1)}(t) + x(t) = 0$), then $S_{n-1} = 0$ and $S_{n-2} = -(a_n(-\partial_s)s + a_{n-1}(-\partial_s)) = 0$, which shows that both $\vartheta_{n-1}$ and $\vartheta_{n-2}$ cannot be estimated since they do not appear in $Q$.

If $t = 0$ is not an ordinary point of (2), then note that a change of independent variable $t = T + \alpha$, $\alpha \in \mathbb{K}^*$, cannot simply be used to solve our problem. Indeed applying the process described above to the new ODE at $T = 0$ will provide the values of $\theta^{(i)}(\alpha)$ but not $\theta^{(i)}(0)$.

Equation (6) is an inhomogeneous linear equation in the $\vartheta_i$'s. If we differentiate $n-1$ times (6) with respect to $s$, we get the following inhomogeneous linear system for the $\vartheta_i$'s

$$ U \Theta = -(P(s, \partial_s) \quad \partial_s P(s, \partial_s) \quad \dots \quad \partial_s^{n-1} P(s, \partial_s))^T \hat{x}(s), $$
where the matrix $U$ is defined by

$$ U = \begin{pmatrix} S_0 & \cdots & S_{n-1} \\ S'_0 & \cdots & S'_{n-1} \\ \vdots & \ddots & \vdots \\ S^{(n-1)}_0 & \cdots & S^{(n-1)}_{n-1} \end{pmatrix}. $$

The polynomials $S_k$'s are defined in (4) and, under the assumption $a_{n0} \neq 0$, we have that $U$ is an invertible upper “anti-triangular” matrix, namely,

$$ U = \begin{pmatrix} S_0 & S_1 & \cdots & S_{n-3} & S_{n-2} & -a_{n0} \\ S'_0 & S'_1 & \cdots & S'_{n-3} & -a_{n0} & 0 \\ S''_0 & S''_1 & \cdots & -2a_{n0} & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ S^{(n-2)}_0 & -(n-2)! a_{n0} & 0 & \cdots & \cdots & 0 \\ -(n-1)! a_{n0} & 0 & \cdots & \cdots & 0 & 0 \end{pmatrix}, $$

and we get

$$ \Theta = \begin{pmatrix} \vartheta_0 \\ \vdots \\ \vartheta_{n-1} \end{pmatrix} = -U^{-1} \begin{pmatrix} P \\ \partial_s P \\ \vdots \\ \partial_s^{n-1} P \end{pmatrix} \hat{x}(s). \quad (7) $$

Note that we can give an explicit expression for the inverse $U^{-1}$ of the matrix $U$ so that there is no need to compute the inverse of a matrix to get the above expression of $\Theta$. For more details, see [Chartouny et al. (2021)]. We then have the following result.

**THEOREM 2.3** ([Chartouny et al. (2021)]). If $t = 0$ is a regular point of an ordinary differential equation (2), then the initial conditions $\{x^k(0)\}_{k=0,...,n-1}$ can be expressed explicitly in terms of $\hat{x}(s)$, its derivatives, and the coefficients $a_{ij}$'s of (2).

Let us set $U^{-1} = (c_{ij})$. Note that the $c_{ij}$'s are polynomial expressions of the $a_{ij}$'s and $a_{n0}^{-1}$. Equation (7) yields

$$ \vartheta_i = - \underbrace{\sum_{j=1}^{n} c_{ij} \partial_s^{j-1} P(s, \partial_s) \hat{x}(s)}_{R_i(s, \hat{x}(s))}, \quad i = 0, \dots, n-1, \quad (8) $$

for the initial conditions $\vartheta_i$'s. Moreover, in [Chartouny et al. (2021)], it has been shown that, for all $i = 0, \dots, n-1$, $\deg(s) = n + i$.

We shall now apply the inverse Laplace transform to get explicit expressions in terms of $x(t)$ and its indefinite repeated integrals (e.g., convolutions). To get formulas for the $\vartheta_i$'s involving indefinite repeated integrals of $x(t)$, we first rewrite the right-hand side of (8) as

$$ \vartheta_i = \frac{N_i}{D_i}, \quad N_i(s, \hat{x}(s)) := \frac{R_i(s, \hat{x}(s))}{s^{n+i+1}}, \quad D_i(s) := \frac{1}{s^{n+i+1}}. $$

Then, we apply the inverse Laplace transform $\mathcal{L}^{-1}$ and since $\vartheta_i$ is a constant and $\mathcal{L}^{-1}$ is a linear transformation, we obtain $\mathcal{L}^{-1}(\vartheta_i D_i) = \mathcal{L}^{-1}(N_i)$ so that

$$ \vartheta_i = \frac{\mathcal{L}^{-1}(N_i)}{\mathcal{L}^{-1}(D_i)}, \quad i = 0, \dots, n-1, \quad \mathcal{L}^{-1}(D_i)(t) = \frac{t^{n+i}}{(n+i)!}. $$
Finally, note that the term in the inverse Laplace transform of the right-hand side of the following equation

$$ \mathcal{L}^{-1}(N_i) = - \sum_{j=1}^{n} \mathcal{L}^{-1} \left( \frac{c_{ij}}{s^{n+i+1}} \partial_s^{j-1} P(s, \partial_s) \hat{x}(s) \right), \quad i = 0, \dots, n-1, $$

is a strictly proper rational function in s, namely, the degree in s of its numerator is strictly less than the degree of its denominator. Hence, using the normal forms of OD operators, we have

$$ \frac{c_{ij}}{s^{n+i+1}} \partial_s^{j-1} P(s, \partial_s) \hat{x}(s) = \sum_{\substack{0 \le k \le n+i+1 \\ 0 \le l \le m+j-1}} \frac{d_{kl}}{s^k} \partial_s^l \hat{x}(s), $$

where all the $d_{kl}$'s are polynomial expressions of the $a_{ij}$'s and $a_{n0}^{-1}$. Using $\mathcal{L} \left( \int_0^t y(\tau) d\tau \right) (s) = s^{-1} \hat{y}(s)$, we get that $\mathcal{L}^{-1}(N_i)$ is a finite sum of repeated indefinite integrals of terms of the form $(-t)^l x(t)$, which can also be expressed as a convolution. The $\vartheta_i$'s can thus be expressed linearly in terms of convolutions of $x$ and depend polynomially on the coefficients $a_{ij}$'s of the ODE (2), $a_{n0}^{-1}$ and $t^{-1}$. We then have the following result.

**THEOREM 2.4** ([Chartouny et al. (2021)]). If $t = 0$ is a regular point of an ordinary differential equation (2), then all the initial conditions $\{x^{(k)}(0)\}_{k=0,...,n-1}$ can explicitly be expressed in terms of convolutions of $x$ and of the coefficients $a_{ij}$'s of the ODE (2).

## 2.2 Estimation of the coefficients of the ODE

The closed-forms for the initial conditions $\{x^k(0)\}_{k=0,...,n-1}$'s obtained in Theorem 2.4 above depend on the constant parameters $a_{ij}$'s of the ODE (2). We now focus on the estimation of the parameters $a_{ij}$'s as repeated indefinite integrals of $x(t)$.

In the frequency domain, with the notations (4), we recall that the ODE defined by (2) is equivalently defined by (5) which can be rewritten as (6). In Section 2.1 above, we have differentiated $(n-1)^{\text{th}}$ times Equation (6) to get $U \Theta = -(P(s, \partial_s) \partial_s P(s, \partial_s) \dots \partial_s^{n-1} P(s, \partial_s))^T \hat{x}(s)$, from which we were able to solve for the $\vartheta_i$'s. From the last row of the matrix $U$, we see that if we differentiate Equation (6) more than $n$ times, then we shall get equations of the form $\partial_s^l P(s, \partial_s) \hat{x}(s) = 0$ for $l \ge n$ that do not depend on the $\vartheta_i$'s. For $l \ge n$, we have

$$ (-\partial_s)^l P(s, \partial_s) \hat{x}(s) = (-\partial_s)^l \sum_{\substack{0 \le i \le n \\ 0 \le j \le m}} a_{ij} (-\partial_s)^j s^i \hat{x}(s) = \sum_{\substack{0 \le i \le n \\ 0 \le j \le m}} (-\partial_s)^{l+j} s^i (\hat{x}(s) a_{ij}) = 0. $$

Hence, if we set $\vec{a} := (a_{00} \ a_{01} \ \dots \ a_{ij} \ \dots \ a_{n(m-1)} \ a_{nm})^T$, we get the following system of linear equations for the $a_{ij}$'s: for all $l \ge n$,

$$ \underbrace{\begin{pmatrix} (-\partial_s)^n & (-\partial_s)^{n+1} & \dots & (-\partial_s)^{n+j} s^i & \dots & (-\partial_s)^{n+m-1} s^n & (-\partial_s)^{n+m} s^n \\ \vdots & \vdots & \dots & \vdots & \dots & \vdots & \vdots \\ (-\partial_s)^l & (-\partial_s)^{l+1} & \dots & (-\partial_s)^{l+j} s^i & \dots & (-\partial_s)^{l+m-1} s^n & (-\partial_s)^{l+m} s^n \end{pmatrix}}_{V(s, \partial_s)} \hat{x}(s) \vec{a} = 0. \quad (9) $$

To obtain explicit formulas for the $a_{ij}$'s in terms of $\hat{x}(s)$ and its derivatives, we then need to solve the homogeneous linear system (9).

Let us first assume that $a_{nm} \ne 0$ and solve the system for the remaining $a_{ij}$'s. Let $p := (n+1)(m+1)-1$ be the number of unknown $a_{ij}$'s and take $l = n+p-1$ so that (9) provides $p$ equations in $p$ unknowns. Let us denote by $V_p(s, \partial_s)$ the $p \times p$ matrix defined by taking all but the last column of $V(s, \partial_s)$, $b(s, \partial_s)$
the last column of $V(s, \partial_s)$, and $\vec{a}_p$ the column vector of size $p$ formed by all but the last entry of $\vec{a}$.
Hence, (9) is equivalent to

$$ (V_p(s, \partial_s) \hat{x}(s)) \vec{a}_p = -a_{nm} b(s, \partial_s) \hat{x}(s), $$

so that if the matrix $V_p(s, \partial_s) \hat{x}(s)$ is invertible, i.e., $\det(V_p(s, \partial_s) \hat{x}(s)) \neq 0$, then we get

$$ \vec{a}_p = -a_{nm} (V_p(s, \partial_s) \hat{x}(s))^{-1} b(s, \partial_s) \hat{x}(s). \quad (10) $$

From (10), we obtain that each $a_{ij}$ can be written as a fraction $n_{ij}(s, \hat{x}(s))/d_{ij}(s, \hat{x}(s))$, where $n_{ij}$ and $d_{ij}$ are polynomials in $s, \hat{x}(s)$ and its derivatives. We now proceed as we did in Section 2.1 to apply the inverse Laplace transform and obtain formulas depending on $x(t)$ and its repeated indefinite integrals. Setting $q_{ij} := \max\{\deg_s n_{ij}, \deg_s d_{ij}\} + 1$ and defining $n'_{ij}(s, \hat{x}(s)) := n_{ij}(s, \hat{x}(s))/s^{q_{ij}}$ and $d'_{ij}(s, \hat{x}(s)) := d_{ij}(s, \hat{x}(s))/s^{q_{ij}}$, we obtain two polynomials in $s^{-1}$ and some derivatives of $\hat{x}(s)$. Applying the inverse Laplace transform to $n'_{ij}$ and $d'_{ij}$, and using the fact that the inverse Laplace transform maps a product to a convolution (see Proposition A.3 in the appendix), the $a_{ij}$'s are then ratios of sums of repeated indefinite integrals of terms of the form $(-t)^\alpha x(t) = \mathcal{L}^{-1}(\partial_s^\alpha \hat{x}(s)) (t)$, i.e., ratios of two convolutions depending only on $x(t)$. The $a_{ij}$'s can thus be expressed linearly in terms of $a_{nm}$ and as ratios of convolutions of $x$.

Now, if $\det(V_p(s, \partial_s) \hat{x}(s)) = 0$, then there exists a non-zero constant vector $\vec{c} := (c_{00} \ c_{01} \ \dots \ c_{0p})^T$ such that $V_p(s, \partial_s) \hat{x}(s) \vec{c} = 0$, and thus, we get equations of the form

$$ (-\partial_s)^k \sum_{(i,j) \in \{(0, \dots, n) \times \{0, \dots, m\} \} \setminus \{(n, m)\}} c_{ij} (-\partial_s)^j s^i \hat{x}(s) = 0, \quad k = n, \dots, l. $$

Using the inverse Laplace transform, the fact that $\mathcal{L}^{-1}(s^i \hat{x}(s)) = x^{(i)}(t) + \sum_{k=0}^{i-1} \delta^{(i-k-1)} x^{(k)}(0)$, where $\delta^{(e)}$ denotes the $e^{th}$ derivative of the Dirac distribution at $t = 0$, and $t^j \delta^{(e)} = 0$ for $j > e$, we obtain that $x$ then satisfies the following ODE

$$ t^n \left( \sum_{(i,j) \in \{(0, \dots, n) \times \{0, \dots, m\} \} \setminus \{(n, m)\}} c_{ij} t^j x^{(i)}(t) \right) = 0. $$

Hence, (10) holds if $x$ is a generic solution of (2), namely, a function which does not satisfy a lower order/degree ODE.

**THEOREM 2.5** ([Chartouny et al. (2021)]). Let us consider the ODE (2) with $a_{nm} \neq 0$. If $x$ is a generic solution of (2), then the coefficients $a_{ij}$'s can be expressed as ratios of convolutions of $x$ and linear expressions of $a_{nm}$.

The explicit expressions of the $a_{ij}$'s defined in Theorem 2.5 can then be substituted into the expressions of the initial conditions $x^{(k)}(0)$'s obtained in Theorem 2.4 to get explicit expressions of the $x^{(k)}(0)$'s as rational expressions of convolutions of $x$.

Finally, by assumption, we have $a_{ij} \in \mathbb{K}$. Hence, if $\mathbb{K} = \mathbb{Q}(\theta_1, \dots, \theta_r)$, i.e., if the $a_{ij}$'s are rational functions of the parameters $\theta_k$'s defining the signal $x$, i.e., $a_{ij} = f_{ij}(\theta_1, \dots, \theta_r)$, where $f_{ij}$ is rational function for $i = 0, \dots, n$ and $j = 0, \dots, m$, then the problem of expressing the $\theta_k$'s in terms of $x$ is reduced to expressing $\theta_k$'s in terms of the $a_{ij}$, i.e., to solving the rational system $f_{ij}(\theta_1, \dots, \theta_r) = a_{ij}$ for $i = 0, \dots, n$ and $j = 0, \dots, m$. This last problem is out of the scope of this paper.
### 3 The algebraic parameter estimation problem for a class of perturbations

In this section, we consider a signal $x$ which satisfies the following ODE

$$ \sum_{i=0}^{n} a_i(t) x^{(i)}(t) = 0, \quad n \ge 1, \quad \forall i = 0, \dots, n, \quad a_i \in \mathbb{K}[t], \quad a_n \ne 0, \qquad (11) $$

where $x^{(i)}$ stands for the $i^{th}$ derivative of $x$ with respect to $t$. We further assume that $t = 0$ is an ordinary point of the differential equation, i.e., $a_n(0) \ne 0$. The signal $x$ is now corrupted by a perturbation $\gamma$ so that we measure $z(t) := x(t) + \gamma(t)$. In what follows, we shall consider the case of a *structured perturbation* of the form

$$ \gamma(t) = c \, t^r \, H(t), \qquad (12) $$

where $c$ is an unknown constant, $r \in \mathbb{Z}_{\ge 0}$ and $H$ is the Heaviside function defined by $H(t) = 1$ for $t \ge 0$ and $0$ elsewhere. These structured perturbations are standard, e.g., in the *disturbance rejection problem* (see, e.g., [Quadrat et al. (2014)] and the references therein). For instance, if $r = 0$, then $\gamma(t) = c H(t)$ is a *step function* which corresponds to a *unknown bias*. If $r = 1$, then $\gamma(t) = c \, t \, H(t)$ is a *ramp function* and if $r = 2$, then $\gamma(t) = c \, t^2 \, H(t)$ is a *parabolic function*.

In Section 3.1, we generalize the approach developed in Section 2 to take into account the perturbation term $\gamma$ defined by (12) and we give the mathematical formulation of the corresponding algebraic parameter estimation problem. To remove the term containing the unknown constant $c$, we have to compute OD operators which annihilate a certain polynomial in $s$. In Section 3.2, we fully characterize the left ideal of $A_1(\mathbb{K})$ defined by all the OD operators with coefficients in $\mathbb{K}[s]$ which annihilate a general polynomial $P$ of $\mathbb{K}[s]$, i.e., the left annihilator $\text{ann}_{A_1(\mathbb{K})}(P)$ of $P$. In Section 3.3, using the above results, we then show how the approach stated in Section 2.1 can be generalized to characterize the initial conditions of the algebraic parameter estimation problem for the class of perturbations (12). Finally, in Section 3.4, we generalize the results given in Section 2.2 to characterize the $a_i$'s of (11) while $x$ is corrupted by the class of perturbation $\gamma$ defined by (12).

#### 3.1 Reformulation of the algebraic parameter estimation problem

We shall study how both the initial conditions $x^{(i)}(0)$, for $i = 0, \dots, n-1$, and the coefficients of the polynomials $a_i$'s can be expressed in terms of the measured function $z$ and its repeated indefinite integrals (e.g., convolutions). To achieve our goal, we first apply the Laplace transform to (11) and if we denote by $\hat{x}(s)$ the Laplace transform of $x$, where $s \in \mathbb{C}$ is the Laplace variable, we obtain

$$ \sum_{i=0}^{n} a_i(-\partial_s) \left( s^i \hat{x}(s) - \sum_{j=0}^{i-1} s^{i-j-1} x^{(j)}(0) \right) = 0, $$

which can be rewritten as

$$ \underbrace{\sum_{i=0}^{n} (a_i(-\partial_s) s^i)}_{R(s, \partial_s)} \hat{x}(s) + \underbrace{\sum_{k=0}^{n-1} \left( \underbrace{\sum_{i=k+1}^{n} -a_i(-\partial_s) s^{i-k-1}}_{S_k(s)} \right)}_{Q(s)} \underbrace{x^{(k)}(0)}_{\partial_k} = 0. $$

In terms of $\hat{z} = \hat{x} + \hat{\gamma}$, the latter equation $R(s, \partial_s) \hat{x}(s) + Q(s) = 0$ yields

$$ R(s, \partial_s) \hat{z}(s) + Q(s) - R(s, \partial_s) \hat{\gamma}(s) = 0. \qquad (13) $$

Setting $l := r + 1$, we have $\hat{\gamma}(s) = c/s^l$. As the constant $c$ is unknown, we shall try to get rid of the last term of (13) by applying a linear differential operator $L(s, \partial_s)$ on the left to Equation (13). This linear differential operator must satisfy $L(s, \partial_s) R(s, \partial_s) \hat{\gamma}(s) = 0$. To compute $L(s, \partial_s)$, we shall first
multiply (13) by a suitable power of s to get polynomial expressions. Let us start by computing $R(s, \partial_s) \hat{\gamma}(s)$. We have

$$R(s, \partial_s) \hat{\gamma}(s) = \sum_{i=0}^{n} a_i(-\partial_s) s^i \frac{1}{s^l} c = c \sum_{i=0}^{n} a_i(-\partial_s) s^{i-l}.$$

If $l \le n$, then $R(s, \partial_s) \hat{\gamma}(s) \in \mathbb{K}[s^{-1}, s]$, where $\mathbb{K}[s^{-1}, s]$ denotes the ring of Laurent polynomials, and $R(s, \partial_s) \hat{\gamma}(s)$ can be decomposed uniquely as the sum of the following two terms

$$\Gamma_{s^{-1}}(s) := c \sum_{i=0}^{l-1} a_i(-\partial_s) \frac{1}{s^{l-i}} \in \mathbb{K}[s^{-1}] \setminus \mathbb{K}, \quad \Gamma_s(s) := c \sum_{i=1}^{n} a_i(-\partial_s) s^{i-l} \in \mathbb{K}[s]. \quad (14)$$

Note that the “constant term” $a_l(-\partial_s)$ 1 in $s^{-1}$ has been chosen to contribute to $\Gamma_s$. But, we could have chosen to add it to $\Gamma_{s^{-1}}$.

Let us consider the “generic case”<fnref n="1" />

$$l \le n, \quad \Gamma_{s^{-1}} \ne 0, \quad \Gamma_s \ne 0. \quad (15)$$

In order to obtain a polynomial, we shall multiply $R(s, \partial_s) \hat{\gamma}(s)$ by a suitable power of s. Let

$$p := \deg_{s^{-1}}(\Gamma_{s^{-1}}) \in \mathbb{Z}_{>0}, \quad p' := \deg_{s}(\Gamma_{s}) \in \mathbb{Z}_{\ge 0}.$$

Note that if we denote $d_i := \deg_t(a_i)$, then we have $p \le \max_{i=0,...,l-1}\{l - i + d_i\}$. Moreover, since $t=0$ is an ordinary point, we have $p' = n-l$.

We thus have $s^p R(s, \partial_s) \hat{\gamma}(s) \in \mathbb{K}[s]$ and multiplying (13) by $s^p$, we get that

$$s^p R(s, \partial_s) \hat{z}(s) + s^p Q(s) - s^p R(s, \partial_s) \hat{\gamma}(s) = 0$$

is an ODE with polynomial coefficients. Setting

$$R'(s, \partial_s) := s^p R(s, \partial_s), \quad S(s) := s^p Q(s) = s^p \sum_{k=0}^{n-1} S_k(s) \vartheta_k, \quad \bar{Q}(s) := \frac{1}{c} s^p R(s, \partial_s) \hat{\gamma}(s),$$

we then get that $\hat{z}$ satisfies

$$R'(s, \partial_s) \hat{z}(s) + S(s) - c \bar{Q}(s) = 0. \quad (16)$$

The last term of (16), which still contains the unknown constant $c$, is now a polynomial in $s$ and its degree $q$ verifies

$$q := \deg_{s}(\bar{Q}) = p + p' = p + n - l \le \max_{i=0,...,l-1} \{l - i + d_i\} + n - l. \quad (17)$$

We can now compute the annihilator of the polynomial $\bar{Q}$ which will allow us to get rid of the term $-c \bar{Q}(s)$ in Equation (16).

**Example 3.1.** Let us consider a sinusoidal signal $x(t) = A \sin(\omega t + \phi)$ satisfying the order $n=2$ ODE $(\frac{d^2}{dt^2} + \omega^2)x(t) = 0$ and corrupted by a perturbation given by a ramp function, i.e., $r=1$ and $\gamma(t) = c t H(t)$, where $c$ is a constant. We thus have $l=2$, $\hat{\gamma}(s) = c/s^2$, and Equation (13) gives

$$(s^2 + \omega^2) \hat{z}(s) - s \vartheta_0 - \vartheta_1 - \frac{(s^2 + \omega^2) c}{s^2} = 0.$$

We are thus in the “generic case” with

$$p = \deg_{s^{-1}}(\Gamma_{s^{-1}}) = \deg_{s^{-1}}\left(\frac{\omega^2}{s^2}\right) = 2, \quad p' = \deg_{s}(\Gamma_{s}) = \deg_{s}(1) = 0.$$

<footnote n="1">¹See the comments after Theorem 3.8 where we discuss the non-generic cases.</footnote>
Equation (16) is thus

$$s^2 (s^2 + \omega^2) \widehat{z}(s) - s^3 \vartheta_0 - s^2 \vartheta_1 - (s^2 + \omega^2) c = 0,$$

and since the constant parameter $c$ is unknown we need to compute the annihilator of the polynomial $\overline{Q}(s) = s^2 + \omega^2$ of degree $q = p + n - l = 2$ to cancel the last term in the left-hand side of the latter equality.

## 3.2 Computation of the annihilator

Let us consider the first Weyl algebra $A := A_1(K') = K'[s] \langle \partial_s | \partial_s s = s \partial_s + 1 \rangle$, where $K'$ is a field of characteristic zero, and let us characterize the annihilator of $P \in K'[s]$, namely, the left ideal of $A$ defined by

$$\mathrm{ann}_A(.P) := \{a \in A \mid a P := a(P) = 0\}.$$

**LEMMA 3.2.** Let $P(s) = \sum_{k=0}^q b_k s^k$, where the $b_k$'s belong to the field $K'$, $q > 0$ and $b_q \neq 0$. Then, we have

$$\mathrm{ann}_A(.P) = \langle -P^{(q)} + P \partial_s^q, -P^{(q)} \partial_s + P^{(1)} \partial_s^q, \dots, -P^{(q)} \partial_s^{q-1} + P^{(q-1)} \partial_s^q, \partial_s^{q+1} \rangle. \quad (18)$$

**PROOF.** The annihilator $\mathrm{ann}_A(.P)$ can be obtained by considering the polynomial relations between the successive derivatives of $P$, i.e., by considering the following $K[s]$-module

$$\ker_{K'[s]}(.L) := \{\lambda := (\lambda_0 \dots \lambda_{q+1}) \in K'[s]^{1 \times (q+2)} \mid \lambda L = 0\}, \quad L := \begin{pmatrix} P \\ P^{(1)} \\ \vdots \\ P^{(q+1)} \end{pmatrix} \in K'[s]^{(q+2) \times 1}.$$

Equivalently, we can compute a generating set of compatibility conditions of the inhomogeneous linear system $L \eta = \zeta$. Let us write $P(s) = \sum_{k=0}^q b_k s^k$, where the $b_k$'s belong to the field $K'$. We first note that

$$P^{(j)}(s) = \sum_{k=j}^{q} \frac{k!}{(k-j)!} b_k s^{k-j} = \sum_{u=0}^{q-j} \frac{(u+j)!}{u!} b_{u+j} s^u.$$

Setting $j = q - v$, we then get $P^{(q-v)}(s) = \sum_{u=0}^v b'_{q-v+u} s^u$, where $b'_{q-v+u} := ((q - v + u)! / u!) b_{q-v+u}$. We thus obtain

$$L \eta = \zeta \iff \left\{ \begin{array}{l} P \eta = \zeta_0, \\ \vdots \\ P^{(q-2)} \eta = (b'_q s^2 + b'_{q-1} s + b'_{q-2}) \eta = \zeta_{q-2}, \\ P^{(q-1)} \eta = (b'_q s + b'_{q-1}) \eta = \zeta_{q-1}, \\ P^{(q)} \eta = b'_q \eta = \zeta_q, \\ P^{(q+1)} \eta = 0 = \zeta_{q+1}. \end{array} \right.$$

Since $b'_q \neq 0$, from the last but one equation, we get $\eta = \zeta_q / b'_q$ and substituting this identity into the rest of the equations, we obtain the following compatibility conditions

$$\frac{P^{(i)}}{P^{(q)}} \zeta_q = \zeta_i \iff P^{(i)} \zeta_q - P^{(q)} \zeta_i = 0, \quad i = 0, \dots, q+1.$$
Hence, we obtain that the rows of the following matrix

$$
M := \begin{pmatrix}
-P^{(q)} & 0 & \dots & 0 & -P & 0 \\
0 & -P^{(q)} & \dots & 0 & -P^{(1)} & 0 \\
0 & 0 & \dots & -P^{(q)} & -P^{(q-1)} & 0 \\
0 & 0 & \dots & \dots & 0 & 1
\end{pmatrix} \in \mathbb{K}'[s]^{(q+1)\times(q+2)}
$$

generates ker$_{K'[s]}(L)$. Then, we have

$$
M \begin{pmatrix} \zeta_0 \\ \vdots \\ \zeta_{q+1} \end{pmatrix} = 0 \iff M \begin{pmatrix} 1 \\ \partial_s \\ \vdots \\ \partial_s^{q+1} \end{pmatrix} P \eta = 0 \iff \begin{pmatrix} -P^{(q)} + P \partial_s^q \\ -P^{(q)} \partial_s + P^{(1)} \partial_s^q \\ \vdots \\ -P^{(q)} \partial_s^{q-1} + P^{(q-1)} \partial_s^q \\ \partial_s^{q+1} \end{pmatrix} P \eta = 0,
$$

which shows (18).

h

For instance, the trivial annihilator $P \partial_s - P^{(1)}$ of $P$ can be expressed as

$$
P \partial_s - P^{(1)} = \frac{1}{P(q)} \left( -P^{(1)} \left( -P^{(q)} + P \partial_s^q \right) + P \left( -P^{(q)} \partial_s + P^{(1)} \partial_s^q \right) \right) \in \mathrm{ann}_A(.P).
$$

According to Lemma 3.2, the left ideal ann<sub>A</sub>(.P) of the Weyl algebra A can thus be generated by
q + 1 elements. However, a famous theorem due to Stafford asserts that every left or right ideal of A
can be generated by two elements [Stafford (1978)] and we can indeed prove the following result.

PROPOSITION 3.3. Let $P(s) = \sum_{k=0}^{q} b_k s^k$, where the $b_k$'s belong to the field $\mathbb{K}', q > 0$ and $b_q \neq 0$. Then, we have

$$
\mathrm{ann}_A(.P) = \left\langle -P^{(q)} + P \partial_s^q, \partial_s^{q+1} \right\rangle .
$$

PROOF. We need to determine two generators for the ideal ann<sub>A</sub>(.P) given by (18). Let us set
G<sub>1</sub> := -P<sup>(q)</sup> + P * ∂<sub>s</sub><sup>q</sup> and G<sub>2</sub> := *∂<sub>s</sub><sup>q+1</sup>*. Then, using P<sup>(q+1)</sup> = 0, we have

$$
\partial_s G_1 - P G_2 = -P^{(q)} \partial_s - P^{(q+1)} + P \partial_s^{q+1} + P^{(1)} \partial_s^q - P \partial_s^{q+1} = -P^{(q)} \partial_s + P^{(1)} \partial_s^q,
$$

which shows that the second generator of ann<sub>A</sub>(.P) given by (18) belongs to the left ideal ⟨G<sub>1</sub>, G<sub>2</sub>⟩
of A generated G<sub>1</sub> and G<sub>2</sub>. By induction, let us suppose that -P<sup>(q)</sup> ∂<sub>s</sub><sup>j</sup> + P<sup>(i)</sup> ∂<sub>s</sub><sup>q</sup> ∈ ⟨G<sub>1</sub>, G<sub>2</sub>⟩ and let us
then prove that -P<sup>(q)</sup> ∂<sub>s</sub><sup>j+1</sup> + P<sup>(i+1)</sup> ∂<sub>s</sub><sup>q</sup> ∈ ⟨G<sub>1</sub>, G<sub>2</sub>⟩. Again, using P<sup>(q+1)</sup> = 0, we then have

$$
\partial_s \left( -P^{(q)} \partial_s^j + P^{(i)} \partial_s^q \right) - P^{(i)} G_2 = -P^{(q)} \partial_s^{j+1} - P^{(q+1)} \partial_s^j + P^{(i)} \partial_s^{q+1} + P^{(i+1)} \partial_s^q - P^{(i)} \partial_s^{q+1}, \\
\qquad = -P^{(q)} \partial_s^{j+1} + P^{(i+1)} \partial_s^q \in \langle G_1, G_2 \rangle,
$$

which proves by induction that all the generators of ann<sub>A</sub>(.P) belong to ⟨G<sub>1</sub>, G<sub>2</sub>⟩. This proves the
result.
h

If $q = 0$, then $P(s) \in \mathbb{K}'$ is a constant and we have $\mathrm{ann}_A(.P) = \langle \partial_s \rangle$. Finally, we note that the first
generator $G_1 = -P^{(q)} + P \partial_s^q$ of $\mathrm{ann}_A(.P)$ depends on the coefficients $b_j$'s of $P$ (and thus of the
coefficients of the polynomials $a_i$'s) contrary to the second one $G_2 = \partial_s^{q+1}$.

*Example 3.4.* In Example 3.1, we found $P(s) = s^2 + \omega^2$. Proposition 3.3 with $\mathbb{K}' = \mathbb{Q}(\omega)$, implies that

$$
\mathrm{ann}_A(.P) = \langle (s^2 + \omega^2) \partial_s^2 - 2, \partial_s^3 \rangle.
$$
### 3.3 Estimation of the initial conditions

We now investigate when the initial conditions $\vartheta_k = x^{(k)}(0)$, $k = 0, \dots, n-1$, can be obtained in terms of convolutions of $z$. Since the perturbation $\gamma(t) = ct^r H(t)$ is not supposed to be known apriori, from (16), we first have to use the elements of the annihilator of $\overline{Q}$ to delete the term $-c\overline{Q}(s)$ containing $c$. Then, we have to study whether or not the $\vartheta_k$'s can be obtained explicitly in terms of derivatives of $\hat{z}$. If so, then using the inverse Laplace transform, we can express the $\vartheta_k$'s as convolutions of the measured signal $z$.

With the notations of Section 3.2, applying the operator $G_1(s, \partial_s) = -\overline{Q}^{(q)} + \overline{Q} \partial_s^q \in \text{ann}_A(\overline{Q})$ to (16), we then get

$$G_1(s, \partial_s) S(s) = -G_1(s, \partial_s) R'(s, \partial_s) \hat{z}(s). \quad (19)$$

Let us now develop the term $G_1(s, \partial_s) S(s)$ to characterize when the initial conditions $\vartheta_k$'s can be expressed in terms of $\hat{z}(s)$. We have

$$G_1(s, \partial_s) S(s, \vartheta) = \sum_{k=0}^{n-1} \underbrace{\left( \left( -\overline{Q}^{(q)} + \overline{Q}(s) \partial_s^q \right) s^p S_k(s) \right)}_{T_k \in \mathbb{K}[s]} \vartheta_k.$$

As in the case without perturbation, a natural idea to express the $\vartheta_k$'s in terms of $\hat{z}$ is to differentiate (19) a certain number of times to obtain a solvable inhomogeneous linear system in the $\vartheta_k$'s. To do that, we first need to study the degrees of the polynomials $T_k$'s.

**LEMMA 3.5.** With the above notations and assumptions, we have

* For all $k \in \{0, \dots, n-1\} \setminus \{r\}$, $\text{deg}_s(T_k) = p + n - k - 1$,

* $\text{deg}_s(T_r) \le p - r - 1 - \min_{i=0, \dots, r} \{v_i - i \mid v_i \ne -\infty\}$.

**PROOF.** We have

$$T_k(s) := \left(-\overline{Q}^{(q)} + \overline{Q}(s) \partial_s^q\right) s^p S_k(s), \quad S_k(s) = \sum_{i=k+1}^{n} a_i (-\partial_s) s^{i-k-1}.$$

Since $t=0$ is an ordinary point of the differential equation (11), we have $S_k(s) = a_n(0)s^{n-k-1} + \dots$, where $\dots$ contains monomials of degree strictly less than $n-k-1$. Moreover, $\overline{Q}(s) = \sum_{k=0}^q b_k s^k$ with $b_q = a_n(0)$ so that $\overline{Q}^{(q)}(s) = q! a_n(0)$, and we further have

$$\partial_s^q s^p S_k(s) = a_n(0) \frac{(p+n-k-1)!}{(p+n-k-1-q)!} s^{p+n-k-1-q} + \dots,$$

where $\dots$ contains monomials of degree strictly less than $p+n-k-1-q$. The leading monomial of $-\overline{Q}^{(q)} s^p S_k(s)$ is thus $-q! a_n(0)^2 s^{p+n-k-1}$ and that of $\overline{Q}(s) \partial_s^q s^p S_k(s)$ is $a_n(0)^2 \frac{(p+n-k-1)!}{(p+n-k-1-q)!} s^{p+n-k-1}$. We then obtain

$$T_k(s) = -\overline{Q}^{(q)} s^p S_k(s) + \overline{Q}(s) \partial_s^q s^p S_k(s), \\ = a_n(0)^2 \left( -q! + \frac{(p+n-k-1)!}{(p+n-k-1-q)!} \right) s^{p+n-k-1} + \dots,$$

where $\dots$ contains monomials of degree strictly less than $p+n-k-1$. We thus have

$$\text{deg}_s(T_k) = p + n - k - 1,$$

except if $\frac{(p+n-k-1)!}{(p+n-k-1-q)!} = q!$, i.e., if $q = p + n - k - 1$ which, from (17), means $k+1=l$, i.e., $k=r$.
It thus remains to determine the degree of $T_r$. Note that, from (14), we have

$$S_r(s) = \sum_{i=r+1}^{n} a_i (-\partial_s) s^{i-r-1} = \sum_{i=l}^{n} a_i (-\partial_s) s^{i-l} = \frac{1}{c} \Gamma_s(s),$$

and as $R(s, \partial_s) \hat{\Upsilon}(s) = \Gamma_{s-1}(s) + \Gamma_s(s)$, we get

$$S_r(s) = \frac{1}{c} (R(s, \partial_s) \hat{\Upsilon}(s) - \Gamma_{s-1}(s)).$$

Consequently, we have

$$T_r(s) = (-\bar{Q}^{(q)} + \bar{Q}(s) \partial_s^q) s^p S_r(s) = -\frac{1}{c} (-\bar{Q}^{(q)} + \bar{Q}(s) \partial_s^q) s^p \Gamma_{s-1}(s),$$

since $(-\bar{Q}^{(q)} + \bar{Q}(s) \partial_s^q) s^p R(s, \partial_s) \hat{\Upsilon}(s) = (-\bar{Q}^{(q)} + \bar{Q}(s) \partial_s^q) \bar{Q}(s) = 0$, and, from (14), we thus have

$$T_r(s) = - \left( -\bar{Q}^{(q)} + \bar{Q}(s) \partial_s^q \right) s^p \sum_{i=0}^{r} a_i (-\partial_s) \frac{1}{s^{r+1-i}}.$$

Denoting $v_i := v_0(a_i)$ the valuation of $a_i$, we thus obtain

$$\deg_s(T_r) \le p - \min_{i=0,...,r} \{r + 1 - i + v_i \mid v_i \ne -\infty\}.$$

h

We now state the following corollary of Lemma 3.5 which implies that (19) provides a triangular inhomogeneous linear system for the $\vartheta_k$'s.

**COROLLARY 3.6.** With the above notations and assumptions, we have

$$\forall k, k' \in \{0, \dots, n-1\}, k \ne k', \deg_s(T_k) \ne \deg_s(T_{k'}).$$

**PROOF.** From Lemma 3.5, it suffices to prove that $\deg_s(T_r) \ne \deg_s(T_j)$ for all $j \in \{0, \dots, n-1\} \setminus \{r\}$. If $r < n-1$, then proving that $\deg_s(T_r) < \deg_s(T_{n-1}) = p$ is enough. By contradiction, let us assume $p \le \deg_s(T_r) \le p - r - 1 - \min_{i=0,...,r} \{v_i - i \mid v_i \ne -\infty\}$. Then, there exists $j \in \{0, \dots, r\}$ such that $p \le p - r - 1 - v_j + j$ and $v_j \ne -\infty$. This yields $v_j \le j - (r+1) < 0$ which is in contradiction with $v_j \ge 0$. Similarly, in the case $r = n-1$, we prove that $\deg_s(T_r = T_{n-1}) < \deg_s(T_{n-2}) = p + 1$. This ends the proof. h

**Example 3.7.** Let us continue Example 3.1. We have $S(s, \vartheta) = -s^3\vartheta_0 - s^2\vartheta_1$ and in Example 3.4, we found $G_1(s, \partial_s) = (s^2 + \omega^2) \partial_s^2 - 2$. This implies

$$G_1(s, \partial_s) S(s, \vartheta) = -2s(2s^2 + 3\omega^2) \vartheta_0 - 2\omega^2 \vartheta_1,$$

so that $T_0(s) = -2s(2s^2 + 3\omega^2)$ and $T_1(s) = -2\omega^2$. We recall that in this example $n = p = 2$, $r = 1$ and we can check that

$$\deg_s(T_0) = 2 + 2 - 0 - 1 = 3 \ne \deg_s(T_1) = 0 \le 2 - 1 - 1 - \min_{i=0,1} \{v_i - i \mid v_i \ne -\infty\} = 0,$$

since $v_0 = 0$ and $v_1 = -\infty$.

Let us finally explain how Corollary 3.6 permits us to compute all the $\vartheta_k$'s. We start by computing $\vartheta_j$, where $j$ is such that $f_j := \deg_s(T_j) > \deg_s(T_k)$ for $k \ne j$. Indeed differentiating $f_j$ times (19), all the $\vartheta_k$'s, for $k \ne j$, disappear and we get a linear equation for $\vartheta_j$. More precisely, if we note $T_j(s) = t_{jf_j} s^{f_j} + \dots$, where $\dots$ contains monomials of degree strictly less than $f_j$, then applying $\partial_s^{f_j}$ to (19), we get

$$\vartheta_j = - \frac{1}{f_j! t_{j f_j}} \partial_s^{f_j} G_1(s, \partial_s) R'(s, \partial_s) \hat{z}(s).$$
Note that, if $j \neq r$, then $t_{jf_j} = a_n(0)^2 \left(-q! + \frac{(p+n-j-1)!}{(p+n-j-1-q)!}\right) \neq 0$. Hence, for $j \neq r$, $\vartheta_j$ depends polynomially on the coefficients $a_{ij}$ of the polynomials $a_i$ defining (11) and on $a_n(0)^{-1}$, i.e., $a_{n0}^{-1}$. Knowing $\vartheta_j$, we can then proceed and compute the next $\vartheta_k$'s using the same technique.

**THEOREM 3.8.** With the above notations, if $t = 0$ is a regular point for the ODE (11) and if we are in the “generic case”, i.e., $l \le n$, $\Gamma_{s-1} \ne 0$, and $\Gamma_s \ne 0$, then all the initial conditions $\{\vartheta_k = x^{(k)}(0)\}_{k=0,...,n-1}$ of (11) can be expressed explicitly in terms of an operator $Q_k(s, \partial_s)$ of $A_1(\mathbb{K})$ applied to $\tilde{z}(s)$, i.e., $\vartheta_k = Q_k(s, \partial_s)\tilde{z}(s)$, for $k = 0, ..., n-1$. Moreover, if $l_k := \deg_s Q_k(s, \partial_s)$, then writing $\vartheta_k = n_k/d_k$, where

$$n_k = \frac{Q_k(s, \partial_s)}{s^{l_k+1}} \tilde{z}(s), \quad d_k = \frac{1}{s^{l_k+1}}, \quad k = 0, \dots, n-1,$$

we then have

$$\vartheta_k = \frac{\mathcal{L}^{-1}(n_k)}{\mathcal{L}^{-1}(d_k)}, \quad \mathcal{L}^{-1}(d_k) = \frac{t^{l_k}}{l_k!}.$$

Hence, all the initial conditions can be explicitly expressed as convolutions of $z$ depending rationally on the coefficients of the polynomial coefficients $a_i$'s of the ODE (11).

Let us make some comments on the assumptions of Theorem 3.8.

(1) If $t = 0$ is not an ordinary point of the ODE (11), then, as in the case without perturbation considered in the previous section, we cannot recover all the initial conditions. For instance, $\vartheta_{n-1} = x^{(n-1)}(0)$ cannot be obtained as it does not appear in the equation after applying the Laplace transform. Indeed, with the above notations, $S_{n-1}(s) = -a_n(-\partial_s)1 = 0$ if $a_n(0) = 0$. For more details, see Example 2.2.

(2) If $\Gamma_{s-1} = 0$, then $T_r(s) = 0$ and we cannot recover the initial condition $\vartheta_r = x^{(r)}(0)$ as it does not appear in the equation after applying the operator $G_1(s, \partial_s)$ to get rid of the unknown constant $c$. All the other initial conditions can be obtained as explained above.

(3) If $l > n$, then we have $\Gamma_s = 0$, and the strategy developed in this section also permits us to obtain all the initial conditions. For more details, see [Chartouny (2021)].

(4) If $\overline{Q}(s) \in \mathbb{K}$ is a constant, i.e., $q = 0$, then we have $\text{ann}_A(\overline{Q}) = \langle \partial_s \rangle$ and the above strategy, where we choose $G_1(s, \partial_s) = -\overline{Q}^{(q)} + \overline{Q} \partial_s^q$ instead of $\partial_s$, allows us to get the initial conditions in a similar manner.

**Example 3.9.** Continuing Example 3.1, we start by computing $\vartheta_0$. Here Equation (19) writes

$$-2 \left(2s^3 + 3s\omega^2\right) \vartheta_0 - 2\omega^2 \vartheta_1 = - \left(s^2(s^2 + \omega^2)\right)^2 \frac{d^2}{ds^2} \tilde{z}(s) + \left(8s^5 + 12\omega^2 s^3 + 4\omega^4 s\right) \frac{d}{ds} \tilde{z}(s) + \left(10s^4 + 12\omega^2 s^2 + 2\omega^4\right) \tilde{z}(s). \quad (20)$$

Differentiating 3 times with respect to $s$ and solving for $\vartheta_0$, we get

$$\vartheta_0 = \frac{1}{24} \omega^4 s^2 \frac{d^5}{ds^5} \tilde{z}(s) + \frac{1}{12} \omega^2 s^4 \frac{d^5}{ds^5} \tilde{z}(s) + \frac{1}{24} s^6 \frac{d^5}{ds^5} \tilde{z}(s) + \frac{5}{12} \omega^4 s \frac{d^4}{ds^4} \tilde{z}(s) + \frac{3}{2} \omega^2 s^3 \frac{d^4}{ds^4} \tilde{z}(s) + \frac{13}{12} s^5 \frac{d^4}{ds^4} \tilde{z}(s) \\ + \frac{5}{6} \omega^4 \frac{d^3}{ds^3} \tilde{z}(s) + 8 \omega^2 s^2 \frac{d^3}{ds^3} \tilde{z}(s) + \frac{55}{6} s^4 \frac{d^3}{ds^3} \tilde{z}(s) + 14 \omega^2 s \frac{d^2}{ds^2} \tilde{z}(s) + 30 s^3 \frac{d^2}{ds^2} \tilde{z}(s) + 6 \omega^2 \frac{d}{ds} \tilde{z}(s) + 35 s^2 \frac{d}{ds} \tilde{z}(s) + 10 s \tilde{z}(s).$$

Then, applying the inverse Laplace transform, after dividing by $s^7$ and introducing the denominator $1/s^7$ as explained in Theorem 3.8, we get

$$\begin{aligned} \vartheta_0 = & 30 \frac{1}{t^6} \left(-\int_0^t \frac{1}{24} \omega^4 (t-\tau_1)^4 \tau_1^5 z(\tau_1) d\tau_1 - \int_0^t \omega^2 (t-\tau_1)^2 \tau_1^5 z(\tau_1) d\tau_1 - \int_0^t \tau_1^5 z(\tau_1) d\tau_1 + \int_0^t \frac{1}{12} \omega^4 (t-\tau_1)^5 \tau_1^4 z(\tau_1) d\tau_1 \right. \\ & + \int_0^t 6 \omega^2 (t-\tau_1)^3 \tau_1^4 z(\tau_1) d\tau_1 + \int_0^t (26t-26\tau_1) \tau_1^4 z(\tau_1) d\tau_1 - \int_0^t \frac{1}{36} \omega^4 (t-\tau_1)^6 \tau_1^3 z(\tau_1) d\tau_1 - \int_0^t 8 \omega^2 (t-\tau_1)^4 \tau_1^3 z(\tau_1) d\tau_1 \\ & - \int_0^t 110 (t-\tau_1)^2 \tau_1^3 z(\tau_1) d\tau_1 + \int_0^t \frac{14}{5} \omega^2 (t-\tau_1)^5 \tau_1^2 z(\tau_1) d\tau_1 + \int_0^t 120 (t-\tau_1)^3 \tau_1^2 z(\tau_1) d\tau_1 - \int_0^t \frac{1}{5} \omega^2 (t-\tau_1)^6 \tau_1 z(\tau_1) d\tau_1 \\ & \left. - \int_0^t 35 (t-\tau_1)^4 \tau_1 z(\tau_1) d\tau_1 + \int_0^t 2 (t-\tau_1)^5 z(\tau_1) d\tau_1 \right). \end{aligned}$$

Now, we can proceed similarly to get $\vartheta_1$ from (20). No differentiation is needed since $\deg_s(T_1) = 0$ and we can replace $\vartheta_0$ by its expression obtained above.
## 3.4 Estimation of the coefficients of the ODE

Theorem 3.8 above provides explicit expressions for the initial conditions {$\theta_k = x^{(k)}(0)$}$_{k=0,...,n-1}$ of (11) in terms of the measured function $z$ and its repeated indefinite integrals (e.g., convolutions). Let the polynomial coefficients of (11) be denoted by $a_i(t) = \sum_{j=v_i}^{d_i} a_{ij} t^j$, for $i = 0, ..., n$, where the $a_{ij}$'s belong to $\mathbb{K}$. The expressions of the $\theta_k$'s can depend on the equation parameters, namely, on the $a_{ij}$'s, some of which could be unknown or unfixed to numerical values. Hence, an important issue studied here is to investigate when the coefficients $a_{ij}$'s of (11) can be found again in terms of convolutions of $z$ which are independent of the value $c$ of the perturbation $\gamma(t)$ (supposed to be unknown).

As noticed at the end of Section 3.2, contrary to the first generator $G_1 = -\bar{Q}^{(q)} + \bar{Q}(s) \partial_s^q$ used in Section 3.3, the second generator $G_2 = \partial_s^{q+1}$ of ann$_A(\bar{Q})$ does not depend on the coefficients $a_{ij}$'s. Hence, applying the operator $G_2$ to (16), i.e., to $s^p R(s, \partial_s) \hat{z}(s) + S(s) - c \bar{Q}(s) = 0$, we then obtain

$$ \partial_s^{q+1} s^p R(s, \partial_s) \hat{z}(s) + \partial_s^{q+1} S(s) = 0, \quad (21) $$

which is a linear expression in the $a_{ij}$'s.

To obtain linear equations for the $a_{ij}$'s, we first remove the $\theta_k$'s from Equation (21). We recall that $S(s) = s^p \sum_{k=0}^{n-1} S_k(s) \theta_k$, where $S_k \in \mathbb{K}[s]$ and $\deg_S(S) = n - k - 1$ so that $\deg_S(S) = p + n - 1$ since we assume that $t = 0$ is an ordinary point of (11). Then, if $q+1 > p+n-1$ which, using (17), yields $l-1 = r < 1$, i.e., $r=0$ and the perturbation is a bias, then we have $\partial_s^{q+1} S(s) = 0$ and Equation (21) reduces to $\partial_s^{q+1} s^p R(s, \partial_s) \hat{z}(s) = 0$. Otherwise, for $r > 0$, we still need to apply $\partial^r$ to (21) to remove $S(s)$ and we get $\partial_s^{r+q+1} s^p R(s, \partial_s) \hat{z}(s) = 0$. Hence, we always obtain

$$ \partial_s^{r+q+1} s^p R(s, \partial_s) \hat{z}(s) = 0. \quad (22) $$

Equation (22) depends linearly on the coefficients $a_{ij}$'s of (11) and on the parameter $r$ defining the perturbation $\gamma(t) = c t^r H(t)$. In practice, it usually makes sense to consider that $r$ is apriori known.

Using $R(s, \partial_s) = \sum_{i=0}^n a_i(-\partial_s) s^i$, where $a_i(-\partial_s) = \sum_{j=v_i}^{d_i} a_{ij}(-\partial_s)^j$, Equation (22) becomes

$$ \sum_{i=0}^{n} \sum_{j=v_i}^{d_i} \left( \partial_s^{r+q+1} s^p (-\partial_s)^j s^i \right) a_{ij} \hat{z}(s) = 0. \quad (23) $$

Hence, if set $g := \sum_{i=0}^n (d_i - v_i + 1) \ge 0$ and

$$ \vec{a} := (a_{0v_0} \quad \dots \quad a_{0d_0} \quad a_{1v_1} \quad \dots \quad a_{1d_1} \quad \dots \quad a_{nv_n} \quad \dots \quad a_{nd_n})^T \in \mathbb{K}^g, $$

then the $g$ coefficients $a_{ij}$'s of (11) satisfy the following linear equation

$$ \left( \partial_s^{r+q+1} s^p (-\partial_s)^{v_0} \quad \dots \quad \partial_s^{r+q+1} s^p (-\partial_s)^{d_n} s^n \right) \vec{a} \hat{z}(s) = 0. $$

Note that (23) can be differentiated, e.g., $w$ times. Hence, the $a_{ij}$'s also satisfy the linear system

$$ \left( \begin{array}{cccc} \partial_s^{r+q+1} s^p (-\partial_s)^{v_0} & \dots & \partial_s^{r+q+1} s^p (-\partial_s)^{d_n} s^n \\ \vdots & \ddots & \vdots \\ \partial_s^{r+q+1+w} s^p (-\partial_s)^{v_0} & \dots & \partial_s^{r+q+1+w} s^p (-\partial_s)^{d_n} s^n \end{array} \right) \vec{a} \hat{z}(s) = 0. \quad (24) $$

We can then repeat what we did in the case without perturbation to obtain $\vec{a}$, i.e., the $a_{ij}$'s in terms of $\hat{z}$ and its derivatives by solving the linear system of equations (24). Finally, applying the inverse Laplace transform as before, we get the following result.
**THEOREM 3.10.** With the above notations, if $t = 0$ is a regular point for the ODE (11), then the coefficients $a_{ij}$'s of (11) can be expressed explicitly as convolutions of the measured function $z$.

**Example 3.11.** Let us continue Example 3.1. The expressions of $\vartheta_0$ and $\vartheta_1$ obtained in Example 3.9 still depend on the parameter $\omega^2$ appearing as a coefficient of the ODE satisfied by the sinusoidal signal. We can then use Theorem 3.10 to estimate $\omega^2$. Since $r = 1$ and $q = 2$, Equation (22) writes

$$24 \tilde{z}(s) + 96 s \frac{d}{ds} \tilde{z}(s) + 12 (s^2 + \omega^2) \frac{d^2}{ds^2} \tilde{z}(s) + 60 s \frac{d^2}{ds^2} \tilde{z}(s) + 8 s (s^2 + \omega^2) \frac{d^3}{ds^3} \tilde{z}(s) + 8 s^3 \frac{d^3}{ds^3} \tilde{z}(s) + s^2 (s^2 + \omega^2) \frac{d^4}{ds^4} \tilde{z}(s) = 0.$$

From the latter equality, we get

$$\omega^2 = - \frac{s^4 \frac{d^4}{ds^4} \tilde{z}(s) + 16 s^3 \frac{d^3}{ds^3} \tilde{z}(s) + 72 s^2 \frac{d^2}{ds^2} \tilde{z}(s) + 96 s \frac{d}{ds} \tilde{z}(s) + 24 \tilde{z}(s)}{s^2 \frac{d^4}{ds^4} \tilde{z}(s) + 8 s \frac{d^3}{ds^3} \tilde{z}(s) + 12 \frac{d^2}{ds^2} \tilde{z}(s)},$$

and applying the inverse Laplace transform, we finally obtain

$$\omega^2 = \frac{-\int_0^t \tau_1^4 z(\tau_1) d\tau_1 - \int_0^t (-16t + 16\tau_1) \tau_1^3 z(\tau_1) d\tau_1 - \int_0^t 36(t - \tau_1)^2 \tau_1^2 z(\tau_1) d\tau_1 + \int_0^t 16(t - \tau_1)^3 \tau_1 z(\tau_1) d\tau_1 - \int_0^t (t - \tau_1)^4 z(\tau_1) d\tau_1}{\int_0^t \frac{1}{2}(t - \tau_1)^2 \tau_1^4 z(\tau_1) d\tau_1 - \int_0^t \frac{4}{3}(t - \tau_1)^3 \tau_1^3 z(\tau_1) d\tau_1 + \int_0^t \frac{1}{2}(t - \tau_1)^4 \tau_1^2 z(\tau_1) d\tau_1}.$$

Finally, combining Theorems 3.8 and 3.10, the initial conditions $\{x^{(k)}(0)\}_{k=0,...,n-1}$ of the ODE (11) can explicitly be expressed as rational functions of convolutions of $x$.

# 4 Implementation in Maple

Maple is a mathematics-based software and services for education, engineering, and research. For instance, it can manipulate mathematical expressions and find symbolic solutions to certain problems, such as those arising from ordinary and partial differential equations.

## 4.1 The NONA package

The NONA package<fnref n="2" /> is dedicated to the effective study of the algebraic parameter estimation problem, introduced by M. Fliess and H. Sira-Ramirez in [Fliess et al. (2003)] and investigated in the former Non-A project-team (Inria Lille-Nord Europe). This Maple package, based on the OREMODULES package [Chyzak et al. (2007)], has been developed by A. Quadrat (Inria Paris, Ouragan project-team). The binary is freely downloadable at https://who.rocq.inria.fr/Alban.Quadrat/Non-A/. Let us shortly describe the main commands of this package.

<table><thead><tr><th>Annihilator</th><th>Compute the annihilator of a given polynomial</th></tr></thead><tbody><tr><td>AnnihilatorOfExpansion</td><td>Compute the annihilator of a finite linear combination of signals defined by ODEs with polynomial coefficients</td></tr><tr><td>ParameterEstimation</td><td>State P, Q and Q&#773; defined in Sections 2 and 3. The option “bias” treats the particular case of a bias perturbation</td></tr><tr><td>ParameterEstimationEq</td><td>State P z(s) + Q + Q&#773; y(s). The option “bias” treats the particular case of a bias perturbation</td></tr></tbody></table>

## 4.2 Development of the NONA package

For some examples, Maple 2020 is not able to calculate the inverse Laplace transform. For instance, if we want to compute the inverse Laplace transform of $Z'(s)^2$ using Maple, i.e., the square of the first derivative of $Z(s)$ with respect to $s$, the output is “Error, (in collect) invalid 1st argument normal”. We encounter the same problem with expressions of the form $Z(s)^2$, $Z(s)^3$, $Z''(s)^4$, $Z(s)Z'(s)^3$. To solve this problem, we have implemented a convolution procedure. This

<footnote n="2">²NONA stands for “non-asymptotic estimation methods”.</footnote>
function takes as input an expression e that can contain powers of functions in s and the parameter t.
The output of the function is the convolution of e. For instance, the command

convolution(Z(s)^2, t);

yields the following output

$$
\int_{0}^{t} Z(t - \tau_{1}) Z(\tau_{1}) d\tau_{1}.
$$

Moreover, if we want to compute the degree in s of an expression that has the following form
e := z(s)^5 s^9 + s^4 + z''(s)s^2, the function degree will display FAIL. This result is expected because
z(s) is a function in s. Therefore, the degree depends on the expression of z. However, in our study,
to calculate the Laplace inverse (as seen in Sections 2 and 3), we need to know the degree in s of
these types of expressions regardless of z(s). Therefore, we implemented a degree function deg to
compute the degree in s of these expressions, i.e., expressions that can contain functions in s as
well as their derivatives. For instance, if we consider the above expression again,

e := z(s)^5*s^9 + s^4 + diff(z,s,s)*s^2;
deg(e,s);

yields 9.

Furthermore, we have implemented a function ExplicitExpression that computes the explicit expression of expression given in the frequency domain.

For instance, if we consider

v_0 := 2*s*diff(x,s,s) + (s^2 + 2*(k + 2)*s)*diff(x,s) + 2*s*x(s);
ExplicitExpression(v_0,s);

we get as output

$$
\frac{2}{t^2} \left( - \int_0^t \tau x(\tau) d\tau + 2 \int_0^t (t-\tau) \tau^2 x(\tau) d\tau + 2 \int_0^t (t-\tau) x(\tau) d\tau - (2+k) \int_0^t (t-\tau)^2 \tau x(\tau) d\tau \right).
$$

We have also implemented a function InvMatrix that directly computes the inverse of the matrix
U defined in Section 2 without using the Maple command MatrixInverse. As input, the function
takes the matrix U and gives U⁻¹ as output. This function enables us to save time in the computation
process since the explicit expression of U⁻¹ can be given (see [Chartouny et al. (2021)]). For instance,
when using the regular function of Maple to compute the inverse of U for n = 9, the function takes
0.183 CPU seconds. On the other hand, the InvMatrix only takes 0.005 CPU seconds. Of course,
the difference highly increases with the size n of the matrix U to inverse.

More importantly, we are collecting different procedures that make automatic the computation
of the initial conditions and the parameters of a signal defined by a linear ODE with polynomial
coefficients affected by a structured perturbation of the form given in Section 3.

Finally, our second contribution to the NonA package is the development of more explicit examples
that have been added to the NonA library of examples (e.g., sums of exponentials, sums of sinusoids,
sums of orthogonal polynomials).

4.3 Worked example

We illustrate the NonA package with Examples 3.1, 3.4, 3.7, 3.9, and 3.11. For more complex and
therefore longer examples (e.g., sums of exponentials, sums of sinusoids, sums of orthogonal
polynomials), see the library of examples of the NonA package.

We first start by loading the following two packages:

> with(OreModules): with(NonA):
We then define the OD operator annihilating the signal $x(t) = a \sin(\omega t + \phi)$ and its order:

> L := dt^2+omega^2; n := 2;

$$L := dt^2 + \omega^2$$

We can then define the following two Weyl algebras

$$A = \mathbb{Q}(\omega, \vartheta_0, \vartheta_1, c)\langle\partial_t, t | \partial_t, t = t \partial_t + 1\rangle, \quad B = \mathbb{Q}(\omega, \vartheta_0, \vartheta_1, c)\langle\partial_s, s | \partial_s, s = s \partial_s + 1\rangle$$

in which we shall make the computations (for the time-domain and for the frequency-domain):

> A := DefineOreAlgebra(diff=[ds,s],polynom=[s],comm=[omega, vartheta[0],

> vartheta[1],c]):

> B := DefineOreAlgebra(diff=[dt,t],polynom=[t],comm=[omega, vartheta[0],

> vartheta[1],c]):

We can now apply the command ParameterEstimation of the NONA package that generates the OD operators and polynomials defining the different terms of Equation (13):

> ParEst := ParameterEstimation(L,B);

$$ParEst := [\omega^2 + s^2, [-s, -1], -\omega^2 - s^2]$$

We state again that we consider the case of a ramp perturbation, i.e., $r = 1$, so that $\hat{\gamma}(s) = c/s^2$. Simply denoting the Laplace transform of $z(t)$ by $z(s)$ (instead of $\hat{z}(s)$ as done in the text), Equation (13) is given by the equation eq defined below:

> P := ParEst[1];

> Q := add(ParEst[2][i]*vartheta[i-1],i=1..n):

> eq := ApplyMatrix(P,z(s),A)+Q-ApplyMatrix(P,c/s^2,A);

$$eq := (\omega^2 + s^2) z(s) - s\vartheta_0 - \vartheta_1 - \frac{(\omega^2 + s^2)c}{s^2}$$

We can define the polynomial $\bar{Q}$:

> Qbar := 1/c*s^2*ApplyMatrix(P,c/s^2,A);

$$Qbar := \omega^2 + s^2$$

and compute its annihilator using the command Annihilator of the NONA package:

> AnnQbar := Annihilator(Qbar,A);

$$AnnQbar := \begin{bmatrix} \omega^2 ds^2 + sds - 2 \\ sds^2 - ds \\ ds^3 \end{bmatrix}$$

We can then define the element $G_1$ used in Example 3.7 and check again that it belongs to $ann_A(\bar{Q})$:

> P1 := simplify((-diff(Qbar,s$2)+Qbar*ds^2));

$$P1 := -2 + (\omega^2 + s^2) ds^2$$

> Factorize(Matrix([P1]),AnnQbar,A);

$$[ 1 \quad s \quad 0 ]$$

We obtain that $G_1$ can be expressed as the first entry of the matrix $AnnQbar$ plus $s$ times the second entry of $AnnQbar$. Moreover, we can also check again that $P1 \bar{Q} = 0$:

> ApplyMatrix(P1,Qbar,A);
Next, we apply $G_1$ to $s^2$ eq to get the equation from which we can compute the initial conditions $\vartheta_0$ and $\vartheta_1$. Note that, as expected, this equation does not depend on the unknown constant $c$.

> EQ1 := ApplyMatrix(P1, s^2*eq, A);

$$EQ1 := s^2 (\omega^2 + s^2)^2 \frac{d^2}{ds^2} z(s) + (4 \omega^4 s + 12 \omega^2 s^3 + 8 s^5) \frac{d}{ds} z(s) + (2 \omega^4 + 12 \omega^2 s^2 + 10 s^4) z(s) \\ -4 s^3 \vartheta_0 - 6 \omega^2 s \vartheta_0 - 2 \omega^2 \vartheta_1$$

Differentiating 3 times (see Examples 3.7 and 3.9 for more details) this equation and solving for $\vartheta_0$ we then get

> vartheta0_s := normal(solve(diff(EQ1, s$3$), vartheta[0]));

$$\text{vartheta0}_s := 1/24 \left( \frac{d^5}{ds^5} z(s) \right) \omega^4 s^2 + 1/12 \left( \frac{d^5}{ds^5} z(s) \right) \omega^2 s^4 + 1/24 \left( \frac{d^5}{ds^5} z(s) \right) s^6 + \frac{5 \left( \frac{d^4}{ds^4} z(s) \right) \omega^4 s}{12} \\ + 3/2 \left( \frac{d^4}{ds^4} z(s) \right) \omega^2 s^3 + \frac{13}{12} \left( \frac{d^4}{ds^4} z(s) \right) s^5 + 5/6 \left( \frac{d^3}{ds^3} z(s) \right) \omega^4 + 8 \left( \frac{d^3}{ds^3} z(s) \right) \omega^2 s^2 \\ + \frac{55 s^4}{6} \frac{d^3}{ds^3} z(s) + 14 \left( \frac{d^2}{ds^2} z(s) \right) \omega^2 s + 30 s^3 \frac{d^2}{ds^2} z(s) + 6 \left( \frac{d}{ds} z(s) \right) \omega^2 + 35 \left( \frac{d}{ds} z(s) \right) s^2 + 10 s z(s)$$

We then apply the command `ExplicitExpression` of the `NonA` package to obtain the following expression of $\vartheta_0$ in the time-domain:

> vartheta0_t := ExplicitExpression(vartheta0_s, s);

$$\text{vartheta0}_t = 30 \frac{1}{6} \left( \int_0^t \frac{1}{1-24 \omega^4 (t-\tau_1)^4 \tau_1^5 z(\tau_1)}{\tau_1 + \int_0^t -\omega^2 (t-\tau_1)^2 \tau_1^5 z(\tau_1)} d\tau_1 + \int_0^t -\tau_1^5 z(\tau_1) d\tau_1 \right. \\ + \int_0^t \frac{1}{12} \omega^4 (t-\tau_1)^5 \tau_1^4 z(\tau_1) d\tau_1 + \int_0^t \frac{6}{26} \omega^2 (t-\tau_1)^3 \tau_1^4 z(\tau_1) d\tau_1 + \int_0^t \frac{6}{26} \omega^2 (t-\tau_1)^4 \tau_1^4 z(\tau_1) d\tau_1 \\ + \int_0^t \frac{1}{1-36} \omega^4 (t-\tau_1)^6 \tau_1^3 z(\tau_1) d\tau_1 + \int_0^t \frac{8}{110} \omega^2 (t-\tau_1)^4 \tau_1^3 z(\tau_1) d\tau_1 + \int_0^t \frac{14}{5} \omega^2 (t-\tau_1)^5 \tau_1^2 z(\tau_1) d\tau_1 \\ \left. + \int_0^t \frac{120}{1-1/5} \omega^2 (t-\tau_1)^6 \tau_1^2 z(\tau_1) d\tau_1 + \int_0^t \frac{1}{35} \omega^2 (t-\tau_1)^4 \tau_1^2 z(\tau_1) d\tau_1 + \int_0^t \frac{1}{2} \omega^2 (t-\tau_1)^5 z(\tau_1) d\tau_1 \right)$$

Using this formula for $\vartheta_0$, we can compute an expression for $\vartheta_1$. We do not print the result here.

> vartheta1_s := normal(subs(vartheta[0]=vartheta0_s, solve(EQ1, vartheta[1]))):

> vartheta1_t := ExplicitExpression(vartheta1_s, s);

The above expressions involved the unknown coefficient $\omega^2$ of the ODE satisfied by $x(t)$. We shall now estimate it. To achieve this task, we first apply $\partial_s^4$ to the equation $s^2$ eq (see Example 3.11):

> P3 := ds^4;

> EQ2 := ApplyMatrix(P3, s^2*eq, A);

$$EQ2 := \left( \frac{d^4}{ds^4} z(s) \right) \omega^2 s^2 + \left( \frac{d^4}{ds^4} z(s) \right) s^4 + 8 \left( \frac{d^3}{ds^3} z(s) \right) \omega^2 s + 16 \left( \frac{d^3}{ds^3} z(s) \right) s^3 \\ + 12 \left( \frac{d^2}{ds^2} z(s) \right) \omega^2 + 72 \left( \frac{d^2}{ds^2} z(s) \right) s^2 + 96 s \frac{d}{ds} z(s) + 24 z(s)$$

and we then solve the result for $\omega^2$ to obtain

> omega2_s := solve(EQ2, omega^2);

$$\text{omega2}_s := - \frac{\left( \frac{d^4}{ds^4} z(s) \right) s^4 + 16 \left( \frac{d^3}{ds^3} z(s) \right) s^3 + 72 \left( \frac{d^2}{ds^2} z(s) \right) s^2 + 96 s \frac{d}{ds} z(s) + 24 z(s)}{\left( \frac{d^4}{ds^4} z(s) \right) s^2 + 8 s \frac{d^3}{ds^3} z(s) + 12 \frac{d^2}{ds^2} z(s)}$$

Then, we use the `ExplicitExpression` command of the `NonA` package again to obtain the following expression of $\omega^2$ in the time domain:

> omega2_t := ExplicitExpression(omega2_s, s);

omega2_t :=

$$\frac{\int_0^t -\tau_1^4 z(\tau_1) d\tau_1 + \int_0^t -(-16 t + 16 \tau_1) \tau_1^3 z(\tau_1) d\tau_1 + \int_0^t -36 (t-\tau_1)^2 \tau_1^2 z(\tau_1) d\tau_1 + \int_0^t 16 (t-\tau_1)^3 \tau_1 z(\tau_1) d\tau_1 + \int_0^t -(t-\tau_1)^4 \tau_1 z(\tau_1) d\tau_1}{\int_0^t 1/2 (t-\tau_1)^2 \tau_1^4 z(\tau_1) d\tau_1 + \int_0^t -4/3 (t-\tau_1)^3 \tau_1^3 z(\tau_1) d\tau_1 + \int_0^t 1/2 (t-\tau_1)^4 \tau_1^2 z(\tau_1) d\tau_1}$$
Finally, if we replace $z(t)$ by $x(t) = a \sin(\omega t + \phi)$ in the formulas obtained for $\omega^2$, $\vartheta_0$, and $\vartheta_1$, then we can check that we find again $\omega^2$, $a \sin(\phi)$, and $\omega a \cos(\phi)$ as expected:

> z := t -> a*\sin(omega*t+phi):

> simplify convert(omega2_t, int);

> simplify convert(vartheta0_t, int);

> simplify convert(vartheta1_t, int);

$\omega^2$

$a \sin(\phi)$

$\omega a \cos(\phi)$

## 5 Conclusion

In this paper, we have studied mathematical and computer algebra aspects of the algebraic estimation problem initiated in [Fliess et al. (2003)] and further developed in [Belkoura et al. (2009), Mboup (2009), Quadrat (2017), Ushirobira et al. (2016)] (see also the references therein). In particular, we have shown how the results obtained in [Chartouny et al. (2021)] on the inverse Cauchy problem for linear ODEs with polynomial coefficients can be extended to handle the case of a signal corrupted by a standard structured perturbation of the form $\gamma(t) = c t^r H(t)$, where $r \in \mathbb{Z}_{\ge 0}$. To our knowledge, the results obtained in this paper are the first general ones that characterize the possibility to estimate the initial conditions and the constant parameters of a signal defined by a linear ODE with polynomial coefficients. These results are implemented in the NOnA package, developed in Maple, dedicated to the effective study of the algebraic estimation problem.

In the future, we shall incorporate constant delays in the class of structured perturbations to cover the main perturbations classically considered in practice. To do that, we shall consider methods and results developed in [Belkoura et al. (2009)]. Moreover, using the general results obtained in this paper based on closed-form solutions, we can now analyze different numerical aspects of the algebraic estimation method developed in [Fliess et al. (2003)] and also study the effects of noise that corrupts the measurement. See [Fliess et al. (2003), Mboup (2009)] for a precise study of these important problems in the case of particular classes of signals.

Finally, based on the recent effective study of rings of ordinary integro-differential operators (see [Quadrat et al. (2020)] and the references therein) and its implementation in Maple, we want to develop a purely time-domain approach to the algebraic estimation problem. Within this new approach, the closed-form solutions for the initial conditions and the constant parameters defining the signal would be obtained directly in the time domain, thus avoiding the transformations forth and back from the time domain to the frequency domain through the (inverse) Laplace transform.

## References

[Becker et al. (1993)] T. Becker, V. Weispfenning, Gröbner Bases. A Computational Approach to Commutative Algebra, Springer, 1993.

[Belkoura et al. (2009)] L. Belkoura, J.-P. Richard, M. Fliess, Parameters estimation of systems with delayed and structured entries, Automatica, 45, 1117–1125 (2009).

[Chartouny et al. (2021)] M. Chartouny, T. Cluzeau, A. Quadrat, On the inverse Cauchy problem for linear ordinary differential equations, Proceedings in Applied Mathematics and Mechanics, Volume 21, Issue 1 (2021).

[Chartouny (2021)] M. Chartouny, Algorithmic Study of the Algebraic Parameter Estimation Problem, Master Thesis, University of Versailles-Saint-Quentin en Yvelines (UVSQ), Mathematics Department, 2021.

[Chyzak et al. (1998)] F. Chyzak, B. Salvy, Non-commutative elimination in Ore algebras proves multivariate identities, J. Symbolic Computation, 26, 187–227 (1998).
[Chyzak et al. (2007)] F. Chyzak, A. Quadrat, D. Robertz, OreModules: A symbolic package for the study of multidimensional linear systems, in *Applications of Time-Delay Systems, Lecture Notes in Control and Information Sciences*, 352, Springer, 233–264 (2007), https://who.rocq.inria.fr/Alban.Quadrat/OreModules.html.

[Cohn (1971)] P. M. Cohn, *Free rings and their relations*, Acad. Press, 1971.

[Coutinho (1995)] S. C. Coutinho, *A Primer of Algebraic D-Modules*, Cambridge University Press, 1995.

[Fliess et al. (2003)] M. Fliess, and H. Sira-Ramírez, An algebraic framework for linear identification, *ESAIM Control Optim. Calc. Variat.*, 9, 151–168 (2003).

[Kailath et al. (2000)] T. Kailath, A. H. Sayed, B. Hassibi, *Linear Estimation*, Prentice Hal, 2000.

[Mboup (2009)] M. Mboup, Parameter estimation for signals described by differential equations, *Applicable Analysis*, 88, 29–52 (2009).

[Quadrat et al. (2014)] A. Quadrat, A. Quadrat, Delay effects in visual tracking problems for an optronic sighting system, *Low-Complexity Controllers for Time-Delay Systems, Advances in Delays and Dynamics* 2, Springer, 2014, 77–92.

[Poor (1994)] H. V. Poor, *An Introduction to Signal Detection and Estimation*, Springer, 1994.

[Quadrat (2017)] A. Quadrat, Towards an effective study of the algebraic parameter estimation problem, *Proceedings IFAC 2017 Workshop Congress*, (Toulouse, 2017). NONA: A symbolic package for the algebraic parameter estimation problem, https://who.rocq.inria.fr/Alban.Quadrat/Non-A.html.

[Quadrat et al. (2020)] A. Quadrat, G. Regensburger, Computing polynomial solutions and annihilators of integro-differential operators with polynomial coefficients, *Algebraic Methods and Symbolic-Numeric Computation in Systems Theory, Advances in Delays and Dynamics*, 9, Springer, 2020, 87–114.

[Schwartz (1966)] L. Schwartz, *Théorie des distributions*, Hermann, 1966.

[Stafford (1978)] J. T. Stafford, Module structure of Weyl algebras, *J. Lond. Math. Soc.*, 18, 429–442 (1978).

[Ushirobira et al. (2016)] R. Ushirobira, W. Perruquetti, M. Mboup, An algebraic continuous time parameter estimation for a sum of sinusoidal waveform signals, *International Journal of Adaptive Control and Signal Processing*, 30, 1689–1713 (2016).

[Van Trees (2004)] H. L. Van Trees, *Detection, Estimation, and Modulation Theory, Part I: Detection, Estimation, and Linear Modulation Theory*, John Wiley & Sons, 2004.

# A Appendix: mathematical prerequisites

## A.1 Laplace transform

The Laplace transform is an integral transform that has been proven to be an important tool for solving ordinary or partial differential equations.

**Definition A.1.** The Laplace transform, denoted $\mathcal{L}$, maps a real-valued integrable Lebesgue measurable function $f$ to a complex-valued function $\mathcal{L}(f)$ defined by:

$$ \mathcal{L}(f)(s) := \int_{0}^{+\infty} e^{-st} f(t) dt. $$

We shall simply denote $\mathcal{L}(f)(s)$ by $\hat{f}(s)$.

If we do not make any hypothesis on $f$, then the integral $\int_{0}^{+\infty} e^{-st} f(t) dt$ does not necessarily exist. In Definition A.1, we assume that $f$ is Lebesgue integrable but the Laplace transform can be extended to *temperate distributions* [Schwartz (1966)]. Hence, if $W$ denotes the *Heaviside distribution* then $\mathcal{L}(W)(s) = s^{-1}$. Similarly, if we denote by $\delta$ the *Dirac distribution*, defined as the derivative of $W$ in the sense of the theory of distributions, i.e., $\delta = \dot{W}$, then we have $\mathcal{L}(\delta) = 1$. See [Schwartz (1966)]. In the following table, where $H$ denotes the *Heaviside function* defined by $H(t) = 1$ for $t \ge 0$ and $0$ elsewhere, we recall a few standard examples of Laplace transforms.

<table><thead><tr><td>function</td><td>Laplace transform</td></tr></thead><tbody><tr><td>$H(t)$</td><td>$1/s$</td></tr><tr><td>$\cos(t) H(t)$</td><td>$s / (s^2 + 1)$</td></tr><tr><td>$\sin(t) H(t)$</td><td>$1 / (s^2 + 1)$</td></tr><tr><td>$e^t H(t)$</td><td>$1 / (s - 1)$</td></tr><tr><td>$t^n H(t)$</td><td>$n! / s^{n+1}$</td></tr></tbody></table>
The following three classical results are used in the present paper.

**PROPOSITION A.2 ([Schwartz (1966)])**. Let $f$ be a function. The Laplace transform satisfies the following identities:

1) $\mathcal{L}(f^{(n)})(s) = s^n \mathcal{L}(f)(s) - \sum_{i=0}^{n-1} s^{n-i-1} f^{(i)}(0)$, where $f^{(n)}$ denotes the $n^{th}$ derivative of $f$,

2) $\mathcal{L}(t^n f(t))(s) = (-\partial_s)^n (\mathcal{L}(f)(s))$, where $\partial_s^n := \frac{d^n}{ds^n}$ denotes the $n^{th}$ iteration of the derivation with respect to $s$,

3) $\mathcal{L}(a_i(t)f(t))(s) = a_i(-\partial_s) \mathcal{L}(f)(s)$, where the $a_i$'s are polynomials in $t$.

**PROPOSITION A.3 ([Schwartz (1966)])**. Let $(f \star g)(t) := \int_0^t f(t-\tau)g(\tau)d\tau$ denotes the convolution product of a function $f$ by a function $g$. The Laplace transform maps a convolution product to a product, i.e., $\mathcal{L}(f \star g) = \hat{f}\hat{g}$. Conversely, the inverse Laplace transform maps a product to a convolution, i.e., $\mathcal{L}^{-1}(\hat{f}\hat{g}) = f \star g$.

## A.2 Weyl algebra

Let $k$ be a field of characteristic 0, $k[s]$ the commutative ring formed by all the polynomials in $s$ with coefficients in $k$, and $\text{End}_k(k[s])$ the ring formed by all the $k$-endomorphisms of $k[s]$, namely, $f \in \text{End}_k(k[s])$ if $f$ is a $k$-linear map from $k[s]$ to $k[s]$.

**Definition A.4 ([Coutinho (1995)])**. Let $A$ be the smallest $k$-sub-algebra of $\text{End}_k(k[s])$ generated by the following two endomorphisms of $\text{End}_k(k[s])$

$$ \begin{array}{ccc} \partial_s : k[s] & \to & k[s] \qquad s : k[s] \to k[s] \\ p : \to & \frac{dp}{ds}, & p : \to s p, \end{array} $$

which satisfy the identity $\partial_s \circ s = s \circ \partial_s + 1$ in $\text{End}_k(k[s])$, where $\circ$ denotes the composition of endomorphisms and $1$ is the identity endomorphism of $k[s]$. The $k$-algebra $A$ is called the first Weyl algebra, also denoted by $A_1(k)$.

To simplify the notations, in what follows, we shall remove the composition symbol $\circ$. For instance, the above identity between endomorphisms of $k[s]$ will be written as $\partial_s s = s \partial_s + 1$.

Let us consider the free associative $k$-algebra $k\langle D, S \rangle$ formed by all finite $k$-linear combinations of words constructed with the two letters $D$ and $S$ (e.g., $SSD$ and $SDS$ are two different words) [Cohn (1971)]. Let $J = \langle DS - SD - 1 \rangle \subset k\langle D, S \rangle$ be the two-sided ideal of $k\langle D, S \rangle$ generated by $DS - SD - 1$, namely, the set of elements defined as finite $k$-linear combinations of words of the form $W_1(DS - SD - 1)W_2$, where $W_1$ and $W_2$ are two words. If we note by $s$ (resp., $\partial$) the residue class of $S$ (resp., $D$) in $k\langle D, S \rangle/J$, then we have $k\langle D, S \rangle/J = k\langle \partial, s | \partial s = s \partial + 1 \rangle = A$.

## A.3 Gröbner bases

In the present paper, we merely use Gröbner bases methods to compute annihilators. For completeness, we briefly recall here the basics of Gröbner bases using the standard commutative setting, i.e., for the case of a polynomial ring in several commuting variables. For more details, see [Becker et al. (1993)] and the references therein.

### A.3.1 Gröbner bases for ideals over a commutative polynomial algebra.
Let $x := x_1, \dots, x_n$ be a collection of variables, $D := k[x]$ the ring of multivariate polynomials in $x_1, \dots, x_n$ with coefficients in the field $k$, $\alpha = (\alpha_1, \dots, \alpha_n)$ and $\beta = (\beta_1, \dots, \beta_n)$.

**Definition A.5.** A monomial order on $D$ is any relation $<_{\alpha}$ on the set $\{x^{\alpha} | \alpha \in \mathbb{N}^n\}$ satisfying:
(1) < is a total order on $\mathbb{N}^n$, i.e., all the elements of the set are comparable to each other,

(2) < is compatible with multiplication in $D$, i.e, if $\alpha < \beta$, then $\alpha + \gamma < \beta + \gamma$ for $\alpha, \beta, \gamma \in \mathbb{N}^n$,

(3) < is a well-ordering, i.e., any nonempty subset of $\mathbb{N}^n$ has a smaller element for <,

(4) < is an admissible term order if < is a total order which is compatible with multiplication in $D$ and for which 1 is the least element.

We implicitly set up an ordering on the variables $x_i$'s in $k[x_1, ..., x_n]: x_n < x_{n-1} < ... < x_1$. With this choice, there are still many ways to define monomial orders. Some of the most important are given in the following example.

**Example A.6.** (1) The lexicographic order on $x$-monomials is defined by $\alpha <_{lex} \beta$ whenever the first nonzero entry of $\beta - \alpha$ is positive. For instance, if we consider $\mathbb{Q}[x_1, x_2]$, then we have

$$1 <_{lex} x_2 <_{lex} x_2^2 <_{lex} x_1 <_{lex} x_1 x_2 <_{lex} x_1^2.$$

(2) The degree reverse lexicographic order on $x$-monomials is defined by $\alpha <_{tdeg} \beta$ whenever $|\alpha| < |\beta|$ or if we have $|\alpha| = |\beta|$, then the last nonzero entry of $\beta - \alpha$ is negative. For instance, if we consider $\mathbb{Q}[x_1, x_2]$, then we have

$$1 <_{tdeg} x_2 <_{tdeg} x_1 <_{tdeg} x_2^2 <_{tdeg} x_1 x_2 <_{tdeg} x_1^2.$$

(3) Let $y := y_1, ..., y_m$. Assume that an admissible monomial order $<_x$ (resp., $<_y$) on $x$-monomials (resp., on $y$-monomials) is given. An elimination order is then defined by

$$u < w < t \iff u <_x w \text{ or } (u = w \text{ and } v <_y t),$$

where $u, w$ (resp., $v, t$) are $x$-monomials (resp., $y$-monomials). An elimination order serves to eliminate the $x_i$'s. The elimination order, which will be used in what follows, is the one induced by the total degree orders on $x$-monomials and $y$-monomials. This is a very common order called lexdeg. For instance, if we consider $\mathbb{Q}[x_1, x_2, x_3]$, $\mathbf{x} = x_1$, $\mathbf{y} = x_3$, $\mathbf{x} <_{lex} \mathbf{y}$ and $\mathbf{x} <_{tdeg}$ and $\mathbf{y} <_{tdeg}$, then we have

$$1 <_{lexdeg} x_3 <_{lexdeg} x_3^2 <_{lexdeg} x_2 <_{lexdeg} x_2 x_3 <_{lexdeg} x_1 <_{lexdeg} x_1 x_3 <_{lexdeg} x_2^2 <_{lexdeg} x_1 x_2 <_{lexdeg} x_1^2.$$

**Definition A.7.** Let $<$ be a monomial order of $D$ and $P \in D$. We can then define

(1) The leading monomial $\text{lm}_<(P)$ of $P$ to be the $<-$maximal monomial that appears in $P$,

(2) The leading coefficient $\text{lc}_<(P)$ of $P$ to be the coefficient of $\text{lm}_<(P)$,

(3) The leading term $\text{lt}_<(P)$ of $P$ to be the product $\text{lc}_<(P) \text{lm}_<(P)$.

**Definition A.8.** Fix a monomial order $<$ on $D$ and let $I \subset D$ be an ideal. A Gröbner basis for $I$ (with respect to $<$) is a finite collection of polynomials $G = \{Q_1, ..., Q_t\} \subset I$ with the property that for every nonzero $P \in I$, $\text{lt}(P)$ is divisible by $\text{lt}(Q_i)$ for some $i$.

Nowadays, Gröbner bases are widely used in computer algebra and many applications. In particular, we use Buchberger's algorithm to compute a Gröbner basis of an ideal. Note that the algorithm initially due to Buchberger has been well improved in the last decades. In Maple, we have an implementation in the Groebner package.

### A.3.2 Gröbner bases for modules
Let us state the definition of a module again.

**Definition A.9.** Let $D$ be a commutative ring. A $D$-module $M$ is an abelian group $(M, +)$ equipped with a scalar multiplication

$$D \times M \to M \\ (d, m) \mapsto d m,$$

which satisfies the following properties: for all $d_1, d_2 \in D$ and for all $m_1, m_2 \in M$
$$ (1) \quad d_1 (m_1 + m_2) = d_1 m_1 + d_1 m_2, $$

$$ (2) \quad (d_1 + d_2) m_1 = d_1 m_1 + d_2 m_1, $$

$$ (3) \quad (d_2 d_1) m_1 = d_2 (d_1 m_1), $$

$$ (4) \quad 1 m_1 = m_1. $$

Note that the definition of a D-module is similar to the one of a vector space but where the scalars belong to a ring D instead of a field.

**Definition A.10.** A D-module $M$ is said to be finitely generated if $M$ admits a finite set of generators, namely there exists a finite set $S := \{m_i\}_{i=1,...,r}$ of elements of $M$ satisfying that for every $m \in M$, there exist $d_i \in D$ for $i = 1, ..., r$ such that $m = \sum_{i=1}^{r} d_i m_i$. Then, $S$ is called a set of generators of $M$.

In what follows, we consider $D$ to be the polynomial algebra $A = k[x_1, ..., x_n]$. Let $\text{Mon}(A)$ be the set of monomials of $A$ and $\{f_j\}_{j=1,...,p}$ the standard basis of the free finitely generated left $A$-module $A^{1 \times p} := \{(\lambda_1, ..., \lambda_p) | \lambda_i \in A, i = 1, ..., p\}$, namely the $k^{th}$ component of $f_j$ is 1 if $k = j$ and 0 otherwise. First, we can extend a monomial order $<$ from $\text{Mon}(A)$ to the set of monomials of the form $u f_j$, where $u \in \text{Mon}(A)$ and $j = 1, ..., p$. This extension is also denoted by $<$ and it has to satisfy the following two conditions:

$$ (1) \quad \forall w \in \text{Mon}(A) : u f_i < v f_j \implies w u f_i < w v f_j. $$

$$ (2) \quad u < v \implies u f_j < v f_j \text{ for } j = 1, ..., p. $$

Without loss of generality, we let $f_p < f_{p-1} < \dots < f_1$. There are two natural extensions of a monomial order to $\text{Mon}(A^{1 \times p})$.

**Definition A.11.** Let $<$ be an admissible monomial order on $\text{Mon}(A)$ and $\{f_j\}_{j=1,...,p}$ the standard basis of the left $A$-module $A^{1 \times p}$.

(1) The term over position order on $\text{Mon}(A^{1 \times p})$ induced by $<$ is defined by

$$ u f_i < v f_j \iff u < v \text{ or } u = v \text{ and } f_i < f_j. $$

(2) The position over term order on $\text{Mon}(A^{1 \times p})$ induced by $<$ is defined by

$$ u f_i < v f_j \iff f_i < f_j \text{ or } f_i = f_j \text{ and } u < v. $$

The term over position order is of more computational value concerning efficiency. The position over term order can be used to eliminate components.

**A.3.3 Computation of a left kernel.** We now give an algorithm that computes the left kernel of a matrix $R \in A^{q \times p}$, namely, the set of row vectors $\lambda \in A^{1 \times q}$ which are such that $\lambda R = 0$. The idea is to consider the inhomogeneous linear system $R \eta = \zeta$ to eliminate the $\eta_i$'s from these equations, and to select the equations in the $\zeta_i$'s only. In other words, computing the left kernel of $R$ is equivalent to computing a generating set of compatibility conditions for the inhomogeneous linear system

$$ R \eta = \zeta. $$

**Algorithm:** Computation of the left kernel of $R \in A^{q \times p}$, i.e., find $S \in A^{r \times q}$ such that:

$$ \ker_A(R) := \{\lambda \in A^{1 \times q} | \lambda R = 0\} = A^{1 \times r} S := \{\mu S | \mu \in A^{1 \times r}\}. $$

* **Input:** $R \in A^{p \times q}$.

* **Output:** A matrix $S \in A^{r \times q}$ such that $\ker_A(R) = A^{1 \times r} S$.

(1) Introduce the indeterminates $\eta_1, \dots, \eta_p, \zeta_1, \dots, \zeta_q$ over $A$ and define the set:

$$ P := \left\{ \sum_{j=1}^{p} R_{ij} \eta_j - \zeta_i \mid i = 1, \dots, q \right\}. $$
(2) Compute a Gröbner basis $G$ of $P$ in the free $A$-module generated by the $\eta_j$'s and the $\zeta_i$'s for $j = 1, \dots, p$ and $i = 1, \dots, q$, namely, $\bigoplus_{j=1}^p A\eta_j \oplus \bigoplus_{i=1}^q A\zeta_i$ with respect to a term order which eliminates the $\eta_j$'s.

(3) Compute $G \cap \bigoplus_{i=1}^q A\zeta_i = \sum_{i=1}^q S_{ki}\zeta_i$, $k = 1, \dots, r$ by selecting the elements of $G$ containing only the $\zeta_i$'s, and return the matrix $S := (S_{ij}) \in A^{r \times q}$.

We finally illustrate how the above algorithm is used in the present paper to compute the annihilator of a polynomial over the Weyl algebra. Let us consider the first Weyl algebra $A := A_1(k) = k[s]\langle\partial_s | \partial_s s = s \partial_s + 1\rangle$ and characterize the annihilator of a polynomial $Q \in k[s]$, namely, the left ideal of $A$ defined by $\text{ann}_A(Q) := \{a \in A | aP := a(P) = 0\}$. If $d := \text{deg}_S(Q)$, then this annihilator can be obtained by considering the polynomial relations between the $d+1$ first derivatives of $Q$, i.e., by considering the following left kernel:

$$ \ker_{k[s]}(\cdot, L) := \{ \lambda := (\lambda_0 \dots \lambda_{d+1}) \in k[s]^{1 \times (d+2)} | \lambda L = 0 \}, \quad L := \begin{pmatrix} Q \\ Q^{(1)} \\ \vdots \\ Q^{(d+1)} \end{pmatrix} \in k[s]^{(d+2) \times 1}. $$

More precisely, if $\ker_{k[s]}(\cdot, L) = \text{im}_{k[s]}(\cdot, S) := k[s]^{1 \times r} S$ with $S \in k[s]^{r \times (d+2)}$ and $r \in \mathbb{N}$, then we have:

$$ \text{ann}_A(Q) = S \begin{pmatrix} 1 \\ \partial_s \\ \vdots \\ \partial_s^{d+1} \end{pmatrix}. $$

*Example A.12.* Let us compute the annihilator of the polynomial $Q(s) = s^2 + \omega^2 \in k[s]$, where $k = \mathbb{Q}[\omega]$. Then, we have $d = 2$. We consider the following matrix:

$$ L := \begin{pmatrix} s^2 + \omega^2 \\ 2s \\ 2 \\ 0 \end{pmatrix} \in k[s]^{4 \times 1}. $$

Using the algorithm described above, we can compute the left kernel $\ker_{k[s]}(\cdot, L)$. To do that, we first consider $\lambda$ and $\mu = (\mu_1, \dots, \mu_4)^T$ and write the entries of $L\lambda - \mu$, namely:

$$ P := \{(s^2 + \omega^2) \lambda - \mu_1, 2s\lambda - \mu_2, 2\lambda - \mu_3, -\mu_4\}. $$

Computing a Gröbner basis for $P$ with respect to a (lexdeg) term order which eliminates $\lambda$ (we can consider $\mu_4 < \dots < \mu_1 < \lambda$), we get $G := \{\mu_4, s\mu_3 - \mu_2, \omega^2\mu_3 + s\mu_2 - 2\mu_1, 2\lambda + \mu_3\}$. Then, in $G$, we can consider the elements which only contain the $\mu_i$'s, i.e., $\{\omega^2\mu_3 + s\mu_2 - 2\mu_1, s\mu_3 - \mu_2, \mu_4\}$. Hence, we have:

$$ \underbrace{\begin{pmatrix} -2 & s & \omega^2 & 0 \\ 0 & -1 & s & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}}_{S} \begin{pmatrix} \mu_1 \\ \mu_2 \\ \mu_3 \\ \mu_4 \end{pmatrix} = \begin{pmatrix} \omega^2 \mu_3 + s \mu_2 - 2\mu_1 \\ s \mu_3 - \mu_2 \\ \mu_4 \end{pmatrix}, $$
which shows that $ker_{k[s]}(.L) = k[s]^{1 \times 3} S$. Finally, we obtain

$$ \mathrm{ann}_{A_1(k)}(.Q) = S \begin{pmatrix} 1 \\ \partial_s \\ \partial_s^2 \\ \partial_s^3 \end{pmatrix} = \begin{pmatrix} \omega^2 \partial_s^2 + s \partial_s - 2 \\ s \partial_s^2 - \partial_s \\ \partial_s^3 \end{pmatrix}. $$

In Maple, we can compute the kernel of $L$ using the OreModules package that contains the SyzygyModule command which computes this kernel as follows:

A := DefineOreAlgebra(diff=[ds,s],polynom=[s],comm=[omega]):

L := Vector(column)(4,[s^2+omega^2, 2*s, 2, 0]):

S:= SyzygyModule(L,A);

$$ S = \begin{bmatrix} -2 & s & \omega^2 & 0 \\ 0 & -1 & s & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} $$

Finally, note that if $k = \mathbb{Q}$, then we can characterize $ker_{k[s]}(.L)$ by simply computing the *Smith normal form* of the univariate polynomial matrix $L$.