A numerical investigation of Brockett’s ensemble optimal control problems

This paper is devoted to the numerical analysis of non-smooth ensemble optimal control problems governed by the Liouville (continuity) equation that have been originally proposed by R.W. Brockett with the purpose of determining an efficient and robust control strategy for dynamical systems. A numerical methodology for solving these problems is presented that is based on a non-smooth Lagrange optimization framework where the optimal controls are characterized as solutions to the related optimality systems. For this purpose, approximation and solution schemes are developed and analysed. Specifically, for the approximation of the Liouville model and its optimization adjoint, a combination of a Kurganov–Tadmor method, a Runge–Kutta scheme, and a Strang splitting method are discussed. The resulting optimality system is solved by a projected semi-smooth Krylov–Newton method. Results of numerical experiments are presented that successfully validate the proposed framework.


Introduction
The ensemble control problems considered in this paper were proposed by Brockett [7][8][9] to design efficient and robust control strategies for steering ensembles of trajectories of dynamical systems in a desired way. For this purpose, the adequate model governing the evolution of the ensembles expressed in terms of a density is the hyperbolic Liouville (continuity) equation. In application, this ensemble may represent the probability density of trajectories of multiple trials of a dynamical system with the initial conditions specified by a distribution function, or the physical density of multiple non-interacting systems (e.g., particles). In both cases, the function that determines the dynamics of these systems appears as the drift coefficient of the Liouville equation. Therefore, the Liouville framework allows to lift the problem of controlling a single trajectory of a finite-dimensional dynamical system to the optimal control problem governed by a partial differential equation (PDE) for a continuum (ensemble) of dynamical systems subject to the same control strategy. We remark that the Liouville equation represents also the fundamental building block for continuity models as the Fokker-Planck equation for stochastic systems and the Boltzmann equation for atomistic models. In particular, it can be used to model particle streaming while neglecting collisions but allowing to consider source terms [19]. Thus, one of the purposes of this work is to present a numerical optimization framework devoted to ensemble optimal control problems that can be applied to similar problems involving continuity equations; see, e.g., [2,13,18,21,31] for different classes of these equations.
The formulation and theoretical investigation of the ensemble control problems considered in this paper are presented in [3], where existence and regularity of solutions for a class of Liouville optimal control problems is discussed in the case where the drift function corresponds to a composition of linear and bilinear (input-affine) control mechanisms for the underlying dynamical system. In this class of problems, the purpose of the controls is to steer the ensemble of trajectories along a given path and to come close to a desired target configuration at a given final time. As in [7,8], these objectives are formulated in terms of minimizing different expectation functionals that include appropriately chosen costs of the controls. In particular, we consider L 2 -, L 1 -and H 1 -costs of the controls, where the L 2 term is a classical regularization term, the L 1 -cost has the purpose to promote minimum action of the control during the time evolution by promoting sparsity [15,37], and the H 1 -cost corresponds to the minimum attention control proposed by Brockett in [9].
The challenges of the numerical investigation presented in this paper are manifold. One of these challenges is that we are a considering a nonlinear control mechanism in the Liouville model where the controls multiply the density function, and this product is subject to spatial differentiation. A further challenge posed is that the numerical approximation of the Liouville equation must guarantee non-negativity of the computed density in addition to the required properties of accuracy and stability. Moreover, turning to the functional structure of the controls' objectives, we notice that ensemble cost functionals are a much less investigated topic, especially in combination with non-smooth costs of the controls. Furthermore, the presence of L 1 costs and box constraints on the values of the controls require further numerical analysis effort due to the resulting lack of Fréchet differentiability of the resulting optimization problem.
On the other hand, our formulation covers and extends Brockett's ensemble optimal control strategy so as to address many possible requirements in applications of this framework. For this purpose, in this work we consider all numerical analysis issues concerning the solution of our general Liouville-based ensemble optimal control problems.
As thoroughly discussed in [3] and illustrated below, the first step in solving our ensemble optimal control problems is the derivation of the corresponding first-order optimality conditions that consist of the controlled Liouville equation, its optimization adjoint (having the structure of a transport equation), and a variational inequality that we may also call (with some abuse of wording) the optimality condition equation. The numerical solution of this optimality system proceeds along two main steps that are the numerical approximation of the equations involved and their solution by a numerical optimization scheme.
For the approximation step, we present a novel formulation and analysis of discretization of the Liouville equation and its optimization adjoint model; the latter is called the adjoint Liouville equation. For the former, we consider the well-known second-order finite-volume Kurganov-Tadmor (KT) discretization scheme for the spatial flux derivatives that results in a generalized monotonic upwind scheme for conservation laws (MUSCL). For the temporal discretization, we use the secondorder strong stability preserving Runge-Kutta (SSPRK2) discretization scheme. Such schemes possess several important properties (such as conservation of probability) that are inherent to the exact solutions of the Liouville equation. In addition, because the solution of the Liouville equation represents a density function, it is crucial that the numerical solution remains non-negative over all times.
We prove that our SSPRK2-KT scheme preserves positivity subject to a restriction on the time-step size. Further, we prove that our scheme is second-order convergent in the L 1 norm. This result is less-known in the context of generic finite-volume schemes for linear conservation laws. For the adjoint Liouville equation, which is a transport equation with a source term, we use a second-order Strang time-splitting scheme combined with the KT spatial discretization scheme, and for the resulting approximation we prove second-order accuracy. Further, we notice that the optimality condition equation is a variational inequality involving an integral for which we use second-order accurate quadratures, and we implement a projection step in the optimization procedure. Notice that, while second-order accuracy for the above three components of the optimality system is separately guaranteed by suitable approximation, we are not able to prove this order of accuracy of the coupled system; this is an issue that remains widely open in the scientific literature, apart of the case of much simpler problems with linear control mechanisms; see, e.g., [6].
The second fundamental step in solving our ensemble optimal control problems is the design of a numerical optimization procedure. For this purpose, one recognizes that the optimality condition equation provides the semi-smooth gradient of the ensemblecost functional along the constraint given by the Liouville model. However, because of the presence of control constraints and the combination of L 2 -, L 1 -and H 1 -costs, the assembling of our gradient is challenging. In particular, by imposing constraints on the value of the control, we are required to implement a H 1 projection of the control update. At this point, we remark that the combination of L 1 -and H 1 -costs and the H 1 projection are less investigated in the literature.
However, this effort is very well justified by our purpose of implementing a state-ofthe-art semi-smooth Krylov-Newton methodology for our new class of PDE optimal control problems. In doing this, we also rely on past experience in [14,15], and the resulting Newton scheme is used to validate our optimal control framework.
In the next Section, we illustrate the formulation of Liouville ensemble optimal control problems, and discuss the chosen control mechanism and the constitutive terms of an ensemble cost functional. Further, in correspondence to our optimization setting, we present the optimality system and discuss the construction of the gradient.
In Sect. 3, we investigate the approximation of the Liouville equation and of adjoint Liouville equation in the optimality system. For the former, we consider a combination of a second order accurate strong stability preserving Runge-Kutta discretization in time and the Kurganov-Tadmor finite volume discretization in space. For the latter, we discuss a scheme that combines the Kurganov-Tadmor discretization and a Strang's splitting technique. For both methods, we present a detailed analysis of stability and accuracy, and in the case of the Liouville equation we prove that our scheme is positive preserving.
Section 4 is devoted to the implementation of our semi-smooth Krylov-Newton method that requires the numerical solution of the Liouville equation and its adjoint and the implementation of the gradient together with a H 1 -projection procedure for the controls.
In Sect. 5, we present results of numerical experiments with our solution methodology that validate our optimal control framework in terms of the ability of the controls to perform the given tasks. For this purpose, we consider the tracking of non-differentiable trajectories and also the case of bimodal distributions. A section of conclusion completes this work.

Notation
In this section, we present our notation that we use throughout the paper.
Given a bounded domain ⊂ R d , the symbol C ∞ c ( ) denotes the space of infinitely often differentiable functions with compact support in . Given k ∈ N, we denote by C k ( ) the space of all k-times continuously differentiable functions defined on , and by C k b ( ) the subspace of C k ( ) formed by functions which are uniformly bounded together with all their derivatives up to the order k. We equip For α ∈ ]0, 1], we denote with C 0,α ( ) the classical Hölder space (Lipschitz space if α = 1), endowed with the norm In particular, For k ∈ N and 1 ≤ p ≤ +∞, we denote with W k, p ( ) the usual Sobolev space of L p functions with all the derivatives up to the order k in L p ; we also set H k ( ) := W k,2 ( ). For 1 ≤ p < +∞, let W −k, p ( ) denote the dual space of W k, p ( ). For any p ∈ [1, +∞], the space L p loc ( ) is the set formed by all functions which belong to L p ( 0 ), for any compact subset 0 of . Furthermore, we make use of the so-called Bochner spaces. Given a Banach space (X , · X ) and a fixed time T > 0, we define for 1 ≤ p < ∞, and a generic representative function φ = φ(x, t), the spaces Given a Banach space X and a sequence n n , we use the notation n n ⊂ X meaning that n ∈ X for all n ∈ N and that this sequence is uniformly bounded in X : there exists some constant M > 0 such that n X ≤ M ∀n ∈ N. Given two Banach spaces (X , · X ) and (Y , · Y ), the space X ∩ Y , endowed with the norm · X ∩Y := · X + · Y , is still a Banach space.
For every p ∈ [1, +∞], we use the notation In addition, given two vectors u and v in R d , we write u ≤ v if the inequality is satisfied component by component by the two vectors: namely, u i ≤ v i for all 1 ≤ i ≤ d.

Formulation of ensemble optimal control problems
Consider a particle whose position at time t is denoted with ξ(t) ∈ R d . Suppose that this particle is subject to a velocity field a(x, t) over R d , where (x, t) ∈ R d × [0, T ], for some final time T > 0; then the particle's trajectory is obtained by integratinġ ξ(t) = a ξ(t), t assuming an initial condition ξ(0) = ξ 0 . Now, suppose we have an infinite number of non-interacting particles subject to the same vector field and being distributed with a smooth initial density ρ |t=0 = ρ 0 ; then the evolution of this material density is modelled by the following Liouville equation with the initial condition at t = 0 given by ρ(x, 0) = ρ 0 (x). Notice that, in this model, the state variable x of the dynamical system defined by a, becomes the space variable in the Liouville equation. We call a the drift function.
We have the same model (2.1) if we consider a unique particle subject to the flow a and having the initial condition ξ 0 chosen based on the probability density ρ 0 . In this case, the Liouville equation governs the evolution of the probability density function ρ of the position of the particle in the interval [0, T ]. However, the significance of the Liouville equation above is not limited to the case where x denotes the space coordinate and a a velocity field. In fact, it applies equally well in the case where x represents the velocity of the particle and a plays the role of acceleration/force. Another possibility is to identify x with the position and velocity of a particle in the phase space, and in this case the Liouville operator corresponds to the streaming part in the transport and Boltzmann equations [10,19].
Clearly, the interpretation of ρ as a probability or material density leads to the requirement that the initial condition for the Liouville model is non negative, ρ 0 ≥ 0. Moreover, we can normalize the total probability or mass requiring that With this conditions, one can show that the evolution of ρ modeled by the Liouville equation (2.1) has the following properties (2. 2) The first property can be proved by the vanishing viscosity method and the maximum principle or solving along characteristics; see, e.g., [20,22]; the second property follows from a simple application of the divergence theorem. We remark that the Liouville equation allows to model the transport of the (material or probability) density also in the case when the drift function is non smooth [2,18,31], and also in the case when it includes a control mechanism.
Therefore, the representation of the ensemble of trajectories in terms of an evolving density and the fact that we can manipulate the drift with a control function to achieve some purposes of the motion of the particles leads to the formulation of the following ensemble optimal control problem min u∈U ad In this problem, the drift a(x, t; u) ∈ R d includes a time-dependent vector-valued control function u, and the purpose of this control is to drive ρ such that the cost functional J is minimized.
As in many application, we consider dynamical systems that are controlled by linear and bilinear control mechanisms as follows where a 0 is a given smooth vector field, which is Lipschitz in x, and b, c ∈ R are constants. Moreover, the function u 1 corresponds to a linear control mechanism, and u 2 represents a bilinear (input-affine) control term. We assume that both functions are defined on the time interval denote the Hadamard product of two vectors; we also use the notation u = (u 1 , u 2 ). Notice that, in the simple case where a 0 = 0, b = c = 1, and ρ 0 is a normalized Gaussian unimodal distribution, then, the Liouville dynamics can be completely described by the first-and second-moment equations that include the controls u 1 and u 2 . To illustrate this fact, consider the following average operator applied to an integrable function g In particular, we have the mean Then, by taking the average of our controlled dynamical system (that is, using the Liouville equation), we obtaiṅ where the control u 1 appears as the main driving force of the mean value of the density, and u 2 determines the evolution of the variance of the density. See [7] for more details on this construction. However, the validity of this setting is very limited by the assumptions above whereas the Liouville framework allows to accommodate a more general drift and to consider multi-modal distributions of the density function. Next, we discuss the meaning of the different terms in our cost functional. The first term in (2.3) has the purpose to model the tracking of a desired trajectory ξ D (t), in the sense that minimizing this term corresponds to having all trajectories of the ensemble of particles being close to ξ D . For this purpose, we choose a function θ(x, t) that, for a fixed t, is required to monotonically increase as a function of the distance |x − ξ D (t)|. Therefore, by minimization of the first term we have that ρ is mainly concentrated on the minimum of θ corresponding to ξ D . We refer to θ as an attracting potential, and focus on the following choice Similarly, the purpose of the second term in J is to model the requirement that the density ρ at final time concentrates on a final position denoted with ξ T ∈ R d . Therefore, we may choose Notice that, because of (2.2), we can augment this functions by adding a constant such that J is bounded from below by zero. However, this shift would not influence the result of the optimization problem. The last term in (2.3) represents the cost of the controls' action. Moreover, it determines the functional space where the control is sought. In our case, the cost function κ is chosen based on the requirement of implementing a slowly varying control function with sparsity along the time horizon. A control that slowly changes in time can be obtained considering the following cost of the control where ν > 0 is a positive weight, and | · | denotes the Euclidean norm. In fact, as ν is taken larger, the emerging optimal control will result in controls with small values of its time derivative, which corresponds to the "minimum attention control" in [9]. In addition, we add a L 2 -cost of the control that should measure the total effort made by u as follows where γ > 0 is a positive weight. Clearly, if ν = γ > 0, then these two terms together correspond to a H 1 (0, T )-cost of the control. As for inducing sparsity of the control function u, we consider the following L 1 -cost of the control where δ ≥ 0. This cost promotes sparsity of the control function, in the sense that, as δ is increased, the u resulting from the minimisation procedure will be zero on larger open intervals in (0, T ); see Figure 2 of [15]. In our framework, the purpose of this control is to promote "minimum action".
Summarizing, we specify the term κ u(t) in (2.3) as follows where γ, ν > 0 and δ ≥ 0. It is clear that, with this setting, the control space U corresponds to a weighted H 1 space given by H 1 T space, endowed with the following weighted H 1 -product The notation = d/dt stands for the weak time derivative. In order to complete the modelling of the control space, we also require that each component of the control function u may only take values in a compact convex set of R. Thus, we define the following set of admissible controls where the inequalities are meant componentwise, and we choose u a = u a 1 , u a 2 and As discussed in detail in [3], the optimal control problem (2.3)-(2.4), with (2.5) and (2.7) and the admissible set of controls U ad , admits a solution. Furthermore, this solution can be characterized as the solution to the so-called first-order optimality condition. This condition can be formulated by first introducing the Liouville controlto-state map G defined by where ρ is the unique solution to the Liouville equation with the chosen control u ∈ U ad and the given initial data ρ 0 ∈ L 2 (R d ). Notice that, in our setting, with ρ 0 ∈ H m (R d ), the Liouville problem admits a unique weak solution ρ ∈ C([0, T ]; H m (R d )); see [3] for all details.
With this map, our optimal control problem can be equivalently formulated as where J (u) := J G(u), u represents the so-called reduced cost functional. Now, if δ = 0, then J (u) is Fréchet differentiable in suitable topologies [3], and a solution u to our optimal control problem necessarily satisfies the first-order optimality condition given by (∇ J (u), v − u) U ≥ 0 for all v ∈ U ad , where ∇ represents the gradient in the control space; see [27,40]. However, since we have δ > 0, our reduced cost functional is only sub-differentiable [17,41], and the formulation of the optimality condition becomes more involved, especially in our case where a H 1 In the following, we give a detailed discussion of the formulation of the optimality condition in terms of an optimality system, where the calculation of ∇ J (u) is achieved by introducing an auxiliary adjoint variable with its own constitutive equation. Further, we use the notion of Clarke's subdifferential that we illustrate in detail in Sect. 4.
The derivation of the optimality system can be conveniently done in the Lagrange framework starting from following Lagrange functional where q and q 0 represents Lagrange multipliers (the adjoint variables). In this framework, the optimality system is obtained by requiring that the Fréchet derivatives of L(ρ, u, q) with respect to each of its arguments are zero. As shown in [3], Fréchet differentiability of L with respect to ρ and q involves a loss of derivatives, which is natural owing to the hyperbolic nature of the underlying PDEs and which requires more regular initial data. Clearly, the derivatives of L with respect to q and q 0 give the Liouville equation and its initial conditions. On the other hand, the Fréchet derivative of L with respect to ρ leads to the following adjoint Liouville equation Notice that this is a transport problem with a given terminal condition and thus evolving backwards in time. However, with our choice of θ and ϕ as Gaussian functions, the adjoint equation admits a solution in the same functional space of the solution of the Liouville equation [3]. We remark that the density ρ does not explicitly appear in this equation, which means that the Liouville equation and its adjoint are decoupled; however, both are driven by the same drift. Next, we discuss the derivative of L with respect to u. For this purpose, we remark that the case δ, γ > 0 and ν = 0 is well-known in the literature [37], whereas the case that includes the minimum attention control term (ν > 0) is novel. For this reason, to ease our discussion and make it more accessible, we first discuss the case ν = 0, for which in U ad given by (2.8), and thereafter we also include the term with ν > 0.
In the case ν = 0, as proved in [3], there exists a λ ∈ ∂ g(u) := δ ∂ u L 1 , the Clarke's subdifferential [16] of the L 1 -cost, such that the following inequality condition must be satisfied at optimality Notice that in this inequality condition, the argument on the left-hand side in the scalar product represents the reduced gradient in L 2 T .
Moreover, because of the presence of constraints for the controls, there exist λ + and λ − , belonging to L ∞ for all m = 1, 2 and all 1 ≤ r ≤ d.
In (2.13), one usually refers to the first equation as the optimality condition equation; the conditions given in the second and third line are the complementarity conditions for the inequality constraints in U ad . Moreover, the last three lines give an equivalent expression for λ ∈ ∂ g(u); see [37]. In our case, λ r m can be understood to be δ sgn(u r m ), where sgn(x) is the sign function. Now, in the case ν > 0, the control is sought in U ad given by (2.8), which requires to construct the reduced gradient in the space Then, assuming that u ∈ U ad ∩ H 1 0 [0, T ]; R 2d , we compute μ by solving the following boundary-value problem which has to be understood in a weak sense. Notice that, in this construction, we have made the modelling choice that the control function u is zero at the beginning and at the end of the time interval. This setting corresponds to having the control switched on at t = 0 and off at t = T . Based on (2.12) and the definition of μ in the H 1 space, we identify the reduced gradient in this space as follows where m = 1, 2 and r = 1 . . . d. Thus, the optimality condition (2.12) becomes The result of this section is that a solution to our ensemble optimal control problems can be characterized as the solution to the following optimality system. (2.20) Notice that, for the sake of better readability, in the following we choose ν = γ > 0 and δ > 0.

Approximation of the Liouville optimality system
In this section, we discuss the spatial and temporal discretization of the Liouville equation and its adjoint in the optimality system. Our aim is to develop an approximation framework that is second-order accurate and preserves the two essential properties of the continuous Liouville model given in (2.2), namely positivity and conservativeness of its solution.
For simplicity of notation, in the following we focus on a two-dimensional problem, i.e. d = 2. Then a = (a 1 , a 2 ) ∈ R 2 . In view of applications to the numerical study of our optimal control problem, we consider a large but bounded convex domain ⊂ R 2 : We also fix a smooth initial density ρ 0 that is (by machine precision) compactly supported in . For θ and ϕ we take Gaussian functions having sufficiently small variance and centred sufficiently far from the boundary of , so that (by machine precision) we can assume that also those functions are compactly supported in . Then, we solve problems (2.17) and (2.18) in × [0, T ], supplemented with homogeneous Dirichlet boundary conditions on ∂ . Notice that, in this setting, it is possible to use the results of [3] to prove existence and uniqueness of smooth enough solutions to (2.17) and (2.18). For this purpose, one extends the functions ρ 0 , θ and ϕ to be zero outside the domain , and the drift function a to a smooth function, which is bounded on R d together with all its space derivatives. We remark that this is always possible, for instance by multiplying a with a smooth compactly supported function χ of the space variable only, such that χ ≡ 1 on a neighbourhood of .
We consider our solutions on a time interval [0, T ], where T > 0 is chosen such that the corresponding solutions ρ to (2.17) and q to (2.18) are still compactly supported in , far away from its boundary ∂ . Observe that this property is true by finite propagation speed, since the (extended) drift is bounded on R d . For = (−B, B) × (−B, B) fixed above, we set a numerical grid that provides a partitioning of in N x × N x , N x > 1, equally-spaced non-overlapping square cells of side length h = 2B/N x . On this partitioning, we consider a cell-centred finite-volume setting, where the nodal points at which the density and adjoint variables are defined are placed at the centres of the square volumes. These nodal points are given by Therefore, the elementary cell is defined as Thus, the computational domain is given by Analogously, the time interval [0, T ] is divided in N t > 1 subintervals of length t and the points t k are given by This defines the time mesh In this setting, the cell average of the density ρ (and so of any integrable function), on the cell with centre (x i 1 , x j 2 ) at time t k , is given bȳ Notice that, in our numerical setting, function values are identified with their cellbased average located at the cell centres. For this reason, our numerical framework aims at determining approximation of theses averages. Specifically, we discuss a discretization scheme that results in values ρ k i, j that approximateρ k i, j . Similarly, we denote with q k i, j the numerical approximation ofq k i, j that is computed as in (3.1). We also consider a piecewise constant approximation to the time-dependent control functions, where we denote with u k+1/2 the value of the control in the time interval [t k , t k+1 ). Further, for the projection of a continuous u to the corresponding approximation space, we set u k+1/2 = u(t k ). For a function f defined on Q h, t , we define the discrete norms · 1,h and · ∞,h as follows:

A Runge-Kutta Kurganov-Tadmor scheme for the Liouville equation
In this section, we discuss a suitable approximation of our controlled Liouville equation in ×[0, T ]. Supposing that ρ 0 has compact support, and because of finite propagation speed, we can choose such that the solution ρ at the boundary ∂ is zero for all times t ∈ [0, T ]. For our purpose, we focus on the finite-volume scheme proposed by Kurganov-Tadmor (KT) in [26] that involves a generalized MUSCL flux. To describe this scheme, we denote the flux in the Liouville equation as a function of ρ with f (ρ) = aρ = a(x, t) ρ(x, t). With this definition, the KT scheme for the Liouville equation in semidiscretized form is given by where f = ( f 1 , f 2 ) = (a 1 ρ, a 2 ρ), and similarly for H x 2 i, j±1/2 (ρ + , ρ − ; t). In this formula, the so-called local speeds V x r (t) are given by , the approximation of ρ at the cell edges is given by intermediate values that approximate the function value from above respectively from below as follows (3.5) We approximate the partial derivatives of ρ using the minmod function as follows: In direction x 1 , this approximation is given by An analogous expression holds in the direction x 2 . The multivariable minmod function for vectors x ∈ R d is given by otherwise.
Next, we discuss the local truncation error of the semi-discrete KT scheme (3.2)-(3.4). Proof The flux H , given in (3.3), is a first-order Rusanov flux [34] that is C 2 with Lipschitz continuous partial derivatives with respect to ρ + , ρ − , in a neighbourhood of ρ i, j . Further, by using a Taylor series expansion, we have the following approximation and a similar result holds in the x 2 direction, except at the points of extrema, which are characterized by (ρ x 1 ) i, j = 0 (see [30,Th. 3.2]). Using the result in [29, Lemma 2.1], we have that the semi-discrete scheme (3.2)-(3.5) is second-order accurate in space except possibly at points of extrema.
For the time discretization of the Liouville equation (2.4), we use a second-order strong stability preserving Runge-Kutta (SSPRK2) method [35] (also known as the Huen's method). The combination of this scheme with the KT discretization of the flux f , given in (3.2), results in a new approximation method that we call the SSPRK2-KT-scheme. This scheme is implemented by the Algorithm 3.1 given below using the following definition where H denotes H x j ·,· , j = 1, 2, given in (3.3) corresponding to the time step t k .

Algorithm 3.1 SSPRK2-KT scheme
Require: ρ 0 i, j , F Ensure: Solve the Liouville equation in ρ as follows 1: Set k = 0 2: while 0 ≤ k < N t do 3: is computed using (3.7). 5: In (t k , t k+1 ), compute ρ (2) i, j = ρ (1) i, j + t F(ρ (1) i, j ) with initial condition ρ (1) i, j , where F(ρ (1) i, j ) is computed using (3.7). 6: Time step update: i, j . 7: end for 8: k = k + 1 9: end while 10: return ρ k Now, we study the properties of the SSPRK2-KT scheme given in Algorithm 3.1. This method has the following strong stability property [23, Lemma 2.1] Further, we have conservativeness of the total probability (or mass) as a consequence of the finite-volume formulation: From the first step of the SSPRK2-KT scheme, given in Algorithm 3.1, summing up over all indices i, j ∈ {1, . . . , N x } and using the fact that the solution has zero flux on the boundary (since it has compact support in ), we get In a similar way, we have Thus, Iterating over k, we have Next, we show that, starting from a non-negative initial density, the solution obtained with the SSPRK2-KT scheme remains non-negative. For this purpose, we define the CFL-number and require the conditions on the components of the drift a given in (2.5) given by Notice that the control u belongs to the set U ad defined in (2.8). Then, for j = 1, 2, we see that so that the aforementioned conditions on the components of the drift a are satisfied under the following CFL condition.
This CFL condition only depends on the components of the vector a 0 , and it does not rest on the unknowns of the problem. With the conditions (3.9), we can prove the following lemma on the positivity of the SSPRK2-KT scheme.

10)
Proof The SSPRK2-KT scheme, given in Algorithm 3.1, comprises of a two-step Euler scheme that results in the computation of ρ (1) and ρ (2) and a final averaging step. To prove positivity of the SSPRK2-KT scheme, it is enough to show that the solutions obtained in each of the two Euler steps are positive. Without loss of generality, we prove that the solution obtained in the first step of the SSPRK2-KT scheme is positive. A similar analysis holds true for the second step. Let ρ k i, j ≥ 0 for fixed 0 ≤ k < N t . We will show that ρ k+1 i, j ≥ 0 for all i, j = 1, . . . , N x . For this purpose, notice that the SSPRK2-KT scheme can be written as follows where all discrete quantities on the right are considered at the timestep t k . We see that the first four terms on the right hand side in (3.11) are always non-negative, provided that ρ ± i±1/2, j , ρ ± i, j±1/2 ≥ 0. The remaining terms are non-negative under the CFL-condition (3.9).
Thus, it remains to show that ρ ± i+1/2, j , ρ ± i, j+1/2 ≥ 0 for all i, j = 1, . . . , N x , where ρ ± i, j is given as in (3.5). We consider each of the expressions for (ρ x 1 ) k i, j in the direction of x 1 given as in (3.6). First, assume that ρ x 1 , which is one of the possible values of the minmod limiter in (3.6).
Then it follows that This is non-negative, since ρ k i, j ≥ 0 for all i, j = 1, . . . , N x by the inductive assump- The other cases for the value of (ρ x ) k i, j = 0 follow analogously. If (ρ x ) k i, j = 0, then ρ ± i+1/2, j = ρ i+1, j ≥ 0 and ρ ± i, j+1/2 = ρ i, j+1 ≥ 0. This proves the lemma. Remark 3. 1 The proof of the above lemma follows similar arguments as in [25, Theorem 2.1]. However, a primary difference is that in [25], the positivity result is proved using one-sided local speeds, exploiting the structure of the hyperbolic equation, whereas in our case, the proof relies on conversion of the intermediate values ρ ± to the cell-average values ρ k i, j and then showing that ρ k i, j ≥ 0 implies ρ k+1 i, j ≥ 0, which seems a much simpler approach.

Remark 3.2
Under the same CFL-like condition (3.9), the proof of Lemma 3.3 can be extended to the case of a SSPRK-KT scheme with a Runge-Kutta method of pth order, p ∈ N, that is given as an average of p Euler steps.

Remark 3.3
For the case where (ρ x 1 ) k i, j = 0, the approximations of ρ at the cell-edges, given by (3.5), are piecewise constant in the cell ω i j h . Thus, the KT scheme given by (3.2)-(3.5), reduces to a linear upwind scheme that is locally first-order accurate, TVD and positive. This is consistent with the Godunov's barrier theorem.
Next, we prove discrete L 1 stability of the SSPRK2-KT scheme without a righthand side.

Lemma 3.4 (Stability) The solution ρ k i, j obtained with the SSPRK2-KT-scheme in
under the CFL condition (3.9).

Proof Using Lemma 3.2, we have
Again, from Lemma 3.3, we have that the solution obtained with the SSPRK2-KT scheme is positive under the CFL condition (3.9). This gives us the following relation which proves the result.
Next, we aim at proving the L 1 convergence of the SSPRK2-KT scheme. For this purpose, we prove the following stability estimate for the discrete solution of the Liouville equation (2.4) with a right-hand side function g(x, t). We remark that in the case g ≡ 0 the solution may not be positive.

Lemma 3.5 Let ρ k i, j be the SSPRK2-KT solution to the Liouville equation (2.4) with a Lipschitz continuous right-hand side g(x, t) and let the CFL condition (3.9) be fulfilled. Then this solution satisfies the following stability estimate
Proof The SSPRK2-KT scheme, given in Algorithm 3.1, for the Liouville equation (2.4) with a right-hand side g(x, t) can be written as where with n = (1) and n = k correspond to the solution ρ (1) and ρ k , respectively, at the time step t k and analogously for ρ n± . Moreover, also the drift is always considered at the time-step t k . The equations in (3.12) can be rewritten in a compact form with a suitable function H as follows Now, the KT flux H , given in (3.3), is a combination of the monotonicity preserving Rusanov flux and the monotonicity preserving MUSCL reconstruction. This leads to the SSPRK2-KT scheme to be monotone preserving under the CFL condition [24]. Thus, H is a monotone non-decreasing function of its arguments. Then the following discrete entropy inequality holds for the specific Kruzkov entropy pair (|ρ|, sgn(ρ)) (see [39,Lemma 2.4]) where 1,k ·,· , 2,k ·,· are the conservative entropy fluxes defined as follows for i, j = 1, . . . , N x Therefore, we have for k = 0, . . . , N t − 1 Summing up over all i, j and because of our assumption on ρ being zero on the boundary, we have which iteratively gives us Next, we consider the local consistency error of our SSPRK2-KT at the point The accuracy result for the KT scheme, given by Lemma 3.1, the MUSCL reconstruction error given in Equation (60)  Define the error at the point ( Notice that ρ k i, j satisfies (3.12) with g ≡ 0 by construction. Further, the exact solution ρ satisfies (3.12) with the consistency error T k i, j as source term. Hence, by taking the difference, the error satisfies (3.12) with the source term given by −T k i, j . From Lemma 3.5, we obtain e k+1 This leads to the following result on the L 1 convergence of the solution obtained using the SSPRK2-KT scheme.

Numerical analysis of the Strang splitting scheme
In this section, we deal with the numerical solution of the adjoint equation (2.11). In this case, we have a terminal condition, and the adjoint problem requires evolution backward in time. For this reason, it is convenient to perform a change of the time variable as follows: With this transformation, we can rewrite (2.11) in the following way To solve this problem, we apply the Strang splitting method [38] by first rewriting (3.15) as follows with q(x, 0) = −ϕ(x), (3.16) which is defined in ×[0, T ]. Furthermore, we assume that ϕ and θ have (by machine precision) compact support for all times inside the interval [0, T ]. See the numerical experiments section for the specific choices for θ and ϕ. Then, we solve problem (3.16), supplemented with homogeneous Dirichlet boundary conditions. We can conveniently illustrate the Strang splitting method applied to (3.16) remaining at the continuous level within one time interval. Let us consider the solution of the adjoint equation (3.16) at time τ k given by q k (x), x ∈ . Then, the first step of our solution scheme is to solve the following equation (3.17) For this purpose, we use the SSPRK2-KT scheme given in Algorithm 3.1. We denote the solution to this problem with q 1 .
In the second step, for each x fixed, we analytically solve the following linear ordinary differential equation (3.18) Let the solution obtained in this step be denoted with q 2 . The analytical solution of (3.18) at the discrete level is given in (3.22) below. Notice that q 1 (x, t) and q 2 (x, t) denote functions of continuous variables. We define q k 1,i, j and q k 2,i, j as their discrete counterparts using the finite-volume approximation strategy.
The last step is to solve (3.17) with the SSPRK2-KT scheme with the initial condition q 2 (·, τ k+1/2 ) in the time interval [τ k+1/2 , τ k+1 ]. This problem is formulated as follows In a numerical setting, the solution obtained in this step is the desired solution of the adjoint equation (3.16), and q k+1 denotes the adjoint variable at time τ k+1 . Notice that, by our numerical approximation strategy for u, the value of u in [τ k , τ k+1 ) is constant. The steps of the Strang splitting scheme are outlined in Algorithm 3.2 below.

Algorithm 3.2 Kurganov-Tadmor-Strang (KTS) scheme
Require: q 0 = −ϕ, F Ensure: Solve adjoint equation in q 1: k = 0 2: while 0 ≤ k < N t do 3: Apply one temporal step of Algorithm 3.1 in (τ k , τ k+1/2 ), with inputs q k i, j , −F, within the time-interval (τ k , τ k+1/2 ). Denote the solution q k+1/2 1,i, j . 5: In (τ k , τ k+1 ), solve (3.18) using exact integration as given in (3.22) 6: Apply one temporal step of Algorithm 3.1 in (τ k+1/2 , τ k+1 ), with inputs q k+1 2,i, j , −F, within the time-interval (τ k+1/2 , τ k+1 ). Denote the solution with q k+1 i, j . 7: end for 8: k = k + 1 9: end while 10: return q k Now, we discuss some properties of the Strang-splitting scheme described in Algorithm 3.2. For this purpose, we denote with q k i, j the numerical solution of (3.16) with the generic right-hand side G, at the grid point (x 1 i , x 2 j , τ k ). We have the following discrete L 1 stability estimate. be the numerical solution obtained from (3.17). Since (3.17) is solved using the SSPRK2-KT scheme, using the entropy inequality computations as in Lemma 3.5, we have Next, denoting the numerical solution as q k+1 2,i, j obtained as the analytical solution from (3.18), using an integrating factor approach in [τ k , τ k+1 ], we have (3.22) where R = div a dτ . Notice that, because the drift and G are given explicitly, and the control u is constant in the sub-interval of integration, this equation can be solved exactly. We exemplify this calculation considering G being constant in [τ k , τ k+1 ) and equal to its value at τ k . In this case, it holds that R(τ k ) = (u 1 2 ) k+1/2 + (u 2 2 ) k+1/2 τ k . In general, without any assumptions on the approximation strategy for u and G, but considering the bilinear structure of our drift function (c.f. (2.5)) and the assumption on a 0 , we can state that there exists an L > 0 such that |div (a(x, τ, u) Thus, we have Hence, by integration we obtain Further, by using (3.21), we have q k+1 (3.19) is solved using the SSPRK2-KT scheme, we have q k+1 which gives the desired result.
Next, we consider the local truncation error of our KTS scheme at the point where H is the SSPRK2-KT operator given in (3.13) and is the exact integration operator for (3.18) at time τ k , defined in (3.22  and with a smooth initial conditionq (·, τ k ). This solution can be written as follows , and S 2 = S 2 ( t) the same operator for [τ k+1/2 , τ k ]. Then the following error estimate holds With this result and the truncation error estimate of the SSPRK2-KT scheme given in 3.6, we have the following result Define the error at the point ( Then e k i, j satisfies (3.16), with the right-hand side being ; see the explanation after Lemma 3.6. Thus, from Lemma 3.7 we obtain e k+1 where M = max k∈{0,...,N t } h 2 i, j |Z k i, j |. This leads to the following result on the L 1 convergence of the solution obtained using the KTS scheme Theorem 3.2 Let q ∈ C 3 be the exact solution of the adjoint equation (3.16), with countably many extrema, and let q 0 ·,· + ϕ(·, ·) 1,h = O(h 2 ). Under the CFL condition (3.9), the solution q k i, j obtained with the KTS scheme, given by Algorithm 3.2, is second-order accurate in the discrete L 1 -norm as follows q k ·,· − q(·, ·, t k ) Remark 3. 4 We remark that results similar to Theorem 3.2 have been obtained in [11,12]. However, in these papers the equation that has been considered is the convection diffusion equation, which is parabolic, whereas we have a hyperbolic transport (adjoint) equation with a source term. Furthermore, we employ a different analysis using an entropy inequality technique for proving the discrete stability estimate that is subsequently used for proving the convergence error estimate.

A projected semi-smooth Krylov-Newton method
In this section, we illustrate a semi-smooth Krylov-Newton (SSKN) method for solving the ensemble optimal control problem (2.3)-(2.4) with the drift given by (2.5) and the cost functional setting specified in Sect. 2. We remark that our SSKN scheme belongs to the class of projected semi-smooth Newton schemes discussed in [41].
In general, a Newton method is an iterative procedure aiming at finding roots of a given function. Its peculiarity is that it may generate a sequence that can converge superlinearly or even quadratically to the solution sought.
In order to explain the Newton method in simple terms, consider the problem to find a root ζ ∈ R N of a map M : R N → R N , as follows where, for the moment, we assume that ζ → M(ζ ) is continuously differentiable. Now, denote with J (ζ ) the Jacobian of M at ζ . The Newton method generates a sequence (ζ ) ∈N by means of the following two steps The steps s 1 -s 2 are performed for = 0, 1, 2, . . ., starting with a given initial guess ζ 0 . Clearly, the Newton sequence is well defined if the Jacobian is invertible at each iterate, and we assume that this is the case in a neighbourhood N of the solution ζ , where also the inverse is uniformly bounded. With these assumptions and requiring that the initial guess ζ 0 ∈ N is sufficiently close to ζ , one can prove that the sequence (ζ ) converges quadratically to the root ζ , that is, ζ +1 − ζ 2 ≤ c ζ − ζ 2 2 , for some constant c > 0, and · 2 denotes the Euclidean norm of a vector in R N . However, in the case where M is only differentiable and provided that the following holds then the Newton sequence converges at least superlinearly, i.e. faster than linearly. The same Newton procedure (4.2) can be applied to find an extremal of the minimization problem min ζ ∈R N f (ζ ), by considering M(ζ ) := ∇ f (ζ ), where ∇ denotes the gradient in R N , and assuming that f : R N → R is twice differentiable. In this case, we have that J (ζ ) = ∇ 2 f (ζ ). Now, in the case of a constrained optimization problem min ζ ∈K f (ζ ), where K ⊂ R N is closed and convex, an extremal ζ of this problem is characterized by the inequality ∇ f (ζ ) · (ζ − ζ ) ≥ 0. However, this inequality can be equivalently written as follows where P K is the projection of R N onto K , and s > 0 is arbitrary but fixed. Therefore the solution of the optimality condition in the form of an inequality can be reformulated as a root problem. However, even if f is continuously differentiable, the function F is not. On the other hand, if ∇ f is locally Lipschitz, then also F is locally Lipschitz continuous. We see that a lack of differentiability of ∇ f (ζ ) or the presence of constraints as above hinder the application of the Newton scheme to solve optimization problems. This situation has motivated a great effort towards the generalization of the notion of differentiability that makes possible to pursue the Newton approach also in nondifferentiable cases; see [16,17,41] for details and further references.
The main assumption for this generalization is the Lipschitz continuity of the map ζ → F(ζ ), in which case Rademacher's theorem [41] states that this map is almost everywhere differentiable. Based on this result, the notion of differentiability has been extended as follows; see [41] for a detailed discussion. Definition 4.1 Assuming F : R N → R N be locally Lipschitz continuous, we have the following generalized Jacobians of F at ζ : (a) the Bouligand subdifferential given by where U nd is the set of points where F fails to be Fréchet differentiable and J (ζ ) denotes the Jacobian of F at ζ ; With this construction, we can apply (4.2) by choosing a generalized Jacobian J ∈ ∂F(ζ ). However, in order to guarantee superlinear convergence of the resulting Newton sequence, the following property of semi-smoothness is required; see [32,41]. Notice that the discussion above has focused on finite dimensional spaces, which is also the case of our numerical optimization problem. However, the subdifferential framework given above has been extended also to maps acting between infinitedimensional Banach spaces [16,17,41]. In particular, we can apply it to our ensemble optimal control problem (2.9), that is, min u∈U ad J (u) := J G(u), u . One can recognize that J (u) is not Fréchet differentiable (see also details in [3]) due to the presence of the L 1 -cost, which however is Lipschitz in u. In fact, in Sect. 2, we have used sub-differential calculus [41] to determine the gradient ∇ J (u), and formulated the first-order optimality condition ( ∇ J (u), v − u) U ≥ 0 for all v ∈ U ad . Therefore we can proceed as in (4.4) and consider to apply the Newton scheme with generalized Jacobian to the equation However, although this procedure is standard with control problems with L 2 -L 1 costs [14,15,41], it becomes very cumbersome in our case with H 1 costs. For this reason, we consider a projected semi-smooth Newton (pSSN) scheme with the following steps with = 0, 1, 2, . . ., and starting with a given initial guess u 0 ∈ U ad . In the following, we discuss the step s 1 and thereafter s 2 . Concerning the step s 1 , we see that the main computational effort in the procedure (4.6) would be the assembly and inversion of the Jacobian J (u ), but this is not possible because of the size of the problem. In fact, in PDE optimization, one implements the action of the Jacobian (reduced Hessian) on a vector and uses a Krylov approach. Thus, we replace the step s 1 in this procedure with the step: Solve by a Krylov method (e.g. minres) to a given tolerance. In this way, we have a projected SSKN scheme.
Next, we illustrate how the action of the Jacobian on the increment u is constructed. For this purpose, we determine the second-order directional derivative of our Lagrange functional (2.10) with respect to u, and this requires to consider the linearizations of the forward and adjoint equations with respect to u. We have the following where L(ρ, u) := ∂ t ρ + div (a(u) ρ) represents the Liouville operator, and * means adjoint. In (4.7), the functionρ is the solution of the following linearized Liouville problem ∂ tρ + div (aρ) = −div (â ρ), withρ |t=0 = 0, (4.8) whereâ = ∂a ∂u u. Equation (4.8) is obtained in the following way. First, define a small variationρ of ρ such that (ρ +ρ) ∈ C T (L 2 (R d )) and (ρ) |t=0 = 0. Then, insert ρ +ρ for ρ in equation (2.4), use the linearity of (2.1) with respect to ρ and take into account that ρ itself solves (2.1).
To complete the discussion of (4.7), we explain how to computeq. It is obtained solving the following linearized adjoint problem, resulting from a linearization procedure similar to that forρ. We have −∂ tq − a · ∇q =â · ∇q, withq |t=T = 0. (4.9) We solve (4.8) and (4.9) with our KTS scheme. For further details on the implementation of the action of the Jacobian on a vector, we refer to, e.g., [5], Chapter 6.3.5.
Specifically, for our case one can verify that (4.7) is explicitly given componentwise by where the components of are solutions to the following boundary-value problem (4.10) We approximate this problem by finite differences for the time derivative, which results in a tridiagonal linear system, and solve it with the Thomas algorithm. Compare this boundary-value problem with (2.14). Notice that, in general in optimal control problems, the reduced Hessian has a favourable spectral structure, in the sense that it is spectrally equivalent to a secondkind Fredholm integral operator, and in this case Krylov solvers can converge in a mesh-independent number of iterations [1]. We remark that our numerical experience is consistent with this statement. Moreover, it is supported by our estimate of the condition number of the Jacobian. For this purpose, we use the power method to estimate the largest eigenvalue of J (u), and the inverse power method to estimate its smallest eigenvalue, and obtain an approximation to the spectral condition number of the Jacobian. (Notice that these methods also do not require the assembly of J (u).) We compute this condition number in correspondence to different mesh sizes and report these values in Table 1. A plot of these values for different mesh sizes is shown in Fig. 1, where we see that the condition number has similar values for different mesh sizes. Next, we discuss step s 2 of (4.6). This step is required to ensure that any control update results in a control function in U ad . To implement the H 1 projection, P U ad , we solve the following optimization problem miñ u∈U ad Since H 1 T is a Hilbert space and U ad is non-empty, closed and convex, we know that there exists a unique projection (see [33], Theorem 4.11). This problem can equivalently be written as follows Notice that, corresponding to this optimization problem, we have the following Lagrange functional with Lagrange multipliers q a , q b , To solve this optimization problem to implement the projection P U ad , we use a gradient descent scheme; see, e.g., [4], Section 2.8.
Next, we summarize the Newton procedure to solve our Liouville ensemble optimal control problem in the following algorithm that determines the reduced gradient at a given u. Notice that δ enters in step 3 of Algorithm 4.1 throughλ. With Algorithm 4.1, we can define our projected semi-smooth Krylov-Newton algorithm as follows.

Algorithm 4.2 Projected semi-smooth Krylov-Newton method
Require: u 0 1: Set = 0, E > tol 2: while E > tol and < max do 3: Compute ∇ J (u ) with Algorithm 4.1 4: Solve J (u ) u = − ∇ J (u ) (we use minres; here we need to solve (4.8), (4.9), (4.10)) 5: Set u +1 = P U ad u + s u , where s is determined by the Armijo linesearch-backtracking scheme. We solve (4.12) within the backtracking scheme. In this algorithm, we use the difference between consecutive iterations of the control as termination criterion, specifically to stop the algorithm, if the difference is less then a threshold tol > 0. Moreover, we define a maximum number of iterations max ∈ N.

Numerical experiments
In this section, we present results of numerical experiments to validate the accuracy of our numerical framework and to demonstrate the ability of the ensemble optimal control in driving the density in order to perform a achieve given objectives.
We have proved second-order accuracy of our SSPRK2-KT scheme for the Liouville equation in Theorem 3.1. In order to validate this estimate, we define a setting that admits an exact solution. Thus, we choose the following control function u(t) = 0.05 t 0.002 0.5 −0.001 , Further, we take the initial condition where v 0 = 1 4 . With this setting, the Liouville problem ∂ t ρ + div (aρ) = 0, with ρ |t=0 = ρ 0 , admits the solution where the mean m(t) = (m 1 (t), m 2 (t)) and the variance v = (v 1 (t), v 2 (t)) are the solutions to (2.6) with the initial conditions m(0) = (0, 0) and v(0) = (1, 1). Now, we use this setting to determine the solution error of our algorithm. For this purpose, we solve the corresponding Liouville problem and report the values of the discrete L 1 norm of the solution error given by In Table 2, the values of e K T corresponding to different grids are presented, and in Fig. 2, we compare the rate of change of these values with that of first-and second-order accuracies. We see that the obtained numerical accuracy lies between these reference rates, becoming closer to second-order by refining the mesh size.
Next, we validate our estimate for the KTS scheme in solving a transport problem with source term (the adjoint problem) as given in Theorem 3.2. We proceed in a whereã(x, t) = −a(x, t), θ =ρ∇a and −ϕ = ρ 0 . Thus, we have the solutionq =ρ.
However, notice that the KTS scheme uses the Strang splitting in order to accommodate the source term −θ . Therefore, the solutionq =ρ is appropriate to independently test the KTS scheme. Thus, we define Hence, we perform a second series of experiments where we compute the values of this norm in correspondence to solutions obtained on different grids. These values are reported in Table 3, and in Fig. 3, we compare the rate of change of e K T S with that of first-and second-order accuracies. Also in this case, we see that the resulting rate of convergence is approximately of second-order.
Next, we validate the ability of our optimization framework to construct controls that steer the ensemble density in order to follow a desired path. For this purpose, we start considering the tracking of a piecewise smooth trajectory with an initial density given by an unimodal distribution. Thereafter we demonstrate that our approach allows to construct control functions that are able to drive the evolution of the density with a bimodal structure.
This ρ 0 ∈ C ∞ ( ) represents a unimodal Gaussian distribution centred in ξ 0 = (−0.5, 0.5), with variance σ = 1 In this experiment, the purpose of the control is to drive the ensemble of trajectories along the following piecewise smooth desired trajectory A plot of this trajectory in is given in Fig. 4a. Correspondingly, our potentials in the objective functional are chosen as follows where x ∈ and t ∈ [0, T ], and the values of C θ , C ϕ , σ θ and σ ϕ are given in Table 4. Now, we specify a setting that facilitates a comparison of our results of ensemble control with a simple dynamics for the trajectory. Specifically, suppose that our desired trajectory is the result of the following dynamicṡ Then we can immediately compute the control u 1 in this equation such that the solution to (5.5) is given by (5.4). This control is plotted in Fig. 5b, and we refer to it as the single-trajectory control, specifically taken u 2 ≡ 0, which corresponds to no change in the variance. Notice that this control is not in our control space U ad (recall its definition (2.8) above), since it is not continuous. Moreover, in its construction, we do not require to satisfy the conditions of its value being zero at initial and final times. In our drift (2.5), we choose a 0 = 0, b = 1 and c = 0, and with this setting we solve our Liouville control problem, taking the numerical values given in Table 4. The resulting control function is depicted in Fig. 5a, which appears similar to the singletrajectory control in Fig. 5b. We see that the former is in H 1 T and is zero at t = 0 and t = T as required, we refer to it as the Liouville control.
Corresponding to the Liouville control, we obtain an evolution of the density with which we compute the function E[x](t) = x ρ(x, t) dx. This function is shown in Fig. 4b. Notice that it closely resembles the desired trajectory.
In our second experiment, we consider the setting a 0 = 0, b = 1 and c = 1, an we take a smooth initial ρ 0 that is given by a bimodal Gaussian distribution as follows The values of the other parameters are specified in Table 5.
In this case, we choose the following desired trajectory We have that ξ D (0) corresponds to the midpoint between the centres of the two Gaussians defining the initial density; see Fig. 6a, where we plot circles around the centres of the two Gaussians with radius of their standard deviation. With this setting, we solve our Liouville optimal control problem and obtain the controls shown in Fig. 6b. The values of C θ , C ϕ , σ θ and σ ϕ together with the values of the numerical parameters are given in Table 5.
Corresponding to these controls, we obtain the evolution of the density depicted in Fig. 7a. Specifically, we plot the shape of the density ρ at all times. One can see that the bimodal density is driven towards the desired trajectory becoming unimodal. The same result is visualized in Fig. 7b from a different perspective.
We would like to conclude this section considering a setting that generalizes our framework. Our purpose is to demonstrate that our control framework is also able to drive a smooth bimodal distribution to follow two trajectories. In this case, we choose Fig. 7 Evolution of the density in the second experiment a potential θ that resembles a double well, so that it provides two basins of attraction corresponding to the two trajectories.
In this experiment, we consider the initial bimodal ρ 0 given in (5.6), and consider the following two desired trajectories T In correspondence to these trajectories, we define θ as follows Similarly, we define ϕ(x) = θ(x, T ). We solve the resulting ensemble control problem with Algorithm 4.2 and obtain the controls depicted in Fig. 8. In correspondence to these controls, we obtain the evolution of the initial bimodal density shown in Fig. 9b. We see that the two initial Gaussians are driven along the desired trajectories ξ D1 and ξ D2 shown in Fig. 9a and merge at the final time.

Conclusion
A numerical optimization framework to solve non-smooth ensemble optimal control problems governed by the Liouville equation was presented. In these problems, the cost functional has the structure of an expected value function for the ensemble of trajectories and includes L 2 , L 1 and H 1 costs of the controls.
The solutions to these Liouville control problems were characterized as solutions to first-order optimality systems consisting of the Liouville equation, its optimization adjoint, and a variational inequality in H 1 . In order to numerically solve the Liouville problem, a strong stability-preserving Runge-Kutta method (SSPRK2) in time and a Kurganov-Tadmor (KT) scheme in space were discussed, proving that the combined SSPRK2-KT scheme provides a second-order accurate solution. Moreover, positivity and conservativeness of the SSPRK2-KT scheme were also shown under a CFL condition. To solve the adjoint Liouville equation, a Strang splitting technique combined with the SSPRK2-KT scheme and exact integration were proposed. The resulting Kurganov-Tadmor-Strang scheme was shown to be second order accurate. The above approximation schemes were used to implement a H 1 -projected semismooth Krylov Newton method to solve the Liouville optimal control problems. Results of numerical experiments were presented that successfully validated the proposed computational framework and the ability of the resulting control functions to drive the ensemble of trajectories in the case of both unimodal and bimodal initial distributions.