CASHOCS: A Computational, Adjoint-Based Shape Optimization and Optimal Control Software

The solution of optimization problems constrained by partial differential equations (PDEs) plays an important role in many areas of science and industry. In this work we present a new software package for Python, called cashocs, which automatically solves such problems in the context of optimal control and shape optimization. Our software implements a discretization of the continuous adjoint approach, which derives the necessary adjoint systems and (shape) derivatives in an automated fashion. As cashocs is based on the finite element software FEniCS, it inherits its simple, high-level user interface. This makes it straightforward to define and solve PDE constrained optimization problems with our software. In this paper we discuss the design and functionalities of cashocs and also demonstrate its straightforward usability and applicability.


Motivation and significance
Shape optimization and optimal control problems constrained by partial differential equations (PDEs) and their numerical solution are important in many areas of science and industry: They are, e.g., used for the optimization of chemical reactors [1], glass cooling processes [2], and semiconductors [3] as well as the optimal design of cooling systems [4], aircrafts [5], and electric machines [6].To solve these problems, the so-called adjoint approach is often employed, which facilitates the computation of (shape) gradients for the problems, which can be used to solve them numerically.However, for complex, coupled, or highly nonlinear problems, such as the ones arising from industrial applications, even the derivation of the necessary equations for the adjoint approach is an extremely involved and error-prone task.Consequently, it is not feasible to carry it out by hand anymore (see, e.g., [7]).For these reasons, there has been a lot of effort recently to automate the tasks for solving PDE constrained optimization problems, resulting in software such as dolfin-adjoint [8], Fireshape [9], NGSolve [10], and our software cashocs.What distinguishes cashocs from these other packages is its novel approach of using automatic differentiation solely to derive the adjoint system and (shape) derivatives, while implementing and automating a discretization of the continuous adjoint approach in all remaining aspects.Particularly, it implements discretizations of continuous, infinite-dimensional optimization algorithms, whereas the other packages use either external optimization libraries for finite-dimensional optimization, or require the user to implement these algorithms themselves.Our approach leads to unique features, such as the possibility of discretizing and solving the state and adjoint systems differently as well as the choice of the scalar product for the computation of the (shape) gradients, and also gives rise to mesh independent behavior, as shown in Section 3.Moreover, cashocs is the only one of these packages that has implemented a remeshing feature for shape optimization problems.

Mathematical Background
Let us begin with stating the general form of the optimization problems our software can solve.Optimal control problems have the form where y is again the state variable, and the set of admissible domains A is used to incorporate additional geometrical constraints.We interpret the PDE constraint e(y, Ω) = 0 in the following weak sense Find y ∈ Y (Ω) such that e(y, Ω), p Z(Ω) * ,Z(Ω) = 0 for all p ∈ Z(Ω).
In particular, this means that the PDE constraint is given on the domain Ω, and it is the shape of this domain that is subjected to optimization.Problems (1) and (2) are prototypes for the kinds of problems that cashocs can solve, and we refer the reader to Section 3 for illustrative examples.As mentioned previously, these kinds of problems are usually solved with the adjoint approach, whose derivation is beyond the scope of this paper.Hence, we refer the reader to [11,12] and [13,14,15] for a discussion and derivation of the adjoint approach for optimal control and shape optimization problems, respectively.

Software Architecture
To solve optimization problems with cashocs, the user has to do the following.First, they have to implement the problem in a Python script, including the definition of the computational mesh, the state system, and the cost functional.To do so, they can use the same syntax as for defining the problem in FEniCS, with only very minor modifications, resulting in a simple, high-level user interface that supports many important types of optimization problems (cf.Section 3).Second, the user has to define a configuration file that specifies the parameters for the solution of the state system and the optimization algorithm, which is loaded into the user script.Then, one can set up an optimization problem using cashocs.OptimalControlProblem or cashocs.ShapeOptimizationProblem, respectively, and solve it with the solve method of the respective class.Internally, our software uses FEniCS to generate and compile C++ code for the finite element assembly of the problems and PETSc [16] is used to solve the arising linear systems, which makes the solution of the problems very efficient.A schematic overview of cashocs' architecture can be seen in Figure 1.

Software Functionalities
The following algorithms are available for shape optimization and optimal control problems in cashocs • the gradient descent method, • nonlinear conjugate gradient methods (NCG) methods, • limited memory BFGS (L-BFGS) methods.
Note, that particularly for shape optimization these algorithms correspond to the state of the art, with the L-BFGS methods being introduced in [14], and the NCG methods in [15].Additionally, the following optimization algorithms are available for optimal control problems only • a truncated Newton method, • a primal-dual active set method.
Note, that for optimal control problems, all methods can also treat box constraints for the control variable using projection techniques.The user can adjust the behavior of these algorithms using the configuration file, where, e.g., the relative and absolute stopping tolerances, maximum number of iterations, and other, algorithm specific, parameters can be modified.
Additional features of cashocs include, among others, the possibility to use different discretizations for state and adjoint systems, the implementation of a Picard iteration for solving coupled systems, the possibility to specify which scalar product is used for the computation of the (shape) gradient, and remeshing for shape optimization problems, which utilizes the mesh generation software Gmsh [17].

Illustrative Examples
To demonstrate our software's simplicity for defining PDE constrained optimization problems as well as its efficiency for solving them, we now investigate two model problems, one for optimal control and one for shape optimization.Note, that a variety of other examples for using cashocs can be found in the tutorial at https://cashocs.readthedocs.io/en/latest/tutorial_index.html.

Optimal Control
As a model optimal control problem we consider the following one from [11] min This optimal control problem has a tracking-type cost functional with a Tikhonov regularization for the control variable.The PDE constraint is given by a Poisson equation with homogeneous Dirichlet boundary conditions, and the control variable u enters the PDE as a right-hand side.The weak formulation of this PDE constraint is given by (4) For this example, let us use Ω = (0, 1) 2 , α = 1e−4, and For the discretization of the domain we use a uniform triangular mesh which divides Ω into n × n squares that are halved to create triangles.To solve this problem with cashocs, we can use the code shown in Listing 1, which we briefly discuss in the following.Note, that as our software is based on FEniCS, we refer the reader to [18, Chapter 1], where the syntax of FEniCS is explained using several descriptive examples.In Listing 1, we begin by importing FEniCS and cashocs in lines 1 and 2. Next, we define the mesh with the UnitSquareMesh function, and set up the volume measure for integration, in lines 5-7.Subsequently, we define a function space of linear Lagrange elements in line 10, and define the functions y, p, and u.These are used to define the weak form of the PDE constraint in line 17, where the function p plays the role of the test function.In lines 19 and 20 the Dirichlet boundary conditions for the Poisson problem are defined.Note, that up until now, we only used commands from FEniCS with the following minor variations.Instead of defining y as a TrialFunction and p as a TestFunction, both are now Function objects.Additionally, instead of defining the (linear) PDE constraint using its left-and right-hand sides, we define it as we would for a nonlinear variational problem in FEniCS, anal-   ogously to the form in ( 1) and ( 4).In lines 23 and 24 we define the desired state, which is used in line 25 to define the cost functional.Again, we have only used FEniCS commands for these operations.To invoke cashocs to solve this problem, all we have to do is loading the configuration file into the script in line 28, initializing the OptimalControlProblem in line 29, and calling its solve method subsequently.In total, we have to add only three additional lines of code to solve the problem.Note, that a minimal configuration file for the code is shown in Listing 2, and for a detailed description of the configuration files for optimal control problems we refer to https://cashocs.readthedocs.io/en/latest/demos/optimal_control/doc_config.html.
A plot of the computed optimal control and state using the Dai-Yuan nonlinear CG method is shown in Figure 2.Moreover, Table 2 shows the amount of iterations the optimization algorithms need to solve this problem for a sequence of finer meshes, using n = 16, 32, 64, and 128 subdivisions.We observe that all algorithms show mesh independent behavior as they basically need the same amount of iterations for convergence regardless of the discretization.

Shape Optimization
As model problem for shape optimization we consider the following one from [15,19] ( For this problem, the PDE constraint is, again, given by a Poisson problem with homogeneous Dirichlet boundary conditions, so that its weak form is given by ( 4) with u replaced by f .We proceed analogously to [15,19] and use as initial guess for the domain Ω 0 the unit circle in R 2 , and for the right-hand side f we use We discretize Ω 0 with a uniform triangular mesh by dividing the circle into n smaller strips, which are then meshed uniformly.This problem can be solved with cashocs using the code provided in Listing 3, which we briefly discuss in the following.As before, we refer to [18, Chapter 1] for a detailed introduction to the syntax of FEniCS, which we also use for the problem definition in cashocs.
from fenics import * import cashocs # define mesh and volume measure n = 64 mesh = UnitDiscMesh .create (MPI.comm_world , n, 1, 2) dx = Measure ('dx ', mesh) # function space of linear Lagrange elements V = FunctionSpace (mesh , 'CG ', 1) # state and adjoint variables u = Function (V) p = Function (V) # right−hand side x = SpatialCoordinate (mesh) The code is very similar to the one in Listing 1 as we again have a Poisson equation as PDE constraint.We start the script by importing FEniCS and cashocs.Then, we define the mesh and volume measure, now using the function UnitDiscMesh, in lines 5-7.For the discretization of the Poisson equation, we again use linear Lagrange elements whose corresponding function space is defined in line 10, and the functions y and p are defined in Table 3: Required iterations to reach the stopping criterion for problem (5).
lines 12 and 13.Thereafter, we define the right-hand side of the Poisson problem, using SpatialCoordinate in lines 16-18, which is then used to define the weak form of the Poisson equation in line 20.As for optimal control problems, the only major differences to traditional FEniCS syntax are that y and p are Function objects, and that the PDE constraint is written in the sense of ( 2) and ( 4).Subsequently, we set up a FEniCS MeshFunction for the boundaries, which is used to define the Dirichlet boundary conditions.Moreover, this is used to define which boundaries are fixed via the configuration file (cf.lines 7-8 of Listing 4).Finally, we define the cost functional in line 28.For solving this problem with cashocs, we proceed analogously to Listing 1, and first load the configuration file, then set up the ShapeOptimizationProblem, and finally call its solve method in lines 31-34.Note, that a minimal configuration file for this problem is shown in Listing 4. A detailed discussion of the configuration files for shape optimization can be found at https://cashocs.readthedocs.io/en/latest/demos/shape_optimization/doc_config.html.A plot of the optimal state on the optimal domain, computed with the Dai-Yuan NCG method in cashocs, is given in Figure 3.Moreover, we also show the number of iterations required by the algorithms on successively finer discretizations of n = 16, 32, 64, and 128 strips for the unit circle in Table 3.As before, we see that the number of iterations basically stays the same regardless of the discretization, which shows that we also have mesh independent behavior for shape optimization problems.

Impact
Our software enables users to treat complex, coupled, and highly nonlinear PDE constrained optimization problems in an automated fashion.The user is only required to define the PDE constraint and cost functional using basically the same syntax as for defining these objects in FEniCS.Thanks to the high-level user interface, the corresponding optimization problem can then be solved by adding only three additional lines of code.Our approach of implementing a discretization of the continuous adjoint approach leads to mesh independent behavior of the optimization algorithms, as shown in Section 3, making our software attractive for science and industry.In fact, cashocs has already been used to treat highly nonlinear optimization problems for parameter identification and optimal control in the context of chemical microreactors in [1].It has also been used in [15] for a numerical benchmark of NCG methods for shape optimization.Moreover, cashocs is used at Fraunhofer ITWM to solve PDE constrained optimization problems for industrial applications.Due to the generality of our software, which can treat lots of important classes of cost functionals and PDE constraints, it can be applied to many relevant problems in science and industry, automating their solution in an efficient and user friendly way.

Conclusions
We have presented cashocs, a software for numerically solving PDE constrained shape optimization and optimal control problems.The software automatically derives the required adjoint systems and (shape) derivatives, and implements a discretization of the continuous adjoint approach.Our software inherits FEniCS' high-level user interface which allows for a straightforward definition and solution of PDE constrained optimization problems.Additionally, the user still retains control over many important parameters for the optimization, ranging from the solution of the PDEs to the optimization algorithm, which allows them to make precise adjustments to the numerical solution of their problems.

Conflict of Interest
We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.
where u ∈ U and y ∈ Y are the control and state variables, U and Y are appropriate Banach spaces, and the set of admissible controls U ad ⊆ U is used to model additional constraints on the control variable.Moreover, J : Y ×U → R is the cost functional and e : Y ×U → Z * is a PDE constraint, which we interpret in the following weak sense Find y ∈ Y such that e(y, u), p Z * ,Z = 0 for all p ∈ Z.Here, Z * denotes the topological dual space of Z, and ϕ, x Z * ,Z denotes the duality pairing of ϕ ∈ Z * and x ∈ Z. Shape optimization problems have the form min y,Ω J (y, Ω) s.t.e(y, Ω) = 0 and Ω ∈ A,
(a) Optimal control u.(b) Optimal state y.

Figure 3 :
Figure 3: Numerical solution of problem (5), optimal state y on the optimal domain Ω, computed with the Dai-Yuan NCG method.

Table 2 :
Required iterations to reach the stopping criterion for problem (3).