Multiple sparsity constrained control node scheduling with application to rebalancing of mobility networks

This paper treats an optimal scheduling problem of control nodes in networked systems. We newly introduce both the L0 and l0 constraints on control inputs to extract a time-varying small number of effective control nodes. As the cost function, we adopt the trace of the controllability Gramian to reduce the required control energy. Since the formulated optimization problem is combinatorial, we introduce a convex relaxation problem for its computational tractability. After a reformulation of the problem into an optimal control problem to which Pontryagin's maximum principle is applicable, we give a sufficient condition under which the relaxed problem gives a solution of the main problem. Finally, the proposed method is applied to a rebalancing problem of a mobility network.

problem finds not only which but also when nodes should become control nodes, and hence it seems more challenging and efficient for achieving high control performance. Indeed, the authors in [11], [12] consider the scheduling problem for discrete-time systems and show its effectiveness over the time-invariant control node selection. Mathematically, all of the aforementioned works consider ℓ 0 constrained optimization problems, in which the number of nodes selected at the same time is constrained. On the other hand, in networked systems it is also important to effectively compress control signals and reduce communication traffic. To achieve this, it is desirable to find the best time duration over which controllers should become active. Then, we considered the L 0 constraint on control inputs and formulated a node scheduling problem for continuous-time systems in [13]. This scheduling problem is furthermore analyzed in [14], which provides an explicit formula of the optimal solutions and shows that the solutions are obtained by a greedy algorithm. However, these two works on continuous-time systems mainly consider the L 0 control cost, and the resulting number of activated control nodes at each time instance (i.e., the ℓ 0 control cost) is not taken into account.
In view of this, this paper newly considers an optimal node scheduling problem under the L 0 and ℓ 0 constraints. By introducing these two constraints, we can find a time-varying small number of control nodes while reducing the support of control inputs. As the network controllability, we adopt the trace of the controllability Gramian. This quantity is closely related to the average energy required to steer the system in all directions in the state space [15]. The formulated problem thus includes a combinatorial structure caused by the L 0 and ℓ 0 norms. To circumvent this, we introduce a convex relaxation problem and establish a condition for the main problem to be exactly solved via the convex optimization. For the analysis, we transform the convex relaxation problem to an optimal control problem to which Pontryagin's maximum principle is applicable. Unlike the previous formulation in [13], [14], our maximum condition in the principle does not boil down to component-wise calculation and the transversality condition is needed to show the equivalence between the main problem and the relaxation problem.
To demonstrate the practicability of the developed method, we address a rebalancing problem of mobility networks with one-way trips, e.g., car-and bike-sharing systems [16], [17]. In such a system, a customer picks up a vehicle in a station and can return it in another station, which enhances the usability of sharing systems but causes a problem of uneven distribution of vehicles. To maintain this system, rebalancing vehicles is required, which should be as infrequent as possible for the reduction of staff cost. This problem is formulated as optimization problems to solve with conventional techniques of optimization and control engineering [18], [19]. Although in the existing research the amount of rebalanced vehicles has been taken into account, the frequency of rebalancing should be really considered to reduce staff cost. In this paper, we show how to apply the developed method to attain infrequent rebalancing, namely sparse rebalancing, in the mobility network systems, and illustrate its effectiveness through numerical examples.
The remainder of this paper is organized as follows: Section II provides mathematical preliminaries. Section III formulates our node scheduling problem. Section IV introduces a convex relaxation problem and gives a sufficient condition for the main problem to boil down to the convex optimization. A numerical example of the proposed node scheduling is also illustrated. Section V extends the proposed method to a rebalancing problem of mobility networks. Section VI offers concluding remarks.

II. MATHEMATICAL PRELIMINARIES
This section reviews notation that will be used throughout the paper.
We denote the set of all positive integers by N and the set of all real numbers by R. Let m ∈ N and Ω ⊂ R. For a vector a = [a1, a2, . . . , am] ⊤ ∈ R m , diag(a) denotes the diagonal matrix whose (i, i)-component is given by ai, and a ∈ Ω m means ai ∈ Ω for all i. The ℓ 0 norm and ℓ 1 norm of a are defined by a ℓ 0 #{i ∈ {1, 2, . . . , m} : ai = 0} and a ℓ 1 m i=1 |ai|, where # returns the number of elements of a set. We denote the Euclidean norm by For any M ∈ R N×N , TrM denotes the trace of M . For any M ∈ R N 1 ×N 2 , M ⊤ denotes the transpose of M . Let S be a closed subset of R m and a ∈ S. A vector ξ ∈ R m is a proximal normal to the set S at the point a if and only if there exists a constant σ ≥ 0 such that ξ ⊤ (b − a) ≤ σ b − a 2 for all b ∈ S. The proximal normal cone to S at a is defined as the set of all such ξ, which is denoted by N P S (a). We denote the limiting normal cone to S at a by N L S (a), i.e., where µL is the Lebesgue measure on R, and ess sup denotes the essential supremum defined by We denote the set of all functions v with v p < ∞ by L p , p ∈ {1, ∞}. We call a vector-valued function with absolutely continuous components arc [20, p. 255].

A. System Description
Let us consider a network model consisting of n nodes and define the overall system bẏ where x(t) = [x1(t), x2(t), . . . , xn(t)] ⊤ ∈ R n is the state vector consisting of n nodes, where xi(t) is the state of the i-th node at time t; u(t) ∈ R m is the exogenous control input that influences the network dynamics; A ∈ R n×n is the dynamics matrix that represents the information flow among nodes; B = [b1, b2, . . . , bm] ∈ R n×m is a constant matrix that represents candidates of control nodes; v(t) ∈ {0, 1} m represents the activation schedule of the control input u(t); T > 0 is the final time of control. In this setting, the control input uj (t), the j-th component of u(t), is able to affect the system through the vector bj at time t if and only if vj (t) = 1, and the nodes that receive the inputs are called control nodes. In other words, control node scheduling problem seeks an optimal variable v(t) over [0, T ] based on a given cost function and some constraints.

B. Main Problem
This paper is interested in the controllability performance as the cost function. The performance is related to the quantity of the required control energy, for which a number of metrics have been proposed; see e.g. [9], [10]. Among them, this paper adopts the trace of the controllability Gramian, which is a metric to approximate the network average controllability [15]. In addition, this paper introduces the L 0 and ℓ 0 constraints on the control input. The L 0 constraint limits the time duration where the control becomes active, and the ℓ 0 constraint limits the number of control nodes at each time. Thus, the main problem of this paper is defined as follows: In this paper, we will show that Problem 1 is exactly solved via a convex optimization problem. Note that given two optimization problems are said to be equivalent if the set of all optimal solutions coincides. Remark 1: Note that, from a property of the trace operator we have Since V (t) 2 = V (t) for v(t) ∈ {0, 1} m , Problem 1 is equivalent to the following problem: where and bj is the j-th column of B. Remark 2: Compared to existing works [13], [14], our formulation considers the component-wise L 0 norm vj L 0 and includes the ℓ 0 norm of inputs, by which we can adjust each L 0 cost of control variables and the number of control nodes. While the optimal solution in the previous works is shown to be constructed from the top slice of the functions f1(t), f2(t), . . . , fm(t) defined by (3) by using a rearrangement [14], this property does not hold for our optimization problem. This will be illustrated in the example section.

IV. ANALYSIS
The convex relaxation problem of Problem 1 is defined as follows: Problem 2: Given A ∈ R n×n , B ∈ R n×m , T > 0, β ∈ {1, 2, · · · , m − 1}, and αj ∈ (0, T ], j = 1, 2, . . . , m, find a function v(t) [v1(t), v2(t), . . . , vm(t)] ⊤ that solves The set of all functions that satisfy the constraints of an optimization problem is called feasible set. Let us denote the feasible set of Problem 1 and Problem 2 by V0 and V1, i.e., 1} m for all t. Then, we first show the discreteness of solutions of Problem 2, which guarantees that the optimal solutions of Problem 2 belong to the set V0. Proof: Then, for each j, the value vj L 1 is equal to the final state yj(T ) of the systemẏj(t) = vj (t) with yj(0) = 0. Hence, Problem 2 is equivalently expressed as follows: This is an optimal control problem to which Pontryagin's maximum principle [20,Theorem 22.2] is applicable. Let the process (y * , v * ) be a local maximizer of the problem (4). Then, it follows from the maximum principle that there exists a constant η equal to 0 or 1 and an arc q : [0, T ] → R m satisfying the following conditions: (i) the nontriviality condition: (ii) the transversality condition: where S {a ∈ R m : aj ≤ αj, ∀j}, (iii) the adjoint equation for almost every t ∈ [0, T ]: where DyH η is the derivative of the function H η at the second variable y, and H η : [0, T ]×R m ×R m ×R m is the Hamiltonian function associated to the problem (4), which is defined by (iv) the maximum condition for almost every t ∈ [0, T ]: (7) and [20,Theorem 6.41] that there exists a constant γ ∈ R m such that q(t) = γ on [0, T ], since our Hamiltonian does not depend on the second variable y. Then, from (6), there exist sequences {ξi} ⊂ R m and {ωi} ⊂ R m such that For any i ∈ N, by definition, there exists a constant σi ≥ 0 such that Fix any ε > 0, j0 ∈ {1, 2, . . . , m}, and i ∈ N. Take w ∈ S such that wj 0 = ω . . , m} and i ∈ N from the arbitrariness of ε > 0, j0, and i. Hence, γj ≤ 0 for all j by (9).
In addition, we have γj which holds for all i > N . From (9), we have γj 0 = 0, and thus γj 0 (y * j 0 (T )−αj 0 ) = 0. Finally, the supremum in (8) is attained by a point in V, since the right hand side is a linear function of v and V is a closed set.
The following theorem is the main result, which shows the equivalence between Problem 1 and Problem 2.
The existence of optimal solutions of Problem 2 is assumed in Theorem 1 and Theorem 2. Although we omit the proof due to the page limitation, we can show the existence in a similar way to [3, Lemma 1].

A. Example
This section gives an example of our node scheduling. We consider a network model (1)  to each node at most 0.4 sec. Note that the functions fi(t) − fj (t) and fj (t) are shown in Fig. 1 and the upper panel of Fig. 2, by which we can confirm the equivalence between Problem 1 and Problem 2 from Theorem 2. Then, we applied CVX [22] in MATLAB, which is a software for convex optimization, to Problem 2. Fig. 2 shows the resulting time series of control nodes on [0, T ]. Certainly, we can see that the set of control nodes depends on the time and satisfies both the L 0 and ℓ 0 constraints. Thus, we can find a finite number of essential nodes at each time and an essential time interval to provide control inputs. For comparison, we also simulated the previous node scheduling proposed in [13], [14], which try to find a function v(t) ∈ R m that solves where α ∈ (0, mT ) is given. Fig. 3 shows the optimal scheduling when α = 1.6, which is equal to the value 4 j=1 αj . Compared to our framework, the problem (16) [13], [14] (bottom) more than 2 control nodes on an interval and tends to select a particular node. Note also that the feasible set of the problem (16) includes the set V0 due to the absence of the ℓ 0 constraint, by which the optimal value is greater than that of Problem 1, where the optimal values of problem (16) and Problem 1 are 2.1287 and 1.8957, respectively. Finally, the values fj(t) at the switching instances of control nodes are all equal, i.e., we have to take the "top slice" of the functions fj. Actually, this property is shown in [14]. On the other hand, our optimal solution can not be obtained by the simple method, as shown in Fig. 2. Indeed, node 4 is chosen as a control node around the final time T , although the value f4(t) is not ranked in the top 2 among f1(t), . . . , f4(t) at the time.

A. Problem Formulation
In this section, we consider a rebalancing problem on the mobility network of a sharing system with one-way trips. This system is first modeled, and the rebalancing problem is then formulated by the main problem tackled in this paper.
The structure of the network is modeled, based on [19], as follows. Let S = {1, 2, . . . , s} be the set of the indexes of s ∈ N stations, where vehicles are parked. Assume that customers take the service according to the Poisson process at a rate gij (t) ≥ R, the number of the demand for travels from station j ∈ S to i ∈ S per time. In addition, assume that the time it takes to travel from station j to i follows the exponential distribution with an average τij(t) > 0. Let uij (t) ≥ 0 be the rate at which the vehicles are rebalanced by the staff from station j to i. Let fij (t) ≥ 0 be the expectation of the number of vehicles traveling from station j to i, and it varies according toḟ where γij (t) = 1/τij (t) corresponds to the arrival rate of vehicles at station i from j per time. Let vi(t) ≥ 0 denote the expectation of the number of vehicles parked at station i, and it varies according tȯ (gji(t) + uji(t)).
See [23] for detailed derivation of the system. Dynamic pricing is applied to this system, which is modeled as follows. Let pij(t) ∈ R be the price for the rent of a vehicle from station j to i. According to an economic model [24], the amount of the vehicles in service depends on the price as follows: whereḡij(t)≥ 0 is the expected demand of the vehicles when pij(t) = 0, and θij (t)> 0 denotes the price elasticity. Assume that the price pij (t) is determined by the amount of the vehicles at stations i and j according to the following rule: wherepij (t)≥ 0 is a standard price, and λij(t) ≥ 0 and λji(t) ≥ 0 are the coefficients to adjust the price according to the amount of vehicles at stations i and j, respectively. Under this rule, as the amount vi(t) (vj (t)) of vehicles at station i (j) becomes larger (smaller), the price pij (t) becomes higher to reduce the amount of vehicles entering station i (leaving station j). Let the standard pricē pij(t) be determined according to the expected demandḡij(t) of vehicles as follows:p From (19), (20), and (21), (17) and (18) are reduced tȯ vi(t)= s j=1

(γij(t)fij (t) + θji(t)(λji(t)vj(t) − λij (t)vi(t))
−uji(t)), respectively. By collecting (22) and (23) for all i ∈ S, we obtain the model of the mobility network aṡ This is the system of dimension n = s 2 with m = s 2 − s inputs for where Ξ= [Ê1Ê2 · · ·Ês],Êi = [e1 · · · ei−1 ei+1 · · · es], , γ13(t), . . . , γs,s−1(t)), Im ∈ R m×m is the identity matrix of dimension m, and ei ∈ R s is the unit vector of dimension s with the ith entry one. For rebalancing, the staff organizes β ∈ N teams to transport vehicles. Each team goes to a station full of vehicles by their management car, and some members of the team transfer vehicles to a vacant station along with the management car driven by other members. Then, the management car picks up the members to go to another station. Assume that the dynamics of staff is sufficiently faster than that of customers and that only one team can rebalance a vehicle on each route between stations to distribute the staff around the area. Then, the number of rebalances at the same time is at most β, which is expressed as u(t) ℓ 0 ≤ β for any t. Without the loss of generality, the possible rebalance number is one through standardization, which is expressed as |uij (t)| ≤ 1 for any i, j ∈ S and any t, i.e., u L ∞ ≤ 1. Additionally, uij (t) ≥ 0 has to be satisfied. See [23] for dynamic models of rebalancing by the staff.
Let x0 ∈ R n be the initial amounts of vehicles in the stations, which are unevenly distributed, and let x d ∈ R n be the desired terminal amounts to attain even distribution. We design control input u(t) to achieve the state x(T ) = x d at the terminal time T from the initial state x(0) = x0 under the dynamics (24). In summary, the balancing problem of the mobility network is formulated as follows.
Problem 3: Given A(t) ∈ R n×n , B(t) ∈ R n×m , T > 0, x0 ∈ R n , x d ∈ R n , and β ∈ {1, 2, · · · , m − 1}, find a control u that solves We denote the feasible set of Problem 3 by U0. We also denote the state-transition matrix of A(t) by Φ(t, τ ). In other words, Φ(t, τ ) is the unique solution of the matrix differential equation where In is the identity matrix. Remark 3: Note that if x d = Φ(T, 0)x0, then the optimal control is trivial (i.e. zero control), and hence we assume x d = Φ(T, 0)x0 throughout the paper. Note also that we assume β < m, since if β = m the optimal control is a standard time-sparse hands-off control discussed in [2], [3].

B. Analysis
We here introduce a convex relaxation problem of Problem 3, where the L 0 and ℓ 0 norms are replaced by the L 1 and ℓ 1 norms, respectively.
Problem 4: Given A(t) ∈ R n×n , B(t) ∈ R n×m , T > 0, x0 ∈ R n , x d ∈ R n , and β ∈ {1, 2, · · · , m − 1}, find a control u that solves We denote the feasible set of Problem 4 by U1. In other words, Note that we have U0 ⊂ U1, since a ℓ 1 ≤ a ℓ 0 for any a ∈ [0, 1] m . Then, we first show the discreteness of the optimal solutions of Problem 4. The property guarantees that the optimal solutions of Problem 4 belong to the set U0, which is illustrated in the proof of Theorem 4. For this, we introduce an assumption on the impulse response Φ(·, ·)B(·) of the system. Assumption 1: For any nonzero ρ ∈ R n and η ∈ {0, 1}, we have for all j, and for all i, j with i = j, where bj(t) ∈ R n is the j-th column vector in matrix B(t). Theorem 3: Under Assumption 1, the optimal solution of Problem 4 is unique and it takes only the values in the set {0, 1} almost everywhere.
Proof: Let the process (x * , u * ) be a local minimizer of Problem 4. Then, it follows from Pontryagin's maximum principle [20,Theorem 22.2] that there exists a constant η equal to 0 or 1 and an arc q : [0, T ] → R n satisfying the following conditions: 1) the nontriviality condition: 2) the adjoint equation for almost every t ∈ [0, T ]: where DxH η is the derivative of the function H η at the second variable x, and H η : [0, T ] × R n × R n × R m → R is the Hamiltonian function associated to Problem 4, which is defined by 3) the maximum condition for almost every t ∈ [0, T ]: where U {u ∈ [0, 1] m : m j=1 uj ≤ β}.
Hence, we have q(t) = Ψ(t, T )q(T ) on [0, T ] from (27), where Ψ(·, ·) is the state-transition matrix of −A(t) ⊤ . Note also that the supremum in (28) is attained by a point in U, since the right hand side is continuous on u and the set U is closed. Hence, we have We here claim that q(T ) = 0. Indeed, if q(T ) = 0, then it follows from the nontriviality condition (26) that η = 1. From (29), we have u * (t) = 0 for almost all t since q(t) = 0 on [0, T ]. This implies x d = x * (T ) = Φ(T, 0)x0. This contradicts the definition of x0 and x d (see Remark 3). Hence, q(T ) = 0. Since Ψ(t, T ) = Φ(T, t) ⊤ for any t from [25], we have for almost all t ∈ [0, T ] from Assumption 1.
The following theorem is the main result, which shows the equivalence between Problem 3 and Problem 4.
We finally provide the existence of optimal controls of Problem 4 without the proof, which is confirmed in a similar way to [3, Lemma 1]. Since we have U0 ⊂ U1, if the set U0 is not empty, then there exists an optimal control of Problem 3 under Assumption 1. Remark 4: Although we consider non-negative controls in Problem 3, with a slight modification we can show results similar to those presented above (i.e., Theorems 3 and 4) for controls with an L ∞ constraint u L ∞ ≤ 1. To be more precise, under a suitable assumption, the optimal control of a corresponding convex relaxation problem is unique and it takes only the values in {0, ±1} almost everywhere, and this shows the equivalence between the original problem and the relaxation problem as illustrated in Theorem 4. In this sense, our optimal control with multiple sparsity includes the existing sparse control in [2], [3] as a special case.

C. Numerical Study
We demonstrate the effectiveness of the developed method through numerical study. Let s = 10 be the number of the stations. The total number of the vehicles is 200, which implies that i∈S xi(t) = 200 for all t. The location (within a radius of about twenty kilometers) and initial number of vehicles of each station are depicted in Fig. 4. The color on the bottom indicates the degree of the road congestion in The ratio γij(t) is determined depending on the congestion and distance between the stations described in Fig. 4. This system is modeled as (24) with A(t) and B(t) given in (25). The number of the service staff is set to β = 10. For this system, the proposed L 0 /ℓ 0 optimization method is applied for sparse rebalance. Figs. 5 and 6 show the simulation results under this setting. The upper panel of Fig. 5 represents the time plots of the components xi(t), the numbers of the vehicles in the stations, which shows that all xi(t) converge to 20 to achieve the rebalance after 4 hours. On the other hand, the lower panel of Fig. 5 represents those when no control input is applied, namely, uij(t) = 0, where the components xi(t) do not agree at that time. These results show that the proposed method accelerates the rebalance speed. The upper panel of Fig. 6 represents the activated time duration (L 0 norm) of the control input uij (t) when the proposed L 0 /ℓ 0 optimization method is applied, while the lower panel shows that by using an L 2 /ℓ 1 optimization method. These figures illustrate that many of the control inputs are equal to zero (inactive) via the L 0 /ℓ 0 optimization while many of the control inputs via the L 2 /ℓ 1 optimization are non-zero (active). Indeed, the L 0 costs of the optimal controls are approximately 32.08 (L 0 /ℓ 0 optimization) and 136.36 (L 2 /ℓ 1 optimization), respectively. Hence, sparse rebalance is successful due to the proposed method.

VI. CONCLUSION
This paper has analyzed an optimal node scheduling that maximizes the trace of the controllability Gramian. This analysis enables us to find an activation schedule of control inputs that steers the system while saving energy. Taking the number of control nodes and the time length of providing inputs into account, our optimization problem newly includes two types of constraints on sparsity. We have shown a sufficient condition under which our sparse optimization problem boils down to a convex optimization problem. This paper assumes the network topology among nodes is given and fixed. Future work includes the design of the time-varying topology and more practical model of the staff dynamics in the mobility system.

Supplementary material
Takuya Ikeda, Kazunori Sakurama, and Kenji Kashima I. MODEL DERIVATION OF THE TRAFFIC SYSTEM Assume that the customers take the service according to the Poisson process, and that the time it takes to travel between stations follows an exponential distribution. The latter assumption is valid because the farther customers are away from their destinations, the more they have chance of wasting time by shopping or getting lost on the way. Then, the traffic system of customers' vehicles is modeled as (17) and (18) in the main paper, that is,ḟ (g ji (t) + u ji (t)).
See Table I for the meaning and units of the variables, where [units] indicates the unit of the number of vehicles. We derive the system consisting in (r1) and (r2) from the assumptions. Here, the following notation is used. For a random variable X(t) ∈ Z, P[X(t) = k] denotes the probability that X(t) = k occurs. The expectation of X(t) is represented by E[X(t)] := k∈Z kP[X(t) = k].
Let F ij (t) ∈ Z [units] be the number of vehicles traveling from station j to i at time t, and let V i (t) ∈ Z [units] be the number of vehicles parking at station i at time t. These numbers vary within time interval [t, t + δ) for δ > 0 according to (G ji (t, δ) + U ji (t, δ)), where G ij (t, δ) ∈ Z [units] is the number of vehicles of customers departing from station j for the destination i, U ij (t, δ) ∈ Z [units] is the number of vehicles for rebalancing, departing from station j for the destination i, and A ij (t, δ) ∈ Z [units] is the number of vehicles arriving at station i from j within interval [t, t + δ).
We assume the following conditions on G ij (t, δ), U ij (t, δ), and A ij (t, δ). First, the number G ij (t, δ) of the customers is a random variable following the Poisson process  for k ∈ Z with a rate g ij (t) ≥ 0 [units·time −1 ], which is the number of the demand for traveling from station j to i per time. From (r5), the expectation of G ij (t, δ) is calculated as follows: Second, the number U ij (t, δ) of rebalanced vehicles from station j to i is determined with a rate u ij (t) ≥ 0 [units·time −1 ], which leads to E[U ij (t, δ)] = u ij (t)δ.
Third, let p ij (t, δ) ∈ (0, 1) be the the probability of the arrival of the vehicle traveling from station j to i within interval [t, t + δ). Then, the number A ij (t, δ) of vehicles arriving at station i from j within interval [t, t + δ) is given as the random variable following the binomial distribution for k, ℓ ∈ Z, which indicates the probability that k vehicles arrive at station i out of ℓ vehicles. From (r8), the conditional expectation of A ij (t, δ) is calculated as