Mourad Baiou
email: mourad.baiou@isima.fr

Aurélien Mombelli
email: aurelien.mombelli@uca.fr

Alain Quilliot
email: alain.quilliot@isima.fr

Monitoring a Fleet of Autonomous Vehicles through A* like Algorithms and Reinforcement Learning

Introduction

Intelligent vehicles, provided with an ability to move with some level of autonomy, recently became a hot spot in the mobility field. Still, determining what can be exactly done with new generations of autonomous or semi-autonomous vehicles able to follow their own way without being physically tied to any kind of track (cable, rail,. . . ) remains an issue. Most people are doubtful about the prospect of seeing such vehicles moving without any external control inside crowded urban areas. Instead they foresee that the use of those vehicles is likely to be restricted to protected areas for specific purposes: relocation of free access vehicles inside large parking areas, rural or urban logistics inside closed areas, pick up and delivery transactions inside warehouses (see [START_REF]amazon prime air[END_REF][START_REF] Ren | Cooperation of multi-robots for disaster rescue[END_REF]), rescue or repair interventions in a context of natural disaster. This point of view raises the general challenge of monitoring a fleet of such vehicles, required to perform internal logistics tasks while safely interacting with other players: workers, machines and standard vehicles. Related decision problems are at the intersection of Robotics and are at the intersection of Robotics and Operations Research. When it comes to the management autonomous vehicle fleets, current trend is to the implementation of a 3-levels supervision architecture:

• The first level, or embedded level, is defined by the monitoring and sensing devices which are embedded inside the vehicles, with the purpose of controlling trajectories in real time and adapting them to the possible presence of obstacles: currently, most effort from the robotics community remains devoted to this embedded level (see [START_REF] Chen | Cooperative intersection management: a survey[END_REF]), which mostly involves optimal control and artificial perception techniques.

• The second one, or middle one, is in charge of small tricky areas, like for instance crossroads (see Figure 2). It sends signals and instructions to the vehicles in order to regulate their transit and avoid them to collide when they get through those areas (see [START_REF] Chen | Cooperative intersection management: a survey[END_REF][START_REF] Philippe | Probability collective algorithm applied to decentralized coordination of autonomous vehicles[END_REF]). • The third one, or global one, refers to the dynamic planning and routing of the fleet, in order to make this fleet achieve some internal logistics requests (see [START_REF] Le-Anh | A review of design and control of automated guided vehicle systems[END_REF][START_REF] Vis | Survey of research in the design and control of AGV systems[END_REF][START_REF] Koes | Heterogeneous multi-robot coordination with spatial and temporal constraints[END_REF]).

A true challenge is about the synchronization of those monitoring levels and about the control of communication processes which will allow them to interact. We deal here with the global control level, while assuming that this level is in charge of vehicle routing and scheduling decisions. At a first glance, one may think into related problem as a kind of PDP: Pick up and Delivery problem (see [START_REF] Berbeglia | Static pick up and delivery problems : a classification scheme and survey[END_REF]), since an elementary task will consists for a vehicle in moving from some origin to some destination, performing some loading or unloading transaction and keeping on. But some specific features impose new challenges:

• The time horizon for autonomous or semi-autonomous vehicles is usually short and decisions have to be taken on line, which means that decisional processes must take into account the communication infrastructure (see [START_REF] Vivaldini | Communication infrastructure in the centralized management system for intelligent warehouses[END_REF]) and the way the global supervisor can be provided, at any time, with a representation of the current state of the system and its short term evolution; • As soon as autonomous vehicles are involved, safety is at stake (see [START_REF] Ryan | Spatial risk modelling of behavioural hotspots: Risk aware paths planning for autonomous vehicles[END_REF][START_REF] Pimenta | Models and algorithms for reliability oriented DARP with autonomous vehicles[END_REF]).

The global supervisor must compute and schedule routes in such a way that not only tasks are going to be efficiently performed, but also that local and embedded supervisors will perform their job more easily.

Taking care of safety requires quantifying the risk induced by the introduction into the system of any additional vehicle. Addressing this issue means turning real time collected traffic data into risk estimators (see [START_REF] Ryan | Spatial risk modelling of behavioural hotspots: Risk aware paths planning for autonomous vehicles[END_REF][START_REF] Zhang | Modeling of workflow congestion and optimization of flow routing in a manufacturing/warehouse facility[END_REF]). We do not do it here. Instead, we focus on the way resulting estimators may be used in order to take safe routing and scheduling decisions. So we assume that, at the time when we are trying to route and schedule a given vehicle V , we are provided with a procedure which, for any arc e = (x, y) or the transit network and any time value t, can compute a rough estimation of the risk related to make V be running on e at time t. Then our goal becomes to schedule the route Γ that V is going to follow, in such a way that its arrival time is minimal and that induced risk estimation remains bounded by some threshold. For the sake of simplicity, we limit ourselves to one vehicle V and one origin/destination move (o, d). Our problem may then be view as the search for a constrained shortest path (see [START_REF] Lozano | On an exact method for the constrained shortest path problem[END_REF]). But two features make it significantly more difficult:

• We must deal with a time dependent network (see [START_REF] Franceschetti | A metaheuristic for the time dependent pollution-routing problem[END_REF][START_REF] Krumke | Relocation in carsharing systems using flows in time-expanded networks[END_REF][START_REF] Park | Time dependent optimal routing in micro-scale emergency situations[END_REF]); • The on line context keeps us from relying on a heavy machinery like those related to mathematical programming.

According to this purpose, we proceed in 3 steps:

• The first one is devoted to the setting of our SPR: Shortest Path under Risk problem and to a discussion about its structural properties. • The second step is devoted to the design of algorithms designed for a static context: almost exact algorithms which adapt well-know A* algorithm for path searching in a large state space (see [START_REF] Nilsson | Artificial Intelligence[END_REF]); local search heuristic algorithms, which estimate the quality of a given route Γ through application of a filtered dynamic programming procedure. In both case, we try several notions of arc traversal decisions, relying on respectively risk versus time, risk versus distance and distance versus time estimations. • The last step deals with the online issue. We turn above mentioned algorithms designed according to a static paradigm into reactive algorithms for on line contexts. According to this prospect, we apply statistical learning and auto-adaptative reinforcement learning techniques, in order to associate, with any current traffic patterns, ad hoc arc traversal decisions.

So the paper is organized as follows: in Section II we formally describe our model and state some structural results. In Section III we describe the global structure of an exact A* algorithm and a local search heuristic, and present the ways the notion of arc traversal decision may be implemented. In Section IV, we address the on line issue and explain how statistical learning techniques may be used in order to turn static tree search or dynamic programming algorithms into fast decision rule based algorithms. Section V is about numerical experiments.

The SPR: Shortest Path under Risk Model

We refer here to a fleet of autonomous vehicles, which evolves throughout the time inside some kind of industrial infrastructure, for instance a warehouse, with the purpose of achieving internal logistic tasks (item storing and retrieving, maintenance, inventory,. . . ). Those tasks have to be performed in a safe way under small standard costs (time, energy,. . . ). It comes that, in order to set SPRC model, we first need to formalize here the safety notion.

Transit Network and Risk Function

We suppose that our fleet of vehicles moves inside a simple almost planar transit network G = (N, A), N denoting the node set and A the arc set. This network G is likely to represent for instance a warehouse (see Figure 3), or any kind of similar industrial or rural restricted area. To any arc e = (x, y) corresponds a length Le and a maximal speed v maxe : an autonomous vehicles traversing e is not allowed to go faster than v maxe while moving along e. We denote by L * the shortest path distance induced by values Le, e ∈ A. We suppose that at time t = 0, when the global supervisor of the fleet must take a decision about a target vehicle V, he is provided with some knowledge about the routes which are followed by the other vehicles, their schedule, and the tasks that they are going to perform. This knowledge allows him to derive a risk estimation Π e (t) function whose meaning comes as follows: For any small value dt, Π e (t).dt is an estimation of the Expected Damage in case V moves at maximal speed v maxe along e between time t and time t + dt.

Obtaining functions Π e is not part of this study: it requires experimental data analysis. But, since the global supervisor must maintain those risk estimation functions all along the process, those functions must be expressed according to a simple format. So we make the assumption that any function Π e is piecewise linear (see Figure 4). As a matter of fact, if additional vehicle V moves across arc e at a speed v less than maximal speed v maxe , induced risk will decrease. We assume that, if V traverses arc e during some interval [t, t + dt] at speed v ≤ v maxe , then related Expected Damage is given by Equation 1:

Risk(v, t) = Φ v v maxe .Π e (t).dt (1)
where Φ is an increasing convex function with values in [0, 1], such that for any u, Φ(u) ≪ u, which mean that Φ(u) is significantly smaller than u. The meaning of condition Φ(u) ≪ u is that, since going slower implies for vehicle V a larger traversal time, Φ(u) ≪ u will also implies that the risk induced by the traversal of e decreases while the speed decreases, even if the traversal time increases. In the next sections, we shall set Φ(u) = u 2 .

It comes that if vehicle V moves across arc e between time T and time T + δ, according to speed function t → v(t), then related Expected Damage is given by Equation 2:

T +δ T Φ v v maxe .Π e (t). (2)
Speed Normalization: We only care here about traversal times of arcs e ∈ A, and not about their true length, in the geometric sense. So we suppose here that, for any arc e, v maxe = 1. According to this we deal with reduced speed values u ∈ [0, 1] and L e means the minimal traversal time for arc e.

Routing Strategies and the SPR Problem

Let us suppose now that origin o and destination d are given, which are both nodes of the transit network G = (N, A). A routing strategy from o to d for additional vehicle V , is going to be defined by a pair (Γ, u), where:

• Γ is a path from o to d in the network G.

• u is a speed function, which, to any time value t ≥ 0, makes correspond the reduced speed u(t) ≤ 1 of the vehicle V . Notice that if we refer to the previously described speed normalization process the true speed v(t) of V is going to depend on the arc e where V is located.

Such a routing strategy (Γ, u) being given, path Γ may be viewed in a standard way as a sequence e 1 , . . . , e n of arcs of G. If we set t 0 = 0 and denote by t i when V arrives to the end-node of e i , then values t i are completely determined by speed function t → u(t). Then we set:

• T ime(Γ, u) = t n = global duration induced by the routing strategy (Γ, u)

• Risk(Γ, u) = i ti ti-1 Φ(u(t))Π e (t)dt = global risk induced by (Γ, u).
Then the SPR: Shortest Path Under Risk comes in a natural way as follows:

{SPR: Shortest Path Under Risk: Given origin o and destination d, together with some threshold R max , compute a routing strategy (Γ, u)such that Risk(Γ, u) ≤ R max and T ime(Γ, u) is the smallest possible}.

Some Structural Results

As it is stated, SPR looks more like an optimal control problem than like a combinatorial one. But, as we are going to show now, we may impose restrictions on speed function u, which are going to make the SPR model get closer to a discrete decision model.

Proposition 1. Optimal solution (Γ, u) of SPR may be chosen in such a way that u is piecewise constant, with breakpoints related to the times t i when vehicle V arrives at the end-nodes of arcs t i , i = 1, . . . , n, and to the breakpoints of function Π e i , i = 1, . . . , n.

Proof. Let us suppose that V is moving along some arc e = e i , and that δ 1 , δ 2 are 2 consecutive breakpoints in above sense. If u(t) is not constant between δ 1 and δ 2 then we may replace u(t) by the mean value u * of function t → u(t) between δ 1 and δ 2 . Time value T ime(Γ, u) remains unchanged, while risk value Risk(Γ, u) decreases because of the convexity of function Φ. So we conclude.

Proposition 2. If optimal SPR trajectory (Γ, u) is such that u(t) ̸ = 1 at some t, then Risk(Γ, u) = R max .
Proof. Let us suppose that path Γ is a sequence e 1 , . . . , e n of arcs of G. We proceed by induction on n.

• First case: n = 1.

Let us suppose above assertion to be false. Breakpoints of e = e 1 , may be written t 0 = 0, t 1 , . . . , t Q = T ime(Γ, u), and we may set:

• q 0 = largest q such that u < 1 between t q and t q+1 ;

• u 0 = related speed; l 0 = distance covered by V at time t q0 . Let us increase u 0 by ϵ > 0, such that u 0 + ϵ ≤ 1 and that induced additional risk taken between t q0 and t q0+1 does not exceed R max -Risk(Γ, u). Then, at time t q0+1 , vehicle V covered a distance l > l 0 . If l < L e , then it keeps on at speed u = 1, and so arrives at the end of e before time t Q , without having exceeded the risk threshold R max . We conclude.

• Second case: n > 1.

Let us suppose above assertion to be false and denote by R 1 the risk taken at the end of arc e, and by t 1 related time value. Induction applied to arcs e 2 , . . . , e n , and risk threshold R max -R 1 implies that the speed of V is equal to 1 all along the arcs e 2 , . . . , e n . Let us denote by τ 0 = 0, τ 1 , . . . , τ Q the breakpoints of e 1 which are between 0 and t 1 and let us set τ Q+1 = t 1 and:

• q 0 = largest q such that u < 1 between τ q and τ q+1 ;

• u 0 = related speed; l 0 = distance covered by V at time t q0+1 . Then we increase u 0 by ϵ > 0, such that u 0 + ϵ ≤ 1 and that induced additional risk taken between τ q0 and τ q0+1 does not exceed (R max -Risk(Γ,u))

2

. While moving at speed u 0 + ϵ along e 1 , vehicle V faces 2 possibilities: either it arrives at the end of e 1 before time τ q0+1 or it may keep on moving from time τ q0+1 on along e 1 at speed u = 1. In any case, it reaches the end of e 1 at some time t 1 -β, β < 0, with an additional risk no larger than (R max -Risk(Γ,u))

2

. So, for any i = 2, . . . , n we compute speed value u i such that moving along e i at speed u i between t i-1 -β and t i-1 does not induce an additional risk more than (R max -Risk(Γ,u))

2n

. So we apply to V the following strategy: move as described above on arc e 1 and next, for any i = 2, . . . , n, move along e i at speed u i between t i-1 -β and t i-1 and next at speed 1 until the end of e i . The additional risk induced by this strategy cannot exceed (R max -Risk(Γ, u)). On another side, this strategy makes vehicle V achieve its trip strictly before time t n . We conclude. Proposition 3. Given an optimal SPR trajectory (Γ, u), with Γ = {e 1 , . . . , e n } and u satisfying Proposition 1. Let us denote by t i the arrival time at the end of arc e i . Then, for any i = 1, . . . , n, and any t in [t i-1 , t i ] such that u = u(t) < 1, the quantity Φ ′ (u(t)).Π eq (t) is independent on t, where Φ ′ (u) denotes the derivative of Φ in u.

Proof. Once again, let us denote by t i time when vehicle V arrives at the end of arc e i . For a given i, we denote by δ 1 , . . . , δ H(i) , the breakpoints of function Π ei which are inside interval ]t i-1 , t i [, by Π iq q related value of Π eq on the interval ]δ j , δ j+1 [, by u 0 , . . . , u q , . . . , u H(q) , the speed values of V when it leaves those breakpoints, and by R q the risk globally taken by V when it moves all along e q . Because of proposition 2, vector (u 0 , . . . , u H(q) ) is an optimal solution of the following convex optimization problem:

• Compute (u 0 , . . . , u H(q) ) such that q u q .(δ q+1 -δ q ) and which minimizes q Φ(u q )Π ei q (δ q+1 -δ q ).

Then, Kuhn-Tucker conditions for the optimality of differentiable convex optimization program tell us that there must exists λ ≥ 0 such that: for any q such that u q < 1, Φ ′ (u q ).Π ei q = λ. As a matter of fact, we see that λ cannot be equal to 0. We conclude.

Remark 1. In case Φ(u) = u 2 , above equality Φ ′ (u q )Π ei q = λ becomes u q Π ei q = λ 2 where u q Π ei q means the instantaneous risk per distance dR dL value at the time when V moves along e i between times δ q and δ q+1 .

A consequence: Risk Versus Distance Reformulation of the SPR Model

Remark 1 leads us to define the Risk versus Time coefficient for arc e i as the value 2Φ ′ (u q )Π ei q involved in Proposition 3. This proposition, combined with Proposition 1, allows us to significantly simplify SPR: We define a risk versus distance strategy as a pair (Γ, λ RD ) where: Let us suppose that we follow a trajectory (Γ, u) which meets Proposition 3, and that we know value λ RD e for any arc e of Γ .Since Φ is supposed to be convex and such that Φ(u) ≪ u, we may state that Φ ′ admits a reciprocal function Φ ′-1 . Then, at any time t when vehicle V is inside arc e, we are able to reconstruct value

• Γ is a path, that
u(t) : Φ ′-1 ( λ RD e 2Πe ), if Φ ′-1 ( λ RD e 2Πe ) < 1 1, otherwise (3)
According to this and Proposition 3, SPR may be rewritten as follows (we extend the notations Risk(Γ, u) and T ime(Γ, u) as Risk(Γ, λ RD ) and T ime(Γ, λ RD ):

Risk versus Distance SPR Reformulation: Compute risk versus distance strategy (Γ, λ RD ) such that Risk(Γ, λ RD ) ≤ R max and T ime(Γ, λ RD ) is the smallest possible.

Discussion about the Complexity

The time dependence of the transit network together with the proximity of the SPR model with Shortest Path Constraint models suggests that SPR is a complex problem. Practical difficulty of SPR may be captured through the following example, which makes appear that if (Γ, u) defines an optimal SPR trajectory, the risk per distance value λ RD e = 2Φ ′ (u(t))Π e may be independent on t on arc e as told in Proposition 3, but cannot be considered as independent on arc e. Path Γ contains 2 arcs, e 1 and e 2 , both with length 1 and maximal speed 2. Function Π e2 is constant and equal to 1. Function Π e1 takes value 2 for 0 ≤ t ≤ 1, and a very large value M (for instance 100) for t > 1. R max = 3 4 ; Function Φis : u → Φ(u) = u 2 .Then we see that vehicle V must go fast all along the arc e 1 , in order to get out of e 1 before this arc becomes very risky. That means that its speed is equal to 1 on e 1 , and that its risk per distance value is equal to 1 2 . Next it puts the brake, in the sense that its speed remains equal to 1 but its risk per distance value decreases to 1 4 . It is easy to check that this routing strategy is the best one, with Risk(Γ, u) = 3 4 and T ime(Γ, u) = 2. Sill, identifying the complexity of SPR is not that simple, since we are dealing with continuous variables. As a matter of fact, complexity also depends on function Φ. We conjecture that:

Conjecture 1. If Φ(u) = u 2 then SPR is in NP time and is NP-Hard.

Algorithms

Our algorithms all rely on notions of state and decision. A state is a 3-uple (i, T, R), where:

• i is a node of G where vehicle V is currently located; • T is the times spent in order to reach i, and R is the amount of risk induced by this process of moving from origin o to node i.

Then a decision will consist in:

• Choosing the arc e = (i, i o ) along which the vehicle is going to move;

• Choosing some parameter λ which is going to determine the speed function u along the arc e.

Previous section 2 suggests the use of risk versus distance coefficient λ RD e as decision parameter λ. But other choices are possible. We restrict ourselves to the case when Φ(u) = u 2 .

Decision Scheme

As told above, a natural approach is to refer to Proposition 3 and consider λ = λ RD as expressing the mean Risk versus Distance coefficient Φ ′ (u)Π e . But another intuitive approach is to consider λ = λ SP as expressing the mean speed of V along e, and deduce this way the arrival time on i o in a straightforward way. Finally, we may also consider that λ = λ RS expresses the mean Risk Speed of V along e, which means the amount of risk vehicle V takes per time unit as it advances along e. We are going to describe here those 3 possibilities, together with the way resulting state (i o , R o , T o ) may be deduce from λ and (i, R, T ).

• First approach: The Risk versus Distance approach.

Since P hi(u) = u 2 , Φ ′ (u(t))Π e (t) = 2u(t)Π e (t) for any t during the traversal of e. It comes that if we fix λ RD the speed value u(t) is given by: u(t) = Inf (1, λ RD /Π e (t)). Resulting state (i o , R o , T o )will be obtained from λ RD and (i, R, T ) through the following iterative process:

Risk_Distance Transition procedure: Let us set t 0 = T , and let us denote by t 1 , . . . , t Q the breakpoints of Π e which are larger than T and by Π e 0, . . . , Π e Q related Π e values. Initialization:

t ← t 0 ; r ← R; L ← 0; q ← 0; N ot Stop ; While N ot Stop do π ← Π e q ; q ← q + 1; δ ← t q -t; u = Inf (1, λ RD /π); If L e > L + uδ then L ← L + uδ; r ← r + Φ(u)πδ; t ← t q ; Else δ ← δ Le-L uδ ; t ← t + δ; r ← r + Φ(u)πδ; L ← L e ; Stop; R o ← r; T o ← t; If R o > R max then Fail else Success;
• Second approach: The Mean Speed approach.

Fixing λ SP means fixing the time T o as: T o = T + Le λ SP . In order to determine the function t → u(t) and the value R o , we solve the following quadratic program:

Mean_Speed Program: Let us denote by t 0 = T , t 1 , . . . , t q = T o the breakpoints of Π e which belong to [T, T o ] and by Π e 1 , . . . , Π e Q related Π e values. Then we must compute speed values u 1 , . . . , u Q ∈ [0, 1] such that: q u q (t q -t q-1 ) = T o -T q u 2 q Π e q (t q -t q-1 ) < R max . This quadratic convex program may be solved through direct application of Kuhn-Tucker 1st order formulas for local optimality. Then we get R o by setting: R o = R + q u 2 q Π e q .(t q -t q-1 ). If R o > R max then the Mean Speed transition related to λ SP yields a Fail result.

• Third approach: The Risk Speed approach.

Since Φ(u) = u 2 we have that at any time t during the traversal of e, related risk speed dR dT (t) is equal to u(t) 2 Π e (t). It comes that if we fix λ RS we get:

u(t) = Inf (1, λ RS Π e (t) 1/2 ).
Resulting state (i o , R o , T o ) will be obtained from λ RD and (i, R, T ) through the same following iterative process a for the Risk versus Distance approach. Several kinds of controls may be applied to above process: one may do a random walk descent. In any case, we need to discuss both Update and Evaluate steps.

A Local

• Update step: It relies on a pre-process which is applied to the transit network G and involves some proximity threshold S_P rox. For any two nodes i, j of G such that L * i,j ≤ S_P rox, we pre-compute a collection P ath i,j of elementary path from i to j. This provides us with an operator Detour, which acts on any path Γ though parameters i, j, γ as follows:

• i, j are nodes of Γ such that i precedes j in Γ ; γ is some path in P ath i,j ;

• Detour(Γ, i, j, γ) replaces the restriction Γ i,j of Γ from i to j by path γ. Since Detour may admit a rather large number of parameters values (i, j, γ), we first identify pairs of nodes (i, j) in Γ , such that the slowdown coefficient

Tj -Ti T * i,j
is large, and pick up such a pair (i, j). Next we choose path γ in P ath i,j under the condition that is not very crowded between time T i and time T j , that means which is such that the sum, for the arcs e of γ of mean Π e (t) values between time T i and time T j is small. • Evaluation step: This evaluation step relies on a dynamic programming procedure DP_Evaluate whose main features come as follows:

• Let us denote by e 1 , . . . , e n the arcs of Γ , and by i 0 , . . . , i n related nodes; So the time space of DP_Evaluate comes in a natural as the set {0, 1, . . . , n} and a state at time q = 0, 1, . . . , n, is a pair (T, R), where T means the time when vehicle V arrives in i q , and R the cumulative risk at this time. Clearly, initial state is (0, 0) and final state should be any pair (T, R) such that R ≤ R max . • Then a decision at time q becomes a value λ, (λ RD , λ SP , λ RS ) in the sense of Section 3.1 and such a decision induces a transition (q, R, T ) → (q +1, R o , T o ) as described in Section 3.1, with cost R o -R. This decision is feasible if it does induce a Fail result.

According to this, Bellman principle may be applied: the algorithm DP_Evaluate scans the time space {0, 1, . . . , n}, and, for any q = 0, 1, . . . , n, computes related state set State[q], according to the following instructions:

• Initialize State[0] as {(0, 0)} and State[q] as N il for any q > 0;

• For q = 1, . . . , n -1 do Generate decision set Λ; For any λ in Λ and any state

(T, R) in State[q] do (I1) Compute (in case λ is feasible) resulting state (T o , R o ); If there does not exist (t 1 , R 1 ) in State[q +1] such that t 1 ≤ T o and R 1 ≤ R o then insert (T o , R o ) into State[q + 1] and remove from State[q + 1] any (t 1 , R 1 ) such that t 1 ≥ T o and R 1 ≥ R o .
In case we are already provided with some feasible SPR solution (Γ * , u * ) with value T * , then we may apply the following filtering rule:

• Lower Bound Based Filtering Rule:

Let (T o , R o ) be the state involved in instruction (I1). If T o + L * i(q+1),d ≥ T * then state (T o , R o
) may be killed: we do not insert it into State[q + 1], since we cannot expect it to be extended into a better solution than current solution (Γ * , u * ).

Remark 2. We turn LS_SPR algorithm into a greedy algorithm by removing the update step and by generating Λ in such a way its cardinality is 1.

A A* Algorithm

A* algorithm [START_REF] Nilsson | Artificial Intelligence[END_REF] was designed in order to deal with path search for robots evolving in very large (possibly infinite) state spaces. It can be adapted to our problem, since solving SRP means searching for a shortest path in a Risk Expanded network, whose nodes are pairs (i, T, R), i ∈ N, 0 ≤ R ≤ RM ax, T ≥ 0, and arcs corresponds to transition ((i, T, R) → decisionλ → (i o , T o , R o )) as described in Section 3.1. In the present case, it will rely on the following data structures:

• An expansion list LE, which contains states (i, T, R), ordered according to increasing optimistic estimation value W . The optimistic estimation value of state (i, T, R) is equal to T + L * i,d and provides us with a lower bound to the best possible value of a SPR solution (Γ, u) which would extend the path which allowed us reaching state (i, T, R).

• A pivot list LP ivot, which contains states (i, T, R), together with optimistic estimation value W , which already appeared as the first element (Head) of LE. There should not exist in LE an element (i, T, R) which is dominated by another element (i,

T 1 , R 1 ) in LP ivot LE, that means which is such that T 1 ≤ T and R 1 ≤ R.
Then A*_SPR algorithm may be described as follows:

A*_SPR Algorithm: Initialize LP ivot as N il and LE as {(o, 0, 0)}; N ot Stop;

While (N otStop) ∧ (LE ̸ = N il) do (i, T, R) ← Head(LE); If i = d
W o = T o + L * i o ,d ; If R o ≤ RM ax and if there does not exist (i o , T 1 , R 1 ) in LP ivot LE such that such that T 1 ≤ T o and R 1 ≤ R o then Insert (i o , T o , R o ) into LE and remove from LE any (i o , T 1 , R 1 ) such that T 1 ≥ T o and R 1 ≥ R o .
Do it in such a way that LE remains ordered according to optimistic estimation values W ;

Remark 3. If we are able to generate all decisions likely to appear inside a given optimal decision sequence, then above algorithm A*_SPR is optimal.

Remark 4. We turn A*_SPR algorithm into a shortest path algorithm by reducing Λ to 1 element.

Discussion: The Decision Set Λ

Both above algorithms rely on an instruction 'Generate λ -decisionsetΛ'. But λ values are continuous ones. So, we must decide about the way we generate a finite λ -decisionsetΛ.

The simplest case is the case when we deal with Risk versus Distance decisions λ RD , since in such a case, propositions Proposition 2 and Proposition 3 suggests us that a mean value for λ RD is going to be given by λ

RD mean = Rmax L * o,d
. Then a natural way to generate Λ is to fix an odd number 2.K + 1 of λ RD values, a geometric step value δ > 0, and to set:

Λ = {λ RD mean } ∧ {(1 + δ) k .λ RD mean , k = 1, . . . , K} ∧ {(1 + δ) -k .λ RD mean , k = 1, . . . , K} (4)
According to this, Λ is determined by K and δ. We may consider K as a flexible parameter. As for the choice of value δ, it becomes determined by K and by the minimal and maximal values λ RD max = (1 + δ) k λ RD mean and λ RD min = (1 + δ) -k λ RD mean which we want to assign to λ RD . If we want to allow the vehicle to move with a speed twice as large as the speed suggested by λ RD mean , then we see that we must choose a value λ RD max = 2λ RD mean . That means that δ is determined by the acceleration coefficient ρ we may want to apply to the vehicle in order to make him possible to reach the end of an arc before some risky situation occurs on this arc. We shall test for instance ρ ∈ [START_REF] Berbeglia | Static pick up and delivery problems : a classification scheme and survey[END_REF][START_REF] Chen | Cooperative intersection management: a survey[END_REF][START_REF] Franceschetti | A metaheuristic for the time dependent pollution-routing problem[END_REF].

In the case of Risk over Time and Mean Speed decisions λ RS and λ SP , we must arbitrarily fix mean values λ SD mean and λ SO mean , and try to learn them throughout the computational process. This opens the way to next section.

Key issue here becomes about the value of ω. Intuitively, ω should be equal to Rmax T o , where T o is the optimal SPR value, and we should be able to learn this value as a function of the main characteristics of SPR instances: most relevant characteristics seem to be risk threshold R max , the length (expected length) L * of path Γ , the mean value ∆ of functions Π e , e ∈ A, and the frequency B of the breakpoints of those functions. We may notice that in case all functions Π e are constant and equal to some value ∆, then speed u is going to be constant and equal to Rmax L * ∆ , and so that time value T o will be equal to

L * u = ∆L * 2
Rmax . This will lead to initialize ω as ω =

R 3 max ∆L * 2 .
In order to refine this initial choice for ω we retrieve, for any instance, the optimal decision sequence {λ 1 , . . . , λ n }, related state sequence s 1 = (T 1 , R 1 ), . . . , s n = (T n , R n ). Then, for any such an instance, we look for the value ω which statistically makes states s i be always among the best ones for the ranking related to quantity ωT + R. As a matter of fact, while performing numerical experiments, we focus on an estimation of optimal value T o as a function of R max , L Lj is large with respect to Rmax L * , or cautious if the converse holds. Then it may happen that our pruning technique yields pairs (T, R) which, taken as a whole, are either too risky or too cautious. In order to control this kind of side-effect, we make ω value become auto-adaptive.

More precisely, we start, as previously explained, from some pre-learned ω value, and make it evolve through Reinforcement Learning, that means throughout an (or several) execution the DP_Evaluate (or A*_SPR or LS_SPR). In order to explain it better, we focus on the case of DP_Evaluate algorithm, while supposing that state threshold S has been fixed and that initial value ω has been computed as described above. So let us suppose that, at some time during the process, we just dealt with arc e i and so computed current state set State[i], while updating value ω. Applying decisions of Λ and filtering resulting states (T, R) through Bellman principle provides us with a state subset State[i + 1] whose size is likely to exceed S. Then we rank states (T, R) of State[i + 1] according to ωT + R values. Ideally, states (T, R) ordered this way should make S best states (T, R) be balanced in the sense that risky states should get along with cautious ones, or, in other words, that the ratio R Rmax should be centered around the ratio L * 0,i L * . If, for instance, those values are centered significantly above this ratio, then we are moving in a too risky way and must make ω decrease. Con-versely, if those best values are centered below this ratio, then we are too careful and must make ω increase.

We implement this principle by performing a kind of statistical analysis of those best values in State[i + 1], in order to derive, from those S best states (T, R), an indicator Risk_Balance, which takes symbolic values {Risky, N ormal, Caref ul} depending on the way the mean Outputs related to the behavior of the procedure DP_Evaluate For every instance, L * value is the length of path Γ in the L sense. We apply DP_Evaluate while testing the role of parameters λ = λ RD , λ SP , λ RS , as well as K and δ. So, for every instance, we compute:

• in Table 2: The risk value R_DP mode , the time value T _DP mode , and the CPU times (in s.) CP U mode , induced by application of DP_Evaluate with λ mode = λ RD , λ SP , λ RS , K = 10, ρ = 4;

• in Table 3: For the specific mode λ RD , related number State of states per node i, together with time value T RD , when K = 1, 3, 5, 7, 10 and ρ = 4;

• in Table 4: For the specific mode λ RD , related number State of states per node i, together with time value T RD , when K = 10 and ρ = 1.5, 2, 3, 4, 8.

CPU times are in seconds. Outputs related to the behavior of A*_SPR and LS_SPR. We test the ability A*_SPR and LS_SPR. to catch optimal solution, and observe the characteristics of resulting path. We rely on λ = λ RD , K = 10 and ρ = 4. For every instance, we compute:

• in Table 5 CPU times are in seconds.

Fig. 1 .

 1 Fig. 1. An Autonomous Vehicle

Fig. 2 .

 2 Fig. 2. Hierarchical Supervision Architecture

Fig. 3 .

 3 Fig. 3. A Warehouse like Transit Network

Fig. 4 .

 4 Fig. 4. A piecewise function Π e

  means a sequence {e 1 , . . . , e n } of arcs, which connects origin node o do destination node d; • λ RD e associates, with any arc e in Γ , Risk versus Distance coefficient λ RD e = 2Φ ′ (u)Π e . In case Φ(u) = u 2 , we notice that this coefficient means the amount of risk per distance unit induced on arc e at any time t such that u(t) < 1, by any trajectory (Γ, u) which satisfies Proposition 3.

Fig. 5 .

 5 Fig. 5. Detour Operator

:

  The time value T _A * , the risk value R_A * , CPU time (in s.) CP U A * , the number Node of visited nodes, the number State A * of generated states, and the deviation Dev A * between the length of resulting path Γ and L * , induced by A*_SPR with λ mode = λ RD , K = 10 and ρ = 4;• in Table6: The time value T L S, CPU time (in s.) CP U LS , the number Trial of trials, the number State of generated states, and the deviation Dev LS between the length of resulting path Γ and L * , which derives from applying LS_SPR with λ mode = λ RD , K = 10 and ρ = 4.

  Search Algorithm Involving Dynamic Programming. Initialize Γ as the shortest path according to L from o to d; N ot Stop; While N ot Stop do First step: Evaluate Γ , and get the arrival time T i of vehicle V in any node i of Γ ; Second step: Update Γ ; Keep the best path solution Γ ever obtained.

	LS_SPR Algorithm:
	This local search heuristic LS_SPR works in 2 steps:

  then Stop; Retrieve the SPR solution Γ related to (i, T, R); Else Remove (i, T, R) from LE and Insert it into LP ivot; Generate λ decision setΛ; (I2) For any arc e = (i, i o ) and any λ in Λ do Compute resulting state (i o , T o , R o ) together with value

  , B and ∆, and next test the ranking of states s i among the state subsets State[i] for ω = Rmax T o . Bounding States through Reinforcement Learning. Still, this way of performing learning may induce distortions. A lack of flexibility in the pruning procedure associated with a non fully well-fitted value ω may yield, for a given node i, a collection State[i] poorly balanced, in the sense that one would expect related values (T, R) to distribute themselves as a wide Pareto set. More precisely, we may qualify a pair (T, R) as risky if

	R
	j≥i+1

*

Table 1 .

 1 . Then our filter&learn Filter_Learn works as follows: Filter_Learn Procedure: Rank states (T, R) of State[i + 1] according to ωT + R values; Select S best pairs (T, R) according to this ranking and compute Risk B alance; If Risk B alance = N ormal then Keep only the S best states in State(i + 1); If Risk B alance = Risky then Split State[i + 1] into 2 subsets S 1 and S 2 with same size: S 1 is made of the best states (S, R) according to our tanking and S 2 = S -S 1 ; Keep only the S 2 best states in S 1 and S 2 in State(i+1); Make ω decrease; If Risk B alance = Cautious thenproceed as in previous case, while making ω increase. We perform numerical experiments with the purpose of studying• the behavior of static DP_Evaluate, LS_SPR, A*_SPR algorithms of Section 3. We pay special attention to the dependence of those algorithms to the choice of the decision mode (Mean Speed, Risk versus Time, Risk versus Distance), to the characteristics of decision set Λ; • the way we may efficiently turn those static algorithms into efficient dynamic algorithms through the use of statistical an reinforcement learning techniques described in Section 4. We generated networks (N, A) as connected symmetric partial grids, which means grids n * n, modified through removal of 25% of its arcs. Those partial grids are summarized through their number |N | of nodes and their number |A| of arcs. Length values L e , e ∈ A, are uniformly distributed between 3 and 10. Function Φ is taken as function u → Φ(u) = u 2 . Function Π e are generated by fixing a time horizon T max , a mean frequency B of break points t e i , and an average value ∆ for value Π e (t): More precisely, values Π e are generated within a finite set {2∆,3∆ 2 , ∆, ∆ 2 , 0}. As for threshold R max , we notice that if functions Π e are constant with value ∆ and if we follow a path Γ with length L diam , where L diam is the diameter of network G, at speed 1 2 = vmax 2 , then the expected risk is L diam ∆2. It comes that we generate R max as a quantity α L diam ∆2, where α is a number between 0.2 and 2. Finally, since an instance is also determined by origin/pair (o, d), we denote by L * the value L * o,d . Table1presents a package of 12 instances with their characteristics. Instances' characteristics

	R Rmax value is lo-

*

Adaptation ot the State Bounding Scheme to A*_SPR Algorithm. In the case of the A*_SPR algorithm, we apply the same principle, with the idea that, in the list LE ∧ LP ivot, the number of states (i, T, R) related to a given i should not exceed some target threshold S. We structure elements of LE and LP ivot according to lists LE[i] and LP ivot

[i]

, where each list LE[i] and LP ivot[i] is a list of states (T, R) in the DP_Evaluate sense and apply the learn&filter process only to the LE[i] lists.

5 Numerical Experiments

Goal:

Table 2 .

 2 Impact of λ mode , with K = 10 and ρ = 4Instance R RD T RD cpu RD R SP T SP cpu SP R RS T RS cpu RS

	1	6.74 104.0	0.48	6.90 117.7	0.73	4.58 112.6	0.47
	2 36.10	37.2	0.62 35.66 38.7	0.52 34.00 38.9	0.78
	3 85.48	41.3	0.72 77.57 44.7	0.86 81.28 41.7	1.02
	4	5.42	99.4	0.51	6.41 105.0	1.57	3.96 102.5	0.49
	5 30.07	41.5	0.97 30.18 45.0	1.99 30.18 43.7	2.03
	6 84.62	34.5	1.10 79.93 37.1	1.89 86.67 36.1	2.22
	7	7.99 362.8	4.02 21.58 390.2	5.66	6.27 370.1	5.28
	8 107.73 144.6	5.88 107.68 177.6	7.27 91.08 164.1	9.05
	9 227.87	98.3	6.61 227.57 99.6	8.47 228.39 98.2	9.26
	10 11.44 410.5	8.42 74.17 423.1 13.07	8.52 388.7 13.59
	11 106.38 143.7	10.16 107.61 167.6 20.06 88.08 157.5 19.35
	12 203.65	87.5	10.55 202.81 91.7 19.53 207.34 86.3 23.54

Table 3 .

 3 Impact of K, with ρ = 4 RD States T RD States T RD States T RD States T RD States

	K	1	3	5	7	10
	instance T 1 111 10.67	108 28.33

Table 5 .

 5 Behavior of A*_SPR, with λ mode = λ RD , K = 10 and ρ = 4

	Instance T _A * R_A *	N ode	cpu A *	States A *	ρ A *
	1	94.1	5.64	14.00	0.41	17.00	4.26
	2	38.9	35.95	5.00	0.20	14.60	0.00
	3	41.4	84.21	5.00	0.04	1.80	0.00
	4	99.5	6.39	12.00	0.62	14.42	0.00
	5	42.9	30.14	9.00	0.63	11.89	0.00
	6	36.9	82.30	5.00	0.06	1.40	0.00
	7 310.4	20.58	97.00	7.80	18.98	10.94
	8 132.3	107.55	44.00	5.53	16.09	5.08
	9	98.3	228.03	17.00	0.42	1.88	0.00
	10 349.2	24.88	96.00	14.16	18.65	1.30
	11 132.1	106.97	45.00	12.81	14.24	0.41
	12	85.5	208.41	17.00	0.81	1.06	0.00

Table 6 .

 6 Behavior of LS_SPR, with λ mode = λ RD , K = 10 and ρ = 4

	Instance T _LS R_LS	T rial	cpu LS	States LS	ρ LS
	1 102.5	5.60	2.00	0.35	31.33	5.09
	2	39.5	34.85	1.00	0.31	47.50	0.00
	3	41.4	84.07	1.00	0.35	49.17	0.00
	4	99.6	4.72	1.00	0.30	35.50	0.00
	5	41.8	29.04	2.00	0.86	47.00	0.00
	6	34.6	84.14	1.00	0.61	49.83	0.00
	7 365.0	6.28	4.00	7.48	46.32	0.00
	8 160.9	107.72	10.00	22.47	58.31	18.50
	9	99.1	227.33	1.00	3.39	63.50	0.00
	10 458.8	9.81	5.00	20.11	60.61	2.48
	11 164.3	88.26	1.00	5.43	60.67	0.00
	12	88.3	205.68	1.00	5.93	66.06	0.00

Speeding Algorithms through Statistical Learning Techniques

We consider here two ways of speeding our algorithms in order fit with a dynamic contexts. The first one impose a small number K of possible decisions and the second one is to bound the number of states (T, R) related to any node i. We do it while focusing on the case λ = λ RD .

Bounding Decisions

Once acceleration parameter ρ has been tuned, controlling the size of decision set Λ means fixing value K. If we set K = 1 (greedy algorithm in the case of DP_Evaluate and shortest path algorithm in the case of A*_SPR), then the choice is about λ RD mean , which, in a first approach, should be equal to Rmax

. If K ̸ = 1, then we apply the following statistical learning process:

• We apply DP_Evaluate to instances which fit parameter ρ, while using some reference decision number K ref .

For any instance I, we retrieve the optimal decision sequence {λ 1 , . . . , λ n }, Every decision λ i is related to some number k i in -K ref , . . . , 0, . . . , K ref .

• Then we compute, for every value k in {-K ref , . . . , 0, . . . , K ref }, the percentage τ (k) of occurrence of k in those decision sequences. • Finally, K being the target decision number, we split the decision range 4 ] and [λ 3 4 , λ max ] in such a way that:

• 1/4 of decisions λ i belong to interval [λ min , λ 1 4];

Then restricted Λ becomes the set {λ

Bounding States

In order to filter state set State[i] related to a given node i and impose a prefixed lower bound S on the cardinality of State[i], several techniques may be applied. One may for instance consider as equivalent 2 states (T, R) and

does not exceed some rounding threshold. We are not going to follow this approach which does not guaranty that we are going to maintain the cardinality of State[i] below imposed threshold S. Instead, we are going to do as if there were existing a natural conversion rate ω which turns risk into time. According to this, we are going to rank pairs (T, R) State[i] according to increasing values ωT + R and keep on with the S best ones according to this ordering, while killing the others.

Comments:

We see that A*_SPR and LS_SPR are very close to each other with λ mode = λ RD . Also, they often use the shortest path but not every time especially when working with small R max values.

Outputs related to characteristics of the solutions (Learning preprocess). We apply the learning devices of Section 4 and test the behavior of DP_Evaluate with no more than 3 possible decisions and 6 possible states for every i = 1, . . . , n. For every instance, we compute in Table 7, related risk value R_F ast_DP , the time value T _F ast_DP , and CPU times CP U _F ast_DP . CPU times are in seconds.

Conclusion

We dealt here with a shortest path problem with risk constraints, which we handled under the prospect of fast, reactive and interactive computational requirements. But in practice, a vehicle is scheduled in order to perform some kind of pick up and delivery trajectory while performing retrieval and storing tasks. It comes that a challenge becomes to adapt previously described models and algorithms to such a more general context. Also, there exist a demand from industrial players to use our models in order to estimate the best-fitted size of an AGV fleet, and the number of autonomous vehicles inside this fleet. We plan addressing those issues in the next months.


Corey Duncan
email: corey.duncan@univ-fcomte.fr

Samir Jemei

Robin Roche

Marie-Cécile Pera

Modeling of hydrogen-based power-to-gas plant operation for evaluation of control strategies

Keywords: gas, optimization, system modeling, techno-economical analysis

This paper presents a simulation model for operation of a power-to-gas plant for the purposes of project feasibility and control strategy analysis. The plant is translated into an optimization problem formulation using a modified unit commitment style technique. A pilot project in France is presented as a case study for simulation model application and results from running the simulation are presented. Using maximizing production and minimizing operational cost as two separate objective functions, the resulting output data are analysed. The results show the influence of wind power, natural gas grid availability and electricity prices on production, highlighting the importance of their inclusion in power-to-gas plant modeling, as well as the consideration of more realistic equipment modeling in feasibility studies.

MOTS-CLES -power-to

Introduction

The first continent to make such a commitment, the European Council recently released the "Fit for 55 Package" which sets an ambitious goal for the European Union (EU) to reach carbon neutrality by 2050, with an initial reduction of 55% of greenhouse gas (GHG) emissions compared to 1990 levels by 2030 [START_REF] Council | Fit for 55[END_REF]. To reach carbon neutrality by mid-century, the International Energy Agency (IEA) estimates that global electricity demand will double by 2050, with renewable energy -namely variable renewable energy sources wind and solar -covering nearly 90% of total electricity generation [START_REF] Iea | Net Zero by 2050 -A Roadmap for the Global Energy Sector[END_REF]. For maximum utilization of these intermittent energy sources at such high penetration levels, flexible and adaptable energy networks must be implemented [START_REF] Lund | Review of energy system flexibility measures to enable high levels of variable renewable electricity[END_REF].

Coupling enables energy that may be curtailed or otherwise wasted to be fed into another energy network or sector through a mediary technology, decoupling the energy source generation from its direct energy consumption [START_REF] Irena | Sector Coupling[END_REF]. Hydrogen (H 2 ) has received increased attention in the past years as the mediary energy carrier as it allows for several sectors and energy networks to be coupled. This process is known as Power-to-Gas or PtG. PtG technologies also involve converting H 2 into methane -also known as synthetic natural gas or SNG when produced in this fashion -via methanation reactors for increased technology applicability and flexibility. Along with the IEA hydrogen tracking report [START_REF] Iea | Hydrogen -Analysis[END_REF] and their Net Zero study [START_REF] Iea | Net Zero by 2050 -A Roadmap for the Global Energy Sector[END_REF], many studies have been done that show the economical and energy efficiency benefits of including such technologies in future energy strategies [START_REF] Henni | Infrastructural coupling of the electricity and gas distribution grid to reduce renewable energy curtailment[END_REF], [START_REF] Frischmuth | Hydrogen sourcing strategies and cross-sectoral flexibility trade-offs in net-neutral energy scenarios for Europe[END_REF]. Further, Frischmuth and Härtel [START_REF] Frischmuth | Hydrogen sourcing strategies and cross-sectoral flexibility trade-offs in net-neutral energy scenarios for Europe[END_REF] demonstrated that a system consisting of both centralized and decentralized PtG systems show lower energy system costs and fewer power generation investments plus other flexibility options. However, designing and operating PtG systems is very complex as the flow of multiple energy vectors in the coupled system must be optimized such that they meet the design criteria as defined by the designer or operator. Further, how this design criteria -or optimization objective function -is evaluated depends upon the desired technique for formulating the optimization problem. The limitations of the equipment as well as internal system and external networks -also known as constraints -must also be considered during optimization. The influence of variable renewable energy, energy markets, energy loads, energy demands and other input data further add to the complexity of modeling a PtG system. The objective of this paper is to present a PtG plant model which can be used to simulate operation over the project lifespan for ideal control strategy investigation and overall project feasibility analysis. A pilot plant will also be presented which is used as a case study of model application. The plant analysis is done by taking input parameters from the user based upon local conditions (i.e. electricity prices, feedstock limitations such as generation profiles), plant configuration (i.e. equipment to be used), end-use application(s) and simulation objective and converting the data into an optimization problem to be solved.

The structure of the rest of the paper will be as follows: a section discussing current state-of-the-art techniques used to model PtG plants will be presented with results of the simulation. Literature knowledge gaps will be shown and how they correlate to the scientific contribution of this work.

State-of-the-Art

McDonagh et al. [START_REF] Mcdonagh | Modelling of a power-to-gas system to predict the levelised cost of energy of an advanced renewable gaseous transport fuel[END_REF] studied synthetic production as a transport fuel, calculating a cost of 107-143 A C/MWh in 2020 when assuming an average electricity price of 35 A C/MWh and 6,500 full-load operational hours. However, this study focuses on a national level (Ireland) and does not discuss cost reduction strategies for individual plants, nor does it use optimization techniques were used for simulation. Breyer et al. [START_REF] Breyer | Power-to-Gas as an Emerging Profitable Business Through Creating an Integrated Value Chain[END_REF] modeled a PtG plant used on-site of a pulp mill in Finland with ancillary market participation and 4,000 full-load operational hours and found the business case to be profitable, with grid services representing 40% of total income. However, this model used an assumed fixed electricity price, did not consider limitations on local gas consumption, and simplified Excel ® spreadsheets for economical calculations. Gorre et al. [START_REF] Gorre | Production costs for synthetic methane in 2030 and 2050 of an optimized Power-to-Gas plant with intermediate hydrogen storage[END_REF] performed an optimization analysis of production costs of SNG in terms of various electricity purchasing and gas selling strategies for a PtG plant in 2030 and 2050, comparing them to the expected SNG prices. The results showed that flexible operation in day-ahead markets with an average price of 30 A C/MWh and 4,000 full-load operational hours can have a positive revenue by 2050. However, the study is generic to Europe and does not consider participation in ancillary services, local technical or physical limitations, governmental support schemes or specific country day-ahead electricity prices. Further the optimization method is not mentioned or presented. Abdelghany et al. [START_REF] Abdelghany | Optimal operations for hydrogen-based energy storage systems in wind farms via model predictive control[END_REF] provide a mixed-integer linear programming method for optimizing energy management in a wind-powered PtG plant for power generation to meet an energy load. Their style includes equipment states and detailed operational constraints, however, its purpose is for real-time energy management not project lifespan analysis and thus uses very short time horizon without overall project feasibility analytics. Thus, from this brief literature review it can be seen that models which include plant operation are either very simplified, ignoring all possible limitations and opportunities of the system, or too detailed of minor plant dynamics and miss the overall picture of project feasibility. The model presented hopes to marry these differences into a model which captures major equipment dynamics, local feedstock or end-use constraints, and lifetime uncertainties to provide tangible results for project feasibility.

Plant Operation Simulation Model: An Optimization Problem

To simulate operation of a PtG plant over the project lifespan, the system can be converted into an optimization problem, using mathematical formulation to represent all the equipment, feedstock parameters, external system limitations, enduse applications and operational objectives. A unit commitment (UC) style model is used to form a mixed-integer linear optimization problem. There are several possible objective functions which can be maximized/minimized based upon simulation objectives. The formulated UC problem of the system is defined in Equations ( 1)-( 3) and extracted from [START_REF] Knueven | On Mixed Integer Programming Formulations for the Unit Commitment Problem[END_REF]:

min g∈G t∈T c g (t) (1)
s.t.

g∈G A g (P g , P g , δ g ) + N (s) = L (2)

(P g , P g , c g ) ∈ Π g ∀g ∈ G (3)
where the objective function ( 1) is the minimization of the summation of the operational cost of each generator c g at each time step t in the simulation period T . Other objective functions could be chosen but only minimization of operational cost is defined here. The Equation (2) matrix A g (P g , P g , δ g ) determine how the generators and interact with system requirements L while N (s) refers to other potential decision variables involving operation of the system. P g , P g and δ g represent the feasible power, maximum power and status of the generator, respectively. Finally, Equation (3) defines the feasible region for each generator and the cost associated to it, with the set Π g representing the description of this feasible operational region for each generator. For reference, generators are defined as equipment which produce either final products or internal system feedstock i.e. electrolysers, reactors, and gas purification systems.

Critical Constraints

There are several constraints which are required to define Π g and A g , however, only the critical constraints which dominantly determine plant operation will be shown here. The proposed model goes beyond standard UC models such that extra operational states of the generators are included which consider startup and idle of the equipment -states which play an important role in real operation of the PtG plant. The modeling strategy presented by Abdelghany et al. [START_REF] Abdelghany | Optimal operations for hydrogen-based energy storage systems in wind farms via model predictive control[END_REF] was modified and applied to a PtG UC model that can incorporate the "CLD" and "STB" states into the already defined "OFF" and "ON" states in the optimization problem formulation. The application of these states and possible transitions between them for plant generators is shown in a sequential graph in Figure 1. Once defined, these states are included in operational constraints to more accurately portray plant operation internally and with external systems.

Beyond internal constraints of equipment such as minimum and maximum production levels, ramping limitations and thermodynamic and kinetic processes, external systems have large governance over plant operation. If a PtG plant is producing hydrogen using exclusively wind farm power, the electrolyser must not consume more power than the farm produces at each time step in the simulation. This can be defined as shown in Equation ( 4):

P ez (t) ≤ P wind (t) ∀t ∈ T (4)
where P ez (t) is power consumption of the electrolyser at time step t and P wind (t) is the wind farm power generation. If the produced H 2 is being injected into the local NG grid or the H 2 is being further converted into SNG in a methanation reactor and and injected in the grid, any local grid capacity limitations must also be considered. This can be expressed as constraints as shown in Equations ( 5) and ( 6):

Q grid H2 (t) ≤ d H2 N G Q N G,avail (t) ∀t ∈ T (5)
Q grid SN G (t) ≤ Q N G,avail (t) ∀t ∈ T (6)
where Q grid H2 (t) and Q grid SN G (t) are the volumetric production of H 2 and SNG at time step t, respectively, Q N G,avail (t) is the volumetric availability in the NG grid at time step t and d H2 N G is the volumetric ratio of H 2 accepted in the local NG grid. As operational cost is to be minimized and the cost of electrolysis is the majority of this cost [START_REF] Van Leeuwen | Report on the costs involved with PtG technologies and their potentials across the EU[END_REF], the price of electricity on the day-ahead market is of great interest. This price changes hourly and can have large variations depending upon the time of day, season or unforeseen events such as nuclear reactor shutdown. The cost of electrolyser operation is thus the most critical cost for simulation consideration and is summarized in Equation ( 7):

c ez (t) = Q ez H20 (t)c H20 + Q H2 (t) η ez c el (t) (7)
where c ez (t) is the total operational cost of the electrolyser at time step t, Q ez H20 (t) is the volume of water consumed at time step t, c H20 is the cost of water, Q H2 (t) is the volumetric production of H 2 at time step t, η ez is the efficiency of the electrolyser and c el (t) is the day-ahead market price of electricity at time step t.

Case Study

A PtG pilot plant named HYCAUNAIS has been used for a case study of the model. Located in Saint-Florentin, France the project is being led by Storengy and has several industrial and public partners [START_REF]Agence de la transition écologique[END_REF]. A diagram of the HYCAUNAIS plant layout and possible configurations is shown in Figure 2. The plant will produce gases with a low environmental footprint for NG grid injection and possibly mobility. Following the combined power production signal of two existing wind farms (1) provided by the wind farm operator and sourcing the electricity via a grid connection (2), a 1 MW el PEM electrolyser (3) will produce hydrogen for a 50 Nm 3 /h biological methanation reactor (4), with intermediate H 2 storage (5) installed between the electrolyser and reactor. Currently, gas from the landfill (6) on-site is being upgraded by a WAGABOX ® unit [START_REF] Frischmuth | Hydrogen sourcing strategies and cross-sectoral flexibility trade-offs in net-neutral energy scenarios for Europe[END_REF], recovering 90 vol% of the bio-CH 4 contained in it. This highly pure bio-CH 4 is being injected into the natural gas (NG) grid [START_REF] Mcdonagh | Modelling of a power-to-gas system to predict the levelised cost of energy of an advanced renewable gaseous transport fuel[END_REF]. The CO 2 stream (including some impurities) normally vented during upgrading will be utilized as the carbon source for methanation, which must be purified [START_REF] Breyer | Power-to-Gas as an Emerging Profitable Business Through Creating an Integrated Value Chain[END_REF] prior to injection into the biological reactor. Intermediate storage of CO 2 [START_REF] Gorre | Production costs for synthetic methane in 2030 and 2050 of an optimized Power-to-Gas plant with intermediate hydrogen storage[END_REF] between the reactor and purifier is also used. Bio-NG captured during purification is combined with the produced SNG from the reactor and injected into the NG grid at the same existing site used by the WAGABOX ® unit. Compressed natural gas (CNG) [START_REF] Abdelghany | Optimal operations for hydrogen-based energy storage systems in wind farms via model predictive control[END_REF] and H 2 [START_REF] Knueven | On Mixed Integer Programming Formulations for the Unit Commitment Problem[END_REF] mobility are also being investigated as secondary applications to improve operational hours of the plant and profitability while NG grid injection is seen as the primary application. In addition to these end-use applications, electrolyser participation in electricity grid ancillary services is also being investigated with the purpose of adding an additional revenue stream while only minorly effecting overall operation. The project boundary is input from the electrical grid and WAGABOX ® effluent to output to NG grid injection as shown in Figure 2.

Based upon the plant description given above, the optimization model presented in section 3. can be used to simulate operation of the HYCAUNAIS pilot plant as a case study of model application.

Results

Case study results of model application will be presented using two objective functions: maximizing SNG production and minimizing operational cost. The objectives were modeled separately to analyse the production behaviour on the the electrolyser and reactor. The simulation was done over a week time horizon at time steps of one hour. Only NG grid injection has been considered in these results for end-use application; mobility is not included. SNG production from the reactor is the primary end-use application of the plant and, hence, was chosen to be the technical objective to maximize. As the maximum operational output of the reactor is desired in this configuration, the only limitations on output gas flow is reactor rated capacity, ramping constraints, and NG grid capacity (referred to as NG consumption). These limitations are demonstrated in the top plot of Figure 3 and their effect on gas flow production in the bottom plot. When NG consumption is below reactor rated capacity, it means the NG grid cannot accept full gas output but only what is currently being consumed locally. Reactor production thus is forced to be limited and eventually terminated until consumption increases again. A similar effect on electrolyser operation can be seen with wind power: the electrolyser can only consume electrical power, and therefore produce H 2 , up to what is available from the wind farm. Electrolyser operation is thus limited several times in the last third of the week. Due to the limited size of the H 2 storage -as seen as the orange dotted line in Figure 3, represented by internal gas pressure -the reactor cannot operate for long before the stock is depleted. Figure 4 shows the critical input data and output production results when the objective function is changed to minimizing operational cost. As price of electricity is now a main concern for optimization (as price of electrolysis is most influential on cost of PtG plants [START_REF] Van Leeuwen | Report on the costs involved with PtG technologies and their potentials across the EU[END_REF]) and NG consumption limits reactor production for a short time, the former has been replaced with the latter in the parameters plot of Figure 4. Further, an additional constraint of 60% full-load operational hours of the reactor in the simulated horizon has been applied, forcing the solver to find the cheapest hours to operate for that minimum time. As the wind power limitations are still a factor in the later hours of the week, the equipment runs full-load during midweek and then chooses cheaper times to run at capacity while still considering wind power in the end of the week. As one could imagine, minimizing operational cost is a much more complex problem to solve when compared to maximizing SNG production as economics of operation increase considered constraints and decision variables. Indeed, solve runtime of maximizing SNG production for the 168 hour (7 days) horizon is 17.8 seconds while it took 6847 seconds (114 minutes) to solve minimizing operational cost, over 380 times slower.

Conclusion and Future Work

A power-to-gas plant operation simulation model has been proposed in this paper to be used for feasibility analysis of potential projects over their lifespan. The model uses optimization to achieve this by converting possible plant configurations into a mathematical problem formulation to be solved. A case study of a pilot plant under construction in France is presented to apply the model for demonstration. By implementing two objective functions separately (maximize reactor production and minimize operational cost), the influence of critical parameters on equipment output was shown to ultimately decide optimal values.

From the results shown, it can be seen that attempting to improve the solver runtime is of great interest in future work. Improving computational efficiency through other problem formulation methods or solvers will be investigated.

Once the formulation method has been chosen, increasing applicability of the model to other PtG scenarios will be done, such as: CNG and H 2 mobility, ancillary services participation, multiple end-use applications, multiple objective functions and optimization of equipment sizing. Further, including forecasting models of variable renewable energy generation, energy markets and energy loads will improve accuracy of model results in future projects.

Figure 1 :

 1 Figure 1: Generator operational states sequential graph. Modified from [11].

Figure 2 :

 2 Figure 2: HYCAUNAIS project plant schematic and process flow.

Figure 3 :

 3 Figure 3: Critical input data (top) and out production results (bottom) are shown when maximizing SNG production is the objective function of simulation optimization. Simulation period is one week.

Figure 4 :

 4 Figure 4: Critical input data (top) and output production results (bottom) are shown when minimizing operational cost is the objective function of simulation optimization. First segment of zero reactor production is due to no NG grid capacity available, which is not shown in the parameter plot to avoid figure clutter. It can be seen in Figure 3. Simulation period is one week.

Acknowledgements

This research was supported by the HYCAUNAIS project (agreement number 1882C0096) with financial support provided by the Agency of Ecological Transition in France (ADEME) under the "Investissements d'Avenir de l'Etat" aid programme, the EIPHI Graduate School (contract number ANR 17-EURE-0002) and the Region Bourgogne Franche-Comté and Bpifrance under "fond régional d'investissement" (FRI).


Member, IEEE Konstantinos Alexandris

Member, IEEE Chia-Yu Chang

Member, IEEE Navid Nikaein

Utility-based Opportunistic Scheduling under Multi-Connectivity with Limited Backhaul Capacity

Keywords: 5G, Multi-connectivity, Backhaul, QoS

Multi-connectivity is a 5G key enabler to fulfill multi-service quality-of-service (QoS) requirements. This work examines an opportunistic resource allocation problem under multi-connectivity and limited backhaul capacity in evolved LTE using two utility functions: proportional fairness (PF) and utility proportional fairness (UPF). We then propose an efficient algorithm to deal with the formulated NP-hard problem. Such algorithm guarantees to produce a solution with an explicit bound on the optimal solution. Finally, simulation results justify that UPF utility function with opportunistic scheduling in multiconnectivity can better satisfy QoS requirements.

I. INTRODUCTION

T O flexibly utilize available radio access technologies (RATs) and spectrum bands, simultaneous multiple connections per end-user perspective can be exploited. Such concept is termed as multi-connectivity [START_REF] Ravanshid | Multi-connectivity functional architectures in 5G[END_REF] and is highlighted by 3GPP as one key enabler toward 5G [START_REF]Multi-connectivity; Overall description; Stage-2[END_REF]. In practice, the multiconnectivity deployment can mix-and-match with several 5G new radio techniques, e.g., millimeter wave [START_REF] Maeder | A scalable and flexible radio access network architecture for fifth generation mobile networks[END_REF], to exploit its benefit. The authors of [START_REF] Tesema | Mobility modeling and performance evaluation of multi-connectivity in 5G intra-frequency networks[END_REF] study the improvement of cell-edge user throughput as well as mobility robustness. In [START_REF] Liu | Dual connectivity in backhaul-limited massive-MIMO HetNets: User association and power allocation[END_REF], the operating expenditure is minimized via jointly considering the user association and power allocation under dual connectivity. Our prior work in [START_REF] Alexandris | Multi-connectivity resource allocation with limited backhaul capacity in evolved LTE[END_REF] stresses a better quality of service (QoS) satisfaction levering multiple connections. Also, it is necessary to consider the non-ideal backhaul condition [START_REF] Michalopoulos | User-plane multi-connectivity aspects in 5G[END_REF], e.g., limited backhaul capacity, to make use of multiple connections.

Nevertheless, current studies of multi-connectivity mostly focus on exploiting its versatility, while less effort is made to allow the scheduling of proper transmission resources over several connections for each user, termed as opportunistic scheduling in [START_REF] Knopp | Information capacity and power control in single-cell multiuser communications[END_REF]. Such opportunistic scheduling approach aims to take the channel quality information into account and can provide several advantages [START_REF] Asadi | A survey on opportunistic scheduling in wireless communications[END_REF], like improving capacity, QoS satisfaction, and power efficiency. An example shown in [START_REF] Grøndalen | Scheduling policies in time and frequency domains for LTE downlink channel: A performance comparison[END_REF] exploits channel variability to allow for better QoS awareness. In summary, the opportunistic scheduling approach is complementary to the multi-connectivity scheme and can further unleash its potential.

To this end, this work focuses on the utility-based opportunistic scheduling in a multi-connectivity deployment with capacity-limited backhaul network. Also, we propose an efficient algorithm for the formulated NP-hard optimization problem maximizing the network utility function. To the best of our knowledge, this is the first work exploring the opportunistic The authors are with Department of Communication Systems, EURECOM, 06410 Sophia-Antipolis, France (e-mail: firstname.lastname@eurecom.fr). This work was supported by the EU's H2020 Framework Programme under Grant 671639 (COHERENT) and Grant 762057 (5G-PICTURE).

scheduling in multi-connectivity and providing an algorithm with an explicit bound on the optimal solution.

II. SYSTEM MODEL

We consider an area served by a set of base stations (BSs)

B = {b 1 , • • • , b |B| } connected to the core network (CN) via the backhaul network to serve a set of user equipment (UEs) U = {u 1 , • • • , u |U | } distributed in this area.

A. Air-interface model

A.1-Carrier frequency: Multi-connectivity deployments can be divided into: (a) Intra-frequency and (b) Interfrequency [START_REF] Maeder | A scalable and flexible radio access network architecture for fifth generation mobile networks[END_REF]. The inter-frequency deployment stands for the case where a UE is connected through multiple carrier frequencies, either in a single or multi RATs. The intra-frequency one refers to transmissions over the same frequency from a single RAT, where the coordinated multipoint (CoMP) processing is required to mitigate interference. This work studies the interfrequency single LTE RAT deployment and denotes the carrier frequency as f j for the j-th BS. It is noted that the current setup can be extended also to scenarios where 5G NR small cells are collocated with an LTE macro cell in non-standalone operation mode [START_REF] Lien | 5G new radio: Waveform, frame structure, multiple access, and initial access[END_REF].

A.2-Physical data rate: The physical data rate that the i-th UE can get from the j-th BS using the k-th sub-channel1 in bps is R j,k,i = W j,k log 2 (1 + SINR j,k,i ), where W j,k is the bandwidth of the k-th sub-channel at the j-th BS in Hz and

SINR j,k,i = RSRP j,k,i bm =bj , fm=fj RSRP m,k,i + W j,k N 0 (1)
is the signal-to-interference-plus-noise-ratio (SINR) of the received signal from the j-th BS to the i-th UE of the kth sub-channel. The reference signal received power (RSRP) RSRP j,k,i = L j,i s j,k,i G j,i P j,k,i includes (a) large-scale fading as L j,i integrating pathloss and shadowing effects from the jth BS to the i-th UE, (b) small-scale fading as s j,k,i comprising the multi-path effects from the j-th BS to the i-th UE over the k-th sub-channel, (c) combined antenna gain of the j-th BS and the i-th UE as G j,i , and (d) transmission power from the j-th BS to the i-th UE of the k-th sub-channel as P j,k,i . Note N 0 is the thermal noise density, and thus the product W j,k N 0 shows the aggregated noise power.

B. Connection and Traffic model B.1-Multi-connectivity:

To identify the feasible connections between a single UE and mutliple BSs, the average SINR over all sub-channels is used as the decision criterion. Specifically, when the average SINR from the j-th BS to the i-th UE SINR j,i = 1 |K| k∈K SINR j,k,i is above a threshold, SINR th , then the connection in between is feasible. Thus, a set E {(b j , u i ) : SINR j,i > SINR th , ∀ b j ∈ B, ∀ u i ∈ U } is formed to represent all feasible connections.

B.2-Backhaul network: For the connected backhaul link from the CN to the j-th BS, we denote its capacity as C j in bps. Further, we consider a homogeneous backhaul capacity for different BSs in this work (i.e., C j = C, ∀j) and leave the heterogeneous case (i.e., C i = C j , ∀i = j) for future work.

B.3-Requested traffic rate: The requested traffic rate Ri of the i-th UE is defined as the sum rate required by a number of user applications running on the UE side in bps. Such requested rate can be provided through multiple BSs and multiple sub-channels to fulfill multi-service QoS requirements.

III. PROBLEM SETUP A. Network utility function

We define x j,k,i ∈ {0, 1} as the binary variable representing the allocation or not of the k-th sub-channel to the i-th UE from the j-th BS. Hence,

y j,k,i = x j,k,i • R j,k,i is the corresponding physical data rate.
First of all, we introduce two applied utility functions [START_REF] Abdel-Hadi | A utility proportional fairness approach for resource allocation in 4G-LTE[END_REF]; a) Proportional Fairness (PF): The utility function is defined as Φ (y) = log (y) that adopts "proportional fair" allocation; b) Utility Proportional Fairness (UPF): Such utility function aims to achieve the UE's requested traffic rate as introduced in Sec.II-B.3 applying the family of the natural logarithm of the sigmoid utility function Φ (y) = log 1/(1 + e -γ(y-β) ) , γ > 0, y ≥ 0, where β ≡ Ri . Then, we form two vectors as x i = x 1,1,i ,...,x 1,|K|,i , x 2,1,i ,...,x 2,|K|,i ,...,x |B|,1,i ,...,x |B|,|K|,i T is the vector form of binary allocation variables and

R i = R 1,1,i ,...,R 1,|K|,i , R 2,1,i ,...,R 2,|K|,i ,...,R |B|,1,i ,...,R |B|,|K|,i
T represents the vectorized physical data rate 2 . Finally, the network utility function of the i-th UE can be expressed as

U (x i ) Φ R T i • x i
, where R T i • x i stands for a user's aggregated data rate over its connected BSs and allocated sub-channels and Φ (•) can be any of the utility functions as presented above.

B. Problem formulation

Optimization Problem: The problem falls into the category of the network utility maximization for resource allocation and is given as follows:

max xi∈{0,1} |B|×|K| ui∈U U (x i ) (2)
s.t. C1:

ui∈U k∈K x j,k,i ≤ B j , ∀ b j ∈ B, (3)
C2: bj ∈B k∈K x j,k,i ≤ M i , ∀ u i ∈ U, (4)
C3:

ui∈U x j,k,i ≤ 1, ∀ b j ∈ B, k ∈ K, (5)
C4:

ui∈U k∈K x j,k,i R j,k,i ≤ C j , ∀ b j ∈ B. (6)
2 The (•) T denotes the transpose operator.

Algorithm 1: Opportunistic scheduling algorithm in multiconnectivity with complexity Objective function: The objective is to allocate the radio resources x j,k,i ∈ {0, 1} for each UE u i ∈ U to maximize the sum of utility functions over all UEs. Both PF/UPF utility functions described in Sec.III-A can be applied. Such objective function considers the channel variation across sub-channels (hence the term opportunistic) and is different from the nonopportunistic scheme, where SINR j,k,i = SINR j,i , ∀ k ∈ K.

O (|B| • |K|) 2 • |U | Input : R j,k,i , Ri , B j , M i , ∀ b j ∈ B, k ∈ K, u i ∈ U Output: S 1 S 0 ← ∅; t ← 0; /*
Constraints: There are four major constraints in the formulated utility maximization problem. C1: It ensures that the number of allocated sub-channels of BS b j to all UEs will not exceed the total number of available sub-channels at b j as B j . C2: It assures the total number of allocated sub-channels to UE u i among all connected BSs will not exceed the maximum allocatable sub-channels at u i as M i , which depends on the user capability defined as the UE-category in 3GPP TS36.306. C3: It ensures that each single sub-channel will be assigned up to one UE simultaneously implying sub-channel exclusivity. C4: It guarantees that the aggregated rate at each BS b j will not exceed the provisioned backhaul capacity C j .

C. Proposed algorithm

The proposed optimization problem is an NP-hard combinatorial problem. To the best of our knowledge, there is not such an algorithm to provide guarantee on the solution given by a bound due to the non-matroid nature of constraint C4, i.e., affine transformation with weights R j,k,i not equal to one [START_REF] Krause | Submodular Function Maximization[END_REF]. In consequence, to provide an explicit bound, we initially exclude the constraint C4 of Eq. ( 6) from the optimization process and we apply it afterwards as discussed below.

Lemma 1. The formulated optimization problem is with submodular and monotone objective function (Eq. (2)) and three matroid constraints (Eq. (3), Eq. (4), and Eq. ( 5)). Theorem 1. Let OPT be the optimal solution of the formulated problem and S be the output of Alg. 1. Then, it holds that f (S ) ≥ 1 4 OPT, where f (•) is the objective function of Eq. (2) rewritten as a set function in Eq. (10) of Appendix A.

Proof. The proof is given in

Proof. Define I = L l=1 I l , where I l is matroid constraints. Even I is generally not a matroid, the output of greedy algorithm of Alg. 1 is guaranteed with f (S ) ≥ 1 L+1 OPT shown in [START_REF] Krause | Submodular Function Maximization[END_REF]. As L = 3, i.e., three matroid constraints, in the formulated problem, we prove Theorem 1.

Subsequently, to deal with the backhaul capacity limitation, i.e., C j , introduced in the pre-excluded constraint C4, we further scale the allocated data rate of the i-th UE as

R i bj ∈B k∈K (x j,k,i ) R j,k,i • min Cj u i ∈U k∈K (x j,k,i ) R j,k,i , 1
, where (x j,k,i ) are the outcomes of the proposed algorithm. In reality, such scaling can be applied at the level of BS or gateway. For instance, each BS can extend legacy rate adaptation approach to be aware of the backhaul condition, while the gateway can leverage the flow shaping approach to indirectly restrict the resource allocation. Last but not least, the continuous relaxation technique is not applied here, since it does not guarantee how far the solution is from the optimal one and violates discrete resource block assignment used in LTE system.

IV. SIMULATION RESULTS

Two major scenarios are simulated 3 using self-developed MATLAB system-level simulator with the proposed algorithm: Scenario A-Empty cell: There is 1 BS that serves 0 UE, i.e, 0% traffic load, while the rest 2 BSs serve non-zero UEs. Such empty cell traffic load can characterize some real-world measurements. For instance, the phantom cell use case [START_REF] Ishii | A novel architecture for LTE-B: C-plane/U-plane split and phantom cell concept[END_REF] aims to deploy small cells that serve fewer UEs than macro cell in order to boost user plane performance via exploiting underutilized radio resources. The 0/z/z notation represents such scenario where z is the number of non-zero served UEs per BS. Scenario B-Loaded cell: All cells serve non-zero UEs and we use z/z/z notation to represent this scenario. The value of z is 2 in our simulations 4 . In both scenarios, backhaul capacity for all BSs is considered to be homogeneous. Then, our considered performance metrics include (i) user satisfaction ratio as the percentage of users that are satisfied with the allocated data rate, i.e., P sat i Prob{R i ≥ Ri } and (ii) unsatisfied normalized error as the normalized Euclidean 3 Simulation parameters applied to UEs, BSs and network planning are mostly taken from 3GPP (TR36.814, TR36.942, TR25.942) and NGMN documents, and some important parameters are listed in TABLE I. 4 Without loss of generality, we focus on the impact of zero and non-zero cell users and any variation of z can cover this case study. distance from the allocated rate to the requested traffic rate formulated as

E i = Ri-Ri Ri , if R i < Ri , 0 , otherwise. (7)
In Fig. 1, empty cell scenario results are shown with several combinations of PF/UPF with single-/multi-connectivity. We can see that using UPF utility function can boost the satisfaction ratio even with single-connectivity; however, some QoS requirements can never be fulfilled even we continuously increase the backhaul capacity. It is due to the lack of utilizing unallocated resources at the empty cell, and thus multi-connectivity shows its advantage. Further, using UPF with multi-connectivity can better fulfill QoS requirements even with fewer backhaul capacity. Then, the results of loaded cell scenario are shown in Fig. 2. Using multiple connections is not enough to fulfill QoS requirements if only PF is adopted. Such observation is due to there are no unallocated resources compared to the empty cell case; hence, we need to efficiently reshuffle allocated resources via using UPF.

Finally, Fig. 3 shows the opportunistic scheduling benefit in loaded cell scenario with UPF and multi-connectivity compared to the non-opportunistic one. Note the non-opportunistic scheduling can optimally rely on the average SINR over all sub-channels, i.e., SINR j,i ; however, its worst case is based on the minimum SINR among all sub-channels, i.e., min k∈K (SINR j,k,i ). In comparison, the opportunistic scheduling can better utilize available sub-channels to provide a smaller unsatisfied normalized error compared to the best or the worst case of the non-opportunistic scheduling.

V. CONCLUSION AND FUTURE WORK

In this work, we formulate the utility-based opportunistic resource allocation problem under multi-connectivity and limited backhaul capacity. The proposed algorithm can efficiently solve the problem with an explicit bound on the optimal solution. The simulation results justify the benefit of UPF utility function with opportunistic scheduling in terms of QoS satisfaction ratio and unsatisfied normalized error. As a future Fig. 3: Unsatisfied error versus backhaul capacity (z/z/z) direction, the efforts will be given on the implementation applicability considering different RAN deployment flavors.

APPENDIX A PROOF OF LEMMA 1

We define the function h The following holds with inequality for every P ⊆ Q ⊆ N :

(x) : {0, 1} |N | → R expressed as h (x) = g a T • x = g   |N | n=1 a n x n   , (8)
n∈Q a n ≥ n∈P a n ⇔ g n∈Q a n ≥ g n∈P a n , (9)
since g is monotone proving the monotonicity of h.

Proposition 1. If a n ∈ R + and g (•) is concave, then the function h (W) = g n∈W a n is submodular. Proof. The proof is based on part (a) of Lemma 2.6.2. in [START_REF] Topkis | Supermodularity and complementarity[END_REF].

Applying Proposition 1, we prove the submodularity of h (•).

The objective function in Eq. ( 2) can be rewritten as the set function f (S) : 2 V → R that is expressed as

f (S) = ui∈U Φ   bj ∈B k∈K R j,k,i 1 S (x j,k,i )   , (10)
where

V = B × K × U, 1 S (•) is the indicator function on the set S = {(j, k, i) ∈ V : x j,k,i = 1}.
Proposition 2. Φ (•) utility functions are monotone and strictly concave.

Proof. The PF utility function is the logarithmic function and thus is monotone and strictly concave. While the UPF takes the natural logarithm of sigmoid function proven as monotone and strictly concave following Lemma III.1. in [START_REF] Abdel-Hadi | A utility proportional fairness approach for resource allocation in 4G-LTE[END_REF].

As each term in the summation of objective function f (•) lying on the family of h (•) in Eq. ( 8), where Φ (•) ≡ g (•) is monotone and concave, N ≡ B × K, a n ≡ R j,k,i and x n ≡ x j,k,i . Thus, we conclude that the objective function f (•) is monotone and submodular as it is the summation of monotone and submodular functions.

Finally, we define the three constraints in Eq. ( 3)-( 5) equivalently to the sets I l , l ∈ {1, 2, 3} in Eq. ( 11) as follows:

I 1 = {S ⊆ 2 V : |S ∩ 2 {j,K,U } | ≤ B j , ∀ b j ∈ B}, (11a)
I 2 = {S ⊆ 2 V : |S ∩ 2 {B,K,i} | ≤ M i , ∀ u i ∈ U}, (11b)
I 3 = {S ⊆ 2 V : |S ∩ 2 {j,k,U } | ≤ 1, ∀ b j ∈ B, k ∈ K}. (11c)
To prove these constraints are matroid (see definition in [START_REF] Krause | Submodular Function Maximization[END_REF]) with V as ground set, two properties shall be fulfilled: (a) For all sets P and Q with P ⊆ Q ⊆ V, if Q ∈ I l , i.e., the sub-channel allocation defined by Q does not violate the constraint I l , then P ∈ I l is held, i.e., the constraint I l is not violated by sub-channel allocation in P, since P ⊆ Q. All three aforementioned constraints possess such property.

(b) For all sets P, Q ∈ I l and |Q| > |P|, i.e., more subchannels are allocated in total defined by Q, ∃ e ∈ Q\P such that P ∪ {e} ∈ I l . Since not all sub-channels are allocated defined by P (otherwise, Q will either violate the constraint I l , i.e., Q / ∈ I l , or has no more allocated sub-channels as P, i.e., |Q| ≤ |P|), there is at least one sub-channel that can be further allocated and will not violate the constraint I l . Hence, such property is held for all three constraints. Thus, we conclude that I l , ∀ l ∈ {1, 2, 3} are matroid.

  Appendix A. Hence, we propose a greedy algorithm in Alg. 1 with f (•) as the objective function and I 1 /I 2 /I 3 as constraints detailed in Eq. (11) of Appendix A. Such algorithm can solve the problem with polynomial time complexity O (|B| • |K|) 2 • |U | with an explicit bound on the optimal solution.

Fig. 1 :Fig. 2 :

 12 Fig. 1: User satisfaction ratio versus backhaul capacity (0/z/z)

  where g (•) is monotone and concave, a = a 1 , • • • , a |N | T and a n ∈ R + . The above function is equivalent to the set function h (W) : 2 N → R, where h (W) = g n∈W a n and W = {n ∈ N : x n = 1}. Definition 1. A set function is characterized as monotone iff h (Q) ≥ h (P) for every P ⊆ Q ⊆ N .

  Parameter initialization */ 2 I ← I 1 I 2 I 3 ; /* Get intersection of three constraints */ 3 while t < bj ∈B B j do ← S t-1 v ; /* Include v in current solution St */

	4	t ← t + 1 ;
	5	v ← arg max q∈I\St-1 f(S t-1 ∪ {q}) ; /* v gets the largest gain */
	6	if v = ∅ then
	8	else
	9	break;
	10 S ← S t ; /* Output final solution S */

7

S t

TABLE I :

 I Simulation parameters

	Parameter	Value
	LTE mode	FDD, SISO, Band 1 (3×20 MHz bandwidth)
	Total PRBs of each BS (Bj )	100 (20 MHz bandwidth)
	Maximum allocated PRBs of UE (Mi)	100
	Number of BSs (|B|)	3
	UE traffic model	Full buffer
	SINR threshold (SINR th )	-12 dB
	BS transmission power (P j,k,i )	46 dBm
	Thermal noise density (N0)	-174 dBm/Hz
	Gamma of UPF utility function (γ)	10/β
	Requested traffic rate ( Ri) distribution	Uniformly distributed between 10 to 100 Mbps

The term sub-channel stands for the minimum radio resource that spans over time and frequency domains, a.k.a. physical resource block (PRB). It is noted that a sub-channel set K can be allocated by each BS.