# Bilinear Optimal Control of an Advection-reaction-diffusion System

Roland Glowinski  Department of Mathematics, University of Houston, 44800
Calhoun Road, Houston, TX 77204, USA<EMAIL_ADDRESS>,  Yongcun Song
Department of Mathematics, The University of Hong Kong, Pok Fu Lam, Hong Kong,
PRC<EMAIL_ADDRESS>,  Xiaoming Yuan  Department of Mathematics, The
University of Hong Kong, Pok Fu Lam, Hong Kong, PRC<EMAIL_ADDRESS>and
Hangrui Yue  Department of Mathematics, The University of Hong Kong, Pok Fu
Lam, Hong Kong, PRC<EMAIL_ADDRESS>To the memory of J. L. Lions
(1928-2001) who suggested investigating problems like (BCP) below.

(Date: January 7, 2021)

###### Abstract.

We consider the bilinear optimal control of an advection-reaction-diffusion
system, where the control arises as the velocity field in the advection term.
Such a problem is generally challenging from both theoretical analysis and
algorithmic design perspectives mainly because the state variable depends
nonlinearly on the control variable and an additional divergence-free
constraint on the control is coupled together with the state equation.
Mathematically, the proof of the existence of optimal solutions is delicate,
and up to now, only some results are known for a few special cases where
additional restrictions are imposed on the space dimension and the regularity
of the control. We prove the existence of optimal controls and derive the
first-order optimality conditions in general settings without any extra
assumption. Computationally, the well-known conjugate gradient (CG) method can
be applied conceptually. However, due to the additional divergence-free
constraint on the control variable and the nonlinear relation between the
state and control variables, it is challenging to compute the gradient and the
optimal stepsize at each CG iteration, and thus nontrivial to implement the CG
method. To address these issues, we advocate a fast inner preconditioned CG
method to ensure the divergence-free constraint and an efficient inexactness
strategy to determine an appropriate stepsize. An easily implementable nested
CG method is thus proposed for solving such a complicated problem. For the
numerical discretization, we combine finite difference methods for the time
discretization and finite element methods for the space discretization.
Efficiency of the proposed nested CG method is promisingly validated by the
results of some preliminary numerical experiments.

###### Key words and phrases:

Bilinear optimal control, advection-reaction-diffusion system, conjugate
gradient method, nested iteration, finite element method, finite difference
method

###### 2020 Mathematics Subject Classification:

49M41,35Q93,49J20

The work of R.G is supported by the US Department of Energy (ORNL) and the
Hong Kong based Kennedy Wong Foundation. The work of X.Y is supported by the
General Research Fund 12302318 from the Hong Kong Research Grants Council.

##  1\. Introduction

###  1.1. Background and Motivation

The optimal control of distributed parameter systems has important
applications in various scientific areas, such as physics, chemistry,
engineering, medicine, and finance. We refer to, e.g. [12, 13, 14, 22, 26,
28], for a few references. In a typical mathematical model of a controlled
distributed parameter system, either boundary or internal locally distributed
controls are usually used; these controls have localized support and are
called additive controls because they arise in the model equations as additive
terms. Optimal control problems with additive controls have received a
significant attention in past decades following the pioneering work of J. L.
Lions [22], and many mathematical and computational tools have been developed,
see e.g., [12, 13, 14, 23, 27, 29]. However, it is worth noting that additive
controls describe the effect of external added sources or forces and they do
not change the principal intrinsic properties of the controlled system. Hence,
they are not suitable to deal with processes whose principal intrinsic
properties should be changed by some control actions. For instance, if we aim
at changing the reaction rate in some chain reaction-type processes from
biomedical, nuclear, and chemical applications, additive controls amount to
controlling the chain reaction by adding into or withdrawing out of a certain
amount of the reactants, which is not realistic. To address this issue, a
natural idea is to use certain catalysts or smart materials to control the
systems, which can be mathematically modeled by optimal control problems with
bilinear controls. We refer to [19] for more detailed discussions.

Bilinear controls, also known as multiplicative controls, enter the model as
coefficients of the corresponding partial differential equations (PDEs). These
bilinear controls can change some main physical characteristics of the system
under investigation, such as a natural frequency response of a beam or the
rate of a chemical reaction. In the literature, bilinear controls of
distributed parameter systems have become an increasingly popular topic and
bilinear optimal control problems constrained by various PDEs, such as
elliptic equations [20], convection-diffusion equations [3], parabolic
equations [18], the Schrödinger equation [16] and the Fokker-Planck equation
[7], have been widely studied both mathematically and computationally.

In particular, bilinear controls play a crucial role in optimal control
problems modeled by advection-reaction-diffusion systems. On one hand, the
control can be the coefficient of the diffusion or the reaction term. For
instance, a system controlled by the so-called catalysts that can accelerate
or slow down various chemical or biological reactions can be modeled by a
bilinear optimal control problem for an advection-reaction-diffusion equation
where the control arises as the coefficient of the reaction term [18]; this
kind of bilinear optimal control problems have been studied in e.g., [3, 4,
18, 19]. On the other hand, the systems can also be controlled by the velocity
field in the advection term, which captures important applications in e.g.,
bioremediation [15], environmental remediation process [21], and mixing
enhancement of different fluids [24]. We note that there is a very limited
research being done on the velocity field controlled bilinear optimal control
problems; and only some special one-dimensional space cases have been studied
in [15, 17, 21] for the existence of an optimal control and the derivation of
first-order optimality conditions. To the best of our knowledge, no work has
been done yet to develop efficient numerical methods for solving multi-
dimensional bilinear optimal control problems controlled by the velocity field
in the advection term. All these facts motivate us to study bilinear optimal
control problems constrained by an advection-reaction-diffusion equation,
where the control enters into the model as the velocity field in the advection
term. Actually, investigating this kind of problems was suggested to one of us
(R. Glowinski), in the late 1990’s, by J. L. Lions (1928-2001).

###  1.2. Model

Let $\Omega$ be a bounded domain of $\mathbb{R}^{d}$ with $d\geq 1$ and let
$\Gamma$ be its boundary. We consider the following bilinear optimal control
problem:

(BCP) $\displaystyle\left\\{\begin{aligned} &\bm{u}\in\mathcal{U},\\\
&J(\bm{u})\leq J(\bm{v}),\forall\bm{v}\in\mathcal{U},\end{aligned}\right.$

with the objective functional $J$ defined by

(1.1)
$J(\bm{v})=\frac{1}{2}\iint_{Q}|\bm{v}|^{2}dxdt+\frac{\alpha_{1}}{2}\iint_{Q}|y-y_{d}|^{2}dxdt+\frac{\alpha_{2}}{2}\int_{\Omega}|y(T)-y_{T}|^{2}dx,$

and $y=y(t;\bm{v})$ the solution of the following advection-reaction-diffusion
equation

(1.2) $\displaystyle\left\\{\begin{aligned} \frac{\partial y}{\partial
t}-\nu\nabla^{2}y+\bm{v}\cdot\nabla y+a_{0}y&=f\quad\text{in}\quad Q,\\\
y&=g\quad\text{on}\quad\Sigma,\\\ y(0)&=\phi.\end{aligned}\right.$

Above and below, $Q=\Omega\times(0,T)$ and $\Sigma=\Gamma\times(0,T)$ with
$0<T<+\infty$; $\alpha_{1}\geq 0,\alpha_{2}\geq 0,\alpha_{1}+\alpha_{2}>0$;
the target functions $y_{d}$ and $y_{T}$ are given in $L^{2}(Q)$ and
$L^{2}(\Omega)$, respectively; the diffusion coefficient $\nu>0$ and the
reaction coefficient $a_{0}$ are assumed to be constants; the functions $f\in
L^{2}(Q)$, $g\in L^{2}(0,T;H^{1/2}(\Gamma))$ and $\phi\in L^{2}(\Omega)$. The
set $\mathcal{U}$ of the admissible controls is defined by

$\mathcal{U}:=\\{\bm{v}|\bm{v}\in[L^{2}(Q)]^{d},\nabla\cdot\bm{v}=0\\}.$

Clearly, the control variable $\bm{v}$ arises in (BCP) as a flow velocity
field in the advection term of (1.2), and the divergence-free constraint
$\nabla\cdot\bm{v}=0$ implies that the flow is incompressible. One can control
the system by changing the flow velocity $\bm{v}$ in order that $y$ and $y(T)$
are good approximations to $y_{d}$ and $y_{T}$, respectively.

###  1.3. Difficulties and Goals

In this paper, we intend to study the bilinear optimal control problem (BCP)
in the general case of $d\geq 2$ both mathematically and computationally.
Precisely, we first study the well-posedness of (1.2), the existence of an
optimal control $\bm{u}$, and its first-order optimality condition. Then,
computationally, we propose an efficient and relatively easy to implement
numerical method to solve (BCP). For this purpose, we advocate combining a
conjugate gradient (CG) method with a finite difference method (for the time
discretization) and a finite element method (for the space discretization) for
the numerical solution of (BCP). Although these numerical approaches have been
well developed in the literature, it is nontrivial to implement them to solve
(BCP) as discussed below, due to the complicated problem settings.

####  1.3.1. Difficulties in Algorithmic Design

Conceptually, a CG method for solving (BCP) can be easily derived following
[14]. However, CG algorithms are challenging to implement numerically for the
following reasons: 1). The state $y$ depends non-linearly on the control
$\bm{v}$ despite the fact that the state equation (1.2) is linear. 2). The
additional divergence-free constraint on the control $\bm{v}$, i.e.,
$\nabla\cdot\bm{v}=0$, is coupled together with the state equation (1.2).

To be more precise, the fact that the state $y$ is a nonlinear function of the
control $\bm{v}$ makes the optimality system a nonlinear problem. Hence,
seeking a suitable stepsize in each CG iteration requires solving an
optimization problem and it can not be as easily computed as in the linear
case [14]. Note that commonly used line search strategies are too expensive to
employ in our settings because they require evaluating the objective
functional value $J(\bm{v})$ repeatedly and every evaluation of $J(\bm{v})$
entails solving the state equation (1.2). The same concern on the
computational cost also applies when the Newton method is employed to solve
the corresponding optimization problem for finding a stepsize. To tackle this
issue, we propose an efficient inexact stepsize strategy which requires
solving only one additional linear parabolic problem and is cheap to implement
as shown in Section 3.

Furthermore, due to the divergence-free constraint $\nabla\cdot\bm{v}=0$, an
extra projection onto the admissible set $\mathcal{U}$ is required to compute
the first-order differential of $J$ at each CG iteration in order that all
iterates of the CG method are feasible. Generally, this projection subproblem
has no closed-form solution and has to be solved iteratively. Here, we
introduce a Lagrange multiplier associated with the constraint
$\nabla\cdot\bm{v}=0$, then the computation of the first-order differential
$DJ(\bm{v})$ of $J$ at $\bm{v}$ is equivalent to solving a Stokes type
problem. Inspired by [9], we advocate employing a preconditioned CG method,
which operates on the space of the Lagrange multiplier, to solve the resulting
Stokes type problem. With an appropriately chosen preconditioner, a fast
convergence of the resulting preconditioned CG method can be expected in
practice (and indeed, has been observed).

####  1.3.2. Difficulties in Numerical Discretization

For the numerical discretization of (BCP), we note that if an implicit finite
difference scheme is used for the time discretization of the state equation
(1.2), a stationary advection-reaction-diffusion equation should be solved at
each time step. To solve this stationary advection-reaction-diffusion
equation, it is well known that standard finite element techniques may lead to
strongly oscillatory solutions unless the mesh-size is sufficiently small with
respect to the ratio between $\nu$ and $\|\bm{v}\|$. In the context of optimal
control problems, to overcome such difficulties, different stabilized finite
element methods have been proposed and analyzed, see e.g., [1, 6]. Different
from the above references, we implement the time discretization by a semi-
implicit finite difference method for simplicity, namely, we use explicit
advection and reaction terms and treat the diffusion term implicitly.
Consequently, only a simple linear elliptic equation is required to be solved
at each time step. We then implement the space discretization of the resulting
elliptic equation at each time step by a standard piecewise linear finite
element method and the resulting linear system is very easy to solve.

Moreover, we recall that the divergence-free constraint $\nabla\cdot\bm{v}=0$
leads to a projection subproblem, which is equivalent to a Stokes type
problem, at each iteration of the CG algorithm. As discussed in [8], to
discretize a Stokes type problem, direct applications of standard finite
element methods always lead to an ill-posed discrete problem. To overcome this
difficulty, one can use different types of element approximations for pressure
and velocity. Inspired by [8, 9], we employ the Bercovier-Pironneau finite
element pair [2] (also known as $P_{1}$-$P_{1}$ iso $P_{2}$ finite element) to
approximate the control $\bm{v}$ and the Lagrange multiplier associated with
the divergence-free constraint. More concretely, we approximate the Lagrange
multiplier by a piecewise linear finite element space which is twice coarser
than the one for the control $\bm{v}$. In this way, the discrete problem is
well-posed and can be solved by a preconditioned CG method. As a byproduct of
the above discretization, the total number of degrees of freedom of the
discrete Lagrange multiplier is only $\frac{1}{d2^{d}}$ of the number of the
discrete control. Hence, the inner preconditioned CG method is implemented in
a lower-dimensional space than that of the state equation (1.2), implying a
computational cost reduction. With the above mentioned discretization schemes,
we can relatively easily obtain the fully discrete version of (BCP) and derive
the discrete analogue of our proposed nested CG method.

###  1.4. Organization

An outline of this paper is as follows. In Section 2, we prove the existence
of optimal controls for (BCP) and derive the associated first-order optimality
conditions. An easily implementable nested CG method is proposed in Section 3
for solving (BCP) numerically. In Section 4, we discuss the numerical
discretization of (BCP) by finite difference and finite element methods. Some
preliminary numerical results are reported in Section 5 to validate the
efficiency of our proposed numerical approach. Finally, some conclusions are
drawn in Section 6.

##  2\. Existence of optimal controls and first-order optimality conditions

In this section, first we present some notation and known results from the
literature that will be used in later analysis. Then, we prove the existence
of optimal controls for (BCP) and derive the associated first-order optimality
conditions. Without loss of generality, we assume that $f=0$ and $g=0$ in
(1.2) for convenience.

###  2.1. Preliminaries

Throughout, we denote by $L^{s}(\Omega)$ and $H^{s}(\Omega)$ the usual Sobolev
spaces for any $s>0$. The space $H_{0}^{s}(\Omega)$ denotes the completion of
$C_{0}^{\infty}(\Omega)$ in $H^{s}(\Omega)$, where $C_{0}^{\infty}(\Omega)$
denotes the space of all infinitely differentiable functions over $\Omega$
with a compact support in $\Omega$. In addition, we shall also use the
following vector-valued function spaces:

$\displaystyle\bm{L}^{2}(\Omega):=[L^{2}(\Omega)]^{d},$
$\displaystyle\bm{L}_{div}^{2}(\Omega):=\\{\bm{v}\in\bm{L}^{2}(\Omega),\nabla\cdot\bm{v}=0~{}\text{in}~{}\Omega\\}.$

Let $X$ be a Banach space with a norm $\|\cdot\|_{X}$, then the space
$L^{2}(0,T;X)$ consists of all measurable functions $z:(0,T)\rightarrow X$
satisfying

$\|z\|_{L^{2}(0,T;X)}:=\left(\int_{0}^{T}\|z(t)\|_{X}^{2}dt\right)^{\frac{1}{2}}<+\infty.$

With the above notation, it is clear that the admissible set $\mathcal{U}$ can
be denoted as $\mathcal{U}:=L^{2}(0,T;\bm{L}_{div}^{2}(\Omega))$. Moreover,
the space $W(0,T)$ consists of all functions $z\in
L^{2}(0,T;H_{0}^{1}(\Omega))$ such that $\frac{\partial z}{\partial t}\in
L^{2}(0,T;H^{-1}(\Omega))$ exists in a weak sense, i.e.

$W(0,T):=\\{z|z\in L^{2}(0,T;H_{0}^{1}(\Omega)),\frac{\partial z}{\partial
t}\in L^{2}(0,T;H^{-1}(\Omega))\\},$

where $H^{-1}(\Omega)(=H_{0}^{1}(\Omega)^{\prime})$ is the dual space of
$H_{0}^{1}(\Omega)$.

Next, we summarize some known results for the advection-reaction-diffusion
equation (1.2) in the literature for the convenience of further analysis.

The variational formulation of the state equation (1.2) reads: find $y\in
W(0,T)$ such that $y(0)=\phi$ and $\forall z\in L^{2}(0,T;H_{0}^{1}(\Omega))$,

(2.1) $\int_{0}^{T}\left\langle\frac{\partial y}{\partial
t},z\right\rangle_{H^{-1}(\Omega),H_{0}^{1}(\Omega)}dt+\nu\iint_{Q}\nabla
y\cdot\nabla zdxdt+\iint_{Q}\bm{v}\cdot\nabla yzdxdt+a_{0}\iint_{Q}yzdxdt=0,$

where
$\left\langle\cdot,\cdot\right\rangle_{H^{-1}(\Omega),H_{0}^{1}(\Omega)}$
denotes the duality pairing between $H^{-1}(\Omega)$ and $H_{0}^{1}(\Omega)$.
The existence and uniqueness of the solution $y\in W(0,T)$ to problem (2.1)
can be proved by standard arguments relying on the Lax-Milgram theorem, we
refer to [22] for the details. Moreover, we can define the control-to-state
operator $S:\mathcal{U}\rightarrow W(0,T)$, which maps $\bm{v}$ to
$y=S(\bm{v})$. Then, the objective functional $J$ in (BCP) can be reformulated
as

$J(\bm{v})=\frac{1}{2}\iint_{Q}|\bm{v}|^{2}dxdt+\frac{\alpha_{1}}{2}\iint_{Q}|S(\bm{v})-y_{d}|^{2}dxdt+\frac{\alpha_{2}}{2}\int_{\Omega}|S(\bm{v})(T)-y_{T}|^{2}dx,$

and the nonlinearity of the solution operator $S$ implies that (BCP) is
nonconvex.

For the solution $y\in W(0,T)$, we have the following estimates.

######  Lemma 2.1.

Let $\bm{v}\in L^{2}(0,T;\bm{L}^{2}_{div}(\Omega))$, then the solution $y\in
W(0,T)$ of the state equation (1.2) satisfies the following estimate:

(2.2) $\|y(t)\|_{L^{2}(\Omega)}^{2}+2\nu\int_{0}^{t}\|\nabla
y(s)\|_{L^{2}(\Omega)}^{2}ds+2a_{0}\int_{0}^{t}\|y(s)\|_{L^{2}(\Omega)}^{2}ds=\|\phi\|_{L^{2}(\Omega)}^{2}.$

###### Proof.

We first multiply the state equation (1.2) by $y(t)$, then applying the
Green’s formula in space yields

(2.3) $\frac{1}{2}\frac{d}{dt}\|y(t)\|_{L^{2}(\Omega)}^{2}=-\nu\|\nabla
y(t)\|_{L^{2}(\Omega)}^{2}-a_{0}\|y(t)\|_{L^{2}(\Omega)}^{2}.$

The desired result (2.2) can be directly obtained by integrating (2.3) over
$[0,t]$. ∎

Above estimate implies that

(2.4) $y~{}\text{is bounded in}~{}L^{2}(0,T;H_{0}^{1}(\Omega)).$

On the other hand,

$\frac{\partial y}{\partial t}=\nu\nabla^{2}y-\bm{v}\cdot\nabla y-a_{0}y,$

and the right hand side is bounded in $L^{2}(0,T;H^{-1}(\Omega))$. Hence,

(2.5) $\frac{\partial y}{\partial t}~{}\text{is bounded
in}~{}L^{2}(0,T;H^{-1}(\Omega)).$

Furthermore, since $\nabla\cdot\bm{v}=0$, it is clear that

$\iint_{Q}\bm{v}\cdot\nabla yzdxdt=\iint_{Q}\nabla
y\cdot(\bm{v}z)dxdt=-\iint_{Q}y\nabla\cdot(\bm{v}z)dxdt=-\iint_{Q}y(\bm{v}\cdot\nabla
z)dxdt,\forall z\in L^{2}(0,T;H_{0}^{1}(\Omega)).$

Hence, the variational formulation (2.1) can be equivalently written as:“ find
$y\in W(0,T)$ such that $y(0)=\phi$ and $\forall z\in
L^{2}(0,T;H_{0}^{1}(\Omega))$,

$\int_{0}^{T}\left\langle\frac{\partial y}{\partial
t},z\right\rangle_{H^{-1}(\Omega),H_{0}^{1}(\Omega)}dt+\nu\iint_{Q}\nabla
y\cdot\nabla zdxdt-\iint_{Q}(\bm{v}\cdot\nabla
z)ydxdt+a_{0}\iint_{Q}yzdxdt=0.$

###  2.2. Existence of Optimal Controls

With above preparations, we prove in this subsection the existence of optimal
controls for (BCP). For this purpose, we first show that the objective
functional $J$ is weakly lower semi-continuous.

######  Lemma 2.2.

The objective functional $J$ given by (1.1) is weakly lower semi-continuous.
That is, if a sequence $\\{\bm{v}_{n}\\}$ converges weakly to $\bar{\bm{v}}$
in $L^{2}(0,T;\bm{L}^{2}_{div}(\Omega))$, we have

$J(\bar{\bm{v}})\leq\underset{n\rightarrow\infty}{\lim\inf}J(\bm{v}_{n}).$

###### Proof.

Let $\\{\bm{v}_{n}\\}$ be a sequence that converges weakly to $\bar{\bm{v}}$
in $L^{2}(0,T;\bm{L}^{2}_{div}(\Omega))$ and $y_{n}:=y(x,t;\bm{v}_{n})$ the
solution of the following variational problem: find $y_{n}\in W(0,T)$ such
that $y_{n}(0)=\phi$ and $\forall z\in L^{2}(0,T;H_{0}^{1}(\Omega))$,

(2.6) $\int_{0}^{T}\left\langle\frac{\partial y_{n}}{\partial
t},z\right\rangle_{H^{-1}(\Omega),H_{0}^{1}(\Omega)}dt+\nu\iint_{Q}\nabla
y_{n}\cdot\nabla zdxdt-\iint_{Q}(\bm{v}_{n}\cdot\nabla
z)y_{n}dxdt+a_{0}\iint_{Q}y_{n}zdxdt=0.$

Moreover, it follows from (2.4) and (2.5) that there exists a subsequence of
$\\{y_{n}\\}$, still denoted by $\\{y_{n}\\}$ for convenience, such that

$y_{n}\rightarrow\bar{y}~{}\text{weakly in}~{}L^{2}(0,T;H_{0}^{1}(\Omega)),$

and

$\frac{\partial y_{n}}{\partial t}\rightarrow\frac{\partial\bar{y}}{\partial
t}~{}\text{weakly in}~{}L^{2}(0,T;H^{-1}(\Omega)).$

Since $\Omega$ is bounded, it follows directly from the compactness property
(also known as Rellich’s Theorem) that

$y_{n}\rightarrow\bar{y}~{}\text{strongly in}~{}L^{2}(0,T;L^{2}(\Omega)).$

Taking $\bm{v}_{n}\rightarrow\bar{\bm{v}}$ weakly in
$L^{2}(0,T;\bm{L}_{div}^{2}(\Omega))$ into account, we can pass the limit in
(2.6) and derive that $\bar{y}(0)=\phi$ and $\forall z\in
L^{2}(0,T;H_{0}^{1}(\Omega))$,

$\int_{0}^{T}\left\langle\frac{\partial\bar{y}}{\partial
t},z\right\rangle_{H^{-1}(\Omega),H_{0}^{1}(\Omega)}dt+\nu\iint_{Q}\nabla\bar{y}\cdot\nabla
zdxdt-\iint_{Q}(\bar{\bm{v}}\cdot\nabla
z)\bar{y}dxdt+a_{0}\iint_{Q}\bar{y}zdxdt=0,$

which implies that $\bar{y}$ is the solution of the state equation (1.2)
associated with $\bar{\bm{v}}$.

Since any norm of a Banach space is weakly lower semi-continuous, we have that

$\displaystyle\underset{n\rightarrow\infty}{\lim\inf}J(\bm{v}_{n})$
$\displaystyle=$
$\displaystyle\underset{n\rightarrow\infty}{\lim\inf}\left(\frac{1}{2}\iint_{Q}|\bm{v}_{n}|^{2}dxdt+\frac{\alpha_{1}}{2}\iint_{Q}|y_{n}-y_{d}|^{2}dxdt+\frac{\alpha_{2}}{2}\int_{\Omega}|y_{n}(T)-y_{T}|^{2}dx\right)$
$\displaystyle\geq$
$\displaystyle\frac{1}{2}\iint_{Q}|\bar{\bm{v}}|^{2}dxdt+\frac{\alpha_{1}}{2}\iint_{Q}|\bar{y}-y_{d}|^{2}dxdt+\frac{\alpha_{2}}{2}\int_{\Omega}|\bar{y}(T)-y_{T}|^{2}dx$
$\displaystyle=$ $\displaystyle J(\bar{\bm{v}}).$

We thus obtain that the objective functional $J$ is weakly lower semi-
continuous and complete the proof. ∎

Now, we are in a position to prove the existence of an optimal solution
$\bm{u}$ to (BCP).

######  Theorem 2.3.

There exists at least one optimal control
$\bm{u}\in\mathcal{U}=L^{2}(0,T;\bm{L}_{div}^{2}(\Omega))$ such that
$J(\bm{u})\leq J(\bm{v}),\forall\bm{v}\in\mathcal{U}$.

###### Proof.

We first observe that $J(\bm{v})\geq 0,\forall\bm{v}\in\mathcal{U}$, then the
infimum of $J(\bm{v})$ exists and we denote it as

$j=\inf_{\bm{v}\in\mathcal{U}}J(\bm{v}),$

and there is a minimizing sequence $\\{\bm{v}_{n}\\}\subset\mathcal{U}$ such
that

$\lim_{n\rightarrow\infty}J(\bm{v}_{n})=j.$

This fact, together with
$\frac{1}{2}\|\bm{v}_{n}\|^{2}_{L^{2}(0,T;\bm{L}^{2}_{div}(\Omega))}\leq
J(\bm{v}_{n})$, implies that $\\{\bm{v}_{n}\\}$ is bounded in
$L^{2}(0,T;\bm{L}^{2}_{div}(\Omega))$. Hence, there exists a subsequence,
still denoted by $\\{\bm{v}_{n}\\}$, that converges weakly to $\bm{u}$ in
$L^{2}(0,T;\bm{L}^{2}_{div}(\Omega))$. It follows from Lemma 2.2 that $J$ is
weakly lower semi-continuous and we thus have

$J(\bm{u})\leq\underset{n\rightarrow\infty}{\lim\inf}J(\bm{v}_{n})=j.$

Since $\bm{u}\in\mathcal{U}$, we must have $J(\bm{u})=j$, and $\bm{u}$ is
therefore an optimal control. ∎

We note that the uniqueness of optimal control $\bm{u}$ cannot be guaranteed
and only a local optimal solution can be pursued because the objective
functional $J$ is nonconvex due to the nonlinear relationship between the
state $y$ and the control $\bm{v}$.

###  2.3. First-order Optimality Conditions

Let $DJ(\bm{v})$ be the first-order differential of $J$ at $\bm{v}$ and
$\bm{u}$ an optimal control of (BCP). It is clear that the first-order
optimality condition of (BCP) reads

$DJ(\bm{u})=0.$

In the sequel of this subsection, we discuss the computation of $DJ(\bm{v})$,
which will play an important role in subsequent sections.

To compute $DJ(\bm{v})$, we employ a formal perturbation analysis as in [14].
First, let $\delta\bm{v}\in\mathcal{U}$ be a perturbation of
$\bm{v}\in\mathcal{U}$, we clearly have

(2.7) $\delta J(\bm{v})=\iint_{Q}DJ(\bm{v})\cdot\delta\bm{v}dxdt,$

and also

(2.8) $\displaystyle{}\begin{aligned} &\delta
J(\bm{v})=\iint_{Q}\bm{v}\cdot\delta\bm{v}dxdt+\alpha_{1}\iint_{Q}(y-y_{d})\delta
ydxdt+\alpha_{2}\int_{\Omega}(y(T)-y_{T})\delta y(T)dx,\end{aligned}$

in which $\delta y$ is the solution of

(2.9) $\displaystyle\left\\{\begin{aligned} \frac{\partial\delta y}{\partial
t}-\nu\nabla^{2}\delta y+\delta\bm{v}\cdot\nabla y+\bm{v}\cdot\nabla\delta
y+a_{0}\delta y&=0\quad\text{in}\quad Q,\\\ \delta
y&=0\quad\text{on}\quad\Sigma,\\\ \delta y(0)&=0.\end{aligned}\right.$

Consider now a function $p$ defined over $\overline{Q}$ (the closure of $Q$);
and assume that $p$ is a differentiable function of $x$ and $t$. Multiplying
both sides of the first equation in (2.9) by $p$ and integrating over $Q$, we
obtain

$\iint_{Q}p\frac{\partial}{\partial t}\delta
ydxdt-\nu\iint_{Q}p\nabla^{2}\delta ydxdt+\iint_{Q}\delta\bm{v}\cdot\nabla
ypdxdt+\iint_{Q}\bm{v}\cdot\nabla\delta ypdxdt+a_{0}\iint_{Q}p\delta ydxdt=0.$

Integration by parts in time and application of Green’s formula in space yield

(2.10) $\displaystyle{}\begin{aligned} \int_{\Omega}p(T)\delta
y(T)dx-\int_{\Omega}p(0)\delta y(0)dx+\iint_{Q}\Big{[}-\frac{\partial
p}{\partial t}-\nu\nabla^{2}p-\bm{v}\cdot\nabla p+a_{0}p\Big{]}\delta ydxdt\\\
+\iint_{Q}\delta\bm{v}\cdot\nabla
ypdxdt-\nu\iint_{\Sigma}(\frac{\partial\delta
y}{\partial\bm{n}}p-\frac{\partial p}{\partial\bm{n}}\delta
y)dxdt+\iint_{\Sigma}p\delta y\bm{v}\cdot\bm{n}dxdt=0.\end{aligned}$

where $\bm{n}$ is the unit outward normal vector at $\Gamma$.

Next, let us assume that the function $p$ is the solution to the following
adjoint system

(2.11) $\displaystyle\qquad\left\\{\begin{aligned} -\frac{\partial p}{\partial
t}-\nu\nabla^{2}p-\bm{v}\cdot\nabla
p+a_{0}p&=\alpha_{1}(y-y_{d})~{}\text{in}~{}Q,\\\
p&=0~{}\qquad\quad\quad\text{on}~{}\Sigma,\\\
p(T)&=\alpha_{2}(y(T)-y_{T}).\end{aligned}\right.$

It follows from (2.8), (2.9), (2.10) and (2.11) that

$\delta J(\bm{v})=\iint_{Q}(\bm{v}-p\nabla y)\cdot\delta\bm{v}dxdt.$

which, together with (2.7), implies that

(2.12) $\left\\{\begin{aligned} &DJ(\bm{v})\in\mathcal{U},\\\
&\iint_{Q}DJ(\bm{v})\cdot\bm{z}dxdt=\iint_{Q}(\bm{v}-p\nabla
y)\cdot\bm{z}dxdt,\forall\bm{z}\in\mathcal{U}.\end{aligned}\right.$

From the discussion above, the first-order optimality condition of (BCP) can
be summarized as follows.

######  Theorem 2.4.

Let $\bm{u}\in\mathcal{U}$ be a solution of (BCP). Then, it satisfies the
following optimality condition

$\iint_{Q}(\bm{u}-p\nabla y)\cdot\bm{z}dxdt=0,\forall\bm{z}\in\mathcal{U},$

where $y$ and $p$ are obtained from $\bm{u}$ via the solutions of the
following two parabolic equations:

(state equation) $\displaystyle\quad\qquad\qquad\qquad\left\\{\begin{aligned}
\frac{\partial y}{\partial t}-\nu\nabla^{2}y+\bm{u}\cdot\nabla
y+a_{0}y&=f\quad\text{in}~{}Q,\\\ y&=g\quad\text{on}~{}\Sigma,\\\
y(0)&=\phi,\end{aligned}\right.$

and

(adjoint equation) $\displaystyle\qquad\qquad\qquad\left\\{\begin{aligned}
-\frac{\partial p}{\partial t}-\nu\nabla^{2}p-\bm{u}\cdot\nabla
p+a_{0}p&=\alpha_{1}(y-y_{d})\quad\text{in}~{}Q,\\\
p&=0\quad\qquad\qquad\text{on}~{}\Sigma,\\\
p(T)&=\alpha_{2}(y(T)-y_{T}).\end{aligned}\right.$

##  3\. An Implementable Nested Conjugate Gradient Method

In this section, we discuss the application of a CG strategy to solve (BCP).
In particular, we elaborate on the computation of the gradient and the
stepsize at each CG iteration, and thus obtain an easily implementable
algorithm.

###  3.1. A Generic Conjugate Gradient Method for (BCP)

Conceptually, implementing the CG method to (BCP), we readily obtain the
following algorithm:

  1. (a)

Given $\bm{u}^{0}\in\mathcal{U}$.

  2. (b)

Compute $\bm{g}^{0}=DJ(\bm{u}^{0})$. If $DJ(\bm{u}^{0})=0$, then
$\bm{u}=\bm{u}^{0}$; otherwise set $\bm{w}^{0}=\bm{g}^{0}$.

  3. For $k\geq 0$, $\bm{u}^{k},\bm{g}^{k}$ and $\bm{w}^{k}$ being known, the last two different from $\bm{0}$, one computes $\bm{u}^{k+1},\bm{g}^{k+1}$ and $\bm{w}^{k+1}$ as follows:

  4. (c)

Compute the stepsize $\rho_{k}$ by solving the following optimization problem

(3.1) $\displaystyle\left\\{\begin{aligned} &\rho_{k}\in\mathbb{R},\\\
&J(\bm{u}^{k}-\rho_{k}\bm{w}^{k})\leq
J(\bm{u}^{k}-\rho\bm{w}^{k}),\forall\rho\in\mathbb{R}.\end{aligned}\right.$

  5. (d)

Update $\bm{u}^{k+1}$ and $\bm{g}^{k+1}$, respectively, by

$\bm{u}^{k+1}=\bm{u}^{k}-\rho_{k}\bm{w}^{k},$

and

$\bm{g}^{k+1}=DJ(\bm{u}^{k+1}).$

  6. If $DJ(\bm{u}^{k+1})=0$, take $\bm{u}=\bm{u}^{k+1}$; otherwise,

  7. (e)

Compute

$\beta_{k}=\frac{\iint_{Q}|\bm{g}^{k+1}|^{2}dxdt}{\iint_{Q}|\bm{g}^{k}|^{2}dxdt},$

and then update

$\bm{w}^{k+1}=\bm{g}^{k+1}+\beta_{k}\bm{w}^{k}.$

  8. Do $k+1\rightarrow k$ and return to (c).

The above iterative method looks very simple, but practically, the
implementation of the CG method (a)–(e) for the solution of (BCP) is
nontrivial. In particular, it is numerically challenging to compute
$DJ(\bm{v})$, $\forall\bm{v}\in\mathcal{U}$ and $\rho_{k}$ as illustrated
below. We shall discuss how to address these two issues in the following part
of this section.

###  3.2. Computation of $DJ(\bm{v})$

It is clear that the implementation of the generic CG method (a)–(e) for the
solution of (BCP) requires the knowledge of $DJ(\bm{v})$ for various
$\bm{v}\in\mathcal{U}$, and this has been conceptually provided in (2.12).
However, it is numerically challenging to compute $DJ(\bm{v})$ by (2.12) due
to the restriction $\nabla\cdot DJ(\bm{v})=0$ which ensures that all iterates
$\bm{u}^{k}$ of the CG method meet the additional divergence-free constraint
$\nabla\cdot\bm{u}^{k}=0$. In this subsection, we show that equation (2.12)
can be reformulated as a saddle point problem by introducing a Lagrange
multiplier associated with the constraint $\nabla\cdot DJ(\bm{v})=0$. Then, a
preconditioned CG method is proposed to solve this saddle point problem.

We first note that equation (2.12) can be equivalently reformulated as

(3.2) $\left\\{\begin{aligned} &DJ(\bm{v})(t)\in\mathbb{S},\\\
&\int_{\Omega}DJ(\bm{v})(t)\cdot\bm{z}dx=\int_{\Omega}(\bm{v}(t)-p(t)\nabla
y(t))\cdot\bm{z}dx,\forall\bm{z}\in\mathbb{S},\end{aligned}\right.$

where

$\mathbb{S}=\\{\bm{z}|\bm{z}\in[L^{2}(\Omega)]^{d},\nabla\cdot\bm{z}=0\\}.$

Clearly, problem (3.2) is a particular case of

(3.3) $\left\\{\begin{aligned} &\bm{g}\in\mathbb{S},\\\
&\int_{\Omega}\bm{g}\cdot\bm{z}dx=\int_{\Omega}\bm{f}\cdot\bm{z}dx,\forall\bm{z}\in\mathbb{S},\end{aligned}\right.$

with $\bm{f}$ given in $[L^{2}(\Omega)]^{d}$.

Introducing a Lagrange multiplier $\lambda\in H_{0}^{1}(\Omega)$ associated
with the constraint $\nabla\cdot\bm{z}=0$, it is clear that problem (3.3) is
equivalent to the following saddle point problem

(3.4) $\left\\{\begin{aligned} &(\bm{g},\lambda)\in[L^{2}(\Omega)]^{d}\times
H_{0}^{1}(\Omega),\\\
&\int_{\Omega}\bm{g}\cdot\bm{z}dx=\int_{\Omega}\bm{f}\cdot\bm{z}dx+\int_{\Omega}\lambda\nabla\cdot\bm{z}dx,\forall\bm{z}\in[L^{2}(\Omega)]^{d},\\\
&\int_{\Omega}\nabla\cdot\bm{g}qdx=0,\forall q\in
H_{0}^{1}(\Omega),\end{aligned}\right.$

which is actually a Stokes type problem.

In order to solve problem (3.4), we advocate a CG method inspired from [9,
10]. For this purpose, one has to specify the inner product to be used over
$H_{0}^{1}(\Omega)$. As discussed in [9], the usual $L^{2}$-inner product,
namely, $\\{q,q^{\prime}\\}\rightarrow\int_{\Omega}qq^{\prime}dx$ leads to a
CG method with poor convergence properties. Indeed, using some arguments
similar to those in [8, 9], we can show that the saddle point problem (3.4)
can be reformulated as a linear variational problem in terms of the Lagrange
multiplier $\lambda$. The corresponding coefficient matrix after space
discretization with mesh size $h$ has a condition number of the order of
$h^{-2}$, which is ill-conditioned especially for small $h$ and makes the CG
method converges fairly slow. Hence, preconditioning is necessary for solving
problem (3.4). As suggested in [9], we choose $-\nabla\cdot\nabla$ as a
preconditioner for problem (3.4), and the corresponding preconditioned CG
method operates in the space $H_{0}^{1}(\Omega)$ equipped with the inner
product $\\{q,q^{\prime}\\}\rightarrow\int_{\Omega}\nabla q\cdot\nabla
q^{\prime}dx$ and the associated norm
$\|q\|_{H_{0}^{1}(\Omega)}=(\int_{\Omega}|\nabla q|^{2}dx)^{1/2},\forall
q,q^{\prime}\in H_{0}^{1}(\Omega)$. The resulting algorithm reads as:

  1. G1

Choose $\lambda^{0}\in H_{0}^{1}(\Omega)$.

  2. G2

Solve

$\left\\{\begin{aligned} &\bm{g}^{0}\in[L^{2}(\Omega)]^{d},\\\
&\int_{\Omega}\bm{g}^{0}\cdot\bm{z}dx=\int_{\Omega}\bm{f}\cdot\bm{z}dx+\int_{\Omega}\lambda^{0}\nabla\cdot\bm{z}dx,\forall\bm{z}\in[L^{2}(\Omega)]^{d},\end{aligned}\right.$

and

$\left\\{\begin{aligned} &r^{0}\in H_{0}^{1}(\Omega),\\\ &\int_{\Omega}\nabla
r^{0}\cdot\nabla qdx=\int_{\Omega}\nabla\cdot\bm{g}^{0}qdx,\forall q\in
H_{0}^{1}(\Omega).\end{aligned}\right.$

If $\frac{\int_{\Omega}|\nabla
r^{0}|^{2}dx}{\max\\{1,\int_{\Omega}|\nabla\lambda^{0}|^{2}dx\\}}\leq
tol_{1}$, take $\lambda=\lambda^{0}$ and $\bm{g}=\bm{g}^{0}$; otherwise set
$w^{0}=r^{0}$. For $k\geq 0$, $\lambda^{k},\bm{g}^{k},r^{k}$ and $w^{k}$ being
known with the last two different from 0, we compute
$\lambda^{k+1},\bm{g}^{k+1},r^{k+1}$ and if necessary $w^{k+1}$, as follows:

  3. G3

Solve

$\left\\{\begin{aligned} &\bar{\bm{g}}^{k}\in[L^{2}(\Omega)]^{d},\\\
&\int_{\Omega}\bar{\bm{g}}^{k}\cdot\bm{z}dx=\int_{\Omega}w^{k}\nabla\cdot\bm{z}dx,\forall\bm{z}\in[L^{2}(\Omega)]^{d},\end{aligned}\right.$

and

$\left\\{\begin{aligned} &\bar{r}^{k}\in H_{0}^{1}(\Omega),\\\
&\int_{\Omega}\nabla\bar{r}^{k}\cdot\nabla
qdx=\int_{\Omega}\nabla\cdot\bar{\bm{g}}^{k}qdx,\forall q\in
H_{0}^{1}(\Omega),\end{aligned}\right.$

and compute the stepsize via

$\eta_{k}=\frac{\int_{\Omega}|\nabla
r^{k}|^{2}dx}{\int_{\Omega}\nabla\bar{r}^{k}\cdot\nabla w^{k}dx}.$

  4. G4

Update $\lambda^{k},\bm{g}^{k}$ and $r^{k}$ via

$\lambda^{k+1}=\lambda^{k}-\eta_{k}w^{k},\bm{g}^{k+1}=\bm{g}^{k}-\eta_{k}\bar{\bm{g}}^{k},~{}\text{and}~{}r^{k+1}=r^{k}-\eta_{k}\bar{r}^{k}.$

If $\frac{\int_{\Omega}|\nabla r^{k+1}|^{2}dx}{\max\\{1,\int_{\Omega}|\nabla
r^{0}|^{2}dx\\}}\leq tol_{1}$, take $\lambda=\lambda^{k+1}$ and
$\bm{g}=\bm{g}^{k+1}$; otherwise,

  5. G5

Compute

$\gamma_{k}=\frac{\int_{\Omega}|\nabla r^{k+1}|^{2}dx}{\int_{\Omega}|\nabla
r^{k}|^{2}dx},$

and update $w^{k}$ via

$w^{k+1}=r^{k+1}+\gamma_{k}w^{k}.$

Do $k+1\rightarrow k$ and return to G3.

Clearly, one only needs to solve two simple linear equations at each iteration
of the preconditioned CG algorithm (G1)-(G5), which implies that the algorithm
is easy and cheap to implement. Moreover, due to the well-chosen
preconditioner $-\nabla\cdot\nabla$, one can expect the above preconditioned
CG algorithm to have a fast convergence; this will be validated by the
numerical experiments reported in Section 5.

###  3.3. Computation of the Stepsize $\rho_{k}$

Another crucial step to implement the CG method (a)–(e) is the computation of
the stepsize $\rho_{k}$. It is the solution of the optimization problem (3.1)
which is numerically expensive to be solved exactly or up to a high accuracy.
For instance, to solve (3.1), one may consider the Newton method applied to
the solution of

$H_{k}^{\prime}(\rho_{k})=0,$

where

$H_{k}(\rho)=J(\bm{u}^{k}-\rho\bm{w}^{k}).$

The Newton method requires the second-order derivative
$H_{k}^{\prime\prime}(\rho)$ which can be computed via an iterated adjoint
technique requiring the solution of _four_ parabolic problems per Newton’s
iteration. Hence, the implementation of the Newton method is numerically
expensive.

The high computational load for solving (3.1) motivates us to implement
certain stepsize rule to determine an approximation of $\rho_{k}$. Here, we
advocate the following procedure to compute an approximate stepsize
$\hat{\rho}_{k}$.

For a given $\bm{w}^{k}\in\mathcal{U}$, we replace the state
$y=S(\bm{u}^{k}-\rho\bm{w}^{k})$ in the objective functional
$J(\bm{u}^{k}-\rho\bm{w}^{k})$ by

$S(\bm{u}^{k})-\rho S^{\prime}(\bm{u}^{k})\bm{w}^{k},$

which is indeed the linearization of the mapping $\rho\mapsto
S(\bm{u}^{k}-\rho\bm{w}^{k})$ at $\rho=0$. We thus obtain the following
quadratic approximation of $H_{k}(\rho)$:

(3.5)
$Q_{k}(\rho):=\frac{1}{2}\iint_{Q}|\bm{u}^{k}-\rho\bm{w}^{k}|^{2}dxdt+\frac{\alpha_{1}}{2}\iint_{Q}|y^{k}-\rho
z^{k}-y_{d}|^{2}dxdt+\frac{\alpha_{2}}{2}\int_{\Omega}|y^{k}(T)-\rho
z^{k}(T)-y_{T}|^{2}dx,$

where $y^{k}=S(\bm{u}^{k})$ is the solution of the state equation (1.2)
associated with $\bm{u}^{k}$, and $z^{k}=S^{\prime}(\bm{u}^{k})\bm{w}^{k}$
satisfies the following linear parabolic problem

(3.6) $\displaystyle\left\\{\begin{aligned} \frac{\partial z^{k}}{\partial
t}-\nu\nabla^{2}z^{k}+\bm{w}^{k}\cdot\nabla y^{k}+\bm{u}^{k}\cdot\nabla
z^{k}+a_{0}z^{k}&=0\quad\text{in}\quad Q,\\\
z^{k}&=0\quad\text{on}\quad\Sigma,\\\ z^{k}(0)&=0.\end{aligned}\right.$

Then, it is easy to show that the equation $Q_{k}^{\prime}(\rho)=0$ admits a
unique solution

(3.7)
$\hat{\rho}_{k}=\frac{\iint_{Q}\bm{g}^{k}\cdot\bm{w}^{k}dxdt}{\iint_{Q}|\bm{w}^{k}|^{2}dxdt+\alpha_{1}\iint_{Q}|z^{k}|^{2}dxdt+\alpha_{2}\int_{\Omega}|z^{k}(T)|^{2}dx},$

and we take $\hat{\rho}_{k}$, which is clearly an approximation of $\rho_{k}$,
as the stepsize in each CG iteration.

Altogether, with the stepsize given by (3.7), every iteration of the resulting
CG algorithm requires solving only _three_ parabolic problems, namely, the
state equation (1.2) forward in time and the associated adjoint equation
(2.11) backward in time for the computation of $\bm{g}^{k}$, and to solving
the linearized parabolic equation (3.6) forward in time for the stepsize
$\hat{\rho}_{k}$. For comparison, if the Newton method is employed to compute
the stepsize $\rho_{k}$ by solving (3.1), at least _six_ parabolic problems
are required to be solved at each iteration of the CG method, which is much
more expensive numerically.

######  Remark 3.1.

To find an appropriate stepsize, a natural idea is to employ some line search
strategies, such as the backtracking strategy based on the Armijo–Goldstein
condition or the Wolf condition, see e.g., [25]. It is worth noting that these
line search strategies require the evaluation of $J(\bm{v})$ repeatedly, which
is numerically expensive because every evaluation of $J(\bm{v})$ for a given
$\bm{v}$ requires solving the state equation (1.2). Moreover, we have
implemented the CG method for solving (BCP) with various line search
strategies and observed from the numerical results that line search strategies
always lead to tiny stepsizes making extremely slow the convergence of the CG
method.

###  3.4. A Nested CG Method for Solving (BCP)

Following Sections 3.2 and 3.3, we advocate the following nested CG method for
solving (BCP):

  1. I.

Given $\bm{u}^{0}\in\mathcal{U}$.

  2. II.

Compute $y^{0}$ and $p^{0}$ by solving the state equation (1.2) and the
adjoint equation (2.11) corresponding to $\bm{u}^{0}$. Then, for a.e.
$t\in(0,T)$, solve

$\left\\{\begin{aligned} &\bm{g}^{0}(t)\in\mathbb{S},\\\
&\int_{\Omega}\bm{g}^{0}(t)\cdot\bm{z}dx=\int_{\Omega}(\bm{u}^{0}(t)-p^{0}(t)\nabla
y^{0}(t))\cdot\bm{z}dx,\forall\bm{z}\in\mathbb{S},\end{aligned}\right.$

by the preconditioned CG algorithm (G1)–(G5); and set $\bm{w}^{0}=\bm{g}^{0}.$

For $k\geq 0$, $\bm{u}^{k},\bm{g}^{k}$ and $\bm{w}^{k}$ being known, the last
two different from $\bm{0}$, one computes $\bm{u}^{k+1},\bm{g}^{k+1}$ and
$\bm{w}^{k+1}$ as follows:

  3. III.

Compute the stepsize $\hat{\rho}_{k}$ by (3.7).

  4. IV.

Update $\bm{u}^{k+1}$ by

$\bm{u}^{k+1}=\bm{u}^{k}-\hat{\rho}_{k}\bm{w}^{k}.$

Compute $y^{k+1}$ and $p^{k+1}$ by solving the state equation (1.2) and the
adjoint equation (2.11) corresponding to $\bm{u}^{k+1}$; and for a.e.
$t\in(0,T)$, solve

$\left\\{\begin{aligned} &\bm{g}^{k+1}(t)\in\mathbb{S},\\\
&\int_{\Omega}\bm{g}^{k+1}(t)\cdot\bm{z}dx=\int_{\Omega}(\bm{u}^{k+1}(t)-p^{k+1}(t)\nabla
y^{k+1}(t))\cdot\bm{z}dx,\forall\bm{z}\in\mathbb{S},\end{aligned}\right.$

by the preconditioned CG algorithm (G1)–(G5).

If $\frac{\iint_{Q}|\bm{g}^{k+1}|^{2}dxdt}{\iint_{Q}|\bm{g}^{0}|^{2}dxdt}\leq
tol$, take $\bm{u}=\bm{u}^{k+1}$; else

  5. V.

Compute

$\beta_{k}=\frac{\iint_{Q}|\bm{g}^{k+1}|^{2}dxdt}{\iint_{Q}|\bm{g}^{k}|^{2}dxdt},~{}\text{and}~{}\bm{w}^{k+1}=\bm{g}^{k+1}+\beta_{k}\bm{w}^{k}.$

Do $k+1\rightarrow k$ and return to III.

##  4\. Space and time discretizations

In this section, we discuss first the numerical discretization of the bilinear
optimal control problem (BCP). We achieve the time discretization by a semi-
implicit finite difference method and the space discretization by a piecewise
linear finite element method. Then, we discuss an implementable nested CG
method for solving the fully discrete bilinear optimal control problem.

###  4.1. Time Discretization of (BCP)

First, we define a time discretization step $\Delta t$ by $\Delta t=T/N$, with
$N$ a positive integer. Then, we approximate the control space
$\mathcal{U}=L^{2}(0,T;\mathbb{S})$ by $\mathcal{U}^{\Delta
t}:=(\mathbb{S})^{N};$ and equip $\mathcal{U}^{\Delta t}$ with the following
inner product

$(\bm{v},\bm{w})_{\Delta t}=\Delta
t\sum^{N}_{n=1}\int_{\Omega}\bm{v}_{n}\cdot\bm{w}_{n}dx,\quad\forall\bm{v}=\\{\bm{v}_{n}\\}^{N}_{n=1},\bm{w}=\\{\bm{w}_{n}\\}^{N}_{n=1}\in\mathcal{U}^{\Delta
t},$

and the norm

$\|\bm{v}\|_{\Delta t}=\left(\Delta
t\sum^{N}_{n=1}\int_{\Omega}|\bm{v}_{n}|^{2}dx\right)^{\frac{1}{2}},\quad\forall\bm{v}=\\{\bm{v}_{n}\\}^{N}_{n=1}\in\mathcal{U}^{\Delta
t}.$

Then, (BCP) is approximated by the following semi-discrete bilinear control
problem (BCP)Δt:

$\displaystyle\text{(BCP)}^{\Delta
t}\qquad\qquad\qquad\qquad\left\\{\begin{aligned} &\bm{u}^{\Delta
t}\in\mathcal{U}^{\Delta t},\\\ &J^{\Delta t}(\bm{u}^{\Delta t})\leq J^{\Delta
t}(\bm{v}),\forall\bm{v}=\\{\bm{v}_{n}\\}_{n=1}^{N}\in\mathcal{U}^{\Delta
t},\end{aligned}\right.$

where the cost functional $J^{\Delta t}$ is defined by

$J^{\Delta t}(\bm{v})=\frac{1}{2}\Delta
t\sum^{N}_{n=1}\int_{\Omega}|\bm{v}_{n}|^{2}dx+\frac{\alpha_{1}}{2}\Delta
t\sum^{N}_{n=1}\int_{\Omega}|y_{n}-y_{d}^{n}|^{2}dx+\frac{\alpha_{2}}{2}\int_{\Omega}|y_{N}-y_{T}|^{2}dx,$

with $\\{y_{n}\\}^{N}_{n=1}$ the solution of the following semi-discrete state
equation: $y_{0}=\phi$; then for $n=1,\ldots,N$, with $y_{n-1}$ being known,
we obtain $y_{n}$ from the solution of the following linear elliptic problem:

(4.1) $\displaystyle\left\\{\begin{aligned} \frac{{y}_{n}-{y}_{n-1}}{\Delta
t}-\nu\nabla^{2}{y}_{n}+\bm{v}_{n}\cdot\nabla{y}_{n-1}+a_{0}{y}_{n-1}&=f_{n}\quad\text{in}\quad\Omega,\\\
y_{n}&=g_{n}\quad\text{on}\quad\Gamma.\end{aligned}\right.$

######  Remark 4.1.

For simplicity, we have chosen a one-step semi-explicit scheme to discretize
system (1.2). This scheme is first-order accurate and reasonably robust, once
combined to an appropriate space discretization. The application of second-
order accurate time discretization schemes to optimal control problems has
been discussed in e.g., [5].

######  Remark 4.2.

At each step of scheme (4.1), we only need to solve a simple linear elliptic
problem to obtain $y_{n}$ from $y_{n-1}$, and there is no particular
difficulty in solving such a problem.

The existence of a solution to the semi-discrete bilinear optimal control
problem (BCP)Δt can be proved in a similar way as what we have done for the
continuous case. Let $\bm{u}^{\Delta t}$ be a solution of (BCP)Δt, then it
verifies the following first-order optimality condition:

$DJ^{\Delta t}(\bm{u}^{\Delta t})=0,$

where $DJ^{\Delta t}(\bm{v})$ is the first-order differential of the
functional $J^{\Delta t}$ at $\bm{v}\in\mathcal{U}^{\Delta t}$.

Proceeding as in the continuous case, we can show that $DJ^{\Delta
t}(\bm{v})=\\{\bm{g}_{n}\\}_{n=1}^{N}\in\mathcal{U}^{\Delta t}$ where

$\left\\{\begin{aligned} &\bm{g}_{n}\in\mathbb{S},\\\
&\int_{\Omega}\bm{g}_{n}\cdot\bm{w}dx=\int_{\Omega}(\bm{v}_{n}-p_{n}\nabla
y_{n-1})\cdot\bm{w}dx,\forall\bm{w}\in\mathbb{S},\end{aligned}\right.$

and the vector-valued function $\\{p_{n}\\}^{N}_{n=1}$ is the solution of the
semi-discrete adjoint system below:

${p}_{N+1}=\alpha_{2}({y}_{N}-y_{T});$

for $n=N$, solve

$\displaystyle\qquad\left\\{\begin{aligned} \frac{{p}_{N}-{p}_{N+1}}{\Delta
t}-\nu\nabla^{2}{p}_{N}&=\alpha_{1}({y}_{N}-y_{d}^{N})&\quad\text{in}\quad\Omega,\\\
p_{N}&=0&\quad\text{on}\quad\Gamma,\end{aligned}\right.$

and for $n=N-1,\cdots,1,$ solve

$\displaystyle\qquad\left\\{\begin{aligned} \frac{{p}_{n}-{p}_{n+1}}{\Delta
t}-\nu\nabla^{2}{p}_{n}-\bm{v}_{n+1}\cdot\nabla{p}_{n+1}+a_{0}{p}_{n+1}&=\alpha_{1}({y}_{n}-y_{d}^{n})&\quad\text{in}\quad\Omega,\\\
p_{n}&=0&\quad\text{on}\quad\Gamma.\end{aligned}\right.$

###  4.2. Space Discretization of (BCP)Δt

In this subsection, we discuss the space discretization of (BCP)Δt, obtaining
thus a full space-time discretization of (BCP). For simplicity, we suppose
from now on that $\Omega$ is a polygonal domain of $\mathbb{R}^{2}$ (or has
been approximated by a family of such domains).

Let $\mathcal{T}_{H}$ be a classical triangulation of $\Omega$, with $H$ the
largest length of the edges of the triangles of $\mathcal{T}_{H}$. From
$\mathcal{T}_{H}$ we construct $\mathcal{T}_{h}$ with $h=H/2$ by joining the
mid-points of the edges of the triangles of $\mathcal{T}_{H}$.

We first consider the finite element space $V_{h}$ defined by

$V_{h}=\\{\varphi_{h}|\varphi_{h}\in
C^{0}(\bar{\Omega});{\varphi_{h}\mid}_{\mathbb{T}}\in
P_{1},\forall\,{\mathbb{T}}\in\mathcal{T}_{h}\\}$

with $P_{1}$ the space of the polynomials of two variables of degree $\leq 1$.
Two useful sub-spaces of $V_{h}$ are

$V_{0h}=\\{\varphi_{h}|\varphi_{h}\in
V_{h},\varphi_{h}\mid_{\Gamma}=0\\}:=V_{h}\cap H_{0}^{1}(\Omega),$

and (assuming that $g(t)\in C^{0}(\Gamma)$)

$\displaystyle V_{gh}(t)=\\{\varphi_{h}|\varphi_{h}\in
V_{h},\varphi_{h}(Q)=g(Q,t),\forall\,Q~{}\text{vertex
of}~{}\mathcal{T}_{h}~{}\text{located on}~{}\Gamma\\}.$

In order to construct the discrete control space, we introduce first

$\Lambda_{H}=\\{\varphi_{H}|\varphi_{H}\in
C^{0}(\bar{\Omega});{\varphi_{H}\mid}_{\mathbb{T}}\in
P_{1},\forall\,{\mathbb{T}}\in\mathcal{T}_{H}\\},~{}\text{and}~{}\Lambda_{0H}=\\{\varphi_{H}|\varphi_{H}\in\Lambda_{H},\varphi_{H}\mid_{\Gamma}=0\\}.$

Then, the discrete control space $\mathcal{U}_{h}^{\Delta t}$ is defined by

$\mathcal{U}_{h}^{\Delta
t}=(\mathbb{S}_{h})^{N},~{}\text{with}~{}\mathbb{S}_{h}=\\{\bm{v}_{h}|\bm{v}_{h}\in
V_{h}\times
V_{h},\int_{\Omega}\nabla\cdot\bm{v}_{h}q_{H}dx\left(=-\int_{\Omega}\bm{v}_{h}\cdot\nabla
q_{H}dx\right)=0,\forall q_{H}\in\Lambda_{0H}\\}.$

With the above finite element spaces, we approximate (BCP) and (BCP)Δt by
(BCP)${}_{h}^{\Delta t}$ defined by

$\displaystyle\text{(BCP)}_{h}^{\Delta
t}\qquad\qquad\qquad\qquad\qquad\qquad\left\\{\begin{aligned}
&\bm{u}_{h}^{\Delta t}\in\mathcal{U}_{h}^{\Delta t},\\\ &J_{h}^{\Delta
t}(\bm{u}_{h}^{\Delta t})\leq J_{h}^{\Delta t}(\bm{v}_{h}^{\Delta
t}),\forall\bm{v}_{h}^{\Delta t}\in\mathcal{U}_{h}^{\Delta
t},\end{aligned}\right.$

where the fully discrete cost functional $J_{h}^{\Delta t}$ is defined by

(4.2) $J_{h}^{\Delta t}(\bm{v}_{h}^{\Delta t})=\frac{1}{2}\Delta
t\sum^{N}_{n=1}\int_{\Omega}|\bm{v}_{n,h}|^{2}dx+\frac{\alpha_{1}}{2}\Delta
t\sum^{N}_{n=1}\int_{\Omega}|y_{n,h}-y_{d}^{n}|^{2}dx+\frac{\alpha_{2}}{2}\int_{\Omega}|y_{N,h}-y_{T}|^{2}dx$

with $\\{y_{n,h}\\}^{N}_{n=1}$ the solution of the following fully discrete
state equation: $y_{0,h}=\phi_{h}\in V_{h}$, where $\phi_{h}$ verifies

$\phi_{h}\in V_{h},\forall\,h>0,~{}\text{and}~{}\lim_{h\rightarrow
0}\phi_{h}=\phi,~{}\text{in}~{}L^{2}(\Omega),$

then, for $n=1,\ldots,N$, with $y_{n-1,h}$ being known, we obtain $y_{n,h}\in
V_{gh}(n\Delta t)$ from the solution of the following linear variational
problem:

(4.3) $\int_{\Omega}\frac{{y}_{n,h}-{y}_{n-1,h}}{\Delta t}\varphi
dx+\nu\int_{\Omega}\nabla{y}_{n,h}\cdot\nabla\varphi
dx+\int_{\Omega}\bm{v}_{n}\cdot\nabla{y}_{n-1,h}\varphi
dx+\int_{\Omega}a_{0}{y}_{n-1,h}\varphi dx=\int_{\Omega}f_{n}\varphi
dx,\forall\varphi\in V_{0h}.$

In the following discussion, the subscript $h$ in all variables will be
omitted for simplicity.

In a similar way as what we have done in the continuous case, one can show
that the first-order differential of $J_{h}^{\Delta t}$ at
$\bm{v}\in\mathcal{U}_{h}^{\Delta t}$ is $DJ_{h}^{\Delta
t}(\bm{v})=\\{\bm{g}_{n}\\}_{n=1}^{N}\in(\mathbb{S}_{h})^{N}$ where

(4.4) $\left\\{\begin{aligned} &\bm{g}_{n}\in\mathbb{S}_{h},\\\
&\int_{\Omega}\bm{g}_{n}\cdot\bm{z}dx=\int_{\Omega}(\bm{v}_{n}-p_{n}\nabla
y_{n-1})\cdot\bm{z}dx,\forall\bm{z}\in\mathbb{S}_{h},\end{aligned}\right.$

and the vector-valued function $\\{p_{n}\\}^{N}_{n=1}$ is the solution of the
following fully discrete adjoint system:

(4.5) ${p}_{N+1}=\alpha_{2}({y}_{N}-y_{T});$

for $n=N$, solve

(4.6) $\displaystyle\qquad\left\\{\begin{aligned} &p_{N}\in V_{0h},\\\
&\int_{\Omega}\frac{{p}_{N}-{p}_{N+1}}{\Delta t}\varphi
dx+\nu\int_{\Omega}\nabla{p}_{N}\cdot\nabla\varphi
dx=\int_{\Omega}\alpha_{1}({y}_{N}-y_{d}^{N})\varphi dx,\forall\varphi\in
V_{0h},\end{aligned}\right.$

then, for $n=N-1,\cdots,1,$, solve

(4.7) $\displaystyle\qquad\left\\{\begin{aligned} &p_{n}\in V_{0h},\\\
&\int_{\Omega}\frac{{p}_{n}-{p}_{n+1}}{\Delta t}\varphi
dx+\nu\int_{\Omega}\nabla{p}_{n}\cdot\nabla\varphi
dx-\int_{\Omega}\bm{v}_{n+1}\cdot\nabla{p}_{n+1}\varphi dx\\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad+a_{0}\int_{\Omega}{p}_{n+1}\varphi
dx=\int_{\Omega}\alpha_{1}({y}_{n}-y_{d}^{n})\varphi dx,\forall\varphi\in
V_{0h}.\end{aligned}\right.$

It is worth mentioning that the so-called discretize-then-optimize strategy is
employed here, which implies that we first discretize (BCP), and to compute
the gradient in a discrete setting, the fully discrete adjoint equation
(4.5)–(4.7) has been derived from the fully discrete cost functional
$J_{h}^{\Delta t}(\bm{v})$ (4.2) and the fully discrete state equation (4.3).
This implies that the fully discrete state equation (4.3) and the fully
discrete adjoint equation (4.5)–(4.7) are strictly in duality. This fact
guarantees that $-DJ_{h}^{\Delta t}(\bm{v})$ is a descent direction of the
fully discrete bilinear optimal control problem (BCP)${}_{h}^{\Delta t}$.

######  Remark 4.3.

A natural alternative has been advocated in the literature: (i) Derive the
adjoint equation to compute the first-order differential of the cost
functional in a continuous setting; (ii) Discretize the state and adjoint
state equations by certain numerical schemes; (iii) Use the resulting discrete
analogs of $y$ and $p$ to compute a discretization of the differential of the
cost functional. The main problem with this optimize-then-discretize approach
is that it may not preserve a strict duality between the discrete state
equation and the discrete adjoint equation. This fact implies in turn that the
resulting discretization of the continuous gradient may not be a gradient of
the discrete optimal control problem. As a consequence, the resulting
algorithm is not a descent algorithm and divergence may take place (see [11]
for a related discussion).

###  4.3. A Nested CG Method for Solving the Fully Discrete Problem
(BCP)${}_{h}^{\Delta t}$

In this subsection, we propose a nested CG method for solving the fully
discrete problem (BCP)${}_{h}^{\Delta t}$. As discussed in Section 3, the
implementation of CG requires the knowledge of $DJ_{h}^{\Delta t}(\bm{v})$ and
an appropriate stepsize. In the following discussion, we address these two
issues by extending the results for the continuous case in Sections 3.2 and
3.3 to the fully discrete settings; and derive the corresponding CG algorithm.

First, it is clear that one can compute $DJ_{h}^{\Delta t}(\bm{v})$ via the
solution of the $N$ linear variational problems encountered in (4.4). For this
purpose, we introduce a Lagrange multiplier $\lambda\in\Lambda_{0H}$
associated with the divergence-free constraint, then problem (4.4) is
equivalent to the following saddle point system

(4.8) $\left\\{\begin{aligned} &(\bm{g}_{n},\lambda)\in(V_{h}\times
V_{h})\times\Lambda_{0H},\\\
&\int_{\Omega}\bm{g}_{n}\cdot\bm{z}dx=\int_{\Omega}(\bm{v}_{n}-p_{n}\nabla
y_{n-1})\cdot\bm{z}dx+\int_{\Omega}\lambda\nabla\cdot\bm{z}dx,\forall\bm{z}\in
V_{h}\times V_{h},\\\ &\int_{\Omega}\nabla\cdot\bm{g}_{n}qdx=0,\forall
q\in\Lambda_{0H}.\end{aligned}\right.$

As discussed in Section 3.2, problem (4.8) can be solved by the following
preconditioned CG algorithm, which is actually a discrete analogue of
(G1)–(G5).

  1. DG1

Choose $\lambda^{0}\in\Lambda_{0H}$.

  2. DG2

Solve

$\left\\{\begin{aligned} &\bm{g}_{n}^{0}\in V_{h}\times V_{h},\\\
&\int_{\Omega}\bm{g}_{n}^{0}\cdot\bm{z}dx=\int_{\Omega}(\bm{v}_{n}-p_{n}\nabla
y_{n-1})\cdot\bm{z}dx+\int_{\Omega}\lambda^{0}\nabla\cdot\bm{z}dx,\forall\bm{z}\in
V_{h}\times V_{h},\end{aligned}\right.$

and

$\left\\{\begin{aligned} &r^{0}\in\Lambda_{0H},\\\ &\int_{\Omega}\nabla
r^{0}\cdot\nabla qdx=\int_{\Omega}\nabla\cdot\bm{g}_{n}^{0}qdx,\forall
q\in\Lambda_{0H}.\end{aligned}\right.$

If $\frac{\int_{\Omega}|\nabla
r^{0}|^{2}dx}{\max\\{1,\int_{\Omega}|\nabla\lambda^{0}|^{2}dx\\}}\leq
tol_{1}$, take $\lambda=\lambda^{0}$ and $\bm{g}_{n}=\bm{g}_{n}^{0}$;
otherwise set $w^{0}=r^{0}$. For $k\geq 0$, $\lambda^{k},\bm{g}_{n}^{k},r^{k}$
and $w^{k}$ being known with the last two different from 0, we define
$\lambda^{k+1},\bm{g}_{n}^{k+1},r^{k+1}$ and if necessary $w^{k+1}$, as
follows:

  3. DG3

Solve

$\left\\{\begin{aligned} &\bar{\bm{g}}_{n}^{k}\in V_{h}\times V_{h},\\\
&\int_{\Omega}\bar{\bm{g}}_{n}^{k}\cdot\bm{z}dx=\int_{\Omega}w^{k}\nabla\cdot\bm{z}dx,\forall\bm{z}\in
V_{h}\times V_{h},\end{aligned}\right.$

and

$\left\\{\begin{aligned} &\bar{r}^{k}\in\Lambda_{0H},\\\
&\int_{\Omega}\nabla\bar{r}^{k}\cdot\nabla
qdx=\int_{\Omega}\nabla\cdot\bar{\bm{g}}_{n}^{k}qdx,\forall
q\in\Lambda_{0H},\end{aligned}\right.$

and compute

$\eta_{k}=\frac{\int_{\Omega}|\nabla
r^{k}|^{2}dx}{\int_{\Omega}\nabla\bar{r}^{k}\cdot\nabla w^{k}dx}.$

  4. DG4

Update $\lambda^{k},\bm{g}_{n}^{k}$ and $r^{k}$ via

$\lambda^{k+1}=\lambda^{k}-\eta_{k}w^{k},\bm{g}_{n}^{k+1}=\bm{g}_{n}^{k}-\eta_{k}\bar{\bm{g}}_{n}^{k},~{}\text{and}~{}r^{k+1}=r^{k}-\eta_{k}\bar{r}^{k}.$

If $\frac{\int_{\Omega}|\nabla r^{k+1}|^{2}dx}{\max\\{1,\int_{\Omega}|\nabla
r^{0}|^{2}dx\\}}\leq tol_{1}$, take $\lambda=\lambda^{k+1}$ and
$\bm{g}_{n}=\bm{g}_{n}^{k+1}$; otherwise,

  5. DG5

Compute

$\gamma_{k}=\frac{\int_{\Omega}|\nabla r^{k+1}|^{2}dx}{\int_{\Omega}|\nabla
r^{k}|^{2}dx},$

and update $w^{k}$ via

$w^{k+1}=r^{k+1}+\gamma_{k}w^{k}.$

Do $k+1\rightarrow k$ and return to DG3.

To find an appropriate stepsize in the CG iteration for the solution of
(BCP)${}_{h}^{\Delta t}$, we note that, for any
$\\{\bm{w}_{n}\\}_{n=1}^{N}\in(\mathbb{S}_{h})^{N}$, the fully discrete
analogue of $Q_{k}(\rho)$ in (3.5) reads as

$Q_{h}^{\Delta t}(\rho)=\frac{1}{2}\Delta
t\sum^{N}_{n=1}\int_{\Omega}|\bm{u}_{n}-\rho\bm{w}_{n}|^{2}dx+\frac{\alpha_{1}}{2}\Delta
t\sum^{N}_{n=1}\int_{\Omega}|y_{n}-\rho
z_{n}-y_{d}^{n}|^{2}dx+\frac{\alpha_{2}}{2}\int_{\Omega}|y_{N}-\rho
z_{N}-y_{T}|^{2}dx,$

where the vector-valued function $\\{z_{n}\\}^{N}_{n=1}$ is obtained as
follows: $z_{0}=0$; then for $n=1,\ldots,N$, with $z_{n-1}$ being known,
$z_{n}$ is obtained from the solution of the linear variational problem

$\left\\{\begin{aligned} &z_{n}\in V_{0h},\\\
&\int_{\Omega}\frac{{z}_{n}-{z}_{n-1}}{\Delta t}\varphi
dx+\nu\int_{\Omega}\nabla{z}_{n}\cdot\nabla\varphi
dx+\int_{\Omega}\bm{w}_{n}\cdot\nabla y_{n}\varphi dx\\\
&\qquad\qquad\qquad\qquad\qquad+\int_{\Omega}\bm{u}_{n}\cdot\nabla{z}_{n-1}\varphi
dx+a_{0}\int_{\Omega}{z}_{n-1}\varphi dx=0,\forall\varphi\in V_{0h}.\\\
\end{aligned}\right.$

As discussed in Section 3.3 for the continuous case, we take the unique
solution of ${Q_{h}^{\Delta t}}^{\prime}(\rho)=0$ as the stepsize in each CG
iteration, that is

(4.9) $\hat{\rho}_{h}^{\Delta t}=\frac{\Delta
t\sum_{n=1}^{N}\int_{\Omega}\bm{g}_{n}\cdot\bm{w}_{n}dx}{\Delta
t\sum_{n=1}^{N}\int_{\Omega}|\bm{w}_{n}|^{2}dxdt+\alpha_{1}\Delta
t\sum_{n=1}^{N}\int_{\Omega}|z_{n}|^{2}dxdt+\alpha_{2}\int_{\Omega}|z_{N}|^{2}dx}.$

Finally, with above preparations, we propose the following nested CG algorithm
for the solution of the fully discrete control problem (BCP)${}_{h}^{\Delta
t}$.

  1. DI.

Given $\bm{u}^{0}:=\\{\bm{u}_{n}^{0}\\}_{n=1}^{N}\in(\mathbb{S}_{h})^{N}$.

  2. DII.

Compute $\\{y_{n}^{0}\\}_{n=0}^{N}$ and $\\{p^{0}_{n}\\}_{n=1}^{N+1}$ by
solving the fully discrete state equation (4.3) and the fully discrete adjoint
equation (4.5)–(4.7) corresponding to $\bm{u}^{0}$. Then, for $n=1,\cdots,N$
solve

$\left\\{\begin{aligned} &\bm{g}_{n}^{0}\in\mathbb{S}_{h},\\\
&\int_{\Omega}\bm{g}_{n}^{0}\cdot\bm{z}dx=\int_{\Omega}(\bm{u}_{n}^{0}-p_{n}^{0}\nabla
y_{n-1}^{0})\cdot\bm{z}dx,\forall\bm{z}\in\mathbb{S}_{h},\end{aligned}\right.$

by the preconditioned CG algorithm (DG1)–(DG5), and set
$\bm{w}^{0}_{n}=\bm{g}_{n}^{0}.$

For $k\geq 0$, $\bm{u}^{k},\bm{g}^{k}$ and $\bm{w}^{k}$ being known, the last
two different from $\bm{0}$, one computes $\bm{u}^{k+1},\bm{g}^{k+1}$ and
$\bm{w}^{k+1}$ as follows:

  3. DIII.

Compute the stepsize $\hat{\rho}_{k}$ by (4.9).

  4. DIV.

Update $\bm{u}^{k+1}$ by

$\bm{u}^{k+1}=\bm{u}^{k}-\hat{\rho}_{k}\bm{w}^{k}.$

Compute $\\{y_{n}^{k+1}\\}_{n=0}^{N}$ and $\\{p_{n}^{k+1}\\}_{n=1}^{N+1}$ by
solving the fully discrete state equation (4.3) and the fully discrete adjoint
equation (4.5)–(4.7) corresponding to $\bm{u}^{k+1}$. Then, for
$n=1,\cdots,N$, solve

(4.10) $\left\\{\begin{aligned} &\bm{g}_{n}^{k+1}\in\mathbb{S}_{h},\\\
&\int_{\Omega}\bm{g}_{n}^{k+1}\cdot\bm{z}dx=\int_{\Omega}(\bm{u}_{n}^{k+1}-p_{n}^{k+1}\nabla
y_{n-1}^{k+1})\cdot\bm{z}dx,\forall\bm{z}\in\mathbb{S}_{h},\end{aligned}\right.$

by the preconditioned CG algorithm (DG1)–(DG5).

If $\frac{\Delta t\sum_{n=1}^{N}\int_{\Omega}|\bm{g}_{n}^{k+1}|^{2}dx}{\Delta
t\sum_{n=1}^{N}\int_{\Omega}|\bm{g}_{n}^{0}|^{2}dx}\leq tol$, take
$\bm{u}=\bm{u}^{k+1}$; else

  5. DV.

Compute

$\beta_{k}=\frac{\Delta
t\sum_{n=1}^{N}\int_{\Omega}|\bm{g}_{n}^{k+1}|^{2}dx}{\Delta
t\sum_{n=1}^{N}\int_{\Omega}|\bm{g}_{n}^{k}|^{2}dx},~{}\text{and}~{}\bm{w}^{k+1}=\bm{g}^{k+1}+\beta_{k}\bm{w}^{k}.$

Do $k+1\rightarrow k$ and return to DIII.

Despite its apparent complexity, the CG algorithm (DI)-(DV) is easy to
implement. Actually, one of the main computational difficulties in the
implementation of the above algorithm seems to be the solution of $N$ linear
systems (4.10), which is time-consuming. However, it is worth noting that the
linear systems (4.10) are separable with respect to different $n$ and they can
be solved in parallel. As a consequent, one can compute the gradient
$\\{\bm{g}^{k}_{n}\\}_{n=1}^{N}$ simultaneously and the computation time can
be reduced significantly.

Moreover, it is clear that the computation of $\\{\bm{g}^{k}_{n}\\}_{n=1}^{N}$
requires the storage of the solutions of (4.3) and (4.5)-(4.7) at all points
in space and time. For large scale problems, especially in three space
dimensions, it will be very memory demanding and maybe even impossible to
store the full sets $\\{y_{n}^{k}\\}_{n=0}^{N}$ and
$\\{p_{n}^{k}\\}_{n=1}^{N+1}$ simultaneously. To tackle this issue, one can
employ the strategy described in e.g., [14, Section 1.12] that can drastically
reduce the storage requirements at the expense of a small CPU increase.

##  5\. Numerical Experiments

In this section, we report some preliminary numerical results validating the
efficiency of the proposed CG algorithm (DI)–(DV) for (BCP). All codes were
written in MATLAB R2016b and numerical experiments were conducted on a Surface
Pro 5 laptop with 64-bit Windows 10.0 operation system, Intel(R) Core(TM)
i7-7660U CPU (2.50 GHz), and 16 GB RAM.

Example 1. We consider the bilinear optimal control problem (BCP) on the
domain $Q=\Omega\times(0,T)$ with $\Omega=(0,1)^{2}$ and $T=1$. In particular,
we take the control $\bm{v}(x,t)$ in a finite-dimensional space, i.e.
$\bm{v}\in L^{2}(0,T;\mathbb{R}^{2})$. In addition, we set $\alpha_{2}=0$ in
(1.1) and consider the following tracking-type bilinear optimal control
problem:

(5.1) $\min_{\bm{v}\in
L^{2}(0,T;\mathbb{R}^{2})}J(\bm{v})=\frac{1}{2}\int_{0}^{T}|\bm{v}(t)|^{2}dt+\frac{\alpha_{1}}{2}\iint_{Q}|y-y_{d}|^{2}dxdt,$

where $|\bm{v}(t)|=\sqrt{\bm{v}_{1}(t)^{2}+\bm{v}_{2}(t)^{2}}$ is the
canonical 2-norm, and $y$ is obtained from $\bm{v}$ via the solution of the
state equation (1.2).

Since the control $\bm{v}$ is considered in a finite-dimensional space, the
divergence-free constraint $\nabla\cdot\bm{v}=0$ is verified automatically. As
a consequence, the first-order differential $DJ(\bm{v})$ can be easily
computed. Indeed, it is easy to show that

(5.2) $DJ(\bm{v})=\left\\{\bm{v}_{i}(t)+\int_{\Omega}y(t)\frac{\partial
p(t)}{\partial
x_{i}}dx\right\\}_{i=1}^{2},~{}\text{a.e.~{}on}~{}(0,T),\forall\bm{v}\in
L^{2}(0,T;\mathbb{R}^{2}),$

where $p(t)$ is the solution of the adjoint equation (2.11). The inner
preconditioned CG algorithm (DG1)-(DG5) for the computation of the gradient
$\\{\bm{g}_{n}\\}_{n=1}^{N}$ is thus avoided.

In order to examine the efficiency of the proposed CG algorithm (DI)–(DV), we
construct an example with a known exact solution. To this end, we set $\nu=1$
and $a_{0}=1$ in (1.2), and

$y=e^{t}(-3\sin(2\pi x_{1})\sin(\pi x_{2})+1.5\sin(\pi x_{1})\sin(2\pi
x_{2})),\quad p=(T-t)\sin\pi x_{1}\sin\pi x_{2}.$

Substituting these two functions into the optimality condition
$DJ(\bm{u}(t))=0$, we have

$\bm{u}=(\bm{u}_{1},\bm{u}_{2})^{\top}=(2e^{t}(T-t),-e^{t}(T-t))^{\top}.$

We further set

$\displaystyle f=\frac{\partial y}{\partial t}-\nabla^{2}y+{\bm{u}}\cdot\nabla
y+y,\quad\phi=-3\sin(2\pi x_{1})\sin(\pi x_{2})+1.5\sin(\pi x_{1})\sin(2\pi
x_{2}),$ $\displaystyle y_{d}=y-\frac{1}{\alpha_{1}}\left(-\frac{\partial
p}{\partial t}-\nabla^{2}p-\bm{u}\cdot\nabla p+p\right),\quad g=0.$

Then, it is easy to verify that $\bm{u}$ is a solution point of the problem
(5.1). We display the solution $\bm{u}$ and the target function $y_{d}$ at
different instants of time in Figure 1 and Figure 2, respectively.

Figure 1. The exact optimal control $\bm{u}$ for Example 1.

Figure 2. The target function $y_{d}$ at $t=0.25,0.5$ and $0.75$ (from left to
right) for Example 1.

The stopping criterion of the CG algorithm (DI)–(DV) is set as

$\frac{\Delta t\sum_{n=1}^{N}|\bm{g}^{k+1}_{n}|^{2}}{\Delta
t\sum_{n=1}^{N}|\bm{g}^{0}_{n}|^{2}}\leq 10^{-5}.$

The initial value is chosen as $\bm{u}^{0}=(0,0)^{\top}$; and we denote by
$\bm{u}^{\Delta t}$ and $y_{h}^{\Delta t}$ the computed control and state,
respectively.

First, we take $h=\frac{1}{2^{i}},i=5,6,7,8$, $\Delta t=\frac{h}{2}$ and
$\alpha_{1}=10^{6}$, and implement the proposed CG algorithm (DI)–(DV) for
solving the problem (5.1). The numerical results reported in Table 1 show that
the CG algorithm converges fairly fast and is robust with respect to different
mesh sizes. We also observe that the target function $y_{d}$ has been reached
within a good accuracy. Similar comments hold for the approximation of the
optimal control $\bm{u}$ and of the state $y$ of problem (5.1). By taking
$h=\frac{1}{2^{7}}$ and $\Delta t=\frac{1}{2^{8}}$, the computed state
$y_{h}^{\Delta t}$ and $y_{h}^{\Delta t}-y_{d}$ at $t=0.25,0.5$ and $0.75$ are
reported in Figures 3, 4 and 5, respectively; and the computed control
$\bm{u}^{\Delta t}$ and error $\bm{u}^{\Delta t}-\bm{u}$ are visualized in
Figure 6.

Table 1. Results of the CG algorithm (DI)–(DV) with different $h$ and $\Delta t$ for Example 1. Mesh sizes | $Iter$ | $\|\bm{u}^{\Delta t}-\bm{u}\|_{L^{2}(0,T;\mathbb{R}^{2})}$ | $\|y_{h}^{\Delta t}-y\|_{L^{2}(Q)}$ | ${\|y_{h}^{\Delta t}-y_{d}\|_{L^{2}(Q)}}/{\|y_{d}\|_{{L^{2}(Q)}}}$
---|---|---|---|---
$h=1/2^{5},\Delta t=1/2^{6}$ | 117 |  2.8820$\times 10^{-2}$ |  1.1569$\times 10^{-2}$ |  3.8433$\times 10^{-3}$
$h=1/2^{6},\Delta t=1/2^{7}$ | 48 |  1.3912$\times 10^{-2}$ |  2.5739$\times 10^{-3}$ |  8.5623$\times 10^{-4}$
$h=1/2^{7},\Delta t=1/2^{8}$ | 48 |  6.9095$\times 10^{-3}$ |  4.8574$\times 10^{-4}$ |  1.6516$\times 10^{-4}$
$h=1/2^{8},\Delta t=1/2^{9}$ | 31 |  3.4845$\times 10^{-3}$ |  6.6231$\times 10^{-5}$ |  2.2196$\times 10^{-5}$

Figure 3. Computed state $y^{\Delta t}_{h}$, error $y^{\Delta t}_{h}-y$ and
$y^{\Delta t}_{h}-y_{d}$ (from left to right) at $t=0.25$ for Example 1.

Figure 4. Computed state $y^{\Delta t}_{h}$, error $y^{\Delta t}_{h}-y$ and
$y^{\Delta t}_{h}-y_{d}$ (from left to right) at $t=0.5$ for Example 1.

Figure 5. Computed state $y^{\Delta t}_{h}$, error $y^{\Delta t}_{h}-y$ and
$y^{\Delta t}_{h}-y_{d}$ (from left to right) at $t=0.75$ for Example 1.

Figure 6. Computed optimal control $\bm{u}^{\Delta t}$ and error
$\bm{u}^{\Delta t}-\bm{u}$ for Example 1.

Furthermore, we tested the proposed CG algorithm (DI)–(DV) with
$h=\frac{1}{2^{6}}$ and $\Delta t=\frac{1}{2^{7}}$ for different penalty
parameter $\alpha_{1}$. The results reported in Table 2 show that the
performance of the proposed CG algorithm is robust with respect to the penalty
parameter, at least for the example being considered. We also observe that as
$\alpha_{1}$ increases, the value of $\frac{\|y_{h}^{\Delta
t}-y_{d}\|_{L^{2}(Q)}}{\|y_{d}\|_{{L^{2}(Q)}}}$ decreases. This implies that,
as expected, the computed state $y_{h}^{\Delta t}$ is closer to the target
function $y_{d}$ when the penalty parameter gets larger.

Table 2. Results of the CG algorithm (DI)–(DV) with different $\alpha_{1}$ for Example 1. $\alpha_{1}$ | $Iter$ | $CPU(s)$ | $\|\bm{u}^{\Delta t}-\bm{u}\|_{L^{2}(0,T;\mathbb{R}^{2})}$ | $\|y_{h}^{\Delta t}-y\|_{L^{2}(Q)}$ | $\frac{\|y_{h}^{\Delta t}-y_{d}\|_{L^{2}(Q)}}{\|y_{d}\|_{{L^{2}(Q)}}}$
---|---|---|---|---|---
$10^{4}$ | 46 | 126.0666 |  1.3872$\times 10^{-2}$ |  2.5739$\times 10^{-3}$ |  8.7666$\times 10^{-4}$
$10^{5}$ | 48 | 126.4185 |  1.3908$\times 10^{-2}$ |  2.5739$\times 10^{-3}$ |  8.6596$\times 10^{-4}$
$10^{6}$ | 48 | 128.2346 |  1.3912$\times 10^{-2}$ |  2.5739$\times 10^{-3}$ |  8.5623$\times 10^{-4}$
$10^{7}$ | 48 | 127.1858 |  1.3912$\times 10^{-2}$ |  2.5739$\times 10^{-3}$ |  8.5612$\times 10^{-4}$
$10^{8}$ | 48 | 124.1160 |  1.3912$\times 10^{-2}$ |  2.5739$\times 10^{-3}$ |  8.5610$\times 10^{-4}$

Example 2. As in Example 1, we consider the bilinear optimal control problem
(BCP) on the domain $Q=\Omega\times(0,T)$ with $\Omega=(0,1)^{2}$ and $T=1$.
Now, we take the control $\bm{v}(x,t)$ in the infinite-dimensional space
$\mathcal{U}=\\{\bm{v}|\bm{v}\in[L^{2}(Q)]^{2},\nabla\cdot\bm{v}=0\\}.$ We set
$\alpha_{2}=0$ in (1.1), $\nu=1$ and $a_{0}=1$ in (1.2), and consider the
following tracking-type bilinear optimal control problem:

(5.3)
$\min_{\bm{v}\in\mathcal{U}}J(\bm{v})=\frac{1}{2}\iint_{Q}|\bm{v}|^{2}dxdt+\frac{\alpha_{1}}{2}\iint_{Q}|y-y_{d}|^{2}dxdt,$

where $y$ is obtained from $\bm{v}$ via the solution of the state equation
(1.2).

First, we let

$\displaystyle y=e^{t}(-3\sin(2\pi x_{1})\sin(\pi x_{2})+1.5\sin(\pi
x_{1})\sin(2\pi x_{2})),$ $\displaystyle p=(T-t)\sin\pi x_{1}\sin\pi
x_{2},~{}\text{and}~{}\bm{u}=P_{\mathcal{U}}(p\nabla y),$

where $P_{\mathcal{U}}(\cdot)$ is the projection onto the set $\mathcal{U}$.

We further set

$\displaystyle f=\frac{\partial y}{\partial t}-\nabla^{2}y+{\bm{u}}\cdot\nabla
y+y,\quad\phi=-3\sin(2\pi x_{1})\sin(\pi x_{2})+1.5\sin(\pi x_{1})\sin(2\pi
x_{2}),$ $\displaystyle y_{d}=y-\frac{1}{\alpha_{1}}\left(-\frac{\partial
p}{\partial t}-\nabla^{2}p-\bm{u}\cdot\nabla p+p\right),\quad g=0.$

Then, it is easy to show that $\bm{u}$ is a solution point of the problem
(5.3). We note that $\bm{u}=P_{\mathcal{U}}(p\nabla y)$ has no analytical
solution and it can only be solved numerically. Here, we solve
$\bm{u}=P_{\mathcal{U}}(p\nabla y)$ by the preconditioned CG algorithm
(DG1)–(DG5) with $h=\frac{1}{2^{9}}$ and $\Delta t=\frac{1}{2^{10}}$, and use
the resulting control $\bm{u}$ as a reference solution for the example we
considered.

Figure 7. The target function $y_{d}$ with $h=\frac{1}{2^{7}}$ and $\Delta
t=\frac{1}{2^{8}}$ at $t=0.25,0.5$ and $0.75$ (from left to right) for Example
2.

The stopping criteria of the outer CG algorithm (DI)–(DV) and the inner
preconditioned CG algorithm (DG1)–(DG5) are respectively set as

$\frac{\Delta t\sum_{n=1}^{N}\int_{\Omega}|\bm{g}_{n}^{k+1}|^{2}dx}{\Delta
t\sum_{n=1}^{N}\int_{\Omega}|\bm{g}_{n}^{0}|^{2}dx}\leq 5\times
10^{-8},~{}\text{and}~{}\frac{\int_{\Omega}|\nabla
r^{k+1}|^{2}dx}{\max\\{1,\int_{\Omega}|\nabla r^{0}|^{2}dx\\}}\leq 10^{-8}.$

The initial values are chosen as $\bm{u}^{0}=(0,0)^{\top}$ and
$\lambda^{0}=0$; and we denote by $\bm{u}_{h}^{\Delta t}$ and $y_{h}^{\Delta
t}$ the computed control and state, respectively.

First, we take $h=\frac{1}{2^{i}},i=6,7,8$, $\Delta t=\frac{h}{2}$,
$\alpha_{1}=10^{6}$, and implement the proposed nested CG algorithm (DI)–(DV)
for solving the problem (5.3). The numerical results reported in Table 3 show
that the CG algorithm converges fast and is robust with respect to different
mesh sizes. In addition, the preconditioned CG algorithm (DG1)–(DG5) converges
within 10 iterations for all cases and thus is efficient for computing the
gradient $\\{\bm{g}_{n}\\}_{n=1}^{N}$. We also observe that the target
function $y_{d}$ has been reached within a good accuracy. Similar comments
hold for the approximation of the optimal control $\bm{u}$ and of the state
$y$ of problem (5.3).

Table 3. Results of the nested CG algorithm (DI)–(DV) with different $h$ and $\Delta t$ for Example 2. Mesh sizes | $Iter_{CG}$ | $MaxIter_{PCG}$ | $\|\bm{u}_{h}^{\Delta t}-\bm{u}\|_{L^{2}(Q)}$ | $\|y_{h}^{\Delta t}-y\|_{L^{2}(Q)}$ | $\frac{\|y_{h}^{\Delta t}-y_{d}\|_{L^{2}(Q)}}{\|y_{d}\|_{{L^{2}(Q)}}}$
---|---|---|---|---|---
$h=1/2^{6},\Delta t=1/2^{7}$ | 443 | 9 |  3.7450$\times 10^{-3}$ |  9.7930$\times 10^{-5}$ |  1.0906$\times 10^{-6}$
$h=1/2^{7},\Delta t=1/2^{8}$ | 410 | 9 |  1.8990$\times 10^{-3}$ |  1.7423$\times 10^{-5}$ |  3.3863$\times 10^{-7}$
$h=1/2^{8},\Delta t=1/2^{9}$ | 405 | 8 |  1.1223$\times 10^{-3}$ |  4.4003$\times 10^{-6}$ |  1.0378$\times 10^{-7}$

Taking $h=\frac{1}{2^{7}}$ and $\Delta t=\frac{1}{2^{8}}$, the computed state
$y_{h}^{\Delta t}$, the error $y_{h}^{\Delta t}-y$ and $y_{h}^{\Delta
t}-y_{d}$ at $t=0.25,0.5,0.75$ are reported in Figures 8, 9 and 10,
respectively; and the computed control $\bm{u}_{h}^{\Delta t}$, the exact
control $\bm{u}$, and the error $\bm{u}_{h}^{\Delta t}-\bm{u}$ at
$t=0.25,0.5,0.75$ are presented in Figures 11, 12 and 13.

Figure 8. Computed state $y^{\Delta t}_{h}$, error $y^{\Delta t}_{h}-y$ and
$y^{\Delta t}_{h}-y_{d}$ with $h=\frac{1}{2^{7}}$ and $\Delta
t=\frac{1}{2^{8}}$ (from left to right) at $t=0.25$ for Example 2.

Figure 9. Computed state $y^{\Delta t}_{h}$, error $y^{\Delta t}_{h}-y$ and
$y^{\Delta t}_{h}-y_{d}$ with $h=\frac{1}{2^{7}}$ and $\Delta
t=\frac{1}{2^{8}}$ (from left to right) at $t=0.5$ for Example 2.

Figure 10. Computed state $y^{\Delta t}_{h}$, error $y^{\Delta t}_{h}-y$ and
$y^{\Delta t}_{h}-y_{d}$ with $h=\frac{1}{2^{7}}$ and $\Delta
t=\frac{1}{2^{8}}$ (from left to right) at $t=0.75$ for Example 2.

Figure 11. Computed control $\bm{u}^{\Delta t}_{h}$ and exact control $\bm{u}$
(left, from top to bottom) and the error $\bm{u}^{\Delta t}_{h}-\bm{u}$
(right) with $h=\frac{1}{2^{7}}$ and $\Delta t=\frac{1}{2^{8}}$ at $t=0.25$
for Example 2.

Figure 12. Computed control $\bm{u}^{\Delta t}_{h}$ and exact control $\bm{u}$
(left, from top to bottom) and the error $\bm{u}^{\Delta t}_{h}-\bm{u}$
(right) with $h=\frac{1}{2^{7}}$ and $\Delta t=\frac{1}{2^{8}}$ at $t=0.5$ for
Example 2.

Figure 13. Computed control $\bm{u}^{\Delta t}_{h}$ and exact control $\bm{u}$
(left, from top to bottom) and the error $\bm{u}^{\Delta t}_{h}-\bm{u}$
(right) with $h=\frac{1}{2^{7}}$ and $\Delta t=\frac{1}{2^{8}}$ at $t=0.75$
for Example 2.

##  6\. Conclusion and Outlook

We studied the bilinear control of an advection-reaction-diffusion system,
where the control variable enters the model as a velocity field of the
advection term. Mathematically, we proved the existence of optimal controls
and derived the associated first-order optimality conditions. Computationally,
the conjugate gradient (CG) method was suggested and its implementation is
nontrivial. In particular, an additional divergence-free constraint on the
control variable leads to a projection subproblem to compute the gradient; and
the computation of a stepsize at each CG iteration requires solving the state
equation repeatedly due to the nonlinear relation between the state and
control variables. To resolve the above issues, we reformulated the gradient
computation as a Stokes-type problem and proposed a fast preconditioned CG
method to solve it. We also proposed an efficient inexactness strategy to
determine the stepsize, which only requires the solution of one linear
parabolic equation. An easily implementable nested CG method was thus
proposed. For the numerical discretization, we employed the standard piecewise
linear finite element method and the Bercovier-Pironneau finite element method
for the space discretizations of the bilinear optimal control and the Stokes-
type problem, respectively, and a semi-implicit finite difference method for
the time discretization. The resulting algorithm was shown to be numerically
efficient by some preliminary numerical experiments.

We focused in this paper on an advection-reaction-diffusion system controlled
by a general form velocity field. In a real physical system, the velocity
field may be determined by some partial differential equations (PDEs), such as
the Navier-Stokes equations. As a result, we meet some bilinear optimal
control problems constrained by coupled PDE systems. Moreover, instead of
(1.1), one can also consider other types of objective functionals in the
bilinear optimal control of an advection-reaction-diffusion system. For
instance, one can incorporate $\iint_{Q}|\nabla\bm{v}|^{2}dxdt$ and
$\iint_{Q}|\frac{\partial\bm{v}}{\partial t}|^{2}dxdt$ into the objective
functional to promote that the optimal velocity field has the least rotation
and is almost steady, respectively, which are essential in e.g., mixing
enhancement for different flows [24]. All these problems are of practical
interest but more challenging from algorithmic design perspectives, and they
have not been well-addressed numerically in the literature. Our current work
has laid a solid foundation for solving these problems and we leave them in
the future.

## References

  * [1] R. Becker and B. Vexler, Optimal control of the convection-diffusion equation using stabilized finite element methods, Numerische Mathematik, 106 (2007), pp. 349–367.
  * [2] M. Bercovier and O. Pironneau, Error estimates for finite element method solution of the Stokes problem in the primitive variables. Numerische Mathematik, 33 (1979), pp. 211–224.
  * [3] A. Borzì, E.-J. Park and M. Vallejos Lass, Multigrid optimization methods for the optimal control of convection-diffusion problems with bilinear control, Journal of Optimization Theory and Applications, 168 (2016), pp. 510–533.
  * [4] P. Cannarsa, G. Floridia and A. Y. Khapalov, Multiplicative controllability for semilinear reaction–diffusion equations with finitely many changes of sign, Journal de Mathématiques Pures et Appliquées 108 (2017), pp. 425–458.
  * [5] C. Carthel, R. Glowinski and J. L. Lions, On exact and approximate boundary controllabilities for the heat equation: a numerical approach, Journal of Optimization Theory and Applications, 82 (1994), pp. 429–484.
  * [6] L. Dede’ and A. Quarteroni, Optimal control and numerical adaptivity for advection-diffusion equations, ESAIM: Mathematical Modelling and Numerical Analysis, 39 (2005), pp. 1019–1040.
  * [7] A. Fleig and R. Guglielmi, Optimal control of the Fokker–Planck equation with space-dependent controls, Journal of Optimization Theory and Applications, 174 (2017), pp. 408–427.
  * [8] R. Glowinski, Ensuring well-posedness by analogy; Stokes problem and boundary control for the wave equation, Journal of Computational Physics, 103 (1992), pp. 189–221.
  * [9] R. Glowinski, Finite Element Methods for Incompressible Viscous Flow, Handbook of Numerical Analysis, Vol. 9, Elsevier, Amsterdam, 2003, pp. 3-1176.
  * [10] R. Glowinski, Variational Methods for the Numerical Solution of Nonlinear Elliptic Problems, Society for Industrial and Applied Mathematics, Philadelphia, 2015.
  * [11] R. Glowinski and J. He, On shape optimization and related issues, In Computational Methods for Optimal Design and Control, J. Borggaard, J. Burns, E. Cliff & S. Schreck (eds.), Birkhäuser, Boston, MA, 1998, pp. 151–179.
  * [12] R. Glowinski and J. L. Lions, Exact and approximate controllability for distributed parameter systems, Part I, Acta Numerica, 3 (1994), pp. 269–378.
  * [13] R. Glowinski and J. L. Lions, Exact and approximate controllability for distributed parameter systems, Part II, Acta Numerica, 4 (1995), pp. 159–328.
  * [14] R. Glowinski, J. L. Lions and J. He, Exact and Approximate Controllability for Distributed Parameter Systems: A Numerical Approach (Encyclopedia of Mathematics and its Applications), Cambridge University Press, 2008.
  * [15] N. Handagama and S. Lenhart, Optimal control of a PDE/ODE system modeling a gas-phase bioreactor, In Mathematical Models in Medical and Health Sciences, M. A. Horn, G. Simonett, and G. Webb (eds.), Vanderbilt University Press, Nashville, TN, 1998.
  * [16] K. Ito and K. Kunisch, Optimal bilinear control of an abstract Schrödinger equation, SIAM Journal on Control and Optimization, 46 (2007), pp. 274–287.
  * [17] H. R. Joshi, Optimal control of the convective velocity coefficient in a parabolic problem, Nonlinear Analysis: Theory, Methods & Applications, 63 (2005), pp. e1383–e1390.
  * [18] A. Y. Khapalov, Controllability of the semilinear parabolic equation governed by a multiplicative control in the reaction term: a qualitative approach, SIAM Journal on Control and Optimization, 41 (2003), pp. 1886–1900.
  * [19] A. Y. Khapalov, Controllability of Partial Differential Equations Governed by Multiplicative Controls, Springer, 2010.
  * [20] A. Kröner and B. Vexler, A priori error estimates for elliptic optimal control problems with a bilinear state equation, Journal of Computational and Applied Mathematics, 230 (2009), pp. 781–802.
  * [21] S. Lenhart, Optimal control of a convective-diffusive fluid problem, Mathematical Models and Methods in Applied Sciences, 5 (1995), pp. 225–237.
  * [22] J. L. Lions, Optimal Control of Systems Governed by Partial Differential Equations (Grundlehren der Mathematischen Wissenschaften), Vol. 170, Springer Berlin, 1971.
  * [23] J. L. Lions, Exact controllability, stabilization and perturbations for distributed systems, SIAM Review, 30 (1988), pp. 1–68.
  * [24] W. Liu. Mixing enhancement by optimal flow advection, SIAM Journal on Control and Optimization, 47 (2008), pp. 624–638.
  * [25] J. Nocedal, and S.J. Wright, Numerical Optimization, Second Edition, Springer, 2006.
  * [26] F. Tröltzsch, Optimal Control of Partial Differential Equations: Theory, Methods, and Applications, Vol. 112, American Mathematical Society, 2010.
  * [27] E. Zuazua, Propagation, observation, and control of waves approximated by finite difference methods, SIAM Review, 47 (2005), pp. 197–243.
  * [28] E. Zuazua, Controllability of Partial Differential Equations, 3rd cycle, Castro Urdiales (Espagne), 2006, pp.311. cel-00392196.
  * [29] E. Zuazua, Controllability and observability of partial differential equations: some results and open problems, Handbook of differential equations: evolutionary equations, Vol. 3, North-Holland, 2007, pp. 527–621.