{'fraction_non_alphanumeric': 0.0499866345896819, 'fraction_numerical': 0.016617660162166977, 'mean_word_length': 4.133089412302767, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Discrete-Time Linear-Quadratic Regulation via Optimal Transport

In this paper, we consider a discrete-time stochastic control problem with uncertain initial and target states. We first discuss the connection between optimal transport and stochastic control problems of this form. Next, we formulate a linear-quadratic regulator problem where the initial and terminal states are distributed according to specified probability densities. A closed-form solution for the optimal transport map in the case of linear-time varying systems is derived, along with an algorithm for computing the optimal map. Two numerical examples pertaining to swarm deployment demonstrate the practical applicability of the model, and performance of the numerical method.


I. INTRODUCTION
The problem of steering the states of a linear system from an initial distribution to a terminal distribution has attracted much interest in recent years [1]- [3]. Applications of such controllers include the density control of swarms [4], [5] and networked dynamical systems [6], as well as opinion dynamics [7].
Optimal mass transport is a mathematical framework for deriving mass-preserving maps between specified distributions that minimize a cost of transport. The optimal transport cost, in some specific contexts called the Wasserstein metric, provides a useful metric on the space of probability distributions. This has been employed in a wide variety of fields, such as economics [8], machine learning [9], [10], computer vision [11], and image processing [12]. The Wasserstein metric also allows one to tractably compute worst-case distributions in optimization problems [13], which have been applied in areas such as state estimation [14], and machine learning [15]. The computation of the Wasserstein metric and corresponding transport map has also attracted much attention, in particular techniques allowing for computational speedup such as entropic regularization and Sinkhorn scaling [16], [17].
The connection of optimal transport to continous-time control began with the seminal reformulation of optimal transport as a PDE-based fluid dynamics optimization problem [18]. In this approach, a velocity field is computed that minimizes the average kinetic energy of a fluid moving from one density to another. Equivalently, this approach can be thought of as a single-integrator particle moving from an initial state with uncertainty described by an initial distribution, to a final state with an uncertainty described by a final distribution. The cases of general linear time-varying (LTV) systems, and general LTV systems driven by noise (so-called Schrödinger bridges) were developed by [19].
The latter paper [19] employs a Lagrangian-based cost function, where the static quadratic cost is replaced with a time-varying cost with dynamical constraints. Such techniques were developed in [20], which dealt with optimal transport with nonholonomic constraints. In a similar problem configuration, the existence and uniqueness of transport maps were determined for linear-quadratic costs by [21]. Other works include distributed optimal transport for swarms of single-integrators [22], [23], Perron-Frobenius operator methods for computing optimal transport over nonlinear systems [24], and a related problem regarding the steering of an LTV systems to a terminal state with specified expected value and covariance [25]- [27].
While much attention has been paid to optimal transport of dynamical systems in continuous-time, there has been a marked lack of works discussing the implementation of such controllers in discrete time, which is a gap in the literature that needs to be addressed before optimal transport techniques can be implemented on digital controllers. One contribution of this paper is to provide a rigorous analysis of the optimal transport problem for linear-quadratic regulation of LTV systems in discrete time.
In the present work, we discuss the theory and implementation of optimal transport for discrete-time linearquadratic regulation for LTV systems. Our contributions are as follows. We formalize a previously-developed method of applying optimal transport methods to control by converting a class of optimal control problems to an optimal transport problem where the cost function is the optimal cost-to-go from an initial state to a terminal state. This formalism is then applied to derive the closed-form solution of the discrete-time LQR problem with state-density constraints. This problem is solved numerically, and the solution is then implemented on an example involving swarm deployment.
The paper is organized as follows. We outline the notation and preliminaries on optimal transport in §II. Our problem statement is outlined in §III, where we discuss formulating optimal transport problems for control systems in terms of value functions. Our results concerning optimal transport for LQR and its numerical computation are in §IV. We present numerical examples and an application to swarm deployment in §V, and conclude the paper in §VI.

II. MATHEMATICAL PRELIMINARIES
In this section, we outline the notation used in the paper, as well as the necessary preliminaries on optimal transport.

A. Notation
The n-dimensional space of real numbers, non-negative real numbers, and positive real numbers are respectively denoted by R n , R n + , and R n ++ . We denote vectors in lowercase x, y, z ∈ R n , and matrices in capital-case A, B, C ∈ R n×m . Inequalities are interpreted component-wise. Symmetric positive-definite and positive semi-definite cones of matrices are respectively denoted as S n ++ and S n + . For Q ∈ S n + , we let x T Qx = x 2 Q . The n × n identity matrix is denoted by I n , or just I if comformable dimensions are assumed. 1 n denotes the length-n vector of all ones, and 0 denotes a matrix of zeros. The identity map is denoted by Id(x) = x. The direct sum of n m × m square matrices matrices {A i } n i=1 is the nm × nm matrix formed by placing A 1 , . . . , A n on the block diagonal. It is denoted by A measure space is a triple (X , A, µ) where X is a set, A is a σ-algebra on X , and µ is measure on (X , A). We write a probability space as (X , µ), where µ(X ) = 1 is a nonnegative measure, and we assume that (X , µ) is equipped with the Borel σ-algebra.
For probability spaces (X , µ 0 ), (Y, µ 1 ), the pushforward map, denoted by µ 1 := T # µ 0 , is defined by the relation for each B ∈ A(Y ). If a random variable x is distributed according to a probability density function ρ(x), then we write x ∼ ρ(x)dx.

B. Optimal Transport
In this section, we summarize four seminal forms of the optimal transport problem, and then specify the form of the optimal transport for our present work. One may consult excellent texts by Villani for a more in-depth discussion of the theory [28], [29].
Consider two probability spaces (X , µ 0 ) and (Y, µ 1 ). A transport map T : X → Y is said to transport µ 0 to µ 1 if T # µ 0 = µ 1 . The Monge optimal transport problem seeks to find an optimal map T that minimizes some cost of transport c(x, T (x)), In general, if one of the measures µ 0 , µ 1 has infinite second moment, then the cost of (OT1) may be infinite. Furthermore, the pushforward constraint of (OT1) makes this problem computationally intractable. Kantorovich formulated a relaxation of (OT1) that obtains the same minimizer under quadratic costs 1 , i.e., c(x, y) = 1 2 x − y 2 2 . The problem considers the set of joint probability distributions π(x, y) on X × Y whose marginals are the initial and target measures, for all Borel sets A ⊆ X and B ⊆ Y. With some abuse of notation, to make variables of operators (e.g., optimization, integration) we may write the above as π(x, ·) = µ 0 (x), π(·, y) = µ 1 (y).
The Kantorovich optimal transport is then given by, For the case of quadratic costs, (OT2) obtains the same minimum as (OT1), and the optimal coupling satisfies π The dual of (OT2) has an explicit interpretation in economic theory of transport pricing [8], but perhaps more importantly, it offers insight into the structure of the optimal map T in the case of quadratic costs. For φ, ψ in the dual space of probability measures, the dual is given by, When c(x, y) = 1 2 x−y 2 2 , the optimal map T * (x) of (OT1) can be written in terms of ψ * from (OT3) as [29], and in particular it can be shown that ( 1 2 x T x + ψ * (x)) is a convex function [30]. Note that in our notation, ψ * (x) refers to the optimal ψ, and not its Fenchel conjugate.
One final formulation of optimal transport we describe here is given by Brenier and Benamou in the form of an optimal control problem in a fluid dynamics setting. Given initial and terminal densities ρ 0 , ρ 1 , one seeks to find a smooth, time-dependent velocity field v(x, t) taking ρ 0 to ρ 1 in unit time, while satisfying the continuity equation. The velocity field minimizes the average kinetic energy of the fluid. The problem is explicitly defined as [18], , the solution to (OT4) is given by a linear interpolation with the optimal map, and so the densities at time t satisfy

III. STOCHASTIC OPTIMAL CONTROL WITH STATE-DENSITY CONSTRAINTS
In this section, we consider an optimal transport approach for the discrete-time linear-quadratic regulator. We present a formal discretization of the continuous-time controllers presented in [19], and extend this to the more general framework of LQR control.
We consider systems with a state z k ∈ R n of the form where the initial condition z 0 has some uncertainty described by a probability density ρ 0 (x) and u k ∈ R m is the control.
Our goal is to translate the system (1) to a terminal state z t f ∼ ρ 1 over a time horizon 0 ≤ k ≤ t f , where ρ 1 captures some desired uncertainty in the terminal state 2 . The control should satisfy some optimality principle under an appropriate cost, and so an optimization problem with dynamics (1) is, where the expectation is with respect to a joint distribution π(z 0 , z t f ), as defined in (OT2). The remark below formalizes a solution technique for problems of the form (2) which was used by [19] to solve continuous-time optimal control problems with control costs. Remark 1: A general method to solve problems of the form (2) is to first solve the deterministic problem to determine a formula C(x, y) for the optimal cost-to-go from x to y. Thus, (2) can be re-written as where the marginal constraints on π encode the relevant distributions on the initial state x and terminal state y. Problem (3) is clearly a Kantorovich optimal transport problem of the form (OT2), where the cost function is now the deterministic value function encoding the cost-to-go from x to y, and the optimal coupling π * encodes a mapping between initial and terminal conditions x and y. The solution to (3), under appropriate assumptions on the cost C(x, y), yields a coupling of the form We solve this problem in the following section.
IV. DERIVATION OF THE OPTIMAL MAP In this section, we outline the solutions to Problem (4), beginning with the simplified case of a cost on the control only. Our main contribution in this section is the moregeneral LQR problem, outlined in IV-B.

A. Discrete-Time Optimal Transport -Control Cost Case
Consider the task of transporting over a time horizon of 0 ≤ k ≤ t f a linear time-varying system with uncertain initial state z 0 characterized by ρ 0 (z), to a final state z t f with an uncertainty characterized by ρ 1 (z), with minimal control cost. We assume that the dynamics z k+1 = A k z k + B k u k are controllable over the interval 0 ≤ k ≤ t f . The problem is posed as Following a similar derivation as the continuous-time case studied in [19], one can consider first solving the deterministic problem, The solution to (P2) is given in closed form as where Substituting the optimal cost (5) into (P1) and applying the coordinate transformation transforms (5) into x −ŷ 2 2 , and so we arrive at the Kantorovich optimal transport problem minπ R n ×R n x −ŷ 2 2 dπ(x,ŷ) s.t.π(x, ·) =ρ 0 (x)dx,π(·,ŷ) =ρ 1 (ŷ)dŷ, with the distributions in the (x,ŷ) coordinates given bŷ Now, suppose thatT is the Monge map that corresponds to the solution of (P3). Then, by using (7), the solution to the original problem (P1) can be written in terms of its Monge map The optimal controls are thus given by,

B. Discrete-Time Optimal Transport -Linear-Quadratic Case
In this subsection, we consider the more general case of a linear-quadratic cost function. The problem is formulated as We proceed using the methodology in Remark 1 by considering the solution to the deterministic problem, We summarize the cost-to-go of (P4) in the following lemma. Note that we utilize a pseudoinverse, present in (12). While at first glance it may seem that this pseudoinverse severely limits the applicability of this lemma, this is not the case. We discuss in Remark 3 (after the proof of the lemma) why the pseudoinverse is well-defined for all controllable LTV systems, and we highlight an example in §V that shows that the pseudoinverse is indeed well-behaved, even in pathological cases.

Lemma 1:
The optimal cost-to-go of (P4) is quadratic in x, y, in that where Q xy is invertible. The optimal control of (P4) is given by where the matrices comprising the optimal cost and control are given by where S 1 , and S 2 are . . .
Remark 2: For an LTI system, these matrices are simply Proof: An analysis in the simpler case of LTI systems with y = 0 may be found in [31]. By setting U (l) = vec({u k }) l k=0 , U = U (t f ), andR = t f −1 k=0 R k , the control cost term can be written as Similarly, by writing z k = Φ(k, 0)x + Υ(k, 0)U (k) andQ = t f −1 k=0 Q k , the state cost term can be written as To eliminate the equality constraints, we can write, and we can thus parameterize U as follows: = Substituting (13) into (10) and (11) yields the total cost as Taking the gradient of (14) with respect to U 1 , setting it to zero, and solving for U * 1 yields, . Substituting this form of U * 1 into (13) yields the optimal control as in (9), and substituting U * 1 into (14) yields the optimal cost as in (8). It can be checked that Q xy is positive definite.
Remark 3: Our technical assumption in the above lemma is that the LTV system is controllable, in the sense that the controllability Gramian is positive definite. The matrix S 2 in (12) may be zero, even in cases when the underlying system is controllable. The pseudoinverse is well-defined in this case due to the elementary property 0 † m×n = 0 n×m . See the example in §V for more details.
△ We can now state and prove the main theorem regarding the solution of Problem (P3).
Theorem 1: Consider the setting of Problem (P3), and the Kantorovich optimal transport problem, where C * (x, y) is given by (8). Then, the optimal coupling π * of (15) is given by where T * # ρ 0 = ρ 1 . Furthermore, the control inputs optimizing Problem (P3) is given by, where the relevant matrices are defined in Lemma 1.

Proof:
We can see that the Kantorovich problem (15) corresponds to Problem (P3). By Theorem 2.2 in [21], the solution to Problem (15) under the cost C * (x, y) from (8) exists, and is of the form where ϕ(x) is a convex function. From the calculation in Lemma 1, the optimal control in (16) follows.

Remark 4:
The solution to Problem (15) yields a Monge map T * that transports x ∼ ρ 0 (x)dx to y := T * (x) ∼ ρ 1 (y)dy, minimizing the expected cost-to-go from x to y. Another interpretation of this map is that it pairs initial and terminal states (x, y) in such a manner that it minimizes the LQR cost averaged over the distribution of initial states. We exploit this interpretation in §V, where we discuss an application to swarm deployment. △ First, we examine the numerical computation of T * (x).

C. Numerical Computation of the Monge Map
In general, the Monge map T (x) is difficult to compute numerically [32], [33]. In fact, the optimal control formulation (OT4) was devised by Brenier and Benamou precisely to numerically compute T (x), and devising fast solvers for this problem is an area of active research [34]. In one dimension, a classical result (used in [19]) determines the Monge map in terms of the cumulative distribution functions of the initial and terminal densities as This can readily be solved to high precision with a bisection algorithm.
For systems with n > 1 states, the situation is more complicated. For example, in the single-integrator system x k+1 = x k + u k , the one-timestep Monge map exists explicitly when the initial and termimal distributions are Gaussian. Suppose N (m 1 , Σ 1 ). Then, the optimal Monge map is a shift and scaling [35], For general distributions, we outline a discretization-based method for computing π * from (OT2), and then generating the image of T * (x) from this approximate π * . Suppose we discretize X = Y := R n into cells , and then define probability mass vectors ρ 0 ∈ R nx , ρ 1 ∈ R ny representing ρ 0 , ρ 1 , as The cost in (OT2) can be written over this discrete space as where x i , y i are representative coordinates of the cell, say their centroids. The marginal constraints can be imposed on π as π1 ny = ρ 0 , π T 1 nx = ρ 1 .
Letting C(x i , y j ) := C ij , and π(x i , y j ) := π ij , we arrive at the linear program To recover a discrete image of the map T , one has to numerically 'un-do' the pushforward operation that π represents. This is done by the element-wise division of π(x i , y j ) by ρ 0 (x i ): .
Note that this definition requires that ρ 0 must be strictly positive over the discrete domain; alternately if ρ 0 (x i ) = 0, then the corresponding row of π must also be 0 from the constraints in (17). In this case, we can define τ (x i , y j ) arbitrarily, since if there is no mass to move from x i , it is irrelevant where that mass should move to. Note that Problem (17) suffers from the 'curse of dimensionality' due to the discretization of R n . Fast approximations of optimal transport are an ongoing area of research, and in the near future one may expect that Problem (17), or approximations of it, may be computationally tractable for large statespaces [16], [36], [37]. The graph of T over {x i } k i=1 can then be determined by applying the map τ to the domain {x i } k i=1 . Suppose [X 1 , . . . , X n ] are appropriately-vectorized coordinates in each of the n directions of the discretized domain in R n . Then, the matrix τ generates the image of T as follows: where T i is the vectorized map over the domain in the ith direction of R n .

V. EXAMPLES
In this section, we provide numerical experiments of the results in §IV. Code (and parameters) for the examples can be accessed at [38]. The runtimes for the computation of the optimal transport maps are 12s and 9s, for examples V-A and V-B respectively, on an Intel Core i7-9700K CPU (3.60GHz).

A. 2D LQR Example on LTI System
First, we provide an example of the numerical implementation of Theorem 1. We implemented the optimal transport method for LQR on a 2-state, 1-input system with matrices Let the states be denoted by z k := (z k ). Our initial states were distributed according to ρ 0 (x) depicted in Fig. 1, and we sought to steer the system to the distribution ρ 1 (y), also depicted in Fig. 1.
The distributions ρ 0 , ρ 1 are supported on a discrete grid on the cube [−1, 1] 2 with a discretization length ∆x = ∆y = 0.0588. The optimal transport map T * was computed by solving (17) discretized on this grid, with a cost matrix C computed using (8) from Lemma 1. In Fig. 2, we show a color map of the image of T * . On the left is the coordinate in the first dimension as a function of the initial (x (1) , x (2) ) = (z (1) 0 , z (2) 0 ), and on the right is the coordinate in the second dimension as a function of the initial (z (1) 0 , z (2) 0 ). One can note that the numerical approximation contains outliers in regions where ρ 0 has little mass.
We ran an experiment with 10,000 i.i.d. random initial conditions sampled from ρ 0 . For each initial condition z 0 , the optimal map computed the corresponding final condition as z t f = T * (z 0 ). The simulation then used the optimal control inputs (16) to guide the system to T * (z 0 ) over a time horizon of 0 ≤ k ≤ 10. Plots of the empirical distributions are shown in Fig. 3 at times t = 0, 3, 7, 10.

B. Swarm Deployment -2D LQR on an LTV System
Consider the task of assigning target positions to n agents whose initial states have an empirical distribution approximating ρ 0 (x)dx, but not randomly instantiated. For example, consider n agents spaced at constant intervals in the cube [−1, 1] 2 , as depicted in the middle-left subfigure of Fig. 5. Clearly, this is an approximation of a uniform distribution. Our target distribution is the logo of the Swiss Federal Institute of Technology, Zürich, discretized over a 35 × 35 pixel domain. A target application could be a swarm of UAVs  providing a background performance act during a university event.
Using LTV discrete single-integrator dynamics A k = Q k = R k = I 2 , 0 ≤ k ≤ 10, we compute the optimal map using (17) with the cost matrix (8) from Lemma 1, depicted in Figure 6. The simulation was again produced over a time horizon of 0 ≤ k ≤ 10. This time, we plot the explicit mapping between points in a grid and their target states as generated by the map T * , as shown in the bottom-right of Figure 6. The dynamics (18) are controllable, in the sense that the controllability Gramian (6) W c (t f , 0) is positive-definite, however the matrix S 2 in (12) is 0 4×2 . Since 0 † 4×2 = 0 2×4 , by (12), this simply means that the control is zero for 6 ≤ k ≤ 10. As the system is controllable, it is steered to the final position by timestep k = 5, as evident in Figure 4.

VI. CONCLUSION
In this paper, we studied the discrete-time linear-quadratic regulator with uncertainties in the initial state, and how optimal transport can be used to guide the system to a final state with an uncertainty specified by a target probability density. We derived the form of the optimal control from the optimal transport map, and discussed numerical implementations of this map. Finally, we provided numerical examples with an application to swarm deployment.
Future work may include online computation of the optimal transport map corresponding to the LQ cost, and studying systems where additional uncertainty comes from being driven by noise of arbitrary distributions.