﻿ Iterative Solution of Mesh Constrained Optimal Control Problems with Two-Level Mesh Approximations of Parabolic State Equation

Journal of Applied Mathematics and Physics
Vol.06 No.01(2018), Article ID:81585,11 pages
10.4236/jamp.2018.61007

Iterative Solution of Mesh Constrained Optimal Control Problems with Two-Level Mesh Approximations of Parabolic State Equation

A. Lapin1,2, E. Laitinen3

1Institute of Computational Mathematics & Information Technology, Kazan Federal University, Kazan, Russia

2Tianjin University of Finance and Economics, Tianjin, China

3Department of Mathematical Sciences, University of Oulu, Oulu, Finland

Received: November 21, 2017; Accepted: January 2, 2018; Published: January 5, 2018

ABSTRACT

We consider a linear-quadratical optimal control problem of a system governed by parabolic equation with distributed in right-hand side control and control and state constraints. We construct a mesh approximation of this problem using different two-level approximations of the state equation, ADI and fractional steps approximations in time among others. Iterative solution methods are investigated for all constructed approximations of the optimal control problem. Their implementation can be carried out in parallel manner.

Keywords:

Parabolic Optimal Control, State Constraints, Finite Difference Method, Constrained Saddle Point Problem, Iterative Method

1. Introduction

Optimal control of time-dependent production processes plays an important role in many real world applications. State constraints in optimal control of systems governed by partial differential equations have to be often included in the mathematical models. For instance, in continuous casting process a need to prevent the cracks in a slab and the solidification at a wrong place leads to the bounds on the temperature variable. Similar demands arise in the processes of crystal growth and cooling of glass melts (see articles [1] [2] [3] [4] [5] and bibliography therein). The introduction of pointwise state constraints yields adjoint variables and multipliers which only admit low regularity complicating both theoretical analysis and the constructing appropriate numerical methods.

There are few results in the area of numerical solution methods for the constrained parabolic optimal control problems. One of the approaches to overcome the difficulty connected with low regularity of the solutions is using Lavrentiev regularization. This approach has been used in [6] [7]. In [8] [9] [10] [11] a priori error estimates for space-time discretizations of linear-quadratic parabolic optimal control problems have been obtained for problems.

The implementations of the iterative methods for parabolic optimal control problems include the solution of the parabolic equation and corresponding adjoint parabolic equation at each iteration and this is the most time consuming part of the algorithms if applying the implicit (backward Euler) approximation of parabolic state equation. On the other hand, easily implementable explicit (forward Euler) approximation of a parabolic equation with a constant step in time requests extremely restrictive constraint for this step. Using explicit approximations of the parabolic equations with special series of non-uniform time steps allows partially avoid this deficiency. Such kind approximation is well-known for the differential equations [12] and they demonstrate the advantage in time of calculations in relation to the implicit schemes. Similar approximation have been used for the continuous casting problem [13] and recently applied to a parabolic state constrained problem [14]. One more approach to construct effective algorithm for parabolic optimal control problem is to use parareal approximation of the state equation [15] [16].

In this article we continue the investigations of [14] [16] [17] [18] [19] [20] on the iterative solution methods for the constrained saddle point problems with applications to optimal control problems. In the cited papers parabolic optimal control problems have been solved by using either backward or forward Euler approximating schemes for the state equation.

The main purpose of this article is to generalize these results for the case when any two-level scheme is used for the approximation of the parabolic state equation, including different splitting (locally one-dimensional) schemes.

To simplify the exposition we restrict ourselves to consider a problem in unit square with distributed control and observation, and to use finite difference schemes to approximate the state equation, while the most of the results can be extended for the case of lowest order finite element method, for a control in Neumann boundary condition, for final observation etc.

2. Formulation and Approximations of the Problem

Consider homogeneous Dirichlet initial-boundary value problem

$\frac{\partial y}{\partial t}-\Delta y=u\text{\hspace{0.17em}}\text{ }\text{in}\text{ }\text{\hspace{0.17em}}{Q}_{T};\text{\hspace{0.17em}}y=0\text{\hspace{0.17em}}\text{ }\text{on}\text{ }\text{\hspace{0.17em}}\Sigma =\partial \Omega ×\left(0,T\right];\text{\hspace{0.17em}}y\left(x,0\right)=0$ (1)

in the cylinder ${Q}_{T}=\Omega ×\left(0,T\right]$ , $\Omega {=\left(0,1\right)}^{2}$ , with lateral surface $\Sigma =\partial \Omega ×\left(0,T\right]$ . We call $y\left(x,t\right)$ and $u\left(x,t\right)$ as state and control functions. It is well-known that for any $u\in {L}_{2}\left({Q}_{T}\right)$ there exists a unique weak solution of problem (1) such that $y\in W\left(0,T\right)={L}^{\infty }\left(0,T;{H}_{0}^{1}\left(\Omega \right)\right)\cap {H}^{1}\left(0,T,{L}^{2}\left(\Omega \right)\right)$ , and the following inequality takes place ([21] [22] [23]):

$\underset{0\le t\le T}{sup}||y\left(t\right)|{|}_{{H}_{0}^{1}\left(\Omega \right)\right)}+||\frac{\partial y}{\partial t}|{|}_{{L}_{2}\left({Q}_{T}\right)}\le c||u|{|}_{{L}_{2}\left({Q}_{T}\right)}.$ (2)

Define objective function

$J\left(y,u\right)=\frac{1}{2}\underset{{Q}_{T}}{\int }{\left(y\left(x,t\right)-{y}_{d}\left(x,t\right)\right)}^{2}dxdt+\frac{1}{2}\underset{{Q}_{T}}{\int }{u}^{2}dxdt$

with a given functions ${y}_{d}\in {L}_{2}\left({Q}_{T}\right)$ , and the sets of constraints for state and control:

$\begin{array}{l}{U}_{ad}=\left\{u\in {L}_{2}\left({Q}_{T}\right):\text{ }|u\left(x,t\right)|\text{\hspace{0.17em}}\le {u}_{max}\text{\hspace{0.17em}}\text{ }a.e.\text{ }\text{\hspace{0.17em}}\left(x,t\right)\in {Q}_{T}\right\},\\ {Y}_{ad}=\left\{y\in W\left(0,T\right):\text{ }{y}_{min}\le y\left(x,t\right)\le {y}_{max}\text{\hspace{0.17em}}\text{ }a.e.\text{\hspace{0.17em}}\text{ }\text{\hspace{0.17em}}\left(x,t\right)\in {Q}_{T}\right\},\end{array}$

Above ${u}_{max}>0$ and $-\infty \le {y}_{min}<{y}_{max}\le \infty$ . We solve the following optimal control problem:

$\begin{array}{c}\underset{\left(y,u\right)\in K}{min}J\left(y,u\right),\\ K=\left\{\left(y,u\right)\in {Y}_{ad}×{U}_{ad}:\text{\hspace{0.17em}}\text{ }\text{state}\text{\hspace{0.17em}}\text{equation}\text{\hspace{0.17em}}\left(1\right)\text{\hspace{0.17em}}\text{is}\text{\hspace{0.17em}}\text{satisfied}\text{ }\right\}.\end{array}$ (3)

Problem (3) has a unique solution (cf., e.g. [14]).

We construct the finite-difference approximations of problem (3) using a uniform in x and t mesh ${\omega }_{x}×{\omega }_{t}$ in ${\stackrel{¯}{Q}}_{T}$ .

Let ${\omega }_{x}$ be the uniform mesh of the meshsize $h$ on $\stackrel{¯}{\Omega }$ , $\partial {\omega }_{x}={\omega }_{x}\cap \partial \Omega$ be the set of the boundary nodes while ${\omega }_{x}^{0}={\omega }_{x}\\partial {\omega }_{x}$ . By ${V}_{h}^{0}$ we denote the space of mesh functions which are defined on ${\omega }_{x}$ and vanish in the boundary nodes $\partial {\omega }_{x}$ , let $\text{dim}{V}_{h}^{0}={N}_{x}$ . The mesh on the time segment we denote by ${\omega }_{t}=\left\{{t}_{j}=j\tau ,\text{ }j=0,1,\dots ,{N}_{t};\text{ }{N}_{t}\tau =T\right\}$ .

We will use the notations $y,u,...$ for the mesh functions from ${V}_{h}^{0}$ and for the vectors of their nodal values as well. We also don’t distinguish the linear operators from ${V}_{h}^{0}$ to ${V}_{h}^{0}$ and corresponding them ${N}_{x}×{N}_{x}$ matrices acting on the vectors of nodal values of mesh functions from ${V}_{h}^{0}$ .

By ${y}_{j}=y\left(x,{t}_{j}\right)\in {R}^{{N}_{x}}$ we denote the values on a time level ${t}_{j}=j\tau \in {\omega }_{t}$ of a mesh function of x and t, and by $||\cdot |{|}_{x}$ -euclidian norm in the space ${R}^{{N}_{x}}$ . We use the following notations: $N={N}_{t}{N}_{x}$ , $\left(.,.\right)$ and $||\cdot ||$ are the inner product and the norm in ${R}^{N}$ . Let ${E}_{x}$ be ${N}_{x}×{N}_{x}$ unit matrix while $E$ be the unit matrix in ${R}^{N×N}$ .

We denote by ${\Lambda }_{1}:\text{\hspace{0.17em}}{V}_{h}^{0}\to {V}_{h}^{0}$ the linear operator such that

${\Lambda }_{1}y\left(x\right)={h}^{-2}\left(2y\left({x}_{1},{x}_{2}\right)-y\left({x}_{1}+h,{x}_{2}\right)-y\left({x}_{1}-h,{x}_{2}\right)\text{\hspace{0.17em}}\text{ }\text{for}\text{ }\text{\hspace{0.17em}}x\in {\omega }_{x}^{0}$

and similarly for ${\Lambda }_{2}y\left(x\right).$ Obviously, $\Lambda ={\Lambda }_{1}+{\Lambda }_{2}$ is the mesh Laplasian with homogeneous Dirichlet boundary conditions. Recall that we use the same notations for corresponding matrices. Matrices ${\Lambda }_{i}$ are symmetric, commute and positive definite with spectrum in a segment $\left[{\xi }_{0},{\xi }_{1}\right]$ , where ${\xi }_{1}$ is of order ${h}^{-2}$ , while ${\xi }_{0}>0$ is bounded below by a constant which doesn’t depend on $h$ .

Let us introduce ${N}_{x}×{N}_{x}$ matrices U, D and R, where U is a symmetric and positive matrix, and approximate state problem (1) by two-level finite difference scheme

$U\frac{{y}_{j}-{y}_{j-1}}{\tau }+D{y}_{j-1}=R{u}_{j},\text{\hspace{0.17em}}j=1,2,\dots ,{N}_{t},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{y}_{0}=0.$ (4)

Below we give several examples of two-level finite difference scheme (4).

1) Finite difference scheme with weights:

$\frac{{y}_{j}-{y}_{j-1}}{\tau }+\Lambda \left(\sigma {y}_{j}+\left(1-\sigma \right){y}_{j-1}\right)={u}_{j},\text{\hspace{0.17em}}j=1,2,\dots ,{N}_{t},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{y}_{0}=0$ (5)

with $\sigma \in \left[0,1\right]$ can be written in the form (4) with $U={E}_{x}+\sigma \tau \Lambda$ , $D=\Lambda$ and $R={E}_{x}$ . Recall that (5) contains forward Euler ( $\sigma =0$ ), backward Euler ( $\sigma =1$ ) and Crank-Nicolson ( $\sigma =1/2$ ) schemes.

Scheme (5) is unconditionally stable for $\sigma \ge 1/2$ and stable in the case $0\le \sigma <1/2$ if $\tau \le 2\left({\xi }_{max}\left(1-2\sigma {\right)\right)}^{-1}$ , where ${\xi }_{max}=2{\xi }_{1}$ is the maximal eigenvalue of $\Lambda$ . The stability estimate is:

$\underset{j=1}{\overset{{N}_{t}}{\sum }}||{y}_{j}|{|}_{x}^{2}\text{\hspace{0.17em}}\le \underset{j=1}{\overset{{N}_{t}}{\sum }}||{u}_{j}|{|}_{x}^{2}.$ (6)

2) Two-level scheme with factorized preconditioner:

$\begin{array}{c}\left({E}_{x}+\sigma \tau {\Lambda }_{1}\right)\left({E}_{x}+\sigma \tau {\Lambda }_{2}\right)\frac{{y}_{j}-{y}_{j-1}}{\tau }+\Lambda {y}_{j-1}={u}_{j},\\ j=1,2,\dots ,{N}_{t},\text{\hspace{0.17em}}{y}_{0}=0.\end{array}$ (7)

Mesh scheme (7) is unconditionally stable for $\sigma \ge 0.5$ with stability estimate (6).

3) Fractional steps scheme:

$\begin{array}{c}0.5\frac{{z}_{1}-{y}_{j-1}}{\tau }+{\Lambda }_{1}{z}_{1}=0.5{u}_{j},\\ 0.5\frac{{z}_{2}-{y}_{j-1}}{\tau }+{\Lambda }_{2}{z}_{2}=0.5{u}_{j},\\ {y}_{j}=0.5{z}_{1}+0.5{z}_{2},\text{\hspace{0.17em}}j=1,2,\dots ,{N}_{t},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{y}_{0}=0.\end{array}$ (8)

System (8) can be written in form (4) with

$U=\left({E}_{x}+2\tau {\Lambda }_{1}\right)\left({E}_{x}+2\tau {\Lambda }_{2}\right),\text{\hspace{0.17em}}D=\Lambda +4\tau {\Lambda }_{1}{\Lambda }_{2},\text{\hspace{0.17em}}R={E}_{x}+\tau \Lambda .$

Scheme (8) is unconditionally stable with stability estimate (6).

Let us further use notation ${y}_{d}:\text{ }{\omega }_{t}\to {V}_{h}^{0}$ for an ${V}_{h}^{0}$ -approximation of the function ${y}_{d}$ . Define a mesh goal function and the sets of the constraints:

${J}_{h}\left(y,u\right)=\frac{1}{2}\underset{j=1}{\overset{{N}_{t}}{\sum }}\tau ||{y}_{j}-{y}_{dj}|{|}_{x}^{2}+\frac{1}{2}\underset{j=1}{\overset{{N}_{t}}{\sum }}\text{ }\text{ }\tau ||{u}_{j}|{|}_{x}^{2},$ (9)

$\begin{array}{c}{U}_{ad}^{h}=\left\{u\in {R}^{N}:|{u}_{i}^{j}|\text{\hspace{0.17em}}\le \stackrel{¯}{u}\text{\hspace{0.17em}}\forall x\in {\omega }_{x},\forall t\in {\omega }_{t}\right\},\\ {Y}^{h}=\left\{y\in {R}^{N}:{y}_{min}\le {y}_{j}\le {y}_{max},\text{\hspace{0.17em}}\forall x\in {\omega }_{x},\forall t\in {\omega }_{t}\right\}.\end{array}$

The mesh optimal control problem reads as follows:

$\begin{array}{c}\underset{\left(y,u\right)\in {K}_{h}}{min}{J}_{h}\left(y,u\right),\\ {K}_{h}=\left\{\left(y,u\right)\in {Y}_{ad}^{h}×{U}_{ad}^{h}:\text{\hspace{0.17em}}\text{ }\text{\hspace{0.17em}}\left(4\right)\text{\hspace{0.17em}}\text{is}\text{\hspace{0.17em}}\text{satisfied}\text{ }\right\}.\end{array}$ (10)

Problem (10) has a unique solution because the set ${K}_{h}$ is a convex compact and the quadratical function ${J}_{h}\left(y,u\right)$ is continuous and strictly convex on ${K}_{h}$ .

3. Saddle Point Problem and Iterative Solution Method

Let us rewrite problem (10) in a “vector-matrix” form. Let the matrices $L,Q\in {R}^{N×N}$ be defined by the equalities:

$\begin{array}{c}{\left(Ly\right)}_{j}=\left\{U\frac{{y}_{1}}{\tau }\text{\hspace{0.17em}}\text{ }\text{for}\text{ }\text{\hspace{0.17em}}j=1,\text{\hspace{0.17em}}\text{\hspace{0.17em}}U\frac{{y}_{j}-{y}_{j-1}}{\tau }+D{y}_{j-1}\text{\hspace{0.17em}}\text{ }\text{for}\text{ }\text{\hspace{0.17em}}j=2,\dots ,{N}_{t}\right\},\\ {\left(Qu\right)}_{j}=\left\{R{u}_{j}\text{\hspace{0.17em}}\text{ }\text{for}\text{ }\text{\hspace{0.17em}}j=1,\dots ,{N}_{t}\right\}.\end{array}$

Denote by $\psi$ and $\phi$ the indicator functions of the sets ${Y}_{0}^{h}$ and ${U}_{ad}^{h}$ , respectively. We obtain the following algebraic form of mesh optimal control problem (10):

$\underset{Ly=u}{min}\left\{I\left(y,u\right)=\frac{1}{2}||y-{y}_{d}|{|}^{2}+\frac{1}{2}||u|{|}^{2}+\psi \left(y\right)+\phi \left(u\right)\right\}.$ (11)

Below we suppose that for a two-level finite difference scheme (4) the following assumption is valid:

$\exists \text{\hspace{0.17em}}\text{const}\text{\hspace{0.17em}}{C}_{stab}:\text{\hspace{0.17em}}||{L}^{-1}Q||\text{\hspace{0.17em}}\le {C}_{stab}.$ (12)

The stability estimate (6) approve that for all cited above particular cases of (4) assumption (12) is true (with constant ${C}_{stab}=1$ ).

Remark 1 More well-known finite difference schemes can be written in form (4) and satisfy assumption (12): different kinds of ADI schemes proposed in [24], [25] and “sequential” variant of fractional steps scheme [26] [27] etc. Also a more general variant of scheme (8) (see [28]) with positive weights ${\alpha }_{1},{\alpha }_{2}$ , such that ${\alpha }_{1}+{\alpha }_{2}=1$ instead of ${\alpha }_{i}=1/2$ as in (8) can be considered.

We construct Lagrange function for problem (11):

$L\left(y,u,\lambda \right)=I\left(y,u\right)+\psi \left(y\right)+\phi \left(u\right)+\left(\lambda ,Ly-Qu\right).$

A saddle point of this Lagrangian satisfies the following system:

$\left(\begin{array}{ccc}E& 0& {L}^{T}\\ 0& E& -{Q}^{T}\\ L& -Q& 0\end{array}\right)\left(\begin{array}{c}y\\ u\\ \lambda \end{array}\right)+\left(\begin{array}{c}\partial \psi \left(y\right)\\ \partial \phi \left(u\right)\\ 0\end{array}\right)∍\left(\begin{array}{c}{y}_{d}\\ 0\\ 0\end{array}\right),$ (13)

where $\partial \psi \left(y\right)$ and $\partial \phi \left(u\right)$ are the subdifferentials of the corresponding functions.

With the notations $z=\left(y,u{\right)}^{T}$ , $f=\left({y}_{d}{,0,0\right)}^{T}$ , $\Psi \left(z\right)=\psi \left(y\right)+\phi \left(u\right)$ and

$A=\left(\begin{array}{cc}E& 0\\ 0& E\end{array}\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}B=\left(\begin{array}{cc}L& -Q\end{array}\right)$

(11) becomes a particular case of minimization problem

$\underset{Bz=0}{min}\left\{\frac{1}{2}\left(Az,z\right)-\left(f,z\right)+\Psi \left(z\right)\right\},$ (14)

while (13)-a particular case of saddle point problem

$\left(\begin{array}{cc}A& {B}^{T}\\ B& 0\end{array}\right)\left(\begin{array}{c}z\\ \eta \end{array}\right)+\left(\begin{array}{c}\partial \Psi \left(z\right)\\ 0\end{array}\right)∍\left(\begin{array}{c}f\\ 0\end{array}\right).$ (15)

We will use the following results from the article [17]:

Proposition 1 Let

$\text{matrix}\text{\hspace{0.17em}}A\in {R}^{m×m}\text{\hspace{0.17em}}\text{issymmetricandpositivedefinite;}$ (16)

$\text{matrix}\text{\hspace{0.17em}}B\in {R}^{s×m}\text{\hspace{0.17em}}\text{hasafullcolumnrank;}$ (17)

$\Psi :{R}^{m}\to \stackrel{¯}{R}\text{\hspace{0.17em}}\text{isaconvex},\text{properandlower}\text{\hspace{0.17em}}\text{semicontinuousfunction;}$ (18)

$\left\{z\in {R}^{m}:Bz=0\right\}\cap \text{int}\text{ }\text{dom}\text{ }\Psi \ne \varnothing .$ (19)

Then

1) Problem (15) has a non-empty set of the solutions $X=\left\{\left(z,\eta \right)\right\}$ , where z is unique solution of (14).

2) Uzawa-type iterative method

$\begin{array}{c}A{z}^{k+1}+\partial \Psi \left({z}^{k+1}\right)∍{B}^{T}{\eta }^{k}+f,\\ \frac{1}{\rho }D\left({\eta }^{k+1}-{\eta }^{k}\right)+B{z}^{k+1}=0,\text{\hspace{0.17em}}\rho >0\end{array}$ (20)

with a symmetric preconditioner $D$ satisfying the inequality

$\left(D\eta ,\eta \right)\frac{\left(1+\epsilon \right)\rho }{2}\left({A}^{-1}{B}^{T}\eta ,{B}^{T}\eta \right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\forall \eta \in {R}^{s},\text{\hspace{0.17em}}\epsilon >0,$ (21)

converges for any initial guess ${\eta }^{0}$ : $\left({z}^{k},{\eta }^{k}\right)\to \left({z}^{*},{\eta }^{*}\right)\in X$ for $k\to \infty$ .

Theorem 1 Problem (13) has a solution $\left(y,u,\lambda \right)$ with unique pair $\left(y,u\right)$ , which coincides with the solution of problem (11).

Proof. Obviously, all assumptions (16)-(18) and (19) of proposition 1 are satisfied for problem (13). In particular, vector $\left({y}_{0},{u}_{0}\right)=\left(0,0\right)$ satisfies assumption (19).

Below we construct for problem (13) the easily implementable preconditioner D which is spectrally equivalent to $B{A}^{-1}{B}^{T}$ with constants independent on $\tau$ and h. As this preconditioner we take $D=L{L}^{T}$ . Then method (20) for problem (13) with this preconditioner reads as follows:

$\begin{array}{c}{y}^{k+1}+\partial \psi \left({y}^{k+1}\right)∍{y}_{d}-{L}^{T}{\lambda }^{k},\\ {u}^{k+1}+\partial \phi \left({u}^{k+1}\right)∍{Q}^{T}{\lambda }^{k},\\ L{L}^{T}\frac{{\lambda }^{k+1}-{\lambda }^{k}}{\rho }=L{y}^{k+1}-Q{u}^{k+1}.\end{array}$ (22)

Theorem 2 Method (22) converges for any initial guess ${\lambda }^{0}$ if

$0<\rho <\frac{2}{1+{C}_{stab}^{2}}.$ (23)

Proof. The matrix $B{A}^{-1}{B}^{T}$ is spectrally equivalent to $D=L{L}^{T}$ , namely,

$L{L}^{T}B{A}^{-1}{B}^{T}\left(1+{C}_{stab}^{2}\right)L{L}^{T},$ (24)

where constant ${C}_{stab}$ is defined in (12). In fact, by direct calculation we find $B{A}^{-1}{B}^{T}=L{L}^{T}+Q{Q}^{T}\ge L{L}^{T}$ . Further, since $||{L}^{-1}Q||\text{\hspace{0.17em}}\le {C}_{stab}$ , then $P{Q}^{T}{L}^{-T}P\le {C}_{stab}$ . This inequality is equivalent to

$Q{Q}^{T}\le {C}_{stab}^{2}L{L}^{T},$ (25)

whence the result.

In virtue of inequality (24) the convergence condition (21) of proposition 1 is true for the parameter $\rho$ from (23).

In the case of optimal control problem without state constraints we can estimate a rate of convergence for iterative method (22) and give an optimal iterative parameter $\rho$ . Namely, the following statement holds:

Theorem 3 Let $\partial \psi =0$ . Then there exists a unique solution $\left(y,u,\lambda \right)$ of saddle point problem (13), and for theoretically optimal iterative parameter

${\rho }_{0}=\frac{1}{{C}_{stab}}$ the following estimate for the rate of convergence of method (22) is

valid:

$||{L}^{T}\left({\lambda }^{k+1}-\lambda \right)||\le {\left(1-\frac{1}{{C}_{stab}}\right)}^{1/2}||{L}^{T}\left({\lambda }^{k}-\lambda \right)||,\text{\hspace{0.17em}}\text{\hspace{0.17em}}k=0,1,\dots$ (26)

Proof. The uniqueness of $\left(y,u\right)$ is proved in theorem 1. The uniqueness of $\lambda$ in the case $\partial \psi =0$ follows from the equation $y+{L}^{T}\lambda ={y}_{d}$ .

Vector $\lambda$ is the solution of the equation

$L{L}^{T}\lambda +Q\circ {\left(E+\partial \phi \right)}^{-1}\circ {Q}^{T}\lambda =L{y}_{d},$ (27)

while iterative method (22) can be written in the form

$L{L}^{T}\frac{{\lambda }^{k+1}-{\lambda }^{k}}{\rho }+L{L}^{T}{\lambda }^{k}+Q\circ {\left(E+\partial \phi \right)}^{-1}\circ {Q}^{T}{\lambda }^{k}=L{y}_{d}.$ (28)

It is well-known (cf., e.g. [29]) that the operator ${\left(E+\partial \phi \right)}^{-1}$ is co-coercive:

$\left(\left(E+\partial \phi {\right)}^{-1}\lambda -{\left(E+\partial \phi \right)}^{-1}\mu ,\lambda -\mu \right)\ge ||{\left(E+\partial \phi \right)}^{-1}\lambda -{\left(E+\partial \phi \right)}^{-1}\mu |{|}^{2}$

Because of this and (25) the operator $P=E+\left({L}^{-1}Q\right)\circ {\left(E+\partial \phi \right)}^{-1}\circ \left({Q}^{T}{L}^{-T}\right)$ satisfies the following properties (strong monotonicity and Lipshitz-continuity):

$\left(P\left(\lambda \right)-P\left(\mu \right),\lambda -\mu \right)\ge ||\lambda -\mu |{|}^{2},$

$\left(P\left(\lambda \right)-P\left(\mu \right),\eta \right)\le {C}_{stab}{\left(P\left(\lambda \right)-P\left(\mu \right),\lambda -\mu \right)}^{1/2}||\eta ||$

The rest of the prove is quite standard (cf. [18]). Namely, with the notations ${\eta }^{k}={L}^{T}{\lambda }^{k}$ , $\eta ={L}^{T}\lambda$ and ${z}^{k}={\eta }^{k}-\eta$ we have the equation

$\frac{{z}^{k+1}-{z}^{k}}{\rho }+P\left({\eta }^{k}\right)-P\left(\eta \right)=0.$

We multiply this equation by $2\rho {z}^{k+1}$ and obtain the equality

$||{z}^{k+1}|{|}^{2}-||{z}^{k}|{|}^{2}+||{z}^{k+1}-{z}^{k}|{|}^{2}+2\rho \left(P\left({\eta }^{k}\right)-P\left(\eta \right),{z}^{k+1}\right)=0,$

Due to the properties of $P$ the following estimate holds:

$\begin{array}{l}2\rho \left(P\left({\eta }^{k}\right)-P\left(\eta \right),{z}^{k+1}\right)=2\rho \left(P\left({\eta }^{k}\right)-P\left(\eta \right),{z}^{k}\right)+2\rho \left(P\left({\eta }^{k}\right)-P\left(\eta \right),{z}^{k+1}-{z}^{k}\right)\\ \ge \left(2\rho -{\rho }^{2}{C}_{stab}^{2}\right)\left(P\left({\eta }^{k}\right)-P\left(\eta \right),{z}^{k}\right)-||{z}^{k+1}-{z}^{k}|{|}^{2}.\end{array}$

Substituting this estimate in the previous equality we get

$||{z}^{k+1}|{|}^{2}\le \left(1-\rho \left(2-\rho {C}_{stab}^{2}\right)\right)||{z}^{k}|{|}^{2},$

whence $||{z}^{k}||\to 0$ if $0<\rho <\frac{2}{{C}_{stab}^{2}}$ and rate of convergence (26) is true for optimal parameter ${\rho }_{0}=\frac{1}{{C}_{stab}^{2}}$ .

Remark 2 Due to the equalities $y={y}_{d}-{L}^{T}\lambda$ , ${y}^{k}={y}_{d}-{L}^{T}{\lambda }^{k-1}$ , we have the following estimate for the rate of convergence of ${y}^{k}$ :

$\begin{array}{l}||{y}^{k}-y||\le {\left(1}^{-}||{y}^{k-1}-y||\le {\left(1}^{-}||{y}^{1}-y||\\ =\left(1-\frac{1}{{C}_{stab}}{\right)}^{\left(k-1\right)/2}||{L}^{T}\left({\lambda }^{0}-\lambda \right)||\text{\hspace{0.17em}}\text{\hspace{0.17em}}k=1,2,\dots \end{array}$

For the sequence of control vectors $\left\{{u}^{k}\right\}$ the estimate is as follows:

$\begin{array}{l}||{u}^{k}-u||=||{\left(}^{E}\left({\lambda }^{k-1}\right)-{\left(}^{E}\left(\lambda \right)||\le ||{L}^{-T}||||{L}^{T}\left({\lambda }^{k-1}-\lambda \right)||\\ \le ||{L}^{-T}||{\left(1}^{-}||{L}^{T}\left({\lambda }^{0}-\lambda \right)||\text{\hspace{0.17em}}\text{\hspace{0.17em}}k=1,2,\dots \end{array}$

Remark 3 For the state constrained problems there are no estimates for rate of convergence for iterative method (22). In this case instead of (27) we have the equation

$L\circ \left(E+\partial \psi \right)\circ {L}^{T}\lambda +Q\circ {\left(E+\partial \phi \right)}^{-1}\circ {Q}^{T}\lambda =L{y}_{d},$

which operator is only co-coercive. The convergence of the iterative methods for such kind of equations have been investigated in [18].

On the other hand, numerous calculations show that in this general case the choice of the iterative parameter $\rho ;{\rho }_{0}$ as in theorem 3 is practical and seems to be close to optimal one.

When implementing method (22) one has to solve the inclusions with respect ${y}^{k+1}$ and ${u}^{k+1}$ , and a system of linear equations with matrix $L{L}^{T}$ . Concerning solving the inclusions we underline that the matrices and the operators in them have diagonal form, so, their solving reduces to easy pointwise projection.

In turn, solving equation with the matrix $L{L}^{T}$ is equivalent to solving direct (with $L$ ) and adjoint (with ${L}^{T}$ ) parabolic mesh schemes. In the case of mesh schemes with factorized preconditioner (7) or fractional steps scheme (8) approximating state equation their solving reduces to solving set of non-coupled “one-dimensional” mesh problems-systems of linear algebraic equations with tridiagonal matrices ${E}_{x}+\text{const}\tau {\Lambda }_{i}$ , $i=1,2$ . Obviously, these systems can be solved by a well-known direct method and parallel.

4. Variants and Generalizations

The effectiveness of the implementation of iterative method (22) is based on two main properties:

・ Preconditioner has factorized form $D=L{L}^{T}$ and is spectrally equivalent to “main” matrix of the problem;

・ Equations with the matrices $L$ and ${L}^{T}$ as in (7) and (8) can be easily implementable.

The first property is ensured by the inequality (12): $||{L}^{-1}Q||\le {C}_{stab}$ . Just this inequality allows us to prove spectral equivalency of the matrix $B{A}^{-1}{B}^{T}=L{L}^{T}+Q{Q}^{T}$ and $D=L{L}^{T}$ with the constants independent on mesh parameters $\tau$ and $h$ . In turn, this inequality is nothing but a stability estimate for a corresponding two-level approximation of a parabolic state equation. Numerous classes of stable two-level finite difference schemes for the parabolic equations can be found in [27].

The second property-easy solution of the equations with matrices $L$ , ${L}^{T}$ in (7) and (8)-is the consequence of their “local one-dimensional” structure. This imposes several limitations to the domains, boundary conditions and using orthogonal meshes. Nevertheless, a lot of different mesh schemes with factorized preconditioner of the form (7), satisfying stability property (12) is known (cf., e.g. [24] [27] and the bibliography therein).

The results of this paper can be extended to the parabolic optimal control problems with other state and control constraints, such as, for example, in [14] and [30].

Acknowledgements

The research was supported by the RFBR grant No. 16-01-00408 and by grants No. 296348 and No. 296928 from Academy of Finland.

Cite this paper

Lapin, A. and Laitinen, E. (2018) Iterative Solution of Mesh Constrained Optimal Control Problems with Two-Level Mesh Approximations of Parabolic State Equation. Journal of Applied Mathematics and Physics, 6, 58-68. https://doi.org/10.4236/jamp.2018.61007

References

1. 1. Dautov, R., Kadyrov, R., Laitinen, E., Lapin, A., Pieskä, J. and Toivonen V. (2003) On 3D Dynamic Control of Secondary Cooling in Continuous Casting Process. Lobachevskii J. Math., 13, 3-13.

2. 2. Gunzburger, M., Ozugurlu, E., Turner, J. and Zhang, H. (2002) Controlling transport Phenomena in the Czochralski Crystal Growth Process. J. Cryst. Growth, 234, 47-62. https://doi.org/10.1016/S0022-0248(01)01635-9

3. 3. Clever, D. and Lang, J. (2008) Optimal Control of Radiative Heat Transfer in Glass Cooling with Restrictions on the Temperature Gradient. Preprint SPP1253-20-01.

4. 4. Pinnau, R. (2007) Analysis of Optimal Boundary Control for Radiative Heat Transfer Modelled by the SP1-System. CMS, 5, 951-969.

5. 5. Hinze, M. and Ziegenbalg, S. (2007) Optimal Control of the Free Boundary in a Two-Phase Stefan Problem. J. Comput. Phys., 223, 657-684. https://doi.org/10.1016/j.jcp.2006.09.030

6. 6. Neitzel, I. and Troltzsch, F. (2009) On Regularization Methods for the Numerical Solution of Parabolic Control Problems with Pointwise State Constraint. ESAIM: Control, Optimisation and Calculus of Variations, 15, 426-452. https://doi.org/10.1051/cocv:2008038

7. 7. Neitzel, I. and Troltzsch, F. (2008) On Convergence of Regularization Methods for Nonlinear Parabolic Optimal Control Problems with Control and State Constraints. Control and Cybernetics, 37, 1013-1043.

8. 8. Meidner, D. and Vexler, B. (2008) A Priori Error Estimates for Space-Time Finite Element Approximation of Parabolic Optimal Control Problems. Part II: Problems with Control Constraints. SIAM J. Control Optim., 47, 1301-1329. https://doi.org/10.1137/070694028

9. 9. Deckelnick, K. and Hinze, M. (2011) Variational Discretization of Parabolic Control Problems in the Presence of Pointwise State Constraints. J. Comp. Math., 29, 1-15.

10. 10. Meidner, D., Rannacher, R. and Vexler, B. (2011) A Priori Error Estimates for Finite Element Discretizations of Parabolic Optimization Problems with Pointwise State Constraints in Time. SIAM J. Control Optim., 49, 1961-1997. https://doi.org/10.1137/100793888

11. 11. Gong, W. and Hinze, M. (2013) Error Estimates for Parabolic Optimal Control Problems with Control and State Constraints. Computational Optimization and Applications, 56, 131-151. https://doi.org/10.1007/s10589-013-9541-z

12. 12. Lebedev, V.I. (1998) Explicit Difference Schemes for Solving Stiff Systems of ODEs and PDEs with Complex Spectrum. Russian J. Numer. Anal. Math. Modelling., 13, 107-116. https://doi.org/10.1515/rnam.1998.13.2.107

13. 13. Kadyrov, R.F., Laitinen, E. and Lapin, A. (2003) Using Explicit Schemes for Control Problems in Continuous Casting Process. Lobachevskii Journal of Mathematics, 13, 25-38.

14. 14. Lapin, A., Laitinen, E. and Lapin, S. (2015) Explicit Algorithms to Solve a Class of State Constrained Parabolic Optimal Control Problems. Russian J. Numer. Analysis Math. Modeling, 30, 351-362. https://doi.org/10.1515/rnam-2015-0032

15. 15. Maday, Y. and Turinici, G. (2002) A Parareal in Time Procedure for the Control of Partial Differential Equations. CRAS, 335, 387-392. https://doi.org/10.1016/S1631-073X(02)02467-6

16. 16. Lapin, A. and Romanenko, A. (2016) Uzawa-Type Iterative Method with Parareal Preconditioner for a Parabolic Optimal Control Problem. IOP Conf. Series: Materials Science and Engineering, 158.

17. 17. Lapin, A. (2010) Preconditioned Uzawa Type Methods for Finite-Dimensional Constrained Saddle Point Problems. Lobachevskii J. Math., 31, 309-322. https://doi.org/10.1134/S1995080210040013

18. 18. Laitinen, E., Lapin, A. and Lapin, S. (2010) On the Iterative Solution of Finite-Dimensional Inclusions with Applications to Optimal Control Problems. Comp. Methods in Appl. Math., 10, 283-301.

19. 19. Laitinen, E. and Lapin, A. (2012) Iterative Solution Methods for a Class of State Constrained Optimal Control Problems. Applied Mathematics, 3, 1862-1867. https://doi.org/10.4236/am.2012.312253

20. 20. Laitinen, E. and Lapin, A. (2013) Iterative Solution Methods for the Large-Scale Constrained Saddle Point Problems. In: Numerical Methods for Differential Equations, Optimization, and Technological Problems, Comp. Meth. Appl. Sc., 27, 19-39. https://doi.org/10.1007/978-94-007-5288-7_2

21. 21. Ladyzhenskaya, O., Solonnikov, V. and Ural’ceva, N. (1968) Linear and Quasilinear Equations of Parabolic Type. Transl. Math. Monographs, 23, AMS, Providence, RI.

22. 22. Lions, J.-L. and Magenes, E. (1972) Non-Homogeneous Boundary Value Problems and Applications. Springer.

23. 23. Quarteroni, A. and Valli, A. (1997) Numerical Approximation of Partial Differential Equations. Springer.

24. 24. Douglas Jr., J., and Gunn, J.E. (1964) A General Formulation of Alternating Direction Methods. Numerishe Mathematik, 6, 428-453. https://doi.org/10.1007/BF01386093

25. 25. D’Yakonov, E.G. (1961) The Method of Alternating Directions in the Solution of Finite Difference Equations. Doki. Akad. Nauk SSSR, 138, 271-274 (Russian).

26. 26. Yanenko, N. (1971) The Method of Fractional Steps. Springer. https://doi.org/10.1007/978-3-642-65108-3

27. 27. Samarsky, A.A. (2001) Theory of Difference Schemes. Marcel Dekker. https://doi.org/10.1201/9780203908518

28. 28. Zaitseva, S.B. and Zlotnik, A.A. (1996) Optimal Error Estimates of a Locally One-Dimensional Method for the Multidimensional Heat Equation. Mathematical Notes, 60, 137-146. https://doi.org/10.1007/BF02305177

29. 29. Barbu, V. (1984) Optimal Control of Variational Inequalities. Pitman Advanced Pub. Program.

30. 30. Lapin, A. and Laitinen, E. (2016) Preconditioned Uzawa-Type Method for a State Constrained Parabolic Optimal Control Problem with Boundary Control. Lobachevskii Journal of Mathematics, 37, 561-569. https://doi.org/10.1134/S1995080216050085