Paolo Mason
email: paolo.mason@centralesupelec.fr

Yacine Chitour
email: yacine.chitour@centralesupelec.fr

Mario Sigalotti
email: mario.sigalotti@inria.fr

On universal classes of Lyapunov functions for linear switched systems *

In this paper we discuss the notion of universality for classes of candidate common Lyapunov functions of linear switched systems. On the one hand, we prove that a family of absolutely homogeneous functions is universal as soon as it approximates arbitrarily well every convex absolutely homogeneous function for the C 0 topology of the unit sphere. On the other hand, we prove several obstructions for a class to be universal, showing, in particular, that families of piecewisepolynomial continuous functions whose construction involves at most l polynomials of degree at most m (for given positive integers l, m) cannot be universal.

Introduction

Common Lyapunov functions constitute the most popular and powerful tool for the stability analysis of switched systems. Roughly speaking, the use of common Lyapunov functions for stability analysis gathers the global behavior 1 of the system and allows to bypass the explicit analysis of single trajectories, which may be extremely complex. Yet, looking for a common Lyapunov function may be a nontrivial task as stability cannot always be checked by means of Lyapunov functions in a simple form, for instance within the class of quadratic forms. Given a family of systems (e.g., the family of all linear switched systems), classes of functions large enough to include a Lyapunov function for each globally asymptotically stable system are called universal [START_REF] Blanchini | A new class of universal Lyapunov functions for the control of uncertain linear systems[END_REF] and a result establishing the existence of such a class is a called a converse Lyapunov theorem. The literature dealing with converse Lyapunov theorems, starting from the works by Massera and Kurzweil in the 1950s (see e.g. [START_REF] Massera | Contributions to stability theory[END_REF][START_REF] Kurzweil | On the inversion of the second theorem of Lyapunov's second theorem on stability of motion[END_REF][START_REF] Wilson | Smoothing derivatives of functions and applications[END_REF][START_REF] Lin | A smooth converse Lyapunov theorem for robust stability[END_REF][START_REF] Hante | Converse Lyapunov theorems for switched systems in Banach and Hilbert spaces[END_REF]) is quite rich. The results concerning the existence of smooth Lyapunov functions for nonlinear systems with global asymptotic stability properties require the development of rather sophisticated techniques. Concerning robust asymptotic stability with respect to a closed invariant set in presence of perturbation terms, converse Lyapunov theorems have been derived in [START_REF] Lin | A smooth converse Lyapunov theorem for robust stability[END_REF]. In the context of switched systems (even in a nonlinear setting) such results establish the equivalence between the global uniform asymptotic stability and the existence of a smooth Lyapunov function. For switched linear systems the construction of a Lyapunov function is much more direct and natural due, essentially, to the homogeneous nature of the system and the equivalence between asymptotic and exponential stability (see e.g. [START_REF] Dayawansa | A converse Lyapunov theorem for a class of dynamical systems which undergo switching[END_REF]). Furthermore, in the linear case and even for the more general class of uncertain systems, it is well-known that the families of piecewise quadratic functions, polyhedral functions, and homogeneous polynomials are universal [START_REF] Molchanov | Lyapunov functions that define necessary and sufficient conditions for absolute stability of nonlinear nonstationary control systems[END_REF][START_REF] Molchanov | Lyapunov functions that define necessary and sufficient conditions for absolute stability of nonlinear nonstationary control systems[END_REF][START_REF] Molchanov | Criteria of asymptotic stability of differential and difference inclusions encountered in control theory[END_REF][START_REF] Blanchini | A new class of universal Lyapunov functions for the control of uncertain linear systems[END_REF]. On the other hand, for every positive integer m, the family of polynomials of degree less or equal than m is not universal even for the simple class of two-dimensional linear switched systems with two modes [START_REF] Mason | Common polynomial lyapunov functions for linear switched systems[END_REF]. Similarly, it is well accepted in the research community (although, to the authors' knowledge, no explicit proof is available) that families of piecewise quadratic and polyhedral functions whose construction involves a uniformly bounded number of quadratic or linear functions cannot be universal. For this reason, all numerical methods investigating the existence of Lyapunov functions within these classes are affected by a certain degree of conservativeness.

The contribution of this paper is twofold. First, we provide a general sufficient condition for a class of functions to be universal (Proposition 3.1), which is a formalization of fundamental ideas already present in [START_REF] Molchanov | Lyapunov functions that define necessary and sufficient conditions for absolute stability of nonlinear nonstationary control systems[END_REF][START_REF] Molchanov | Lyapunov functions that define necessary and sufficient conditions for absolute stability of nonlinear nonstationary control systems[END_REF][START_REF] Molchanov | Criteria of asymptotic stability of differential and difference inclusions encountered in control theory[END_REF][START_REF] Blanchini | A new class of universal Lyapunov functions for the control of uncertain linear systems[END_REF]. As a corollary, we recover the universal classes of functions obtained in these references. We next derive the main results of this paper, which provide some necessary conditions for the universality of classes of functions. The first one, Theorem 4.2, is an abstract result which applies to families of real analytic functions. The fact that polynomials with a uniform bound m on their degree do not form a universal class [START_REF] Mason | Common polynomial lyapunov functions for linear switched systems[END_REF] follows as a simple consequence of this result. Finally, Theorem 4.5 states that families of piecewise-polynomial continuous functions whose construction involves at most l polynomials of degree at most m (for given positive integers l, m) cannot be universal.

Universal classes of common Lyapunov functions

We consider linear switched systems of the form

ẋ(t) = A(t)x(t), t ≥ 0, x ∈ R n , (Σ M )
where the switching law A is an arbitrary function belonging to the space S arb (M ) of measurable functions taking values on a bounded subset M of the set of n × n matrices, denoted by M n (R). We use Φ A (t, s) to denote the fundamental matrix from s to t for (Σ M ) associated with the switching law A so that every solution of (Σ M ) can be written as x(t) = Φ A (t, 0)x(0). We are interested in the following uniform stability properties.

Definition 2.1. The switched system (Σ M ) is said to be

• uniformly stable if there exists C > 0 such that, for every switching law A and t ≥ 0, Φ A (t, 0) ≤ C;

• uniformly exponentially stable if there exist C, γ > 0 such that, for every switching law A and t ≥ 0, Φ A (t, 0) ≤ Ce -γt .

Stability in the previous senses may be assessed through common Lyapunov functions, defined below. Definition 2.2. We say that a continuous function V : R n -→ R + is a nonstrict common Lyapunov function for (Σ M ) if it is positive definite, that is, V (0) = 0 and V (x) > 0 for every x = 0, and V is non-increasing along trajectories of (Σ M ). If, moreover, V is strictly decreasing along nonzero trajectories of (Σ M ), we say that V is a common Lyapunov function for (Σ M ).

Remark 2.3. If V is a (possibly nonstrict) common Lyapunov function and ϕ : R + → R + is continuous, strictly increasing, and such that ϕ(0) = 0, then ϕ • V is also a (nonstrict) common Lyapunov function. In particular, the positive multiple of a common Lyapunov function is a common Lyapunov function and if there exists a absolutely homogeneous common Lyapunov function, then for every α > 0 there exists a absolutely homogeneous common Lyapunov function of degree α.

We state here the version for linear switched systems of the classical direct Lyapunov theorem.

Theorem 2.4. A linear switched system (Σ M ) admitting a nonstrict common Lyapunov function is uniformly stable. If there exists a strict common Lyapunov function for (Σ M ), then the latter is uniformly exponentially stable.

In case the strict common Lyapunov function V : R n → R + in the above theorem is of class C 1 on R n \ {0}, then a standard test for checking the strict decrease of V along non-trivial trajectories of (Σ M ) goes as follows:

∇V (x) M x < 0, M ∈ M , x ∈ R n \ {0}. (1)
Definition 2.5. A set P of functions from R n to R is a universal class of Lyapunov functions if for every bounded set M ⊂ M n (R) such that (Σ M ) is uniformly exponentially stable there exists a common Lyapunov function for (Σ M ) in P.

An equivalent formulation of the universality of a class P is that the converse Lyapunov theorem holds true within P.

As mentioned in the introduction, the construction of a common Lyapunov function for switched linear systems can be easily obtained. For instance, a locally Lipschitz continuous Lyapunov function may be defined as

V (x) = sup A∈S arb (M ) +∞ 0 Φ A (t, 0)x dt,
and may be regularised outside the origin by convolution with a smooth function (see e.g. [START_REF] Dayawansa | A converse Lyapunov theorem for a class of dynamical systems which undergo switching[END_REF][START_REF] Mason | Common polynomial lyapunov functions for linear switched systems[END_REF] for more details). Recalling that a function V : R n → R is absolutely homogeneous of degree α if V (λx) = |λ| α V (x) for every x ∈ R n and λ ∈ R, this classical construction leads to the following result.

Proposition 2.6. Let α ≥ 1 and P be the class of convex absolutely homogeneous functions of degree α on R n that are positive and smooth on R n \ {0}. Then P is universal. Moreover, for every bounded set M ⊂ M n (R) such that (Σ M ) is uniformly exponentially stable, there exists ε > 0 and a common Lyapunov function

V ∈ P for (Σ M ) such that ∇V (x) M x ≤ -x α for every x ∈ R n \ {0} and M ∈ M .
Note that a globally smooth Lyapunov function may be constructed by classical regularisation techniques developed in a nonlinear setting (see e.g. [START_REF] Kurzweil | On the inversion of the second theorem of Lyapunov's second theorem on stability of motion[END_REF][START_REF] Lin | A smooth converse Lyapunov theorem for robust stability[END_REF]), at the price of losing homogeneity.

Similar to Proposition 2.6, the following straightforward converse Lyapunov result links the uniform stability of (Σ M ) with the existence of a nonstrict common Lyapunov function.

Proposition 2.7. Assume that M ⊂ M n (R) is bounded and (Σ M ) is uniformly stable. Then the convex and absolutely homogeneous function of degree one

V (x) = sup t≥0,A∈S arb (M ) Φ A (t, 0)x . (2)
is a nonstrict common Lyapunov function for (Σ M ).

Due to Proposition 2.6, the continuous differentiability outside the origin of the common Lyapunov function is not a restrictive assumption when checking the uniform exponential stability of a linear switched system. On the other hand, the uniform stability of (Σ M ) does not always imply the existence of a C 1 nonstrict common Lyapunov function (see e.g. [START_REF] Chitour | Geometric and asymptotic properties associated with linear switched systems[END_REF]Example 3]). Furthermore, even in case of uniform exponential stability, it may be useful to provide a criterion to ensure the existence of a common Lyapunov function in a class of non-differentiable functions, such as piecewise linear or piecewise quadratic ones. For these reasons we introduce below a criterion which generalizes Equation ( 1) and characterizes the family of (possibly nonstrict) common Lyapunov functions in a nonsmooth setting. We refer to [START_REF] Bacciotti | Stability and stabilization of discontinuous systems and nonsmooth Lyapunov functions[END_REF]Proposition1] for a similar result in the context of differential inclusions.

We need the following preliminary result, which expresses the variation of a convex function V along a trajectory in terms of the subdifferential of V . Recall that the subdifferential ∂V (x) at a point x ∈ R n is defined as

∂V (x) = {l ∈ R n | l (y -x) ≤ V (y) -V (x), ∀y ∈ R n }.
The proof of the lemma is provided for completeness. Proof. As V is convex, V is Lipschitz and the composition V •ϕ is absolutely continuous. Hence for almost every t ∈ I the derivatives of both ϕ and V • ϕ are well-defined. By definition of subdifferential, for every t, s ∈ I and l ∈ ∂V (ϕ(t)) we have

l (ϕ(s) -ϕ(t)) ≤ V (ϕ(s)) -V (ϕ(t)).
We deduce that

d dt V (ϕ(t)) = lim s→t + V (ϕ(s)) -V (ϕ(t)) s -t ≥ l lim s→t + ϕ(s) -ϕ(t) s -t = l φ(t)
holds true for almost every t ∈ I and for every l ∈ ∂V (ϕ(t)). Similarly, taking the limit as s → t -, we obtain that d dt V (ϕ(t)) ≤ l φ(t) for almost every t ∈ I and for every l ∈ ∂V (ϕ(t)). This concludes the proof of the lemma.

Here follows an adaptation to the nonsmooth setting of the characterization of common Lyapunov function.

Proposition 2.9. Let M be a bounded subset of M n (R) and V : R n → R + be a convex positive definite function. Then V is a nonstrict common Lyapunov function for (Σ M ) if and only if

l M x ≤ 0, ∀x ∈ R n \ {0}, ∀l ∈ ∂V (x), ∀M ∈ M . (3)
Moreover, if the inequality in (3) is strict then V is a common Lyapunov function for (Σ M ).

Proof. The second part of the proposition and the if implication in the first part directly follow from Lemma 2.8. We are left to show that if V is a nonstrict common Lyapunov function for (Σ M ), then the inequality (3) holds true. By contradiction, suppose that there exist x ∈ R n , l ∈ ∂V (x), and M ∈ M such that l M x > 0. By [START_REF] Rockafellar | Convex Analysis[END_REF]Theorem 25.6] one may find a differentiability point y of V such that the pair (y, ∇V (y)) is arbitrarily close to (x, l). In particular we may assume ∇V (y) M y > 0, that is, V is increasing at t = 0 along the trajectory t → e tM y, leading to a contradiction.

Positive universality results

Given a switched linear system (Σ M ), the family P identified by Proposition 2.6 is unsuitable for numerical investigation of the existence of a Lyapunov function. With this goal in mind, interesting candidate classes P are those parametric families of functions for which the property of being positive definite and strictly decreasing along all admissible dynamics can be translated into numerically verifiable algebraic relations or inequalities (e.g., linear matrix inequalities). It is well-known that piecewise-quadratic, polynomial, and polyhedral functions represent examples of such families [START_REF] Molchanov | Lyapunov functions that define necessary and sufficient conditions for absolute stability of nonlinear nonstationary control systems[END_REF][START_REF] Molchanov | Lyapunov functions that define necessary and sufficient conditions for absolute stability of nonlinear nonstationary control systems[END_REF][START_REF] Molchanov | Criteria of asymptotic stability of differential and difference inclusions encountered in control theory[END_REF][START_REF] Blanchini | A new class of universal Lyapunov functions for the control of uncertain linear systems[END_REF]. We next provide a general sufficient condition for a class P to be universal. Roughly speaking, we exploit the fact, specific to convex functions defined on compact sets, that being close in the uniform norm is equivalent to possessing "close" subdifferentials. Proposition 3.1. Let P be a subset of the family of convex absolutely homogeneous functions of degree one from R n to R + . Assume that for every convex absolutely homogeneous function V : R n → R + of degree one and every δ > 0 there exists a function W in P such that W (x) -V (x) ≤ δ for every x in the unit sphere S n-1 of R n . Then P is a universal class of Lyapunov functions.

Proof. Let (Σ M ) be uniformly exponentially stable. Let V be the absolutely homogeneous of degree one common Lyapunov function provided by Proposition 2.6. In order to prove the proposition, it is enough to show that any convex absolutely homogeneous function close enough to V on S n-1 in uniform norm is itself a Lyapunov function for (Σ M ).

We proceed by contradiction: we assume that there exists a sequence of convex absolutely homogeneous functions (W k ) k∈N converging uniformly to V on S n-1 as k goes to infinity and such that each W k is not strictly decreasing along at least one trajectory of the system. In particular the derivative of W k along such a trajectory is nonnegative on a set of times of positive measure. By Lemma 2.8 and absolute homogeneity of W k , we deduce that there exist x k ∈ S n-1 and M k ∈ M such that, for every fixed l k ∈ ∂W k (x k ), one has l k M k x k ≥ 0. By compactness, we may assume that x k tends to x ∈ S n-1 as k goes to infinity. Then, by [14, Theorem 24.5], l k converges to ∇V (x), so that lim k→∞ ∇V (x) Mk x = lim k→∞ l k M k x k ≥ 0. However, it follows by the choice of V and Proposition 2.6 that ∇V (x) Mk x ≤ -1, yielding a contradiction.

Remark 3.2. By the absolute homogeneity property, the statement of Proposition 3.1 could be equivalently reformulated by fixing δ = 1.

As an application of the previous result, two classical examples of universal classes of Lyapunov functions (cf. [START_REF] Blanchini | A new class of universal Lyapunov functions for the control of uncertain linear systems[END_REF][START_REF] Molchanov | Lyapunov functions that define necessary and sufficient conditions for absolute stability of nonlinear nonstationary control systems[END_REF][START_REF] Molchanov | Lyapunov functions that define necessary and sufficient conditions for absolute stability of nonlinear nonstationary control systems[END_REF][START_REF] Molchanov | Criteria of asymptotic stability of differential and difference inclusions encountered in control theory[END_REF]) are recalled in the following corollary. Proof. Let V be a convex absolutely homogeneous function of degree one. Let (x i ) i∈N be a dense sequence in S n-1 , and l i ∈ ∂V (x i ). We consider the increasing sequence of absolutely homogeneous functions of degree one defined by

W i (x) = max j=1,...,i |l j x|.
Observe that each W i is convex and W i (x) ≤ V (x) for every x ∈ R n . Indeed

|l j x| = max{l j x, l j (-x)} = max{l j (x -x j ), l j (-x -x j )} + l j x j ≤ V (x) -V (x j ) + l j x j = V (x),
for every positive integer j, by definition of subgradient and since V (x) = V (-x) and V (x j ) = l j x j . We deduce that W i (x k ) = V (x k ) for every k ∈ N and i ≥ k, hence lim i→∞ W i (x k ) = V (x k ) and we can apply [14, Theorem 10.8] to conclude that the sequence of functions W i converges to V uniformly on S n-1 . By applying Proposition 3.1 we get that the family of polyhedral functions is a universal class of Lyapunov functions.

Let us now consider the absolutely homogeneous functions of degree one

Z i (x) = i j=1 |l j x| 2i 1 2i
.

The function Z i is convex since it is the composition of the 2i-norm on R i , i.e., y 2i = i j=1 y 2i j 1 2i for y ∈ R i , with the linear function from R n to R i mapping x to (l 1 x, . . . , l i x) . Moreover it is immediate to see that

W i (x) ≤ Z i (x) ≤ i 1 2i W i (x)
, and in particular Z i tends to W i uniformly on S n-1 as d goes to infinity. By applying again Proposition 3.1, it follows that the family of homogeneous sums of squares is a universal class of Lyapunov functions.

Remark 3.4. According to Remark 2.3, the first part of Corollary 3.3 remains valid if one replaces the piecewise linear functions max k=1,...,N |l k x| with the functions max k=1,...,N |l k x| q = max k=1,...,N |l k x| q , for any given q > 1. In particular, for q = 2, we have that the family of piecewise quadratic functions is a universal class of Lyapunov functions.

Remark 3.5. The proof of Proposition 3.1 relies on the fact that, whenever a linear switched system is uniformly exponentially stable, there exists a common Lyapunov function which is convex and homogeneous. In the classical construction, convexity and homogeneity are direct consequences of the convexity and homogeneity of the map x 0 → Φ A (t, 0)x 0 for given t ≥ 0 and A ∈ S arb (M ). Proposition 3.1, and hence Corollary 3.3, can then be extended to any class of nonlinear switched systems

ẋ(t) = f σ(t) (x(t)), σ(t) ∈ Σ,
with Σ ⊂ R m a bounded set of parameters, satisfying the following conditions:

• f σ is 1-homogeneous for every σ ∈ Σ, that is f σ (λx) = λf σ (x) for every x ∈ R n and λ ∈ R, • for every R > 0, {f σ | B(0,R) | σ ∈ Σ} is a compact subset of C(B(0, R), R n ),
• the class of admissible switching laws is L ∞ (R + , Σ),

• denoting by x(t, x 0 , σ(•)) the solution at time t of the system starting at x 0 and corresponding to the switching law σ ∈ L ∞ (R + , Σ), the function

x 0 → x(t, x 0 , σ(•)) is convex.

Negative universality results

Next, we provide restrictions on the classes of functions which may be candidate to be universal. For this purpose, we introduce the following technical result.

Lemma 4.1. Let M 1 , M 2 be bounded subsets of M n (R) and assume that (Σ M 1 ) is uniformly stable. For ν > 0, denote by M ν 2 the set of matrices of the form M -νId n for M ∈ M 2 , where Id n is the n × n identity matrix. Set

M = M 1 ∪ M ν 2 .
Then, the switched system (Σ M ) is uniformly stable for ν > 0 large enough.

Proof. Proposition 2.7 guarantees the existence of a convex nonstrict Lyapunov function V for (Σ M 1 ), absolutely homogeneous of degree one. Intuitively speaking, the lemma follows from the fact that, for λ large enough, the vectors M x, with M ∈ M λ 2 and x ∈ R n , point towards the interior of the sublevel set V -1 ([0, V (x)]). Let us formalize this idea. Since l x = V (x) whenever l ∈ ∂V (x), by boundedness of M 2 and of ∪ x∈V -1 (1) ∂V (x) [14, Theorem 24.7], for λ > 0 large enough one has l (M -λId n )x = l M x -λl x = l M x -λ < 0 for every M ∈ M 2 , x ∈ V -1 (1), and l ∈ ∂V (x). The result is then an immediate consequence of Proposition 2.9. Theorem 4.2. Let n ≥ 2 and P be a compact subset of the space of analytic functions from R n to R endowed with the topology of uniform convergence on bounded sets. Assume that P does not contain the zero function. Then P cannot be a universal class of Lyapunov functions.

Proof. We start by showing the theorem in the case n = 2. We proceed by contradiction, assuming that every uniformly exponentially stable switched system in the case where M consists of two matrices in M 2 (R) admits a Lyapunov function in P. We consider a switched system corresponding to

M 0 = {M 0 1 , M 0 2 } ⊂ M 2 (R)
, where M 0 1 , M 0 2 are Hurwitz, the corresponding trajectories rotate clockwise around the origin, the system is uniformly stable, but not attractive, and starting from every initial nonzero condition there exists a unique periodic trajectory, with four switches per period. The existence of such a system is obtained in [2, Theorem 1], where it corresponds to the case S4, R = 1. In particular, there exist t 1 , t 2 > 0 such that e t 1 M 0 1 e t 2 M 0 2 has an eigenvalue equal to -1, corresponding to an eigenvector x 0 . Set T = t 1 + t 2 and consider the switched systems associated with

M ε = {M ε 1 , M ε 2 } ⊂ M 2 (R), where M ε i = M 0 i -εId 2 for i = 1, 2. For ε ≥ 0 we consider the T -periodic switching sequence A ε (•) which takes values M ε 1 for t ∈ [0, t 1 ) and M ε 2 for t ∈ [t 1 , T ). Since λ(M ε ) = -ε, system (Σ M ε
) is uniformly exponentially stable for ε > 0. Hence, by assumption, it admits a Lyapunov function V ε (•) in P. Since the latter is compact, there exists a sequence (ε k ) k∈N converging to zero such that (V ε k ) k∈N converges to some V ∈ P. Moreover, for every t ≥ 0,

V (Φ A 0 (t, 0)x 0 ) = lim k→∞ V ε k (Φ A ε k (t, 0)x 0 ) ≤ lim k→∞ V ε k (x 0 ) = V (x 0 ).
Since V (Φ A 0 (2T, 0)x 0 ) = V (x 0 ) we deduce that V is constant along the trajectory Φ A 0 (•, 0)x 0 . The function t → V (e tM 0 1 x 0 ) is analytic, being the composition of analytic functions, and it is constantly equal to V (x 0 ) for t ∈ [0, t 1 ]. By analyticity, t → V (e tM 0 1 x 0 ) is constant for all t > 0 and therefore it must be identically equal to 0 since V (lim t→∞ e tM 0 1 x 0 ) = V (0) = 0. Since every nonzero point of R 2 may be written as µe tM 0 1 x 0 for some positive µ and t, we deduce that V must be identically zero, contradicting the assumptions on P.

We are left to prove the result for n > 2. For this purpose we consider

M 1 = { M 0 1 , M 0 2 } with M 0 i = M 0 i 0 0 -Id n-2 ,
where the matrices M 0 i are defined as above. Let λ > 0 and M λ 2 be given by Lemma 4.1 with M 2 = {M ∈ so(n) | M ≤ 1}, where so(n) denotes the space of skew-symmetric n × n matrices. Define M 0 = M 1 ∪ M λ 2 and, for ε > 0, consider the switched system corresponding to M ε = M 0 -εId n . It is clear that (Σ M ε ) is uniformly exponentially stable for every ε > 0.

Letting Π 1,2 be the (x 1 , x 2 ) plane, i.e., Π 1,2 = {x ∈ R n | x 3 = • • • = x n = 0}, we notice that, starting from every x ∈ Π 1,2 , there exists a periodic trajectory of (Σ M 0 ) lying on Π 1,2 . The restrictions of functions in P to Π 1,2 form a compact set of analytic functions on this plane. As in the case n = 2, we prove by contradiction that P is not universal. Assume that there exists a sequence (V ε k (•)) k∈N of Lyapunov functions in P for (Σ M ε k ) converging to V ∈ P. We can show as before that V is equal to 0 on Π 1,2 . Because of the choice of M 2 and by construction of M 0 , every 1-dimensional linear subspace of R n may be reached in finite time from Π 1,2 via a trajectory of (Σ M 0 ). Since V is non-increasing along such a trajectory, we deduce that V ≡ 0 on R n , obtaining a contradiction.

Remark 4.3. The assumption that the zero function is not in P cannot be removed from the hypotheses of Theorem 4.2. Indeed, consider the subset of polynomial functions made of the zero polynomial and, for every N ≥ 0, the polynomials of degree N whose maximum of the absolute value of the coefficients is upper bounded by a positive constant c N , where the latter is chosen in such a way that the supremum on the ball of radius N of the polynomial is less than or equal to 1/(N + 1). Since the class P contains a multiple of any polynomial, it is universal by Corollary 3.3. It is also compact since, given a sequence in P, it has a subsequence with either degree going to infinity or constant degree. In the former case, the subsequence converges to zero for the topology of uniform convergence on bounded sets, while in the latter one the coefficients are uniformly bounded and hence the sequence admits a further converging subsequence.

As a consequence of the previous result we obtain a partial counterpart to Corollary 3.3 for homogeneous polynomial functions. Namely, we recover that, if we impose a uniform bound on the degree, such functions do not form a universal class of Lyapunov functions, as already established in [START_REF] Mason | Common polynomial lyapunov functions for linear switched systems[END_REF]. Corollary 4.4. For every n ≥ 2 and every positive integer m the set of polynomial functions of degree at most m from R n to R is not a universal class of Lyapunov functions.

Proof. Without loss of generality we normalize the polynomial functions assuming (for instance) that the maximum of the modulus of their coefficients is equal to 1. In this case the coefficients of the polynomials belong to a compact subset of an Euclidean space and, since polynomials depend continuously (in the topology of uniform convergence on bounded sets) on such coefficients, we are reduced to a compact family of nonzero polynomials. The corollary then follows as an immediate consequence of the previous theorem. most m lines through the origin. Hence, it follows that either Q ε jk ≡ 0 (i.e., P ε j ≡ P ε k ) or Q ε jk vanishes at most 2d times on the unit circle. Moreover, for every ε > 0, the integer N = 2d l-1 2 + 1 is a strict upper bound for the total number of zeroes of Q ε jk for j, k ∈ {1, . . . , l}. Partitioning the circle into N arcs C 1 , . . . , C N of equal length, for every ε > 0 there exists an arc C nε which contains no zero of the nontrivial polynomials Q ε jk in its interior. Denote by A nε the closed middle third of C nε .

We next claim that for every ε > 0 there exists ν ε > 0 large enough such that the restriction of V ε to the dilated arc ν ε A nε coincides with the restriction to the same arc of one of the polynomials P ε 1 , . . . , P ε l . By definition of the function V ε and taking into account its continuity, it is enough to prove that, for every ε > 0 there exists ν ε > 0 large enough such that, in the interior of the arc ν ε A nε , one has for every j, k ∈ {1, . . . , l} that either P ε j ≡ P ε k or P ε j -P ε k is never vanishing. To see that, it is enough to prove that if Q ε jk does not vanish on A nε then P ε j -P ε k does not vanish on ν ε A nε for ν ε > 0 large enough independent of j, k. In that case, one has, for ν > 0 large enough, x ∈ S 1 and j, k ∈ {1, . . . , l},

(P ε j -P ε k )(νx) = ν d Q ε jk (x) + o(1) ,
where d is the positive degree of Q ε jk and o(1) is a function of x and ν tending to 0 as ν tends to infinity uniformly with respect to x ∈ S 1 and j, k ∈ {1, . . . , l}. This concludes the proof of the claim.

Since the arcs A 1 , . . . , A N do not depend on ε, there exist one of them, denoted by A, and sequences (ε m ) m∈N , (ν m ) m∈N in (0, +∞), and (V m ) m∈N in P 2 d such that lim m→∞ ε m = 0 and V εm = V m on νA for every m ∈ N and every ν ≥ ν m . Let Vm ∈ P 2 d be the homogeneous term of maximal degree of V m . Notice that

Vm (x) = lim ν→+∞ ν -dm V m (νx), ∀x ∈ R 2 ,
where d m denotes the degree of V m . Up to normalizing V εm , we may assume that the maximum of the moduli of the coefficients of the polynomial Vj is equal to 1. Thus, up to extracting a subsequence, Vm converges uniformly on compact sets to some nonzero V ∈ P 2 d . Similarly to the proof of Theorem 4.2, we can construct a periodic trajectory t → Φ A 0 (t, 0)x starting at x in the interior of the arc A, with A 0 (•) piecewise constant taking values in M 0 . Consider the switching laws

Lemma 2 . 8 .

 28 Let V : R n → R be a convex function and ϕ : I → R n be an absolutely continuous function, with I ⊆ R an open interval. Then V • ϕ is absolutely continuous and it holds d dt V (ϕ(t)) = l φ(t), ∀l ∈ ∂V (ϕ(t)), for a.e. t ∈ I.

Corollary 3 . 3 .

 33 The family of polyhedral functions {max k=1,...,N |l k x| | l k ∈ R n , N ∈ N} and that of homogeneous sums of squares { N k=1 (l k x) 2d | l k ∈ R n , d, N ∈ N} are universal classes of Lyapunov functions.

This research was partially supported by the iCODE institute, research project of the Idex Paris-Saclay.

The conclusion of Corollary 4.4 can be proved to hold true for functions involving maxima and minima within a finite family of polynomials such as the class of polyhedral functions V of the form

with l 1 , . . . , l N ∈ R n , where N is fixed. This partial counterpart to Corollary 3.3 is a consequence of the following more general result. Theorem 4.5. Let n ≥ 2 and P n d be the family of polynomial functions in R n of degree at most d and l be a positive integer. Consider the family

Then, P n d,l is not universal.

Proof. We first claim that if P n d,l is universal, the same is true for

If V ∈ P n d,l is a common Lyapunov function for (Σ M ), then V : R 2 x → V (x, 0) is a common Lyapunov function for (Σ M ) and V ∈ P 2 d,l . We are left to prove that P 2 d,l is not universal. Consider the switched systems (Σ M ε ) introduced in the proof of Theorem 4.2, which are uniformly exponentially stable for ε > 0, and only uniformly stable for ε = 0. Assume by contradiction that P 2 d,l is universal and, in particular, that for every ε > 0 there exists a Lyapunov function V ε ∈ P 2 d,l for (Σ M ε ). By definition of P 2 d,l , for every ε > 0 there exist l polynomials P ε 1 , . . . , P ε l of degree at most d such that V ε (x) ∈ {P ε 1 (x), . . . , P ε l (x)}. Given j, k ∈ {1, . . . , l}, we investigate the set of zeroes of the polynomial Q ε jk defined as the homogeneous polynomial corresponding to the terms of maximal degree of P ε j -P ε k . For this purpose, recall that every homogeneous polynomial Q of positive degree m may be factorized as

, where α k , β k ∈ C for k = 1, . . . , m, so that its zeroes correspond to the intersection of the unit circle S 1 with at A ε (•) = A 0 (•) -εId 2 taking values in M ε . For every t ≥ 0 such that Φ A 0 (t, 0)x ∈ A and for every m ∈ N, we have

and therefore

We then deduce that t → V (Φ A 0 (t, 0)x 0 ) is constant on {t ≥ 0 | Φ A 0 (t, 0)x 0 ∈ A}. Moreover Φ A 0 (t, 0) = e tM 0 i for some i = 1, 2, for t small enough, and by repeating the argument in the proof of Theorem 4.2 we obtain that V ≡ 0, yielding a contradiction.


Thi Thuy

Nga Nguyen

Olivier Brun

Balakrishna J Prabhu

Learning resource allocation algorithms for cellular networks

. On numerical examples with realistic mobility patterns, we show that the learning algorithm yields an approximate yet satisfactory solution with much less computation time.

Introduction

In cellular wireless networks, a central scheduling problem is to choose one among several concurrent users (for example, mobile phones) to which the scheduler (henceforth also referred to as the base station) must send data to. This scheduling decision is taken every time-slot which is of the order of 2 ms [START_REF] Margolies | Exploiting mobility in proportional fair cellular scheduling: Measurements and algorithms[END_REF] and is based on what are called as channel conditions of the users. Roughly, the channel condition of a user determines the data rate at which the base station can communicate with this user. In wireless networks, these conditions can vary randomly on short as well as on long time scales. Also called fading and shadowing, these random variations are a consequence of the interference patterns induced by the different obstacles (building, trees, etc.) in the path of the radio waves used for wireless communications [START_REF] Tse | Fundamentals of Wireless Communication[END_REF]. However, the decision is not as easy as scheduling the user with the best channel condition as such a policy could lead to unfairness between users. Imagine a user with a direct line of sight path with the base station and another who is inside a building or an underground metro station. Quite possibly, the former user will always have a better channel which would starve the latter user of any communication.

To avoid these unfair allocations, a typical solution is to define a utility which is a usually a concave function of the throughput1 of the users and then compute the scheduling decision as the one that maximizes the sum of the utilities of the users. For example, the utility function could be the logarithm of the average data rates of the users. These solutions fall under the umbrella of the network utility maximization problem [START_REF] Yi | Stochastic network utility maximisation-a tribute to kelly's paper published in this journal a decade ago[END_REF].

In an ideal scenario, the base station will solve this utility maximization problem every 2 ms. However, there are two practical issues that make this infeasible. First, the throughput of a user depends upon future channel conditions which are unknown to the scheduler. Second, the integer scheduling problem is known to be NP-complete [START_REF] Margolies | Exploiting mobility in proportional fair cellular scheduling: Measurements and algorithms[END_REF]. Solving such an optimization problem over a time horizon of seconds (thousands of time-slots) and with hundreds of users can be unrealistic in time-slots of 2 ms (which are actually becoming shorter as the technology progresses). To overcome these practical issues, several heuristics have been proposed that are based on an estimation of the future data rates [START_REF] Margolies | Exploiting mobility in proportional fair cellular scheduling: Measurements and algorithms[END_REF][START_REF] Nguyen | An algorithm for improved proportional-fair utility for vehicular users[END_REF][START_REF] Nguyen | Joint downlink power control and channel allocation based on a partial view of future channel conditions[END_REF]. The heuristics (called STO1 and STO2) in [START_REF] Nguyen | An algorithm for improved proportional-fair utility for vehicular users[END_REF][START_REF] Nguyen | Joint downlink power control and channel allocation based on a partial view of future channel conditions[END_REF] use the estimated future data rates as an input to a relaxed version of the original problem restricted to a shorter time horizon thereby reducing the dimensionality of the problem. These heuristics are thus much faster to solve than the original problem but they still require solving rather frequently a large-scale concave optimization problem which can be time consuming. It was shown numerically in [START_REF] Nguyen | An algorithm for improved proportional-fair utility for vehicular users[END_REF][START_REF] Nguyen | Joint downlink power control and channel allocation based on a partial view of future channel conditions[END_REF] that STO1 and STO2 performed better than the one in [START_REF] Margolies | Exploiting mobility in proportional fair cellular scheduling: Measurements and algorithms[END_REF] as well as the popular PF algorithm [START_REF] Kushner | Convergence of proportional-fair sharing algorithms under general conditions[END_REF] that does not use future estimated rates. Therefore, in this paper, we shall address the problem of improving the speed of these heuristics without compromising on their superior total utility.

Contributions

We propose a machine learning based solution to speed up the operations of STO12 . The key idea is to use a Deep Feedforward Neural Network (DFNN) to approximate the output of the relaxed optimization problem in the heuristics. It will be shown on numerical examples that once the DFNN is trained, it takes much less time to generate a reasonable accurate solution compared to using the specialized Python package CVXPY [START_REF] Diamond | CVXPY: A Python-embedded modeling language for convex optimization[END_REF] that uses the solver MOSEK [START_REF] Aps | MOSEK Optimizer API for Python manual[END_REF] for solving of a concave optimization problem. We compare different DFNN architectures and different loss functions to find the most appropriate ones for our problem. We then compare the behavior of the learning algorithm with STO1 and with other existing algorithms as well. The comparison shall be done on scenarios created by SUMO [START_REF] Lopez | Microscopic traffic simulation using sumo[END_REF] which can generate realistic mobility patterns on road networks with vehicular mobility. Based on numerical results, the learning algorithm is shown to perform close to STO1 with much less computation time.

Related works

The theory of approximation for DFNN has been studied in many papers. Motivated by Komogorov's superposition theorem [START_REF] Kolmogorov | On the representation of continuous functions of many variables by superposition of continuous functions of one variable and addition[END_REF] in 1957, many approximation results have proven the approximation capabilities of feedforward neural networks for the class of continuous functions such as [START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF], [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF], [START_REF] Montanelli | Error bounds for deep relu networks using the kolmogorov-arnold superposition theorem[END_REF]. In his theorem, Komogorov proved that any continuous function can be represented as a superposition of continuous functions of one variable. In [START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF] (1989), Cybenko proved that any multivariate continuous function with support in a hypercube can be uniformly approximated by a linear finite combinations of compositions of a sigmoidal functions and a set of affine functions. This representation is in fact a feed-forward neural networks with sigmoidal activation functions. Independently with the work of Cybenko, Hornik [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF] (1989) also proved a similar result. Two years later, Hornik [START_REF] Hornik | Approximation capabilities of multilayer feedforward networks[END_REF] showed that multi-layer feed-forward neural networks with arbitrary bounded and non-constant activation function can approximate arbitrary well real-valued continuous functions on compact subsets of R n as long as sufficiently many hidden layers are available. The word "deep" in "deep learning" thus simply means many layers.

Learning an algorithm to produce an approximate algorithm in order to reduce computation time has been proposed in several recent research papers [START_REF] Gregor | Learning fast approximations of sparse coding[END_REF], [START_REF] Sun | Learning to optimize: Training deep neural networks for wireless resource management[END_REF]. In [START_REF] Gregor | Learning fast approximations of sparse coding[END_REF], the authors consider a Sparse Coding problem which is used for extracting features from raw data. The problem is that Sparse Coding is often too slow for real-time processing in several applications such as pattern recognition. The authors propose a method using a non linear, feed-forward function to learn Sparse Coding to produce an approximate algorithm with 10 times less computation.

Learning an algorithm for wireless resource management has been proposed in [START_REF] Sun | Learning to optimize: Training deep neural networks for wireless resource management[END_REF]. In that work, the authors used DFNN to learn an algorithm for the interference channel power control problem. They obtain an almost real time algorithm, since passing the input through a DFNN to get the output only requires a small number of simple operations as compared to an iterative optimization algorithm. They show that, by choosing an appropriate initialization, the initial power control algorithm performs a continuous mapping which can be efficiently learnt.

In this paper, we use DFNNs for learning a channel allocation algorithm maximizing the proportional fairness between vehicular users. The proposed method is however potentially applicable to other convex optimization problems.

Organization

In Section 2, we recall the resource allocation problem in the case of a single Base Station (BS). We also remind the reader of the STO1 algorithm and state the learning problem we address in this paper. In Section 3, we formally define the input-output relationship for the DFNN model. Numerical results are presented in Section 4. We compare the computing times of the DFNNbased prediction algorithm against those of the original algorithm in Section 5 to evaluate the reduction in computing times. Finally, in Section 6 we discuss several research directions that can be followed in future work.

Problem Formulation

Consider the following downlink discrete-time channel allocation problem for a single Base Station (BS) (see [START_REF] Kushner | Convergence of proportional-fair sharing algorithms under general conditions[END_REF] and references therein):

               maximize O(α) = K i=1 log   T j=1 α ij r ij   subject to K i=1 α ij ≤ 1, j = 1, . . . , T ; α ij ∈ {0
, 1}, j = 1, . . . , T, i = 1, . . . , K.

(I)

Here r ij ≥ 0 is the data rate of user i in time-slot j, and α ij is the corresponding allocation in this slot. The constraints impose that the BS can choose at most one user in each time-slot. The objective of the BS is to maximize the sum of the individual user utilities which are defined as the logarithm of the total throughput of the user over a time horizon of T . As mentioned in Section 1, solving this problem is not practical because the data rates r i,j become known to the scheduler only in slot j. Further, users arrive and leave and it is not possible to know in advance which users will be present in the network in the future. Hence, the algorithms proposed use either no information on the future rate (see [START_REF] Kushner | Convergence of proportional-fair sharing algorithms under general conditions[END_REF] and references therein) or an estimation of the future rates [START_REF] Margolies | Exploiting mobility in proportional fair cellular scheduling: Measurements and algorithms[END_REF][START_REF] Nguyen | An algorithm for improved proportional-fair utility for vehicular users[END_REF].

A more fine-grained scheduling problem on the downlink involves joint power control and channel allocation which allows the BS to vary the transmit power in addition to choosing how much of the channel (or bandwidth) it can allocate to the users [START_REF] Huang | Downlink scheduling and resource allocation for ofdm systems[END_REF], and is defined as:

                         maximize K i=1 log   T j=1 x ij log 1 + p ij γ ij x ij   subject to K i=1 x ij ≤ 1, ∀j; x ij ≥ 0, ∀i, j 1 T j i p ij ≤ P ; i p ij ≤ P max , ∀j. (P)
Here p ij (resp. x ij ) is the power (resp. fraction of bandwidth) allocated by the BS to user i in slot j. The parameter γ ij represents the channel condition of the user i in slot j and x ij log 1 +

p ij γ ij x ij
is the Shannon rate obtained by this user when it is allocated power p ij and fraction of bandwidth x ij . The utility of the user is again the logarithm of its total throughput and the objective of the BS is to maximize the sum utility. Note that there are constraints on the power allocation which involve both the maximum power that can be expended in a time slot as well as the average power spend over the whole horizon. For this joint power and channel allocation problem, an adapted version of STO1 was proposed in [START_REF] Nguyen | Joint downlink power control and channel allocation based on a partial view of future channel conditions[END_REF].

Remark 1 (Joint power control and channel allocation). For brevity, we shall explain the STO1 algorithm and its associated machine learning solution only for the channel allocation problem (I). A remark shall be made wherever the treatment of this problem differs from that of the joint power control and channel allocation problem (P).

The STO1 algorithm

STO1 is a sequential algorithm which computes the allocation on time-slot j based on the current and the past data rates, and the estimated future data rates. It operates on two time-scales (shown in Figure 1) in order to reduce its complexity. Define big-slots as a certain number (order of hundreds) of time-slots (or small-slots) that reflect the duration over which the distance of each user to the BS does not change much, from that the data rate of each user in each time-slot can change a certain amount but the sum of its data rate over big-slot does not change much by law of large number.

For example, while the scheduling time-slots are 2 ms in length, one big slot can be equal to 100 time-slots, and one user can move at most by a few meters in a big slot. Big slots is defined to reduce dimension of the original optimization problem as explained below. It shall be assumed that estimations are available for users' future positions at the granularity of bigslots, and that the mean of future rates based on the future positions are estimated.

Before giving a formal definition of the STO1 algorithm, we give an intuitive explanation. In each small-slot j, STO1 does two steps. In the first step, it solves the problem (I) but with several restrictions: (i) the horizon is shortened to J big-slots; (ii) the future allocations are computed only on the aggregated level of big-slots; and (iii) the integer constraints on α ij are relaxed. In the second step, the fractional allocation for the current small-slot is projected onto the set of the feasible integral allocations. The first step reduces the number of variables and hence the dimensionality of the problem as the future allocations are computed only for big-slots. Note that the second step is optional and is relevant only to problems with integral constraints.

A more formal definition is as follows. Let δ be the size of the small-slot, and let ∆ be the size of the big-slot in absolute time units. Let m = ∆/δ be the number of small slots in a big-slot (see Figure 1). Denote by rij the mean rate in slot j for user i. At each small-slot t, with a slight abuse of notation, we shall denote by ρi,0 = (m-(t mod m))+t j=t+1 rij the total rate for user i in the remaining channel allocation slots of the current big-slot τ = 0, where t mod m denotes the remainder when dividing t by m. We also define ᾱi,0 as the corresponding allocation for the current big-slot τ = 0.

Denote by ρiτ = τ m+A j=(τ -1)m+A+1 rij where A = t m + 1 m, is the total average data rate that user i will get in the future big-slot τ (τ = 1, 2, ..., J -1), where big slot τ starts after the current big-slot, and J is the short time horizon in term of big slots over which we can estimate the mean future rate. We also define ᾱiτ as the corresponding allocation for user i in future big slot τ . These allocations ᾱiτ can be interpreted as the fraction of small slots that user i will be allocated in the big-slot τ .

Note that this definition is slightly different from the definition in [START_REF] Nguyen | An algorithm for improved proportional-fair utility for vehicular users[END_REF]. The differences are as follows: in [START_REF] Nguyen | An algorithm for improved proportional-fair utility for vehicular users[END_REF], there is no current big slot and the future big-slot starts just after the current small slots, but it does not change much. The above definition of two time slot types corresponds in fact to the ones introduced in [START_REF] Nguyen | Joint downlink power control and channel allocation based on a partial view of future channel conditions[END_REF].

Denote by a i (t) = t j=1 α ij r ij the total throughput allocated to user i up to time slot t, and let K(t) be the number of users inside the coverage range of the BS at time t.

The algorithm STO1 contains two steps which are as follows:

• Step 1solve the following optimization problem over a short-term horizon of J big-slots:

                   maximize K(t) i=1 log a i (t -1) + α it r it + J τ =1
ᾱiτ ρiτ subject to

K(t) i=1 α it = 1, K(t) i=1 ᾱiτ = 1, τ = 0, . . . , J -1, α it , ᾱiτ ∈ [0, 1], , τ = 0, . . . , J -1, i = 1, . . . , K(t).
(STO1-Opt) The decision variables in Problem (STO1-Opt) are the channel allocations in the current small slot, α it , and the channel allocations in the current and future big-slots, ᾱiτ . Since the future allocations are only computed on the time-scale of big-slots, there is a reduction by factor m in the number of variables in (STO1-Opt).

• Step 2obtain a feasible integral allocation from the fractional one.

For example, set α i * j = 1 with i * = argmax i α ij and set it to 0 for the other users.

Remark 2. (STO1-Opt) is solved thanks to the python package CVXPY [START_REF] Diamond | CVXPY: A Python-embedded modeling language for convex optimization[END_REF] and the solver MOSEK [START_REF] Aps | MOSEK Optimizer API for Python manual[END_REF]. In [START_REF] Nguyen | An algorithm for improved proportional-fair utility for vehicular users[END_REF], this optimization problem was solved using a projected gradient algorithm, since it allows to iteratively solve a convex optimization problem when the feasible set is a simplex or a Cartesian product of simplices, no matter how complex the objective function is as long as it is smooth and convex. The advantage of CVXPY [START_REF] Diamond | CVXPY: A Python-embedded modeling language for convex optimization[END_REF] is that it can be used to generalize this idea to other convex optimization problems, with more complex constraints such as power control constraints or others QoS constraints (e.g. delay constraints).

Learning STO1 with DFNN

As presented in [START_REF] Nguyen | An algorithm for improved proportional-fair utility for vehicular users[END_REF] and [START_REF] Nguyen | Joint downlink power control and channel allocation based on a partial view of future channel conditions[END_REF], STO1 and STO2 are better than the other existing algorithms. However they have to solve an optimization problem with a large number of variables and constraints frequently (every small slot for STO1 or every big slot for STO2), so even if their performance is good, their computations are heavy. When the system is large and requires many more QoS contraints, they may not be able to run in real time. In addition, even when they are able to run in real time, it is good to reduce the computation time without reducing too much the quality of the allocation. Therefore in this paper, our objective is to learn the STO1 algorithm using Deep Feedforward Neural Networks to obtain a new algorithm that behaves like STO1 but with a significantly reduced computation time. In other words, we want to learn the input-output relationship of STO1, by approximating the input-output mapping of STO1 with a DFNN (see [START_REF] Goodfellow | Deep Learning[END_REF] for background material on supervised learning with DFNN). After getting the approximation function (that is, the DFNN), the output can be computed by feeding the DFNN with the input value, instead of solving an optimization problem. This simpler method is expected to work faster than the original algorithm.

The main idea is to train a DFNN to predict an approximate solution to (STO1-Opt) instead of using a specialized convex optimization package. Obviously, the same idea could be used for other problems in order to obtain an approximate method which performs almost as well as the original algorithm but requires much less computing time. In short, the approach is as follows:

1. Design a well-performing algorithm based on the best available information;

2. Learn it by an approximate algorithm which behaves closely to the original one and has less computation.

The basic idea can be summarized as follows. Supervised learning is a learning task that amounts to learning an input-output relationship from a bunch of examples of input-output pairs which is called training data. The relationship can be derived from those examples by analyzing the training data and an inferred function will then be produced. A learning algorithm is considered as good if it can generalize the input-output relationship, i.e, it is able to determine outputs of unseen inputs with small error. So the objective here is to find the function representing as well as possible the input-output relationship.

Mathematically, let F be the true function that represents the relation between input and output. However, if the true function is unknown, our task is to find an approximate function F that is inferred from examples taken from training data (x n , y n ) N n=1 . The STO1 algorithm can be seen as a function F that maps an input x n ∈ X (a problem instance) to an output y n ∈ Y (a channel allocation), where X and Y denote the input and output spaces of STO1, respectively. Unfortunately, STO1 is too complicated to get the exact formula of F . Therefore here we want to approximate it by another function F : X → Y which is in the form of a DFNN (an example of DFNN is illustrated in Figure . 2). Our objective is thus to find F such that it minimizes the empirical risk

R erm (F ) = 1 N N n=1 l(y n , F (x n )),
where l(•) is a loss function which measures how far y n is from F (x n ). Note that l(•) is not necessarily a mathematical distance but it is similar to a distance in the sense that it is small if the difference between y n and f (x n ) is small.

Output: Allocation

Input: state of the system A DFNN is composed of many linear functions (sum of matrix multiplications and bias vectors) and non-linear functions (relu, sigmoid, softmax, etc.), and inside linear functions there are many parameters. So finding a good DFNN function means finding a good architecture, that is: the way the linear and non-linear functions are combined, the linear and non-linear functions in each layer, the size (number of units in each hidden layer) and then their parameters. Finding a good architecture is in general not an easy task [START_REF] Sun | On the depth of deep neural networks: A theoretical view[END_REF]. In this paper, we shall empirically compare some architectures through experiments presented in Section 4.2. After fixing the architecture, we have to find appropriate parameters by minimizing the empirical risk defined above.

System Setup for learning

Recall that, in STO1, we have two types of time scales: one is the big slot ∆ and the other one is the small slot δ. In STO1, we solve problem (STO1-Opt) with variables of size (1 + J) * K every small slot, where J is the time horizon in terms of big slots and K is the number of users in the system. The size of the allocation vector for each user is equal to 1 + J since it contains the allocation for the current small slot, α it , and the average allocation ᾱiτ for the subsequent J big slots (including the current big slot). The input of STO1 is a data rate vector of size (1 + J) * K and the total allocated throughput for the K users. The output of STO1 is the current allocation vector (α it ) i=1,...,K which is of size K, since we shall only use current allocation for making decision.

As it is defined at the moment, STO1 is not well suited to be modelled as a learning problem for the two reasons stated below.

Firstly, since K can vary over time, the dimension of the input vector will also vary. To circumvent this problem and to properly define STO1 as a function, we have to fix the size of the state. To do that, we extend the real state of the system by adding some pseudo users. Let us assume that there are at most K M users inside the system. We will then add K M -K pseudo users, where K is the number of real users in the system at time t. We will actually learn an extended version of STO1 which is STO1 when we restrict it to K users. There are many ways to extend STO1, but here we try to define an extended version that preserves as much as possible the continuity of STO1. When we mention "learning STO1", it means "learning the extended function" of STO1.

Secondly, the output of STO1 as defined above is the solution of an optimization problem. So in fact STO1 is a set-valued mapping since the solution need not be unique. But by using the CVXPY package to solve the convex optimization problem (STO1-Opt), we agree with the way it determines one of the solutions. This makes STO1 becomes a function (instead of set-valued mapping).

Remark 3 (Joint power control and channel allocation). For the joint power control and channel allocation problem, the state needs to be augmented by the remaining total power. The output of the DFNN will now give the transmit power to each user as well as the fraction of the channel it gets allocated.

State

We define a state as a matrix of size (2+J)×K M , where K M is the maximum number of users in the system. There are thus 2 + J rows, and each row has K M elements. The interpretation is as follows:

• The first row gives the current rates of the K users. We fill in the K positions on the left hand side with these current rates, and the remaining K M -K positions are filled with -1.

• The next J rows (from 2, .., J + 1) give the average rates of the users in the next J big slots. For pseudo users (the K M -K columns on the right hand side), we use the value (-1) • (∆/δ -(t mod ∆/δ)) for the current big slot and the value (-1) • ∆/δ for the other big slots.

• The last row gives the total allocated throughput of the K users. For pseudo users, we use a large enough value which is significantly greater than the total allocated throughput of real users.

By observing how STO1 works, we remark, as expected, that STO1 gives priority to users with a low allocated throughput and a high current rate. Therefore, the way we define the state (that is, by using negative values for the current and future rates of pseudo users, and extremely large values for their allocated throughput) is intended to help the model ignore quickly the pseudo users. Remark 4. Remark that there are K real users in the system at present time. Therefore, the K places of the real users in the first row which give the current rates of those users have to be strictly positive. The future rates (from the second row to the (J + 1)-th row) can be zero.

Target

We remind the reader that we want to learn only the current allocation, not the future allocation. Therefore, the target will be a vector of size K M , where the first K positions represent the fractional allocation α it of the K users as computed by STO1, and the last positions are filled with zero. Since in the optimization problem, the sum of allocation should be equal to 1, when there is no user in the system (all positions correspond to pseudo users), the allocation vector will be set to (1/K M , ..., 1/K M ) by convention.

Figure 3 illustrates the input and output of the DFNN model as described above.

Loss, DFNN architecture, initial parameters and optimizer

We will try several different loss functions and architectures and compare them in the numerical section. The initial parameters (weights) of the DFNN will be chosen as proposed in [START_REF] Lecun | Efficient backprop[END_REF], which allows the initial parameters to be not too big and not too small. The optimizer is Adam, which was first introduced in [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] and is a stochastic first-order gradient-based algorithm.

The convergence of Adam is proven in [START_REF] Reddi | On the convergence of adam and beyond[END_REF].

Numerical Comparisons

In this section, we do simulations to evaluate the influence of many factors on the behavior of the DFNN model (loss functions, architecture of the DFNN). We use the keras library [START_REF] Chollet | [END_REF] to implement our code.

There are actually a lot of factors that can have an impact on the behavior of the learning procedure such as the initial learning rate, the learning rate decay, the optimizer, the initial weight, the number of parameters, the activation functions in layers... Here we are not able to justify all our choices, but we focus on the factors which have the most significant impact on the learning algorithm in our opinion. The initial learning rate is chosen equal to 0.0015 and after each epoch, this learning rate decays by a factor 0.998.

An Unified Data Generator for Comparison

To support the comparisons in this section, the data (both for training and validation) is generated as follows. The number of users is generated randomly from 0 to K M = 10. The sojourn time of each user is generated in (0, 400) seconds. This value could of course be increased, but here in order to reduce the learning time and be able to make many comparisons, we consider only small scenarios. At each learning epoch, the model will go through 1600 samples (that is, input-output pairs). The transmission rate in each small slot is generated randomly between 0 and 5 * δ/∆. The rate we use for evaluating in SUMO scenarios are given by

r(x) = η 1 + κ e -d(x,BS)/σ , (1)
where d(x, BS) is the distance from position x to the BS, and η represents the noise level. For the SUMO scenarios in Section 4.4, we use κ = 3, σ = 100 and η ∼ Uniform(0.7, 1.3). The others parameters are equal to J = 10, ∆ = 1 s, δ = 2 ms.

Comparison of different DFNN architectures

In this part, we will consider 4 different architectures of the DFNN model and compare their performances. For the 4 models, the activation function used in hidden layers is the relu function, whereas the output layer uses the softmax function since we want the sum of the allocations to be equal to 1.

In this comparison, we use the same loss function for all models, the huber loss [START_REF]Huber loss[END_REF].

Model 1

The first model used in this section contains 2 layers which are 1 hidden layer and 1 output layer. The hidden layer contains 500 units, and in total the model has 67, 510 parameters. The architecture of this model is illustrated in Figure 4.

Model 2

As the first model, the second model contains 2 layers: 1 hidden layer and 1 output layer. However, the hidden layer contains 1000 units, and in total the model has 135, 010 parameters. We take the same number of layers as in model 1 (but with more units in hidden layers) in order to compare whether it is better to have more parameters.

Model 3

As the two previous models, the third model contains 2 layers (1 hidden layer and 1 output layer). The hidden layer contains 100 units, and we have 13, 510 parameters in total. We take the same number of layers as in model 1 (but fewer units in the hidden layer) to compare whether it is better to have fewer parameters.

Model 4

The last model contains 10 layers which are 9 hidden layers and 1 output layer. Each hidden layer contains 82 units, and in total the model contains 67, 496 parameters. We take a model that has almost the same number of parameters as Model 1, to compare whether it is better to have more layers or fewer layers.

Remark 5 (Joint power control and channel allocation). For the joint power control and channel allocation problem, we still compare the four above models except that the output layer of each model will be modified since it includes not only the channel allocation but also the power.

Figure 5a illustrates the loss of the 4 models on training and validation data. Figure 5b plots loss and absolute error of the 4 models on the same axis on validation set. The same quantities but for the problem of joint power control and channel allocation are shown in 6.

From these figures, we observe that for the model without power control:

• Having almost the same number of parameters, Model 1 with fewer layers is better than Model 4. • Having the same layers, Model 1 and Model 2 with more parameters are better than Model 3.

• Model 1 and Model 2 behave similarly and have the same number of layers. However Model 1 has less parameters than Model 2 so it is less costly from a computational point of view. Therefore from now on we shall use Model 1 for other comparisons in the sequel.

For the model with power control:

• Having almost the same number of parameters, Model 1 with few layers is slightly better than Model 4, but the difference is quite small in this case.

• Having the same layers, models 1, 2 and 3 are almost the same but Model 3 has fewer parameters.

Comparisons of different loss functions

To compare the quality of the learning model obtained using different loss functions, we use the same model, that is Model 1. Figure 7 presents the results obtained with the huber loss [START_REF]Huber loss[END_REF], with the sum of binary cross-entropy [START_REF]Cross entropy[END_REF] and dice loss (which equals to 1-dice coefficient [START_REF]Sorensen dice coefficient[END_REF]). The second loss function is denoted by bce dice loss. From the figures, the bce dice loss function is better in this case. So for the next comparisons, we shall use Model 1 and bce dice loss.

Performance Evaluation on SUMO scenarios

In this section, we shall use a mobility simulation software for comparisons of the algorithms, that is Simulation of Urban Mobility application (SUMO) [START_REF] Lopez | Microscopic traffic simulation using sumo[END_REF]. SUMO is an open source software designed for simulating mobility of moving users (vehicles, bus, truck, bicycle, pedestrian, ...) in large road traffic networks. It allows to import maps of different cities and simulate realistic mobility traces. This application is used to simulate the complex moving dynamic systems in several specific regions of Toulouse city to compare our heuristics against existing algorithms in realistic scenarios. The performance evaluation of the heuristic is done in two steps: firstly, SUMO is used for generating the mobility traces of vehicles; finally, these traces are then fed to a Python script which implements the algorithms and computes the value of the objective function of those algorithms. We shall compare the learning-based allocation scheme with STO1 and other existing algorithms on two different scenarios created with SUMO. The first scenario contains 244 users and lasts 61.7 minutes. The map of this scenario is shown in Figure 8. The results obtained with the different learning schemes on this scenario are shown in Figure 10. The second scenario contains 214 users and lasts 62.4 minutes. The map of this scenario is shown in Figure 9. The data for BS location can be found on the website3 of the French Frequency Agency (ANFR), which manages all radio frequencies in France. The results obtained with the different learning schemes on this scenario are shown in Figure 11. We also simulate two existing algorithms, (PS) 2 S [START_REF] Margolies | Exploiting mobility in proportional fair cellular scheduling: Measurements and algorithms[END_REF] and PF [START_REF] Kelly | Charging and rate control for elastic traffic[END_REF] (which are also used in [START_REF] Nguyen | An algorithm for improved proportional-fair utility for vehicular users[END_REF] for comparisons) in order to show that the approximation algorithm performs better than the existing algorithms. As mentioned above, for the learning algorithm, we use Model 1 and bce dice loss.

When the number of learning epochs is large enough, the learning-based scheme performs well compared to STO1 and other algorithms.

Computing times

The computing time of STO1 depends on the convex optimization solver used, whereas the learning algorithm has only to feed the DFNN model with the input matrix. We consider the same setting as in Section 4.2, that is K = 10 (there are 10 users in the system) and the short term horizon is J = 10 seconds. For these values, the average computing time of Mosek is around 43.7 ms, whereas the prediction with the DFNN model (Model 1, which contains 67,510 parameters) takes only 0.65 ms on average. When adding power control, Mosek solves the optimization problem in around 113.4 ms, while the prediction of the DFNN model (almost the same with Model 1 but the output layer contains power vector in addition, which contains 73,020 parameters) takes 0.68 ms on average. These computing times are averaged over 10000 samples, all are measured on a machine using GPU (graphics processing unit) which allows computing many calculations in parallel.

From the above measurements, we can conclude that the computing time with a solver can vary widely with the number of QoS constraints for the same network (number of users, time horizon). In contrast, the prediction time of the learning-based algorithm with DFNN is almost insensitive to such changes.

Summary and Discussion

We have proposed to use DFNN for learning the channel allocation obtained with one of the heuristics (STO1) introduced in [START_REF] Nguyen | An algorithm for improved proportional-fair utility for vehicular users[END_REF] and [START_REF] Nguyen | Joint downlink power control and channel allocation based on a partial view of future channel conditions[END_REF]. Numerical results on SUMO scenarios show that the learning-based method yields approximate yet satisfactory channel allocations with much less computation time as long as there are enough learning epochs. The state of the DFNN is defined in such a way that the model is not restricted to a particular scenario, that is, it can learn the channel allocation for a general network.

There are several directions of research that can be investigated to improve the learning algorithm, such as a better generator of data, a better loss function, a better architecture of the DFNN model, and other things such as the optimizer, the learning rate, etc.

Figure 1 :

 1 Figure 1: The different types of time slots in the STO1 algorithm.

Figure 2 :

 2 Figure 2: An example of Deep Feed-forward Neural Network.

Figure 3 :

 3 Figure 3: Input and output of the DFNN model.

Figure 4 :

 4 Figure 4: Model 1 architecture.

( a )

 a Loss on training set and validation set of each model (b) Plot on same axis for loss and absolute error on validation set of all the four models

Figure 5 :

 5 Figure 5: Comparison of the 4 DFNN models. 14

Figure 6 :

 6 Figure 6: Comparison of the 4 DFNN models for the joint power control and channel allocation problem.

( a )

 a Absolute error on training and validation set of the two losses (b) Plot on same axis of the absolute error on validation of the two losses

Figure 7 :

 7 Figure 7: Comparison of Loss functions.

Figure 8 :

 8 Figure 8: The Carmes borough in Toulouse, with one BS (Free Mobile type LTE1800). The actual size is 200m × 400m.

Figure 9 :

 9 Figure 9: Duroux, one BS type LTE1800, operator SFR. The actual size is around 350 × 500m.

Figure 10 :

 10 Figure 10: Comparisons of evaluated on Carmes scenario created by SUMO.

Figure 11 :

 11 Figure 11: Comparisons of evaluated on Duroux scenario created by SUMO.

We use data rate and the throughput to denote two different but related quantities. The throughput only takes into account the data rate of the time-slots in which a user is served. It is thus no more than the total data rate of this user.

STO2 is similar to STO1 but solves the optimization problem less frequently. In this paper we shall focus on STO1 but the ideas developed here can be applied to STO2 as well.

https://data.anfr.fr/anfr/portail

Acknowledgements

This work was partially funded by a contract with Continental Digital Services France.