Driving bifurcating parametrized nonlinear PDEs by optimal control strategies: application to Navier-Stokes equations with model order reduction

This work deals with optimal control problems as a strategy to drive bifurcating solution of nonlinear parametrized partial differential equations towards a desired branch. Indeed, for these governing equations, multiple solution configurations can arise from the same parametric instance. We thus aim at describing how optimal control allows to change the solution profile and the stability of state solution branches. First of all, a general framework for nonlinear optimal control problem is presented in order to reconstruct each branch of optimal solutions, discussing in detail the stability properties of the obtained controlled solutions. Then, we apply the proposed framework to several optimal control problems governed by bifurcating Navier-Stokes equations in a sudden-expansion channel, describing the qualitative and quantitative effect of the control over a pitchfork bifurcation, and commenting in detail the stability eigenvalue analysis of the controlled state. Finally, we propose reduced order modeling as a tool to efficiently and reliably solve parametric stability analysis of such optimal control systems, which can be challenging to perform with standard discretization techniques such as Finite Element Method.


Introduction
Nonlinear parametrized partial differential equations (PDE(µ)s) play a fundamental role in several fields, ranging from Continuum Mechanics to Quantum Mechanics, passing through Fluid Dynamics. When compared to linear PDE(µ)s with smooth parametric dependence, the most striking difference is that in the linear case a stable solution evolves in a continuous (and thus, unique) manner when the parameter changes slightly. In contrast, in the nonlinear case the solution for a given parameter µ of a given PDE(µ) may not be unique. Indeed, the model can suddenly change its behavior, together with the stability properties of its solutions. Models with such a feature are called bifurcation problems [3,16,52]. Examples of models which are characterized by the non-uniqueness of the solution are the Von Kármán plate model for buckling [59,8,12,44], the Gross-Pitaevskii equation for Bose-Einstein condensates [38,31,18,43] and Navier-Stokes equations in a channel [48,46,45].
The critical point at which the system loses its original features is called bifurcation point, and it is usually denoted by µ * . Many different bifurcating configurations can emerge from these points, and an intuitive way to visualize them is through a plot of a scalar output of the solution against the parameter value for which the solution has been computed, e.g. the so called bifurcation diagram. As an example, a branch of solutions can give rise at µ * to two further symmetric branches of solutions which coexist with the pre-existing one, and switch their stability properties with it. When such situation occurs, the model is said to undergo a pitchfork bifurcation phenomenon [52]. In particular, we are interested in an application in Fluid Dynamics, where we consider sudden-expansion channel flows, which are motivated by many practical scenarios. Flow profiles are described by Navier-Stokes equations, parametrized by the viscosity value µ. Here we consider a simplified version for a model of a cardiac disease, called mitral valve regurgitation, that may cause either a symmetric regurgitant jet or a wall-hugging, non-symmetric regurgitant jet. The latter phenomenon, which can be clinically detected through echocardiography, is called the Coanda effect [57], and expresses the tendency of a fluid jet to be attracted to a nearby surface. This represents an issue from the medical point of view, because the wall-hugging jet might lead to inaccurate echocardiography measurements. It is therefore of the utmost practical interest to try to drive the system towards the branch of symmetric solutions, which are more favorable for the measurement process.
Towards this goal, parametrized optimal control problems (OCP(µ)s) governed by PDE(µ)s might be employed. Indeed, OCP(µ)s can be interpreted as an input-output system which achieves an observable configuration [13,24,28,36,58]. They have been exploited in several applications in different scientific fields, see e.g. [35] for an overview. The main goal of this work is to use OCP(µ)s to drive bifurcating state profiles towards a different desired state, which might possibly belong to another state solution branch.
In such OCP(µ)s, a parametric study of the state solution is necessary in order to understand the behavior of the system. The complexity of this task can be tackled using a combination of existing methodologies, which allow the complete reconstruction of the aforementioned bifurcation diagram. However, the discretization through standard techniques can lead to a possibly huge system to be solved for many values of the parameters. In this work we exploit Galerkin Finite Element (FE) approach, which can be challenging when several instances of parametrized problem have to be studied, most of all because the optimal control setting requires additional equations to be solved on top of the original PDE(µ)s. For this reason we also propose Reduced Order Modeling (ROM) as a tool to overcome this issue, see [27] as an introductory reference. Namely, we build a lowdimensional space from FE high-fidelity optimal solutions and we perform a Galerkin projection in a reduced space, usually much smaller than the considered high-fidelity one. The main ingredients are the reduced space construction exploiting Proper Orthogonal Decomposition (POD) algorithm [6,15,17] and then, using a Galerkin projection, solving a reduced problem in the lower dimensional reduced space, for every parameter instance.
The main novelty of this work is the formulation, numerical simulation and subsequent model reduction of OCP(µ)s governed by bifurcating parametrized Navier-Stokes equations, together with the analysis of the properties of their eigenvalues. We develop several test cases corresponding to different control actions in order to understand and validate our findings. At the best of our knowledge, a thorough description of the mathematical analysis and widespread exploitation of OCP(µ)s for bifurcating nonlinear PDE(µ)s have not been extensively addressed in literature. This work could then pave the way towards application of OCP(µ)s as a tool to drive the behavior of complex bifurcation problems to a desired branch of solutions.
The work is outlined as follows: in Section 1 we describe the structure of OCP(µ)s for general nonlinear bifurcating systems both at the continuous and discretized level, exploiting a branch-wise approach through a FE approximation. Furthermore, we introduce the basic notions of stability analysis and its connection to eigenvalue problems for the uncontrolled and the controlled systems. In Section 2 we introduce the uncontrolled sudden-expansion channel problem, while in Section 3 we build several OCP(µ)s assigning different roles to the control variable. There we also perform a global eigenvalue analysis in order to understand the main features of the achieved optimal solutions. We present two different boundary control problems and two distributed ones, which give very different solution configurations. The ROM strategy is described in Section 4: the presented approach is tested for all the numerical simulations performed at the FE level. Conclusions follow in Section 5.

Nonlinear Parametrized Optimal Control Problems and Bifurcating systems
In this Section we introduce a generic nonlinear OCP(µ)s. We will focus on the minimization of quadratic cost functional under nonlinear PDE(µ) constraint for Hilbert spaces, following the Lagrangian approach [24,28]. In Section 1.1 and 1.2 we provide existence results and optimality conditions for nonlinear OCP(µ)s in their continuous and discretized version, respectively. Then, Section 1.3 will describe the spectral properties of the optimization system at hand.

Problem Formulation.
Optimal control is a mathematical tool which aims at modifying the natural behavior of a system. Let us suppose to have a state PDE(µ) (1) G(y; µ) = f, with state variable y := y(µ) ∈ Y, i.e. G : Y×P → Y * where Y is a Hilbert space, f ∈ Y * is a forcing term, P ⊂ R P is a parameter space of dimension P ≥ 1, while G(y; µ) = E n (y; µ) + E (y; µ) is the state operator, with E ∈ L(Y, Y * ) and E n representing the linear and nonlinear contributions, respectively. Here, we call L(·, ·) the space of linear continuous functions between two spaces. We now want y to be the most similar to a known solution profile y d := y d (µ) ∈ Y obs ⊇ Y. To this end, a new variable is introduced in the equation, the control variable u := u(µ) ∈ U, with U another possibly different Hilbert space. Let us define the controlled equation E : Y × U × P → Y * as where C ∈ L(U, Y * ) is the control operator describing the action of the variable u on the system 1 . In other words, we are trying to change the behavior of the state PDE(µ) through C(u). The OCP(µ) reads: given a µ ∈ P, find the pair (y, u) ∈ Y × U which solves (i) U is convex, bounded and closed; (ii) Y is convex and closed; (iii) for every µ ∈ P, the controlled system E(y, u; µ) = 0 has a bounded solution map u ∈ U → y(u) ∈ Y; (iv) for a given µ ∈ P, the map (y, u, µ) ∈ Y × U × P → E(y, u; µ) ∈ Y * is weakly continuous with respect to (w.r.t.) the first two arguments; (v) for a given y d ∈ Y obs , the objective functional J(y, u; y d ) is weakly lower semicontinuous w.r.t. y and u. We now discuss the Lagrangian structure and the necessary first order optimality conditions. First of all, let z := z(µ) ∈ Y * * = Y be an arbitrary variable called adjoint variable. Let us call X = (y, u, z) ∈ X := Y × U × Y and let us build the Lagrangian functional L : X × Y obs × P → R as where ·, · YY * is the duality pairing of Y and Y * . The introduction of the adjoint variable allows to treat problem (2) in an unconstrained fashion finding the stationary point of (4). We remark that we consider z in the same space of the state variable for a proper definition of the discretized problem: we will clarify the reason in Section 1.2. Moreover, the variable X will inherit the parameter dependence by definition, i.e X := X(µ). Furthermore, let us assume that the following hold: y and u; (viii) given µ ∈ P, the controlled system E(y, u; µ) = 0 has a unique solution y = y(u) ∈ Y for all u ∈ U; (ix) given µ ∈ P, D y E(y, u; µ) ∈ L(Y, Y * ) has a bounded inverse for all control variables u. 1 Parametrized control operators are also possible, with a straightforward extension of the methodology presented herein.
The Fréchet derivative w.r.t. a variable will be indicated as D , as already done in (ix). Assuming (y, u) ∈ Y × U to be a solution to (2) for a given µ ∈ P, thanks to hypotheses (vi) -(ix), there exists an adjoint variable z ∈ Y such that the following variational system is satisfied [28]: or, in strong form, where D y E(y, u; µ) * ∈ L(Y, Y * ) is the adjoint operator of the Fréchet linearization of E(y, u; µ) w.r.t. the state variable, while C * ∈ L(Y, U * ) is the adjoint of the control operator. We will refer to problem (5) as optimality system, in weak or strong form, respectively. Moreover, writing the latter in compact form, it reads: given µ ∈ P, find X ∈ X such that In the nonlinear case, even considering only the state equation, the solution for a given parameter µ may not be unique. Therefore, the local well-posedness of the problem (6), which strongly relies on its local invertibility assumptions (viii) and (ix), can fail due to the singularity of the state equation. Therefore, we can talk about solution branches, i.e. multiple solution behaviors for a given value of the parameter. We denote by k the number of branches, and by X i , i = 1, . . . , k the set of each solution on the i-th branch. We call solution ensemble the set of all the solution branches X i : In the next Section, we will discuss the FE approximation of a solution for a fixed value of the parameter to the nonlinear OCP(µ), restricting ourselves to the well-posed setting.
1.2. The FE Approximation. We are interested in the numerical approximation of the solution ensemble (7) of the nonlinear OCP(µ) in (6), defined over an open and bounded regular domain Ω ⊂ R d . Our aim is to discretize the problem at hand, in order to investigate its qualitative changes w.r.t. the values of the parameter. We remark that, even considering a single branch, the fulfillment of the well-posedness conditions can fail at some critical point µ * . Hence, in the following, we assume µ = µ * and X(µ) ∈ X i for some i ∈ {1, . . . , k}, thus we call X i a non-singular branch. Furthermore, we assume the nonlinearity to be at most quadratic in the state variable, guided by the numerical results we are going to present in the following Sections. However, the structure and the methodology do not change when dealing with nonlinearities of higher order. To approximate the system in (6), first of all we define the triangulation T N T of Ω, where K is an element of T N T and N T is the number of cells. Then, let us consider the discrete spaces Y Ny = Y ∩ K ry and U Nu = U ∩ K ru , where K r = {v ∈ C 0 (Ω) : v| K ∈ P r , ∀ K ∈ T N T }, and P r is the space of all the polynomials of degree at most equal to r. Let us consider the FE function The FE approximation of the parametrized problem (6) reads: given µ ∈ P, find X N := X N (µ) ∈ X N such that We now want to make the algebraic structure of the system (8) explicit. After the FE discretization, we can define y, u and z as the column vectors which entries are given by the FE coefficients of the state, control and adjoint variables in their approximated spaces, respectively. In the same fashion, we call y d the column vector of the FE coefficients representing the desired state profile. Let us focus on the structure of the optimality problem. At the FE level, applying our controlled state to the FE basis, we can derive the matrices E n + E − C and the forcing term vector f. Moreover, we define the mass matrices M y and M u for state/adjoint variables and control, respectively. We still need to understand the algebraic structure of D y E(y, u; µ). The Fréchet derivative of the controlled state equation w.r.t. the state y will be E n [y] + E . In other words, the linear state structure is preserved, the nonlinear operator is linearized in E n [y] and the control contribution disappears. Then, the global matrix formulation of the optimization system (8) is which in compact form reads: where R(X; µ) represents the global residual of the optimality system. To solve system (10), we rely on Newton's method and we solve until a residual based convergence criterion is satisfied. Since the matrix E n [y] T still depends on y, the Jacobian matrix will be of the following nature: does not depend anymore on the state, but only on the j−th value of the adjoint variable. Now, one can write We remark that, in the considered numerical settings, A is symmetric and, thus, we will always refer to (13) as a saddle point structure. To guarantee the solvability of the system we consider A an invertible matrix. Furthermore, we need the following Brezzi inf-sup condition to be verified: The assumption z ∈ Y, will guarantee the fulfillment of the inf-sup stability condition in the FE approximation [40,41]. For general nonlinear problems, A may possibly be different from A T , however the well-posedness results can be extended, and the interested reader may refer to [9,14].
In the next Section, we will describe the branch-wise procedure implemented to reconstruct the bifurcation diagram.

Bifurcation and stability analysis.
For nonlinear PDE(µ)s, the assumptions (vii)-(ix) ensure the applicability of the well-known Implicit Function Theorem [19,3]. Indeed, under those assumptions, one expects that when the parameter changes slightly, a stable solution evolves continuously in a unique manner. When such conditions fail to be fulfilled, the model may undergo bifurcation phenomena. In particular, we consider the case in which the state equation has a singularity at some parameter value µ * , and for the sake of clarity, we will assume f = 0 (or, equivalently, include the forcing term f in the expression of G). Indeed, from the mathematical perspective, we can give a precise definition of such points [3]. Definition 1.1. A parameter value µ * ∈ P is a bifurcation point for (1) from the solution y * := y(µ * ), if there exists a sequence (y n , µ n ) ∈ Y × P, with y n = y * , such that (i) G(y n ; µ n ) = 0, and (ii) (y n , µ n ) → (y * , µ * ).
Thus, the bifurcation phenomena is a paradigm for non-uniqueness in nonlinear analysis, and a necessary condition is the failure of the Implicit Function Theorem. Proposition 1.1. A necessary condition for µ * to be a bifurcation point for G is that the partial derivative D y G(y * ; µ * ) is not invertible.
Moreover, given the existence of many possible configurations for the system, a natural question is to understand which one inherits the stability of the unique solution, when it exists. To perform a stability analysis, one of the most common and widely-studied method is the spectral analysis of the problem, which consists in the investigation of the eigenvalues of the system. This technique, which has its roots in the ordinary differential equation (ODEs) theory [52,34,33], allows to understand the stability property of a solution to (1) by means of the sign of the spectrum of the model operator.
In particular, for a general nonlinear uncontrolled problem, one linearizes the equation around the solution under investigation,ŷ = y(μ), and then solves the eigenvalue problem given by (16) D y G(ŷ;μ)y e = ρμy e , where the pair (ρμ, y e ) represents respectively the eigenvalues and the eigenvector of D y G atŷ for eachμ fixed. This analysis provides us information about the physical stability of the problem. Indeed, the stability analysis is strongly bound to the investigation of the solution features after a small perturbation. If the perturbation is small enough and the dynamics of the system remains in a neighborhood of the solution, then it will be called a stable solution. Thus, it is fundamental to observe that, in connection with ODEs stability theory, a positive eigenvalue gives an exponentially divergent behavior, while a negative one produces only small oscillations around the solution. Therefore, it is clear that in order to have a stable solution, all eigenvalues must have negative real parts.
Dealing with the controlled problem (6) makes the analysis more involved. Indeed, the adjoint variable does not have physical meaning, making a straightforward application of the considerations above not applicable. Due to the high indefiniteness of the saddle point matrix (13), a standard sign-analysis is no longer possible, see [9,10,11] as references on the topic. Nevertheless, we can consider the eigenvalue problem for the system of the optimality conditions, in order to investigate a posteriori the spectral property of (6), as: (17) D X G(X;μ)X e = σμX e , whereX = X(μ) is the solution of which we are investigating the stability property and (σμ, X e ) is the eigenpair formed by theμ-dependent eigenvalues σμ and eigenvectors X e . We will refer to (16) as the state eigenvalue problem and to (17) as the global eigenvalue problem.
As we said, the main issue with bifurcating system is the lack of the invertibility for the Fréchet derivative of the operator G due to Proposition 1.1. The latter, being equivalent to the injectivity and surjectivity of G, can be rewritten in terms of the continuous Babuška inf-sup stability: there exists an inf-sup constantβ Ba > 0 such that It is clear that the inclusion property X N ⊂ X is only a necessary but not sufficient condition for (18) and (19) to hold at the discrete level. Hence, an additional assumption has to be required for the discrete Babuška inf-sup stability of G: there exists a constantβ N Ba > 0 such that We remark that the continuous condition on the surjectivity (19) is no longer needed for the discrete inf-sup stability (20). In fact, while the assumption (20) corresponds to the non singularity of the matrix Jac, a discrete counterpart of the assumption (19) would require its surjectivity, or equivalently the injectivity of the transpose matrix Jac T , which being square would be the same as requiring (20). We also highlight that both continuous and discrete inf-sup conditions are satisfied as long as µ = µ * .
We can finally present the branch-wise procedure we developed in order to deal with the numerical computation of multiple branches of solutions. Such an approach requires the combination of different methodologies for the approximation of the bifurcation diagram and the analysis of its stability properties. In order to keep the presentation simple, we consider that the first component µ of the parameter µ ∈ P ⊂ R P is the one responsible for the bifurcation behavior of the model, and in order to follow a single branch we consider that all the P − 1 remaining parameters (thus the global physical/geometrical configuration) are unchanged on the branch. Even though in this work we will only deal with bifurcation phenomena with co-dimension one, the following methodology can be adapted to the multi-parameter and/or co-dimension > 1 case, see e.g. [42], by carefully choosing the technique to follow the branching behavior.
Algorithm 1 summarizes how to reconstruct each branch X i of solutions. More precisely, we combine, respectively: • Newton's method, as the nonlinear solver, • Galerkin FE method, as the discretization phase, • simple continuation method, as the bifurcation path tracer, • generalized eigenvalue problems, as the stability detectors. At the very beginning, one has to chose the branch to approximate, and the most preferable way to "guide" the nonlinear solver to the desired configuration is through the initial guess. Thus, in order to reconstruct a branch we consider the discrete version of the parameter space P K = [µ 1 , . . . , µ K ] ⊂ P of cardinality K. We can take P K as an ordered set, with the natural ordering induced by the first parameter component. Such ordering serves to assign the solution obtained for a given parameter µ j−1 as the initial guess for the nonlinear solver at next iteration for µ j . This allows us to follow the bifurcation behavior of the model. We choose the simplest variant of the continuation methods [2], where the parametric set is fixed a priori, since it works well with pitchfork like bifurcation [44,43]. A more involved methodology have to be implemented when dealing with e.g. turning points or secondary bifurcations. The next step is the actual discretization of the problem by means of the Newton-Kantorovich method [19] combined with the Galerkin FE method. The initial guess for the former is set to the solution for the previous parameter value, while the latter projects the problem into a finite dimensional space, obtaining a linear system that we repeatedly solve until a convergence criterion is satisfied (here we chose a threshold tolerance for the norm of the global residual (10) at the i−th iteration of the Newton's method). In Algorithm 1, we call Jac y (ŷ,μ) the Jacobian matrix referred to the state equation (1) and with V and V y the scalar product matrices of the global optimization variable and of the state variable, respectively. Finally, having computed a solution X j of the problem (6) for the parameter µ j , we can investigate its stability properties solving the two generalized eigenproblems, to recover the physical stability and the spectral properties, respectively for the state and global eigenvalue problems. Remark 1.1. Moreover, we highlight that the choice of the initial guess is fundamental, but it is not always sufficient to recover the full bifurcation diagram. In such cases, different techniques were proposed, e.g. manipulating the set P K through predictor-corrector continuation methods, which involves pseudo-arclength strategy and homotopy [2,52]. Finally, when it is difficult to choose a proper guess, one can rely on: • the discretized version of analytic expressions, resembling the main properties of the sought solution [43]; • a deflation method, which requires only one initial guess, and discovers the full diagram preventing the convergence to already discovered solutions, helping the solver to find new branches [45,18]; • the eigenvectors of the global eigenvalue problem, that have been used in [44] to obtain the direction of the bifurcation branch in a neighborhood of the bifurcation points.
Algorithm 1 A pseudo-code for the reconstruction of a branch 1: X 0 = X guess Initial guess 2: for µ j ∈ P K do Continuation loop 3: end while 8: Jac y (y j ; µ j )y e = ρ µ j V y y e State eigenproblem 9: Jac(X j ; µ j )X e = σ µ j VX e Global eigenproblem 10: end for

Bifurcations for Navier-Stokes Equations: the Coanda Effect
In this Section we analyze a bifurcating phenomenon deriving from Navier-Stokes equations in a sudden-expansion channel flow problem. Indeed, consider the channel geometry depicted in Figure  1. A fluid characterized by a high viscosity presents a jet which is symmetric w.r.t. the horizontal axis. Furthermore, a pair of vortices, called Moffatt eddies [39], form downstream of the expansion. Lowering the viscosity, the inertial effects of the fluid become more important and the two symmetric recirculation regions break the symmetry. Indeed, as the length of the recirculation zones increases, one can observe a non-uniform decrease of the pressure along the vertical axis. Thus, when we reach the aforementioned critical value, one recirculation zone expands whereas the other shrinks, giving rise to an asymmetric jet. This phenomenon is called the Coanda effect and has been extensively studied in literature within different contexts [57,48,32,46,25,45,42].
From the mathematical point of view, this translates to a PDE(µ) which, decreasing the viscosity µ below a certain critical value µ * , admits the existence of more solutions for the same value of µ ∈ P. During the study of the solution for different viscosity values, we expect the system to show two qualitatively different configurations: • a physically unstable configuration with a symmetric jet flow, the symmetric solution, • a physically stable configuration with a wall-hugging jet, the asymmetric solution.
These solutions, depicted in Figure 3, coexist for parameter values below the critical one µ * and belong to different branches that intersect in the bifurcation point, forming the so-called pitchfork bifurcation.
In the next Sections we introduce the mathematical formulation of Navier-Stokes equations describing the flow in a channel. This will serve us to highlight the bifurcated behavior of the system and its stability properties and it will be fundamental to understand how different controls affect the original system.

Navier-Stokes problem as the state equation.
Here we consider a simplified setting with a two-dimensional planar straight channel with a narrow inlet and a sudden expansion, depicted in Figure 1, which represents a simplification of the left atrium and the mitral valve, respectively. We define Γ in = {0} × [2.5, 5] and Γ out = {50} × [0, 7.5], where inflow and outflow boundary conditions are imposed, respectively. We indicate with Γ wall , the boundaries representing the walls, in this case The steady and incompressible Navier-Stokes equations for a viscous flow in Ω read as: is the velocity of the fluid, p is its pressure normalized over a constant density and µ represents the kinematic viscosity. We supplement the system (21) with proper boundary conditions: a stress free boundary condition on the velocity at the outlet, Γ out with outer normal n, a no-slip (homogeneous) Dirichlet boundary condition on Γ wall and a non-homogeneous Dirichlet boundary conditions v in at the inlet Γ in given by v in ( For later convenience, we introduce the dimensionless Reynolds number, which represents the ratio between inertial and viscous forces, and is given by Re = U h/µ, where U and h are characteristic velocity (i.e., maximum inlet velocity, U = 31.25) and characteristic length of the domain (i.e., length of the inlet section, h = 2.5), respectively. In the following we will consider µ as the parameter. Fixed the domain Ω, the flow regime varies as we consider different values for the viscosity µ in P ⊂ R. As we said in the introduction to this Section, this model exhibits a bifurcating behavior. Indeed we have the existence and uniqueness of the solution only above a certain critical value for the viscosity, that for this test case corresponds to µ * ≈ 0.96. Such value has been found in different works and numerical contexts for this benchmark [45,32,42]. It can be obtained either "a posteriori" by looking at the behaviour of the flow while varying the viscosity, or "a priori" by investigating the change of sign of the leading eigenvalue w.r.t. the parameter. To investigate the loss of uniqueness in a neighborhood of this pitchfork bifurcation, we set the parameter space as P = [0.5, 2.0], such that the first critical point µ * is included. These values for the viscosity correspond to Re in the interval [39.0, 156.0].
be the function spaces for velocity. Furthermore, let Q = L 2 (Ω) be the function space for pressure. The weak formulation of (21) reads as: given µ ∈ P, find v ∈ V in and p ∈ Q such that Figure 2. Uncontrolled Navier-Stokes system: Bifurcation diagram.
We can rewrite the formulation of (22) in an equivalent way as: given µ ∈ P, find v ∈ V in and p ∈ Q such that having introduced the following bilinear and trilinear forms for all v,v, ψ ∈ V and p ∈ Q,

Numerical approximation of the problem.
We can now discuss the numerical approximation of the Navier-Stokes equation, that will be the state equation of the control problem in the next Sections. We consider a mesh on the domain Ω with N T = 2785 cells and N y = 24301 degrees of freedom associated to a Taylor-Hood P 2 -P 1 discretization of V × Q. This choice is motivated by the well-known stability results of the Taylor-Hood Finite Element pair [50]. In order to plot the bifurcation diagram, we choose an output value that results in a symmetry indicator of the approximated solution. This function is given by the value of the vertical component of the velocity in a point of the channel, i.e. v x2 evaluated for (x 1 , x 2 ) = (14, 4), or the nearest node to that value since the mesh is unstructured. We remark that this output provides a merely graphical intuition about the symmetry breaking, and can be chosen to be efficiently computed by means of the analysis in Section 4. In Figure 2 we plot the bifurcation diagram with all the solution branches found for the system (23) in the viscosity range chosen. The numerical approximation clearly shows that a supercritical pitchfork bifurcation occurs around the critical viscosity value µ * ≈ 0.96. It is evident that we have a unique solution for all µ > µ * , thus when the fluid behaves like a Stokes one, while we find three qualitatively different solutions increasing the Reynolds number. The bifurcation point µ * is also the one responsible for the change in stability properties of the model. Indeed, the unique symmetric solution remains stable until it encounters the critical value µ * , where it becomes unstable. Moreover, this feature is inherited by the bifurcating solutions, which evolve as a physically stable branch.
Some representative solutions for the lower and middle branch are presented in Figure 3. Velocity and pressure fields belonging to different branches, for the same viscosity value µ = 0.5, present qualitatively dissimilar behavior. Indeed, the pressure for the lower branch decreases near the bottom-left corner of the expansion, causing the velocity to deflect, hugging the lower wall. Finally,  thanks to the no-slip boundary condition the flux goes back to mid line, ending with a non-axissymmetric outflow.
The stability analysis is performed through the eigenvalue analysis depicted in Section 1.3, where Algorithm 1 has been applied exclusively to the state equation (1). In particular, we analyzed the behavior of the first N eig = 100 eigenvalues of (16), by means of the Krylov-Schur algorithm, varying the viscosity of the fluid. Such eigenvalues are plotted in Figure 4 for the stable lower branch (left panel) and unstable middle branch (right panel). Note that since the Navier-Stokes operator is not symmetric, we have the presence of both real and complex eigenvalues. Here we are just interested in pitchfork bifurcation, thus in the zoom we follow the behavior of the biggest real eigenvalue and its sign. When investigating the stability of the lower branch, all the eigenvalues of the Navier-Stokes system, linearized around this stable solution, have negative real part. From the consideration of Section 1.3, we can assert the stability of the wall-hugging branch. Indeed, the zoom in the left plot of Figure 4 shows no crossing of negative-real part eigenvalues. On the contrary, the close up in the right plot, which corresponds to the symmetric flow, shows the sign change of the biggest eigenvalue, thus characterizing a physically unstable solution.

Steering bifurcating governing equations towards desired branches by means of Optimal Control
In this Section we focus on several OCP(µ)s governed by Navier-Stokes equations (21) in the geometrical configuration of a contraction-expansion channel. We aim at understanding how different control problems can affect the solution behavior discussed in Section 2 for the uncontrolled case, especially when bifurcation phenomena are taken into account. This leads us to analyze the controlled systems, trying to reach state profiles which are different from the expected uncontrolled solution. Our goal is to investigate and better understand the role that optimal control plays as an attractor towards a desired configuration. We thus follow the general procedure described in Section 1.1, and discuss the specific case of optimal control for the Coanda effect. Nonetheless, the procedure adopted here is general and can be used in wide variety of applications. For all the applications, we will simulate the physical phenomenon over the domain Ω shown in Figure 1. Moreover, for the OCPs structure, we will require the velocity solution v ∈ V to be the most is a line near the end of the channel. This structure allows the control to change the solution at the outflow following a prescribed convenient configuration. During the rest of this work, we will employ two velocity solution profiles, which are showed in Figure 5: we will denote them as the symmetric desired profile (or target) for Figure 5a and the asymmetric desired profile (or target) for Figure 5b.
The first is the result of a Stokes system over Ω for µ = 1 with the same boundary conditions of the Navier-Stokes uncontrolled equations (21). On the contrary, the latter is the physically stable solution of (21) for µ = 0.49. While the former choice aims at controlling the system towards a globally symmetric configuration with a weaker outgoing flux, the latter is set to achieve the opposite goal. The purpose of steering the bifurcating behavior is summarized in the minimization of the functional where U := (L 2 (Ω u )) 2 with Ω u ⊂ Ω: indeed, the control action can be performed even over a portion of the boundary ∂Ω. We will refer to Ω u as the control domain. Within this study we will analyze how the choice of Ω u , combined with different values of the penalization parameter α, affects the solution behavior of the system, compared to the uncontrolled Navier-Stokes state equation. We remark that hypothesis (i)-(ix) are verified in this specific OCP(µ) setting, see e.g. [28]. Then, after a general introduction in Section 3.1, we will provide the analysis of different optimal control systems, as follows. Section 3.2. A weak control is built by controlling a Neumann boundary and the optimality system slightly affects the usual bifurcating nature of the uncontrolled Navier-Stokes equations. Section 3.3. A strong control effect can be observed over the classical bifurcating behavior of the uncontrolled solution by acting on the forcing term. Section 3.4. The penalization parameter α is analyzed while acting at the end of the inlet channel, and we discuss how changing α results in different orders of magnitude for the optimal control. Section 3.5. We show how imposing different boundary flux conditions completely changes the known behavior of the starting system. Finally, in Section 3.6, remarks and comparisons on the spectral analysis of the four test cases are presented.

OCP(µ)s governed by Navier-Stokes equations.
We recast OCP(µ)s constrained to Navier-Stokes equations in the algebraic formulation presented in Section 1.1.
The steady and incompressible controlled Navier-Stokes equations in a given domain Ω are: accompanied by some boundary conditions. The control operator C : U → V * can represent an external forcing term or a boundary term. If C is defined in the whole domain we will say that the control is distributed, while if it is defined in a portion of the internal domain, we will deal with localized control. Furthermore, we will refer to Neumann control and Dirichlet control, if the control acts as Neumann or Dirichlet boundary conditions, respectively. The weak formulation of (26) reads: given µ ∈ P, find v ∈ V, p ∈ Q and u ∈ U such that where a(·, ·; µ), b(·, ·) and s(·, ·, ·) have been already defined in (24) while c : U × Y → R is a bilinear form associated to the operator C. First of all, to derive the optimality conditions, we need the adjoint variables w ∈ V and q ∈ Q for velocity and pressure, respectively. Let X = ((v, p), u, (w, q)) ∈ X := Y×U×Y be an optimal solution, where Y := V×Q. The Lagrangian functional for this specific problem is The optimality system built through Fréchet differentiation is given by: where the first two equations form the adjoint system, the last two form the state system (27), while differentiating w.r.t. the variable u leads to the optimality equation. In particular, the adjoint system combined with the optimality equation has the following form: where m : V × V → R and r : U × U → R terms come from the Fréchet derivative of (25) w.r.t. the velocity and control, respectively. They represent the L 2 scalar product in Γ obs and Ω u , respectively. Furthermore, we remark that s(ϕ, v, w) + s(v, ϕ, w) is the linearization around v of the trilinear form s(v, v, ϕ), by definition. Therefore, the strong formulation for (30) reads: where I Ωu and I Ω obs are the indicator functions of the control and observation domains, respectively. The detailed derivation of the optimality system is addressed in several works, see e.g. [28,21,23]. The global optimization problem reads: given µ ∈ P, find X = ((v, p), u, (w, q)) ∈ X such that (26) and (31) are verified. We remark that, if we call y := (v, p) and z := (w, q), we recover the global algebraic formulation presented in Section 1.1 and the saddle point structure is preserved. Indeed, let us suppose that we apply the Taylor-Hood approximation P 2 -P 1 for state y and adjoint variable z. Furthermore, we discretize the space U with FE using P 2 polynomials. Recalling the notation of Section 1.2, we define the quantities where v, p, w, q are the column vectors of FE coefficients for state and adjoint, velocities and pressures respectively, while M v is the mass velocity matrix and C v derives by the bilinear form c(·, ·). Furthermore, the linearized state equation structure can be now expressed as where K is the stiffness matrix associated to the bilinear form a(·, ·; µ), D is the continuity equation matrix coming from b(·, ·) and S[v j ] is the algebraic formulation of s(v, ·, ·) + s(·, v, ·) evaluated at the FE velocity basis functions. It remains to understand the specific structure of D y (E n [y] T )[z j ] defined in (12). To this end, we define s ad (v, w, ϕ) as the adjoint operator of the linearized trilinear form s(v, v, ·) around the state velocity v. Applying s ad (v, ·, ·) to the basis functions of V Nv will result in S[v j ] T . In the Jacobian matrix evaluation, a linearization of s ad (w, v, ϕ) is performed not only in w, but also w.r.t. the variable v. This process will lead to is given by the form s ad (w, ·, ·) applied to the FE velocity basis. Then, the Jacobian reads where X is the FE coefficient vector of the optimal solution and As already specified in Section 1.2, we assume that for µ = µ * the saddle point (35) is well-posed. Moreover, we highlight that we are dealing with a nested saddle point structure: indeed, for the state equation (33) we require that, for a given µ = µ * and fixed v j , the matrix K + S[v j ] is invertible and that Brezzi inf-sup condition holds, i.e.
This is indeed the case for the Taylor-Hood discretization introduced in Section 2.
In the next subsections we analyze how the controlled problem behaves, comparing its properties with the ones of the uncontrolled system presented in Section 2. We will use the word natural optimal branch to describe the branch that is obtained by running Algorithm 1 with a trivial initial guess. This branch may consist of either symmetric or asymmetric configurations, depending on the test case. Further branches may exist, but are much harder to compute in practice and require very tailored initial guesses that can be provided by running Algorithm 1 in a neighborhood of µ * , and will be named non-natural optimal branches 2 . We interpret the concept of natural optimality as a numerical stability property of the optimal control system. Indeed, as already specified in Section 1.1, for OCP(µ)s it makes no sense to talk about the physical stability of the global optimal solution.
In fact, the system is "artificially" built by adding non-physical adjoint variables, with the aim of changing the system behavior.
3.2. Neumann Control: weak steering. The first test case we present is a Neumann control over the boundary Γ out , where homogeneous Dirichlet conditions are applied to Γ wall := Γ 0 ∪ Γ D . More specifically, in this case, the optimality conditions read: given µ ∈ P find X ∈ X such that The desired velocity v d will always be of the symmetric type for this specific example. In other words, we are studying which is the best choice for Neumann boundary condition, to reach the exiting symmetric profile shown in Figure 5a. We study the behavior of the controlled solution varying α = 1, 0.1, 0.001, 0.0001, where the greater is the value of α the lower is the strength of the control. In Figure 6 we show some representative solutions for α = 0.01 and µ = 0.5, for state velocity and pressure variables. In this case, the natural optimal branch is composed by asymmetric solutions (Figure 6, top), while there is a further non-natural optimal branch made up by symmetric solutions (Figure 6, bottom). Results obtained following the natural optimal and non-natural optimal branches are shown in Figures 7 and 8, respectively. Therefore, we conclude that the Neumann control affects weakly the system, as it is not able to steer the system towards the desired symmetric configuration after bifurcation has occurred, thus not changing drastically the features already observed for the uncontrolled state equations (see Figure 2). The left plot of Figure 7 depicts the velocity profile magnitude over Γ obs for the highest value of the Reynolds number when following the natural optimal branch. Even though the obtained velocity (marked by an orange line) is indeed different from the desired profile (denoted by a blue line), especially for what concerns peak values, we observe that the Neumann control straightens the flux near the end of the channel (compare the orange line to the green line, which represents the uncontrolled asymmetric profile), even when high Reynolds numbers are considered. The resulting profile is similar to the uncontrolled symmetric velocity (red line), even though full symmetry is not achieved. The action of the control variable is shown in the right plot of Figure 7 when changing the parameter µ following the natural optimal branch: the control is stronger for µ < µ * (i.e., when the wall-hugging phenomenon occurs and straightening in necessary), while it remains low in magnitude for µ > µ * . Similarly, the left plot of Figure 8 shows the velocity profile magnitude over Γ obs for the highest value of the Reynolds number when following the non-natural optimal branch. In this case, the controlled symmetric profile (orange line) coincides with the uncontrolled symmetric profile (red line). Furthermore, the right plot of Figure 8 shows that the control variable around the critical µ * (e.g., µ = 1 and µ = 0.95) is asymmetric to counteract the stable wall-hugging physically driven behavior of the uncontrolled system. We further remark that, compared to the natural optimal branch, the control variable of the non-natural optimal branch is much lower in magnitude. Table 1 shows the value of the cost functional (25) for several values of µ (rows) and α (columns), following either the natural optimal or non-natural optimal branches. The first column also shows the value of the uncontrolled functional, i.e. (25) evaluated for the uncontrolled velocity v of the equation (21) and zero control. The main observation is that decreasing the value of α results in (c) (d) Figure 6. Neumann Control: optimal solutions with α = 0.01 and µ = 0.5, belonging to the natural optimal (panels (a) and (b) for state velocity and pressure, respectively) and the non-natural optimal (panels (c) and (d)) branches. lower cost functional values, since a lower value of α allows stronger control to take place and drive the velocity to the desired configuration. In all cases, the non-natural branch presents lower values of the functional compared to the natural branch; this has to be expected, as the cost functional measures deviation from a symmetric target, and the non-natural branch is clearly closer to the target being made of symmetric solution (compare e.g. for µ = 0.05 and α = 0.001 the left panels of Figures 7-8). However, the natural branch is the one for which the control procedure is influencing the most the cost functional values: for instance, for µ = 0.05 and α = 0.001, the cost functional is decreased by 6% on the non-natural branch and by 55% on the natural one w.r.t. the corresponding uncontrolled configuration. Again, this has to be expected from the previous discussion of Figures  7-8, which shows larger impact of the control procedure to straighten the solution on the natural branch. Finally, large values of µ have negligible cost functionals, as the target velocity almost coincides with the uncontrolled velocity. From such an analysis we deduce that, when bifurcating phenomena occur, a configuration can perform better than another one, and finding all the solution branches can be of great importance to understand the solution that best recover the desired profile.
Concerning the stability of the solution, we performed the eigenvalues analysis described in Algorithm 1. We can derive several insights from the Figure 9, which represents the global eigenvalue

Distributed Control: strong steering.
This Section deals with a distributed control in Ω u ≡ Ω, thus the control variable u acts as an external forcing term on the whole domain. Here we consider again Γ wall = Γ 0 ∪ Γ D . Given µ ∈ P, the optimal solution X ∈ X satisfies the following system: First of all, we underline that in distributed OCP(µ)s the action of the control is usually stronger and it affects deeply the original system. To show this, we will steer the system towards either symmetric or asymmetric desired profiles v d : • Symmetric target: the aim of this setting is to steer the solution of (39) to a symmetric profile. We plot two representative control solutions in Figures 10a and 10b, obtained for µ = 2 and µ = 0.5 when following the natural optimal branch, which is composed of symmetric solutions. The stronger action of the control allows the controlled velocity profile to be more diffusive compared to the uncontrolled symmetric profile, as represented in the left plot of Figure 11, corresponding to the observed slice of the velocity solution for µ = 0.5: in this case the controlled velocity (orange line) and the symmetric target (blue line) almost coincide. The right plot of Figure 11 shows that a slightly asymmetric control is required only near the critical value µ * (also compare to Figures 10a and 10b for the cases µ = 2 and µ = 0.5). Furthermore, the control action is clearly higher when the Re value increases. Indeed, for µ = 2 the control exclusively acts in the proximity of Γ obs with a maximum magnitude of 1.8·10 −4 , while for µ = 0.5 it reaches a value of 1.6 of magnitude. A further non-natural optimal branch exists, and is made of symmetric solutions, but is hardly reachable by numerical continuation methods unless tailored guesses are provided restarting Algorithm 1 in a small neighborhood of µ * . • Asymmetric target: in this case, we desire to recover the asymmetric target for all µ ∈ P.
We plot two representative control solutions in Figures 10c and 10d, obtained for µ = 2 and µ = 0.5 when following the natural optimal branch, which is made of asymmetric solutions. The action of the control is also visible in the left plot of Figure 12, obtained for µ = 2: indeed, we see how the flux over Γ obs is pushed towards the domain wall (orange line), in contrast to the symmetric profile of the uncontrolled velocity (green line). Namely, also in this case, the distributed control is able to drive the solution towards the desired state. In order to do so, the control variable has to be large when µ > µ * , i.e. when the uncontrolled configuration on Γ obs would lead to a symmetric profile. Indeed, in Figure 10c the maximum control value reaches 7 for µ = 2 in the upper part of the domain. In contrast Figure 10d it lowers to 10 −11 for µ = 0.5 when the stable asymmetric velocity solution does not need to be controlled by an external forcing term. This is confirmed in the right plot of Figure 12 for several values of µ. Also in this case a non-natural optimal branch (featuring symmetric solutions) continues to exist, but it is numerically difficult to reach.
We plot the eigenvalues for α = 1 in (σ µ ) = [−0.01, 0.01] and for α = 0.1 in (σ µ ) = [−0.005, 0.005]. For this test case, the predominance of positive eigenvalues is visible also for large values of the penalization parameter. The smaller is α, the more negative eigenvalues are lowered. For all the α taken into account, the shears phenomenon does not appear: for α = 1 a small trace of the shears structure is still visible (highlighted in blue) in Figures 14a and 14c where the bottom part of the shear is pushed away from (σ µ ) = 0. Instead, for α = 0.1 the shears structure is completely lost: Figures 14c and 14d show that only one eigenvalue (representing the top of the shears, and marked in blue) approaches (σ µ ) = 0 without crossing it.
We finally notice that the point µ * * , where the upper shears curve is the closest to the axis (σ µ ) = 0, allows to obtain further information on the bifurcating phenomenon. From Figure 14, µ * * ≈ 0.96 for the symmetric target, regardless of α, while requiring an asymmetric target leads to µ * * ∈ [1.0, 1.2] with a mild dependence on α.
Thanks to the aid of Figure 13, which shows the bifurcation diagram for the controlled solution with asymmetric target (a similar plot can be obtained for the symmetric target as well, but is here omitted because the lines almost overlap), we can state that optimal control is not only able to steer the state solution towards a desired branch, but may also affect the location of the bifurcation point. The role of the penalization parameter α will be clarified in the next Section and it will result into a completely new optimal solution behavior in Section 3.5.

Channel Control: the α effect.
This Section aims at describing how the value of the penalization parameter α can affect the natural convergence towards a symmetric target over Γ obs . Towards this goal, we analyzed the action of a control variable defined at the end of inlet channel, i.e. Ω u = Γ ch , as depicted in Figure 1. The boundary Γ wall is, once again, Γ 0 ∪ Γ D . Within this setting, the problem reads: given µ ∈ P find the optimal solution X ∈ X such that the following holds The optimal control acts as a forcing term capable to change the way the flow enters in the expansion channel. In Figure 15 we show the adjoint velocity and pressure profiles obtained for µ = 0.5 and two different penalization values, namely α = 1 and α = 0.01. In the first case, following Algorithm 1, the natural optimal branch presents a wall-hugging behavior, while for smaller values of α the  control variable is able to drive the velocity towards a straight flux (see the left panels of Figures  16 and 17). Therefore, for large values of α the natural optimal branch is composed by asymmetric solutions (i.e., far away from the target), while for smaller values of α the natural optimal branch is made of symmetric solutions. From the right plots of Figures 16 and 17, the control is very sensitive close to µ * and this is shown by its asymmetric configuration both for the wall-hugging solution and the straight one. For α = 1, 0.1, 0.01, we were able to detect two solutions using different initial guesses in the continuation method, showing symmetric and asymmetric features coexisting for some values of µ < µ * . The smaller was α, the more difficult was to recover the non-natural branch. For example, when α = 0.001, the action of the control variable drives the wall-hugging phenomenon towards a straight flux so strongly that we were not able to really reconstruct the whole optimal non-natural branch. Indeed, either the Newton's solver did not converge (this happens also for α = 0.1 and µ = 0.5, compare Table 3) or converged to the natural branch consisting in symmetric features.
As usual, the role of α is highlighted in reducing objective functional, as Table 3 shows. As already specified in Sections 3.2 and 3.3, the straight configuration is lowering much more the functional than the asymmetric solution, due to its similarity with the symmetric v d , which is the fixed target for this test case. In this case, the role of α is crucial in order to reach a solution which represents better the desired state. Indeed, the control was able to steer the solution to the symmetric profile
for α = 0.1, 0.01, 0.001. From the functional point of view, we do not have a notable decrease, as can be observed in Table 3, where the value of (25) is presented for different values of µ and α w.r.t. the uncontrolled problem solution. Yet, acting at the end of the inlet channel still allows to drive the optimal solution towards a natural convergence to the symmetric v d , but the parabolic profile on Γ obs is not reached (the functional decreases only of a 10% for µ = 0.5 and α = 0.001 w.r.t. the uncontrolled symmetric solution). Figure 18 shows the eigenvalues of the global eigenproblem in the range (σ µ ) = [−0.01, 0.01] when following the natural optimal branch. For α = 1, we can observe the shears phenomenon, which disappears for other values of the penalization parameter. Lowering the value of α, leads to a positive-dominated eigenvalues ensemble. Furthermore, a clustering around the value of α can be observed in plots 18c and 18d. In the next Section, very peculiar features have been observed as well, while changing the value of the penalization parameter α in a Dirichlet control.
3.5. Dirichlet Control: flux action. In this final example, we propose a Dirichlet control over the boundary Ω u ≡ Γ D . We fix the symmetric configuration as desired state v d over the line Γ obs , while we set Γ wall = Γ 0 . In other words, we are trying to control a Dirichlet boundary condition in order to lead the controlled solution towards the symmetric profile. The problem to be solved reads: Figure 18. Channel Control: spectral analysis with α = 1, 0.1, 0.01, 0.001 following the natural branch.
given µ ∈ P find X ∈ X such that Allowing the flux to freely enter or exit from the boundary Γ D drastically changes the optimal solution behavior. Since we are asking for a symmetric desired profile, the main action of the control is to straighten the flow: this behavior can be observed from Figure 20a and the left plot of Figure  19. Indeed, even for large values of α, the velocity profile reaches the symmetric configuration, while for lower values of the penalization parameter, the velocity on Γ obs is parabolic. This feature is highlighted also from the functional values in Table 4, where the functional (25) is shown for several µ (rows) and α (columns) w.r.t. the uncontrolled stable and unstable solutions. The cost functional largely decreases for smaller values of α, e.g. α = 0.001: for example, focusing on µ = 0.5, Figure 19. Dirichlet Control. Left: comparison of velocity profiles in the controlled and uncontrolled cases for α = 1, 0.01, µ = 0.5 on Γ obs w.r.t. the symmetric desired profile when following the natural optimal branch. Right: representation of control variable evolution for α = 0.1, 0.01, 0.001, 0.001 and µ = 0.5 at x1 = 10 when following the natural optimal branch. the functional only lowers of 18% for α = 0.01, in contrast to almost 82% for α = 0.001. Within the setting of α = 0.001, the system manifests an interesting and unexpected profile, shown in Figure 20b. The flux presents an asymmetric configuration for low value of µ. Namely, for low α a bifurcating solution appears as depicted in Figure 20b. The asymmetric behavior is due to the control variable which not only allows the flow to exit from Γ D (in order to avoid the asymmetric recirculation of the wall-hugging solution), but it is adding flux near the channel, in order to achieve the straight configuration and the parabolic velocity profile given by the symmetric target velocity over the observation domain, as it is represented in the right plot of Figure 19.    Figure 21 provides the eigenvalue analysis, where we show some close-ups starting with (σ µ ) = [−0.001, 0.001] for α = 1 in the top-right image, and restricting the vertical interval due to the order of the lowered α in the remaining images. As already noticed in Section 3.3, the strongest is the control, the more the eigenvalues are positive. Furthermore, as it already happened in the distributed control case with asymmetric target in Section 3.3, the value µ * seems to be not relevant anymore as an indication of the bifurcation point. Recalling the definition of µ * * in Section 3.3 as the value of the parameter for which the top curve of the shears (marked in blue) is approaching (σ µ ) = 0 from above, Figure 21 shows that such a curve is moving away from (σ µ ) = 0 as α decreases, and thus no such point µ * * exists. The results of the previous Sections have shown that the top shear structure is typically associated to the wall-hugging bifurcation, and that µ * * provides an indication of the bifurcation: we are thus lead to believe that the standard bifurcating configuration observed in the uncontrolled case, consisting of a branch of symmetric solutions and a branch of wall-hugging ones, is not present here, with the latter branch disappearing. However, the system seems to be featuring a different bifurcation, presented in Figure 20b. Indeed, we can see an eigenvalue crossing the line (σ µ ) = 0 for the global eigenproblem in Figure 21d for α = 0.001. Therefore, if we plot the eigenvalues of the state eigenproblem of Algorithm 1 (Figure 22), we see how the symmetric profile does never cross the origin, while the asymmetric solution in Figure 20b for α = 0.001 does. In the setting with the modified BCs, the physical stable solution behavior is a feature of straight profile. Moreover, from Figure 22, we can clearly observe a couple of complex and conjugate eigenvalue crossing the imaginary axis. This is, in fact, a paradigm for Hopf bifurcation [48,52] and represents another evidence of how deeply the system changed its inner features.
Remark 3.1 (Lagrange multipliers). From a numerical point of view, we employed Lagrange multipliers to solve the optimality system (41). The condition v = u on Γ D has been weakly imposed in integral form. This will result in extra terms in the adjoint equations. Furthermore, we weakly impose also the boundary condition w = 0 on Γ D with another multiplier. The reason of this latter decision will be explained in Section 4.
3.6. Comparative Eigenvalue Analysis. In this Section, we sum up all the observations and results derived from the global eigenvalue analysis over the four test cases. Therefore, we now list the similarities between these, especially for what concerns the variation against different values of the penalization parameter α: • the eigenvalues cluster around the value of α. This behavior is well represented in Figures  9c, 9d, 18c and 18d. These eigenvalues come from the optimality equation; • the predominance of positive eigenvalues over the negative ones. In all the performed spectral analysis we have observed that the control action lowers the negative eigenvalues. The stronger is the control, the greater is the number of positive eigenvalues, as it is represented in Figures 14b and 21; • the shears effect for low controlled system. The shears configuration is characteristic of control problems which do not highly change the uncontrolled system solution. It is the case of Neumann control in Figure 9a and of the channel control for α = 1 as shown in Figure  18a. For the other cases, the smaller is α the least visible is this eigenvalue configuration: in some cases, the structure is completely broken; • the µ * * identification. It is clear that the shears (or their top curve, if shears are broken) approach to the real line in the point µ * * for which the bifurcating phenomenon of the controlled system occurs. This is the same situation we found for the uncontrolled problem, in which the path of the eigenvalue identifies the value of bifurcation parameter µ * . Moreover, this situation is often preserved regardless of α. Indeed, the positive shears eigenvalue is still present in Figures 14b, 18b, 18c and even in the Dirichlet optimal control, as shown in Figure 21. In some cases, a shift of the µ * * compared to µ * has been observed.
Since the structure of the spectral analysis is highly influenced by the control strength, we tried to perform an eigenvalue analysis dealing with only state and adjoint equations. For all the test cases, shears occurs. The shears structure is symmetric when the solution shows the wall-hugging property, while it is slightly asymmetric when the state flow is straight. We guess that this behavior is due to the different reaction of state and adjoint blocks to the bifurcating phenomena. Indeed, for the symmetric flux, the behavior of the state equation has to be preserved for all µ, while the adjoint problem, which is strictly linked to the control variable, puts more effort in rebalancing the flux, resulting in an asymmetric contribution that causes the shears to be slightly asymmetric.
The spectral analysis of a nonlinear system is a indispensable tool to understand bifurcation phenomena, eventually. Under the point of view of computational costs, it is a very tough task, most of all for nonlinear OCP(µ)s. Indeed, FE discretization leads to huge systems to be solved for a wide sample of parameter µ ∈ P. In the next Section we propose ROMs as a suitable approach to overcome this issue.

ROMs for Nonlinear OCP(µ)s
This Section introduces ROM approximation techniques for nonlinear OCP(µ)s. The proposed reduced strategy is independent from the governing state equation. Indeed, we refer to [54,55,56,60] for previous contributions to ROM for nonlinear OCP(µ)s, extending them to standard techniques proposed in [25,26,43,44,45,47], for bifurcating systems. Section 4.1 introduces the basics ideas of ROM approach and the standard assumptions to guarantee its efficiency and applicability. Moreover, in Section 4.2, we will describe the used reduction strategy, relying on POD-Galerkin basis construction, see [6,15,17,27] as introductory references, combined with aggregated spaces techniques, following the linear OCP(µ)s fashion, as already presented in [4,5,20,22,29,30,40,41,49]. Finally, numerical results are shown in Section 4.3.

General Reduction Strategy.
In Section 3, we analyzed how optimal flow control can modify an expected behavior in bifurcating systems. In such problems, a parametric study of the state solution is necessary in order to understand the solution properties. The FE method cannot be exploitable when several instances of parametrized problem have to be studied. Indeed, the computational time required can be too costly, especially in an optimal control setting, where the state equation is associated to an optimization system. ROM techniques aim at building a reduced surrogate of the FE approximation. The main idea of ROM is to spend computational resources to build a new (smaller) model starting from FE solutions. Such a reduced system, although having a considerably lower dimension, is built in a way that it does not lose accuracy and can be used to analyze several parametric configurations in a versatile low-dimensional framework. We now briefly introduce ROM ideas in the OCP(µ)s setting. At the continuous level, we have already defined the solution branch X i , which represents how an optimal solution X(µ) = (y(µ), u(µ), z(µ)) of (6) changes w.r.t. the parameter µ ∈ P. For the sake of clarity, in this Section we will underline the parameter dependence of our solution branch, since it is of primary importance to understand the ROM basics and notions. From now on we will refer to FE as high-fidelity approximation. Indeed, we suppose that the FE discretization reliably represents the solution branch X i as its discretized counterpart X N i . The ROM aims at representing the high-fidelity X N i through the construction of basis derived from snapshots, i.e. properly chosen solutions of (8). The reduced function spaces are contained in the FE spaces: a standard Galerkin projection is performed, resulting in an efficient low-dimensional solution which is still accurate w.r.t. the high-fidelity model. We exploited a branch-wise reduction, namely, for every bifurcating solution branch X i , we build a different ROM. Of course this is the best approach form the accuracy standpoint, but in [44,43] a different context a global approach was pursued, where the loss in accuracy is balanced with the construction of a single ROM. We suppose to have already constructed the reduced spaces Y N ⊂ Y Ny ⊂ Y and U N ⊂ U Nu ⊂ U, the former for state and adjoint variables, the latter for control, respectively (description of the construction of reduced spaces is postponed to Section 4.2). The function space dimension N is usually much lower than N . We will refer to this stage as offline phase: here, apart from the basis construction, all the parameter independent quantities are assembled and stored. After this reduced spaces building process, one can solve the following low-dimensional problem in an online phase given µ ∈ P, find X N (µ) = (y N (µ), u N (µ), z N (µ)) ∈ X N := Y N × U N × Y N such that it holds: Namely, at each new parametric instance µ ∈ P, the system (42), which inherits the features from the high-fidelity dynamics (8), is assembled and solved. Also in this case, we can deal with the nonlinearity applying a Newton's method, as in the FE setting. Moreover, the stability analysis can be performed as described in Algorithm 1, employing the Galerkin projection in the reduced space. It is clear that, it is convenient to exploit the ROM only if one does not have to build from scratch the reduced model for any parametric instance. For this reason, the system (42) is assumed to be affinely decomposed, i.e. all the variational forms involved can be written as the product of µ−independent forms and µ−dependent functions [27]. When the affine dependency assumption is verified, the online phase does not depend on N , and can be performed in a small amount of time. Conversely, the offline process is performed only once and can take advantage of High Performance Computing (HPC) resources.
Remark 4.1. Our test cases deal with Navier-Stokes governing equations (21), then, it has at most quadratically nonlinear terms and the affine decomposition is not fulfilled. One can employ hyperreduction techniques as the Empirical Interpolation Method (EIM) to recover it. We refer the interested reader to [7] or [27,Chapter 5].
In the next Section, we will focus on the ROM offline and online phase, showing the strategy employed to build the reduced function spaces.

Offline and online stages.
In this Section, we present how to build a reduced space for OCP(µ)s. We exploit here the POD algorithm. Thanks to this algorithm, N max snapshots are sampled and then compressed in order to generate function spaces of dimension N < N max .
It is well known that optimization governed by PDEs(µ) constraints leads to the solution of a saddle point system [9,13,28,53], as already specified in Section 1.2. In order to guarantee the well-posedness of such a structure, the matrix B of system (13) must verify the inf-sup stability condition (15) for every µ ∈ P. In the FE approximation, the above-mentioned relation holds since state and adjoint spaces are equally discretized. However, the inf-sup stability must hold at the reduced level too, since the relation is provable if the reduced spaces for state and adjoint variables coincide. The standard POD construction process leads to the reduced function spaces for state and adjoint which may be different even under the assumption of the same starting FE spaces. To overcome this issue, the basis are usually manipulated in order to stabilize the system. Indeed, we apply aggregated spaces technique, as already done in several papers about ROM for OCP(µ)s, see [4,5,20,22,29,30,40,41,49] as references. The strategy aims at building a common space for state and adjoint which is able to describe both state and adjoint variables. Let us suppose to have applied a standard POD for all the involved variables and to have defined the following spaces Y N = span {χ y n , χ z n , n = 1, . . . , N } and U N = span{χ u n , n = 1, . . . , N }, with where Z y ≡ Z z = [χ y 1 | · · · |χ y N |χ z 1 | · · · |χ z N ] ∈ R Ny×2N and Z u = [χ u 1 | · · · |χ u N ] ∈ R Nu×N are the reduced basis matrices for each variable and Z spans the global space X N . We want to solve the optimality system in a low dimensional framework at each parametric instance. To this end, we employ a Galerkin projection into the reduced spaces and the system (6) will be (43) G (ZX N ; µ), and F N := Z T F.
The system (43) is nonlinear, thus we apply Newton's method and we iteratively obtain from the FE approximation, the Fréchet derivative inherits the saddle point structure, i.e.
We now have all the ingredients to define a reduced Brezzi inf-sup condition as follows (46) β Br,N (µ) := inf If µ = µ * , relation (46) is verified thanks to the aggregated space definition. We remark that this technique is increasing the dimension of the global reduced system. However, it is usually still much smaller then N . For the sake of simplicity and for a consistent construction of state and adjoint space, we always choose N max and N equal for all the involved variables.
Remark 4.2 (Supremizer Stabilization). Dealing with Navier-Stokes governing equations, one has to take care not only of the global inf-sup condition (46), but also with the state equation infsup condition. Indeed, Navier-Stokes problem can be recast as saddle point itself, which results in a nested saddle point structure when a Navier-Stokes problem is used as state equation of an optimal control problem. The aggregated space techniques has to be accompanied by supremizer stabilization for the reduced velocity space. This approach [51] consists in defining a supremizer operator T µ : P Np → V Nv as where b(·, ·) is the bilinear form representing the continuity equation defined in Section 2.1. Then, we enrich the reduced velocity space through the pressure supremizers as follows: where χ Tp n and χ Tq n are the basis supremizers obtained by state and adjoint pressure snapshots, respectively. Enlarging in this way the reduced space for velocity will guarantee inf-sup stability for the Navier-Stokes equation. This approach, i.e. supremizer stabilization combined with aggregated spaces, is the key for the well-posedness of the whole optimality system (3). This will lead to a reduced system of dimension 13N , which is still convenient compared to the global FE approximation dimension.   [42,45,32]. The former choice is due to the presence of two multiplier variables, which increase the global ROM dimension (from 13N to 15N ), and jeopardize the robustness of the reduced nonlinear solver. We remark that the final reduced systems are still much smaller than their high-fidelity counterparts, which involve from 50 to 70 thousands of degrees of freedom, depending on the type of control imposed, boundary or distributed, respectively. We perform an online phase solving (42) for 151 equidistant value of µ in the same parameter space P. The performance has been tested through separate error analysis for each variable. The reliability of the ROM approach has been evaluated through • an average error over the parameter space against an increasing value of the reduced spaces dimension N from one up to N ; • a µ-dependent error computed for the value N . The two error analyses highlight different features of the reduced system that we are going to discuss in the following. Indeed, the average error gives us information about how the reliability changes as the behavior of the solution changes. The straight profile appears to be always the best approximated, due to its Stokes-like (symmetric) nature for all the value µ ∈ P. This is the case of Neumann and Channel control, which average error is depicted in Figures 24a and 28a. Their asymmetric counterparts, Figures 24b and 28b, show how representing the two different features of the solution, a Stokes-like one for µ > µ * and a wall-hugging for the lower values of µ, using the same value of N is more difficult than the symmetric one. Nonetheless, the provided accuracy for basis size N is satisfactory for many practical applications in both target cases. This argument applies to the control and adjoint variables, yet the state ones are the best described by ROM for all the test cases. Because of the optimality equation, the adjoint variables feel the direct influence of the control, which is the most challenging one to be approximated by the reduced model due to its high variability in µ. Indeed, the control variable presents a sort of on-off behavior which drastically affects the efficiency of reduced representation.
For example, if we deal with Stokes target v d , the control is off for high values of viscosity, but, when µ ∼ µ * , it starts to grow in magnitude and to change drastically its features. This is represented in Figures 25b and 31b, where the higher values for the control error and, thus, for the adjoint variable, is for higher values of µ. Since the control magnitude is essentially zero, for Channel and Dirichlet test cases with low Reynolds number (high viscosity), we chose to plot the absolute errors instead of the relative ones, in order to prevent division by zero. This is not  the case of Distributed control, see for example Figure 26b, which presents good errors decay for all the variables, since its strong action causes the control magnitude to always be a meaningful normalization factor. The most challenging case to be approximated is the Dirichlet one for α = 0.001. It is the most complex dynamics, where a new bifurcation appears. Indeed, it can reach an error of almost 10 −3 for the controlled state. Even though this result is worse than the ones obtained for the other test cases, where average error are ranging between 10 −5 and 10 −8 , the accuracy provided by the ROM is  still acceptable for many practical purposes. We remark that this performance is strictly correlated to the more complex features of the Dirichlet problem. This is evident also from the µ-dependent errors in Figure 31a and 31b. In both the pictures, we see a great increment of the error for high Re. Although the phenomenon appears also for the other test cases, see Figures 27b, 29a and 29b, it is not as strong as the Dirichlet case. Furthermore, the µ-dependent error gives an a posteriori information about the bifurcation point. Indeed, in order to have good accuracy property, ROM approach requires regularity on the parametric dependence of the solution instead of the spatial one. This means that reduced errors will generally exhibit a peak at µ * . In fact, an increasing value of the error can be seen around µ * ∼ 0.96 for Neumann, Distributed and Channel control as one can observe from Figures 25a, 27a and 29a for the state variable, respectively.
This feature can be very useful when there is no previous knowledge about bifurcating behaviors. In this sense, ROM is not only an useful approach to solve in a fast way very complicated time consuming systems, but also to detect parameters which can be related to the bifurcating nature of the problem at hand, since their instances will be the worst approximated. Namely, the ROMs confirm "a posteriori" the location of the bifurcation points. We conclude this analysis by noticing that the same considerations hold for the Dirichlet control, but this time at the left end of the parametric domain P where such phenomenon is clearly influenced by the new configuration observed in Figure 20b.

Conclusions
In this work we proposed a first attempt to steer bifurcation phenomena arising from nonlinear PDE(µ) through an optimal control formulation. First of all, we built a general framework which can be applied to general nonlinear OCP(µ)s and we proposed a global stability analysis trough the solution of the eigenvalue problem associated to the optimization system. We tested the stability analysis over four control problems governed by bifurcating Navier-Stokes equations in a suddenexpansion channel. We studied how the control can affect the classical stable wall-hugging solution of the uncontrolled state equation and we proposed some observations on different configurations and features of the optimal solution based on the eigenvalue systems. Furthermore, we employed ROM for all the test cases, proposing it as a strategy to solve the parametrized analysis of the optimization system in a low-dimensional setting, while confirming the applicability of reduction strategies to complex nonlinear models. We believe that this work is a first step towards a better comprehension of the very complex action of optimal control over a nonlinear bifurcating system. We believe that the content of this work could pave the way to many improvements on the topic. Among them, one could be the deeper analysis of the Dirichlet test case, which gives the more unexpected and, consequently, difficult to interpret, results.
with Applications in Computational Fluid Dynamics". We also acknowledge the PRIN 2017 "Numerical Analysis for Full and Reduced Order Methods for the efficient and accurate solution of complex systems governed by Partial Differential Equations" (NA-FROM-PDEs) and the INDAM-GNCS project "Tecniche Numeriche Avanzate per Applicazioni Industriali". The computations in this work have been performed with RBniCS [1] library, developed at SISSA mathLab, which is an implementation in FEniCS [37] of several reduced order modelling techniques; we acknowledge developers and contributors to both libraries.