Mean-field optimal control for biological pattern formation
MMEAN-FIELD OPTIMAL CONTROL FOR BIOLOGICAL PATTERNFORMATION
MARTIN BURGER, LISA MARIA KREUSSER, AND CLAUDIA TOTZECK
Abstract.
We propose a mean-field optimal control problem for the parameter iden-tification of a given pattern. The cost functional is based on the Wasserstein distancebetween the probability measures of the modeled and the desired patterns. The first-order optimality conditions corresponding to the optimal control problem are derivedusing a Lagrangian approach on the mean-field level. Based on these conditions we pro-pose a gradient descent method to identify relevant parameters such as angle of rotationand force scaling which may be spatially inhomogeneous. We discretize the first-orderoptimality conditions in order to employ the algorithm on the particle level. Moreover,we prove a rate for the convergence of the controls as the number of particles used forthe discretization tends to infinity. Numerical results for the spatially homogeneous casedemonstrate the feasibility of the approach. Introduction
In the past years, interacting particle or agent systems have been widely used to modelcollective behavior in biology, sociology and economics. Among the many examples ofapplications are biological phenomena such as animal herding or flocking [7, 8, 13], cellmovement [18], as well as sociological and economical processes like opinion formation [15],pedestrian flow dynamics [6, 27], price formation [4], robotics [26] and data science [22].Most of these models start from interacting particle systems and encode ‘first principles’ ofbiological, social or economical interactions, inspired by Newtonian physics with its classicaldynamical systems of first or second order. A key feature is the formation of global patternseven if the agents only interact with each other at a local scale. These patterns includeconsensus, polarisation and clustering. The qualitative results on pattern formation hasbeen achieved by intensive studies of the limit for infinitely many particles.Interacting particle models have been extended to include control actions. The controlof the dynamics has recently become an active research area [1, 9, 12, 17, 24, 25, 28]. Theimpact of the control on pattern formation has been studied both on the level of agentsas well as in the mean-field limit, and has been applied successfully to a wide range ofapplications including traffic flow [20] and herds of animals [7, 8].While optimal control has mainly been studied for isotropic interacting particle modelssuch as the Cucker-Smale model [24], this paper focuses on control actions for a class ofagent-based models with anisotropic interaction forces. For a large number N of interactingcells with positions x j = x j ( t ) ∈ R , j = 1 , . . . , N , at time t and interaction forces G , weconsider models of the form d x j d t = 1 N N X k =1 k = j G ( x j − x k , T ( x j )) , (1)equipped with initial data x j (0) = x in j , j = 1 , . . . , N , for given x in j ∈ R , j = 1 , . . . , N .Here, T denotes an underlying tensor field influencing the interaction force F in additionto the distance vector x j − x k . This tensor field is responsible for anisotropic patternformation and (1) can be regarded as a prototype for understanding complex phenomena innature. An example of this class of models is the K¨ucken-Champod model [5, 23] describingthe formation of fingerprint patterns based on the interaction of N cells. The interactingparticle model (1) has been studied in [5, 10, 11, 14]. In particular, the particles align inline patterns. For spatially homogeneous T , the stationary solution is given by straightline patterns, while for general T , more complex line patterns are observed in numericalsimulations as stationary solutions. a r X i v : . [ m a t h . O C ] S e p M. BURGER, L. M. KREUSSER, AND C. TOTZECK
Fingerprint identification algorithms are of great importance in forensic science and areincreasingly used in biometric applications. Due to data protection and privacy regulations,many of these algorithms are usually developed based on realistic synthetic fingerprint im-ages which motivates the simulation of realistic patterns based on biological models. Usingparticle models like (1) fingerprint patterns can be produced as stationary solution wherethe adjustment of one of the model parameters is related to the distances between the fin-gerprint lines. The motivation of this paper is to produce patterns with desired features,such as certain angles of rotation or scalings. We therefore propose to consider an optimalcontrol problem constrained by (1). This allows us to estimate the model parameters fora given desirable stationary pattern. In this framework, we focus on the estimation of theangle of rotation of the pattern and the strength of the interaction force at each point of thegiven domain. The spatially inhomogeneous control problem considered here is an inverseproblem and can be regarded as the first step towards modeling complex fingerprint patternswith specific features in the future.The cost functional used for the parameter identification measures the Wasserstein dis-tance of the current states and the desired states and penalizes parameter settings that aretoo far from some given reference states. The algorithm proposed for the numerical sim-ulations is based on projected gradient descent methods, where the gradient is computedusing an adjoint approach. The first-order optimality conditions, required for evaluating thegradient, are derived with the help of a Lagrangian ansatz, similar to the one in [7]. Choos-ing the first-optimize-then-discretize approach pays off in the numerical implementation. Infact, we can use different time discretizations for the forward and the adjoint solver. Thisallows us to save a lot of effort, as the forward solver is implemented with an explicit Eulerscheme and the linear and stiff adjoint system is solved implicitly.The main novelties of our approach are the modeling of a spatially inhomogeneous mean-field optimal control problem with periodic boundary conditions, and the treatment of theWasserstein cost. The latter is, in particular, challenging from the numerical point of view.We show some simulation results on the particle level to demonstrate the feasibility of ourapproach and show a rate for the convergence of the controls as the number of cells tendsto infinity.This paper is organised as follows. In Section 2, we introduce the mean-field model andthe considered spatially inhomogeneous optimal control problem in the macroscopic setting.The first-order optimality conditions are derived in Section 3. In Section 4, we derive thediscrete optimal control problem and the first-order optimality conditions by discretizing theforward and adjoint models. We also show the convergence of the discrete optimal controlproblem. The derived discretizations of the problem form the basis for our numerical schemesand we describe the resulting algorithm in Section 5. Numerical simulation results for thespatially homogeneous case are shown in Section 6 before we conclude the paper in Section 7.2.
Description of the problem
In this section, we introduce an interacting particle model on the torus that includescontrol actions. Starting from (1), we introduce the considered interaction forces first, thenwe pass to the continuum model and formulate the optimal control problem.2.1.
Interaction forces.
For formulating an optimal control problem, we introduce spa-tially inhomogeneous interaction forces depending on a control variable u . We introducea Hilbert space H ( R ) of real-valued functions on R and require that it is continuouslyembedded in L ∞ ( R ) and C ( R ). An example for H ( R ) is given by the Sobolev space H ( R ) which is continuously embedded in L ∞ ( R ) and C ( R ) in two dimensions. We de-fine the space of controls as U = H ( R ) × H ( R ) and consider controls u = ( θ, η ) ∈ U with θ ∈ H ( R ) and η ∈ H ( R ). To introduce a control in (1), we replace the force G dependingon T in (1) by a force F depending on a control u ∈ U . This results in an interacting particlemodel of the form d x j d t = 1 N N X k =1 k = j F ( x j − x k , u ( x j )) . (2) EAN-FIELD OPTIMAL CONTROL FOR BIOLOGICAL PATTERN FORMATION 3
A typical aspect of aggregation models is the competition of social interactions (repulsionand attraction) between particles and thus we assume that the force F is of the form F ( d = d ( x, y ) , u ( x )) = F A ( d, u ( x )) + F R ( d, u ( x )) , (3)where d = d ( x, y ) = x − y . Here, F R ( d, u ( x )) denotes the repulsion force that a particle atlocation y exerts on particle at location x subject to the control parameter u ( x ), and F A isthe attraction force a particle at location y exerts on particle at location x , again subject tothe control parameter u ( x ). We assume that F R and F A are of the form F R ( d = d ( x, y ) , u ( x )) = η ( x ) f R ( η ( x ) | d | ) d (4)and F A ( d = d ( x, y ) , u ( x )) = η ( x ) f A ( η ( x ) | d | ) R θ ( x ) (cid:18) χ (cid:19) R Tθ ( x ) d, (5)respectively, where the spatially inhomogeneous control parameter u ( x ) is defined as u ( x ) =( θ ( x ) , η ( x )). Here, we consider χ ∈ [0 ,
1] and radially symmetric coefficient functions f R and f A , where, again, d = d ( x, y ) = x − y ∈ R . The rotation matrix R θ ( x ) is defined as R θ ( x ) = (cid:18) cos( θ ( x )) − sin( θ ( x ))sin( θ ( x )) cos( θ ( x )) (cid:19) . The unusual form of the attraction force F A is motivated by [5, 10, 11, 14, 23] where thedirection of the interaction force depends on a spatially homogeneous or inhomogeneoustensor field T = T ( x ) with T ( x ) := χs ( x ) ⊗ s ( x ) + l ( x ) ⊗ l ( x ) ∈ R , for orthonormal vector fields s = s ( x ) and l = l ( x ) ∈ R . Writing s, l in polar coordinatesresults in s ( x ) = ( − sin( θ ( x )) , cos( θ ( x ))) , l ( x ) = (cos( θ ( x )) , sin( θ ( x )))and the tensor field T is given by T ( x ) = R θ ( x ) (cid:18) χ (cid:19) R Tθ ( x ) . This expression occurs in the definition of the attraction force F A in (5).The parameter χ introduces an anisotropy to the force F if χ <
1. The force F along l ( x ) = (cos( θ ( x )) , sin( θ ( x ))) is independent of χ and we assume that F is short-range re-pulsive, long-range attractive along l . This implies that F is also short-range repulsive,long-range attractive along s ( x ) = ( − sin( θ ( x )) , cos( θ ( x ))) for χ = 1, while for χ = 0 thetotal force F along s is purely repulsive.For the forces F R and F A , we make the following assumption: Assumption 1.
We assume that the force coefficients f R and f A are bounded with f R ( | d | ) → , f A ( | d | ) → as | d | → . . (6) Further, f R , f A are continuously differentiable, implying that the partial derivatives of F R and F A are bounded, i.e. sup d ∈ R (cid:12)(cid:12)(cid:12)(cid:12) ∂F R ∂d (cid:12)(cid:12)(cid:12)(cid:12) < + ∞ , sup d ∈ R (cid:12)(cid:12)(cid:12)(cid:12) ∂F A ∂d (cid:12)(cid:12)(cid:12)(cid:12) < + ∞ and sup η ∈ [ η min ,η max ] (cid:12)(cid:12)(cid:12)(cid:12) ∂F R ∂η (cid:12)(cid:12)(cid:12)(cid:12) < + ∞ , sup η ∈ [ η min ,η max ] (cid:12)(cid:12)(cid:12)(cid:12) ∂F A ∂η (cid:12)(cid:12)(cid:12)(cid:12) < + ∞ . Note that these conditions are satisfied for the exponentially decaying force coefficients in[5, 10, 11, 14, 23]. The rather unusual assumption (6) is considered to guarantee physicallyrelevant forces through periodic extension on the torus which will be introduced in thefollowing.Motivated by the numerical simulations in [5, 14], the aim of this work is to study (2) onthe torus and we consider Ω = T ⊂ R . Since the torus T can be associated with the unitsquare [0 , with periodic boundary conditions, it is useful to consider periodically defined M. BURGER, L. M. KREUSSER, AND C. TOTZECK forces for the associated discretized problems. Hence, we assume that ¯ F : R × U → R isthe periodic extension of some force F , defined by¯ F ( d + k, u ( x )) = F ( d, u ( x )) , x ∈ R , d = d ( x, y ) ∈ [ − . , . , k ∈ Z , (7)for any u ∈ U , see [11] for more details. As the solutions are very sensitive to the scalingparameter η, we restrict the space of controls U to the space of admissible controls U ad ⊂ U ,defined as U ad = { u = ( θ, η ) ∈ U : η ∈ [ η min , η max ] a.s. } (8)for 0 < η min < η max . The force F and its periodic extension ¯ F are Lipschitz continuouswith respect to their first variable with Lipschitz constant C d , i.e. | F ( d, u ( x )) − F ( ¯ d, u ( x )) | ≤ C d | d − ¯ d | for all u ∈ U ad where the Lipschitz constant C d is independent of u = ( θ, η ) due to theboundedness of η . Clearly, the partial derivatives with respect to θ are bounded and hence,due to Assumption 1, there exists a Lipschitz constant C u such that | F ( d, u ( x )) − F ( d, ¯ u ( x )) | ≤ C u | u ( x ) − ¯ u ( x ) | for all d ∈ R .2.2. State model.
Before formulating the state model on Ω = T , we introduce our nota-tion. Let P (Ω) denote the space of Borel probability measures on Ω. By P (Ω), we denotethe space of Borel probability measures on Ω with finite second moments, endowed withthe 2-Wasserstein distance W , and by P ac (Ω) ⊂ P (Ω) we denote the space of Borel prob-ability measures with finite second moments which are absolutely continuous with respectto the Lebesgue measure. For any map g : Ω → R and any measure µ : Ω → [0 , + ∞ ] wedenote by g µ : R → [0 , + ∞ ] the pushforward measure, defined by g µ ( B ) = µ ( g − ( B ))for any Borel set B ⊂ R . In the following, we consider the restriction of U to Ω, given by U = H (Ω) × H (Ω) and U ad as defined in (8).Next, we formulate the state model. For this, we consider a control u = ( θ, η ) ∈ U , i.e. θ ( x ) ∈ R and η ( x ) ∈ R for all x ∈ R , and introduce the interacting particle modeld x j d t = 1 N N X k =1 k = j ¯ F ( x j − x k , u ( x j ))(9)for the periodic interaction force ¯ F on the unit torus Ω = T ⊂ R . The associated contin-uum model on Ω is given by ∂ t ρ ( t, x ) + ∇ x · (cid:2) ρ ( t, x )( ¯ F ( · , u ( x )) ∗ ρ ( t, · ))( x ) (cid:3) = 0 in [0 , T ] × Ω(10)for any
T > ρ (0 , · ) = ρ in in Ωfor some given probability ρ in ∈ P (Ω).For ¯ F and u given, (10) has a unique global (weak measure) solution ρ ∈ C ([0 , T ] , P (Ω))for any T >
T > ρ t := ρ ( t, · ) forthe solution of (10) at time t ∈ [0 , T ]. We assume that T > T is close to the corresponding steady state ρ ∞ of (10)satisfying ρ ∞ ( x )( ¯ F ( · , u ( x )) ∗ ρ ∞ )( x ) = 0 , (cid:0) ¯ F ( · , u ( x )) ∗ ρ ∞ (cid:1) ( x ) = Z Ω ¯ F ( x − y, u ( x )) ρ ∞ ( y ) d y. (11)Provided χ ≥ θ ( x ) for x ∈ Ω. Since s ( x ) rotatesanticlockwise as θ ( x ) increases with s ( x ) = (0 ,
1) for θ ( x ) = 0, θ ( x ) is the angle between thedirection of the stationary pattern at x and the vertical axis. Note that rotations θ ( x ) + kπ for any k ∈ Z result in the same direction of the pattern at x , implying that it is suffi-cient to consider θ ( x ) ∈ [0 , π ) for x ∈ Ω. This can also be seen by the fact that we have F A ( d, ( θ ( x ) , η ( x ))) = F A ( d, ( θ ( x ) + π, η ( x ))) for any θ ( x ) ∈ [0 , π ) and hence also for itsperiodic extension. EAN-FIELD OPTIMAL CONTROL FOR BIOLOGICAL PATTERN FORMATION 5
The solution of the mean-field PDE (10) depends on the choice of initial data. In par-ticular, this implies that the solution to (11) is not unique. For concentrated initial data,a single vertical line along s is expected as stationary solution, while for other initial data,more complex patterns may arise as stationary solution. The distance between those linescan be controlled by rescaling the total force which is controlled by the positive function η . Remark 2.1.
The solution to (11) is not unique in general. To see this, note that anysolution ρ ∞ of (11) implies that ρ ∞ + c for any constant c ∈ R is also a solution to (11) since we have R Ω ¯ F ( x − y, u ( x )) d y = 0 for any x ∈ Ω . To guarantee a unique solution, weconsider an approximation of a stationary solution to (11) as the state problem, given by thesolution to the anisotropic aggregation equation (10) after a sufficiently large time T > forspecific initial data. Remark 2.2.
The term ( ¯ F ( · , u ( x )) ∗ ρ )( x ) for ρ ∈ P (Ω) in the state problem (11) canbe regarded as a macroscopic velocity field. We denote the space of Lipschitz continuousfunctions on Ω by Lip(Ω) and we write h· , ·i for the scalar product on R . The velocity field v : P (Ω) × U ad → Lip(Ω) , v ( ρ, u )( x ) = ( ¯ F ( · , u ( x )) ∗ ρ )( x ) satisfies h v ( ρ, u )( x ) − v ( ρ, u )( y ) , x − y i = (cid:28)Z Ω ( ¯ F ( x − z, u ( x )) − ¯ F ( y − z, u ( y ))) ρ ( z ) d z, x − y (cid:29) (12a) ≤ C d | x − y | , x, y ∈ Ω(12b) for all ( ρ, u ) ∈ P (Ω) × U ad where the Lipschitz constant C d > is independent of ( ρ, u ) . For any function g on Ω , we write k g k ∞ = sup x ∈ Ω | g ( x ) | for the supremum norm and wehave for any ( ρ, u ) , (¯ ρ, ¯ u ) ∈ P (Ω) × U ad : k v ( ρ, u ) − v (¯ ρ, ¯ u ) k ∞ ≤ k v ( ρ, u ) − v (¯ ρ, u ) k ∞ + k v (¯ ρ, u ) − v (¯ ρ, ¯ u ) k ∞ . Indeed, we have k v ( ρ, u ) − v (¯ ρ, u ) k ∞ ≤ C v W ( ρ, ¯ ρ ) for some constant C v > and k v (¯ ρ, u ) − v (¯ ρ, ¯ u ) k ∞ = sup x ∈ Ω (cid:12)(cid:12)(cid:12)(cid:12)Z Ω ( ¯ F ( x − z, u ( x )) − ¯ F ( x − z, ¯ u ( x )))¯ ρ ( z ) d z (cid:12)(cid:12)(cid:12)(cid:12) ≤ C u k u − ¯ u k ∞ due to the continuous embedding of H (Ω) in L ∞ (Ω) . This implies k v ( ρ, u ) − v (¯ ρ, ¯ u ) k ∞ ≤ C v W ( ρ, ¯ ρ ) + C u k u − ¯ u k ∞ . (13) For weak solutions ρ, ¯ ρ of (10) with initial data ρ in , ¯ ρ in and controls u, ¯ u , respectively, onecan show that there exists positive constants a, b such that W ( ρ t , ¯ ρ t ) ≤ (cid:0) W ( ρ in , ¯ ρ in ) + b k u − ¯ u k ∞ (cid:1) exp( at )(14) for all t ∈ [0 , T ] where we used again the continuous embedding of H (Ω) in L ∞ (Ω) . Theproof is based on computing the derivative of W ( ρ t , ¯ ρ t ) with respect to t , estimates usingthe inequalities (12) , (13) , and the Gronwall inequality, see [7] for more details. Optimal control problem in the macroscopic setting.
With the state problemdefined in (11), we can now formulate the associated optimal control problem. For ¯ F and u given, we identify the solution of the state problem as the weak solution of the macroscopicproblem (10) at time T > ρ in the following. Given the distributionfunction ρ des ∈ P (Ω) of a desired fingerprint pattern, the task at hand is to find the forcefield ¯ F characterized by u = ( θ, η ) corresponding to the given density ρ des . This task can beregarded as an inverse problem which is mathematically formulated as an optimal controlproblem with a PDE constraint. We define the cost functional J ( ρ T , u ) = 12 W ( ρ T , ρ des ) + λ k θ − θ ref k H (Ω) + λ k η − η ref k H (Ω) (15)where λ , λ > θ ref , η ref ∈ H (Ω) are given reference values. Tosummarize, we obtain: Problem 1.
Find u ∈ U ad such that ( ρ T ( u ) , u ) = arg min ρ T ,u J ( ρ T , u ) subject to (10) . M. BURGER, L. M. KREUSSER, AND C. TOTZECK
Considering the well-defined solution operator S : U → P (Ω) with Su = ρ T associatedwith the unique weak solution of (10), the optimization problem is then given bymin ˜ J ( u )where ˜ J is defined by ˜ J ( u ) := J ( Su, u ) and is called the reduced cost functional.3.
First-order optimality conditions in the macroscopic setting
The main objective of this section is to derive the first-order optimality conditions (FOC)for the optimal control problem by using the Lagrangian approach based on Wassersteincalculus. The arguments of this section are similar to the ones in [7]. Since there exists aunique global (weak measure) solution ρ ∈ C ([0 , T ] , P (Ω)) of the mean-field PDE (10) forany u ∈ U , we define the state operator(16) E ( ρ, u )[ ϕ ] := h ϕ T , ρ T i − h ϕ , ρ in i − Z T h ∂ t ϕ ( t, x ) + ( ¯ F ( · , u ( x )) ∗ ρ t )( x ) · ∇ ϕ ( t, x ) , ρ ( t, x ) i d t with E ( ρ, u ) = 0 for all ϕ ∈ A := C c ([0 , T ] × Ω). Further, let J ( ρ ) := 12 W ( ρ T , ρ des ) , J ( u ) := λ k θ − θ ref k H (Ω) , J ( u ) := λ k η − η ref k H (Ω) . (17)Hence, we have J ( ρ, u ) = J ( ρ ) + J ( u ) + J ( u ). Since the functional J is not handyfor deriving the first-order optimality conditions, we consider the extended Lagrangian I defined bymin ( ρ,u ) I ( ρ, u ) = min u n J ( u )+ J ( u )+min ρ sup ϕ ∈A (cid:8) J ( ρ )+ E ( ρ, u )[ ϕ ] (cid:9)o = min u n J ( u )+ J ( u )+ κ ( u ) o with κ ( u ) = min ρ sup ϕ ∈A n J ( ρ ) + E ( ρ, u )[ ϕ ] o . Note that E has the property sup ϕ ∈A E ( ρ, u )[ ϕ ] ≥ ϕ ≡ E ( ρ, u )[0] = 0for every ( ρ, u ) . Therefore, if E ( ρ, u )[ ϕ ] > ϕ , the linearity in ϕ of E yields E ( ρ, u )[ aϕ ] = aE ( ρ, u )[ ϕ ] for every a > ϕ E ( ρ, u )[ ϕ ] = + ∞ . We say a pair ( ρ, u ) ∈ C ([0 , T ] , P (Ω)) × U is admissible if E ( ρ, u )[ ϕ ] = 0 for all ϕ ∈ A . Inthe following, we derive a necessary condition for an admissible pair ( ρ, u ) to be a stationarypoint. For this, we consider perturbations of the admissible pair ( ρ, u ), given by an admissibleperturbation u δ = u + δh ∈ U of u for δ ≥ h = ( θ, η ) ∈ U , and the associated ρ δ ∈ P (Ω) satisfying E ( ρ δ , u δ )[ ϕ ] = 0 for all ϕ ∈ A . These conditions result in thefollowing assumption: Assumption 2.
Let ( ρ, u ) be an admissible pair and suppose that h ∈ U . Suppose ¯ δ > isgiven such that for any ≤ δ ≤ ¯ δ the perturbation u δ := u + δh satisfies • u δ is an admissible control, i.e. u δ ∈ U . • There exists ρ δ ∈ C ([0 , T ] , P (Ω)) such that E ( ρ δ , u δ )[ ϕ ] = 0 for all ϕ ∈ A . For an admissible pair ( ρ, u ) satisfying Assumption 2 we have κ ( u δ ) = min ρ sup ϕ ∈A n J ( ρ ) + E ( ρ, u δ )[ ϕ ] o = J ( ρ δ )= J ( ρ δ ) − J ( ρ ) + min ρ sup ϕ ∈A n J ( ρ ) + E ( ρ, u )[ ϕ ] o = J ( ρ δ ) − J ( ρ ) + κ ( u ) , and the directional derivative of G := J + J + κ at u along h is given bylim δ → G ( u δ ) − G ( u ) δ = lim δ → [ J ( ρ δ ) − J ( ρ )] + [ J ( u δ ) − J ( u )] + [ J ( u δ ) − J ( u )] δ , which depends on the relationship between ρ δ and ρ . Remark 3.1.
Using estimate (14) , we obtain W ( ρ δt , ρ t ) ≤ δ √ b k h k ∞ exp (cid:18) aT (cid:19) EAN-FIELD OPTIMAL CONTROL FOR BIOLOGICAL PATTERN FORMATION 7 for any t ∈ [0 , T ] . Hence, for t ∈ [0 , T ] given, the curve [0 , ∞ ) δ ρ δt ∈ P (Ω) whichstarts from ρ t at δ = 0 is absolutely continuous with respect to the 2-Wasserstein distanceand there exists ψ t ∈ L ( ρ t ; Ω) satisfying [2, Proposition 8.4.6](18) lim δ → W ( ρ δt , ( id + δψ t ) ρ t ) δ = 0 . Furthermore, W (( id + δψ t ) ρ t , ρ t ) ≤ Z Z Ω × Ω | x + δψ t ( x ) − x | d ρ t ( x ) = δ Z Z Ω × Ω | ψ t ( x ) | d ρ t ( x ) . In particular, we have that lim sup δ → W ( ρ δt , ρ t ) δ = lim sup δ → W (( id + δψ t ) ρ t , ρ t ) δ ≤ sZ Z Ω × Ω | ψ t | d ρ t ( x ) . In the following we establish an explicit relationship between the perturbations ψ t and h as in [7, Thm 3.4]. Note that we omit the proofs of the following lemma and the theorembecause they are very similar to the ones in the reference. Lemma 3.1.
Let ( ρ, u ) be an admissible pair and let δ > sufficiently small, h ∈ U and u δ = u + δh such that (i) u δ ∈ U , and (ii) there exists ρ δ ∈ C ([0 , T ] , P ac2 (Ω)) satisfying E ( ρ δ , u δ ) = 0 .Suppose that ψ ∈ C ((0 , T ) × Ω) with ψ ≡ satisfies ∂ t ψ t + Dψ t v ( ρ t , u ) = K ( ρ t , u )[ ψ t , h ] for ρ t -almost every x ∈ Ω , (19) where ( t, x )
7→ K ( ρ t , u )[ ψ t , h ]( x ) is a bounded Borel map satisfying lim δ → Z T Z Ω (cid:12)(cid:12)(cid:12)(cid:12) v ( ν δt , u δ ) ◦ ( id + δψ t )( x ) − v ( ρ t , u )( x ) δ − K ( ρ t , u )[ ψ t , h ]( x ) (cid:12)(cid:12)(cid:12)(cid:12) d ρ t ( x ) d t = 0(20) with ν δt := ( id + δψ t ) ρ t . Then, (18) holds for this ψ , i.e. lim δ → W ( µ δt , ( id + δψ t ) µ t ) δ = 0 . Remark 3.2.
Note that for any h ∈ U we obtain by Taylor expansion: (21) K ( ρ, u )[ ψ, h ] = Dv ( ρ, u ) ψ − Z Ω ( ∇ x ¯ F )( · − y, u ) ψ ( y ) d ρ ( y ) − ∇ u v ( ρ, u ) h + O ( δ ) , and (19) may be written as ∂ t ψ t + { ψ t , v ( ρ t , u ) } = − Z Ω ( ∇ x ¯ F )( · − y, u ) ψ t ( y ) d ρ t ( y ) − ∇ u v ( ρ t , u ) h, where {· , ·} denotes the Lie bracket given by { ξ, ψ } = ( Dξ ) ψ − ( Dψ ) ξ for vector fields ξ, ψ . The existence of ψ ∈ C b ((0 , T ) × Ω) satisfying the assumptions of Lemma 3.1 is providedin the following statement.
Theorem 3.1.
Suppose that the assumptions of Lemma 3.1 hold. For the velocity field v : P (Ω) × U → Lip loc (Ω) given by v ( ρ, u )( x ) = ( F ( · , u ( x )) ∗ ρ )( x ) there exists ψ ∈ C b ((0 , T ) × Ω) with ψ = 0 satisfying ∂ t ψ t + Dψ t v ( ρ t , u ) = K ( ρ t , u )[Ψ t , h ] for ρ t dt -almost every ( t, x ) ∈ (0 , T ) × Ω , where K is given by (21) . Now, we are able to state the first-order necessary condition for ( ρ, u ) to be a stationarypoint. For this, the Gˆateaux-derivatives of J , J and J in (17) are required. By [2,Proposition 8.5.2], we have δ ρ J ( ρ T )( x ) = t ρ des ρ T ( x ) − x (22)where t ρ des ρ T ( x ) is the unique optimal transport map from ρ T to ρ des . For u = ( θ, η ) ∈ U and h = ( h θ , h η ) ∈ U , we haved J ( u )[ h ] + d J ( u )[ h ] = λ h h θ , θ − θ ref i H (Ω) + λ h h η , η − η ref i H (Ω) (23) M. BURGER, L. M. KREUSSER, AND C. TOTZECK for the inner product h· , ·i H (Ω) on H (Ω) since the Gˆateaux-derivatives are given byd k θ − θ ref k H (Ω) [ h θ ] = (cid:28) h θ , θ − θ ref k θ − θ ref k H (Ω) (cid:29) H (Ω) for θ − θ ref = 0 , d k η − η ref k H (Ω) [ h η ] = (cid:28) h η , η − η ref k η − η ref k H (Ω) (cid:29) H (Ω) for η − η ref = 0 . Theorem 3.2.
Let (¯ ρ, ¯ u ) be an optimal pair, J , J , J be Gˆateaux-differentiable. Supposethat h ∈ U and assume that there exists ¯ δ > such that for any ≤ δ ≤ ¯ δ it holds ¯ u + δh ∈ U .Then, d J (¯ u )[ h ] + d J (¯ u )[ h ] + Z Ω h δ ρ J (¯ ρ T ) , ψ T i d¯ ρ T = lim δ → G (¯ u + δh ) − G (¯ u ) δ = 0(24) where J , J , J are defined in (22) , (23) and t ψ t ∈ L (¯ ρ t ) satisfies (19) with ψ = 0 . The optimality condition (24) together with (22) and (23) is implicit. In order to formu-late an optimization algorithm, we derive an explicit form of the optimality conditions bycomputing the adjoint. Thus, to derive the adjoint-based first-order optimality system, webegin with the dual problem corresponding to (19). It can be obtained by testing (19) witha family of vector-valued measures ( m t ) t ∈ (0 ,T ) with m t ∈ P (Ω; R ), resulting in Z T Z Ω (cid:16) ∂ t ψ t + Dψ t v (¯ ρ t , ¯ u ) − K (¯ ρ t , ¯ u )[ ψ t , h ] (cid:17) · d m t d t = 0or equivalently Z T Z Ω (cid:18) ∂ t ψ t + { ψ t , v (¯ ρ t , ¯ u ) } + Z Ω ( ∇ x ¯ F )( · − y, ¯ u ) ψ t ( y ) d¯ ρ t ( y ) + ∇ u v (¯ ρ t , ¯ u ) h (cid:19) · d m t d t = 0 . Integrating by parts and using ψ = 0, we obtain Z T h ∂ t m t + ∇ · ( v (¯ ρ t , ¯ u ) ⊗ m t ) + ∇ v (¯ ρ t , ¯ u ) m t − ¯ ρ t Z Ω ( ∇ x ¯ F )( y − · , ¯ u ) d m t ( y ) , ψ t i d t = Z Ω ψ T · d m T + Z T Z Ω ∇ u v (¯ ρ t , ¯ u ) h · d m t d t. We choose ¯ m t such that the dual problem ∂ t ¯ m t + ∇ · ( v (¯ ρ t , ¯ u ) ⊗ ¯ m t ) + ∇ v (¯ ρ t , ¯ u ) ¯ m t − ¯ ρ t Z Ω ( ∇ x ¯ F )( y − · , ¯ u ) d ¯ m t ( y ) = 0subject to the terminal condition ¯ m T = ¯ ρ T δ ρ J (¯ ρ T ). By using the optimality condition(24), we obtain − d J (¯ u )[ h ] − d J (¯ u )[ h ] + Z T Z Ω ∇ u v (¯ ρ t , ¯ u ) h · d ¯ m t d t = 0for any h ∈ U . The optimality condition can be summarised as: Theorem 3.3.
A minimizing pair (¯ ρ, ¯ u ) satisfies the state problem ∂ t ¯ ρ t + ∇ · (¯ ρ t v (¯ ρ t , ¯ u )) = 0 subject to the initial condition ¯ ρ (0 , · ) = ¯ ρ = ρ in , and the optimality condition d J (¯ u ) + d J (¯ u ) = Z T Z Ω ∇ u v (¯ ρ t , ¯ u ) d ¯ m t d t, (25) where ¯ m t solves the adjoint equation given by ∂ t ¯ m t + ∇ · ( v (¯ ρ t , ¯ u ) ⊗ ¯ m t ) + ∇ v (¯ ρ t , ¯ u ) ¯ m t − ¯ ρ t Z Ω ( ∇ x ¯ F )( y − · , ¯ u ) d ¯ m t ( y ) = 0 subject to the terminal condition ¯ m T = ¯ ρ T δ ρ J (¯ ρ T ) . Discretization of the optimality conditions
In this section, we formally derive the adjoints and the optimality conditions for thediscrete optimal control problem. In addition, we introduce the discretised reduced costfunctional and its gradient which are required for applying adjoint-based descent methodsto solve the discretised control problems numerically.
EAN-FIELD OPTIMAL CONTROL FOR BIOLOGICAL PATTERN FORMATION 9
Discrete adjoint system.
We derive the discrete system in order to formulate analgorithm for numerical simulations. Note that the domain Ω = T can be regarded as theunit square [0 , with periodic boundary conditions. First, we make an ansatz, similar to[7], to obtain an equation for the vector-valued adjoint variable. Indeed, assuming | ¯ m t | (cid:28) ¯ ρ t for every t ∈ [0 , T ] yields the existence of a vector field ¯ ξ t : Ω → R such that ¯ m t = ¯ ξ t ¯ ρ t where ¯ ρ t is the weak solution to the state equation (10). The dual problem allows us toobtain an equation for ¯ ξ t given by ∂ t ¯ ξ t + D ¯ ξ t v (¯ ρ t , ¯ u ) = −∇ v (¯ ρ t , ¯ u ) ¯ ξ t + Z Ω ( ∇ x ¯ F )( y − · , ¯ u ) ξ t ( y ) d¯ ρ t ( y ) , implying that we can write the adjoint equation as ∂ t ¯ ξ t + ∇ ( v (¯ ρ t , ¯ u ) · ¯ ξ t ) = Z Ω ( ∇ x ¯ F )( y − · , ¯ u ) ¯ ξ t ( y ) d¯ ρ t ( y ) . (26)This PDE is equipped with the terminal condition ¯ ξ T = δ ρ J (¯ ρ T ) where δ ρ J (¯ ρ T )( x ) = t ρ des ρ T ( x ) − x by (22). Note that (26) has a unique solution ξ ∈ C ([0 , T ] × Ω) such that ξ ( t ) ∈ Lip b (Ω) for t ∈ [0 , T ]. We omit the proof, as it is very similar to [7, Thm 3.10]. Using¯ m t = ¯ ξ t ¯ ρ t , (25) can be written asd J (¯ u ) + d J (¯ u ) = Z T Z Ω ∇ u v (¯ ρ t , ¯ u ) ¯ ξ t d¯ ρ t d t = Z T Z Ω Z Ω ∇ u ¯ F ( x − y, ¯ u ( x )) ¯ ξ t ( x ) d¯ ρ t ( y ) d¯ ρ t ( x ) d t. (27)Now, we formally derive the associated optimality conditions on the particle level. Forthis, we consider the equation of characteristics for the mean-field PDE (10), given byd x d t = v ( ρ N , ¯ u )( x ) = Z Ω ¯ F ( x − y, ¯ u ( x )) d ρ N ( y )with the empirical measure ρ N ( t, x ) = 1 N N X j =1 δ ( x − x j ( t ))of particle positions x i ∈ R , i = 1 , . . . , N . The periodicity of the boundary conditions inthe macroscopic setting is mimicked in the discrete setting by considering periodic forcesdefined by (7). We emphasize that this is crucial, as otherwise the same two particles mayinteract various times. For a given control u ∈ U , this leads to the particle systemd x i d t = 1 N N X j =1 ¯ F ( x i − x j , u ( x i )) , i = 1 , . . . , N, (28)subject to the initial conditions x i (0) = X i , i = 1 , . . . , N, (29)where X i ∈ Ω , i = 1 , . . . , N, are given independent realizations of random variables withlaw( X i ) = ρ . Note that (28) is identical to (9).We introduce adjoint variables ξ i = ξ ( x i ) ∈ R , i = 1 , . . . , N , and the discrete adjointsystem readsdd t ξ i = 1 N N X j =1 ∇ x ¯ F ( x i − x j , u ( x i )) ξ i − N N X j =1 ∇ x ¯ F ( x i − x j , u ( x i )) ξ j (30)with terminal condition ξ i ( T ) = δ ρ J ( x i ( T )) = t ρ des ρ T ( x i ( T )) − x i ( T ) . First-order optimality conditions.
We introduce the cost functional in the discretesetting. For this, we replace the probability measure ρ des by the empirical measure ρ N des = 1 N N X i =1 δ x des i , where x des i ∈ Ω , i = 1 , . . . , N , are assumed to be given. Then, the discrete analogue J N : R N × U → R to the cost functional (15) in the continuum setting is given by J N (( x ( T ) , . . . , x N ( T )) , u )= N W N N X i =1 δ x i ( T ) , ρ N des ! + λ k θ − θ ref k ∞ + λ k η − η ref k ∞ (31)where ( θ ref , η ref ) ∈ U ad is assumed to be given reference data. Remark 4.1.
Note that we rescale the first term of the cost functional by N . This isnecessary to have well-balanced terms in the discrete Lagrangian below, see (32) . Indeed, as N → ∞ we obtain infinity many terms in the dual accounting for the constraint, whereasthe Wasserstein distance has fixed order. This scaling is also reported in [3, 7] . Denoting x N = ( x , . . . , x N ), the discrete optimal control problem is given by Problem 2.
For N ∈ N fixed, find u N ∈ U ad such that ( x N ( T ) , u N ) = arg min x N ( T ) ,u N J N ( x N ( T ) , u N ) subject to (28) , where x N = x N ( u N ) . We define the state operator e N by e N ( x, u ) = dd t x ( t ) − N (cid:16)P Nj =1 ¯ F ( x i ( t ) − x j ( t ) , u ( x i ( t ))) (cid:17) Ni =1 x (0) − X ! and the discrete state system (28)–(29) can be rewritten as e N ( x, u ) = 0. Its weak formula-tion is given by h e N ( x, u ) , ( ξ, ζ ) i = Z T dd t x ( t ) − N N X j =1 ¯ F ( x i ( t ) − x j ( t ) , u ( x i ( t ))) Ni =1 · ξ ( t ) d t + ( x (0) − X ) · ζ = N X i =1 Z T dd t x i ( t ) − N N X j =1 ¯ F ( x i ( t ) − x j ( t ) , u ( x i ( t ))) · ξ i ( t ) d t + ( x (0) − X ) · ζ where h· , ·i denotes the scalar product on R . For Lagrange multipliers ( ξ, ζ ), the Lagrangiancorresponding to Problem 2 reads(32) L ( x, u, η, ζ ) = J N (( x ( T ) , . . . , x N ( T )) , u ) + h e N ( x, u ) , ( ξ, ζ ) i for N ∈ N fixed. We have d u J N ( x, u )[ h ] = d J ( u )[ h ] + d J ( u )[ h ]and h d u e N ( x, u )[ h ] , ( ξ, ζ ) i = − N N X i =1 N X j =1 Z T ∇ u ¯ F ( x i ( t ) − x j ( t ) , u ( x i ( t ))) h ( x i ( t )) · ξ i ( t ) d t, implyingd J ( u )[ h ] + d J ( u )[ h ] − N N X i =1 N X j =1 Z T ∇ u ¯ F ( x i ( t ) − x j ( t ) , u ( x i ( t ))) h ( x i ( t )) · ξ i ( t ) d t = 0for any h ∈ U . The first-order optimality condition in the rescaled (see Remark 4.1) micro-scopic setting reads: Theorem 4.1.
Let N ∈ N be given and let ( x N , u N ) be an optimal pair. The optimalitycondition corresponding to Problem 2 reads d J ( u N ) + d J ( u N ) = 1 N N X i =1 N X j =1 Z T ∇ u ¯ F ( x i ( t ) − x j ( t ) , u N ( x i ( t ))) ξ i ( t ) d t (33) where ξ satisfies (30) with terminal condition ξ i ( T ) = N δ ρ J ( x i ( T )) = N ( t ρ des ρ T ( x i ( T )) − x i ( T )) . EAN-FIELD OPTIMAL CONTROL FOR BIOLOGICAL PATTERN FORMATION 11
Gradient of the reduced cost functional.
In this section, we introduce the discre-tised reduced cost functional and its gradient. Motivated by the control-to-state operator S and the reduced cost functional ˜ J in the mean-field setting, introduced in Section 2.3, weconsider the control-to-state operator S N : U → R N ; u ( x , . . . , x N ) where ( x , . . . , x N )satisfies the forward particle system (28) on [0 , T ] with initial conditions (29). The reducedcost functional ˜ J N in the discrete setting is then given by ˜ J N ( u ) := J N ( S N ( u ) , u ) where J N is defined in (31). Since the force ¯ F satisfies Assumption 1, we can implicitly obtaind S N ( u ) via 0 = d x e N ( S N ( u ) , u )[d S N ( u )] + d u e N ( S N ( u ) , u ) , i.e. d S N ( u )[ h ] = − (d x e N ( S N ( u ) , u )) − d u e N ( S N ( u ) , u )[ h ] . The Gˆateaux derivative ˜ J N in the direction h = ( h θ , h η ) ∈ U is obtained fromd ˜ J N ( u )[ h ] = h d x J N ( S N ( u ) , u ) , d S N ( u )[ h ] i + h d u J N ( S N ( u ) , u ) , h i = h d u J N ( S N ( u ) , u ) − (d u e N ( S N ( u ) , u )) ∗ (d x e N ( S N ( u ) , u )) −∗ d x J N ( S N ( u ) , u ) , h i . Defining the adjoint variable ξ = ( ξ k ) Nk =1 by(d x e N ( S N ( u ) , u )) ∗ [ ξ ] = − d x J N ( S N ( u ) , u )yields d ˜ J N ( u )[ h ] = h d u J N ( S N ( u ) , u ) + (d u e N ( S N ( u ) , u )) ∗ ξ, h i = d J ( u )[ h ] + d J ( u )[ h ] − N N X i =1 N X j =1 Z T ∇ u ¯ F ( x i ( t ) − x j ( t ) , u ( x i ( t ))) h ( x i ( t )) · ξ i ( t ) d t. Using the variational lemma we can identify the gradient as ∇ ˜ J N ( u ) = d J ( u ) + d J ( u ) − N N X i =1 N X j =1 Z T ∇ u ¯ F ( x i ( t ) − x j ( t ) , u ( x i ( t ))) ξ i ( t ) d t. (34)With the gradient of the reduced cost functional (34) at hand, we have everything requiredto state the gradient descent algorithm used for the computation of optimal controls.4.4. Convergence of the discrete optimal control problem.
As a first step towardsthe convergence of the discrete optimal control problem, we consider a stability estimate ofthe solutions to the discrete and the continuous adjoint problems (26) and (30). Similarlyas in [7], one can show the following stability estimate for the adjoint solution:
Lemma 4.1.
Let x N = ( x , . . . , x N ) be the solution to the forward particle system (28) with initial condition (29) and given control u N ∈ U ad . Let ρ N ( t, · ) denote the empiricalmeasure corresponding to x N ( t ) for any t ∈ [0 , T ] . Let ρ ∈ C ([0 , T ] , P (Ω)) be the solutionto the mean-field state problem (10) for given control ¯ u ∈ U ad . Let ξ N = ( ξ , . . . , ξ N ) denote the solution to the discrete adjoint system (30) for the pair ( x N , u N ) and supposethat ¯ ξ ∈ C ([0 , T ] , Lip b (Ω)) satisfies (26) for ( ρ, ¯ u ) . Then, there exist positive constants a and b , independent of N ∈ N such that sup t ∈ [0 ,T ] N N X i =1 | ξ i ( t ) − ¯ ξ t ( x i ( t )) | ≤ b exp( aT ) Z T W ( ρ N ( s, · ) , ρ ( s, · )) + k u N − ¯ u k ∞ d s. Denoting by ρ N (0 , · ) the empirical measure which corresponds to the particle locationsin x N ( t ) at time t and using (14), we obtainsup t ∈ [0 ,T ] N N X i =1 | ξ i ( t ) − ¯ ξ t ( x i ( t )) | ≤ C T (cid:0) W ( ρ N (0 , · ) , ρ in ) + k u − ¯ u k ∞ (cid:1) (35)for some constant C T , depending on T > N ∈ N . Theorem 4.2.
Let (¯ ρ, ¯ u ) and ( x N , u N ) be the optimal pairs for Problem 1 and Problem 2with initial data ρ in ∈ P (Ω) and X N = ( X , . . . , X N ) ∈ Ω N , respectively. Let ρ N (0 , · ) denotethe empirical measure corresponding to the initial configuration X N . Then, there exists aconstant c > depending only on T , ¯ ξ and ¯ F such that for λ , λ > c in J , J , definedin (23) , it holds k u N − ¯ u k ∞ ≤ c min { λ , λ } − c W ( ρ N (0 , · ) , ρ in ) . Proof.
Let ξ N = ( ξ , . . . , ξ N ) denote the solution to the discrete adjoint system (30) forthe pair ( x N , u N ) and suppose that ¯ ξ ∈ C ([0 , T ] × Ω) satisfies (26) for ( ρ, u ). We denotethe empirical measure by ρ Nt which corresponds to the particle locations in x N ( t ) at time t . Considering the optimality conditions (27) and (33) in the macroscopic and microscopicsetting, we have for h N := u N − ¯ u ∈ U (d J ( u N ) − d J (¯ u ))[ h N ] + (d J ( u N ) − d J (¯ u ))[ h N ]= 1 N N X i =1 N X j =1 Z T ∇ u ¯ F ( x i ( t ) − x j ( t ) , u N ( x i ( t ))) h N ( x i ( t )) · ξ i ( t ) d t − Z T Z Ω Z Ω ∇ u ¯ F ( x − y, ¯ u ( x )) h N ( x ) · ¯ ξ t ( x ) d¯ ρ t ( y ) d¯ ρ t ( x ) d t = 1 N N X i =1 Z T Z Ω ∇ u ¯ F ( x i ( t ) − y, u N ( x i ( t ))) h N ( x i ( t )) d ρ N ( y ) · ξ i ( t ) d t − Z T Z Ω Z Ω ∇ u ¯ F ( x − y, ¯ u ( x )) h N ( x ) d¯ ρ t ( y ) · ¯ ξ t ( x ) d¯ ρ t ( x ) d t = 1 N N X i =1 Z T Z Ω ∇ u ¯ F ( x i ( t ) − y, u N ( x i ( t ))) h N ( x i ( t )) d ρ Nt ( y ) · ( ξ i ( t ) − ¯ ξ t ( x i ( t ))) d t + Z T Z Ω Z Ω ∇ u ¯ F ( x − y, u N ( x )) h N ( x ) · ¯ ξ t ( x ) d ρ Nt ( x ) d ρ Nt ( y ) d t − Z T Z Ω Z Ω ∇ u ¯ F ( x − y, ¯ u ( x )) h N ( x ) · ¯ ξ t ( x ) d¯ ρ t ( x ) d¯ ρ t ( y ) d t. Using (35), the first term can be estimated by C T (cid:0) W ( ρ N (0 , · ) , ρ in ) + k u N − ¯ u k ∞ (cid:1) for some constant C T depending of T and independent of N . Denoting the optimal couplingbetween ρ Nt and ¯ ρ t by π t , the remaining terms can be rewritten as Z T Z Ω Z Ω × Ω (cid:0) ∇ u ¯ F ( x − y, u N ( x )) h N ( x ) · ¯ ξ t ( x ) − ∇ u ¯ F ( x − y, u N ( x )) h N ( x ) · ¯ ξ t ( x ) (cid:1) d π t ( x, x ) d ρ Nt ( y ) d t + Z T Z Ω Z Ω (cid:2) ∇ u ¯ F ( x − y, u N ( x )) − ∇ u ¯ F ( x − y, ¯ u ( x )) (cid:3) h N ( x ) · ¯ ξ t ( x ) d¯ ρ t ( x ) d ρ Nt ( y ) d t + Z T Z Ω × Ω Z Ω (cid:2) ∇ u ¯ F ( x − y, ¯ u ( x )) − ∇ u ¯ F ( x − y , ¯ u ( x )) (cid:3) h N ( x ) · ¯ ξ t ( x ) d¯ ρ t ( x ) d π t ( y, y ) d t ≤ c T k∇ u ¯ F k ∞ sup t ∈ [0 ,T ] Lip( ¯ ξ t ) + k∇ x ∇ u ¯ F k ∞ sup t ∈ [0 ,T ] k ¯ ξ t k ∞ ! k u N − ¯ u k ∞ W ( ρ Nt , ¯ ρ t )+ k∇ u ∇ u ¯ F k ∞ sup t ∈ [0 ,T ] k ¯ ξ t k ∞ k u N − ¯ u k ∞ + k∇ x ∇ u ¯ F k ∞ sup t ∈ [0 ,T ] k ¯ ξ t k ∞ k u N − ¯ u k ∞ W ( ρ Nt , ¯ ρ t ) ! for some constant c T . This yields(d J ( u N ) − d J (¯ u ))[ h N ] + (d J ( u N ) − d J (¯ u ))[ h N ] ≤ c (cid:0) k u N − ¯ u k ∞ + W ( ρ Nt , ¯ ρ t ) (cid:1) EAN-FIELD OPTIMAL CONTROL FOR BIOLOGICAL PATTERN FORMATION 13 where c depends on T , ¯ ξ and ¯ F . For u N = ( θ N , η N ), ¯ u = (¯ θ, ¯ η ) and h N = ( θ N − ¯ θ, η N − ¯ η ),we have(d J ( u N ) − d J (¯ u ))[ h N ] + (d J ( u N ) − d J (¯ u ))[ h N ] = λ k θ N − ¯ θ k H (Ω) + λ k η − ¯ η k H (Ω) ≥ min { λ , λ }k u N − ¯ u k ∞ by (23) and the continuous embedding of H (Ω) in L ∞ (Ω). This implies(min { λ , λ } − c ) k u N − ¯ u k ∞ ≤ c W ( ρ Nt , ¯ ρ t ) . Choosing λ , λ > c we obtain the desired inequality. (cid:3) Numerical schemes
In this section, we introduce the numerical schemes, used for solving the forward andadjoint initial value problems including the terminal condition of the adjoint problem, aswell as the optimal control problem.5.1.
Forward and adjoint initial value problems.
For solving the discrete forwardproblem (28) with initial condition (29) on the unit square [0 , we apply the simpleexplicit Euler scheme. The discrete adjoint problem (30) is linear and stiff, it is thereforesolved implicitly. Note that the force ¯ F is defined periodically by (7) in both problems.Further, periodic boundary conditions guarantee that the particle positions x i cannot leavethe domain, i.e. x i ∈ [0 , for i = 1 , . . . , N . Remark 5.1.
We emphasize that the first-optimize then discretize approach allows us tochoose different discretizations for the forward and the adjoint solver. This is a huge advan-tage, as otherwise the computational effort increases tremendously due to very small stepsizes or a complicated implementation of the forward solver.
Terminal condition of the adjoint problem.
The main challenge of the implemen-tation of the particle optimization is the evaluation of the terminal condition of the adjoints,given by ξ k ( T ) := ¯ ξ T ( x k ) = N δ ρ J (¯ ρ T )( x k ) = N t ρ des ρ T ( x k ( T )) − x k ( T ) . (36)We realize it with the help of the Python Optimal transport library [16].While ρ T and ρ des are probability densities in the macroscopic setting, we consider theassociated empirical measures ρ N ( T, · ) = N P Nk =1 δ x k ( T ) and ρ N des = N P Nk =1 δ x des k , respec-tively. Let a,b denote the sample weights for the 1D histograms corresponding to x k ( T ) and x des k for k = 1 , . . . , N, i.e. a k = N = b k for k = 1 , . . . , N . Instead of the usual ground costmatrix, we use a ground matrix M that accounts for the periodic boundary conditions. Wecompute the earth mover’s distance (EMD) using the function G0 = ot.emd(a,b,M) where G0 solves Kantorovich’s optimal transport problem: G0 = arg min G0 ∈ U(a,b) h G0,M i for U(a,b) = { G ∈ R N,N : G1 N =a, G T N =b } , where N is the vector of ones of length N . In particular, each entry G0 ij of the couplingmatrix G0 describes the amount of mass flowing from the mass found at x i ( T ) towards x des j .Having G0 at hand, we compute the vectors connecting the x k ( T ) with the desired po-sition t ρ des ρ T ( x k ( T )) for all k = 1 , . . . , N. To do this, we assemble position vectors ~x, ~y ∈ R N containing the first and the second component of the positions of the particles, respectively.We proceed analogously with the positions of the particles of the desired distribution. Thenwe construct the matrices X = ( ~x, ~x, . . . , ~x ) , X = ( ~x des , ~x des , . . . , ~x des ) T ,Y = ( ~y, ~y, . . . , ~y ) , Y = ( ~y des , ~y des , . . . , ~y des ) T . The distances between each particle of the current and the desired distribution, again ac-counting for the periodic boundary conditions, are contained in the matrices M X = X − X , M Y = Y − Y . For the right-hand side of the adjoints we use r x = sum( M X ∗ ( G0 > , r y = sum( M Y ∗ ( G0 > , where ∗ denotes componentwise multiplication and entries of G0 are only considered for G0 > Optimal control problem.
As mentioned before, the simulation results are verysensitive to the value of η. We therefore restrict the domain to η ∈ [ η min , η max ] . Thus, givena control u N , we update the control by v N via a projected gradient decent [21], where thestep size τ is determined via line search, i.e. we consider v N = P U ad (cid:0) u N + τ ∇ ˜ J N ( u N ) (cid:1) , where ∇ ˜ J N ( u N ) denotes the gradient of the reduced cost functional in (34) and P U ad is theprojection onto U ad . Using the solvers for the state system (28), the adjoint system (30)and steepest descent to update the control u N , we obtain Algorithm 1. Algorithm 1:
Optimal Control Algorithm
Data:
Initial data x N (0) = ( x (0) , . . . , x N (0)) with x i (0) ∈ [0 , given by (29) for i = 1 , . . . , N ; simulation time T >
0; desired values x N des ; other parametervalues; Result:
Control u N , state x N ( T ) and optimal function value J N ( x N ( T ) , u N )initialization;solve state problem (28) for x N ;solve adjoint problem (30) for ξ N ;evaluate gradient of the reduced cost functional ∇ ˜ J N ( u N ) in (34); while stopping criterion not satisfied do perform a line search to update the control u N , compute x N and evaluate J N ;solve adjoint problem (30) for ξ N ;evaluate gradient of the reduced cost functional ∇ ˜ J N ( u N ) in (34); end Numerical results
In this section we discuss numerical results, obtained with the particle algorithm intro-duced in Algorithm 1. Since we are mainly interested in the Wasserstein distance in thecost functional, we restrict ourselves to the simple case of spatially homogeneous controlparameters, i.e., u = ( θ, η ) ∈ R × R , for the numerical simulations. Note that the norms in J and J reduce to the standard norms in R . First, we set the force coefficients and the parameter values. Then we describe howartificial data is obtained for the parameter estimation, as we have no real data available.Finally, we show simulation and convergence results.6.1.
Parameter values for the results.
For the numerical examples we set N = 1200and choose the force coefficients f R and f A of the repulsion and attraction forces (4) and(5) as f R ( η | d | ) = ( αη | d | + β ) e − e R η | d | , f A ( η | d | ) = − γη | d | e − e A η | d | , resulting in the total forces F R ( d = d ( x, y ) , u ) = ηf R ( η | d | ) d, F A ( d = d ( x, y ) , u ) = ηf A ( η | d | ) R θ (cid:18) χ (cid:19) R Tθ d, with parameters α = 270 , β = 0 . , γ = 35 , e R = 100 , e A = 95 as in [5]. Moreover, we set dt = 2 and T = 10000 leading to 5000 time steps for each solve of the forward, adjointand gradient computation. The optimization iteration is stopped when the relative gradientsatisfies the condition k∇ J k k k∇ J k < (cid:15) stop , where ∇ J corresponds to the first gradient of the computation and ∇ J k denotes the gradientof the current iteration in the optimization procedure. We choose (cid:15) stop = 0 .
05 for allsimulations. The reference values in the cost function are set to θ ref = 0 . π and η ref = 1, and EAN-FIELD OPTIMAL CONTROL FOR BIOLOGICAL PATTERN FORMATION 15 the scaling factors in the cost function are λ = 1 e − and λ = 1 e − . The parameter valuefor the anisotropy is χ = 0 . . The admissible interval for η is given by [ η min , η max ] = [0 . , . . Artificial data.
As we do not have any real data available, we compute some arti-ficial data to validate our approach. We therefore choose the parameters θ data and η data ,and consider some random initial condition x (0) which consists of uniformly distributedpositions in [0 , . The initial positions for the optimization is another sample of uniformlydistributed positions in [0 , . This induces some noise and therefore, we do not expect thatthe algorithm fits the data perfectly. We use different values for the data parameters andgive details for every simulation together with the corresponding simulation results.6.3.
Simulation results.
Let θ , η denote the initialization values for the optimizationprocedure and let θ data , η data be the values for the artificial data. The simulation resultsshown below correspond to the following parameters:P1) data: θ = 0 . π, θ data = 0 . π, η = 0 . , η data = 1 . θ opt = 0 . π, η opt = 1 . θ = 0 . π, θ data = 0 . π, η = 0 . , η data = 0 . θ opt = 0 . π, η opt = 0 . θ = 0 . π, θ data = 0 . π, η = 0 . , η data = 0 . θ opt = − . π, η opt = 0 . λ and λ are chosen very small, the total cost is mainlydriven by J which corresponds to the Wasserstein distance of the discrete densities. Thisis the desired behaviour.It is interesting to see that for P
3) the optimal angle is approximating the shifted referenceangle, i.e., θ opt ≈ θ art − π . This occurs as we allow θ ∈ R and have very small regularizationparameters λ and λ . We also see the difference in the plots of the cost functional J .Indeed, note that for P
3) the value of J at the end of the optimization is about one orderof magnitude larger then the ones corresponding to P
1) and P . Conclusion
We proposed a mean-field optimal control ansatz to identify parameters underlying given,artificially generated patterns. The state system is an agent-based model with anisotropicinteraction forces that lives on the torus. The identification algorithm used gradient infor-mation that is computed with the help of the first order optimality conditions. The costfunctional penalizes the Wasserstein distance of the data pattern and the modelled pat-tern resulting from the state system for large times. Numerical results on the particle leveldemonstrate the performance of the proposed method. These results can be seen as a firststep towards the modelling of complex fingerprint patterns with specific features in futurework. acknowlegements
MB has been partially supported by the German Science Foundation (DFG) through CRCTR 154 ”Mathematical Modelling, Simulation and Optimization Using the Example of GasNetworks”. MB and LMK acknowledge support from the European Union Horizon 2020 re-search and innovation programmes under the Marie Sk lodowska-Curie grant agreement No.777826 (NoMADS). LMK acknowledges support from the European Union Horizon 2020research and innovation programmes under the Marie Sk lodowska-Curie grant agreementNo. 691070 (CHiPS), the EPSRC grant EP/L016516/1, the German National AcademicFoundation (Studienstiftung des Deutschen Volkes), the Cantab Capital Institute for theMathematics of Information and Magdalene College, Cambridge (Nevile Research Fellow-ship). CT was partly supported by the European Social Fund and by the Ministry Of Science,
Figure 1.
P1) On the top-left we see the state at T = 10000 for the initialvalues θ = 0 . π and η = 0 . . The plot on the top-right reports the statesof the artificial data at T = 10000 in black and the state corresponding tothe optimized values θ opt = 0 . π, η opt = 1 . References [1] G. Albi, Y.-P. Choi, M. Fornasier, and D. Kalise. Mean field control hierarchy.
Applied Mathematics &Optimization , 76(1):93–135, 2017.[2] L. Ambrosio, N. Gigli, and G. Savare.
Gradient Flows: In Metric Spaces and in the Space of ProbabilityMeasures . Lectures in Mathematics. ETH Z¨urich. Birkh¨auser Basel, 2005.[3] M. Bongini, M. Fornasier, F. Rossi, and F. Solombrino. Mean-field pontryagin maximum principle.
Opt.Theo. Appl. , 175(1):1–38, 2017.[4] M. Burger, L. Caffarelli, P. A. Markowich, and M.-T. Wolfram. On a boltzmann-type price forma-tion model.
Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences ,469(2157):20130126, 2013.
EAN-FIELD OPTIMAL CONTROL FOR BIOLOGICAL PATTERN FORMATION 17
Figure 2.
P2) On the top-left we see the state at T = 10000 for the initialvalues θ = 0 . π and η = 0 . . The plot on the top-right reports the statesof the artificial data at T = 10000 in black and the state corresponding tothe optimized values θ opt = 0 . π, η opt = 0 . [5] M. Burger, B. D¨uring, L. M. Kreusser, P. A. Markowich, and C.-B. Sch¨onlieb. Pattern formation of anonlocal, anisotropic interaction model. Math. Models Methods Appl. Sci. , 28(03):409–451, 2018.[6] M. Burger, M. Di Francesco, P. A. Markowich, and M.-T. Wolfram. Mean field games with nonlinearmobilities in pedestrian dynamics.
Discrete & Continuous Dynamical Systems - B , 19(5):1311–1333,2014.[7] M. Burger, R. Pinnau, C. Totzeck, and O. Tse. Mean-field optimal control and optimality conditions inthe space of probability measures. Preprint, arXiv:1902.05339, 2019.[8] M. Burger, R. Pinnau, C. Totzeck, O. Tse, and A. Roth. Instantaneous control of interacting particlesystems in the mean-field limit.
Journal of Computational Physics , 405:109181, 2020.[9] J. A. Carrillo, Y.-P. Choi, C. Totzeck, and O. Tse. An analytical framework for a consensus-based globaloptimization method.
Math. Mod. Meth. Appl. Sci. , 28(6), 2018.[10] J. A. Carrillo, B. D¨uring, L. M. Kreusser, and C.-B. Sch¨onlieb. Equilibria of an anisotropic nonlocalinteraction equation: Analysis and numerics. Preprint, arXiv:1912.09337, 2019.[11] J. A. Carrillo, B. D¨uring, L. M. Kreusser, and C.-B. Sch¨onlieb. Stability analysis of line patterns of ananisotropic interaction model.
SIAM J. Appl. Dyn. Syst. , 18(4):1798–1845, 2019.
Figure 3.
P3) On the top-left we see the state at T = 10000 for the initialvalues θ = 0 . π and η = 0 . . The plot on the top-right reports the statesof the artificial data at T = 10000 in black and the state correspondingto the optimized values θ opt = − . π, η opt = 0 . [12] P. Degond, M. Herty, and J. G. Liu. Meanfield games and model predictive control. Communicationsin Mathematical Sciences , 15:1403–1422, 2017.[13] M. R. D’Orsogna, Y. L. Chuang, A. L. Bertozzi, and L. S. Chayes. Self-propelled particles with soft-coreinteractions: Patterns, stability, and collapse.
Phys. Rev. Lett. , 96:104302, Mar 2006.[14] B. D¨uring, C. Gottschlich, S. Huckemann, L. M. Kreusser, and C.-B. Sch¨onlieb. An anisotropic inter-action model for simulating fingerprints.
J. Math. Biol. , 78(7):2171–2206, 2019.[15] B. D¨uring, P. Markowich, J.-F. Pietschmann, and M.-T. Wolfram. Boltzmann and Fokker–Planck equa-tions modelling opinion formation in the presence of strong leaders.
Proceedings of the Royal SocietyA: Mathematical, Physical and Engineering Sciences , 465(2112):3687–3708, 2009.[16] R. Flamary and N. Courty. Pot: Python optimal transport library. https://pythonot.github.io , 2017.[17] M. Fornasier and F. Solombrino. Mean-field optimal control.
ESAIM: Control, Optimisation and Cal-culus of Variations , 20(4):1123–1152, 2014.[18] A. Gerisch and M. A. J. Chaplain. Mathematical modelling of cancer cell invasion of tissue: local andnon-local models and the effect of adhesion.
J. Theoret. Biol. , 250:684–704, 2008.
EAN-FIELD OPTIMAL CONTROL FOR BIOLOGICAL PATTERN FORMATION 19 [19] F. Golse.
Macroscopic and Large Scale Phenomena: Coarse Graining, Mean Field Limits and Ergodic-ity , chapter On the Dynamics of Large Particle Systems in the Mean Field Limit, pages 1–144. SpringerInternational Publishing, Cham, 2016.[20] M. Herty, C. Kirchner, and A. Klar. Instantaneous control for traffic flow.
Mathematical Methods in theApplied Sciences , 30(2):153–169, 2007.[21] M. Hinze, R. Pinnau, M. Ulbrich, and S. Ulbrich.
Optimization with PDE Constraints . Springer, 2009.[22] L. M. Kreusser and M.-T. Wolfram. On anisotropic diffusion equations for label propagation. arXiv:2007.12516 , 2020.[23] M. K¨ucken and C. Champod. Merkel cells and the individuality of friction ridge skin.
J. Theoret. Biol. ,317:229 – 237, 2013.[24] B. Piccoli, F. Rossi, and E. Tr´elat. Control to flocking of the kinetic cucker–smale model.
SIAM Journalon Mathematical Analysis , 47(6):4685–4719, 2015.[25] R. Pinnau, C. Totzeck, O. Tse, and S. Martin. A consensus-based model for global optimization and itsmean-field limit.
Math. Mod. Meth. Appl. Sci. , 27(1), 2017.[26] J. P. Taylor-King, B. Franz, C. A. Yates, and R. Erban. Mathematical modelling of turning delays inswarm robotics.
IMA Journal of Applied Mathematics , 80(5):1454–1474, 2015.[27] C. Totzeck. An anisotropic interaction model with collision avoidance. arXiv:1912.04234 , 2019.[28] C. Totzeck and M.-T. Wolfram. Consensus-based global optimization with personal best. arXiv:2005.07084arXiv:2005.07084