An error estimator for spectral method approximation of flow control with state constraint

1.   Introduction

There have been lots of theoretical studies on the optimal control of PDEs with state constraints, which form a foundation for its numerical approximation. Existence and uniqueness of the solutions, Lagrange multipliers, optimality conditions, and important regularity results were derived for control problems with state constraints in the pioneer work ^[1]. More discussion on these topics can be found in ^[2,3]. In respect of numerical methods, the finite element approximation of state constrained control problems has been widely investigated, and we don't try to give a detailed introduction here, one can find more works ^[4,5], including pointwise constraints, integral constraint and so on. At the same time, many numerical strategies were developed to provide efficient approximation for these control problems. Bergounioux and Kunisch ^[6] used the primal-dual active set algorithm to solve state-constrained problems. A semi-smooth Newton method was proposed to compute state-constrained control problems by De los Reyes and Kunisch ^[7]. Gong and Yan ^[8] established a mixed variational scheme for control problems with pointwise state constrains, and a direct numerical algorithm was adopted without the optimality conditions.

In recent years, spectral method has been used to approximate control problems. Despite the restriction on the higher regularity of the approximated solutions, spectral method can provide fast convergence rate and high-order accuracy with a smaller number of unknowns, which is significant to successful applications of control problems. The spectral method was considered to solve the control problems with integral state constraint, and a priori error estimates was established in ^[9]. Chen et al. ^[10,11] derived both the a priori and a posteriori error estimates for optimal control with control constraints. However, they only provided the upper bound estimation for the a posteriori error indicators, and numerical tests didn't illustrate the performance of estimators. To investigate the efficiency of adaptive strategy in the spectral method, we have to establish some successful a posteriori error estimators as in the finite element framework. However, there are not much work on these aspects to the best of our knowledge. In this work, we derive a posteriori error estimator, and prove that it can be constructed as the upper and lower bounds of approximation error.

The plan of the article is as follows. Spectral approximation of the control problem is presented, and optimality conditions of the exact problem and discretized problem are provided in the next section. We establish the a posteriori error estimator, and construct it as the upper and lower bounds of the approximation error in section 3. Numerical example confirms the theoretical result, and shows the behaviour of the indicator in section 4.

In this paper, we let $\Omega = (-1, 1)^2$ and denote ${\boldsymbol U} = L^2(\Omega)^2, {\boldsymbol Y} = H^1_0(\Omega)^2, Q = L_0^2(\Omega) = \{ q\in L^2(\Omega) \ \mid \ \int_{\Omega} q = 0\}$, where $H^1_0(\Omega)$ and $H^m(\Omega)$ with $m$ being a positive integer are usual Sobolev spaces on $\Omega$. Let $C$ denote a positive constant independent of $N$, the order of the spectral method.

2.   Spectral approximation and optimality conditions

In this section, we state the Galerkin spectral approximation and optimality conditions for the control problem with state constraint. The model under consideration is as follows: find $({\boldsymbol y}, r, {\boldsymbol u}) \in {\boldsymbol Y} \times Q \times {\boldsymbol U}$ such that

where $\boldsymbol y_0, \boldsymbol f\in L^2(\Omega)^2$, $G_{ad} = \{\boldsymbol v \in \boldsymbol U \ \mid \ \parallel {\boldsymbol v}\parallel_{0, \Omega} \ \leq d\}$, and $d$, $\upsilon$ are positive constants. It is necessary to introduce a weak formula of the state equations for spectral approximation of the optimal control. Let

Then there exist two positive constants $\sigma$ and $\delta$ such that for any ${\boldsymbol u}, \ {\boldsymbol v} \in {\boldsymbol Y}, \ q\in Q$,

By Poincare inequality, we can know that there exists a constant $\gamma > 0$ such that

Furthermore, there exists a constant $\beta > 0$ (see, e.g., ^[12]) such that

Then the standard weak formula for the state equation can be presented as: find $(\boldsymbol {y(u)}, r(\boldsymbol u))\in {\boldsymbol Y}\times Q$ such that

The problem (2.4) is well-posed by Babuška-Brezzi theorem and (2.2)-(2.3), and the control problem (2.1) can be restated as follows (OCP): find $({\boldsymbol y}, r, {\boldsymbol u}) \in {\boldsymbol Y} \times Q \times {\boldsymbol U}$ such that

We can prove existence and uniqueness of the solutions for the control problem $(\rm OCP)$ by reformulating it to a control-constrained problem. In fact, the control problem (\rm OCP) can be equivalent to the control-constrained problem

where constraint set $\mathscr{U}_{ad} = \{\boldsymbol u\in \boldsymbol U: \|G\boldsymbol u\|_{0, \Omega}\leq d \}$ with the operator $G:\boldsymbol u\rightarrow \boldsymbol y(\boldsymbol u)$. It is clear that $\mathscr{U}_{ad}$ is a closed and convex subset in $\boldsymbol U$, and the control problem (2.5) has a unique solution by the standard theorem ^[13].

2.1.   Spectral method approximation

Considering the spectral approximation of the control problem, we introduce the finite dimensional spaces ${\boldsymbol Y}_{N}$, $M_{N-2}$, and ${\boldsymbol U}_N$ to approximate spaces ${\boldsymbol Y}$, $Q$, and ${\boldsymbol U}$ respectively, where ${\boldsymbol Y}_{N} = (Q^0_N)^2$ with $Q^0_{N} = \{x_N \in Q_N: x_N| _{\partial\Omega} = 0\}$, $M_{N-2} = Q_{N-2}\cap L^2_0(\Omega)$, ${\boldsymbol U}_N = (Q_N)^2$, and $Q_N$ denotes the space of all algebraic polynomials of degree less than or equal to $N$ with respect to each single variable $x_i, \ (i = 1, 2)$. Letting $L_n(x)$ be the $n$th-degree Legendre polynomial and $\phi_k(x) = \frac{1}{\sqrt{4k+6}}(L_k(x)-L_{k+2}(x))$ (see ^[14], for example), then it holds

We now approximate the state equations as follows

It can be known from ^[12] that there exists a constant $\beta_N = O(N^{-\frac{1}{2}})$ satisfying

It follows from Babuška-Brezzi's theory that problem (2.8) is well-posed, and the Legendre Galerkin spectral approximation of $(\rm OCP)$ can be stated as ${{\text{(OCP)}}_{\text{N}}}:$

where $G_{ad}^N = \boldsymbol Y_N \cap G_{ad}$.

2.2.   Optimality conditions

We first derive optimality conditions for the exact problem (OCP) by the techniques and refined results developed in ^[2], though we can complete the derivation by other strategies. Then the similar conclusion is presented for the discretized problem ${{\text{(OCP)}}_{\text{N}}}$.

Lemma 2.1. The triplet $({\boldsymbol y}, r, {\boldsymbol u})$ is the solution of control problem $({\rm OCP})$ if and only if there is a $(\boldsymbol y^*, r^*, \lambda) \in \boldsymbol Y \times Q \times \mathbb{R}$ such that

where

Proof. We denote $G:\boldsymbol u\rightarrow \boldsymbol y$ the operator to solve the state equation (2.4), $D(G(\boldsymbol u))$ the Gâteaux derivative of $G$ at $\boldsymbol u$. Furthermore, let $U = Z = K = L^2(\Omega)^2$, $C = G_{ad}$ respectively and $\boldsymbol u_0 = -\boldsymbol f$ satisfy the Slater type condition in Theorem 5.2 in ^[2], then there exists $\boldsymbol \mu \in L^2(\Omega)^2$ such that $(\boldsymbol u, \boldsymbol \mu)$ satisfying

where $(G(\boldsymbol u), \boldsymbol u)\in \boldsymbol Y\times \boldsymbol U$ satisfies the control problem (OCP).

In the next analysis, it can be derived that

where $\lambda$ satisfies (2.12). In fact, we can complete the proof by two cases: $\|\boldsymbol y\|_{0, \Omega} < d$ and $\| \boldsymbol y\|_{0, \Omega} = d$.

In the first case $\| \boldsymbol y\|_{0, \Omega} < d$, the following inequality holds by (2.13)

Thus we have for all $\boldsymbol \omega\in \boldsymbol U$ and $\|\boldsymbol \omega\|_{0, \Omega} = 1$

due to $(d-\|\boldsymbol y \|_{0, \Omega})\boldsymbol \omega+\boldsymbol y\in G_{ad}$. Similarly,

By (2.15) and (2.16)

In the second case $\| \boldsymbol y\|_{0, \Omega} = d$, it follows from (2.13) that

Then, we have

which implies that

and

Then, the identity (2.14) follows from (2.17)–(2.19).

It can be derived from (2.4) that for any $\boldsymbol v \in L^2(\Omega)^2$

By $\boldsymbol y-\boldsymbol y_0+\lambda\boldsymbol y\in L^2(\Omega)^2$, we can introduce the co-state equation as follows

Letting $\boldsymbol w = \boldsymbol y^*$ in (2.20) and $\boldsymbol q = \boldsymbol y'(\boldsymbol u)\boldsymbol v$ in (2.21), we have that for all $\boldsymbol v \in L^2(\Omega)^2$

Then (2.11)-(2.12) follows from (2.4), (2.13), (2.19), (2.21)-(2.22). Furthermore, lemma 2.1 can be proved as soon as the uniqueness of the solution is derived for (2.11). In fact, letting both $(\boldsymbol y_1, r_1, \boldsymbol u_1, \boldsymbol y_1^*, r^*_1, \lambda)$ and $(\boldsymbol y_2, r_2, \boldsymbol u_2, \boldsymbol y_2^*, r^*_2, \lambda)$ satisfy (2.11), then we have

It follows that

By (2.13)-(2.14), it holds that for all $\boldsymbol v\in G_{ad}$

It follows from (2.23) and (2.24) that

which implies $\boldsymbol u_1 = \boldsymbol u_2, \boldsymbol y_1 = \boldsymbol y_2$. Furthermore, it can be deduced that $r_1 = r_2, \boldsymbol y^*_1 = \boldsymbol y^*_2, r^*_1 = r^*_2, \lambda_1 = \lambda_2$. Then we can complete the proof of lemma 2.1 as we argue above.

Similarly, we can derive optimality conditions for the discretized control problem (2.10), and obtain the following result.

Lemma 2.2. The control problem ${{\text{(OCP)}}_{\text{N}}}$ has a unique solution, and the triple $({\boldsymbol y}_N, r_N, {\boldsymbol u}_N) \in{\boldsymbol Y}_{N}\times M_{N-2}\times{\boldsymbol U}_{N}$ is the solution of ${{\text{(OCP)}}_{\text{N}}}$ if and only if there is a $({\boldsymbol y}^*_N, r^*_N, \lambda_N) \in {\boldsymbol Y}_N\times M_{N-2}\times \mathbb{R}$ such that

where

Remark 2.1. (see, e.g., ^[10]). The optimal control $\boldsymbol u \in H^2(\Omega)^2$ if the initial data functions $\boldsymbol y_0, \boldsymbol f\in L^2(\Omega)^2$ by (2.11).

3.   A posteriori error estimator

In this section, the a posteriori error estimates are derived, and an error indicator is established for the spectral approximation of the control problem. We first recall two important spectral projection operators and more details can be found in ^[15].

Lemma 3.1. Let $P^0_{1, N}: H_0^1(\Omega)^2\rightarrow (Q^0_N)^2$ be the projection operator, satisfying for any ${\boldsymbol{w}}\in H^1_0(\Omega)^2$

If ${\boldsymbol{w}}\in (H^m(\Omega)\cap H^1_0(\Omega))^2$ with $m\geq 1$, then

Lemma 3.2. Define $P_N: L^2(\Omega)\rightarrow Q_N$ being the $L^2$ orthogonal projection operator, which satisfies for any $r\in L^2(\Omega)$

If $r\in H^m(\Omega)$ with $m\geq 0$, then

The following lemma is helpful for further analysis.

Lemma 3.3. Let $(\boldsymbol y_N, \boldsymbol u_N, \boldsymbol y^*_N, \lambda_N)$ be the solution of (2.25), then

Proof. By $(\boldsymbol y_N, \boldsymbol u_N, \boldsymbol y^*_N, \lambda_N)$ satisfying the optimality conditions (2.25), we have

associating with (2.25a) and (2.25e), it holds that

The discussion on estimating $\|\boldsymbol y^*_N\|_{1, \Omega}$ can be divided into two cases: $\|\boldsymbol y_N\|_{0, \Omega} < d$ and $\|\boldsymbol y_N\|_{0, \Omega} = d$.

The first case $\|\boldsymbol y_N\|_{0, \Omega} < d$, we have $\lambda_N = 0$. Letting $\boldsymbol q_N = \boldsymbol y^*_N$ in (2.25c), we have

The second case $\|\boldsymbol y_N\|_{0, \Omega} = d$, let $\boldsymbol q_N = \boldsymbol y^*_N-\frac{1}{d^2}(\boldsymbol y^*_N, \boldsymbol y_N)\boldsymbol y_N$ in (2.25c). Then it holds

It follows that

It follows from (3.3) and (3.4) that

Similarly, we can estimate $|\lambda_N|$ by two cases: $\|\boldsymbol y_N\|_{0, \Omega} < d$ and $\|\boldsymbol y_N\|_{0, \Omega} = d$.

If $\|\boldsymbol y_N\|_{0, \Omega} < d$, then $|\lambda_N|\leq C$.

If $\|\boldsymbol y_N\|_{0, \Omega} = d,$ letting $\boldsymbol q_N = \boldsymbol y_N$ in (2.25c), we derive

Therefore, it can be derived that

Then we can complete the proof by (3.2), (3.5), and (3.6).

3.1.   Upper error bound

In this subsection, we establish the a posteriori error estimator and prove that it provides an upper bound for the discretization errors. It is convenient to introduce an auxiliary system: find $(\boldsymbol y(\boldsymbol u_N), \boldsymbol r(\boldsymbol u_N), \boldsymbol y^*(\boldsymbol u_N), \boldsymbol r^*(\boldsymbol u_N))$ such that

By utilizing lemma 3.1–3.3, we can get the following lemmas.

Lemma 3.4. Let $(\boldsymbol y, r, \boldsymbol u, \boldsymbol y^*, r^*, \lambda)$ and $(\boldsymbol y_N, r_N, \boldsymbol u_N, \boldsymbol y^*_N, r^*_N, \lambda_N)$ be the solutions of (2.11) and (2.25) respectively. Then we have that

Proof. By (2.11) and (3.7), we have

The following analysis is divided into three cases.

If $\| \boldsymbol y\|_{0, \Omega} = d$, let $\boldsymbol q = \boldsymbol y$ in (3.9d). Then we have

Let $\boldsymbol w = \boldsymbol y^*-\boldsymbol y^*(\boldsymbol u_N)$ in (2.11a), we have

Thus it can be deduced that

If $\| \boldsymbol y_N\|_{0, \Omega} = d$, let $\boldsymbol q = \boldsymbol y(\boldsymbol u_N)$ in (3.9d) to obtain that

Let $\boldsymbol w = \boldsymbol y^*-\boldsymbol y^*(\boldsymbol u_N)$ in (3.7a) to derive that

It follows that

If $\|\boldsymbol y\|_{0, \Omega} < d$ and $\|\boldsymbol y_N\|_{0, \Omega} < d$, it follows from (2.12) and (2.26) that $\lambda = \lambda_N = 0$, which implies the identity

Letting $\boldsymbol w = \boldsymbol y-\boldsymbol y(\boldsymbol u_N)$ in (3.9a) and $\phi = r-r(\boldsymbol u_N)$ in (3.9b), we have that

which implies that

Then, we can derive from (2.11e), (2.25e), (3.10)–(3.13), and lemma 3.3 that

which implies the lemma 3.4.

Lemma 3.5. Let $(\boldsymbol y, r, \boldsymbol u, \boldsymbol y^*, r^*, \lambda)$ be the solution of (2.11), and $(\boldsymbol y_N, r_N, \boldsymbol u_N, \boldsymbol y^*_N, r^*_N, \lambda_N)$ be the solution of (2.25). Then we have that

and

Proof. By (2.25) and (3.7), we have

Noting that

It follows that

Similarly, we can derive from (2.25) and (3.7) that

and

Then, we can complete the proof by (3.16)–(3.19) and Cauchy's inequality with $\epsilon$.

The main result of this subsection can now be stated as follows.

Theorem 3.1. Let $(\boldsymbol y, r, \boldsymbol u, \boldsymbol y^*, r^*, \lambda)$ and $(\boldsymbol y_N, r_N, \boldsymbol u_N, \boldsymbol y^*_N, r^*_N, \lambda_N)$ be the solutions of (2.11) and (2.25) respectively. Then we have that

where the total approximation error $\|e\|$ is defined by

the estimator $\eta$ is given below

and $\theta$ is presented as follows

Proof. By (2.11e), (2.25e), and (3.9), we have

It is clear that by (2.12) and (2.26)

which implies that

It follows from (3.24) and (3.26) that

associating with lemma 3.1 and lemma 3.3, we have

It can be derived from (3.9), (3.13), and lemma 3.3 that

associating with (3.13), (3.27), and lemma 3.4, it can be derived that

Then the theorem follows from lemma 3.4, lemma 3.5, (3.13), and (3.28)-(3.29).

3.2.   Lower error bound

A lower bound for the error $\|e\|$ is established to investigate the sharpness of the indicator in this subsection. We first need some polynomial inverse estimates presented in the following lemma, which can be found in ^[16].

Lemma 3.6. Let $\alpha, \beta \in \mathbb{R}$ satisfy $-1 < \alpha < \beta$ and $\delta\in [0, 1]$. Then there exist constant $C_1, C_2 = C_2(\alpha, \beta)$, and $C_3 = C_3(\delta)$ such that for all $z_N\in Q_N$

where $\Phi_\Omega(x)$ is the distance function defined by

In the subsequent analysis, we establish the lower bound for discretization error $\|e\|$ by using lemma 3.6. In fact, letting $\alpha = 0, \beta = \gamma\in(\frac{1}{2}, 1]$ in (3.30b), we can derive that

Denote $\boldsymbol \vartheta = (\nu\Delta \boldsymbol y_N- \nabla r_N+\boldsymbol u_N+P_N\boldsymbol f)\Phi^\gamma_\Omega$, then we have

Furthermore, it follows from (3.30b) and (3.30c) that

Thus, we have

It follows from (3.31)–(3.33) that

It is clear that

Then, we have

Similarly, letting $\alpha = 0, \beta = \gamma\in(\frac{1}{2}, 1]$ in (3.30b), it holds that

Let $\boldsymbol \vartheta^* = (\nu\Delta \boldsymbol y^*_N+\nabla r^*_N+(1+\lambda_N)\boldsymbol y_N-P_N\boldsymbol y_0)\Phi^\gamma_\Omega$ to obtain that

Furthermore, it can be derived from (3.30b) and (3.30c) that

which implies the inequality

It follows from (3.36)–(3.38) that

Similar to $\eta_2$, we have

Then the estimation of lower error bound can be stated as following theorem.

Theorem 3.2. Let $(\boldsymbol y, r, \boldsymbol u, \boldsymbol y^*, r^*, \lambda)$ and $(\boldsymbol y_N, r_N, \boldsymbol u_N, \boldsymbol y^*_N, r^*_N, \lambda_N)$ be the solutions of (2.11) and (2.25) respectively. Then we have that

where $e$ and $\theta$ are defined in theorem 3.1.

Proof. The theorem follows from (3.34)-(3.35) and (3.39)-(3.40).

4.   Numerical experiment

In this section, we carry out a numerical example for the control problem (OCP) to investigate whether the indicator $\eta$ tends to zero at the same rate as the error $\|e\|$. We are interested in the following model: find $({\boldsymbol y}, r, {\boldsymbol u}) \in {\boldsymbol Y} \times Q \times {\boldsymbol U}$ such that

where $G_{ad} = \{\boldsymbol v \in \boldsymbol U \ \mid \ \parallel {\boldsymbol v}\parallel_{0, \Omega} \ \leq d\}$. The data of this example are as follows

and $d = \sqrt{6}\pi$. The exact solutions are given by

We solve the discrete system (2.25) via Arrow-Hurwicz algorithm (see, for example, ^[17] and ^[18]).

Arrow-Hurwicz Algorithm: we describe the main steps of the algorithm as follows.

● Step 1: Let $k = 0$, choose a step size $\rho > 0$, give the initial values $\lambda_N^0$ and $\boldsymbol u_N^0$.

● Step 2: Let $l = 0$ and $\boldsymbol u_N^{k, 0} = \boldsymbol u_N^k$.

● Step 3: Solve the state equations

and the co-state equations

● Step 4: Let $\boldsymbol u_N^{k, l+1} = \boldsymbol u_N^{k, l}-P_N(\boldsymbol y_N^{*k, l} +\alpha\boldsymbol u_N^{k, l})$. If

let $l = l+1$ and then we turn to Step 3.

● Step 5: $\lambda_N^{k+1} = max\{0, \lambda_N^k+\rho(\|\boldsymbol y_N^{k, l}\|_{0, \Omega}-d)\}.$

● Step 6: Stop if $\mid \lambda_N^{k+1}-\lambda_N^k\mid \ < {\rm Tol_\lambda}$ and output

else, let $\boldsymbol u_N^{k+1} = \boldsymbol u_N^{k, l+1}$, $k = k+1$, and turn to Step 2.

Denote ${\boldsymbol y}_N = (y_{N1}, y_{N2})$, ${\boldsymbol y}_N^* = (y^*_{N1}, y^*_{N2})$, and ${\boldsymbol u}_N = (u_{N1}, u_{N2})$, we have the following expressions by the property of (2.7)

The numerical results are presented in Table 1. The table shows that the errors decrease rapidly with a relatively small number of unknowns, which is important in a number of applications. We further plot the error $\|e\|$ and the error indicator $\eta$ versus the ploynomial degree $N$ in Figure 1, where the longitudinal axis is in logarithmic scale. It can be observed from the figure that the two curves are very close to each other, which implies that the indicator is nearly equivalent to the error $\|e\|$. Moreover, it seems that the error $\|e\|$ and the indicator $\eta$ decay with the rate $\frac{10^4}{10^{0.625N}}$, which shows that the proposed method for the control problem is very efficient and the spectral accuracy is achieved. Compared with the case of finite element method, it can be seen that the indicator developed in this work can be implemented more simply and can provide successful estimation for the errors with less computational load, which is helpful for developing the hp adaptive spectral element method for the optimal control problems.

5.   Conclusions

In this paper, the upper and lower bounds of approximation error are provided with the help of the a posteriori error indicator. The illustrative numerical experiment shows the performance of the error estimator. In our future work, we hope to extend these results to adaptive method in the hp spectral element framework, which will be compared with the adaptive hp finite element method for optimal control problems.

Acknowledgments

The authors are grateful to the referees for their useful comments and suggestions, which lead to improvements of the presentation. The authors were supported by the State Key Program of National Natural Science Foundation of China (11931003), the National Natural Science Foundation of China (Grant No. 11601466 and 41974133), the Nanhu Scholars Program for Young Scholars of XYNU, the Postgraduate Education Reform and Quality Improvement Project of Henan Province (Grant No. YJS2022ZX34), the National Natural Science Foundation of Henan (Grant No. 202300410343), the Training Plan of Young Key Teachers in Universities of Henan (Grant No. 2018GGJS096), and the Scientific Research Fund Project for Young Scholars of XYNU (Grant No. 2020-QN-050).

Conflict of interest

The authors declare there is no conflicts of interest.