Relaxation approach for optimization of free boundary problems

We consider an optimal control problem (OCP) constrained by a free boundary problem (FBP). FBPs have various applications such as in fluid dynamics, flow in porous media or finance. For this work we study a model FBP given by a Poisson equation in the bulk and a Young‐Laplace equation accounting for surface tension on the free boundary. Transforming this coupled system to a reference domain allows to avoid dealing with shape derivatives. However, this results in highly nonlinear partial differential equation (PDE) coefficients, which makes the OCP rather difficult to handle. Therefore, we present a new relaxation approach by introducing the free boundary as a new control variable, which transforms the original problem into a sequence of simpler optimization problems without free boundary. In this paper, we formally derive the adjoint systems and show numerically that a solution of the original problem can be indeed asymptotically approximated in this way.


INTRODUCTION
Free boundary problems (FBPs) occur in various fields of application.These include, for example, solidification, technical textiles, or thin film manufacturing processes [1][2][3].Phase separation and optimal control problems (OCPs) constrained by the Stefan problem have been widely studied [2,[4][5][6].Many of these problems have in common that there is a relation between the surface tension and the curvature of the boundary [7,8].
In this paper, we study a model FBP that is very similar to the one considered by Saavedra and Scott [9].Analytical results such as existence and uniqueness of solutions, existence of a minimizer and first and second order optimality conditions are shown by Antil et al. [10,11].Throughout this work we assume that the free boundary Γ  (dashed line) is represented as a graph of a function , cf. Figure 1 .The model consists of two equations: a Poisson equation, which describes the quantity in the interior of the domain, and a Young-Laplace equation, which connects the surface tension to the curvature of the free boundary.
Free boundaries make the optimization highly complex and costly, particularly the derivation of the necessary optimality conditions and efficient numerical handling.When solving OCPs with FBPs, besides the state variables, the domain is also an unknown.To avoid having to deal with shape derivatives, we transform the FBP to a fixed domain.This is done at the expense of highly nonlinear PDE coefficients.
F I G U R E 1 Ω  is the physical domain of the original FBP with boundary Ω  = Σ ∪ Γ  .The fixed part of the boundary is denoted by Σ and Γ  is free.We assume that the top boundary is given in graph form through .By transforming the domain Ω  onto a reference domain Ω, we obtain a fixed top boundary Γ.The complete boundary of Ω is given by Ω = Γ ∪ Σ.All analysis and computations are done on Ω.
In this paper, we propose a relaxation-based approach that allows solving a sequence of simpler optimization problems instead of the optimization problem with the coupled constraint.Thereby we use the free boundary as a control, which is an idea already presented by Volkov [3,12,13].
This paper is organized as follows.In the remainder of this section we formulate the model FBP on a reference domain so that we do not have to deal with shape derivatives.We also state the nonlinear PDE constrained OCP and the relaxed OCP.Section 2 is devoted to the analysis of the well-posedness of the relaxed OCP and the formal derivation of the adjoint systems.The numerical solution of the OCPs based on a Broyden-Fletcher-Goldfarb-Shanno (BFGS) method is covered in Section 3. Here, we also show that the sequence of minimizers converges towards a minimizer of the original optimization problem as desired.We conclude in Section 4 with a summary and outlook.

Coupled state problem
We consider a model FBP on a domain Ω  ⊂ ℝ 2 similar to the problem analyzed by Antil et al. [10].The difference is that we do distributed control, whereas their control variable  acts on the boundary.The top boundary Γ  of the physical domain Ω  is free and parametrized by , compare Figure 1.Throughout this work we require that  is Lipschitz continuous with constant 1, that is, |d  1 | ≤ 1 a.e. = (0, 1), which ensures that || ≤ 1∕2.The governing equations are given by a Poisson equation in the bulk and a Young-Laplace equation accounting for surface tension on the free boundary.In order not to deal with changes in the domain when we consider OCPs for the FBP, we transform Ω  to Ω = (0, 1) 2 and Γ  to Γ = (0, 1) × {1}, compare with ref. [10].Consequently, the coefficients in our PDE are highly nonlinear.Further, we identify the function spaces on Γ and  = (0, 1).Since we are only interested in small perturbations of the flat case  = 0, we assume like in ref. [10] that the curvature can be linearized and our control is scaled accordingly.This results in the following transformed PDE problem where  > 0 is some fixed constant,  is the outward pointing unit normal vector along Ω and  is the Nemytskii operator defined by

Graph as state variable
Our goal is to find an optimal control  ∈   ⊆  in the set of admissible controls, so that the solution pair (, ) ∈  ×  of the FBP approximates a given boundary   ∶  → ℝ and potential   ∶ Ω → ℝ.Hence, we consider the following tracking-type cost functional for fixed parameters   ,   ,   ≥ 0 weighting the objective and the control costs, respectively.Note, that √ 1 +  is the square root of the functional determinant of the transformation of Ω  to Ω. Treating the graph as a state variable the minimization problem reads min ∈   0 (, , ), where   is a set of admissible controls.

Graph as control variable
The PDE constraint in ( 0 ) is a highly nonlinear and numerically challenging due to the coupling with the free boundary.
To simplify the problem we incorporate the Equation (1c) for  into the cost function; see also refs.[3,12].By this relaxation the graph becomes an extra control variable.We introduce a parameter  > 0, relax the residual of the equation for  on the top boundary Γ and solve the resulting minimization problem.The relaxation is done in an  2 -sense, which yields with state variable  and control variables  and .Here,  and  are the function spaces from above with suitable admissible subsets   and   specified later.Note, that the state constraint in ( 0 ) becomes a control constraint, which is contained in   .

ANALYSIS
In the following we assume (A1) the set of admissible controls   or   ×   , respectively, is a weakly closed and convex subset of the control space, (A2) the desired potential   is in  2 (Ω), and the desired boundary   in  3 () ∩  1 0 ().
It is, for example, reasonable to consider the unconstrained case   =  and a small ball around  = 0 for   ⊆ .This choice of   is consistent with the small perturbations assumption and linearization of the curvature in (1c).

Well-posedness of (𝑷 𝝀 )
Due to the highly nonlinear coefficients in the PDE constraint of (  ), the following theorem provides the necessary existence and uniqueness of a solution .
The proof follows standard arguments; see, for example [17].

Formal Lagrange principle
Due to limited space, we only present the formal adjoint equations to ( 0 ) and (  ).Assuming enough regularity to ensure Fréchet differentiability, we apply the Lagrangian multiplier technique, where  denotes the Lagrange multiplier for the equation constraint and  is the corresponding one for the additional boundary condition; see, for example [18] or [17].
The adjoint equations to ( 0 ) in strong form can be derived analogously as in ref. [10], where also the definition of  1 [] and  2 [] can be found.They read and The optimality condition of the Lagrange function  0 for problem ( 0 ) is simply ∇   0 () = −(   + ).
The adjoint equation to (  ) is Consequently, to identify the descent direction, we need to compute an  3 ()-representative  of the derivative by solving the weak form of a 6ℎ order equation.The optimality condition of the Lagrange function   of (  ) with respect to  is as before, so the complete one is given by ) .

Setting
To solve the optimization problems ( 0 ) and (  ) numerically, we use the BFGS algorithm together with the Armijo rule [18].We choose   and   as the solution pair of (1) for ) inspired by a two-dimensional scaled density of the normal distribution with  = 0.5,  = 1,  = 150, and  = 20.The derivative in  2 direction of this function  is quite large, which causes an actual deformation of the top boundary and gives an interesting optimization setting.It further ensures that   and   are reachable.In the following the surface tension coefficient is fixed to  = 1.Recall that the cost function consists of two tracking type terms for   and   and a regularization term for .To get enough regularity of  in the analysis we needed that the term  −   is measured in the  3 ()-norm.However, for the numerical analysis we reduce these regularity assumptions and measure the difference only in  2 ().Since the case, where  acts as a state variable, automatically yields a solution with  2 -regularity, we want to ensure the same for the case where  acts as a control variable.This is realized by computing suitable gradients and details are given below in (6).The parameters in the cost functional (2) are chosen as   = 1.0,   = 100 and   = 10 −2 .In the relaxed problem (  ) we initially choose  = 10 −2 .The asymptotic solution for  → ∞ is computed via a homotopy approach.In our implementation  is increased by a fixed factor   = 10, until the residual of the equation  1 in the cost function of (  ) is smaller than a given tolerance chosen as   = 10 −6 .The critical point found for the maximal  is compared to the critical point computed by ( 0 ).
The PDE constraints are solved with a finite element method (FEM)in the framework of FENICS [19,20].For this purpose, the function spaces  and  are discretized with continuous galerkin (CG) elements of first order.The function space  is discretized with CG elements of second order [20].
The optimization procedure is initialized with  = 0 and  = 0 and stopped when the termination criterion is fulfilled, which is composed of two conditions.Either the relative gradient norm is smaller than a given tolerance, or the relative change of the controls is no longer large.This second criterion has been shown to be reasonable by numerical experiments.The tolerances for both criteria are chosen as  = 10 −3 .
Due to the discretization of  with second order CG elements, we have to be careful when implementing the  2 -scalar product.To compute a gradient representant in  2 () numerically we need to solve the following equation with , ,  as before.The notation ⟨⋅, ⋅⟩  2 (), denotes the numerical implementation of the  2 -scalar product.Since we deal with CG elements, the derivative is a piecewise linear function by convention with possible discontinuities.Therefore we have to implement terms that give a continuous solution without any jumps at the interface of two elements.This is done with a Nitsche approach, for which we need the following two notations.Let  be a piecewise continuous function on a mesh  .On a face  between two cells  1 and  2 , we define the average operator and the jump operator as respectively, where   are the outer normal vectors of the two cells.The numerical implementation of the  2 -scalar product then reads where the differential d  denotes the integral along the inner boundaries of all elements.More details can be found in ref. [21].

Results
Repeatedly solving problem (  ) for increasing values of  eventually leads to a critical point, where the second equation of the coupled PDE problem (1) is also satisfied in the weak sense.The evolution of the cost functional can be seen in Figure 2. The jumps result from the term  1 in the cost function and enlarging .Between two jumps we see the evolution of the cost functional while computing a critical point of (  ) for fixed .Next to the composite cost functional, one sees the decrease in the term  1 , which essentially measures the square of the residual of the second equation of the coupled PDE.
The relation between the increase in  and the decrease in  1 is somewhat better than linear, what we would expect.Note, however, that in ref. [22] it is shown for a different problem that the order of convergence can be better for solutions with higher regularity, which we might observe here as well.In our setting  was increased from 10 −2 to 10 3 .
It remains to illustrate that the final minimizer of (  ) converges towards the minimizer ( 0 ) for finer discretizations.For this purpose, the difference of the minimizers are calculated for ,  and  and plotted against the parameter  for different grid sizes in Figure 3.
In Figure 4 different (interim) solutions for  are shown, including the final critical point of (  ), and various intermediate solutions for  from 10 −2 to 10 3 .Next to it a zoom into the interesting part at the center of the interval is given.In both figures the critical point of ( 0 ) is shown as well.

CONCLUSION
We presented a highly nonlinear optimization problem for a FBP model that was transformed onto a reference domain and used a penalty method to relax the condition for the free boundary.First, the free boundary was treated as a state variable, while in the second approach it acts as a control variable.We stated results concerning the existence and uniqueness of solutions with  2 (Ω) regularity to the state equation in the relaxed case.Further, we discussed the well-posedness of the OCP for the relaxed problem, which required a careful regularity analysis of the co-normal derivative involving the highly nonlinear Nemytskii operator [].For numerical simulations, we formally calculated derivatives using the Lagrangian approach.Finally, we examined numerically the asymptotic behavior of the two OCPs for a reachable setting, underlining the feasibility of our approach.

A C K N O W L E D G M E N T S
Open access funding enabled and organized by Projekt DEAL.

2 3 4
Optimization results of control problem.(A) Cost functional with increased  at jumps.(B) Corresponding mismatch  1 after each outer iteration with fixed .Relative difference between critical point of state problem and critical points for ,  and  as  increases in regularized case.Asymptotic behavior of  for ℎ = 1 64 .(A) Evolution of critical points of (  ) for  from 10 −2 to 10 3 with critical point  * of problem ( 0 ).(B) Zoom into left figure.