Reduced Order Optimal Control of the Convective FitzHugh-Nagumo Equation

TL;DR

This paper compares three model order reduction methods—POD, DEIM, and DMD—for optimal control of the convective FitzHugh-Nagumo equations, focusing on accuracy and computational efficiency in blood coagulation modeling.

Contribution

It introduces a comparative analysis of POD, DEIM, and DMD methods for reduced order optimal control of convective FHN equations, highlighting their accuracy and speed.

Findings

POD provides the highest accuracy.

POD-DMD offers the fastest computation.

All methods are compared against full order solutions.

Abstract

In this paper, we compare three model order reduction methods: the proper orthogonal decomposition (POD), discrete empirical interpolation method (DEIM) and dynamic mode decomposition (DMD) for the optimal control of the convective FitzHugh-Nagumo (FHN) equations. The convective FHN equations consists of the semi-linear activator and the linear inhibitor equations, modeling blood coagulation in moving excitable media. The semilinear activator equation leads to a non-convex optimal control problem (OCP). The most commonly used method in reduced optimal control is POD. We use DEIM and DMD to approximate efficiently the nonlinear terms in reduced order models. We compare the accuracy and computational times of three reduced-order optimal control solutions with the full order discontinuous Galerkin finite element solution of the convection dominated FHN equations with terminal controls. Numerical results show that POD is the most accurate whereas POD-DMD is the fastest.

Tables

Table 1. Table 1: Relative L 2 superscript 𝐿 2 L^{2} -errors between full and reduced order states and control at the final time.

	$\Vert u_h-u_h^r\Vert$	$\Vert v_h-v_h^r\Vert$	$\Vert f_h-f_h^r\Vert$
POD	4.644e-03	6.167e-02	6.614e-01
POD–DEIM	4.988e-03	6.511e-02	7.123e-01
POD–DMD	5.787e-03	7.670e-02	1.061e+00

Table 2. Table 2: The values of the optimal objective functionals J 𝐽 J , nonlinear CG iterations and line searches, average number of Newton iterations per time step, CPU times and speedup factors.

	$J$	$\#$ CG	$\#$ Line	$\#$Newton	CPU	Speedup
		iterations	searches	iterations	times	factors
FOM	2.376e-03	184	710	3.22	15496.7	-
POD	2.396e-03	85	317	3.19	828.3	18.7
POD–DEIM	2.392e-03	88	332	3.20	386.2	40.1
POD–DMD	2.396e-03	99	370	1	8.0	1937.1

Equations

$$\underset{f\in {\mathcal{F}}_{ad}}{\min }J(u,v,f)=$$

$$+\frac{\upsilon }{2}\underset{0}{\overset{T}{\int }}\underset{\Omega }{\int }f(x,t{)}^{2}dxdt,$$

$$u_t(x,t)-d_u\Delta u(x,t)+\mathbf{b}(x)\cdot \nabla u(x,t)+g(u(x,t))+v(x,t)$$

$$v_t(x,t)-d_v\Delta v(x,t)+\mathbf{b}(x)\cdot \nabla v(x,t)+ϵ(v(x,t)-c_{3}u(x,t))$$

$${\partial }_{\text{n}}u(x,t)=0,{\partial }_{\text{n}}v(x,t)=0$$

$$u(x,t)=u_D(x,t),v(x,t)=v_D(x,t)$$

$$f\in {\mathcal{F}}_{ad}:=\{f\in L^{\infty }(Q):f_l\le f(x,t)\le f_{r}for\; a.e(x,t)\in Q\},$$

$$g(u)=c_{1}u(u-c_2)(u-1),$$

$$b_1(x)=a{x}_2(H-x_2),b_{max}=\frac{1}{4}a{H}^2,a>0,b_2(x)=0.$$

$$u_{0}(x)=\left\{\begin{array}[]{ll}\bar{u},&\hbox{if}\;\;x_{a}\leq x_{1}\leq x_{b},\\ 0,&\hbox{otherwise},\end{array}\right.\quad v_{0}(x)=0.$$

$$-p_t(x,t)-d_u\Delta p(x,t)-\mathbf{b}(x)\cdot \nabla p(x,t)+g^{\prime }(u)p(x,t)-ϵc_{3}q(x,t)$$

$$-q_t(x,t)-d_v\Delta q(x,t)-\mathbf{b}(x)\cdot \nabla q(x,t)+ϵq(x,t)+p(x,t)$$

$$d_u{\partial }_{\text{n}}p(x,t)+(\mathbf{b}\cdot \mathbf{n})p(x,t)=0,$$

$$p(x,t)=0,$$

$$p(x,T)=u(x,T)-u_T(x),$$

$$\displaystyle\int\limits_{0}^{T}\int\limits_{\Omega}\big{(}p(x,,t)+\upsilon f(x,t)\big{)}\big{(}\nu(x,t)-f(x,t)\big{)}\;dx\;dt\geq 0,\qquad\forall\nu\in\mathcal{F}_{ad},$$

$$f(x,t)={\mathbb{P}}_{[f_l,f_r]}\left(-\frac{1}{\upsilon }p(x,t)\right),$$

$$W_h=\left\{w\in L^2(\Omega ):w{\mid }_K\in {\mathbb{P}}^l(K),\forall K\in {\mathcal{T}}_h\right\}.$$

$${\Gamma }^{-}=\left\{x\in \partial \Omega :\mathbf{b}(x)\cdot \mathbf{n}(x)<0\right\},{\Gamma }^{+}=\partial \Omega \setminus {\Gamma }^{-}.$$

$$\partial K^{-}=\left\{x\in \partial K:\mathbf{b}(x)\cdot {\mathbf{n}}_K(x)<0\right\},\partial K^{+}=\partial K\setminus \partial K^{-},$$

$$\mbox{$\left[\!\left[w\right]\!\right]$}=w|_{E}\bm{n}_{K}+w^{e}|_{E}\bm{n}_{K^{e}},\quad\mbox{$\left\{\!\!\left\{w\right\}\!\!\right\}$}=\frac{1}{2}\big{(}w|_{E}+w^{e}|_{E}\big{)}.$$

$$\mbox{$\left[\!\left[\nabla w\right]\!\right]$}=\nabla w|_{E}\cdot\bm{n}_{K}+\nabla w^{e}|_{E}\cdot\bm{n}_{K^{e}},\quad\mbox{$\left\{\!\!\left\{\nabla w\right\}\!\!\right\}$}=\frac{1}{2}\big{(}\nabla w|_{E}+\nabla w^{e}|_{E}\big{)}.$$

$$⟨{\partial }_t{u}_h,w⟩+a_{h,u}(u_h,w)+g_h(u_h,w)+⟨v_h,w⟩$$

$$⟨u_h(x,0),w(x)⟩$$

$$⟨{\partial }_t{v}_h,w⟩+a_{h,v}(v_h,w)+ϵ⟨v_h,w⟩-ϵc_3⟨u_h,w⟩$$

$$⟨v_h(x,0),w(x)⟩$$

$$-⟨{\partial }_t{p}_h,w⟩+a_{h,p}(p_h,w)+⟨g^{\prime }(u_h)p_h,w⟩-ϵc_3⟨q_h,w⟩$$

$$⟨p_h(x,T),w(x)⟩$$

$$-⟨{\partial }_t{q}_h,w⟩+a_{h,q}(q_h,w)+ϵ⟨q_h,w⟩+⟨p_h,w⟩$$

$$⟨q_h(x,T),w(x)⟩$$

$$a_{h,i}(\tilde{w},w)=$$

$$\displaystyle+\sum\limits_{E\in\mathcal{E}^{0}_{h}\cup\mathcal{E}^{D}_{h}}\frac{d_{i}\gamma}{h_{E}}\int\limits_{E}\mbox{$\left[\!\left[\tilde{w}\right]\!\right]$}\cdot\mbox{$\left[\!\left[w\right]\!\right]$}\;ds+\sum\limits_{K\in\mathcal{T}_{h}}\int\limits_{K}\bm{b}\cdot\nabla\tilde{w}w\;dx$$

$$+\sum _{K\in {\mathcal{T}}_h}\underset{\partial K^{-}\backslash \partial \Omega }{\int }\mathbf{b}\cdot \mathbf{n}({\tilde{w}}^e-\tilde{w})wds-\sum _{K\in {\mathcal{T}}_h}\underset{\partial K^{-}\cap {\Gamma }^{-}}{\int }\mathbf{b}\cdot \mathbf{n}\tilde{w}wds,$$

$${\ell }_{h,i}(t,w)$$

$$g_h(u,w)$$

$$\sum\limits_{E\in\mathcal{E}^{N}_{h}}\int\limits_{E}\big{(}\bm{b}\cdot\bm{n}\big{)}p_{h}w\;ds\quad\hbox{and}\quad\sum\limits_{E\in\mathcal{E}^{N}_{h}}\int\limits_{E}\big{(}\bm{b}\cdot\bm{n}\big{)}q_{h}w\;ds,$$

$$u_h(x,t)=\sum _{i=1}^{N_K}\sum _{j=1}^{N_l}{u}_j^i(t){\varphi }_j^i(x),v_h(x,t)=\sum _{i=1}^{N_K}\sum _{j=1}^{N_l}{v}_j^i(t){\varphi }_j^i(x),f_h(x,t)=\sum _{i=1}^{N_K}\sum _{j=1}^{N_l}{f}_j^i(t){\varphi }_j^i(x).$$

$$u_{h,T}(x)=\sum _{i=1}^{N_K}\sum _{j=1}^{N_l}{u}_{T,j}^i{\varphi }_j^i(x),v_{h,T}(x)=\sum _{i=1}^{N_K}\sum _{j=1}^{N_l}{v}_{T,j}^i{\varphi }_j^i(x).$$

Full text

AMS subject classifications. 35K58, 65M60, 49J20.

††journal: Computers and Mathematics with Applications

1 Introduction

Optimal control problems (OCPs) governed by semilinear partial differential equations (PDEs) with wave-type solutions have been recently investigated for the Schlögl equation with traveling waves fronts and for the Nagumo equation with spiral waves [1], for the FitzHugh-Nagumo (FHN) equations with spiral and traveling waves [2, 3, 4], lambda-omega systems with spiral waves [5], Allee equation [6], optimal harvesting [7] and the convective FHN equations [8] with traveling wave solutions. PDE-constrained optimization is also applied to parameter identification in pattern formation. Optimal parameters leading to Turing patterns in semi-linear reaction-diffusion equations like the Schnakenberg and Gierer-Meinhardt equations are identified in [9, 10, 11]. Controlling traveling wave solutions or patterns of semi-linear PDEs is computationally challenging. For stable solutions of these OCPs, at each iteration step of the nonlinear optimization algorithms, large scale linear systems should be solved. Efficient and stable solutions of this kind of PDE constrained optimization problems are the focus of the above mentioned studies.

For space discretization of PDEs, commonly used methods are finite differences, finite elements, finite volume, and spectral methods. One of the most stable and accurate methods is the discontinuous Galerkin (dG) method for convection dominated PDEs like the FHN equations. In this paper, we use the symmetric interior penalty Galerkin (SIPG) method for space discretization [12, 13], which is the most relevant and popular dG method. The dG methods possess better properties for convective problems compared with continuous finite element methods, but they require a much larger number of degrees of freedom. Recently OCPs governed by linear stationary and time-dependent convection-diffusion-reaction equations [14, 15, 16] and semilinear steady state equations [17] are solved efficiently by the several dG methods. OCPs with PDE constraints are discretized usually following two approaches: the discretize-then-optimize approach, where the optimality conditions are imposed in the discrete setting, or the optimize-then-discretize approach in which the optimality conditions are derived in variational form and then discretization holds. There is no preferred approach [18, 19], here we follow the optimize-then-discretize approach: the optimality conditions are discretized in space by the SIPG method and the backward Euler method in time. We employ the projected nonlinear conjugate gradient (CG) method [20] for solving the resulting nonlinear discrete optimization problem.

The optimal control of instationary PDEs requires numerous evaluations of the optimality system, which lead to large scale optimization problems. In the last decade, different types of reduced order modeling (ROM) methods were developed to approximate these by low order or surrogate models. The resulting small optimization problems can be solved efficiently with less computational effort. For an overview of ROMs for OCPs for linear and semilinear PDEs, we refer the reader to [18, 19]. Despite its heuristic nature, currently, proper orthogonal decomposition (POD) is the preferred model reduction technique for linear and semilinear OCPs. The POD basis functions are constructed from the snapshots, solutions of space-time discretized PDE at pre-specified time-instances, then the reduced order OCP is solved applying the Galerkin projection. ROMs using the POD-Galerkin projection for OCPs for linear PDE constraints are investigated in [21, 22, 23] and for semi-linear PDE constraints in [24, 25, 26]. With increasing number of POD basis functions, more accurate solutions are obtained. But the success of POD depends on the type of the problem; for semi-linear PDEs with traveling or spiral wave type solutions like the Schlögl and FHN equations [1, 4, 27] require more POD modes than for the OCPs with semi-linear parabolic equations [24]. Because the reduced nonlinear term still depends on the dimension of the full order model (FOM), i.e., high dimensional finite element discretized model, hyper reduction techniques are developed. The empirical interpolation method (EIM) [28] and the discrete empirical interpolation method (DEIM) [29] are commonly used as hyperreduction techniques to reduce the computational cost for evaluation of the nonlinear terms. More recently, as an alternative to the DEIM, the dynamic mode decomposition (DMD) [30] is used to approximate the nonlinear terms of the PDEs. DMD was first introduced by Schmid [31], Rowley [32] and it is based on the linear approximation of the infinite dimensional Koopman operator [33, 34].

We compare the accuracy and computational efficiency of POD, POD with DEIM and POD with DMD for the optimal control of the convective FHN equations modeling blood coagulation and bio-reactors [35, 36]. The FHN equations are the simplest and most widely used model for describing the complex spatio-temporal behavior of traveling waves in excited media. In contrast to the classical FHN equations [37, 38] with a semi-linear PDE and with a linear ordinary differential equation (ODE), the convective FHN equations consists of a semi-linear PDE with the monotone cubic nonlinear term for the activator and a linear PDE for the inhibitor, modeling excitable systems in moving media. For semi-linear OCPs like the Burger’s equation, quadratic nonlinear terms can be reduced to the linear terms [39]. But for higher order polynomial or general nonlinearities, hyperreduction techniques should be used. As far as we know, reduced order OCPs using POD–DEIM and POD–DMD are not yet investigated for semilinear PDE constraint OCPs. Among these three methods, POD–DMD is the fastest, because after collecting snapshots, the non-linearity disappears and the reduced model becomes a coupled system of linear ODEs. Also, the reduced non-convex optimization problem for POD and POD–DEIM becomes convex for POD–DMD. The POD is the most accurate but, the slowest. The POD–DMD is less accurate than the POD but fastest.

In the sequel, we describe the OCP constrained by the convective FHN equations in Section 2. The fully discrete optimal control system is given in Section 3, based on the first order optimality conditions. In Section 4, the reduced order OCPs (POD, POD–DEIM, and POD–DMD) are constructed. In Section 5, we compare the three ROMs with respect to accuracy and computational efficiency for a test problem with terminal control. The paper ends with concluding remarks in Section 6.

2 Optimal control problem

We consider the following OCP:

[TABLE]

governed by the convective FHN equations

[TABLE]

where $d_{u}$, $d_{v}$ are the diffusion coefficients, $c_{3}$, $\epsilon$ are real constants, $x=(x_{1},x_{2})^{T}$ is the spatial element, $u_{T}(x)$ and $v_{T}(x)$ are given desired terminal states, and $u_{D}(x,t)$ and $v_{D}(x,t)$ are given functions. The control function $f(x,t)$ satisfy the point-wise box constraints

[TABLE]

for some real numbers $-\infty<f_{l}<f_{r}<+\infty$. The cubic polynomial $g(u)$ is given by

[TABLE]

where $c_{i}$, $i=1,2$, are non-negative real numbers. The non-linearity $g(u)$ is monostable for $0<c_{1}<20$ and $c_{2}=0.02$ [35], whereas for the Schlögl equation [1], the classical FHN equations [2], and the diffusive FHN equations [40] the bistable cubic non-linearity is considered. The velocity field $\bm{b}(x)=(b_{1}(x),b_{2}(x))^{T}$ is divergence-free, and for the maximum wave speed $b_{\max}$, it is given by the parabolic profile along the $x_{1}$-direction

[TABLE]

In (1) and (2), $T>0$ and $Q:=\Omega\times(0,T]$ denote the final time and the time-space domain, respectively. The space domain $\Omega=(0,L)\times(0,H)$ is bounded, Lipschitz in $\mathbb{R}^{2}$, and $\Sigma=\partial\Omega\times(0,T]$ is the lateral surface. The Dirichlet and Neumann boundaries are $\Sigma^{D}:=\Gamma_{D}\times(0,T]$ and $\Sigma^{N}:=\Gamma_{N}\times(0,T]$, where $\Gamma_{D}=\{x=\{0,L\},\;0<y<H\}$ and $\Gamma_{N}=\{0<x<L,\;y=\{0,H\}\}$. The outward unit normal vector is denoted by $\bm{n}(x)$ and the outward normal derivative on $\partial\Omega$ is defined as $\partial_{\emph{n}}u=\nabla u\cdot\bm{n}$.

The convective FHN equations (2) consists of an activator equation for $u(x,t)$ and an inhibitor equation for $v(x,t)$ in an excitable medium. In blood coagulation process, the activator variable $u(x,t)$ describes the concentrations of thrombin in the excitation, and $v(x,t)$ the inhibition of this excitation and recovery of the medium variable, which is activated factor XI [36]. The complex process of coagulation consists of cascadic enzymatic reactions and feedback loops [41]. These reactions enable generation of auto-catalytic thrombin outside of the damaged region. The most important property of the blood coagulation process is the formation of auto-waves with the velocity independent initial conditions [35, 36, 41]. The flow propagates within the impermeable channel walls (Neumann boundary conditions). Different types of standing and triggering waves occur depending on the constants of the convective FHN equations [35]. The waves become more curved when the flow velocity $b_{max}$ increases. For a detailed discussion of the complex spatio-temporal wave phenomena occurring in the uncontrolled convective FHN equations, we refer the reader to [35, 42].

The waves are initiated at $t=0$ inside an initial narrow excitation region $\{x_{a}\leq x_{1}\leq x_{b},\;0\leq x_{2}\leq H\}$

[TABLE]

The interval from $x_{a}$ to $x_{b}$ is the damaged region of the vessel wall. In blood coagulation, for fixed thrombin concentration the initial condition (3) mimics the activation process [35, 36, 41]. In this study, we take a rectangle as $x_{a}=0$, $x_{b}=0.1$, and we set $\bar{u}=0.1$.

Optimal control aims that the activator $u(x,t)$ and the inhibitor $v(x,t)$ follow desired terminal states $u_{T}(x)$ and $v_{T}(x)$ as close as possible in the $L_{2}$ norm by minimizing the cost functional $J(u,v,f)$. The inhibitor $v(x,t)$ has only some auxiliary character as for the FHN equations in [2, 3, 4], therefore only the activator $u(x,t)$ is to be controlled by $f(x,t)$, which seems to be experimentally feasible [43]. The parameter $\upsilon>0$ in the optimal cost functional $J(u,v,f)$ denotes the penalization or the Tikhonov regularization parameter.

Remark 2.1**.**

We have considered the OCP (1) with terminal observation functionals to show the applicability of reduced order modeling. Reduced solutions of OCPs with observation functionals, distributed in space–time for the FHN equations [8] can also be handled using the same approach.

3 Optimality system and discretization

The local optimal solutions are constructed using the first-order necessary optimality conditions. In [8], we consider a more general OCP constrained by the FHN equations with distributed, terminal and sparse controls. Here we consider the OCP (1) only with the terminal controls and without sparse controls. As in [8], we follow the optimize-then-discretize approach, first by deriving the first order optimality conditions in the variational form and then discretize the OCP. The second order necessary and sufficient optimality conditions for semi-linear parabolic equations are investigated in [2, 44]. The first-order necessary optimality conditions [44, Chapter 1], lead to the following adjoint system

[TABLE]

with the mixed boundary and terminal time conditions

[TABLE]

The flow in the adjoint system (4) is in the opposite direction of the flow in the state system (2). In view of the box constraints, the optimal solution $f(x,t)$ has to satisfy the following variational inequality

[TABLE]

which leads to the point-wise projection formula

[TABLE]

with the projection operator $\mathbb{P}_{[a,b]}(\phi)=\max\{a,\min\{\phi,b\}\}$.

3.1 Space discretization of the optimality system

The optimality system consisting of the state equation (2) and the adjoint equation (4) is discretized in space using the SIPG method [13]. We denote the family of meshes consisting of non-overlapping triangular elements $K$ by $\{\mathcal{T}_{h}\}_{h}$. The diameter of an element $K$ and the length of an edge $E$ are denoted by $h_{K}$ and $h_{E}$, respectively. We assume that the mesh is regular in the following sense: for different triangles $K_{i},\;K_{j}\in\{\mathcal{T}_{h}\}_{h},i\neq j$, the intersection $K_{i}\cap K_{j}$ is either empty or a vertex or an edge, i.e., hanging nodes are not allowed.

The discrete test, state, and control spaces are defined by the space of piecewise discontinuous finite element functions

[TABLE]

Here, the set of polynomials on $K\in\mathcal{T}_{h}$ of degree at most $l$ is denoted by $\mathbb{P}^{l}(K)$. In numerical examples, we use linear discontinuous finite elements, i.e., $l=1$.

The set $\mathcal{E}_{h}$ of all edges are split into the sets of interior, Dirichlet and Neumann boundary edges, denoted respectively by $\mathcal{E}_{h}^{0},\,\mathcal{E}_{h}^{D},\,\mathcal{E}_{h}^{N}$, so that $\mathcal{E}_{h}=\mathcal{E}_{h}^{0}\cup\mathcal{E}_{h}^{D}\cup\mathcal{E}_{h}^{N}$. The inflow and outflow boundaries are defined by

[TABLE]

Similarly, the inflow and outflow boundaries of an element $K\in\mathcal{T}_{h}$ are defined by

[TABLE]

where $\bm{n}_{K}(x)$ is the exterior unit normal vector on the boundary $\partial K$ of an element $K$. Let the common edge for two elements $K$ and $K^{e}$ be $E$. Then there are two traces of a function $w\in W_{h}$ along the edge $E$, denoted by $w|_{E}$ from the interior of $K$ and $w^{e}|_{E}$ from the interior of $K^{e}$. Accordingly, the jump and average of $w$ across the edge $E$ are defined by

[TABLE]

Similarly, the jump and average of a vector-valued function are given by

[TABLE]

We set \mbox{\left[!\left[w\right]!\right]}=w|_{E}\bm{n} and \mbox{\left{!!\left{\nabla w\right}!!\right}}=w|_{E} on a boundary edge $E\in K\cap\partial\Omega$. Then, the state equations (2) and the adjoint equations (4) are discretized using the SIPG method leading to the system: $\forall w\in W_{h}$ and a.e. $t\in(0,T]$

[TABLE]

[TABLE]

where $\left\langle{\cdot,\cdot}\right\rangle$ denotes the usual $L^{2}$-inner product. The bilinear forms $a_{h,u}$ and $a_{h,v}$ in the discrete state system (5) are given for $i=u,v$ by

[TABLE]

for any $\tilde{w},\;w\in W_{h}$, and $\gamma$ is called the penalty parameter [13]. For $i=u,v$, the linear right-hand side vectors, and the nonlinear form are defined by

[TABLE]

The bilinear forms $a_{h,p}$ and $a_{h,q}$ in the adjoint system (6) are defined similar to the bilinear forms $a_{h,u}$ and $a_{h,v}$ for states, respectively, but with the negative velocity field. They contain, respectively, the additional terms

[TABLE]

which come from the mixed boundary conditions. The SIPG discrete states and the control in (5) are given in the form

[TABLE]

In a similar way, the discrete desired terminal states (just $L^{2}$-projections) are given by

[TABLE]

The number of (triangular) elements are denoted by $N_{K}$, and the local dimension of each element is given by $N_{l}$. $\phi_{j}^{\,i}$ denotes the $j$-th finite element basis function defined on the $i$-th triangle. Setting dG degrees of freedom $N:=N_{K}\times N_{l}$, the corresponding time dependent unknown coefficient vectors can be written as

[TABLE]

On the other hand, the known coefficient vectors of the discrete desired terminal states are given by

[TABLE]

which are computed from the projection identities

[TABLE]

Then, the SIPG semi-discretized system of the state equations (5) leads to the following system of ODEs

[TABLE]

and the semi-discrete objective functional takes form

[TABLE]

Here $\bm{M}$ denotes the mass matrix, and $\bm{S}_{u}$ and $\bm{S}_{v}$ are the stiffness matrices corresponding to the bilinear forms $a_{h,u}$ and $a_{h,v}$, respectively. The nonlinear vector $\bm{g}(\bm{u})$ corresponds to the nonlinear form $g_{h}$, and the time-dependent vectors ${\bm{\ell}}_{u}:={\bm{\ell}}_{u}(t)$ and ${\bm{\ell}}_{v}:={\bm{\ell}}_{v}(t)$ are the vectors corresponding to the linear forms $\ell_{h,u}$ and $\ell_{h,v}$, respectively. By a similar setting, the SIPG semi-discretized system of the adjoint equations (6) are given as the following system of ODEs

[TABLE]

where the state $u(x,t)$ dependent matrix $\bm{R}_{g}(\bm{u})$ is related to the term $\left\langle{g^{\prime}(u_{h})p_{h},w}\right\rangle$ in the adjoint system (6), and $\bm{S}_{p}$ and $\bm{S}_{q}$ are the stiffness matrices corresponding to the bilinear forms $a_{h,p}$ and $a_{h,q}$, respectively.

3.2 Fully discrete optimality system

The time interval $[0,T]$ is partitioned uniformly as $0=t_{0}<t_{1}<\cdots<t_{N_{T}}=T$, with the step-size $\Delta t=T/N_{T}$. For $n=0,1,\ldots,N_{T}$, we denote by $\bm{u}_{n}$, $\bm{v}_{n}$, $\bm{f}_{n}$, $\bm{p}_{n}$ and $\bm{q}_{n}$ the approximate solution vectors to the semi-discrete solutions $\bm{u}(t)$, $\bm{v}(t)$, $\bm{f}(t)$, $\bm{p}(t)$ and $\bm{q}(t)$ at the time instance $t=t_{n}$, respectively. Also, we set $\bm{\ell}_{u}^{n}:=\bm{\ell}_{u}(t_{n})$ and $\bm{\ell}_{v}^{n}:=\bm{\ell}_{v}(t_{n})$. Then, applying the backward Euler method to the semi-discrete state equation (7) and to the semi-discrete adjoint equation (8), we obtain the following fully discrete state system

[TABLE]

for $n=1,2,\ldots,N_{T}$, and the following fully discrete adjoint system

[TABLE]

for $n=N_{T},\ldots,2,1$, where $\bm{u}_{T}$ and $\bm{v}_{T}$ are the known coefficient vectors of the desired terminal states. We note that the fully discrete state system (9) is solved forward in time, whereas the fully discrete adjoint system (10) is solved backward in time. The equations (9) and (10) are semi-explicit in time. Because the mass matrix ${\mathbf{M}}$ is symmetric positive definite and therefore invertible, the discrete solutions of (9) and (10) exist and are unique. The resulting semilinear discrete OCP can be solved by several optimization algorithms. We use the projected nonlinear CG method [20], which is applied for solving the OCP governed by Schlögl and FHN equations in [1, 2, 3, 4] and by the convective FHN equations in [8]. The details of the implementation of the projected nonlinear CG algorithm can be found in these papers.

4 Reduced order optimal control

We construct ROMs for the state equation (2) and the adjoint equation (4). For the construction of the POD modes to the state variables $u$ and $v$, the POD method uses the snapshot matrices generated by the coefficient vectors $\bm{u}$ and $\bm{v}$ of the optimal state solution vectors from the full order optimal control problem, and similarly, the snapshot matrix generated by the coefficient vector $\bm{f}$ of the optimal control solution is used to construct the POD modes to the control variable $f$. For the construction of the POD modes to the adjoint variables $p$ and $q$, we use the same POD modes generated by the snapshots of the states $u$ and $v$, motivated by the error analysis in [19]. This approach might be not the best option, but the construction of different POD modes using the snapshot matrices of the adjoint states does not improve the suboptimal solutions much [24] and requires more computational work. It is not known a-priorly how far the optimal solution of the reduced POD problem is from the exact one unless the snapshots are generating a sufficiently large state space. Different approaches are discussed in [18] for the proper choice of the POD basis in the reduced OCP, which affects the accuracy of the reduced solutions. Here we construct the POD basis for the reduced OCP from the snapshots of controls of the FOM.

4.1 POD Galerkin discretization

The reduced order system for the state equations (2) of lower dimensions $k_{u}$ and $k_{v}$ is formed by the Galerkin projection of the equations onto the subspaces

[TABLE]

For any time $t$, the reduced order states $u_{h}^{r}(x,t)\in W_{u,h}^{r}$ and $v_{h}^{r}(x,t)\in W_{v,h}^{r}$ stand for the approximations of the full order states $u_{h}(x,t)$ and $v_{h}(x,t)$, respectively, and they are given by

[TABLE]

where $\bm{u}^{r}(t):=(u^{r}_{1}(t),\ldots,u^{r}_{k_{u}}(t))^{T}$ and $\bm{v}^{r}(t):=(v^{r}_{1}(t),\ldots,v^{r}_{k_{v}}(t))^{T}$ are the reduced coefficient vectors, and $\{\psi_{u,i}\}_{i=1}^{k_{u}}$ and $\{\psi_{v,i}\}_{i=1}^{k_{v}}$ are the $L^{2}$-orthogonal reduced basis functions. Similarly, the full order control $f_{h}(x,t)$ is approximated by the reduced order control $f_{h}^{r}(x,t)$ from the $k_{f}$ dimensional subspace $W_{f,h}^{r}$

[TABLE]

where $\bm{f}^{r}(t):=(f^{r}_{1}(t),\ldots,f^{r}_{k_{f}}(t))^{T}$ is the reduced coefficient vector for the control, and $\{\psi_{f,i}\}_{i=1}^{k_{f}}$ are the $L^{2}$-orthogonal reduced basis functions. All the reduced basis functions are linear combinations of dG basis functions

[TABLE]

where $\Psi_{\cdot,\cdot,i}$ are the coefficient vectors of the $i$-th reduced basis functions. For $w\in\{u,v,f\}$, the following POD matrices are constructed

[TABLE]

where the columns $\Psi_{w,i}$ are the POD modes. The POD modes are computed through the snapshot matrices $\mathcal{U}=[\bm{u}_{1},\ldots,\bm{u}_{N_{T}}]$, $\mathcal{V}=[\bm{v}_{1},\ldots,\bm{v}_{N_{T}}]$ and $\mathcal{F}=[\bm{f}_{1},\ldots,\bm{f}_{N_{T}}]$ in $\mathbb{R}^{N\times N_{T}}$ [40, 45], respectively, where the columns of the snapshot matrices are simply the coefficient vectors of the full order states and control from the FOM (9) at the time instances $t_{i}$, $i=1,\ldots,N_{T}$. Then, for $w\in\{u,v,f\}$, the reduced basis functions $\{\psi_{w,i}\}$, $i=1,2,\ldots,k_{w}$, are given by the solution of the following minimization problem

[TABLE]

where $\delta_{ij}$ is the Kronecker delta. Again for $w\in\{u,v,f\}$, and for $\mathcal{W}\in\{\mathcal{U},\mathcal{V},\mathcal{F}\}$, the above minimization problem is equivalent to the eigenvalue problems

[TABLE]

for the coefficient vectors $\Psi_{w,\cdot,i}$ of the reduced basis functions $\psi_{w,i}$. Setting $\widehat{\mathcal{W}}=R\mathcal{W}$ ($R^{T}$ is the Cholesky factor of the mass matrix ${\bm{M}}$), we obtain the equivalent formulation of (11) as

[TABLE]

where $\widehat{\Psi}_{w,\cdot,i}=R\Psi_{w,\cdot,i}$. Because the singular value decomposition (SVD) is more stable and efficient than the eigenvalue decomposition, we apply the SVD to compute the first $k_{w}$ left singular vectors, $\widehat{\Psi}_{w,\cdot,i}:=\zeta_{w,i}$, of the matrix $\widehat{\mathcal{W}}$

[TABLE]

where the diagonal matrix $\Sigma_{w}$ contains the singular values $\sigma_{w,i}$, and $\zeta_{w,i}$ are the columns of the orthogonal matrix $\zeta_{w}$. Using the fact that $\widehat{\Psi}_{w,\cdot,i}=R\Psi_{w,\cdot,i}$, the coefficient vectors $\Psi_{w,\cdot,i}$ of the reduced basis functions are computed as

[TABLE]

By the above settings, the following relations can be found between the coefficient vectors $\bm{u}$, $\bm{v}$ and $\bm{f}$ of the FOM solutions, and the reduced coefficient vectors $\bm{u}^{r}$, $\bm{v}^{r}$ and $\bm{f}^{r}$

[TABLE]

Then the reduced system for the state equations takes the form

[TABLE]

and using the $L^{2}$-orthogonality of the reduced basis functions, the semi-discrete reduced order objective functional takes form

[TABLE]

which is absolutely cheaper to calculate, comparing to the full order objective functional $J_{h}$. Using the same POD modes generated by the snapshot matrix of the states, we can obtain the reduced adjoint system in a similar way

[TABLE]

In the ROMs (12) and (13), the reduced matrices are defined by

[TABLE]

The sub-optimality systems (12) and (13) are solved by backward Euler method in time.

4.2 Discrete empirical interpolation method

In this part, we rewrite the reduced state equation (12) resulting from the application of the DEIM [29] to the cubic nonlinear term $\bm{g}(\bm{\Psi}_{u}\bm{u}^{r})$. The DEIM aims to approximate the nonlinear vector $\bm{g}(\bm{\Psi}_{u}\bm{u}^{r})$ by projecting it onto a subspace of the space generated by the non-linear functions and spanned by a basis of dimension $m\ll N$

[TABLE]

where $\bm{s}(t)\in\mathbb{R}^{m}$ is the corresponding coefficient vector, and the DEIM basis matrix $\bm{W}\in\mathbb{R}^{N\times m}$ is calculated by the application of the POD to the nonlinear snapshot matrix $\mathcal{G}=[\bm{g}(\bm{u}_{1})\cdots\bm{g}(\bm{u}_{N_{T}})]\in\mathbb{R}^{N\times N_{T}}$. Since the system (14) is overdetermined, we find a projection matrix $\bm{P}$ such that $\bm{P}=[\bm{e}_{\mathfrak{p}_{1}}\cdots\bm{e}_{\mathfrak{p}_{m}}]\in\mathbb{R}^{N\times m}$ where $\bm{e}_{\mathfrak{p}_{i}}$ is the $i$-th column of the identity matrix $\bm{I}\in\mathbb{R}^{N\times N}$. After elimination of the coefficient vector $\bm{s}(t)$ in (14), we arrive at the following ROM for the state equation with the DEIM approximation as follows

[TABLE]

where the matrix $\bm{Q}:=\bm{\Psi}_{u}^{T}\bm{W}(\bm{P}^{T}\bm{W})^{-1}$ is computed once in the off-line stage. The computation of the reduced nonlinear vector $\bm{g}^{r}(\bm{u}_{r}):=\bm{P}^{T}\bm{g}(\bm{\Psi}_{u}\bm{u}^{r})\in\mathbb{R}^{m}$ and its Jacobian matrix require only $m\ll N$ and $m\times N_{l}$ integral evaluations, respectively, while they need $N$ and $N\times N_{l}$ integral evaluations when DEIM approximation is not used.

4.3 Dynamic mode decomposition

The DMD extracts dynamically relevant spatio-temporal information content from numerical or experimental data sets [34]. Without explicit knowledge of the dynamical system, the DMD algorithm determines eigenvalues, eigenmodes, and spatial structures for each mode. The snapshot data is decomposed in spatio-temporal modes both by POD and DMD that correlates the spatial features of the data. Additionally, DMD associates the snapshot data to the temporal Fourier modes. DMD is a special case of the Koopman operator [33] approximating nonlinear systems via an associated infinite dimensional system. The connection between the DMD and Koopman operator was established in [32, 31] . The Koopman operator ${\mathcal{K}}$ acts on a set of scalar observable functions $g:{\mathcal{M}}\rightarrow{\mathbb{C}}$

[TABLE]

for the nonlinear dynamical system

[TABLE]

where ${\bm{y}}\in{\mathcal{M}}$, and ${\mathcal{M}}$ is an n-dimensional manifold. The DMD determines the Koopman eigenvalues and modes directly from the data, when linear observables are considered as state space, $g({\bm{y}})={\bm{y}}$, directly from the data.

We consider a snapshot matrix and a time shifted version of the snapshot matrix

[TABLE]

The DMD involves the decomposition of the unknown best-fit linear operator ${\mathbf{A}}\in\mathbb{R}^{n\times n}$ relating the matrices above

[TABLE]

The DMD modes are computed by the exact DMD algorithm [46].

In [30], the DMD is proposed for producing low-rank approximations of the nonlinearities in PDEs as an alternative to DEIM. Here, after building the POD basis functions of rank $k_{u}$, we collect snapshots ${\mathcal{G}},{\mathcal{G}}^{\prime}$ for the non-linear term ${\bm{g}}(\bm{u})$ from the FOM (7) at $N_{T+1}$ time instances and divide them into sets as in (16)

[TABLE]

The snapshots matrices satisfy

[TABLE]

where the unknown matrix ${\bm{A}}_{\mathcal{G}}$ is given as the solution of the minimization problem

[TABLE]

where $\|\cdot\|_{F}$ denotes the Frobenius norm, ${\bm{A}}_{\mathcal{G}}={\mathcal{G}}^{\prime}{\mathcal{G}}^{{\dagger}}$, and ${\dagger}$ is the Moore-Penrose pseudoinverse. The Algorithm 1 gives the exact DMD procedure [46] to find the DMD basis matrix $\bm{\Psi}^{\text{DMD}}=[\bm{\Psi}^{\text{DMD}}_{1}\cdots\bm{\Psi}^{\text{DMD}}_{\widetilde{m}}]\in\mathbb{R}^{N\times\widetilde{m}}$. After calculating the DMD basis function $\bm{\Psi}^{\text{DMD}}$ by the Algorithm 1, the nonlinear vector ${\bm{g}}(\bm{u})$ in the full order state equations (7) is approximated as the following time-dependent vector

[TABLE]

where $\alpha=[\alpha_{1},\ldots,\alpha_{\widetilde{m}}]^{T}$ is the initial amplitudes given by $\alpha=(\bm{\Psi}^{\text{DMD}})^{{\dagger}}{\bm{g}}(\bm{u}_{1})$, and $\omega=[\omega_{1},\ldots,\omega_{\widetilde{m}}]^{T}$ includes the eigenvalues $\lambda_{j}$ as $\omega_{j}=\log{(\lambda_{j})}/\Delta t,j=1,\ldots\widetilde{m}$.

Inserting the identity (17) into the reduced state equations (12), we obtain the following linear ROM for the activator equation

[TABLE]

Although the dimension of the reduced system (18) with $\widetilde{m}$ DMD modes is small as the reduced system (15) with $m$ DEIM modes, the POD–DMD reduced state equation (18) is linear, and the OCP problem becomes convex and Newton’s iterations are no more needed. Therefore, the POD–DMD is much faster than the POD and the POD–DEIM.

5 Numerical results

We consider the OCP (1)-(2) with a fast wave speed $V_{max}=64$, and with the parameters and initial condition given by

[TABLE]

[TABLE]

The desired terminal states are chosen as

[TABLE]

where $u_{\mathrm{nat}}$ and $v_{\mathrm{nat}}$ stand for the solutions of the uncontrolled convective FHN equations. The box constraint for the control is given by

[TABLE]

The space domain is taken as a rectangle with $H=5$ and $L=100$. We use uniform step size in space, $\Delta x_{1}=\Delta x_{2}=0.5$, and in time, $\Delta t=0.05$ for the final time $T=1$. Stopping criteria for the OCP is the relative error $|J_{old}-J|/|J_{old}|\leq 10^{-3}$ on the objective functional. For $w\in\{u,v,f\}$, the number of POD modes are determined according to the following energy criteria (relative information content (RIC) )

[TABLE]

which represents the energy captured by the first $k_{w}$ POD modes over all $d_{w}$ POD modes, $d_{w}$ is the rank of the related snapshot matrix, and $\sigma_{w,i}$ is the corresponding singular value of the $i$-th mode. In the computations, the number of POD modes are taken as the same, $k:=\max\{k_{u},k_{v},k_{f}\}$, so that the energy criteria $\min_{k_{w}}\epsilon(k_{w})\geq\%99.99$ is satisfied for $w\in\{u,v,f\}$. Moreover, for $w\in\{u,v,f\}$, we define the following relative $L^{2}$-error

[TABLE]

to compute the errors between full and reduced order solutions at a given time $t$.

The singular values of the snapshot matrices of the states $u$, $v$, of the control $f$, and of the non-linearity $g(u)$ in Figure 1 decay slowly almost at the same rate, which is typical for convection dominated problems. According to the RIC for $\%99.99$ of energy captured, we have chosen $k=9$ POD modes for the states $u$, $v$ and the control $f$.

The FOM state solutions at the final time in Figures 2-3, top, show the same wave type characteristics as in [8]. The error plots between full and reduced order state solutions in Figures 2-3, bottom, agree for POD, POD–DEIM and POD–DMD with $m=14$ and $\widetilde{m}=18$ DEIM and DMD modes, respectively. In Figure 4, we give the full and reduced order control profiles at the final time. We see that there is a bit difference between the full and reduced order profiles, but the reduced order profiles are similar in between. The relative $L^{2}$-errors between the full and reduced order solutions of the states and the control at the final time are shown in Table 1. Both states $u$ and $v$ are approximated by the reduced models with the same accuracy, whereas the errors in control are larger for POD–DMD due to linearization of the reduced OCP.

Because the Galerkin projected POD–DMD system is linear, the resulting OCP is convex. Therefore, no Newton iteration is recorded for the POD–DMD in Table 2, and instead of nonlinear CG, we use linear CG method. Comparing the speedup factors (the ratio of the CPU times (in seconds) required for evaluation of FOM to the CPU times required for evaluation of ROM) in Table 2, the computational efficiency of the POD–DMD approach is clearly visible. The POD reduced sub-optimal solutions are the most accurate, the POD–DEIM and POD–DMD solutions are less accurate, whereas the POD solutions oscillate at higher modes, as shown for the relative FOM-ROM errors in the Frobenius norm in Figure 5, and in Table 1 as well . The value of the objective function $J$ obtained by full and reduced order solutions are almost the same with increasing number of POD modes in Figure 6, left. The ROMs require significantly less CPU time compared to the FOM. In terms of computational cost, the POD–DMD is the fastest, Figure 6, right. For a small number of POD modes, the CPU times oscillate until $k\approx 6$, and after they do not change much. Therefore it is not necessary to use a larger number of POD modes considering the CPU times, and the accuracy of the reduced order solutions and reduced order objective functional.

6 Conclusions

In this study, the computational efficiency and accuracy of three ROM methods are investigated for solving semilinear convection dominated OCP problem. We derived the ROMs for the state and adjoint equations using the same POD basis functions generated by the snapshot matrices of the states only. Two hyperreduction techniques, the DEIM and DMD are applied to reduce the computational cost arising from the nonlinear term in the activator equation. Among the three ROM techniques the POD without DEIM/DMD is the most accurate as expected. The POD–DEIM and POD–DMD errors are close, but the POD–DMD is the fastest due to the fact that the reduced state equations are no more nonlinear and the OCP is convex. In a future study, we plan to compare these ROM techniques for model predictive control and feedback control problems using compressive POD and DMD.

Acknowledgements The authors would like to thank the reviewers for the comments and suggestions that helped to improve the manuscript.