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      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  . In     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  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  and recently applied to a parabolic state constrained problem . One more approach to construct effective algorithm for parabolic optimal control problem is to use parareal approximation of the state equation  .
In this article we continue the investigations of       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
in the cylinder , , with lateral surface . We call and as state and control functions. It is well-known that for any there exists a unique weak solution of problem (1) such that , and the following inequality takes place (  ):
Define objective function
with a given functions , and the sets of constraints for state and control:
Above and . We solve the following optimal control problem:
Problem (3) has a unique solution (cf., e.g. ).
We construct the finite-difference approximations of problem (3) using a uniform in x and t mesh in .
Let be the uniform mesh of the meshsize on , be the set of the boundary nodes while . By we denote the space of mesh functions which are defined on and vanish in the boundary nodes , let . The mesh on the time segment we denote by .
We will use the notations for the mesh functions from and for the vectors of their nodal values as well. We also don’t distinguish the linear operators from to and corresponding them matrices acting on the vectors of nodal values of mesh functions from .
By we denote the values on a time level of a mesh function of x and t, and by -euclidian norm in the space . We use the following notations: , and are the inner product and the norm in . Let be unit matrix while be the unit matrix in .
We denote by the linear operator such that
and similarly for Obviously, is the mesh Laplasian with homogeneous Dirichlet boundary conditions. Recall that we use the same notations for corresponding matrices. Matrices are symmetric, commute and positive definite with spectrum in a segment , where is of order , while is bounded below by a constant which doesn’t depend on .
Let us introduce matrices U, D and R, where U is a symmetric and positive matrix, and approximate state problem (1) by two-level finite difference scheme
Below we give several examples of two-level finite difference scheme (4).
1) Finite difference scheme with weights:
with can be written in the form (4) with , and . Recall that (5) contains forward Euler ( ), backward Euler ( ) and Crank-Nicolson ( ) schemes.
Scheme (5) is unconditionally stable for and stable in the case if , where is the maximal eigenvalue of . The stability estimate is:
2) Two-level scheme with factorized preconditioner:
Mesh scheme (7) is unconditionally stable for with stability estimate (6).
3) Fractional steps scheme:
System (8) can be written in form (4) with
Scheme (8) is unconditionally stable with stability estimate (6).
Let us further use notation for an -approximation of the function . Define a mesh goal function and the sets of the constraints:
The mesh optimal control problem reads as follows:
Problem (10) has a unique solution because the set is a convex compact and the quadratical function is continuous and strictly convex on .
3. Saddle Point Problem and Iterative Solution Method
Let us rewrite problem (10) in a “vector-matrix” form. Let the matrices be defined by the equalities:
Denote by and the indicator functions of the sets and , respectively. We obtain the following algebraic form of mesh optimal control problem (10):
Below we suppose that for a two-level finite difference scheme (4) the following assumption is valid:
The stability estimate (6) approve that for all cited above particular cases of (4) assumption (12) is true (with constant ).
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 ,  and “sequential” variant of fractional steps scheme   etc. Also a more general variant of scheme (8) (see ) with positive weights , such that instead of as in (8) can be considered.
We construct Lagrange function for problem (11):
A saddle point of this Lagrangian satisfies the following system:
where and are the subdifferentials of the corresponding functions.
With the notations , , and
(11) becomes a particular case of minimization problem
while (13)-a particular case of saddle point problem
We will use the following results from the article :
Proposition 1 Let
1) Problem (15) has a non-empty set of the solutions , where z is unique solution of (14).
2) Uzawa-type iterative method
with a symmetric preconditioner satisfying the inequality
converges for any initial guess : for .
Theorem 1 Problem (13) has a solution with unique pair , 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 satisfies assumption (19).
Below we construct for problem (13) the easily implementable preconditioner D which is spectrally equivalent to with constants independent on and h. As this preconditioner we take . Then method (20) for problem (13) with this preconditioner reads as follows:
Theorem 2 Method (22) converges for any initial guess if
Proof. The matrix is spectrally equivalent to , namely,
where constant is defined in (12). In fact, by direct calculation we find . Further, since , then . This inequality is equivalent to
whence the result.
In virtue of inequality (24) the convergence condition (21) of proposition 1 is true for the parameter
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 . Namely, the following statement holds:
Theorem 3 Let . Then there exists a unique solution of saddle point problem (13), and for theoretically optimal iterative parameter
the following estimate for the rate of convergence of method (22) is
Proof. The uniqueness of is proved in theorem 1. The uniqueness of in the case follows from the equation .
Vector is the solution of the equation
while iterative method (22) can be written in the form
It is well-known (cf., e.g. ) that the operator is co-coercive:
Because of this and (25) the operator satisfies the following properties (strong monotonicity and Lipshitz-continuity):
The rest of the prove is quite standard (cf. ). Namely, with the notations , and we have the equation
We multiply this equation by and obtain the equality
Due to the properties of the following estimate holds:
Substituting this estimate in the previous equality we get
and rate of convergence (26) is true for optimal parameter
Remark 2 Due to the equalities , , we have the following estimate for the rate of convergence of :
For the sequence of control vectors the estimate is as follows:
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
which operator is only co-coercive. The convergence of the iterative methods for such kind of equations have been investigated in .
On the other hand, numerous calculations show that in this general case the choice of the iterative parameter 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 and , and a system of linear equations with matrix . 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 is equivalent to solving direct (with ) and adjoint (with ) 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 , . 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 and is spectrally equivalent to “main” matrix of the problem;
・ Equations with the matrices and as in (7) and (8) can be easily implementable.
The first property is ensured by the inequality (12): . Just this inequality allows us to prove spectral equivalency of the matrix and with the constants independent on mesh parameters and . 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 .
The second property-easy solution of the equations with matrices , 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.   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  and .
The research was supported by the RFBR grant No. 16-01-00408 and by grants No. 296348 and No. 296928 from Academy of Finland.