A penalization scheme for the numerical solution of optimal control problems with control complementarity constraints

We suggest a simple penalty method for the numerical solution of optimal control problems with control complementarity constraints which is based on the famous Fischer–Burmeister function. The distinct advantage of this approach is the resulting smoothness of the penalty term and the unconstrainedness of the associated penalized surrogate problems. A numerical example from elliptic PDE control is presented as well.


Introduction
Mathematical programs with complementarity constraints (MPCCs) can be used to model numerous real-world applications where mechanical, chemical, or economical equilibrium systems have to be optimized. Noting that MPCCs are inherently irregular, huge effort has been put into the development of suitable problem-tailored optimality conditions, constraint qualifications, and numerical methods during the last decades. However, most of these contributions address the finite-dimensional setting while the infinite-dimensional scenario has received comparatively little attention; nevertheless, a growing number of underlying applications demands some respective progress. In this work, we consider optimal control problems with pointwise control complementarity constraints, where a pair of distributed controls should have the property that each control is non-negative while their product vanishes. We describe a numerical solution approach which is based on the penalization of the full complementarity constraints (including the non-negativity that is usually treated separately) using the well-known Fischer-Burmeister function. This leads to a sequence of smooth and unconstrained problems that can be solved using standard methods.

Problem statement
For a bounded Lipschitz domain Ω ⊂ R d , we consider the optimal control problem with control complementarity constraints where S : H 1 (Ω) × H 1 (Ω) → H is a linear and continuous operator which maps each pair of controls to some element of the observation Hilbert space H. Classically, S is the composition of the solution operator for a given ODE or PDE and a suitable embedding into H. Furthermore, y d ∈ H is a given desired state, and α 1 ≥ 0, α 2 ≥ 0, as well as ε > 0 are fixed. Emphasizing that the controls are chosen from the Sobolev space H 1 (Ω) which is compactly embedded into the Lebesgue space L 2 (Ω), the feasible set of (OC 4 ) is weakly sequentially closed. Furthermore, the objective functional of (OC 4 ) is convex, continuous, and (due to ε > 0) H 1 -coercive. Thus, the existence of an optimal solution to (OC 4 ) follows by standard arguments. Let us note that similar arguments do not apply to the situation where the controls are chosen from L 2 (Ω) while ε := 0 holds. In this case, the associated feasible set is not weakly sequentially closed.

A penalty method
For the construction of our solution method, we exploit the well-known Fischer-Burmeister function ϕ : [3]. The penalty term Φ : , v(x))dx for all u, v ∈ H 1 (Ω).

of 2
Section 16: Optimization Using technical arguments based on the smoothness of the squared Fischer-Burmeister function provided in [1], one can show that Φ is continuously Fréchet differentiable. Next, we choose a sequence {σ k } k∈N of positive penalty parameters which tend to ∞ as k → ∞ and consider the family of unconstrained optimization problems Standard arguments can be applied in order to show that (P k ) possesses a global minimizer for each k ∈ N. The following theorem summarizes the qualitative properties of the associated penalty method. Theorem 3.1 For each k ∈ N, let (ū k ,v k ) ∈ H 1 (Ω) × H 1 (Ω) be a global minimizer of (P k ). Then, {(ū k ,v k )} k∈N possesses a convergent subsequence whose limit point is a global minimizer of (OC 4 ).
Noting that Φ is smooth, the local minimizers of (P k ) can be characterized with the aid of associated first-order optimality conditions. On the other hand, one has to keep in mind that (P k ) is not convex, which means that these optimality conditions are, in general, not sufficient. Particularly, in numerical practice one cannot ensure that the penalized problems (P k ) are solved to global optimality. In this regard, one has to check whether the actual output of the suggested penalty method is reasonable. One possibility to do so is to test whether the computed point satisfies a strong stationarity-type optimality condition which can be shown to hold at each local minimizer of (OC 4 ), see [1].
Here, χ A : D → R denotes the characteristic function of A ⊂ D which equals 1 on A and 0 otherwise. Above, we use D u := Ω × (0, 1 4 ) and D v := Ω × ( 3 4 , 1). The other data in (OC 4 ) is chosen to be y d := 2χ Dd + 1 where we used 1)], α 1 = α 2 := 0, and ε := 10 −5 . For the numerical solution of the problem, we choose σ k := 10 k for each k ∈ N. Next, a discretization using piecewise affine finite elements for state and control is performed on (P k ) and (PDE) in order to respect the overall H 1 -setting. Afterwards, the associated discrete first-order system is derived. Noting that the (discretized) derivative of the penalty term Φ is nonsmooth but Lipschitz continuous, this system can be solved with the aid of a damped semismooth Newton method. Putting things together, the overall algorithm, which is based on a first-discretize-then-optimize strategy, possesses an outer loop where the penalty parameter is increased and an inner loop where the penalized and discretized surrogate problems are solved with a Newton-type method. In order to start the algorithm, we compute the solution of the program (OC 4 ) where the controls are only demanded to be non-negative but not necessarily complementary, see [2] for details on the solution of such problems. These initial controls as well as the output of our algorithm are visualized in Figure 1. The above mentioned test indicates that the computed point is indeed strongly stationary for the associated problem (OC 4 ) up to a numerical tolerance.