An optimal control problem for the continuity equation is considered. The aim of a controller is to maximize the total mass within a target set at a given time moment. An iterative numerical algorithm for solving this problem is presented. Copyright ⃝c by the paper's authors. Copying permitted for private and academic purposes. In: Yu. G. Evtushenko, M. Yu. Khachay, O. V. Khamisov, Yu. A. Kochetov, V.U. Malkova, M.A. Posypkin (eds.): Proceedings of the OPTIMA-2017 Conference, Petrovac, Montenegro, 02-Oct-2017, published at http://ceur-ws.org
Consider a mass distributed on Rn that drifts along a controlled vector field v = v(t, x, u). The aim of the controller is to bring as much mass as possible to a target set A by a time moment T .
Let us give the precise mathematical statement of the problem. Suppose that ρ = ρ(t, x) is the density of the distribution and u = u(t) is a strategy of the controller. Then, ρ evolves in time according to the continuity equation where ρ0 denotes the initial density. Our aim is to find a control u that maximizes the following integral Typically, u belongs to a set U of admissible controls. Here we take the following one:
U = fu( ) is measurable, u(t) 2 U a.e. t 2 [0, T ]g , where U is a compact subset of Rm.
In this paper we propose an iterative method for solving problem (1)–(3), which is based on the needle linearization algorithm for classical optimal control problems [Srochko, 2000]. Given an initial guess u0, the algorithm produces a sequence of controls uk with the property J [uk+1] J [uk], for all k 2 N.
A J [u] =
ρ(T, x) dx . (3)
A different approach for numerical solution of (1)–(3) was proposed by S. Roy and A. Borz`ı in [Roy, 2017]. The authors used a specific discretization of (1) to produce a finite dimensional optimization problem. It seems difficult to compare the efficiency of both algorithms, because one was tested for 2D and the other for 1D problems.
Finally, let us remark that problem (1)–(3) is equivalent to the following optimal control problem for an ensemble of dynamical systems:
Maximize ∫ {x : x=y(T )}
In what follows, Φs;t denotes the flow of a time-dependent vector field w = w(t, x), i.e., Φs;t(x) = y(t), where y( ) is a solution to the Cauchy problem {y˙(t) = w (t, y(t)) ,
y(s) = x.
Given a set A Rn and a time interval [0, T ], we use the symbol At for the image of A under the map ΦT;t, i.e., At = ΦT;t(A). The Lebesgue measure on R is denoted by L1.
The map v : [0, T ]
Rn
U ! Rn is continuous.
The map x 7! v(t, x, u) is twice continuously differentiable, for all t 2 [0, T ] and u 2 U . There exist positive constants L, C such that jv(t, x, u) for all t 2 [0, T ], u 2 U , and x, x′ 2 Rn.
v(t, x′, u)j
L x x′j and jv(t, x, u)j j
C (1 + jxj), The target set A Rn is a compact tubular neighborhood, i.e., A is a compact set that can be expressed as a union of closed n-dimensional balls of a certain positive radius r.
In addition, to guarantee the existence of an optimal control
Description 2. For each t, find 3. Let 4. For each ε 2 (0, T ], find 5. Construct uw;I(") by (5). 6. Find 7. Let uk+1 = uw;I(" ).
{∫ w(t) 2 argmin ρ(t, x) v (t, x, ω) nAt (x) dσ(x) : ω 2 U } . g(t) = ∫
v (t, x, w(t))] nAt (x) dσ(x) .
{∫ I(ε) 2 argmax g(t) dt : ι
}
[0, T ] is measurable and L1(ι) = ε .
ε∗ 2 argmax {J [uw;I(")] : ε 2 (0, T ]}.
(4)
(5)
(6)
(7)
(8)
(9)
Theorem 2.1
∫ ρ(t, x) v(t, x, ω) nAt (x) dσ(x) .
Here At = Φt;T (A), where Φ is the phase flow of the vector field (t, x) 7! v (t, x, u(t)), nAt (x) is the measure theoretic outer unit normal to At at x, σ is the (n 1)-dimensional Hausdorff measure.
The precise definitions of the measure theoretic unit normal and the Hausdorff measure can be found, e.g., in [Evans, 1992]. We remark that whenever ∂A is an (n 1)-dimensional surface, nAt coincides with the usual outer unit normal to At and σ coincides with the usual (n 1)-dimensional volume form.
Let I(ε) [0, T ] be a measurable set of Lebesgue measure ε. Given two controls u and w, we consider their mixture uw;I(")(t) = {w(t), t 2 I(ε),
u(t), otherwise.
The proof of Theorem 2.1 gives, as a byproduct, the following increment formula ∫
∫ which will be used in the next section. 3
In this section we describe the algorithm, prove the improvement property J [uk+1] implementation.
J [uk], and discuss a possible 1. Let uk be a current guess. For each t, compute the set ∂At and ρ(t, ) on ∂At.
The algorithm produces an infinite sequence of admissible controls. Of course, any its implementation should contain obvious modifications that would cause the algorithm to stop after a finite number of iterations. Note that it may happen that problems (8) and (9) admit no solution. In this case I(ε) and ε∗ must be taken so that the values of the corresponding cost functions lie near the supremums.
Since the integral from the right-hand side is positive for all small ε, we conclude that J [uk+1] = J [uw;I(" )] > J [uk], as desired. In this step we must compute ρ(t, x) for all t and x satisfying x 2 ∂At. Recall that
ρ(t, x) =
ρ0(y) det DΦ0;t(y) ,
where y = Φt;0(x). d dt
[ ] (det DΦ0;t(y)) = (det DΦ0;t(y)) tr DΦ0;t(y)−1 ddt DΦ0;t(y) .
d dt d dt
DΦ0;t(y) = Dxv (t, Φ0;t(y), u(t)) DΦ0;t(y).
(det DΦ0;t(y)) = (det DΦ0;t(y)) div v (t, Φ0;t(y), u(t)) .
Thus, computing of ρ(t, x) requires solving two Cauchy problems, one for finding Φ0;t(y) and one for finding det DΦ0;t(y).
In general, the optimization problem (7) is nonlinear, which makes it difficult. On the other hand, in many cases
∫
A ρ(T, x) dx = must be computed. To that end, we must know the whole set A0, while on the other steps of the algorithm we deal only with the boundaries of At. It is interesting to note that, under the additional assumption that the target set A
Rn is contractible and its boundary ∂A is an (n
the knowledge of ∂A0 is enough for computing the cost.
Indeed, since the target A = AT is contractible, the set A0 is contractible as well. Any differential form on a contractible set is exact [Tu, 2011]. Hence ρ0 dx1 ^ ^ dxn = da, for some (n 1)-dimensional differential form α. Now the Stokes theorem gives: ∫
A0 ρ0 dx1 ^ ^ dxn = ∫
α.
.
4
4.1
This section describes several toy problems, which we used for testing the algorithm.
Consider a boat floating in the middle of a river at night. Since it is dark, the boatmen cannot see any landmarks, and therefore are unsure about the boat’s position. They want to reach a river island at a certain time with highest probability. How should they act?
the island is a unit circle centered at x0, the initial position of the boat is described by the density function
where u 2 R2 is a component of the boat’s velocity due to rowing. Here juj umax.
Parameters for the computation: σ = 1, α = β = 0.5, umax = 0.75, x0 = ( 3, 0), T = 12. 4.2
Here we want to stop a moving pendulum whose initial position is uncertain. In this case we have (11) , v1(x) = (1)
0 .
x˙ = v0(x) + u v1(x),
where u 2 [ umax, umax] is an external force. The initial position of the pendulum is given by (11). The target is a unit circle centered at (π/2, 0).
Parameters for the computation: σ = 1, umax = 0.5, x0 = (π/2, 0), T = 6. 4.3
Consider a herd of sheep located near the origin. The sheep are effected by a vector field v0(x) pushing them away from the origin. To prevent this we can turn on repellers, which are located at the following positions ( xk =
R cos 2π(k m 1) , R sin 2π(k m 1) ) ,
m ∑ ukvk(x), k=1 where uk is an intensity of k-th repeller. The control u = (u1, . . . , um) belongs to the simplex
{ U = (u1, . . . , um) : m ∑ uk = 1, uk 2 [0, 1], k = 1, . . . , m k=1 }
.
x v0(x) = α √1 + jx x0 x0j2 , where x0 is a certain point not far from the origin, and vk(x) = β e−|x−xk|4 (x xk), k = 1, . . . , m.
Suppose that the initial distribution is given by (11), the target is an ellipse centered at x0 whose major and minor semi-axes are a and b.
Parameters for the computation: σ = 1, x0 = (0, 0), T = 3, m = 6, a = 2, b = 1.2.
Remark 4.1 The answer to the minimization problem
ω 2 U, arising in the second step of the algorithm, is very simple. Let j be such that then an optimal solution is given by ω¯ = (0, . . . , 0, 1, 0, . . . , 0), where 1 is located at the j-th position. In particular, this means that at every time moment t only one repeller is turned on. Hence instead of repellers, we may think of a dog that jumps from one place to another.
The work was supported by the Russian Science Foundation, grant No 17-11-01093.