Keywords: PDE-constrained optimal control, preconditioning, feedback stabilization, conjugate gradient methods

1 Introduction↩︎

Optimal control of differential equations is a well-established field of mathematical analysis, numerical methods and scientific computing. Problems subject to ordinary or partial differential equations (ODEs or PDEs) have been thoroughly analyzed and numerically approached [1]–[3]. The fundamental structure in these systems is that pairs of optimization variables \((x,u)\) satisfying the PDE may be characterized via a control-to-state map. This structure is useful analytically but may also be used to derive efficient numerical algorithms as the control-to-state map immediately provides access to a projection onto the feasible set. Many PDE-constrained optimal control problems may be viewed as a realization of the abstract linear-quadratic optimal control problem (OCP) \[\label{eq:OCP} \min_{(x,u)\in \mathcal{X}\times \mathcal{U}} \frac{1}{2}\|\mathcal{C}x-y_\mathrm{ref}\|_\mathcal{Y}^2 + \frac{\alpha}{2} \|u\|_\mathcal{U}^2 \quad \mathrm{s.t.} \quad \mathcal{A}x-\mathcal{B}u = f,\tag{1}\] where

1. \(\mathcal{X}\) is a Banach space and \(\mathcal{Y},\mathcal{U},\mathcal{P}\) are Hilbert spaces,

2. \(\mathcal{A}\in L(\mathcal{X},\mathcal{P}^*)\) is continuously invertible and encodes some ODE or PDE,

3. \(\mathcal{B}\in L(\mathcal{U},\mathcal{P}^*)\) is an input operator and \(\mathcal{C}\in L(\mathcal{X},\mathcal{Y})\) is an observation operator,

4. \(y_\mathrm{ref} \in \mathcal{Y}\) is a reference output, \(f \in \mathcal{P}^*\) is a forcing term, and \(\alpha > 0\) is a regularization parameter.

Here, the constraint could correspond either to a static problem that is, an elliptic PDE for which \(\mathcal{A}\) is usually a second-order (elliptic) operator or a time-dependent problem such as a parabolic or hyperbolic partial differential equation, where the constraint operator \(\mathcal{A}\) involves space and time derivatives. Such linear-quadratic problems occur, e.g., in in each step of an SQP-method applied to nonlinear problems.

One of the simplest approaches to solve 1 is to eliminate the PDE-constraint using invertibility of \(\mathcal{A}\), that is, leveraging the control-to-state map \(u \mapsto \mathcal{A}^{-1}\mathcal{B}u\). Hence, setting \(x=\mathcal{A}^{-1}(\mathcal{B}u + f)\) one may eliminate the state from the optimal control problem and consider the reduced or condensed problem. Then, standard approaches such as gradient methods may be applied [1], [2]. Alternatively, one may keep the constraint and introduce a Lagrange multiplier and derive suitable first-order optimality conditions in function space. Then, these could be approached by linear or nonlinear solvers, leading to an algorithm in state, control and adjoint variables such as conjugate gradient methods endowed with suitably designed constraint preconditioners [4]–[6]. In both of these approaches, however, the condition number of the systems to be solved strongly depends on the norm of the solution operator \(\|\mathcal{A}^{-1}\|_{\mathcal{P}^*\to \mathcal{X}}\), where \(\| \cdot \|_{\mathcal{P}^* \to \mathcal{X}}\) denotes the standard operator norm in \(L(\mathcal{P}^*,\mathcal{X})\). This operator norm, however, can be prohibitively large, e.g., when considering unstable time-dependent problems on large time horizons or elliptic problems with small ellipticity constant. While usually this is treated as given by the problem, we show in this work how to shape the constraint operator \(\mathcal{A}\) by an equivalence transformation such that the transformed problem is subject to an operator \(\widetilde{\mathcal{A}} = \mathcal{A}+ \mathcal{B}\mathcal{K}\), where \(\mathcal{K}\in L(\mathcal{X},\mathcal{U})\) is a feedback operator leading to the property \(\|{\widetilde{\mathcal{A}}}^{-1}\|_{\mathcal{P}^*\to \mathcal{X}} \ll \|{\mathcal{A}}^{-1}\|_{\mathcal{P}^*\to \mathcal{X}}\). Feedback transformations and feedback stabilization are classical and well-understood tools in mathematical systems theory [7]. While they are typically studied from a control-theoretic perspective, we show that they can also be leveraged in optimal control to drastically improve the conditioning and performance of numerical solvers.

The central tools are feedback stabilizations that are borrowed from mathematical systems theory where they are usually used to stabilize dynamic control systems; here, we showcase how they may also be leveraged for reformulations in optimal control in view of efficiency and robustness of numerical methods.

To illustrate the guiding idea of this work, we introduce a state feedback operator \(\mathcal{K}:\mathcal{X}\to \mathcal{U}\) and substitute \[\label{eq:feedback} u=\mathcal{K}x+v,\tag{2}\] where \(x\in \mathcal{X}\) is the state and \(v\in \mathcal{U}\) is a new control variable. In this new control variable, the OCP 1 reads \[\label{eq:OCP2} \min_{(x,v)\in \mathcal{X}\times \mathcal{U}} \frac{1}{2}\|\mathcal{C}x-y_\mathrm{ref}\|_\mathcal{Y}^2 + \frac{\alpha}{2} \|\mathcal{K}x+v\|_\mathcal{U}^2 \quad \mathrm{s.t.} \quad (\mathcal{A}-\mathcal{B}\mathcal{K})x - \mathcal{B}v = f.\tag{3}\] This problem is clearly equivalent to the problem 1 as the transformation 2 can easily be inverted via \(v = u - \mathcal{K}x\), hence it defines a bijection (independently of \(x\in \mathcal{X}\)). However, as illustrated in Figure 1, a suitably chosen feedback operator \(\mathcal{K}\) leads to strongly improved convergence behavior of iterative numerical optimal control solvers applied to 3 instead of 1 .

In the following we abbreviate the norm of the stabilized control-to-observation map as follows: \[\label{eq:norm-c-to-obs} \sigma_\mathcal{K}:= \|\mathcal{C}(\mathcal{A}-\mathcal{B}\mathcal{K})^{-1}\mathcal{B}\|_{\mathcal{U}\to \mathcal{Y}}\tag{4}\] so that the unfeedbacked cases is recovered by \(\sigma_0=\|\mathcal{C}\mathcal{A}^{-1}\mathcal{B}\|_{\mathcal{U}\to \mathcal{Y}}\) which relates to the original control-to-observation map of 1 .

Code availability. The code for all numerical experiments presented in this work is available in the repository

2 Iterative methods for PDE-constrained optimal control↩︎

In this part, we introduce two well-established methods in numerical optimal control that particularly benefit from our suggested feedback stabilization, namely a gradient method for the reduced problem as introduced in Subsection 2.1 and a preconditioned conjugate gradient method for the optimality system as recalled in Subsection 2.2.

In the context of iterative solvers for the optimality system corresponding to 3 it is crucial to design preconditioners that are tailored to the underlying problem formulation. The framework of operator preconditioning for such systems (cf. [6], [8], [9]) requires bounding the energy of the operator inner product by the preconditioner inner product.

While the most descriptive measure for convergence of Krylov methods is the distribution of the eigenvalues [10] it is often difficult to obtain robust bounds for the eigenvalues in context of optimal control, see [11] for the heat equation with full observation. A measure that is also very informative for the convergence behavior is the condition number [12] giving rise to relative energy error estimates of the form \[\frac{\|e_k\|_{A}}{\|e_0 \|_{A}} \le 2 \left(\frac{\sqrt{\kappa(A)} - 1}{\sqrt{\kappa(A)} + 1} \right)^k\] with \(e_k\) being the error at step \(k\in \mathbb{N}\) of the conjugate gradient method [13]. Hence, an estimate of the condition number \(\kappa\) can be central in the convergence analysis. Moreover, independence of the condition number w.r.t.problem parameters allows for robust convergence rates in the considered system class. The condition number of a symmetric linear operator \(A: \mathcal{U}\to\mathcal{U}^*\) with respect to a preconditioner \(Q:\mathcal{U}\to \mathcal{U}^*\) can be computed as follows. Assume that there are \(0 <m \le M\), such that \[m \langle v,Qv\rangle_{\mathcal{U},\mathcal{U}^*} \le \langle v,Av\rangle_{\mathcal{U},\mathcal{U}^*} \le M \langle v, Q v\rangle_{\mathcal{U},\mathcal{U}^*} \quad \forall v\in \mathcal{U}.\] If the constants in these estimates are sharp, then the condition number of \(A\) relative to \(Q\) is given as \[\kappa(A) := \frac{M}{m}.\]

2.1 Gradient method for the reduced problem↩︎

Due to continuous invertibility of the state operator associated with the constraint, we may eliminate the state via \(x = \mathcal{A}^{-1}\mathcal{B}u + \mathcal{A}^{-1}f\) such that we may equivalently solve the unconstrained problem \[\label{eq:redprob} \min_{u\in \mathcal{U}} h(u) := \frac{1}{2}\|\mathcal{C}(\mathcal{A}^{-1}\mathcal{B}u + \mathcal{A}^{-1}f) - y_\mathrm{ref}\|_\mathcal{Y}^2 + \frac{\alpha}{2} \|u\|_\mathcal{U}^2.\tag{5}\] Due to strict convexity of the quadratic function \(h:\mathcal{U}\to \mathbb{R}\), the solution of this minimization problem amounts to solving the linear system \[\label{eq:cg} 0 = h'(u^*) = \left((\mathcal{C}\mathcal{A}^{-1}\mathcal{B})^*\mathcal{R}_y\mathcal{C}\mathcal{A}^{-1}\mathcal{B}+ \alpha \mathcal{R}_u\right)u^* + (\mathcal{C}\mathcal{A}^{-1}\mathcal{B})^*\mathcal{R}_y(\mathcal{C}\mathcal{A}^{-1}f - y_\mathrm{ref}).\tag{6}\] Here, \(\mathcal{R}_y: \mathcal{Y}\to \mathcal{Y}^*\) and \(\mathcal{R}_u: \mathcal{U}\to \mathcal{U}^*\) denote the Riesz isomorphisms corresponding to the Hilbert input and output spaces, respectively. A gradient direction can be computed by solving the Riesz equation on \(\mathcal{U}\): \[\nabla h(u) = \mathcal{R}_u^{-1}h'(u).\] As the system 6 is governed by a self-adjoint operator \[\mathcal{M}_\mathrm{red} : \mathcal{U}\to \mathcal{U}^*, \qquad \mathcal{M}_\mathrm{red}:= (\mathcal{C}\mathcal{A}^{-1}\mathcal{B})^*\mathcal{R}_y\mathcal{C}\mathcal{A}^{-1}\mathcal{B}+ \alpha \mathcal{R}_u,\] a suitable method to solve it is the conjugate gradient method, where in each iteration the evaluation of the governing operator \(\mathcal{M}_\mathrm{red}\) involves one solution of the state equation (due to \(\mathcal{A}^{-1}\)) and one solution of the adjoint equation (due to \(\mathcal{A}^{-*}\)). To assess the convergence behavior of a conjugate gradient method applied to 6 , we estimate the relevant condition number.

For the self-adjoint operator defined above, we obtain the following well-known bound on the condition number, see e.g. [14]. Note that the proof follows also as a particular case of Proposition 4.

This result shows that the bound on the condition number strongly depends on the inverse operator \(\|\mathcal{A}^{-1}\|_{\mathcal{P}^* \to \mathcal{X}}\). In Section 4 we will show for various applications of elliptic, parabolic and hyperbolic equations that this operator norm may become prohibitively large and that the feedback transformation, that is, in essence replacing \(\mathcal{A}\) by \(\mathcal{A}-\mathcal{B}\mathcal{K}\) as in 3 may alleviate this effect.

2.2 Preconditioned projected conjugate gradients↩︎

The second approach used to solve the optimal control problem 1 builds upon the first-order optimality condition of the constrained problem without eliminating the state.

To this end, we recall the optimality conditions of the abstract optimal control problem 1 , see, e.g., [15]. Therein, the central assumption is that the constraint operator \([\mathcal{A},-\mathcal{B}]\) has closed range, which in our setting follows from the fact that \(\mathcal{A}\) is an isomorphism, hence surjective such that \([\mathcal{A},-\mathcal{B}]\) is surjective (and thus has closed range) as well. Let \((x,u)\in \mathcal{X}\times \mathcal{U}\) be optimal for 1 . Then, there is an adjoint state \(p\in \mathcal{P}^*\) such that \[\label{eq:optsys} \mathcal{M}_\mathrm{opt}\begin{bmatrix} x\\u\\p \end{bmatrix}= \begin{bmatrix} \mathcal{C}^*\mathcal{R}_y\mathcal{C}& 0 & \mathcal{A}^*\\ 0 & \alpha \mathcal{R}_u& -\mathcal{B}^*\\ \mathcal{A}& -\mathcal{B}&0 \end{bmatrix}\begin{bmatrix} x\\u\\p \end{bmatrix} = \begin{bmatrix} \mathcal{C}^*y_\mathrm{ref}\\ 0\\ f \end{bmatrix}.\tag{7}\] While this system is symmetric, it is indefinite due to the presence of the zero block in the bottom right. Hence, to apply a conjugate gradient method, one may incorporate a constraint preconditioner \[\label{eq:const95prec} \mathcal{Q}= \begin{bmatrix} 0 & 0& \mathcal{A}^*\\ 0& \alpha\mathcal{R}_u& -\mathcal{B}^*\\ \mathcal{A}& -\mathcal{B}& 0 \end{bmatrix}\tag{8}\] which is invertible due to invertibility of \(\mathcal{A}\) (cf. [4]–[6]). Due to its triangular structure, application of its inverse decouples into one adjoint solve, the solution of the Riesz equation \(\alpha \mathcal{R}_uw = v\) and one state solve.

When endowed with this preconditioner, it is guaranteed that the state and control iterates of the conjugate gradient method evolve in the subspace \(\ker [\mathcal{A},-\mathcal{B}]\times \mathcal{P}\) on which the governing operator \(\mathcal{M}_\mathrm{opt} : \mathcal{X}\times\mathcal{U}\times\mathcal{P}\to (\mathcal{X}\times\mathcal{U}\times\mathcal{P})^*\) as defined in 7 is positive definite.

When using a constraint preconditioner, the speed of convergence of PPCG is governed by the condition number of the top left \(2\times2\) block of \(\mathcal{M}_\mathrm{opt}\) with respect to the top left \(2\times 2\) block of \(\mathcal{Q}\) on \(\ker [\mathcal{A},-\mathcal{B}] = \{(\mathcal{A}^{-1}\mathcal{B}u, u)\,|\,u\in \mathcal{U}\}\).

One may use this insight to actually show the equivalence of the two condition number estimates of Proposition 1 and Proposition 2. In particular, we note that for \(z:= (x,u,p)\in \mathcal{Z}:= \mathcal{X}\times \mathcal{U}\times \mathcal{P}\) with \((x,u)\in \ker [\mathcal{A},-\mathcal{B}]\), \[\begin{align} \label{eq:itsthesame} \langle z,\mathcal{M}_\mathrm{opt}z\rangle_{\mathcal{Z},\mathcal{Z}^*} = \|\mathcal{C}x\|^2_\mathcal{Y}+ \alpha\|u\|^2_\mathcal{U}= \|\mathcal{C}\mathcal{A}^{-1}\mathcal{B}u\|^2_\mathcal{Y}+ \alpha\|u\|^2_\mathcal{U}= \langle u, \mathcal{M}_\mathrm{red}u\rangle_\mathcal{U}. \end{align}\tag{9}\] In this sense, the scalar product induced by the restriction of \(\mathcal{M}_\mathrm{opt}\) to \(\ker [\mathcal{A},-\mathcal{B}] \times \mathcal{P}\) (as done for the PPCG method in Proposition 2) coincides with the scalar product induced by \(\mathcal{M}_\mathrm{red}\). Hence, we immediately obtain the following result.

2.3 Illustrative scalar example↩︎

In this part, we briefly provide a very simple scalar problem to illustrate the main idea of the suggested approach. For more complex problems, we refer to Section 4 where we illustrate for a wide range of partial differential equations (such as elliptic, hyperbolic and parabolic problems) how to suitably choose the feedback operator to obtain a small solution operator norm in the transformed problem. For now, let us consider the discrete-time problem \[\begin{align} \label{eq:discopt} \begin{aligned} \frac{1}{2}\sum_{k=1}^{N} |x_k - 5|^2 + \frac{\alpha}{2} \sum_{k=0}^{N-1} |u_k|^2. \quad \mathrm{s.t.}\quad x_{k+1} &= ax_k + u_k, \;k=0,\ldots N-1, \\[-.3cm] x_0 &= 0. \end{aligned} \end{align}\tag{10}\] Here, \(\mathcal{U}= \mathcal{X}= \mathcal{P}=\mathcal{Y}= \mathbb{R}^N\) and \(0\leq a\in \mathbb{R}\) is a parameter to be specified. It is easily seen that the control-to-state map \(\mathcal{A}^{-1}\) corresponds to solving the forward iteration of the discrete time constraint, that is, \[\begin{align} \label{eq:gradient:calA} { \left[\begin{smallmatrix} x_{1}\\x_{2}\\\vdots\\x_{N} \end{smallmatrix}\right]} = { \left[\begin{smallmatrix} u_0\\ a x_0 + u_1\\ \vdots \\ a x_{N-1} + u_{N-1} \end{smallmatrix}\right]} = { \left[\begin{smallmatrix} 1 & 0 & 0 & \ldots & 0 \\ a & 1 & 0 & \ldots & 0 \\ a^2 & a & 1 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a^{N-1} & \ldots & \ldots & a & 1 \end{smallmatrix}\right]} { \left[\begin{smallmatrix} u_0\\ u_1\\\vdots \\ u_{N-1} \end{smallmatrix}\right]} =: T_N u \end{align}\tag{11}\] and input resp.output operators are just identity matrices \(\mathcal{C}= \mathcal{B}= I \in \mathbb{R}^{N\times N}\). Hence, one may easily show that (see Appendix 5.1)

Hence, if the original problem in the constraint is unstable (\(a>1\)), the operator norm deteriorates for larger horizons. This reflects the instability of the uncontrolled problem when \(a>1\). In this case, \(\mathcal{A}^{-1}\) and its adjoint significantly amplify any kind of errors, which relates to the condition number of the governing matrix derived in Proposition 1 or Proposition 2.

However, we may easily stabilize the constraint by introducing a feedback in the sense of 2 via \(u=kx+v\) with \(k=-a+0.5\) such that \(\widetilde{\mathcal{A}}\) is governed by the stable equation \[x_{k+1} = ax_k + u = ax_k + (-a+0.5)x_k +v = 0.5x_k + v.\] Consequently, the associated control-to-state map denoted by \(\widetilde{\mathcal{A}}^{-1}\) satisfies 11 with \(a= 0.5\) such that \(\|{\widetilde{\mathcal{A}}}^{-1}\| \sim \max\{1,|0.5|^{N-1}\} = 1\). In particular, the solution operator norm is bounded uniformly in the time horizon \(N\).

We briefly illustrate the effect of the feedback for the gradient method introduced in Subsection 2.1 in Figure 2. We depict the relative error of the control over the iterations for different parameters \((a,\alpha)\) of the problem 10 with horizon \(N=100\). We observe that the transformation leads to a better convergence even for the stable problem in the top row of Figure 2. This effect gets more drastic, when going to an unstable problem (bottom row) or when decreasing the Tikhonov parameter \(\alpha\). Importantly, as expected due to the feedback, the convergence behavior for the transformed problem is very robust w.r.t.the stability parameter \(a\) when comparing the red lines of the top and bottom row of Figure 2. This clearly is not the case for the original problem. Last, we note that, when increasing the instability parameter to \(a=1.5\) or increasing the time horizon, the gradient method for the original problem terminates with an overflow error, while the one for the feedback-transformed problem still converges.

3 Solution of the feedback-transformed problem↩︎

In the previous section, we have introduced two numerical algorithms and illustrated how they might benefit from feedback stabilizations. However, the algorithms suggested in Subsections 2.1 and 2.2 particularly leveraged the fact that there are no mixed terms involving both \(x\) and \(u\) in the cost functional. For the gradient method of the reduced problems, this implies that the governing operator \(\mathcal{M}_\mathrm{red}\) of 6 consists of only two summands of which only one involves state and adjoint solves. For the PPCG method, as a consequence of the absence of coupling terms in the cost functional, the top left \(2\times 2\) block of \(\mathcal{M}_{\mathrm{opt}}\) was block-diagonal.

In this part, we thus illustrate how to apply these methods to the feedback-transformed problem. In particular, we provide bounds on the condition number of the transformed problem. Before analyzing the algorithms, however, we briefly analyze the feedbacked optimal control problem. To this end we note that in view of optimality conditions (see 7 ), it is usually very helpful if the constraint operator has closed range. Moreover, we provide a relation of the optimal solutions of the original problem 1 and its feedback-transformed counterpart 3 .

We now state the main assumption on the choice of the feedback.

We briefly comment on this assumption. First, we note that Assumption 1 can always be satisfied by choosing \(\mathcal{K}=0\). In this case, however, the transformed problem is identical to the original problem. When choosing a feedback operator \(0\neq \mathcal{K}\in L(\mathcal{X},\mathcal{U})\), Assumption 1[ass:iso] is to ensure that the methods are applicable as both build upon solving the constraint for the state. Further, Assumption 1[ass:smallness] ensures smallness of the feedback w.r.t.the cost. If the observation operator \(\mathcal{C}: \mathcal{X}\to \mathcal{Y}\) has closed range, this may be characterized geometrically in terms of a kernel inclusion. Precisely, it is easily checked that if \(\operatorname{ran}(\mathcal{C})\) is closed, then Assumption 1[ass:smallness] is equivalent to \(\ker \mathcal{C}\subset \ker \mathcal{K}.\)

In the following Subsections 3.1 and 3.2 we show how to use the discussed numerical methods for the feedback-transformed problems and provide bounds on the condition number.

3.1 Gradient method for the reduced problem↩︎

In this subsection, we analyze the gradient method for the condensed problem when applied to the feedback-transformed optimal control problem 3 .

To this end, we denote \(\widetilde{\mathcal{A}} = \mathcal{A}-\mathcal{B}\mathcal{K}\) such that eliminating the state by the control in 3 via \(x = \widetilde{\mathcal{A}}^{-1} (\mathcal{B}v + f)\) we obtain the reduced transformed problem \[\min_{v\in \mathcal{U}} \tilde{h}(v) := \frac{1}{2}\|\mathcal{C}({\widetilde{\mathcal{A}}}^{-1}\mathcal{B}v + {\widetilde{\mathcal{A}}}^{-1}f) - y_\mathrm{ref}\|_\mathcal{Y}^2 + \frac{\alpha}{2} \|\mathcal{K}({\widetilde{\mathcal{A}}}^{-1}\mathcal{B}v + {\widetilde{\mathcal{A}}}^{-1}f)+v\|_\mathcal{U}^2.\] Expanding the squares, this cost functional is (up to constant terms) given by \[\begin{align} \tilde{h}(v) &= \frac{1}{2} \|\mathcal{C}{\widetilde{\mathcal{A}}}^{-1}\mathcal{B}v\|_\mathcal{Y}^2 + \langle \mathcal{C}{\widetilde{\mathcal{A}}}^{-1}\mathcal{B}v,\mathcal{C}{\widetilde{\mathcal{A}}}^{-1}f-y_\mathrm{ref} \rangle_{\mathcal{Y}} \\ & \qquad + \frac{\alpha}{2}\| \mathcal{K}{\widetilde{\mathcal{A}}}^{-1}\mathcal{B}v + v\|_\mathcal{U}^2 + \alpha\langle \mathcal{K}{\widetilde{\mathcal{A}}}^{-1}\mathcal{B}v + v,\mathcal{K}\widetilde{\mathcal{A}}^{-1}f\rangle \end{align}\] leading to the optimality condition \[\begin{align} \label{eq:cg2} \begin{aligned} 0 = \nabla \tilde{h}(v^*) &= \left((\mathcal{C}\widetilde{\mathcal{A}}^{-1}\mathcal{B})^*\mathcal{R}_y\mathcal{C}\widetilde{\mathcal{A}}^{-1}\mathcal{B}+ \alpha (\mathcal{K}\widetilde{\mathcal{A}}^{-1}\mathcal{B}+ I)^*\mathcal{R}_u(\mathcal{K}\widetilde{\mathcal{A}}^{-1}\mathcal{B}+ I)\right)v^* \\ & \qquad + (\mathcal{C}\mathcal{A}^{-1}\mathcal{B})^*\mathcal{R}_y(\mathcal{C}\mathcal{A}^{-1}f - y_\mathrm{ref}) + \alpha(\mathcal{K}\widetilde{\mathcal{A}}^{-1}\mathcal{B}+ I)^*\mathcal{R}_u\mathcal{K}\widetilde{\mathcal{A}}^{-1}f. \end{aligned} \end{align}\tag{13}\] Hence, the transformed governing operator is \(\widetilde{\mathcal{M}}_\mathrm{red} : \mathcal{U}\to \mathcal{U}^*\) defined by \[\widetilde{\mathcal{M}}_\mathrm{red} := (\mathcal{C}\widetilde{\mathcal{A}}^{-1}\mathcal{B})^*\mathcal{R}_y\mathcal{C}\widetilde{\mathcal{A}}^{-1}\mathcal{B}+ \alpha (\mathcal{K}\widetilde{\mathcal{A}}^{-1}\mathcal{B}+ I)^*\mathcal{R}_u(\mathcal{K}\widetilde{\mathcal{A}}^{-1}\mathcal{B}+ I).\]

Before discussing possible improvements of the condition number, we briefly discuss the complexity of applying \(\widetilde{\mathcal{M}}_\mathrm{red}\) in comparison to the original operator \(\mathcal{M}_\mathrm{red}\). We may organize the application of \(\widetilde{\mathcal{M}}_\mathrm{red}\) to \(v\in \mathcal{U}\) as follows: First we compute the state and the original control via \(y = \widetilde{\mathcal{A}}^{-1}\mathcal{B}v\) and \(u= \mathcal{R}_u(\mathcal{K}y+v)\). Then we may express the evaluation as \(\widetilde{\mathcal{M}}_\mathrm{red}v = \mathcal{B}^*\widetilde{\mathcal{A}}^{-*}(\mathcal{C}^*\mathcal{R}_y\mathcal{C}y+\alpha \mathcal{K}^* u) + \alpha u\). In this way, just as in the original problem governed by \(\mathcal{M}_{\mathrm{red}}\), one state and one adjoint solve are sufficient to compute a gradient step. The only additional cost, compared to the case \(\mathcal{K}=0\), is the application of the feedback operator and its adjoint.

We now prove a condition number estimate for the feedbacked reduced optimality system.

We briefly mention that while one may leverage the feedback to reduce a possible instability manifesting in the solution operator norm, the feedback bound \(\delta_\mathcal{K}\) also appears explicitly in the estimates (e.g. ?? ) such that too large feedbacks might worsen the condition number. Hence, it is in practice crucial to balance the feedback size and the solution operator norm. If the solution operator norm \(\sigma_\mathcal{K}\) may be linked with the relative bound on the feedback \(\delta_\mathcal{K}\) of Assumption 1[ass:smallness], then one may optimize the upper bound ?? w.r.t.the feedback, see Figure 3.

We will discuss these examples and how to obtain the optimal feedback strength \(\delta_\mathcal{K}^*\) that minimizes condition number bound ?? in detail in Section 4. However, we already see here that the theoretical optimization of the condition number leads to the best performance of the numerical method, that is, the smallest error after a fixed number of iterations.

Next, we analyze the PPCG method for the feedbacked problem.

3.2 Preconditioned projected conjugate gradients↩︎

To apply the PPCG method as introduced in Subsection 2.2 to the transformed problem 3 , we first deduce the optimality system for the transformed problem. To this end, observe by Proposition 3 that the transformed constraint operator has closed range such that we may again apply the abstract optimality conditions from [15].

Let \((x,v) \in \mathcal{X}\times \mathcal{U}\) be optimal for 3 . Then, there is \(p\in \mathcal{P}^*\) such that \[\label{eq:optsys2} \begin{bmatrix} \mathcal{C}^*\mathcal{R}_y\mathcal{C}+ \alpha \mathcal{K}^*\mathcal{R}_u\mathcal{K}& \alpha \mathcal{K}^*\mathcal{R}_u& (\mathcal{A}- \mathcal{B}\mathcal{K})^*\\ \alpha \mathcal{R}_u\mathcal{K}& \alpha \mathcal{R}_u& -\mathcal{B}^*\\ (\mathcal{A}- \mathcal{B}\mathcal{K})& -\mathcal{B}&0 \end{bmatrix}\begin{bmatrix} x\\v\\p \end{bmatrix} = \begin{bmatrix} \mathcal{C}^*y_\mathrm{ref}\\ 0\\ f \end{bmatrix}.\tag{15}\] We briefly comment on the relation between the optimality systems of 1 and 3 .

In the transformed case, a suitable constraint preconditioner for 15 is given by \[\widetilde{\mathcal{Q}}= \begin{bmatrix} 0 & 0 & (\mathcal{A}-\mathcal{B}\mathcal{K})^*\\ 0 & \alpha\mathcal{R}_u& -\mathcal{B}^*\\ \mathcal{A}-\mathcal{B}\mathcal{K}& -\mathcal{B}& 0 \end{bmatrix}.\]

We briefly comment on the estimates of Propositions 4 and 6. Therein, we observe a dependence of the condition number on the feedbacked solution operator \(\widetilde{\mathcal{A}}^{-1} = (\mathcal{A}- \mathcal{B}\mathcal{K})^{-1}\) instead of the original solution operator \(\mathcal{A}^{-1}\). Hence, if the feedback is chosen suitably to reduce the solution operator norm, the condition number may be drastically improved. This improvement may come from stabilization of a time-dependent problem as done for the illustrating scalar example of Subsection 2.3 and as will be done for general ODEs, heat and wave equations in the subsequent Section 4. For stationary problems, such as elliptic partial differential equations, we will show how to enlarge the coercivity constant of the associated bilinear form.

4 Applications↩︎

In this part, we provide several important applications of the above developed abstract approach. In Subsection 4.1 we discuss the case of ODE constraints to illustrate the main mechanism of the feedback stabilization in view of linear system theory. Then, in Subsection 4.2 we move to elliptic PDEs, where feedbacks are used to shape the coercivity of the governing main operator. Last, we move again to time-dependent problems and apply our method to parabolic equations in Subsection 4.3 and to hyperbolic equations in Subsection 4.4.

4.1 Finite-dimensional control systems↩︎

In this part, we consider linear finite-dimensional control systems \[\label{eq:fin95dim95cs} \dot{x}(t) = A x(t) + Bu(t) \quad t\in [0,T], \qquad x(0) = x_0,\tag{16}\] where \(T>0\) is a time horizon, \(A \in \mathbb{R}^{n \times n}\) denotes the system matrix, \(B \in \mathbb{R}^{m \times n}\) is a control input matrix, \(u: [0, T) \to \mathbb{R}^m\) is a control function and \(x: [0, T) \to \mathbb{R}^n\) is the associated state. For brevity, we refer to the control system 16 as the pair \((A,B)\). For every \(x_0 \in \mathbb{R}^n\) and \(u \in L^2(0,T;\mathbb{R}^m)\), 16 admits a unique mild solution \(x \in H^1(0,T;\mathbb{R}^n) \subset C([0,T];\mathbb{R}^n)\) given by the variation-of-constants formula \[\label{eq:var95of95const} x(t) = e^{tA}x_0 + \int_0^t e^{(t-s)A} B u(s)\,\mathrm{d}s, \quad t \in [0,T].\tag{17}\] Moreover, \(x\) satisfies 16 almost everywhere. Note that, due to [16], for any matrix \(A\in \mathbb{R}^{n\times n}\) there exist \(\omega \in \mathbb{R}\) and \(M \geq 1\) such that \[\label{eq:Gronwall} \left\|e^{t A} \right\|_{\mathbb{R}^n\to\mathbb{R}^n} \leq M e^{\omega t} , \qquad t \geq 0.\tag{18}\] We briefly comment on the constants involved in this estimate in view of spectral properties of the matrix \(A\).

To rewrite the control system 16 as an abstract constraint as in 1 , we introduce the operators \(\mathcal{A}: H^1(0,T; \mathbb{R}^n) \to L^2(0,T; \mathbb{R}^n) \times \mathbb{R}^n\) and \(\mathcal{B}: L^2(0,T; \mathbb{R}^m) \to L^2(0,T; \mathbb{R}^n) \times \mathbb{R}^n\) where \[\label{eq:defnA} \mathcal{A}x = { \left[\begin{smallmatrix} \tfrac{\mathrm{d}}{\mathrm{d}t}x-Ax \\ x(0) \end{smallmatrix}\right]} , \qquad \mathcal{B}u = { \left[\begin{smallmatrix} Bu \vphantom{\tfrac{\mathrm{d}}{\mathrm{d}t}x-Ax} \\ 0 \vphantom{x(0)} \end{smallmatrix}\right]} .\tag{19}\] Consequently, 16 is abstractly written as \[\label{eq:abstracode} \mathcal{A}x - \mathcal{B}u = { \left[\begin{smallmatrix} 0 \\ x_0 \end{smallmatrix}\right]} .\tag{20}\] To include state feedbacks into the abstract operator setting, we define the bounded multiplication operator \(\mathcal{K}: L^2(0,T; \mathbb{R}^n) \to L^2(0,T; \mathbb{R}^m)\) via \[\label{eq:multop} (\mathcal{K}x)(t) = Kx(t)\tag{21}\] for a given state feedback matrix \(K\in \mathbb{R}^{n \times m}\). Then, the closed-loop dynamics associated with the static feedback \(u(t) = Kx(t) + v(t)\) for all \(t\in [0,T]\) reads \[\label{eq:fin95dim95closed95loop} \dot{x}(t) = (A + BK)x(t) + Bv(t) , \qquad x(0) = x_0.\tag{22}\] with unique solution given by the variations-of-constants formula \[x(t) = e^{t(A+BK)}x_0 + \int_0^t e^{(t-s)(A+BK)} B v(s)\,\mathrm{d}s. \quad t \in [0,T].\] Hence, if one can suitably shape the spectrum of \(A\) with the feedback matrix \(K\) such that \(\omega_0(A+BK) < 0\), then 18 holds also with any \(\omega_0(A+BK) < \omega < 0\). This setting is termed exponential stabilizability in the realm of mathematical system theory [16], as the system 22 is subject to an exponentially stable matrix.

Correspondingly, the feedbacked version of 20 may be obtained by replacing \(u\) with \(v\) and \(\mathcal{A}\) with \[\label{eq:feedbackedcalA} \widetilde{\mathcal{A}}x = (\mathcal{A}- \mathcal{B}\mathcal{K}) x = { \left[\begin{smallmatrix} \tfrac{\mathrm{d}}{\mathrm{d}t}x-(A+BK) x \\ x(0) \end{smallmatrix}\right]}\tag{23}\]

The following result yields a bound on the operator norm in dependence of the stability parameter \(\omega\). The proof is presented in Appendix 5.2.

The following result is immediate from Proposition 8.

4.2 Elliptic PDEs↩︎

In this part, we study an application to problems governed by elliptic partial differential equations. To this end, let \(\Omega \subset \mathbb{R}^n\) be a bounded Lipschitz domain. We consider the following distributed control problem governed by a linear elliptic convection–diffusion–reaction equation: \[\label{eq:elliptic} - \operatorname{div}(\kappa(\omega)\nabla x(\omega)) + b(\omega)^\top \nabla x(\omega) + c(\omega) x(\omega) = \chi_{\Omega_c}(\omega) u(\omega) \qquad \text{for }\omega \in \Omega,\tag{24}\] supplemented with suitable boundary conditions to be specified later. Here, \(x:\Omega\to \mathbb{R}\) is the state variable and \(u:\Omega \to \mathbb{R}\) denotes the distributed control acting on a subdomain which enters the equation via a characteristic function \(\chi_{\Omega_c}\) of the control domain \(\Omega_c\subset \Omega\). Throughout, we assume for the coefficients the standard regularity assumptions, that is \(\kappa\in L^\infty(\Omega;\mathbb{R}^{n\times n}),b\in L^\infty(\Omega;\mathbb{R}^n), c\in L^\infty(\Omega)\). Further, we assume that the matrix field is uniformly elliptic, that is, there is \(\underline{\kappa}>0\) such that \(\|v\|_{\mathbb{R}^n}^2 \leq v^\top \kappa(\omega)v\) for all \(v\in \mathbb{R}^n\) and \(\omega \in \Omega\). Let \(\mathcal{X}\subset H^1(\Omega)\), depending on the imposed boundary conditions. For instance, in the case of homogeneous Dirichlet boundary conditions one chooses \(\mathcal{X}= H_0^1(\Omega)\) whereas for Neumann or Robin boundary conditions one typically takes \(\mathcal{X}= H^1(\Omega)\). Define the bilinear form \(a : \mathcal{X}\times \mathcal{X}\to \mathbb{R}\) by \[\begin{align} \label{eq:secondorder95op} a(x,y)\mathrel{\vcenter{:}}= \langle \kappa\nabla x, \nabla y \rangle_{L^2(\Omega;\mathbb{R}^n)} + \langle b^\top \nabla x, y \rangle_{L^2(\Omega)} + c \langle x, y \rangle_{L^2(\Omega)} + a_{\partial \Omega}(x,y), \end{align}\tag{25}\] where the boundary contribution \(a_{\partial \Omega}\) depends on the imposed boundary conditions. For example, \(a_{\partial \Omega} \equiv 0\) in the Dirichlet and Neumann cases, whereas for Robin boundary conditions \(\partial x = \kappa x\) on \(\partial \Omega\) one has \(a_{\partial \Omega}(x,y) = \int_{\partial \Omega} \kappa x y \, ds\).

To ensure well-posedness, we assume throughout that the bilinear form \(a\) is continuous and coercive on \(\mathcal{X}\).

Our setting covers, for example, the following standard situations:

1. \(b=0\), \(c=0\) and homogeneous Dirichlet boundary conditions. In this case the coercivity constant is usually given by the Poincaré inequality.

2. \(b=0\) and \(c>0\) and Neumann boundary conditions. Here, coercivity is obtained through the reaction term.

3. convection–diffusion-type operators with \(c_0 - \frac{1}{2} \operatorname{div}(b) \ge 0\), assuming that \(b\) admits a weak divergence.

Continuity of \(a\) ensures the existence of a bounded linear operator \(\mathcal{A}: \mathcal{X}\to \mathcal{X}^*\) with \[\label{eq:opdefinedelliptic} \langle \mathcal{A}x,y \rangle_{\mathcal{X}^*,\mathcal{X}}= a(x,y) , \qquad x,y \in \mathcal{X}.\tag{26}\] Define the operator \(\mathcal{B}: L^2(\Omega) \to \mathcal{X}^*\) with \(\mathcal{B}u = \mathcal{R}_{L^2(\Omega)} \chi_{\Omega_c}\). Then the problem 24 may abstractly be written by \[\mathcal{A}x - \mathcal{B}u=0.\] The Lax-Milgram Theorem, [17], yields that \(\mathcal{A}\) is boundedly invertible. Further, one straightforwardly may compute for \(x\in \mathcal{X},f\in \mathcal{X}^*\), \[\label{eq:coerc1} \mathcal{A}x = f \quad \Longrightarrow \quad m\|x\|^2_\mathcal{X}\leq \langle \mathcal{A}x,x\rangle_{\mathcal{X}^*,\mathcal{X}} = \langle f,x\rangle_{\mathcal{X}^*,\mathcal{X}} \leq \|f\|_{\mathcal{X}^*}\|x\|_{\mathcal{X}}\tag{27}\] and hence \[\begin{align} \label{eq:coerc2} \| \mathcal{A}^{-1} \|_{\mathcal{X}^*\to \mathcal{X}} \leq \frac{1}{m}. \end{align}\tag{28}\] However, this also shows that if the coercivity constant \(m\) tends to \(0\), then the operator norm of the inverse \(\mathcal{A}^{-1}\) may become unbounded. Thus, as \(m \to 0\), the problem becomes increasingly ill-conditioned, which may lead to significant difficulties in numerical approximations. Hence, one may introduce a feedback operator to improve the coercivity constant of the feedbacked state operator, that is, to obtain \[\begin{align} \langle (\mathcal{A}- \mathcal{B}\mathcal{K}) x,x\rangle_{\mathcal{X}^*,\mathcal{X}} \geq \widetilde{m} \|x\|^2_{\mathcal{X}} \quad \forall x\in \mathcal{X}\quad \mathrm{with} \quad \widetilde{m} \gg m. \end{align}\] such that the feedbacked operator satisfies \[\begin{align} \| \widetilde{\mathcal{A}}^{-1} \|_{\mathcal{X}^*\to \mathcal{X}} = \| (\mathcal{A}- \mathcal{B}\mathcal{K})^{-1} \|_{\mathcal{X}^*\to \mathcal{X}} \leq \frac{1}{\widetilde{m}} \ll \frac{1}{m}. \end{align}\] We now illustrate this approach for a concrete example.

Distributed control of Neumann problem. We consider the Poisson problem with homogeneous Neumann boundary condition \[\begin{align} \label{eq:Neumann} \begin{aligned} -\Delta x + \gamma x & = u & \qquad \text{on } & \Omega, \\ \partial_N x & = 0 & \qquad \text{on } & \partial\Omega. \end{aligned} \end{align}\tag{29}\] with \(\gamma > 0\) and where \(\partial_N\) denotes the normal derivative of \(x\). Note that this corresponds to 24 by choosing \(\kappa \equiv I, b\equiv 0, c \equiv \gamma\) and \(\Omega_c = \Omega\). The bilinear form associated with 29 is given by \(a: H^1(\Omega) \times H^1(\Omega) \to \mathbb{R}\) with \[a(x,y) = \langle \nabla x, \nabla y \rangle_{L^2(\Omega;\mathbb{R}^n)} + \gamma \langle x,y \rangle_{L^2(\Omega)}\] We note that \[\begin{align} \label{eq:coerc95ex} a(x,x) = \|\nabla x\|^2_{L^2(\Omega;\mathbb{R}^d)} + \gamma\|x\|^2_{L^2(\Omega)} \ge \min\{\gamma,1\}\|x\|^2_{H^1(\Omega)}. \end{align}\tag{30}\] For constant functions \(z\in H^1(\Omega)\), we get \(a(z,z) = \gamma \|x\|^2_{L^2(\Omega)} = \gamma \|x\|^2_{H^1(\Omega)}\) such that the coercivity constant in 30 is optimal in the case \(\gamma < 1\) in the sense that it can not be chosen larger. In view of 27 and 28 , this yields that the associated operator \(\mathcal{A}: H^1(\Omega) \to H^1(\Omega)^*\) defined via 26 satisfies the solution operator norm \[\label{eq:elliptic95original} \| \mathcal{A}^{-1} \|_{H^{1}(\Omega)^*\to H^{1}(\Omega)} = \frac{1}{\gamma}, \qquad \mathrm{if} \;\gamma < 1.\tag{31}\] Hence, if \(\gamma>0\) is small, the solution operator norm becomes very large; in the limit \(\gamma \to 0\), \(\mathcal{A}\) even lacks to be invertible due to a kernel consisting of the subspace of constant functions. We now show how to improve this solution operator norm by means of a feedback transformation. To this end, let \(\mathcal{K}: L^2(\Omega) \to L^2(\Omega)\) with \[\label{eq:feedback95elliptic} \mathcal{K}x = -(-\gamma + 1)x = (\gamma-1) x\tag{32}\] Substituting \(u=\mathcal{K}x +v\), the control problem 29 is equivalently rewritten by \[\begin{align} \label{eq:ellipticfeedbacked} 0 = -\Delta x +\gamma x - u = -\Delta x + \gamma x - \mathcal{K}x -v = -\Delta x + x -v \qquad \text{on } \Omega . \end{align}\tag{33}\] This new control problem is now driven by \(\widetilde{\mathcal{A}}: H^1(\Omega) \to H^{1}(\Omega)^*\) with \[\langle \widetilde{\mathcal{A}}x, y\rangle_{H^1(\Omega)^*,H^1(\Omega)} = \langle (\mathcal{A}- \mathcal{B}\mathcal{K})x, y\rangle_{H^1(\Omega)^*,H^1(\Omega)} = \langle \nabla x, \nabla y \rangle_{L^2(\Omega;\mathbb{R}^n)} + \langle x,y \rangle_{L^2(\Omega)},\] which clearly satisfies \(\langle \widetilde{\mathcal{A}}x,x\rangle_{H^1(\Omega)^*,H^1(\Omega)} = \|x\|^2_{H^1(\Omega)}\) such that \[\label{eq:elliptic95feedbacked} \|\widetilde{\mathcal{A}}^{-1}\|_{H^1(\Omega)^*\to H^1(\Omega)} = 1,\tag{34}\] that is, the feedbacked solution operator is bounded independently of \(\gamma\). Even the limiting case \(\gamma = 0\), in which \(\mathcal{A}\) is surjective but not injective may thus be treated using a gradient or PPCG method, as we only require the feedbacked control-to-state map to solve 3 .

To apply Propositions 4 and 6, we remain to verify Assumption 1. Assumption 1[ass:iso], i.e., continuous invertibility of \(\widetilde{\mathcal{A}}\) follows directly from 34 . Assumption 1[ass:smallness] is clearly satisfied with \(\delta_\mathcal{K}= (1-\gamma)>0\). Hence, both Proposition 4 and 6 are applicable.

In Figure 4, we show the results for the condensed conjugate gradient method as outlined and analyzed in Subsection 3.1. We solve the problem 1 subject to the elliptic PDE 29 on the unit square \([-1,1]^2\subset \mathbb{R}^2\), \(\mathcal{Y}= L^2(\Omega)\), \(\mathcal{C}= I\), \(y_\mathrm{ref}(\omega_1,\omega_2) = \sin(2\pi\omega_1) + \sin(2\pi \omega_2)\) for varying regularization parameters \(\alpha\) in the cost functional and \(\gamma > 0\) in the PDE. As a comparison, we solve the feedbacked counterpart 3 using the feedback as defined in 32 . Combining the deduced condition number estimates of Proposition 1 for the original problem and Proposition 6 for its feedbacked version and the solution operator estimates 31 and 34 , we see that the condition number for the original problem deteriorates for \(\gamma \to 0\), while it remains uniformly (in \(\gamma\)) bounded for the feedbacked problem. This may clearly be observed in Figure 4. When decreasing the parameter \(\gamma\) (left to right), we observe that the convergence behavior of the original problem becomes worse and worse (for both choices of the regularization parameter depicted in the upper and lower row). In contrast, the convergence behavior of the feedbacked problem is visibly independent of \(\gamma\). This clearly reflects the uniform bound on the solution operator 34 .

In Figure 5 we show the results for the PPCG method applied to the same problem. We observe the same behavior as for the gradient method: For decreasing parameters \(\gamma > 0\), the convergence of the PPCG method deteriorates, while for the feedbacked problem it remains unchanged. This again reflects the theoretical findings corresponding to the condition number bounds of Propositions 2 and 6 together with the solution operator norms 31 and 34 .

Optimization of the condition number ?? . Here, we briefly discuss the extreme case in which \(\gamma = 0\). In this case, \(\mathcal{A}\) is not invertible. However, one may easily show that the feedbacked problem is invertible, as \[\sigma_{-\delta I} = \|(\mathcal{A}- \mathcal{B}\mathcal{K})^{-1}\|_{H^1(\Omega)^*\to H^1(\Omega)} = \max\{\delta^{-1},1\}\] when using the feedback \(\mathcal{K}= -\delta I\) for some \(\delta>0\). Further, in this case, we may minimize the expression ?? over \(\delta\) and obtain \(\delta_{\min} = \frac{1}{\sqrt{2\alpha}}\) with associated minimal value \(\kappa_{\min}=9.\)

In Figure 6, we illustrate the performance of the feedback transformation for different choices of \(\alpha\). Therein, in the respective left figures plot the error of the PPCG iteration after 5 iterations over different choices of \(\delta\) and the right figures plot the upper bound in ?? . We clearly observe a non-monotone behavior of the performance with respect to the chosen feedback. Further, the optimal performance matches quite well the optimization of the condition number estimate (depicted with a vertical red line), despite being only an upper bound.

4.3 Parabolic PDEs↩︎

In this part, we provide an application to parabolic problems such as the heat equation, the advection-reaction-diffusion equation, Oseen equations or beam equation with viscous damping. To this end, let \(V\) be a separable and reflexive Banach space that is densely and continuously embedded into a Hilbert space \(H\) giving rise to a Gelfand triple \(V\hookrightarrow H \hookrightarrow V^*\), where \(\|v\|_H \le \|v\|_V \quad \forall v\in V\). We assume that the main operator \(A\in L(V,V^*)\) satisfies the generalized ellipticity or Gårding inequality \[\begin{align} \label{eq:garding} \exists c \in \mathbb{R}, m > 0, \beta \ge 0: \quad \langle Av,v\rangle_{V^*,V} + c \|v\|^2_H \geq \|v\|_{m,\beta}^2 := m \|v\|^2_V+\beta\|v\|^2_H. \end{align}\tag{35}\] Observe that we have some freedom in the choice of \(m\), \(c\), and \(\beta\), which will turn out to be useful later on. The use of these weighed norms will allow us to obtain quantitative results on the norm of the solution operator.

For example, \(A\) could be a second-order operator in weak form induced by a bilinear form 25 , possibly with feed-back. Throughout this section, we abbreviate \(\|A\|_{V\to V^*}\) by \(\|A\|\) and write \(\dot{x}\) for the weak time-derivative of \(x\). As a solution space, we consider the space \[\begin{array}{rcl} W(0,T) & \mathrel{\vcenter{:}}= & \left\lbrace x \in L^2(0,T; V) \, \middle| \, \dot{x} \in L^2(0,T; V^*) \right\rbrace. \end{array}\] Note that we have \(W(0,T) \hookrightarrow C([0,T]; H)\) continuously [18] such that point evaluations, e.g., to prescribe initial values, may be defined.

A parabolic partial differential equation is then given by the initial value problem \[\begin{align} \label{eq:parab} \mathrm{Find}\;x\in W(0,T):\;\dot{x} + Ax = g, \qquad x(0)=x^0, \end{align}\tag{36}\] where \(g\in L^2(0,T;V^*)\) is a source term and \(x^0\in H\) is an initial datum. To put this into perspective w.r.t.the abstract optimal control problem 1 , let \(B\in L(U,V^*)\) be a control operator and \(u\in L^2(0,T;U)\) a control function. Then, for given \(l \in L^2(0,T;V^*), x^0\in H\), we may set \[\begin{align} \mathcal{A} x := \begin{bmatrix} \dot{x} + Ax \\ x(0) \end{bmatrix}, \quad \mathcal{B}u = \begin{bmatrix} Bu\\ 0 \end{bmatrix}, \quad f = \begin{bmatrix} g\\ x^0 \end{bmatrix} \end{align}\] such that \(\mathcal{X}=W(0,T)\), \(U = L^2(0,T;U)\), \(P^*=L^2(0,T;V^*)\times H\). It is well-known by parabolic theory [1], [18], [19], that under the above assumptions and for given \(u\in L^2(0,T;U)\), the problem \(\mathcal{A}x + \mathcal{B}u = f\) has a unique solution \(x\in W(0,T)\). For a concise proof that \(\mathcal{A}\) is indeed an isomorphism, we refer to [15]. In this work, we are particularly interested in norms of the solution operator and suitably feedbacked counterparts. To alleviate notation, we define the following scaled norms \[\|x\|_{L^2_{c,m,\beta}} := \|e^{-c t}\|x(t)\|_{m,\beta}\|_{L^2(0,T)}, \qquad \|f\|_{L^2_{c,m,\beta},*}:= \|e^{-c t}\|f(t)\|_{m,\beta,*}\|_{L^2(0,T)},\] where \(\|\cdot\|_{m,\beta}\) is defined in 35 and \(\|\ell\|_{m,\beta,*} := \sup_{\|v\|_{m,\beta}\le 1} |\ell(v)|\).

We briefly comment on the consequence of the estimate ?? that directly implies an estimate on the standard norm due to \[\begin{align} \|e^{-ct}x\|_{L_2(0,T;V)} &\le \tfrac{1}{\sqrt{m}}\|x\|_{L^2_{c,m,\beta}}, \quad \|e^{-ct}x\|_{_{L_2(0,T;H)}} \le \tfrac{1}{\sqrt{m+\beta}} \|x\|_{L^2_{c,m,\beta}} \\ \|g\|_{L^2_{c,m,\beta},*} &\le \tfrac{1}{\sqrt{m}}\|e^{-ct}g\|_{L_2(0,T;V^*)}, \quad \|g\|_{L^2_{c,m,\beta},*} \le \tfrac{1}{\sqrt{m+\beta}} \|e^{-ct}g\|_{L_2(0,T;H)} \end{align}\] for \(x\in L^2(0,T;V)\) and \(f\in L^2(0,T;H)\).

Let us consider now the case of distributed control and observation, i.e., \(\mathcal{U}= L^2(0,T;H)\) and \(B\in L(H,H)\). As an immediate result, we obtain the following.

Next, we provide a concrete example to illustrate the results.

Heat equation with Neumann boundary conditions. We consider the very same setup as in 29 , but now allow for \(\gamma \in \mathbb{R}\), that is, also negative values. Here we have that 35 holds for \(A=-\Delta + \gamma I\) (associated with \(\mathcal{A}\)) with \(c=0\) and \(m=\min\{\gamma,1\}\) if \(\gamma >0\) and with \(c= - \gamma + 1>0\) and \(m=1\) for \(\gamma \leq 0\). In the former case, the uncontrolled system is exponentially stable; in the latter it is stable (if \(\gamma = 0\)) or exponentially unstable (if \(\gamma <0\)). For the governing operator of \(\widetilde{\mathcal{A}}\), i.e. \(A-BK\), 35 always holds with \(c=0\) by the choice of the feedback due to 33 .

In Figure 7, we illustrate the results for the condensed gradient method. For \(\gamma>0\) (left), the problem is exponentially stable. However, for \(\gamma < 0\), the (uncontrolled) problem corresponding to \(\mathcal{A}\) is exponentially unstable. Hence, when increasing the horizon \(T=10\) to \(T=30\) (from top to bottom), the convergence behavior significantly worsens for the original problem. This reflects the solution operator norm of Proposition 10 depending exponentially on the time horizon. and its impact on the condition number in Proposition 1. For the feedbacked problem, the convergence behavior again is robust w.r.t.an increased time horizon. Hereby we note, however, that the iteration numbers increase (in contrast to the elliptic case of Subsection 4.2), as an increased horizon leads to a higher number of time discretization points.

In Figure 8 we report the iteration numbers for the PPCG method for the fixed regularization parameter \(\alpha = 10^{-3}\) here. However, we see again that while the feedbacked version is very robust with respect to increasing the horizon length, the original problem faces severe numerical difficulties in particular for the unstable case \(\gamma = -0.5\) depicted in the lower row.

Optimization of the condition number ?? . We consider the heat equation with Neumann boundary conditions and with \(\mathcal{C}=\mathcal{B}= I\). When choosing the feedback \(\mathcal{K}= -\delta I\) for \(\delta > 0\) and by a similar argumentation as for the finite-dimensional case proven in Proposition 8 replacing the matrix exponential by the operator semigroup, one can derive the solution operator norm

If \(\delta \to 0\) then the first term approaches \(T\) (which may be seen using L’Hôpital’s rule) that is typical for stable (but not asymptotically stable) systems due to integration of a constant function over \([0,T]\); a similar bound will be observed for the wave equation in the next subsection. For any \(\delta > 0\), however, this term is bounded uniformly in \(T\) and for \(\delta \to \infty\) it approaches one. Again, we may plug this relation into ?? and optimize over \(\delta\). This is illustrated in Figure 9.

4.4 Hyperbolic PDEs↩︎

Let \(\Omega \subset \mathbb{R}^n\) be a bounded open domain. Consider the controlled wave equation with homogeneous Dirichlet boundary conditions and prescribed initial data \[\label{eq:wave} \begin{array}{rcll} x_{tt}(t,\omega) - \Delta x(t,\omega) &=& u(t, \omega) & \qquad \text{on } [0,T] \times \Omega \\ x(t,\omega) &=&0 & \qquad \text{on } [0,T] \times \partial \Omega \\ x(0,\omega) &=&x_0(\omega) & \qquad \text{on } \Omega \\ x_t(0,\omega) &=&x_1(\omega) & \qquad \text{on } \Omega . \end{array}\tag{38}\] Define the operator \(A: L^2(\Omega) \supset H^2(\Omega) \cap H_0^1(\Omega) \to L^2(\Omega)\) with \(A x = - \Delta x\). Note that \(A\) is self-adjoint and uniformly positive and there exists an operator square root \(A^{1/2} : L^2(\Omega) \supset H_0^1(\Omega) \to L^2(\Omega)\) where \(H_0^1(\Omega)\) is endowed with the energy inner product \(\langle \cdot,\cdot \rangle_{H_0^1(\Omega)} =\langle A^{1/2} \cdot,A^{1/2}\cdot \rangle_{L^2(\Omega)}\). Let \(H^{-1}(\Omega)\) be the dual space of \(H_0^1(\Omega)\) with respect to the pivot space \(L^2(\Omega)\). Precisely, \((H_0^1(\Omega), L^2(\Omega), H^{-1}(\Omega))\) defines a Gelfand triple and there exists a unique bounded extension of \(A\) from \(H_0^1(\Omega)\) to \(H^{-1}(\Omega)\) that coincides with the Riesz isomorphism \(\mathcal{R}_{H_0^1(\Omega)}\), again denoted by \(A: H_0^1(\Omega) \to H^{-1}(\Omega)\). For a Banach space \(Z\), denote by \(D'(0,T;Z)\) the space of \(Z\)-valued distributions on \((0,T)\), i.e., \(D'(0,T;Z)\mathrel{\vcenter{:}}= L(C_c^\infty(0,T),Z)\). Since every continuous \(H_0^1(\Omega)\)-valued function defines an \(H_0^1(\Omega)\)-valued distribution, the second distributional time derivative \(\ddot x\in D'(0,T;H_0^1(\Omega))\) is well-defined for every \(x\in C([0,T];H_0^1(\Omega))\). Due to the continuous embedding \(H_0^1(\Omega)\hookrightarrow H^{-1}(\Omega)\), we may equally regard \(\ddot x\in D'(0,T;H^{-1}(\Omega)).\) Because \(A\in L(H_0^1(\Omega),H^{-1}(\Omega))\), the pointwise mapping \(t\mapsto Ax(t)\) belongs to \(C([0,T];H^{-1}(\Omega))\subset D'(0,T;H^{-1}(\Omega))\). Consequently, \[\tfrac{\mathrm d^2}{\mathrm dt^2}+A: C([0,T];H_0^1(\Omega)) \rightarrow D'(0,T;H^{-1}(\Omega))\] is well-defined. Define the space and norm \[\begin{align} E(0,T;A) &\mathrel{\vcenter{:}}= \left\lbrace x \in C([0,T];H_0^1(\Omega)) \cap C^1([0,T];L^2(\Omega)) : \ddot x + A x \in L^2(0,T;L^2(\Omega)) \right\rbrace\\ \| x \|_{E(0,T;A)} &\mathrel{\vcenter{:}}= \| x \|_{C([0,T];H_0^1(\Omega))} + \| \dot{x} \|_{C([0,T];L^2(\Omega))} + \| \ddot x + A x\|_{L^2(0,T;L^2(\Omega))}. \end{align}\] which forms a Banach space. Indeed, completeness follows by identifying \(E(0,T;A)\) with the corresponding first-order energy space associated with the abstract wave operator \({ \left[\begin{smallmatrix} 0 & I \\ -A & 0 \end{smallmatrix}\right]}\). Let \(\mathcal{P}\mathrel{\vcenter{:}}= L^2(0,T;L^2(\Omega)) \times H_0^1(\Omega) \times L^2(\Omega)\) and define \(\mathcal{A}: E(0,T;A) \to \mathcal{P}\) and \(\mathcal{B}: L^2(0,T;L^2(\Omega)) \to \mathcal{P}\) by \[\mathcal{A}x = { \left(\begin{smallmatrix} \ddot x + A x \\ x(0) \\ \dot{x}(0) \end{smallmatrix}\right)} , \qquad \mathcal{B}u = { \left(\begin{smallmatrix} u \vphantom{A} \\ 0 \vphantom{(0)} \\ 0 \vphantom{(0)} \end{smallmatrix}\right)} .\] Introducing the first-order state variable \((x,v) \mathrel{\vcenter{:}}= (x, \dot{x})\) (displacement-velocity formulation) transforms the second-order equation into an abstract Cauchy problem generated by the associated wave operator. Consequently, \(\mathcal{A}\) is an isomorphism, see e.g. [16]. We rewrite the hyperbolic control problem 38 by \[\label{eq:abstract95wave} \mathcal{A}x - \mathcal{B}u = { \left[\begin{smallmatrix} 0 \\ x_0 \\ x_1 \end{smallmatrix}\right]} .\tag{39}\] For the cost functional, we introduce position and velocity observation \[\mathcal{C}={ \left[\begin{smallmatrix} I \\ \frac{\mathrm{d}}{\mathrm{d}t} \end{smallmatrix}\right]} : E(0,T;A) \to L^2(0,T;H_0^1(\Omega) \times L^2(\Omega)).\]

The following result shows a solution operator bound for the uncontrolled problem. In particular, we observe a dependence on the time horizon \(T\) which (in constrast to the previously discussed parabolic case) is not exponential but linear; this is due to the energy conservation property of the wave equation and essentially stems from integrating a constant norm (energy) over the time horizon \(T\). The proof is presented in Appendix 5.3.

The previous theorem shows that \(\mathcal{C}\mathcal{A}^{-1}\) is not uniformly bounded in the time horizon. We now show how to suitably design a feedback stabilization \(\widetilde{\mathcal{A}}= \mathcal{A}- \mathcal{B}\mathcal{K}\) such that \(\mathcal{C}\widetilde{\mathcal{A}}^{-1}\) is uniformly bounded in the time horizon. As the observation includes both position and velocity, in view of Assumption 1[ass:smallness] we may use the feedback operator \[\mathcal{K}: E(0,T;A) \to L^2(0,T;L^2(\Omega)), \qquad \mathcal{K}x = -(x + \dot{x})\] which induces linear friction in the model. The resulting closed loop system reads \(\ddot x + \dot{x} - \Delta x + x =0\) and in displacement-velocity formulation, \[\begin{bmatrix} \dot{x} \\ \dot{v} \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ \Delta -1& - 1 \end{bmatrix} \begin{bmatrix} x \\v \end{bmatrix}.\] The next result shows that using this feedback operator, the transformed system is exponentially stable and in particular, the feedbacked constraint operator is uniformly boundedly invertible in the time horizon. The proof is presented in Appendix 5.3.

We illustrate the results for the gradient method in Figure 10. As the feedbacked system has a condition number bounded uniformly in the time horizon (due to Theorem 12) and the original system has a bound depending on \(T\) (see Theorem 11), we observe that when increasing the time horizon (left to right), the gradient method applied to the feedbacked problem more and more outperforms the gradient method for to the original problem.

The results for the PPCG method are shown in Figure 11 for \(\alpha = 10^{-3}\). Here, we observe a similar behavior, that is, while both formulations perform equally well for small time horizons, the difference becomes more and more significant as we increase the time horizon.

5 Supplementary proofs↩︎

5.1 Proof of 12↩︎

We first consider the case \(a\neq 1\). To obtain a lower bound, we apply the Toeplitz matrix \(T_N\) to the first canonical basis vector, which gives the lower bound \(\|T_N\|_2 \geq \|T_N e_1\|_2 = \left(\sum_{k=0}^{N-1} a^{2k}\right)^{1/2}\). From \(\|T_N\|_2 \leq \sqrt{\|T_N\|_1\|T_N\|_\infty}\) and \(\|T_N\|_1 = \|T_N\|_\infty = \frac{1-a^N}{1-a}\) we get an upper bound for the spectral norm. Putting both together for \(a\neq 1\) yields \(\sqrt{\tfrac{1-a^{2N}}{1-a^2}} \;\leq\; \|T_N\|_2 \;\leq\; \tfrac{1-a^N}{1-a}\) that is, the claim. We now consider the case \(a=1\). For the upper bound, we may estimate \(\|T_N\|_2 \leq \|T_N\|_F = \sqrt{N^2} = N\). For the lower bound, we choose \(\mathbf{1}:=(1,\ldots,1)^\top\in\mathbb{R}^N\) for which we have \(T_N\mathbf{1} = (1,2,\ldots,N)^\top,\) leads to \(\|T_N\|_2 \ge \tfrac{\|T_N\mathbf{1}\|_2}{\|\mathbf{1}\|_2} = \left( \tfrac{\sum_{k=1}^N k^2}{N} \right)^{1/2} = \left( \tfrac{(N+1)(2N+1)}{6} \right)^{1/2}.\) As a result \(\|T_N\|_2 \ge \frac{N}{\sqrt3}\) which together with the upper bound derived above yields the result.

5.2 Proof of Proposition 8↩︎

Let \(\gamma_0 : H^1(0,T; \mathbb{R}^n) \to \mathbb{R}^n\) be the boundary evaluation map, i.e., \(\gamma_0(x) = x(0)\). It is well-known that \(\gamma_0\) is a bounded operator, see [20]. Hence, it is enough to prove boundedness of \(\tfrac{\mathrm{d}}{\mathrm{d}t}-A: H^1(0,T; \mathbb{R}^n) \to L^2(0,T; \mathbb{R}^n)\) which is immediate from \[\| \tfrac{\mathrm{d}}{\mathrm{d}t}x-Ax \|_{L^2(0,T; \mathbb{R}^n)}^2 \leq \| \tfrac{\mathrm{d}}{\mathrm{d}t}x \|_{L^2(0,T; \mathbb{R}^n)}^2 + \| Ax \|_{L^2(0,T; \mathbb{R}^n)}^2 \leq \operatorname{max}(1, \|A\|^2) \| x \|_{H^1(0,T; \mathbb{R}^n)}^2\] for all \(x \in H^1(0,T; \mathbb{R}^n)\). We proceed to show bijectivity of \(\mathcal{A}\). To this end, let \(x \in \ker \mathcal{A}\). In view of formula 17 , we obtain \(x(t) = e^{tA}x_0 = e^{tA}0 =0,\) which proves \(x=0\). Furthermore, let \((f, x_0) \in L^2(0,T; \mathbb{R}^n) \times \mathbb{R}^n\) and denote by \(x \in H^1(0,T; \mathbb{R}^n)\) the unique mild solution of \[\dot{x}(t) = A x(t) + f(t) , \qquad x(0) = x_0.\] Then, \(\mathcal{A}x = { \left[\begin{smallmatrix} f \\ x_0 \end{smallmatrix}\right]}\) and hence, \(\operatorname{ran}(\mathcal{A})= L^2(0,T; \mathbb{R}^n) \times \mathbb{R}^n\). Hence, there exists \(\mathcal{A}^{-1}: L^2(0,T; \mathbb{R}^n) \times \mathbb{R}^n \to H^1(0,T; \mathbb{R}^n)\) given by \[\mathcal{A}^{-1} { \left[\begin{smallmatrix} f \\ x_0 \end{smallmatrix}\right]} = e^{\cdot A}x_0 + \int_0^\cdot e^{(\cdot-s)A} f(s)\,\mathrm{d}s \eqqcolon x(\cdot).\] We remain to show boundedness of \(\mathcal{A}^{-1}\). First, assume that \(\omega \neq 0\). Due to 18 we observe \[\begin{align} \int_0^{T} \left\|e^{t A}x_0 \right\|_{\mathbb{R}^n}^2 \, \mathrm{d}t \leq M^2 \left\|x_0 \right\|_{\mathbb{R}^n}^2 \int_0^{T} e^{2t\omega} \, \mathrm{d}t \leq M^2\frac{e^{2T\omega}-1}{2\omega} \left\|x_0 \right\|_{\mathbb{R}^n}^2 . \end{align}\] Moreover, tedious but straightforward computations show that \[\begin{align} \int_0^T \!\!\int_0^t \|e^{(t-s)A}f(s)\|^2 \,\mathrm ds\,\mathrm dt &\le M^2 \int_0^T \!\!\int_0^t e^{2(t-s)\omega}\|f(s)\|^2 \,\mathrm ds\,\mathrm dt \\ &= M^2 \int_0^T \|f(s)\|^2 \left(\int_s^T e^{2(t-s)\omega}\,\mathrm dt\right) \,\mathrm ds \\ &= M^2 \int_0^T \tfrac{e^{2(T-s)\omega}-1}{2\omega} \|f(s)\|^2 \,\mathrm ds \le M^2 \tfrac{e^{2T\omega}-1}{2\omega} \|f\|_{L^2(0,T;\mathbb{R}^n)}^2 . \end{align}\] Combining these two estimates yields the \(L^2\)-estimate \[\| x \|_{L^2(0,T; \mathbb{R}^n)}^2 \leq M^2 \tfrac{e^{2T\omega} -1}{\omega} \left( \| x_0 \|_{\mathbb{R}^n}^2 +\| f\|_{L^2(0,T;\mathbb{R}^n)}^2 \right).\] For the full \(H^1\)-estimate, we compute \[\begin{align} \| \mathcal{A}^{-1} { \left[\begin{smallmatrix} f \\ x_0 \end{smallmatrix}\right]} \|_{H^1(0,T; \mathbb{R}^n)}^2 &= \| x \|_{L^2(0,T; \mathbb{R}^n)}^2 + \| \dot{x} \|_{L^2(0,T; \mathbb{R}^n)}^2 \\ &= \| x \|_{L^2(0,T; \mathbb{R}^n)}^2 + \| A x + f\|_{L^2(0,T; \mathbb{R}^n)}^2 \\ & \leq \left(1 + 2 \|A\|^2\right) \| x \|_{L^2(0,T; \mathbb{R}^n)}^2 + 2 \| f\|_{L^2(0,T; \mathbb{R}^n)}^2 \\ &\leq \left( \tfrac{M^2}{\omega} \left(1 + 2 \|A\|^2\right) (e^{2T\omega}-1)+2\right) \left(\left\|x_0 \right\|_{\mathbb{R}^n}^2 + \| f \|_{L^2(0,T; \mathbb{R}^n)}^2\right) . \end{align}\] The case \(\omega = 0\) is derived analogously.

\(\square\)

5.3 Proofs of Theorem 11 and Theorem 12↩︎

In this proof, we abbreviate the Gelfand triple \((H_0^1(\Omega), L^2(\Omega), H^{-1}(\Omega))\) by \((H_{1/2}, H,H_{-1/2})\). The operator \(\mathcal{A}\) is bounded by construction of the solution space \(E(0,T;A)\) and continuity of the evaluation operators \(x \mapsto x(0)\) and \(x \mapsto \dot{x}(0)\). Let \(D \in L(H)\) fulfill \(\langle Dv,v\rangle_H \geq 0\) for all \(v \in H\). Let \(G: X \mathrel{\vcenter{:}}= H_{1/2} \times H \supset \operatorname{dom}(G) \to X\) be defined by \[\operatorname{dom}(G) = \left\lbrace (x,v) \in H_{1/2} \times H_{1/2} : Ax \in H \right\rbrace, \qquad G { \left[\begin{smallmatrix} x \\ v \end{smallmatrix}\right]} = { \left[\begin{smallmatrix} v \\ -Ax - D v \end{smallmatrix}\right]}\] It is well-known, that is \(G\) maximally dissipative and hence the generator of a contraction semigroup \((T(t))_{t \geq 0}\) on \(X\). Observe that Theorem 11 corresponds to the choice \(D=0\), whereas Theorem 12 is obtained for \(D=I\). For \((x(0),v(0))=(x_0,v_0) \in X\), the unique mild solution of the inhomogeneous Cauchy problem \[{ \left[\begin{smallmatrix} \dot{x}(t) \\ \dot{v}(t) \end{smallmatrix}\right]} = G { \left[\begin{smallmatrix} x(t) \\ v(t) \end{smallmatrix}\right]} + { \left[\begin{smallmatrix} f_1(t) \\ f_2(t) \end{smallmatrix}\right]} , \;\text{is given by } { \left[\begin{smallmatrix} x(t) \\ v(t) \end{smallmatrix}\right]} \mathrel{\vcenter{:}}= T(t){ \left[\begin{smallmatrix} x_0 \vphantom{(0)} \\ x_1 \vphantom{(0)} \end{smallmatrix}\right]} + \int_0^t T(t-s){ \left[\begin{smallmatrix} f_1(s) \\ f_2(s) \end{smallmatrix}\right]} \, \mathrm{d}s.\] Since the above equation is the first order (in time) formulation of 38 (with \(f_1=0,f_2 =f\)), we obtain \[\label{eq:solution95hyperb} x(t) = \pi_1 T(t){ \left[\begin{smallmatrix} x_0 \vphantom{(0)} \\ x_1 \vphantom{(0)} \end{smallmatrix}\right]} + \int_0^t \pi_1 T(t-s){ \left[\begin{smallmatrix} 0 \\ f(s) \end{smallmatrix}\right]} \, \mathrm{d}s\tag{40}\] where \(\pi_1 : X \to H_{1/2}\) denotes the projection onto the first component. Note that this map fulfills \(\| \pi_1 \|_{X \to H_{1/2}}=1\). Let \((f,x_0,x_1) \in L^2(0,T;H) \times H_{1/2} \times H\) be arbitrary. A combination of the above and [21] yields that the unique solution \(x \in E(0,T;A)\) of \[\mathcal{A}x = { \left[\begin{smallmatrix} f \\ x_0 \\ x_1 \end{smallmatrix}\right]} , \qquad \widetilde{\mathcal{A}}x\mathrel{\vcenter{:}}= (\mathcal{A}- \mathcal{B}\mathcal{K})x={ \left[\begin{smallmatrix} f \\ x_0 \\ x_1 \end{smallmatrix}\right]}\] is given by 40 , respectively. Hence, both \(\mathcal{A}\) and \(\widetilde{\mathcal{A}}=\mathcal{A}- \mathcal{B}\mathcal{K}\) are invertible and straightforward computations show that the inverses are continuous. We now show the upper bounds on \(\| \mathcal{C}\mathcal{A}^{-1} \|, \| \mathcal{C}\widetilde{\mathcal{A}}^{-1} \|\), respectively. It is well-known, that if \(D=0\), then \(G\) is skew-adjoint and hence the generator of a unitary group \((T(t))_{t \geq 0}\) on \(H_{1/2} \times H\). In this case, \[\begin{align} \left|\left|\mathcal{C}\mathcal{A}^{-1}{ \left[\begin{smallmatrix} f \\ x_0 \\ x_1 \end{smallmatrix}\right]} \right|\right|_{L^2(0,T;H_{1/2} \times H)}^2 &= \int_0^T \|x(t)\|_{H_{1/2}}^2 +\|\dot{x}(t)\|_{H}^2 \, \mathrm{d}t \\ &= \int_0^T \left\| T(t){ \left[\begin{smallmatrix} x_0 \\ x_1 \end{smallmatrix}\right]} + \int_0^t T(t-s){ \left[\begin{smallmatrix} 0 \\ f(s) \end{smallmatrix}\right]} \, \mathrm{d}s\right\|_{X}^2 \, \mathrm{d}t \\ &\leq \int_0^T 2 \| { \left[\begin{smallmatrix} x_0 \\ x_1 \end{smallmatrix}\right]} \|_X^2 \, \mathrm{d}t + \int_0^T t \int_0^t 2 \| f(s) \|_H^2 \, \mathrm{d} s \, \mathrm{d}t \\ &\leq \operatorname{max}(2T,T^2) \left( \| { \left[\begin{smallmatrix} x_0 \\ x_1 \end{smallmatrix}\right]} \|_X^2 + \| f \|_{L^2(0,T;H)}^2 \right) \end{align}\] and hence ?? . Conversely, if \(D=I\), then \(G\) is the generator of an exponentially stable semigroup \((T(t))_{t \geq 0}\) on \(H_{1/2} \times H\), see e.g. [22], [23]. Hence, there exist constants \(M \geq 0\), \(\omega>0\) such that \(\|T(t)\|_{H_{1/2}\times H \to H_{1/2}\times H} \leq Me^{-\omega t}\). Therefore, \[\left\| \int_0^t T(t-s) { \left[\begin{smallmatrix} 0 \\ f(s) \end{smallmatrix}\right]} \,\mathrm{d}s \right\|_X \leq M \int_0^t e^{-\omega(t-s)} \|f(s)\|_H\,\mathrm{d}s.\] Squaring and applying the weighted Cauchy-Schwarz inequality yields \[\begin{align} \left( \int_0^t e^{-\omega(t-s)} \|f(s)\|_H\,\mathrm{d}s \right)^2 &= \left(\int_0^te^{-\omega(t-s)/2}e^{-\omega(t-s)/2} \|f(s)\|_H \,\mathrm{d}s \right)^2\\ &\leq \left( \int_0^t e^{-\omega(t-s)} \,\mathrm{d}s \right) \left(\int_0^t e^{-\omega(t-s)}\|f(s)\|_H^2 \,\mathrm{d}s\right). \end{align}\] Since \(\int_0^t e^{-\omega(t-s)}\,\mathrm{d}s = \frac{1-e^{-\omega t}}{\omega}\leq \frac{1}{\omega},\) it follows that \[\left\| \int_0^t T(t-s) { \left[\begin{smallmatrix} 0 \\ f(s) \end{smallmatrix}\right]} \,\mathrm{d}s \right\|_X^2 \leq \frac{M^2}{\omega} \int_0^t e^{-\omega(t-s)} \|f(s)\|_H^2 \,\mathrm{d}s.\] Consequently, \[\begin{align} \left|\left|\mathcal{C}\widetilde{\mathcal{A}}^{-1}{ \left[\begin{smallmatrix} f \\ x_0 \\ x_1 \end{smallmatrix}\right]} \right|\right|^2 &= \int_0^T \|x(t)\|_{H_{1/2}}^2 +\|\dot{x}(t)\|_{H}^2 \, \mathrm{d}t \\ &= \int_0^T \left\| T(t){ \left[\begin{smallmatrix} x_0 \\ x_1 \end{smallmatrix}\right]} + \int_0^t T(t-s){ \left[\begin{smallmatrix} 0 \\ f(s) \end{smallmatrix}\right]} \, \mathrm{d}s\right\|_{X}^2 \, \mathrm{d}t \\ &\leq \int_0^T 2M^2e^{-2\omega t} \| { \left[\begin{smallmatrix} x_0 \\ x_1 \end{smallmatrix}\right]} \|_X^2 \, \mathrm{d}t + \frac{2M^2}{\omega} \int_0^T \int_0^t e^{-\omega(t-s)} \|f(s)\|_H^2 \,\mathrm{d}s\, \mathrm{d}t \\ &\leq M^2 \tfrac{1-e^{-2\omega T}}{\omega} \| { \left[\begin{smallmatrix} x_0 \\ x_1 \end{smallmatrix}\right]} \|_X^2 + 2M^2 \int_0^T \| f(s) \|_H^2 \left( \int_s^T e^{-\omega (t-s)} \, \mathrm{d}t \right) \, \mathrm{d} s . \end{align}\] Since \(\int_s^T e^{-\omega (t-s)}\, \mathrm{d} t = \int_0^{T-s} e^{-\omega r}\, \mathrm{d} r \leq \int_0^\infty e^{-\omega r}\, \mathrm{d} r = \frac{1}{\omega}\) and \(T \mapsto \tfrac{1-e^{-2\omega T}}{\omega}\) is monotonically decreasing,?? follows due to \[\left|\left|\mathcal{C}\widetilde{\mathcal{A}}^{-1}{ \left[\begin{smallmatrix} f \\ x_0 \\ x_1 \end{smallmatrix}\right]} \right|\right|_{L^2(0,T;H_{1/2} \times H)}^2 \leq \max \left(\tfrac{M^2}{\omega},\tfrac{2M^2}{\omega} \right) \left| \left| { \left[\begin{smallmatrix} f \\ x_0 \\ x_1 \end{smallmatrix}\right]} \right| \right|_{L^2(0,T; H) \times H^{1/2} \times H}^2\]

\(\square\)

References↩︎