Model order reduction for optimality systems through empirical gramians

In the present article, optimal control problems for linear parabolic partial differential equations (PDEs) with time-dependent coefficient functions are considered. One of the common approach in literature is to derive the first-order sufficient optimality system and to apply a finite element (FE) discretization. This leads to a specific linear but high-dimensional time variant (LTV) dynamical system. To reduce the size of the LTV system, we apply a tailored reduced order modeling technique based on empirical gramians and derived directly from the first-order optimality system. For testing purpose, we focus on two specific examples: a multiobjective optimization and a closed-loop optimal control problem. Our proposed methodology results to be better performing than a standard proper orthogonal decomposition (POD) approach for the above mentioned examples.

1 Introduction

Optimization problems constrained by time-dependent partial differential equations (PDEs) arise in many fields of application in engineering and across all sciences. Examples of such problems include optimal (material) design or optimal control of processes and inverse problems, where parameters of a PDE model are unknown and need to be estimated from measurements. The numerical solution of such problems is very challenging as the underlying PDEs have to be solved repeatedly within outer optimization algorithms, and the dimension of the parameters that need to be optimized might be very high or even infinite dimensional, especially when one is interested in multiobjective optimization or closed-loop optimal control. In a classical approach, the underlying PDE (forward problem) is approximated by a high dimensional full order model (FOM) that results from discretization. For the spatial discretization, often a finite element (FE) method is used, leading to high dimensional dynamical systems. Hence, the complexity of the optimization problem directly depends on the numbers of degrees of freedom (DOF) of the FOM. Mesh adaptivity has been advised to minimize the number of DOFs (see, e.g., [1, 2]).

A more recent approach is the usage of model order reduction (MOR) methods in order to replace the FOM by a surrogate reduced order model (ROM) of possibly very low dimension. MOR is a very active research field that has shown tremendous development in recent years, both from a theoretical and application point of view. For an introduction and overview, we refer to the monographs and collections [3–6]. In the context of optimal control, ROM is utilized [7–11]. In multiobjective optimization (cf., e.g., [12, 13]) and model predictive control (cf., e.g., [14, 15]), the situation is even more complex because many optimization problems have to be solved for varying data functions.

An appropriate ROM for the optimal control problem governed by high dimensional dynamical systems has to guarantee not only a sufficiently accurate surrogate model for the dynamical system but also, in particular, provide an acceptable suboptimal control. For that reason, MOR should inheritate properties such as controllability and observability of the ROM for the dynamical system (see [9, 16–19]). In the context of linear time-invariant (LTI) systems, it is well known that balanced truncation is an MOR strategy that ensures controllability and observability of the reduced LTI system under reasonable conditions [7, 9]. However, for linear time-variant (LTV) systems, balanced truncation cannot be applied directly. Here, the empirical gramian approach offers a promising extension [20], where an empirical controllability gramian and an empirical observability gramian are computed by simulations of the dynamical system for different impulse controls and different initial conditions, respectively. Then, the MOR can be derived from a suitable singular value decomposition; cf. Section 3. Recalling that for LTI systems, this approach is equivalent to balanced truncation (see [20, 21] for more details). In this study, we consider optimal control problems governed by linear parabolic PDE constraints with transient inhomogeinities or coefficients. After an FE discretization, we derive an LTV system of the form as follows:

$\begin{array}{ll}ẏ(t)=A(t)y(t)+B(t)u(t)+f(t)\text{ for }t\in (0,T],y(0)=y_{\circ },& (\mathrm{1a})\end{array}$

$\begin{array}{ll}z(t)=C(t)y(t)\text{ for }t\in (0,T].& (\mathrm{1b})\end{array}$

Let us mention that the empirical gramian approach can be also used for general non-linear dynamical systems. However, to get a low-dimensional ROM for a non-linear dynamical system, an efficient realization of the non-linear term is required, for instance, by applying the (discrete) empirical interpolation method (cf. [22–24]) or missing point evaluation (cf. [25]). For this reason, in this study, we focus on LTV systems and quadratic cost functionals. The optimal solutions can be characterized by a linear first-order optimality system due to convexity. This first-order optimality system can be observed as a extended coupled linear LTV system in the state and the adjoint variable. The new contribution of the present study is the development of a ROM, which is tailored to the optimization problem by utilizing empirical gramians computed from solutions of the coupled LTV system that describes the first-order necessary optimality conditions. The parameters of the state equation (different convection functions) and the cost functional (different desired states) lead to different inputs for the simulations required to compute the empirical gramians. It turns out that the obtained ROM is more reliable and robust regarding changes in the parameters than a multiple snapshot proper orthogonal decomposition (POD) method [11]. Furthermore, our MOR approach is certified by a-posteriori error estimates for the controls.

The study is organized as follows: In Section 2 we recall the empirical gramian framework and explain how this framework can be utilized to get a ROM for (1). Our empirical gramian approach is stated in Section 3. Based on the first-order optimality system, an MOR is computed. In Section 4, we first recall the multiple snapshot POD method. Moreover, the multiobjective optimal control problem and the closed-loop control are discussed. We presented some conclusions in Section 5. Finally, the a-posteriori error estimate is briefly presented in the Appendix.

2 Empirical gramians

In this section, we explain the concepts of empirical gramians and how they are used to perform model order reduction for linear time-variant input—output systems of the form (1). We mainly follow the study by [26] and suppose that T = ∞ holds. Throughout this study, we make use of the following hypothesis.

Assumption 1. For given control $u\in L^{\infty }(0,\infty ;{ℝ}^{n_u})$, the state and the output of system (1) satisfy $y\in L^{\infty }(0,\infty ;{ℝ}^{n_y})$ and $z\in L^{\infty }(0,\infty ;{ℝ}^{n_z})$, respectively.

Remark 1. Assumption 1 is a standard assumption in the context of dynamical systems. As also noted in the study mentioned in [20], this assumption is generally satisfied for stable linear systems and exponentially stable non-linear systems.      ◇

First, gramian-based methods are introduced for linear time-invariant (LTI) systems when

$\begin{array}{ll}A(t)=A,B(t)=B,C(t)=C\text{and }f(t)=0\text{ for }t>0& (2)\end{array}$

hold with constant matrices $A\in {ℝ}^{n_y\times n_y}$, $B\in {ℝ}^{n_y\times n_u}$, and $C\in {ℝ}^{n_z\times n_y}$ (see [27]). We recall the following definition from the study mentioned in [18] (Definitions 3.1 and 3.4).

Definition 2. The LTI system

$\begin{array}{ll}\begin{array}{c}ẏ(t)=Ay(t)+Bu(t)\text{ for }t>0,x(0)=y_{\circ },\\ z(t)=Cy(t)\text{ for }t>0\hfill \end{array}& (3)\end{array}$

is called controllable if, for any initial state $y_{\circ }\in {ℝ}^{n_y}$, final time T > 0 and final state $y_T\in {ℝ}^{n_y}$; there exists a piecewise continuous input $u:\left[0,T\right]\to {ℝ}^{n_u}$ such that the solution of (3) satisfies y(T) = y_T. Otherwise, the LTI system is said to be uncontrollable. The LTI system (3) is said to be observable, if, for any T > 0, the initial condition y_° can be determined from the knowledge of the input u(t) and the output y(t) for t ∈ [0, T]. Otherwise, (3) is called unobservable.

Suppose that the LTI system is stable, i.e., all eigenvalues of A have strictly negative real part; cf. [18] (Definition 3.1). In that case, the linear controllability gramian is defined as the symmetric matrix

$\begin{array}{l}{L}_c:=\mathcal{C}{\mathcal{C}}^{*}={\displaystyle {\int }_0^{\infty }}{e}^{At}B{B}^{\top }{e}^{A^{\top }t}\text{d}t={\displaystyle {\int }_0^{\infty }}(e^{At}B)(e^{At}B{)}^{\top }\text{d}t\in {ℝ}^{n_y\times n_y},\end{array}$

where $\mathcal{C}:L^2(0,\infty )\to {ℝ}^{n_y}$, $u\mapsto {\int }_0^{\infty }{e}^{At}Bu(t)\mathrm{\text{d}}t$, is the (linear) controllability operator. Throughout the study, the symbol “⊤” stands for the transpose of a vector or matrix. Since the LTI system is stable, the operator $\mathcal{C}$ is bounded and its (Hilbert space) adjoint ${\mathcal{C}}^{*}:{ℝ}^{n_y}\to L^2(0,\infty )$ satisfies

$\begin{array}{l}{\langle {\mathcal{C}}^{*}v,\phi \rangle }_{L^2(0,T)}={\langle v,\mathcal{C}\phi \rangle }_{{ℝ}^{n_y}}\text{ for }v\in {ℝ}^{n_y}\text{and}\phi \in L^2(0,T).\end{array}$

Moreover, the linear observability gramian is given as the symmetric matrix

$\begin{array}{l}{L}_o:={\mathcal{O}}^{*}\mathcal{O}={\displaystyle {\int }_0^{\infty }}{e}^{A^{\top }t}{C}^{\top }C{e}^{At}\text{d}t=\\ {\displaystyle {\int }_0^{\infty }}(e^{A^{\top }t}{C}^{\top })(e^{A^{\top }t}{C}^{\top }{)}^{\top }\text{d}t\in {ℝ}^{n_y\times n_y},\end{array}$

where $\mathcal{O}:{ℝ}^{n_y}\to L^2(0,\infty ;{ℝ}^{n_z})$, x ↦ Ce^Atx is the (linear and bounded) observability operator. Recalling that L_c and L_o have the following properties (cf., e.g., [18, 28]; Section 3.8):

Lemma 3. Let us consider the LTI system (3). We assume that all eigenvalues of A have strictly negative real part. Then, the controllability gramian L_c is a positive semidefinite solution to the Lyapunov equation

$\begin{array}{ll}A{L}_c+L_c{A}^{\top }+B{B}^{\top }=0.& (4)\end{array}$

If $\mathrm{\text{Im}}(\mathcal{C})={ℝ}^{n_y}$, (4) admits a unique solution L_c which has full rank n_y and is positive definite. Moreover, the observability gramian L_o is a positive semidefinite solution to the Lyapunov equation

$\begin{array}{ll}{A}^{\top }{L}_o+L_{o}A+C^{\top }C=0.& (5)\end{array}$

If $ker(\mathcal{O})=\left\{0\right\}$, system (5) admits a unique solution L_o. Moreover, L_o is positive definite.

The matrices L_c and L_o contain essential information which states that the dynamical system is controllable and observable, respectively. This is utilized in model order reduction approaches such as balanced truncation, where states which are not or only less controllable and observable are neglected; cf., e.g., [9, 29, 30]. To combine both information, the so-called cross gramian matrix has been introduced; cf., e.g., [31–33]. Here, we have to assume that the input and output dimensions are the same. Thus, we have n_u = n_z and define the cross gramian matrix

$\begin{array}{l}{L}_x:={\displaystyle {\int }_0^{\infty }}{e}^{At}BC{e}^{At}\text{d}t={\displaystyle {\int }_0^{\infty }}(e^{At}B)(e^{A^{\top }t}{C}^{\top }{)}^{\top }\text{d}t\in {ℝ}^{n_y\times n_y}.\end{array}$

Unfortunately, the previous relationships cannot be used in the time variant case (or even in the more general non-linear case). In the case of linear time variant systems, both gramians are functions of two variables: the initial time moment and the final time moment. Unlike linear time-invariant systems, the difference in time moments alone does not uniquely characterize the gramian matrices. If one changes the initial time moment while maintaining the interval length for both gramians, the resulting solutions will also differ. This discrepancy arises because the system description evolves over time, leading to changes in system characteristics. Consequently, the shifted system will differ from the original one. Therefore, in the literature, various methods exist for computing gramians in the context of linear time variant systems. A general extension is a data-driven approach, where tentative candidates are constructed by subsequent simulations of the model. That is the reason why such candidates are called empirical gramians. This concept is introduced in the study by [30] and extended in the study by [20]. It turns out that we get a data-driven approach because the empirical gramians can be computed from measured or computed trajectories of the dynamical system.

Definition 4. For given k ∈ ℕ let $I_k\in {ℝ}^{k\times k}$ be the identity matrix, we denote by $\mathcal{T}$^k an arbitrary set of n_{$\mathcal{T}$} orthogonal matrices, i.e.,

$\begin{array}{l}{\mathcal{T}}^k:=\{T_1,\dots ,T_{n_{\mathcal{T}}}|T_i\in {ℝ}^{k\times k},T_i^{\top }{T}_i=I_k,i=1,\dots ,n_{\mathcal{T}}\}\\ \subset {ℝ}^{k\times k}\end{array}$

with $\mathcal{M}$ an arbitrary set of s positive constants, i.e.,

$\begin{array}{l}\mathcal{M}:=\{c_1,\dots ,c_{n_{\mathcal{M}}}|c_i\in ℝ,c_i>0\text{ for }i=1,\dots ,n_{\mathcal{M}}\}\subset ℝ,\end{array}$

and with $\mathcal{E}$^k ⊂ ℝ^k the set of standard unit vectors $e_1,\dots ,e_k\in {ℝ}^k$. Furthermore, given n ∈ ℕ and a function w ∈ L^∞(0, ∞; ℝⁿ), we define the mean $\bar{w}\in {ℝ}^n$ as

$\begin{array}{l}\bar{w}:={\displaystyle \underset{T\to \infty }{lim}}\frac{1}{T}{\displaystyle {\int }_0^T}w(t)\text{d}t\in {ℝ}^n\end{array}$

provided the limit exists.

Definition 5 (see [20, Definition 4]). Suppose that the sets ${\mathcal{T}}^{n_u}$, $\mathcal{M}$, and ${\mathcal{E}}^{n_u}$ are given as in Definition 4 for k = n_u. For (1) the empirical controllability gramian ${\hat{L}}_c\in {ℝ}^{n_y\times n_y}$ is defined as follows:

$\begin{array}{l}{\hat{L}}_c:={\displaystyle \sum _{i=1}^{n_{\mathcal{T}}}}{\displaystyle \sum _{j=1}^{n_{\mathcal{M}}}}{\displaystyle \sum _{l=1}^{n_u}}\frac{1}{c_j^2{n}_{\mathcal{T}}{n}_{\mathcal{M}}}{\displaystyle {\int }_0^{\infty }}{Y}^{ijl}\left(t\right)\text{ d}t\text{ with }{Y}^{ijl}\left(t\right)=\\ (y^{ijl}\left(t\right)-{ȳ}^{ijl})(y^{ijl}\left(t\right)-{ȳ}^{ijl}{)}^{\top },\end{array}$

where $y^{ijl}(t)\in {ℝ}^{n_y}$ solves (1) corresponding to the impulse input u(t) = c_jT_ie_lδ(t) and ȳ^ijl(t) stands for its mean.

Remark 6. In our numerical example, we utilize the following weighted inner product in the state space ${ℝ}^{n_y}$ induced by the symmetric positive definite matrix $W\in {ℝ}^{n_y\times n_y}$:

$\begin{array}{ll}{\langle y,ỹ\rangle }_W=y^{\top }Wỹ,{\Vert y\Vert }_W=\sqrt{{\langle y,y\rangle }_W}\text{ for }y,ỹ\in {ℝ}^{n_y}.& (6)\end{array}$

In Definiton 5, the matrices Y^ijl(t) have to be replaced by (y^ijl(t) − ȳ^ijl)(y^ijl(t) − ȳ^ijl)^⊤W.      ◇

Definition 7 (see [20, Definition 6]). Let the sets ${\mathcal{T}}^{n_y}$, $\mathcal{M}$, and ${\mathcal{E}}^{n_y}$ be defined as in Definition 4 for the choice k = n_u. For (1) the empirical observability gramian, ${\hat{L}}_o\in {ℝ}^{n_y\times n_y}$ is defined as

$\begin{array}{l}{\hat{L}}_o:={\displaystyle \sum _{i=1}^{n_{\mathcal{T}}}}{\displaystyle \sum _{j=1}^{n_{\mathcal{M}}}}\frac{1}{c_j^2{n}_{\mathcal{T}}{n}_{\mathcal{M}}}{\displaystyle {\int }_0^{\infty }}{T}_i{Z}^{ij}\left(t\right)T_i^{\top }\text{d}t\text{with }{Z}_{\nu \mu }^{ij}\left(t\right)=\\ (z^{ij\nu }\left(t\right)-{\bar{z}}^{ij\nu }{)}^{\top }(z^{ij\mu }\left(t\right)-{\bar{z}}^{ij\mu }),\end{array}$

where $z^{ij\nu }(t)\in {ℝ}^{n_z}$, ν = 1, …, n_y, is the output of (1) corresponding to the initial condition y_° = c_jT_ie_ν, and ${\bar{z}}^{ij\nu }$ denotes its mean.

Remark 8. 1) In our numerical examples, we have n_y = n_z so that the output space is also supplied by the weighted inner product 〈·, ·〉_W introduced in (6). Then, we replace the matrix elements $Z_{\nu \mu }^{ij}(t)$ by ${(z^{ij\nu }(t)-{\bar{z}}^{ij\nu })}^{\top }W(z^{ij\mu }(t)-{\bar{z}}^{ij\mu })$.

2) Recalling that the empirical gramians reconstruct the true gramians in the case of LTI systems; cf. [20] (Lemmas 5 and 7). In all the other cases, they can be observed as an extension of the gramian notion.

3) In [26] (Section 3.1.3), also an empirical cross gramian is introduced. Since in our numerical examples, the results based on empirical cross gramians are not as good as the ones using a combination of empirical controllability and observability gramians, we skip the empirical cross gramian here.     ◇

Given the symmetric empirical gramians ${\hat{L}}_c$ and ${\hat{L}}_o$, it is possible to define a reduced order model through a balancing transformation which is based on the singular value decomposition (SVD). We follow the approach proposed in the study by [21]. First, we compute the SVD of the symmetric matrices ${\hat{L}}_c$ and ${\hat{L}}_o$:

$\begin{array}{ll}{\hat{L}}_c=U_c{\Sigma }_c{U}_c^{\top }\in {ℝ}^{n_y\times n_y}\text{and }{\hat{L}}_o=U_o{\Sigma }_o{U}_o^{\top }\in {ℝ}^{n_y\times n_y}& (\mathrm{7a})\end{array}$

with orthogonal matrices U_c, $U_o\in {ℝ}^{n_y\times n_y}$ and diagonal matrices Σ_c, ${\Sigma }_o\in {ℝ}^{n_y\times n_y}$ containing the non-negative singular values in descending order. Then, we can compute the matrices ${\hat{L}}_c^{1/2}$ and ${\hat{L}}_o^{1/2}$ as well as their product

$\begin{array}{ll}{\hat{L}}_{co}^{1/2}:={\hat{L}}_c^{1/2}{\hat{L}}_o^{1/2}=U_c{\Sigma }_c^{1/2}{U}_c^{\top }{U}_o{\Sigma }_o^{1/2}{U}_o^{\top }\in {ℝ}^{n_y\times n_y}.& (\mathrm{7b})\end{array}$

Finally, we derive the SVD of ${\hat{L}}_{co}^{1/2}$:

$\begin{array}{ll}{\hat{L}}_{co}^{1/2}=U\Sigma V^{\top }& (\mathrm{7c})\end{array}$

with orthogonal matrices U, $V\in {ℝ}^{n_y\times n_y}$ and a diagonal matrix $\Sigma \in {ℝ}^{n_y\times n_y}$ containing the non-negative singular values σ_i, i = 1, …, n_y, of ${\hat{L}}_{co}^{1/2}$ in descending order.

Remark 9. Recalling that for LTI systems, the matrix L_co = L_cL_o possesses an eigenvalue decomposition (cf. [18, p. 77]):

$\begin{array}{l}{L}_c{L}_o=Ũ\left[\begin{array}{cc}{\tilde{\Sigma }}^2& 0\\ 0& 0\end{array}\right]{Ũ}^{-1}\end{array}$

which is equivalent to the fact that $L_{co}^{1/2}$ has the following SVD

$\begin{array}{l}{L}_c^{1/2}{L}_o^{1/2}=[{Û}_1|{Û}_2]\left[\begin{array}{cc}\tilde{\Sigma }& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{\widehat{V}}_1^{\top }\\ {\widehat{V}}_2^{\top }\end{array}\right],\end{array}$

see [21] (Theorem 1). Then, the LTI system can be transformed into a minimal realization, i.e., the transformed LTI system is balanced, and its controllability and observability gramians are both equal to Σ.      ◇

Next, we turn to the reduced order modeling (ROM). Here, we utilize the approach from balanced truncation. To build a reduced-order model, we truncate the sufficiently small singular values, taking only ℓ ≪ n_y so that σ_i < ε holds for i = ℓ + 1, …, n_y. For that purpose, we partition the SVD as follows:

$\begin{array}{l}{\hat{L}}_c^{1/2}{\hat{L}}_o^{1/2}=U\Sigma V^{\top }=[U_1|U_2]\left[\begin{array}{cc}{\Sigma }_1& 0\\ 0& {\Sigma }_2\end{array}\right][V_1|V_2{]}^{\top }=\\ \left(\begin{array}{cc}{U}_1{\Sigma }_1{V}_1^{\top }& 0\\ 0& U_2{\Sigma }_2{V}_2^{\top }\end{array}\right),\end{array}$

where $U_1\in {ℝ}^{n_y\times \ell }$, $U_2\in {ℝ}^{n_y\times (n_y-\ell )}$, and ${\Sigma }_1\in {ℝ}^{\ell \times \ell }$, ${\Sigma }_2\in {ℝ}^{(n_y-\ell )\times (n_y-\ell )}$. We suppose that the matrix $W^{\ell }=V_1^{\top }{U}_1\in {ℝ}^{\ell \times \ell }$ is invertible. Now the reduced-order model for (1) is derived in a standard way (cf., e.g., [9]). For given large terminal time T > 0 and t ∈ [0, T], we approximate $y(t)\in {ℝ}^{n_y}$ by $U_1{y}^{\ell }(t)$, where y^ℓ(t) ∈ ℝ^ℓ solves together with $z^{\ell }(t)\in {ℝ}^{n_z}$

$\begin{array}{l}{W}^{\ell }{ẏ}^{\ell }(t)=A^{\ell }(t)y^{\ell }(t)+B^{\ell }(t)u(t)+\\ V_1^{\top }f(t)\text{ for }t\in (0,T],y^{\ell }(0)=V_1^{\top }{y}_{\circ },& (\mathrm{8a})\end{array}$

$\begin{array}{ll}{z}^{\ell }(t)=C(t)U_1{y}^{\ell }(t)\text{ for }t\in \left[0,T\right].& (\mathrm{8b})\end{array}$

with $A^{\ell }(t)=V_1^{\top }A(t)U_1\in {ℝ}^{\ell \times \ell }$, $B^{\ell }(t)=V_1^{\top }B(t)\in {ℝ}^{\ell \times n_u}$, $f^{\ell }(t)=V_1^{\top }f(t)$ and $y_{\circ }^{\ell }=V_1^{\top }{y}_{\circ }$.

Remark 10. 1) Equation (8a) is an ℓ-dimensional nonlinear system of differential equation. Since we assume ℓ ≪ n_y, (8) is called a low-dimensional or reduced order model for (1).

2) Note that for an LTI system this approach is equivalent to balanced truncation; see, e.g., [20, 21] for more details.

3) In our numerical examples the matrices are given as $W^{\ell }=V_1^{\top }W{U}_1$, $A^{\ell }(t)=V_1^{\top }WA(t)U_1\in {ℝ}^{\ell \times \ell }$, $B^{\ell }(t)=V_1^{\top }WB(t)\in {ℝ}^{\ell \times n_u}$, $f^{\ell }(t)=V_1^{\top }Wf(t)$ and $y_{\circ }^{\ell }=V_1^{\top }W{y}_{\circ }$; cf. Remarks 6 and 8.     ◇

3 Model reduction for linear time-variant optimality systems

Notices that (1) is an LTV system and therefore to reduce it one cannot apply directly balanced truncation. Here we follow a different approach. We are interested in controlling (1) in an optimal way. Thus, we define the following objective function

$\begin{array}{ll}\begin{array}{l}\mathcal{J}(y,u)=\frac{1}{2}{\displaystyle {\int }_0^T}{\Vert C(t)y(t)-z_{\text{d}}(t)\Vert }_Q^2\text{d}t+\frac{1}{2}{\displaystyle {\int }_0^T}{\Vert u(t)\Vert }_R^2\text{d}t\\ =\frac{1}{2}{\displaystyle {\int }_0^T}{\Vert z(t)-z_{\text{d}}(t)\Vert }_Q^2\text{d}t+\frac{1}{2}{\displaystyle {\int }_0^T}{\Vert u(t)\Vert }_R^2\text{d}t,\end{array}& (9)\end{array}$

where $z_{\text{d}}\in L^2(0,T;{ℝ}^{n_z})$ is a desired state, $Q\in {ℝ}^{n_z\times n_z}$ and $R\in {ℝ}^{n_u\times n_u}$ are symmetric positive definite matrices, $\Vert ·|{|}_Q={\langle ·,·\rangle }_Q^{1/2}$ and $\Vert ·|{|}_R={\langle ·,·\rangle }_R^{1/2}$ hold. Now the optimal control problems reads

$\begin{array}{l}min\mathcal{J}(y,u)\text{subject to (s.t.) }(y,u)\in H^1(0,T;{ℝ}^{n_y})\times (\text{P})\\ L^2(0,T;{ℝ}^{n_u})\text{ satisfies }(1).\end{array}$

The Lagrange functional associated with P is given by

$\begin{array}{r}ℒ(y,u,p,p_{°})=\mathcal{J}(y,u)+{\displaystyle {\int }_0^T\langle \dot{y}(}t)-A(t)y(t)-B(t)u(t)-\\ f(t),p(t){\rangle }_W\text{ d}t+{\langle y(0)-y_{\circ },p_{\circ }\rangle }_W,\end{array}$

where the weighted inner product has been defined in (6). Due to convexity, a first-order sufficient optimality system can be derived from stationarity conditions of the Lagrangian (see, e.g., [34])

$\begin{array}{ll}ẏ(t)=A(t)y(t)+B(t)R^{-1}B{(t)}^{\top }Wp(t)+f(t),t\in (0,T],y(0)=y_{\circ },& (\mathrm{10a})\end{array}$

$\begin{array}{l}-ṗ(t)=W^{-1}(A{(t)}^{\top }Wp(t)+\\ C{(t)}^{\top }Q(z_{\text{d}}(t)-C(t)y(t))),t\in [0,T),p(T)=0& (\mathrm{10b})\end{array}$

and

$\begin{array}{ll}u(t)=R^{-1}B{(t)}^{\top }Wp(t)\text{ for }t\in (0,T).& (11)\end{array}$

Note that utilizing (11) the control u is replaced in (10a). System (10) can be rewritten to get the form of a dynamical system, although formally it is not. We introduce the transformed variable q(t): = p(T − t). For brevity, we set

$\begin{array}{l}{\tilde{z}}_{\text{d}}(t)=z_{\text{d}}(T-t)\in {ℝ}^{n_z},Ã(t)=A(T-t)\in {ℝ}^{n_y\times n_y},\\ \hfill \tilde{C}(t)=C(T-t)\in {ℝ}^{n_z\times n_y}\end{array}$

for t ∈ [0, T]. This allows us to express (10) as

$\begin{array}{ll}\begin{array}{c}ẏ(t)=A(t)y(t)+B(t)R^{-1}B{(t)}^{\top }Wq(T-t)+f(t),t\in (0,T],y(0)=y_{\circ },\\ \stackrel{\circ }{q}(t)=W^{-1}(Ã{(t)}^{\top }Wq(t)+\tilde{C}{(t)}^{\top }Q({\tilde{z}}_{\text{d}}(t)-\hfill \\ \tilde{C}(t)y(T-t))),t\in [0,T),q(0)=0.\hfill \end{array}& (12)\end{array}$

For t ∈ (0, T) let

$\begin{array}{l}x(t)=\left[\begin{array}{c}y(t)\\ q(t)\end{array}\right]\in {ℝ}^{2n_y},\mathcal{F}(t)=\left[\begin{array}{c}f(t)\\ \tilde{C}{(t)}^{\top }Q{\tilde{z}}_{\text{d}}(t)\end{array}\right]\in {ℝ}^{2n_y},\\ x_{\circ }=\left[\begin{array}{c}{y}_{\circ }\\ 0\end{array}\right]\in {ℝ}^{2n_y}.\end{array}$

Then, (12) can be written as

$\begin{array}{ll}ẋ(t)=\mathcal{A}(t)x(t)+\tilde{\mathcal{A}}(t)x(T-t)+\mathcal{F}(t)\text{ for }t\in (0,T),x(0)=x_{\circ }& (13)\end{array}$

with the two (2n_y) × (2n_y)-matrices

$\begin{array}{l}\mathcal{A}(t)=\left[\begin{array}{cc}A(t)& 0\\ 0& W^{-1}Ã{(t)}^{\top }W\end{array}\right],\tilde{\mathcal{A}}(t)=\\ \left[\begin{array}{cc}0& B(t)R^{-1}B{(t)}^{\top }W\\ -W^{-1}\tilde{C}{(t)}^{\top }Q\tilde{C}(t)& 0\end{array}\right].\end{array}$

Although (13) has not a canonical form, it can be seen as a dynamical system, because the solution x evolves from 0 to the final time T. Thus, we can apply a model-order reduction scheme based on empirical gramians as introduced in Section 2. For that purpose we interprete (13) as an input-output system, where we consider the desired states z_d as input and the control $ũ(t)=R^{-1}B{(t)}^{\top }W\widehat{x}(T-t)$ as output with $\widehat{x}(t)={(x_{n_y+1}(t),\dots ,x_{2n_y}(t))}^{\top }$. Note that $\widehat{x}(t)$ is nothing else than the adjoint variable (modulus a time transformation) for P. Hence the output ũ is the optimal control sought. Applying the empirical gramian approach explained in Section 2, we can construct a reduced-order model for (13) which should achieve satisfactory approximation performances while the target z_d varies. What we aim in fact is not to construct a reduced order model for (1), but a way to cheaply compute the optimal solution of problem P for different values of its underlying parameters, such as the right-hand side f or the target z_d. This technique is particularly useful in the context of multiobjective optimization or model predictive control (MPC), where many optimization problems with changing parameters must be solved; cf. Section 4.

Remark 11. Let us define the right-hand side

$\begin{array}{l}F(t,x)=\left(\begin{array}{cc}A(t)& B(t)R^{-1}B{(t)}^{\top }W\\ W^{-1}C{(t)}^{\top }QC(t)& -W^{-1}A{(t)}^{\top }W\end{array}\right)x+\\ \left(\begin{array}{c}f(t)\\ -W^{-1}C{(t)}^{\top }Q{z}_{\text{d}}(t)\end{array}\right)\end{array}$

for t ≥ 0 and $x=(y,p)\in {ℝ}^{n_y}\times {ℝ}^{n_y}$. Then we can write the differential equations (10) as follows

$\begin{array}{ll}ẋ(t)=F(t,x(t))\text{ for }t>0.& (14)\end{array}$

However, we cannot pose an initial condition for the second component p of x (due to the terminal condition for the dual variable). Suppose that Φ_t(x_°) denotes the solution to (14) at time t ≥ 0 for an arbitrarily chosen initial condition $x_{\circ }\in {ℝ}^{2n_y}$. Then, Φ_t is a symplectic function [35, Lemma 1], and we can introduce symplectic Koopman operators to get accurate approximations of certain optimal nonlinear feedback laws; cf. [35] for more details. Here we follow a different approach.      ◇

4 Two applications

In this section we illustrate our proposed strategy for two different examples, which are highlighting the power of the empirical gramians in comparison to a standard approach based on Proper Orthogonal Decomposition (POD). In both examples we consider optimization problems with parabolic PDE constraints. In Section 4.1 we introduce the POD method, which will serve as comparison for our results. In Section 4.2 we apply our proposed methodology to a multiobjective optimization framework, while in Section 4.3 we test it in a feedback control context.

From now on, let T > 0 and Ω ⊂ ℝ^d (d ∈ {1, 2, 3}) with boundary Γ = ∂Ω. We set Q = (0, T) × Ω and Σ = (0, T) × ∂Ω. Furthermore, we define the function spaces V = H¹(Ω), H = L²(Ω), $\mathcal{Y}$ = W(0, T) = L²(0, T; V) ∩ H¹(0, T; V′) and $\mathcal{H}$ = $\mathcal{U}$ = L²(0, T; H). For more details on Lebesgue and Sobolev spaces we refer the reader to [36], for instance.

4.1 The POD method

The (discrete) POD method is based on constructing a low-dimensional subspace that can resemble the information carried out by a given set of vectors ${\{y_j^k\}}_{j=1}^n$, 1 ≤ k ≤ ℘, (the so-called snapshots) belonging to the Euclidean space ℝ^m; cf., e.g., [37] and [11] (Section 2.1). Let

$\begin{array}{l}\mathcal{V}=\text{span}\{y_j^k|1\le j\le n\text{ and }1\le k\le \wp \}\subset {ℝ}^m\end{array}$

be the space spanned by the snapshots with dimension n_{$\mathcal{V}$} = dim$\mathcal{V}$ ≤ min{n℘, m}. To avoid trivial cases we assume n_{$\mathcal{V}$} ≥ 0. For ℓ ≤ n_{$\mathcal{V}$} the POD method generates pairwise orthonormal functions ${\left\{{\psi }_i\right\}}_{i=1}^{\ell }$ such that all $y_j^k$, 1 ≤ j ≤ n and 1 ≤ k ≤ ℘, can be represented with sufficient accuracy by a linear combination of the ψ_i's. This is done by a minimization of the mean square error between the $y_j^k$'s and their corresponding ℓ-th partial Fourier sum:

$\begin{array}{ll}\{\begin{array}{l}min{\displaystyle \sum _{k=1}^{\wp }}{\displaystyle \sum _{j=1}^n}{\alpha }_j\Vert {\displaystyle y_j^k}-{\displaystyle \sum _{i=1}^{\ell }}{\langle {\displaystyle y_j^k},{\psi }_i\rangle }_{\tilde{W}}{\psi }_i{\displaystyle {\Vert }_{\tilde{W}}^2}\\ \text{s.t.}{\left\{{\psi }_i\right\}}_{i=1}^{\ell }\subset {ℝ}^m\text{ and }{\langle {\psi }_i,{\psi }_j\rangle }_{\tilde{W}}={\delta }_{ij},1\le i,j\le \ell ,\end{array}& (15)\end{array}$

where the α_j's are positive weighting parameters for j = 1, …, n. In fact, the weights α_j are chosen to resemble a trapezoidal rule for the temporal integration. Moreover, $\tilde{W}\in {ℝ}^{m\times m}$ defines the weighted inner product ${\langle ·,·\rangle }_{\tilde{W}}$ and associated norm $\Vert ·|{|}_{\tilde{W}}$; cf. (6). The symbol δ_ij denotes the Kronecker symbol satisfying δ_ii = 1 and δ_ij = 0 for i ≠ j. An optimal solution to (15) is denoted as a POD basis of rank ℓ. It can be proven that such a solution is characterized by the eigenvalue problem

$\begin{array}{ll}\mathcal{R}{\psi }_i={\lambda }_i{\psi }_i\text{ for }1\le i\le \ell ,& (16)\end{array}$

where λ₁ ≥ … ≥ λ_ℓ ≥ … ≥ λ_{n$\mathcal{V}$} > 0 denote the eigenvalues of the linear, compact, nonnegative and self-adjoint operator

$\mathcal{R}:{ℝ}^m\to {ℝ}^m,\mathcal{R}\psi ={\displaystyle \sum _{k=1}^{\wp }}{\displaystyle \sum _{j=1}^n}{\alpha }_j{\langle \psi ,y_j^k\rangle }_{\tilde{W}}{y}_j^k\text{ for }\psi \in {ℝ}^m;$

cf., e.g., [11] (Lemma 2.2). We refer to [11] (Remark 2.11) for an explicit form of the operator $\mathcal{R}$ in the case ℘ = 2. Recall that for a solution ${\left\{{\psi }_i\right\}}_{i=1}^{\ell }$ to (16) the following approximation error formula holds true:

${\displaystyle \sum _{k=1}^{\wp }}{\displaystyle \sum _{j=1}^n}{\alpha }_j\Vert y_j^k-{\displaystyle \sum _{i=1}^{\ell }}{\langle y_j^k,{\psi }_i\rangle }_{\tilde{W}}{\psi }_i{\Vert }_{\tilde{W}}^2={\displaystyle \sum _{i=\ell +1}^{n_{\mathcal{V}}}}{\lambda }_i;$

cf. [11] (Theorem 2.7). In our numerical tests the (temporal) snapshots ${\{y_j^k\}}_{j=1}^n$ will generally come from simulations of (13) for ℘ different choices of the inputs. The index n is related to the discrete time steps 0 = t₁ < … < t_n = T and ℘ to the number of inputs used. We define the POD matrix $\text{Ψ}:=\left[{\psi }_1|\dots |{\psi }_{\ell }\right]\in {ℝ}^{m\times \ell }$ and derive – analogous to (8a) – a POD-based reduced order model for (13) by choosing U₁ = Ψ and $V_1=\tilde{W}\text{Ψ}$ which satisfy $U_1^{\top }{V}_1={\text{Ψ}}^{\top }\tilde{W}\text{Ψ}=I_{\ell }$.

4.2 Multiobjective optimization problem

There is hardly ever a situation, where only one goal is of interest at a time. When carrying out a purchase for example, we want to pay a low price while getting a high quality product. In the same manner, multiple objectives are present in most technical applications such as fast and energy efficient driving or designing light and stable constructions. This dilemma leads to the field of multiobjective optimization (cf., e.g., [12]), where the aim is to minimize all relevant objectives J = (J₁, …, J_k) simultaneously.

While we are usually satisfied with one (global or even local) optimal solution in the single-objective setting, there generally exists an infinite number of optimal compromises in the situation where multiple objectives are present since the different objectives are in conflict with each other. Let us recall the concept of Pareto optimality:

1) A point ū dominates a point u, if J(ū) ≤ J(u) componentwise and J_i(ū) < J_i(u) for at least one i ∈ {1, …, k};

2) A feasible point ū is called Pareto optimal if there exists no feasible u dominating ū;

3) The set of Pareto optimal points is called the Pareto set, its image under J the Pareto front.

In multiobjective optimization, the goal can be therefore to compute the Pareto set, i.e., the set of points ū for which no u exists that is superior in all objective components, cf. Figure 1.

FIGURE 1

Figure 1. Sketch (red lines) of a Pareto set (A) and Pareto front (B) for min{J(u)|u ∈ ℝ²} with J:ℝ² → ℝ², $J_1(u)={(u_1+1)}^2+{(u_2+1)}^2$ and $J_2(u)={(u_1-1)}^2+{(u_2-1)}^2$. The goal is to simultaneously minimize J₁ and J₂ represented by the blue and the orange paraboloid, respectively.

4.2.1 Problem formulation

Now we turn to our PDE-constrained multiobjective optimization problem (MOP). Given the canonical embedding $\mathcal{C}:\mathcal{Y}\to \mathcal{H}$ and two desired states ${𝔷}_{\text{d}}^1$, ${𝔷}_{\text{d}}^2\in \mathcal{H}$ we consider the following strictly convex MOP:

$\begin{array}{l}\begin{array}{l}minJ(𝔲)=\left[\begin{array}{c}{J}_1(𝔲)\\ J_2(𝔲)\end{array}\right]=\frac{1}{2}\left[\begin{array}{l}{\displaystyle {\Vert \mathcal{C}{𝔶}_{𝔲}-{𝔷}_{\text{d}}^1\Vert }_{\mathcal{H}}^2}+\sigma {\Vert 𝔲\Vert }_{\mathcal{U}}^2\\ {\displaystyle {\Vert \mathcal{C}{𝔶}_{𝔲}-{𝔷}_{\text{d}}^2\Vert }_{\mathcal{H}}^2}+\sigma {\Vert 𝔲\Vert }_{\mathcal{U}}^2\hfill \end{array}\right]\\ \text{s.t.}𝔲\in \mathcal{U}\text{ and }𝔶={𝔶}_{𝔲}\in \mathcal{Y}\text{ solves the linear parabolic problem}\hfill \\ \begin{array}{l}\frac{\partial {𝔶}_{𝔲}}{\partial t}(t,\text{x})-\Delta {𝔶}_{𝔲}(t,\text{x})+𝔳(\text{x})·\nabla {𝔶}_{𝔲}(t,\text{x})=𝔲(t,\text{x})+\\ 𝔣(t,\text{x})\text{ for }(t,\text{x})\in Q,\\ \frac{\partial {𝔶}_{𝔲}}{\partial n}(t,\text{s})=0\text{ for }(t,\text{s})\in \Sigma ,{𝔶}_{𝔲}(0,\text{x})={𝔶}_{\circ }(\text{x})\text{ for }\text{x}\in \Omega ,\end{array}\end{array}\\ \hfill (\widehat{P})\end{array}$

where σ > 0, 𝔳 ∈ C(Ω; ℝ^d), 𝔣 ∈ $\mathcal{H}$ and 𝔶_° ∈ H hold. Recall that a unique (weak) solution 𝔶_𝔲 ∈ $\mathcal{Y}$ follows from standard results for any control 𝔲 ∈ $\mathcal{U}$; cf., e.g., [38, 39].

Due to convexity it can be proven that a sufficient first-order optimality condition for $\widehat{P}$ is given as follows: If there exist a weight $\bar{\alpha }\in \left[0,1\right]$ satisfying together with an element $\overline{𝔲}\in \mathcal{U}$ the conditions

$\begin{array}{ll}\bar{\alpha }\nabla J_1(\overline{𝔲})+(1-\bar{\alpha })\nabla J_2(\overline{𝔲})=0& (17)\end{array}$

then $\overline{𝔲}$ is a Pareto optimal point for $\widehat{P}$; see [12] (Theorem 3.27) and [13] (Section 2.1), for instance.

Remark 12. For fixed weight $\bar{\alpha }\in \left[0,1\right]$ we observe that (17) is a sufficient first-order optimality condition for the (strictly) convex problem

$\begin{array}{ll}min{J}^{\alpha }(𝔲)=\alpha J_1(𝔲)+(1-\alpha )J_2(𝔲)\text{s.t. }𝔲\in \mathcal{U}.& (18)\end{array}$

with the choice $\alpha =\bar{\alpha }$. This motivates the well-known weighted sum method; cf. [40]. By varying the weight α ∈ [0, 1] the strictly convex problem (18) is solved. This requires many optimization solves, and model order reduction offers the possibility to speed-up the solution process.      ◇

For arbitrarily chosen α ∈ [0, 1] we set

$\begin{array}{l}{\mathcal{C}}^{\alpha }=\left(\begin{array}{c}\alpha \mathcal{C}\\ (1-\alpha )\mathcal{C}\end{array}\right):\mathcal{Y}\to \mathcal{H}\times \mathcal{H}\text{and }{𝔷}_{\text{d}}^{\alpha }=\left(\begin{array}{c}\alpha {𝔷}_{\text{d}}^1\\ (1-\alpha ){𝔷}_{\text{d}}^2\end{array}\right)\\ \in \mathcal{H}\times \mathcal{H}.\end{array}$

Following Remark 12, we consider the (strictly) convex scalar optimization problem

$\begin{array}{l}\begin{array}{l}min{J}_{\alpha }(𝔲)=\frac{1}{2}{\Vert {\mathcal{C}}^{\alpha }{𝔶}_{𝔲}-{𝔷}_{\text{d}}^{\alpha }\Vert }_{\mathcal{H}\times \mathcal{H}}^2+\frac{\sigma }{2}{\Vert 𝔲\Vert }_{\mathcal{U}}^2\\ \text{s.t.}u\in \mathcal{U}\text{ and }{𝔶}_{𝔲}\in \mathcal{Y}\text{solves the linear parabolic problem}\\ \frac{\partial {𝔶}_{𝔲}}{\partial t}(t,\text{x})-\Delta {𝔶}_{𝔲}(t,\text{x})+\\ 𝔳(\text{x})·\nabla {𝔶}_{𝔲}(t,\text{x})=𝔲(t,\text{x})+𝔣(t,\text{x})\text{ for }(t,\text{x})\in Q,\\ \frac{\partial {𝔶}_{𝔲}}{\partial n}(t,\text{s})=0\text{ for }(t,\text{s})\in \Sigma ,{𝔶}_{𝔲}(0,\text{x})={𝔶}_{\circ }(\text{x})\text{ for }\text{x}\in \Omega .\end{array}\\ \hfill ({\widehat{P}}^{\alpha })\end{array}$

Next, we discretize ${\widehat{P}}^{\alpha }$ in space using piecewise linear finite elements (FE). Let m = n_y = n_u, n_z = 2m be valid and φ₁, …, φ_m ∈ V denote the FE ansatz functions. The FE subspace is $V^{\text{fe}}=\mathrm{\text{span}}\left\{{\phi }_1,\dots ,{\phi }_m\right\}\subset V$. For (t, x) ∈ Ω_T we approximate the state and control as

$\begin{array}{l}𝔶(t,\text{x})\approx {𝔶}^{\text{fe}}(t,\text{x})={\displaystyle \sum _{i=1}^m}{y}_i(t){\phi }_i(\text{x}),𝔲(t,\text{x})\approx {𝔲}^{\text{fe}}(t,\text{x})=\\ {\displaystyle \sum _{i=1}^m}{u}_i(t){\phi }_i(\text{x}),\end{array}$

respectively. Moreover, let

$\begin{array}{l}{𝔷}_{\text{d}}^j(t,\text{x})\approx {\displaystyle \sum _{i=1}^m}{z}_{\text{d}i}^j(t){\phi }_i(\text{x})(j=1,2),𝔣(t,\text{x})\\ \approx {\displaystyle \sum _{i=1}^m}{f}_i(t){\phi }_i(\text{x}),{𝔶}_{\circ }(\text{x})\approx {\displaystyle \sum _{i=1}^m}{y}_{\circ i}{\phi }_i(\text{x})\end{array}$

for (t, x) ∈ Q. Finally, we introduce the m × m matrices

$\begin{array}{l}M=(({\langle {\phi }_j,{\phi }_i\rangle }_H)),D=(({\langle {\phi }_j,v·\nabla {\phi }_i\rangle }_H)),\\ S=(({\langle \nabla {\phi }_j,\nabla {\phi }_i\rangle }_{H^d})),A=S+D,B=I_m,R=\sigma M,W=\tilde{W}=M\\ C=\left(\begin{array}{c}\alpha I_m\\ (1-\alpha )I_m\end{array}\right),z_{\text{d}}^{\alpha }(t)=\left(\begin{array}{c}\alpha z_{\text{d}}^1(t)\\ (1-\alpha )z_{\text{d}}^2(t)\end{array}\right),Q^{\alpha }=\left(\begin{array}{cc}\alpha M& 0\\ 0& (1-\alpha )M\end{array}\right).\end{array}$

Now problem ${\widehat{P}}_{\text{fe}}^{\alpha }$ is approximated by the following FE optimization problem

$\begin{array}{l}\begin{array}{l}min{\displaystyle J_{\text{fe}}^{\alpha }}(u)={\displaystyle {\int }_0^T}\frac{1}{2}{\Vert C^{\alpha }y(t)-z_{\text{d}}^{\alpha }(t)\Vert }_{Q^{\alpha }}^2+\frac{1}{2}{\Vert u(t)\Vert }_R^2\text{d}t\\ \text{s.t.}u\in {\mathcal{U}}^{\text{fe}}\text{ and }y\in {\mathcal{Y}}^{\text{fe}}\text{ solves the initial value problem}\\ Mẏ(t)=Ay(t)+Mu(t)+Mf(t)\text{ for }t\in (0,T],y(0)=y_{\circ },\end{array}\\ \hfill ({\widehat{P}}_{\text{fe}}^{\alpha })\end{array}$

where $\mathcal{U}$^fe = L²(0, T; ℝ^m) and $\mathcal{Y}$^fe = H¹(0, T; ℝ^m). Notice that in this semidiscrete setting, our problem is equivalent to (1a) of Section 3. For (10) we derive

$\begin{array}{ll}ẏ(t)=M^{-1}Ay(t)+\frac{1}{\sigma }p(t)+f(t),t\in (0,T],y(0)=y_{\circ },& (\mathrm{19a})\end{array}$

$\begin{array}{l}-ṗ(t)=M^{-1}{A}^{\top }p(t)+\alpha z_d^1(t)+\\ (1-\alpha )z_d^2(t)-y(t),t\in [0,T),p(T)=0,& (\mathrm{19b})\end{array}$

where we have used u(t) = R⁻¹B(t)^⊤Wp(t) = p(t)/σ.

In particular, a solution to ${\widehat{P}}_{\text{fe}}^{\alpha }$ can be computed by solving the sufficient first-order optimality system (13). For the specific case of ${\widehat{P}}_{\text{fe}}^{\alpha }$ we have

$\begin{array}{ll}ẋ(t)=\mathcal{A}x(t)+\tilde{\mathcal{A}}x(T-t)+{\mathcal{F}}^{\alpha }(t)\text{ for }t\in (0,T),x(0)=x_{\circ }& (20)\end{array}$

with

$\begin{array}{l}\mathcal{A}=\left[\begin{array}{cc}{M}^{-1}A& 0\\ 0& M^{-1}{A}^{\top }\end{array}\right],\tilde{\mathcal{A}}=\left[\begin{array}{cc}0& \frac{1}{\sigma }{I}_m\\ -I_m& 0\end{array}\right],\\ {\mathcal{F}}^{\alpha }(t)=\left[\begin{array}{c}f(t)\\ {\tilde{z}}_{\text{d}}^{\alpha }(t)\end{array}\right],x_{\circ }=\left[\begin{array}{c}{y}_{\circ }\\ 0\end{array}\right]\end{array}$

and ${\tilde{z}}_{\text{d}}^{\alpha }(t)=\alpha z_{\text{d}}^1(T-t)+(1-\alpha )z_{\text{d}}^2(T-t)$.

In conclusion, to compute the Pareto front of the multiobjective optimization problem ${\widehat{P}}_{\text{fe}}^{\alpha }$ with a weighted sum approach we have to solve several optimality systems of the form (20) with changing parameter α. For the gramian based MOR approach we proceed as described in Algorithm 1.

ALGORITHM 1

Algorithm 1. Gramian-based MOP using the weighted sum approach.

Remark 13.(1) Steps 1–4 of Algorithm 1 we call training or offline phase. For the gramian orthogonal matrices set $\mathcal{T}$^2m, respectively the positive constant set $\mathcal{M}$ we choose n_{$\mathcal{T}$} = n_{$\mathcal{M}$} = 3, where $I_{2m}\in {\mathcal{T}}^{2n_y}$ and 1 ∈ $\mathcal{E}$. The other two orthogonal matrices and scalars were chosen randomly. Furthermore, we restrict the set $\mathcal{E}$^2m to ten random unit vectors to save computational effort. Altogether, we compute 3 × 3 × 10 = 90 trajectories to (20) to get the empirical controllability gramian. For ẑ_d we choose a random input instead of δ-impulse inputs; cf., e.g., [26].

(2) Note that our gramian-based approach ensures sufficiently good approximation qualities for varying desired state ${\tilde{z}}_{\text{d}}^{\alpha }$ for α ∈ [0, 1]. This is required in the weighted sum method.

(3) For solving the dynamical systems (20) and (21) over a fixed time horizon [0, T] we use the implicit Euler scheme with the equidistant time grid t_j = (j − 1)Δt and Δt = T/(n − 1) for j = 1, …, n with step size Δt > 0. One can show that the resulting linear system is uniquely solvable.

(4) Since the empirical gramians are trained with respect to the desired state, the resulting basis functions are flexible while using different targets. Therefore, they are also more robust for a varying parameter α, which is pretty important while computing the Pareto front.     ◇

In our numerical experiments we compare the results that we got from Algorithm 1 with the ones obtained from a POD reduced-order modeling approach. The snapshots used to generate the POD basis are the solutions of (20) while computing the controllabilty and observability gramians in Algorithm 1. Clearly, this is done only for fair comparison purposes among the techniques. In general, building a POD model capable of approximating the solution of problems with different targets would require to run many simulations varying targets and controls. Unfortunately this is possible only if the targets are known a-priori. Otherwise the only option remaining is choosing random quantities, which would require exactly the same effort of the gramian approach. As we can see, the big difference is that the gramian approach is capable to exploit the randomness to construct a good approximation of the full order model, while the POD seriously struggles if the generated snapshots are not close to the optimal solution of the problem. Furthermore, notice that if we train our POD model by computing snapshots for varying targets and controls we would get a huge amount of snapshots, which will lead to a costly SVD, while the gramian have the advantage to lead to perform SVDs on matrix which have the same size of the problem. The costly SVD for the POD model is directly noticeable at the CPU time of the offline phase; cf. Table 2.

4.2.2 Numerical results

In the following numerical test the chosen parameter and problem required functions are summarized in Table 1.

TABLE 1

Table 1. Domain, optimization and PDE data for the numerical implementation of $\widehat{P}$.

First of all, already in Figure 2 one can notice the better accuracy of the proposed approach based on empirical gramians for two different number of reduced order models. To better quantify the accuracy of the reduced-order solutions we compare the absolute errors

$\begin{array}{l}t\mapsto \text{er}{\text{r}}_{\text{abs}}^J=\Vert (\begin{array}{c}{J}_1\left(u_{\alpha }^{\text{fe}}\right)\\ J_2\left(u_{\alpha }^{\text{fe}}\right)\end{array})-(\begin{array}{c}{J}_1\left(u_{\alpha }^{\ell }\right)\\ J_2\left(u_{\alpha }^{\ell }\right)\end{array}){\Vert }_2,\\ t\mapsto \text{er}{\text{r}}_{\text{abs}}^y={\Vert y_{\alpha }^{\text{fe}}\left(t\right)-y_{\alpha }^{\ell }\left(t\right)\Vert }_M,t\mapsto \text{er}{\text{r}}_{\text{abs}}^u={\Vert u_{\alpha }^{\text{fe}}\left(t\right)-u_{\alpha }^{\ell }\left(t\right)\Vert }_M,\end{array}$

and also the relative errors

$\begin{array}{l}t\mapsto \text{er}{\text{r}}_{\text{rel}}^J=\frac{\text{er}{\text{r}}_{\text{abs}}^J}{\Vert (\begin{array}{c}{J}_1(u_{\alpha }^{\text{fe}})\\ J_2(u_{\alpha }^{\text{fe}})\end{array}){\Vert }_2},t\mapsto \text{er}{\text{r}}_{\text{rel}}^y=\frac{\text{er}{\text{r}}_{\text{abs}}^y}{\Vert y_{\alpha }^{\text{fe}}(t){\Vert }_M},\\ t\mapsto \text{er}{\text{r}}_{\text{rel}}^u=\frac{\text{er}{\text{r}}_{\text{abs}}^u}{\Vert u_{\alpha }^{\text{fe}}(t)|{|}_M}.\end{array}$

FIGURE 2

Figure 2. Pareto front for the full-order model, for the empirical gramian-based model and the classical POD model with ℓ = 30 (left plot) and ℓ = 50 (right plot) basis functions.

These errors are reported in Figure 3 for the gramian and POD approaches. As one can see, the empirical gramians perform better than the POD method. Since the desired states ${𝔷}_{\text{d}}^1$ and ${𝔷}_{\text{d}}^2$ are not part of the training set, it is way harder for the POD-based reduced order model to reconstruct the optimal solution. Also switching the desired states ${𝔷}_{\text{d}}^1$ and ${𝔷}_{\text{d}}^2$ afterwards one get similar approximation results. This shows the advantage of flexibility of the gramian-based approach. With significant reduction of time (due to cheaper SVDs), one can prepare a basis which is more suitable for varying targets and thus different optimal control problem. This is a large advantage if one considers that often the POD basis must be tailored for each problem at a time. In Table 2, one can see the significant speed-up of the offline phase with the empirical gramians method with respect to the POD one. We picked here the cases ℓ = 30 and ℓ = 50. As said, for the gramian approach we can add all snapshots up in two matrices, whose sizes coincide with the dimension of the problem. For the POD basis, instead, we need a way larger matrix, which depends on the number of snapshots. If one generates many snapshots, as in this case, the SVD of such a matrix gets numerical costly. The reason is that for multiple snapshots the matrix that has to be decomposed by SVD gets larger; cf. [11] (Remark 2.12).

FIGURE 3

Figure 3. Errors for the reduced order solutions with empirical gramians (solid line) compared to POD approximimation (dashed line). The parameter ℓ represents the number of basis functions.

TABLE 2

Table 2. CPU times to solve the multiobjective optimization problem $\widehat{P}$ and corresponding speed-ups for the respective MOR techniques for ℓ basis functions.

In the brackets there are the respective values including the training time. Notice that the offline time for the gramian-based MOR is significantly lower than for the POD method which leeds to better speed-up factors. The reason for this is that the required SVD for the POD basis computation has to be computed for a much larger data matrix including all snapshots.

4.3 Model predictive control

In the second example we show the potential of the empirical gramians in a model predictive control (MPC) framework. MPC is an optimization based method for the feedback control of dynamical systems where the time horizon usually tends to infinity; see, e.g., [14] for a general introduction and [41] for parabolic PDEs. For large terminal time T we are considering problems of the form

$\begin{array}{l}\begin{array}{l}{\displaystyle \underset{(y,u)}{min}}\mathcal{J}(y,u)=\frac{1}{2}{\displaystyle {\int }_0^T}{\Vert C(t)y(t)-z_{\text{d}}(t)\Vert }_Q^2\text{d}t+\frac{1}{2}{\displaystyle {\int }_0^T}{\Vert u(t)\Vert }_R^2\text{d}t\\ =\frac{1}{2}{\displaystyle {\int }_0^T}{\Vert z(t)-z_{\text{d}}(t)\Vert }_Q^2\text{d}t+\frac{1}{2}{\displaystyle {\int }_0^T}{\Vert u(t)\Vert }_R^2\text{d}t,\\ \text{s.t.}(y,u)\in H^1(0,T;{ℝ}^{n_y})\times L^2(0,T;{ℝ}^{n_u})\text{ solves }\\ ẏ(t)=A(t)y(t)+B(t)u(t)+f(t)\text{ for }t\in (0,T],y(0)=y_{\circ }.\end{array}\\ \hfill (\text{MPC})\end{array}$

MPC is based on iterative, finite-horizon optimization of dynamical control problems. Let first fix a time discretization with initial time step t₀ and t_k = t₀ + kΔt for Δt > 0 and k ∈ ℕ. Fixed a time shift T_pred = NΔt for N > 0, N ∈ ℕ, we compute an optimal control u_k such that u_k(t_i) is minimizing the cost function for t_i in the discrete time interval [t_k, t_k + T_pred], i.e., i = k, …, k + N. Now, only the optimal control over the discrete time interval [t_k, t_k + T_f] (T_f = MΔt with 0 < M ≤ N) is stored and the previous procedure is repeated starting from the new resulting current state at time t_k + T_f. This yields to a new optimal control u_k+1 which will contribute to another piece of the final reconstructed suboptimal control for the whole infinite horizon. In conclusion we are solving iteratively optimal control problems, where the finite prediction horizon keeps being shifted forward. The parameter T_pred is called MPC prediction horizon and T_f is the MPC feedback time step. The reason why we are building a feedback control relies on the fact that each new computed optimal control for a given time window depends on the current state of the system. We present our MPC method in Algorithm 2.

ALGORITHM 2

Algorithm 2. Model predictive control (MPC).

4.3.1 Problem formulation

In our numerical example we study the following problem:

$\begin{array}{l}\begin{array}{l}minJ(u)=\frac{1}{2}{\displaystyle {\Vert \mathcal{C}{𝔶}_u-{𝔷}_d\Vert }_{\mathcal{H}}^2}+\frac{\sigma }{2}{\Vert u\Vert }_{\mathcal{U}}^2\\ \text{s.t.}𝔲\in \mathcal{U}\text{ and }{𝔶}_{𝔲}\in \mathcal{Y}\text{ solves the linear parabolic problem }\\ \frac{\partial {𝔶}_{𝔲}}{\partial t}(t,\text{x})-\Delta {𝔶}_{𝔲}(t,\text{x})+\alpha (t)𝔳(\text{x})·\nabla {𝔶}_{𝔲}(t,\text{x})={\chi }_{\omega }(\text{x})𝔲(t,\text{x})+\\ 𝔣(t,\text{x})\text{ for }(t,\text{x})\in Q,\\ \frac{\partial {𝔶}_{𝔲}}{\partial n}(t,\text{s})=0\text{ for }(t,\text{s})\in \Sigma ,{𝔶}_{𝔲}(0,\text{x})={𝔶}_{\circ }(\text{x})\text{ for }\text{x}\in \Omega ,\end{array}\\ \hfill (\tilde{P})\end{array}$

where the spaces $\mathcal{Y}$, $\mathcal{U}$ and $\mathcal{H}$, the linear embedding operator $\mathcal{C}$ and the data σ, 𝔶_°, 𝔳 and 𝔣 have been introduced at the beginning of Sections 4, 4.2.1. Moreover, let α ∈ L^∞(0, T) be scalar-valued parameter function andχ_ω the characteristic function of a nonzero subset ω ⊂ Ω. Recall that a unique (weak) solution 𝔶_𝔲 ∈ $\mathcal{Y}$ for any control 𝔲 ∈ $\mathcal{U}$ follows from standard results; see, e.g., [38, 39].

Next, we proceed analog to Section 4.2 to derive the semidiscrete problem with weighted norms. The semidiscrete first-order necessary and sufficient optimality system of $\tilde{P}$ is equivalent to solving the system

$\begin{array}{ll}ẋ(t)={\mathcal{A}}^{\alpha }(t)x(t)+\tilde{\mathcal{A}}x(T-t)+\mathcal{F}(t)\text{ for }t\in (0,T),x(0)=x_{\circ }\in {ℝ}^{2m},& (22)\end{array}$

for x = (y, q), q(t) = p(T − t), with y and p solution of state and adjoint equation of $\tilde{P}$ respectively, and

$\begin{array}{l}{\mathcal{A}}^{\alpha }(t)=\left[\begin{array}{cc}{M}^{-1}{A}^{\alpha }(t)& 0\\ 0& M^{-1}{Ã}^{\alpha }{(t)}^{\top }\end{array}\right],{\tilde{\mathcal{A}}}^{\alpha }(t)=\left[\begin{array}{cc}0& \frac{1}{\sigma }{M}^{-1}B{M}^{-1}{B}^{\top }\\ -I_m& 0\end{array}\right],\\ B=(({\langle {\chi }_{\omega }{\phi }_j,{\phi }_i\rangle }_H)),\mathcal{F}(t)=\left[\begin{array}{c}f(t)\\ {\tilde{z}}_d(t)\end{array}\right],x_{\circ }=\left[\begin{array}{c}{y}_{\circ }\\ 0\end{array}\right],\\ A^{\alpha }(t)=S+\alpha (t)D,{Ã}^{\alpha }(t)=S+\alpha (T-t)D.\end{array}$

Moreover, we have u(t) = M⁻¹B^⊤p(t)/σ.

Now, the possible approach for applying model-order reduction to $\tilde{P}$ would be to apply the POD method; cf., e.g., [11]. One would run simulations of state and adjoint equations varying the control u in order to construct a POD basis. The more the inputs u for the training are closer to the unknown optimal control, the better the POD approximation will be. Clearly this shows the limitation of this approach, since it is not really possible to guess a-priorly what the optimal control will be. Anyway, since we are in an MPC framework, one could think of paying the price of solving occasionally FE models to generate (and update) the POD basis and perform the rest of the runs with the resulting reduced-order model. This strategy has also the advantage that can be controlled by the use of an a-posteriori error estimator to decide when trigger the update. For this MPC-POD method we refer to [15, 42]. With respect to this method, the gramian-based approach offers additional advantages:

1) The possibility of training the controllability gramian (and thus the reduced-order model) choosing as input the advection coefficient α(t). In such a way, the resulting basis functions will be less sensitive to perturbations of this advection coefficient.

2) Similarly, we can train the observability gramian by choosing as output z(t) = M⁻¹B^⊤p(t)/σ, which nothing else that the first-order optimality condition of the perturbed optimal control problem that has (22) as optimality system. In such a way, the reduced-order model will be capable to reconstruct the optimal control with high accuracy.

As also pointed out in [15], in fact, the MPC-POD method suffers in presence of advection already for a small Peclet number and any strong deviation from the original POD snapshots can also lead to inaccuracy. Although the a-posteriori error estimator mitigates these problems, the continue triggering of the update of the basis has also a negative impact in the overall time performances of the algorithm. The gramian based MPC can overcome this by preparing a reduced-order model slightly sensible to perturbations of advection and optimal control. The basis functions are capable to keep track of the changes in this case, thanks to the randomized training performed at the begin. For a more detailed description of the gramian based MPC we refer to Algorithm 3, for the MPC-POD algorithm we refer to [15] (Algorithm 11).

ALGORITHM 3

Algorithm 3. Gramian-based MPC.

4.3.2 Numerical results

For the following numerical tests we choose the same domain and PDE parameter than in Section 4.1.1 apart from the end time point, which is here set to T = 250. The additional advection coefficient function α ∈ L^∞(0, T; ℝ) is chosen randomly in a specified range, which will be varied to defined different test cases. We report the other chosen parameters in Table 3.

TABLE 3

Table 3. Optimization and MPC data for the numerical implementation of $\tilde{P}$.

The numerical results are reported in Tables 4, 5. The quantity $\text{er}{\text{r}}_{\text{abs}}^y$, $\text{er}{\text{r}}_{\text{abs}}^u$, $\text{er}{\text{r}}_{\text{rel}}^y$ and $\text{er}{\text{r}}_{\text{rel}}^u$ are the absolute errors and the relative errors in approximating full-order optimal state ȳ and control ū, respectively. More precisely, being ȳ^ℓ and ū^ℓ a reduced-order model solution, we define

$\begin{array}{l}er{r}_{abs}^y={\left({\int }_0^T\Vert \overline{y}(t)-{\overline{y}}^{\ell }(t){\Vert }_M^2\text{ d}t\right)}^{1/2},er{r}_{rel}^y=\frac{er{r}_{abs}^y}{{({\int }_0^T\Vert \overline{y}(t){\Vert }_M^2\text{ d}t)}^{1/2}},\\ er{r}_{abs}^u={\left({\int }_0^T\Vert \overline{u}(t)-{\overline{u}}^{\ell }(t){\Vert }_M^2\text{ d}t\right)}^{1/2},er{r}_{rel}^u=\frac{er{r}_{abs}^u}{{({\int }_0^T\Vert \overline{u}(t){\Vert }_M^2\text{ d}t)}^{1/2}},\end{array}$

TABLE 4

Table 4. Error in approximating the optimal state ȳ and control ū solution of the discretized full order model $\tilde{P}$ with randomly chosen advection coefficient function α in the reported range and ℓ ∈ {10, 30, 50} basis functions.

TABLE 5

Table 5. Average CPU times to solve the MPC framework for $\tilde{P}$ and corresponding speed-ups for the respective MOR techniques after 1,000 simulations with randomly chosen advection coefficients function α.

where the time integrals are realized numerical by applying a standard trapezoidal approximation. As said, since we train the gramians with respect to the advection coefficient α, we obtain basis functions which are robust with respect to changes of the advection term. This can be seen in Table 4, where the relative error increases significantly only for the largest range of [0, 10].

The MPC-POD method, instead, starts loosing accuracy already at the stage [0, 5]. Furthermore, the relative error for the control resulting from Algorithm 3 is from two to four order of magnitude smaller with respect to the MPC-POD one. In particular, in the last case of α(t) ∈ [0, 10] the MPC-POD method is not practically usable. Therefore, we can say that Algorithm 3 is able to recover the optimal solution with a better approximation than the MPC-POD method. Let us mention, that for a fair comparison, we used the same snapshots for the two approaches. The problem for POD is that the resulting basis has no capability to track the information carried by perturbing the advection, while the gramian can. Anyway, the deterioration of accuracy for increasing advection emerges also in the context of the gramian basis, although does not affect its usability in our numerical tests. From Table 5, one can also see another disadvantage of the POD method.

Since we are using the same number of snapshots for the two techniques, the POD matrix has the dimension of N_x × L, where N_x is the of FE nodes and L the number of snapshots. Clearly, if L ≫ N_x, the time to compute a POD basis, i.e., performing a SVD on the matrix N_x × L is significantly larger than the two SVDs required for the N_x × N_x gramians. This reflects the difference on the overall time speed-up between MPC-POD (1.7) and Algorithm 3 (4.8). We can then conclude that our proposed gramian based approach has the double advantage of robustness with respect to variation of the advection field and overall time speed-up in comparison to the MPC-POD method.

4.3.3 A-Posteriori error analysis

To validate our numerical approach we control the MOR error by an a-posteriori error estimate. The idea is that the error of the difference between the (unknown) optimal control and its suboptimal approximation computed by the MOR-based discretization can be estimated without knowing the optimal solution (𝔶_α, 𝔲_α) to $\tilde{P}$. The techniques are based on well-known Galerkin-type estimates. In the Appendix we have summarized the results. For more details we refer to [8, 10, 11], for instance.

From Figure 4 it can be seen that in each step of the MPC algorithm the error in the computed suboptimal controls are smaller than 0.1. The impact of the truncation value ℓ on the error is thoroughly examined in [43].

FIGURE 4

Figure 4. Reduced MPC with one reduced gramain (left) and POD (right) basis with reduced basis rank ℓ = 50 and range(α) = [0.1]. One every column are the absolute L²-spatial-difference in state and control over the time, i.e., $t\mapsto \Vert ȳ(t)-{ȳ}^{\ell }(t)|{|}_M$ and $t\mapsto \Vert ū(t)-{ū}^{\ell }(t)|{|}_M$ (top) and the particular error estimates for every MPC step.

5 Conclusions

In the present paper empirical gramian have been used to derive reduced order models tailored for the first-order optimality system which consists of the coupled state and dual equation. Here, data of the optimization problem serve as inputs for the computation of the gramian matrices. The validation by a-posteriori error analysis shows that we get reliable reduced order models which turn out to perform better than reduced order models based on a standard POD, which utilizes exactly the same snapshots for the ROM. The reason of such performances is the possibility of including controllability and observability knowledge through the empirical gramian approach. This leads to the fact that the gramian-based ROM is more appropriate for the optimization purposes. Finally, as future perspective, the proposed approach can be applied to nonlinear dynamical systems as well. Here, the accuracy could be checked by hierarchical a-posteriori error analysis; see, e.g., [44], which extends the work in [45]. To get an efficient numerical reduced order realization for the nonlinear dynamical system, one would need to apply the (discrete) empirical interpolation method (cf. [22–24]) or missing point evaluation (cf. [25]).

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

The authors acknowledge funding by the Independent Research Grant Efficient model order reduction for model predictive control of non-linear input-output dynamical systems of the Zukunftskolleg at the University of Konstanz and by the Deutsche Forschungsgemeinschaft (DFG) for the project Localized Reduced Basis Methods for PDE-Constrained Parameter Optimization under contract VO 1658/6-1.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fams.2023.1144142/full#supplementary-material

References

1. Becker R, Kapp H, Rannacher R. Adaptive finite element methods for optimal control of partial differential equations: basic concept. SIAM J Control Optimiz. (2000) 39:113–32. doi: 10.1137/S0363012999351097

CrossRef Full Text | Google Scholar

2. Liu W, Yan N. A posteriori error estimates for distributed convex optimal control problems. Adv Comp Math. (2001) 15:285–309. doi: 10.1023/A:1014239012739

CrossRef Full Text | Google Scholar

3. Holmes P, Lumley J, Berkooz G, Rowley C. Turbulence, Coherent Structures, Dynamical Systems and Symmetry. 2nd ed. Cambridge: Cambridge University Press (2013).

Google Scholar

4. Hesthaven J, Rozza G, Stamm B. Certified Reduced Basis Methods for Parametrized Partial Differential Equations. SpringerBriefs in Mathematics. Cham: Springer (2016).

Google Scholar

5. Quarteroni A, Manzoni A, Negri F. Reduced Basis Methods for Partial Differential Equations: An Introduction. UNITEXT-La Matematica per il 3+2. Cham: Springer (2016).

Google Scholar

6. Benner P, Ohlberger M, Cohen A, Willcox K. Model Reduction and Approximation. Philadelphia, PA: Society for Industrial and Applied Mathematics (2017).

Google Scholar

7. Benner P, Mehrmann V, Sorensen D. Dimension reduction of large-scale systems: Proceedings of a Workshop held in Oberwolfach, Germany, October 19–25, 2003. In: Lecture Notes in Computational Science and Engineering. Heidelberg: Springer (2006).

Google Scholar

8. Hinze M, Volkwein S. Error estimates for abstract linear-quadratic optimal control problems using proper orthogonal decomposition. Comput Optim Appl. (2008) 39:319–45. doi: 10.1007/s10589-007-9058-4

CrossRef Full Text | Google Scholar

9. Antoulas AC. Approximation of Large-Scale Dynamical Systems. Advances in Design and Control. Philadelphia, PA: SIAM (2009).

Google Scholar

10. Tröltzsch F, Volkwein S. Pod a-posteriori error estimates for linear-quadratic optimal control problems. Comput Optim Appl. (2009) 44:83–115. doi: 10.1007/s10589-008-9224-3

CrossRef Full Text | Google Scholar

11. Gubisch M, Volkwein S. Proper Orthogonal Decomposition for Linear-Quadratic Optimal Control, Chapter 1. Philadelphia, PA: SIAM (2017). p. 3–63.

Google Scholar

12. Ehrgott M. Multicriteria Optimization. 2nd ed. Berlin: Springer (2005).

Google Scholar

13. Banholzer S. ROM-Based Multiobjective Optimization with PDE Constraints (Ph.D. thesis). University of Konstanz (2021). Available online at: http://nbnresolving.de/urn:nbn:de:bsz:352-2-1g98y1ic7inp29 (accessed November 28, 2023).

Google Scholar

14. Grüune L, Pannek J. Nonlinear Model Predictive Control: Theory and Algorithms. 2nd ed. London: Springer (2016).

Google Scholar

15. Mechelli L. POD-based State-Constrained Economic Model Predictive Control of Convection-Diffusion Phenomena (Ph.D. thesis), University of Konstanz (2019). Available online at: http://nbn-resolving.de/urn:nbn:de:bsz:352-2-2zoi8n9sxknm1 (accessed November 28, 2023).

Google Scholar

16. Bellman R. Introduction to the Mathematical Theory of Control Processes. Vol. II: Nonlinear Processes. New York, NY: Academic Press (1971).

Google Scholar

17. Isidori A. Nonlinear Control Systems: An Introduction. 2nd ed. Berlin; Heidelberg: Springer-Verlag (1989).

Google Scholar

18. Zhou K, Doyle J, Glover K. Robust and Optimal Control. Upper Saddle River, NJ: Prentice-Hall (1996).

Google Scholar

19. Hernández-Lerma O, Laura-Guarachi L, Mendoza-Palacios S, González-Sánchez D. An Introduction to Optimal Control Theory. The Dynamic Programming Approach. Texts in Applied Mathematics. Cham: Springer (2023).

Google Scholar

20. Lall S, Marsden JE, Glavaški S. Empirical model reduction of controlled nonlinear systems. In: IFAC Proceedings Volumes, 32:2598–2603. 14th IFAC World Congress 1999, Beijing, Chia, 5–9 July. Oxford: Pergamon (1999).

Google Scholar

21. Garcia J, Basilio J. Computation of reduced-order models of multivariable systems by balanced truncation. Int J Syst Sci. (2002) 33:847–54. doi: 10.1080/0020772021000017308

CrossRef Full Text | Google Scholar

22. Barrault M, Maday Y, Nguyen N, Patera A. An empirical interpolation method: application to efficient reduced-basis discretization of partial differential equations. Comp Rendus Math. (2004) 339:667–72. doi: 10.1016/j.crma.2004.08.006

CrossRef Full Text | Google Scholar

23. Chaturantabut S, Sorensen D. Nonlinear model reduction via discrete empirical interpolation. SIAM J Sci Comp. (2010) 32:2737–64. doi: 10.1137/090766498

CrossRef Full Text | Google Scholar

24. Chaturantabut S, Sorensen D. A state space estimate for POD-DEIM nonlinear model reduction. SIAM J Numer Anal. (2012) 50:46–63. doi: 10.1137/110822724

CrossRef Full Text | Google Scholar

25. Astrid P, Weiland S, Willcox K, Backx T. Missing point estimation in models described by proper orthogonal decomposition. IEEE Trans Automat Contr. (2008) 53:2237–51. doi: 10.1109/TAC.2008.2006102

CrossRef Full Text | Google Scholar

26. Himpe C. emgr-The empirical gramian framework. Algorithms. (2018) 11:91. doi: 10.3390/a11070091

CrossRef Full Text | Google Scholar

27. Kalman R. Mathematical description of linear dynamical systems. SIAM J. Control Optim. (1963) 1:182–92. doi: 10.1137/0301010

CrossRef Full Text | Google Scholar

28. Kailath T. Linear Systems. Englewood Cliffs, NJ: Prentice-Hall (1980).

Google Scholar

29. Baur U, Benner P. Cross-gramian based model reduction for datasparse systems. Electron Transact Numer Anal. (2008) 31:256–70. doi: 10.1137/070711578

CrossRef Full Text | Google Scholar

30. Moore B. Principal component analysis in linear systems: controllability, observability, and model reduction. IEEE Trans Automat Contr. (1981) 26:17–32. doi: 10.1109/TAC.1981.1102568

CrossRef Full Text | Google Scholar

31. Aldhaheri R. Model order reduction via real schur-formdecomposition. Int J Control. (1991) 53:709–16. doi: 10.1080/00207179108953642

CrossRef Full Text | Google Scholar

32. Gallivan K, Vandendorpe A, Van Dooren P. Sylvester equations and projection-based model reduction. J Comput Appl Math. (2004) 162:213–29. doi: 10.1016/j.cam.2003.08.026

CrossRef Full Text | Google Scholar

33. Sorensen D, Antoulas A. The sylvester equation and approximate balanced reduction. Linear Algebra Appl. (2002) 351/2:671–700. doi: 10.1016/S0024-3795(02)00283-5

CrossRef Full Text | Google Scholar

34. Betts JT. Practical methods for optimal control and estimation using nonlinear programming. In: Advances in Design and Control. 2nd, ed. Philadelphia, PA: SIAM (2010).

Google Scholar

35. Villanueva M, Jones C, Houska B. Towards global optimal control via koopman lifts. Automatica. (2021) 132:109610. doi: 10.1016/j.automatica.2021.109610

CrossRef Full Text | Google Scholar

36. Evans L. Partial Differential Equations. Graduate Studies in Mathematics. American Mathematical Society. Providence, RI: American Mathematical Society (2010).

Google Scholar

37. Kunisch K, Volkwein S. Galerkin proper orthogonal decomposition methods for parabolic problems. Numer Math. (2001) 90:117–48. doi: 10.1007/s002110100282

CrossRef Full Text | Google Scholar

38. Dautray R, Lions J-L. Mathematical Analysis and Numerical Methods for Science and Technology. Vol. 5: Evolution Problems I. Berlin: Springer (2000).

Google Scholar

39. Pazy A. Semigoups of Linear Operators and Applications to Partial Differential Equation. Berlin: Springer (1983).

Google Scholar

40. Zadeh L. Optimality and non-scalar-valued performance criteria. IEEE Trans Automat Contr. (1963) 8:59–60. doi: 10.1109/TAC.1963.1105511

CrossRef Full Text | Google Scholar

41. Azmi B, Kunisch K. A hybrid finite-dimensional rhc for stabilization of time-varying parabolic equations. SIAM J Control Optimiz. (2019) 57:3496–526. doi: 10.1137/19M1239787

CrossRef Full Text | Google Scholar

42. Mechelli L, Volkwein S. POD-based economic model predictive control for heat-convection phenomena. In: Radu FA, Kumar K, Berre I, Nordbotten JM, Pop IS, , editors. Numerical Mathematics and Advanced Applications ENUMATH 2017. Cham: Springer International Publishing (2019). p. 663–671.

Google Scholar

43. Rohleff J. Data-Driven Model-Order Reduction for Model Predictive Control (Master thesis). University of Konstanz (2023).

Google Scholar

44. Azmi B, Petrocchi A, Volkwein S. Adaptive parameter optimization for an elliptic-parabolic system using the reduced-basis method with hierarchical a-posteriori error analysis. arXiv. (2023)

Google Scholar

45. Hain S, Ohlberger M, Radic M, Urban K. A hierarchical a posteriori error estimator for the reduced basis method. Adv Comp Math. (2019) 45:2191–214. doi: 10.1007/s10444-019-09675-z

CrossRef Full Text | Google Scholar