domain: spi sde

text:
Xuan-Linh Dang
email: xuan-linh.dang@supelec.fr

Marc Petit

Paul Codani
email: paul.codani@mpsa.com

Keywords: Energy Management System, Micro-grid, Power Transformer, Smart Charging, Vehicle-To-Grid

An analysis of the operating conditions of a power transformer with the introduction of photovoltaic (PV) panels and Electric Vehicles (EV) in an eco-district is conducted. The study focuses on: (a) determining the optimal transformer sizing when there is neither PV nor EV; and (b) characterizing the overloading periods of the transformer induced by the penetration of PV sources and EVs in the district. We compare two EV integration strategies: they are either considered as unmanaged loads, or controlled by an Energy Management System (EMS) with Vehicle-to-Grid (V2G) capabilities. Transformer operating curves and rated power are deduced from technical standards supplied by Schneider Electric. The results show that EVs and PV introduction has a substantial impact on overloading periods. However, the latters are seriously mitigated both in terms of energy flows and peak power with the implementation of the EMS.

I. INTRODUCTION

In order to cope with the objectives of reduction in both air pollution and CO 2 emissions, the penetration of Renewable Energy Sources (RES) (such as photovoltaic (PV) panels, wind turbines, etc.) in power systems has substantially increased during the last few years. Similarly, plug-in Electric Vehicles (EV) moved by electric motors and powered by electrochemical batteries represent a promising solution with respect to these goals. With the upcoming decrease in battery costs, and the deployment of charging stations, EV sales are expected to increase within the next few years.

However, the increasing penetration of these two new units brings up concerns regarding their impacts on the electrical grid security. On one hand, RES are asynchronous and intermittent by nature, and distributed at the distribution level. They could trigger local congestion, frequency-and voltagerelated problems, as well as system wide balancing issues [START_REF] Sharma | System inertial frequency response estimation and impact of renewable resources in ERCOT Interconnection[END_REF]- [START_REF] O'sullivan | Studying the Maximum Instantaneous Non-Synchronous Generation in an Island SystemFrequency Stability Challenges in Ireland[END_REF]. On the other hand, if not managed properly, the massive This research is a part of the Energy Positive IT 2.0 project and benefit from the support of the Conseil General de lEssonne.

Paul CODANI is preparing a PhD with the Armand Peugeot Chair dedicated to hybrid technologies and economy of electro mobility, and is funded by PSA Peugeot Citroen. In this paper, we are concerned with the overloading periods of the transformer of an eco-district induced by the introduction of PV panels and EVs. Initially, such a district, pictured in Fig. 1, only has residential and commercial facilities. It is connected to the distribution grid through an oil-immersed power transformer.

First, we propose an approach to determine the optimal rated power (in kVA) of the transformer when there is neither PV nor EV penetration. Then, we assess the impacts of the introduction of PV units and EVs on the operating conditions of the transformer. At this point, EVs are considered as dumb loads, i.e. they abide by a charge-as-plugged strategy. Finally, we conduct the same analysis considering that EVs have bidirectional capabilities, i.e. they can inject power back to the grid, and that their charging patterns are controlled by an Energy Management System (EMS), whose aim is to mitigate the transformer overloading periods.

As underlined in [START_REF] Igualada | Optimal Energy Management for a Residential Microgrid Including a Vehicle-to-Grid System[END_REF], very little research has been conducted on the coupling of microgrid energy management with Vehicle-to-grid (V2G) strategies. The main contributions of this paper are:

• an approach for an optimal sizing of the rated power of a distribution transformer, based on its load profile over a whole year and on technical guidelines from a transformer manufacturer;

• a method to analyze the operating mode of the power transformer during the overloading periods; • an optimal Energy Management System for the microgrid of the eco-district connected to the distribution grid. It includes residential and commercial consumption units; and three-different types of EV fleets; • a V2G system taking into account different types of EV fleets with their users behaviors and a smart charging strategy for each EV fleet. The paper is organized as follows. Section II provides a description of the system understudy and all its components. In section III, the transformer is characterized including the optimal sizing and the definitions of the overloading assessment criteria and the EMS is described. Section IV features the simulation use cases and the associated results. Finally, section V is the conclusion.

II. ECO-DISTRICT ENERGY SYSTEMS DESCRIPTION

A. Residential consumption modeling

The residential consumption model is based on the study of Ian and Murray [START_REF] Richardson | Domestic electricity use: A high-resolution energy demand model[END_REF], in which a stochastic model is proposed to simulate the habitants behaviors. More specifically, the duration of electricity use in households is highly dependent on the timing of the occupants activities. Thus, the use of electric appliances, lighting, heating and domestic hot water within a residential household are also taken into in the model.

B. Commercial building modeling

A load model is built to simulate the electric consumption of the commercial building, taking into account the difference between working days and weekend; and typical load shape.

From commercial building data, it is observed that most buildings have a base load during the night time. In the early morning, the buildings electric demand increases for cooling/warming depending on the outdoor air temperature. This will trigger a short-lived load spike called morning start-up. Later, as people come to work in the building, the electricity consumption increases with the arrival of occupants. Its value varies daily, seasonally and yearly. After working time, the power consumption reduces quite significantly, due to the fact that electric appliances are successively turned off.

C. Electric Vehicle fleet modeling

All EVs are assumed to be full-electric vehicles. They have a 22kWh battery with a state of charge (SOC) in the range [20%, 90%]. These values are based on the optimal operation range of the battery. The EV fleet is shared into three different sub-fleets:

• EV fleet A: people living in the eco-district. These EV are typically plugged-in during the night, and leave in the morning to go to work and come back at late afternoon; • EV fleet B: people working in the eco-district. These EVs typically arrive in the morning in the eco-district, and leave in the late afternoon. They charge at work during the day;

• EV fleet C: a company fleet, for instance belonging to the mail services. These EVs are typically used from early morning to noon for the first-round deliveries, and from the early afternoon to the middle of the afternoon for the second-round deliveries.

Several driver behaviors are considered. First, the drivers of fleet B are divided into two categories: (a) the fleet of drivers having an ordinary behavior, that is charging at home during the night and participating in the EMS during the day; and (b) the fleet of drivers who aim at charging the most as they can at work, thus charging only the required energy for their morning trip during the night at home, and not participating in the EMS during the day.

Then, range anxiety is taken into account for all the drivers: when they inform their EV controller with their future driving needs (see section III-C2), we assume that they overestimate their future trip distances and that they always provide the EV controller with their yearly maximum distance traveled.

Based on [START_REF] Cgdd | Livre Vert sur les infrastructures de recharge ouvertes au public pour les vehicules decarbones[END_REF], the Electric Vehicle Supply Equipment (EVSE) powers are distributed according to table I. Home charging is mainly done at low power, while working charging stations are more equally distributed (although fast charging is still marginal). Finally, the company fleet has 22kW EVSEs only. The EV fleets model is dynamic and stochastic. EV fleets average distance trip (D), departure time (T d ), daily number trips (N) and seasonal energy consumption (E) are obtained from PSA internal data and from the CROME project results and from French government agencies [START_REF]Cross-border Mobility for EVs[END_REF], [START_REF]Les vehicules electriques en perspectives[END_REF]

D. PV production data

The photovoltaic generation data has been measured in the area of Paris over one full year.

III. TRANSFORMER AND EMS CHARACTERIZATIONS

In this section, all active and apparent powers are respectively provided in kW and kVA.

A. Optimal transformer sizing

In this part, the approach used to size the transformer is presented. It is based on the technical guidelines provided by a transformer manufacturer [START_REF] Schneider Electric | Guide de la distribution electrique basse tension et HTA[END_REF]. To optimize the rated power (in kVA) of the transformer, it is necessary to analyze its load profile at this step, induced only by the residential and commercial loads within a whole year. The net instantaneous power, P sub (t), at the substation during the j th day is calculated according to (1):

P j sub (t) = P j resi (t) + P j com (t) (1)
where P j resi (t) and P j com (t) are respectively the consumption of residential households and commercial buildings.

Then, the yearly maximum substation power, P sub max , is deduced from the daily maximum substation power as described by (2):

P sub max = max j max t∈j (P j sub (t)) with 1 ≤ j ≤ 365 (2)
From this value, the optimal rated power is determined according to the algorithm below: for all i such that 0 ≤ i ≤ n do The optimal apparent power of the transformer is determined based on a power factor of 0.928 in order to avoid penalties when supplying energy to the consumers [START_REF] Schneider Electric | Guide de la distribution electrique basse tension et HTA[END_REF]:

S opt = P opt size 0.928 (3)
To take into account the growth of the demand, a 25% security factor is added to get the rating power S rated . Then, we choose the commercial transformer with the closest power. Finally, the active rated power P rated is calculated with a low power factor (0.8) to be conservative.

S rated = 1.25 × S opt (4)

B. Transformer Operating Conditions Analysis

In this scenario, we introduce PV panels on the roof of the commercial facilities, and EVs which charge according to a "charge-as-plugged" strategy. In other words, the EV fleet is not controlled by any external operator, and drivers plug-in as soon as possible (reflecting a normal and intuitive behavior).

Based on the hypothesis detailed in section II-C, the consumption of all EVs, P EV (t), is calculated. Its value is computed on a minute-basis. It differs significantly from one day to the others and from one season to the other. Then, the instantaneous net power is described by ( 5):

(5) P sub (t) = P P V (t) -(P resi (t) + P com (t) + P EV (t))

To characterize the overloading periods, the average daily substation power and the cumulative durations of the overloading periods are also calculated according to the equations provided in the algorithm. Moreover, the energy E ex is determined by [START_REF] Clement-Nyns | The impact of vehicle-togrid on the distribution grid[END_REF] when the substation power is overloaded:

(6) E ex = n i=1 t k tj (|P i sub (t)|-P rated )dt ,
with |P i sub (t)|> P rated with [t j , t k ] the overloading time interval during i th the day.

C. Energy Management System

1) District Energy Management System: At each time interval, the power flow, P f low (t), between the micro-grid and the distribution grid (excluding the available EVs for the EMS) is calculated as following:

P f low (t) = P P V (t) -(P resi (t) + P com (t) + P EV nonEM S (t)) (7)
where P EV nonEM S (t) denotes the consumption of EVs that are not available for the EMS most of the time, because they need to charge at full power for transportation needs. Then, based on the power flexibility that the EVs can provide (see [START_REF]Cross-border Mobility for EVs[END_REF]), the EMS calculates the power P EV EM S (t) that has to be supplied by the available EV fleets. It depends on the current balance between supply and demand in the district as shown in Fig. 2. The mathematical formulation is described by [START_REF] Richardson | Domestic electricity use: A high-resolution energy demand model[END_REF]: where P charg (t) and P disch (t) denote respectively the charging and discharging power flexibilities that can be provided by the EVs.

P EV EM S (t) =      max(P charg (t), -P f low (t) -P rated ), cases A or B max(P charg (t), |P f low (t)|-P rated ), case C max(P disch (t), |P f low (t)|-P rated ), case D (8)
The objective of this EV power management is to maximize EVs charging power during the off-peak periods (cases A, B and C in Fig. 2) either from the PV (case A, B) or from the distribution grid (case C); and to minimize the discharging energy from available EVs when the demand has reached the peak load or when the generation is too low to supply to the total demand (case D in Fig. 2). Finally, the net power flow at the substation level is given by ( 9):

P sub (t) = P f low (t) + P EV EM S (t) (9)
2) EV fleet smart charging regime: Once the required power from the available EVs P EV EM S (t) has been determined by the EMS, it must be dispatched among those vehicles. If EV charging is ordered, we choose to first charge the EVs that have the biggest difference (in absolute value) between their current SOC and their required energy for next trip (SOC i req (t d )). On the contrary, we consecutively discharge the EVs that have the highest SOC.

To respect the future driving needs, each EV is able to decide for each moment whether to participate in the EMS or to charge quickly for transportation. Therefore, individual EVs provide the EMS with their available charging powers over the next timeframe as described by [START_REF]Cross-border Mobility for EVs[END_REF]:

P i charg (t) = -min(P i EV SE , SOC i max -SOC i (t) ∆t ) P i disch (t) = min(P i EV SE , SOC i (t) -SOC i req (t + ∆t) ∆t ) (10)
Where P i charg (t) and P i disch (t) are respectively the available charging and discharging powers of i th EV for the next time interval, P i EV SE is the EVSE power; SOC i (t), SOC i max SOC i req (t) denote respectively the current, the maximum and the required state of charge for the i th EV.

It is noticeable that we assume EVs to have a good estimation of their future trip distances and departure times.

This could be achieved in practice by either implementing a machine learning algorithm, the EV controller being able to estimate future driving needs depending on the previous ones, or by having the EV user directly informing the EV controller with its future driving needs (through, for instance, a smartphone application).

IV. RESULTS AND DISCUSSION

A. Simulation use cases

We consider that 800 people are living in the eco-district, and that 1000 persons are working there. Then, it is assumed that 20% of the living people and 10% of the working people own an electric vehicle (fleets A and B). Additionally, the surface area of the PV panels is supposed to be 3000 m 2 and there are 10 electric vehicles in the company fleet (the fleet C).

B. Optimal Transformer sizing

According to the steps detailed in III-A, we compute the rated active power, P rated = 252kW , in the case without PV panels nor EV.

Considering this rated power, we compare the overloading periods with the requirements from [START_REF] Schneider Electric | Guide de la distribution electrique basse tension et HTA[END_REF]. Results are provided in Fig. 3. The sizing method is efficient, and no overloading period excesses the recommended limits.

C. Transformer operating conditions with the introduction of PV and EVs

Then, according to the methodology explained previously, we successively add the PV panels and the EVs in the ecodistrict, first without any smart charging strategy, then with the EMS described in section III-C. Fig. 4 shows the monotonically decreasing maximal daily substation power (in absolute value) without EV and with EV (with and without EMS). It is noticeable that the EMS allows a 25% reduction of the substation peak power in comparison with the unmanaged charging scenario.

Finally, we compare the transformer overloading occurrences under the two scenarios (with and without EMS). Results are provided in Fig. 5. Operating conditions of the transformer during overloading occurrences, with respect to the limitations provided in [START_REF] Schneider Electric | Guide de la distribution electrique basse tension et HTA[END_REF] Because of the introduction of PV panels and unmanaged EVs in the eco-district, there are more than 30 occurrences where transformer overloading capability is exceeded (Fig. 5a). Conversely, with the implementation of the EMS, all overloading periods are within the acceptable limits(Fig. 5b). The results are corroborated with those provided in table III, which compares the two strategies in terms of absolute maximal daily substation power (P max global ), exceeded exchanged energy (E ex , (6)), cumulative overloading durations and average substation power (P sub ave ) during the overloading periods.

V. CONCLUSION

In this paper, the rated power of a distribution transformer is computed to be compliant with the thermal limits. Then the operating conditions with the introduction of renewable sources and electric vehicles are also studied.

The flexibility provided by the EVs is used to reduce the risks of large overloading in the transformer. The proposed EMS allows to increase the margin regarding the thermal limits. This margin can be used to accept more EVs without risks of critical overloading or to postpone transformer reinforcement. An additional study will allow to define the optimal relationship between the PV surface and the number of EV to reach the transformer thermal limits.

Fig. 1 .

 1 Fig. 1. System Overview

E

  (kWh/km) Winter: cw = 0.18 // Summer: cs = 0.13 a A.M. values, b P.M. values

1 :

 1 Build the overloading limitation curve: P lim sub (t) = f (t d ) 2: α ← 0 {Initialization of α} 3: sized ← false 4: while sized do 5: P size ← α × P sub max 6:

Fig. 2 .

 2 Fig. 2. Energy Management System strategy in the district

Fig. 3 .

 3 Fig. 3. Operating conditions of the transformer (no EV & no PV)

Fig. 4 .

 4 Fig. 4. Daily maximal substation power during the whole year with the different scenarios

  as described in table II. Then, trips distances and departure times are distributed according to Gaussian distributions with mean (µ) and standard deviations (σ) values provided in table II.


text:
Xuan-Linh Dang
email: xuan-linh.dang@centralesupelec.fr

Marc Petit

Paul Codani
email: paul.codani@mpsa.com

Energy optimization in an Eco-district with Electric Vehicles smart charging

Keywords: Energy Management System, Micro-grid, Smart Charging, Vehicle-To-Grid

The energy modeling of an eco-district, composed of residential and commercial facilities, is proposed. First, an optimal strategy aiming at sizing the rated power of the substation transformer is introduced. Then, considering the introduction of photovoltaic (PV) panels and electric vehicles (EV) in the district, the new operating conditions of the transformer are detailed. In order to reduce (a) the peak power of the transformer and (b) the cumulative durations of the transformer overloading periods, an Energy Management System (EMS) is implemented. The latter uses the district EVs, which have Vehicle-to-grid (V2G) capabilities, as flexibility providing units. An economics layer is added in order to compute the electricity costs considering the French pricing. This approach is used to deduce the optimal contracted power with the local system operator. Results show a reduction of 70% in overloading durations, 71% in average overloading power, and 17% in electricity costs with the EMS.

I. INTRODUCTION

Urbanization, population growth and climate change put increasing pressure on electricity consumption and CO 2 emission reductions. These challenges have fostered the development of Renewable Energy Sources (RES) (such as photovoltaic (PV) panels, wind turbines, etc.) in power systems. Their penetration rate has substantially increased during the last few years. In France, even though a slight decrease in growth rate was observed in 2013, the development of RES sources is expected to continue [START_REF] Rte | Bilan previsionnel de l'equilibre offre-demande d'electricite en France[END_REF]. Similarly, plug-in Electric Vehicles (EV) appear as promising solutions to deal with the aforementioned issues. EV sales went up by 60% in Europe in 2014, and by around 20% in France for privately-owned EVs [START_REF] Avere | No Title[END_REF]. The expected trends in battery cost reduction and charging station deployment should strengthen the EV sales in the near future.

However, the increasing penetration of these two new units brings up concerns regarding their impacts on the electrical This research is a part of the Energy Positive IT 2.0 project and benefit from the support of the Conseil General de l'Essonne Paul Codani benefits from the support of the Chair "PSA Peugeot Citroen Automobile: Hybrid technologies and Economy of Electromobility", so-called Armand Peugeot Chair led by Ecole CentraleSupelec, ESSEC and sponsored by PEUGEOT CITROEN Automobile. grid security. On one hand, RES are asynchronous and intermittent by nature, and distributed at the distribution level. They could trigger local congestion, frequency and voltagerelated problems, as well as system wide balancing issues [START_REF] Sharma | System inertial frequency response estimation and impact of renewable resources in ERCOT Interconnection[END_REF]- [START_REF] O'sullivan | Studying the Maximum Instantaneous Non-Synchronous Generation in an Island SystemFrequency Stability Challenges in Ireland[END_REF]. On the other hand, if not managed properly, the massive introduction of plug-in vehicles could jeopardize grid security [START_REF] Darabi | Aggregated Impact of Plug-in Hybrid Electric Vehicles on Electricity Demand Profile[END_REF]- [START_REF] Clement-Nyns | The impact of vehicle-togrid on the distribution grid[END_REF].

In this paper, the authors focus on the electric and economics impacts induced by the introduction of PV and EVs in an eco-district, pictured in Fig. 1. As underlined in [START_REF] Igualada | Optimal Energy Management for a Residential Microgrid Including a Vehicle-to-Grid System[END_REF], very little research has been conducted on the coupling of microgrid energy management with Vehicle-to-grid (V2G) strategies.

This paper builds on [START_REF] Dang | Transformer Operating conditions under introduction of PV and EVs in an Eco-district[END_REF]. In this previous work, the authors proposed an approach to determine the optimal rated power (in kVA) of the transformer when there was neither PV nor EV penetration. Then, considering the introduction of PV and EVs, the operating conditions of the transformer were characterized with and without the implementation of an Energy Management System using EVs as flexibility providing units.

In this paper, we go a step further and propose the following novel contributions:

• an economics analysis of the electricity costs of the ecodistrict, considering the French price scale fixing; • a method to optimally choose the contracted power of the transformer with respect to the associated costs;

• a sensitivity analysis on the EV penetration rate in the district. The authors would like to draw the reader's attention on the fact that we will distinguish two characteristic powers for the substation transformer:

• the rated power P rated , which is a physical feature of the transformer. This power is used in the energy analysis of the transformer operating conditions. • the contracted power P c , which has no physical meaning; it is the power that has been contracted between the ecodistrict energy operator and the system operator (SO). This power is used for the economics considerations, as peak power values are evaluated by the SO with respect to this value. The paper is organized as follows. The system description is recalled in section II. The transformer sizing method, and the EMS description are then provided in section III. In section IV, the electricity French pricing scales for ecodistrict-type units are explained. Finally, section V and VI provide respectively the results and the conclusion.

II. ECO-DISTRICT ENERGY SYSTEMS DESCRIPTION

A. Residential consumption modeling

The residential consumption model is based on the study of Ian and Murray [START_REF] Richardson | Domestic electricity use: A high-resolution energy demand model[END_REF], in which a stochastic model is proposed to simulate the habitants behaviors. More specifically, the duration of electricity use in households is highly dependent on the timing of the occupants activities. Thus, the use of electric appliances, lighting, heating and domestic hot water within a residential household are also taken into in the model.

B. Commercial building modeling

A load model is built to simulate the electric consumption of the commercial building, taking into account the difference between working days and weekend; and typical load shape.

From commercial building data, it is observed that most buildings have a base load during the night time. In the early morning, the buildings electric demand increases for cooling/warming depending on the outdoor air temperature. This will trigger a short-lived load spike called morning start-up. Later, as people come to work in the building, the electricity consumption increases with the arrival of occupants. Its value varies daily, seasonally and yearly. After working time, the power consumption reduces quite significantly, due to the fact that electric appliances are successively turned off.

C. Electric Vehicle fleet modeling

All EVs are assumed to be full-electric vehicles. They have a 22kWh battery with a state of charge (SOC) in the range [20%, 90%]. These values are based on the optimal operation range of the battery. The EV fleet is shared into three different sub-fleets:

• EV fleet A: people living in the eco-district. These EV are typically plugged-in during the night, and leave in the morning to go to work and come back at late afternoon;

• EV fleet B: people working in the eco-district. These EVs typically arrive in the morning in the eco-district, and leave in the late afternoon. They charge at work during the day; • EV fleet C: a company fleet, for instance belonging to the mail services. These EVs are typically used from early morning to noon for the first-round deliveries, and from the early afternoon to the middle of the afternoon for the second-round deliveries. Several driver behaviors are considered. First, the drivers of fleet B are divided into two categories: (a) the fleet of drivers having an ordinary behavior, that is charging at home during the night and participating in the EMS during the day; and (b) the fleet of drivers who aim at charging the most as they can at work, thus charging only the required energy for their morning trip during the night at home, and not participating in the EMS during the day.

Then, range anxiety is taken into account for all the drivers: when they inform their EV controller with their future driving needs (see section III-B2), we assume that they overestimate their future trip distances and that they always provide the EV controller with their yearly maximum distance traveled.

Based on [START_REF] Cgdd | Livre Vert sur les infrastructures de recharge ouvertes au public pour les vehicules decarbones[END_REF], the Electric Vehicle Supply Equipment (EVSE) powers are distributed according to table I. Home charging is mainly done at low power, while working charging stations are more equally distributed (although fast charging is still marginal). Finally, the company fleet has 22kW EVSEs only.

TABLE I EVSE REPARTITIONS

EVSE power plug

Fleet A Fleet B Fleet C

Slow (a) -3kW 93% 35% 0% Slow (b) -7kW 7% 34% 0% Intermediate charging -22kW 0% 29% 100% Fast charging -43kW 0% 2% 0%
The EV fleets model is dynamic and stochastic. EV fleets average distance trip (D), departure time (T d ), daily number trips (N) and seasonal energy consumption (E) are obtained from PSA internal data and from the CROME project results and from French government agencies [START_REF]Cross-border Mobility for EVs[END_REF], [START_REF]Les vehicules electriques en perspectives[END_REF] as described in TABLE II. Then, trips distances and departure times are distributed according to Gaussian distributions with mean (µ) and standard deviations (σ) values provided in TABLE II.

D. PV production data

The photovoltaic generation data has been measured in the area of Paris over one full year.

III. TRANSFORMER SIZING AND EMS CHARACTERIZATION

A. Transformer sizing

Only the residential and commercial consumptions are considered when sizing the rated power of the transformer, P rated . With respect to thermal limitations provided by a transformer manufacturer [START_REF] Schneider Electric | Guide de la distribution electrique basse tension et HTA[END_REF], we design an algorithm which will determine the optimal rated power of the transformer, while still allowing some overloading occurrences. The final choice for the transformer depends on the typical rated power proposed by the manufacturers. In the present case S rated = 315 kVA is chosen, with a 0.8 power factor, then P rated = 252kW. The transformer operating conditions after the sizing are presented in Fig. 2. More details on the sizing algorithm can be found in [START_REF] Dang | Transformer Operating conditions under introduction of PV and EVs in an Eco-district[END_REF]. B. EMS description 1) District Energy Management System: At each time interval, the power flow, P f low (t), between the micro-grid and the distribution grid (excluding the available EVs for the EMS) is calculated as following:

P f low (t) = P P V (t) -P resi (t) + P com (t) + P EV nonEM S (t) (1)
where P EV nonEM S (t) denotes the consumption of EVs that are not available for the EMS most of the time, because they need to charge at full power for transportation needs. Then, based on the power flexibility that the EVs can provide (see (4)), the EMS calculates the power P EV EM S (t) that has to be supplied by the available EV fleets. It depends on the current balance between supply and demand in the district as shown in Fig. 3. The mathematical formulation is described by (2):

P EV EM S (t) =      max(P charg (t), -P f low (t) -P rated ), cases A or B max(P charg (t), |P f low (t)|-P rated ), case C min(P disch (t), |P f low (t)|-P rated ), case D (2)
where P charg (t) and P disch (t) denote respectively the charging and discharging power flexibilities that can be provided by the EVs.

The objective of this EV power management is to maximize EVs charging power during the off-peak periods (cases A, B and C in Fig. 3) either from the PV (case A, B) or from the distribution grid (case C); and to minimize the discharging energy from available EVs when the demand has reached the peak load or when the generation is too low to supply to the total demand (case D in Fig. 3). Finally, the net power flow at the substation level is given by (3):

P sub (t) = P f low (t) + P EV EM S (t) (3)
2) EV fleet smart charging regime: Once the required power from the available EVs P EV EM S (t) has been determined by the EMS, it must be dispatched among those vehicles. If EV charging is ordered, we choose to first charge the EVs that have the biggest difference (in absolute value) between their current SOC and their required energy for next trip (SOC i req (t d )). On the contrary, we consecutively discharge the EVs that have the highest SOC.

To respect the future driving needs, each EV is able to decide for each moment whether to participate in the EMS or to charge quickly for transportation. Therefore, individual EVs provide the EMS with their available charging powers over the next timeframe as described by (4):

P i charg (t) = -min P i EV SE , SOC i max -SOC i (t) ∆t P i disch (t) = min P i EV SE , SOC i (t) -SOC i req (t + ∆t) ∆t (4)
Where P i charg (t) and P i disch (t) are respectively the available charging and discharging powers of i th EV for the next time interval, P i EV SE is the EVSE power; SOC i (t), SOC i max SOC i req (t) denote respectively the current, the maximum and the required state of charge for the i th EV.

It is noticeable that we assume EVs to have a good estimation of their future trip distances and departure times. This could be achieved in practice by either implementing a machine learning algorithm, the EV controller being able to estimate future driving needs depending on the previous ones, or by having the EV user directly informing the EV controller with its future driving needs (through, for instance, a smartphone application).

IV. ECONOMICS CONSIDERATIONS

We consider the French electricity price scale fixing to compute the electricity costs in the eco-district. All the details on the French pricing are provided in [START_REF] Erdf | Tarif d'Utilisation du Reseau Public de Distribution d'Electricite (TURPE)[END_REF]. Considering the value of the transformer rated power (P rated = 252kW , see section III-A), the district is connected at the medium voltage level (so-called HTA level in France). The cost structure is twofold:

• an energy component C E which is related to the electricity consumption/production (kWh) • a capacity component C P which is related to the contracted power C P (kW)

These different costs are all dependent on the season and time during which the energy exchanges occur. We subscribe to a pricing structure with five different temporal classes; they are described in TABLE III. We also provide the power and energy coefficients that will be used in the power and energy components calculation. For each temporal class i, a contracted power P ci has to be chosen, while still abiding by the constraint expressed in equation ( 5):

∀i ∈ {1; 4}, P ci ≤ P ci+1 (5)
It is noticeable that we do not take into account the metering and customer management costs, which are other components included in the cost structure but would not be impacted by the energy management strategy.

A. Capacity component calculation

The capacity component C P is itself divided into two subcomponents: the yearly capacity subscription cost C sub , and the monthly overloading costs C over . Both components are computed according to equation ( 6):

C sub = a 2 ×   k 1 × P c1 + 5 i=2 k i × (P ci -P ci-1 )   C over = 0.15 × k i × a 2 × ∆P 2 (6)
with the constant a 2 = 9.24 (e/kW/year), P c1→5 and k 1→5 respectively the contracted powers and power coefficients for the various temporal classes, and ∆P the difference between the substation power and the contracted one (only when the substation power exceeds the contracted one).

The results presented in section V-B showing the optimal contracted power result from an optimization algorithm, which minimizes the total capacity costs C P while abiding by constraint (5).

B. Energy component calculation

The energy cost component C E is computed according to equation [START_REF] Green | The impact of plug-in hybrid electric vehicles on distribution networks: A review and outlook[END_REF]:

C E = 5 i=1 d i × E i (7)
with d 1→5 the energy coefficients (see TABLE III) and E i the net energy exchanged during temporal class i. It is noticeable that we assume grid parity in the energy cost calculation, i.e. energy is sold at the same price as it is bought.

C. Capacity costs optimization

We formulate an optimization problem, which objective function is the total capacity costs C P :

min Pc 1→5 C P (8)
The problem is subject to the constraint expressed in equation [START_REF] O'sullivan | Studying the Maximum Instantaneous Non-Synchronous Generation in an Island SystemFrequency Stability Challenges in Ireland[END_REF]. The power and energy flows between the district and the grid are given as parameters; they are computed for the baseline scenario, once with the EMS implemented and once without.

V. RESULTS AND DISCUSSIONS

A. Baseline simulation scenario

The baseline scenario is the following. We consider that 800 people are living in the eco-district, and that 1000 persons are working there. Then, it is assumed that 20% of the living people and 10% of the working people own an electric vehicle (fleets A and B). Additionally, the surface area of the PV panels is supposed to be 3000 m 2 and there are 10 electric vehicles in the company fleet (the fleet C). The total number of EVs per fleet is provided in TABLE VII.

B. Economics results

The results from the optimization problem described in section IV-C are shown in TABLE IV. Without any EMS, the transformer peak powers occur during the critical peak tariff periods. As a consequence, and because of constraint ( 5), all the contracted power values are set by the critical peak contracted power (P c1 ). On the other hand, thanks to the EMS, it is possible to have lower hightariff contracted powers than low-tariff contracted powers. Moreover, the absolute values of the contracted powers are significantly reduced. These results have an impact on the electricity costs as shown in TABLE V. As the contracted power values are much lower with the EMS, the capacity costs are reduced by 24%. Overall, the savings in electricity payment are reduced by 17%.

C. Transformer overloading 1) Baseline scenario:

We analyze the transformer operating conditions under the baseline scenario, with and without the EMS. Fig. 4 shows the monotonically decreasing maximal daily substation power (in absolute value) without EV and with EV (with and without EMS). It is noticeable that the EMS allows a 25% reduction of the substation peak power in comparison with the unmanaged charging scenario. Results show a drastic improvement in the transformer operating conditions with the implementation of the EMS, with a reduction of 70% in the cumulative overloading durations, and a reduction of 90% in the energy exchanged during overloading periods.

2) EV penetration rate: sensitivity analysis: Finally, we aim at conducting a sensitivity analysis on the EV penetration rate, and assessing the impact of the latter on the transformer operating conditions. This analysis could provide insights on the ability of the EMS to accept more ore less EV charging periods, depending on the EV sales' evolution. We consider three scenarios in terms of EV numbers in each fleet. They are summarized in TABLE VII. Fig. 5a and5b show respectively the distribution function of the overloading average powers and overloading durations.

The EMS performs well in terms of overloading average powers. The EV penetration rate only has a small impact on this parameter. On the other side, the impacts on the overloading durations is more significant. As indications, the transformer rated power that would be required, if there was no EMS, are provided in TABLE VIII (with the EMS, in the main scenario, there is no need to change the rated power (P rated = 252kW )).

These results could be used in future work to compute savings in investment costs. VI. CONCLUSION In this paper, the energy modeling of an eco-district is presented. First, only residential and commercial consumptions are considered, and an optimal sizing for the substation transformer is determined. Then, PV panels and EVs are introduced in the district. In order to limit the impacts both in terms of transformer operating conditions and electricity costs, and Energy Management System is proposed. Results show a reduction of 70% in overloading durations, 71% in average overloading power, and 17% in electricity costs with the EMS.

Overall, yearly electricity savings amount to 3370e(TABLE V), considering 270 flexible EVs (baseline scenario A). Thus, a direct remuneration of the EV end users does not appear as a very suitable business model, as drivers would only earn 12e/year. Policy recommendations, such as changes in electricity pricing, driven by the environmental gains of this solution could foster the developments of such EMS. In a next work, the impact of single phase 7kW chargers in residential parkings will be analysed because there seems to be a consensus for installing them. Their cost is barely more expensive than a 3 kW charge.

Fig. 1 .

 1 Fig. 1. System Overview

Fig. 2 .

 2 Fig. 2. Operating conditions of the transformer (no EV & no PV)

Fig. 3 .

 3 Fig. 3. Energy Management System strategy in the district

Fig. 4 .

 4 Fig. 4. Daily maximal substation power during the whole year with the different scenarios

Fig. 5 .

 5 Fig.5. Operating conditions of the transformer during overloading occurrences, with respect to the limitations provided in[START_REF] Schneider Electric | Guide de la distribution electrique basse tension et HTA[END_REF]

TABLE III TEMPORAL

 III CLASSES CONSIDERED FOR ELECTRICITY PRICING[START_REF] Erdf | Tarif d'Utilisation du Reseau Public de Distribution d'Electricite (TURPE)[END_REF]

	Class number	Corresponding	Power coeff	Energy coeff
		period of time	k i	d i (ce/kWh)
	1	critical peak	100%	2.98
	2	on-peak winter	92%	2.56
	3	off-peak winter	55%	1.53
	4	on-peak summer	40%	1.30
	5	off-peak summer	12%	0.87

Table

  VI details the comparison of the energy efficiency with and without EMS, in terms of absolute maximal daily substation power (P max global ), exceeded exchanged energy (E ex ), cumulative overloading durations and average substation power (P sub ave ) during the overloading periods. The energy E ex is determined by (9) when the substation power is overloaded:

		n	t k
	E ex =	i=1	tj	|P i sub (t)|-P rated dt ,	(9)
	with |P i sub (t)|> P rated

TABLE VIII TRANSFORMER

 VIII RATED POWER THAT WOULD BE REQUIRED WITHOUT EMS

	Scenario	Apparent power	Active power
		S rated (kVA)	P rated (kW)
	Pessimistic Scenario (P)	400	320
	Baseline Scenario (A)	500	400
	Optimistic Scenario (O)	630	504


Céline Gicquel

Laurence A Wolsey

Michel Minoux

On Discrete Lot Sizing and Scheduling on Identical Parallel Machines

Keywords:

We consider the multi-item discrete lot-sizing and scheduling problem on identical parallel machines. Based on the fact that the machines are identical, we introduce aggregate integer variables instead of individual variables for each machine. For the problem with start-up costs, we show that the inequalities based on a unit flow formulation for each machine can be replaced by a single integer flow formulation without any change in the resulting LP bound. For the resulting integer lot-sizing with start-ups subproblem, we show how inequalities for the unit demand case can be generalized and how an approximate version of the extended formulation of Eppen and Martin can be constructed. The results of some computational experiments carried out to compare the effectiveness of the various mixed-integer programming formulations are presented.

Introduction

In the present paper, we consider the multi-item discrete lot-sizing and scheduling problem on identical parallel machines.

The case with one machine, known as the discrete lot-sizing and scheduling problem (DLSP) has the particularity that at most one item is produced per period, that production of that item (if any) is at full capacity and that start-up or changeover costs have to be incurred when switching production from one item to another. A natural subproblem of this model is the single item discrete lot-sizing problem with start-ups, where there are linear production and storage costs and fixed production and start-up costs. For this Van Eijl and van Hoesel [START_REF] Van Eijl | On the discrete lot-sizing and scheduling problem with Wagner-Whitin costs[END_REF] have presented valid inequalities and van Hoesel and Kolen [START_REF] Van Hoesel | A linear description of the discrete lot-sizing and scheduling problem[END_REF] have provided a tight and compact extended formulation.

A different single-item parallel resource model with start-ups was used by Lasdon and Terjung [START_REF] Lasdon | Optimization theory for large systems[END_REF][START_REF] Lasdon | An efficient algorithm for multi-item scheduling[END_REF]. Here the demands are integers, the amount produced is integral and there are linear production, storage and start-up costs. Eppen and Martin [START_REF] Eppen | Solving multi-item capacitated lot-sizing problems using variable redefinition[END_REF] have given a tight extended formulation and Vanderbeck and Wolsey [START_REF] Vanderbeck | Valid inequalities for the Lasdon-Terjung production model[END_REF] have derived several families of valid inequalities. The multi-item problem with several machines has also been studied recently by Gicquel et al. [START_REF] Gicquel | A tight MIP formulation for the discrete lot sizing and scheduling problem with parallel resources[END_REF].

The purpose of this paper is threefold. We first introduce in Section 2 two formulations of the problem, using either machine-specific binary variables or aggregate integer variables, and provide equivalence results between these two formulations. In Section 3 we propose two ways of strengthening the second formulation involving the aggregate integer variables. On the one hand, we develop an approximate extended formulation based on the extended formulation proposed in [START_REF] Eppen | Solving multi-item capacitated lot-sizing problems using variable redefinition[END_REF], and on the other we derive a new family of valid inequalities, that can be seen as an extension to the parallel-machine case of the valid inequalities proposed in [START_REF] Van Eijl | On the discrete lot-sizing and scheduling problem with Wagner-Whitin costs[END_REF] for the single-machine problem. Finally in Section 4 we provide some computational results carried out to compare the various formulations of the problem. Our results indicate that it is possible to solve instances with 10 or more items, up to 10 machines and up to 50 periods to optimality, and that for instances with up to 150 periods, solutions within 5% of optimality can be obtained in less than 10 minutes.

Problem Formulations and Equivalence Results

We first give a formulation for the problem with start-up costs and machine specific 0-1 variables. There are I items, T periods and K identical machines. As the machines are identical and produce just one item at full capacity in a period, the demand is measured in units of full capacity. The data consists of the demands d i t for item i in period t, the (machine-independent) production cost p i t of producing item i on one of the machines in period t, the (machine-independent) start-up cost q i t of starting to produce item i on a machine in period t and h i t the storage cost for item i at the end of period t. Note that one can assume without loss of generality that d i t ∈ {0, 1, . . . , K}. Also an additional item, item 0, is introduced to represent the case when a machine is idle.

The variables are: -y ik t = 1 if item i ∈ {0, . . . , I} is produced on machine k in period t, and y ik t = 0 otherwise -z ik t = 1 if production of item i on machine k is started in period t, and z ik t = 0 otherwise, -s i t is the stock of item i at the end of period t.

The corresponding formulation is:

v BIN = min I i=1 K k=1 T t=1 (p i t y ik t + q i t z ik t ) + I i=1 T t=1 h i t s i t ( 1
)
s i t-1 + K k=1 y ik t = d i t + s i t ∀ i, t (2)
z ik t ≥ y ik t -y ik t-1 ∀ i, k, t (3)
z ik t ≤ y ik t ∀ i, k, t (4)
y ik t-1 + z ik t ≤ 1 - j:j =i (y jk t -z jk t ) ∀ i, k, t (5)
I i=0

y ik t = 1 ∀ k, t (6)
s ∈ R IT + , y, z ∈ {0, 1} (I+1)KT (7)
Note that inequality ( 5) is a strengthened version of the more standard inequality

y ik t-1 +z ik t ≤ 1, see Constantino [2]. Let Q yz be the polyhedron {(y, z) ∈ R (I+1)KT + × R (I+1)KT + : (3) -(6)} and X yz = Q yz ∩ Z (I+1)KT × Z (I+1)KT .
We now introduce the aggregate variables

Y i t = K k=1 y ik t and Z i t = K k=1 z ik t ( 8
)
and the resulting aggregated formulation

v IN T = min I i=1 T t=1 (p i t Y i t + q i t Z i t + h i t s i t ) (9)
s i t-1 + Y i t = d i t + s i t ∀ i, t (10)
Z i t ≥ Y i t -Y i t-1 ∀ i, t (11)
Z i t ≤ Y i t ∀ i, t (12)
Y i t-1 + Z i t ≤ K - j:j =i (Y j t -Z j t ) ∀ i, t (13)
I i=0 Y i t = K ∀ t (14) s ∈ R IT + , y, z ∈ Z (I+1)T + . ( 15
)
Note that the slack variable in [START_REF] Van Eijl | On the discrete lot-sizing and scheduling problem with Wagner-Whitin costs[END_REF],

Ω i t-1 = Z i t -Y i t + Y i t-1
, can be seen as the number of machines turned off in period t -1.

As above, let

Q Y Z be the polyhedron {(Y, Z) ∈ R (I+1)T + × R (I+1)T + : (11) -(14)} and X Y Z = Q Y Z ∩ Z (I+1)T × Z (I+1)T .
It is easily verified that if (y, z) ∈ Q yz and Y, Z are defined as in [START_REF] Lasdon | An efficient algorithm for multi-item scheduling[END_REF], then (Y, Z) ∈ Q Y Z . Below we establish the converse.

Theorem i) If (Y, Z) ∈ Q Y Z /X Y Z , then there exists (y, z) ∈ Q yz /X yz satisfying (8), ii) Q yz = conv(X yz ) and Q Y Z = conv(X Y Z ).

Proof

The proof involves an extended formulation for Q yz and Q Y Z . As in Belvaux and Wolsey [START_REF] Belvaux | Modelling practical lot-sizing problems as mixed-integer programs[END_REF], we define the changeover variables as w ijk t = 1 if item i is produced on machine k in period t -1 and item j is produced in period t and w ijk t = 0 otherwise, and the polyhedron P yzw :

i w ijk t = y jk t ∀ j, k, t (16)
j w ijk t = y ik t-1 ∀ i, k, t (17)
i:i =j

w ijk t = z jk t ∀ j, k, t (18)
i

y ik 1 = 1 ∀ k (19) w ∈ R (I+1) 2 KT + , y, z ∈ R (I+1)KT + . ( 20
)
We also introduce the aggregate changeover variables

W ij t = K k=1 w ijk t (21)
leading to the aggregated polyhedron P Y ZW :

i W ij t = Y j t ∀ j, t (22)
j W ij t = Y i t-1 ∀ i, t (23)
i:i =j

W ij t = Z j t ∀ j, t (24)
i Y i 1 = K (25) W ∈ R (I+1) 2 T + , Y, Z ∈ R (I+1)T + . ( 26
)
Constantino [2] has shown that Q yz is the projection of P yzw onto the space of the y, z variables. The argument uses the max-flow/min cut theorem applied to ( 16)-( 17) viewed as a transportation network. As P yzw is an integral polytope, this implies that Q yz is also an integral polytope. An identical argument holds for P XY Z and Q Y Z . It remains to show that if (Y, Z, W ) ∈ P Y ZW , then there exists a decomposed machine specific solution (y, z, w) ∈ P yzw satisfying the aggregation equations ( 8) and (21). However ( 22),( 23),( 26) The result follows as now any (Y, Z) ∈ Q Y Z corresponds to a (Y, Z, W ) ∈ P Y ZW which, as just shown, corresponds to a (y, z, w) ∈ P yzw and hence to a (y, z) ∈ Q yz , and integer solutions also lead to integer solutions.

It follows that the aggregated and disaggregated formulations are equivalent both as integer and linear programs.

Let v BIN LP , v IN T LP
denote the values of the linear programming relaxations of ( 1)-( 7) and ( 9)-( 15) respectively.

Corallary v IN T = v BIN and v IN T LP = v BIN LP .
The potential benefits from use of the aggregated formulation are i) its smaller size and ii) the fact that the most obvious problems of symmetry of the solutions of the disaggregated model are removed. Note also that whereas complete enumeration of the y vectors in the disaggregated model gives (2 K ) (I+1)T distinct vectors, the number of Y vectors in the aggregated model is (K + 1) (I+1)T .

Strengthening the Aggregate Formulation

Here we consider ways to strengthen the formulation ( 9)-(15), and in particular the single item set X LT :

s t-1 + Y t = d t + s t ∀ t (27) Z t ≥ Y t -Y t-1 ∀ t (28) Z t ≤ Y t ∀ t (29) Y, Z ∈ {0, 1, • • • , K} T (30)
treated by Lasdon-Terjung [START_REF] Lasdon | An efficient algorithm for multi-item scheduling[END_REF].

From now on we will use the notation

d ut ≡ t j=u d j , Y ut ≡ t j=u Y j , etc.
, with a superscript i to denote the item when appropriate.

Two ways to strengthen this formulation have been proposed. Eppen and Martin [START_REF] Eppen | Solving multi-item capacitated lot-sizing problems using variable redefinition[END_REF] gave an extended formulation for conv(X LT ) with O(K 2 T d 1T ) variables and O(KT d 1T ) constraints. Vanderbeck and Wolsey [START_REF] Vanderbeck | Valid inequalities for the Lasdon-Terjung production model[END_REF] gave three families of valid inequalities along with separation heuristics.

Here we give an approximate version of the extended formulation, and also show how a family of inequalities developed by Van Eijl and Van Hoesel [START_REF] Van Eijl | On the discrete lot-sizing and scheduling problem with Wagner-Whitin costs[END_REF] can be adapted to give valid inequalities for this model.

An Approximate Extended Formulation

We modify the extended formulation given by Eppen and Martin [START_REF] Eppen | Solving multi-item capacitated lot-sizing problems using variable redefinition[END_REF]. The approximation is obtained by aggregating nodes and arcs whenever the stock level exceeds some user specified value ∆.

We now describe the construction, given an approximation parameter ∆ and the fact Y t is bounded by the number of machines K. One constructs a directed graph D = (V, A) with nodes (t, α, β) where t represents the time period, α represents Y 1t and β represents Y t for d 1t ≤ α ≤ d 1t + ∆ and nodes (t, α) for d 1t + ∆ < α ≤ d 1T . The arcs (t, α, β, γ) depart from node (t, α, β) and arrive at either node (t+1, α +γ, γ) or node (t+1, α +γ) with corresponding variables u(t, α, β, γ). The second set of arcs depart from node (t, α) and arrive at either node (t + 1, α + γ, γ) or node (t + 1, α + γ) with corresponding variables v(t, α, γ). See Figure 2.

Thus the arc variables are defined as follows:

-u(t, α, β, γ) = 1 if Y 1t = α, Y t = β, Y t+1 = γ, and 0 otherwise with d 1t ≤ α ≤ min(d 1t + ∆, d 1T ), d 1T ≥ α + Y t+1 ≥ d 1,t+1 , Y t , Y t+1 ∈ {0, • • • , K}, and -v(t, α, γ) = 1 if Y 1t = α, Y t+1 = γ, and 0 otherwise with d 1T ≥ α > d 1t + ∆, d 1T ≥ α + Y t+1 ≥ d 1,t+1 , Y t , Y t+1 ∈ {0, • • • , K}.
The resulting "approximate dynamic programming" formulation is:

u(0, 0, 0) = 1 (31) γ u(t, α, β, γ) = γ u(t -1, α -β, γ, β) + v(t -1, α -β, β) ∀nodes (t, α, β) (32) γ v(t, α, γ) = β γ u(t -1, α -β, γ, β) + v(t -1, α -β, β) ∀nodes (t, α) (33)
Y t = α,γ γ β u(t -1, α, β, γ) + v(t -1, α, γ) ∀ t (34) Z t ≥ α,β,γ (γ -β) + u(t, α, β, γ) ∀ t (35) Z t ≥ Y t -Y t-1 ∀ t (36) Z t ≤ Y t ∀ t (37) Y, Z ∈ [0, K] T , u(t, α, β, γ) ∈ R 1 + ∀ t, α, β, γ, v(t, α, β) ∈ R 1 + ∀ t, α, β, (38)
where D = max t d t . Note that when ∆ = d 1T , the v variables and the equations (33) disappear and one obtains the original extended formulation from [START_REF] Eppen | Solving multi-item capacitated lot-sizing problems using variable redefinition[END_REF] with equality in (35).

Figure 2: The variables in the approximate extended formulation

A Family of Valid Inequalities

We show that the valid inequalities of Van Eijl and Van Hoesel [START_REF] Van Eijl | On the discrete lot-sizing and scheduling problem with Wagner-Whitin costs[END_REF] derived for the case where Y t ∈ {0, 1} and d t ∈ {0, 1} for all t are also valid for X LT . Specifically we break the integer demands d t into demands of one unit. Thus given an interval [t, l], we say that the q th unit of demand occurs in period t q if d t,tq-1 < q and d t,tq ≥ q. Thus t ≤ t

1 ≤ • • • ≤ t q ≤ • • • ≤ t p ≤ l.

Proposition

The inequality

s t-1 + p q=1 (Y t+q-1 + Z t+q,tq ) ≥ p (39)
is valid for X LT , where d tl = p and t + q < t q for q = 1, . . . , p.

Proof

Let α(t, p) be the expression on the left of the inequality. We use induction.

Initial
+ Y t,t+p-1 = s t-1 + Y t,tp+1 ≥ d t,tp+1 ≥ p + 1,
where the first inequality follows from the flow conservation equations ( 27) and the second from the definition of t p+1 . However α(t, p + 1) ≥ s t-1 + Y t,t+p-1 and the claim follows.

Observe that the inequality (39) is dominated by

s t-1 + Y t,tp ≥ d t,tp if t + p ≥ t p .

Numerical results

We now present the results of some computational experiments carried out to compare the effectiveness of the various mixed-integer programming formulations discussed in sections 2 and 3 in solving the problem under study.

Problem instance generation

We randomly generated instances of the problem using a procedure adapted from that described in [START_REF] Salomon | Solving the discrete lotsizing and scheduling problem with sequence dependent set-up costs and set-up times using the travelling salesman problem with time windows[END_REF] for the multi-item DLSP with a single machine. More precisely, the various instances tested differ with respect to the following characteristics:

-Problem dimension: The problem dimension is represented by the number of items I, the number of periods T and the number of machines K. We define 4 sets of instances:

- -Production capacity utilization: Production capacity utilization ρ is defined as the ratio between the total cumulated demand ( N i=1 d i 1T ) and the total cumulated available capacity (K × T ). ρ was varied between 0.75 and 0.95, in steps of 0.05. -Demand pattern: Integer demands d i t ∈ {1, ...K} for each item have been randomly generated according to the following procedure:

1. We randomly select an item i * from a discrete uniform DU (1, N ) distribution and generate a value d i * T from a discrete uniform DU (1, K) distribution. 2. For each item i, except item i * , we randomly select a period t i from a discrete uniform DU (1, T ) distribution and generate a value d i ti from a discrete uniform DU (1, K) distribution. 3. For each entry in a N × T matrix, except for the entries corresponding to the (i, t) combinations for which we set d i t > 0 in steps 1 or 2, we randomly generate a number α it from a discrete uniform DU (1, N T ) distribution. 4. While the total cumulated demand (

N i=1 d i 1T
) does not exceed ρKT , we consider the entries (i, t) one by one in the increasing order of the corresponding value α it and generate a value d i t from a discrete uniform DU (1, K) distribution. 5. When the total cumulated demand reaches ρKT , we examine whether the corresponding instance is feasible by checking that

N i=1 d i 1t ≤ Kt for all t.
If the instance is infeasible, we repeat steps 1 to 4.

For each possible combination of problem dimension and production capacity utilization, 5 instances were generated, resulting in a total of 4 × 5 × 5 = 100 instances.

Computational experiments

All tests were run on a Pentium 4 (2.8 Ghz) with 505 Mb of RAM, running under Windows XP. We used a standard MILP software (CPLEX 11.1) with the solver default settings to solve the problems with one of the following formulations: -AGG: Aggregate formulation ( 9)-( 15). -AGG1: Aggregate formulation ( 9)-( 15) strengthened by the valid inequalities proposed in [START_REF] Vanderbeck | Valid inequalities for the Lasdon-Terjung production model[END_REF] for each item. These valid inequalities are added to the formulation at the root node of the Branch & Bound tree. The corresponding cutting-plane generation algorithm is based on the heuristic separation routine described by the authors of [START_REF] Vanderbeck | Valid inequalities for the Lasdon-Terjung production model[END_REF].

-AGG2: Aggregate formulation ( 9)-(15) strengthened by the valid inequalities (39) proposed in section 3 for each item. These valid inequalities are added to the formulation at the root node of the Branch & Bound tree using a standard cutting-plane generation procedure. -AGG12: Aggregate formulation ( 9)-(15) first strengthened by the valid inequalities (39) proposed in section 3 and then by the valid inequalities proposed in [START_REF] Vanderbeck | Valid inequalities for the Lasdon-Terjung production model[END_REF].

-EXT: Extended formulation proposed by Eppen and Martin in [START_REF] Eppen | Solving multi-item capacitated lot-sizing problems using variable redefinition[END_REF] for each item. We use the barrier algorithm embedded in CPLEX to solve the initial linear relaxation of the problem. -APP0.5: Approximate extended formulation (31)-(38) with ∆ i = 0.5d i,1,T . We use the barrier algorithm embedded in CPLEX to solve the initial linear relaxation of the problem.

Tables 1 to 4 display the computational results. For each set of 25 instances and each formulation, we provide: -Var. and Const.: the number of variables and constraints.

-#VI1 (resp. #VI2 ): the average number of valid inequalities from the family proposed in [START_REF] Vanderbeck | Valid inequalities for the Lasdon-Terjung production model[END_REF] (resp. proposed in section 3 of the present paper) added to the formulation by the cutting-plane generation procedure.

-#Feas (resp. #Opt): the number of instances out of the corresponding 25 instances for which a feasible solution (resp. a guaranteed optimal solution) could be found within 30 minutes of computation.

-Gap 0 : the integrality gap, i.e. the relative difference between the lower bound provided by the linear relaxation of the problem and the value of an optimal solution. For formulations AGG1, AGG2 and AGG12, we consider the lower bound obtained after the cutting-plane generation procedure has stopped. -#N odes: the average number of nodes explored by the Branch & Bound procedure before a guaranteed optimal solution is found or the computation time limit is reached.

-CP U IP : the average computation time in seconds required to find a guaranteed optimal solution. If one could not be found, we use the computation time limit of 1800 seconds.

-Gap: the average relative gap value obtained after 30 minutes of computation between the best integer solution and the best lower bound found.

Comparison of formulations AGG1, AGG2 and AGG12

We first discuss the results obtained with formulations AGG1 and AGG2 (see Tables 1234). These results show that the proposed valid inequalities (39) are more efficient at strengthening the aggregate formulation than the ones previously proposed in [START_REF] Vanderbeck | Valid inequalities for the Lasdon-Terjung production model[END_REF]. Namely, as can be seen e.g. for set A instances: -the reduction of the initial integrality gap (Gap 0 ) is greater with formulation AGG2 than with formulation AGG1 (2% vs 13%).

-the number of cutting-planes generated at the root node of the Branch & Bound tree is lower with formulation AGG2 than with formulation AGG1 (1951 vs 8463).

-the average computation time is significantly lower with formulation AGG2 than with formulation AGG1 (17s vs 594s). Moreover, results from Tables 1 to 3 also suggest that, in general, combining both families of valid inequalities (formulation AGG12) does not lead to shorter computation times. A noticeable exception can be found for the case where a large number of resources is involved (see Table 4). Specifically for set D instances, the average computation time is reduced from 711s to 636s by using AGG12.

Comparison of formulations AGG2, EXT and APP0.5

We now compare the strengthened aggregate formulation AGG2 with the extended formulation EXT. Results from Table 1 show that formulation EXT is more efficient than formulation AGG2 at solving small instances. Namely, the average computation time is reduced from 17s with formulation AGG2 to 6s with formulation EXT for set A instances.

However, formulation AGG2 is more efficient at solving large instances (see Tables 234). Indeed, the number of instances for which a feasible solution could be found within the computation time limit is higher with formulation AGG2 than with formulation EXT (75 vs 56 instances). This can be explained by the fact that, even if formulation EXT provides lower bounds of good quality, its size significantly increases with the problem dimensions, especially with the number of periods and/or resources. As a consequence, solving its linear relaxation is computationally demanding and only a limited number of nodes can be explored before the time limit is reached.

Moreover, the use of the approximate formulation APP0.5 does not improve the results obtained with formulation EXT with respect to the number of instances for which a feasible solution could be found within the computation time limit. This may be explained by the fact that aggregating part of the nodes and arcs in the corresponding network leads to a reduction in the size of the formulation but also to a decrease in the quality of the lower bounds provided by the linear relaxation of the problem. This is due to the fact that equalities Z t = α,γ (γ -β) + u(t, α, β, γ) used in the extended formulation proposed by [START_REF] Eppen | Solving multi-item capacitated lot-sizing problems using variable redefinition[END_REF] to evaluate startup costs have to be replaced by the weaker inequalities (35)-(36) in the approximate extended formulation.

Mathematical programming-based heuristic

Finally, we devised a simple heuristic as a first attempt to combine the advantages of formulations EXT and AGG2, namely the quality of the lower bound provided by EXT and the quality of the feasible solutions found by AGG2. More precisely, we implemented the following procedure:

1. We solve the linear relaxation of the problem using extended formulation EXT or approximate extended formulation APP0.5. 2. We look for the variables Y it with an integer non-negative value in the optimal continuous solution and fix them to this value.

3. We solve the resulting problem using formulation AGG2.

I = 10, T = 150, K = 2
Table 5 show the computational results obtained with this heuristic for the instances of sets B and D. We display:

-Formulation: the formulation (EXT or APP0.5) used to solve the initial linear relaxation.

-NVfix : the average number of variables Y it fixed at step 2 of the heuristic.

-#Feas: the number of instances out of the corresponding 25 instances for which a feasible solution could be found within the computation time limit of 600s.

-CP U tot : the average total computation time in seconds before the heuristic stops or before the computation time limit is reached.

-Gap heur : the average relative gap between the lower bound provided by the linear relaxation of formulation EXT or APP0.5 and the best integer solution found.

Results from table 5 show that we are able to obtain with the proposed heuristic better integer solutions within a shorter computation time for set B instances. Namely, the average gap after 10 minutes of computation (Gap heur =1.1%) is smaller than the gap obtained with formulation AGG2 after 30 minutes of computation (Gap=2.2%). However, no improvement was obtained for set D instances.

Finally we observe that the inequalities and reformulations presented above can be used to treat problems with sequence-dependent changeover costs in which constraints [START_REF] Van Eijl | On the discrete lot-sizing and scheduling problem with Wagner-Whitin costs[END_REF]

Figure 1 :

 1 Figure 1: An Integer K-flow (with K=3)

  set A (small instances): I = 10, T = 50, K = 2 -set B (instances with a large number of periods): I = 10, T = 150, K = 2 -set C (instances with a large number of items): I = 25, T = 50, K = 2 -set D (instances with a large number of machines): I = 10, T = 50, K = 10 -Costs: For each item, inventory holding costs (resp. startup costs) have been randomly generated from a discrete uniform DU (5, 10) distribution (resp. from a discrete uniform DU (100, 200) distribution).

  step. From flow conservation and the definition of t 1 , s t-1 + Y t,t1 ≥ 1. So either s t-1 ≥ 1, or Y t,t1 ≥ 1. However Y t,t1 ≥ 1 if and only if Y t + Z t+1,t1 ≥ 1, and the claim follows.

Inductive step. We assume that α(t, p) ≥ p is valid. There are two cases. Case 1. Y t+p + Z t+p+1,tp+1 ≥ 1. Adding this inequality to the inequality α(t, p) ≥ p gives α(t, p + 1) ≥ p + 1. Case 2. Y t+p + Z t+p+1,tp+1 = 0. This implies that Y u = 0 for u = t + p, • • • , t p+1 , which in turn implies that s t-1

Table 1 :

 1 Results for set A instances: I = 10, T = 50, K = 2

		Var.	Const. VI1	VI2 Feas/Opt Gap 0 (%) Nodes CPU IP (s) Gap(%)
	AGG	1550	2050			25/0	45	4.5e6	1800	14
	AGG1	1550	2050	8463		25/22	13	6981	594	0.4
	AGG2	1550	2050		1951	25/25	2	25	17	0
	AGG12	1550	2050	4130 1951	25/25	1.6	118	49	0
	EXT	16449	7666			25/25	0.3	1	6	0
	APP0.5 13096	6574			25/25	1	5	15	0
		Var.	Const. VI1	VI2	Feas/Opt Gap 0 (%) Nodes CPU IP (s) Gap(%)
	AGG1 *	4650	6150
	AGG2	4650	6150		17546	25/5	3	1287	1541	2.2
	AGG12 *	4650	6150
	EXT	155790 59014			16/10	1.5	20	1491	0.5
	APP0.5	129812 50501			16/10	2	10	1353	0.6

* Formulations AGG1 and AGG12 could not be solved due to exceeded memory capacity limitations.

Table 2 :

 2 for set B instances:

Table 3 :

 3 -(15) are replaced by (22)-(26).

		Var.	Const.	VI1	VI2 Feas/Opt Gap 0 (%) Nodes CPU IP (s) Gap(%)
	AGG1	3850	5200	15159		25/14	7	3360	817	0.8
	AGG2	3850	5200		2070	25/25	0.7	128	17	0
	AGG12	3850	5200	6432 2070	25/25	0.5	159	77	0
	EXT	16777 11233			25/25	0.3	8	17	0
	APP0.5 12822	9327			25/25	7	37	28	0

Results for set C instances: I = 25, T = 50, K = 2

Table 4 :

 4 Results for set D instances: I = 10, T = 50, K = 10 Formulation NVfix Feas CPU tot (s) Gap heur (%)

	set B	EXT	159	25	298	1.1
	set D	APP0.5	235	14	687	0.5

Table 5 :

 5 Results obtained with the heuristic for sets B and D instances[2] M. Constantino, A polyhedral approach to production planning models: startup costs and times and lower bounds on production. PhD thesis (1995), Université catholique de Louvain, Belgique.