## 1 Introduction

A striking property of nonlinear partial differential equations (PDEs) is that they may support multiple solutions, not trivially related via the nullspace of a linear operator. These solutions are captured in a bifurcation diagram, recording the solutions to an equation

(1) |

as a parameter is varied, where is a Banach space on a bounded domain , , and . When a system is to be designed, consideration of multiple solutions can be of primary importance: for example, one may wish to ensure that an aircraft remains in a high-lift rather than low-lift regime in takeoff [kamenetskiy2014], or may wish to exploit bistability in the design of a liquid crystal display that only consumes power when an image is changed [jones2012].

In this work, we consider the mathematical *control* of the
bifurcation structure of a physical system. Both bifurcation analysis
and PDE-constrained optimization are now mature fields, but their
combination remains unexplored. As a first step in this direction, we
formulate an algorithm that modifies the domain on which a PDE
is posed in order to advance or delay a single specified bifurcation
point. Several engineering problems may be understood as such a task.
For instance, in wing design one typically aims to avoid instabilities
due to the buckling of the wing structure. In this example, the
bifurcation parameter is the scaling of the load applied to
the wing. To avoid buckling, one can modify the design of the structure
to delay the first bifurcation point (with respect to
), and thus remove potentially unstable solutions. A similar
approach might be taken to the design of aircraft
stiffeners [xia2020nonlinear] and other components. On the
contrary, in other applications, inducing new branches might be
desirable. For instance, this is the case in the design of new
mechanical metamaterials that exploit instabilities to rapidly switch
between configurations, as in the snapping of a Venus
flytrap [medina2020navigating, vidoli2008, forterre2005].

We formulate the task of controlling an isolated bifurcation point as a PDE-constrained shape optimization problem. The essential idea of the formulation is to characterize the bifurcation point as the solution of an augmented system of equations. There are many choices for how to describe bifurcation points in this way; see e.g.[seydel2009practical, §5.4]. In our case we choose the system proposed by Seydel [seydel1979] and Moore & Spence [moore1980calculation], referred to as the Moore–Spence system. We then propose a numerical algorithm for solving the constrained optimization problem. This is a rather challenging task: none of the state equation Eq. 1, the Moore–Spence system, or the outer shape optimization problem have unique solutions.

Much of the literature on PDE-constrained optimization requires that the control-to-state map is single-valued [troltzsch2010, hinze2009, Al07, DeZo11, HaMae03, SoZo92]. In fact, in many analyses some smoothness of the map is further required. However, the case where the implicit function theorem cannot be applied has also been considered, in both finite and infinite dimensions [heinkenschloss2014, dennis1997, heinkenschloss2002]. Some experiments involving multiple solutions of optimal control problems, arising from multiple solutions of the underlying control-to-state map, were recently presented in [pichi2020driving]. In our work, we address the challenges posed by the co-existence of multiple solutions by careful selection of the initial guesses employed and their updates, which occur while the optimization is carried out, to exploit local uniqueness of isolated bifurcation points. Our algorithm proves successful in several numerical experiments, and allows us to find a new shape for which a branch of solutions arises at a specific value .

While this paper focuses on controlling simple pitchfork and fold bifurcation points with respect to the shape of the domain, we expect that the ideas developed here also apply to other settings, such as the control of the bifurcation diagram with respect to a parameter as another parameter in some (possibly infinite-dimensional) parameter space is varied, or the control of other kinds of bifurcations such as Hopf bifurcations.

The paper is organized as follows. We first review the characterization we use of bifurcation points using augmented systems of partial differential equations in Section 2. Then, in Section 3, we describe an algorithm for performing PDE-constrained shape optimization constrained by the augmented system. In Section 4, we combine these numerical techniques in an algorithm that modifies a shape to delay or advance a given bifurcation point. Finally, in Section 5, we apply this technique to delay the bifurcations of three numerical examples involving the Allen–Cahn equation, the Navier–Stokes equations, and the compressible hyperelasticity equation.

## 2 Characterization of singular points

Bifurcation diagrams of problems of the form may contain bifurcation points (see [seydel2009practical]), which we define as follows.

###### Definition 1 (Bifurcation point)

Assume that the operator does not have a symmetry group. For , define and, for and , let . A bifurcation point (also called branch point) with respect to is a solution pair to such that, for any , the cardinality of changes when passes .

###### Remark 1

The nonlinear operator may have a symmetry group such that for all group elements and , . Here represents the action of the element on the solution . Such a group is called a Lie group when it has a structure of a smooth manifold and when the group multiplication and inversion operations are smooth maps [olver2000applications]. In the presence of such symmetries, one is typically interested in points at which the number of distinct group orbits varies, which are sets of the form

We illustrate some different situations with bifurcation points arising in the bifurcation diagram of a nonlinear PDE in Fig. 1. Fig. 1

(a) features a fold bifurcation with a turning point. Panels (b) and (c) illustrate two potential scenarios for the intersection of two branches. In a transcritical bifurcation (panel (b)), the trivial branch (indicated by a black line) may lose stability when one of the state eigenvalues crosses zero at the bifurcation point, which results in an exchange of stability to the other branch, indicated by a red line. A subcritical pitchfork bifurcation (c) might also feature a loss of stability at the bifurcation point, where the trivial branch gives birth to a secondary stable branch (indicated in red).

According to Definition 1, bifurcation points are characterized by the change in the number of solutions as varies in a small interval around . This implies that the Fréchet derivative with respect to the variable at the bifurcation point cannot be an isomorphism (otherwise, by the implicit function theorem, there would exist locally a unique solution curve , but this is not compatible with Definition 1). Therefore, we are interested in the computation of singular points, i.e., solution pairs to Eq. 1 such that does not have a bounded inverse.

The idea of Seydel and Moore & Spence for finding singular solutions to Eq. 1 is to solve the following extended system (2): find such that equationparentequation

(2a) | ||||

(2b) | ||||

(2c) |

where is a functional chosen to normalize . Equations (2b) and (2c) imply that

is a non-zero eigenfunction associated to the zero eigenvalue of

.###### Remark 2

Several choices of normalization conditions in the Moore–Spence system have been considered [cliffe1986use, winters1988onset]. In this work we set

###### Remark 3

In this work, we only consider simple singular points, those for which .
This includes nondegenerate turning points and simple symmetry-breaking pitchfork bifurcations [moore1980calculation, werner1984computation]. We further assume that the solutions
of (2) are isolated.
Bifurcation points of higher multiplicity () also satisfy (2) *a fortiori*,
but may not be isolated.
It is possible to design augmented systems that are selective, in the sense that they admit certain kinds of bifurcation points as solutions but not others
(e.g. [moore1980numerical, weber1981numerical]). These typically have the advantage that the solutions are isolated, but the augmented systems become larger and more complicated.

###### Remark 4

Different generalizations have also been studied to compute Hopf bifurcations [strogatz2018nonlinear], i.e. points at which a periodic solution arises, in contrast to the stationary solutions considered in this work [griewank1983calculation, jackson1987finite, jepson1981numerical, roose1985direct].

Since we are interested in changing the shape of the domain to control bifurcation diagrams, we now present a shape optimization technique constrained by partial differential equations.

## 3 PDE-constrained shape optimization

In this section, we give a brief introduction to PDE-constrained shape optimization. A more detailed exposition of the mathematical and implementation techniques employed here is given in [paganini2018higher, paganini2020fireshape].

A shape optimization problem is an optimization problem of the following form: find a domain that minimizes a shape functional , where denotes a set of admissible domains. If evaluating the shape functional on a domain requires solving a PDE stated on , the shape optimization problem is said to be PDE-constrained.

There are competing approaches to define the set of admissible domains , such as the phase-field method [GaHeHiKa15, BlGaFaHaSt14], the level-set method [La18, AlJoTo02] and the moving mesh method [AlPa06, paganini2018higher], among others. In this work, we consider the moving mesh method, because its neat integration with standard finite element software allows the automation of several key steps in shape optimization [paganini2020fireshape]. The main idea of the moving mesh method is to construct the set of admissible domains by applying diffeomorphisms to an initial domain , that is,

where is the group of bi-Lipschitz isomorphisms on , that is and are Lipschitz continuous. In this setting, shape optimization is carried out by constructing a geometric transformation so that minimizes . This approach is called moving mesh method because, in its lowest order discretization, this method is equivalent to replacing the initial domain with a finite element mesh and optimizing the coordinates of the mesh nodes while retaining the mesh connectivity [paganini2018higher].

The transformation is commonly constructed iteratively. First, one initializes with the identity transformation, that is, for every . Then one constructs the subsequent iterates using the composition formula

(3) |

where is the geometry update. These geometry updates can be computed in a standard fashion using derivatives of the function . As an example, in the next paragraph we present a simple steepest descent shape optimization algorithm. More advanced optimization strategies can be adapted to shape optimization in a similar way.

To construct the sequence using a steepest descent algorithm, we compute the update as the Riesz representative of the shape derivative evaluated at [paganini2018higher, PaHi16]. This means that satisfies the equation

(4) |

where is a suitable Hilbert space with inner product , and

(5) |

Once has been computed, the new transformation is obtained using formula (3), possibly scaling by a sufficiently small step length to ensure that is bijective [Al07, Lemma 6.13].

###### Remark 5

There are competing requirements on the choice of the Hilbert space [paganini2018higher]. On the one hand, it is desirable that the Riesz representative is sufficiently regular to ensure that is bi-Lipschitz. This suggest using an -type space, which is embedded in [Ev10, Sect. 5.6.3]. On the other hand, it is desirable that (4) can be solved efficiently using standard finite elements, which suggest using an -type space. In this work, we follow the latter option and employ an inner product based on the linear elasticity equations. This is a standard choice in shape optimization and, although it is not guaranteed that at the continuous level, its finite element approximation (based on standard Lagrangian finite elements) belongs to , and is thus suitable to define the transformation .

In our numerical simulations, we employ the shape optimization software Fireshape [paganini2020fireshape], which is based on the moving mesh method. This choice offers numerous convenient features, including automated shape differentiation in the Unified Form Language [ham2019automated], automated adjoint computation via pyadjoint [farrell2013automated, dokken2020automatic], integration with the finite element software Firedrake [rathgeber2016firedrake] for efficient linear solvers (e.g. block preconditioners and geometric multigrid) from PETSc [balay2020petsc], and access to numerous optimization algorithms via the Rapid Optimization Library [ridzal2017rapid].

## 4 Optimization of isolated bifurcation points

In light of Section 2 and Section 3, we formulate the problem of finding a shape of the domain such that the operator admits a bifurcation point at a target bifurcation parameter as follows:

(6) |

Compared to classical PDE-constrained shape optimization, the problem (6) presents an additional difficulty: the state constraint admits multiple solutions by design. More precisely, unless the bifurcation diagram of is trivial, there are values of for which the problem to find such that admits multiple solutions . Additionally, the bifurcation diagram of typically includes multiple bifurcation points (in the sense of Definition 1). For these reasons, to solve (6), we must select not only an initial domain , but also an initial isolated bifurcation point that we wish to control. This selection can be performed by first computing the bifurcation diagram of the initial domain and selecting via visual inspection. Since most algorithms for bifurcation analysis compute bifurcation points approximately, we then refine the approximation of using the following procedure. We select a pair from the computed bifurcation diagram that both satisfies and that is sufficiently close . Then, we solve the following eigenvalue problem:

(7) |

to determine the first (a typical value used is ) smallest (single) eigenvalues in magnitude and corresponding eigenfunctions (normalized to with respect to the norm ). The eigenvalue problem Eq. 7 is solved using a Krylov–Schur algorithm with a shift-and-invert spectral transformation [stewart2002krylov], implemented in the SLEPc library [hernandez2005slepc]. Finally, we form initial guesses and solve the Moore–Spence system with Newton’s method. This generates a set of (at most) distinct solutions , from which we finally select the tuple minimizing the difference between and in the -norm as the starting point for shape optimization, that is, as .

###### Remark 6

Using this procedure ensures that the initial state is feasible. In practice, selecting is important, because situations like the one illustrated in Fig. 2 might arise; the eigenfunction associated with the lowest eigenvalue for may not be a good guess for the lowest eigenfunction at , since the eigenvalues may swap order.

Once and have been determined, we use Fireshape [paganini2020fireshape] to solve (6) iteratively with the trust-region optimization algorithm [conn2000trust] and inner solver L-BFGS [liu1989limited] implemented in the ROL optimization library [ridzal2017rapid]. To ensure that the computational mesh is not tangled, we reject optimization steps for which the minimum of the determinant of the Jacobian of the transformation defined in (3) becomes negative. Additionally, we need to ensure that, after a domain update, Newton’s method converges to a solution to the Moore–Spence system that lives on the same branch as the one computed at the previous iteration, i.e. is close enough to a solution obtained before. To enforce this, we reject domain updates that do not satisfy the inequality

(8) |

where denotes the pushforward operator. For example, if , then .

## 5 Numerical examples

In this section, we apply the shape optimization strategy described in Section 4 to three test cases based on three different differential equations : the Allen–Cahn equation, the Navier–Stokes equations, and the hyperelasticity equations. The bifurcation diagrams displayed in this section are computed with deflated continuation [farrell2016computation, farrell2015deflation], complemented with arclength continuation. These techniques have been successfully applied to a wide range of physical problems such as the deformation of a hyperelastic beam [farrell2018a], liquid crystals [emerson2018computing, xia2021a], Bose–Einstein condensates [boulle2020deflation, charalampidis2020bifurcation, charalampidis2018computing], and fluid dynamics [boulle2021bifurcation]. However, our optimization strategy is not tied to a specific algorithm for bifurcation analysis, and other options may be used.

### 5.1 Allen–Cahn equation

As first test case, we consider the cubic-quintic Allen–Cahn equation, which can be used to model phase separation [allen1979microscopic, uecker2014pde2path]. In strong form, the Allen–Cahn equation is given by

(9) |

with homogeneous Dirichlet conditions. In this example, we consider two-dimensional domains . Therefore, we can set .

We consider two different initial domains: a square with rounded corners and edge length 2 (see Fig. 5(a)), and a unit disk . The computational meshes of and are both generated using Gmsh [geuzaine2009gmsh] and are composed of and triangles, respectively. The Allen–Cahn equation is discretized using piecewise linear continuous Lagrangian finite elements.

In (9), several bifurcations occur as varies. In the left panels of Figs. 4 and 3, we show the bifurcation diagrams computed with deflated continuation [farrell2016computation] for when the computational domain is and , respectively. We notice that several branches arise in both cases from the ground state . We can order these branches by numbering them from the left. We observe that the solutions associated to higher branches feature more complicated patterns than the states belonging to the first branches (see Figs. 6 and 5). Another interesting observation is the presence a fold bifurcation on each branch, due to the cubic-quintic terms [uecker2014pde2path].

We now apply the optimization method described in Section 4 to delay the birth of the first branch from the ground state in the Allen–Cahn equation to a target bifurcation parameter of . In Figs. 6 and 5, we display the evolution of the optimization functional with respect to the trust-region iterations. We also show some snapshots of the solution and domain for certain values of . We observe that the optimization algorithm reduces the value of the objective functional to less than (Figs. 6 and 5) in about 20 optimization steps.

To verify that the returned solutions satisfy the original objective, we recompute the bifurcation diagrams (using deflated and arclength continuation) on the optimized shapes. The right panels of Figs. 4 and 3 illustrate the resulting bifurcation diagrams using the optimized shapes, starting from and , respectively. In both cases, we observe that the bifurcation structure is controlled in the desired way, and that the first branch arises at the desirable value of the bifurcation parameter. As expected, the new shape also affects the birth of the later branches. For instance, the second branch of the diagram with originally arises at and gets delayed to in Fig. 4, while the higher branches are still not visible in the range .

This example also showcases the complexity and nonlinearity of the problem studied in this work. None of the problems at any level have a unique solution: the PDE may have multiple solutions for a given bifurcation parameter ; the Moore–Spence system has multiple solutions due to the presence of several bifurcation points; and the shape optimization problem admits multiple solutions, as found from different initial shapes. The shapes displayed in Fig. 5(a) and Fig. 6(a) are both solutions to the same shape optimization problem but initialized using and .

Finally, we demonstrate that the same approach can be employed to control arbitrary branches. We repeat the previous experiment and shape optimize to postpone the birth of the third branch (which arises near , see the left panel of Fig. 3) to . In Fig. 5(b), we observe that the shapes obtained by this procedure, the value of falls below machine precision in roughly 26 optimization steps. Similarly, in Fig. 6 (b) we show the optimization history when delaying the birth of the fifth bifurcation of from (roughly) to . In this case, the optimized shape resembles a (smoothened) hexagon with edges at the extrema of the function (indicated by blue and red colours).

### 5.2 Navier–Stokes equations

In the second test case, we consider the functional associated to the nondimensionalized steady incompressible Navier–Stokes equations. These read

equationparentequation

(10a) | ||||

(10b) |

where is the fluid velocity, , is the pressure, and Re is the Reynolds number, which in this example plays the role of bifurcation parameter. In particular, we consider an initial geometry formed by the union of two rectangles union of two rectangles,

, and represents a channel subject to sudden expansion. The inflow is modelled with a Poiseuille flow at the leftmost boundary inflow boundary condition, whereas an outflow boundary condition is imposed on the rightmost boundary. The remaining parts of the boundary carry no-slip boundary conditions. The velocity and pressure are discretized with the inf-sup stable Taylor–Hood finite element, i.e. the vector-valued piecewise quadratic continuous Lagrangian element for velocity, and piecewise linear continuous Lagrangian element for pressure.

It is well known that this problem admits multiple solutions depending on the values of the bifurcation parameter Re [fearn1990]. The bifurcation diagram for the initial mesh is shown in Fig. 7(a). The birth of the first bifurcation from the symmetric ground state occurs at . If the mesh were -symmetric about , the bifurcation diagram would be symmetric and connected, with branches originating at pitchfork bifurcations from the symmetric ground state. Since symmetry is broken by the use of the unstructured mesh, the pitchforks are perturbed into disconnected fold bifurcations. This phenomenon is very difficult to avoid in practice (especially on complex industrial geometries), but deflation deals with it robustly.

Using the shape optimization approach described in Section 4, we modify the shape of the inlet to delay the birth of the first bifurcation from the ground state solution. Since we are interested in changing the shape of the inlet, we fix all boundaries during the shape optimization algorithm except for the following segments:

Fig. 8(a) displays the velocity magnitude of the solution on the initial mesh at the bifurcation point located at . In Fig. 8(b)-(d), we report the shape of the inlet channel after delaying the bifurcation point to , respectively. We observe that the symmetry of the solution illustrated in Fig. 8(a) is broken as we delay the birth of the bifurcation. The shape of the inlet channel becomes bent and forces the fluid to flow downward as the value of the bifurcation point increases from to

. In each case, the shape optimization algorithm required between 20 and 30 steps to minimize the loss function to nearly machine precision.

Finally, we apply deflated continuation to compute the bifurcation diagram on the bent domain shown in Fig. 8(c), for which the first bifurcation point happens at . The resulting bifurcation diagram is displayed in Fig. 7(b). We observe that the first bifurcation branch effectively arises at through a fold bifurcation and is disconnected from the original solution branch.

### 5.3 Hyperelastic beam

We consider the shape optimization of a hyperelastic beam under compression. The potential energy of the beam is given by

(11) |

where is the initial domain, is the unknown displacement, is a constant body force per unit area, and is a Neo–Hookean strain energy density given by

(12) |

Here and denote the Lamé parameters (chosen to yield a Young’s modulus and Poisson ratio ), , , , and . Finally, we complement the potential energy (11) with the Dirichlet boundary conditions

where the variable load plays the role of bifurcation parameter. On the remaining part of the boundary we impose natural traction-free boundary conditions.

The presence of the bifurcation parameter in the Dirichlet boundary conditions presents an additional difficulty. Dirichlet boundary conditions such as this are typically enforced in the definition of the trial space, by seeking a solution in

where is the given boundary data and is the boundary on which to apply the conditions. However, this approach is not feasible when deriving the Moore–Spence system (2) symbolically using the automatic differentiation facilities of UFL [alnaes2014], as these require having the bifurcation parameter as a variable in the weak form of the equations, rather than in the spaces employed.

To circumvent this difficulty, we choose to seek a solution in

where , and to impose the Dirichlet boundary condition on the right boundary using Nitsche’s method [nitsche1971variationsprinzip, sime2020automatic] (as in Xia et al. [xia2020nonlinear]). More precisely, we add the following terms

(13) |

to the weak formulation arising from the Fréchet differentiation of the hyperelastic energy (11). In formula (13), is the Nitsche penalization parameter, denotes the outward normal vector, and

is the first Piola–Kirchhoff stress tensor, calculated as the derivative of the stored energy density function:

where . In the numerical experiments, we choose as in Xia et al. (justified in [xia2020nonlinear, Appendix A]). The final weak form is to find such that

(14) |

for all . We discretize Eq. 14 using piecewise linear continuous Lagrangian finite elements.

We observe several fold bifurcations as the bifurcation parameter varies in the interval on the initial domain. The bifurcation diagram of this problem is illustrated in Fig. 9(a) using the diagnostic (i.e. the second component of the displacement) to distinguish the different bifurcation branches. In this problem the -symmetry is broken by the nonzero gravitational body force , causing the structure to prefer to buckle downwards rather than upwards. We seek to delay the birth of the first fold bifurcation arising at to (see Fig. 9(a) and the corresponding displacement illustrated in Fig. 10(a)). As in Section 5.2, the shape optimization algorithm required between 20 and 30 steps to minimize the loss function to nearly machine precision in each case.

In Fig. 10(b)-(d), we report the shapes returned by
the shape optimization procedure and the magnitude of the displacement at the
bifurcation points (respectively arising at
). We observe that the shape of the beam is bent
downward, which makes it difficult to buckle and therefore delays the location of the first fold bifurcation. The bifurcation
diagrams of the hyperelasticity equations solved on these new shapes via
deflated continuation are displayed in Fig. 9(c)-(d). We
find that the bifurcation structure of the problem changes as the location of
the fold bifurcation is delayed and as the shape of the beam is bent. In
Fig. 9(d), we see that the branch that we chose to
control does indeed arise at as desired, but that another fold bifurcation
now exists for a smaller bifurcation parameter. This is because the
numerical technique described in this paper controls a *specific* branch, but
the other bifurcation points are also modified as the shape of the domain changes
during the optimization procedure.

## 6 Conclusions

We have presented a robust computational technique to control the location of a specified bifurcation point using shape optimization. The key idea is to describe the bifurcation point using the Moore–Spence system. We applied this procedure to vary bifurcation points in the bifurcation diagrams of the cubic-quintic Allen–Cahn equation, the incompressible Navier–Stokes equations, and a hyperelasticity equation. We expect that this numerical technique could have a wide range of applications in physics and engineering to control the stability of solutions in areas such as the design of new mechanical metamaterials, structural components, flow devices, and liquid crystals.

Comments

There are no comments yet.