Houssam Eddine Zahaf

Richard Olejnik

Giuseppe Lipari

Abou El

Hassen Benyamina

Modelling the Energy Consumption of Soft Real-Time Tasks on Heterogeneous Computing Architectures

The problem of reducing the energy consumption of embedded processors is of paramount importance for mobile devices powered by batteries. Many of these devices embed heterogeneous processors like the ARM bitLITTLE for adapting to the performance requirements and energy consumption profiles of modern applications. However, before proposing algorithms for optimally managing the energy consumption of a device, it is necessary to build a realistic model of the performance and power profile of an application.

In this paper we address this problem by proposing a simple model for the execution time of soft real-time applications and the energy consumption of the hardware platform. We also propose a simple benchmarking methodology for obtaining the parameters of the model from measures. We identify in the memory access pattern of software tasks a key factor that influences both the performance and the energy consumption. We then show how to apply our methodology on the ODROID development board, and we present some preliminary results of our study. We conclude the paper by discussing current research and future directions.

INTRODUCTION

Many modern electronics embedded devices are powered by batteries, hence it is of paramount importance to reduce their power consumption as much it is possible, so to prolong their autonomy.

There are two contrasting objectives: reducing the energy consumption of the device (for lowering the "cost" of recharging but also for reducing heating and prolonging the autonomy of battery powered embedded systems), while at the same time optimising the Quality of Service provided to the user. Many of the applications running on model smart-ACM ISBN 978-1-4503-2138-9. DOI: 10.1145/1235 phones and tablets have soft real-time requirements, because the quality of service experienced by the users depends on the response time of the computations. For example, users want to watch audio/video content streamed from Internet at the maximum quality on their tablets while at the same time prolonging the duration of the battery. A delayed visualisation of some video frames, or even worse, the delayed decoding of an audio sample, may greatly compromise the user experience.

In order to propose efficient techniques for optimising the energy consumption of a complex device consisting of heterogeneous computing cores without compromising the response time of the running applications, it is of foremost importance to build a model of the system. Therefore, we need to identify the main factors that impact on the energy consumption and on the response time of an application.

The complexity of the problem forces us to make a choice between accuracy and complexity. An exact model of the hardware/software architecture would have the same complexity as the hardware/software architecture itself, and so it would be of no use in practice. We need then to abstract away the many details and focus on the main parameters that can influence energy consumption and execution time.

This paper is a first step toward a model of the energy consumption and of the execution time of a task running on a ARM big/LITTLE architecture. For the hardware architecture, our model considers a few parameters as the type of the processor being used (big or LITTLE), the operating frequency, and the memory architecture (size of the cache). For the execution time of a soft real-time task we consider the number of cache misses, the processor on which it is running and the operating frequency. We show how to compute the model parameters with a set of measurement, and we show that our model captures the behaviour of the application in a sufficiently precise manner.

The proposed model will be used in the future to build scheduling and allocation algorithms that aim at minimising energy consumption without compromising on the Quality of Service.

RELATED WORK

In the research literature, many works have focused on measuring execution time of executing tasks, both in the real-time scheduling domain and in High Performance Com-puting (HPC). However, not so many papers focus on benchmarking both execution time and energy consumption in mobile embedded systems.

Erich et al., in [START_REF] Strohmaier | Architecture independent performance characterization and benchmarking for scientific applications[END_REF], proposed an approach for designing benchmarks for scientific computing. They focused basically on task granularity, and on temporal and spatial locality. The main result of their investigation is that one of the most important factors influencing task performances is the memory access. However, they only addressed big processors (like Intel Xeon Phi, . . .) designed for HPC, and did not consider low-power embedded processors, nor heterogeneous architectures. In particular, here we are interested in the trade-off between power consumption and performance.

Kambel et al. [START_REF] Kimball Dunn | Development of smart phone usability benchmarking tasks[END_REF] investigate benchmarking for mobile phones. They assume that classic benchmarking methodologies are not suitable for the target devices. However, their work focuses on a subset of applications, notably application heavily relying on GUI interaction, and low-intensive computation applications (e-mail clients, chats, etc.), where the emphasis is on networking.

Gal-On and Levy [START_REF] Gal-On | Measuring multicore performance[END_REF] discussed the effectiveness of current benchmarks for modern processors. They focus on parallelism, heterogeneity, memory access impacts and on performances. The report does not present any experiments and discusses in a general way the possible impact of these aspects.

In the area of Real-Time scheduling, benchmarking is focused on worst-case execution time rather then on performance and energy consumption. In addition, benchmarks are used to compare scheduling and worst-case execution time algorithms, rather than on build a model of the task behaviour. Here instead we are interested in average execution time and its relationship to average memory access patterns, and how this influence the power consumption of the architecture.

SYSTEM MODEL AND ASSUMPTIONS

The objective of this paper is to build a model of the execution time and of the energy consumption of soft real-time task executing on a heterogeneous multi-core architecture. In this section we describe the model of a soft real-time task and of the hardware architecture.

A soft real-time task is a software process executing periodically. An example of C code that uses the POSIX RT API for implementing a periodic real-time task is reported in Figure 1.

From a mathematical point of view, a periodic real-time task is a sequence of jobs Ji = (ai, ci, di). Each job corresponds to the execution of the code denoted as periodic code in the figure. The arrival time ai is the activation time of the job (corresponding to the values of variable next); the computation time ci corresponds to the execution time of the jobs, that is the processing time necessary to complete the periodic code; di is the job's deadline that is the absolute time within which the job should complete its execution in order to satisfy the performance requirements of the user.

Per a periodic task, we have that ai = ai-1 + T , where T is the task's period. Very often the job deadline is equal to the next arrival time of the job (each job should complete before the next one is activated). Hence di = ai + Ti = ai+1.

The execution time of a job depends on many factors: 1) the number of hardware instructions to be executed, 2) how many clock cycles per each instruction 3) the operating In turn, 1) depends on the input data and on the internal state of the task; 2) depends on the hardware architecture and its internal state. Notably, the pipeline state and the cache state may have a great influence on the execution time of the task. Typically, we can control the processor operating frequency, to a certain extent. In particular, it is possible to set the frequency so to reduce energy consumption without degrading too much the performance of the task. For example, under the assumption that the system executes one single task, and that we have an exact model of the execution time of a job, we could set the processor frequency so that ci = di -ai = T , that is the task finishes always exactly at its deadline.

This technique is called Dynamic Voltage and Frequency Scaling (DVFS), and will be described in the following. Unfortunately, it is impossible to have a precise model of the execution time of a task, so we have to derive approximate models.

In the following, we restrict our attention to tasks with a weak dependence from input and state variables; in other words, we assume that the number of processor instructions to execute in each jobs is approximately constant. The variability due to input data and state variables is very dependent on the application class, and it will be the focus of future works.

In the next section we describe the hardware architecture used for our experiments, and a set of benchmarks for measuring execution time and power consumption of a specific task.

Odroid architecture description

The ODROID XU3 [START_REF]Odroid xu3 board description[END_REF] The frequency of little cores can be calibrated homogeneously from 200Mhz to 1400Mhz in discrete steps of 100Mhz (13 modes). Respectively, the frequency of big cores can be calibrated from 200Mhz to 2000Mhz, in discrete steps of 100Mhz (19 modes). In this paper, core operational frequency and state are set offline and do not change dynamically. Therefore, for the moment we are not interested in the core state and frequency changing costs.

Odroid bench design

In order to build a simple task and architecture models, we implemented a simple bench of periodic matrix multiplication of size S × S. The matrix multiplication example allows us to change the data size and can be easily parallelized. Also, the matrix multiplication implementation allows us to control the spatial and temporal locality which make the cache behaviour easily predictable. To measure the power consumption, we used 3 out of the 4 available embedded sensors: the big cores sensor; the little cores sensor; the memory sensor. The external power sensor allows us to have an external view on the power consumption of all components of the ODROID board.

The thread has a higher priority than all other threads in the system. In Figures 3,5, 6, 7, 8, 10, the task is implemented in one periodic thread locked on a big core or a little core using POSIX thread affinity. The frequencies are set using the CPUFreq tool. Each experiment is run 500 times, the presented results are the average values. After each experimentation, L1 and L2 caches are reset.

Benchmark results

Figure 3 presents the average execution time of the thread locked on a big core (red squares) and on a little core (blue circles) as a function of the frequency. When big and little cores operate at the same frequency, big cores perform approximately 3 times faster than little cores. This can be explained by the more complex architecture of ARM A15 compared to ARM A7 (Figure ( 4)), a big core executes in average more instructions than a little core per one cycle.

Figure 5 presents the average power consumption of big (red squares) and little (blue circles) cores as functions of the frequency. A little core consumes 3 -5 times less power than a big core operating at the same frequency. This is due may be to the complex architecture of big cores compared to little cores.

Figure 6 presents the average consumed energy of one instance of matrix multiplication thread (the cumulative power by one instance of the thread) for big (red squares) and little (blue circles) cores as functions of the frequency respectivly. Even if it finishes the executions in a shorter time, a big core still consumes more energy than a little one that operates on the same frequency. Even more, a big core set on the lowest frequency (200Mhz) still consumes more than a little core set at the maximum frequency (1400Mhz).

Figures 8 and7 present power consumption of the whole Odroid Board as function of frequency. The power consumption of all componenets except the big and little cores is constant.

Figure (9) presents the power consumption ratio of big (red squares) and little (blue circle) cores as functions of frequency. When the frequency is low, the board consumes (40-60) times more than little cores, and (6-9) more than little core. However, this ratio decreases substantially when frequency is high (only 8 time for little cores and 1.7 for big cores).

Figure (10) presents the memory power consumption as function of the frequency. The memory consumption is higher When the thread is locked on a little core than when the same thread is locked on a big core. This is probably due to the fact that Little cores have only 1 4 of big core L2 cache, Then, more cache-misses occur causing more energy consumption. This assumption can not be verified because

ARCHITECTURE MODELING

To flexibly model heterogeneous architectures as big.LITTLE, we introduce the concept of core group. A multicore architecture A of m cores is compound of G core group. All cores (Pj) of the same group Gg have the same micro-architecture and frequency characteristics (minimal frequency (f g min ), maximal frequency (f g max ), operating frequency (f g op )). Each group has a set of discrete frequencies modes and its own energy calculation coefficients (discussed in Section 4.1).

Cores are indexed alternatively, core j belongs to the group g = (j%G). For a processor with 4 cores and 2 groups, cores 0, 2 belong to group 0 and cores 1, 3 belong to group 1. We choose this representation to allow to compute the core group just by having its index j. Moreover, by using the concept of core group we can address homogeneous architecture (e.g. SMP) by setting the number of groups G = 1, and heterogeneous architectures like ARM big.LITTLE architecture by setting the number of groups to 2 or more. A = {Gg, g ∈ {1 . . . G}} Gg = ({Pj, (j mod G) = g, j < m}, fmin, fmax, fop), g ∈ {1 . . . G}

Energy model

Increasing a core's frequency involves switching its tran- sistors more rapidly, and transistors that are switched more rapidly dissipate more power. The power dissipated due to switching is called dynamic power. In order to reduce dynamic power, we calibrate the operating frequency for big and little cores. Even transistors that aren't switching will still leak current during idle periods. This leakage current constantly dissipates power. The amount of power dissipated per unit of area due to leakage current is called static power. In order to reduce the static power consumption, we have the ability to set the cores in a deep state power using the ARM Wait For Interrupt and Wait For Event hardware instructions.

The overall power consumption is the sum of the dynamic and static power (total power consumption).Jing Mie et al. [START_REF] Mei | Energy-aware preemptive scheduling algorithm for sporadic tasks on DVS platform[END_REF] defined the dynamic power dissipated by a CMOS circuit as the product of a constant coefficient ξ that depends on the technology, by the square of the voltage, and by the frequency.

E = ξ × V 2 × f (1)
They also defined the frequency as the ratio of the difference between the actual voltage V and V T h raised to the power of a, where V th is the threshold voltage by the product of a constant K and the logic depth L d . a and K are constants that depend on the technology.

f = (V -V T h ) a K × L d (2)
By combining the two equations, the dynamic power E can be expressed as the product of a constant Const by the

E(f ) = Const × f λ (3)
In order to evaluate the total power consumption in the rest of this work, we used a polynomial regression (we used the online regression tool available at http://www.xuru.org/ rt/PR.asp) of the power values measured for big and little cores (described in Figure [START_REF] Mei | Energy-aware preemptive scheduling algorithm for sporadic tasks on DVS platform[END_REF] as a function of variable frequency. Figure (11) presents the results of the polynomial regression of degree 3 for big cores (asterisks) and little cores (squares).

Figure (11) shows that the match between the values and the regression curve for both big and little cores is rather satisfactory. The regression polynomial in Equation ( 4) is for little cores and the one of Equation ( 5) is for big cores.

P ower(f ) l = 0.07477f 3 + 0.045706f + 0.008425

P ower(f ) b = 1.0561341f 3 + 1.376928f 2 + 2.452f + 1.0216 (5) We assume that the execution time C of a thread depends on the micro-architecture of the core group on which the thread is allocated and on its operating frequency. Thus, to model the execution time of a task we obviously need to compute model coefficients for every core group, in our case for big and LITTLE cores.

TASK MODEL

Only a part of the execution time depends on the core operating frequency. For example, when a thread suspends itself waiting for operating on external devices, the suspension time depends on the response time of the external device which may not be linked to the core operating frequency. As an relevant example, consider the time to access main memory in case of a cache miss: the response time depends on the bus contention (due to other cores/devices) and on the response time of the memory circuitry.

The component of the execution time that directly depends on the frequency is denoted by ct(f ), whereas the part that does not depends on the frequency is denoted by it. For simplicity, in this paper we consider this second part as a constant which only depends on the thread average characteristics (e.g. number of memory requests and average number of cache misses). Thus, the execution time of a thread on a core of frequency f is modelled by a semi-linear 0. 6)

C g (f ) = ct • fmax × fmax f + it (6)
To estimate the two components of the execution time, we measure the thread execution time under different conditions (different inputs, differents matrix sizes, frequencies, groups). Figure 12 presents a non-linear regression of the execution time values (Figure 3) as function of frequency for matrix multiplication of size 300x300. Equation ( 7) presents the semi-linear function of the regression for little cores (resp. Equation 8 for big cores).

C l = 4.27 f + 0.32 (7) C b = 1.45 f + 0.07 (8)
We observe that, in this specific case the it time for little cores is 4 time bigger than the it time of little cores.

We assume that the it component of the execution time depends on the number of cache misses experienced by the task, which in turns depends not only on the size of the task data set, but also on the interference by the other threads in the system. However, we cannot currently measure the number of cache misses experienced by a task on the bigLIT-TLE processor: due to a software error in the Linux device driver for the ODROID platform, we cannot currently collect this precious data. We are currently actively working at solving the problem, and we hope to report experiments directly linking the number of cache misses to the execution time of the thread.

CONCLUSIONS AND FUTURE WORK

In this paper we presented some preliminary results of a study for modelling the execution time and the power consumption of soft real-time tasks running on a heterogeneous processor architecture, namely the bigLITTLE architecture.

We have shown how to compute the parameters of the power consumption model of the architecture, and how to model the execution time of a thread as a function of the core frequency and the memory access. We found that the number of cache misses has a big impact on the execution time and hence on the power consumption of the memory.

In the future we plan to carefully investigate the relationship between data set size, number of cache misses, execution time and power consumption. In fact, accessing main memory slows down the task and increases the energy consumed by the main memory itself; moreover DVFS techniques do not influence the component of the execution time due to main memory access.

It is clear that this kind of study is preliminary to the proposition of power management heuristics. In addition, we are currently studying the use of machine learning techniques to identify the parameters of a task at run-time.

Figure 1 :

 1 Figure 1: Pseudo-code of a soft-real time periodic task

Figure 2 :

 2 Figure 2: ODROID XU3 board architecture

Figure 3 :Figure 5 :

 35 Figure 3: Execution Time with different frequencies on big and little cores

Figure 4 :Figure 6 :

 46 Figure 4: Pipeline architecture of A7 (left) and A15 (right) cores.

Figure 7 :

 7 Figure 7: Little Core power consumption: Internal vs External measures

Figure 9 :

 9 Figure 8: Big Core power consumption: Internal vs External measures

Figure 10 :

 10 Figure 10: Memory power consumption: Big VS Little cores

Figure 11 :Figure 12 :

 1112 Figure 11: Regression and real values for big and little power consumption: cores consumption

  board (Figure[START_REF] Kimball Dunn | Development of smart phone usability benchmarking tasks[END_REF]) is compound of a samsung Exynos 5422, a Mali GPU, a RAM memory and I/O peripherals. The Samsung Exynos 5422 is an ARM big.LITTLE multicore architecture. It is composed of 8 cores: 4 big cores (ARM cortex A15 ) and 4 LITTLE cores (ARM Cortex A7 ). The ODROID XU3 board embeds 4 power sensors: a sensor for big cores, a sensor for little cores, a sensor for the GPU and the last one for the memory. An external energy sensor (Odroid SmartPower) is plugged to the power supply of the board in order to measure the overall power consumption.Each core of the Exynos 5422 has a private L1 cache of 2 × 32Kb. Little cores share 512Kb of L2 cache. Big cores

					Central Memory
		L2 Cache		L2 Cache
			Big's cluster	Little's cluster
		L1 Cache	L1 Cache	L1 Cache	L1 Cache
		ARM	ARM	ARM	ARM
	Mali	Cortex A15	Cortex A15	Cortex A7	Cortex A7
	GPU	L1 Cache	L1 Cache	L1 Cache	L1 Cache
		ARM	ARM	ARM	ARM
		Cortex A15	Cortex A15	Cortex A7	Cortex A7


Alain Billionnet

Marie-Christine Costa

Pierre-Louis Poirion
email: pierre-louis.poirion@ensta-paristech.fr

2-Stage Robust MILP with continuous recourse variables

come

Introduction

This paper deals with robust mixed-integer linear programming (MILP) to study problems with uncertain data. This is a possible alternative to two-stage stochastic linear programming introduced by Dantzig in [START_REF] Dantzig | Linear programming under uncertainty[END_REF]. In this framework the uncertain data of the problem are modeled by random variables, and the decision-maker looks for an optimal solution with respect to the expected objective value. He makes decisions in two stages: first before discovering the actual value taken by the random variables, second once uncertainty has been revealed. However, this approach requires to know the underlying probability distribution of the data, which is, in many cases, not available; furthermore the size of the resulting optimization model increases in such a way that the stochastic optimization problem is often not tractable. Robust optimization is a recent approach that does not rely on a prerequisite precise probability model but on mild assumptions on the uncertainties involved in the problem, as bounds or reference values of the uncertain data. It looks for a solution that remains satisfactory for all realizations of the data (i.e. for worst scenarios). It was first explored by Soyster [START_REF] Soyster | Convex programming with set-inclusive constraints and applications to inexact linear programming[END_REF] who proposed a linear optimization model for data given in a convex set. However this is an over conservative approach that leads to optimal solutions too far from the one of the nominal problem. Robust adjustable optimization models have been proposed and studied to address this conservatism. More precisely, a lot of recent published works cover robust linear programming with row-wise uncertainty for continuous variables [START_REF] Ben-Tal | Adjustable robust solutions of uncertain linear programs[END_REF][START_REF] Ben-Tal | Robust solutions of uncertain linear programs[END_REF][START_REF] Bertsimas | Robust linear optimization under general norms[END_REF][START_REF] Bertsimas | The price of robustness[END_REF] or discrete variables [START_REF] Atamturk | Strong formulations of robust mixed 0-1 programming[END_REF][START_REF] Kouvelis | Robust Discrete Optimization and its Applications[END_REF] and, even more recently, column-wise right-hand side uncertainty [START_REF] Bertsimas | Theory and applications of robust optimization[END_REF][START_REF] Minoux | On 2-stage robust LP with RHS uncertainty: complexity results and applications[END_REF][START_REF] Thiele | Robust linear optimization with recourse[END_REF].

In [START_REF] Gabrel | Robust location transportation problems under uncertain demands[END_REF], Gabrel et al. propose a solution based on the approaches given in these last papers to solve a location transportation problem. We first show that their solution can be applied to any linear program with mixed-integer first stage variables and continuous recourse variables. We will also see that problems with left-hand side uncertainty can be solved in the same way. Then, we show that the method can also be used for an affine definition of the uncertainty set which is more general than the one used in [START_REF] Gabrel | Robust location transportation problems under uncertain demands[END_REF][START_REF] Minoux | On 2-stage robust LP with RHS uncertainty: complexity results and applications[END_REF][START_REF] Thiele | Robust linear optimization with recourse[END_REF].

To the extent of our knowledge, in all works published until now, the authors always assumed that the problem satisfies a "full recourse property" (see Section 2) which cannot be always satisfied for real problems: we show that, when this property is not verified, we can modify the objective function in order to use the previous approach to solve the problem. Some tests on a simple production problem prove the feasibility of our approach.

We focus here on a linear robust problem with right-hand and left hand-side uncertainty, mixed-integer first-stage variables and continuous second-stage variables. In Section 2 we present the general problem. For the sake of clarity, we first study the robust problem with right-hand side uncertainty and full recourse property with a specific definition of the uncertainty set. In Section 3 we show how to modelize and solve the recourse problem. In Section 4 we present the solution for the robust problem. In Section 5 we show that our results can be applied in case of left-hand side uncertainty and we extend our results to other definitions of the uncertainty set. Finally, in Section 6 we study the cases where the full recourse property is not verified. To facilitate reading the article, we present in Annex 1 and 2 two proofs of results given in Section 6.1 .

A mixed-integer linear robust problem

We consider applications requiring decision-making under uncertainty which can be modeled as a two-stage mixed-integer linear program with recourse. The set of variables is partitioned into two distinct sets: the x variables, called decision variables, concern the decisions to be taken in the first stage, before knowing the realization of the uncertain events; the second stage variables y, called recourse variables, will be fixed only after the uncertainty has been revealed.

We focus here on robust mixed-integer linear problems when the constraints coefficients are uncertain, as well on the right-hand side as on the left-hand side. In addition, we restrict our study to the case where the recourse variables y are continuous variables while the decision variables x are mixed-integer variables.

The deterministic problem can be formulated as the following MILP (in this paper we will omit the transpose sign tr when there is no possible confusion):

P min x,y αx + βy Ax + By ≥ d Cx ≥ b x i ∈ N, i = 1, ..., p 1 , x i ∈ R + , i = (p 1 + 1), ..., p, y ∈ R q + . (1) (2) (3)
where

A ∈ Q T ×p , B ∈ Q T ×q , d ∈ Q T , C ∈ Q n×p , b ∈ Q n , α ∈ Q p + , β ∈ Q q + ,
and Q is the set of rational numbers.

We assume that there exists (x, y) such that (1)-( 3) are satisfied and we say that a solution x is feasible if x satisfies constraints (2) and (3). The uncertain coefficients are those of d (right-hand side), A and B (left-hand side).

We assume that the program P satisfies the property P, called "full recourse property": for any feasible values of the decision variables (here x) and for any possible value of A,B and d, there exist values of the recourse variables (here y) such that (1) is satisfied, that is such that there exists a feasible solution of P . Let us notice that the property P is always satisfied if there is a column of B whose all terms are positive. The hypothesis that P satisfies the property P cannot be always satisfied for real problems: we show in Section 6 that we can extend our results when P is not satisfied.

We suppose that d belongs to a given set D which defines the set of possible scenarios and for the sake of clarity we assume at first that the uncertainty concerns only the right-hand side d of (1). We show in Section 5.1 that our results can be extended when the matrices A and B are also uncertain.

Our robustness objective is to find a feasible solution x, y of P that minimizes the total cost involved by the worst possible scenario of d in connection with x. We can state the robust problem as the following mathematical program: + . Given a mathematical program π, we denote by v(π) the value of an optimal solution. The robust program can then be rewritten as:

P R min x αx + max d∈D min y βy By ≥ d -Ax y ∈ R q + Cx ≥ b x i ∈ N, i = 1, ..., p 1 , x i ∈ R + , i = (p 1 + 1), ..., p.
P R min x αx + v(R(x)) Cx ≥ b x i ∈ N, i = 1, ..., p 1 , x i ∈ R + , i = (p 1 + 1), ..., p.
Let us define more precisely the uncertainty. Following the idea proposed by Bertsimas and Sim [START_REF] Bertsimas | The price of robustness[END_REF] and Minoux [START_REF] Minoux | On 2-stage robust LP with RHS uncertainty: complexity results and applications[END_REF], we suppose that each coefficient d t , t = 1, ..., T belongs to an interval [ dt -∆ t , dt + ∆ t ] where dt is a given value and where ∆ t ≥ 0 is a given bound of the uncertainty of d. The uncertainty set D is therefore given by:

D = {d : d t ∈ [ dt -∆ t , dt + ∆ t ], ∀t = 1, ..., T }.
For a fixed x, the worst scenario is obtained for d t = dt + ∆ t , t = 1, ..., T. Indeed, By ≥ d + ∆ -Ax implies By ≥ d -Ax for all d ∈ D. Thus this uncertainty definition brings the robust problem back to a deterministic one. It provides a high "protection" against uncertainty, but it is very conservative in practice and leads to very expensive solutions. To avoid overprotecting the system, we impose, as in [START_REF] Thiele | Robust linear optimization with recourse[END_REF], the constraint

T t=1 ∆t>0 | d t -dt ∆ t | ≤ δ,
where δ is a positive integer which bounds the total scaled deviation of d from its nominal value d. Notice that there always exists a worst scenario with d t ≥ dt , ∀t, hence we can redefine the uncertainty set D as

D = {d : d t = dt + δ t ∆ t , 0 ≤ δ t ≤ 1, ∀t = 1, ..., T, T t=1 δ t ≤ δ}.

The recourse problem

To solve the recourse "max min" problem for given values of the decision variables, the minimization linear sub-program is transformed in a maximization program by considering its dual. But that leads to a quadratic objective function. We show that whatever the coefficients in the recourse problem, the quadratic terms can be written as products of a 0-1 variable and a continuous but bounded variable, which allows a linearization of these products. Let x be a feasible solution and let d ∈ D, we define the following linear program

R(x, d) min y βy By ≥ d -Ax y ∈ R q
+ . We notice that R(x, d) has a finite solution for all feasible x and for all possible scenario d since P satisfies P, and since βy ≥ 0 for any feasible solution y of R(x, d). Thus by the strong duality theorem, we have

v(R(x, d)) = v(DR(x, d)), where DR(x, d) is the dual program of R(x, d): DR(x, d) max λ (d -Ax)λ λB ≤ β λ ∈ R T + . (4) (5)
Then, for any feasible x, v(R(x)) = max d∈D v(DR(x, d)), and we can reformulate DR(x) as

max δ: T t=1 δt≤ δ 0≤δt≤1, t=1,...,T max λ: λB≤β λ∈R T + . T t=1 [( dt + δ t ∆ t -(Ax) t )λ t ]
where for a vector (u), we denote by (u) t the t-th coordinate of (u).

DR(x) can be written as follows:

DR(x) max λ,δ T t=1 [( dt -(Ax) t )λ t + ∆ t δ t λ t ] λB ≤ β T t=1 δ t ≤ δ 0 ≤ δ t ≤ 1 t = 1, ..., T λ ∈ R T + . (4) (6) (7)
(5)

However, this bilinear program with linear constraints is not concave. Therefore computing the optimal solution of DR(x) written as above is not an easy task.

We now prove that we can solve DR(x) by solving an equivalent mixed-integer linear program. To prove this claim, we need the following proposition:

Proposition 1. There is an optimal solution λ * , δ * of DR(x) such that δ * t ∈ {0, 1}, 1 ≤ t ≤ T.
Proof. For any fixed λ, there is an optimal solution, (λ, δ * ), of DR(x), where δ * is an extreme point of the polyhedron defined by ( 6) and [START_REF] Bertsimas | The price of robustness[END_REF], that is to say, a point such that δ * t ∈ {0, 1}, 1 ≤ t ≤ T, since δ is an integer. More precisely δ * t = 1 for indices corresponding to the δ largest ∆ t λ t . Therefore we can assume that there is an optimal solution of DR(x), such that λ t δ t belongs to {0, λ t }. To linearize λ t δ t , we now prove that we can restrict ourselves to the case where λ t is bounded by a constant Λ, for all t. Proposition 2. There exists Λ > 0 which can be calculated from B and β such that the conditions λ t ≤ Λ, t = 1, ..., T, can be added to DR(x) without loss of generality.

Proof. Let x be feasible. Let us rewrite DR(x) with the slack variables λ t ≥ 0, t = 1, ..., T. The constraints (4) become: B tr λ + λ = β. Let (λ * , λ * , δ * ) be an optimal solution of DR(x), we can assume w.l.o.g. that (λ * , λ * ) is an optimal basic solution of DR(x) when δ is set to δ * . Therefore, there exists a basic matrix E = (e ij ) of (B tr I T ) and basic vectors

λ * E , λ * E such that: (λ * E λ * E ) tr = E -1 β.
Let ê be an upper bound on the absolute value of the coefficients of E -1 for all basic matrices E of (B tr I T ), and let β = max i=1,...,q β i , we have λ * t ≤ ê βq, t = 1, ..., T. Therefore there exists an optimal solution (λ * , δ * ) of DR(x) such that λ * t is bounded by Λ = ê βq for any t = 1, ..., T.

We can now linearize DR(x) by substituting the new variables ν t to the products λ t δ t and by adding the constraints:

ν t ≤ λ t , ν t ≤ Λδ t , ν t ≥ λ t -Λ(1 -δ t ), ν t ≥ 0.
DR(x) is equivalent to the following mixed-integer linear program:

LDR(x) max λ,δ,ν T t=1 [( dt -(Ax) t )λ t + ∆ t ν t ] λB ≤ β T t=1 δ t ≤ δ ν t ≤ λ t , t = 1, ..., T ν t ≤ Λδ t , t = 1, ..., T λ, ν ∈ R T + δ t ∈ {0, 1}, t = 1, ..., T.
Notice that the linearization constraints, ν t ≥ λ t -Λ(1 -δ t ), t = 1, ..., T , can be omitted since the coefficients of ν t in the objective function to maximize are positive.

Solving the robust problem

In order to solve the robust problem P R, we will first reformulate it as a linear program and then use a constraint generation algorithm. In the previous section, we proved that the recourse problem is equivalent to the linear program LDR(x). Thus the robust problem can be reformulated as:

P R min x αx + v(LDR(x)) Cx ≥ b x i ∈ N, i = 1, ..., p 1 , x i ∈ R + , i = (p 1 + 1), ..., p.
Let P Q be the polyhedron defined by the constraints of LDR(x) where we replace δ t ∈ {0, 1} by 0 ≤ δ t ≤ 1, and let (P Q ) I = conv(P Q ∩ {δ ∈ N m }), be the convex hull of the feasible solution of LDR(x). Notice that this convex hull does not depend on x. (P Q ) I is a polyhedron, thus we have

LDR(x) max λ,δ,ν T t=1 [( dt -(Ax) t )λ t + ∆ t ν t ]   λ δ ν   ∈ (P Q ) I ,
Let S = {(λ s , δ s , ν s ) 1≤s≤S }, be the set of extreme points of (P Q ) I . For any feasible x, there is s ∈ {1, ..., S} such that (λ s , δ s , ν s ) is an optimal solution of Algorithm 1 Constraint generation algorithm 1: (λ 0 , δ 0 , ν 0 ) = (0, 0, 0). Set L ← -∞, U ← +∞, k ← 1. 2: Solve the master problem :

P R k min x,z αx + z z ≥ T t=1 ( dt -(Ax)t)λ s t + ∆tν s t , 0 ≤ s ≤ k -1 Cx ≥ b x i ∈ N, i = 1, ..., p 1 , x i ∈ R + , i = (p 1 + 1), ..., p z ∈ R Let (x k , z k ) be the obtained solution. L ← αx k + z k . 3: Solve LDR(x k ). Let (λ k , δ k , ν k ) be the optimal solution. U ← min{U, αx k + v(DR(x k ))}. if U = L, then return (x k , z k ) else go to 4. 4: Add the constraint z ≥ T t=1 ( dt -(Ax)t)λ k t + ∆tν k t ,
to the master problem P R k , k ← k + 1 and go to 2.

LDR(x).

Thus the robust problem can be reformulated as the linear program:

P R min x,z αx + z z ≥ T t=1 [( dt -(Ax) t )λ s t + ∆ t ν s t ], 1 ≤ s ≤ S Cx ≥ b x i ∈ N, i = 1, ..., p 1 , x i ∈ R + , i = (p 1 + 1), ..., p, z ∈ R (8)
However, due to the potentially tremendous number of constraints, we solve P R by a constraint generation algorithm as in [START_REF] Thiele | Robust linear optimization with recourse[END_REF] or [START_REF] Gabrel | Robust location transportation problems under uncertain demands[END_REF]. Initially, we consider a subset S 0 of S; at a step k, we consider a subset S k of S and we solve a relaxed program P R k of P R, called master problem, which consists in solving P R with the subset of constraints [START_REF] Dantzig | Linear programming under uncertainty[END_REF] corresponding to S k . The obtained solution in denoted by (x k , z k ). Then we solve DR(x k ), called slave problem, to check if (x k , z k ) is optimal. If not, then a new constraint is added, i.e. an extreme point is added to S k (See Algorithm 1). On the basis that the number of extreme points of (P Q ) I is finite, one can prove that this algorithm converges in a finite number of steps.

5 Some generalizations

Left-hand side uncertainty

In the previous sections, we assumed that the uncertainty concerned only the righthand side d of constraints [START_REF] Atamturk | Strong formulations of robust mixed 0-1 programming[END_REF]. We now prove that our approach can be generalized to the case where the constraint coefficients (A = (A ti ) 1≤t≤T, 1≤i≤p ), are also likely to be uncertain. As before, we assume that each coefficient A ti belongs to an interval [ Āti -Γ ti , Āti +Γ ti ], where Āti is a given value and where Γ ti is a given bound of the uncertainty of A ti . Furthermore, in order to avoid overprotecting the system, we assume that the total scaled deviation of the uncertainty of the i-th column of A, A i = (A ti , t = 1, ..., T ), is bounded. Similarly to D, the uncertainty set A i of A i is defined as :

A i = {A i : A ti = Āti -γ ti Γ ti , 0 ≤ γ ti ≤ 1, t = 1, ..., T, T t=1 γ ti ≤ γi },
where γi is a given integer.

The robust problem can thus be formulated as:

P R min x αx + max A i ∈A i , ∀i=1,...,p d∈D min y βy By ≥ d -(A 1 , ..., A p )x y ∈ R q + Cx ≥ b x i ∈ N, i = 1, ..., p 1 , x i ∈ R + , i = (p 1 + 1), ..., p.
And the recourse problem becomes:

DR (x) max λ,δ,γ T t=1 [( dt -p i=1 Āti x i )λ t + ∆ t δ t λ t + p i=1 Γ ti x i γ ti λ t ] λB ≤ β T t=1 δ t ≤ δ 0 ≤ δ t ≤ 1, t = 1, .., T λ ∈ R T + T t=1 γ ti ≤ γi , i = 1, ..., p 0 ≤ γ ti ≤ 1, i = 1, .., p, t = 1, ..., T.
We can then linearize the quadratic terms (δ t λ t and γ ti λ t ), to obtain a mixedinteger linear recourse problem and then solve the robust problem as we did in the previous sections.

Let us now consider that matrix B has uncertain coefficients. As before we assume that, for t = 1, ..., T and i = 1, ..., q, B ti ∈ [ Bti -Φ ti , Bti + Φ ti ] and that column i of B is in the set {B i :

B ti = bti -φ ti Φ ti , 0 ≤ φ ti ≤ 1, t = 1, ..., T, T t=1 φ ti ≤ φi }.
The dual of the recourse problem is similar to the program DR(x) defined in Section 3: we add to DR(x) variables φ and constraints T t=1 φ ti ≤ φi , i = 1, ..., q and 0 ≤ φ ti ≤ 1, i = 1, ..., q, t = 1, ..., T , and we replace constraints (4) by:

T t=1 ( Bti -φ ti Φ ti )λ t ≤ β i , i = 1, ..., q (4 )
.

Now there are quadratic terms (φ ti λ t ) in the constraints. To linearized these terms, we must verify that the propositions 1 and 2 can be extended, i.e. there is an optimal solution such that φ ∈ {0, 1} T q and λ is bounded. The proofs given in Section 3 can easily be extended. Let λ be fixed. If constraints (4 ) and T t=1 φ ti ≤ φi , i = 1, ..., q are verified, a feasible integer solution can be obtained by setting φ ti to 1 if i corresponds to one of the φi larger values of Φ ti λ t and to 0 otherwise.

For Proposition 2, by using the same arguments as in Section 3, for any fixed φ we can bound λ t by Λ(φ). Since φ ti ∈ {0, 1} and T t=1 φ ti ≤ φi , there are 2 qT different values of φ and we can bound λ t by max φ Λ(φ).

Generalization to other uncertainty sets

In the previous sections, we assumed that uncertain coefficients could be written as d t = dt +δ t ∆ t ∀t, where δ t expresses the uncertainty on d t and satisfies T t=1 δ t ≤ δ. Now we generalize our results when the vector d can be written as d = d + Dδ, where the vector d and the matrix D are given and where δ belongs to a bounded polyhedron D whose extreme points (d 1 , ..., d S ) are known. Notice that this definition of the uncertainty covers the one given by Babonneau et al. in [START_REF] Babonneau | Robust capacity expansion solutions for telecommunication networks with uncertain demands[END_REF] Let us rewrite the recourse problem:

DR (x) max λ,δ ( d + Dδ -Ax)λ λB ≤ β δ ∈ D λ ∈ R T + .
Let v 1 , ..., v S ∈ [0, 1] be variables such that δ = S s=1 d s v s and S s=1 v s = 1. We can rewrite the recourse problem as

DR (x) max λ,v ( d -Ax)λ + S s=1 (v s (Dd s )λ) λB ≤ β S s=1 v s = 1 0 ≤ v s ≤ 1, s = 1, ..., S λ ∈ R T + .
Using the same argument as in Proposition 1, we can prove that there exists an optimal solution (λ * , v * ) of DR (x) such that either v s * = 1 or v s * = 0, s = 1, ..., S. Therefore we can linearize the quadratic terms v s λ t , for all s and for all t, as we did in Section 3, to obtain a mixed-integer linear recourse problem, and finally we can solve the robust problem by using Algorithm 1.

6 The problem without the full recourse property

Solving the problem

In the previous sections, we assumed that the deterministic problem satisfied the property P. We now prove that we can extend our results to the case where we only assume that the robust problem P R has a solution, i.e. there is x such that for all d in D there is a feasible solution y to R(x, d) and so, there is M such that v(P R) ≤ M . In addition, the method detects if the problem has no solutions. For the sake of clarity, we give here the proof for the case where all the decision variables are integer. The complete proof for mixed integer decision variables is given in Annex 1. In Annex 2, we propose a general theoretical value of M ; nevertheless, a better specific value of M can be calculated for each problem, as it is the case for the application presented in Section 6.2.

First we show how to obtain a new MILP, denoted P ε such that the robust associated problem has the same optimal solution as the initial robust problem if there is one, and P ε satisfies P. To obtain P ε , we add new recourse variables w t , t = 1, ..., T. As in Sections 2, 3 and 4, for the sake of clarity and w.l.o.g., we consider only right-hand side uncertainty.

Let ε be a given strictly positive value, we define the following MILP:

P ε min x,y,w αx + βy + M ε T t=1 w t Ax + By + w ≥ d Cx ≥ b x i ∈ N, i = 1, ..., p y ∈ R q + , w ∈ R T + (1ε) (2) (3) (4ε)
Since the variables w t , t = 1, ..., T, are not bounded, P ε satisfies property P.

We denote by P R ε , the robust problem associated to P ε ,

P R ε min x αx + max d∈D min y,w βy + M ε T t=1 w t By + w ≥ d -Ax y ∈ R q + , w ∈ R T + Cx ≥ b x i ∈ N, i = 1, ..., p,
and by R ε (x), R ε (x, d) and DR ε (x) the associated subproblems as those defined in Section 3. Notice that since all the inputs, A, B, C, b, d, α, β, ∆, of P R have rational coefficients, we can reduce P R ε and all the corresponding subproblems to programs where all the inputs are integer. Therefore we assume from now that all the inputs are integer.

Proposition 3. v(P R ε ) satisfies 0 ≤ v(P R ε ) ≤ v(P R) ≤ M.
Proof. Let (x, ŷ) be an optimal solution of P R, By hypothesis,

v(P R) ≤ M thus v(R(x)) ≤ M. Let d be a scenario in D. Since v(R(x)) = max d∈D v(R(x, d)),
we have v(R(x, d)) ≤ M. Let ȳ be an optimal solution of R(x, d), we notice that (y, w) = (ȳ, 0) is a feasible solution of (R ε (x, d)) with the same cost. Thus

v(R ε (x, d)) ≤ v(R(x, d)), for any d ∈ D, which implies v(R ε (x)) ≤ v(R(x)), and
0 ≤ v(P R ε ) ≤ v(P R) ≤ M.
Let (x * , d * , y * , w * ) be an optimal solution of P R ε : d * is the worst scenario for x * . Notice that Proposition 1 is valid for DR ε (x * ). Therefore d * t = dt or d * t = dt + ∆ t , and d * t is an integer. From Proposition 3, we have

αx * + βy * + M ε T t=1 w * t ≤ M.
Since αx * + βy * ≥ 0, we have T t=1 w * t ≤ ε, and thus w * t ≤ ε, ∀t = 1, .., T. We now prove that if (y * , w * ) is a basic optimal solution of (R ε (x * , d * )), then for ε small enough, we have w * t = 0, ∀t; and then (x * , y * ) is feasible for P R, and therefore from Proposition 3, v(P R) = v(P R ε ).

Let us rewrite P R ε with the positive slack variables σ = (σ t , t = 1, ..., T ): the constraint (1ε) becomes Ax + By + w -σ = d. Let (x * , d * , y * , w * , σ * ) be an optimal solution where (y * , w * , σ * ) is a basic optimal solution of the program

R ε (x * , d * ) min y,w,σ βy + M ε T t=1 w t By + w -σ = d * -Ax * y ∈ R q + , σ, w ∈ R T + .
which is equivalent to

R ε (x * , d * ) min y,w,σ βy + M ε T t=1 w t B I T -I T   y w σ   = d * -Ax * y ∈ R q + , σ, w ∈ R T + . Let L = B I T -I T = (l ij ) ∈ Z T ×(q+2T
) and l M = max i,j |l ij |. We notice that L has rank T .

Proposition 4. Let ε < 1 (l M ) T T T /2 . If (x * , d * , y * , w *
) is an optimal solution of P R ε then w * = 0 and (x * , d * , y * ) is an optimal solution of P R.

Proof. First, assume that (y * , w * ) is a basic optimal solution of R ε (x * , d * ).

There exists a basic matrix E ∈ Z T ×T of L and basic vectors

y * E , w * E , σ * E such that: E (y * E w * E σ * E ) tr = d * -Ax * .
The matrix E = (e kj ) is invertible, and

E -1 = 1 det(E) adj(E), where adj(E) is the adjugate matrix of E. Therefore (y * E w * E σ * E ) tr = 1 det(E) adj(E)(d * -Ax * ), and 0 ≤ w * t = 1 det(E) (adj(E)(d * - Ax * )) t ≤ ε, where w *
t is a basic variable and where t is the associated index in

(y * E , w * E , σ * E ). Thus |adj(E)(d * -Ax * )) t | ≤ ε| det(E)|. (9)
Hadamard's inequality states that

| det(E)| ≤ rank(E) j=1 rank(E) k=1 e 2 kj .
Thus we have

| det(E)| ≤ (l M ) T T T /2 . If ε < 1 (l M ) T T T /2 , then according to (9), |(adj(E)(d * -Ax * )) t | < 1. Since |(adj(E)(d * -Ax * )) t | ∈ N, we have |(adj(E)(d * -Ax * )) t | = 0,
therefore w t = 0 for any basic variable w t and thus w t = 0 for all t = 1, ..., T.

Thus, with ε <

1 (l M ) T T T /2 , if (y * , w * ) is a basic optimal solution of R ε (x * , d * ) then w * = 0.
Now, assume that (y * , w * ) is a non basic optimal solution of R ε (x * , d * ). We can write (y * , w * ) as the sum of a convex combination of basic optimal solutions (satisfying w s * = 0 for all s) and a positive combination of extreme rays of the constraints polyhedron of R ε (x * , d * ). Since all the optimal solutions of R ε (x * , d * ) satisfy w s * ≤ ε for all s = 1, ..., S, no extreme ray in the decomposition can have w s t * = 0 for some (s, t) with s ∈ {1, ..., S} and t ∈ {1, ..., T }.

Finally, for any optimal solution (x * , y * , w * , s * ) of P R ε we have w * = 0, and then (x * , d * , y * ) is a feasible solution of P R. Then from Proposition 3, v(P R) = v(P R ε ) and so (x * , d * , y * ) is an optimal solution of P R.

To summarize, we solve P R by solving P R ε : if the optimal solution (x * , y * , w * ) of P R ε verifies w * = 0 then (x * , y * ) is an optimal solution of P R. Otherwise w * t > ε for some t, v(P R ε ) > M and P R has no solution.

Application to a production problem

In this section, we test our approach on a production problem. A company decides to install factories to manufacture several products. There are p possible factory locations and a factory at site i produces a quantity a it of each product t, t = 1, ..., T , for a total production cost equal to α i , i = 1, ..., p. The demand of a product t is uncertain and is denoted as before by d t = dt + δ t ∆ t , t = 1, ..., T . If the production is not sufficient, the company must buy the missing product at a high cost to meet the demand but there is a given bound to product t available on the market: β t is the unit purchasing cost for product t, and K t the maximal quantity of t which can be bought, t = 1, ..., T . The decision variables are denoted by x i , x i = 1 if a factory is installed at site i, x i = 0 otherwise, i = 1, ..., p. The recourse variables are denoted by y t , t = 1, ..., T , y t is the (possible) lacking quantity of product t to purchase. We don't take into account the selling price of products and profit in our model. The problem to solve is the following:

P R min x p i=1 α i x i + max d∈D min y T t=1 β t y t y t ≥ d t -p i=1 a it x i ∀t = 1, ..., T y t ≤ K t ∀t = 1, ..., T y ∈ R T + x i ∈ {0, 1} i = 1, ..., p.
We assume that if there are factories on the p sites then the demand is always met, i.e. p i=1 a it + K t ≥ dt + ∆ t for all t = 1, ..., T . In that way, the problem always has a solution. In addition, the data verify a it + K t < dt , for at least one (i, t) ∈ {1, ..., p} × {1, ..., T }. In this way the full recourse property is not verified since for some values of the decision variables such that x î = 1 and x i = 0 for all i = î there is no solution to the recourse problem whatever the values of δ t .

To ensure that w t = 0 for all t at the end of the algorithm, we set M = p i=1 α i + T j=1 β j K j and ε = 1 . So M is greater than the worst possible cost and ε is equal to the absolute value of any determinant of a submatrix of the constraint matrix in the recourse problem which is totally unimodular in this case (we don't need to use Hadamard's inequality in inequation ( 9)).

The data are generated in the following way: For a given p and a given T , we generate randomly the coefficients α, β, d, K, between 0 and 100. Then for each t = 1, ..., T we generate randomly ∆ t between 0 and dt . Then, we compute randomly a it to ensure that a it + K t < dt , for at least one (i, t) ∈ {1, ..., p} × {1, ..., T }. Finally we generate randomly δt between 0 and T.

As might be expected, applying Algorithm 1 without the artificial variables (i.e. with M = 0) gives no solution: the program returns "Dual infeasible due to empty colomn". Eight values of (p, T ) are tested. For each value, we generate five instances and we give average results in Table 1. The column 'CPU time' gives the time in seconds to obtain an optimal solution on a Bi-pro Intel Nehalem XEON 5570 at 2.93 GHz with 24 Go of RAM. The next column gives the total number of iterations in Algorithm 1. The last column precises the number of iterations where the recourse property is not satisfied, that is the number of iterations where the solution returned by the recourse problem is not a feasible solution (i.e. w > 0). The most difficult instances are those where the number of sites and products are similar while it is easy to solve instances with up to 1500 sites (resp. products) and 100 products (resp. sites). The CPU time seems to be unrelated to the ratio between the number of iterations and the number of unfeasible solutions in the recourse problem.

Finally, we notice that our method to solve robust problems which does not satisfy the full recourse property is practicable and allow large scale instances to be solved.

Annex 1: problem including continuous decision variables and not verifying the full recourse property

Now we consider the problem with mixed integer decision variables: x i ∈ N, i = 1, ..., p 1 , x i ∈ R + , i = (p 1 + 1), ..., p. We assume that P R has a solution and so there is M such that v(P R) ≤ M . Following the ideas of Section 6.1 for integer decision variables, we are going to show that there is ε such that we can solve P R by solving P R ε .

From Proposition 1 the worst scenario for any feasible x is an extreme point of D. Let {d 1 , ..., d S }, be the set of extreme points of D; all these points are integers and S = T δ + 1. We can rewrite the problem P R ε as the following MILP where y s and w s correspond to the recourse variables associated to scenario s:

P R ε min x,z y 1 ,...,y S w 1 ,...,w S α 1 x 1 + α 2 x 2 + z z ≥ βy s + M ε T t=1 w s t , s = 1, ..., S A 1 x 1 + A 2 x 2 + By s + w s ≥ d s , s = 1, ..., S C 1 x 1 + C 2 x 2 ≥ b x 1 ∈ N p 1 , x 2 ∈ R p-p 1 + y s ∈ R q + , w s ∈ R T + , s = 1, ..., S (10)
where

x = x 1 x 2 , A = (A 1 , A 2 ) and C = (C 1 , C 2 )
. Let (x * , z * , y 1 * , ..., y S * , w 1 * , ..., w S * ) be an optimal solution of this program; there is ŝ ∈ {1, ..., S} such that constraint (10) is saturated:

z * = βy ŝ * + M ε T t=1 w ŝ *
t and ŝ is the worst scenario associated to x * . The scheme of the proof is the following: first we prove that for ε small enough we have w s = 0 for all s = 1, ..., S in any optimal solution of P R ε , then we deduce that w = 0 in any optimal solution of P R ε and finally, from an optimal solution of P R ε we obtain an optimal solution of P R. (ST +n) (ST +n)/2 then any optimal solution of P R ε verifies w s = 0, for all s = 1, ..., S.

Proof. Let (x 1 * , x 2 * , z * , y 1 * , ..., y S * , w 1 * , ..., w S * ) be an optimal solution of P R ε . From constraints [START_REF] Kouvelis | Robust Discrete Optimization and its Applications[END_REF] and Proposition 3, which is true for mixed integer variables, we have for all s = 1, ..., S: x 2 ∈ R p-p 1 + y s ∈ R q + , w s ∈ R T + , s = 1, ..., S which verifies v(D ε (x 1 * )) = v(P R ε (x 1 * )) = v(P R ε ) . Let (λ , x 2 , z , y 1 , ..., y S , w 1 , ..., w S ) be an optimal solution of D ε (x 1 * ), then (x 2 , z , y 1 , ..., y S , w 1 , ..., w S ) is an optimal solution of D ε (x 1 * , λ ) which is the minimization part of D ε (x 1 * ) obtained by fixing λ to λ .

Considering the slack variable e s , s = 1, ..., S and f , we can rewrite D ε (x 1 * , λ ) as x 2 ∈ R p-p 1 + y s ∈ R q + , w s , e s ∈ R T + , s = 1, ..., S.

Now we are going to prove that if ε < 1 l (ST +n) M (ST +n) (ST +n)/2 , every optimal solution u = (x 2 , z , y 1 , ..., y S , w 1 , ..., w S , e 1 , ...., e S , f ), of D ε (x 1 * , λ ) satisfies w s = 0, s = 1, ..., S. Then we shall be able to conclude the proof for P R ε .

We can rewrite the constraint of D ε (x 1 * , λ ) as Lu = v, where L is the following integer matrix of rank ST + n,

For

  any feasible x, we define the following linear program R(x) called "Recourse Program":

Proposition 5 .

 5 If ε < 1 l (ST +n) M

α 1 x 1 *≤ α 1 x 1 * + α 2 x 2 *A 2 x 2 +

 11222 + α 2 x 2 * + βy s * + M ε + z * ≤ M.Therefore, since α and β are positive, M ε T t=1 w s * t ≤ M andw s * t ≤ ε, ∀t, s.(14)Let us notice that (x 2 * , z * , y 1 * , ..., y S * , w 1 * , ..., w S * ) is an optimal solution of P R ε (x 1 * ). Dualizing Constraints[START_REF] Kouvelis | Robust Discrete Optimization and its Applications[END_REF] in P R ε (x 1 * ), we obtainD ε (x 1 * ) max λ≥0 min x 2 ,z y 1 ,...,y S w 1 ,...,w S α 2 x 2 + z + By s + w s ≥ d s -A 1 x 1 * , s = 1, ..., S C 2 x 2 ≥ b -C 1 x 1 *

D ε (x 1 *A 2 x 2 +

 122 , λ ) min x 2 ,z y 1 ,...,y S w 1 ,...,w S e 1 ,...,e S ,fα 2 x 2 + z + By s + w s -e s = d s -A 1 x 1 * , s = 1, ..., S C 2 x 2 -f = b -C 1 x 1 *

  . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 1 :

 1 Average CPU time and number of iterations for p sites and T products

	p	T	CPU time (s) # iterations # "non feasible" iterations
	100 100	10	32	18
	100 1500	754	95	65
	300 500	735	87	44
	500 300	317	56	39
	500 500	2362	117	67
	500 800	3480	114	80
	600 700	4989	139	98
	800 100	118	46	25
	1500 100	714	89	41

u tr = (x 2 , y 1 , w 1 , e 1 , ..., y S , w S , e S , f )

First, assume that u is a basic optimal solution of D ε (x 1 * , λ ). There is a (ST + n) × (ST + n)-basic matrix E = (e ij ) such that

Then, for any (t, s), t ∈ {1, ..., T } and s ∈ {1, ..., S}, there is τ such that w s t = 1 det(E) (adj(E) v) τ . If (x 2 , z , y 1 , ..., y S , w 1 , ..., w S ) is an optimal solution of D ε (x 1 * , λ ), then (x 1 * , x 2 , z , y 1 , ..., y S , w 1 , ..., w S ) is an optimal solution of P R ε and, from (14) we have w s t ≤ ε, and so

As in Section 6.1, let l M = max i,j |l ij |. According to Hadamard's inequality, we have Assume now that u is a non basic optimal solution of D ε (x 1 * , λ ). By using the same arguments as at the end of proof of Proposition 4 we can prove that no extreme ray in the u decomposition can have w s t = 0 for some (s, t) with s ∈ {1, ..., S} and t ∈ {1, ..., T }.

Therefore, in any optimal solution of D ε (x 1 * , λ ) we have:

Finally, consider an optimal solution (x 1 * , x 2 * , z * , y 1 * , ..., y S * , w 1 * , ..., w S * ) of P R ε . There is λ * such that (λ * , x 2 * , z * , y 1 * , ..., y S * , w 1 * , ..., w S * ) is an optimal solution of D ε (x 1 * ) and thus (x 2 * , z * , y 1 * , ..., y S * , w 1 * , ..., w S * ) is an optimal solution of D ε (x 1 * , λ * ), and thus w s * = 0 for all s = 1, ..., S and the proof is completed.

) is an optimal solution of P R ε then w * = 0 and (x * , d * , y * ) is an optimal solution of P R.

Proof. Assume that ε < 1 l (ST +n) M (ST +n) (ST +n)/2 . From Proposition 5, w s = 0 for all s = 1, ..., S in any optimal solution of P R . Then w = 0 in any optimal solution of P R ε . Indeed, assume there is an optimal solution (x * , d * , y * , w * ) of P R ε such that w * > 0; then there is an optimal solution of P R ε with x = x * and w s > 0 for s corresponding to the worst scenario associated to x * : a contradiction. Now, let (x * , d * , y * , w * ) be an optimal solution of P R ε . Since w * = 0, (x * , d * , y * ) is a solution of P R with value v(P R ε ). Since, from Proposition 3, v(P R ε ) ≤ v(P R), (x * , d * , y * ) is an optimal solution of P R.

Annex 2: a general bound M of P R

We give here a general bound for P R, i.e. a value M which is greater than an optimal value of the robust problem if this problem has a solution. This theoretical bound can be easily improved according to each specific problem.

Assume that P R has at least one feasible x. We show how to compute M such that v(P R) ≤ M . As in Annex 1, we can rewrite P R as the following linear program :

Adding the slack variables, e s , f s , s = 1, ...S and g, we can rewrite the constraints of the program above as Lu = v with

, where we assume w.l.o.g. that all the coefficients are integer (see Section 6.1), u tr = (x, z, y 1 , e 1 , f 1 , ..., y S , e S , f S , g) and v tr = (0, ..., 0, d 1 , ..., d S , b).

Consider the relaxed problem ( P R) where the variables, denoted by ū, are in R + . Let ū * be an extreme point of ( P R). Then there exists a basic (S + ST + n) × (S + ST + n)-matrix L B of L and a basic vector ū * B such that

Since all the coefficients of L are integers,

Furthermore, since all the coefficients of adj(L B ) are determinants of (S + ST + n-1)×(S +ST +n-1)-matrices, we have, according to Hadamard's inequality: v ∞ .

However according to theorem 5.6 from [?, chap. 5.1], we can choose a feasible u for the robust problem such that u -ū * ∞ ≤ (S + ST + n)Ξ(L), where Ξ(L) is the maximum over all the sub-determinants of L, and which can be bound according to Hadamard 's inequality. Then v(P R) ≤ M = (p α ∞ + 1) u ∞ ≤ (p α ∞ + 1)( ū * ∞ + u -ū * ∞ ) .