T Ehret

A De Truchis

M Mazzolini

J.-M Morel

G Facciolo

AUTOMATIC METHANE PLUME QUANTIFICATION USING SENTINEL-2 TIME SERIES

Keywords: Methane, quantification, automatic, monitoring, satellite, sentinel-2

Methane emissions monitoring is essential to control methane pollution. In this paper, we propose an automatic practical methodology using time series to estimate the quantity of methane in a given plume using a multispectral satellite like Sentinel-2. Sentinel-2 proposes a low revisit time, a good spatial resolution and a low acquisition cost. Contrary to previous methods, the proposed approach does not require a manual selection of an optimal reference image. We compared its performance on an oil-and-gas site in Kazakhstan. This is the first step toward an automatic global monitoring system for methane plume detection and quantification with these satellites.

INTRODUCTION

Recently, the detection of large and frequent methane emissions has raised concerns in the ability of natural gas to effectively reduce greenhouse gas emissions as a substitute to coal [START_REF] Riley M Duren | California's methane super-emitters[END_REF][START_REF] Daniel | Intermittency of large methane emitters in the permian basin[END_REF]. Over a 20-year horizon, a methane molecule has a global warming potential about 85 times larger than carbon dioxide [START_REF]Anthropogenic and Natural Radiative Forcing[END_REF]. However, it has also been shown that a large part of methane emissions could be controlled or avoided [START_REF] David R Lyon | Aerial surveys of elevated hydrocarbon emissions from oil and gas production sites[END_REF]. This show that methane monitoring is now essential to reduce methane pollution.

Airborne campaigns [START_REF] Del | An evaluation of airborne swir imaging spectrometers for ch4 mapping: Implications of band positioning, spectral sampling and noise[END_REF] are often used to perform monitoring. While these campaigns have very good spatial resolutions and low detection thresholds, they suffer from a limited spatial coverage. Indeed, these campaigns can only cover a limited area and are costly to perform so they often have very limited revisit.

In order to perform a more global monitoring of greenhouse gas emissions produced by human activities, several satellites, such as GOSAT and TROPOMI [START_REF] Veefkind | Tropomi on the esa sentinel-5 precursor: A gmes mission for global observations of the atmospheric composition for climate, air quality and ozone layer applications[END_REF], have been placed in orbit over the past ten years. Contrary to airborne campaigns, these satellites cover the entire Earth daily and were used successfully to detect unintended very large emissions [START_REF] Pandey | Satellite observations reveal extreme methane leakage from a natural gas well blowout[END_REF]. However, it is difficult to detect small emissions due to its very low spatial resolution (5-7km). An alternative would be high resolution hyperspectral satellites such as PRISMA [START_REF] Daniel | Potential of next-generation imaging spectrometers to detect and quantify methane point sources from space[END_REF] and GHGSat [START_REF] Mckeever | First methane sensing results from ghgsat's commercial constellation[END_REF]. However, the tasking nature, the relatively small field of view and the cost of data makes them unsuitable for a recurrent global monitoring.

The Sentinel-2 mission presents many advantages: it is recurrent with a low revisit time, it covers the entire planet and has a low acquisition cost. Sentinel-2 images comprise 12 bands, from visible to SWIR, with spatial resolutions from 10m per pixel to 60m per pixel. Although it was not initially designed for methane detection, it has been shown by Varon et al. [START_REF] Daniel | High-frequency monitoring of anomalous methane point sources with multispectral sentinel-2 satellite observations[END_REF] that some the bands are sensitive enough to methane to allow detection of moderate CH 4 emissions.

In this work, we propose an automatic methodology to quantify methane plumes using a multi-spectral satellites such as Sentinel-2 or Landsat-8. Starting from a timeseries of Sentinel-2 images, the background image without any methane plume is estimated by doing a temporal regression. The quantity of methane corresponding to the plume is then estimated for each pixel using an approximated atmospheric model.

AUTOMATIC METHANE QUANTIFICATION FROM MULTI-SPECTRAL TIME SERIES

We focus on the quantification of isolated excess concentrations of methane in the atmosphere. As shown in Fig. 1, bands B11 and B12 of Sentinel-2A are impacted by the presence of methane in the atmosphere. This means that it is possible to detect the presence of methane based on the dimming caused by the concentration of the gas. The Beer-Lambert law states that for a light source with intensity I 0 and a wavelength λ

I = I 0 e -N i=0 Ai(λ)li , (1)
where the light goes through N gases defined by their absorption A i (λ) and equivalent optical path length l i defined as the product of the actual optical path and the concentration of the i th gas. In our case, the light source is the sun and I 0 corresponds to its intensity in the SWIR spectrum. It is usually assumed that I 0 is constant. Taking into account that the sensor of a satellite integrates over a band of wavelengths described by a sensitivity function s, the intensity of the light seen by a space-borne sensor becomes

I = I 0 s(λ)α(λ)e -γ N i=0 Ai(λ)li dλ, (2)
where the two passes through the atmosphere are taken into account in γ, which is a function of both the sun azimuth angle θ sun and the satellite view angle θ sat (both defined in Fig. 2 such as γ = 1 cos(θsun) + 1 cos(θsat) ). The albedo α(λ) models the reflection on the ground.

In practice, the atmosphere model can be well approximated with a simple fixed "pure methane atmosphere". This means that N = 1, l atm the quantity of methane in the atmosphere is considered fixed and the only absorption coefficient is A CH4 . We discuss this approximation in Section 4. In the presence of a methane emission, characterized by l em , the intensity of the light seen by the sensors then becomes

I em = I 0 s(λ)α(λ)e -γA CH4 (λ)(latm+lem) dλ. (3)
The problem is that α is usually unknown. We propose to estimate its contribution by estimating the observation without methane, also called background observation, I bg . When we assume that methane emissions are anomalous events, it is to be expected that most observations in a time series should not contain excess methane. So, if we suppose that α is rather stable in time, the time series can be used to estimate a methane free background model that can be compared with the current observation. Here, we compute the background for a given date as its linear regression over the previous dates. If we denote by I t the observation at time t, then the Fig. 2. Observation model for a satellite like Sentinel-2 with a passive sensor. regression computes the optimal weights w i that solve

min {wi} I t - t-1 i=0 w i I i 2 . (4)
Then the background is obtained as the linear combination

I bg = t-1 i=0 w i I i .
The estimation is performed on the logarithm of the images so as to limit the impact of abnormal high values in the SWIR bands. To further improve the background subtraction we combine this estimation with a band ratio that exploits the correlation between SWIR bands, similarly to the multiple-band single-pass (MBSP) from Varon et al. [START_REF] Daniel | High-frequency monitoring of anomalous methane point sources with multispectral sentinel-2 satellite observations[END_REF].

Since γ is known for each acquisition, this ratio only depends on the methane excess. Therefore, it is possible to estimate the value of l em as the solution of a simple optimization problem arg min lem

I em I bg -B12 s(λ)e -γA CH4 (λ)(lem+latm) dλ B12 s(λ)e -γA CH4 (λ)latm dλ 2 2
.

(5) This quantification scheme can also be adapted when using band ratio.

EXPERIMENTS

We consider an area of interest of size approximately 5x5 km 2 in Kazakhstan located on an oil-and-gas location. The size of this area is large enough so that the background estimation is not impacted too much by the presence of methane in the input images. We collected the L1C Sentinel-2 time series corresponding to this location and, using a cloud detection algorithm such as the one proposed by Dagobert et al. [START_REF] Dagobert | Cloud Detection by Luminance and Inter-band Parallax Analysis for Pushbroom Satellite Imagers[END_REF], we removed all images with a cloud coverage of more than 15% of the image pixels. We then manually detected the plume shown in this section. We compared the quality of the reconstruction using a different amount of images for the background estimation. From Fig. 3, we can see that using a longer time series allows for a better background reconstruction. Indeed, after subtraction of the estimated background, the plume becomes much more visible and therefore its signal-to-noise ratio (SNR) with respect to the background is improved.

In Fig. 4 we also compared the proposed method to different approaches proposed by Varon et al. [START_REF] Daniel | High-frequency monitoring of anomalous methane point sources with multispectral sentinel-2 satellite observations[END_REF]. Some of these methods require a manual selection of a reference image. So, in order to have a fair comparison, we selected the reference image that gave the best SNR out of all available for each method. Nevertheless, our automatic method still outperforms all of them.

EFFECT OF MODEL SIMPLIFICATIONS IN THE METHANE QUANTIFICATION.

Atmosphere model. We study how different assumptions on the atmospheric model ( 5) affect the methane quantification with a numerical example. Computations are done without the logarithm to simplify the notations, this however does not change the quantification. We also assume that the sun and viewing angles are zero, i.e. γ = 2 for the numerical applications. For that we suppose that our observation yields a given attenuation

δR = u(x, y, t) û(x, y, t) = 0.95, (6)
we then estimate the concentration from the attenuation using three different models.

For the first model, let us suppose that the atmosphere doesn't absorb any light, i.e.

∀f, e -2 l i=0 Ai(f )li = 1.

In this case, the observed attenuation δR is Table 1. Estimated excess methane using different atmospheric models with increasing complexity.

δR = u(x, y, t) û(x, y, t) = s(f )e -2A CH4 (f )lem df s(f )df . ( 8
For the second model, we consider the atmospheric methane by assuming an atmosphere purely composed of methane and excess methane. The atmosphere consists of l atm cm of methane, which leads to

u(x, y, t) û(x, y, t) = s(f )e -2A CH4 (f )(latm+lem) df s(f )e -2A CH4 (f )latm df . (9)
This case can be brought back to the case without atmosphere by considering a second sensibility s ′ such that s ′ (f ) = s(f )e -2A(f )latm . For the third and last model, we consider a complete, but fixed, atmospheric model as in ( 5) that also includes atmospheric methane (l atm = 1.3cm ≈ 1600ppb) generated using [START_REF] Steven | A new software tool for computing Earth's atmospheric transmission of near-and far-infrared radiation[END_REF]. Table 1 summarizes the estimated concentration of methane using the different models. We observe that modeling the atmospheric methane cannot be neglected as it induces a notable bias. On the other hand, the simulation of the rest of the atmospheric gases can be negleted.

Sentinel-2A vs Sentinel-2B. The Sentinel-2 constellation is comprised of two satellites: Sentinel-2A and Sentinel-2B. While they both have very similar sensitivities, they are not perfectly identical. We show here that it is important to take this difference into consideration during the quantification step. For that we consider the same attenuation given in [START_REF] Veefkind | Tropomi on the esa sentinel-5 precursor: A gmes mission for global observations of the atmospheric composition for climate, air quality and ozone layer applications[END_REF]. We then estimate the concentration for both configurations of satellites. Using the sensitivity s A of Sentinel-2A yields an estimated excess of CH 4 of 2393ppb. Doing the MBMP [START_REF] Daniel | High-frequency monitoring of anomalous methane point sources with multispectral sentinel-2 satellite observations[END_REF], SNR 20.25

Ours, SNR 58.18

Fig. 4. Comparison between multiple methane excess quantization methods. Our automatic method provides a much better background reconstruction allowing to discriminate and quantify the methane plume better.

same estimation with the sensitivity s B of Sentinel-2B yields 3194ppb. This show that it is very important to use the correct sensitivity during the estimation process.

Choosing a fixed methane atmosphere. As mentioned previously, we choose a fixed atmospheric quantity of methane of l atm = 1.6 ≈ 1800ppb. However, it has been shown that the quantity of methane in the atmosphere is location dependent. We estimated the error due to this approximation by estimating l em for different values of l atm . We saw that for l atm between 1700 and 1900 ppb, the error on l em was at most of the order of 30ppb. This justify selecting a single fixed value for the quantity of methane in the atmosphere independently from the location.

CONCLUSION

We have presented an automatic method for methane quantification using a multispectral satellite like Sentinel-2. Although in this paper we only presented results using Sentinel-2, we have had a similar success using Landsat-8 imagery. The next step is to use controlled releases and a source emission method like IME [START_REF] Daniel | Quantifying methane point sources from fine-scale satellite observations of atmospheric methane plumes[END_REF] to measure the total error of our quantification process with an actual ground truth. This will also help us define the precision of the quantification process based on the methane excess emitted from the calibrated source. Having an automatic pipeline paves the way for an automatic detection and quantification pipeline that could provide a global automatic monitoring of the Earth with a short revisit time for methane pollution.

  It was partly financed by Office of Naval research grant N00014-17-1-2552, DGA Astrid project « filmer la Terre » n o ANR-17-ASTR-0013-01, MENRT.

Fig. 1 .

 1 Fig.1. Methane transmittance spectrum for 1cm of methane, in gray, and Sentinel-2A spectral sensitivity for all its bands, one color for each band. Bands B11 (1568-1659 nm, in orange on the right) and B12 (2114-2289 nm, in green on the right) are impacted by the presence of methane in the atmosphere thus allowing its measurement.

Fig. 3 .

 3 Fig.3. Impact of the number of images used during the background estimation. A longer time series improves the SNR of the extracted plume in the residual thus helps its quantification. The SNR was estimated using a plume mask that was manually annotated.


Dinil Sebastien Soudan

Mon Divakaran

Eitan Altman

Vicat-Blanc Primet

Pascale Vicat-Blanc

Primet

EN INFORMATIQUE ET EN AUTOMATIQUE Extending Routing Games to Flows over Time

Keywords:

Routing game presents an interesting framework to analyse the practical problem of source routing in the Internet. It is particularly useful in quantifying the inefficiency of selfish user behavior that results in any transportation network without any central authority. This game assumes that the only user criteria for decision making is path cost. In this work, we take a step further, and model a routing game where user decision is based not only on path but also on time. We show that, under convex cost functions, this new routing game over time can be mapped to the classical routing game, thereby presenting a model that can exploit well-established results in the subject. Using a simple example, we demonstrate the usefulness of the model, and motivate the need for resource coordination to minimize inefficiency or cost.

Introduction

Routing games is a game-theoretical subject that sheds light on an important practical problem in the Internet (or for that matter, any transportation network) -the benefits and implications of routing traffic without any central authority. The 'routing' here refers to source routing, where an individual user is able to decide its route, or path, from the source to the destination.

The literature on routing games is wealthy, briefed in Section 2. While the previous works consider path selection as the only decision variable for a user, we explore, along with path selection, another important variable -time preference. We call this game, routing game over time, RGoT in short. The users in this game face a cost on the selected path which depends on the time as well as the total amount of traffic sent on this path.

The motivation for including the time criteria for decision making comes from the need to meet time constraints during data transfers. This has lead researchers to explore in-advance reservation, for example [START_REF] Burchard | Networks with advance reservations: Applications, architecture, and performance[END_REF] and [START_REF] Chen | Scheduling deadline-constrained bulk data transfers to minimize network congestion[END_REF]. We propose a model which enables to quantify the inefficiency -price of anarchy -of realizing such a system as a result of selfish decisions, with no centralized coordination. The model for RGoT in presented in Section 3. Therein, we define the concept of equilibrium and optimal allocation. We show that under convex cost function, allocation can be described using step functions. In Section 4 we prove that RGoT can be solved by solving an instance of classical routing game, and map some of the important results known for this kind of game. Finally, in Section 5, using a simple example we illustrate and motivate the need for coordination.

Related work

Since the introduction of the concept of equilibrium in transportation network by J. G. Wardrop in 1952, this kind of Nash equilibrium has been widely studied and extended. Among these, we note the introduction of price of anarchy and the results on classes of cost functions by Roughgarden et al. [START_REF] Roughgarden | Selfish routing[END_REF]. The dynamics of convergence to equilibrium has been investigated in [START_REF] Fischer | On the evolution of selfish routing[END_REF]. The concepts of non-atomic and atomic routing games are detailed in [START_REF] Nisan | Algorithmic Game Theory[END_REF]Ch. 18].

On the other side, while concepts of network flows and flows over time was introduced by Ford and Fulkerson, network flows have been extensively studied as summarized in [START_REF] Cormen | Introduction to Algorithms[END_REF]Ch. 26], development on the version considering time is recent [START_REF] Hall | Multicommodity flows over time: Efficient algorithms and complexity[END_REF]. But since it is used in several works on bandwidth and flow allocation, it is interesting and important to extend non-atomic routing games to routing games that take time into consideration. This paper presents this extension.

Routing games over time

In this section, we define the model for RGoT, and then extend the concepts of equilibrium flows and optimal allocations to this model.

Model

This section describes the RGoT model used in this paper. It is inspired from the non-atomic routing game model and from the model of flow over time. This model of routing game assumes infinite number of users. Each user is small, and as such doesn't contribute significantly to the cost. Request as presented in the following definition refers to a kind of request; meaning that the users that form a particular request have same constraints in terms of source, destination, and paths; and total aggregate v r demand (volume) (s r , d r )

source/destination pair t s r start time (assumed to be a rational) t d r deadline (assumed to be a rational) P r set of unique labels for s r -d r paths of G = (V, E) volume is known. We now proceed to give some important definitions. The notations used are listed in Table 1 for a given request r.

Definition 3.1 (Request). In the directed graph G, r (s r , d r , v r , t s r , t d r , P r ), is a request between vertices s r and d r , dates t s r , t d r , allowed to transfer over one or more paths from P r and with a total volume of v r . Every path in P r is unique, in the sense that if two requests use the same path (as set of edges) their labels will be different. In addition, we let P i P r .

These requests induce a traffic on the network which is characterized by the rate traffic going on each path. This is called allocation vector. Definition 3.2 (Allocation vector). If f p is the rate for request r over the path p ∈ P r , the allocation vector f is defined as being the vector

[f p ] p∈P with f p : R + → R + t →fp(t)
. In addition, we define link allocation as:

f e : R + → R + t →fe(t)
where f e (t) p∈P:e∈p f p (t)

This allocation vector gives the rates over time with which each request will send on its paths to achieve transfer of the specified volume during the time window. If it does so, it is said to be a feasible allocation. Link allocation gives the total rate seen on a link at a given date.

Definition 3.3 (Feasible allocation).

∀r ∈ R,

t d r t s r p∈Pr f p (t)dt = v r
Next, we define per-edge cost. It is assumed to be a piecewise constant function with regard to the time.

Definition 3.4 (Per-edge cost function). For an edge e, c(e, b, t) is the cost function for transferring at rate b at time t over link e. For the cost function, we furthermore assume: given I c , a finite partition of time (R + ) with rational breaks, c has the following structure: c(e, b, t)

J∈Ic t∈J c J,e (b) and c J,e (.) is assumed to be non-negative, continuous and non-decreasing.

Using the definitions of allocation vector and cost, allocation f incurs cost c(e, f e (t), t) at time t on edge e. From the assumption on the form of cost, cost on this interval can be split on intervals I c as done below. But before doing this, we define RGoT. In the remaining, we focus on the interval:

T min r∈R {t s r }, max r∈R {t d r }
All time breakst s r , t d r , and changes of cost function -are rationals. Hence, we can divide time axis in intervals of same lengths, such that all these breaks come at the boundaries by taking this length as the least common multiplier of denominators of interval lengths. This partitioning of time is used in Def. 3.7 to define the set of time intervals. It is illustrated in Fig. 1(a). This is an artifact, used later, to reduce this game to normal non-atomic routing games.

c p (f ) = J∈I e∈p t∈J c J,e (f e (t))dt (1)
Having defined the model, next sections define some specific feasible flows which arise from selfish user behaviors, equilibrium flows, and from social interest, optimal allocations.

Equilibrium Flows

As said before, users behave selfishly on their infinitely small fraction of volume belonging to one request. This leads to allocations which form a subset of feasible allocations, called equilibrium flows. We show that they can be defined using step functions with steps on I when cost is convex. The definition of equilibrium that follows, is similar to (characterization of) Wardrop equilibrium.

Definition 3.11 (Equilibrium). f is an equilibrium flow in (G, R, c) if:
1. f is a feasible allocation vector, and; 2. for every r ∈ R, for every interval J and J included in [t s r , t d r ] such that |J| = | J| and any p, p ∈ P r where t∈J f p dt > 0, c J,p (f ) ≤ c J,p (f ).

If f is an equilibrium flow, for every r ∈ R and any p ∈ P r , c p (f ) ≤ c p ( f ) with f such that, f p has been replaced by fp , and

t d r t s r f p (t)dt = t d r t s
r fp (t)dt, as it is a feasible flow and following second point of Def. 3.11 with I = J = T (f p is supposed to be null out of [t s r , t d r ] as it doesn't contribute to feasibility).

We proceed to show that step functions are suitable for equilibrium flows, when the cost is convex.

Proposition 3.12. If cost function c J,e (b) is convex, for any link allocation f e , there is one step function constant on each J ∈ I that transfers the same volume on each time interval for a cost at least as good.

Proof. We proceed in two steps.

(1) f e (t) is not better than constant function on J: Let J be an interval without arrivals/departures/cost-breaks, we suppose f e (t) is integrable,

J f e (t)dt = f J,e |J| and c J,e (b) be a convex function. Let µ(.) be the Lebesgue measure, µ(J) = |J| and η(.) µ(.)/µ(J). Obviously η(J) = 1 and:

J c J,e (f e (t))µ(dt) = µ(J) J c J,e (f e ( t µ(J)
))η(dt)

As η(J) = 1, using Jensen's inequality:

µ(J)c J,e ( J f e ( t µ(J) )η(dt)) ≤ µ(J) J c J,e (f e ( t µ(J) ))η(dt)
(from definition of η(.) and above equation)

⇒ µ(J)c J,e ( J f e (t) µ(J) µ(dt)) ≤ J c J,e (f e (t))µ(dt) (from definition of f J,e ) ⇒ µ(J)c J,e ( f J,e |J| µ(J) ) ≤ J c J,e (f e (t))µ(dt) ⇒ |J|c J,e (f J,e ) ≤ J c J,e (f e (t))µ(dt) (2)
(2) says that the cost of a function constant on J and equal to f J,e has a cost as good as any other function which transfer the same volume on J.

(2) For any link allocation f e (.), there is one step function as "good": First part of this proof can be applied on each intervals and f e (.) be replaced by a step function with a cost at least as good.

Proposition 3.13. If cost function c J,e (b) is convex, for any feasible allocation f there is f made of step functions which: [START_REF] Burchard | Networks with advance reservations: Applications, architecture, and performance[END_REF] for all p ∈ P has a cost c p ( f ) not worse than c p (f );

(2) has steps on I.

Proof. We start with one f . For each p ∈ P, we define fp as the step function which transfers the same volume as f p on each intervals of I but using a constant rate. This allocation is also feasible as it transfers the same volumes during the same intervals and thus volume constraints are satisfied. Furthermore, its steps are on I. fe are step functions, as sum of step functions fp . Using Eq. ( 1) and Prop. 3.12, c p ( f ) ≤ c p (f ).

Next we show that for a convex cost function, for any equilibrium flow, there is a corresponding equilibrium flow made of step function which is as good in terms of cost. If we come back to infinitely small users making requests, this basically means, they have no incentive of using anything other than a step function constant on intervals of T or I. Proposition 3.14. Equilibrium flows can be taken as step functions.

Proof. Equilibrium flows are feasible allocation and can thus be taken as step function for a cost not higher, as proved by Prop. 3.13.

It can be seen that with cost functions that are not non-negative, convex or non-decreasing in f e , the link allocation, this proposition does not apply. This is because, we can't exploit the structure of time intervals to define the step functions since it might be cheaper to reduce the duration of transfer while increasing the utilized rate, eventually without limits. As an example, this kind of allocation (equilibrium flow and later optimal allocation) as step functions of time do not apply for sub-linear cost as it is always better to group the utilization of the resources.

Social Cost and Optimal Allocation

Having defined equilibrium allocation that will result from selfish decisions of users realizing requests of R, we now define optimal allocation that minimizes the total cost charged to serve the requests. The total cost is referred as social cost. By dividing it by t p f p , average cost can be obtained. As in equilibrium flow, if the social cost is convex, we can take optimal allocation as step function as stated in Prop. 3.20. To begin with, we observe that convexity of cost function implies convexity of social cost.

C(f )

Remark 3.17. x → x.c J,e (x) is convex on R + . Since this is a product of two convex, non-decreasing and positive functions on R + . Proposition 3.18. If cost function c J,e (b) is convex, positive and non-decreasing, for any link allocation f e , there is one link allocation as step function, constant on each J ∈ I, that transfers the same volume on each time interval for a social cost at least as good.

Proof. As stated in Remark 3.17, x → x.c J,e (x) is convex on R + . Allocation are positive functions. Then, using Jensen's inequality (or Hermite-Hadamard inequality) we show that there is a step function with steps in I which is as good as any feasible allocation and transfer the same volume on each time interval. Proof. We start with one f . For each p ∈ P, we define fp as the step function which transfers the same volume as f p on each interval of I, but using a constant rate. This allocation is also feasible as it transfers the same volumes during the same intervals and thus volume constraints are satisfied, steps are on I. fe are step functions as sum of step functions fp . Using Eq. (3) and Prop. 3.18, C( f ) ≤ C(f ).

We conclude this section showing that there is an optimal allocation among the set of feasible allocations made of step functions. Proposition 3.20. Optimal allocation can be taken as step function with steps on I.

Proof. Optimal allocations are feasible allocation and can thus be taken as step function with step on I for a cost not higher as proved by Prop. 3.19.

In the next section, we exploit the structure of step function and revisit previously introduced definitions to establish equivalence between step allocations and generic allocations under convex costs.

Discrete Allocations -Reduction to Step Function

As users and regulator try to minimize their costs (equilibrium flow) or social cost (optimal allocation), all the feasible allocations of interest are step functions and their steps are on I (Prop. 3.13 and 3.19).

Hence, interesting feasible allocations are step functions with steps on I, they can be described by discrete value giving their value on the steps.

∀r ∈ R, ∀p ∈ P r , ∀J ∈ I\I r , f J,p = 0
This can be justified, as it would not contribute to the total volume, but will add to the cost; hence rendering useless from a social cost and individual cost point of view.

Remark 3.23. As f is supposed to be a step function on I, it can be reconstructed from the matrix [f J,p ] by:

∀p ∈ P, f p (t) = J∈I J (t)f J,p (4)
Similarly we define discrete link allocation and show that continuous-time link allocation can be reconstructed from the discrete version.

Proposition 3.25.

f e (t) = J∈I J (t)f J,e (6)
Proof.

f e (t) = p∈P:e∈p f p (t) (from Def. 3.2) = p∈P:e∈p J∈I J (t)f J,p (from Eq. ( 4
)
)

f e (t) = J∈I J (t)f J,e (from Def. 3.24)
From above definition and definition of feasible flow, we get:

Proposition 3.26. f is feasible iff: ∀r ∈ R, J∈Ir p∈Pr |J|f J,p = v r (7)
Proof. From Def. 3.

Reduction of Cost Function, Equilibrium and Optimal

In this section, we prove that, under convex cost function, both continuous and discrete formulations are equivalent.

Proposition 3.27. If cost function c J,e (b) is convex, for any f (supposed to be a vector of step functions constant on each J), for all p ∈ P, and for any J in I,

c J,p (f ) = e∈p |J|c J,e (f J,e ) (8)
Proof. This comes from Def. (3.9) applied to step function f and interval J of I and [f J,e ] defined using equations ( 6) and (5).

Proposition 3.28. If cost function c J,e (b) is convex, for any f (supposed to be a vector of step functions constant on each J), for all p ∈ P,

c p (f ) = J∈I e∈p |J|c J,e (f J,e ) (9)
Proof. This is obtained by applying Eq. ( 1) to the step function f and [f J,e ] defined using [f J,p ] as given by Eq. ( 6) and ( 5).

Proof. Same as previous proof but using Eq. (3).

Proposition 3.30. f is an equilibrium flow in (G, R, c) iff:

1. f is a feasible allocation vector, and;

2. for every r ∈ R, for every J, J ∈ I r , c J,p (f ) ≤ c J,p (f ) with p, p in P r where J f J,p > 0.

INRIA Proof. We prove the equivalence in two parts:

⇒ f is supposed to be made of step functions with steps on I. First point is straightforward. By replacing [START_REF] Hall | Multicommodity flows over time: Efficient algorithms and complexity[END_REF] in second point of Def. 3.11 and applying it on J and J in I r we obtain above mentioned conditions.

⇐ First point is again straightforward. Regarding second point, condition 2) in Def. 3.11 is different is the sense that it is for any pair of sub-interval of [t s r , t d r ]. For any such sub-interval and by linearity of the integration over time, cost can be expressed as a sum of cost c J,p with J in I r or part of a such interval. From this, we get two partitions of J and J. They can be different but we can obtain a one to one mapping of sub-intervals of same length from J to J. On each of these pairs of sub-intervals, costs are independent of time and f P is constant. Applying condition 2) of Prop. 3.30 on the underlying intervals in I r and scaling the results with respect to the length of the sub-interval, once summed, proves condition 2) of Def. 3.11.

We summarize this section by the following proposition: Proposition 3.31. f is an optimal allocation iff [f J,p ] is feasible and C(f ) as given in (10) is minimized.

Proof. Trivial.

With this, we have shown the equivalence of the continuous and discrete versions of the routing game under convex cost functions. Next, we move on to the time-expanded network.

Reduction to Non-atomic Routing Game

Since we have defined the RGoT and exhibited the step function structure under convex costs, we can now reduce the game to non-atomic routing game, and then exploit the results already known on this class of games.

Time-expanded Routing Game

To use the results from the classical non-atomic routing game, which we refer as Routing Game in Time-Expanded Network, we now present that the two games can be mapped, and we show how to do so. Definition 4.1 (Time-expanded (TE) network). For a RGoT (G, R, c), Ḡ is the time-expanded graph obtained from G by duplicating it one for each interval of I and adding one pair of virtual source and virtual sink for each request r in R. These extra sources/sinks are connected to the actual source/sink of each duplicate of G in time intervals of I r . We note ē = (J, e).

All edges of Ḡ, except those connecting the actual and virtual sources/sinks, are of the form ē = (J, e). Ḡ has |V ||I| + 2|R| vertices and |E||I| + 2 r∈R |I r | edges. Definition 4.2 (Path in TE network). For every request r in R, we define Pr as the set of paths in Ḡ: for every path p in P r and every interval J in I r , Pr contains a path p = (J, p) which contains the same edges in Ḡ as p in G in addition to the two extra edges that connect the virtual source and virtual sink of r to the nodes in Ḡ, which represents, for each interval their actual source/destination in G. We let, P = r∈R Pr Observe that Pr has |P r ||I r | paths.

To illustrate, the time-expanded graph of requests presented in Fig. 1 is as in Fig. 2. The network of Fig. 1(b) can be seen replicated in Fig. 2 for the intervals of a request. If two requests (e.g., r 1 and r 2 ) have overlapping interval (I 3 ), the graph is shared by paths of the same requests ( P1 and P2 ). Virtual source/sink nodes connect paths of one request. The edges of Ḡ that don't have their counterparts in G are only used by one request and have a null cost in the time-expanded routing games. Hence, they won't contribute to the social cost.

We now define the requests that are used in the new graph: R, the allocations and feasible allocations in these games.

f ē = f Je,e
where f Je,e is the one defined in Def. ⇐⇒ f is feasible

In the following two sections, we map the cost functions, equilibrium and optimal allocation from RGoT to routing games on the TE graph defined in this section.

Cost Function in TE Network

Here we show the mapping of cost between RGoT and time-expanded routing games. Definition 4.11 (Costs). For all ē = (J, e), c ē(.) |J|c J,e (.). We call c the function that associates c ē(.) to ē. For any f and for all p, c p( f ) ē∈p c ē(f ē).

Proposition 4.12. For all p in P,

c p (f ) = p∈I×{p} c p( f ) f p (t)
is still supposed to be null out of I r for each p ∈ P r .

Proof. Transition between first and second line is obtained by recalling that per-edge cost over edges that are not of the form ē = (J, e) is null.

c p (f ) =
c p( f ) = c J,p (f ) (from Prop. 3.27) Definition 4.14 (Social cost). C( f ) = ē f ēc ē(f ē) Proposition 4.15. C( f ) = C(f ) Proof. C(f ) =

Equilibrium and Optimal Allocations

The two definitions given here are from non-atomic routing games. The two propositions basically state that our RGoT can be solved by recasting it to this well-known game and using established results, as will be shown in next section. Definition 4.16. f is an equilibrium flow in ( Ḡ, R, c) if:

1. f is a feasible allocation vector, and; 2. for every r ∈ R, c p( f ) ≤ c p( f ) with p, p in Pr where f p > 0.

Proposition 4.17. f is an equilibrium flow in ( Ḡ, R, c) iff f is an equilibrium flow in (G, R, c).

Proof. First point comes from Prop. 4.10. Second point is obtained by replacing c J,p in Prop. 3.30 as established in Prop. 4.13. Definition 4.18. A feasible flow of ( Ḡ, R, c) is optimal if it minimizes C( f ) over other feasible flows.

Proposition 4.19. A feasible flow f of ( Ḡ, R, c) is optimal iff f is optimal in (G, R, c).
Proof. The proof comes from the equivalence established in propositions 4.10 and 4.15.

Results -Existence and PoA

Next theorem on the existence of equilibrium flows for non-atomic routing games comes from the literature [START_REF] Nisan | Algorithmic Game Theory[END_REF].

Theorem 4.20 ( [5, Theorem 18.8]). Under the assumption of existence of feasible allocation(s), for a non-atomic instance ( Ḡ, R, c) there exists at least one equilibrium, and for any two equilibria f and f , cē ( fe ) = cē ( f e ).

From Theorem 4.20, it follows that, C( f ) = C( f ); i.e., all equilibrium have the same social cost. We also derive from the theorem that our routing game over time have equilibria since the time-expansion of the problem always exists and this always has at least one equilibrium.

Corollary 4.21. From Theorem 4.20 and reduction of the game made in previous sections (game, costs, equilibrium and optimal), it follows that (G, R, c) has at least one equilibrium, and they all have the same social cost.

Existence of optimal allocation is trivial, as there is always some (possibly one) feasible allocation.

Definition 4.22 (PoA). Price of Anarchy (PoA) is the ratio of the cost of worst equilibrium to the cost of optimal allocation.

Since all equilibria in routing game have the same cost, PoA is the ratio of the cost of any equilibrium to the optimal social cost. Furthermore, since cost and social cost are same for RGoT and its TE counterpart, PoA of ( Ḡ, R, c) is equal to the PoA of the corresponding (G, R, c).

Application to Scheduling

In this section, we illustrate the mapping of PoA between the two forms of game, and its implications, for example, on the benefit of doing centralized scheduling to transfer data. We consider a simple example of a single link network with time-slot preference and congestion. First, we describe the cost functions associated with this model, and then discuss equilibrium flows, optimal allocations and PoA.

Lateness and Time-slot Preference as Cost

In this model, f p (and thus f e ) is the expected rate that players plan to have. As long as the link is not oversubscribed, they won't be late with respect to their deadlines. If it is congested, lack of bandwidth increases lateness -cost -as a linear function of this over-subscription. Figures 3(b) and 3(a) shows the cost model for congestion. Another component of this cost represents users' preference of one time-slot over the other. We illustrate this cost on the following example. This RGoT uses one single link and one request of total volume v, but the time interval [0, T ], where 0 < T , is divided in two equal parts and is assigned user preferences: Using time-expansion as defined earlier, we get the routing game with two links (1 and 2) sharing same source and destination; first of these links has cost c1 , and second c2 as defined presently: with f e (t) = t∈[0,T /2) f1 + t∈[T /2,T ] f2 in RGoT. Fig. 4(a) show the cost functions on the timeexpanded network. It includes time preference and congestion aversion for the two links of Ḡ. Note that, the slope s equals d 2 -d 1 v-u . Let r = v/T . If the cheapest link is not congested, i.e., r ≤ u, as its cost is less than that of the other link, both equilibrium flow and optimal allocation exclusively use this path/link. This results in a PoA of 1 as social cost of both are the same. But when the cheaper link starts to get congested, it might remain less expensive than the second link up to a point. Beyond this point, it is less expensive from a social cost point of view to route the flow on the two links. But the equilibrium flow still uses only one link as its cost c p (.) remains lesser than the other. Fig. 4(b) illustrates this situation with two alternative feasible flow and their social cost as shaded. The sum of the areas of the two light gray rectangle are less than the one of the dark rectangle (which is the social cost of equilibrium). More formally: Proposition 5.1. For max{u, v -u} < r and r < min{v, 2u}, equilibrium flow uses exclusively link 1 and P oA > 1.

cē (f ē) = d ē if f ē ∈ [0, u] d ē + s(f ē -u) else cē (f ē) u d 2 d 1 v f ē (
Proof. Provided r is between max{u, v -u} and min{v, 2u}, equilibrium flows only use link 1 and the social cost of such an allocation is:

C eq = (d 1 + (r -u)s)r
If we consider a different feasible allocation which allocates only u on first link and the remaining on second link, its social cost is:

C = u d 1 + (r -u)d 2
We now demonstrate that this allocation has a social cost strictly less than the equilibrium flow:

C < C eq ⇐⇒ u d 1 + (r -u) d 2 < (d 1 + (r -u)s)r ⇐⇒ (r -u) d 2 < (r -u) d 1 + (r -u)s r ⇐⇒ d 2 -d 1 s < r ⇐⇒ v -u < r • u ≤ r ≤ 2u: x * 1 (r) = max{min{r, v/2}, u} and x * 2 (r) = r -max{min{r, v/2}, u} • 2u ≤ r: x * 1 (r) = max{min{r, (r + (v -u)/2))/2}, u} and x * 2 (r) = r -max{min{r, (r + (v -u)/2))/2}}
The optimal allocations as a function of r are shown on Fig. 5 by the emphasized lines (plain and dashed) for different values of v/2.

Equilibrium flows Regarding equilibrium flows:

• r < v: Since c 1 (r) < c 2 (r), x 1 (r) = r and x 2 (r) = 0. • v ≤ r < u + v: x 1 (r) = v and x 2 (r) = r -v is the equilibrium flow. • u + v ≤ r: x 1 (r) = (v -u)/2 + r/2 and x 2 (r) = (u -v)/2 + r/2 is the equilibrium flow.
PoA For the specific case where u < v/2 < 2u, we have:

• r ≤ v/2: since x * 1 (r) = r, P oA = 1. • v/2 ≤ r ≤ u + v/2: in this region, x * 1 (r) = v/2 and C * = (d 1 -d 2 + (v/2 -u)s) v/2 + r d 2 while x 1 (r) = r and C eq = r (d 1 + (r -u)s). P oA = C eq C * = r (d 1 + (r -u)s) (d 1 -d 2 + (v/2 -u)s)v/2 + r d 2 = r (d 1 + (r -u)s) r d 2 -(v/2) 2 s
PoA is increasing as long as s < (2r +v)d/v 2 since r is less than v if s is less than 2(u+v)d/v 2 maximum PoA is reached for r = u + v/2 and is:

P oA = (2u + v)(-2d + s v -2s u) -4d u -2d v + s v 2
In the other case, PoA will first increase and then decrease before leaving the interval. Maximum PoA is reached for r = s v 2 /d-v 2 and its value is:

P oA = (s v -d) 2 d 2
• u + v/2 < r: Since v/2 > u, r > 2u and then the system is clearly overloaded (both links congested) which is not really interesting in this scenario.

Conclusion

In this example, we see that as long as none of the links are congested, PoA remains 1. This basically means that selfish behavior leads to optimal allocation. Therefore, under the presented model no centralized control is needed to achieve optimal allocation. This model, following Wardrop equilibrium concept, assumes that users have complete and up-to-date information to decide which path to use based on the cost. The cost depends on the link allocation reached so far. Thus, it assumes that users know the current link allocation. This information is not needed for the time-slot preference part of the cost, but for the congestion aversion part. Hence, implementing this scheduling mechanism without coordinator mandates the congestion information be available to users, that too, in advance.

In both cases coordination is essential: either to take the decision, or to publish the information.

In the case where the system is near congestion, as studied above, PoA is greater than 1. Here, coordination can improve the social cost by sacrificing some of the flow. This can not be achieved by selfish users unwilling to spontaneously take a more costly path.

From [START_REF] Roughgarden | The price of anarchy is independent of the network topology[END_REF], we know that using affine cost functions it is possible to construct instance of the nonatomic routing game with PoA up to 4/3. This is value is attained in two-link settings. When the class of function quadratic functions, this bound becomes 3 √ 3/(3 √ 3 -2). Actually, the more non-linear the cost functions are, the higher is the worst PoA. Using polynomial cost functions of degree p, PoA tends to infinity as p does. In the next section, we see how this inefficiency can be mitigated through coordination or modification of the cost function.

Reducing Social Cost

Since PoA can be large, it is worth using optimal allocation to minimize social cost. Minimization of social cost can be achieved by a centralized system using optimal allocation or the game can be modified so that equilibrium of the new game is an optimal allocation for the original game. This corresponds to a reduction of the PoA in the new instance.

Cost Sharing and Resource Coordination Consider a resources manager (RM) who proposes some commodities for a cost which depends on the total demand. Users know only about their own individual decisions. The cost function is likely to be convex, positive and non-decreasing as it ultimately results from the cost of congestion. A resource coordinator (RC) aggregating the requests and demanding the optimal allocation for the aggregate can do better than selfish user behaviors. It follows that RC can improve the social cost, i.e., total cost. RC can then share the benefit among all entities: RC and clients. This makes this configuration profitable for both clients and RC.

Pigouvian Tax Another alternative is to modify the game so that cost perceived by users deter them from ultimately reaching an allocation other than the optimal of original game. For this purpose the related theorem is: Theorem 5.2 (Theorem 18.27 of [START_REF] Nisan | Algorithmic Game Theory[END_REF]). Let f * be an optimal flow for ( Ḡ, R, c) and let τ ē = f * ē c′ ē( f * ē ) denote the marginal cost tax for edge ē with respect to f * . Then f * is an equilibrium flow for ( Ḡ, R, c + τ ) It basically says that, equilibrium of the modified game is optimal allocation of the original game. From the realization point of view, this solution still has the problem of sharing the information, and also in addition making the users value the modified cost.

Concluding Remarks

This paper presented a new model for RGoT. It shows that under convex cost this can be solved using the results from non-atomic routing games. Based on an example, we have seen, selfish user behavior increases the social cost. Just because they prefer one time interval to the other, users selfishly raise the cost of cheapest links instead of sacrificing to keep it and social cost low. This motivates the need to coordinate. Possible extension of this work is to use atomic or splittable flow model instead of non-atomic model in order to consider the perception users have of their own impact.

INRIA

Definition 3 . 5 (

 35 RGoT). (G, R, c), where G is a network, R a set of requests and c a cost function, defines an RGoT.

1

 1 Time c J,. (.) c J,. (.) (a) Timeline and its discrete versions. (b) Example of requests with their alternate paths.

Figure 1 :

 1 Figure 1: Example of requests with their intervals.

Definition 3 . 6 (

 36 Time intervals). For a game (G, R, c) with I c being the set of time interval used to define cost function c, we define T as the set of time intervals defined by all the start time t s r and deadline t d r of request i ∈ R and start and end of intervals I c included in T. T is a partition of T. Definition 3.7 (Same length time intervals). We define I as the refined set of intervals of T such that all intervals of I have the same length (called |J| in the remaining). This set is obtained by subdividing I. Definition 3.8 (Intervals of a request). I r {J ∈ I|J ⊂ [t s r , t d r ]} In the remainder of this work, c J,e has been extended to J in I by using the function c I,e of the interval I ∈ I c that contains J. Definition 3.9 (Per-path cost on an interval). For any time interval J, c J,p (f ) t∈J e∈p c(e, f e (t), t)dt Definition 3.10 (Per-path cost).

  c p (f ) c T,p (f ) = t∈T e∈p c(e, f e (t), t)dt (Def. 3.9) = t∈T e∈p J∈Ic 1 t∈J c J,e (f e (t))dt (def. of c())

Definition 3 .

 3 15 (Social cost).

  t e f e (t)c(e, f e (t), t)dt = t e J∈Ic

1

  t∈J f e (t)c J,e (f e (t))dt (def. of c())=J∈Ic t∈J e f e (t)c J,e (f e (t))dtC(f ) = J∈I e t∈J f e (t)c J,e (f e (t))dt(3)Optimal allocationin the social sense -is: Definition 3.16 (Optimal allocation). A feasible flow of (G, R, c) is optimal if it minimizes C(f ) over other feasible flows.

Proposition 3 .

 3 19. If cost function c J,e (b) is convex, positive and non-decreasing, for any feasible allocation f , there is one f made of step functions which: (1) has a costs C( f ) not worse than C(f ); and (2) has steps on I.

Definition 3 .

 3 24 (f J,e ).f J,e p∈P:e∈p f J,p

Proposition 3 .

 3 29. If cost function c J,e (b) is convex, for any f (supposed to be a vector of step functions constant on each J), C(f ) = J∈I e |J|f J,e c J,e (f J,e )

2 Figure 2 :

 22 Figure 2: Time-expanded network.

Proposition 4 . 3 .

 43 ē = (J e , e) ∈ p = (J p , p) iff e ∈ p J e = J p Proof. By construction of Ḡ in Def. 4.1.

Definition 4 . 4 (

 44 Requests in the TE network). For every request r in R, we define: r (s r , dr , Pr , vr vr |J| ) and the set R as the set of requests r.

Definition 4 . 5 (Remark 4 . 6 .

 4546 Allocations in the TE network). Allocations in the TEnetwork: f [f p] [f J,p ]Link allocation: f ē p∈ P:ē∈p f p For each edge of Ḡ which doesn't have a ē = (J, e) form, there is no need to define f ē. In the remaining, we use the notation ē = (J, e) to refer to edges that contribute to the costs and define the counterpart in G and time intervals.

Proposition 4 . 7 .

 47 For any ē = (J e , e), f ē = f Je,e . Proof. By Def. 4.5, f ē = p∈ P:ē∈p f p and by Def. 4.5, f p = f J,p . Using Prop. 4.3 and definition of p and ē, we get: for any ē = (J e , e),

  J∈I e∈p |J|c J,e (f J,e ) (from Prop. 3.28) = J∈I ē∈p c ē(f ē) (from Def. 4.5, 4.11 and Prop. 4.3) c p (f ) = p∈I×{p} c p( f ) (from Def. 4.11) Proposition 4.13. For all p in P and J in I, c p( f ) = c J,p (f ) with p = (J, p).

  f ē) (from Def. 4.11) = e∈p |J|c J,e (f J,e ) (from Def. 4.5, 4.11 and Prop. 4.3)

  e c J,e (f J,e ) (by Prop. 3.29) = ē f ēc ē(f ē) (by Prop. 4.7 and Def. 4.11) C(f ) = C( f ) (by Def. 4.14)

  Cost model for congestion on a link of capacity ue. c e (f e ) u e f e (b) Cost of congestion ce(fe) for edge e.

Figure 3 :

 3 Figure 3: Cost model for congestion.

d 1 and d 2

 2 with 0 ≤ d 1 < d 2 . The cost function adopted for this link of capacity u is then:c(e, f e , t) = d 1 + c e (f e ) if t < T /2 d 2 + c e (f e ) else with: c e (f e ) = 0 if f e ∈ [0, u] fe-u ue.α(e,R)else Let s = 1 ue.α(e,R) , s > 0. For each time interval, the cost is convex, positive and non-decreasing in link allocation.

Figure 4 :

 4 Figure 4: Time preference and congestion aversion.

Table 1 :

 1 Table of notations.

INRIA

  Definition 3.21 (f J,p ). For a feasible allocation f , we define the |I| × |P| matrix [f J,p ] (J,p)∈I×P , Remark 3.22. It can be noted that in Def. 3.21, f J,p for J not used by a request r using path p is supposed to be null. More formally:

	where f J,p interval J.	t∈J	f p (t) |J|	dt is the constant rate of a step function which transfers as much as f on

  3, f is feasible iff:

		t d r
	∀r,	r p∈Pr t s	f p (t)dt = v r
	⇐⇒ ∀r,	t d r r p∈Pr J∈I t s	J (t)f J,p dt = v r (from Remark 3.23)
			t d r
	⇐⇒ ∀r,	p∈Pr J∈I	r t s	J (t)f J,p dt = v r
	⇐⇒ ∀r,	p∈Pr J∈Ir	|J|f J,p = v r (from Remark 3.22)

  3.24.

	Remark 4.8. Using Prop. 4.7 and Prop. 3.25, we get:
	f e (t) =
	Proof.
	f is feasible
	⇐⇒ ∀r ∈ R,	J∈Ir p∈Pr	|J|f J,p = v r (from Prop. 3.26)
	⇐⇒ ∀r ∈ R,	(J,p)∈Ir×Pr	f J,p =	v r |J|	(as |J| is uniform)
	⇐⇒ ∀r ∈ R,	p∈ Pr	f

ē=(J,e)∈I×{e} J (t)f ē Using Remark 3.23 and Def. 4.5, we also have:

f p (t) = p=(J,p)∈I×{p} J (t)f p

Flow allocation f can thus be reconstructed. Definition 4.9 (Feasible allocation). f is a feasible allocation for R iff: ∀r ∈ R, p∈ Pr f p = vr This is the definition of feasibility in normal routing game. Proposition 4.10. f is feasible iff f is feasible. p = vr (Def. of Ḡ, R and f p)

RR n°6931

Acknowledgement

This work was done in the framework of the INRIA and Alcatel-Lucent Bell Labs Joint Research Lab on Self Organized Networks.

This text is also

which holds by assumption.

As the cost of this feasible allocation is strictly less than the one of equilibrium flow, PoA is strictly more than 1.

Equilibrium and Optimal Allocation

Define x 1 (r) and x 2 (r) as the part of the traffic r sent on links 1 and 2 for equilibrium flow. Similarly, x * 1 (r) and x * 2 (r) are the fraction for optimal allocation. We have x 1 (r) + x 2 (r) = r and x * 1 (r) + x * 2 (r) = r. C eq (r) denotes the social cost of equilibrium under demand r and C * (r) the social cost of optimal allocation.

Optimal allocations Since costs are defined in pieces, social cost are also piecewise. In the following, x denotes x 1 and x 2 is r -x.

It follows that the social cost of a feasible allocation C(x, r) = x c 1 (x) + (r -x)c 2 (r -x) has the following expression in the different regions (a),(b),(c) and (d) of Fig. 5:

, which is strictly positive when x > v/2, and strictly negative when x < v/2.

))/2 and strictly positive above.

As cost functions are continuous, so are optimal and equilibrium allocations as functions of r.

We conclude from previous list that optimal allocations are:

• r ≤ u: x * 1 (r) = r and x * 2 (r) = 0 INRIA


Ngoc Long Nguyen

Jérémy Anger

Axel Davy

Pablo Arias

Gabriele Facciolo

SELF-SUPERVISED PUSH-FRAME SUPER-RESOLUTION WITH DETAIL-PRESERVING CONTROL AND OUTLIER DETECTION

Keywords: push-frame burst satellite imaging, multiimage super-resolution, deep-learning, self-supervised learning, robust estimators

Self-supervised training enables the application of deeplearning based methods for multi-image super-resolution of satellite imagery. In this work we propose two improvements on the self-supervised Deep-Shift-and-Add (DSA) method introduced by Nguyen et al. First, we demonstrate how the self-supervised loss of DSA can be extended to provide the image interpreter with a spatially varying parameter to control the trade-off between detail preservation and noise removal at test time. Second, we endow the DSA architecture with a mechanism that enables the network to be robust to outliers produced for example by dead pixels, reflections or registration errors.

INTRODUCTION

Multi-image super-resolution (MISR) by numerical methods has recently been adopted as way to increase the resolution of push-frame satellites [START_REF] Murthy | Skysat-1: very high-resolution imagery from a small satellite[END_REF][START_REF] Anger | Fast and accurate multi-frame super-resolution of satellite images[END_REF]. By leveraging high framerate low-resolution acquisitions, low-cost constellations can be effective competitors to traditional high-cost satellites.

Satellite images are usually noisy. This noise is often considered as an unwanted artifact and is subject to removal. In doing so, it is inevitable that some details will also be lost. In this work, our goal is to perform robust joint super-resolution and denoising from a burst of satellite images, with control over the trade-off between noise removal and detail preservation. We focus on push-frame satellite sensors such as the SkySat constellation from Planet.

Most MISR methods for satellite images are still based on classic model-based techniques [START_REF] Murthy | Skysat-1: very high-resolution imagery from a small satellite[END_REF][START_REF] Anger | Fast and accurate multi-frame super-resolution of satellite images[END_REF][START_REF] Anger | Parallax estimation for push-frame satellite imagery: application to super-resolution and 3d surface modeling from skysat products[END_REF]. For example, This work was supported by a grant from Région Île-de-France. It was also partly financed by IDEX Paris-Saclay IDI 2016, ANR-11-IDEX-0003-02, Office of Naval research grant N00014-17-1-2552, DGA Astrid project « filmer la Terre » n o ANR-17-ASTR-0013-01, MENRT. This work was performed using HPC resources from GENCI-IDRIS (grant 2021-AD011012458) and from the "Mésocentre" computing center of Centrale-Supélec and ENS Paris-Saclay supported by CNRS and Région Île-de-France (http://mesocentre.centralesupelec.fr/). We thank Planet for providing the L1A SkySat images. Anger et al. [START_REF] Anger | Fast and accurate multi-frame super-resolution of satellite images[END_REF] propose solving a least-squares problem that fits spline polynomials to the observed samples (we denote it ACT, as its early trigonometric polynomial formulation [START_REF] Feichtinger | Efficient numerical methods in non-uniform sampling theory[END_REF]).

In [START_REF] Farsiu | Robust shift and add approach to superresolution[END_REF], Farsiu et al. propose a variational method that extends the classic shift-and-add MISR to be robust to outliers, which amounts to applying pixel-wise medians.

On the other hand, training supervised deep-learning based MISR methods is challenging because it is difficult to obtain real training satellite data with ground truth. To overcome this problem, Nguyen et al. [START_REF] Nguyen | Self-supervised multi-image super-resolution for push-frame satellite images[END_REF] propose Deep Shiftand-Add (DSA), a neural network for MISR from satellite image bursts that can be trained in a self-supervised manner, i.e. it does not require ground truth HR data. This work was recently extended in [START_REF] Nguyen | Self-supervised super-resolution for multi-exposure push-frame satellites[END_REF] to support multi-exposure bursts. Notwithstanding its advantages, DSA does not handle outliers. For satellite imagery, moving objects misregistration, reflections or dead pixels are examples of such outliers.

Besides, like most deep-learning based image restoration algorithms, DSA tends to smooth out textures or details that have a low contrast relative to the noise level. This is a known problem for restoration methods based on minimizing a distortion measure (such as the MSE or L 1 loss) [START_REF] Blau | The perception-distortion tradeoff[END_REF]. As these results have a bad perceptual quality, several works combine distortion losses with adversarial losses that aim at reducing the distance between the distribution of restored images and that of real HR images [START_REF] Blau | The perception-distortion tradeoff[END_REF][START_REF] Ledig | Photo-realistic single image super-resolution using a generative adversarial network[END_REF]. However, this requires the net- work to "invent" plausible information in the regions where the original content cannot be recovered. This behavior is desired for applications where the goal is to create an aesthetically pleasing natural looking image, but it is unacceptable in cases where critical decisions are made based on the data.

In this paper we propose two contributions to address these limitations of DSA.

1. A data-fitting term that allows controlling the amount of detail in the solution with a spatially varying map that is given to the network as input at test time. In this way a user can control this noise reduction and detail preservation trade-off depending on the application. This is inspired by a common photography trick for recovering low contrast details lost to image denoising, which consists in adding back to the output a fraction of the noisy input image.

2.

We also propose a robust version of DSA, thanks to the introduction of outlier masks computed by the network itself (see Fig. 1). Our proposed framework contains DSA, ACT and ACT-robust (a robust variant of the ACT method formulated as an L1 fitting) as particular cases (without the need to solve a computationally complex optimization problem at test time) and yields state-of-the-art results.

PROPOSED METHOD

Our network (Fig. 2) is built upon the Deep Shift-and-Add (DSA) architecture [START_REF] Nguyen | Self-supervised multi-image super-resolution for push-frame satellite images[END_REF] with four major modules: Motion Estimator, Encoder, Feature Shift-and-Add block (FS&A) and Decoder. The self-supervised DSA loss [START_REF] Nguyen | Self-supervised multi-image super-resolution for push-frame satellite images[END_REF] drives the network to produce a super-resolved image I HR such that when subsampled, it coincides with the reference frame

I LR 0 DSA loss = Π 2 ( I HR ) -I LR 0 1 , (1)
where I HR is the network output and Π 2 is the subsampling operator. Since the reference frame is withheld from fusion during training, the network cannot learn to reproduce the noise in the reference. Thus the training converges to produce a noise-free high-resolution images. This self-supervised loss is based on the minimization of a distortion measure with respect to a target. It is a known fact that these type of losses tend to smooth fine details whose magnitude is comparable with that of the noise [START_REF] Blau | The perception-distortion tradeoff[END_REF].

We propose a loss that permits to control the trade-off between noise reduction and detail preservation. For that we incorporate a multi-frame data fitting term controlled by a spatial map D (Sec. 2.1). In addition we introduce in the DSA network architecture MaskNet, a new trainable module (Fig. 2), whose purpose is to produce outlier masks O that indicate the presence of outliers (Sec. 2.2).

Noise reduction -detail preservation trade-off

The self-supervised DSA loss imposes a data-driven prior that favors smooth reconstruction in regions with lower contrast. The cost of producing a noiseless image is that some details might be lost. To counteract this, we add the following loss, which corresponds to ACT [START_REF] Anger | Fast and accurate multi-frame super-resolution of satellite images[END_REF][START_REF] Feichtinger | Efficient numerical methods in non-uniform sampling theory[END_REF] when p = 2:

LL p loss = 1 T T i=1 Π 2 Warp( I HR , 2F i→0 ) -I LR i p p , (2)
with 1 ≤ p ≤ 2 (in this work we only consider p = 1 and p = 2) and where Warp is an operator that warps its input according to the estimated motion field F i→0 (see [START_REF] Nguyen | Self-supervised multi-image super-resolution for push-frame satellite images[END_REF] for details). This loss corresponds to the likelihood of the data under a generalized Gaussian noise model [START_REF] Farsiu | Robust shift and add approach to superresolution[END_REF].

To understand the rationale behind this loss, consider a case in which the images can be aligned with an integer shift (so that no interpolation is needed). Then the minimizer of the LL p loss is obtained by aligning and aggregating the LR images on the HR grid. For p = 2 the aggregation is the average, and for p = 1 is the median. For Gaussian noise, these solutions are unbiased estimators of the high-resolution image. Thus no details are lost and the noise is reduced via the temporal aggregation. Of course some noise will remain, with variance depending on the number of images in the sequence. Note that these solutions do not require any datadriven learning: i.e. they do not depend on any priors learned from the data; they only depend on the set of LR images. This is because the target frames I LR i are all part of the input. Since the least-squares (LS) solution (p = 2) is sensitive to outliers [START_REF] Farsiu | Robust shift and add approach to superresolution[END_REF], we propose using the p = 1, which we call least absolute value (LAV) loss. Complete training loss. Most of the time, our priority is to produce a noise-free HR image. Nevertheless keeping details might be preferred when we have few images or when we want to detect very high-frequency objects such as crosswalks, solar panels, etc. To control the trade-off between noise removal and detail conservation, we introduce a noise-detail map (denoted D) as a parameter to balance the losses ( 1) and ( 2). This map is spatially varying and takes values between 0 and 1. Values closer to 1 indicate that we want to keep details (and noise), whereas small values imply that the corresponding region should be denoised. The training loss is defined as the balance between the DSA and the LAV loss per pixel loss = Π 2 ( I HR ) -

I LR 0 • (1 -Π 2 (D)) 1 + 1 T T i=1 Π 2 (Warp( I HR , 2F i→0 )) -I LR i • Π 2 (D) 1 , (3)
where I HR = Net(I LR i=1,...,T , D) is the network output and "•" denotes the element-wise multiplication. To simplify, we assume that D is smooth and that the images are coarsely prealigned so that the D does not have to be warped in the loss.

Outlier handling

In DSA [START_REF] Nguyen | Self-supervised multi-image super-resolution for push-frame satellite images[END_REF], the features computed by the encoder are averaged by the Feature Shift-and-Add module. Because of this averaging, outliers have a strong impact that the decoder cannot entirely mitigate. For this reason, we propose removing them from the averaging by incorporating a submodule MaskNet to the DSA architecture to predict outlier masks. We take inspiration from a video denoising application [START_REF] Maggioni | Efficient Multi-Stage Video Denoising with Recurrent Spatio-Temporal Fusion[END_REF] where a similar mask predicting network is used for removing misaligned areas in a recursive frame fusion method. We define outliers as regions that are inconsistent with the majority of frames in the sequence, and masks allow to exclude them from fusion. To estimate such masks, we first approximate a low-resolution outlier-free image using a temporal median of the LR frames aligned to the reference, which we denote by M LR . Then the absolute difference [START_REF] Maggioni | Efficient Multi-Stage Video Denoising with Recurrent Spatio-Temporal Fusion[END_REF] between the warped median image and each image is used as input for MaskNet

O i = MaskNet | Warp(M LR , F i→0 ) -I LR i | . (4)
We also impose the smoothness of the produced masks by adding a TV regularization term in the loss. The outlier masks are then used as weights in the weighted FS&A block where ϵ is a threshold to avoid division by 0, {J LR i } are the features computed by the Encoder, {F i→0 } are the optical flows estimated by the Motion Estimator, and the SPMC module [START_REF] Nguyen | Self-supervised multi-image super-resolution for push-frame satellite images[END_REF][START_REF] Tao | Detail-revealing deep video super-resolution[END_REF] maps the LR features onto a common HR grid. The outliers will be assigned negligible weights in the outlier masks so that they do not contribute in the fusion.

J HR = T i=1 SPMC(J LR i • O i , F i→0 ) T i=1 Max(SPMC(O i , F i→0 ), ϵ) , (5)

EXPERIMENTS

In our experiments, we first demonstrate the trade-off between denoising and detail restoration using the noise-detail map. Then we justify our choice of architecture and losses in order to handle outliers.

Examining detail preservation map D

In order to train the network, we prepare sets of training input data {I LR i=0,...,T ; D}. The random spatially-varying noisedetail maps D are generated first by thresholding a filtered Gaussian noise image (σ = 40) (the filter itself is a Gaussian filter with σ = 28), then the resulting binary image is then smoothed with a small Gaussian filter (σ = 3).

Fig. 3 illustrates the trade-off between noise reduction and detail preservation when we change D for the cases of 4 and 14 input frames. As expected, with D = 0 the network removes noise while smoothing out the textures as it cannot distinguish high frequencies from noise. On the other hand, with D = 1 the output is noisier but it better preserves the high frequency details. The difference between the two behaviors is particularly noticeable when using few input frames. Moreover, we can use a spatially varying map to reduce noise in uniform regions and preserve details in textured regions as shown in Fig. 3. Since the map D is a network input provided at test time, it lends itself to applications where a user interactively edits the map.

Robustness to outliers

To train MaskNet we add synthetic outliers in the LR images during training. To this aim, we first generate random blobs in an image and then substitute the pixels of these regions with data from a different stack. As the LL p loss optimization is not data-driven, it is strongly affected by outliers. Here, we justify two key features in our framework that enable its robustness to outliers: The MaskNet and the L1 norm in the LAV loss.

The authors of [START_REF] Farsiu | Robust shift and add approach to superresolution[END_REF] show that the optimal solution of the LS (resp. LAV) problem is the pixelwise average (resp. median) of the LR images. Consequently, the LS problem is not robust to outliers. We experimentally observed (Fig. 4) that using the L2 norm, MaskNet learns to produce only constant maps and the network produces artifacts in the SR result. Conversely, with the L1 norm, the MaskNet is able to detect outliers in the LR images and impose a negligible weight to these regions.

SR on synthetic data with outliers. Table 1 highlights the usefulness of the MaskNet when the image stacks contain outliers. We can notice that with few or many frames, both for D = 0 or D = 1, MaskNet helps to increase significantly the PSNR by 0.3 -1.5dB.

SR on real satellite data with moving objects. Moving objects that are not correctly aligned can be considered as out-liers. Fig. 1 illustrates how our architecture with and without MaskNet handles moving objects. As expected, the motion estimator of L1B [START_REF] Murthy | Skysat-1: very high-resolution imagery from a small satellite[END_REF] and DSA predicts smooth optical flows and ignores small moving objects. Consequently, without MaskNet we observe a blur trait on the highway. On the other hand, when we use MaskNet, the network is able to filter out the motion of the car, leading to a better reconstruction.

CONCLUSION

We extended the self-supervised DSA method [START_REF] Nguyen | Self-supervised multi-image super-resolution for push-frame satellite images[END_REF] by providing a spatially varying parameter to control the trade-off between detail preservation and noise removal at test time. In addition we endow the DSA architecture with a mechanism that enables the network to be robust to outliers produced for example by dead pixels, reflections or registration errors. All within a self-supervised framework. These improvements lead to state-of-the-art results.

In this work we assume that outliers are everything not consistently visible in most of the frames. However, if we want to preserve the content of the reference image some modifications might be necessary. This will be the subject of future work.

Fig. 1 :

 1 Fig. 1: Super-resolution from a sequence of 15 real lowresolution SkySat L1A frames. (a) L1B from Planet, (b) DSA [3], (c) Our improvement with an additional CNN to detect the outliers.

Fig. 2 :

 2 Fig.2: Overview of our method, which builds upon the Deep Shift-and-Add (DSA) architecture[START_REF] Nguyen | Self-supervised multi-image super-resolution for push-frame satellite images[END_REF].

Fig. 3 :

 3 Fig. 3: Super-resolution from a stack of 4 (first row) and 14 (second row) noisy synthetic images. From left to right: Reconstruction with D = 0 (without detail preservation), with D = 1, and with a mixed D map.

6 LS

 6 loss, with Mask LAV loss, w/o Mask LAV loss, with Mask

Fig. 4 :

 4 Fig. 4: Effect of the architecture and the loss in robustness to outliers. First line: 7 LR with outliers. Second line: Reconstruction with the LS or LAV loss, with and without MaskNet.

Table 1 :

 1 Average PSNR on the synthetic test set with outliers.

		D = 1	D = 0	(D = 1) + Mask (D = 0) + Mask
	T = 4	31.10	36.87	32.65	37.16
	T = 14	36.52	40.45	37.25	40.77