An adaptive semismooth Newton method for approximately solving control-constrained elliptic optimal control problems

Introduction

Lately, there has been significant progress in the context of optimal control of partial differential equations. Optimization with partial differential equation constraints is widely used in many fields, such as the natural sciences (Aniţa et al., 2011; Cantrell et al., 2009; Christofides et al., 2008) and industry (Capasso et al., 2000; Field and Komkov, 1992; Hu, 2012; Manchanda et al., 2017). Thus, the development, analysis and computer implementation of efficient optimization techniques is required for such optimization problems. Because of these properties, the solution of such problems has become increasingly popular. Also, in the past years, some studies have been published in the context of optimization with partial differential equation constraints (Borzì and Schulz, 2012; Hinze et al., 2008; Lions, 1971; Tröltzsch, 2010).

The goal of this paper is to prove the global convergence property for a nonmonotone semismooth Newton (NSN) method. Using this property, we can exploit an adaptive nonmonotone line search method for the approximate solution of elliptic optimal control problems with control constraints.

Let U , Y be reflexive Banach spaces and W be a Banach space. Consider the infinite dimensional optimization problem given by

min ( y , u ) ∈ Y × U J ( y , u )

(1)

subject to

E ( y , u ) = 0

(2)

u ∈ U ad

(3)

where J : Y × U → R is a quadratic cost functional and E : Y × U → W is an elliptic partial differential equation that is called the state equation. Also, the variables u and y are, respectively, called the control and the state variables, and the set of admissible controls U ad is given by

U ad = { u ∈ U | u a ≤ u ≤ u b } with u a , u b ∈ U

There are two main classes of numerical methods for solving optimal control problems: direct methods and indirect methods. In direct methods, the optimal control problem is transformed to a finite-dimensional optimization problem through discretization and then the obtained finite-dimensional optimization problem is solved using a suitable nonlinear optimization algorithm. In indirect methods, the optimality system is derived first and then the obtained optimality system is discretized using a suitable discretization method. In this work, an efficient direct method is used to control-constrained elliptic optimal control problems. To discretize such problems, a finite-difference scheme is used to discretize the state equation. Also, a Newton–Cotes rule is applied to discretize the cost functional. After discretization of the problem (equations (1) to (3)), the obtained finite-dimensional problem is given by

min ( y ^ , u ^ ) ∈ R ( n y + n u ) J ( y ^ , u ^ )

(4)

subject to

E ( y ^ , u ^ ) = 0

(5)

u ^ ∈ U ^ ad

(6)

where

J ( y ^ , u ^ ) : R ( n y + n u ) → R

E ( y ^ , u ^ ) : R ( n y + n u ) → R n y

and

U ^ ad = { u ^ ∈ R n u | u ¯ ≤ u ^ ≤ u ¯ } with u ¯ , u ¯ ∈ R n u

Here, n y and n u denote the number of state and control variables, respectively.

There are various numerical methods for solving the problem (equations (4) to (6)), which include penalty methods, augmented Lagrangian methods, active-set methods, interior point methods, gradient projection methods, trust region techniques or some combination of these methods (Bazaraa et al., 2006; Nocedal and Wright, 2006).

In this work, we prove the global convergence of a nonmonotone semismooth Newton (NSN) method developed by Bonettini and Tinti (2007) to exploit an adaptive nonmonotone line search method. Then, we implement this adaptive method to find the approximate solution of the control problem (equations (4) to (6)). In nonmonotone line search methods, some growth in the function value is permitted. Many researchers pointed out that utilizing nonmonotone schemes can improve the likelihood of finding a global optimum, as well as the convergence rate in cases where a monotone scheme is forced to creep along the bottom of a narrow curved valley (Bonettini and Tinti, 2007; Shen et al., 2012; Su and Pu, 2009; Su and Yu, 2009; Ulbrich and Ulbrich, 2003). Encouraging numerical results have been reported by combining nonmonotone techniques with a Newton-type direction for highly nonlinear and ill-conditioned optimization problems. Many partial differential equation constrained optimization problems have nonlinear and ill-posed structures. Because of this property, using the nonmonotone schemes is more reliable for these problems.

This discussion motivates us to apply a nonmonotone version of the semismooth Newton (NSN) method in the line search framework to find approximate solutions of control-constrained elliptic optimal control problems.

The structure of this paper is outlined as follows. The following section is devoted to discretization techniques for control-constrained elliptic optimal control problems. In the next section, the adaptive NSN method is described for the obtained discretized form. Numerical examples and conclusions are given and discussed in the final sections.

Control-constrained elliptic optimal control problems and their discretization procedure

Consider the following control-constrained elliptic optimal control problem

min ( y , u ) ∈ Y × U ad J ( y , u ) = 1 2 ∫ Ω ( ( y ( x ) − z ( x ) ) 2 + ν u 2 ( x ) ) d x

(7)

− Δ y − uy = f in Ω

(8)

y = 0 on ∂ Ω

(9)

where Ω = ( 0 , 1 ) 2 ⊂ R 2 , with boundary ∂ Ω , y ∈ Y ⊂ L 2 ( Ω ) and u ∈ U ad ⊂ L 2 ( Ω ) . Here, f ∈ L 2 ( Ω ) and z ( x ) ∈ L 2 ( Ω ) is the desired state and ν > 0 is a regularization parameter.

A finite-difference discretization scheme is used to discretize equations (8) and (9). By defining the mesh size steps h 1 = 1 / n and h 2 = 1 / m , with n , m ∈ N , we get the mesh

Ω h 1 h 2 = { ( x 1 i , x 2 j ) : x 1 i = i h 1 , x 2 j = j h 2 , i = 0 , 1 , … , n , j = 0 , 1 , … , m }

Let y ij represent an approximation to y ( x 1 i , x 2 j ) . To discretize equations (8) and (9), x 1 - and x 2 -derivatives are replaced with centred finite differences. Then equations (8) and (9) can be approximated in the following way

1 h 1 2 ( − y i − 1 , j + 2 y ij − y i + 1 , j ) ] + 1 h 2 2 ( − y i , j − 1 + 2 y ij − y i , j + 1 ) − u ij y ij = f ij

(10)

for i = 1 , … , n − 1 : j = 1 , … , m − 1 . For simplicity, we will consider the special case where h 1 = h 2 ≡ h . Therefore, equation (10) can be rewritten as follows

1 h 2 ( − y i − 1 , j − y i + 1 , j + 4 y ij − y i , j − 1 − y i , j + 1 ) − u ij y ij = f ij

(11)

for i , j … n − 1 . From equation (11) and by defining the vectors

y → = ( y [ 1 ] y [ 2 ] ⋮ y [ n − 1 ] ) , where y [ j ] = ( y 1 j y 2 j ⋮ y n − 1 , j ) u → = ( u [ 1 ] u [ 2 ] ⋮ u [ n − 1 ] ) , where u [ j ] = ( u 1 j u 2 j ⋮ u n − 1 , j )

and

f → = ( f [ 1 ] f [ 2 ] ⋮ f [ n − 1 ] ) , where f [ j ] = ( f 1 j f 2 j ⋮ f n − 1 , j )

the matrix form of equations (8) and (9) can be written as

A h y → − U h y → = f →

(12)

where

U h = diag ( u → ) = [ u 1 j u 2 j ⋱ u n − 2 , n − 1 u n − 1 , n − 1 ]

Here, A h is an ( n − 1 ) × ( n − 1 ) block tridiagonal matrix and has the form

A h = 1 h 2 ( T − I − I T − I − I T − I ⋱ ⋱ ⋱ − I T )

where each block T or I is itself an ( n − 1 ) × ( n − 1 ) matrix

T = ( 4 − 1 − 1 4 − 1 − 1 4 − 1 ⋱ ⋱ ⋱ − 1 4 )

and I is the ( n − 1 ) × ( n − 1 ) identity matrix. Also, the trapezoidal rule known as the Newton–Cotes integration formula is used to discretize the cost functional as follows

1 2 ∫ Ω ( ( y ( x ) − z ( x ) ) 2 + ν u 2 ( x ) ) d x = h 2 2 ∑ i = 1 n − 1 ∑ j = 1 n − 1 ( ( y ij − z ij ) 2 + ν u ij 2 ) + ν h 2 8 ( u 00 2 + u 0 n 2 + u n 0 2 + u nn 2 ) ︸ * * + ν h 2 4 ( ∑ i = 1 n − 1 u i 0 2 + ∑ i = 1 n − 1 u in 2 + ∑ j = 1 n − 1 u 0 j 2 + ∑ j = 1 n − 1 u nj 2 ) ︸ * * *

(13)

where h is defined as before, u ij = u ( x 1 i , x 2 j ) , u ¯ ij = u a ( x 1 i , x 2 j ) , u ¯ ij = u b ( x 1 i , x 2 j ) and z ij = z ( x 1 i , x 2 j ) . Furthermore, the square of terms u 00 , u 0 n , u n 0 , u nn , u i 0 , u in , u 0 j and u nj appear in (**) and (***) but they do not appear in equation (11). Consequently, the optimal values of these terms are zero and equation (13) can be written in the following matrix form

( ( y → − z → ) T B h ( y → − z → ) + ν u → T B h u → )

where

z → = ( z [ 1 ] z [ 2 ] ⋮ z [ n − 1 ] ) , where z [ j ] = ( z 1 j z 2 j ⋮ z n − 1 , j )

and

B h = h 2 2 ( I I I ⋱ I )

which is an ( n − 1 ) × ( n − 1 ) block diagonal matrix in which each block I is itself an ( n − 1 ) × ( n − 1 ) identity matrix.

Thus, the fully discretized problem is given as follows

min ( y → , u → ) ∈ R n y × U ^ ad J ( y → , u → ) = ( ( y → − z → ) T B h ( y → − z → ) + ν u → T B h u → )

(14)

where y → is the solution of

A h y → − U h y → = f →

(15)

and

u ¯ ≤ u → ≤ u ¯ →

(16)

where

u ¯ → = ( u ¯ [ 1 ] u ¯ [ 2 ] ⋮ u ¯ [ n − 1 ] ) , where u ¯ [ j ] = ( u ¯ 1 j u ¯ 2 j ⋮ u ¯ n − 1 , j )

u ¯ → = ( u ¯ [ 1 ] u ¯ [ 2 ] ⋮ u ¯ [ n − 1 ] ) , where u ¯ [ j ] = ( u ¯ 1 j u ¯ 2 j ⋮ u ¯ n − 1 , j )

NSN method to solve the obtained bound constrained optimization problem

This section deals with how the NSN method introduced by Bonettini and Tinti (2007) can be applied to solve the bound constrained optimization problem (equations (14) to (16)). To describe the NSN algorithm model, we need first to prove the existence of an optimal solution to the general control problem (equation (4) to (6)). The feasible set is denoted by

F : = { ( y ^ , u ^ ) : ( y ^ , u ^ ) ∈ R n y × U ^ ad , E ( y ^ , u ^ ) = 0 }

Now, consider the following assumptions:

Definition 1. A state-control pair ( y * , u * ) ∈ F is called optimal for equations (4) to (6), if

J ( y * , u * ) ≤ J ( y ^ , u ^ ) ∀ ( y ^ , u ^ ) ∈ F

Now, we introduce a theorem for the existence of an optimal solution ( y * , u * ) .

Theorem 1. Let assumptions (i) to (iv) hold. If n y , n u ∈ N such that n y , n u < ∞ , and the level set L θ : = { ( y ^ , u ^ ) ∈ F | J ( y ^ , u ^ ) ≤ θ } is a nonempty, closed and bounded set for θ < ∞ , the problem (equations (4) to (6)) has an optimal solution ( y * , u * ) .

Proof. Since J ≥ 0 and assumptions (i) to (iv) hold, the infimum

J * : = inf ( y ^ , u ^ ) ∈ F J ( y ^ , u ^ )

exists and hence a minimizing sequence { ( y ^ k , u ^ k ) } ⊆ F can be found with

lim k → ∞ J ( y ^ k , u ^ k ) = J *

For all k greater than some number K < ∞ , all ( y ^ k , u ^ k ) have to lie in L θ with θ > J * . Since L θ is closed and bounded, it is compact. Thus, there exists a convergent sub-sequence { ( y ^ k p , u ^ k p ) } ⊆ F with ( y ^ k p , u ^ k p ) → ( y * , u * ) as p → ∞ . Since assumptions (i) to (iv) hold, therefore J is continuous. Thus, ( y * , u * ) is a minimum, because

J * = lim p → ∞ J ( y ^ k p , u ^ k p ) = J ( y * , u * )

Under the assumptions (i) to (iv) and regarding equations (4) to (6), the implicit function theorem (Griva et al., 2009) allows us to consider y ^ = y ^ ( u ^ ) such that E ( y ^ ( u ^ ) , u ^ ) = 0 . Consequently, the reduced problem can be formulated as follows

min u ^ ∈ R n u J ^ ( u ^ ) subject to u ^ ∈ U ^ ad

(17)

where J ^ ( u ^ ) = J ( y ^ ( u ^ ) , u ^ ) is the reduced objective function and its reduced gradient is denoted ∇ u ^ J ^ ( u ^ ) . The implicit function theorem helps us to obtain ∇ u ^ J ^ ( u ^ ) . The implicit function theorem based on the assumptions (i) to (iv) guarantees the differentiability of y ^ ( u ^ ) . The Jacobian of y ^ ( u ^ ) is the solution of

E y ^ ( y ^ ( u ^ ) , u ^ ) y ^ u ^ ( u ^ ) = − E u ^ ( y ^ ( u ^ ) , u ^ )

where y ^ u ^ ( u ^ ) is the Jacobian of y ^ ( u ^ ) and E u ^ ( y ^ ( u ^ ) , u ^ ) is the partial Jacobian of E with respect to u ^ at ( y ^ , u ^ ) . Since y ^ ( u ^ ) is differentiable, the function J ^ is differentiable and its gradient is given by

∇ u ^ J ^ ( u ^ ) = y ^ u ^ T ( u ^ ) ∇ y ^ J ( y ^ ( u ^ ) , u ^ ) + ∇ u ^ J ( y ^ ( u ^ ) , u ^ ) = − E u ^ T ( y ^ ( u ^ ) , u ^ ) E y ^ − T ( y ^ ( u ^ ) , u ^ ) ∇ y ^ J ( y ^ ( u ^ ) , u ^ ) + ∇ u ^ J ( y ^ ( u ^ ) , u ^ ) = − E u ^ T ( y ^ ( u ^ ) , u ^ ) p ^ ( u ^ ) + ∇ u ^ J ( y ^ ( u ^ ) , u ^ )

(18)

where p ^ ( u ^ ) is the solution of

E y ^ T ( y ^ ( u ^ ) , u ^ ) p ^ ( u ^ ) = ∇ y ^ J ( y ^ ( u ^ ) , u ^ )

Theorem 2. Let the assumptions (i) to (iv) hold. If u * is a local optimal solution of the reduced problem (equation (17)), then ∇ u ^ J ^ ( u * ) T ( u ^ − u * ) ≥ 0 for all u ^ ∈ U ^ ad .

Proof. The convexity of U ^ ad yields

γ u ^ + ( 1 − γ ) u * = u * + γ ( u ^ − u * ) ∈ U ^ ad

for all u ^ ∈ U ^ ad and γ ∈ [ 0 , 1 ] . As u * is a local optimal solution, one has

J ^ ( u * ) ≤ J ^ ( u * + γ ( u ^ − u * ) )

for γ sufficiently small. Dividing by γ > 0 and taking the limit as γ → 0 yields

∇ u ^ J ^ ( u * ) T ( u ^ − u * ) = lim γ → 0 J ^ ( u * + γ ( u ^ − u * ) ) − J ^ ( u * ) γ ≥ 0

From equation (18), the adjoint representation of the reduced gradient ∇ u → J ^ ( u → ) for the bound constrained optimization problem (equations (14) to (16)) is given by

∇ u → J ^ ( u → ) = ( 2 ν B h u → + Y h p → )

(19)

where p → ∈ R n y is the Lagrange multiplier or the adjoint state for the equality constraint (equation (15)) and solves the adjoint equation

( A h − U h ) p → = 2 B h ( y → − z → )

Combining Theorem 2 and equation (19) leads to the following inequality

( 2 ν B h u → ∗ + Y h p → ) T u → ∗ ≤ ( 2 ν B h u → ∗ + Y h p → ) T u → , ∀ u → ∈ U ^ ad

which implies that u → * is a solution of the minimization problem

min u → ∈ R n u ( 2 ν B h u → * + Y h p → ) T u → subject to u → ∈ U ^ ad

(20)

where Y h = diag ( y → ) . From equation (20), it also follows that

u → * = { u ¯ → if ( 2 ν B h u → * + Y h p → ) < 0 ∈ ( u ¯ , u ¯ → ) if ( 2 ν B h u → * + Y h p → ) = 0 u ¯ if ( 2 ν B h u → * + Y h p → ) > 0

(21)

The Lagrangian function for the bound constrained optimization problem (equation (14) to (16)) is given by

L ( y → , u → , p → , λ ¯ , λ ¯ ) = ( y → − z → ) T B h ( y → − z → ) + ν u → T B h u → − p → T ( A h y → − U h y → − f → ) + λ → ¯ T ( u ¯ − u → ) + λ ¯ → T ( u → − u ¯ → )

where p → ∈ R n y is defined as before. Also, λ → ¯ ∈ R n u and λ ¯ → ∈ R n u represent the Lagrange multipliers for the inequality constraints on the control variable. Now, the following proposition determines the so-called Karush–Kuhn–Tucker conditions as the necessary optimality conditions (Bazaraa et al., 2006; Nocedal and Wright, 2006).

Proposition 1. Let ( y → * , u → * ) be an optimal solution for the control problem (equation (14) to (16)). Then, there exist Lagrange multipliers p → ∈ R n y and λ → ¯ , λ ¯ → ∈ R n u such that

{ ∇ y → L ( y → * , u → * , p → , λ → ¯ , λ ¯ → ) = 0 ⇒ ( A h − U h * ) p → − 2 B h ( y → * − z → ) = 0 ∇ u → L ( y → * , u → * , p → , λ → ¯ , λ ¯ → ) = 0 ⇒ − Y h * p → + λ → ¯ − λ ¯ → − 2 ν B h u → * = 0 ∇ p → L ( y → * , u → * , p → , λ → ¯ , λ ¯ → ) = 0 ⇒ ( A h − U h * ) y → * − f → = 0 λ → ¯ ≥ 0 , λ ¯ → ≥ 0 λ → ¯ T ( u ¯ − u → * ) = λ ¯ → T ( u → * − u ¯ → ) = 0 u ¯ ≤ u → * ≤ u ¯ →

(22)

where U h * = diag ( u → * ) and Y h * = diag ( y → * ) .

From equation (21) and the second optimality condition in equation (22), the last three optimality conditions in equation (22) can be reformulated as

λ → = max ( 0 , λ → + c ( u → * − u ¯ → ) ) + min ( 0 , λ → + c ( u → * − u → ¯ ) )

(23)

where c > 0 and λ → = λ → ¯ − λ ¯ → with λ → ¯ = max ( 0 , λ → ) and λ ¯ → = − min ( 0 , λ → ) .

The first three optimality conditions in equations (22) and (23) for ( y → , u → ) lead to the following system to be solved

H ( y → , u → , p → , λ → ) = ( ( A h − U h ) y → − f → ( A h − U h ) p → − 2 B h ( y → − z → ) − Y h p → + λ → − 2 ν B h u → λ → − max ( 0 , λ → + c ( u → − u ¯ → ) ) − min ( 0 , λ → + c ( u → − u → ¯ ) ) ) = 0

(24)

Since the max and min functions are not differentiable, the classical Newton method cannot be applied. For this reason, the semismooth Newton method can be used to solve such systems. The concept of the semismooth Newton method has been introduced first by Pang and Qi (1993), Qi (1993) and Qi and Sun (1993) and then developed in many papers and books (Chen et al., 2000; Chi et al., 2017; Hinze et al., 2008; Ito and Kunisch, 2008; Kröner et al., 2011). The semismooth Newton iteration can be stated regarding an active set strategy and the generalized derivatives of the max and min operators.

Let us set c = α and define the active sets as

A + = { x ij ∈ Ω h : λ → + α ( u → − u ¯ → ) > 0 } A − = { x ij ∈ Ω h : λ → + α ( u → − u → ¯ ) < 0 }

where Ω h = { x ij = ( x 1 i , x 2 j ) : x 1 i = ih , x 2 j = jh , i , j = 1 , … , n − 1 } . Also, I = Ω h \ ( A + ∪ A − ) denotes the inactive set. The generalized derivatives of the max operator, max ( 0 , λ → + α ( u → − u ¯ → ) ) , and min operator, min ( 0 , λ → + α ( u → − u → ¯ ) ) , can be defined as

G max ( λ → + α ( u → − u ¯ → ) ) = χ A + G min ( λ → + α ( u → − u → ¯ ) ) = χ A −

where χ A + and χ A − denote the characteristic functions of A + and A − , respectively.

The last two equations in equation (24) imply that

u → = { u ¯ → in A + − Y h p → ν h 2 in I u ¯ in A −

(25)

The semismooth Newton step to equation (24) is determined by the solution of the following system and the update of equation (25)

G ( y → , p → , λ → ) . d → = − H 1 ( y → , p → , λ → )

(26)

where

G ( y → , p → , λ → ) = ( ( A h − U h ) 0 0 − 2 B h ( A h − U h ) 0 0 0 C − D ) d → = ( d → y d → p d → λ )

and

H 1 ( y → , p → , λ → ) = ( ( A h − U h ) y → − f → ( A h − U h ) p → − 2 B h ( y → − z → ) λ → − max ( 0 , λ → + c ( u → − u ¯ → ) ) + min ( 0 , λ → + c ( u → − u ¯ ) ) )

Also, C is an ( n − 1 ) 2 × ( n − 1 ) 2 identity matrix and D = diag ( χ A ) where χ A denotes the characteristic function of A = A + ∪ A − .

Once the direction d → k is computed in iteration k, a nonmonotone line search procedure introduced as follows is used to compute the step size α k . The nonmonotone globalization technique is useful in difficult nonlinear problems, because it may help in escaping from steep-sided valleys and may improve both the possibility of finding the global optimum and the rate of convergence. Bonettini and Tinti (2007) introduced a NSN method that generates a sequence { t → k } such that

∥ H 1 ( t → k ) + G ( t → k ) . d ∥ ≤ η k ∥ H 1 ( t → l ( k ) ) ∥

(27)

and

∥ H 1 ( t → k + α k d k ) ∥ ≤ ( 1 − α k β ( 1 − η k ) ) ∥ H 1 ( t → l ( k ) ) ∥

(28)

where t → k = ( y → k , p → k , λ → k ) , η k ∈ [ 0 , 1 ) , β ∈ ( 0 , 1 ) and

∥ H 1 ( t → l ( k ) ) ∥ = max 0 ≤ j ≤ m ( k ) ∥ H 1 ( t → k − j ) ∥ , k = 0 , 1 , 2 , …

(29)

where m ( 0 ) = 0 and 0 ≤ m ( k ) ≤ min { m ( k − 1 ) + 1 , M } with M ≥ 0 . Here, the fixed parameter M is the maximum of the memory element m ( k ) .

Now, by attention to these descriptions, we establish the global convergence of this NSN method in the following theorem.

Theorem 3. Assume that at each iterate k, it is possible to find a vector d k such that equation (27) is satisfied. If ∑ k ≥ 0 ( 1 − η k ) is divergent, then H 1 ( t → k ) → 0 . Furthermore, if t → * k is a limit point of { t → k } , then H 1 ( t → * k ) = 0 and if G ( t → * k ) is nonsingular, then the sequence { t → k } converges to t → * k .

Proof. By equation (28), we have