{'fraction_non_alphanumeric': 0.0368682128971712, 'fraction_numerical': 0.007909907494302185, 'mean_word_length': 4.301705756929637, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On optimal control of free boundary problems of obstacle type

A numerical study of an optimal control formulation for a shape optimization problem governed by an elliptic variational inequality is performed. The shape optimization problem is reformulated as a boundary control problem in a fixed domain. The discretized optimal control problem is a non-smooth and non-convex mathematical programing problem. The performance of the standard BFGS quasi-Newton method and the BFGS method with the inexact line search are tested.


Introduction
In the present paper a numerical method for shape optimization for a class of free boundary problems is proposed. The shape optimization problem is reformulated as a boundary control problem following the ideas discussed in [11]. The optimal shape is obtained as a level set. Thus, in numerical realization using finite elements, there is no need for remeshing during optimization iterations.
As a concrete free boundary (state) problem we consider the contact problem for the Poisson problem. In this case the free boundary is defined by the a-priori unknown contact zone. As a shape optimization problem we consider the problem of finding a shape having minimum area such that the contact zone includes a given subdomain [2].
It is well-known that in control problems for contact problems the control-to-state mapping is not smooth, in general. Despite this, many nonlinear programming codes for smooth optimization can be used to obtain reasonable approximate solutions to those control problems.
The paper is organized as follows. In Section 2 we present a shape optimization problem and an alternate boundary control formulation of it. Section 3 is devoted to disceretization and algebraic sensitivity analysis. Finally, a numerical study on the performance of two nonlinear programming codes applied to the discretized control problem is done in Section 4.

Setting of the problem
We consider an abstract shape optimization problem given formally as follows: Here α is the optimization (control) parameter, A represents a free boundary problem (in abstract operator form) whose solution consists both of the domain Ω α where the PDE is posed as well as the solution y α to that PDE.
Problem (1) is more complicated than a classical shape optimization problem. In the latter the state problem is solved on a fixed domain, in the sense that Ω α is considered to be given instead of being one of the unknowns in the state problem.
Next we introduce the model problem to be studied in the rest of this paper. Let Ω ⊂ R 2 a domain, f ∈ L 2 (Ω), ψ ∈ C 1 (Ω). We consider the following classical free boundary problem: Problem (2) models e.g. the contact between elastic membrane (represented by Ω) and a rigid obstacle defined by ψ. The solution y Ω to (2) then gives the vertical displacement of the membrane under the vertical load represented by f . The set Z Ω is the contact zone between the membrane and the obstacle. Problem (2) can be formulated as the following variational inequality ( [6]): y Ω ∈ K 0 (Ω) : where This formulation does not include explicitly the contact zone Z Ω as an unknown.
Let Ω ⊂ R 2 and Ω 0 ⊂ R 2 be given bounded domains such that Ω 0 Ω. We introduce a system of bounded domains Let y Ω be the unique solution to (3) and let C(y Ω ) = {x ∈ Ω | y Ω (x) = ψ(x) a.a. x} denote the corresponding contact region. We consider the following shape optimization problem introduced in We are thus looking for a domain Ω (representing the membrane) having minimum area such that the contact region "covers" the given domain Ω 0 .
The state constraint Ω 0 ⊂ C(y Ω ) can be relaxed by introducing a penalty term resulting in the shape optimization problem One of the main difficulties in the theoretical and numerical treatment of the problem (5) consists in the variable character of the domain Ω on which the state problem (3) is given. Theoretical and numerical aspects of shape optimization using "moving" domains are discussed e.g. in refs [12], [13], [8], [7].
Alternatively, there has been significant amount of interest in fixed domain formulations of shape optimization problems (see e.g. [5]). The most common approaches to get rid of moving domains are the following: • By scaling the domain such that it becomes fixed, the optimization parameter appears as a coefficient in the state problem.
• The domain Ω is represented by a level set Ω = {x ∈ Ω | Φ(x) < 0}, where Φ : Ω → R is an unknown level set function to be determined.
• The state problem is modified by adding control variable to it. An optimal shape is defined implicitly by a level set (but without separate level set function).

Boundary control approach
In this paper we utilize the boundary control approach discussed in [11]. In what follows we assume that f ≤ 0 and ψ < 0. We define the set of admissible controls as follows where u min , u max > 0 are given constants. Next, to each u ∈ U we associate y(u), the solution to the variational inequality Let H : R → R denote the Heaviside step function. We consider the following boundary control problem posed in a fixed domain: subject to (6).
As u > 0 on ∂ Ω and f ≤ 0, ψ < 0 in Ω it follows that for suitable choice of u the set of points where y(u) is strictly negative is non-empty. The first term in the cost functional J ε gives the area of this set while the latter term adds a penalty if the contact zone C(y(u)) does not cover Ω 0 . Assume that there exists an optimal pair (u ε , y ε ) := (u ε , y(u ε )) for (7). We can now consider Ω ε := {x ∈ Ω | y ε (x) < 0} as an approximate solution to the original shape optimization problem (4).
The control approach clearly makes sense if Ω 0 is simple enough star-like domain. If, however, Ω 0 is e.g. multiply connected the set Ω ε might not approximately solve (5).

Approximation and numerical realization
To realize (7) numerically, we must discretize both the control and state variables. In what follows, we assume that Ω is the disk B(0, R). We discretize the state problem (6) by using piecewise linear triangular elements. Instead of the exact Heaviside function, we use the smoothed one H β (x) := 1 2 tanh( x β ) + 1 2 , β > 0. The use of piecewise linear triangular elements implies that the most obvious way to discretetize the control u would be to use piecewise linear and continuous discretization in the same finite element mesh where the state variable is discretized. However, this approach has two well-known drawbacks. Firstly, the number of optimization variables is very large whenever dense meshes are used. Secondly, the piecewise linear approximation (without a suitable regularization term in the cost function) is prone to spurious oscillations.
Instead, we look for a differentiable and periodic function u a : [0, 2π] → R that is fully defined by a vector of parameters a = (a 1 , a 2 , ..., a n ). Examples of such functions are Bezier functions, cubic Hermite, and cubic spline interpolating polynomials, e.g.
Here we shall use shape preserving periodic piecewise cubic Hermite interpolation polynomial (see [4], [1]) to represent the control function. The parametrized control function solves the following interpolation problem The advantage of this kind of parametrization is that the number of discrete optimization variables is small but at the same time the control is a smooth function without excessively prone to wild oscillations. Moreover, the interpolant u a does not overshoot the data, so u a ∈ U if u min ≤ a i ≤ u max ∀i.
Sensitivity analysis for the discrete state problem  (6) can be expressed as an equivalent quadratic programming problem where K and f are the stiffness matrix and the force vector, respectively. The discretization of the boundary control problem (7) leads to the nonlinear programming problem minimize a J ε (a) subject to (8) and u min ≤ a k ≤ u max , k = 1, ..., n.
To be able to use descent type optimization methods we need to evaluate the gradient of J ε with respect to the discrete control variable a.
Let I c (a) denote the set of contact nodes, i.e. the solution q = q(a) satisfies q i = ψ i , i ∈ I c (a). Let a be the current approximation of the optimal control. If we assume that I c ( a) is known a priori and it is invariant under small perturbations of a, then J is differentiable and the gradient ∇ a J ε ( a) can be obtained using the standard adjoint equation technique. However, if the contact set changes due to arbitrary small perturbation, the mapping a → J ε (a) is not differentiable at a and the mechanical application of the adjoint technique gives at most an element from the subdifferential.

Additional source of nonshmoothness
Consider piecewise cubic Hermite interpolation of the data (x 1 , y 1 ), (x 2 , y 2 ), ...., (x n , y n ). On each subinterval [x k , x k+1 ] the interpolant can be expressed in terms of a local variable s = x − x k as follows The interpolant satisfies the following interpolation conditions In the classical piecewise cubic Hermite interpolation, the slope parameters d k , k = 1, ..., n are a priori given constants.
In shape-preserving piecewise cubic interpolation, the slopes are not supplied by the user. Instead, they are computed algorithmically as a part of the interpolating process in the following way: Let δ k = (y k+1 − y k )/(x k+1 − x k ) and h k = x k+1 − x k . The slopes d k are determined as follows. If sign(δ k ) · sign(δ k−1 ) < 0 we set d k = 0. If sign(δ k ) · sign(δ k−1 ) > 0 and the two intervals have the same length, then d k is the harmonic mean of δ k−1 and δ k : If the intervals have different lengths, then d k is the following weighted harmonic mean The advantage of shape-preserving piecewise cubic interpolation over cubic spline interpolation is that the interpolant does not overshoot the data.
Consider next the case where the vector of values y = (y 1 , ..., y n ) is an optimization parameter. In case of cubic spline interpolation, the vector of slopes is obtained by solving a tridiagonal system Ad = b, where the vector b is a linear function of y. Thus the value of the interpolant at given fixed point x is a smooth function of y.
With shape-preserving piecewise cubic interpolation the situation is different. The mapping y → d is clearly nonlinear. In addition, it is nonsmooth. A simple example shown in Figure 1 demonstrates this.

Numerical examples
In this section we present two numerical examples. All finite element computations are implemented in Matlab [10]. The quadratic programming problem (8) is solved using Matlab's quadprog function. The box constraints u min ≤ a i ≤ u max are implemented by adding a smooth penalty function  The optimal control problem is solved by using two different optimizers: • Function fminunc from Matlab's Optimization Toolbox with 'quasi-newton' option that implements the standard BFGS method.
For the options of fminunc that control the optimality tolerance we use their default values. Also we use the default options of HANSO. Especially, we do not use the option available in HANSO to continue the optimization using the gradient sampling method as this option is very expensive in terms of the number of function evaluations.
Example 1. Here we use an unstructured finite element mesh with nominal mesh size h = 0.05. The mesh is constructed in such way that ∂Ω 0 coincides element edges.
We ran both optimizers using n = 30 and a (0) i =2, i=1, ..., n until they were not able to improve the control. The evolution of the best objective value versus computational work is depicted for both methods in Figure 2. The implicitly defined optimized domains Ω ε (with the corresponding contours of y ε ) and the optimal boundary controls are plotted in Figures 3  and 4.
One can observe that in this example HANSO 2.2 finds better control than fminunc. However, HANSO 2.2 spends considerable amount of computational work with only very marginal improvement in the cost function.
Example 2. In this example the mesh is generated in the same way as in Example 1 except the nominal mesh size is h = 0.025. The evolution of the best objective value versus computational work is depicted for both methods in Figure 5. The implicitly defined optimized domains Ω ε (with the corresponding contours of y ε ) and the optimal boundary controls are plotted in Figures 6 and 7. This time both optimizers end up to almost same final objective value using essentially same amount of computational work. However, the optimized controls and the corresponding domains Ω ε differ considerably.
The results of Examples 1 and 2 show that the method works reasonably well. The obtained shapes are approximately feasible (due to penalization) and the objective function value was substantially reduced. However, there is no guarantee that any of the obtained domains is even a local minimizer.
From the results, we can also deduce that the boundary control problem is ill-conditioned. The effect of ill-conditioning can be seen by comparing the optimized controls obtained by using the two different optimizers. Comparing the convergence history we see that significant differences in the details of the boundary control have only a minor effect on the value of the objective function. This is due to the fact that the solution of the Poisson problem is generally smoother than the data. Similar stiffness appears e.g. in identification problem related to the Bernoulli free boundary problem studied in [14].

Conclusions
In this paper we have considered some computational aspects of a shape optimization problem with the state constraint given by a free boundary/obstacle problem. The free boundary problem is formulated as a quadratic programming problem which is then solved using the state-of-the-art tools. To solve the optimal shape design problem, a boundary control approach was used. Its main advantage is that there is no need to consider moving domains. This is advantageous especially from the computational point of view.
A well-known feature of optimal control problems governed by obstacle type problems is that the control-to-state mapping is not smooth in general. However, in discrete setting, it is piecewise smooth. In those points where it is smooth, the gradient of the objective function can be evaluated in a straightforward way using the adjoint approach.
Numerical examples show that the location of the boundary of the contact zone can be adjusted by changing the boundary control. However, it seems that the problem is illconditioned in the sense that relatively large changes in the boundary control have only a little effect on the location of the contact zone boundary. We did not make a detailed tradeoff study between accuracy and oscillations. That will be a topic in further studies. One possible remedy (and a topic for further studies) would be to consider distributed control instead of boundary control as was also done in [11].