Optimal Boundary Control Problem for the Stokes Equation in Cardiovascular Applications

1. Introduction

Optimal control problems of systems governed by partial differential equations (PDEs) have garnered significant attention in recent decades. The Stokes equations, a simplified version of the Navier-Stokes equations, are pivotal in describing the flow behavior of incompressible fluids under low Reynolds number conditions where viscous forces dominate inertial forces. These problems may be applied in different science and engineering fields such as atmospheric and ocean sciences, aeronautics, and industrial design. In [15,25], several relevant topics related to both theoretical and numerical aspects may be found, as a result of, at least, a decade of research by different authors. More recently, the application of flow control to computational blood flow modeling has also been investigated. Examples can be seen in [4,14], for the stationary case, or [7] for time-dependent equations.

From the applications point of view, the Stokes and Navier-Stokes equations have been complemented with Dirichlet boundary conditions, mixed with imposed stress conditions, also known as traction boundary conditions. Nevertheless, as it was indicated in [8], for the purpose of considering the coupling of the fluid with a model representing the blood vessel wall, imposing the so-called total stress vector, instead of the usual stress vector, should be preferred.

Optimal control problems, in general, aim to determine control inputs that optimize a given performance criterion while ensuring that the system dynamics adhere to certain constraints. When these systems are governed by partial differential equations (PDEs) like the Stokes equations, the complexity of the control problems increases significantly. In particular, boundary optimal control problems involve finding controls applied on the boundary of the domain to achieve desired outcomes within the fluid flow.

Mixed boundary conditions, combining Dirichlet and stress (Neumann) conditions, add another layer of complexity to these control problems. While Dirichlet conditions, which specify the value of the function on the boundary, are well-understood and lend themselves to classical regularity results, the introduction of stress conditions complicates the analysis. Mixed boundary conditions often lead to a loss of regularity that poses significant challenges in both the theoretical and numerical study of these problems. In this sense, following the approach from [1,2], they employ mixed boundary conditions, including the homogeneous Dirichlet condition for the pressure on a portion of the boundary of the flow domain; but at the same time, they formulate the optimal control problem without using the curl operator and boundary conditions associated with this operator. However, in [1] they take into account the dependence of both the viscosity and the thermal conductivity coefficient on the temperature. This expands the range of possible applications of the results. But, in our case, we consider a phisical model based on the Stokes equations where the viscosity is a real positive constant and we do not consider viscosity to be dependent on temperature. This is justified by the specific medical applications we are focusing on, where the viscosity of blood can be considered constant. The temperature variations in these scenarios are minimal and do not significantly affect the viscosity. Therefore, assuming constant viscosity simplifies the model without compromising the accuracy of the results within the application we want to study.

In this work, we analyse the boundary optimal control problem associated with the time dependent Stokes equation under mixed Dirichlet and total stress boundary conditions. The boundary of the domain is assumed to be composed of three distinct components. A fixed total-stress condition is assumed on one part, a homogeneous Dirichlet condition on another and, finally, a Dirichlet-type control is assumed to act on the remaining third component. This configuration is inspired by the applications addressed in [7,14]. There, a Dirichlet control was applied at the inlet boundary, stress conditions were considered on the outlet boundaries, whereas homogeneous Dirichlet conditions were associated with the physical no-slip assumption. In those works, the stress boundary condition is to be understood as a prescribed force per unit area, computed as the normal component of the Cauchy stress tensor.

The application of the Stokes equations in cardiovascular problems is particularly noteworthy. In modeling blood flow in microcirculation and small arteries, where the flow is typically slow and viscosity effects are pronounced, the Stokes equations provide a more appropriate model than the full Navier-Stokes equations. In such contexts, optimal control techniques are used to simulate and control blood flow, aiming to improve the understanding and treatment of cardiovascular diseases.

Recent studies have applied these principles to computational blood flow modeling. For instance, Dirichlet boundary conditions are used to model the no-slip condition on vessel walls, while stress conditions are applied at the interfaces where blood interacts with medical devices or at the outlet boundaries of the vascular network, see [13]. These mixed boundary conditions mimic physiological conditions more accurately but also introduce additional analytical challenges due to the loss of regularity.

The coupling of fluid flow with models representing the artery walls, often referred to as fluid-structure interaction (FSI) problems, further complicates the scenario. Here, the use of total stress boundary conditions, which account for both pressure and viscous stress, becomes crucial. These conditions are essential for accurately modeling the interaction between blood flow and arterial walls but require advanced mathematical techniques to handle the resulting optimal control problems.

Boundary control problems, under the stationary assumption, have long been studied by different authors. Concerning mixed boundary conditions of Dirichlet-Neumann (stress vector) type, in [11] a Neumann-type control problem is studied. There, the goal was the drag minimization, under state constraints, on a two-dimensional exterior problem. Following some of those ideas, in [13], the existence of an optimal boundary control was established under different types of cost functionals, motivated by some of the aforementioned applications. In [23], the authors used a Lagrange multipliers approach to address the control of the Dirichlet boundary, still in the stationary case.

Concerning time-dependent problems, Dirichlet controls have been considered in [18,21,22] for bounded domains in 2D. In [10], the authors considered a Dirichlet control acting on an interior boundary of an unbounded 3D domain. The case of mixed boundary conditions, particularly of the case of Dirichlet-Total Stress remains, up to our knowledge, to be analysed.

Our work aims to provide a deeper understanding of mixed boundary optimal control problems for the time-dependent Stokes equations, particularly in the context of cardiovascular applications. We draw inspiration from established methods in the literature, such as those used for Dirichlet controls in bounded domains and mixed Dirichlet-Neumann conditions. However, we extend these approaches to handle the more complex mixed boundary conditions of Dirichlet and Total Stress. This extension is non-trivial and requires careful analysis to ensure the existence and regularity of solutions, addressing the unique challenges posed by the loss of classical regularity, see [24].

In this sense, our domain $\Omega$ is an idealized domain representing a partially obstructed vessel capillary. The domain $\Omega$ is represented in Figure 1. In this case, the left side is the inlet boundary (where we imposed the control) while the right side corresponds to outlet boundary (where we imposed the total stress condition). The remaining boundaries represent the vessel wall where we imposed homogeneous Dirichlet conditions.

Moreover, this work aims to serve as a foundational study for future research in the control of the Navier-Stokes equations and fluid-structure interaction problems. By establishing a robust theoretical and numerical framework for the Stokes equation, we try to pave the way for addressing these more complex and realistic problems, ultimately contributing to the advancement of optimal control techniques in cardiovascular applications.

The plan of this paper reads as follows. In Section 2 we present an analysis of the direct problem, i.e., the Stokes equation with mixed boundary conditions. This section provides the mathematical setting for the analysis of the optimal control problem, which will be the subject of Section 3. In the later, existence of solution and first order conditions are studied. In Section 4, we provide numerical results and we finish by summarizing some conclusions and perspectives.

2. Unsteady Stokes Equation

In this section, we consider the direct problem, that is, the non-stationary Stokes equation with mixed boundary conditions of Dirichlet and total stress types.

Despite their relevance when modeling viscous incompressible flows, the Stokes equation, endowed with such type of mixed boundary conditions, have been seldom analysed from a theoretical standpoint. An exception are the results found in [3,20], for the homogeneous case. As explained by the authors there, one of the advantages of this configuration is the possibility of deriving an energy inequality, contrary to the configuration based on the ubiquitously used stress (Neumann) boundary conditions, for which no control over the energy flux can be achieved.

We can describe our problem as follows. Given functions $f,g,h$ and the initial condition $u_0$, we want to solve system

$$\left\{\begin{array}{cc}{\displaystyle {\partial }_{t}u-\nabla ·\left(2\nu D\left(u\right)\right)+\nabla p=f}\hfill & (0,T)\times \Omega ,\hfill \\ \nabla ·u=0\hfill & (0,T)\times \Omega ,\hfill \\ u=g\hfill & (0,T)\times {\Gamma }_{in},\hfill \\ u=0\hfill & (0,T)\times {\Gamma }_{\omega },\hfill \\ \left(2\nu D\left(u\right)-pI\right)·\overline{n}=h,u\times \overline{n}=0\hfill & (0,T)\times {\Gamma }_{out},\hfill \\ u(0,x)=u_0\left(x\right)\hfill & \Omega \hfill \end{array}\right.$$

(1)

where the unknowns are the fluid velocity u and the pressure p, $\overline{n}$ is the exterior unit normal to ${\Gamma }_{out}$ and $\partial \Omega =\overline{{\Gamma }_{in}}\cup \overline{{\Gamma }_{out}}\cup \overline{{\Gamma }_{\omega }}$ is the boundary of an open bounded domain $\Omega \subset {\mathbb{R}}^N$ with $N=2,3$.

Now, we provide a brief nomenclature below that defines the physical parameters and symbols used throughout this paper. This nomenclature is try to help readers fully understand the equations and the model presented in the study and to apply the study to the medical cases.

• $\nu :$ Kinematic viscosity of the fluid.
• ${\displaystyle D\left(u\right)=\frac{1}{2}\left(\nabla u+\nabla u^T\right)}$ represents the strain rate tensor.
• $f:$ external force term applied to the fluid.
• g: control velocity condition at the inlet ${\Gamma }_{in}.$
• $\overline{n}$: unit normal vector to the boundary.
• h: boundary condition related to stress on ${\Gamma }_{out}$.
• $u_0$ initial velocity.

Under the above assumptions, we can say that system (1) models an incompressible Newtonian fluid at a constant temperature.

The non-homogeneous Dirichlet condition for the boundary ${\Gamma }_{in}$ corresponds to what will later be considered as the control action. Also, there is another boundary component, ${\Gamma }_{out}$, where we apply a prescribed total stress condition. Finally, in the remaining part of the boundary ${\Gamma }_{\omega }$ we impose a homogeneous Dirichlet condition.

Therefore, whenever required, we denote by ${\Gamma }_D={\Gamma }_{in}\cup {\Gamma }_{\omega }$ and ${\Gamma }_N={\Gamma }_{out}$ as a reference to Dirichlet (essential) and Neumann (natural) boundary conditions, respectively.

2.1. Functional Spaces and Preliminary Results

In this section we introduce some functional spaces along with several properties and preliminary results that will be required later on, for our subsequent analysis.

In the sequel, let $\Omega$ be a connected open set of ${\mathbb{R}}^N$$N=2,3$, with a locally Lipschitz boundary $\partial \Omega$ composed by three smooth open subsets, mutually disjoint, denoted by ${\Gamma }_{in}$, ${\Gamma }_{\omega },$${\Gamma }_{out}$, such that:

$$\begin{array}{c}{\Gamma }_{in},{\Gamma }_{\omega },{\Gamma }_{out}areofclass{\mathcal{C}}^{1,1}\hfill \\ \partial \Omega =\overline{{\Gamma }_{in}}\cup \overline{{\Gamma }_{\omega }}\cup \overline{{\Gamma }_{out}}.\hfill \end{array}$$

Following the classical analysis for the Stokes equation (see [26]), let us consider:

$$\mathcal{D}\left(\overline{\Omega }\right)=\{\phi \in {\mathcal{C}}^{\infty }{\left(\overline{\Omega }\right)}^2:\nabla ·\phi =0in\Omega ,supp\phi \cap {\Gamma }_{\omega }=\varnothing ,\phi \times \overline{n}=0in{\Gamma }_{out}\}$$

and

$$\mathcal{E}\left(\overline{\Omega }\right)=\{\phi \in {\mathcal{C}}^{\infty }{\left(\overline{\Omega }\right)}^2:\nabla ·\phi =0in\Omega ,supp\phi \cap {\Gamma }_D=\varnothing ,\phi \times \overline{n}=0in{\Gamma }_{out}\}.$$

For convenience, we introduce the classical space V and H for Stokes equations, but we adapt the definitions to the boundary conditions:

$$\mathbb{V}:=\{u\in H^1{(\Omega )}^2:\nabla ·u=0in\Omega ,{tr}_{\omega }\left(u\right)=0,u\times \overline{n}=0in{\Gamma }_{out}\}$$

$$\left(ortheclousureof\mathcal{D}\left(\overline{\Omega }\right)in{H}^1{(\Omega )}^2\right),$$

$$V:=\{u\in H^1{(\Omega )}^2:\nabla ·u=0in\Omega ,{tr}_D\left(u\right)=0,u\times \overline{n}=0in{\Gamma }_{out}\}$$

$$\left(ortheclousureof\mathcal{E}\left(\overline{\Omega }\right)in{H}^1{(\Omega )}^2\right)$$

where ${tr}_K$ is the trace operator to a subset ${\Gamma }_K.$

Also, we consider:

$$\mathbb{H}:={\overline{\mathcal{D}\left(\overline{\Omega }\right)}}^{L^2},H:={\overline{\mathcal{E}\left(\overline{\Omega }\right)}}^{L^2}$$

and

$$\mathbb{W}(0,T)=\{u\in L^2(0,T;\mathbb{V}),u^{\prime }\left(t\right)\in L^2(0,T;V^{*})\},$$

$$W(0,T)=\{u\in L^2(0,T;V),u^{\prime }\left(t\right)\in L^2(0,T;V^{*})\},$$

where ${\mathbb{V}}^{*}$ is the dual space of $\mathbb{V}$ and $V^{*}$ the dual of the $V.$

Lemma 1. 

The space $\mathbb{V}$ is dense in $\mathbb{H}$ and the space V is also dense in $H.$

Proof:

Thanks to the Poincaré inequality (it remains true when we have that the functional space $\mathbb{V}$ and V is null on a part of boundary), we have that $i_1:\mathbb{V}\hookrightarrow \mathbb{H}$ and $i_2:V\hookrightarrow H$ are continuous embeddings. Then, the conclusion follows from the classical results found in [26], Chapter 1.2.

□

Note that, from the definition of the functional spaces, and thanks to Riesz’s Theorem, we also have the following chain of continuous injections:

$$\mathbb{V}\hookrightarrow \mathbb{H}\hookrightarrow {\mathbb{V}}^{*}V\hookrightarrow H\hookrightarrow V^{*}.$$

In addition, by using Lemma 1.1 (Chapter III) of [26], we can ensure that $\mathbb{W}(0,T)\hookrightarrow {\mathcal{C}}^0([0,T];\mathbb{H})$ and $W(0,T)\hookrightarrow {\mathcal{C}}^0([0,T];H).$

In the sequel, we consider $|·|$ to be the usual $L^2(\Omega )$ norm, and we denote by $[·,·]$ the inner product in $L^2(\Omega )$. Also, we consider ${|·|}_1$ and ${|·|}_0$ to be the norma in $H^1(\Omega )$ and $H_0^1(\Omega )$ respectively, defined by

$${\displaystyle {\left|u\right|}_1={\left(\right|u|}_{L^2}^2+{\left|D\left(u\right)\right|}_{L^2}^2{)}^{\frac{1}{2}}{and\left|u\right|}_0={\left|D\left(u\right)\right|}_{L^2}.}$$

Besides, we denote by $(·,·)$ the usual $L^2(0,T;L^2(\Omega ))$ inner product. When convenient, we represent $L^2(0,T;V)$ and $L^2(0,T;\mathbb{V})$ by $L^2\left(V\right)$ and $L^2\left(\mathbb{V}\right)$, respectively, for example:

$${\displaystyle {\parallel f\parallel }_{L^2(0,T;\mathbb{V})}={\left({\int }_0^T{\left|D\left(f\right)\right|}_{L^2}^{2}dt\right)}^{\frac{1}{2}}.}$$

2.2. Existence and Uniqueness Result

This section is devoted to the existence of a solution for the Stokes equation (1). To this purpose, we start by introducing the weak formulation as well as the definition of a weak solution of these equations.

Definition 1. 

Given $f\in L^2((0,T)\times \Omega ),g\in L^2((0,T)\times {\Gamma }_{in}),h\in L^2((0,T)\times {\Gamma }_{out})$ and $u_0\in \mathbb{H},$ we say that u is a weak solution of the Stokes Equation (1) if

$$\left\{\begin{array}{c}u\in \mathbb{W}(0,T),\hfill \\ {\displaystyle {\langle u_t,v\rangle }_{L^2\left(V^{*}\right),L^2\left(V\right)}+2\nu (D\left(u\right),D\left(v\right))=}\hfill \\ (f,v)+{(h,v)}_{L^2(0,T;L^2\left({\Gamma }_{out}\right))}\forall v\in L^2\left(V\right),\hfill \\ u=g\mathrm{in}{L}^2((0,T)\times {\Gamma }_{in}),\hfill \\ u\left(0\right)=u_0\mathrm{in}\mathbb{H},\hfill \end{array}\right.$$

(2)

Notice that, the pressure is not appearing in the weak formulation. Nevertheless, as in the cases with a single Dirichlet condition is considered, we can recover the pressure using De Rham’s - type results ([26]). We could also consider a mixed weak formulation that includes pressure in the problem’s formulation. In this case, the solution space would change and would no longer include the condition of zero divergence. Accordingly, this equation would also be incorporated into the weak formulation as an additional equation, allowing us to solve for the pressure as well.

Because we are going to use the Dirichlet boundary condition as a control variable, let us assume some additional requirements for g. We recall the following result

Lemma 2. 

There exists a bounded extension operator $E:H_0^1\left({\Gamma }_{in}\right)\mapsto \mathbb{V}$, where

$$H_0^1\left({\Gamma }_{in}\right)=\{w\in L^2\left({\Gamma }_{in}\right):{\nabla }_{{\Gamma }_{in}}w\in L^2\left({\Gamma }_{in}\right)andw=0in\partial {\Gamma }_{in}\}$$

with ${\nabla }_{{\Gamma }_{in}}w$ the tangential gradient.

Note that, the tangential gradient ${\nabla }_{\Gamma }$ of f at a point $x\in \Gamma$ is defined as the tangential component of the gradient of the trace of f, i.e., the projection of the gradient of the trace of f onto the tangent plane to $\Gamma$ at x, that is,

$${\displaystyle {\nabla }_{\Gamma }f=\nabla \tilde{f}-(\nabla \tilde{f}·\overline{n})\overline{n}=\nabla \tilde{f}-\frac{\partial \tilde{f}}{\partial n}\overline{n}on\Gamma ,}$$

wher $\tilde{f}$ is the trace of f on $\Gamma$.

The proof of this lemma follows the proof of the Corollary 3.5 of [11].

In order to extend the E operator to the time domain, see that

Lemma 3. 

Let a real number $T>0$, then, there exists a bounded extension operator $E_T:L^2(0,T;H_0^1\left({\Gamma }_{in}\right))\mapsto L^2(0,T;\mathbb{V})$, with

$$E_T\left(v\left(t\right)\right)=E\left(v\left(t\right)\right)a.e.t\in [0,T].$$

Proof:

We use the definition of the operator $E_T$ to prove that it is bounded. Let $v\in L^2(0,T;H_0^1\left({\Gamma }_{in}\right))$,

$${\displaystyle \parallel E_T\left(v\right){\parallel }_{L^2(0,T;\mathbb{V})}={\left({\int }_0^T{\parallel E\left(v\left(t\right)\right)\parallel }_{\mathbb{V}}^{2}dt\right)}^{1/2}\le {\left({\int }_0^T{C}^2{\parallel v\left(t\right)\parallel }_{H_0^1\left({\Gamma }_{in}\right)}^{2}dt\right)}^{1/2})}$$

where $C>0$ doesn’t depend on t by the definition of bounded operator. So, we have that:

$$\parallel E_T{\left(v\right)\parallel }_{L^2(0,T;\mathbb{V})}\le C{\parallel v\parallel }_{L^2(0,T;H_0^1\left({\Gamma }_{in}\right))}.$$

Therefore, $E_T$ is well defined and is a bounded operator.

□

The existence of such extension operators allows us to define the following set:

$$\begin{array}{c}\mathcal{G}=\{g=t{r}_{(0,T)\times {\Gamma }_{in}}\left(g_{ext}\right),g_{ext}\in E_T\left(L^2(0,T;H_0^1\left({\Gamma }_{in}\right))\right)\cap \mathbb{W}(0,T),\\ {\displaystyle {\int }_{{\Gamma }_{in}}{g}_{ext}\overline{n}dx+{\int }_{{\Gamma }_{out}}{g}_{ext}\overline{n}dx=0\}}\end{array}$$

which characterizes the Dirichlet boundary conditions that we are going to consider.

At this point, we can state the following existence result

Theorem 1. 

Let $\Omega \subset {\mathbb{R}}^d$, $d=2,3$ be an open boundary set with a locally Lipschitz boundary, $g\in \mathcal{G}$, $h\in L^2((0,T)\times {\Gamma }_{out})$, $f\in L^2((0,T)\times \Omega )$ and $u_0\in \mathbb{H}$ with $u_0-g_{ext}\left(0\right)\in H.$ Then, there exists a unique weak solution $u\in \mathbb{W}(0,T)$ of the Stokes system (2), satisfying:

$$\begin{array}{c}{\displaystyle \underset{t\in [0,T]}{sup}{\left|u\left(t\right)\right|}_{\mathbb{V}}+{\parallel u\parallel }_{L^2\left(\mathbb{V}\right)}\le {C(\parallel f\parallel }_{L^2(0,T;\Omega )}+{\parallel h\parallel }_{L^2(0,T;L^2\left({\Gamma }_{out}\right))}}\hfill \\ +\parallel {\left(g_{ext}\right)}_t{\parallel }_{L^2\left(V^{*}\right)}+\parallel g_{ext}{\parallel }_{L^2\left(V\right)}+|u_0|+|g_{ext}\left(0\right)\left|\right)\mathrm{for} \mathrm{a} \mathrm{constant} C>0.\hfill \end{array}$$

(3)

Proof:

We start by remarking that there exists a function $g_{ext}\in E_T\left(L^2(0,T;H_0^1\left({\Gamma }_{in}\right))\right)\cap \mathbb{W}(0,T)$ with $g=t{r}_{(0,T)\times {\Gamma }_{in}}\left(g_{ext}\right)$ and $u_0\in \mathbb{H}$ with $u_0-g_{ext}\left(0\right)\in H$.

Then, $g_{ext}\left(0\right)\in H$ and for all $\psi \in L^2\left(V\right)$ we have that

$${\langle {\left(g_{ext}\right)}_t,\psi \rangle }_{L^2\left(V^{*}\right),L^2\left(V\right)}+2\nu (D\left(g_{ext}\right),D\left(\psi \right))\in L^2\left(V^{*}\right).$$

Next, consider the following homogeneous problem:

$$\left\{\begin{array}{c}\tilde{u}\in W(0,T),\hfill \\ {\displaystyle {\langle {\tilde{u}}_t,v\rangle }_{L^2\left(V^{*}\right),L^2\left(V\right)}+2\nu (D\left(\tilde{u}\right),D\left(v\right))=(f,v)+{(h,v)}_{L^2(0,T;L^2\left({\Gamma }_{out}\right))}-}\hfill \\ -{\langle {\left(g_{ext}\right)}_t,v\rangle }_{L^2\left(V^{*}\right),L^2\left(V\right)}-2\nu (D\left(g_{ext}\right),D\left(v\right))\forall v\in L^2\left(V\right),\hfill \\ \tilde{u}=0(0,T)\times {\Gamma }_{in},\hfill \\ \tilde{u}\left(0\right)=u_0-g_{ext}\left(0\right)\in H\Omega .\hfill \end{array}\right.$$

We use the compacity method and Galerkin approximation to prove the existence of the solution, as in Theorem 1.1 of Chapter III of [26].

So, we can ensure that

$$\begin{array}{c}{\displaystyle \underset{t\in [0,T]}{sup}|\tilde{u}\left(t\right)|+\parallel \tilde{u}{\parallel }_{L^2\left(V\right)}\le {C(\parallel f\parallel }_{L^2(0,T;\Omega )}+{\parallel h\parallel }_{L^2(0,T;L^2\left({\Gamma }_{out}\right))}}\hfill \\ +\parallel {\left(g_{ext}\right)}_t{\parallel }_{L^2\left(V^{*}\right)}+\parallel g_{ext}{\parallel }_{L^2\left(V\right)}+|u_0-g_{ext}\left(0\right)|)\mathrm{for} \mathrm{a} \mathrm{constant} C>0.\hfill \end{array}$$

Therefore, taking $u=\tilde{u}+g_{ext}$, by linearity, u is the unique solution of (1) and satifies:

$$\begin{array}{c}{\displaystyle \underset{t\in [0,T]}{sup}{\left|u\left(t\right)\right|}_{\mathbb{V}}+{\parallel u\parallel }_{L^2\left(\mathbb{V}\right)}\le {C(\parallel f\parallel }_{L^2(0,T;\Omega )}+{\parallel h\parallel }_{L^2(0,T;L^2\left({\Gamma }_{out}\right))}}\hfill \\ +\parallel {\left(g_{ext}\right)}_t{\parallel }_{L^2\left(V^{*}\right)}+\parallel g_{ext}{\parallel }_{L^2\left(V\right)}+|u_0|+|g_{ext}\left(0\right)\left|\right)\mathrm{for} \mathrm{a} \mathrm{constant} C>0.\hfill \end{array}$$

(4)

□

Note that, the solution $u\in \mathbb{W}(0,T)$ then u is equal a.e. to a continuous function $[0,T]\mapsto \mathbb{H}$, see [26]. This fully clarifies the initial condition in time.

Consequently, these results indicate that mathematically, there is a unique weak solution, which physically implies that the proposed problem is well-posed and can be solved. In other words, for any desired initial data, the Stokes equation can be used to approximate blood flow.

3. Optimal Control Problem

This section is devoted to the study of the optimal control problem associated to the Stokes Equation (2).

3.1. Existence of an Optimal Pair

Consider the quadratic functional defined by:

$$\begin{array}{c}{\displaystyle J(g,u):=\frac{{\alpha }_1}{2}\underset{(0,T)\times \Omega }{\int \int }|u-u_d{|}^{2}dxdt+\frac{{\alpha }_2}{2}{\int }_{\Omega }{|u\left(T\right)-{\tilde{u}}_d|}^{2}dx}\hfill \\ {\displaystyle +\frac{{\alpha }_3}{2}\underset{(0,T)\times {\Gamma }_{in}}{\int \int }{\left|g\right|}^{2}dxdt+\frac{{\alpha }_4}{2}\underset{(0,T)\times {\Gamma }_{in}}{\int \int }{\left|{\nabla }_{{\Gamma }_{in}}g\right|}^{2}dxdt,}\hfill \end{array}$$

(5)

where ${\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4$ are positive given constants and $u_d\in L^2((0,T)\times \Omega )$, ${\tilde{u}}_d\in L^2(\Omega )$ are given functions and ${\nabla }_{{\Gamma }_{in}}g$ is the tangential gradient in ${\Gamma }_{in}$.

Consider also the admissible control set:

$${\mathcal{G}}_{ad}:=\{g\in \mathcal{G},g_{ext}\left(0\right)\in H\}\subset L^2(0,T;H^{1/2}(\partial \Omega ))$$

and the feasible set

$${\mathcal{F}}_{ad}:=\{(g,u)\in {\mathcal{G}}_{ad}\times \mathbb{W}(0,T):uistheweaksolutionof\left(\right)associatedtog\}.$$

The problem that we intend to solve in this section is the following:

Find $(\widehat{g},\widehat{u})\in {\mathcal{F}}_{ad}$ such that

$$J(\widehat{g},\widehat{u})\le J(g,u),\forall (g,u)\in {\mathcal{F}}_{ad}$$

(6)

The existence of an optimal control $(\widehat{g},\widehat{u})$ can then be ensured as follows.

Theorem 2. 

Let $\Omega \subset {\mathbb{R}}^d$, $d=2,3$, be an open boundary set with a locally Lipschitz boundary, $h\in L^2((0,T)\times {\Gamma }_{out})$, $f\in L^2((0,T)\times \Omega )$ and $u_0\in \mathbb{H}$. Then, there exists a pair $(\widehat{g},\widehat{u})\in {\mathcal{F}}_{ad}$ solution of (6).

Proof:

Firstly, notice that if we consider for instance $g\in {\mathcal{G}}_{ad}$ to be null, thanks to Theorem 1, we can immediately see that the feasible set ${\mathcal{F}}_{ad}$ is non-empty. Since $J(g,u)\ge 0$ for all $(g,u)\in {\mathcal{F}}_{ad}$, which is non-empty, there exists a minimizing sequence $(g_n,u_n)\subset {\mathcal{F}}_{ad}$ such that

$$\underset{n\to \infty }{lim}J(g_n,u_n)=\overline{J}=\underset{(g,u)\in {\mathcal{F}}_{ad}}{inf}J(g,u).$$

Also,

$$\underset{n\to \infty }{lim}J(g_n,u_n)\ge {\parallel g_n\parallel }_{L^2(0,T,H_0^1\left({\Gamma }_{in}\right))}.$$

Therefore, $g_n$ is uniformly bounded in $L^2(0,T;H_0^1\left({\Gamma }_{in}\right))$.

Besides, there exists ${\left(g_{ext}\right)}_n\subset E_T\left(L^2(0,T;H_0^1\left({\Gamma }_{in}\right))\right)\cap \mathbb{W}(0,T)$ with $E_T\left(g_n\right)={\left(g_{ext}\right)}_n$ the bounded extension operator defined in Lemma 3. Since $g_n⟶gweaklyin{L}^2((0,T)\times {\Gamma }_{in})$ and ${\nabla }_{{\Gamma }_{in}}{g}_n⟶{\nabla }_{{\Gamma }_{in}}gweaklyin{L}^2((0,T)\times {\Gamma }_{in})$, then ${\left(g_{ext}\right)}_n$ in uniformly bounded in $L^2(0,T,\mathbb{V})$.

Taking into account that $u_n$ are the corresponding solutions of the weak Stokes equation, by using the estimates, we conclude that $u_n$ is an uniformly bounded sequence in $L^2(0,T;\mathbb{V})\cap C^0([0,T];\mathbb{H}).$

Hence, we can extract at least a subsequence, still renamed as $(g_n,u_n)$, which weakly converges to a certain pair $(g,u)$.

Moreover, there exists a pair $(g,u)$ such that, at least for a subsequence,

$$\begin{array}{c}{g}_n⟶gweaklyin{L}^2((0,T)\times {\Gamma }_{in}),\hfill \\ u_n⟶uweaklyin{L}^2(0,T;\mathbb{V}),\hfill \\ u_n⟶uweakly-*in{L}^{\infty }(0,T;\mathbb{H}),\hfill \\ u_n\left(0\right)\to u_{0}in\mathbb{H},\hfill \\ u_n\left(T\right)\to u\left(T\right)in\mathbb{H}.\hfill \end{array}$$

Inserting $\left(u_n\right)$ into (2) we obtain:

$$\left\{\begin{array}{c}{\displaystyle {\langle {\left(u_n\right)}_t,v\rangle }_{L^2\left(V^{*}\right),L^2\left(V\right)}+2\nu (D\left(u_n\right),D\left(v\right))=(f,v)+{(h,v)}_{L^2(0,T;L^2\left({\Gamma }_{out}\right))}\forall v\in L^2\left(V\right),}\hfill \\ u_n=g_n\mathrm{in}{L}^2((0,T)\times {\Gamma }_{in}),\hfill \\ u_n\left(0\right)=u_0\mathrm{in}\mathbb{H}.\hfill \end{array}\right.$$

We can, therefore, conclude that $(g,u)$ belongs to ${\mathcal{F}}_{ad}$. Finally, as J is coercive and lower semicontinuous with respect to the weak convergence, we conclude that

$$\overline{J}=\underset{n\to \infty }{lim}J(g_n,u_n)\ge J(g,u)\ge \overline{J}$$

and $(g,u)$ is a minimizer.

□

3.2. Differentiability and Characterisation Results

Differentiability is an important issue in the study of optimal control problems, as it allows the derivation of the first order optimality conditions which provide further characterisation of the optimal pair.

With this aim, we will follow an approach that is inspired in [21], for the case of a single boundary condition.

We start by introducing two auxiliary systems which are related, respectively, to the linearized and the adjoint state equations which will later be used for the optimality conditions:

$$\left\{\begin{array}{cc}{\displaystyle {\left(w_1\right)}_t-\nabla ·\left(2\nu D\left(w_1\right)\right)+\nabla q_1=f_1}\hfill & \Omega \times (0,T),\hfill \\ \nabla ·w_1=0\hfill & \Omega \times (0,T),\hfill \\ w_1=g_1\hfill & {\Gamma }_{in}\times (0,T),\hfill \\ w_1=0\hfill & {\Gamma }_{\omega }\times (0,T)\hfill \\ \left(2\nu D\left(w_1\right)-q_{1}I\right)·\overline{n}=0\hfill & {\Gamma }_{out}\times (0,T),\hfill \\ w_1(0,x)=0\hfill & \Omega \hfill \end{array}\right.$$

(7)

and

$$\left\{\begin{array}{cc}{\displaystyle -{\left(w_2\right)}_t-\nabla ·\left(2\nu D\left(w_2\right)\right)+\nabla q_2=f_2}\hfill & \Omega \times (0,T),\hfill \\ \nabla ·w_2=0\hfill & \Omega \times (0,T),\hfill \\ w_2=0\hfill & {\Gamma }_{in}\times (0,T),\hfill \\ w_2=0\hfill & {\Gamma }_{\omega }\times (0,T)\hfill \\ \left(2\nu D\left(w_2\right)-q_{2}I\right)·\overline{n}=0\hfill & {\Gamma }_{out}\times (0,T),\hfill \\ w_2(T,x)=d_2\hfill & \Omega .\hfill \end{array}\right.$$

(8)

We can give the following regularity result:

Proposition 1. 

Let $u\in \mathbb{W}$, $f_1\in L^2(0,T;V^{*})$, $f_2\in L^2(0,T;V^{*})$, $g_1\in G_{ad}$ and $d_2\in H.$ Then, the linear system (7) admits a unique weak solution $w_1\in \mathbb{W}(0,T)$ and the system (8) has a unique weak solution $w_2\in L^2(0,T;V)$.

The weak formulations of (7) and (8), although associated to linear problems, include transport type terms that can be recast into the frame of Theorem 1.

Let us now write problem (6) as:

$$\begin{array}{c}minJ(u_L+w,(u_L+w)\left(T\right),g)\hfill \\ \mathrm{subject} \mathrm{to}:e(u_L,w,g)=0\hfill \end{array}$$

where $e:\mathbb{W}(0,T)\times {\mathcal{G}}_{ad}\mapsto L^2\left(V^{*}\right)\times \mathbb{H}$.

In this case, $e=(e_1,e_2)$ with

$$\begin{array}{c}{\displaystyle \langle e_1(\tilde{u},g_{ext}),v\rangle ={\langle {\tilde{u}}_t,v\rangle }_{L^2\left(V^{*}\right),L^2\left(V\right)}+2\nu (D\left(\tilde{u}\right),D\left(v\right))-(f,v)}\hfill \\ -{(h,v)}_{L^2(0,T;L^2\left({\Gamma }_{out}\right))}+{\langle {\left(g_{ext}\right)}_t,v\rangle }_{L^2\left(V^{*}\right),L^2\left(V\right)}+2\nu (D\left(g_{ext}\right),D\left(v\right))\hfill \end{array}$$

$$e_2(\tilde{u},g_{ext})=\tilde{u}\left(0\right)-u_0+g_{ext}\left(0\right),$$

with $u=\tilde{u}+g_{ext}$ and $g=t{r}_{(0,T)\times {\Gamma }_{in}}{g}_{ext}.$

Proposition 2. 

The mapping $e:=(e_1,e_2)$ is twice continuously Fréchet differentiable with local Lipschitz continuous second derivative.

Note that, the differentiability of the linear operators $e_1,e_2$ is immediate. So, we can then conclude that the application $e_1$ and $e_2$ are continuous and twice continuously Fréchet differentiable, and the second derivative is Lipschitz continuous.

We are now in conditions to provide a necessary condition that characterizes the optimal solution, known to exist thanks to Theorem (2).

In the sequel, we denote by $\left(\mathcal{A}\right)$ the list of assumptions:

$$\left(\mathcal{A}\right):=\left\{\begin{array}{c}\Omega \subset {\mathbb{R}}^d,d=2,3,\mathrm{be} \mathrm{an} \mathrm{open} \mathrm{boundary} \mathrm{set} \mathrm{of} \mathrm{class} {\mathcal{C}}^1,\hfill \\ g\in \mathcal{G} \mathrm{with} g=t{r}_{(0,T)\times {\Gamma }_{in}}\left(g_{ext}\right),\hfill \\ h\in L^2((0,T)\times {\Gamma }_{out}),\hfill \\ f\in L^2((0,T)\times \Omega )\hfill \\ u_0\in \mathbb{H}with{u}_0-\tilde{g}\left(0\right)\in H.\hfill \end{array}\right.$$

Theorem 3. 

Let $\left(\mathcal{A}\right)$ be satisfied and $(\widehat{g},\widehat{u})$ a minimum provided by Theorem 2, then there exists an adjoint state $\widehat{\phi }$ and an adjoint pressure $\widehat{q}$ such that $(\widehat{g},\widehat{u},\widehat{\phi },\widehat{q})$ verify:

$$\left\{\begin{array}{c}\widehat{u}\in \mathbb{W}(0,T),\hfill \\ {\displaystyle {\langle {\widehat{u}}_t,v\rangle }_{L^2\left(V^{*}\right),L^2\left(V\right)}+2\nu (D\left(\widehat{u}\right),D\left(v\right))=}\hfill \\ (f,v)+{(h,v)}_{L^2(0,T;L^2\left({\Gamma }_{out}\right))}\forall v\in L^2\left(V\right),\hfill \\ \widehat{u}=\widehat{g}\mathrm{in}{L}^2((0,T)\times {\Gamma }_{in}),\hfill \\ \widehat{u}\left(0\right)=u_0\mathrm{in}\mathbb{H},\hfill \end{array}\right.$$

(9a)

and

$$\left\{\begin{array}{c}\widehat{\phi }\in W(0,T),\hfill \\ {\displaystyle -{\langle {\widehat{\phi }}_t,v\rangle }_{L^2\left(V^{*}\right),L^2\left(V\right)}+2\nu (D\left(\widehat{\phi }\right),D\left(v\right))}\hfill \\ {\displaystyle ={\alpha }_1\underset{(0,T)\times \Omega }{\int \int }(\widehat{u}-u_d)vdxdt\forall v\in L^2\left(V\right),}\hfill \\ \widehat{\phi }\left(T\right)={\alpha }_2(\widehat{u}\left(T\right)-{\tilde{u}}_d),\hfill \end{array}\right.$$

(9b)

and

$${\alpha }_3\left(\widehat{g},\tilde{g}\right)+{\alpha }_4\left({\nabla }_{{\Gamma }_{in}}\widehat{g},{\nabla }_{{\Gamma }_{in}}\tilde{g}\right)=-\left(\left((-2\nu D\left(\widehat{\phi }\right)+\widehat{q}I)\overline{n}\right),\tilde{g}\right)\forall \tilde{g}\in L^2(0,T;H_0^1\left({\Gamma }_{in}\right))$$

(9c)

Proof:

Note that, as we have a unique solution each control g, we can identify the reduced cost functional, i.e, we can write that

$${\displaystyle J\left(g\right)=J(g,u)=\frac{{\alpha }_1}{2}\underset{(0,T)\times \Omega }{\int \int }|u-u_d{|}^{2}dxdt+\frac{{\alpha }_2}{2}{\int }_{\Omega }{|u\left(T\right)-{\tilde{u}}_d|}^{2}dx}$$

$$+{\displaystyle \frac{{\alpha }_3}{2}}\underset{(0,T)\times {\Gamma }_{in}}{\int \int }{\left|g\right|}^{2}dxdt+{\displaystyle \frac{{\alpha }_4}{2}}\underset{(0,T)\times {\Gamma }_{in}}{\int \int }{\left|{\nabla }_{{\Gamma }_{in}}g\right|}^{2}dxdt.$$

Remark that the above reasoning is possible because we have a uniqueness result for the Stokes equations. In case of not having a uniqueness of solution, as in the 3D case of Navier-Stokes equations, this classical reasoning could not be used and we would have to resort to more advanced techniques to the proof. Specifically, we could try to follow the reasoning already used in other works, such as [5,6], in which a uniqueness of solution was not assumed.

It is easy to check that J is Gâteaux differentiable, $J^{\prime }\left(g\right):{\mathcal{G}}_{ad}\mapsto L^2((0,T)\times {\Gamma }_{in})$ and

$${\displaystyle {〈J^{\prime }\left(g\right),\tilde{g}〉}_{L^2(0,T;H^{1/2}(\partial \Omega ))}=\underset{\epsilon \to 0}{lim}\frac{J(g+\epsilon \tilde{g})-J\left(g\right)}{\epsilon }}$$

$$={\alpha }_1\underset{(0,T)\times \Omega }{\int \int }(u-u_d)vdxdt+{\alpha }_2{\int }_{\Omega }(u\left(T\right)-{\tilde{u}}_d)v\left(T\right)dx$$

$$+{\alpha }_3\underset{(0,T)\times {\Gamma }_{in}}{\int \int }g\tilde{g}dxdt+{\alpha }_4\underset{(0,T)\times {\Gamma }_{in}}{\int \int }{\nabla }_{{\Gamma }_{in}}g{\nabla }_{{\Gamma }_{in}}\tilde{g}dxdt,$$

with v solution of:

$$\left\{\begin{array}{c}v\in W(0,T),\hfill \\ {\displaystyle {\langle v_t,v_1\rangle }_{L^2\left(V^{*}\right),L^2\left(V\right)}+2\nu (D\left(v\right),D\left(v_1\right))=0\forall v_1\in L^2\left(V\right),}\hfill \\ v=\tilde{g}\mathrm{in}{L}^2((0,T)\times {\Gamma }_{in}),\hfill \\ v\left(0\right)=0.\hfill \end{array}\right.$$

Therefore, integrating by parts we obtain:

$${〈J^{\prime }\left(g\right),\tilde{g}〉}_{L^2(0,T;H_0^1\left({\Gamma }_{in}\right))}=\underset{(0,T)\times {\Gamma }_{in}}{\int \int }(((-2\nu D\left(\phi \right)+qI)\overline{n}+{\alpha }_{3}g)\tilde{g}+{\alpha }_4{\nabla }_{{\Gamma }_{in}}g{\nabla }_{{\Gamma }_{in}}\tilde{g})dxdt$$

for all $\tilde{g}\in L^2(0,T;H_0^1\left({\Gamma }_{in}\right))$ with $\phi$ the adjoint state associated (9b) to the weak solution to Stokes equations (9a) with $\widehat{g}=g.$

Note that our problem satisfies the conditions of Theorem 9.6 of [24]:

• J and e are Fréchet differentiable, thanks to Proposition 2.
• For each $g\in {\mathcal{G}}_{ad}$, the state equation $e(u,g)=0$ defines a unique control-to-state map $g\mapsto u\left(g\right)$, that is $e\left(u\right(g),g)=0$ for any $g\in {\mathcal{G}}_{ad}$, thanks to the uniqueness result proved in Theorem 1. Moreover, the control-to-state map is Gateaux differentiable, thanks to Proposition 2.
• The partial derivative $e_u^{*}(u\left(g\right),g)$ is a continuous isomorphism. In fact, thanks to Proposition 2, $e_u(u\left(g\right),g):W(0,T)\mapsto L^2\left(V^{*}\right)\times \mathbb{H}$ where $e_{u}v=({\left(e_1\right)}_{u}v,{\left(e_2\right)}_{u}v)$ with

  $$\begin{array}{c}{\displaystyle \langle {\left(e_1\right)}_{u}v,v_1\rangle ={\langle {\tilde{v}}_t,v_1\rangle }_{L^2\left(V^{*}\right),L^2\left(V\right)}+2\nu (D\left(\tilde{v}\right),D\left(v_1\right))}\hfill \\ {\langle {\left(g_{ext}\right)}_t,v_1\rangle }_{L^2\left(V^{*}\right),L^2\left(V\right)}+2\nu (D\left(g_{ext}\right),D\left(v_1\right)),\hfill \end{array}$$

  $${\left(e_2\right)}_{u}v=\tilde{v}\left(0\right)+g_{ext}\left(0\right),$$

  with $v=\tilde{v}+g_{ext}$ and $g=t{r}_{(0,T)\times {\Gamma }_{in}}{g}_{ext}.$
  As a consequence of Proposition 1, the mapping $e_u(u\left(g\right),g):W(0,T)\mapsto L^2\left(V^{*}\right)\times \mathbb{H}$ is a homeomorphism. Hence by Proposition 2 and the Implicit Function Theorem the first derivative $u^{\prime }\left(g\right)\tilde{g}$ at $\tilde{g}\in \mathcal{G}$ can be expressed as

  $$u^{\prime }\left(g\right)\tilde{g}=-e_u^{-1}(u\left(g\right),g)e_g(u\left(g\right),g)\tilde{g}\in W(0,T).$$

Hence, in view of that result, we can ensure that, if $\widehat{g}$ is a minimum control provided by Theorem 2 associated to the state u, then there exists an adjoint pair of state and pressure variables satisfying the optimality system.

□

If we want to consider strong solution, then the optimality system is:

$$\left\{\begin{array}{cc}{\displaystyle {\partial }_t\widehat{u}-\nabla ·\left(2\nu D\left(\widehat{u}\right)\right)+\nabla \stackrel{++}{=}=f}\hfill & (0,T)\times \Omega ,\hfill \\ {\displaystyle \nabla ·\widehat{u}=0}\hfill & (0,T)\times \Omega ,\hfill \\ {\displaystyle \widehat{u}=\widehat{g}}\hfill & (0,T)\times {\Gamma }_{in},\hfill \\ {\displaystyle \widehat{u}=0}\hfill & (0,T)\times {\Gamma }_{\omega },\hfill \\ {\displaystyle \left(2\nu D\left(\widehat{u}\right)-\stackrel{++}{=}I\right)·\overline{n}=h}\hfill & (0,T)\times {\Gamma }_{out},\hfill \\ {\displaystyle \widehat{u}(0,x)=u_0\left(x\right)}\hfill & \Omega \hfill \end{array}\right.$$

(10a)

and

$$\left\{\begin{array}{cc}{\displaystyle -{\partial }_t\widehat{\phi }-\nabla ·\left(2\nu D\left(\widehat{\phi }\right)\right)+\nabla \widehat{q}={\alpha }_1(\widehat{u}-u_d)}\hfill & (0,T)\times \Omega ,\hfill \\ {\displaystyle \nabla ·\widehat{\phi }=0}\hfill & (0,T)\times \Omega ,\hfill \\ {\displaystyle \widehat{\phi }=0}\hfill & (0,T)\times {\Gamma }_{in},\hfill \\ {\displaystyle \widehat{\phi }=0}\hfill & (0,T)\times {\Gamma }_{\omega },\hfill \\ {\displaystyle \left(2\nu D\left(\widehat{\phi }\right)-\widehat{q}I\right)·\overline{n}=0}\hfill & (0,T)\times {\Gamma }_{out},\hfill \\ {\displaystyle \widehat{\phi }(T,x)={\alpha }_2(\widehat{u}\left(T\right)-{\tilde{u}}_d)}\hfill & \Omega \hfill \end{array}\right.$$

(10b)

and

$${\alpha }_3\widehat{g}-{\alpha }_4{\nabla }_{{\Gamma }_{in}}\left({\nabla }_{{\Gamma }_{in}}\widehat{g}\right)={\left.-\left((-2\nu D\left(\widehat{\phi }\right)+\widehat{q}I)\overline{n}\right)\right|}_{(0,T)\times {\Gamma }_{in}.}$$

(10c)

These results allow us to conclude that the boundary control problem has a unique solution and we can compute it by using the optimallity system. From an application perspective, this means that for any desired data, it is possible to control the blood flow by adjusting only the fluid inlet profile, thereby directing the flow to achieve a desired and optimal movement for proper functioning.

4. Algorithm and Numerical Simulations

In the previous sections, we have conducted a theoretical study that has allowed us to characterize the optimal control. In particular, we have obtained an adjoint variable that enables us to calculate the gradient of the reduced cost. The gradient calculation is the main ingredient of the so-called descent methods. Thus, classic descent algorithms such as the gradient method will be proposed.

ALG: Optimal Step Gradient Method

(a) 

Choose $g^0\in L^2((0,T)\times {\Gamma }_{in})$ and $u^0\in L^2(0,T;\mathbb{V}).$

(b) 

Then, for given $n\ge 0$ and $g^n\in L^2((0,T)\times {\Gamma }_{in})$, compute the solution $(u^{n+1},p^{n+1})$ to (11)

$$\left\{\begin{array}{cc}{\displaystyle {\partial }_t{u}^{n+1}-\nabla ·\left(2\nu D\left(u^{n+1}\right)\right)+\nabla p^{n+1}=f^{n+1}}\hfill & (0,T)\times \Omega ,\hfill \\ \nabla ·u^{n+1}=0\hfill & (0,T)\times \Omega ,\hfill \\ u^{n+1}=g^n\hfill & (0,T)\times {\Gamma }_{in},\hfill \\ u^{n+1}=0\hfill & (0,T)\times {\Gamma }_{\omega },\hfill \\ \left(2\nu D\left(u^{n+1}\right)-p^{n+1}I\right)·\overline{n}=h^{n+1}\hfill & (0,T)\times {\Gamma }_{out},\hfill \\ u^{n+1}(0,x)=u_0\left(x\right)\hfill & \Omega \hfill \end{array}\right.$$

(11)

and the solution $({\phi }^{n+1},q^{n+1})$ to (12)

$$\left\{\begin{array}{cc}{\displaystyle -{\partial }_t{\phi }^{n+1}-\nabla ·\left(2\nu D\left({\phi }^{n+1}\right)\right)+\nabla q^{n+1}={\alpha }_1(u^{n+1}-u_d)}\hfill & (0,T)\times \Omega ,\hfill \\ \nabla ·{\phi }^{n+1}=0\hfill & (0,T)\times \Omega ,\hfill \\ {\phi }^{n+1}=0\hfill & (0,T)\times {\Gamma }_{in},\hfill \\ {\phi }^{n+1}=0\hfill & (0,T)\times {\Gamma }_{\omega },\hfill \\ \left(2\nu D\left({\phi }^{n+1}\right)-q^{n+1}I\right)·\overline{n}=0\hfill & (0,T)\times {\Gamma }_{out},\hfill \\ {\phi }^{n+1}(T,x)={\alpha }_2(u^{n+1}\left(T\right)-{\tilde{u}}_d)\hfill & \Omega \hfill \end{array}\right.$$

(12)

and set

$$g^{n+1}=g^n-{\rho }^n{d}^{n+1},$$

(13)

where

$$d^{n+1}={\left.(-2\nu D\left({\phi }^{n+1}\right)+q^{n+1}I)\right|}_{(0,T)\times {\Gamma }_{in}}+{\alpha }_3{g}^n-{\alpha }_4{\partial }_{yy}^2{g}^n$$

(14)

and

$${\displaystyle {\rho }^n=arg\left(\underset{\rho >0}{min}J(g^n-\rho d^{n+1})\right).}$$

(15)

Note that, we are considering that ${\nabla }_{{\Gamma }_{in}}$ is ${\partial }_y$ due to our ${\Gamma }_{in}$ is a vertical line.

Let us now provide further details on the numerical approximation used to solve (11) and (12). We start by introducing the variational form of the equations and then, we discretize our systems in time, by using the backward implicit Euler finite difference. Finally, we use the finite element method to discretize in space. Consequently, the discrete problems solved at each iteration are:

$$\left\{\begin{array}{c}{\displaystyle \left[\frac{u_h^{k+1}-u_h^k}{dt},v\right]+4\nu \left[D\left(u_h^{k+1}\right),D\left(v\right)\right]-\left[p_h^{k+1},\nabla ·v\right]=\left[f^{k+1},v\right]+{\int }_{{\Gamma }_{out}}{h}^{k+1}vdx}\hfill \\ \left[\nabla ·u_h^{k+1},pp\right]=0\hfill \\ u_h^{k+1}=g^n(t^{k+1},x_h)x_h\in {\Gamma }_{in},\hfill \\ u_h^{k+1}=0x_h\in {\Gamma }_{\omega }\hfill \\ u_h^0=u_0\left(x_h\right)x_h\in \Omega ,\hfill \end{array}\right.$$

with $v,pp$ the corresponding test functions and $u_h^{k+1}=u(t^{k+1},x_h),p_h^{k+1}=p(t^{k+1},x_h)$ (being u, p the $n+1$ iterate of the algorithm) for $k=0,1,\cdots ,N-1$. Similarly, we have:

$$\left\{\begin{array}{c}{\displaystyle -\left[\frac{{\phi }_h^{k+1}-{\phi }_h^k}{dt},v\right]+2\nu \left[D\left({\phi }_h^k\right),D\left(v\right)\right]-\left[q_h^k,\nabla ·v\right]={\alpha }_1\left[u_h^k-u_{id}(t^k,x_h),v\right]}\hfill \\ \left[\nabla ·{\phi }_h^k,pp\right]=0\hfill \\ {\phi }_h^k=0x_h\in \partial \Omega \hfill \\ {\phi }_h^N={\alpha }_2(u_h^{k+1}-{\tilde{u}}_d\left(x_h\right))x_h\in \Omega .\hfill \end{array}\right.$$

On the other hand, in the numerical experiment, we use an explicit approximation formula to compute the optimal step $\rho$. For a given search direction d, we introduce the linearization of the mapping $\rho \mapsto u(g-\rho d)=u\left(g\right)-\rho u^{\prime }\left(g\right)d.$ And the optimal step is given by the expression $I^{\prime }\left(\rho \right)=0$ with $I\left(\rho \right)=J(u\left(g\right)-\rho u^{\prime }\left(g\right)d,g-\rho d).$ So the unique root of the $I^{\prime }\left(\rho \right)$ is given by:

$${\displaystyle {\rho }^{*}:=\frac{(J^{\prime }\left(g\right),d)}{{\displaystyle {\alpha }_1\underset{(0,T)\times \Omega }{\int \int }{\left|v\right|}^2+{\alpha }_2{\int }_{\Omega }{\left|v\left(T\right)\right|}^2+{\alpha }_3\underset{(0,T)\times {\Gamma }_{in}}{\int \int }{\left|d\right|}^2+{\alpha }_4\underset{(0,T)\times {\Gamma }_{in}}{\int \int }{\left|{\partial }_{y}d\right|}^2}},}$$

where v is the solution to the homogeneous Stokes system with $v\left(0\right)=0.$

• Test: 2D stenotic vessel capillary with noise data

Our test addressed an idealized domain representing a partially obstructed vessel capillary, and we aimed to approximate the solution to a noisy solution in order to simulate the noise obtained from medical tests

The domain $\Omega$ is represented in Figure 2. In this case, the left side is the inlet boundary (where we imposed the control) while the right side corresponds to outlet boundary (where we imposed the total stress condition). The remaining boundaries represent the vessel wall where we imposed homogeneous Dirichlet conditions.

In order to solve numerically the systems appearing in ALG, we had to fix a mesh, the finite element approximation spaces, and a time step. For the space discretization, we used a mixed finite element formulation with continuous piece-wise P2 and P1 functions for the velocity field and the pressure, respectively. For details, see [12]. The simulations for the current test, as well as for the following ones, have been performed with the FreeFem++ package (see [19]).

We take $T=1$, $dt=0.1$, $f=0$, considering stationary total stress $h=(2x-40,10-2y)$. We want to drive the solution close to desired states, $u_d(t,x)$ and ${\tilde{u}}_d\left(x\right)=u_d(x,T)$, where $u_d$ is a noisy perturbation of the solution $U_d$ of (16), and $T=1$:

$$\left\{\begin{array}{cc}{\displaystyle {\partial }_{t}u-\nabla ·\left(2\nu D\left(u\right)\right)+\nabla p=0}\hfill & (0,T)\times \Omega ,\hfill \\ \nabla ·u=0\hfill & (0,T)\times \Omega ,\hfill \\ u=(1.5t·y·(2-y),0)\hfill & (0,T)\times {\Gamma }_{in},\hfill \\ u=0\hfill & (0,T)\times {\Gamma }_{\omega },\hfill \\ \left(2\nu D\left(u\right)-pI\right)·\overline{n}=(0,0)\hfill & (0,T)\times {\Gamma }_{out},\hfill \\ u(0,x)=(0,0)\hfill & \Omega \hfill \end{array}\right.$$

(16)

We considered, as intial guess for the control,

$${\displaystyle g^0=1.5sin\left(\frac{\pi }{2}t\right)·y·(2-y).}$$

We will now present relative error for different experiments in which we have varied the Reynolds number and the noise rate of the data. It is important to note that in this case, we have chosen small Reynolds numbers (on the order of ${10}^{-2}$ and ${10}^{-3}$) because these correspond to the Reynolds number of blood in capillaries, which is where these equations can be applied. In the following table, see Figure 3, the relative errors between the calculated solution and the desired noiseless solution will be shown:

$${\displaystyle \frac{\parallel u-U_d{\parallel }_{L^2((0,T)\times \Omega )}}{\parallel U_d{\parallel }_{L^2((0,T)\times \Omega )}}.}$$

As can be observed in Figure 3, the algorithm performs quite well and successfully approximates the solution even when the noise levels in the data exceed $50\%.$ It is also notable that, as expected, when we decrease the value of ${\alpha }_1$, the relative errors increase, although they remain small. All of this suggests that the algorithm works efficiently and quickly converges to a solution that is very close to the desired one.

Next, we will present the solution obtained for a specific case at three different times throughout the interval $[0,1]$, see Figure 4, Figure 5 and Figure 6.

As can be seen in the images, see Figure 4, Figure 5 and Figure 6, the solution obtained with the algorithm is very close to the desired noiseless solution, even when we have high levels of noise in the data, as is the case here. Specifically, it should be noted that the noise is introduced into the Stokes equation due to the adjoint state, but as we can see, the algorithm is able to suppress that noise and find the noiseless solution. This aspect is very important from the perspective of application and the use of medical data, as data obtained from medical tests often have a high noise rate. Therefore, it is a very positive aspect that the algorithm is capable of overcoming this problem.

5. Conclusions

In this work, we have provided necessary conditions for a boundary control problem associated with the time-dependent Stokes equations, under boundary conditions of mixed type. In a certain sense, it can be seen as an extension of the work performed in [13], to the the time-dependent case under mixed boundary conditions. We could see that an initial guess, associated with a noise can be improved until decreasing the relative fitting error to less than $5\%.$ This is the case, even when high noise rate are considered, as well as nonphysical data.

One of the major contributions of this work is to lay the groundwork for extending it to the time dependent case. It is important to note that this extension is by no means trivial, and we are currently working on it. One of the key aspects to consider, if we want to extend the work to the Navier-Stokes equations in the two-dimensional case us that we don’t have a linarity equation so we have to use another technique to get the results. Moreover, in the three-dimensional case, we do not have uniqueness of solution nor good regularity even with Dirichlet-type conditions. Therefore, in the three-dimensional case, it will be necessary to work under the assumption that a unique solution exists and assuming sufficient regularity to dprove similar results to those presented in this work. Additionally, since gradient-type algorithms are known to become computationally expensive, parameterized approaches are now being considered for 3D simulations. Future extensions may consider with Navier-Stokes equations, as well as multiple boundary controls.