# Parallelized POD-based Suboptimal Economic Model Predictive Control of a
State-Constrained Boussinesq approximation

Julian Andrej111Lawrence Livermore National Laboratory, 7000 East Ave.,
Livermore, 94550, USA. Lars Grüne222University of Bayreuth, Universitätsstraße
30, 95447 Bayreuth, Germany. Luca Mechelli333University of Konstanz,
Universitätsstraße 10, 78464 Konstanz, Germany. Thomas Meurer444Chair of
Automatic Control, Kiel University, Kaiserstraße 2, 24143, Kiel, Germany.
Simon Pirkelmann† and Stefan Volkwein‡

###### Abstract

Motivated by an energy efficient building application, we want to optimize a
quadratic cost functional subject to the Boussinesq approximation of the
Navier-Stokes equations and to bilateral state and control constraints. Since
the computation of such an optimal solution is numerically costly, we design
an efficient strategy to compute a sub-optimal (but applicationally
acceptable) solution with signifcantly reduced computational effort. We employ
an economic Model Predictive Control (MPC) strategy to obtain a feedback
control. The MPC sub-problems are based on a linear-quadratic optimal control
problem subjected to mixed control and state constraints and a convection-
diffusion equation, reduced with proper orthogonal decomposition. To solve
each sub-problem, we apply a primal-dual active set strategy. The method can
be fully parallelized, which enables the solution of large problems with real-
world parameters.

Keywords: Boussinesq approximation, Model predictive control, State
constraints, Proper orthogonal decomposition, Parallel computing.

##  1 Introduction

Coupled heat-convection phenomena like those occurring in heating, cooling and
air-conditioning (HVAC) of residual buildings can be accurately modeled by the
incompressible Navier-Stokes equations; see [2, 38, 40], for instance. These
equations describe how the physical quantities temperature, air velocities and
pressure are connected and how they influence each other. While the underlying
physical principles are completely described by the equations, working with
them directly for the purpose of control of HVAC processes has turned out to
be challenging due to the high computational complexity. Instead, simplified
models of varying granularity are used in control applications.

One common simplification that still retains much of the relevant physical
relations in the HVAC regime is the Boussinesq approximation [38, 40]. The
general idea is to reduce the coupling between heat and airflow by only
considering bouyancy induced changes of the fluid velocity. More specifically,
we only consider density variations along the gravity field which lead to,
e.g., the rise of warmer air. Although we obtain a simplified model, this
simplification is not enough to prevent high computational cost for computing
its solution. Many degrees of freedom (DOFs) are in fact necessary to
accurately capture the involved dynamics, especially when dealing with real-
world parameters; see, e.g., [5].

One viable approach is to apply model order reduction (MOR). This field was
extensively developed during the last two decades and several techniques were
proposed. Among all, in case of computational fluid dynamics, proper
orthogonal decomposition (POD) seems the one to guarantee best performances
for nonlinear dynamical systems; see, e.g., [6, 24, 37]. In the particular
case of the Boussinesq approximation, the POD-based MOR approach is also
certified by an a-posteriori error analysis [32, 33], that guarantees the
validity of approximation for the computed reduced order model. Due to the
request of computing many PDE solutions, MOR techniques are naturally applied
in the field of optimal control. In [22], for example, a POD-based model
predictive control (MPC) algorithm is considered to approximate a feedback
solution to the Boussinesq equation.

MPC is a well-established model-based control method, sometimes also called
Receding Horizon Control; cf. [16, 34]. The key idea is to decompose an
infinite time horizon optimal control problem into several finite horizon sub-
problems and, therefore, to optimize predictions of the model behavior in
relatively short time. After a sub-problem is solved, the initial part of the
computed optimal control is stored and used as feedback control, the time
horizon is shifted and it is possible to update the problem parameters before
the procedure is iterated. This technique is particularly useful in HVAC
applications. It permits, in fact, to merge predictions, such as weather
forecast, with system measurements [30, 43]. Let us mention that POD-based MPC
is considered also in [10, 12, 27, 29] for different models. Furthermore, we
are particularly interested in economic MPC; cf. [4, 9, 11, 17, 19, 31]. In
this particular branch of MPC, more general cost functionals than those
employing tracking type costs are allowed. In our specific application, the
cost functional depends only on the control variables, without any dependence
on the state variable. The focus is then to keep the system trajectory inside
pre-defined bounds with less control effort as possible and it deals perfectly
with the usual requests of HVAC applications.

This work aspires to unify the above mentioned contributions in a reliable
method, in order to compute a suboptimal solution to an optimal control
problem subject to the Boussinesq approximation and bilateral state and
control constraints. In particular, starting from [5, 27, 31], we propose a
POD-based economic MPC scheme which considers linear-quadratic subproblems
governed by convection-diffusion equations with a virtual control approach
[26, 27] and possible updates of the reduced order model through parallel
solution of the Boussinesq approximation [5]. Furthermore a primal-dual active
set strategy (PDASS) [21, 25] is considered to solve each sub-problem. The
goal is to have a robust and reliable method that speeds up the computational
time.

In Section 2, we present the optimal control problem. Section 3 describes the
virtual control approach and the PDASS for linear-quadratic optimal control
problems, while Section 4 contains a brief explanation of the POD method and
the a-posteriori error estimate. The POD-based economic MPC algorithm is
explained in Section 5 and further details about the implementation are given
in Section 6. Finally, we test our algorithm for a simplified problem in
Section 7.

##  2 The setting of the optimal control problem

The Boussinesq approximation of the Navier-Stokes equations is described by
the following PDE:

$\displaystyle\frac{\partial{\mathbf{v}}}{\partial
t}+{\mathbf{v}}\cdot\nabla{\mathbf{v}}$ $\displaystyle=-\frac{1}{\rho}\nabla
p+\nu\Delta{\mathbf{v}}-\textbf{g}\alpha({y}-\tilde{y})$
$\displaystyle\text{in }Q:=(0,T)\times\Omega,$ (1a)
$\displaystyle\nabla\cdot{\mathbf{v}}$ $\displaystyle=0$
$\displaystyle\text{in }Q,$ (1b) $\displaystyle\mathbf{v}(0)$
$\displaystyle=\mathbf{v}_{\circ},$ $\displaystyle\text{in }\Omega,$ (1c)
$\displaystyle\frac{\partial y}{\partial t}+{\mathbf{v}}\cdot\nabla y$
$\displaystyle=\frac{\kappa}{\rho c_{p}}\,\Delta y$ $\displaystyle\text{in
}Q,$ (1d) $\displaystyle y(0)$ $\displaystyle=y_{\circ}$
$\displaystyle\text{in }\Omega,$ (1e) where $\mathbf{v}:Q\to\mathbb{R}^{2}$ is
the _air velocity_ , $p:Q\to\mathbb{R}$ stands for the _air pressure_ , and
$y:Q\to\mathbb{R}$ is the _temperature_ , $\Omega\subset\mathbb{R}^{2}$
denotes an open and bounded domain with sufficiently smooth boundary, $T>0$ is
the time horizon, $\textbf{g}\in\mathbb{R}^{2}$ is the _gravitational
acceleration_ , $\mathbf{v}_{\circ}$, $y_{\circ}$ are given initial conditions
and $\tilde{y}\in\mathbb{R}$ is the reference temperature at which we measure:

  * •

the (constant) density $\rho>0$,

  * •

the (constant) kinematic viscosity $\nu>0$,

  * •

the (constant) coefficient of thermal expansion $\alpha>0$,

  * •

the (constant) thermal conductivity $\kappa>0$,

  * •

the (constant) isobaric specific heat $c_{p}>0$

of the studied fluid (air in our case).

Figure 1: Domain $\Omega$ in case of underfloor heating system.

In addition, we consider time-varying boundary conditions of the following
form

$\displaystyle v$ $\displaystyle=0$ $\displaystyle\text{ on
}(0,T)\times\Gamma,$ (1f) $\displaystyle-\kappa\frac{\partial
y}{\partial{\bm{n}}}$ $\displaystyle=\gamma(y-y_{\text{out}})$
$\displaystyle\text{ on }(0,T)\times\Gamma\setminus\Gamma_{\mathsf{c}},$ (1g)
$\displaystyle-\kappa\frac{\partial y}{\partial{\bm{n}}}$
$\displaystyle=\gamma_{\mathsf{c}}(y-u)$ $\displaystyle\text{ on
}(0,T)\times\Gamma_{\mathsf{c}},$ (1h)

where $\Gamma=\partial\Omega$ and $\gamma,\gamma_{c}>0$ are fixed parameters.
The function $u\in\mathcal{U}=L^{2}(0,T)$ is a control input placed at a
predefined part of the domain’s boundary, namely $\Gamma_{\mathsf{c}}$, and it
represents a device for controlling the temperature in the room, e.g. a
radiator or an underfloor heating system; cf. Figure 1. The quantity
$y_{\text{out}}\in L^{2}(0,T)$ represents the outside temperature.

We are interested in the weak solution to (1). For more details of the
straightforward variational formulation for (1) we refer to [2, 38] and [22]
for more details.

The goal is to minimize the control effort while keeping the temperature $y$
inside desired bounds, i.e. we are interested in minimizing

$\mathcal{J}(u)=\frac{1}{2}\|u\|^{2}_{\mathcal{U}}$ (2a) subject to (s.t.) the
weak formulation of (1) and the inequality constraints $\displaystyle
u_{\mathsf{a}}(t)\leq u(t)\leq u_{\mathsf{b}}(t)$ $\displaystyle\text{ a.e. in
}(0,T),$ (2b) $\displaystyle y_{\mathsf{a}}(t,{\bm{x}})\leq y(t,{\bm{x}})\leq
y_{\mathsf{b}}(t,{\bm{x}})$ $\displaystyle\text{ a.e. in }Q.$ (2c)

Here $y$ is the (weak) solution of (1),
$u_{\mathsf{a}},u_{\mathsf{b}}\in\mathcal{U}$ given control constraints and
$y_{\mathsf{a}},y_{\mathsf{b}}\in C(Q)$ given state constraints. It follows
from [38, Chap. III, §3] that (1) admits a unique weak solution

$\displaystyle(\mathbf{v},y)\in\mathcal{Y}_{\mathsf{B}}=W(0,T;\tilde{V})\times
W(0,T;V),$

where $V=H^{1}(\Omega)$ and the Hilbert space $\tilde{V}$ is the closure of
the space $\\{v\in C^{\infty}_{0}(\Omega)^{2}\,|\,\nabla\cdot v=0\\}$ with
respect to the norm $\|\cdot\|_{V\times V}$. Moreover, we have
$W(0,T;\tilde{V})=L^{2}(0,T;\tilde{V})\cap H^{1}(0,T;\tilde{V}^{\prime})$, and
the set $W(0,T;V)$ is defined analogously; see [8, Chap. XVIII] for more
details. Thus, considering the control-to-state map

$\displaystyle\mathcal{S}_{\mathsf{B}}:\mathcal{U}\to\mathcal{Y}_{\mathsf{B}},\quad
u\mapsto\mathcal{S}_{\mathsf{B}}(u)=(\textbf{v}(u),y(u))\in\mathcal{Y}_{\mathsf{B}}$

and the admissible set

$\displaystyle\mathcal{U}_{\mathsf{B},\mathsf{ad}}=\big{\\{}u\in\mathcal{U}\,\big{|}\,u\text{
and }y=\mathcal{S}_{\mathsf{B}}(u)\text{ satisfy \eqref{control_constraints}
and \eqref{state_constraints}, respectively}\big{\\}},$

problem (2) can be rewritten as the following purely control-constrained
optimization problem

$\min\mathcal{J}(u)\quad\text{s.t.}\quad
u\in\mathcal{U}_{\mathsf{B},\mathsf{ad}}.$ ($\mathbf{\hat{P}}_{\mathsf{B}}$)

We call ($\mathbf{\hat{P}}_{\mathsf{B}}$) the _reduced problem_ because – in
contrast to (2) – it is an optimization problem in the control variable only.
Note that ($\mathbf{\hat{P}}_{\mathsf{B}}$) admits a global solution
$\bar{u}_{\mathsf{B}}\in\mathcal{U}_{\mathsf{B},\mathsf{ad}}$ provided
$\mathcal{U}_{\mathsf{B},\mathsf{ad}}$ has non-empty interior; cf., e.g.,
[23]. Throughout the work a bar indicates optimality.

While the Boussinesq approximation already reduces numerical effort for
simulation of the fluid processes compared to the full compressible Navier-
Stokes equations, the nonlinearity induced by the two-way coupling of
temperature and velocity still makes it computationally infeasible, when it is
used, e.g., as the prediction model in an MPC scheme. This is because each
step of the MPC algorithm typically requires many evaluations of the
prediction model, which would be too costly for an online application of the
method.

For that reason we elaborate a simplified strategy that aims at decreasing the
computational cost significantly. We introduce several simplifications from a
model and optimal control perspectives. At first, we notice that the buoyancy
effect in the Boussinesq model has a big impact only for large time horizons,
therefore in the MPC open-loop problem, one can consider to fix the current
velocity field $\mathbf{v}$ and solve only the linear convection-diffusion
equation (1d).

However, we will still carry out occasional simulations of the Boussinesq
approximation in order to generate the source for so-called snapshots of the
system state. From those we can then construct a local (with respect to time)
reduced-order model of the convection-diffusion equation (1d) system using
POD. This reduced model summarizes information about the essential dynamics
currently going on in the system and is much smaller than the original system
which allows us to use it for fast computation of an approximated solution of
optimal control problems within the MPC algorithm.

One may ask, why not performing POD directly for the coupled system (1), as in
[32, 33]. Although this is a viable approach, the use of the velocity field v
as datum of the problem and not as unknown improves and simplifies further the
POD approximation. This has the double advantage of using less basis functions
and easier (mathematically and computationally) a-posteriori error estimator
for the problem, which may result in tighter estimates compared to the one for
the fully coupled system.

The last challenge is caused by the state constraints. These do not guarantee
the existence of regular ($L^{2}$) Lagrange multipliers (cf. e.g. [42]) ,
introducing additional computational effort. To overcome this issue, we employ
a virtual control approach together with a primal-dual active set strategy.

##  3 Virtual control approach and PDASS for the linear-quadratic problem

As mentioned in Section 2, as first simplification we assume that the velocity
field v is given and we focus on a model based only on the heat equation
(1d)-(1e) with convection together with the boundary conditions (1g)-(1h). It
is well-known that this equation admits a unique weak solution
$y\in\mathcal{Y}=W(0,T;V)$ for any given control $u\in\mathcal{U}$; cf., e.g.,
[8]. Let $\mathcal{S}:\mathcal{U}\to\mathcal{Y}$ be the control-to-state map
for the heat equation, then we define the admissible set of controls

$\mathcal{U}_{\mathsf{ad}}=\big{\\{}u\in\mathcal{U}\,\big{|}\,u\text{ and
}y=\mathcal{S}(u)\text{ satisfy \eqref{control_constraints} and
\eqref{state_constraints}, respectively}\big{\\}}.$

We are, thus, interested in solving the strictly convex optimal control
problem

$\min\mathcal{J}(u)\quad\text{s.t.}\quad u\in\mathcal{U}_{\mathsf{ad}},$
($\mathbf{\hat{P}}$)

which admits a unique solution $\bar{u}$ if the set
$\mathcal{U}_{\mathsf{ad}}$ has non-empty interior; see, e.g., [23].

Note that the difference between ($\mathbf{\hat{P}}_{\mathsf{B}}$) and
($\mathbf{\hat{P}}$) relies on the different control-to-state map, i.e. on the
different underlying state model. As said, we expect that
$\bar{u}_{\mathsf{B}}$ and $\bar{u}$ will be close for short time horizon $T$,
since the buoyancy effect influences the temperature trajectory $y$ only on
larger time scales. This will allow computing a good enough approximation of
the velocity field v for (1d), to be used as datum of the problem. We will
explain in Section 6 how we compute this approximation. What is important to
clarify now is how to compute the solution to ($\mathbf{\hat{P}}$). As
mentioned, the presence of pointwise state constraints leads to non-regular
Lagrange multipliers; cf. [7, 35, 36]. To overcome this issue, we apply a
virtual control approach as in [26, 27, 28]. Let $w\in\mathcal{W}=L^{2}(Q)$ be
an additional (artificial) control. For given $\varepsilon>0$ we relax the
bilateral state constraints (2c) as follows:

$y_{\mathsf{a}}\leq\mathcal{S}(u)+\varepsilon w\leq
y_{\mathsf{b}}\quad\text{a.e. in }Q$ (3)

and introduce the artificial constraints

$w_{\mathsf{a}}\leq w(t,{\bm{x}})\leq w_{\mathsf{b}}\quad\text{a.e. in }Q$ (4)

with fixed (safeguard) scalars $w_{\mathsf{a}}\ll 0$ and $w_{\mathsf{b}}\gg
0$, i.e., sufficiently small and large, respectively. Furthermore, we define
the admissible set of controls

$\mathcal{Z}^{\varepsilon}_{\mathsf{ad}}=\big{\\{}(u,w)\in\mathcal{U}\times\mathcal{W}\,\big{|}\,u,w\text{
and }\mathcal{S}(u)\text{ satisfy
}\eqref{control_constraints},\eqref{relaxed_state_constraints},\eqref{artificial_constraints}\big{\\}}.$

Now we are interested in solving

$\min
J(u,w)\quad\text{s.t.}\quad(u,w)\in\mathcal{Z}^{\varepsilon}_{\mathsf{ad}},$
($\mathbf{\hat{P}}^{\varepsilon}$)

where

$J(u,w)=\mathcal{J}(u)+\frac{\sigma}{2}\,{\|w\|}^{2}_{\mathcal{W}}$

with $\sigma>0$. Note that ($\mathbf{\hat{P}}^{\varepsilon}$) admits a unique
solution
$(\bar{u}^{\varepsilon},\bar{w}^{\varepsilon})\in\mathcal{U}\times\mathcal{W}$
for each $\varepsilon>0$ provided $\mathcal{Z}^{\varepsilon}_{\mathsf{ad}}$
has non-empty interior (cf. [23, 27]).

######  Remark 1.

We point out that the artificial bounds $w_{\mathsf{a}}$, $w_{\mathsf{b}}$ are
necessary to guarantee $\bar{u}^{\varepsilon}\to\bar{u}$ in $\mathcal{U}$ as
$\varepsilon\to 0$ with an order of convergence $O(\sqrt{\varepsilon})$; cf.
[27, Section 1.3.3]. They can be omitted if $\varepsilon$ is taken small but
fixed. In this case, the existence of a unique solution is still guaranteed
[28]. Moreover, we can choose $w_{\mathsf{a}}$ and $w_{\mathsf{b}}$ as large
in absolute value as one desires, because they do not have any physical
meaning [27]. $\Diamond$

Let us define the non-symmetric continuous weakly-coercive bilinear form
$a(t;\cdot\,,\cdot):V\times V\to\mathbb{R}$ as

$a(t;\varphi,\psi)=\int_{\Omega}{\textstyle\frac{\kappa}{\rho
c_{p}}}\nabla\varphi\cdot\nabla\psi+(\textbf{v}(t,\cdot)\cdot\nabla\varphi)\psi\,\mathrm{d}{\bm{x}}+\gamma\int_{\Gamma\setminus\Gamma_{c}}\varphi\psi\,\mathrm{d}{\bm{s}}+\gamma_{c}\int_{\Gamma_{c}}\varphi\psi\,\mathrm{d}{\bm{s}}.$

Then the first-order necessary optimality conditions for
($\mathbf{\hat{P}}^{\varepsilon}$) are as follows (cf. [27]):

######  Theorem 1.

Let $\mathcal{Z}^{\varepsilon}_{\mathsf{ad}}$ have a non-empty interior.
Suppose that the pair
$(\bar{u}^{\varepsilon},\bar{w}^{\varepsilon})\in\mathcal{Z}^{\varepsilon}_{\mathsf{ad}}$
is the solution to ($\mathbf{\hat{P}}^{\varepsilon}$) with associated optimal
state $\bar{y}^{\varepsilon}=\mathcal{S}(\bar{u}^{\varepsilon})$, i.e.
satisfying

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\,{\langle\bar{y}^{\varepsilon}(t),\varphi\rangle}_{H}+a(t;\bar{y}^{\varepsilon}(t),\varphi)$
$\displaystyle=\gamma
y_{\text{out}}(t)\int_{\Gamma\setminus\Gamma_{c}}\varphi\,\mathrm{d}{\bm{s}}+\gamma_{c}u(t)\int_{\Gamma_{c}}\varphi\,\mathrm{d}{\bm{s}},$
(5a) $\displaystyle\bar{y}^{\varepsilon}(0)$ $\displaystyle=y_{\circ}$ for all
$\varphi\in V$ and a.e. in $[0,T]$, where $H=L^{2}(\Omega)$. Then, there exist
unique Lagrange multipliers $\bar{p}^{\varepsilon}\in\mathcal{Y}$,
$\bar{\beta}^{\varepsilon},\bar{\vartheta}^{\varepsilon}\in\mathcal{W}$ and
$\bar{\alpha}^{\varepsilon}\in\mathcal{U}$ satisfying the dual equation
$\displaystyle-\frac{\mathrm{d}}{\mathrm{d}t}\,{\langle\bar{p}^{\varepsilon}(t),\varphi\rangle}_{H}+a(t;\varphi,\bar{p}^{\varepsilon}(t))$
$\displaystyle=-{\langle\bar{\beta}^{\varepsilon}(t),\varphi\rangle}_{H},$
(5b) $\displaystyle\bar{p}^{\varepsilon}(T)$
$\displaystyle=-\bar{\beta}^{\varepsilon}(T)$ for all $\varphi\in V$ and a.e.
in $[0,T]$ and the optimality system
$\displaystyle\bar{u}^{\varepsilon}-\gamma_{c}\int_{\Gamma_{c}}\bar{p}^{\varepsilon}\mathrm{d}s+\bar{\alpha}^{\varepsilon}$
$\displaystyle=0$ $\displaystyle\text{in }\mathcal{U},$ (5c)
$\displaystyle\sigma\bar{w}^{\varepsilon}+\varepsilon\bar{\beta}^{\varepsilon}+\bar{\vartheta}^{\varepsilon}$
$\displaystyle=0$ $\displaystyle\text{in }\mathcal{W}.$ (5d) Moreover,
$\displaystyle\bar{\beta}^{\varepsilon}$
$\displaystyle=\max\big{\\{}0,\bar{\beta}^{\varepsilon}+\eta(\bar{y}^{\varepsilon}+\varepsilon\bar{w}^{\varepsilon}-y_{\mathsf{b}})\big{\\}}+\min\big{\\{}0,\bar{\beta}^{\varepsilon}+\eta(\bar{y}^{\varepsilon}+\varepsilon\bar{w}^{\varepsilon}-y_{\mathsf{a}}\big{\\}},$
(5e) $\displaystyle\bar{\alpha}^{\varepsilon}$
$\displaystyle=\max\big{\\{}0,\bar{\alpha}^{\varepsilon}+\eta_{u}(\bar{u}^{\varepsilon}-u_{\mathsf{b}})\big{\\}}+\min\big{\\{}0,\bar{\alpha}^{\varepsilon}+\eta_{u}(\bar{u}^{\varepsilon}-u_{\mathsf{b}})\big{\\}},$
(5f) $\displaystyle\bar{\vartheta}^{\varepsilon}$
$\displaystyle=\max\big{\\{}0,\bar{\vartheta}^{\varepsilon}+\eta_{w}(\bar{w}^{\varepsilon}-w_{\mathsf{b}})\big{\\}}+\min\big{\\{}0,\bar{\vartheta}^{\varepsilon}+\eta_{w}(\bar{w}^{\varepsilon}-w_{\mathsf{a}})\big{\\}}$
(5g)

for arbitrarily chosen $\eta,\eta_{u},\eta_{w}>0$, where the max- and min-
operations are interpreted componentwise in the pointwise everywhere sense.

######  Remark 2.

The well-defined bounded solution operator
$\mathcal{A}:\mathcal{W}\to\mathcal{Y}$ is defined as follows: for given
$\beta\in\mathcal{W}$ the function $p=\mathcal{A}\beta$ is the unique solution
to

$\displaystyle-\frac{\mathrm{d}}{\mathrm{d}t}\,{\langle
p(t),\varphi\rangle}_{H}+a(t;\varphi,p(t))$
$\displaystyle=-{\langle\beta(t),\varphi\rangle}_{H}$
$\displaystyle\forall\varphi\in V\text{ a.e. in }[0,T),$ $\displaystyle p(T)$
$\displaystyle=-\beta(T)$

for given $\beta\in\mathcal{W}$; cf. Remark 2.3 in [28]. Then,
$\bar{p}^{\varepsilon}=\mathcal{A}(\bar{\beta}^{\varepsilon})$ solves (5b).
$\Diamond$

Letting
$\nu=(u,w,\theta)\in\mathcal{N}=\mathcal{U}\times\mathcal{W}\times\mathcal{W}$,
we define the active sets corresponding to (5) as

$\displaystyle{\mathscr{A}_{{\mathsf{a}}}^{\mathcal{U}}}(\nu)$
$\displaystyle=\big{\\{}t\in(0,T)\,\big{|}\,\alpha(\nu)+u-u_{\mathsf{a}}<0\text{
a.e.}\big{\\}},$ (6a)
$\displaystyle{\mathscr{A}_{{\mathsf{b}}}^{\mathcal{U}}}(\nu)$
$\displaystyle=\big{\\{}t\in(0,T)\,\big{|}\,\alpha(\nu)+u-u_{\mathsf{b}}>0\text{
a.e.}\big{\\}},$ $\displaystyle{\mathscr{A}_{\mathsf{a}}^{\mathcal{W}}}(\nu)$
$\displaystyle=\bigg{\\{}(t,{\bm{x}})\in
Q\,\big{|}\,\beta(\nu)+\frac{\sigma}{\varepsilon^{2}}\big{(}y(\nu)+\varepsilon
w-y_{\mathsf{a}}\big{)}<0\text{ a.e.}\bigg{\\}},$
$\displaystyle{\mathscr{A}_{\mathsf{b}}^{\mathcal{W}}}(\nu)$
$\displaystyle=\bigg{\\{}(t,{\bm{x}})\in
Q\,\big{|}\,\beta(\nu)+\frac{\sigma}{\varepsilon^{2}}\big{(}y(\nu)+\varepsilon
w-y_{\mathsf{b}}\big{)}>0\text{ a.e.}\bigg{\\}},$
$\displaystyle{\mathscr{A}_{\mathsf{a},2}^{\mathcal{W}}}(\nu)$
$\displaystyle=\big{\\{}(t,{\bm{x}})\in Q\,\big{|}\,w-w_{\mathsf{a}}<0\text{
a.e.}\big{\\}},$
$\displaystyle{\mathscr{A}_{\mathsf{b},2}^{\mathcal{W}}}(\nu)$
$\displaystyle=\big{\\{}(t,{\bm{x}})\in Q\,\big{|}\,w-w_{\mathsf{b}}>0\text{
a.e.}\big{\\}}.$ Similarly, the associated inactive sets are
$\displaystyle{\mathscr{I}^{\mathcal{U}}}(\nu)$
$\displaystyle=(0,T)\setminus\big{(}{\mathscr{A}_{{\mathsf{a}}}^{\mathcal{U}}}(\nu)\cup{\mathscr{A}_{{\mathsf{b}}}^{\mathcal{U}}}(\nu)\big{)},$
(6b) $\displaystyle{\mathscr{I}^{\mathcal{W}}}(\nu)$
$\displaystyle=Q\setminus\big{(}{\mathscr{A}_{\mathsf{a}}^{\mathcal{W}}}(\nu)\cup{\mathscr{A}_{\mathsf{b}}^{\mathcal{W}}}(\nu)\big{)},\quad{\mathscr{I}^{\mathcal{W}}_{2}}(\nu)=Q\setminus\big{(}{\mathscr{A}_{\mathsf{a},2}^{\mathcal{W}}}(\nu)\cup{\mathscr{A}_{\mathsf{b},2}^{\mathcal{W}}}(\nu)\big{)}.$

One can solve (5) applying a primal-dual active set strategy or, equivalently,
a semi-smooth Newton method; see, e.g., [21, 25, 27]. Doing so, at each
iteration $k$ of the semi-smooth Newton method we have to solve the following
system

$\displaystyle\gamma_{c}\int_{\Gamma_{c}}p^{k+1}\,\mathrm{d}{\bm{s}}-u^{k+1}$
$\displaystyle=0$ $\displaystyle\text{in
}{\mathscr{I}^{\mathcal{U}}}(\nu^{k}),$ (7a) $\displaystyle u^{k+1}$
$\displaystyle=u_{\mathsf{a}}$ $\displaystyle\text{in
}{\mathscr{A}_{{\mathsf{a}}}^{\mathcal{U}}}(\nu^{k}),$ (7b) $\displaystyle
u^{k+1}$ $\displaystyle=u_{\mathsf{b}}$ $\displaystyle\text{in
}{\mathscr{A}_{{\mathsf{b}}}^{\mathcal{U}}}(\nu^{k}),$ (7c) $\displaystyle
w^{k+1}$ $\displaystyle=0$ $\displaystyle\text{in
}{\mathscr{I}^{\mathcal{W}}}(\nu^{k}),$ (7d) $\displaystyle
y^{k+1}+\varepsilon\,w^{k+1}$ $\displaystyle=y_{\mathsf{a}}$
$\displaystyle\text{in }{\mathscr{A}_{\mathsf{a}}^{\mathcal{W}}}(\nu^{k}),$
(7e) $\displaystyle y^{k+1}+\varepsilon\,w^{k+1}$
$\displaystyle=y_{\mathsf{b}}$ $\displaystyle\text{in
}{\mathscr{A}_{\mathsf{b}}^{\mathcal{W}}}(\nu^{k}),$ (7f)
$\displaystyle\theta^{k+1}$ $\displaystyle=0$ $\displaystyle\text{in
}{\mathscr{I}^{\mathcal{W}}_{2}}(\nu^{k}),$ (7g) $\displaystyle\sigma
w^{k+1}-\varepsilon\theta^{k+1}$ $\displaystyle=\sigma w_{\mathsf{a}}$
$\displaystyle\text{in }{\mathscr{A}_{\mathsf{a},2}^{\mathcal{W}}}(\nu^{k}),$
(7h) $\displaystyle\sigma w^{k+1}-\varepsilon\theta^{k+1}$
$\displaystyle=\sigma w_{\mathsf{b}}$ $\displaystyle\text{in
}{\mathscr{A}_{\mathsf{b},2}^{\mathcal{W}}}(\nu^{k}).$ (7i)

Taking into account (5), conditions (7) lead to the following coupled
equations:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}{\langle
y^{k+1}(t),\varphi\rangle}_{H}+a(y^{k+1}(t),\varphi)-\gamma_{c}\mathcal{H}^{k}(t;p^{k+1}(t))\int_{\Gamma_{c}}\varphi\,\mathrm{d}{\bm{s}}$
(8a) $\displaystyle\hskip 156.10345pt=\gamma
y_{\text{out}}(t)\int_{\Gamma\setminus\Gamma_{c}}\varphi\,\mathrm{d}{\bm{s}}+\gamma_{c}r^{k}(t)\int_{\Gamma_{c}}\varphi\,\mathrm{d}{\bm{s}},$
$\displaystyle-\frac{\mathrm{d}}{\mathrm{d}t}\,{\langle
p^{k+1}(t),\varphi\rangle}_{H}+a(t;\varphi,p^{k+1}(t))+\frac{\sigma}{\varepsilon^{2}}\,\left\langle\mathcal{G}^{k,\varepsilon}(t;y^{k+1}(t)),\varphi\right\rangle_{H}$
(8b) $\displaystyle\hskip 273.18271pt=\frac{\sigma}{\varepsilon^{2}}\,{\langle
r^{k,\varepsilon}(t),\varphi\rangle}_{H}$ for all $\varphi\in V$ and a.e. in
$(0,T)$, together with the conditions $\displaystyle y^{k+1}(0)$
$\displaystyle=y_{\circ},$ (8c) $\displaystyle
p^{k+1}(T)+\frac{\sigma}{\varepsilon^{2}}\,\mathcal{G}^{k,\varepsilon}(T;y^{k+1}(T))$
$\displaystyle=\frac{\sigma}{\varepsilon^{2}}\,r^{k,\varepsilon}(T),$ (8d)

where

$\displaystyle\mathcal{G}^{k,\varepsilon}(t;y^{k+1}(t))$
$\displaystyle=y^{k+1}(t)\sum_{i=1}^{6}\chi_{\mathscr{A}_{i}(\nu^{k})}(t)+\frac{1}{\varepsilon}y^{k+1}(t)\sum_{i=3}^{6}\chi_{\mathscr{A}_{i}(\nu^{k})}(t)$
$\displaystyle\mathcal{H}^{k}(t;p^{k+1}(t))$
$\displaystyle=\chi_{{\mathscr{I}^{\mathcal{U}}}(z^{k})}(t)\gamma_{c}\int_{\Gamma_{c}}p^{k+1}(t)\mathrm{d}s$
$\displaystyle r^{k}(t)$
$\displaystyle=\chi_{{\mathscr{A}_{{\mathsf{a}}}^{\mathcal{U}}}(z^{k})}(t)u_{\mathsf{a}}(t)+\chi_{{\mathscr{A}_{{\mathsf{b}}}^{\mathcal{U}}}(z^{k})}(t)u_{\mathsf{b}}(t)$
$\displaystyle r^{k,\varepsilon}(t)$
$\displaystyle=y_{\mathsf{a}}\chi_{\mathscr{A}_{1}^{\mathcal{W}}(\nu^{k})}(t)+y_{\mathsf{b}}\chi_{\mathscr{A}_{2}^{\mathcal{W}}(\nu^{k})}(t)-w_{\mathsf{a}}\big{(}\chi_{\mathscr{A}_{3}^{\mathcal{W}}(\nu^{k})}(t)+\chi_{\mathscr{A}_{5}^{\mathcal{W}}(\nu^{k})}(t)\big{)}$
$\displaystyle\quad-
w_{\mathsf{b}}\big{(}\chi_{\mathscr{A}_{4}^{\mathcal{W}}(\nu^{k})}(t)+\chi_{\mathscr{A}_{6}^{\mathcal{W}}(\nu^{k})}(t)\big{)}$
$\displaystyle\quad+\frac{\varepsilon+1}{\varepsilon}y_{\mathsf{a}}\big{(}\chi_{\mathscr{A}_{3}^{\mathcal{W}}(\nu^{k})}(t)+\chi_{\mathscr{A}_{4}^{\mathcal{W}}(\nu^{k})}(t)\big{)}$
$\displaystyle\quad+\frac{\varepsilon+1}{\varepsilon}y_{\mathsf{b}}\big{(}\chi_{\mathscr{A}_{5}^{\mathcal{W}}(\nu^{k})}(t)+\chi_{\mathscr{A}_{6}^{\mathcal{W}}(\nu^{k})}(t)\big{)},$

with $\chi_{\bullet}(t)$ being the indicator functions of the sets

$\displaystyle\mathscr{A}_{1}^{\mathcal{W}}(\nu^{k})$
$\displaystyle={\mathscr{A}_{\mathsf{a}}^{\mathcal{W}}}(\nu^{k})\cap{\mathscr{I}^{\mathcal{W}}_{2}}(\nu^{k}),$
$\displaystyle\mathscr{A}_{2}^{\mathcal{W}}(\nu^{k})$
$\displaystyle={\mathscr{A}_{\mathsf{b}}^{\mathcal{W}}}(\nu^{k})\cap{\mathscr{I}^{\mathcal{W}}_{2}}(\nu^{k}),$
$\displaystyle\mathscr{A}_{3}^{\mathcal{W}}(\nu^{k})$
$\displaystyle={\mathscr{A}_{\mathsf{a}}^{\mathcal{W}}}(\nu^{k})\cap{\mathscr{A}_{\mathsf{a},2}^{\mathcal{W}}}(\nu^{k}),$
$\displaystyle\mathscr{A}_{4}^{\mathcal{W}}(\nu^{k})$
$\displaystyle={\mathscr{A}_{\mathsf{a}}^{\mathcal{W}}}(\nu^{k})\cap{\mathscr{A}_{\mathsf{b},2}^{\mathcal{W}}}(\nu^{k}),$
$\displaystyle\mathscr{A}_{5}^{\mathcal{W}}(\nu^{k})$
$\displaystyle={\mathscr{A}_{\mathsf{b}}^{\mathcal{W}}}(\nu^{k})\cap{\mathscr{A}_{\mathsf{a},2}^{\mathcal{W}}}(\nu^{k}),$
$\displaystyle\mathscr{A}_{6}^{\mathcal{W}}(\nu^{k})$
$\displaystyle={\mathscr{A}_{\mathsf{b}}^{\mathcal{W}}}(\nu^{k})\cap{\mathscr{A}_{\mathsf{b},2}^{\mathcal{W}}}(\nu^{k}),$

respectively. Summarizing, the coupled system of equations (8) can be
formulated only in the variables $y^{k+1}$ and $p^{k+1}$:

$\left(\begin{array}[]{cc}\mathcal{A}_{11}^{k}&\mathcal{A}_{12}^{k}\\\\[4.30554pt]
\mathcal{A}_{21}^{k}&\mathcal{A}_{22}^{k}\end{array}\right)\left(\begin{array}[]{c}y^{k+1}\\\\[4.30554pt]
p^{k+1}\end{array}\right)=\left(\begin{array}[]{c}\mathcal{Q}_{1}(\nu^{k};u_{\mathsf{a}},u_{\mathsf{b}},\gamma_{c},\gamma,y_{\text{out}})\\\\[4.30554pt]
\mathcal{Q}_{2}(\nu^{k};y_{\mathsf{a}},y_{\mathsf{b}},w_{\mathsf{a}},w_{\mathsf{b}},\varepsilon,\sigma)\end{array}\right),$
(9)

where the operators $\mathcal{A}_{11}^{k}$ and $\mathcal{A}_{22}^{k}$ are of
the form

$\mathcal{A}_{11}^{k}=\mathcal{C}+\tilde{\mathcal{A}}_{11}^{k},\quad\mathcal{A}_{22}^{k}=\mathcal{C}^{\star}+\tilde{\mathcal{A}}_{22}^{k}$

and $\mathcal{C}$ stands for the heat-convection operator, which does not
depend on $k$. A discretization of (9) leads to a discretized system of the
form

$\left(\begin{array}[]{cc}\mathrm{A}_{11}^{k}&\mathrm{A}_{12}^{k}\\\\[4.30554pt]
\mathrm{A}_{21}^{k}&\mathrm{A}_{22}^{k}\end{array}\right)\left(\begin{array}[]{c}\mathrm{y}^{k+1}\\\\[4.30554pt]
\mathrm{p}^{k+1}\end{array}\right)=\left(\begin{array}[]{c}\mathrm{Q}_{1}^{k}\\\\[4.30554pt]
\mathrm{Q}_{2}^{k}\end{array}\right)$ (10)

with

$\mathrm{A}_{11}^{k}=\mathrm{C}+\tilde{\mathrm{A}}_{11}^{k},\quad\mathrm{A}_{22}^{k}=\mathrm{C}^{\top}+\tilde{\mathrm{A}}_{22}^{k},$

where $\mathrm{C}$ stands for the discretized heat-convection operator, which
is again independent of $k$. For further details, we refer to [27]. We
summarize the PDASS in Algorithm 1.

Algorithm 1 (Primal-dual active set strategy)

1: Choose starting value $\nu^{0}=(u^{0},w^{0},\vartheta^{0})\in\mathbb{N}$,
set $k=0$ and flag = false;

2: Determine $y^{0}=\mathcal{S}u^{0}$ and $p^{0}=-\mathcal{A}(\sigma
w^{0}+\vartheta^{0})/\varepsilon$.

3: Determine ${\mathscr{A}_{{\mathsf{a}}}^{\mathcal{U}}}(\nu^{0})$,
${\mathscr{A}_{{\mathsf{b}}}^{\mathcal{U}}}(\nu^{0})$,
${\mathscr{I}^{\mathcal{U}}}(\nu^{0})$ for $i=1,\ldots,m$,
${\mathscr{A}_{\mathsf{a}}^{\mathcal{W}}}(\nu^{0})$,
${\mathscr{A}_{\mathsf{b}}^{\mathcal{W}}}(\nu^{0})$,
${\mathscr{I}^{\mathcal{W}}}(\nu^{0})$ and
${\mathscr{A}_{\mathsf{a},2}^{\mathcal{W}}}(\nu^{0})$,
${\mathscr{A}_{\mathsf{a},2}^{\mathcal{W}}}(\nu^{0})$ and
${\mathscr{I}^{\mathcal{W}}_{2}}(\nu^{0})$ from (6);

4: repeat

5: Compute the solution $(y^{k+1},p^{k+1})$ by solving (10);

6: Compute $\nu^{k+1}=(u^{k+1},w^{k+1},\vartheta^{k+1})\in\mathbb{N}$ from (7)
and set $k=k+1$;

7: Determine ${\mathscr{A}_{{\mathsf{a}}}^{\mathcal{U}}}(\nu^{k})$,
${\mathscr{A}_{{\mathsf{b}}}^{\mathcal{U}}}(\nu^{k})$,
${\mathscr{I}^{\mathcal{U}}}(\nu^{k})$ for $i=1,\ldots,m$,
${\mathscr{A}_{\mathsf{a}}^{\mathcal{W}}}(\nu^{k})$,
${\mathscr{A}_{\mathsf{b}}^{\mathcal{W}}}(\nu^{k})$,
${\mathscr{I}^{\mathcal{W}}}(\nu^{k})$ and
${\mathscr{A}_{\mathsf{a},2}^{\mathcal{W}}}(\nu^{k})$,
${\mathscr{A}_{\mathsf{a},2}^{\mathcal{W}}}(\nu^{k})$ and
${\mathscr{I}^{\mathcal{W}}_{2}}(\nu^{k})$ from (6);

8: if
$\mathscr{A}_{{\mathsf{a}}i}^{\mathcal{U}}(\nu^{k})=\mathscr{A}_{{\mathsf{a}}i}^{\mathcal{U}}(\nu^{k-1})$
and
$\mathscr{A}_{{\mathsf{b}}i}^{\mathcal{U}}(\nu^{k})=\mathscr{A}_{{\mathsf{b}}i}^{\mathcal{U}}(\nu^{k-1})$
for $i=1,\dots,m$ then

9: if
${\mathscr{A}_{\mathsf{a}}^{\mathcal{W}}}(\nu^{k})={\mathscr{A}_{\mathsf{a}}^{\mathcal{W}}}(\nu^{k-1})$
and
${\mathscr{A}_{\mathsf{b}}^{\mathcal{W}}}(\nu^{k})={\mathscr{A}_{\mathsf{b}}^{\mathcal{W}}}(\nu^{k-1})$
then

10: if
${\mathscr{A}_{\mathsf{a},2}^{\mathcal{W}}}(\nu^{k})={\mathscr{A}_{\mathsf{a},2}^{\mathcal{W}}}(\nu^{k-1})$
and
${\mathscr{A}_{\mathsf{b},2}^{\mathcal{W}}}(\nu^{k})={\mathscr{A}_{\mathsf{b},2}^{\mathcal{W}}}(\nu^{k-1})$
then

11: Set flag = true;

12: end if

13: end if

14: end if

15: until flag = true;

######  Remark 3.

The discrete linear system (10) can be obtained discretizing the primal and
dual equations (8), for example, with piecewise linear finite elements (FE) in
space and the implicit Euler scheme in time. Let $V_{h}\subset V$ be the FE
space with $\dim V_{h}=N_{\bm{x}}$ and $N_{t}$ be the number of time steps.
Notice that the matrix in (10) has the size
$(2N_{t}N_{\bm{x}})\times(2N_{t}N_{\bm{x}})$. Clearly, its dimension may cause
problems regarding memory consumption and computational time. In our
application, in fact, $N_{\bm{x}}$ is not small, since a coarse FE grid can
not capture the complexity of the dynamics involved. Therefore, we apply POD-
based reduced-order modeling to approximate the solution of (10) by a reduced-
order system of dimension $2\ell N_{t}$ with $\ell\ll N_{\bm{x}}$, gaining
computational time but paying in approximation; cf. Section 4. The number of
time steps may not be small either, since the time horizon $T$ could be
extremely large (e.g. one month) and the time step small (e.g. a minute). This
is not a disadvantage for us, since because we apply an MPC method, we
iteratively solve optimal control subproblems on a smaller time horizon,
computing a feedback control that approximates the optimal control
$\bar{u}^{\varepsilon}$; cf. Section 5. Finally, note that the matrix in (10)
has a sparse block structure that allows to store only the FE matrices, e.g.,
the mass and stiffness matrices, which have dimension $N_{\bm{x}}\times
N_{\bm{x}}$. These blocks will then be projected into the POD space. Doing
that we can save memory significantly, in particular when the number of FE
nodes is large. Note that, in comparison to the Boussinesq approximation, we
gain even more, since in general computing the solution of this model requires
an additional picewise quadratic FE approximation of the velocity space, which
would definitely become infeasible in combination with the PDASS.

##  4 POD and a-posteriori error estimate

Model order reduction and, in particular, POD is a well established field of
numerical analysis, which was intensively developed in the past 20 years. The
key idea is to construct a low-dimensional subspace $V_{\ell}\subset
V_{h}=\mathrm{span}\,\\{\varphi_{1},\ldots,\varphi_{N_{\bm{x}}}\\}$, spanned
by the so-called POD basis functions. For the sake of completeness, in the
first part of this section, we describe briefly the POD method in the discrete
case, but we focus on one snapshot ensemble for brevity. For more details we
refer the reader to [20], for instance. Let $t_{j}=j\Delta t$, $j=0,\ldots,n$,
be the time discretization for a suitable time interval with constant time
step $\Delta t>0$. By $\\{y_{j}\\}_{j=0}^{n}\subset V_{h}$ we denote a given
discrete trajectory, also called snapshots ensemble. We define the snapshot
subspace $\mathscr{V}_{n}=\mathrm{span}\,\\{y_{0},\ldots,y_{n}\\}\subset
V_{h}$. Let $d_{n}=\dim\mathscr{V}_{n}\leq n+1$. For any
$\ell\in\\{1,\ldots,d_{n}\\}$ we are interested in finding an orthonormal set
$\\{\psi_{i}\\}_{i=1}^{\ell}\subset\mathscr{V}_{n}$, which is the solution of
the following minimization problem

$\min\sum_{j=0}^{n}\alpha_{j}\Big{\|}y_{j}-\sum_{i=1}^{\ell}{\langle
y_{j},\psi_{i}\rangle}_{V}\,\psi_{i}\Big{\|}_{V}^{2}\quad\text{s.t.}\quad{\langle\psi_{i},\psi_{j}\rangle}_{V}=\delta_{ij},$
(11)

where $\delta_{ij}$ is the Kronecker delta and $\alpha_{j}>0$ are given
(trapezoidal) weights. We set
$V^{\ell}=\mathrm{span}\,\\{\psi_{1},\ldots,\psi_{\ell}\\}\subset V_{h}$.
Obviously, if $\dim V_{\ell}=\ell\ll m$ holds we gain a computational speed-
up. From $y_{j}\in V_{h}$ ($0\leq j\leq n$) we infer that there exists a
snapshot matrix $Y\in\mathbb{R}^{N_{\bm{x}}\times(n+1)}$, such that

$\displaystyle
y_{j}({\bm{x}})=\sum_{i=1}^{N_{\bm{x}}}Y_{ij}\varphi_{i}({\bm{x}})\quad\text{for
}{\bm{x}}\in\Omega.$

Introducing the weighting matrices
$D=\mathrm{diag}\,(\alpha_{0},\ldots,\alpha_{n})\in\mathbb{R}^{(n+1)\times(n+1)}$
and $W=((W_{ij}))\in\mathbb{R}^{N_{\bm{x}}\times N_{\bm{x}}}$ with
$W_{ij}=\langle\varphi_{j},\varphi_{i}\rangle_{V}$, a solution
$\\{\psi_{j}\\}_{j=1}^{\ell}$ to (11) can be computed as follows:

  * •

Solve the symmetric eigenvalue problem

$\displaystyle\big{(}D^{1/2}Y^{\top}WYD^{1/2}\big{)}\bm{\phi}_{j}=\lambda_{j}\bm{\phi}_{j}\quad\text{for
}j=1,\ldots,\ell$

with $\lambda_{1}\geq\ldots\geq\lambda_{\ell}>0$.

  * •

Define the vectors
$\bm{\psi}_{j}=YD^{1/2}\bm{\phi}_{j}/\sqrt{\lambda_{j}}\in\mathbb{R}^{N_{\bm{x}}}$
for $1\leq j\leq\ell$ and the matrix
$\Psi=[\bm{\psi}_{1}|\ldots|\bm{\psi}_{\ell}]\in\mathbb{R}^{N_{\bm{x}}\times\ell}$.

  * •

Set

$\displaystyle\psi_{j}({\bm{x}})=\sum_{i=1}^{N_{\bm{x}}}\Psi_{ij}\varphi_{i}({\bm{x}})\quad\text{for
}j=1,\ldots,\ell\text{ and }{\bm{x}}\in\Omega.$

Using a Galerkin ansatz, the reduced-order version of (8) can be derived
replacing the test functions $\varphi$ by the POD basis
$\left\\{\psi_{i}\right\\}_{i=1}^{\ell}$. We indicate with
$\mathcal{S}^{\ell}$ the POD solution operator corresponding to the reduced-
order model. We can then construct the discretized system

$\left(\begin{array}[]{cc}\mathrm{A}_{11}^{k,\ell}&\mathrm{A}_{12}^{k,\ell}\\\\[4.30554pt]
\mathrm{A}_{21}^{k,\ell}&\mathrm{A}_{22}^{k,\ell}\end{array}\right)\left(\begin{array}[]{c}\mathrm{y}^{k+1,\ell}\\\\[4.30554pt]
\mathrm{p}^{k+1,\ell}\end{array}\right)=\left(\begin{array}[]{c}\mathrm{Q}_{1}^{k,\ell}\\\\[4.30554pt]
\mathrm{Q}_{2}^{k,\ell}\end{array}\right),$ (12)

of which the matrix belongs to $\mathbb{R}^{(2\ell N_{t})\times(2\ell
N_{t})}$. From (12), one can easily implement a POD-based PDASS following the
structure of Algorithm 1. Now, a first question is how to check that the
chosen number of POD basis $\ell$ is good enough to guarantee a small
approximation error in reconstructing the snapshots. For this, the following
error formula holds

$\sum_{j=0}^{n}\alpha_{j}\bigg{\|}y_{j}-\sum_{i=1}^{\ell}{\langle
y_{j},\psi_{i}\rangle}_{V}\psi_{i}\bigg{\|}_{V}^{2}=\sum_{i=\ell+1}^{d_{n}}\lambda_{i}.$
(13)

More importantly, to check that the constructed POD basis is able to
reconstruct the full-order optimal solution $\bar{u}^{\varepsilon}$, we need
an a-posteriori error estimator. Let $\bar{u}^{\varepsilon,\ell}$ be the
solution of the POD-based PDASS, then the following result from [27] holds:

######  Proposition 1.

Let
$(\bar{u}^{\varepsilon},\bar{w}^{\varepsilon})\in\mathcal{Z}^{\varepsilon}_{\mathsf{ad}}$
be the optimal solution to ($\mathbf{\hat{P}}^{\varepsilon}$). Suppose that
$z^{\mathsf{ap}}=(u^{\mathsf{ap}},w^{\mathsf{ap}})\in\mathcal{Z}^{\varepsilon}_{\mathsf{ad}}$
is given arbitrarily. Then, there exists a perturbation
$\zeta=(\zeta^{u},\zeta^{w})\in\mathcal{U}\times\mathcal{W}$, which is
independent of $(\bar{u}^{\varepsilon},\bar{w}^{\varepsilon})$, so that

$\|\bar{u}^{\varepsilon}-u^{\mathsf{ap}}\|_{\mathcal{U}}+\|\bar{w}^{\varepsilon}-w^{\mathsf{ap}}\|_{\mathcal{W}}\leq\frac{1}{\sigma^{\mathsf{ap}}}\,{\|\mathcal{T}^{\star}\zeta\|}_{\mathcal{U}},$
(14)

where $\sigma^{\mathsf{ap}}:=\min\\{\sigma,1\\}>0$ and $\mathcal{T}^{\star}$
is defined as

$\mathcal{T}^{\star}=\left(\begin{array}[]{cc}\mathcal{I}_{\mathcal{U}}&\mathcal{S}^{\star}\\\
0&\varepsilon\mathcal{I}_{\mathcal{W}}\end{array}\right):\mathcal{U}\times\mathcal{W}\to\mathcal{U}\times\mathcal{W},$

where $\mathcal{I}_{\mathcal{U}}$ and $\mathcal{I}_{\mathcal{W}}$ are the
identities on $\mathcal{U}$ and $\mathcal{W}$, respectively. The perturbation
$\zeta$ is computed as follows: Let
$\xi=(\xi^{u},\xi^{w})\in\mathcal{U}\times\mathcal{W}$ be given as the
solution of the linear system
$\mathcal{T}^{\star}\xi=\nabla\hat{J}(z^{\mathsf{ap}})$, i.e.,

$\left(\begin{array}[]{cc}\mathcal{I}_{\mathcal{U}}&\mathcal{S}^{\star}\\\
0&\varepsilon\mathcal{I}_{\mathcal{W}}\end{array}\right)\left(\begin{array}[]{c}\xi^{u}\\\
\xi^{w}\end{array}\right)=\left(\begin{array}[]{c}u^{\mathsf{ap}}\\\ \sigma
w^{\mathsf{ap}}\end{array}\right),$ (15)

then

$\zeta^{u}(t)=\left\\{\begin{array}[]{ll}-\min\\{0,\xi^{u}(t)\\}&\text{for
}t\in{\mathscr{A}_{{\mathsf{a}}}^{\mathcal{U}}}(z^{\mathsf{ap}}),\\\\[2.84526pt]
-\max\\{0,\xi^{u}(t)\\}&\text{for
}t\in{\mathscr{A}_{{\mathsf{b}}}^{\mathcal{U}}}(z^{\mathsf{ap}}),\\\\[2.84526pt]
-\xi^{u}(t)&\text{for
}t\in{\mathscr{I}^{\mathcal{U}}}(z^{\mathsf{ap}})\end{array}\right.$ (16a) and
$\zeta^{w}(t,{\bm{x}})=\left\\{\begin{array}[]{ll}-\min\\{0,\xi^{w}(t,{\bm{x}})\\}&\text{for
}(t,{\bm{x}})\in{\mathscr{A}_{\mathsf{a}}^{\mathcal{W}}}(z^{\mathsf{ap}}),\\\\[2.84526pt]
-\max\\{0,\xi^{w}(t,{\bm{x}})\\}&\text{for
}(t,{\bm{x}})\in{\mathscr{A}_{\mathsf{b}}^{\mathcal{W}}}(z^{\mathsf{ap}}),\\\\[2.84526pt]
-\xi^{w}(t,{\bm{x}})&\text{for
}(t,{\bm{x}})\in{\mathscr{I}^{\mathcal{W}}}(z^{\mathsf{ap}}).\end{array}\right.$
(16b)

######  Remark 4.

Note that we can easily decouple the equations in system (15) by computing
$\xi^{w}=\sigma_{w}w^{\mathsf{ap}}/\varepsilon$ from the second equation. Then
$\xi^{u}=-\mathcal{S}^{\star}\xi^{w}+u^{\mathsf{ap}}=-\sigma\mathcal{S}^{\star}w^{\mathsf{ap}}/\varepsilon+u^{\mathsf{ap}}$.
Moreover, we recall that
$\mathcal{S}^{\star}:W_{0}(0,T)^{\prime}\to\mathcal{U}$ denotes the dual
operator of the linear solution operator $\mathcal{S}$; see [41, Lemma 2.4].
Finally, we remark that in our numerical realization $z^{\mathsf{ap}}$ is
given by the POD suboptimal solution pair
$(\bar{u}^{\varepsilon,\ell},\bar{w}^{\varepsilon,\ell})$. Thus, (14) can be
utilized as an a-posteriori error estimate. For further details regarding the
a-posteriori error estimator, we refer to [27]. $\Diamond$

##  5 Economic MPC

MPC is a well-established method for computing a closed-loop control for a
dynamical system on an infinite (or large) time horizon. MPC splits the
solution of such problems into the consecutive solution of problems on a small
finite time horizon, the so called prediction horizon. The goal is to reduce
the complexity of the problem in time, while at the same time being able to
react to model uncertainties, disturbances, and possible parameter changes. We
are interested in economic MPC [4, 11]. This term comprises all MPC schemes in
which the cost function does not merely penalize the distance to a pre-defined
target. Here we do not have such a pre-defined target. Rather, the goal is to
minimize the control effort, while satisfying certain control and state
constraints; see, e.g., [17, 19, 31]. Applying MPC to the virtual control
problem ($\mathbf{\hat{P}}^{\varepsilon}$) has several advantages, which will
be more clear after having introduced the method. Let $N\Delta t$,
$N\in\mathbb{N}$, be the prediction horizon. We define the Hilbert spaces

$\mathcal{U}_{N}^{j}=L^{2}(t_{j},t_{j+N}),\quad\mathcal{W}_{N}^{j}=L^{2}((t_{j},t_{j+N})\times\Omega).$

for $j=0,\ldots,N_{t}$ and with the notation $t_{j}=j\Delta t$ for $j>N_{t}$.
Then the MPC cost functional is

$J^{j}_{N}(u,w)=\frac{1}{2}\|u\|^{2}_{\mathcal{U}_{N}^{j}}+\frac{\sigma}{2}\|w\|^{2}_{\mathcal{W}_{N}^{j}}$

and the MPC admissible sets are

$\displaystyle\mathcal{Z}^{\varepsilon,j}_{N,\mathsf{ad}}=\big{\\{}$
$\displaystyle(u,w)\in\mathcal{U}_{N}^{j}\times\mathcal{W}_{N}^{j}:u,w,\mathcal{S}u+\varepsilon
w\text{ satisfy }$ $\displaystyle u_{\mathsf{a}}(t)\leq u(t)\leq
u_{\mathsf{b}}(t)\text{ a.e. in }(t_{j},t_{j+N}],$ $\displaystyle
w_{\mathsf{a}}\leq w(t,{\bm{x}})\leq
w_{\mathsf{b}},\,y_{\mathsf{a}}\leq\mathcal{S}u+\varepsilon w\leq
y_{\mathsf{b}}\text{ a.e. in }(t_{j},t_{j+N}]\times\Omega,$
$\displaystyle\text{ respectively}\big{\\}},$

for $j=0,\ldots,N_{t}$. Furthermore, note that the fully discretized version
of (5a) can be rewritten as the following discrete dynamical system

$\displaystyle y_{j+1}=f(j,y_{j},u_{j})\text{ for }j=0,\ldots,N_{t},\quad
y_{0}=y_{\circ},$ (17)

where the function $f$ comprises all the formulas that come from the
discretization and the model, which we avoid to write for brevity. The MPC
method is then summarized in Algorithm 2.

Algorithm 2 Model predictive control

1: Choose an initial state $y_{0}$, an MPC prediction horizon $N$ and the
regularization parameter $\varepsilon>0$;

2: for each time instant $j=0,1,2,\ldots,N_{t}$ do

3: Measure the current state of the system $y^{j}$ and the model parameters;

4: Solve the optimal control problem  $\min
J_{N}^{j}(u,w)\quad\text{s.t.}\quad(u,w)\in\mathcal{Z}^{\varepsilon,j}_{N,\mathsf{ad}}$
(18) to obtain the open-loop optimal control sequence
$\bar{u}^{\varepsilon}_{N}$;

5: Apply the first element of $\bar{u}^{\varepsilon}_{N}$ as a control to the
system (17) during the next sampling period, i.e. use the feedback law
$\mu_{N}(y_{j})=\bar{u}^{\varepsilon}_{N}(0)$.

6: end for

To guarantee well-posedness and existence of the optimal control to each open-
loop optimal control problem in Algorithm 2, we assume that the admissible
sets $\mathcal{Z}^{\varepsilon,j}_{N,\mathsf{ad}}$ have non-empty interior for
all $j=0,\ldots,N_{t}$.

Algorithm 2 produces a sequence of optimal controls in discrete time and uses
the first element of each optimal control sequence as a feedback control
value. This results in a closed-loop sub-optimal trajectory, which solves (17)
for $u_{j}=\mu_{N}(y_{j})$, $j=0,\ldots,N_{t}$. Note that in step 3 of the MPC
algorithm it is possible to measure (and update) the parameters of the model.
This is particularly advantageous for our application, since usually we do not
have, e.g., a precise weather forecast for the entire long-time horizon $T$.
In this way, we can update the outside temperature while the algorithm runs.
To determine $y^{j}$ typically an estimator is needed to process available
sensor data properly to reconstruct the spatial-temporal distribution. In
addition, we can take advantage for the velocity field v. We mentioned already
that this will be provided from an approximation of the Boussinesq velocity
fields. We thus gain the possibility of updating this approximation while the
simulation runs, improving the accuracy of our simplified model.

The optimal control problems (18) are solved on a much smaller time horizon
$(t_{j},t_{j+N})$ since $N\ll N_{t}$. This guarantees a good computational
speed-up. For general optimal control problems it is not necessarily the case
that the MPC closed-loop trajectory is approximately optimal on the long time
horizon $N_{t}$. However, in recent years structural conditions were
discovered that allow to conclude this property for $N$ large enough; cf. [16,
Chapter 7] and [11, 19]. The most important among these conditions are the
turnpike property and the closely related strict dissipativity property [39,
14]. The turnpike property can be seen as a similarity condition for optimal
trajectories on different time horizons and we will check numerically that
this property holds for the optimal control problem considered in this paper;
cf. Section 7.1. It is also possible to check the quality of the MPC horizon
and trajectory with a-posteriori error estimators, such as in [13, 15], but we
do not utilize them in this work.

What remains to be explained is how to successfully combine MPC and POD. The
challenge is to guarantee a small approximation error for the reduced-order
model while the MPC algorithm proceeds. This is clearly a big challenge if one
thinks about the fact that the problem parameters are not only time-variant,
but that the parameter values at a certain time instant in the future may also
change when new measurement information becomes available. When this happens,
it means that we need to update the POD model. Luckily, the a-posteriori error
estimate from Proposition 1 is a useful tool to check whether this should be
done or not.

Algorithm 3 MPC-POD

1: Choose an initial state $y_{0}$, an MPC prediction horizon $N$, the
regularization parameter $\varepsilon>0$, and tolerances $\tau_{1}$ for the
a-priori estimate (13) and $\tau_{2}$ for the a-posteriori estimate;

2: Set $t_{0}=0$, $y_{0}(0)=y_{\circ}$, and flag = true;

3: for $j=0,1,2,\dots,N_{t}$ do

4: Measure the current state $y_{j}$ of the system at time $t_{j}$;

5: if flag = true then

6: Update the model parameters;

7: Compute primal and dual snapshots;

8: Compute a POD basis $\\{\psi_{i}\\}_{i=1}^{\ell}$ of rank $\ell$ according
to $\tau_{1}$;

9: Set flag = false;

10: end if

11: Solve the MPC-POD optimal control problem  $\min J^{j}_{N}(u,w)\text{ s.t.
}(u,w)\in\mathcal{Z}^{\varepsilon,\ell,j}_{N,\mathsf{ad}}$ (19) to obtain the
open-loop suboptimal control sequence $\bar{u}^{\varepsilon,\ell}_{N}$;

12: Apply the first element of $\bar{u}^{\varepsilon,\ell}_{N}$ as a control
to the system (17) during the next sampling period, i.e. use the feedback law
$\mu_{N}(y_{j})=\bar{u}^{\varepsilon,\ell}_{N}(0)$;

13: Compute the a-posteriori error estimate
$e={\|\mathcal{T}^{\star}\zeta\|}_{\mathcal{U}\times\mathcal{W}}$ from (14);

14: if $e>\tau_{2}$ then

15: Set flag = true;

16: end if

17: end for

We report the MPC-POD algorithm in Algorithm 3. For further details, we refer
to [27, 29]. The admissible sets in (19) are defined as

$\displaystyle\mathcal{Z}^{\varepsilon,\ell,j}_{N,\mathsf{ad}}=\big{\\{}$
$\displaystyle(u,w)\in\mathcal{U}_{N}^{j}\times\mathcal{W}_{N}^{j}:u,w,\mathcal{S}^{\ell}u+\varepsilon
w\text{ satisfy }$ $\displaystyle u_{\mathsf{a}}(t)\leq u(t)\leq
u_{\mathsf{b}}(t)\text{ a.e. in }(t_{j},t_{j+N}],$ $\displaystyle
w_{\mathsf{a}}\leq w(t,{\bm{x}})\leq
w_{\mathsf{b}},\,y_{\mathsf{a}}\leq\mathcal{S}^{\ell}u+\varepsilon w\leq
y_{\mathsf{b}}\text{ a.e. in }(t_{j},t_{j+N}]\times\Omega,$
$\displaystyle\text{ respectively}\big{\\}},$

where we recall that $\mathcal{S}^{\ell}$ is the POD solution operator for the
reduced-order state equation. Obviously, to guarantee existence and uniqueness
of the suboptimal controls solution to (19), we assume that the sets
$\mathcal{Z}^{\varepsilon,\ell,j}_{N,\mathsf{ad}}$ have non-empty interior for
all $j=0,\ldots,N_{t}$.

The only aspect to be explained is how to compute the snapshots in Algorithm
3. We will explain this in the next section. Let us mention that other
possibilities are explored in [10, 12, 27, 29].

##  6 Computational details and additional parallelization issues

In this section, we clarify some aspects of our simplification strategy to
compute an an approximate solution to ($\mathbf{\hat{P}}_{\mathsf{B}}$). We
give details about the implementation of the algorithm, mentioning which parts
can be performed in parallel.

Looking at Algorithm 3, there are still some unclear aspects. At first, it is
not clear how to compute the POD snapshots in Step 7. At iteration $j=0$, we
initialize our model by solving the Boussinesq approximation (1) with a given
initial control $u_{0}$ for the time horizon $[0,T+N\Delta t]$ using the FE
method for the spatial and the implicit Euler method for the temporal
discretization. The intervall $[0,T+N\Delta t]$ is discretized by
$t_{j}=j\Delta t$, $j=0,\ldots,m$ with $m=N_{t}+N$. Exploiting the power of
parallel computing we obtain, therefore, the solution sequence
$\\{(\textbf{v}_{j,h}^{0},y_{j,h}^{0})\\}_{j=0}^{m}$. Then, we compute the
sequence $\\{p_{j,h}^{0}\\}_{j=0}^{m}$ by solving in parallel

$\displaystyle\frac{1}{\Delta t}\,{\langle
p_{j,h}^{0}-p_{j+1,h}^{0},\varphi_{h}\rangle}_{H}+a(t_{j};\varphi_{h},p_{j,h}^{0})+\frac{\sigma}{\varepsilon^{2}}\,{\langle\mathcal{G}^{0,\varepsilon}(t_{j};y_{j,h}^{0}),\varphi_{h}\rangle}_{H}$
$\displaystyle\qquad=\frac{\sigma}{\varepsilon^{2}}\,{\langle
r^{0,\varepsilon}(t_{j}),\varphi_{h}\rangle}_{H}\quad\text{for all
}\varphi_{h}\in V_{h}\text{ and }j=m-1,\ldots,1,$
$\displaystyle\bigg{\langle}p_{m,h}^{0}+\frac{\sigma}{\varepsilon^{2}}\,\mathcal{G}^{0,\varepsilon}(t_{m};y_{m,h}^{0}),\varphi_{h}\bigg{\rangle}_{H}=\frac{\sigma}{\varepsilon^{2}}\,{\langle
r^{0,\varepsilon}(t_{m}),\varphi_{h}\rangle}_{H}\quad\text{for all
}\varphi_{h}\in V_{h}$

(compare (8b) and (8d)), where the terms $\mathcal{G}^{0,\varepsilon}$ and
$r^{0,\varepsilon}$ involve the computation of initial active and inactive
sets (6a)-(6b), which can be also performed in parallel. The snapshots
$\\{y_{j,h}^{0}\\}_{j=0}^{m}$ and $\\{p_{j,h}^{0}\\}_{j=0}^{m}$ are used to
compute the POD basis in step 8. Note that the discrete POD eigenvalue problem
has dimension $2m$, which may be large. We discuss later about this issue. For
subsequent basis updates at an MPC step $j$, triggered by the a-posteriori
error estimator, we will compute our snapshots with the same procedure, but
using the control $\bar{u}^{\varepsilon,\ell}_{N}$ computed at the previous
iteration $j-1$. Note that we get also an update for the velocity field v.

We use the velocity field v (computed at step 7) as fixed one for the open-
loop optimal control problem (19), solved with POD-PDASS. In such a way, every
time the basis and the velocity field are updated, the MPC can correct its
approximation. We remark that the POD-PDASS matrix for problem (19) has the
size $(2N\ell)\times(2N\ell)$ (and is thus small), therefore there is no need
of parallelization for solving the linear system (12). We anyway compute the
active and inactive sets (6a)-(6b) in parallel, since they depend on $N_{x}$
and not on $\ell$.

With the current information, one can already run Algorithm 3. In addition, we
use a couple of improvements, which we describe further. A first disadvantage
of the POD basis updates is that any time we need to compute the snapshots, we
need to solve the Boussinesq equation for $t\in[t_{j},T+N\Delta t]$. Although
the measure of this time domain decreases for increasing $j$, computing the
solution of the Boussinesq approximation (1) for this large time horizon can
be computationally demanding even using parallel computing. This computational
effort can be reduced, considering to generate snapshots in a smaller time
interval $[t_{j},t_{j}+(N+M)\Delta t]$, where $M>0$ is chosen to be small. In
such a way, we have an updated velocity field v (and snapshots) for the
subsequent $M$ iterations of the algorithm. If a POD basis update is performed
in one of these steps, we have saved computational time. If not, we need to
solve (1) for the remaining time horizon up to the end (or for other $M$
steps, for example). Having only $2(N+M)$ snapshots is also beneficial for
computing the POD basis, since we have to solve a smaller eigenvalue problem.

A second disadvantage of the POD basis update is that overwriting the past
information while computing a new basis may not be always a good choice. The
old snapshots, in fact, may still carry useful information for the future
computations. Yet, adding only snapshots over snapshots will increase the
number of POD basis needed to achieve the a-priori tolerance $\tau_{1}$. We
employ, therefore, the snapshots selection strategy proposed in [27]. We
report it in Algorithm 4. This strategy consists in keeping only those old
snapshots that represent dynamics close to the new ones. This procedure is
motivated by the fact that the MPC algorithm gives a forecast of the future
dynamics (new computed snapshots). Therefore, it makes sense to use them as
starting point and add eventually previously computed information, if this
describes similar trajectories.

Algorithm 4 Snapshots selection from [27].

0: Snapshots previously computed and stored in a list $\mathsf{L}$ and
tolerances $0<\rho<\varrho\ll 1$.

1: Compute the new $2(N+M)$ snapshots and store them in a list $\mathsf{S}$;

2: for $i\leq$ length($\mathsf{L}$) do

3: for $j\leq$ $2(N+M)$ do

4: if
$(1-\varrho)\|\mathsf{S}[j]\|\|\mathsf{L}[i]\|\leq|\langle\mathsf{S}[j],\mathsf{L}[i]\rangle|\leq(1-\rho)\|\mathsf{S}[j]\|\|\mathsf{L}[i]\|$
then

5: Add $\mathsf{L}[i]$ to $\mathsf{S}$;

6: break

7: end if

8: end for

9: end for

Furthermore, in step 12 of Algorithm 3, one can solve in parallel one time
step of the Boussinesq model (1) applying the feedback law, instead of solving
(5a). In such a way, we further improve the MPC approximation. This operation
is not costly if one considers that it is done only for one time-step. At
last, note that the computation of the a-posteriori error estimate (14)
requires full-order solves of state and adjoint equations. This can also be
carried out in parallel. To conclude, we resume every step of Algorithm 3,
which is performed in parallel:

  * •

Computation of POD snapshots: this involves solving the Boussinesq
approximation (1) and the adjoint equation (8b);

  * •

Computation of the active and inactive sets (6a)-(6b) of the POD-PDASS;

  * •

Computation of the initial guess for the next MPC iteration $j$: this involves
solving one time step of (1);

  * •

The POD a-posteriori error estimator (14).

##  7 Numerical tests

The numerical tests were performed on the compute server of the Chair of
Applied Mathematics at the University of Bayreuth, a computer with 32 physical
cores (64 logical) Intel(R) Xeon(R) CPU E7-2830 @ 2.13GHz and a total RAM of
512GB. For the numerical realization of the Boussinesq approximation, we refer
to [5].

###  7.1 Numerical verification of the turnpike property

In this subsection we check numerically that the turnpike property holds for
(18), because this is the most important ingredient for proving that MPC
provides nearly optimal solutions, see [19]. Before stating it in a
mathematical formalism, we need to introduce the concept of optimal operation
for an infinite horizon optimal control problem. Consider the discrete-time
time-varying control system (17). In this subsection, we denote with
$y_{u}(\cdot\,;t_{j},\tilde{y}_{\circ})$ the solution of (17) for a given
control $u\in\mathcal{U}_{N}^{j}$ starting in $\tilde{y}_{\circ}\in H$ at time
$t_{j}$. Considering the stage cost function $l:\mathbb{R}\times
H\to\mathbb{R}$, $l(u,w)=(|u|^{2}+\sigma\,\|w\|_{H}^{2})/2$, the cost
functional $J_{N}^{j}$ can be written as

$J_{N}^{j}(u,w)=\sum_{i=0}^{N}l(u(t_{j+i}),w(t_{j+i})).$

We recall that the finite horizon optimal control problem (18) is

$\min_{(u,w)\in\mathcal{Z}^{\varepsilon,j}_{N,\mathsf{ad}}}J_{N}^{j}(u,w)$

with (locally) optimal solution
$(\bar{u}^{\varepsilon}_{N},\bar{w}^{\varepsilon}_{N})\in\mathcal{Z}^{\varepsilon,j}_{N,\mathsf{ad}}$.
By $\mathcal{Z}^{\varepsilon,j}_{\infty,\mathsf{ad}}$, $J_{\infty}^{j}$,
$\bar{u}^{\varepsilon}_{\infty}$ and $\bar{w}^{\varepsilon}_{\infty}$ we
denote the extension of the above quantities to infinite time horizon problem.
Note that $\mathcal{Z}^{\varepsilon,j}_{N,\mathsf{ad}}$ changes if the
starting value $\tilde{y}_{\circ}\in H$ of the PDE is different. Therefore, in
what follows we indicate this dependency with
$\mathcal{Z}^{\varepsilon,j}_{N,\mathsf{ad}}(\tilde{y}_{\circ})$.

######  Definition 1.

Let $\tilde{y}_{\circ}\in H$ be an arbitrary feasible starting value at time
$t_{j}$. Consider a control sequence
$(\bar{u},\bar{w})\in\mathcal{Z}^{\varepsilon,j}_{\infty,\mathsf{ad}}(\tilde{y}_{\circ})$
with corresponding trajectory $\bar{y}$. We say that the system (17) is
optimally operated at $(\bar{y},\bar{u},\bar{w})$ if

$\liminf_{I\to\infty}\sum_{i=0}^{I-1}l(u(t_{j+i}),w(t_{j+i}))-l(\bar{u}(t_{j+i}),\bar{w}(t_{j+i}))\geq
0$

for all feasible $\hat{y}\in H$ and all
$(u,w)\in\mathcal{Z}^{\varepsilon,j}_{\infty,\mathsf{ad}}(\hat{y})$. The
trajectory $\bar{y}$ is called optimal operation trajectory.

We remark that there may exist more than one solution to
$\eqref{eq:discrete_dynamical_system}$ that satisfies Definition 1. Moreover,
it is well-known that, in general, it is not easy to compute an optimal
operation trajectory explicitly. Still, one can expect that the MPC trajectory
converges to this trajectory for sufficiently large MPC horizon $N$. To ensure
this behavior one has to assume structural conditions such as the turnpike
property. This property demands that the solution to the finite (and infinite)
time horizon optimal control problem is close for most of the time to an
optimal operation trajectory $\bar{y}$, according to the following definition.

Figure 2: Time-varying turnpike property for the finite horizon optimal
control problem.

######  Definition 2.

Consider $(\bar{y},\bar{u},\bar{w})$ at which system (17) is optimally
operated. The finite horizon optimal control problem (18) has the turnpike
property at $(\bar{y},\bar{u},\bar{w})$ if there exist $\eta\in\mathcal{L}$†††
$\mathcal{L}=\\{\eta\,:\,\mathbb{R}^{+}_{0}\to\mathbb{R}^{+}_{0}\big{|}\eta\text{
is continuous and strictly decreasing with }\lim_{s\to\infty}\eta(s)=0\\}.$
such that for each $j\in\mathbb{N}_{0}$, each feasible initial guess $y_{j}$,
each optimal trajectory $y_{\bar{u}^{\varepsilon}_{N}}(\cdot\,;t_{j},y_{j})$
and all $N,P\in\mathbb{N}$ there is a set
$\mathcal{Q}(t_{j},y_{j},P,N)\subseteq\\{0,\ldots,N\\}$ with at most $P$
elements and

$\displaystyle\|y_{\bar{u}^{\varepsilon}_{N}}(t_{j+M};t_{j},y^{j})-\bar{y}(t_{j+M})\|_{H}+|\bar{u}^{\varepsilon}_{N}(t_{j+M})-\bar{u}(t_{j+M})|$
$\displaystyle\quad+\|\bar{w}^{\varepsilon}_{N}(t_{j+M})-\bar{w}(t_{j+M})\|_{H}\leq\eta(P)$

for all $M\in\\{0,\ldots,N\\}\setminus\mathcal{Q}(t_{j},y_{j},P,N)$.

The definition of the turnpike property seems technically elaborated but it is
easy to illustrate graphically; cf. Figure 2. This property demands that for
each optimal open-loop trajectory there can only be a finite number $P$ of
time steps, $P$ independent from $N$, at which the trajectory is far from the
optimal operation trajectory. Moreover, it permits to define a bound on the
distance for the remaining close points. Contrary to the computation of the
optimal operation trajectory, it is easy to find numerical evidence that a
given optimal control problem has the turnpike property [18].

Figure 3:  Advection field v for the turnpike test at time instances
$t=1.0,2.0,3.0$.

For this test, we choose $\Omega=(0,1)^{2}$ discretized with a structured
quadrilateral mesh of $N_{\bm{x}}=1769$ nodes,
$\Gamma_{c}=(0,1)\times\\{0\\}$, $\rho=C_{p}=\alpha=1.0$, $\nu=0.1$,
$\kappa=0.025$ and $\textbf{g}=(0,-9.8)$. Moreover, we consider $\gamma=0.1$
and $\gamma_{c}=10^{5}$ for the boundary conditions, so that the action of the
control on the boundary can be seen (numerically) as the one of a Dirichlet
boundary condition. We solve at first (1) with fixed control $u(t)=22.5$,
outside temperature $y_{\text{out}}=\max(13,16-t)$ and initial guess
$(\textbf{v}_{0},y_{0})(t,{\bm{x}})=(0,20)$ for all ${\bm{x}}\in\Omega$ and
$t\in[0,T]$, in order to generate the advection field v to use in (18); cf.
Figure 3. Then, we solve (18) with $u_{\mathsf{a}}(t)=0$,
$u_{\mathsf{b}}(t)=10^{6}$, $w_{\mathsf{b}}=-w_{\mathsf{a}}=10^{9}$,
$\varepsilon=0.025$, $\sigma=1.0$, $y_{\mathsf{a}}(t)=17.5+\min(t,2)$ and
$y_{\mathsf{b}}(t)=23$ for all $t\in[0,T]$. We use a time step $\Delta t=0.01$
and compute the open-loop problem solution for $j=0$ (i.e. $t_{0}=0$) for
different horizons $N$, thus $T=N\Delta t$ in each test. As shown in Figure
4(a), the $L^{2}$ norm of the solution to (18) is the same for a larger number
of time steps as $N$ increases. Moreover, the solutions show a similar
behavior.

(a) $\|y(t)\|_{H}$, with $y$ solution to (18) for different horizons $N$. (b)
$\|y(t)\|_{H}$, with $y$ solution to (18) for different initial guesses and
fixed horizon $N=360$.

Figure 4: Turnpike test results.

This provides clear evidence for the turnpike behavior. This evidence is
further strengthened by the results from Figure 4(b), where the open-loop
trajectories starting from different initial guess are converging to the same
region as the time passes. For more information about how to check the
conditions for near optimal performance of MPC numerically, we refer to [18].

###  7.2 MPC-POD algorithm for the optimal control problem (2)

(a) Minimum (blue), maximum (red) and average (purple) temperature in the room
and state constraints $y_{\mathsf{a}}$ and $y_{\mathsf{b}}$ (black).  (b)
Number of active points with respect to the state constraints.

Figure 5: Results for Algorithm 3 with $\gamma=0.1$ and $\varepsilon=0.025$.

In this section we test Algorithm 3 with the same data of Section 7.1,
choosing an MPC horizon $N=300$. This horizon is a good comprise between
computational time and turnpike property; cf. Figure 4(a). We fix also $M=60$,
the a-priori tolerance $\tau_{1}=10^{-8}\sum_{i=1}^{\ell}\lambda_{i}$ and the
a posteriori error estimator tolerance $\tau_{2}=3.5$. Although $\tau_{2}$
seems large, one has to consider that the control will have a value close to
the state constraints range, roughly between 20 and 30, due to the fact that
the control boundary condition is (numerically) close to a Dirichlet one.
Moreover, the MPC horizon $N=300$ corresponds to a time horizon for the open-
loop problem equal to three. Therefore, one can expect an average $L^{2}$ norm
for the control of $43$ ($=25\sqrt{3}$). Hence, fixing $\tau_{2}=3.5$ equals
asking a relative error approximately of $8\%$. We expect, moreover, an
overestimation factor for the a-posteriori error estimator, since
$(\sigma^{\mathsf{ap}})^{-1}\|\mathcal{T}^{\star}\zeta\|_{\mathcal{U}}$ in
(14) gives an upper-bound also for the error in reconstructing the artificial
control, which might attain significantly greater values according to the
chosen $\varepsilon$. In addition, a too small $\tau_{2}$ will trigger the
basis updates too often and thus worsen the performance of the algorithm.

Figure 5(a) shows the minimum, maximum and average value of the optimal MPC
solution together with the desired bounds. One can notice that the average
temperature is kept inside the bounds, where instead the minimum and the
maximum one exceed them. This is typical of the $L^{2}$ regularization effect
of the virtual control approach and can be mitigated taking a small
$\varepsilon$. This effect is accentuated also by the POD method, which
provides a sub-optimal solution as well-known from the literature [32, 41]. As
a result, the average number of active points at each time step amounts to
$3.5\%$ of the total number of grid points. This can be seen in Figure 5(b).
Note that the number of points below the lower bound $y_{\mathsf{a}}$
increases as the time passes. The reason is that the upper bound
$y_{\mathsf{b}}$ prevents the control at the boundary from growing
significantly to counteract the decreasing outside temperature. This together
with the sub-optimality of the POD solution will prevent the control from
reacting faster and stronger to this external inputs. A possible solution is
to consider local constraints in the domain $\Omega$ [31] and/or local POD
basis functions [1, 3]. Nevertheless, the fact that the average temperature is
inside the bounds and only a small portion of discretization points are active
make Algorithm 3 a viable approach in HVAC of residual buildings.

(a) Minimum (blue), maximum (red) and average (purple) temperature in the room
and state constraints $y_{\mathsf{a}}$ and $y_{\mathsf{b}}$ (black).  (b)
Number of active points with respect to the state constraints.

Figure 6: Results for Algorithm 3 with $\gamma=0.01$ and $\varepsilon=0.0075$.

For smaller $\gamma=0.01$ and $\varepsilon=0.0075$, we can significantly
improve the results in terms of violation of the state constraints; cf. Figure
6. In this test, in fact, the influence of the outside temperature is lower
and the method is capable to reduce significantly the number of active points;
see Figure 6(b). In particular, the minimum temperature in the room never lies
below the given lower bound $y_{\mathsf{a}}$ (Figure 6(a)). The maximum
temperature, instead, is slightly above the upper bound $y_{\mathsf{b}}$. This
can be justified from the chosen value of the relaxation parameter
$\varepsilon$ and also from the fact that we are using the POD method, thus
the reduced-order model approximation introduce a further discretization
error, which may not be canceled even when passing to the limit
$\varepsilon\to 0$. This violation is anyway negligible for HVAC of residual
buildings, since the maximum value of the temperature attained all over the
time-space domain is $23.07$. Furthermore, we already discuss the conflict
(and a possible solution) between the state constraints and the control. This
can be clearly seen here, since the 33 active points are exactly the grid
points belonging to $\Gamma_{c}$.

##  8 Conclusion

We presented a POD-based economic MPC algorithm to handle HVAC of residual
buildings. The strategy consists in computing a sub-optimal solution to an
optimal control problem subjected to the Boussinesq approximation of Navier-
Stokes and bilateral state and control constraints. The goal was to have a
good compromise between robustness and computational cost. The first is
achieved via the MPC and a goal-oriented a-posteriori error estimate, while
the second is obtained through POD and parallel computing. The numerical
results confirm the expectations. While some limitations of the strategy show
up, particularly the emergence of several discretization points over time in
which the state constraints are violated, the resulting control strategy can
be regarded accurate enough for typical HVAC applications. The peformance can
be further improved, for example, by considering localized MOR and localized
constraints. This will be the focus of future work.

## Acknowledgments

The authors gratefully acknowledge support of DFG grants 274852524, 274852737
and 274853298. This work was performed under the auspices of the U.S.
Department of Energy by Lawrence Livermore National Laboratory under Contract
DE-AC52-07NA27344 (LLNL-JRNL-818244). This document was prepared as an account
of work sponsored by an agency of the United States government. Neither the
United States government nor Lawrence Livermore National Security, LLC, nor
any of their employees makes any warranty, expressed or implied, or assumes
any legal liability or responsibility for the accuracy, completeness, or
usefulness of any information, apparatus, product, or process disclosed, or
represents that its use would not infringe privately owned rights. Reference
herein to any specific commercial product, process, or service by trade name,
trademark, manufacturer, or otherwise does not necessarily constitute or imply
its endorsement, recommendation, or favoring by the United States government
or Lawrence Livermore National Security, LLC. The views and opinions of
authors expressed herein do not necessarily state or reflect those of the
United States government or Lawrence Livermore National Security, LLC, and
shall not be used for advertising or product endorsement purposes.

## References

  * [1] A. Abdulle and P. Henning.  A reduced basis localized orthogonal decomposition.  Journal of Computational Physics, 295:379–401, 2015.
  * [2] F. Abergel and R. Temam.  On Some Control Problem in Fluid Mechanics. Theoretical and Computational Fluid Dynamics.  Springer, 1990.
  * [3] F. Albrecht, B. Haasdonk, S. Kaulmann, and M. Ohlberger.  The localized reduced basis multiscale method.  In Proceedings of Algoritmy 2012, Conference on Scientific Computing, Vysoke Tatry, Podbanske, September 9-14, 2012, pages 393–403. Slovak University of Technology in Bratislava, Publishing House of STU, 2012.
  * [4] R. Amrit, J. B. Rawlings, and D. Angeli.  Economic optimization using model predictive control with a terminal cost.  Annual Rev. Control, 35:178–186, 2011.
  * [5] J. Andrej.  Modeling and optimal control of multiphysics problems using the finite element method. PhD thesis, Kiel University, 2019.
  * [6] F. Ballarin, A. Manzoni, A. Quarteroni, and G. Rozza.  Supremizer stabilization of POD–Galerkin approximation of parametrized steady incompressible Navier–Stokes equations.  International Journal for Numerical Methods in Engineering, 102(5):1136–1161, 2015.
  * [7] E. Casas.  Pontryagin’s principle for state-constrained boundary control problems of semilinear parabolic equations.  SIAM Journal on Control and Optimization, 35(4):1297–1327, 1997\.
  * [8] R. Dautray and P.-L. Lions.  Mathematical Analysis and Numerical Methods for Science and Technology. Volume 5: Evolution Problems I.  Springer, 2000.
  * [9] M. Ellis, J. Liu, and P. D. Christofides.  Economic Model Predictive Control: Theory, Formulations and Chemical Process Applications.  Springer Berlin, 2017.
  * [10] G. Fabrini, M. Falcone, and S. Volkwein.  Coupling MPC and HJB for the computation of POD-based feedback laws.  In F. A. Radu, K. Kumar, I. Berre, J. M. Nordbotten, and I. S. Pop, editors, Numerical Mathematics and Advanced Applications ENUMATH 2017, pages 941–949, Cham, 2019. Springer International Publishing.
  * [11] T. Faulwasser, L. Grüne, and M. A. Müller.  Economic nonlinear model predictive control.  Foundations and Trends® in Systems and Control, 5(1):1–98, 2018.
  * [12] J. Ghiglieri and S. Ulbrich.  Optimal flow control based on POD and MPC and an application to the cancellation of Tollmien–Schlichting waves.  Optimization Methods and Software, 29(5):1042–1074, 2014.
  * [13] L. Grüne, L. Mechelli, S. Pirkelmann, and S. Volkwein.  Performance estimates for economic model predictive control and their application in proper orthogonal decomposition- based implementations.  To appear in Mathematical Control and Related Fields., Preprint available at http://nbn-resolving.de/urn:nbn:de:bsz:352-2-ywxqb2eru5d19, 2020.
  * [14] L. Grüne and M. A. Müller.  On the relation between strict dissipativity and the turnpike property.  Syst. Contr. Lett., 90:45–53, 2016.
  * [15] L. Grüne and J. Pannek.  Practical NMPC suboptimality estimates along trajectories.  Systems & Control Letters, 58:161–168, 2009.
  * [16] L. Grüne and J. Pannek.  Nonlinear Model Predictive Control. Theory and Algorithms. 2nd Edition. Communications and Control Engineering. Springer International Publishing, 2017.
  * [17] L. Grüne and S. Pirkelmann.  Closed-loop performance analysis for economic model predictive control of time-varying systems.  In 2017 IEEE 56th Annual Conference on Decision and Control (CDC), page 5563–5569, 2017.
  * [18] L. Grüne and S. Pirkelmann.  Numerical verification of turnpike and continuity properties for time-varying pdes.  IFAC-PapersOnLine, 52(2):7 – 12, 2019.  3rd IFAC Workshop on Control of Systems Governed by Partial Differential Equations CPDE 2019.
  * [19] L. Grüne and S. Pirkelmann.  Economic model predictive control for time-varying system: Performance and stability results.  Optimal Control Applications and Methods, 41(1):42–64, 2020.
  * [20] M. Gubisch and S. Volkwein.  Proper orthogonal decomposition for linear-quadratic optimal control.  In P. Benner, A. Cohen, M. Ohlberger, and K. Willcox, editors, Model Reduction and Approximation: Theory and Algorithms, pages 5–66. SIAM, Philadelphia, PA, 2017.
  * [21] M. Hintermüller, K. Ito, and K. Kunisch.  The primal-dual active set strategy as a semismooth Newton method.  SIAM Journal on Optimization, 13(3):865–888, 2002.
  * [22] M. Hinze and U. Matthes.  Optimal and model predictive control of the Boussinesq approximation.  In K. Kunisch, J. Sprekels, G. Leugering, and F. Tröltzsch, editors, Control of Coupled Partial Differential Equations, pages 149–174, Basel, 2007. Birkhäuser Basel.
  * [23] M. Hinze, R. Pinnau, M. Ulbrich, and S. Ulbrich.  Optimization with PDE Constraints.  Springer Netherlands, 2009.
  * [24] P. Holmes, J. Lumley, G. Berkooz, and C. Rowley.  Turbulence, Coherent Structures, Dynamical Systems and Symmetry. 2nd Edition. Cambridge Monographs on Mechanics. Cambridge University Press, 2012.
  * [25] K. Ito and K. Kunisch.  Semi-smooth Newton methods for variational inequalities of the first kind.  ESAIM: Mathematical Modelling and Numerical Analysis, 37(1):41–62, 2003.
  * [26] K. Krumbiegel and A. Rösch.  A virtual control concept for state constrained optimal control problems.  Computational Optimization and Applications, 43:213–233, 2009\.
  * [27] L. Mechelli.  POD-based State-Constrained Economic Model Predictive Control of Convection-Diffusion Phenomena.  PhD thesis, University of Konstanz, 2019.
  * [28] L. Mechelli and S. Volkwein.  POD-based economic optimal control of heat-convection phenomena.  In M. Falcone, R. Ferretti, L. Grüne, and W. M. McEneaney, editors, Numerical Methods for Optimal Control Problems, pages 63–87, Cham, 2018. Springer International Publishing.
  * [29] L. Mechelli and S. Volkwein.  POD-based economic model predictive control for heat-convection phenomena.  In F. A. Radu, K. Kumar, I. Berre, J. M. Nordbotten, and I. S. Pop, editors, Numerical Mathematics and Advanced Applications ENUMATH 2017, pages 663–671, Cham, 2019. Springer International Publishing.
  * [30] F. Oldewurtel, A. Parisio, C. N. Jones, D. Gyalistras, M. Gwerder, V. Stauch, B. Lehmann, and M. Morari.  Use of model predictive control and weather forecasts for energy efficient building climate control.  Energy and Buildings, 45:15–27, 2012.
  * [31] S. Pirkelmann.  Economic Model Predictive Control and Time-Varying Systems.  PhD thesis, University of Bayreuth, 2020.
  * [32] S. S. Ravindran.  A reduced-order approach for optimal control of fluids using proper orthogonal decomposition.  International Journal for Numerical Methods in Fluids, 34(5):425–448, 2000.
  * [33] S. S. Ravindran.  Error analysis for Galerkin POD approximation of the nonstationary boussinesq equations.  Numerical Methods for Partial Differential Equations, 27(6):1639–1665, 2011.
  * [34] J. Rawlings, D. Mayne, and M. Diehl.  Model Predictive Control: Theory and Design. Second edition.  Nob Hill Publishing, 2017.
  * [35] J. P. Raymond.  Nonlinear boundary control of semilinear parabolic problems with pointwise state constraints.  Discrete & Continuous Dynamical Systems - A, 3:341–370, 1997.
  * [36] J. C. G. Reyes, P. Merino, J. Rehberg, and F. Tröltzsch.  Optimality conditions for state-constrained PDE control problems with time-dependent controls.  Control and Cybernetics, 37:5–38, 2008.
  * [37] G. Stabile and G. Rozza.  Finite volume POD-Galerkin stabilised reduced order methods for the parametrised incompressible Navier–Stokes equations.  Computers & Fluids, 173:273 – 284, 2018.
  * [38] R. Temam.  Navier-Stokes Equations, Theory and Numerical Analysis.  North-Holland, 1984.
  * [39] E. Trélat and E. Zuazua.  The turnpike property in finite-dimensional nonlinear optimal control.  J. Differential Equations, 258(1):81–114, 2015.
  * [40] D. J. Tritton.  Physical Fluid Dynamics.  Springer Netherlands, 1977.
  * [41] F. Tröltzsch and S. Volkwein.  POD a-posteriori error estimates for linear-quadratic optimal control problems.  Computational Optimization and Applications, 44:83–115, 2009.
  * [42] F. Trötzsch.  Regular Lagrange multipliers for control problems with mixed pointwise control- state constraints.  SIAM Journal on Optimization, 22:616–634, 2013.
  * [43] J. Široký, F. Oldewurtel, J. Cigler, and S. Prívara.  Experimental analysis of model predictive control for an energy efficient building heating system.  Applied Energy, 88:3079–3087, 2011.