# LPV system identification with unknown scheduling variable

Xavier Bombois, Guillaume Mercère, Kévin Colin

► To cite this version:

Xavier Bombois, Guillaume Mercère, Kévin Colin. LPV system identification with unknown scheduling variable. 2026. hal-05505584

HAL Id: hal-05505584

https://hal.science/hal-05505584v1

Preprint submitted on 11 Feb 2026

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
# LPV system identification with unknown scheduling variable

X. Bombois<sup>a,b</sup>, G. Mercère<sup>c</sup>, K. Colin<sup>d</sup>

<sup>a</sup> Laboratoire Ampère, Ecole Centrale de Lyon, Université de Lyon, Ecully, France

<sup>b</sup> Centre National de la Recherche Scientifique (CNRS), France

<sup>c</sup> Université de Poitiers, ISAE-ENSMA Poitiers, LIAS, Poitiers, France

<sup>d</sup> Université de Lorraine, CNRS, CRAN, F-54000 Nancy, France

## ARTICLE HISTORY

Compiled February 11, 2026

## ABSTRACT

A Linear Parameter Varying (LPV) system is a system whose parameters depend on an exogenous variable, the so-called scheduling variable. LPV systems are generally controlled using LPV controllers i.e., controller whose parameters also depend on the scheduling variable. In order to be able to use an LPV controller for an LPV system, it is generally assumed that the scheduling variable of the LPV system is either fully known or measurable. Here, we consider the case where the scheduling variable has to be identified based on data collected on the system. We do that in the case of LPV systems that are successively operated at different constant values of the scheduling variable. We show how to identify the maximum likelihood estimate of the unknown constant scheduling variable and we derive the statistical properties of this estimate. Moreover, we also develop an optimal experimental design framework in order to optimally design the experimental conditions of the identification experiment yielding the estimate of the scheduling variable.

## KEYWORDS

LPV systems, System identification, Optimal Experiment Design

## 1. Introduction

Due to its flexibility and the ability to represent many real-life dynamical systems, the Linear Parameter Varying (LPV) framework has known a lot of interest in the last decades, leading to efficient control design methods (Apkarian & Gahinet, 1995; Apkarian, Gahinet, & Biannic, 1994; Dinh, Scorletti, Fromion, & Magarotto, 2005; Scorletti & Ghaoui, 1998) and efficient identification methods (Toth, 2010). LPV systems are systems whose parameters vary as a function of a (possibly time-varying) scheduling variable $p(t)$ (Toth, 2010).

In the LPV literature, the scheduling variable $p(t)$ can take different forms (Toth, 2010). In the so-called quasi-LPV framework, the scheduling variable $p(t)$ is a (possibly complex) function of the internal variables of the system i.e., its input, its output and/or its state vector. A opposed to this, in the classical definition of an LPV system, the scheduling variable $p(t)$ is an external variable (e.g., the temperature, the altitude, the payload, ...). In some cases, this external variable $p(t)$ is time-varying. In

CONTACT X. Bombois. Email: xavier.bombois@ec-lyon.fr
other cases, the external variable remains constant ($p(t) = p$) for a (relatively) long period of time. In this paper, we will restrict attention to this particular type of LPV systems which, in the literature, are often referred as *frozen* LPV systems or Linear Parameter Invariant (LPI) systems (see, e.g., (Dinh et al., 2005)). The value **p** of the constant scheduling variable is in turn often called *operating point* in the literature (Toth, 2010) and it is important to note that, at a given operating point **p**, the LPV system reduces to a Linear Time Invariant (LTI) system whose dynamics depend on **p**.

*Frozen* LPV systems are particularly adapted to, e.g., describe batch-processed systems. In batch-processed engineering, the dynamics of the to-be-described system can depend on a particular property **p** of the raw material that has to be treated/transformed in a given batch. This property **p** can be different from batch to batch, but it obviously remains constant for the duration of a given batch. An example is the batch-processed rubber milling system described in (Hamdan, Coirault, Mercère, & Dairay, 2024) where **p** is the viscosity of the to-be-treated rubber, a property which is known to vary from batch to batch.

In the LPV control literature, control design techniques have been specifically devised for these *frozen* LPV systems (see, e.g., (Dinh et al., 2005)). Like for other LPV systems, the resulting controller is also a controller whose parameters are, such as the system itself, a function of the scheduling variable **p** i.e., $K(\mathbf{p})$. As shown in, e.g., (Bombois, Ghosh, Scorletti, & Huillery, 2021), such a LPV controller can be designed based on a model of the LPV system. Determining such an LPV model is the objective of LPV system identification. Among the LPV identification methods (Toth, 2010), the so-called *local* approach is particularly adapted for the identification of *frozen* LPV systems (see, e.g., (Ghosh, Bombois, Huillery, Scorletti, & Mercere, 2018; Liu, Gross, Pfeiffer, & Werner, 2014; Steinbuch, van de Molengraft, Heuberger, & der Voort, 2003; Wassink, van de Wal, Scherer, & Bosgra, 2005)).

In the local approach, a model of the LPV system is obtained through the interpolation of a number of so-called *local models* obtained at different operating points $\mathbf{p}_i$ ($i = 1, ..., M$). The local model at the operating point $\mathbf{p}_i$ is a model of the LTI dynamics of the LPV system at $\mathbf{p}_i$ and it can therefore be derived using classical LTI identification (e.g., prediction error identification (Ljung, 1999)) using data collected on the LPV system operated at $\mathbf{p}_i$. In (Ghosh et al., 2018), it is shown how the interpolation yielding the LPV model can be achieved in an optimal way.

LPV identification techniques derive the LPV model under the assumption that the scheduling variable is either fully known or measurable (see, e.g., (Piga, Cox, Toth, & Laurain, 2015; Toth, 2010)). Similarly, when an LPV controller is available, its implementation requires knowledge of the value of the scheduling variable corresponding to the current batch. In certain industrial applications, if it is possible, in the commissioning phase, to operate the system at a certain number of known operating points in order to identify the LPV model, it is not possible to obtain a measurement of the scheduling variable at the beginning of each batch in order to operate the LPV controller (i.e., the controller designed based on the identified LPV model). Consequently, in order to operate this LPV controller, the unknown scheduling variable (say $\mathbf{p}^0$) will have to be identified at the beginning of each batch. A good example of this situation is the batch-processed rubber milling system described in (Hamdan et al., 2024).
This paper addresses precisely the problem described above i.e., the identification of the unknown scheduling variable $\mathbf{p}^0$ at the beginning of a batch. More precisely, we develop a maximum likelihood approach to identify an estimate of $\mathbf{p}^0$. This estimate will be obtained by optimally combining the LPV model derived in the commissioning phase and a local model at $\mathbf{p}^0$ which is identified based on data collected on the LPV system operated at $\mathbf{p}^0$. We also give an expression for the variance of this estimate of $\mathbf{p}^0$. Finally, assuming that the LPV model (obtained in the commissioning phase) is given, we provide optimal experiment design tools for the identification experiment yielding the local model at $\mathbf{p}^0$. These optimal experiment tools allow to determine the experimental conditions of this experiment in such a way that the resulting estimate of $\mathbf{p}^0$ has the smallest possible variance.

The maximum likelihood approach to estimate $\mathbf{p}^0$ which is presented in this paper extends, to a more realistic situation, the approach in (Hamdan et al., 2024) where it is assumed that noise-free data can be used in the identification of both the LPV model and the local model at $\mathbf{p}^0$.

As indicated above, the approach presented in this paper to estimate the unknown scheduling variable is adapted to the case of frozen LPV systems. Consequently, if the system is a quasi-LPV system where $p(t)$ is assumed to be a complex function of the internal variables of the system, the approach presented in the present paper is no longer appropriate and other techniques have then to be used instead, e.g., the techniques in the papers (Rehmer & Kroll, 2022; Verhoek, Beintema, Haesaert, Schoukens, & Tóth, 2022; Yang, Liu, & Han, 2018). In these papers, the relation between the internal variables of the system and the time-varying scheduling variable $p(t)$ is modeled using a black-box model and the parameters of this black-box model are identified based on data collected on the to-be-modeled system. Due to the inherent complexity of the problem, the papers (Rehmer & Kroll, 2022; Verhoek et al., 2022; Yang et al., 2018) do not provide optimal experiment design tools to tune the uncertainty of the identified model.

*Notations.* Discrete-time signals will be denoted $x(t)$ where $t$ is the sample number. The variable $z$ will denote both the Z-transform variable and the shift operator. For a matrix $A$, $A^T$, $\lambda_{max}(A)$ and $det(A)$ are, respectively, its transpose, its largest eigenvalue and its determinant. With some abuse, we will use 0 both for zero and for zero matrices. The matrix

$$ \begin{pmatrix} X_1 & 0 & 0 \\ 0 & \ddots & 0 \\ 0 & 0 & X_n \end{pmatrix} $$

will be denoted $diag(X_1, ..., X_n)$ if the elements $X_i$ ($i = 1, ..., n$) are scalar quantities, while it will be denoted $bdiag(X_1, ..., X_n)$ if the elements $X_i$ ($i = 1, ..., n$) are matrices. In addition, $I_n$ represents the identity matrix of dimension $n \times n$ and $\otimes$, the Kronecker product. Finally, $E$ denotes the expected value operator.
## 2. Problem description

### 2.1. True LPV system

In this paper, we consider a Single-Input Single-Output (SISO) LPV system i.e. a SISO system whose parameters vary with time as a function of an exogenous variable, the so-called scheduling variable $p(t)$. We assume that $p(t)$ is a scalar signal that can vary in a given scheduling interval $\mathcal{P} = [p_{min}, p_{max}]$. The considered SISO LPV system is as follows:

$$ 
\begin{aligned}
y(t) &= \check{y}(t) + \tilde{v}(t) \\
\text{with } A(z, \zeta^0(p(t))) \check{y}(t) &= B(z, \zeta^0(p(t))) u(t) \quad (1) \\
\text{and } D(z, \zeta^0(p(t))) \tilde{v}(t) &= C(z, \zeta^0(p(t))) e(t)
\end{aligned}
$$

with $u(t)$ the input signal, $y(t)$ the measured output and $e(t)$ a white noise signal of variance $\sigma_e^2$. In (1), $A(z, \zeta^0)$, $B(z, \zeta^0)$, $C(z, \zeta^0)$ and $D(z, \zeta^0)$ are polynomials in the delay operator $z^{-1}$ and having arbitrary orders: $A(z, \zeta^0) = 1 + a_1^0 z^{-1} + \dots + a_{n_a}^0 z^{-n_a}$, $B(z, \zeta^0) = z^{-n_k}(b_0^0 + b_1^0 z^{-1} + \dots + b_{n_b}^0 z^{-n_b})$, $C(z, \zeta^0) = 1 + c_1^0 z^{-1} + \dots + c_{n_c}^0 z^{-n_c}$ and $D(z, \zeta^0) = 1 + d_1^0 z^{-1} + \dots + d_{n_d}^0 z^{-n_d}$. The vector $\zeta^0$ in (1) is a column vector of dimension $n \triangleq n_a + n_b + 1 + n_c + n_d$ containing the coefficients of the four polynomials describing the LPV system i.e. $a_1^0, \dots, b_0^0, \dots, d_{n_d}^0$. All the entries $\zeta_i^0$ ($i = 1, \dots, n$) of $\zeta^0$ depend on the scheduling variable $p(t)$ as follows:

$$ \zeta_i^0(p(t)) = \sum_{j=0}^{n_p} \zeta_{i,j}^0 f_j(p(t)) \quad i = 1, \dots, n \quad (2) $$

where $n_p$ is an arbitrary order and $f_j(p(t))$ ($j = 0, \dots, n_p$) correspond to a set of unisolvent basis functions, e.g. $f_j(p(t)) = p^j(t)$ (see e.g. (Dankers, Toth, Heuberger, Bombois, & den Hof, 2011)). In (2), we suppose that $n_p$ is the same for each entry $\zeta_i^0$ ($i = 1, \dots, n$) in order to simplify the notations (this is nevertheless not a necessity). In the sequel, we will assume that the structure of the LPV system is known i.e. that $n_a, n_b, n_k, n_c$ and $n_d$ are known in (1) and that $n_p$ is known in (2).

Let us define by $\theta^0 = (\zeta_{1,0}^0, \zeta_{1,1}^0, \dots, \zeta_{n,n_p}^0)^T$, the vector of dimension $k \triangleq (n_p + 1)n$ which contains all the time-invariant coefficients of the polynomial expansions in (2). The relation between $\theta^0$ and the time-varying parameter vector $\zeta^0(p(t))$ can thus be expressed as follows:

$$ \zeta^0(p(t)) = T(p(t)) \theta^0, \quad (3) $$

for a given matrix $T(p(t))$ of dimension $n \times k$ which is uniquely function of $f_j(p(t))$ ($j = 0, \dots, n_p$). More precisely, using (2), we have that $T(p(t)) = I_n \otimes (f_0(p(t)), f_1(p(t)), \dots, f_{n_p}(p(t)))$. From the latter expression, it is clear that $T(p(t))$ is a full row rank matrix.

When the scheduling variable is equal to a constant value $\mathbf{p} \in \mathcal{P}$ i.e.,

$$ p(t) = \mathbf{p}, \quad (4) $$
for a sufficiently long interval of time, the scheduling variable is then called an operating point (Toth, 2010) and the dynamics of the LPV system (1) can be represented by the following LTI system:

$$y(t) = G(z, \zeta^0(\mathbf{p})) u(t) + H(z, \zeta^0(\mathbf{p})) e(t), \quad (5)$$

where $G(z, \zeta) = B(z, \zeta)/A(z, \zeta)$ and $H(z, \zeta) = C(z, \zeta)/D(z, \zeta)$ and where $\zeta^0(\mathbf{p}) \in \mathbb{R}^n$ is a time-invariant parameter vector defined as:

$$\zeta^0(\mathbf{p}) = T(\mathbf{p})\theta^0. \quad (6)$$

As already mentioned above, $T(\mathbf{p})$ is a full row rank matrix for any value of $\mathbf{p}$. Let us also make the following very mild identifiability assumption about (6):

**Assumption 2.1.** The mapping (6) and the scheduling interval $\mathcal{P}$ are such that

$$T(\mathbf{p}_1) \theta^0 \neq T(\mathbf{p}_2) \theta^0$$

for any arbitrary scalars $\mathbf{p}_1$ and $\mathbf{p}_2$ in $\mathcal{P}$ such that $\mathbf{p}_1 \neq \mathbf{p}_2$.

In this paper, we will restrict attention to LPV systems as described by equations (4)-(6) i.e., frozen LPV systems. As indicated in the introduction, this type of LPV systems are particularly adapted to, e.g., describe batch-processed systems. In batch-processed engineering, the dynamics of the to-be-described system can depend on a particular property $\mathbf{p}$ of the raw material that has to be treated/transformed in a given batch. This property $\mathbf{p}$ can be different from batch to batch, but it obviously remains constant for the duration of a given batch. An example is the batch-processed rubber milling system described in (Hamdan et al., 2024) where $\mathbf{p}$ is the viscosity of the to-be-treated rubber, a property which is known to vary from batch to batch.

Within this frozen LPV framework, a local LPV identification procedure such as the one in (Ghosh et al., 2018) allows to determine an estimate $\hat{\theta}_l$ of $\theta^0$ and, based on that estimate $\hat{\theta}_l$, an LPV controller $K(\mathbf{p})$ can be easily designed (Dinh et al., 2005). This LPV controller $K(\mathbf{p})$ is a controller whose parameters depend on the operating point $\mathbf{p}$. If we would know or be able to measure the current value $\mathbf{p}^0$ of the operating point at the beginning of each batch, we would then also be able to determine the particular controller $K(\mathbf{p}^0)$ that has to be applied during this batch. However, in this paper, we consider a situation where knowing or measuring the constant scheduling variable is only possible in the commissioning phase of the LPV system i.e., the phase leading to the estimate $\hat{\theta}_l$ and to the LPV controller $K(\mathbf{p})$. We are thus in a situation where knowing or measuring the constant operating point $\mathbf{p}^0$ is no longer possible in the operation phase. Consequently, at the beginning of each batch, the unknown operating point $\mathbf{p}^0$ will need to be identified based on data collected on the system.

If it is clear that, with classical prediction error identification, we can easily obtain an estimate $\hat{\zeta}_f$ of the parameter vector $\zeta^0(\mathbf{p}^0) = T(\mathbf{p}^0)\theta^0$ describing the LPV system at the operating point $\mathbf{p}^0$, it is also clear that, since $\theta^0$ is unknown, we cannot derive an estimate of $\mathbf{p}^0$ based on $\hat{\zeta}_f$ alone. An estimate of $\mathbf{p}^0$ will therefore only be
derivable by combining the estimate $\hat{\zeta}_f$ (obtained in the operation phase) with the estimate $\hat{\theta}_l$ (obtained in the commissioning phase). Before showing how this can be done in practice, let us first present in details how these estimates $\hat{\zeta}_f$ and $\hat{\theta}_l$ can be determined and also how their respective uncertainty can be evaluated.

## 2.2. Identification of an estimate $\hat{\zeta}_f$ of the vector of coefficients $\zeta^0(\mathbf{p}^0) = T(\mathbf{p}^0)\theta^0$ in the operation phase

In the operation phase, at the beginning of each batch, the LPV system will be described by (5) operated at an unknown scheduling variable $\mathbf{p}^0$ i.e.,

$$y(t) = G(z, \zeta^0(\mathbf{p}^0))) u(t) + H(z, \zeta^0(\mathbf{p}^0))) e(t) \quad (7)$$

$$\text{with } \zeta^0(\mathbf{p}^0) = T(\mathbf{p}^0)\theta^0. \quad (8)$$

In this system description, both $\mathbf{p}^0 \in \mathbb{R}$ and $\theta^0 \in \mathbb{R}^k$ are unknown. For further reference, let us define the vector $\eta^0$ which gathers these two unknown quantities:

$$\eta^0 = \begin{pmatrix} \mathbf{p}^0 \\ \theta^0 \end{pmatrix} \in \mathbb{R}^{k+1} \quad (9)$$

Let us now show how we can identify an estimate $\hat{\zeta}_f$ of the unknown parameter vector (8). Since $\zeta^0(\mathbf{p}^0)$ is the parameter vector describing the LTI system (7), identifying an estimate $\hat{\zeta}_f$ of $\zeta^0(\mathbf{p}^0)$ boils down to a standard (open-loop) identification experiment for a LTI system (Ljung, 1999). Let us describe this identification experiment. If we apply a sufficiently informative input signal $u_f(t)$ (e.g., a (filtered) white noise) to the system (7) for $t = 1, ..., N_f$ and collect the corresponding output $y_f(t)$, we can gather the data set $Z^{N_f} = \{u_f(t), y_f(t) | t = 1, ..., N_f\}$. Based on this data set $Z^{N_f}$ and a full order model structure $M = \{G(z, \zeta) \ H(z, \zeta) | \zeta \in \mathbb{R}^n\}$ for the LTI representation (7), we can determine an estimate $\hat{\zeta}_f$ of $\zeta^0(\mathbf{p}^0)$ using prediction error identification (Ljung, 1999):

$$\hat{\zeta}_f = \arg \min_{\zeta} \frac{1}{N_f} \sum_{t=1}^{N_f} \epsilon_f^2(t, \zeta) \quad (10)$$

with $\epsilon_f(t, \zeta) = H^{-1}(z, \zeta)(y_f(t) - G(z, \zeta)u_f(t))$. This estimate has the following asymptotic distribution:

$$\hat{\zeta}_f \sim \mathcal{N}(\zeta^0(\mathbf{p}^0), P_{\zeta_f}) \quad (11)$$

where the covariance matrix $P_{\zeta_f}$ has the following structure

$$P_{\zeta_f} = \frac{1}{N_f} \mathcal{I}_{\zeta_f}^{-1} \quad (12)$$
with $\mathcal{I}_{\zeta_f}$ a strictly positive-definite matrix given by

$$ \mathcal{I}_{\zeta_f} = \left( \frac{1}{\sigma_e^2} \frac{1}{2\pi} \int_{-\pi}^{\pi} F_u(e^{j\omega}, \zeta^0(\mathbf{p}^0)) F_u^*(e^{j\omega}, \zeta^0(\mathbf{p}^0)) \Phi_{u_f}(\omega) d\omega \right) + \dots \\ \dots + \left( \frac{1}{2\pi} \int_{-\pi}^{\pi} F_e(e^{j\omega}, \zeta^0(\mathbf{p}^0)) F_e^*(e^{j\omega}, \zeta^0(\mathbf{p}^0)) d\omega \right) > 0 \quad (13) $$

with $\Phi_{u_f}(\omega)$ the power spectrum of the excitation signal $u_f(t)$ and with $F_u(z, \zeta) = \frac{1}{H(z, \zeta)} \frac{\partial G(z, \zeta)}{\partial \zeta}$ and $F_e(z, \zeta) = \frac{1}{H(z, \zeta)} \frac{\partial H(z, \zeta)}{\partial \zeta}$. For the sequel, it is important to note that $\mathcal{I}_{\zeta_f}$ does not depend on $N_f$ and that, since the power of the excitation signal $u_f$ will always be limited, we have that $\lambda_{max}(\mathcal{I}_{\zeta_f}) < \infty$. Finally, note also that, since $P_{\zeta_f} \to 0$ when $N_f \to \infty$, $\hat{\zeta}_f$ converges to $\zeta^0(\mathbf{p}^0) = T(\mathbf{p}^0)\theta^0$ with probability 1 when $N_f \to \infty$.

**Remark.** It is shown in (Ljung, 1999) that an estimate of the covariance matrix $P_{\zeta_f}$ (see (12)-(13)) can be deduced based on the data $Z^{N_f}$ and $\hat{\zeta}_f$. Moreover, even if the statistical property (11) is in fact only valid when $N_f \to \infty$, (11) is also widely used for finite $N_f$ (Ljung, 1999). In other words, the estimate of $P_{\zeta_f}$ that can be deduced from the data $Z^{N_f}$ and $\hat{\zeta}_f$ is used as a measure of the uncertainty of $\hat{\zeta}_f$.

## 2.3. Deriving an estimate $\hat{\theta}_l$ of $\theta^0$ using classical local LPV identification in the commissioning phase

As indicated in the introduction, different local LPV identification techniques can be used to derive an estimate $\hat{\theta}_l$ of $\theta^0$. Let us here focus on the technique presented in (Ghosh et al., 2018) since it has the advantage to not only deliver an estimate $\hat{\theta}_l$ of $\theta^0$, but also an expression for the covariance matrix $P_{\theta,l}$ of $\hat{\theta}_l$. Like in all other local LPV identification techniques, the technique in (Ghosh et al., 2018) is based on the identification of an estimate $\hat{\zeta}_{\mathbf{p}_i}$ of the constant parameter vector $\zeta^0(\mathbf{p}_i) = T(\mathbf{p}_i)\theta^0$ i.e., the vector of coefficients of the transfer functions describing the LPV system (5) operated at a fixed and known operating point $\mathbf{p}_i$. The estimate $\hat{\zeta}_{\mathbf{p}_i}$ of $\zeta^0(\mathbf{p}_i) = T(\mathbf{p}_i)\theta^0$ can be determined similarly as in the previous subsection i.e., using prediction error identification with a set of $N_l$ input-output data collected on (5) with $\mathbf{p} = \mathbf{p}_i$. As indicated in the previous subsection, this estimate $\hat{\zeta}_{\mathbf{p}_i}$ is then (asymptotically) normally distributed around $\zeta^0(\mathbf{p}_i) = T(\mathbf{p}_i)\theta^0$ with a covariance matrix $P_{\zeta, \mathbf{p}_i}$ that can be easily estimated.

If this procedure is repeated for a number $M$ of different operating points $\mathbf{p}_i$ ($i = 1, ..., M$), we have thus:

$$ \hat{\zeta}_{\mathbf{p}_i} \sim \mathcal{N}(\zeta^0(\mathbf{p}_i), P_{\zeta, \mathbf{p}_i}) \quad i = 1, ..., M \quad (14) $$

Using the information (14) collected in the $M$ local identification experiments of duration $N_l$ and assuming that $M \ge n_p + 1$ with $n_p$ as defined in (2) (Ghosh et al., 2018), we can show that an efficient estimate<fnref n="1" /> $\hat{\theta}_l$ of $\theta^0$ is given by:

<footnote n="1">¹Using similar arguments as in Section 3, we can indeed show that, given the information (14), $\hat{\theta}_l$ is the maximum likelihood estimate of $\theta^0$.</footnote>
$$ \hat{\theta}_l = \arg \min_{\theta} V_l(\theta) \quad \text{with} \qquad (15) $$

$$ V_l(\theta) = \frac{1}{2} \sum_{i=1}^{M} \left( T(\mathbf{p}_i)\theta - \hat{\zeta}_{\mathbf{p}_i} \right)^T P_{\zeta, \mathbf{p}_i}^{-1} \left( T(\mathbf{p}_i)\theta - \hat{\zeta}_{\mathbf{p}_i} \right) \qquad (16) $$

A closed-form solution of the weighted least-squares problem (15)-(16) is:

$$ \hat{\theta}_l = \left( \sum_{i=1}^{M} T^T(\mathbf{p}_i) P_{\zeta, \mathbf{p}_i}^{-1} T(\mathbf{p}_i) \right)^{-1} \left( \sum_{i=1}^{M} T^T(\mathbf{p}_i) P_{\zeta, \mathbf{p}_i}^{-1} \hat{\zeta}_{\mathbf{p}_i} \right) \qquad (17) $$

and it can be shown (Ghosh et al., 2018) that, asymptotically in $N_l$, this estimate $\hat{\theta}_l$ has the following properties:

$$ \hat{\theta}_l \sim \mathcal{N}(\theta^0, P_{\theta,l}) \quad \text{with} \qquad (18) $$

$$ P_{\theta,l} = \left( \sum_{i=1}^{M} T^T(\mathbf{p}_i) P_{\zeta, \mathbf{p}_i}^{-1} T(\mathbf{p}_i) \right)^{-1} . \qquad (19) $$

As was the case in the previous subsection, we will also use (18) for finite $N_l$ and since (19) can be easily evaluated, we will use $P_{\theta,l}$ as the measure of uncertainty for the estimate $\hat{\theta}_l$.

As shown in (12)-(13), each of the covariance matrices $P_{\zeta, \mathbf{p}_i}$ appearing in (19) has a structure $P_{\zeta, \mathbf{p}_i} = \frac{1}{N_l} \mathcal{I}_{\zeta, \mathbf{p}_i}^{-1}$ with $\mathcal{I}_{\zeta, \mathbf{p}_i}^{-1}$ a matrix which is not function of $N_l$ and which has also the following other properties: $\mathcal{I}_{\zeta, \mathbf{p}_i} > 0$ and $\lambda_{max}(\mathcal{I}_{\zeta, \mathbf{p}_i}) < \infty$. The covariance matrix $P_{\theta,l}$ given in (19) can thus be rewritten as:

$$ P_{\theta,l} = \frac{1}{N_l} \left( \underbrace{\sum_{i=1}^{M} T^T(\mathbf{p}_i) \mathcal{I}_{\zeta, \mathbf{p}_i} T(\mathbf{p}_i)}_{=\mathcal{I}_{\theta,l}} \right)^{-1} \qquad (20) $$

Since $M \ge n_p + 1$, it is shown in (Ghosh et al., 2018) that $\mathcal{I}_{\theta,l}$ is strictly positive definite. Moreover, $\mathcal{I}_{\theta,l}$ is not function of $N_l$ and is such that $\lambda_{max}(\mathcal{I}_{\theta,l}) < \infty$. Using (20), it is also clear that, when $N_l \to \infty$, $P_{\theta,l}$ converges to the zero matrix and thus that $\hat{\theta}_l$ converges to $\theta^0$ with probability one.
<figure><img src="image_2.png" /><figcaption>Figure 1.: Schematic view of the proposed methodology</figcaption></figure>

**Remark.** Note that the original paper (Ghosh et al., 2018) does not give the interpretation that $\hat{\theta}_l$ in (17) is the solution of the weighted least-squares problem (15)-(16).

### 3. Maximum likelihood estimate of $\eta^0$

Let us sum up the previous section. We consider an LPV system (5) whose parameter vector $\zeta^0$ depends on the unknown scheduling variable $\mathbf{p}^0$ via the expression (8) in which $\theta^0$ is an unknown parameter vector. Our objective is to obtain an accurate estimate of the unknown scheduling variable $\mathbf{p}^0$. We assume for this purpose that, in the commissioning phase, $M$ identification experiments of duration $N_l$ have delivered an estimate $\hat{\theta}_l$ for $\theta^0$ and that $\hat{\theta}_l$ has the distribution (18). We also assume that, in the operation phase, an identification experiment of duration $N_f$ has delivered an estimate $\hat{\zeta}_f$ of $\zeta^0(\mathbf{p}^0) = T(\mathbf{p}^0)\theta^0$ which has the distribution given in (11). We have thus an estimate $\hat{\theta}_l$ of $\theta^0$ and an estimate $\hat{\zeta}_f$ of $\zeta^0(\mathbf{p}^0) = T(\mathbf{p}^0)\theta^0$ and these two estimates $\hat{\theta}_l$ and $\hat{\zeta}_f$ are statistically independent. In this section, we will show how we can derive in an optimal way the desired estimate of $\mathbf{p}^0$ from $\hat{\theta}_l$ and $\hat{\zeta}_f$. For this purpose, we will derive the maximum likelihood estimate of $\hat{\eta}_{comb} \triangleq (\hat{\mathbf{p}}_{comb} \quad \hat{\theta}_{comb})$ of $\eta^0 = (\mathbf{p}^0 \quad \theta^0)$ wherein $\hat{\mathbf{p}}_{comb}$ (resp. $\hat{\theta}_{comb}$) is therefore also the maximum likelihood estimate of $\mathbf{p}^0$ (resp. $\theta^0$). It is important to notice that the maximum likelihood approach ensures that, among the possible unbiased estimates that can be deduced from $\hat{\theta}_l$ and $\hat{\zeta}_f$, $\hat{\mathbf{p}}_{comb}$ (resp. $\hat{\theta}_{comb}$) will be the minimum variance estimate of $\mathbf{p}^0$ (resp. $\theta^0$). A schematic view of the proposed methodology is presented in Figure 1.
Using the distributions (18) and (11) and the fact that $\hat{\zeta}_f$ and $\hat{\theta}_l$ are independent, we have

$$ \left( \begin{array}{c} \hat{\zeta}_f \\ \hat{\theta}_l \end{array} \right) - \left( \begin{array}{c} T(\mathbf{p}^0)\theta^0 \\ \theta^0 \end{array} \right) \sim \mathcal{N} \left( 0, \underbrace{\begin{pmatrix} P_{\zeta_f} & 0 \\ 0 & P_{\theta,l} \end{pmatrix}}_{=\Sigma} \right). \quad (21) $$

Consequently, the maximum likelihood estimate of $\eta^0$ is obtained by maximizing, over $\eta = (\mathbf{p}, \theta^T)^T \in \mathbb{R}^{k+1}$, the log-likelihood function $L(\eta)$ corresponding to (21). This log-likelihood function $L(\eta)$ is given by

$$ L(\eta) = -V(\eta) - \log \left( \sqrt{(2\pi)^{k+1}} \det(\Sigma) \right), \quad (22) $$

with

$$ \begin{aligned} V(\eta) &= \frac{1}{2} \left( \left( \begin{array}{c} T(\mathbf{p})\theta \\ \theta \end{array} \right) - \left( \begin{array}{c} \hat{\zeta}_f \\ \hat{\theta}_l \end{array} \right) \right)^T \Sigma^{-1} \left( \left( \begin{array}{c} T(\mathbf{p})\theta \\ \theta \end{array} \right) - \left( \begin{array}{c} \hat{\zeta}_f \\ \hat{\theta}_l \end{array} \right) \right) \\ &= \frac{1}{2} \left( \left( T(\mathbf{p})\theta - \hat{\zeta}_f \right)^T P_{\zeta_f}^{-1} \left( T(\mathbf{p})\theta - \hat{\zeta}_f \right) + \left( \theta - \hat{\theta}_l \right)^T P_{\theta,l}^{-1} \left( \theta - \hat{\theta}_l \right) \right). \end{aligned} \quad (23) $$

Observing that the term $\log(\sqrt{(2\pi)^{k+1}} \det(\Sigma))$ is not a function of $\eta$, the maximum likelihood estimate $\hat{\eta}_{comb}$ of $\eta^0$ is also the solution of the following weighted least-squares problem:

$$ \hat{\eta}_{comb} = \arg \min_{\eta} V(\eta) \quad (24) $$

Let us first note that, if we use similar arguments as the ones above, we can show that the estimate $\hat{\theta}_l$, which is the solution of the weighted least-squares problem (15)-(16), is the maximum likelihood estimate of $\theta^0$ given the information (14). As opposed to (15)-(16), the weighted least-squares problem (24) has no closed-form solution. The optimization problem (24) is indeed a nonlinear optimization problem due to the term $T(\mathbf{p})\theta$ in (23) with $\mathbf{p}$ and $\theta$ being both decision variables. In Appendix A, we discuss the methods that can be used to compute its solution $\hat{\eta}_{comb}$.

Recall that the distributions (18) and (11) are only attained asymptotically i.e., when $N_l, N_f \to \infty$. Consequently, $\hat{\eta}_{comb}$ is thus in fact the maximum likelihood estimate of $\eta^0$ only if both $N_l$ and $N_f$ are sufficiently large. We will use the same asymptotic assumption to show in the following proposition that $\hat{\eta}_{comb}$ is normally distributed around $\eta^0$ and to give an expression of its covariance matrix $P_{\eta,comb}$. Note that the expression of this covariance matrix cannot be derived in a similar manner
as the one of $\hat{\theta}_l$ in (Ghosh et al., 2018), since the weighted least-squares problem (24) has no closed-form solution.

**Proposition 3.1.** Consider the estimate $\hat{\eta}_{comb}$ (see (24)) of the true parameter vector $\eta^0$ (see (9)) describing the LPV system (7) operated at the (unknown) operating point $\mathbf{p}^0$. This estimate combines the estimate $\hat{\theta}_l$ for $\theta^0$ (see (18)) and the estimate $\hat{\zeta}_f$ of $\zeta^0(\mathbf{p}^0) = T(\mathbf{p}^0)\theta^0$ (see (11)). Recall also the duration $N_l$ and $N_f$ of the identification experiments yielding these two estimates. Then, asymptotically in $N_l$ and in $N_f$, the estimate $\hat{\eta}_{comb}$ is normally distributed around $\eta^0$ with a covariance matrix $P_{\eta,comb}$ given by:

$$P_{\eta,comb} = (W^T(\eta^0) P_{\zeta_f}^{-1} W(\eta^0) + bdiag(0, P_{\theta,l}^{-1}))^{-1} \quad (25)$$

$$\text{where } W(\eta) = (D(\mathbf{p})\theta \quad T(\mathbf{p})) \text{ with } D(\mathbf{p}) = \frac{dT(\mathbf{p})}{d\mathbf{p}}. \quad (26)$$

An estimate of $P_{\eta,comb}$ can be obtained by replacing, in (25), the unknown $\eta^0$ by $\hat{\eta}_{comb}$.

**Proof.** See Appendix B. □

Proposition 3.1 shows that the covariance matrix $P_{\eta,comb}$ can be easily estimated using $\hat{\eta}_{comb}$. Consequently, as was also the case for $P_{\zeta_f}$ and $P_{\theta,l}$, this estimate of $P_{\eta,comb}$ can be used as a measure of the uncertainty of the estimate $\hat{\eta}_{comb}$ (and this even in the case where $N_l$ and $N_f$ are finite). As shown in (25), $P_{\eta,comb}$ depends both on $P_{\zeta_f}$ and $P_{\theta,l}$ which are measures of the uncertainty of, respectively, $\hat{\zeta}_f$ and $\hat{\theta}_l$. More precisely, a smaller $P_{\eta,comb}$ (i.e., a smaller uncertainty on $\hat{\eta}_{comb}$) can be obtained by decreasing $P_{\zeta_f}$ and/or by decreasing $P_{\theta,l}$. The latter can, e.g., be achieved by increasing the duration of the identification experiments yielding $\hat{\zeta}_f$ and $\hat{\theta}_l$. It can also be achieved by increasing the power and/or tuning the power spectrum of the excitation signals involved in these experiments.

When $N_l$ and $N_f$ are finite, the matrix $P_{\eta,comb}$ given in (25) is always a strictly positive-definite matrix (see Appendix D). As also shown in this appendix, when both $N_f$ and $N_l$ converge to infinity, the covariance matrix $P_{\eta,comb}$ converges to the zero matrix. This is not surprising since, in that case, the estimates $\hat{\theta}_l$ and $\hat{\zeta}_f$ are equal to $\theta^0$ and $\zeta^0(\mathbf{p}^0) = T(\mathbf{p}^0)\theta^0$, respectively, and $\eta^0$ can thus be exactly retrieved. On the contrary, when $N_l$ remains finite and $N_f \to \infty$, $P_{\eta}$ will be a positive semi-definite matrix different of the zero matrix. This fact is also not surprising since $\eta^0$ cannot be exactly retrieved from $\zeta^0(\mathbf{p}^0) = T(\mathbf{p}^0)\theta^0$ and an uncertain estimate $\hat{\theta}_l$ of $\theta^0$. Obviously, the same can be said in the case where $N_f$ remains finite and $N_l \to \infty$. More details can be found in Appendix D.

**Remark.** Given the information $\hat{\theta}_l$ and $\hat{\zeta}_f$, the smallest covariance matrix for an unbiased estimator of $\eta^0$ can be related to the log-likelihood function $L(\eta)$ given in (22). More precisely, it is equal to the inverse of $J$ with:

$$J \triangleq E \left( - \frac{d^2 L(\eta)}{d\eta^2} \bigg|_{\eta=\eta^0} \right) = E \left( \frac{d^2 V(\eta)}{d\eta^2} \bigg|_{\eta=\eta^0} \right) = E(\mathcal{H}(\eta_0)), \quad (27)$$
where $\mathcal{H}(\eta)$ is the Hessian of $V(\eta)$. Using the expression of $\mathcal{H}(\eta)$ given in Appendix A.2, we find $\mathcal{H}(\eta_0) = P_{\eta,comb}^{-1} + bdiag(A(T(\mathbf{p}^0)\theta^0 - \hat{\zeta}_f), 0)$ with $A$ a given deterministic matrix. Using now (11), it is clear that $E(\mathcal{H}(\eta_0)) = P_{\eta,comb}^{-1}$. The observation that $P_{\eta,comb} = J^{-1}$ is not a surprising observation since $\hat{\eta}_{comb}$ is the maximum likelihood estimate (and thus the minimum variance estimate).

## 4. Accuracy of $\hat{\mathbf{p}}_{comb}$ and of $\hat{\theta}_{comb}$

Let us summarize the results. We suppose that we have an initial estimate $\hat{\theta}_l$ of $\theta^0$ distributed as (18) and an estimate $\hat{\zeta}_f$ of $\zeta^0(\mathbf{p}^0) = T(\mathbf{p}^0)\theta_0$ distributed as (11). This information is combined via the criterion (24) to yield an unbiased estimate $\hat{\eta}_{comb} \triangleq \begin{pmatrix} \hat{\mathbf{p}}_{comb} \\ \hat{\theta}_{comb} \end{pmatrix}$ of $\eta^0 = \begin{pmatrix} \mathbf{p}^0 \\ \theta^0 \end{pmatrix}$ with a covariance matrix $P_{\eta,comb}$ given in (25).

By decomposing the vector $\hat{\eta}_{comb} \triangleq \begin{pmatrix} \hat{\mathbf{p}}_{comb} \\ \hat{\theta}_{comb} \end{pmatrix}$, we obtain an estimate $\hat{\mathbf{p}}_{comb}$ of the unknown operating point $\mathbf{p}^0$ and we also obtain a novel estimate $\hat{\theta}_{comb}$ of $\theta^0$. In this section, we will analyze, in details, the respective uncertainty of these estimates $\hat{\mathbf{p}}_{comb}$ and $\hat{\theta}_{comb}$ by investigating the structure of (25).

Let us for this purpose decompose $P_{\eta,comb}$ as follows:

$$ P_{\eta,comb} = \begin{pmatrix} P_{p,comb} & P_{p\theta,comb} \\ P_{p\theta,comb}^T & P_{\theta,comb} \end{pmatrix}, \qquad (28) $$

with $P_{p,comb} \in \mathbb{R}$ and $P_{\theta,comb} \in \mathbb{R}^{k \times k}$. Since $\hat{\eta}_{comb} \sim \mathcal{N}(\eta^0, P_{\eta,comb})$, the estimates $\hat{\mathbf{p}}_{comb}$ and $\hat{\theta}_{comb}$ are therefore distributed as follows:

$$ \hat{\mathbf{p}}_{comb} \sim \mathcal{N}(\mathbf{p}^0, P_{p,comb}), \qquad (29) $$

$$ \hat{\theta}_{comb} \sim \mathcal{N}(\theta^0, P_{\theta,comb}), \qquad (30) $$

which show that $P_{p,comb}$ is a measure of the uncertainty of $\hat{\mathbf{p}}_{comb}$ and $P_{\theta,comb}$ a measure of the uncertainty of $\hat{\theta}_{comb}$. Using the matrix inversion formula (Zhou & Doyle, 1998), we can derive the following expressions for $P_{p,comb}$ and $P_{\theta,comb}$:

$$ P_{p,comb} = (\mathcal{R}_{\eta,11} - \mathcal{R}_{\eta,12} \mathcal{R}_{\eta,22}^{-1} \mathcal{R}_{\eta,12}^T)^{-1} \qquad (31) $$

$$ P_{\theta,comb} = (\mathcal{R}_{\eta,22} - \mathcal{R}_{\eta,12}^T \mathcal{R}_{\eta,11}^{-1} \mathcal{R}_{\eta,12})^{-1} \qquad (32) $$
where $\mathcal{R}_{\eta,ij}$ ($i = 1,2$ and $j = 1,2$) are the block-elements of the inverse $P_{\eta,comb}^{-1}$ of the symmetric matrix $P_{\eta,comb}$. More precisely, $\mathcal{R}_{\eta,11} = (D(\mathbf{p}^0)\theta^0)^T P_{\zeta_f}^{-1} (D(\mathbf{p}^0)\theta^0)$, $\mathcal{R}_{\eta,12} = (D(\mathbf{p}^0)\theta^0)^T P_{\zeta_f}^{-1} T(\mathbf{p}^0)$ and $\mathcal{R}_{\eta,22} = T^T(\mathbf{p}^0) P_{\zeta_f}^{-1} T(\mathbf{p}^0) + P_{\theta,l}^{-1}$.

Since, in our framework, we have already an estimate $\hat{\theta}_l$ of $\theta^0$, we are more interested in the accuracy of the estimate $\hat{\mathbf{p}}_{comb}$ of $\mathbf{p}^0$ than in the one of the novel estimate $\hat{\theta}_{comb}$ of $\theta^0$. Before analyzing in details the accuracy of the estimate $\hat{\mathbf{p}}_{comb}$ of $\mathbf{p}^0$, let us thus just give the following relation about the covariance matrix $P_{\theta,comb}$ of $\hat{\theta}_{comb}$ (see (32)):

$$0 \le P_{\theta,comb} \le P_{\theta,l} \qquad (33)$$

Consequently, if the initial estimate $\hat{\theta}_l$ has a nonzero<fnref n="2" /> uncertainty $P_{\theta,l} > 0$, the estimate $\hat{\theta}_{comb}$ derived via (24) is a more accurate estimate than the initial estimate $\hat{\theta}_l$. This is not surprising since $\hat{\theta}_{comb}$ combines the information about $\theta^0$ in both $\hat{\theta}_l$ and $\hat{\zeta}_f$. See Appendix E for a formal proof of (33).

In the following proposition, we now analyze in details the variance $P_{p,comb}$ of $\hat{\mathbf{p}}_{comb}$ (see (31)).

**Proposition 4.1.** Consider the estimate $\hat{\eta}_{comb} \triangleq \begin{pmatrix} \hat{\mathbf{p}}_{comb} \\ \hat{\theta}_{comb} \end{pmatrix}$ of $\eta^0 = \begin{pmatrix} \mathbf{p}^0 \\ \theta^0 \end{pmatrix}$ with a covariance matrix $P_{\eta,comb}$ whose expression as a function of the covariance matrices $P_{\theta,l}$ of $\hat{\theta}_l$ and $P_{\zeta_f}$ of $\hat{\zeta}_f$ is given in (25). Consider also the variance $P_{p,comb}$ of $\hat{\mathbf{p}}_{comb}$ (see (31)). Define

$$\alpha_{min,1}(P_{\theta,l}) = \frac{1}{(D(\mathbf{p}^0)\theta^0)^T (T(\mathbf{p}^0)P_{\theta,l}T^T(\mathbf{p}^0))^{-1} D(\mathbf{p}^0)\theta^0} \ge 0 \qquad (34)$$

$$\alpha_{min,2}(P_{\zeta_f}) = \frac{1}{(D(\mathbf{p}^0)\theta^0)^T P_{\zeta_f}^{-1} D(\mathbf{p}^0)\theta^0} \ge 0 \qquad (35)$$

and note that $\alpha_{min,1}(P_{\theta,l})$ (resp. $\alpha_{min,2}(P_{\zeta_f})$) is a scalar quantity depending on $P_{\theta,l}$ (resp. $P_{\zeta_f}$). Define also $P_{p,des}$ as a desired value for $P_{p,comb}$.

Then, the variance $P_{p,comb}$ of $\hat{\mathbf{p}}_{comb}$ has the following properties:

i. Consider the case where the estimate $\hat{\theta}_l$ and its nonzero uncertainty $P_{\theta,l} > 0$ are given. In that particular case, whatever the uncertainty $P_{\zeta_f}$ of $\hat{\zeta}_f$ (i.e., whatever the experimental conditions yielding $\hat{\zeta}_f$), the variance $P_{p,comb}$ of $\hat{\mathbf{p}}_{comb}$ will always be larger than (or equal to) the quantity $\alpha_{min,1}(P_{\theta,l})$ i.e., $P_{p,comb} \ge \alpha_{min,1}(P_{\theta,l})$. Moreover, $P_{p,comb} = \alpha_{min,1}(P_{\theta,l})$ only when $P_{\zeta_f} \to 0$ (i.e., when $N_f \to \infty$). Observe also that the scalar quantity $\alpha_{min,1}(P_{\theta,l})$

<footnote n="2">²If the covariance matrix $P_{\theta,l}$ of the initial estimate $\hat{\theta}_l$ is already equal to the zero matrix, the relation (33) shows that the covariance matrix $P_{\theta,comb}$ is also equal to the zero matrix.</footnote>
increases when $P_{\theta,l}$ increases. Consequently, we can only achieve the desired variance $P_{p,des}$ if the given uncertainty $P_{\theta,l}$ of $\hat{\theta}_l$ is sufficiently small to ensure $\alpha_{min,1}(P_{\theta,l}) < P_{p,des}$.

ii. Consider the case where the estimate $\hat{\zeta}_f$ and its nonzero uncertainty $P_{\zeta_f} > 0$ are given. In that particular case, whatever the uncertainty $P_{\theta,l}$ of $\hat{\theta}_l$ (i.e., whatever the experimental conditions yielding $\hat{\theta}_l$), the variance $P_{p,comb}$ of $\hat{\mathbf{p}}_{comb}$ will always be larger than (or equal to) the quantity $\alpha_{min,2}(P_{\zeta_f})$ i.e., $P_{p,comb} \ge \alpha_{min,2}(P_{\zeta_f})$. Moreover, $P_{p,comb} = \alpha_{min,2}(P_{\zeta_f})$ only when $P_{\theta,l} \to 0$ (i.e., when $N_l \to \infty$). Observe also that the scalar quantity $\alpha_{min,2}(P_{\zeta_f})$ increases when $P_{\zeta_f}$ increases. Consequently, we can only achieve the desired variance $P_{p,des}$ if the given uncertainty $P_{\zeta_f}$ of $\hat{\zeta}_f$ is sufficiently small to ensure $\alpha_{min,2}(P_{\zeta_f}) < P_{p,des}$.

iii. When both $P_{\theta,l} \ne 0$ and $P_{\zeta_f} \ne 0$, we have that $P_{p,comb} > \alpha_{min,1}(P_{\theta,l})$ and $P_{p,comb} > \alpha_{min,2}(P_{\zeta_f})$.

iv. $P_{p,comb} \to 0$ only when both $P_{\theta,l} \to 0$ and $P_{\zeta_f} \to 0$ (i.e., when both $N_l, N_f \to \infty$)

*Proof.* See Appendix F. $\square$

From this proposition, it is clear that the uncertainty of $\hat{\zeta}_f$ and the uncertainty of $\hat{\theta}_l$ limits the achievable accuracy of $\hat{\mathbf{p}}_{comb}$.

# 5. Optimal experiment design to tune the variance $P_{p,comb}$ of $\hat{\mathbf{p}}_{comb}$

## 5.1. Optimal experiment design as a convex optimization problem

In the previous section, we have discussed the statistical properties of the estimate $\hat{\eta}_{comb}$ and, in particular, of its first entry. This first entry is the estimate $\hat{\mathbf{p}}_{comb}$ of the unknown operating point $\mathbf{p}^0$. We have observed that the accuracy of $\hat{\mathbf{p}}_{comb}$ depends on the covariance matrices $P_{\theta,l}$ and $P_{\zeta_f}$ of $\hat{\theta}_l$ and $\hat{\zeta}_f$. Since both $\hat{\theta}_l$ and $\hat{\zeta}_f$ follow from an identification procedure, it is possible to tune their covariance matrices $P_{\theta,l}$ and $P_{\zeta_f}$ in order to tune the accuracy of $\hat{\mathbf{p}}_{comb}$. Since the optimal design of an identification experiment always requires some information about the to-be-identified system (Bombois, Scorletti, Gevers, Van den Hof, & Hildebrand, 2006; Jansson & Hjalmarsson, 2005), we will assume in the sequel that the initial local LPV procedure yielding $\hat{\theta}_l$ has been performed<fnref n="3" /> and we will therefore only design in an optimal way the identification experiment performed on the LPV system operated at $\mathbf{p}^0$ i.e., the identification experiment yielding $\hat{\zeta}_f$. In other words, we will optimally design the covariance matrix $P_{\zeta_f} = \frac{1}{N_f} \mathcal{I}_{\zeta_f}$ with $\mathcal{I}_{\zeta_f}$ given in (13). We will do that with the objective of obtaining the smallest possible variance $P_{p,comb}$ for $\hat{\mathbf{p}}_{comb}$. Having a small variance $P_{p,comb}$ is indeed crucial since a possible use of this estimate $\hat{\mathbf{p}}_{comb}$ would be to replace the unknown $\mathbf{p}^0$ in an LPV controller $K(\mathbf{p})$ designed based on $\hat{\theta}_l$ and it is obvious that, in such a situation, the estimate $\hat{\mathbf{p}}_{comb}$ of $\mathbf{p}^0$ must be accurate to guarantee the control performance of $K(\hat{\mathbf{p}}_{comb})$ on the LPV system operated at $\mathbf{p}^0$.

Recall that, as shown in (28), $P_{p,comb}$ is the (1,1)-entry of the covariance matrix $P_{\eta,comb}$ of $\hat{\eta}_{comb}$. As shown in (25), $P_{\eta,comb}$ depends on a number of variables i.e., on

<footnote n="3">³The optimal design of the local LPV identification procedure yielding $\hat{\theta}_l$ is described in (Ghosh et al., 2018).</footnote>
the fixed $P_{\theta,l}$, on the true parameter vector $\eta_0 = (\mathbf{p}^0, (\theta^0)^T)^T$ and on the covariance matrix $P_{\zeta_f}$ that we wish to design. Given (12) and (13), we see that $P_{\zeta_f}$ is function of $\eta^0$ and of the experimental conditions under which $\hat{\zeta}_f$ will be identified (i.e., the duration $N_f$ of the experiment and the spectrum $\Phi_{u_f}$ of the excitation signal $u_f$). Since $P_{p,comb}$ is the (1,1)-entry of $P_{\eta,comb}$, it is clear that $P_{p,comb}$ depends on the same variables. In order to stress this dependence, we will use the following notations for $P_{\eta,comb}$ and $P_{p,comb}$ in the sequel: $P_{\eta,comb}(\Phi_{u_f}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l})$ and $P_{p,comb}(\Phi_{u_f}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l})$.

Since we suppose that $\hat{\theta}_l$ and its covariance matrix $P_{\theta,l}$ are fixed, we know from Proposition 4.1 that the minimal achievable variance $P_{p,comb}(\Phi_{u_f}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l})$ for $\hat{\mathbf{p}}_{comb}$ is given by $\alpha_{min,1}(P_{\theta,l})$. In practice, this minimal variance $\alpha_{min,1}(P_{\theta,l})$ will not be achievable due the inherent limitations on the duration $N_f$ of the experiment and on the power of the excitation signal $u_f$ ($P_{\zeta_f}$ will indeed be different from zero due to these limitations). The objective of this section will therefore be to take into account these inherent limitations on $N_f$ and $u_f$ and to determine the experimental conditions yielding the smallest value of $P_{p,comb}(\Phi_{u_f}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l})$. Here, these limitations will be that $N_f$ is fixed and that the power of $u_f$ cannot be larger than a given threshold $P_{u,max}$. The optimal experiment design problem can therefore be formulated as follows:

$$ \arg \min_{\Phi_{u_f}, \gamma>0} \gamma \qquad (36) $$

$$ \text{subject to } P_{p,comb}(\Phi_{u_f}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l}) \leq \gamma \qquad (37) $$

$$ \text{and to } \frac{1}{2\pi} \int_{-\pi}^{\pi} \Phi_{u_f}(\omega) \, d\omega \leq P_{u,max} \qquad (38) $$

In order to convexify this optimization problem (where (37) is not affine in the decision variable $\Phi_{u_f}$), we will replace (37) by an equivalent constraint linear in the inverse of $P_{\eta,comb}(\Phi_{u_f}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l})$. From (25), (12) and (13), it is indeed clear that the inverse of $P_{\eta,comb}$ is an affine function of $\Phi_{u_f}$. This equivalent reformulation is given in the following proposition.

**Proposition 5.1.** Consider the covariance matrix $P_{\eta,comb}(\Phi_{u_f}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l}) > 0$ and its (1,1)-entry $P_{p,comb}(\Phi_{u_f}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l})$ (see (25) and (28)). Then, the solution $\gamma^{opt}$ of the optimization problem (36)-(38) is equal to $\frac{1}{\beta^{opt}}$ where $\beta^{opt}$ is the solution of the following convex optimization problem:

$$ \arg \max_{\Phi_{u_f}, \beta>0} \beta \qquad (39) $$

subject to (38) and to

$$ P_{\eta,comb}^{-1}(\Phi_{u_f}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l}) \geq \beta e_1 e_1^T \qquad (40) $$

where $e_1 = (1, 0, ..., 0)^T \in \mathbb{R}^{k+1}$. Moreover, the solution $\Phi_{u_f}^{opt}$ of (39)-(40) is also the
solution of (36)-(38).

**Proof.** Let us first note that, due to (28), we have that

$$P_{p,comb}(\Phi_{u_f}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l}) = e_1^T P_{\eta,comb}(\Phi_{u_f}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l}) e_1 \quad (41)$$

with the vector $e_1$ defined in the statement of the theorem. Using the facts that $P_{\eta,comb}(\Phi_{u_f}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l}) > 0$ and that $\gamma > 0$, we can use the lemma in Appendix G to show that the initial accuracy constraint (37) is equivalent to (40) with $\beta = \frac{1}{\gamma}$. Using the change of variables $\beta = \frac{1}{\gamma}$, it is thus clear that the optimization problem (39)-(40) is equivalent to the optimization problem (36)-(38).

The optimization problem (39)-(40) is moreover convex since the objective function and the two constraints of the optimization problem (39)-(40) are affine in the decision variables $\beta$ and $\Phi_{u_f}$ ($P_{\eta,comb}^{-1}(\Phi_{u_f}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l})$ being affine in $\Phi_{u_f}$). $\square$

If we choose a linear parametrization for the to-be-determined spectrum $\Phi_{u_f}$, the convex optimization problem (39)-(40) can be efficiently solved using LMI optimization (Boyd, El Ghaoui, Feron, & Balakrishnan, 1994). Here, we will use the following parametrization (Jansson & Hjalmarsson, 2005) that corresponds to the spectrum of a signal that is generated by a white noise filtered by an arbitrary FIR filter of degree $L$ ($L$ is a user-chosen parameter):

$$\Phi_{u_f}(\omega) = c_0 + 2 \sum_{i=1}^{L} c_i \cos(i\omega). \quad (42)$$

The positivity of (42) for all $\omega$ can be imposed by a Linear Matrix Inequality (LMI) constraint on the coefficients $c_i$ ($i = 0, \dots, L$) (Jansson & Hjalmarsson, 2005). These coefficients will become the decision variables of the optimization problem (39)-(40).

Solving the optimization problem (39)-(40) for fixed $N_f$ and $P_{\theta,l}$ yields $\beta^{opt}$ and $\Phi_{u_f}^{opt}$. According to Proposition 5.1, if we identify $\hat{\zeta}_f$ with an identification experiment of duration $N_f$ and an excitation spectrum $\Phi_{u_f}^{opt}$, the estimate $\hat{\mathbf{p}}_{comb}$ deduced based on that $\hat{\zeta}_f$ and the fixed $\hat{\theta}_l$ will then have a variance:

$$P_{p,comb}^{opt} \triangleq P_{p,comb}(\Phi_{u_f}^{opt}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l}) = \frac{1}{\beta^{opt}}$$

and this variance $P_{p,comb}^{opt}$ is the smallest variance for $\hat{\mathbf{p}}_{comb}$ that can be achieved given the fixed duration $N_f$ and the excitation power constraint (38). Note also that the optimization problem can of course be solved for different acceptable values of $N_f$ in order to fine tune the optimal variance $P_{p,comb}^{opt}$: the larger $N_f$, the smaller the optimal variance $P_{p,comb}^{opt}$ will be.
## 5.2. Dealing with the chicken-and-egg problem

In the optimal experiment design problem (39)-(40), the accuracy constraint (40) depends on the unknown quantities $\mathbf{p}^0$ and $\theta^0$. This is an occurrence of the so-called chicken and egg problem encountered in almost all optimal experiment design problems (Bombois et al., 2006; Jansson & Hjalmarsson, 2005). The chicken and egg problem is generally circumvented by replacing the unknown quantities by initial estimates for these quantities (Bombois et al., 2006; Jansson & Hjalmarsson, 2005). In the situation considered in this paper, we have an initial estimate $\hat{\theta}_l$ of $\theta^0$. However, before performing the identification yielding $\hat{\zeta}_f$, the only information about the unknown operating point $\mathbf{p}^0$ is that it lies in the interval $\mathcal{P}$ in which $\mathbf{p}$ can vary. We will therefore aim at determining a power spectrum $\Phi_{uf}$ which will lead to an acceptable $P_{\eta,comb}$ for all possible operating points $\mathbf{p} \in \mathcal{P}$. This leads to the following robust optimal experiment design problem:

$$ \arg\max_{\Phi_{uf}, \beta>0} \beta \qquad (43) $$

subject to (38) and to

$$ P_{\eta,comb}^{-1}(\Phi_{uf}, N_f, \mathbf{p}, \hat{\theta}_l, P_{\theta,l}) \geq \beta e_1 e_1^T \quad \forall \mathbf{p} \in \mathcal{P} \qquad (44) $$

Note that, in this robust optimal experiment design problem, the unknown $\theta^0$ has been replaced by its available initial estimate $\hat{\theta}_l$ and that, instead of one unique accuracy constraint, we have an accuracy constraint for all $\mathbf{p} \in \mathcal{P}$ i.e., we have an infinite number of accuracy constraints (see (44)). Due to their inherent complexity, robust optimization problem such as (43)-(44) are commonly addressed using the so-called gridding approach (Rojas, Welsh, Goodwin, & Feuer, 2007) where the compact set $\mathcal{P}$ is replaced by a gridded version of $\mathcal{P}$. In our case, this gridding approach is easy to implement since $\mathcal{P}$ is a scalar interval. Let thus replace (43)-(44) where we have a constraint for each $\mathbf{p} \in \mathcal{P}$, by the optimization problem (45)-(46) where we have a constraint for $n_{grid}$ points $\mathbf{p}_j$ ($j = 1, ..., n_{grid}$) in $\mathcal{P}$:

$$ \arg\max_{\Phi_{uf}, \beta>0} \beta \qquad (45) $$

subject to (38) and to

$$ P_{\eta,comb}^{-1}(\Phi_{uf}, N_f, \mathbf{p}_j, \hat{\theta}_l, P_{\theta,l}) \geq \beta e_1 e_1^T \quad j = 1, ..., n_{grid} \qquad (46) $$

Let us assume that the $n_{grid}$ points $\mathbf{p}_j$ sufficiently cover<fnref n="4" /> the interval $\mathcal{P}$, then the solution $\beta^{opt}$ and $\Phi_{uf}^{opt}$ of (45)-(46), will have the following property:

$$ P_{p,comb}(\Phi_{uf}^{opt}, N_f, \mathbf{p}, \hat{\theta}_l, P_{\theta,l}) \geq \frac{1}{\beta^{opt}} \quad \forall \mathbf{p} \in \mathcal{P}. $$

<footnote n="4">⁴By using the solution $\beta^{opt}$ and $\Phi_{uf}^{opt}$ of (45)-(46), we can easily check if the $n_{grid}$ points $\mathbf{p}_j$ sufficiently cover the interval $\mathcal{P}$. Indeed, we can compute $P_{\eta,comb}(\Phi_{uf}^{opt}, N_f, \bar{\mathbf{p}}, \hat{\theta}_l, P_{\theta,l})$ for all values $\bar{\mathbf{p}}$ in an extremely fine grid of $\mathcal{P}$ (i.e., a grid containing (much) more grid points than $n_{grid}$) and verify whether it is satisfied that, for all these $\bar{\mathbf{p}}$, $P_{\eta,comb}(\Phi_{uf}^{opt}, N_f, \bar{\mathbf{p}}, \hat{\theta}_l, P_{\theta,l}) \geq \beta^{opt} e_1 e_1^T$. In the eventuality that this constraint is not satisfied for some of these $\bar{\mathbf{p}}$, then (45)-(46) should be re-solved using a grid containing the initial $n_{grid}$ points and the $\bar{\mathbf{p}}$ for which the constraint was not satisfied.</footnote>
For the unknown $\mathbf{p}^0$, we thus also have:

$$P_{p,comb}(\Phi_{u_f}^{opt}, N_f, \mathbf{p}^0, \hat{\theta}_l, P_{\theta,l}) \geq \frac{1}{\beta^{opt}}.$$

The above expressions of the variance $P_{p,comb}$ are evaluated at the estimate $\hat{\theta}_l$ of $\theta^0$. If we, in addition, assume that

$$P_{p,comb}(\Phi_{u_f}^{opt}, N_f, \mathbf{p}^0, \hat{\theta}_l, P_{\theta,l}) \approx P_{p,comb}(\Phi_{u_f}^{opt}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l}),$$

we can then say that the estimate $\hat{\zeta}_f$ of $\zeta^0(\mathbf{p}^0) = T(\mathbf{p}^0)\theta^0$ obtained using an experiment of duration $N_f$ and with an excitation spectrum $\Phi_{u_f}^{opt}$ on the LPV system operated at $\mathbf{p}^0$ will yield, when combined with $\hat{\theta}_l$, an estimate $\hat{p}_{comb}$ of the unknown $\mathbf{p}^0$ having a variance $P_{p,comb}(\Phi_{u_f}^{opt}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l})$ which will be most probably larger than $\frac{1}{\beta^{opt}}$. Here also, the obtained variance can be tuned by solving the optimization problem for different values of acceptable $N_f$.

**Remark.** The identification experiment with excitation signal $u_f$ which yields the estimate $\hat{\zeta}_f$ of $\zeta^0(\mathbf{p}^0)$ has to be performed at the beginning of each batch in the operation phase i.e., when the value of $\mathbf{p}^0$ changes (see the discussion in Section 2). Note that an advantage of the robust optimal experiment problem (45)-(46) is that the design of the optimal spectrum $\Phi_{u_f}^{opt}$ of the excitation signal $u_f$ has only to be performed once. Indeed, this spectrum is appropriate whatever the value of the scheduling variable $\mathbf{p}^0$ may be.

## 6. Numerical illustration

We consider a true LPV system which can be represented by the following first order ARX system at each operating point $\mathbf{p} \in \mathcal{P} = [0 \ 3]$:

$$y(t) = \frac{b^0(\mathbf{p})z^{-1}}{1 + a^0(\mathbf{p})z^{-1}}u(t) + \frac{1}{1 + a^0(\mathbf{p})z^{-1}}e(t) \quad (47)$$

with $e(t)$ a Gaussian white noise of variance $\sigma_e^2 = 0.01$ and with:

$$\zeta^0(\mathbf{p}) \triangleq \begin{pmatrix} a^0(\mathbf{p}) \\ b^0(\mathbf{p}) \end{pmatrix} = \underbrace{\begin{pmatrix} 1 & \mathbf{p} & 0 & 0 & 0 \\ 0 & 0 & 1 & \mathbf{p} & \mathbf{p}^2 \end{pmatrix}}_{=T(\mathbf{p})} \underbrace{\begin{pmatrix} -0.5 \\ -0.1 \\ 1 \\ 0.2 \\ 0.05 \end{pmatrix}}_{=\theta^0} \quad (48)$$

In this description, we have<fnref n="5" /> thus $n = 2$ and $k = 5$. We will also suppose in the sequel

<footnote n="5">⁵Here, the degree $n_p$ (see (2)) of the polynomial expansion in $\mathbf{p}$ is different for $a^0(\mathbf{p})$ and for $b^0(\mathbf{p})$: $n_p = 1$ for $a^0(\mathbf{p})$ and $n_p = 2$ for $b^0(\mathbf{p})$. The resulting matrix $T(\mathbf{p})$ is nevertheless also full row rank.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2.: Magnitude plot of $G(e^{j\omega}, \zeta^0(\mathbf{p}))$ for $\mathbf{p} = 0, 0.5, 1, 1.5, 2, 2.5$ and $3$. The static gain of $G(z, \zeta^0(\mathbf{p}))$ is growing for increasing value of $\mathbf{p}$. The orange curve (the second from the bottom) thus corresponds to $\mathbf{p} = 0.5$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3.: Cumulative distribution $f_X(r) = \Pr(\chi^2(6) \le r)$ of a $\chi^2$-distribution with 6 degrees of freedom (red curve) and sampled cumulative distribution of $d(\hat{\eta}_{comb}, \eta^0) \triangleq (\hat{\eta}_{comb} - \eta^0)^T P_{\eta,comb}^{-1} (\hat{\eta}_{comb} - \eta^0)$ following from 1000 Monte-Carlo simulations (blue curve).</figcaption></figure>
that the maximal power of the excitation signal $u(t)$ is $P_{u,max} = 0.1$. The magnitude plots of $G(z, \zeta^0(\mathbf{p}))$ for different $\mathbf{p} \in \mathcal{P}$ are represented in Figure 2. Let us consider this LPV system at the operating point $\mathbf{p}^0 = 0.5$ where $\zeta^0(\mathbf{p}^0) = (-0.55, 1.1125)^T$ and whose magnitude plot is represented in orange in Figure 2 (the curve with a static gain of 8 dB). We will use criterion (24) to determine an estimate $\hat{\eta}_{comb}$ of the parameter vector $\eta^0 = (\mathbf{p}^0 (\theta^0)^T)^T$ describing this system. Let us first illustrate that the distribution $\hat{\eta}_{comb} \sim \mathcal{N}(\eta^0, P_{\eta,comb})$ which has been proven for $N_f, N_l \to \infty$ in Proposition 3.1 can be also used for finite $N_l$ and $N_f$. For this purpose, we will repeat the following procedure 1000 times with different realization of the noise perturbation (Monte-Carlo simulations).

**Procedure for the Monte-Carlo simulations.** An initial estimate $\hat{\theta}_l$ of $\theta^0$ as well as its variance matrix $P_{\theta,l}$ is obtained using the technique described in Section 2.3. More precisely, we perform local identification experiments on the LPV system (47) at $M=3$ different operating points in $\mathcal{P}$ (i.e. $\mathbf{p}_1 = 0.1, \mathbf{p}_2 = 2$ and $\mathbf{p}_3 = 3$) using Gaussian white noise excitation signals $u_i(t)$ ($i=1,...,3$) of variance 0.1 and of duration $N_l = 100$. We then perform an identification experiment of duration $N_f = 100$ on the LPV system operated at $\mathbf{p}^0 = 0.5$ with a Gaussian white noise excitation $u_f(t)$ of variance 0.1 yielding an estimate $\hat{\zeta}_f$ of $\zeta^0(\mathbf{p}^0)$ and its covariance matrix $P_{\zeta_f}$. Using $\hat{\zeta}_f$ and $\hat{\theta}_l$, we determine $\hat{\eta}_{comb} \in \mathbb{R}^6$ by solving (24) and we compute the covariance matrix $P_{\eta,comb}$ as indicated in Proposition 3.1 i.e., using (25) in which we replace the unknown $W(\eta^0)$ by $W(\hat{\eta}_{comb})$. Using $\hat{\eta}_{comb}$ and $P_{\eta,comb}$, we finally compute the scalar quantity $d(\hat{\eta}_{comb}, \eta^0) \triangleq (\hat{\eta}_{comb} - \eta^0)^T P_{\eta,comb}^{-1} (\hat{\eta}_{comb} - \eta^0)$.

**Verification of the statistical properties.** If $\hat{\eta}_{comb} \sim \mathcal{N}(\eta^0, P_{\eta,comb})$, then we also have that $d(\hat{\eta}_{comb}, \eta^0) \triangleq (\hat{\eta}_{comb} - \eta^0)^T P_{\eta,comb}^{-1} (\hat{\eta}_{comb} - \eta^0)$ is distributed as a chi-squared distribution $\chi^2(k+1)$ with $k+1=6$ degrees of freedom ($k+1=6$ is the dimension of $\eta^0$). Since we have repeated the above procedure 1000 times, we have 1000 realizations of $d(\hat{\eta}_{comb}, \eta^0)$ and let us verify if these 1000 realizations follows a $\chi^2$-distribution with 6 degrees of freedom. For this purpose, let us compare, in Figure 3, the curve $f_\chi(r) = \text{Pr}(\chi^2(6) \le r)$ (i.e. the so-called cumulative distribution) and the curve $f_{mc}(r)$ defined as the proportion of the 1000 realizations of $d(\hat{\eta}_{comb}, \eta^0)$ which satisfies $d(\hat{\eta}_{comb}, \eta^0) \le r$. As shown in Figure 3, these two curves are very close to each other showing that $\hat{\eta}_{comb}$ is close to be distributed as $\mathcal{N}(\eta^0, P_{\eta,comb})$ even if $N_l = N_f = 100 \ll \infty$.

**Variance** $P_{p,comb}(\Phi_{u_f}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l})$ **of the estimate** $\hat{\mathbf{p}}_{comb}$ **of p**$^0$. An initial estimate $\hat{\theta}_l$ of $\theta^0$ as well as its variance matrix $P_{\theta,l}$ is now obtained using a local LPV procedure on the LPV system in (47) at the same three different operating points as above (i.e. $\mathbf{p}_1 = 0.1, \mathbf{p}_2 = 2$ and $\mathbf{p}_3 = 3$). The three local identification experiments have been performed with Gaussian white noise excitation signals $u_i(t)$ ($i=1,...,3$) of variance 0.1 and of duration $N_l = 5000$. The obtained $\hat{\theta}_l$ is $\hat{\theta}_l = (-0.5094, -0.0968, 1.0037, 0.2003, 0.0490)^T$ and the covariance matrix $P_{\theta,l}$ is
$$ P_{\theta,l} = 10^{-4} \begin{pmatrix} 0.0992 & -0.0359 & 0.0027 & 0.0003 & -0.0005 \\ -0.0359 & 0.0144 & -0.0010 & -0.0001 & 0.0002 \\ 0.0027 & -0.0010 & 0.2440 & -0.2584 & 0.0605 \\ 0.0003 & -0.0001 & -0.2584 & 0.7868 & -0.2492 \\ -0.0005 & 0.0002 & 0.0605 & -0.2492 & 0.0838 \end{pmatrix} \quad (49) $$

Let us now suppose that the experiment yielding $\hat{\zeta}_f$ will have a duration $N_f = 100$ and that the excitation signal will be a white noise of variance equal to $P_{u,max} = 0.1$ i.e., $\Phi_{uf}(\omega) = 0.1 \forall \omega$. Using this information, we can compute $P_{\zeta_f}$ using (13) and (12) and we can then subsequently compute $P_{p,comb}(\Phi_{uf}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l})$ using (25), (28) and (49). This yields here $P_{p,comb}(\Phi_{uf}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l}) = 0.0124$. In order to assess this variance, we also compute the lower bounds on the variance given in Proposition 4.1: $\alpha_{min,1}(P_{\theta,l}) = 1.78 10^{-4}$ and $\alpha_{min,2}(P_{\zeta_f}) = 0.0122$. We observe that $P_{p,comb}(\Phi_{uf}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l}) = 0.0124$ is indeed larger than these two lower bounds. Moreover, since $P_{p,comb}(\Phi_{uf}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l}) = 0.0124$ is very close to $\alpha_{min,2}(P_{\zeta_f}) = 0.0122$ for the obtained $P_{\zeta_f}$, there is very little to gain in further reducing the covariance matrix $P_{\theta,l}$. Unlike this, by inspecting $P_{p,comb} = 0.0124$ and $\alpha_{min,1}(P_{\theta,l}) = 1.78 10^{-4}$, we can still strongly reduce $P_{p,comb}$ by reducing the covariance matrix $P_{\zeta_f}$ of $\hat{\zeta}_f$ via an appropriate design of the experimental conditions (i.e., $N_f$ and $\Phi_{uf}$).

**Design of the experimental conditions $N_f$ and $\Phi_{uf}$**. Let us first consider that it is acceptable to increase the duration of the experiment to $N_f = 10000$. Using $N_f = 10000$ and $\Phi_{uf}(\omega) = 0.1 \forall \omega$, we then obtain $P_{p,comb}(\Phi_{uf}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l}) = 3 10^{-4}$ which is much closer<fnref n="6" /> to the lower bound $\alpha_{min,1}(P_{\theta,l}) = 1.78 10^{-4}$. If it is not possible to increase the number of data $N_f$, another approach to efficiently tune $P_{\zeta_f}$ (while keeping $N_f = 100$) is to optimize the shape of the spectrum $\Phi_{uf}(\omega)$. For this purpose, we choose $L = 5$ in (42) and we solve the optimal experiment design problem (39)-(40) using that parametrization of the spectrum $\Phi_{uf}(\omega)$ and posing $P_{u,max} = 0.1$. The optimal spectrum in this case is represented in Figure 4 (blue dashdot curve) and this spectrum corresponds, such as it was the case for the white noise, to a signal $u_f$ of power $P_{u,max} = 0.1$. For this spectrum and $N_f = 100$, we obtain $P_{p,comb}(\Phi_{uf}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l}) = 0.005$. In other words, without changing the duration of the experiment and by just optimizing the frequency content of the signal $u_f$ of power $P_{u,max}$, we have been able to tune the covariance matrix $P_{\zeta_f}$ in such a way that the variance $P_{p,comb}$ of $\hat{\mathbf{p}}_{comb}$ is reduced by 57% ($P_{p,comb}$ was indeed equal to 0.0124 for the white noise $u_f$ of power $P_{u,max}$ and of duration $N_f = 100$).

**Dealing with the chicken-and-egg problem.** In the previous paragraph, we considered the idealized case where the unknown $\mathbf{p}^0$ and $\theta^0$ could be used to determine the optimal spectrum in Figure 4. We will now consider the modified optimal experiment design problem in Section 5.2 where $\theta^0$ is replaced by $\hat{\theta}_l$ and where $\mathbf{p}^0$ is replaced by a grid of the interval $\mathcal{P}$ containing $n_{grid} = 60$ points. Using here also $L = 5$ in (42), we solve the optimization problem (45)-(46). This leads to an optimal spectrum $\Phi_{uf}(\omega)$ represented in Figure 4 (red solid curve) and which can

<footnote n="6">⁶If we can both have $N_f = 10000$ and $\Phi_{uf}(\omega) = 10 \forall \omega$, we then further reduce $P_{\zeta_f}$ and obtain $P_{p,comb}(\Phi_{uf}, N_f, \mathbf{p}^0, \theta^0, P_{\theta,l}) = 1.79 10^{-4}$ which is then almost equal to $\alpha_{min,1}(P_{\theta,l}) = 1.78 10^{-4}$.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 4: Spectrum $\Phi_{u_f}(\omega)$ of a white noise signal of power $P_{u,max} = 0.1$ (black dashed) and optimal spectrum $\Phi_{u_f}(\omega)$ corresponding to the parametrization (42) with $L = 5$ and to the optimization problem (39)-(40) (blue dashdot). The optimal spectrum corresponding to the parametrization (42) with $L = 5$ and to the optimization problem (45)-(46) is given by the red solid curve.</figcaption></figure>
almost not be distinguished from the optimal spectrum in the idealized case (see the blue dashdot curve in Figure 4). If we compute the inverse of the optimal cost function $\beta^{opt}$, we obtain 0.008 which means that, independently of the value of the current operating point $\mathbf{p}^0 \in \mathcal{P}$, we will be able to achieve a variance $P_{p,comb}$ of 0.008. Moreover, not surprisingly since the two frequency-dependent spectra in Figure 4 can almost not be distinguished, the covariance matrix $P_{p,comb}$ at $\mathbf{p}^0 = 0.5$ is also equal to 0.005 (such in the idealized case). So, in this particular example, we arrive at the same result in the idealized case of (39)-(40) and in the realistic case of (45)-(46).

In order to verify the validity of the result, we have performed 1000 Monte-Carlo simulations where $\hat{\mathbf{p}}_{comb}$ was evaluated using criterion (24) with the initial estimate $\hat{\theta}_l$ and its covariance matrix $P_{\theta,l}$ given above and with $\hat{\zeta}_f$ and $P_{\zeta_f}$ determined using an experiment on the LPV system at operating point $\mathbf{p}^0 = 0.5$ with an excitation $u_f(t)$ of duration $N_f = 100$ and having the spectrum given by the red solid curve in Figure 4. The difference between the 1000 Monte-Carlo simulations is the realization of the white noise sequence $e$ during that experiment. This procedure has delivered 1000 estimates $\hat{\mathbf{p}}_{comb}^i$ of $\mathbf{p}^0$ ($i = 1, ..., 1000$) and we have observed that the sampled variance $\frac{1}{1000} \sum_{i=1}^{1000} (\hat{\mathbf{p}}_{comb}^i - \mathbf{p}^0)^2$ of these 1000 estimates is equal to 0.0047 and is thus extremely close to the variance of 0.005 predicted by our experimental design procedure.

# 7. Concluding remarks

In this paper, we consider the case of frozen LPV systems for which, in the operation phase, the constant operating point cannot be measured and has therefore to be identified using data collected on the system. We develop a maximum likelihood approach for this purpose and we also give optimal experiment design tools to optimize the identification of the estimate of the unknown operating point. In the future, we will aim at extending the results of this paper to the case where the scheduling variable is time-varying.

# Disclosure statement

There is no conflict of interest related to this paper.

# References

Apkarian, P., & Gahinet, P. (1995). A convex characterization of gain-scheduled $H_\infty$ controllers. *IEEE Transactions on Automatic Control* **40**(5), pp. 853-864.

Apkarian, P., Gahinet, P., & Biannic, J. (1994). Gain-scheduled $H_\infty$ control of missile via linear matrix inequalities. *Journal of Guidance, Control and Dynamics*, vol 28 no 3.

Bombois, X., Ghosh, D., Scorletti, G., & Huillery, J. (2021). LPV system identification for control using the local approach. *International Journal of Control*, **94**(2), 390-410.

Bombois, X., Scorletti, G., Gevers, M., Van den Hof, P., & Hildebrand, R. (2006). Least costly identification experiment for control. *Automatica*, **42**(10), 1651-1662.

Boyd, S., El Ghaoui, L., Feron, E., & Balakrishnan, V. (1994). *Linear matrix inequalities in systems and control theory* (Vol. 15). Philadelphia: SIAM.
Dankers, A., Toth, R., Heuberger, P., Bombois, X., & den Hof, P. V. (2011). Informative data and identifiability for LPV-ARX prediction error identification. In *Proc. iee conference on decision and control* (p. 799-804). Orlando.

Dinh, M., Scorletti, G., Fromion, V., & Magarotto, E. (2005). Parameter dependent $H_\infty$ control by finite dimensional lmi optimization. *Int. J. Robust Nonlinear Control*, **15**, 383-406.

Ghosh, D., Bombois, X., Huillery, J., Scorletti, G., & Mercere, G. (2018). Optimal identification experiment design for LPV systems using the local approach. *Automatica*, **87**, 258-266.

Hamdan, T., Coirault, P., Mercère, G., & Dairay, T. (2024). Data enabled predictive control of LPV systems. *Control Engineering Practice*, **149**, 105969.

Jansson, H., & Hjalmarsson, H. (2005). Input design via LMIs admitting frequency-wise model specifications in confidence regions. *IEEE Transactions on Automatic Control*, **50**(10), 1534-1549.

Liu, Q., Gross, J., Pfeiffer, S., & Werner, H. (2014). A local approach for the LPV identification of an actuated beam using piezoelectric actuators and sensors. *Mechatronics*, *Vol. 24*, pp. 289-297.

Ljung, L. (1999). *System identification: Theory for the user, 2nd edition*. Englewood Cliffs, NJ: Prentice-Hall.

Piga, D., Cox, P., Toth, R., & Laurain, V. (2015). LPV system identification under noise corrupted scheduling and output signal observations. *Automatica*, **53**, 329-338.

Rehmer, A., & Kroll, A. (2022). On affine quasi-LPV system identification with unknown state-scheduling using (deep) recurrent neural networks. In *26th international conference on system theory, control and computing* (pp. 446-451).

Rojas, C., Welsh, J., Goodwin, G., & Feuer, A. (2007). Robust optimal experiment design for system identification. *Automatica*, **43**(6), 993-1008.

Ruhe, A., & Wedin, P. (1980). Algorithms for separable nonlinear least squares problems. *SIAM review*, **22**(3), 318-337.

Scorletti, G., & Ghaoui, L. E. (1998). Improved LMI conditions for gain scheduling and related problems. *International Journal of Robust and Nonlinear Control*, **8**(10), 845-877.

Steinbuch, M., van de Molengraft, M., Heuberger, P., & der Voort, A. V. (2003). Experimental modelling and LPV control of a motion system. In *Proceedings of the american control conference*, pp. 1374-1379.

Toth, R. (2010). *Modeling and identification of linear parameter varying systems*. Berlin: Springer.

Verhoek, C., Beintema, G., Haesaert, S., Schoukens, M., & Tóth, R. (2022). Deep-learning-based identification of LPV models for nonlinear systems. In *61st IEEE conference on decision and control* (pp. 3274-3280).

Wassink, M. G., van de Wal, M., Scherer, C., & Bosgra, O. (2005). LPV control for a wafer stage: beyond the theoretical solution. *Control Engineering Practice*, *Vol 13*, pp. 231-245.

Yang, X., Liu, X., & Han, B. (2018). LPV model identification with an unknown scheduling variable in the presence of missing observations—a robust global approach. *IET Control Theory & Applications*, **12**(10), 1465-1473.

Zhou, K., & Doyle, J. (1998). *Essentials of robust control*. Prentice Hall, Upper Saddle River, New Jersey.

## Appendix A. Optimization problem (24)

### A.1. General considerations

The optimization problem (24) is a nonlinear optimization problem due to the product $T(\mathbf{p})\theta$ and the fact that both $\mathbf{p}$ and $\theta$ are decision variables. However, solving this nonlinear optimization problem using a gradient-based technique is here facilitated by
the fact that the gradient and the Hessian of $V(\eta)$ are computable (see Appendix A.2) and by the fact that the gradient-based algorithm can easily be initialized at a reasonable value $\eta_{init} = (\mathbf{p}_{init}, \theta_{init}^T)^T$ where $\theta_{init}$ is chosen as the available initial estimate $\hat{\theta}_l$ of $\theta^0$ and where

$$ \mathbf{p}_{init} = \arg\min_{\mathbf{p}} \mathcal{V}(\mathbf{p}) \qquad (A1) $$

$$ \text{with } \mathcal{V}(\mathbf{p}) = V \left( \eta = (\mathbf{p}, \hat{\theta}_l^T)^T \right) = \frac{1}{2} \left( \left( T(\mathbf{p})\hat{\theta}_l - \hat{\zeta}_f \right)^T P_{\zeta_f}^{-1} \left( T(\mathbf{p})\hat{\theta}_l - \hat{\zeta}_f \right) \right). \quad (A2) $$

The nonlinear optimization problem (A1)-(A2) has one scalar decision variable $\mathbf{p}$ and can thus, e.g., be addressed using line search in the scheduling interval $\mathcal{P}$. If the mapping (3) is a polynomial mapping (i.e., $f_j(p(t)) = p^j(t)$ in (2)), the optimization problem (A1)-(A2) can also be solved as the solution of a polynomial equation (see Appendix A.3).

**Remark.** An alternative to the approach presented in the previous paragraph follows from the fact that (24) is a so-called separable least-squares problem (Ruhe & Wedin, 1980). Indeed, if we fix the value of $\mathbf{p}$ to $\mathbf{p}^*$, we observe that $\theta(\mathbf{p}^*) = \arg\min_{\theta} V(\eta = (\mathbf{p}^*, \theta^T)^T)$ has a closed-form solution:

$$ \theta(\mathbf{p}^*) = \left( T^T(\mathbf{p}^*) P_{\zeta_f}^{-1} T(\mathbf{p}^*) + P_{\theta,l}^{-1} \right)^{-1} \left( T^T(\mathbf{p}^*) P_{\zeta_f}^{-1} \hat{\zeta}_f + P_{\theta,l}^{-1} \hat{\theta}_l \right) $$

We then have that $\hat{\eta}_{comb}$ (see (24)) is given by $\hat{\eta}_{comb} = \begin{pmatrix} \hat{\mathbf{p}}_{comb} \\ \theta(\hat{\mathbf{p}}_{comb}) \end{pmatrix}$ with $\hat{\mathbf{p}}_{comb}$ the value of $\mathbf{p}$ minimizing $V(\eta = \begin{pmatrix} \mathbf{p} \\ \theta(\mathbf{p}) \end{pmatrix})$ i.e., the cost function (23) obtained by replacing $\theta$ by $\theta(\mathbf{p})$. The latter optimization problem can, e.g., be addressed using a line search in $\mathcal{P}$.

## A.2. Gradient and Hessian of $V(\eta)$

The gradient $\Gamma(\eta)$ of the cost function $V(\eta)$ where $\eta = (\mathbf{p}, \theta^T)^T \in \mathbb{R}^{k+1}$ is given by:

$$ \Gamma(\eta) \triangleq \frac{dV(\eta)}{d\eta} = \begin{pmatrix} \frac{dV(\eta)}{d\mathbf{p}} \\ \frac{dV(\eta)}{d\theta} \end{pmatrix} = \begin{pmatrix} \theta^T D^T(\mathbf{p}) P_{\zeta_f}^{-1} (T(\mathbf{p})\theta - \hat{\zeta}_f) \\ T^T(\mathbf{p}) P_{\zeta_f}^{-1} (T(\mathbf{p})\theta - \hat{\zeta}_f) + P_{\theta,l}^{-1} (\theta - \hat{\theta}_l) \end{pmatrix}. \quad (A3) $$

with $D(\mathbf{p}) = \frac{dT(\mathbf{p})}{d\mathbf{p}}$. Using (26), $\Gamma(\eta)$ can be rewritten in the following compact form:

$$ \Gamma(\eta) = W^T(\eta) P_{\zeta_f}^{-1} (T(\mathbf{p})\theta - \hat{\zeta}_f) + \begin{pmatrix} 0 \\ P_{\theta,l}^{-1} (\theta - \hat{\theta}_l) \end{pmatrix} \quad (A4) $$
The Hessian $\mathcal{H}(\eta)$ of $V(\eta)$ is given by:

$$ \mathcal{H}(\eta) \triangleq \frac{d^2 V(\eta)}{d\eta^2} = \begin{pmatrix} \mathcal{H}_{pp}(\eta) & \mathcal{H}_{p\theta}(\eta) \\ \mathcal{H}_{p\theta}^T(\eta) & \mathcal{H}_{\theta\theta}(\eta) \end{pmatrix} \qquad (A5) $$

$$ \mathcal{H}_{pp}(\eta) = \theta^T \frac{d D^T(\mathbf{p})}{d\mathbf{p}} P_{\zeta_f}^{-1} (T(\mathbf{p})\theta - \hat{\zeta}_f) + \theta^T D^T(\mathbf{p}) P_{\zeta_f}^{-1} D(\mathbf{p})\theta $$

$$ \mathcal{H}_{p\theta}(\eta) = \theta^T D^T(\mathbf{p}) P_{\zeta_f}^{-1} T(\mathbf{p}) $$

$$ \mathcal{H}_{\theta\theta}(\eta) = T^T(\mathbf{p}) P_{\zeta_f}^{-1} T(\mathbf{p}) + P_{\theta,l}^{-1} $$

### A.3. Solution of the optimization problem (A1) when the mapping (3) is a polynomial mapping

The gradient of $\mathcal{V}(\mathbf{p})$ (see (A2)) is $\hat{\theta}_l^T D^T(\mathbf{p}) P_{\zeta_f}^{-1} (T(\mathbf{p})\hat{\theta}_l - \hat{\zeta}_f)$. The solution $\mathbf{p}_{init}$ of (A1) will therefore be a value $\mathbf{p}$ such that

$$ \hat{\theta}_l^T D^T(\mathbf{p}) P_{\zeta_f}^{-1} (T(\mathbf{p})\hat{\theta}_l - \hat{\zeta}_f) = 0 \qquad (A6) $$

Using now the fact that $f_j(p(t)) = p^j(t)$ in (2), we have that

$$ T(\mathbf{p})\hat{\theta}_l = \sum_{j=0}^{n_p} \alpha_j p^j $$

$$ D(\mathbf{p})\hat{\theta}_l = \sum_{j=0}^{n_p-1} \beta_j p^j $$

for known vectors $\alpha_j \in \mathbb{R}^n$ ($j = 0, ..., n_p$) and for $\beta_j = (j+1) \alpha_{j+1}$ ($j = 0, ..., n_p - 1$). Using these expressions, we see that (A6) boils down to the following polynomial equation:

$$ \sum_{j=0}^{n_p} \sum_{i=0}^{n_p-1} (\beta_i^T P_{\zeta_f}^{-1} \alpha_j) p^{i+j} = \sum_{j=0}^{n_p-1} (\beta_j^T P_{\zeta_f}^{-1} \hat{\zeta}_f) p^j $$

If this polynomial equation of degree $2n_p - 1$ has one single real root, this real root is then the solution $\mathbf{p}_{init}$ of the optimization problem (A1). If it has more than one
real roots, $\mathbf{p}_{init}$ is the real root corresponding to the smallest value of the cost function $V(\mathbf{p})$.

## Appendix B. Proof of Proposition 3.1

For the simplicity of the presentation, let us first normalize the cost function $V(\eta)$ by dividing it by $N_f$. Using (23), (12) and (20), we obtain:

$$ 
\begin{aligned}
V_n(\eta) &= \frac{1}{N_f} V(\eta) \\
&= \frac{1}{2} \left( \left( (T(\mathbf{p})\theta - \hat{\zeta}_f)^T \mathcal{I}_{\zeta_f} (T(\mathbf{p})\theta - \hat{\zeta}_f) + \rho (\theta - \hat{\theta}_l)^T \mathcal{I}_{\theta,l} (\theta - \hat{\theta}_l) \right) \right) \quad (B1)
\end{aligned}
$$

with $\mathcal{I}_{\zeta_f}$ and $\mathcal{I}_{\theta,l}$ defined in (13) and (20) and with $\rho \stackrel{\Delta}{=} \frac{N_l}{N_f}$. Note that, when $N_l, N_f \to \infty$, the scalar quantity $\rho = \frac{N_l}{N_f}$ remains strictly positive ($\rho > 0$) and finite ($\rho < \infty$). Moreover, since the matrices $\mathcal{I}_{\zeta_f}$ (resp. $\mathcal{I}_{\theta,l}$) are independent of the value of $N_f$ (resp. $N_l$) and since $\lambda_{max}(\mathcal{I}_{\zeta_f}) < \infty$ and $\lambda_{max}(\mathcal{I}_{\theta,l}) < \infty$, the normalized criterion $V_n(\theta)$ will therefore remain finite even when $N_l, N_f \to \infty$ (which is the situation considered in Proposition 3.1).

It is clear that the criterion (24) can be equivalently rewritten as:

$$ \hat{\eta}_{comb} = \arg \min_{\eta} V_n(\eta) \quad (B2) $$

To prove Proposition 3.1, we will use similar arguments as in Sections 9.2 and 9.3 of (Ljung, 1999), but adapt them to the situation of this proposition.

Due to (B2), we have that $\Gamma_n(\hat{\eta}_{comb}) = 0$ where $\Gamma_n(\eta)$ is the gradient of $V_n(\eta)$. This expression can be expanded around $\eta^0$ using Taylor series:

$$ \Gamma_n(\eta^0) + \mathcal{H}_n(\tilde{\eta}) (\hat{\eta}_{comb} - \eta^0) = 0 \quad (B3) $$

with $\tilde{\eta}$ in the segment between $\hat{\eta}_{comb}$ and $\eta^0$ and with $\mathcal{H}_n(\tilde{\eta})$ the Hessian of $V_n(\eta)$. Under the asymptotic assumption made in Proposition 3.1, the expression (B3) will allow us to derive the statistical properties of $\hat{\eta}_{comb}$. For this purpose, let us first specify the consequences of this asymptotic assumption:

(1) in the asymptotic case, we have that $\hat{\theta}_l$ and $\hat{\zeta}_f$ present the normal distribution given respectively in (18) and (11),

(2) as shown in Appendix C, the matrix $\mathcal{H}_n(\tilde{\eta})$ (which is a random variable for finite $N_l$ and $N_f$) asymptotically converges to the deterministic and strictly positive definite matrix $X$:
$$X = W^T(\eta^0) \mathcal{I}_{\zeta_f} W(\eta^0) + \rho b \mathrm{diag}(0, \mathcal{I}_{\theta_l}) \quad (B4)$$

Consequently, in the asymptotic situation, (B3) can be rewritten as<fnref n="7" />:

$$\hat{\eta}_{comb} - \eta^0 = -X^{-1} \Gamma_n(\eta^0) \quad (B5)$$

and we can use (A4) to find the following expression for the term $\Gamma_n(\eta^0)$ in (B5):

$$\begin{align}
\Gamma_n(\eta^0) &= \frac{1}{N_f} \Gamma(\eta^0) = \frac{1}{N_f} \left( W^T(\eta^0) P_{\zeta_f}^{-1} (T(\mathbf{p}^0)\theta^0 - \hat{\zeta}_f) + \begin{pmatrix} 0 \\ P_{\theta_l}^{-1}(\theta^0 - \hat{\theta}_l) \end{pmatrix} \right) \tag{B6} \\
&= W^T(\eta^0) \mathcal{I}_{\zeta_f} (T(\mathbf{p}^0)\theta^0 - \hat{\zeta}_f) + \rho \begin{pmatrix} 0 \\ \mathcal{I}_{\theta_l}(\theta^0 - \hat{\theta}_l) \end{pmatrix} \nonumber
\end{align}$$

with $W(\eta^0)$ defined in (26).

Let us observe that, in the asymptotic case, $\Gamma_n(\eta^0)$ (see (B6)) is normally distributed since (B6) is a linear combination of the normally distributed variables $\hat{\zeta}_f$ and $\hat{\theta}_l$ and that, by virtue of (18) and (11), the expected value of $\Gamma_n(\eta^0)$ is equal to zero i.e., $E \Gamma_n(\eta^0) = 0$. Combining these facts and the fact that $X$ is deterministic, (B5) yields:

$$\hat{\eta}_{comb} \sim \mathcal{N}(\eta^0, P)$$

with $P = X^{-1}QX^{-1}$ where $Q = E \Gamma_n(\eta^0)\Gamma_n^T(\eta^0)$. We have now to prove that $P = X^{-1}QX^{-1}$ is equal to $P_{\eta,comb}$ (see (25)).

From (B4), it is clear that $X = \frac{1}{N_f}P_{\eta,comb}^{-1}$. Using the fact that $\hat{\zeta}_f$ and $\hat{\theta}_l$ are independent and using also the distributions (18) and (11), we also easily derive that $Q = E \Gamma_n(\eta^0)\Gamma_n^T(\eta^0) = \frac{1}{N_f^2}P_{\eta,comb}^{-1}$. We have thus indeed that $P = X^{-1}QX^{-1} = P_{\eta,comb}$.

## Appendix C. Asymptotic expression for $H_n(\tilde{\eta})$ in (B3)

The matrix $H_n(\tilde{\eta})$ is the Hessian $H_n(\eta)$ of $V_n(\eta)$ evaluated at a value $\tilde{\eta}$ in the segment between $\hat{\eta}_{comb}$ and $\eta^0$. Consequently, in order to determine its asymptotic expression, we need to first determine the asymptotic expression of the cost function $V_n(\eta)$ (see (B1)). For infinitely large value of $N_f$ and $N_f$, the estimates $\hat{\theta}_l$ and $\hat{\zeta}_f$ will converge respectively to $\theta^0$ and $\zeta^0(\mathbf{p}^0) = T(\mathbf{p}^0)\theta^0$. This therefore yields the following asymptotic expression for $V_n(\eta)$:

<footnote n="7">⁷$X$ is invertible since it is strictly positive definite (see Appendix C).</footnote>
$$ \bar{V}_n(\eta) = \frac{1}{2} \left( (T(\mathbf{p})\theta - \zeta^0(\mathbf{p}^0))^T \mathcal{I}_{\zeta_f} (T(\mathbf{p})\theta - \zeta^0(\mathbf{p}^0)) + \rho (\theta - \theta^0)^T \mathcal{I}_{\theta,l} (\theta - \theta^0) \right) \quad (C1) $$

with $\mathcal{I}_{\zeta_f} > 0$ and $\mathcal{I}_{\theta,l} > 0$. Observe that, by virtue of Assumption 2.1, the unique minimizer of $\bar{V}_n(\eta)$ is $\eta^0$. Consequently, for infinitely large $N_l$ and $N_f$, $\hat{\eta}_{comb}$ converges to $\eta^0$.

Since $V_n(\eta)$ asymptotically converges to $\bar{V}_n(\eta)$, the Hessian $\mathcal{H}_n(\eta)$ of $V_n(\eta)$ asymptotically converges to the Hessian $\bar{\mathcal{H}}_n(\eta)$ of $\bar{V}_n(\eta)$. Moreover, since $\hat{\eta}_{comb}$ asymptotically converges to $\eta^0$, $\tilde{\eta}$ also converges to $\eta^0$. Combining these two facts, we see that the asymptotic expression of $\mathcal{H}_n(\tilde{\eta})$ is $\bar{\mathcal{H}}_n(\eta^0)$. The Hessian $\bar{\mathcal{H}}_n(\eta)$ of $\bar{V}_n(\eta)$ has the following expression:

$$ \bar{\mathcal{H}}_n(\eta) = \begin{pmatrix} \bar{\mathcal{H}}_{pp}(\eta) & \bar{\mathcal{H}}_{p\theta}(\eta) \\ \bar{\mathcal{H}}_{p\theta}^T(\eta) & \bar{\mathcal{H}}_{\theta\theta}(\eta) \end{pmatrix} \quad (C2) $$

$$ \text{with } \bar{\mathcal{H}}_{pp}(\eta) = \theta^T \frac{dD^T(\mathbf{p})}{d\mathbf{p}} \mathcal{I}_{\zeta_f} (T(\mathbf{p})\theta - \zeta^0(\mathbf{p}^0)) + \theta^T D^T(\mathbf{p}) \mathcal{I}_{\zeta_f} D(\mathbf{p})\theta $$

$$ \bar{\mathcal{H}}_{p\theta}(\eta) = \theta^T D^T(\mathbf{p}) \mathcal{I}_{\zeta_f} T(\mathbf{p}) $$

$$ \bar{\mathcal{H}}_{\theta\theta}(\eta) = T^T(\mathbf{p}) \mathcal{I}_{\zeta_f} T(\mathbf{p}) + \rho \mathcal{I}_{\theta,l} $$

and $\bar{\mathcal{H}}_n(\eta^0)$ is thus given by:

$$ \bar{\mathcal{H}}_n(\eta^0) = W^T(\eta^0) \mathcal{I}_{\zeta_f} W(\eta^0) + \rho b \mathrm{diag}(0, \mathcal{I}_{\theta_l}) \quad (C3) $$

with $W(\eta^0)$ as defined in (26). It is clear that the matrix $\bar{\mathcal{H}}_n(\eta^0)$ given in (C3) is equal to the matrix $X$ given in (B4). This matrix $\bar{\mathcal{H}}_n(\eta^0) = X$ is also clearly deterministic and it is furthermore strictly positive definite. To show the latter, let us consider an arbitrary vector $x \in \mathbb{R}^{k+1}$ and let us decompose it as $x = (x_1, x_2^T)^T$ with $x_1 \in \mathbb{R}$ and $x_2 \in \mathbb{R}^k$. Then, we have that:

$$ x^T X x = x^T (W^T(\eta^0) \mathcal{I}_{\zeta_f} W(\eta^0)) x + \rho x_2^T \mathcal{I}_{\theta_l} x_2 $$

Since $\mathcal{I}_{\zeta_f} > 0$, $\mathcal{I}_{\theta_l} > 0$ and $\rho > 0$, we have that $x^T(W^T(\eta^0) \mathcal{I}_{\zeta_f} W(\eta^0))x \ge 0$ for all vectors $x$ and $\rho x_2^T \mathcal{I}_{\theta_l} x_2 > 0$ for all vectors $x_2 \ne 0$. In other words, $x^T X x$ is strictly larger than zero for all vectors $x$ such that $x_2 \ne 0$. When $x_2 = 0$, $x^T X x$ reduces to $x_1^2(D(\mathbf{p}^0)\theta^0)^T \mathcal{I}_{\zeta_f} D(\mathbf{p}^0)\theta^0$ and this scalar is also strictly positive since $\mathcal{I}_{\zeta_f} > 0$ and since $D(\mathbf{p}^0)\theta^0 \ne 0$ (the latter being a consequence of Assumption 2.1). We have thus proven that $x^T X x > 0$ for all $x \in \mathbb{R}^{k+1}$ and thus that $X > 0$.
**Remark.** Recall that $\lambda_{max}(\mathcal{I}_{\zeta_f}) < \infty$ and $\lambda_{max}(\mathcal{I}_{\theta_l}) < \infty$ and also that $\rho$ remains finite when $N_f, N_l \to \infty$. Consequently, we have also that $\lambda_{max}(X) < \infty$ when $N_f, N_l \to \infty$.

## Appendix D. Analysis of the expression (25) for $P_{\eta,comb}$

Using the expressions (20) and (12), the expression (25) for $P_{\eta,comb} \in \mathbb{R}^{(k+1) \times (k+1)}$ can be rewritten as:

$$P_{\eta,comb}^{-1} = N_f \underbrace{W^T(\eta^0) \mathcal{I}_{\zeta_f} W(\eta^0)}_{=A_f} + N_l \underbrace{bdiag(0, \mathcal{I}_{\theta_l})}_{=A_l} \quad (D1)$$

where $\mathcal{I}_{\zeta_f} \in \mathbb{R}^{n \times n}$ and $\mathcal{I}_{\theta_l} \in \mathbb{R}^{k \times k}$ are both strictly positive-definite. Let us also recall that $\mathcal{I}_{\zeta_f}$ and $\mathcal{I}_{\theta_l}$ are not function of $N_f$ and $N_l$ and that $\lambda_{max}(\mathcal{I}_{\zeta_f}) < \infty$ and $\lambda_{max}(\mathcal{I}_{\theta_l}) < \infty$.

Following the same reasoning as in Appendix C, it can be proven that (D1) is always strictly positive-definite. Moreover, when both $N_l$ and $N_f$ are finite, we have that $\lambda_{max}(P_{\eta,comb}^{-1}) < \infty$. Consequently, $P_{\eta,comb} > 0$ when both $N_l$ and $N_f$ are finite.

Let us now observe that, in (D1), the matrix $A_l \ge 0$ is of rank $k$ and that $W(\eta) \in \mathbb{R}^{n \times (k+1)}$ (see (26)) is full row rank since $T(p)$ is full row rank. Consequently, using the Sylvester rank formula, $A_f \ge 0$ is of rank $n$.

Let us now analyze the case where $N_f \to \infty$ and $N_l$ is finite. In this case, since $A_f$ is of rank $n$, the matrix $P_{\eta,comb}^{-1}$ will have $n$ infinite eigenvalues and $k+1-n$ finite eigenvalues and thus $P_{\eta,comb}$ will be positive semi-definite with $n$ zero eigenvalues.

Let us now analyze the case where $N_l \to \infty$ and $N_f$ is finite. In this case, since $A_l$ is of rank $k$, the matrix $P_{\eta,comb}^{-1}$ will have $k$ infinite eigenvalues and one finite eigenvalue and thus $P_{\eta,comb}$ will be positive semi-definite with $k$ zero eigenvalues.

Finally, let us analyze the case where both $N_l, N_f \to \infty$. For this purpose, let us further develop (D1) as follows:

$$P_{\eta,comb} = \frac{1}{N_f} \left( \underbrace{W^T(\eta^0) \mathcal{I}_{\zeta_f} W(\eta^0) + \rho bdiag(0, \mathcal{I}_{\theta_l})}_{=X} \right)^{-1}$$

where the matrix $X$ has the properties that $X > 0$ and $\lambda_{max}(X) < \infty$ (see Appendix C). Since $\rho = \frac{N_l}{N_f}$ remains finite when both $N_l, N_f \to \infty$, it is clear that $P_{\eta,comb}$ converges to the zero matrix when both $N_l, N_f \to \infty$.
Appendix E. Proof of the relation (33)

Let us first consider the case where the initial estimate $\hat{\theta}_l$ has a nonzero covariance matrix $P_{\theta,l}$. In this case, $P_{\theta,l}$ is thus a strictly positive definite matrix (see (20)).

Let us introduce the following shorthand notations $v = D(\mathbf{p}^0)\theta^0$ and $T = T(\mathbf{p}^0)$ and let us recall for the sequel that $T = T(\mathbf{p}^0)$ is full row rank and that, by virtue of Assumption 2.1, $v = D(\mathbf{p}^0)\theta^0 \neq 0$. Using these shorthand notations, the expression (32) of $P_{\theta,comb}$ yields:

$$ P_{\theta,comb}^{-1} = P_{\theta,l}^{-1} + T^T P_{\zeta_f}^{-1} T - T^T P_{\zeta_f}^{-1} v \frac{1}{v^T P_{\zeta_f}^{-1} v} v^T P_{\zeta_f}^{-1} T \quad (E1) $$

Let us also consider the expression $P_{\zeta_f} = \frac{1}{N_f} \mathcal{I}_{\zeta_f}^{-1}$ with $\mathcal{I}_{\zeta_f} > 0$ (see (12)) and let us introduce the notations: $\mathcal{I}_{\zeta_f} = \tilde{B}^T \tilde{B}$ with $\tilde{B} > 0$, $\tilde{X} = \tilde{B}T$, $B = \sqrt{N_f} \tilde{B} > 0$, $X = BT$ and finally $\phi = Bv \neq 0$. Then (E1) can be equivalently rewritten as:

$$ P_{\theta,comb}^{-1} = P_{\theta,l}^{-1} + X^T \left( I_n - \phi (\phi^T \phi)^{-1} \phi^T \right) X \quad (E2) $$

Since $\phi$ is a column vector of dimension $n$, we have that $I_n - \phi (\phi^T \phi)^{-1} \phi^T \geq 0$, which also implies that $A \geq 0$ and thus that $P_{\theta,comb}^{-1} \geq P_{\theta,l}^{-1}$. Since $P_{\theta,l}^{-1} > 0$, the latter is equivalent to $P_{\theta,comb} \leq P_{\theta,l}$ (see Appendix G).

We have thus proven that the relation (33) holds when $P_{\theta,l} \neq 0$. Note also that the two matrices $P_{\theta,comb}$ and $P_{\theta,l}$ are different since $A \neq 0$. Combining this fact with $P_{\theta,comb} \leq P_{\theta,l}$, we see that $\hat{\theta}_{comb}$ is indeed a more accurate estimate than $\hat{\theta}_l$.

Let us now consider the case where $P_{\theta,l}$ is equal to the zero matrix i.e., the case where $N_l \to \infty$ in the expression (20). We will show below that $P_{\theta,comb}$ is also equal to the zero matrix when $N_l \to \infty$ and thus that the relation (33) also holds when $P_{\theta,l} = 0$.

For this purpose, let us use (20) to rewrite (E2) as follows:

$$ \begin{aligned} P_{\theta,comb} &= (P_{\theta,l}^{-1} + A)^{-1} \\ &= \frac{1}{N_l} \left( \mathcal{I}_{\theta,l} + \frac{1}{N_l} A \right)^{-1} \end{aligned} \quad (E3) $$

where $\mathcal{I}_{\theta,l} > 0$ is independent of $N_l$. Let us now use the notations $\tilde{X} = \tilde{B}T$ and $\tilde{\phi} = \tilde{B}v$ to rewrite $A$ as follows
$$A = N_f \underbrace{\tilde{X}^T (I_n - \tilde{\phi} (\tilde{\phi}^T \tilde{\phi})^{-1} \tilde{\phi}^T) \tilde{X}}_{=\tilde{A}} \quad (E4)$$

where $\tilde{A} \ge 0$ is independent of $N_f$. We have thus:

$$P_{\theta, comb} = \frac{1}{N_l} \left( \mathcal{I}_{\theta_l} + \frac{N_f}{N_l} \tilde{A} \right)^{-1}$$

where $\frac{N_f}{N_l}$ is a finite positive scalar when both $N_l, N_f \to \infty$, while $\frac{N_f}{N_l} = 0$ when $N_l \to \infty$ and $N_f$ is finite. Consequently, since $\mathcal{I}_{\theta_l} > 0$ and $\tilde{A} \ge 0$, we have that, independently of the value of $N_f$, $P_{\theta, comb}$ is equal to the zero matrix when $N_l \to \infty$, which was the desired result.

## Appendix F. Proof of Proposition 4.1

Let us here also introduce the following shorthand notations $v = D(\mathbf{p}^0)\theta^0$ and $T = T(\mathbf{p}^0)$ and let us recall for the sequel that $T = T(\mathbf{p}^0)$ is full row rank and that, by virtue of Assumption 2.1, $v = D(\mathbf{p}^0)\theta^0 \ne 0$. Using these shorthand notations, the expression (31) of $P_{p, comb}$ becomes:

$$P_{p, comb}^{-1} = v^T P_{\zeta_f}^{-1} v - v^T P_{\zeta_f}^{-1} T (T^T P_{\zeta_f}^{-1} T + P_{\theta, l}^{-1})^{-1} T^T P_{\zeta_f}^{-1} v \quad (F1)$$

Let us now consider the expression $P_{\zeta_f} = \frac{1}{N_f} \mathcal{I}_{\zeta_f}^{-1}$ with $\mathcal{I}_{\zeta_f} > 0$ (see (12)) and let us use the following notations: $\mathcal{I}_{\zeta_f} = \tilde{B}^T \tilde{B}$ with $\tilde{B} > 0$, $\tilde{X} = \tilde{B}T$, $B = \sqrt{N_f} \tilde{B}$, $X = BT = \sqrt{N_f} \tilde{B}T$. Then (F1) can be equivalently rewritten as:

$$P_{p, comb}^{-1} = v^T B^T \left( I_n - X \left( P_{\theta, l}^{-1} + X^T X \right)^{-1} X^T \right) B v \quad (F2)$$

Let us first consider item (i) of Proposition 4.1 (where we assume $P_{\theta, l} \ne 0$ and thus $P_{\theta, l} > 0$) and let us observe that, by virtue of the Woodbury matrix identity, the matrix $I_n - X (P_{\theta, l}^{-1} + X^T X)^{-1} X^T$ in (F2) is equal to $(I_n + X P_{\theta, l} X^T)^{-1}$. Consequently, (F2) is equivalent to:

$$P_{p, comb}^{-1} = v^T B^T (I_n + X P_{\theta, l} X^T)^{-1} B v \quad (F3)$$

Using the expressions $B = \sqrt{N_f} \tilde{B}$ and $X = BT = \sqrt{N_f} \tilde{B}T$, (F3) becomes:
$$P_{p,comb}^{-1} = v^T \tilde{B}^T \left( \frac{I_n}{N_f} + \tilde{B} T P_{\theta,l} T^T \tilde{B}^T \right)^{-1} \tilde{B} v \quad (F4)$$

where $\tilde{B} T P_{\theta,l} T^T \tilde{B}^T > 0$ since $\tilde{B} > 0$, $P_{\theta,l} > 0$ and $T$ is full row-rank. Since $\tilde{B} > 0$, let us also observe the following equalities:

$$v^T \tilde{B}^T \left( \tilde{B} T P_{\theta,l} T^T \tilde{B}^T \right)^{-1} \tilde{B} v = v^T (T P_{\theta,l} T^T)^{-1} v = \frac{1}{\alpha_{min,1}(P_{\theta,l})} \quad (F5)$$

Let us now analyze two cases. The first case is the case where $N_f$ is finite (i.e, $N_f < \infty$ and thus $P_{\zeta_f} \neq 0$ and $P_{\zeta_f} > 0$) and the second case is the case $N_f \to \infty$ (i.e, $P_{\zeta_f} = 0$).

In the first case (finite $N_f$ and $P_{\zeta_f} \neq 0$), we have that $\frac{I_n}{N_f} + \tilde{B} T P_{\theta,l} T^T \tilde{B}^T > \tilde{B} T P_{\theta,l} T^T \tilde{B}^T$ and thus that $(\frac{I_n}{N_f} + \tilde{B} T P_{\theta,l} T^T \tilde{B}^T)^{-1} < (\tilde{B} T P_{\theta,l} T^T \tilde{B}^T)^{-1}$ (see Appendix G). Consequently, using (F4) and (F5), we obtain that:

$$P_{p,comb} > \alpha_{min,1}(P_{\theta,l}) \text{ when } P_{\zeta_f} \neq 0. \quad (F6)$$

In the second case ($N_f \to \infty$ and $P_{\zeta_f} = 0$), we have that $\frac{I_n}{N_f} + \tilde{B} T P_{\theta,l} T^T \tilde{B}^T = \tilde{B} T P_{\theta,l} T^T \tilde{B}^T$ and thus that $P_{p,comb} = \alpha_{min,1}(P_{\theta,l})$. This concludes the proof of item (i) of the proposition.

Let us now consider item (ii) of Proposition 4.1 (where we assume $P_{\zeta_f} \neq 0$ and thus $P_{\zeta_f} > 0$). Using the fact that $P_{\theta,l}^{-1} = \frac{1}{\pi_l} \mathcal{I}_{\theta,l}$ with $\pi_l \stackrel{\Delta}{=} \frac{1}{N_l}$ (see (20)), we can equivalently rewrite (F2) as:

$$P_{p,comb}^{-1} = v^T B^T \left( I_n - \pi_l \underbrace{X (\mathcal{I}_{\theta_l} + \pi_l X^T X)^{-1} X^T}_{=L} \right) B v \quad (F7)$$

where $L > 0$ since $\pi_l \ge 0$, $\mathcal{I}_{\theta_l} > 0$ and $X = BT$ is full row rank. Let us also observe that $\alpha_{min,2}(P_{\zeta_f})$ can be equivalently rewritten as:

$$\alpha_{min,2}(P_{\zeta_f}) = \frac{1}{v^T B^T B v} \quad (F8)$$

Like for item (i), let us here also analyze two cases. The first case is the case where $N_l$ is finite (i.e, $\pi_l > 0$ and also $P_{\theta,l} \neq 0$ and $P_{\theta,l} > 0$) and the second case is the case $N_l \to \infty$ (i.e, $\pi_l = 0$ and also $P_{\theta,l} = 0$).

In the first case ($P_{\theta,l} \neq 0$ and $\pi_l > 0$), we have that $I_n - \pi_l L < I_n$. Consequently, using (F7), we have that $P_{p,comb}^{-1} < v^T B^T B v$ and thus that (see (F8)):
$$P_{p,comb} > \alpha_{min,2}(P_{\zeta_f}) \text{ when } P_{\theta,l} \neq 0. \qquad (F9)$$

In the second case ($N_l \to \infty$, $\pi_l = 0$ and $P_{\theta,l} = 0$), we have that $I_n - \pi_l L = I_n$ and thus that $P_{p,comb} = \alpha_{min,2}(P_{\zeta_f})$. This concludes the proof of item (ii) of the proposition.

We have also proven item (iii) since this item is equivalent to the combination of (F6) and of (F9). Finally, item (iv) is a consequence of the fact that $P_{\eta,comb} \to 0$ (and thus $P_{p,comb} \to 0$) when both $P_{\theta,l} \to 0$ and $P_{\zeta_f} \to 0$ (see Appendix D). Item (iv) can also be deduced from items (i) and (ii). These two items have indeed as consequence that, when both $P_{\theta,l} \to 0$ and $P_{\zeta_f} \to 0$, $P_{p,comb} = \alpha_{min,1}(P_{\theta,l} = 0) = \alpha_{min,2}(P_{\zeta_f} = 0) = 0$. That $\alpha_{min,1}(P_{\theta,l} = 0) = \alpha_{min,2}(P_{\zeta_f} = 0) = 0$ can be proven as follows. Using the expressions (12) for $P_{\zeta_f}$ and the expression (20) for $P_{\theta,l}$, we have that:

$$\alpha_{min,1}(P_{\theta,l} = \frac{1}{N_l} \mathcal{I}_{\theta_l}^{-1}) = \frac{1}{N_l} \frac{1}{(D(\mathbf{p}^0)\theta^0)^T (T(\mathbf{p}^0)\mathcal{I}_{\theta_l}^{-1} T^T(\mathbf{p}^0))^{-1} D(\mathbf{p}^0)\theta^0}$$

$$\alpha_{min,2}(P_{\zeta_f} = \frac{1}{N_f} \mathcal{I}_{\zeta_f}^{-1}) = \frac{1}{N_f} \frac{1}{(D(\mathbf{p}^0)\theta^0)^T \mathcal{I}_{\zeta_f} D(\mathbf{p}^0)\theta^0}$$

from which we observe that $\alpha_{min,1}(P_{\theta,l})$ converges to zero when $N_l \to \infty$ (i.e., when $P_{\theta,l} = 0$) and that $\alpha_{min,2}(P_{\zeta_f})$ converges to zero when $N_f \to \infty$ (i.e., when $P_{\zeta_f} = 0$).

## Appendix G. Lemma based on the Schur complements

Let us develop the following result which is a consequence of the Schur complement (Boyd et al., 1994).

**Lemma G.1.** Consider two strictly positive-definite matrices $A \in \mathbb{R}^{k \times k}$ and $B \in \mathbb{R}^{n \times n}$. Consider also an arbitrary matrix $C \in \mathbb{R}^{n \times k}$. Then, the three matrix inequalities (G1), (G2) and (G3) are equivalent:

$$B - CAC^T \geq 0 \qquad (G1)$$

$$\begin{pmatrix} B & C \\ C^T & A^{-1} \end{pmatrix} \geq 0 \qquad (G2)$$

$$A^{-1} - C^T B^{-1} C \geq 0. \qquad (G3)$$

*Proof.* Since $A > 0$, the fact that (G1) is equivalent to (G2) is the classical Schur complement result (Boyd et al., 1994). Let us also note that (G2) is equivalent to:
$$ \begin{pmatrix} A^{-1} & C^T \\ C & B \end{pmatrix} \geq 0 \qquad (G4) $$

Since $B > 0$, another application of the Schur complement result on (G4) shows that (G3) is indeed equivalent to (G4) and thus also to (G1) and (G2). □

**Remark.** Lemma G.1 also holds when the sign $\ge$ is replaced by $>$ in (G1), (G2) and (G3).