# Flow recovery from distal pressure in linearized hemodynamics: an optimal control approach

Sébastien Imperiale, Jessica Manganotti, Philippe Moireau

► To cite this version:

Sébastien Imperiale, Jessica Manganotti, Philippe Moireau. Flow recovery from distal pressure in linearized hemodynamics: an optimal control approach. Inverse Problems, 2023, 39 (7), pp.075004. 10.1088/1361-6420/acd274. hal-04178851

HAL Id: hal-04178851

https://inria.hal.science/hal-04178851

Submitted on 8 Aug 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Flow recovery from distal pressure in linearized hemodynamics: An optimal control approach

Sébastien Imperiale, Jessica Manganotti, Philippe Moireau

Inria & Ecole Polytechnique, CNRS, Institut Polytechnique de Paris,
1 rue Honoré d'Estienne d'Orves, 91128, Palaiseau, France

E-mail: sebastien.imperiale@inria.fr

April 2023

**Abstract.** The goal of this work is to derive a reliable – stable and accurate – inverse problem strategy for reconstructing cardiac output – blood flow entering the ascending aorta – from pressure measurements at a distal site of the arterial tree, assumed here to be the descending aorta. We assume that a reduced one-dimensional model of the aorta can be linearized around its steady state, resulting in a wave system with absorbing boundary condition at the outlet. Using this model, we attempt to reconstruct the inlet flow from a pressure measurement at the distal outlet. First, we investigate the observability of the problem and prove that the inversion of the input-output operator for the flow and pressure in the space of time-periodic solutions is ill-posed of degree one. We then develop a variational approach where we minimize the discrepancy between measurements and a simulated state and penalize the error with respect to a periodic state. It is shown that the penalty strategy is convergent and provides an efficient solution for the minimization. Numerical results illustrate the robustness of our approach to noise and the potential of our method to reconstruct inlet flow from real pressure recordings during anesthesia.

# 1. Introduction

During general anesthesia, close monitoring of the patient's hemodynamic state is required to assure a positive surgery outcome [17]. Cardiac output and central aortic pressure – located at the heart's outlet (see Figure 1) – are two markers known to provide important information about the cardiovascular system [7]. Their measurement is achieved via invasive technique and the data can be very noisy. However, the recording at the periphery is affected by lower noise levels and requires less invasive techniques.

For this reason, studies have been performed using available peripheral pressure data to access the hemodynamic state at central sites. The most common strategy is to use a transfer function (TF), usually based on an autoregressive exogenous model [1, 16, 20, 29]. This type of TF has been shown to be more efficient than the one derived by standard Fourier techniques [13]. However, the TF method does not usually
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Qualitative representation of the cardiac output and central aortic pressure curve and their location within the aortic vessel.</figcaption></figure>

account for inter-subject hemodynamic variability [25] and it is not able to provide a good prediction for some important pressure characteristics [13].

To overcome these limitations, we develop an inverse problem (IP) strategy based on a reduced-order (RO) model of the circulation and, in particular, on a one-dimensional (1D) model, since it is able to capture wave propagation effects [30, 32]. This RO model is suitable for representing the larger arteries of the circulation [24, 36], a single vessel [5, 22], or multiple interbranched vessels forming a network [12, 15]. Moreover, they can be coupled with data assimilation techniques [10] to build personalized models and move towards patient-specific simulations [1, 6, 20, 27]. The IP solution leads to the estimation of the physical parameters of the problems reflecting the hemodynamic characteristics of the patient [18, 23, 26] or to the reconstruction of physiological markers, such as fluid velocity in [33].

In this work, we derive a variational IP strategy and apply it to a stable RO model. For the latter, we rely on a stable formulation of a 1D model representing the upper thoracic aorta [22]. Our data assimilation technique is able to estimate the cardiac output as well as the pressure and blood flow at every location of the thoracic aorta, using the distal pressure in the vessel. The introduction of a model provides the opportunity to adapt the parameters to different clinical conditions and to take advantage of the ability of 1D models to describe wave propagation. In this work, we rely on a linearized formulation for blood flow propagation because we focus on the mathematical analysis of the inverse problem and, in particular, prove an observability result. We believe that such an analysis is a fundamental prerequisite for developing an inversion strategy for the nonlinear blood flow model. However, despite linearizing the
original nonlinear model [22] around a steady state, we are able to obtain good estimates using synthetic noisy data and real clinical data registered in vivo during operations as observations for the inverse problem.

In Section 2, we present the choice of configuration for the aortic model and its mathematical and numerical non-linear formulation, as well as the linearization around a reference vessel radius. Section 3 introduces the analysis of the forward problem that leads to important properties such as the existence and uniqueness of a solution. In Section 4 we focus on the inverse problem and we study the observability and we introduce the optimal control approach. Finally, in Section 5, we provide a description of the discretized optimization and we show the results obtained with both synthetic and real in vivo observations.

## 2. Problem Settings

### 2.1. A non-linear hemodynamics model

The model used in this work derives from the non-linear 1D equations that describe the blood flow in a generic arterial segment [14, 31]. Typically, for all $t > 0$ and $s \in (0, L)$ the blood flow $Q(t, s)$ and the cross section $A(t, s)$ are given by

$$ egin{cases} 
\partial_t A + \partial_s Q = 0, \\
\partial_t Q + \alpha_{vp} \partial_s \frac{Q^2}{A} + \frac{A}{\varrho} \partial_s P(A) + k_\nu \frac{Q}{A} = 0, \\
P(A) = P_{ext} + \Psi_e(A) + \Psi_v(A),
\end{cases} 
\quad (1) $$

where $\alpha_{vp}$ is a coefficient for the velocity profile, $\varrho$ the density of blood and $k_\nu$ a dissipative parameter, while $P(A)$ is the pressure inside the vessel, linked to cross-section variations through $\Psi_e$ and $\Psi_v$ that are respectively the elastic and viscoelastic terms of the tube law, and $p_{ext}$ is the external pressure. The model employed in this paper is described in [22], the formulation is based on an appropriate change of variables so that the energy density of the system becomes quadratic with respect to the new unknowns $-\varphi$ and $v$ that represent respectively the square root of the kinetic energy density and the square root of the potential elastic density energy, that are defined as

$$ v := Ru; $$

$$ \varphi := \phi(R) := \begin{cases} \sqrt{\Psi(\pi R^2)} & R \ge R_0, \\ -\sqrt{\Psi(\pi R^2)} & 0 \le R < R_0, \end{cases} \quad \text{with} \quad \Psi(\pi R^2) = \frac{\beta}{A_0} \left[ \frac{2}{3}\pi^{\frac{3}{2}}R^3 - \sqrt{A_0\pi R^2 + \frac{1}{3}A_0^3} \right], $$

where $R = \sqrt{\frac{A}{\pi}}$ is the radius of the vessel, $u = \frac{Q}{A}$ the blood velocity, $\beta$ is a parameter that is linked to the Young modulus of the vessel and $A_0$ the initial cross-section of the
arterial segment. The formulation reads, for all $t \ge 0$ and $s \in (0, L)$

$$ \begin{cases} 2\partial_t\varphi + \pi\xi(R)\partial_s(Rv) = 0, \\ \pi\varrho\partial_t v + \frac{\pi\varrho}{2}\left(2\frac{v}{R}\partial_s v + v\partial_s\frac{v}{R}\right) + \varrho k_\nu \frac{v}{R^2} + \pi R \partial_s(\xi(R)\varphi) = 0, \end{cases} \quad (2) $$

where $\xi(R)$ is a smooth function of $R$. At the inlet and outlet of the segment we use the following boundary conditions,

$$ \left\{ \begin{aligned} \pi R \xi(R) \varphi &= \frac{\pi \varrho}{2} \frac{v^2}{R} - \pi R P_{tot}, & \text{at } s = 0, \\ \pi R \xi(R) \varphi &= -\frac{\pi \varrho}{2} \frac{v^2}{R} + \pi R P_{tot}, & \text{at } s = L, \end{aligned} \right. \quad (3) $$

with $P_{tot}$ being the total pressure. Moreover, we have

$$ \left\{ \begin{aligned} \pi R v &= \nu_{in}, & \text{at } s = 0, \\ P_{tot} &= R_c(\pi R v) + P_c, & \text{at } s = L, \\ C_c \dot{P}_c &= \pi R v - R_{per}^{-1} P_c, & \text{at } s = L, \end{aligned} \right. \quad (4) $$

where $\nu_{in}$ is the input blood flow, $P_c$ the pressure at the outlet of the vessel and $R_c, R_{per}$ and $C_c$ the outlet Windkessel parameters. In (2, 3, 4) the main unknowns are $(v, \varphi, P_c)$ whereas $P_{tot}$ at $s = 0$ acts as a Lagrange multiplier for the constraint $\pi R v = \nu_{in}$. The problem is completed by periodic conditions in time, i.e., $v(\cdot, 0) = v(\cdot, T)$, $\varphi(\cdot, 0) = \varphi(\cdot, T)$ and $P_c(0) = P_c(T)$ where $T$ is a heart-beat period.

## 2.2. A linearized hemodynamics model

Our mathematical analysis is carried out on the linearized version of System (2) around a reference area and zero blood flow, which corresponds to a linearization around $\varphi = 0$ and $v = 0$ – and this is a benefit of the formulation (2) proposed in [22]. The linearized problem associated with the small perturbations $(\hat{\varphi}, \hat{v})$ is a system of standard 1D dissipative wave equations. For all $t > 0$ and $s \in (0, L)$

$$ \left\{ \begin{aligned} 2\partial_t\hat{\varphi} + \pi R_0 \xi(R_0)\partial_s\hat{v} = 0, \\ \pi\varrho\partial_t\hat{v} + \varrho k_\nu \frac{\hat{v}}{R_0^2} + \pi R_0 \xi(R_0)\partial_s\hat{\varphi} = 0, \end{aligned} \right. \quad (5) $$

where $R_0 = \phi^{-1}(0)$. The boundary conditions for (5) are deduced from (3) and (4) and read,

$$ \left\{ \begin{aligned} \pi R_0 \hat{v} &= \nu_{in}, & \text{at } s = 0, \\ \xi(R_0) \hat{\varphi} &= \pi R_c R_0 \hat{v} + \hat{P}_c & \text{at } s = L, \\ C_c \dot{P}_c &= \pi R_0 \hat{v} - R_{per}^{-1} \hat{P}_c & \text{at } s = L. \end{aligned} \right. \quad (6) $$

To ease the presentation, we introduce a set of variables $(a, v, P)$ such that (5, 6) have a simpler form. More precisely we set,

$$ a(s, t) = \sqrt{2} \hat{\varphi}(\alpha s, \beta t), \quad v(s, t) = \sqrt{\pi \varrho} \hat{v}(\alpha s, \beta t), \quad P(t) = \gamma \hat{P}_c(\beta t), $$
and rewrite (5) in the domain $(0, 1) = (0, L/\alpha)$. Then, choosing $\alpha = L$ gives, for all $t > 0$ and $s \in (0, 1)$

$$ \left\{ \begin{aligned} \partial_t a + \frac{\beta \sqrt{\pi} R_0 \xi(R_0)}{L \sqrt{2\varrho}} \partial_s v &= 0, \\ \partial_t v + \frac{\beta k_\nu}{\pi R_0^2} v + \frac{\beta \sqrt{\pi} R_0 \xi(R_0)}{L \sqrt{2\varrho}} \partial_s a &= 0. \end{aligned} \right. \qquad (7) $$

and, for all $t > 0$

$$ \left\{ \begin{aligned} v(0, t) &= \frac{\sqrt{\varrho}}{R_0 \sqrt{\pi}} \nu_{in}(\beta t) =: \nu(t) \\ a(1, t) &= \frac{\sqrt{2\pi} R_c R_0}{\xi(R_0) \sqrt{\varrho}} v(1, t) + \frac{\sqrt{2}}{\gamma \xi(R_0)} P(t), \\ \dot{P}(t) &= \beta \gamma \frac{\sqrt{\pi} R_0}{C_c \sqrt{\varrho}} v(1, t) - \frac{\beta R_{per}^{-1}}{C_c} P(t). \end{aligned} \right. \qquad (8) $$

In order to simplify the re-scaled equation we choose the parameter $\beta$ such that

$$ \frac{\beta \sqrt{\pi} R_0 \xi(R_0)}{L \sqrt{2\varrho}} = 1 \Rightarrow \beta = \frac{L \sqrt{2\varrho}}{\sqrt{\pi} R_0 \xi(R_0)}, $$

and $\gamma$ such that

$$ \frac{\sqrt{2}}{\gamma \xi(R_0)} = \beta \gamma \frac{\sqrt{\pi} R_0}{C_c \sqrt{\varrho}} \Rightarrow \gamma^2 = \frac{C_c}{L}, $$

for symmetry reasons. More precisely, the coupling terms of the second and third equations of (8) are now identical hence simplifying the future energy balance.

### 2.3. Statement of the problem

Following the developments of the previous section we consider the following 1D wave propagation problem: find $a(s, t)$ and $v(s, t)$ in $(0, 1) \times [0, T]$ solution of

$$ \left\{ \begin{aligned} \partial_t a(s,t) + \partial_s v(s,t) &= 0, \\ \partial_t v(s,t) + \partial_s a(s,t) + k v(s,t) &= 0, \end{aligned} \right. \qquad (9) $$

where $k$ accounts for dissipative effects. This system is completed with boundary conditions that involves an input $\nu(t)$ at $s = 0$ as well as a new unknown $P(t) \in \mathbb{R}$ at $s = 1$.

$$ v(0,t) = \nu(t), \quad \left\{ \begin{aligned} a(1,t) &= k_r v(1,t) + c_r P(t), \\ \dot{P}(t) + RP(t) - c_r v(1,t) &= 0. \end{aligned} \right. \qquad (10) $$

Above, $k_r$ and $R$ account for dissipation while $c_r$ is a coupling coefficient. They satisfy

$$ k > 0, \quad k_r > 0, \quad R > 0, \quad c_r \neq 0. \qquad (11) $$

To close the set of equations we consider periodic initial conditions

$$ (a(0), v(0), P(0)) = (a(T), v(T), P(T)), \qquad (12) $$
for some given time period $T$.

We can now define the inverse problem that is the main focus of this paper. Given a noisy measurement of the output pressure $P(t), t \in [0, T]$, can we propose a stable – with respect to the measurement noise – reconstruction of the input flow $v(t)$?

### 3. The forward problem

#### 3.1. Semigroup formalism

Let us first assume that the problem (9, 10) is completed with the initial conditions

$$ (a(s,0), v(s,0), P(s,0)) = (a_0(s), v_0(s), P_0(s)). \qquad (13) $$

Now, to write (9, 10, 13) in an abstract form, we first introduce the state space

$$ \mathcal{Z} = L^2(0, 1) \times L^2(0, 1) \times \mathbb{R}, $$

that is an Hilbert space equipped with the scalar product: for all $z = (a, v, P) \in \mathcal{Z}$ and $q = (b, w, Q) \in \mathcal{Z}$,

$$ (z, q)_{\mathcal{Z}} = (v, w)_{L^2(0,1)} + (v, w)_{L^2(0,1)} + PQ. $$

We then define the operator $(A, \mathcal{D}(A))$ by

$$ A = \begin{pmatrix} 0 & -\partial_s & 0 \\ -\partial_s & -k & 0 \\ 0 & c_r \gamma_r & -R \end{pmatrix} \qquad (14) $$

with $\gamma_r$ the trace operator in $s = 1$ and

$$ \mathcal{D}(A) := \{ (a, v, P) \in H^1(0, 1) \times H^1(0, 1) \times \mathbb{R} \mid a(1) = k_r v(1) + c_r P \}, (15) $$

where $H^1(0, 1) = \{w \in H^1(0, 1) \mid w(0) = 0\}$.

**Proposition 3.1** *The operator $(A, \mathcal{D}(A))$ is maximal dissipative.*

**Proof:** Proofs of such result are rather classical, however the maximal property is subject to technicalities to invert a right hand side in $\mathcal{Z}$. To this end, we propose a proof based on T-coercivity [8] in Appendix A ■

We here recall that as $A$ is maximal dissipative, it is the generator of $C^0$-semigroup of contraction $\Phi: \mathbb{R}^+ \ni t \mapsto e^{tA} \in \mathcal{L}(\mathcal{Z})$, allowing to define solutions of

$$ \left\{ \begin{array}{ll} \dot{z} = Az, & \text{in } (0, T), \\ z(0) = z_0, & \end{array} \right. \qquad (16) $$

that corresponds to solutions of (9, 10, 13). If $z_0 \in \mathcal{D}(A)$, then $z \in C^1([0, T]; \mathcal{Z}) \cap C^0([0, T]; \mathcal{D}(A))$ defined by

$$ z(t) = e^{tA}z_0 \qquad (17) $$
will be called a strict solution of the homogeneous problem. If $z_0 \in \mathcal{Z}$, then $z \in C^0([0, T]; \mathcal{Z})$ defined by (17) will be called a mild solution of the homogeneous problem.

In a second step, we need to characterize the adjoint operator of $A$ that is well defined since $\mathcal{D}(A)$ is dense in $\mathcal{Z}$. Given the space

$$ \mathcal{D}(A^*) = \left\{ q \in \mathcal{Z} \left| \sup_{z \in \mathcal{D}(A), z \neq 0} \frac{|(Az, q)_z|}{\|z\|_{\mathcal{Z}}} < +\infty \right. \right\}, \quad (18) $$

the adjoint operator $A^*$ is the unique operator such that $(Az, q)_z = (z, A^*q)_z$, for all $z \in \mathcal{D}(A)$ and all $q \in \mathcal{D}(A^*)$.

**Lemma 3.2** We have

$$ A^* = \begin{pmatrix} 0 & \partial_s & 0 \\ \partial_s & -k & 0 \\ 0 & -c_r \gamma_r & -R \end{pmatrix} \qquad (19) $$

with

$$ \mathcal{D}(A^*) = \{(b, w, Q) \in H^1(0, 1) \times H^1(0, 1) \times \mathbb{R} \mid b(1) = -k_r w(1) + c_r Q\} $$

**Proof:** Choosing $z = (a, v, P) \in \mathcal{D}(A)$ and $q = (b, w, Q) \in [H^1(0, 1) \times H^1(0, 1) \times \mathbb{R}]$, we obtain, using (Appendix A), that $(Az, q)_z = ((z, q)) + \langle z, q \rangle$ with

$$ ((z, q)) = (b', v)_{L^2(0,1)} + (w', a)_{L^2(0,1)} - k(w, v)_{L^2(0,1)} - RPQ, $$

and, since $z \in \mathcal{D}(A)$,

$$ \langle z, q \rangle = -b(1)v(1) - w(1)a(1) + w(0)a(0) + c_r v(1)Q \\ = -v(1)(b(1) + k_r w(1) - c_r Q) - c_r w(1)P + w(0)a(0). $$

Since $((\cdot, q))$ is obviously a continuous linear form on $\mathcal{Z}$ for all $q \in [H^1(0, 1) \times H^1(0, 1) \times \mathbb{R}]$, the space $\mathcal{D}(A^*)$ should be chosen as an appropriate subspace such that $\langle \cdot, q \rangle$ is also a continuous linear form on $\mathcal{Z}$. It can then be proven that this implies

$$ w(0) = 0, \quad b(1) = -k_r w(1) + c_r Q, $$

hence fixing $\mathcal{D}(A^*)$. Then we get that $\langle z, q \rangle = -c_r w(1)P$ which implies that for all $(z, q) \in \mathcal{Z} \times \mathcal{D}(A^*)$,

$$ (z, A^*q)_z = ((z, q)) - c_r w(1)P, \qquad (20) $$

which defines $A^*$ as expected.

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

Note that using the adjoint definition, we have that any mild solution $z \in C^0([0, T]; \mathcal{Z})$ of the homogeneous problem (16) also belongs to $C^1([0, T]; \mathcal{D}(A^*)')$, see [35, Thm. 4.1.6].
## 3.2. Non-homogeneous problem

To deal with the non-homogeneous problem, that is the problem when considering $\nu$ non-zero, we first need to introduce a lifting of $\nu$. To do so, let us consider the solution $(a_\nu, v_\nu, P_\nu) \in H^1(0, 1) \times H^1(0, 1) \times \mathbb{R}$ of

$$ 
\left\{
\begin{array}{ll}
v'_\nu(s) = 0, & s \in (0, 1) \\
a'_\nu(s) + kv_\nu(s) = 0, & s \in (0, 1) \\
a_\nu(1) = k_r v_\nu(1) + c_r P_\nu, & \\
RP_\nu = c_r v_\nu(1), & \\
v_\nu(0) = \nu. &
\end{array}
\right.
\qquad (21) $$

namely

$$ 
\left\{
\begin{array}{ll}
a_\nu(s) = \left(k(1-s) + k_r + \frac{c_r}{R}\right) \nu, & s \in (0, 1) \\
v_\nu(s) = \nu, & s \in (0, 1) \\
P_\nu = \frac{c_r}{R}\nu. &
\end{array}
\right.
\qquad (22) $$

This allows to introduce the Dirichlet mapping

$$ D_\ell: \left| \begin{array}{lcl} \mathbb{R} & \rightarrow & \mathcal{Z} \\ \nu & \mapsto & z = (a_\nu, v_\nu, P_\nu) \end{array} \right. \qquad (23) $$

In the spirit of [3], giving $\nu \in H^1(0, T)$, a solution $z = (a, v, P)$ of (9,10,13) can then be defined from

$$ \tilde{z} = z - D_{\ell}\nu, \qquad (24) $$

such that $\tilde{z} = (\tilde{a}, \tilde{v}, \tilde{P})$ is solution of

$$ 
\left\{
\begin{array}{ll}
\partial_t \tilde{a}(s,t) + \partial_s \tilde{v}(s,t) = -\partial_t a_\nu(s,t), & (s,t) \in (0,1) \times (0,T) \\
\partial_t \tilde{v}(s,t) + \partial_s \tilde{a}(s,t) + k \tilde{v}(s,t) = -\partial_t v_\nu(s,t), & (s,t) \in (0,1) \times (0,T) \\
\dot{\tilde{P}}(t) + R\tilde{P}(t) - c_r \tilde{v}(1,t) = -\dot{P}_\nu(t), & t \in (0,T)
\end{array}
\right.
$$

with the boundary conditions and the initial conditions

$$ 
\left\{
\begin{array}{ll}
\tilde{v}(0,t) = 0, & t \in (0, T) \\
\tilde{a}(1,t) = k_r \tilde{v}(1,t) + c_r \tilde{P}(t), & t \in (0, T)
\end{array}
\right.
\qquad
\left\{
\begin{array}{ll}
\tilde{a}(s,0) = a_0 - a_{\nu(0)}, & s \in (0,1) \\
\tilde{v}(s,0) = v_0 - v_{\nu(0)}, & s \in (0,1) \\
\tilde{P}(0) = P_0 - P_{\nu(0)}.
\end{array}
\right.
$$

Namely, $\tilde{z}$ is solution of

$$ 
\left\{
\begin{array}{ll}
\dot{\tilde{z}} = A\tilde{z} + D_\ell\dot{v}(t), & t \in (0, T), \\
\tilde{z}(0) = z_0 - D_\ell\nu(0), &
\end{array}
\right.
\qquad (25) $$

Following [3], we have used in (25) the following abuse of notation: the operator $A$ is in fact an extension of $(A, D(A))$ as an operator in $D(A^*)'$ with a dense domain that is identified with $\mathcal{Z}$, hence the two equations of (25) should be understood in $D(A^*)'$.

The standard semi-group theory can be used to show existence and uniqueness result for (25). Namely, when $z_0 - D_\ell\nu(0) \in \mathcal{Z}$ (note that it suffices that $z_0 \in \mathcal{Z}$) and
$\nu \in H^2(0, T)$, there exists one and only one solution $\tilde{z} \in C^0([0, T]; \mathcal{Z}) \cap C^1([0, T]; \mathcal{D}(A^*)')$ obtained from the Duhamel formulae, namely for all $t \in [0, T]$,

$$ \begin{aligned} \tilde{z}(t) &= e^{tA}\tilde{z}_0 + \int_0^t e^{(t-s)A} D_{\ell} \dot{\nu}(s) \, ds \\ &= e^{tA}z_0 - e^{tA}D_{\ell}\nu(0) + [e^{(t-s)A}D_{\ell}\nu(s)]_0^t - \int_0^t e^{(t-s)A}AD_{\ell}\nu(s) \, ds \\ &= e^{tA}z_0 - D_{\ell}\nu(t) - \int_0^t e^{(t-s)A}AD_{\ell}\nu(s) \, ds. \end{aligned} $$

Therefore, when $z_0 \in \mathcal{Z}$ and $\nu \in H^2(0, T)$, there exists a unique solution of (9, 10, 13), it satisfies,

$$ z = (a, v, P) \in C^0([0, T]; \mathcal{Z}) \cap C^1([0, T]; \mathcal{D}(A^*)') $$

and is expressed by the Duhamel formulae

$$ \forall t \in [0, T], \quad z(t) = e^{tA}z_0 - \int_0^t e^{(t-s)A} AD_{\ell}\nu(s) \, ds, \qquad (26) $$

hence is solution of the dynamics

$$ \left\{ \begin{array}{ll} \dot{z} = Az + B\nu, & \text{in } (0, T) \\ z(0) = z_0, & \end{array} \right. \qquad (27) $$

where $B = -AD_{\ell} \in \mathcal{L}(\mathbb{R}, \mathcal{D}(A^*)')$ should be understood in the following sense: for all $\nu \in \mathbb{R}$ and $q = (b, w, Q) \in \mathcal{D}(A^*)$

$$ \begin{aligned} -\langle B\nu, q \rangle_{\mathcal{D}(A^*)', \mathcal{D}(A^*)} &= \langle AD_{\ell}\nu, q \rangle_{\mathcal{D}(A^*)', \mathcal{D}(A^*)} \\ &= (D_{\ell}\nu, A^*q)_{\mathcal{Z}} \\ &= (b', v_{\nu})_{\mathbb{L}^2(0,1)} + (w', a_{\nu})_{\mathbb{L}^2(0,1)} - k(w, v_{\nu})_{\mathbb{L}^2(0,1)} \\ &\quad - c_r w(1)P_{\nu} - RP_{\nu}Q \\ &= b(1)v_{\nu}(1) - b(0)v_{\nu}(0) + w(1)a_{\nu}(1) - w(0)a_{\nu}(0) \\ &\quad - c_r w(1)P_{\nu} - RP_{\nu}Q. \end{aligned} $$

where, to obtain the last line, we have used integration by parts and (21). Using again (21), this can be further simplified to get

$$ \langle B\nu, q \rangle_{\mathcal{D}(A^*)', \mathcal{D}(A^*)} = b(0)\nu, $$

giving additionally $B^* = (\delta_0, 0, 0) \in \mathcal{L}(\mathcal{D}(A^*), \mathbb{R})$ where $\delta_0$ is the trace operator in $s = 0$. As $B \in \mathcal{L}(\mathbb{R}, \mathcal{D}(A^*)')$, we can introduce the operator

$$ \Lambda \in \mathcal{L}(L^2(0, T)), C^0([0, T]; \mathcal{D}(A^*)') $$

defined by

$$ (\Lambda\nu)(t) = - \int_0^t e^{(t-s)A} B\nu(s) \, ds. \qquad (28) $$

**Proposition 3.3** The operator $B$ is an admissible control operator for the semi-group associated with $A$, that is to say,

$$ \Lambda \in \mathcal{L}(L^2(0, T)), C^0([0, T]; \mathcal{Z})). \qquad (29) $$
**Proof:** Following [35, Theorem 4.4.3], we are going to prove that $B^*$ is an admissible observation operator for $\Phi^*$, the semigroup associated with $A^*$. Let us denote by $q_0 \in \mathcal{D}(A^*)$, and $q \in C^1([0, T]; \mathcal{Z}) \cap C^0([0, T]; \mathcal{D}(A^*))$ the solution of $\dot{q} = A^*q$, $q(0) = q_0$. The adjoint operator $A^*$ is also maximal dissipative, hence $q = (b, w, Q) \in C^1([0, T; \mathcal{Z}]) \cap C^0([0, T; \mathcal{D}(A^*)])$ is solution to

$$ egin{cases} 
\partial_t b(s,t) - \partial_s w(s,t) = 0, & (s,t) \in (0,1) \times (0,T) \\
\partial_t w(s,t) - \partial_s b(s,t) + k w(s,t) = 0, & (s,t) \in (0,1) \times (0,T) \\
\dot{Q}(t) + RQ(t) + c_r w(r,t) = 0, & t \in (0,T)
\end{cases} 
\quad (30) $$

with boundary conditions

$$ egin{cases} 
w(0,t) = 0, & t \in (0,T) \\ 
b(1,t) = -k_r w(1,t) + c_r Q(t). & t \in (0,T) 
\end{cases} 
\quad (31) $$

Let us introduce a multiplier $m \in C^1([0, 1])$ to be specified later. We multiply the first equation of (30) by $m(s)w(s,t)$ and the second one by $m(s)b(s,t)$. We integrate with respect to time and space, after summation, we get

$$ -\frac{1}{2} \int_0^T \int_0^1 m \partial_s (b^2 + w^2) \, ds \, dt = - \int_0^T \int_0^1 m \partial_t (b w) \, ds \, dt + \int_0^T \int_0^1 m b w \, ds \, dt. $$

Using by parts integration and Young's inequality we arrive at

$$ -\frac{1}{2} \int_0^T [m(b^2 + w^2)]_0^1 dt \leq \mathcal{F}_m(T) + \mathcal{F}_m(0) + \int_0^T \mathcal{F}_m(t) + \mathcal{F}_{|m'|}(t) dt $$

where, for any smooth function $n \in C^0([0, 1])$ (here $n = m$ or $n = |m'|$),

$$ \mathcal{F}_n(t) = \frac{1}{2} \int_0^1 n(s) (b^2(s,t) + w^2(s,t)) \, ds. $$

Now choosing $m(s) = 1 - s$ and denoting $\mathcal{F}(t)$ the functional $\mathcal{F}_m(t)$ with $m \equiv 1$ we obtain that

$$ \frac{1}{2} \int_0^T b^2(0,t) \, dt \leq \mathcal{F}(T) + \mathcal{F}(0) + 2 \int_0^T \mathcal{F}(t) \, dt. $$

Since one can easily shows by standard energy estimates that $\mathcal{F}(t) \leq \frac{1}{2}\|q_0\|_Z^2$ for all $t > 0$, we finally obtain

$$ \int_0^T b^2(0,t) \, dt \leq 2(1+T) \|q_0\|_Z^2. $$

Moreover, since $b(0,t) = B^*\Phi^*(t)q_0$ we have shown that, for all $T > 0$, there exists a constant $c_{st}(T)$ such that

$$ \forall q_0 \in \mathcal{D}(A^*), \quad \int_0^T |B^*\Phi^*(t)q_0|^2 \, dt \leq c_{st}(T) \|q_0\|_Z^2 $$

ensuring that $B^*$ is an admissible control operator for $\Phi^*$, hence $B$ is an admissible observation operator for $\Phi$. ■
We finally obtain the existence result for our problem (9, 10, 13) by the mean of the abstract dynamics (27). When $\nu \in L^2(0, T)$, and $z_0 \in \mathcal{Z}$, we have from [35, Proposition 4.2.5] the existence of a solution

$$z = (a, v, P) \in C^0([0, T]; \mathcal{Z})$$

given by (26), i.e.,

$$\forall t \in [0, T], \quad z(t) = e^{tA}z_0 + (\Lambda\nu)(t). \qquad (32)$$

When $\nu \in H^1(0, T)$, and $z_0 \in \mathcal{Z}$ satisfies,

$$Az_0 + B\nu(0) \in \mathcal{Z},$$

then, the solution of (27) satisfies

$$z = (a, v, P) \in C^1([0, T]; \mathcal{Z}). \qquad (33)$$

**Corollary 3.4** When $\nu \in H^1(0, T)$, $z_0 = (a_0, v_0, P_0) \in H^1(0, 1) \times H^1(0, 1) \times \mathbb{R}$ and

$$v_0(0) = \nu(0), \quad a_0(1) = k_r v_0(1) + c_r P_0, \qquad (34)$$

then the unique solution of (27) satisfies

$$z = (a, v, P) \in C^1([0, T]; \mathcal{Z}) \cap C^0([0, T]; H^1(0, 1) \times H^1(0, 1) \times \mathbb{R}),$$

and the boundary conditions

$$v(0, t) = \nu(t), \quad a(1, t) = k_r v(1, t) + c_r P(t).$$

This solution is called the strict solution of the non-homogeneous problem (27)

**Proof:** The first step is to apply [35, Proposition 4.2.10] to prove that the unique solution of (27) belongs to $C^1([0, T]; \mathcal{Z})$. To do so, one should check that $Az_0 + B\nu(0)$ – a priori an element of $\mathcal{D}(A^*)'$ – can be identified with an element of $\mathcal{Z}$. We have, for all $q = (b, w, Q) \in \mathcal{D}(A^*)$, using (20) and (28),

$$\langle Az_0 + B\nu(0), q \rangle_{\mathcal{D}(A^*)', \mathcal{D}(A^*)} = (z_0, A^*q)_\mathcal{Z} + b(0)\nu(0) \\ = ((z, q)) - c_r w(1)P_0 + b(0)\nu(0).$$

Integrating by parts yields

$$\langle Az_0 + B\nu(0), q \rangle_{\mathcal{D}(A^*)', \mathcal{D}(A^*)} = -(b, v_0')_{L^2(0,1)} - (w, a_0')_{L^2(0,1)} - k(w, v)_{L^2(0,1)} \\ + b(1)v_0(1) - b(0)v_0(0) + w(1)a_0(1) \\ - RPQ - c_r w(1)P_0 + b(0)\nu(0),$$

using now that $b(1) = -k_r w(1) + c_r Q$ – since $q \in \mathcal{D}(A^*)$ – and (34), one can simplify the equation above,

$$\langle Az_0 + B\nu(0), q \rangle_{\mathcal{D}(A^*)', \mathcal{D}(A^*)} = -(b, v_0')_{L^2(0,1)} - (w, a_0')_{L^2(0,1)} - k(w, v)_{L^2(0,1)} \\ - RPQ + c_r Qv_0(1).$$
This shows that $Az_0+B\nu(0)$ can be identified with an element of $\mathcal{Z}$, hence the regularity $C^1([0, T]; \mathcal{Z})$ of the solution. This last regularity property implies in particular that $Az+B\nu$ can be identified with an element of $C^0([0, T]; \mathcal{Z})$, $\dot{z}$ to be specific, and with the same computations as before, but using a function $q \in \mathcal{D}(A^*)$ with $(b, w)$ with compact support and $Q = 0$,

$$ \langle Az + B\nu, q \rangle_{\mathcal{D}(A^*)', \mathcal{D}(A^*)} = (\dot{z}, q)_\mathcal{Z} = (b, \dot{v})_{L^2(0,1)} + (w, \dot{a})_{L^2(0,1)} \\ = -(b, v')_{L^2(0,1)} - (w, a')_{L^2(0,1)} - k(w, v)_{L^2(0,1)}. $$

This equality implies that (9) are satisfied in a strong sense, hence $(a, v)$ belong to $C^0([0, T]; H^1(0, 1) \times H^1(0, 1))$. Now, for any, $q \in \mathcal{D}(A^*)$ with $Q = 0$, we have

$$ (\dot{z}, q)_\mathcal{Z} = ((z, q)) - c_r w(1)P + b(0)\nu(t). $$

Integration by part yields, after simplifications,

$$ 0 = -k_r w(1)v(1, t) - b(0)v(0, t) + w(1)a(1, t) - c_r w(1)P(t) + b(0)\nu(t). $$

Choosing $Q = 0$ and $w(1) = 0$ we deduce that $v(0, t) = \nu(t)$. Then, choosing again $Q = 0$ we deduce that $a(1, t) = k_r v(1, t) + c_r P(t)$. ■

### 3.3. Periodic solutions

We now expect to prove an additional fundamental property for the homogeneous problem (16), namely that the associated semigroup satisfies for all $t \in [0, T]$, $\|\Phi(t)\|_{\mathcal{L}(\mathcal{Z})} < \rho < 1$, which will later be used to prove the existence of periodic solutions for our model (9, 10, 13) or equivalently for the abstract dynamic (27).

The energy functional. Let $m \in C^1([0, 1])$ be a smooth function, We denote by $\mathcal{E}(t)$ the energy of system (9), defined by

$$ \forall t > 0, \quad \mathcal{E}(t) := \frac{1}{2} \int_0^1 (a^2(s, t) + v^2(s, t)) \, ds + \frac{1}{2} P^2(t). $$

We have the so-called energy identity.

**Lemma 3.5** Any strict solution $z$ in the sense of Corollary 3.4 satisfies for all $t \in [0, T]$,

$$ \begin{aligned} \mathcal{E}(t) + k \int_0^t \int_0^1 v^2(s, \tau) \, ds \, d\tau + \int_0^t k_r v^2(1, \tau) \, d\tau + \int_0^t RP^2(\tau) \, d\tau \\ = \mathcal{E}(0) + \int_0^t \nu(\tau) a(0, \tau) \, d\tau. \end{aligned} \quad (35) $$

**Proof:** We multiply the first equation of (9) by $a(s, t)$ and the second one by $v(s, t)$. We integrate with respect to time and space, after summation, and integration by part
we get

$$ \begin{aligned} & \frac{1}{2} \int_{0}^{1} a^{2}(s,t) + v^{2}(s,t) \, ds + \int_{0}^{t} \int_{0}^{1} k v^{2}(s, \tau) \, ds \, d\tau + \int_{0}^{t} v(1, \tau) a(1, \tau) \, d\tau \\ & = \frac{1}{2} \int_{0}^{1} a^{2}(s, 0) + v^{2}(s, 0) \, ds + \int_{0}^{t} v(0, \tau) a(0, \tau) \, d\tau. \end{aligned} \quad (36) $$

Then, we multiply the second equation of (10) by $P(t)$ and integrate with respect to time to obtain

$$ \frac{1}{2}(P^2(t) - P^2(0)) + \int_{0}^{t} (R P^2(\tau) - c_r v(1, \tau) P(\tau)) \, d\tau = 0. \quad (37) $$

Substituting (37) in (36) gives (35). ■

Note that the previous lemma, justifies the classical energy balance of (16), namely, we have, in absence of input source $v$,

$$ \forall t > 0, \quad \mathcal{E}(t) \le \mathcal{E}(0), \quad (38) $$

consistent with the dissipative aspect of $(A, \mathcal{D}(A))$. In fact we have an even stronger property.

**Theorem 3.6** Assume that $\nu = 0$, then for every $t_0 > 0$ there exists $\rho < 1$ such that, any mild solution $z = (a, v, P) \in C^0([0, T]; \mathcal{Z})$ of (16) satisfies,

$$ \forall t \ge t_0 > 0, \quad \mathcal{E}(t) \le \rho \mathcal{E}(0). \quad (39) $$

**Proof:** We first show that (39) is true for strict solutions in $C^1([0, T]; \mathcal{Z}) \cap C^0([0, T]; \mathcal{D}(A))$. The result is then extended to mild solution by density arguments. The proof of (39) for strict solutions is obtained by contradiction. If (39) does not hold there exists a sequence of solutions $\{(a_n, v_n, P_n)\}_{n>0}$ such that the corresponding energy functional $\mathcal{E}_n(t)$ satisfies for some $t_1 > 0$,

$$ \mathcal{E}_n(0) = 1 \quad \text{and} \quad \lim_{n \to +\infty} \mathcal{E}_n(t_1) = 1. \quad (40) $$

Note that the energy relation (35) holds since strict solutions of the homogeneous problem are in particular strict solutions of the non-homogeneous problem with $\nu = 0$. Therefore, since $\mathcal{E}_n$ is a decaying function of time, we deduce that, thanks to the energy relation (35), that for all $\tau \le t_1$,

$$ \lim_{n \to +\infty} \left( \int_{0}^{\tau} \int_{0}^{1} v_n^2(s,t) \, ds \, dt + \int_{0}^{\tau} v_n^2(1,t) \, dt + \int_{0}^{\tau} P_n^2(t) \, dt \right) = 0. \quad (41) $$

Moreover, using the boundary conditions (10), we also have, for all $\tau \le t_1$,

$$ \lim_{n \to +\infty} \int_{0}^{\tau} a_n^2(1, t) \, dt = 0. \quad (42) $$

Let us then rely on a multiplier strategy by multiplying the first equation of (9) by $stv(s,t)$ and the second one by $sta(s,t)$ and integrating with respect to time and space, we get after summation

$$ \int_{0}^{\tau} \int_{0}^{1} \left( \partial_{t}(stav) + st \partial_{s} \frac{v^{2}}{2} + st \partial_{s} \frac{a^{2}}{2} + s(kt - 1)av \right) ds dt = 0. $$
Integrating by parts and rearranging terms shows that, for all strict solutions of (9) with $\nu = 0$ and all $\tau > 0$,

$$ \frac{1}{2} \int_{0}^{\tau} \int_{0}^{1} t a^2 \, ds \, dt = \frac{1}{2} \int_{0}^{\tau} (v^2(1,t) + a^2(1,t)) \, dt. \\ + \int_{0}^{1} s \tau a v(s, \tau) \, ds + \frac{1}{2} \int_{0}^{\tau} \int_{0}^{1} (2(kt-1) s a v - t v^2) \, ds \, dt. \quad (43) $$

Thanks to (43) applied to our sequence of solutions $\{(a_n, v_n, P_n)\}_{n>0}$ and thanks to (41) and (42) we deduce that, for all $\tau \le t_1$,

$$ \liminf_{n \to +\infty} \int_{0}^{\tau} \int_{0}^{1} t a_n^2 \, ds \, dt = 0. $$

This implies that

$$ \liminf_{n \to +\infty} \int_{t_1/2}^{t_1} \int_{0}^{1} a_n^2 \, ds \, dt = 0, $$

hence, with (41) and (42),

$$ \liminf_{n \to +\infty} \int_{t_1/2}^{t_1} \mathcal{E}_n(t) \, dt = 0. \quad (44) $$

However, $\lim_{n \to +\infty} \mathcal{E}_n(t_1) = 1$ and the decay property of the energy functional $\mathcal{E}_n(t)$, implies that

$$ \liminf_{n \to +\infty} \int_{t_1/2}^{t_1} \mathcal{E}_n(t) \, dt \ge \liminf_{n \to +\infty} \mathcal{E}_n(t_1) \int_{t_1/2}^{t_1} \, dt = \frac{t_1}{2}. $$

which contradicts (44) and ends the proof. ■

*Existence and uniqueness of a periodic solution.* We conclude this section by obtaining that $\Phi$ satisfies, thanks to Theorem 3.6, for $T > 0$

$$ \exists \rho < 1, \quad \forall z_0 \in \mathcal{Z}, \quad \|e^{TA}z_0\|_{\mathcal{Z}} \le \rho \|z_0\|_{\mathcal{Z}}. \quad (45) $$

Then, considering the non-homogeneous dynamics (27), we recall the Duhamel formulae that

$$ z(T) = e^{TA}z(0) + (\Lambda\nu)(T), $$

such that seeking a T-periodic solution imposes

$$ (\mathrm{Id} - e^{TA})z(0) = (\Lambda\nu)(T). $$

Using the strict contraction property (45), we have that $(\mathrm{Id} - e^{TA})$ is invertible with a bounded inverse giving by the Neumann series

$$ \forall z \in \mathcal{Z}, \quad (\mathrm{Id} - e^{TA})^{-1}z = \sum_{k \ge 0} e^{kTA}z. $$
Therefore, giving $\nu \in L^2(0, T)$, there exists one and only one T-periodic solution of (27)

$$z \in C^0([0, T]; \mathcal{Z}) \cap C^1([0, T]; \mathcal{D}(A^*)'), \quad z(0) = z(T),$$

given for $t \in [0, T]$ by

$$z(t) = e^{tA}(\mathrm{Id} - e^{TA})^{-1}(\Lambda\nu)(T) + (\Lambda\nu)(t). \quad (46)$$

This solution is called the mild periodic solution of (27) When $\nu$ is smoother one can show that solutions are also smooth. To state this results we introduce the space

$$H_{\#}^{1}(0, T) := \{u \in H^{1}(0, T) | u(T) = u(0)\}.$$

**Theorem 3.7** For $T > 0$ given, assume that $\nu \in H_{\#}^{1}(0, T)$, then the unique T-periodic solution of (27) satisfies

$$z = (a, v, P) \in C^1([0, T]; \mathcal{Z}) \cap C^0([0, T]; H^1(0, 1) \times H^1(0, 1) \times \mathbb{R}),$$

and the boundary conditions

$$v(0, t) = \nu(t), \quad a(1, t) = k_r v(1, t) + c_r P(t). \quad (47)$$

Such solution is called a strict periodic solution of (27).

**Proof:** It suffices to show that $z(0)$ satisfies $Az(0) + B\nu(0) \in \mathcal{Z}$ where $z(0)$ is given by

$$Az(0) = -A(\mathrm{Id} - e^{TA})^{-1}(\Lambda\nu)(T) = -(\mathrm{Id} - e^{TA})^{-1}A(\Lambda\nu)(T).$$

One can show, using the periodicity of $\nu$, that

$$\begin{align*} (\Lambda\nu)(T) &= \int_{0}^{T} e^{(T-s)A} AD_{\ell}\nu(s) \, ds \\ &= -\int_{0}^{T} e^{(T-s)A} D_{\ell}\dot{\nu}(s) \, ds + \left[e^{(T-s)A} D_{\ell}\nu(T)\right]_{0}^{T}, \\ &= -\int_{0}^{T} e^{(T-s)A} D_{\ell}\dot{\nu}(s) \, ds + (\mathrm{Id} - e^{TA}) D_{\ell}\nu(0). \end{align*}$$

The computations above are valid in $\mathcal{D}(A^*)'$. Moreover, in the final inequality, each terms belongs to $\mathcal{Z}$. Therefore, We can apply the operator $A$ – seen as an operator from $\mathcal{Z}$ to $\mathcal{D}(A^*)'$ – to get

$$\begin{align*} A(\Lambda\nu)(T) &= -\int_{0}^{T} e^{(T-s)A} AD_{\ell}\dot{\nu}(s) \, ds + (\mathrm{Id} - e^{TA})AD_{\ell}\nu(0) \\ &= (\Lambda\dot{\nu})(T) + (\mathrm{Id} - e^{TA})B\nu(0). \end{align*}$$

We have shown that

$$Az(0) + B\nu(0) = -(\mathrm{Id} - e^{TA})^{-1}(\Lambda\dot{\nu})(T).$$

Thanks to the admissibility of $B$, the operator $\Lambda$ is continuous from $L^2(0, T)$ to $C^0([0, T]; \mathcal{Z})$. Therefore, the last term in the previous identity is in fact an element of $\mathcal{Z}$.

<figure><img src="image_4.png" /></figure>
# 4. The inverse problem

## 4.1. The input-output operator

The observed quantity is the output pressure $P$. We therefore introduce the bounded observation operator $C \in \mathcal{L}(\mathcal{Z}, \mathbb{R})$.

$$ C: \left| \begin{array}{l} \mathcal{Z} \to \mathbb{R} \\ (a, v, P) \mapsto P \end{array} \right. $$

and the input-output linear operator

$$ \Psi_T: \left| \begin{array}{l} L^2(0, T) \to L^2(0, T) \\ \nu \mapsto [(0, T) \ni t \mapsto C z_{|\nu}(t) \in \mathbb{R}] \end{array} \right. $$

where $z_{|\nu}$ is solution of

$$ \left\{ \begin{array}{l} \dot{z}_{|\nu} = Az_{|\nu} + B\nu, \quad \text{in } (0, T) \\ z_{|\nu}(0) = z_{|\nu}(T) \end{array} \right. \qquad (48) $$

Our objective is then to invert $\Psi_T$ which, from (46), is explicitly given by

$$ (\Psi_T \nu)(t) = C e^{tA} (\mathrm{Id} - e^{tA})^{-1} (\Lambda \nu)(T) + C(\Lambda \nu)(t). $$

As we prove that $B$ is admissible, (29) and the fact that $(\mathrm{Id} - e^{AT})^{-1}$ is bounded implies that $\Psi_T \in \mathcal{L}(L^2(0, T))$.

## 4.2. Observability

We must now quantify what can be reconstructed from the measurements through stability estimates similar to observability inequality in control theory. This will be done with two similar results where the source term is controlled by the output pressure, with two choices of norms.

**Theorem 4.1** There exists a constant $c_{st}$ such that the strict periodic solution of (27) satisfies

$$ \int_{0}^{T} |\nu(t)|^{2} \mathrm{d}t \leq c_{st} \int_{0}^{T} (|P(t)|^{2} + |\dot{P}(t)|^{2}) \mathrm{d}t. $$

**Proof:** By multiplying the first equation of (9) by $e^{\alpha s}v(s, t)$ and the second one by $e^{\alpha s}a(s, t)$ and integrating with respect to time and space, after summation, we get

$$ \int_{0}^{T} \int_{0}^{1} \left( \partial_{t} (e^{\alpha s} a v) + e^{\alpha s} \partial_{s} \frac{v^{2}}{2} + e^{\alpha s} \partial_{s} \frac{a^{2}}{2} + k e^{\alpha s} a v \right) \mathrm{d}s \mathrm{d}t = 0. $$

Integrating by parts, using the periodicity, and rearranging terms gives

$$ \begin{aligned} & \frac{1}{2} \int_{0}^{T} (v^{2}(0, t) + a^{2}(0, t)) \mathrm{d}t - \frac{1}{2} \int_{0}^{T} (v^{2}(1, t) + a^{2}(1, t)) e^{\alpha} \mathrm{d}t \\ & = \frac{1}{2} \int_{0}^{T} \int_{0}^{1} e^{\alpha s} (2k a v - \alpha v^{2} - \alpha a^{2}) \mathrm{d}s \mathrm{d}\tau. \end{aligned} \qquad (49) $$
We remark that the right hand side of the equation above is negative if $\alpha = k$, so that, using the boundary condition $v(0, t) = \nu(t)$, we can write

$$ \frac{1}{2} \int_{0}^{T} \nu^{2}(t) \, dt \leq \frac{e^{k}}{2} \int_{0}^{T} (\nu^{2}(1, t) + a^{2}(1, t)) \, dt. $$

Moreover, using (47) and the equations satisfied by $P$, we get that there exists a constant $c_{st}$ such that

$$ \frac{1}{2} \int_{0}^{T} \nu^{2}(t) \, dt \leq c_{st} e^{k} \int_{0}^{T} \left[ \left( \frac{1}{c_{r}^{2}} + \frac{k_{r}^{2}}{c_{r}^{2}} \right) (\dot{P} + RP)^{2} + c_{r}^{2} P^{2} + 2k_{r}P(\dot{P} + RP) \right] dt, $$

hence, applying the Young inequality, we obtain that there exists another constant $c_{st}$ such that

$$ \frac{1}{2} \int_{0}^{T} |\nu(t)|^{2} \, dt \leq c_{st} \int_{0}^{T} (|\dot{P}(t)|^{2} + |P(t)|^{2}) \, dt. $$

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

In fact, Theorem 4.1 can be completed with a similar observability inequality but in a weaker norm. In order to state the second result, we need to introduce the following notation for the time average of functions: for every $f \in L^2(0, T)$, we define

$$ \langle f \rangle_T := \frac{1}{T} \int_0^T f(t) \, dt. $$

We also need the following decomposition Lemma.

**Lemma 4.2** For every $f \in L^2(0, T)$ there exists a unique $g = \mathcal{I}_T(f) \in H_\#^1(0, T)$ denoted such that

$$ f = \dot{g} + \langle f \rangle_T \quad \text{and} \quad \langle g \rangle_T = 0. \tag{50} $$

Moreover, there exists a constant $c_{st}$ such that, for any $f \in L^2(0, T)$,

$$ \| \mathcal{I}_T(f) \|_{L^2(0,T)} \leq c_{st} \| f \|_{L^2(0,T)}. \tag{51} $$

**Proof:** To prove the existence, observe that the function $g(t)$ defined by

$$ g(t) = g(0) + \int_{0}^{t} f(\tau) \, d\tau - t \langle f \rangle_{T}, \tag{52} $$

satisfies (50) as soon as

$$ g(0) = - \left\langle \int_{0}^{t} f(\tau) \, d\tau \right\rangle_{T} + \langle t \rangle_{T} \langle f \rangle_{T}. \tag{53} $$

Since the decomposition is linear, uniqueness is proven by showing that $f = 0$ implies $g = 0$. When $f = 0$ we obtain $\partial_t g = 0$ from (52) and since $\langle g \rangle = 0$ we indeed obtain $g = 0$. The estimates (51) is a simple consequence of the explicit formulae (52) and (53).

<figure><img src="image_6.png" /></figure>
For smoother periodic functions an interpolation result can be obtained.

**Theorem 4.3** Every $f \in H_{\#}^1(0, T)$ satisfies

$$ \|f - \langle f \rangle_T\|_{L^2(0,T)}^2 \le \|\dot{f}\|_{L^2(0,T)} \|\mathcal{I}_T(f)\|_{L^2(0,T)}. $$

**Proof:** Let $f \in H_{\#}^1(0, T)$ and denote $\bar{f} = f - \langle f \rangle_T$. We have

$$ \begin{aligned} \int_0^T |\bar{f}(t)|^2 dt &= \left| \int_0^T \bar{f}(t) \left( \frac{d}{dt}(\mathcal{I}_T(\bar{f}))(t) + \underbrace{\langle \bar{f} \rangle_T}_{0} \right) dt \right| \\ &= \left| \left[ \mathcal{I}_T(\bar{f})(t) \bar{f}(t) \right]_0^T - \int_0^T \dot{\bar{f}}(t) \mathcal{I}_T(\bar{f})(t) dt \right| \\ &\le \|\dot{\bar{f}}\|_{L^2(0,T)} \|\mathcal{I}_T(\bar{f})\|_{L^2(0,T)} \\ &\le \|\dot{f}\|_{L^2(0,T)} \|\mathcal{I}_T(f)\|_{L^2(0,T)}. \end{aligned} $$

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

**Theorem 4.4** There exists a constant $c_{st}$ such that, for any $\nu \in L^2(0, T)$, the mild periodic solution of (27) satisfies

$$ \langle \nu \rangle_T^2 + \int_0^T |\mathcal{I}_T(\nu)(t)|^2 dt \le c_{st} \int_0^T |P(t)|^2 dt. \quad (54) $$

**Proof:** We denote by $\tilde{z} = (\tilde{a}, \tilde{v}, \tilde{P})$ the strict periodic solution associated with the source term $\mathcal{I}_T(\nu) \in H_{\#}^1(0, T)$ as well as $\bar{z} = (\bar{a}, \bar{v}, \bar{P})$ the strict periodic solution of associated with the constant source term $\langle \nu \rangle_T$. Applying Theorem 4.1 for each solution $\tilde{z}$ and $\bar{z}$ we have

$$ \langle \nu \rangle_T^2 + \int_0^T |\mathcal{I}_T(\nu)(t)|^2 dt \le c_{st} \int_0^T \left( |\tilde{P}|^2 + |\dot{\tilde{P}}|^2 + \frac{1}{T} |\bar{P}|^2 + \frac{1}{T} |\dot{\bar{P}}|^2 \right) dt. $$

It remains to estimate the right hand side of the previous equation with the norm of $P$. This can be done using the linearity of the problem,

$$ \nu = \frac{d}{dt} \mathcal{I}_T(\nu) + \langle \nu \rangle_T \quad \Rightarrow \quad z = \dot{\tilde{z}} + \bar{z} \quad \Rightarrow \quad P = \dot{\tilde{P}} + \bar{P}. \quad (55) $$

This last equation implies $\bar{P}(t) = \langle P \rangle_T$ hence $\dot{\bar{P}} = 0$. The function $\mathcal{I}_T(\nu)$ having zero average one can easily show that this is also the case for $\tilde{P}$. And therefore by Poincaré-Wirtinger inequality,

$$ \begin{aligned} \int_0^T |\tilde{P}|^2 dt &\le c_{pw}^2 \int_0^T |\dot{\tilde{P}}|^2 dt \\ &\Rightarrow \langle \nu \rangle_T^2 + \int_0^T |\mathcal{I}_T(\nu)(t)|^2 dt \le c_{st} (1 + c_{pw}^2) \int_0^T |\dot{\tilde{P}}|^2 dt + c_{st} \langle P \rangle_T^2. \end{aligned} $$
Estimate (54) can then be deduced from (55) using standard techniques.

As a first direct consequence of Theorem 4.4 and the linearity of $\Psi_T$, we have that $\Psi_T \in \mathcal{L}(L^2(0, T))$ is injective. Indeed (54) is nothing but

$$ \langle \nu \rangle_T^2 + \int_0^T |\mathcal{I}_T(\nu)(t)|^2 \, dt \le c_{st} \|\Psi_T \nu\|_{L^2(0,T)}^2. \quad (56) $$

### 4.3. Least-square approach

We now introduce $y_\delta$ which takes into account the observation corresponding to an unknown noisy solution of System (48). From the previous analysis, we may now decide to invert $\Psi_T$ to find the unknown solution using a least squares minimization including a generalized Tikhonov regularization [11, 28] of the form

$$ \min_{\nu \in H_{\#}^1(0,T)} \left\{ \mathcal{J}_{\kappa}(\nu) = \frac{1}{2} \int_0^T \left[ |y_{\delta}(t) - Cz_{|\nu}(t)|^2 + \kappa(|\dot{\nu}(t)|^2 + |\nu|^2) \right] dt \right\}, (57) $$

with $\kappa > 0$ and where $z_{|\nu}$ is constrained to be solution of (48). For further computations it is useful to rewrite the functional $\mathcal{J}_{\kappa}: H_{\#}^1(0, T) \mapsto \mathbb{R}$, as follows,

$$ \mathcal{J}_{\kappa}(\nu) = \frac{\kappa}{2} \|\nu\|_{H^1(0,T)}^2 + \frac{1}{2} \|y_{\delta} - \Psi_T \nu\|_{L^2(0,T)}^2, \quad (58) $$

which shows – since $\Psi_T \in \mathcal{L}(L^2(0, T))$ is continuous – that $\mathcal{J}_{\kappa}$ is continuous and quadratic.

**Theorem 4.5** For all $\kappa > 0$, $\mathcal{J}_{\kappa}$ is strongly convex, hence admits one unique minimiser $\bar{\nu} \in H_{\#}^1(0, T)$.

**Proof:** The functional $\mathcal{J}_{\kappa}$ is differentiable in the sense of Fréchet with

$$ \mathcal{J}_{\kappa}(\nu_2) - \mathcal{J}_{\kappa}(\nu_1) = \langle D\mathcal{J}_{\kappa}(\nu_1), \nu_2 - \nu_1 \rangle + \sigma(\nu_2 - \nu_1, \nu_2 - \nu_1) $$

where, here, $\langle \cdot, \cdot \rangle$ stands for the duality product in $H_{\#}^1(0, T)$ and where

$$ D\mathcal{J}_{\kappa} \in \mathcal{L}(H_{\#}^1(0, T), H_{\#}^1(0, T))' $$

is the Fréchet derivative and $\sigma$ is a continuous bilinear form in $H_{\#}^1(0, T) \times H_{\#}^1(0, T)$ given by

$$ \sigma(\nu, \mu) = (\Psi_T \nu, \Psi_T \mu)_{L^2(0,T)} + \kappa(\nu, \mu)_{H^1(0,T)}. $$

As $\sigma$ is coercive in $H_{\#}^1(0, T)$, we deduce that $\mathcal{J}_{\kappa}$ is a strongly convex function, hence the existence of one unique minimizer.
**Remark 4.6** The previous results holds even if we replace the H¹-norm by the semi norm $||\dot{\nu}||_{L^2(0,T)}$. Indeed in this case, from (56), we know that there exists a constant $c_{st}$ such that

$$ \sigma(\nu, \bar{\nu}) \geq c_{st} (\langle \nu \rangle_T^2 + \|\dot{\nu}\|_{L^2(0,T)}^2) $$

From Poincaré-Wirtinger inequality, $\sigma$ is coercive in $H_\#^1(0,T)$. Therefore, it will be possible to consider the penalty $\kappa \|\dot{\nu}\|_{L^2(0,T)}^2 + \kappa' \|\nu\|_{L^2(0,T)}^2$ with a vanishing second contribution $\kappa' \to 0$. This will penalize large variations of $\nu$ without penalizing too much large values of $\nu$.

We are now able to prove the convergence of our least-square methods when seeking to recover a regular enough input flow from noisy measurement.

**Theorem 4.7** Given a regular target input flow $\tilde{\nu} \in H^1(0,T)$ such that $||\tilde{\nu}||_{H^1(0,T)} \leq M$ and some observations $y_\delta \in L^2(0,T)$ such that $\|y_\delta - \Psi_T \tilde{\nu}\|_{L^2(0,T)} \leq \delta ||\Psi_T \tilde{\nu}||_{L^2(0,T)}$, then by choosing $\kappa = \delta^2 N^2 M^{-2}$, where $N = ||\Psi_T \tilde{\nu}||_{L^2(0,T)}$, there exists a constant $c_{st}$ independent of $M$, $N$ and $\delta$ such that

$$ ||\tilde{\nu} - \bar{\nu}||_{L^2(0,T)} \leq c_{st}(\sqrt{M} + \sqrt{\delta})\sqrt{N\delta}, \quad (59) $$

where $\bar{\nu} = \arg\min_\nu \mathcal{J}_\kappa$ is the estimated input flow in $H_\#^1(0,T)$.

**Proof:** On the one hand, the error system $\tilde{z} = z_{|\tilde{\nu}} - z_{|\bar{\nu}}$ is solution of

$$ \begin{cases} \dot{\tilde{z}} = A\tilde{z} + B(\tilde{\nu} - \bar{\nu}), & t > 0 \\ \tilde{z}(0) = \tilde{z}(T) \end{cases} $$

hence, using Theorem 4.4, there exists a constant $c_{st}$ such that

$$ |\langle \tilde{\nu} \rangle_T - \langle \bar{\nu} \rangle_T|^2 + \int_0^T |\mathcal{I}_T(\tilde{\nu} - \bar{\nu})(t)|^2 dt \leq c_{st} \int_0^T |\check{P}(t) - \bar{P}(t)|^2 dt. $$

On the other hand, Theorem 4.3 gives that

$$ ||\tilde{\nu} - \bar{\nu}||_{L^2(0,T)} \leq ||\dot{\tilde{\nu}} - \dot{\bar{\nu}}||_{L^2(0,T)}^{\frac{1}{2}} ||\mathcal{I}_T(\tilde{\nu} - \bar{\nu})||_{L^2(0,T)}^{\frac{1}{2}} + |\langle \tilde{\nu} \rangle_T - \langle \bar{\nu} \rangle_T| $$

Combining these two inequality, we obtain that there exists a constant $c_{st}$ such that

$$ ||\tilde{\nu} - \bar{\nu}||_{L^2(0,T)}^2 \leq c_{st} \left( ||\dot{\tilde{\nu}}||_{L^2(0,T)} + ||\dot{\bar{\nu}}||_{L^2(0,T)} \right) \left( \int_0^T |\check{P}(t) - \bar{P}(t)|^2 dt \right)^{\frac{1}{2}} + c_{st} \int_0^T |\check{P}(t) - \bar{P}(t)|^2 dt. $$

Moreover, we have that

$$ ||\dot{\tilde{\nu}}||_{L^2(0,T)}^2 \leq \frac{2}{\kappa} \mathcal{J}_{\kappa}(\bar{\nu}) \leq \frac{2}{\kappa} \mathcal{J}_{\kappa}(\tilde{\nu}) \leq M^2 + \frac{\delta^2 ||\tilde{y}||_{L^2(0,T)}^2}{\kappa} = M^2 + \frac{\delta^2 N^2}{\kappa}, $$

leading when $\kappa = \delta^2 N^2 M^{-2}$ to

$$ ||\dot{\tilde{\nu}}||_{L^2(0,T)} + ||\dot{\bar{\nu}}||_{L^2(0,T)} \leq (1 + \sqrt{2})M. $$
Likewise, since $\check{P}(t) = \Psi_T \check{\nu}$ and $\bar{P}(t) = \Psi_T \bar{\nu}$, we have that

$$ \begin{aligned} \int_0^T |\check{P} - \bar{P}|^2 dt &\le 2 \int_0^T |\Psi_T \check{\nu} - y_\delta|^2 dt + 2 \int_0^T |\Psi_T \bar{\nu} - y_\delta|^2 dt \\ &\le 2\delta^2 N^2 + 4 \mathcal{J}_\kappa(\bar{\nu}) \le 2\delta^2 N^2 + 4 \mathcal{J}_\kappa(\check{\nu}) \le 4\delta^2 N^2 + 2\kappa M, \end{aligned} $$

leading when $\kappa = \delta^2 N^2 M^{-2}$ to

$$ \|\check{P}(t) - \bar{P}(t)\|_{L^2(0,T)} \le \sqrt{6\delta N}. $$

This ultimately gives a constant $c_{st}$ independent of $M$, $N$ and $\delta$ such that

$$ \|\check{\nu} - \bar{\nu}\|_{L^2(0,T)}^2 \le c_{st}(M+N)N\delta, $$

and concludes the proof.

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

### 4.4. Penalization of the periodicity constraint

Minimizing (58) is made difficult by the dynamics constraint with an additional periodicity constraint since (58) could also be understood as the minimization of a functional $\mathcal{J}_\kappa : \mathcal{Z} \times H_\#^1(0, T) \mapsto \mathbb{R}$ given by

$$ \mathcal{J}_\kappa(\zeta, \nu) = \frac{1}{2} \int_0^T \left[ |y_\delta(t) - Cz_{|\zeta,\nu}(t)|^2 + \kappa(|\dot{\nu}(t)|^2 + |\nu|^2) \right] dt, \quad (60) $$

over the convex set

$$ \mathcal{K} = \left\{ (\zeta, \nu) \in \mathcal{Z} \times H_\#^1(0, T) \mid z_{|\zeta,\nu}(T) = z_{|\zeta,\nu}(0) \right\} \quad (61) $$

with $z_{|\zeta,\nu}$ solution of the dynamics

$$ \left\{ \begin{array}{l} \dot{z}_{|\zeta,\nu} = Az_{|\zeta,\nu} + B\nu, \quad t \in [0, T] \\ z_{|\zeta,\nu}(0) = \zeta \end{array} \right. \qquad (62) $$

or equivalently, $z_{|\zeta,\nu} \in C^0([0, T]; \mathcal{Z})$ satisfying

$$ z_{|\zeta,\nu}(t) = e^{tA}\zeta + (\Lambda\nu)(t). $$

In order to practically minimize the functional (60), we can therefore rely on a penalized approach with the minimization of $\mathcal{J}_{\kappa,\epsilon} : \mathcal{Z} \times H_\#^1(0, T) \mapsto \mathbb{R}$ given by

$$ \begin{aligned} \mathcal{J}_{\kappa, \epsilon}(\zeta, \nu) ={}& \epsilon^{-1} \|z_{|\zeta, \nu}(T) - z_{|\zeta, \nu}(0)\|_{\mathcal{Z}}^2 \\ &+ \frac{1}{2} \int_0^T \left[ |y_\delta(t) - Cz_{|\zeta, \nu}(t)|^2 + \kappa(|\dot{\nu}(t)|^2 + |\nu|^2) \right] dt, \end{aligned} \quad (63) $$

over its domain of definition. To simplify further computations we introduce the bounded operator

$$ {}^b\Psi_T : \begin{array}{l} \mathcal{Z} \times L^2(0, T) \rightarrow L^2(0, T) \\ (\zeta, \nu) \mapsto Cz_{|\zeta, \nu}. \end{array} $$
and rewrite the functional $\mathcal{J}_{\kappa, \epsilon}$ as

$$ \begin{aligned} \mathcal{J}_{\kappa, \epsilon}(\zeta, \nu) ={}& \epsilon^{-1} \|(\mathrm{Id} - e^{TA})\zeta + (\Lambda\nu)(T)\|_{\mathcal{Z}}^2 \\ & + \frac{\kappa}{2} \|\nu\|_{\mathrm{H}^1(0,T)}^2 + \frac{1}{2} \left\| y_{\delta} - {}^{\flat}\Psi_T \begin{pmatrix} \zeta \\ \nu \end{pmatrix} \right\|_{\mathrm{L}^2(0,T)}^2 . \end{aligned} \quad (64) $$

The following result is proved in Appendix B.

**Theorem 4.8** For all $\kappa > 0$, there exists a unique minimiser $\bar{\xi} = (\bar{\zeta}, \bar{\nu}) \in \mathcal{Z} \times \mathrm{H}_{\#}^1(0, T)$ of the strongly-convex functional $\mathcal{J}_{\kappa, \epsilon}$.

We now consider a sequence $(\epsilon_n)_{n \in \mathbb{N}}$ strictly decreasing and converging to 0, leading to the sequence of optimisers

$$ \mathcal{Z} \times \mathrm{H}_{\#}^1(0, T) \ni \bar{\xi}_n = (\bar{\zeta}_n, \bar{\nu}_n) = \arg \min_{(\zeta, \nu)} \mathcal{J}_{\kappa, \epsilon_n}(\zeta, \nu). $$

Finally, we provide in Appendix B, the proof of the following result of convergence for the minimizers sequence hence justifying our penalization strategy.

**Theorem 4.9** The sequence $(\bar{\zeta}_n, \bar{\nu}_n)_{n \in \mathbb{N}}$ minimizer of $\mathcal{J}_{\kappa, \epsilon_n}$ converges to the minimizer $(\bar{\zeta}, \bar{\nu})$ of $\mathcal{J}_{\kappa}$ over $\mathcal{K}$ for the norm in $\mathcal{Z} \times \mathrm{H}_{\#}^1(0, T)$, and we have

$$ \lim_{n \in \mathbb{N}} \mathcal{J}_{\kappa, \epsilon_n}(\bar{\zeta}_n, \bar{\nu}_n) = \mathcal{J}_{\kappa}(\bar{\zeta}, \bar{\nu}). $$

### 4.5. An optimal control approach

We are now in position to write the optimality system associated to each penalized criterion $\mathcal{J}_{\kappa, \epsilon}$ and the limit system associated with $\mathcal{J}_{\kappa}$. Let us now introduce for all $z \in \mathrm{L}^2((0, T); \mathcal{Z})$ and $y_{\delta} \in \mathrm{L}^2(0, T)$ the adjoint dynamics [2]

$$ \left\{ \begin{array}{ll} \dot{q}_{\epsilon} + A^{*}q_{\epsilon} = -\gamma C^{*}(y_{\delta} - Cz), & \text{in } (0, T) \\ q_{\epsilon}(T) = -\epsilon^{-1}(z(T) - z(0)) & \end{array} \right. \qquad (65) $$

which is also well posed as it is considered backward in time with the semi-group generator $-A^*$. Namely, we have $q \in C^0[0, T]$, $\mathcal{Z}$ given by the Duhamel formula

$$ q_{\epsilon}(t) = \epsilon^{-1} e^{(T-t)A^{*}} (z(0) - z(T)) + \gamma \int_{t}^{T} e^{(s-t)A^{*}} C^{*}(y_{\delta}(s) - Cz(s)) \mathrm{d}s. $$

The adjoint variable allows to easily compute the Fréchet derivatives with respect to $\zeta$ and $\nu$. We find for a given $(\zeta, \nu) \in \mathcal{Z} \times \mathrm{H}_{\#}^1(0, T)$

$$ \forall \eta \in \mathcal{Z}, \quad \langle \mathrm{D}_{\zeta} \mathcal{J}_{\kappa, \epsilon}(\zeta, \nu), \eta \rangle = \epsilon^{-1}(\zeta - z_{|\zeta, \nu}(T), \eta) + (q_{\epsilon}|_{\zeta, \nu}(0), \eta)_{\mathcal{Z}}, \quad (66) $$

and

$$ \forall \mu \in \mathrm{H}_{\#}^1(0, T), $$

$$ \langle \mathrm{D}_{\nu} \mathcal{J}_{\kappa, \epsilon}(\zeta, \nu), \mu \rangle = \kappa(\nu, \mu)_{\mathrm{H}^1(0,T)} - \int_{0}^{T} (q_{\epsilon}|_{\zeta, \nu}(t), B\mu(t))_{\mathcal{Z}} \mathrm{d}t, \quad (67) $$
where $q_{\epsilon|\zeta,\nu}$ is the adjoint variable associated with $z_{|\zeta,\nu}$. We obtain the Euler equation associated with the minimization

$$ \forall (\eta, \mu) \in \mathcal{Z} \times H_{\sharp}^{1}(0, T), \quad \epsilon^{-1}(\bar{\zeta}_{\epsilon} - \bar{z}_{\epsilon}(T), \eta)_{\mathcal{Z}} \\ - (\bar{q}_{\epsilon}(0), \eta)_{\mathcal{Z}} + \kappa(\bar{\nu}_{\epsilon}, \mu)_{H^{1}(0,T)} - \int_{0}^{T} (\bar{q}_{\epsilon}(t), B\mu(t))_{\mathcal{Z}} dt = 0, \quad (68) $$

where $\bar{q}_{\epsilon}$ is the adjoint variable associated with the optimal trajectory $\bar{z}_{\epsilon} = z_{|\zeta_{\epsilon}, \bar{\nu}_{\epsilon}}$ and the available measurements $y_{\delta}$. We then introduce $\Delta_{\sharp}$ the laplacian in time with periodic boundary condition. From Lax-Millgram theorem $\mathrm{Id} - \Delta_{\sharp}$ is invertible from $L^2(0, T)$ to $H_{\sharp}^{1}(0, T)$. This leads to the so-called two-ends problem defining the optimal dynamics of the estimator

$$ \left\{ \begin{array}{ll} \dot{\bar{z}}_{\epsilon} = A\bar{z}_{\epsilon} + \kappa^{-1}B(\mathrm{Id} - \Delta_{\sharp})^{-1}B^{*}\bar{q}_{\epsilon}, & \text{in } (0, T) \\ \dot{\bar{q}}_{\epsilon} + A^{*}\bar{q}_{\epsilon} = -C^{*}(y_{\delta} - C\bar{z}_{\epsilon}), & \text{in } (0, T) \\ \bar{z}_{\epsilon}(0) = \epsilon\bar{q}_{\epsilon}(0) + \bar{z}_{\epsilon}(T) \\ \bar{q}_{\epsilon}(T) = \bar{q}_{\epsilon}(0) \end{array} \right. \quad (69) $$

We point out that as $(\mathrm{Id} + \Delta_{\sharp})^{-1}$ is a non-local term with respect to time. Finally from (68) and from $\bar{\zeta}_{\epsilon}$ converging to $\bar{\zeta}$ when $\epsilon$ goes to 0 that the optimality system naturally becomes at the limit

$$ \left\{ \begin{array}{ll} \dot{\bar{z}} = A\bar{z} + \kappa^{-1}B(\mathrm{Id} - \Delta_{\sharp})^{-1}B^{*}\bar{q}, & \text{in } (0, T) \\ \dot{\bar{q}} + A^{*}\bar{q} = -C^{*}(y_{\delta} - C\bar{z}), & \text{in } (0, T) \\ \bar{z}(0) = \bar{z}(T) \\ \bar{q}(T) = \bar{q}(0) \end{array} \right. \quad (70) $$

# 5. Numerical methods

Based on the preceding analysis, we propose an adapted numerical reconstruction strategy based on the following principle: Discretization of the forward problem and control-based formulation of the inverse problem redeveloped at the discrete level. As it is well known that the interaction between control and discretization can be complicated [38], the proposed strategy will ensure that our method is numerically well suited.

## 5.1. Model discretization

Concerning the spatial discretization, we rely on a finite element discretization based on the variational formulation associated with problem (9, 10, 13). For any test function in $q = (b, w, Q) \in \mathcal{D}(A)$ we obtain from (9, 10) that the strict solution $z = (a, v, P)$ satisfies at each time $t$,

$$ (\dot{a}, b)_{L^2(0,1)} + (v', b)_{L^2(0,1)} + (\dot{v}, w)_{L^2(0,1)} + (a', w)_{L^2(0,1)} \\ + k(v, w)_{L^2(0,1)} + (\dot{P} + RP - c_r v(1, \cdot)) Q = 0. $$

The second term is integrated by parts and boundary terms are substituted either using the property that $q \in \mathcal{D}(A)$ or the boundary condition $v(0, \cdot) = \nu$.

$$ (\dot{a}, b)_{L^2(0,1)} + (\dot{v}, w)_{L^2(0,1)} + \dot{P}Q $$
$$+(a', w)_{L^2(0,1)} - (v, b')_{L^2(0,1)} + k(v, w)_{L^2(0,1)} \\ + RPQ + k_r v(1, \cdot)w(1) = b(0)\nu.$$

Then the formulation above is discretized using a conforming approach. A finite-dimensional space $V_h \subset D(A)$ is constructed using a space-grid $0 \le x_i = ih \le 1$ with $h = N_x^{-1}$, $0 \le i \le N_x$ and Lagrangian finite element of order $k = 1$,

$$V_h = \left\{ (a_h, v_h, P_h) \in H^1(0,1) \times H_l^1(0,1) \times \mathbb{R}, \right. \\ \left. \qquad a_h(1) = k_r v_h(1) + c_r P_h, \ a_h|_{[x_i, x_{i+1}]} \in \mathbb{P}_k, \ v_h|_{[x_i, x_{i+1}]} \in \mathbb{P}_k \right\}. \quad (71)$$

The space $V_h$ is equipped with the scalar product $(\cdot, \cdot)_h$ that is an approximation of the scalar product in $\mathcal{Z}$ using the trapezoidal quadrature formula. We introduce the operator $A_h \in \mathcal{L}(V_h)$ such that for all functions $z_h = (a_h, v_h, P_h) \in V_h$ and $q_h = (b_h, w_h, Q_h) \in V_h$

$$(A_h z_h, q_h)_h = (v_h, b'_h)_{L^2(0,1)} - (a'_h, w_h)_{L^2(0,1)} \\ -k(v_h, w_h)_{L^2(0,1)} - k_r v_h(1)w_h(1) - RPQ, \quad (72)$$

as well as the operator $B_h \in \mathcal{L}(\mathbb{R}, V_h)$ defined by

$$(B_h \nu, q_h)_h = b_h(0) \nu.$$

The semi-discrete problem is then equivalent to

$$\dot{z}_h = A_h z + B_h \nu. \quad (73)$$

For the time-discretization, we use a regular time-grid $0 \le t_n = n\tau \le T$ with $\tau = N_T^{-1}T$, $0 \le n \le N_T$. We choose a backward-Euler approach for three reasons. First, the scheme is unconditionally stable. Second, the additional numerical dissipation has a strong effects on spurious high frequency components of the solutions that otherwise could have compromised the reconstruction (see for instance the review [38]).

After time discretization, (73) is replaced by

$$\frac{z_h^{n+1} - z_h^n}{\tau} = A_h z^{n+1} + B_h \nu^{n+1}. \quad (74)$$

that can be rewritten

$$z_h^{n+1} = \Phi_{h,\tau} z^n + B_{h,\tau} \nu^{n+1}$$

with $\Phi_{h,\tau} = (\mathrm{Id} - \tau A_h)^{-1}$ and $B_{h,\tau} = \tau(\mathrm{Id} - \tau A_h)^{-1}B_h$. Since $\{\nu^n\}$ is an approximation of a periodic function we impose

$$\nu^{N_T} = \nu^0.$$

## 5.2. Optimization based on generalized inverse computation

We first emphasize that we assume that the measurements are compatible with the chosen time discretization, i.e., we have $(y_\delta^n)_{0 \le n \le N}$. Note that if the time step of the model is smaller than the time samples of the measurements, it is always possible to re-sample by interpolation up to additional measurement errors with controlled effect
[9]. Our goal is to minimize the discrete counterpart of the penalized criterion $\mathcal{J}_{\kappa, \epsilon}$ introduced in (63). We define the functional $\mathcal{J}_{\kappa, \epsilon}^{N_T}$ in the finite dimensional space $V_h \times \mathbb{R}^{N_T}$ as.

$$ \mathcal{J}_{\kappa, \epsilon}^{N_T}(\zeta_h, \nu_h) = \epsilon^{-1} \left\| z_{h|\zeta,\nu}^{N_T} - z_{h|\zeta,\nu}^{0} \right\|_{\mathcal{Z}}^2 + \frac{\tau}{2} \sum_{n=0}^{N_T-1} \left[ |y_{\delta}^{n} - C z_{h|\zeta,\nu}^{n}|^2 + \kappa \left( |\nu^n|^2 + \left| \frac{\nu^{n+1} - \nu^n}{\tau} \right|^2 \right) \right], (75) $$

subject to

$$ \left\{ \begin{array}{l} z_{h|\zeta,\nu}^{n+1} = \Phi_{h,\tau} z_{h|\zeta,\nu}^{n} + B_{h,\tau} \nu^{n+1}, \quad 0 \le n \le N_T - 1 \\ z_{h|\zeta,\nu}^{0} = \zeta_h. \end{array} \right. \qquad (76) $$

The minimization of such function leads to a discretization of the optimality system (69), see Appendix C for detail computation.

Minimization of $\mathcal{J}_{\kappa, \epsilon}^{N_T}(\zeta_h, \nu_h)$ can be performed by inverting the underlying linear system associated with the Euler equation of this quadratic problem. This task is made possible here because our initial problem is one-dimensional in space and therefore has a limited number of degrees of freedom after space and time discretization. Developing the solution of (73) using a discrete Duhamel formula

$$ z_{h|\zeta,\nu}^{n} = \Phi_{h,\tau}^{n} \zeta_{h} + \sum_{k=1}^{n} \Phi_{h,\tau}^{n-k} B_{h,\tau} \nu^{k}, $$

the minimization is unconstrained by making the initial condition dependency explicit, namely

$$ \mathcal{J}_{\kappa, \epsilon}^{N_T}(\zeta_h, (\nu_n)_{n \in [1, N_T]}) = \epsilon^{-1} \left\| (\Phi_{h,\tau}^{N_T} - \text{Id})\zeta_h + \sum_{n=1}^{N_T} \Phi_{h,\tau}^{N_T-n} B_{h,\tau} \nu^n \right\|_{\mathcal{Z}}^2 \\ + \frac{\tau}{2} \sum_{n=0}^{N_T-1} \left[ \begin{aligned} & |y_{\delta}^{n} - C \Phi_{h,\tau}^{n} \zeta_h + \sum_{k=1}^{n} C \Phi_{h,\tau}^{n-k} B_{h,\tau} \nu^k|^2 \\ & + \frac{\kappa}{2} \left( |\nu^n|^2 + \left| \frac{\nu^{n+1} - \nu^n}{\tau} \right|^2 \right) \end{aligned} \right]. \quad (77) $$

In (75) and (77), we have that $\nu^0 = \nu^{N_T}$, such that $\nu^{N_T}$ becomes a non independent variable. Introducing

$$ X = \begin{pmatrix} \zeta_h \\ \nu^0 \\ \vdots \\ \nu^{N_T-1} \end{pmatrix} \in V_h \times \mathbb{R}^{N_T} \quad \text{and} \quad Y = \begin{pmatrix} y_{\delta}^0 \\ \vdots \\ y_{\delta}^{N_T-1} \end{pmatrix} $$

we rewrite

$$ \mathcal{J}_{\kappa, \epsilon}^{N_T}(X) = \frac{\epsilon^{-1}}{2} \|\Psi_{\sharp} X\|_{\mathcal{Z}}^2 + \frac{1}{2} \|Y - \Psi_C X - \Psi_B X\|_{\mathbb{R}^{N_T}}^2 \\ + \frac{\kappa}{2} \left( X, \begin{pmatrix} 0 & 0 \\ 0 & K_{\sharp} \end{pmatrix} X \right)_{\mathbb{R}^{N_T}}, $$
where

$$ \Psi_{\#} = ( \Phi_{h,\tau}^{N_T} - \mathrm{Id} \quad \Phi_{h,\tau}^{N_T-1} B_{h,\tau} \quad \dots \quad B_{h,\tau} ) \in \mathcal{L}(\mathcal{V}_h \times \mathbb{R}^{N_T}), $$

while

$$ \Psi_C = \begin{pmatrix} C_h \\ \vdots \\ C_h \Phi_{h,\tau}^{N_T-1} \end{pmatrix} \in \mathcal{L}(\mathcal{V}_h), $$

and

$$ \Psi_B = \begin{pmatrix} C_h B_{h,\tau} & 0 \\ \vdots & \vdots \\ C_h \Phi_{h,\tau}^{N_T-1} B_{h,\tau} & 0 \end{pmatrix} \in \mathcal{L}(\mathbb{R}^{N_T}). $$

Therefore, minimizing $\mathcal{J}_{\kappa, \epsilon}^{N_T}$ can be performed by solving

$$ \left[ \epsilon^{-1} \Psi_{\#}^* \Psi_{\#} + \tau \begin{pmatrix} \Psi_C^* \Psi_C & \Psi_C^* \Psi_B \\ \Psi_B^* \Psi_C & \Psi_B^* \Psi_B + \kappa K_{\#} \end{pmatrix} \right] X = \tau \begin{pmatrix} \Psi_C & \Psi_B \end{pmatrix}^* Y, \quad (78) $$

to the prize of being able to store the necessary dense matrices.

It should be noted that this system is solved with a direct solver in the cases presented below. This is possible due to the limited number of degrees of freedom in the time and space discretization, which is here of the order $10^3$. For problems with a larger number of degrees of freedom, a gradient descent approach should be used instead, see Appendix C.

## 5.3. Numerical results

To illustrate our method, we propose to apply the resulting algorithm to the test case shown in Figure 2, which corresponds to the third benchmark case presented in [4]. It consists of a one-dimensional cylindrical vessel with a homogeneous circular cross-section representing the upper thoracic aorta, with a given blood flow as inlet and a three-element Windkessel model as outlet boundary condition. The model parameters were taken from [37] and are reproduced in Table 1. Our objective is to estimate the inlet blood flow using some distal pressure measurements ($P_{obs}$). We present in Section 5.4 the results obtained with synthetic observations for different amounts of artificial noise added to the measurements. Finally in Section 5.5, we exploit this configuration to invert real in vivo recording of aortic pressure.

## 5.4. Synthetic results

We perform the forward problem giving an inlet blood flow as input and obtain the state along the vessel. Then, we use the output distal pressure as observation for the inverse problem in order to estimate the inlet flow. The "physical" parameters of the model are the same as in the direct simulation.
<figure><img src="image_9.png" /><figcaption>**Figure 2.** Configuration of the model including a cylindrical homogeneous vessel representing the upper thoracic aorta and a three-elements Windkessel that takes into account the resistive and compliant behavior of the distal vessels. The inlet blood flow is the input for the forward model.</figcaption></figure>

<table><caption>**Table 1.** Model parameters of the upper thoracic aorta (Readapted from [37]).</caption><thead><tr><td>Property</td><td>Value</td></tr></thead><tbody><tr><td>Length, L</td><td>24.137 cm</td></tr><tr><td>Radius at diastolic pressure, R<sub>0</sub></td><td>1.2 cm</td></tr><tr><td>Initial flow velocity, u(s, 0)</td><td>0 m s<sup>-1</sup></td></tr><tr><td>Initial pressure, P(s, 0)</td><td>0 Pa</td></tr><tr><td>Wall thickness, h<sub>0</sub></td><td>1.2 mm</td></tr><tr><td>Blood density, ρ</td><td>1060 Kg m<sup>-3</sup></td></tr><tr><td>Friction parameter, K<sub>r</sub></td><td>276.46 mPa s</td></tr><tr><td>Velocity profile, α<sub>vp</sub></td><td>1</td></tr><tr><td>Young's modulus, E</td><td>400.0 kPa</td></tr><tr><td>Diastolic pressure, P<sub>D</sub></td><td>9.46 kPa</td></tr><tr><td>External pressure, P<sub>ext</sub></td><td>0 Pa</td></tr><tr><td>Windkessel resistance, R<sub>c</sub></td><td>1.1752 &middot; 10<sup>7</sup> Pa s m<sup>-3</sup></td></tr><tr><td>Windkessel compliance, C<sub>c</sub></td><td>1.0163 &middot; 10<sup>-8</sup> m<sup>3</sup> Pa<sup>-1</sup></td></tr><tr><td>Windkessel resistance, R<sub>per</sub></td><td>1.1167 &middot; 10<sup>8</sup> Pa s m<sup>-3</sup></td></tr></tbody></table>

*Noise and regularization* In order to observe the behavior of our reconstruction strategy in presence of noise, very likely when real clinical data are involved, we propose to corrupt the synthetic pressure data $(y_n)_{1 \le n \le N_T}$ of sampling rate $T/N_T$ with some random
artificial noise. We define the discretized noise data as

$$y_n^\delta = y_n + \frac{\delta}{\sqrt{T}} \|(y_n)_{1 \le n \le N_T}\|_{L^2(0,T)} \chi_n, \quad 1 \le n \le N_T$$

where $\chi_n$ are independent and identically distributed gaussian variables in $N(0, 1)$ and with the slight abuse of notation the $L^2(0, T)$ norm for such discrete signal is defined by

$$\|(y_n)_{1 \le n \le N_T}\|_{L^2(0,T)}^2 = \frac{T}{N_T} \sum_{n=1}^{N_T} y_n^2$$

One can note that with such choice, we have

$$\frac{\|(y_n^\delta - y_n)_{1 \le n \le N_T}\|_{L^2(0,T)}^2}{\|(y_n)_{1 \le n \le N_T}\|_{L^2}^2} = \frac{\delta^2}{N_T} \sum_{i=1}^{N_T} \chi_n^2,$$

which, by property of the empirical variance, gives

$$\lim_{N_T \to \infty} \frac{\|(y_n^\delta - y_n)_{1 \le n \le N_T}\|_{L^2(0,T)}^2}{\|(y_n)_{1 \le n \le N_T}\|_{L^2(0,T)}^2} = \delta^2.$$

As a consequence, our choice of noise corruption is a good numerical approximation of the deterministic measurement error required in Theorem 4.7. We test the technique with different levels of noise, in particular we choose $\delta = \{0, 0.005, 0.01, 0.02, 0.03\}$ corresponding to realistic cases. For each different $\delta$, we need to adapt the values of $\kappa$ to adjust the weight of the Tikhonov regularization. Using Theorem 4.7, we obtain the optimal theoretical values for $\kappa$, as reported in Table 2. However, another way to find the right balance between regularization and agreement with the measurements is through the L-curve criterion [16], where we plot the norm of the solution – $\|\nu\|_{L^2(0,T)}^2$ – versus the norm of the discrepancy between observations and computed state – $\|y_\delta - Cz_{|\zeta,\nu}|_{L^2(0,T)}^2$ – for different values of $\kappa$ on a log-log plot. The resulting 'L-shape' curve gives the optimal value for the regularization parameter, located in the corner of the 'L' [16]. In practice, we sample $\kappa$ regularly in log-scale between $10^{-8}$ and 1 and we find the ideal regularization parameters, as reported in Table 2, to produce the L-curve. Note that both methods were meaningless when dealing with $\delta = 0$, in this particular case the optimal value of $\kappa$ was manually optimized. In Figure 3 we show the results obtained for the different levels of noise, comparing the input estimations obtained choosing the optimal value of $\kappa$ by using Theorem 4.7 with the ones obtained relying on the L-criterion [16]. In the fourth column, we picture the L-curves obtained for each value of $\delta$. One can observe that there are almost no difference in the reconstructed output pressure between the two choices of $\kappa$. Note that the estimated input flow obtained with the theoretical value of $\kappa$ shows, as expected, a lower dampening of the oscillations caused by noise while the L-curve based $\kappa$, being slightly higher for every value of $\delta$, leads to a smoother curve. Nonetheless both choices seem adapted and work quite well to obtain a coherent input flow.
<figure><img src="image_10.png" /><figcaption>**Figure 3.** The estimated input flow obtained with the theorem-based (see Theorem 4.7) choice (black curve in the second column) and with the L-curve-based choice (red curve in the third column) of the regularization parameter $\kappa$ are compared with the expected flow (blue curve in the second and third column), being the input used for the forward problem. This comparison is done for an increasing level of noise, from the top to the button, and the pressure observation (grey line in the first column) is plotted together with the computed pressure with both choices of $\kappa$ (black and red curve for the theorem-based and L-curve based $\kappa$, respectively). The L-curves obtained for each value of $\delta$ are shown in the right column.</figcaption></figure>
<table><caption>**Table 2.** Optimal values for the regularization parameter $\kappa$ obtained applying Theorem 4.7, in the second column, and using the approach involving the analysis of the L-curve, in the third column.</caption><thead><tr><th>δ</th><th>Theorem based κ</th><th>L-curve based κ</th></tr></thead><tbody><tr><td>0.005</td><td>9.55 10<sup>-5</sup></td><td>10<sup>-3</sup></td></tr><tr><td>0.01</td><td>3.82 10<sup>-4</sup></td><td>5 10<sup>-3</sup></td></tr><tr><td>0.02</td><td>1.5 10<sup>-3</sup></td><td>10<sup>-2</sup></td></tr><tr><td>0.03</td><td>3.4 10<sup>-3</sup></td><td>10<sup>-2</sup></td></tr></tbody></table>

**Sampling mismatch** When it comes to real-world applications, measurements have their own sampling time and this may be far from the $\tau$ used for numerical simulations, imposing to either exploit the data only when they are available or interpolating the data to be compatible with the simulation time-step. In this work, as we are dealing with periodic signals, we choose to rely on a Fourier interpolation. We test two sampled version of the original synthetic observations – $\tau = 8ms$ and $\tau = 5ms$ – being the sampling time often found in catheter-measured arterial pressure. On these sample data, we apply two levels of noise following the approach explained in the previous paragraph, $\delta = \{0.005, 0.02\}$, considered representative of the pressure data usually registered. Finally, we interpolate these data to have discretisation time-step $\tau = 0.001ms$. The results shown in Figure 4 are obtained by setting the L-curve based $\kappa$ values reported in Table 2 for the corresponding $\delta$. The results demonstrate the validity of our model-based inversion approach for realistic catheter measurements.

### 5.5. Real data results

Finally, we decided to evaluate our strategy with real clinical data, although the proposed method is only a first step in which we linearized the original blood flow model. The data, an *in vivo* measurement of aortic pressure, were recorded and anonymized by F. Vallée (M.D. and D.r), anesthetist at Lariboisière Hospital in Paris (APHP). These data are measured with an intra-aortic catheter equipped with a pressure transducer capable of recording blood pressure and velocity over time. The pressure measured directly over the renal artery is used as an observation and the blood flow coming from the left ventricle is reconstructed. The configuration of the IP procedure is shown in Figure 5, and the data injected into the model and the estimated source are shown in the boxes on the right and left sides of the vessel, respectively. Unfortunately, when the catheter is close to the aortic valve, the data are very noisy and, especially for velocity measurements, difficult to analyze from a physiological point of view. For this reason, the results shown in Figure 6 show the comparison between the calculated blood velocity and that recorded by the catheter near the pressure observation site rather than at the inlet. From Figure 6, we can also see that the computed pressure is in agreement with the data and that the computed blood velocity displays the same trend of the *in vivo* curve. The latter exhibits a high degree of noise, making it difficult to perform a deeper analysis of the discrepancy between the calculated condition and the measurements.
<figure><img src="image_11.png" /><figcaption>**Figure 4.** Synthetic pressure data (grey line in the left column) sampled with $\tau = 5$ ms (first and third rows) and $\tau = 8$ ms (second and fourth rows) and with an additional noise of magnitude $\delta = 0.005$ (first and second rows) and $\delta = 0.2$ (third and fourth rows) are interpolated using a MATLAB function based on Fourier interpolation. These data used as observations are compared with the computed output pressure (red line in the left column) and the estimated input flow (red line in the right column) is plotted against the expected estimates (blue line in the right column).</figcaption></figure>
We see that the positive peak of the velocity as well as the trend of the curve are well reproduced. Some differences are observable in the baseline value and in the computed negative peak, which are not present in the data.

<figure><img src="image_12.png" /><figcaption>**Figure 5.** Qualitative representation of an aortic vessel with its main branching vessel, from the heart up to the femoral arteries bifurcation. The location of the measurement site is pointed out, on the right, and the distal pressure data used as observation are shown in the box. The estimation spot is also indicated, on the left side of the vessel, and the source estimated by using our method is reported.</figcaption></figure>

# 6. Conclusions and Perspectives

In this paper, we present a strategy for an inverse problem based on optimal control theory to recover the inlet flow and initial conditions for a reduced-order model of blood circulation based on a linearized formulation. The forward and inverse problems are analyzed from the mathematical point of view, with a particular emphasis on obtaining a fundamental observability condition. This work sets the stage for a new approach to the clinical problem of deriving inaccessible information about cardiac output from accessible peripheral pressure measurements. In practice, we describe an inversion strategy based on a variational approach in which we penalize the periodicity required for the observability condition. Our tests in the presence of noise, accounting for sampling error, and using real *in vivo* data are promising but should be extended to
<figure><img src="image_13.png" /><figcaption>**Figure 6.** Comparison between the blood velocity computed (in red) and expected (in blue) at a distal location close to the observation site, on the left, and between the pressure computed (in red) and used as input observation (in grey), on the right.</figcaption></figure>

evaluate the potential of the estimate for clinical applications. In addition, preliminary parameter estimation is essential to enable patient-specific application. In the clinical context, we could envision such parameter calibration being performed initially using a parameter identification technique such as [6] before the catheter is pulled down during the cardiac procedure. A natural evolution of this work is to propose a strategy adapted to the nonlinear model, possibly by iterating over a linearized configuration like the one studied in this work. A second perspective would be to consider a more complex configuration involving upper limb vessels. This would pave the way of using radial pressure measurements as noninvasive observations to reconstruct cardiac output.

## Acknowledgement

We are very grateful to Fabrice Vallée for providing us with the data that support the findings.

## Data availability statement

The data that support the findings of this study are available upon reasonable request from the authors.

## References

[1] Christopher J Arthurs, Nan Xiao, Philippe Moireau, Tobias Schaeffter, and C Alberto Figueroa. A flexible framework for sequential estimation of model parameters in computational hemodynamics. *Advanced modeling and simulation in engineering sciences*, 7(1):1–37, 2020.

[2] A. Bensoussan. *Filtrage optimal des systèmes linéaires*. Dunod, 1971.
[3] A. Bensoussan, M. C. Delfour, G. Da Prato, and S. K. Mitter. *Representation and Control of Infinite Dimensional Systems*. Birkhauser Verlag, Boston, second edition edition, 2007.

[4] Etienne Boileau, Perumal Nithiarasu, Pablo J. Blanco, Lucas O. Müller, Fredrik Eikeland Fossan, Leif Rune Hellevik, Wouter P. Donders, Wouter Huberts, Marie Willemet, and Jordi Alastruey. A benchmark study of numerical schemes for one-dimensional arterial blood flow modelling. *International Journal for Numerical Methods in Biomedical Engineering*, 31(10):e02732, July 2015.

[5] E. Bollache, N. Kachenoura, A. Redheuil, F. Frouin, E. Mousseaux, P. Recho, and D. Lucor. Descending aorta subject-specific one-dimensional model validated against in vivo data. *Journal of Biomechanics*, 47(2):424-431, 2014.

[6] A. Caiazzo, Federica Caforio, Gino Montecinos, Lucas O. Muller, Pablo J. Blanco, and Eluterio F. Toro. Assessment of reduced-order unscented kalman filter for parameter identification in 1-dimensional blood flow models using experimental data. *International Journal for Numerical Methods in Biomedical Engineering*, 33(8):e2843, jan 2017.

[7] Chen-Huan Chen, Erez Nevo, Barry Fetics, Peter H. Pak, Frank C.P. Yin, W. Lowell Maughan, and David A. Kass. Estimation of central aortic pressure waveform by mathematical transformation of radial tonometry pressure. *Circulation*, 95(7):1827-1836, April 1997.

[8] Patrick Ciarlet. T-coercivity: Application to the discretization of helmholtz-like problems. *Computers & Mathematics with Applications*, 64(1):22-34, 2012.

[9] N. Cîndea, A. Imperiale, and P. Moireau. Data assimilation of time under-sampled measurements using observers, the wave-like equation example. *ESAIM. Control, Optimisation and Calculus of Variations*. *European Series in Applied and Industrial Mathematics*, 21(3):635-669, 2015.

[10] Marta D'Elia, Mauro Perego, and Alessandro Veneziani. A variational data assimilation procedure for the incompressible navier-stokes equations in hemodynamics. *Journal of Scientific Computing*, 52(2):340-359, 2012.

[11] H. W. Engl, M. Hanke, and A. Neubauer. *Regularization of inverse problems*, volume 375 of *Mathematics and its Applications*. Kluwer Academic Publishers Group, Dordrecht, 1996.

[12] Sally Epstein, Marie Willemet, Phil J. Chowienczyk, and Jordi Alastruey. Reducing the number of parameters in 1d arterial blood flow modeling: less is more for patient-specific simulations. *American Journal of Physiology-Heart and Circulatory Physiology*, 309(1):H222-H234, jul 2015.

[13] B. Fetics, E. Nevo, Chen-Huan Chen, and D.A. Kass. Parametric model derivation of transfer function for noninvasive estimation of aortic pressure by radial tonometry. *IEEE Transactions on Biomedical Engineering*, 46(6):698-706, June 1999.

[14] Luca Formaggia, Jean-Frédéric Gerbeau, Fabio Nobile, and Alfio Quarteroni. On the Coupling of 3D and 1D Navier-Stokes Equations for Flow Problems in Compliant Vessels. Research Report RR-3862, INRIA, 2000. Projet M3N.

[15] Fredrik E Fossan, Jorge Mariscal-Harana, Jordi Alastruey, and Leif R Hellevik. Optimization of topological complexity for one-dimensional arterial blood flow models. *Journal of the Royal Society Interface*, 15(149):20180546, 2018.

[16] Per Christian Hansen. *The L-Curve and Its Use in the Numerical Treatment of Inverse Problems*, volume 4, pages 119-142. 01 2001.

[17] Jona Joachim, Fabrice Vallée, Arthur Le Gall, Joaquim Matéo, Stéphanie Lenck, Sandrine Millasseau, Emmanuel Houdart, Alexandre Mebazaa, and Etienne Gayat. Velocity-pressure loops for continuous assessment of ventricular afterload: influence of pressure measurement site. *Journal of Clinical Monitoring and Computing*, 32(5):833-840, November 2017.

[18] Rajnesh Lal, Bijan Mohammadi, and Franck Nicoud. Data assimilation for identification of cardiovascular network characteristics. *International journal for numerical methods in biomedical engineering*, 33(5):e2824, 2017.

[19] L. Ljung. *System Identification: Theory for the User*. Prentice Hall information and system sciences series. Prentice Hall PTR, 1999.

[20] Damiano Lombardi. Inverse problems in 1d hemodynamics on systemic networks: A sequential
approach. *International journal for numerical methods in biomedical engineering*, 30(2):160-179, 2014.

[21] D. Luenberger and Yinyu Ye. *Linear and nonlinear programming*. Springer, 2011.

[22] Jessica Manganotti, Federica Caforio, François Kimmig, Philippe Moireau, and Sebastien Imperiale. Coupling reduced-order blood flow and cardiac models through energy-consistent strategies: modeling and discretization. *Advanced Modeling and Simulation in Engineering Sciences*, 8(1), September 2021.

[23] Vincent Martin, Francois Clément, Astrid Decoene, and Jean-Frédéric Gerbeau. Parameter identification for a one-dimensional blood flow model. In *ESAIM: Proceedings*, volume 14, pages 174-200. EDP Sciences, 2005.

[24] Koen S. Matthys, Jordi Alastruey, Joaquim Peiró, Ashraf W. Khir, Patrick Segers, Pascal R. Verdonck, Kim H. Parker, and Spencer J. Sherwin. Pulse wave propagation in a model human arterial network: Assessment of 1-d numerical simulations against in vitro measurements. *Journal of Biomechanics*, 40(15):3476-3486, jan 2007.

[25] Sandrine C. Millasseau, Sundip J. Patel, Simon R. Redwood, James M. Ritter, and Philip J. Chowienczyk. Pressure wave reflection assessed from the peripheral pulse. *Hypertension*, 41(5):1016-1020, May 2003.

[26] Philippe Moireau, Cristobal Bertoglio, Nan Xiao, C Alberto Figueroa, CA Taylor, Dominique Chapelle, and J-F Gerbeau. Sequential identification of boundary support parameters in a fluid-structure vascular model using patient image data. *Biomechanics and modeling in mechanobiology*, 12(3):475-496, 2013.

[27] Lucas O Müller, Alfonso Caiazzo, and Pablo Javier Blanco. Reduced-order unscented kalman filter with observations in the frequency domain: application to computational hemodynamics. *IEEE Transactions on Biomedical Engineering*, 66(5):1269-1276, 2018.

[28] Tikhonov A. N. Solution of incorrectly formulated problems and the regularization method. *Soviet Math.*, 4:1035-1038, 1963.

[29] Alfredo L. Pauca, Michael F. O'Rourke, and Neal D. Kon. Prospective evaluation of a method for estimating ascending aortic pressure from the radial artery pressure waveform. *Hypertension*, 38(4):932-937, October 2001.

[30] J Peiro, SJ Sherwin, KH Parker, V Franke, Luca Formaggia, Daniele Lamponi, and Alfio Quarteroni. Numerical simulation of arterial pulse propagation using one-dimensional models. *Wall-fluid interactions in physiological flows*, pages 1-36, 2003.

[31] Alfio Quarteroni and Luca Formaggia. Mathematical modelling and numerical simulation of the cardiovascular system. In *Computational Models for the Human Body*, volume 12 of *Handbook of Numerical Analysis*, pages 3-127. Elsevier, 2004.

[32] Philippe Reymond, Fabrice Merenda, Fabienne Perren, Daniel Rüfenacht, and Nikos Stergiopulos. Validation of a one-dimensional model of the systemic arterial tree. *American Journal of Physiology-Heart and Circulatory Physiology*, 297(1):H208-H222, jul 2009.

[33] B Sixou, L. Boissel, and M. Sigovan. Vascular blood flow reconstruction from tomographic projections with the adjoint method and receding optimal control strategy. *Journal of Physics: Conference Series*, 904:012001, oct 2017.

[34] Wim J. Stok, Berend E. Westerhof, Ilja Guelen, and John M. Karemaker. Aortic pressure wave reconstruction during exercise is improved by adaptive filtering: a pilot study. *Medical & Biological Engineering & Computing*, 49(8):909-916, 2011.

[35] M. Tucsnak and G. Weiss. *Observation and control for operator semigroups*. Birkhäuser Advanced Texts: Basler Lehrbücher. Birkhäuser Verlag, Basel, 2009.

[36] Marie Willemet and Jordi Alastruey. Arterial pressure and flow wave analysis using time-domain 1-d hemodynamics. *Annals of biomedical engineering*, 43(1):190-206, 2015.

[37] Nan Xiao, Jordi Alastruey, and C. Alberto Figueroa. A systematic comparison between 1-d and 3-d hemodynamics in compliant arterial models. *International Journal for Numerical Methods in Biomedical Engineering*, 30(2):204-231, sep 2013.
[38] Enrique Zuazua. Propagation, Observation, and Control of Waves Approximated by Finite Difference Methods. SIAM Review, 47(2):197–243, 2005.

**Appendix A. The operator $(A, \mathcal{D}(A))$ is maximal dissipative.**

Let us recall that

$$A = \begin{pmatrix} 0 & -\partial_s & 0 \\ -\partial_s & -k & 0 \\ 0 & c_r \gamma_r & -R \end{pmatrix}$$

with $\gamma_r$ the trace operator in $s = 1$ and

$$\mathcal{D}(A) := \{ (a, v, P) \in \mathbb{H}^1(0, 1) \times \mathbb{H}_\ell^1(0, 1) \times \mathbb{R} \mid a(1) = k_r v(1) + c_r P \},$$

where $\mathbb{H}_\ell^1(0, 1) = \{w \in \mathbb{H}^1(0, 1) \mid w(0) = 0\}$.

We first prove that $A$ is dissipative. Indeed, for all $z = (a, v, P) \in \mathcal{D}(A)$, we have

$$\begin{aligned} (z, Az)_\mathcal{Z} &= \int_0^1 \left[ -a'(s)v(s) - v'(s)a(s) - kv(s)^2 \right] ds + c_r v(1)P - RP^2 \\ &= -a(1)v(1) - \int_0^1 kv(s)^2 ds + c_r v(1)P - RP^2 \\ &= -k_r v(1)^2 - k \|v\|_{\mathbb{L}^2(0,1)}^2 - RP^2 \le 0. \end{aligned}$$

Then, we need to prove that $A$ is maximal. To do so we consider $\lambda = 0$ and we first show that the bilinear form $\mathcal{D}(A) \times \mathcal{D}(A) \ni (z, q) \mapsto -(Az, q)$ is T-coercive on $\mathcal{D}(A)$ equipped with the norm $\|z\|_{\mathcal{D}(A)}^2 = \|a'\|_{\mathbb{L}^2(0,1)} + \|v'\|_{\mathbb{L}^2(0,1)} + P^2$ for $z = (a, v, P) \in \mathcal{D}(A)$ equivalent to the domain graph norm. Note in particular that if $\|a'\|_{\mathbb{L}^2(0,1)} + \|v'\|_{\mathbb{L}^2(0,1)} + P^2 = 0$ for $(a, v, P) \in \mathcal{D}(A)$, then $P = 0$, $v = 0$, hence $a' = 0$ and $a(1) = k_r v(1) + c_r P = 0$ implies that $a = 0$. For a given $(\alpha, \beta)$ two strictly positive scalar parameters, we define the transformation

$$T_\beta^\alpha : \mathcal{D}(A) \ni \begin{pmatrix} a \\ v \\ P \end{pmatrix} \mapsto \begin{pmatrix} a + v' \\ \alpha v + \beta a' \\ P \end{pmatrix}$$

which is clearly bounded in $\mathcal{D}(A)$. This transformation is bijective from $\mathcal{D}(A)$ to $\mathcal{D}(A)$. Indeed, let us consider $(g, h, I) \in \mathcal{D}(A)$ such that

$$\left\{ \begin{array}{ll} a + v' = g, & \text{in } (0, 1), \\ \alpha v + \beta a' = h, & \text{in } (0, 1), \\ P = I. & \end{array} \right.$$

where we seek a solution $(a, v, P) \in \mathcal{D}(A)$. We must find a solution in $\mathbb{H}^1(0, 1)$ of

$$\left\{ \begin{array}{ll} \alpha v - \beta v'' = \beta g' - h, & \text{in } (0, 1), \\ v(0) = 0 \\ v'(1) + k_r v(1) = g(1) - c_R I & \end{array} \right.$$
which is given by a direct application of Lax-Milgram Theorem. Then, we reconstruct $a(x) = \int_{0}^{x} \beta^{-1}(h(s) - \alpha v(s)) ds$. The solution $(a, v, P) \in \mathcal{D}(A)$ is unique by unicity in Lax-Milgram Theorem, hence $T_{\beta}^{\alpha}$ is bijective from $\mathcal{D}(A)$ to $\mathcal{D}(A)$.

Then, we compute for all $z = (a, v, P) \in \mathcal{D}(A)$

$$ -(Az, T_{\beta}^{\alpha}z)_{\mathcal{Z}} = \int_{0}^{1} \left[ (a(s) + v'(s))v'(s) + (\alpha v(s) + \beta a'(s))a'(s) \right] ds \\ + \int_{0}^{1} k(\alpha v(s) + \beta a'(s))v(s) ds - c_{r}v(1)P + RP^{2}, $$

By integrating by part the first term and using $a(1) = k_{r}v(1) + c_{r}P = 0$, we get

$$ -(Az, T_{\beta}^{\alpha}z)_{\mathcal{Z}} = \int_{0}^{1} \left[ -a'(s)v(s) + v'(s)v'(s) + (\alpha v(s) + \beta a'(s))a'(s) \right] ds \\ + \int_{0}^{1} k(\alpha v(s) + \beta a'(s))v(s) ds + k_{r}v(1)^{2} + RP^{2}, $$

Choosing $\alpha > 0, \beta > 0$ such that $\alpha + k\beta = 1$, we obtain

$$ -(Az, T_{\beta}^{\alpha}z)_{\mathcal{Z}} = \int_{0}^{1} \left[ v'(s)v'(s) + \beta a'(s)a'(s) + k\alpha v(s)v(s) \right] ds \\ + k_{r}v(1)^{2} + RP^{2} \\ \geq \inf(\alpha, \beta) \|z\|_{\mathcal{D}(A)}^{2}. $$

Therefore from [8, Theorem 1], for all $f \in \mathcal{Z} \subset \mathcal{D}(A)'$, there exists $z \in \mathcal{Z}$ such that $-Az = f$ which implies that $\lambda \mathrm{Id} - A$ is bijective from $\mathcal{D}(A)$ to $\mathcal{Z}$ for $\lambda = 0$, namely $A$ is maximal.

## Appendix B. Justifying the penalization strategy proposed in Section 4.4

**Proof of Theorem 4.8.** The functional $\mathcal{J}_{\kappa, \epsilon}$ is continuous and quadratic. Further, $\mathcal{J}_{\kappa, \epsilon}$ is differentiable in the sense of Fréchet with

$$ \mathcal{J}_{\kappa, \epsilon}(\zeta_2, \nu_2) - \mathcal{J}_{\kappa}(\zeta_1, \nu_1) = \left\langle D \mathcal{J}_{\kappa}(\zeta_1, \nu_1), \begin{pmatrix} \zeta_2 - \zeta_1 \\ \nu_2 - \nu_1 \end{pmatrix} \right\rangle \\ + {}^{\flat}\sigma((\zeta_2 - \zeta_1, \nu_2 - \nu_1), (\zeta_2 - \zeta_1, \nu_2 - \nu_1)) $$

where $D\mathcal{J}_{\kappa, \epsilon} \in \mathcal{L}(\mathcal{Z} \times H_{\sharp}^{1}(0, T), (\mathcal{Z} \times H_{\sharp}^{1}(0, T))')$ is the Fréchet derivative and ${}^{\flat}\sigma$ is a bilinear form continuous in $(\mathcal{Z} \times H_{\sharp}^{1}(0, T)) \times (\mathcal{Z} \times H_{\sharp}^{1}(0, T))$ given by

$$ {}^{\flat}\sigma((\zeta, \nu), (\eta, \mu)) = \kappa(\nu, \mu)_{H^{1}(0,T)} + \left( {}^{\flat}\Psi_{T}\begin{pmatrix} \zeta \\ \nu \end{pmatrix}, {}^{\flat}\Psi_{T}\begin{pmatrix} \eta \\ \mu \end{pmatrix} \right)_{L^{2}(0,T)} \\ ((\mathrm{Id} - e^{TA})\zeta + (\Lambda\nu)(T), (\mathrm{Id} - e^{TA})\eta + (\Lambda\mu)(T))_{L^{2}(0,T)} $$

Therefore, ${}^{\flat}\sigma$ is positive. Further, let us consider $(\zeta, \nu)$ such that ${}^{\flat}\sigma((\zeta, \nu), (\zeta, \nu)) = 0$. We then have $\nu = 0$, and $(\mathrm{Id} - e^{TA})\zeta = 0$. As $(\mathrm{Id} - e^{TA})$ is invertible, $\zeta = 0$. Therefore,
$^b\sigma$ is definite and $\mathcal{J}_{\kappa,\epsilon}$ is convex. Finally we show that $^b\sigma$ is coercive by contradiction.
Then, let assume there exists a sequence $(\zeta_n, \nu_n)_{n\in\mathbb{N}}$ such that

$$ \lim_{n \to \infty} {}^b\sigma((\zeta_n, \nu_n), (\zeta_n, \nu_n)) = 0 \text{ whereas } \|\zeta_n\|_{\mathcal{Z}}^2 + \|\nu_n\|_{\text{H}^1(0,T)}^2 = 1. $$

From the definition of $^b\sigma$, we have

$$ \lim_{n \to \infty} \|\nu_n\|_{\text{H}^1(0,T)}^2 = 0. $$

Then, we find that

$$ \lim_{n \to \infty} \|(\text{Id} - e^{TA})\zeta_n\|_{\mathcal{Z}} = \lim_{n \to \infty} \|\Lambda(T)\nu_n\|_{\mathcal{Z}}. $$

As we have already noticed that there exists a constant $c_{st}$ such that

$$ \|\Lambda(T)\nu_n\|_{\mathcal{Z}} \le c_{st} \|\nu_n\|_{\text{L}^2(0,T)}, $$

and $(\text{Id} - e^{TA})$ has a bounded inverse, we conclude that

$$ \lim_{n \to \infty} \|\zeta_n\|_{\mathcal{Z}} = 0 \text{ and } \lim_{n \to \infty} \|\zeta_n\|_{\mathcal{Z}} = 1, $$

which is absurd. Therefore, $^b\sigma$ is coercive and $\mathcal{J}_{\kappa,\epsilon}$ is a strongly-convex functional. ■

**Proof of Theorem 4.9.** The proof of this theorem relies on two preliminary propositions. To state this preliminary results we introduce

$$ \mathcal{P}(\zeta, \nu) = \|z_{|\zeta,\nu}(T) - z_{|\zeta,\nu}(0)\|_{\mathcal{Z}}^2, $$

hence

$$ \mathcal{J}_{\kappa,\epsilon}(\zeta, \nu) = \epsilon^{-1} \mathcal{P}(\zeta, \nu) + \mathcal{J}_{\kappa}(\zeta, \nu). $$

By a direct application of classical results in convex optimization, one can show the following propositions.

**Proposition** We have the following bounds

$$ \mathcal{J}_{\kappa, \epsilon_n}(\bar{\zeta}_n, \bar{\nu}_n) \le \mathcal{J}_{\kappa, \epsilon_{n+1}}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}), \qquad (B.1) $$

and

$$ \mathcal{P}(\bar{\zeta}_n, \bar{\nu}_n) \ge \mathcal{P}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}). \qquad (B.2) $$

**Proof:** First, we have

$$ \begin{aligned} \mathcal{J}_{\kappa, \epsilon_n}(\bar{\zeta}_n, \bar{\nu}_n) &\le \mathcal{J}_{\kappa, \epsilon_n}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}) = \epsilon_n^{-1} \mathcal{P}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}) + \mathcal{J}_{\kappa}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}) \\ &\le \epsilon_{n+1}^{-1} \mathcal{P}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}) + \mathcal{J}_{\kappa}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}) = \mathcal{J}_{\kappa, \epsilon_{n+1}}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}). \end{aligned} $$

Second, we remark that

$$ \epsilon_n^{-1} \mathcal{P}(\bar{\zeta}_n, \bar{\nu}_n) + \mathcal{J}_{\kappa}(\bar{\zeta}_n, \bar{\nu}_n) \le \epsilon_n^{-1} \mathcal{P}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}) + \mathcal{J}_{\kappa}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}) $$

while

$$ \epsilon_{n+1}^{-1} \mathcal{P}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}) + \mathcal{J}_{\kappa}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}) \le \epsilon_{n+1}^{-1} \mathcal{P}(\bar{\zeta}_n, \bar{\nu}_n) + \mathcal{J}_{\kappa}(\bar{\zeta}_n, \bar{\nu}_n), $$
Adding the last two inequalities leads to

$$ (\epsilon_{n+1}^{-1} - \epsilon_n^{-1}) \mathcal{P}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}) \leq (\epsilon_{n+1}^{-1} - \epsilon_n^{-1}) \mathcal{P}(\bar{\zeta}_n, \bar{\nu}_n), $$

and gives (B.2).

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

**Proposition** We have

$$ \mathcal{J}_{\kappa}(\bar{\zeta}_n, \bar{\nu}_n) \leq \mathcal{J}_{\kappa}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}), \qquad (B.3) $$

and

$$ \mathcal{J}_{\kappa}(\bar{\zeta}, \bar{\nu}) \geq \mathcal{J}_{\kappa, \epsilon_n}(\bar{\zeta}_n, \bar{\nu}_n) \geq \mathcal{J}_{\kappa}(\bar{\zeta}_n, \bar{\nu}_n), \qquad (B.4) $$

where $(\bar{\zeta}, \bar{\nu}) = \arg \min_{(\zeta, \nu) \in \mathcal{K}} \mathcal{J}_{\kappa}$ with $\mathcal{K}$ defined by (61).

**Proof:** From (B.1) and (B.2) we deduce that

$$ \begin{aligned} \mathcal{J}_{\kappa}(\bar{\zeta}_n, \bar{\nu}_n) &= \mathcal{J}_{\kappa, \epsilon_n}(\bar{\zeta}_n, \bar{\nu}_n) - \epsilon_n^{-1} \mathcal{P}(\bar{\zeta}_n, \bar{\nu}_n) \\ &\leq \mathcal{J}_{\kappa, \epsilon_n}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}) - \epsilon_n^{-1} \mathcal{P}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}) \\ &\leq \mathcal{J}_{\kappa, \epsilon_n}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}) - \epsilon_{n+1}^{-1} \mathcal{P}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}) = \mathcal{J}_{\kappa}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1}). \end{aligned} $$

Then, we write

$$ \mathcal{J}_{\kappa}(\bar{\zeta}, \bar{\nu}) = \epsilon_n^{-1} \mathcal{P}(\bar{\zeta}, \bar{\nu}) + \mathcal{J}_{\kappa}(\bar{\zeta}, \bar{\nu}) \geq \mathcal{J}_{\kappa, \epsilon_n}(\bar{\zeta}_n, \bar{\nu}_n) \geq \mathcal{J}_{\kappa}(\bar{\zeta}_n, \bar{\nu}_n). $$

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

From (B.4), the sequence $(\mathcal{J}_{\kappa, \epsilon_n}(\bar{\zeta}_n, \bar{\nu}_n))$ is bounded. Therefore, $(\bar{\nu}_n)_{n \in \mathbb{N}}$ is bounded in $H^1(0, T)$ hence we can extract a subsequence weakly converging in $H_\#^1(0, T)$ and strongly converging in $L^2(0, T)$. Further, From (B.2) the sequence $\mathcal{P}_{\epsilon_{n+1}}(\bar{\zeta}_{n+1}, \bar{\nu}_{n+1})$ is decreasing, hence bounded. As $\epsilon_n^{-1}$ tends to infinity, this implies that

$$ \lim_{n \to \infty} \|(Id - e^{TA})\bar{\zeta}_n + (\Lambda \bar{\nu}_n)(T)\|_{\mathcal{Z}}^2 = 0. \qquad (B.5) $$

Therefore, we find that

$$ \lim_{n \to \infty} \|(Id - e^{TA})\bar{\zeta}_n\|_{\mathcal{Z}} = \lim_{n \to \infty} \|\Lambda(T)\bar{\nu}_n\|_{\mathcal{Z}}. \qquad (B.6) $$

By admissibility

$$ \|\Lambda(T)\bar{\nu}_n\|_{\mathcal{Z}} \leq c_{st} \|\bar{\nu}_n\|_{L^2(0,T)}, $$

we have that $((Id - \Phi(T))\bar{\zeta}_n)_{n \in \mathbb{N}}$ converges weakly and in a norm in $\mathcal{Z}$, hence is a strongly convergent sequence in $\mathcal{Z}$. By consequence, $\bar{\zeta}_n$ strongly converges in $\mathcal{Z}$. To summarize, we have proven that $(\bar{\zeta}_n, \bar{\nu}_n)_{n \in \mathbb{N}}$ is strongly convergent in $\mathcal{Z} \times L^2(0, T)$ and we denote the limit $(\bar{\zeta}_\infty, \bar{\nu}_\infty)$.
The sequence $(\bar{\zeta}_n, \bar{\nu}_n)_{n \in \mathbb{N}}$ being in particular weakly convergent in $\mathcal{Z} \times H_\#^1(0, T)$, $\mathcal{J}_\kappa$ being convex and lower-semicontinuous (hence weakly lower-semicontinuous), we deduce that

$$ \liminf_{n \to \infty} \mathcal{J}_\kappa(\bar{\zeta}_n, \bar{\nu}_n) \ge \mathcal{J}_\kappa(\bar{\zeta}_\infty, \bar{\nu}_\infty). \quad (B.7) $$

The associated trajectory $z_{|\bar{\zeta}_\infty, \bar{\nu}_\infty}$ is periodic due to (B.5), hence $(\bar{\zeta}_\infty, \bar{\nu}_\infty) \in \mathcal{K}$. As a consequence of (B.7) and (B.4), we obtain that

$$ \mathcal{J}_\kappa(\bar{\zeta}_\infty, \bar{\nu}_\infty) \le \mathcal{J}_\kappa(\bar{\zeta}, \bar{\nu}) \Rightarrow (\bar{\zeta}_\infty, \bar{\nu}_\infty) = (\bar{\zeta}, \bar{\nu}). $$

To conclude we need to show that the sequence $\bar{\nu}_n$ strongly converges in $H^1(0, T)$. From (B.4) we deduce that

$$ \liminf_{n \to \infty} \mathcal{J}_{\kappa, \epsilon_n}(\bar{\zeta}_n, \bar{\nu}_n) = \limsup_{n \to \infty} \mathcal{J}_{\kappa, \epsilon_n}(\bar{\zeta}_n, \bar{\nu}_n) = \mathcal{J}_{\kappa}(\bar{\zeta}, \bar{\nu}), $$

which shows that the limit of $\mathcal{J}_{\kappa, \epsilon_n}(\bar{\zeta}_n, \bar{\nu}_n)$ is unique and given by $\mathcal{J}_\kappa(\bar{\zeta}, \bar{\nu})$. Moreover

$$ \begin{aligned} \mathcal{J}_\kappa(\bar{\zeta}, \bar{\nu}) &= \lim_{n \to \infty} \mathcal{J}_{\kappa, \epsilon_n}(\bar{\zeta}_n, \bar{\nu}_n) = \liminf_{n \to \infty} (\epsilon_n^{-1} \mathcal{P}(\bar{\zeta}_n, \bar{\nu}_n) + \mathcal{J}_\kappa(\bar{\zeta}_n, \bar{\nu}_n)) \\ &\ge \liminf_{n \to \infty} \epsilon_n^{-1} \mathcal{P}(\bar{\zeta}_n, \bar{\nu}_n) + \liminf_{n \to \infty} \mathcal{J}_\kappa(\bar{\zeta}_n, \bar{\nu}_n) \\ &\ge \liminf_{n \to \infty} \epsilon_n^{-1} \mathcal{P}(\bar{\zeta}_n, \bar{\nu}_n) + \mathcal{J}_\kappa(\bar{\zeta}, \bar{\nu}) \end{aligned} $$

hence $\liminf \mathcal{P}(\bar{\zeta}_n, \bar{\nu}_n) = 0$ and since $(\mathcal{P}(\bar{\zeta}_n, \bar{\nu}_n))_{n \ge 0}$ is a decaying sequence, it converges and

$$ \lim_{n \to \infty} \mathcal{P}(\bar{\zeta}_n, \bar{\nu}_n) = 0. \quad (B.8) $$

We conclude by observing that the convergence of $(\bar{\zeta}_n, \bar{\nu}_n)_{n \in \mathbb{N}}$ in $\mathcal{Z} \times L^2(0, T)$ implies

$$ \lim_{n \to \infty} \|y_\delta - {}^b\Psi_T \bar{\nu}_n\|_{L^2(0,T)}^2 = \|y_\delta - {}^b\Psi_T \bar{\nu}\|_{L^2(0,T)}^2, $$

hence, with this equation, (B.8), and the property that $\mathcal{J}_{\kappa, \epsilon_n}(\bar{\zeta}_n, \bar{\nu}_n)$ converges towards $\mathcal{J}_\kappa(\bar{\zeta}, \bar{\nu})$ we obtain that

$$ \lim_{n \to \infty} \|\bar{\nu}_n\|_{H^1(0,T)} = \|\bar{\nu}\|_{H^1(0,T)}. $$

The convergence of the norm implies here that $\bar{\nu}_n$ strongly converge in $H^1(0, T)$. ■

## Appendix C. A discrete-time optimal control approach

*Optimality system at the discrete level* – The minimization of the discretized $\mathcal{J}_{\kappa, \epsilon}^N$ under the finite dimensional constraint (76) can be derived using the associated Lagrangian

$$ \begin{aligned} \mathcal{L}((z_h^n)_{0 \le n \le N_T}, (q_h^n)_{0 \le n \le N}, (\nu^n)_{1 \le n \le N_T}) ={}& \epsilon^{-1} \|z_h^{N_T} - z_h^0\|_{\mathcal{Z}}^2 \\ & + \frac{\tau}{2} \sum_{n=0}^{N_T-1} \left[ |y_\delta^n - Cz_h^n|^2 + \kappa \left( |\nu^n|^2 + \left| \frac{\nu^{n+1} - \nu^n}{\tau} \right|^2 \right) \right] \\ & + \sum_{n=0}^{N_T-1} (q_h^{n+1}, z_h^{n+1} - \Phi_{h,\tau} z_h^n - B_{h,\tau} \nu^{n+1})_{\mathcal{Z}}. \end{aligned} \quad (C.1) $$
From Kuhn-Tucker Theorem in finite dimension, we now seek for the saddle-point of the Lagrangian. Writing the optimality condition leads to, for $0 < n < N_T$,

$$ 
\begin{aligned}
\forall \eta_h \in \mathcal{V}_h, \quad & 0 = \langle \mathrm{D}_{z^n} \mathcal{L}((\bar{z}_h^n)_{0\le n\le N}, (\bar{q}_h^n)_{0\le n\le N}, (\bar{\nu}^n)_{1\le n\le N}), \eta_h \rangle \\
& = -\tau(y_\delta^n - C\bar{z}_h^n, C\eta) + (\bar{q}_h^n, \eta_h) - (\bar{q}_h^{n+1}, \Phi_{h,\tau}\eta_h),
\end{aligned} 
\quad (C.2) $$

namely

$$ \bar{q}_h^n = \Phi_{h,\tau}^* \bar{q}_h^{n+1} + \tau C^*(y_\delta^n - C\bar{z}_h^n), \quad 0 < n < N_T. \quad (C.3) $$

Then, from

$$ 
\begin{aligned}
\forall \eta \in \mathcal{V}_h, \quad & 0 = \langle \mathrm{D}_{z^n} \mathcal{L}((\bar{z}_h^n)_{0\le n\le N}, (\bar{q}_h^n)_{0\le n\le N}, (\bar{\nu}^n)_{1\le n\le N}), \eta_h \rangle \\
& = \epsilon^{-1}(\bar{z}_h^{N_T} - \bar{z}_h^0, \eta_h)_Z + (\bar{q}_h^{N_T}, \eta_h)_Z,
\end{aligned} 
\quad (C.4) $$

we deduce

$$ \bar{q}_h^{N_T} = -\epsilon^{-1}(\bar{z}_h^{N_T} - \bar{z}_h^0). $$

Further, we have also

$$ 
\begin{aligned}
\forall \eta \in \mathcal{V}_h, \quad & 0 = \langle \mathrm{D}_{z^0} \mathcal{L}((\bar{z}_h^n)_{0\le n\le N}, (\bar{q}_h^n)_{0\le n\le N}, (\bar{\nu}^n)_{1\le n\le N}), \eta \rangle \\
& = -\epsilon^{-1}(\bar{z}^{N_T} - \bar{z}^0, \eta_h)_Z - \tau(y_\delta^0 - C\bar{z}^0, C\eta_h)_Z - (\bar{q}^1, \Phi\eta_h)_Z.
\end{aligned} 
\quad (C.5) $$

Therefore, by defining $\bar{q}^0 = \Phi_{h,\tau}^* \bar{q}_h^1 + \tau C^*(y_\delta^0 - C\bar{z}_h^0)$, namely considering (C.3) also for $n=0$, we find that

$$ \bar{q}_h^0 = \epsilon^{-1}(\bar{z}_h^{N_T} - \bar{z}_h^0). \quad (C.6) $$

Finally, by using an Abel transform with $\nu^0 = \nu^{N_T}$, we have the classical identity

$$ 
\begin{aligned}
\sum_{n=0}^{N_T-1} \left( \frac{\nu^{n+1} - \nu^n}{\tau} \right) \left( \frac{\mu^{n+1} - \mu^n}{\tau} \right) &= \sum_{n=1}^{N_T} \left( \frac{\nu^{n+1} - 2\nu^n + \nu^{n-1}}{\tau^2} \right) \mu^n \\
&= \sum_{n=1}^{N_T} \nu^n \left( \frac{\mu^{n+1} - 2\mu^n + \mu^{n-1}}{\tau^2} \right),
\end{aligned} $$

which leads to the partial derivative

$$ 
\begin{aligned}
0 &= \partial_{\nu^n} \mathcal{L}((\bar{z}_h^n)_{0\le n\le N_T}, (\bar{q}_h^n)_{0\le n\le N_T}, (\bar{\nu}^n)_{1\le n\le N}), \\
&= \kappa \left( \nu^n + \frac{\nu^{n+1} - 2\nu^n + \nu^n}{\tau^2} \right) - B_{h,\tau}^* \bar{q}^n.
\end{aligned} 
\quad (C.7) $$

We consequently introduce the linear operator $K_\sharp \in \mathcal{M}_{N_T}(\mathbb{R})$ such that for any vector $(\mu_n)_{n \in [1, N_T]}$ and $(\nu_n)_{n \in [1, N_T]}$ such that $\mu_0 = \mu_{N_T}$ and $\nu_0 = \nu_{N_T}$.

$$ ((\mu_n)_{n \in [1, N_T]}, K_\sharp(\nu_n)_{n \in [1, N_T]})_{\ell^2([1, N_T])} = - \sum_{n=1}^{N_T} \left( \frac{-\nu^{n+1} + 2\nu^n - \nu^n}{\tau^2} \right) \mu^n. $$

We ultimately obtain the optimality system

$$ 
\left\{
\begin{array}{l}
\bar{z}_h^{n+1} = \Phi_{h,\tau} \bar{z}_h^n + B_{h,\tau} ((\mathrm{Id} + K_\sharp)^{-1} (B_{h,\tau} \bar{q}_h^n)_{n \in [1,N]})^{n+1}, \quad 0 \le n \le N_T - 1 \\
\bar{q}_h^n = \Phi_{h,\tau}^* \bar{q}_h^{n+1} + \tau C^*(y_\delta^n - C\bar{z}_h^n), \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad 0 \le n \le N_T - 1 \\
\bar{z}_h^0 = \bar{z}_h^{N_T} + \epsilon \bar{q}_h^0, \\
\bar{q}_h^{N_T} = -\epsilon^{-1}(\bar{z}_h^{N_T} - \bar{z}_h^0),
\end{array}
\right. 
\quad (C.8) $$
A gradient descent approach – Solving the optimality system (C.8) implies to solve a fully coupled space-time problem which can be overcome, for large-dimensional discretization, using an iterative procedure based on a descent-gradient minimization of the associated functional $\mathcal{J}_{\kappa, \epsilon}^{N_T}$. We easily obtain from the previous computation on the Lagrangian that the gradient of $\mathcal{J}_{\kappa, \epsilon}^{N_T}$ is given by

$$ \nabla_{\zeta} \mathcal{J}_{\kappa, \epsilon}^{N_T} = \epsilon^{-1}(\bar{z}_h^0 - \bar{z}_h^{N_T}) - q_{h|\zeta,\nu}^0 $$

and

$$ \nabla_{\nu^n} \mathcal{J}_{\kappa, \epsilon}^{N_T} = \kappa \left( \nu^n + \frac{\nu^{n+1} - 2\nu^n + \nu^n}{\tau^2} \right) - B_{h,\tau}^* q_{h|\zeta,\nu}^n $$

with $(q_{h|\zeta,\nu}^n)_{1\le n \le N_T}$ computed from a given trajectory $(z_{h|\zeta,\nu}^n)$ by following adjoint equation

$$ \left\{ \begin{aligned} q_{h|\zeta,\nu}^n &= \Phi_{h,\tau}^* \bar{q}_h^{n+1} + \tau C^* (y_\delta^n - C \bar{z}_{h|\zeta,\nu}^n) \\ \bar{q}_{h|\zeta,\nu}^{N_T} &= -\epsilon^{-1} (z_{h|\zeta,\nu}^{N_T} - z_{h|\zeta,\nu}^0) \end{aligned} \right. \qquad (C.9) $$

Therefore, the solution of two-end problem (C.8) is obtained as the limit of the iterative back and forward scheme: given a small enough relaxation sequence $(\rho^j)$, we compute at each iteration $j \in \mathbb{N}$

$$ \left\{ \begin{aligned} z_h^{n+1, j+1} &= \Phi_{h, \tau} z_h^{N_T} + B_{h, \tau} (1 - \rho^j) \nu^{n+1} \\ & \quad + \rho^j ((\mathrm{Id} + K_\#)^{-1} (B_{h, \tau} q_h^{n,j})_{n \in [1, N]})^{n+1}, \quad n \ge 0 \quad (C.10) \\ z_h^{0, j+1} &= (1 - \rho^j) z_h^{0, j} + \rho^j (\bar{z}_h^{N_T} + \epsilon q_{h|\zeta,\nu}^{N_T}) \end{aligned} \right. $$

followed by

$$ \left\{ \begin{aligned} \bar{q}_h^{n,j+1} &= \Phi_{h,\tau}^* \bar{q}_h^{n+1,j+1} + \tau C^* (y_\delta^n - C \bar{z}_h^{n,j+1}) \\ q_h^{N_T,j+1} &= -\epsilon^{-1} (z_{h|\zeta,\nu}^{N_T,j+1} - z_h^{0,j+1}) \end{aligned} \right. \qquad (C.11) $$

This can be implemented using the following algorithm:

<caption>**Algorithm 1:** Gradient descent algorithm</caption>

Choose the initial guess $z_h^{0,0}$ and $(\nu^{n,0})_{1\le n \le N_T}$;
while $z_h^{0,j}$ and $(\nu^{n,j})_{1\le n \le N_T}$ are not converged do
Fix $\rho^j$ from the selected descent strategy;
Compute (C.10);
Compute (C.11);

end

However, it is well-known that solving a penalized optimisation problem using a gradient descent approach is unfavorable as, by nature, the underlying linear system to be solved faces ill-conditioning [21].