Christophe Rapine

Bernard Penz

Céline Gicquel

Ayse Akbalik

Christrophe Rapine

Capacity acquisition for the single-item lot sizing problem under energy constraints

Keywords: Energy-efficient manufacturing, lot sizing, capacity acquisition, energy constraints, polynomial time algorithm

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Due to the rarefaction of natural resources, the hardening of environmental public legislation and the growing customer awareness for green products, industrial companies are increasingly under pressure to reduce the environmental impact of their industrial operations. In particular, the availability and affordability of energy is becoming a critical parameter affecting the whole life cycle of an industrial product, including its production phase [START_REF] Salonitis | Energy efficient manufacturing from machine tools to manufacturing systems[END_REF]. This translates into a growing attention of both researchers and practitioners on energy-efficient manufacturing.

Energy efficiency in manufacturing can be achieved by designing more ecoefficient products and by carrying out technological improvements of the production processes, but this usually requires large long-term financial investments. In contrast, energy-efficient production planning, by adjusting managerial parameters of the existing production processes, can be seen as a relatively inexpensive organizational measure whose benefit can be visible in the short-term (see [START_REF] Biel | Systematic literature review of decision support models for energy-efficient production planning[END_REF] and [START_REF] Gahm | Energy-efficient scheduling in manufacturing companies: review and research framework[END_REF]). Energy-efficient production planning is defined by [START_REF] Biel | Systematic literature review of decision support models for energy-efficient production planning[END_REF] as the computation of production plans which takes into account not only economic objectives (such as production costs or completion time minimization) but also energy-related objectives (e.g. energy cost minimization) or energy-related constraints (e.g. compliance with a maximum contracted consumption). It is particularly relevant for energy-intensive production processes as can be seen from the numerous case studies reported in the literature such as tile curing [ Özdamar and Birbil, 1999], melting shops and hot rolling mills for steel production [START_REF] Castro | Resource-task network formulations for industrial demand side management of a steel plant[END_REF], Tang et al., 2012], pulp and paper mills [Santos andAlmada-Lobo, 2012, Waldemarsson et al., 2013] or thermal cracking for ethylene production [Zhao et al., 2016]. Note that energy efficiency is also becoming a crucial element in the short-term planning of server farms activities [START_REF] Chen | Managing server energy and operationnal costs in hosting centers[END_REF], Mitrani, 2013].

Biel and Glock [START_REF] Biel | Systematic literature review of decision support models for energy-efficient production planning[END_REF] provide a recent literature review on decision support models for energy-efficient production planning. They show that most papers focus either on long-term aggregate planning or on short-term scheduling problems but that energy-efficient dynamic lot-sizing problems have been scarcely studied. Moreover, most papers make use of energy-related objectives aiming at minimizing either the energy cost or the energy consumption but implicitly assume that the energy is available in an unlimited quantity. However, due to the development of renewable energy sources and the difficulty to store electric energy on a large scale, energy providers have to increasingly rely on measures such as demand side management programs to shift demand for energy to off-peak periods and balance the electric grid [START_REF] Paulus | The potential of demand-side management in energy-intensive industries for electricity markets in germany[END_REF], Schultz et al., 2015, Shoreh et al., 2016]. For the manufacturers, this means in practice that electric energy cannot be considered anymore as a limitless resource and that they should seek not only to minimize their energy consumption but also to synchronize it with the available energy supply. Note that production scheduling under power availability restrictions has been investigated among others in [START_REF] Artigues | The energy scheduling problem: Industrial case-study and constraint propagation techniques[END_REF], Bruzzone et al., 2012, Mitra et al., 2012]. In this work, we consider production planning at a more aggregate decision level. We namely study a dynamic lot-sizing problem and introduce an energy constraint setting a limit on the total quantity of energy consumed in each time period rather than on the maximum instantaneous power demand.

More precisely, we investigate the single-item multi-resource Continuous Setup Lot sizing Problem (CSLP) with energy-related constraints. The CSLP is a small-bucket lot-sizing problem introduced by Karmakar and Schrage [START_REF] Karmarkar | The deterministic dynamic product cycling problem[END_REF]]. In the CSLP, a production lot on a machine consists of a set of consecutive planning periods in which the item is produced, and start-up costs are incurred at the first period of a production lot to take into account the costs related to startup operations. Note that the term "Continuous Setup" refers here to the fact that the quantity produced on a machine in a given period can take any value between 0 and the machine capacity. This contrasts with another small-bucket in Section 7.

Literature review

The classical lot-sizing problem aims at deciding, for a finite horizon subdivided into discrete periods, how much to produce in each period in order to satisfy a deterministic time-varying demand, with the objective to reach the best possible trade-off between production costs and inventory holding costs. We refer the reader to the seminal paper of Wagner and Whitin [Wagner and Whitin, 1958] and to the survey of Brahimi et al. citeBrahimi for a general introduction to this field. This work is mostly related to a small-bucket variant of the lotsizing problem known as the Continuous Setup Lot-sizing Problem (CSLP). Industrial applications of the CSLP to plan production on a set of unrelated parallel machines can be found in a variety of environments such as tile production [de Matta and Guignard, 1994], injection molding for healthcare products [START_REF] Dastidar | Scheduling injection molding operations with multiple resource constraints and sequence dependent setup times and costs[END_REF], packaging lines for dairy products [START_REF] Marinelli | Capacitated lot sizing and scheduling with parallel machines and shared buffers: A case study in a packaging company[END_REF] or glass melting furnaces [START_REF] Almada-Lobo | Multiple machine continuous setup lotsizing with sequence-dependent setups[END_REF]. In case of identical parallel machines, symmetry issues can arise during the solution phase due to the existence of many alternative optimal solutions obtained by renumbering the machines (see e.g. [Jans, 2009]). In order to avoid that, in small bucket models, it is useful to use an aggregate representation of the machines and to introduce integer rather than binary start-up variables. This has been exploited in a variety of real-world cases such as [Jans andDegraeve, 2004, Lasdon andTerjung, 1971] for tire manufacturing, [Sawik, 2009] for electronic devices production and [Kaczmarczyk, 2011] for mobile phone production. Note that the modeling of all the industrial applications mentioned in this paragraph does not include fixed setup costs to be incurred to keep a machine producing a given item once it has been started up for it, but only start-up costs to be incurred each time the setup state of the machine is changed from one item to the other.

As can be seen e.g. in the literature review provided by Jans and Degraeve [START_REF] Jans | Modeling industrial lot sizing problems: a review[END_REF], the research on lot-sizing has mainly focused on extending the basic problem in order to better model operational aspects or to integrate it as a substructure into larger tactical planning problems. However, it seems that only a limited effort has been done towards taking into account the growing pressure to establish sustainability in manufacturing and to comply with stricter environmental standards. In particular, as shown by the literature review provided by [START_REF] Biel | Systematic literature review of decision support models for energy-efficient production planning[END_REF], research on integrating energy efficiency in the lot-sizing problem appears to be only in its early stages.

A first situation in which energy should be explicitly taken into account is found when the energy needed by the production system is not only obtained from an external energy provider but is also produced internally, for instance by converting by-products of the process into steam or electrical energy. In this case, one has to simultaneously plan the activities of the pro-duction system with the activities of the energy system in order to coordinate them. Such problems are considered among others by Santos and Almada-Lobo [START_REF] Santos | Integrated pulp and paper mill planning and scheduling[END_REF], Waldemarsson et al. [Waldemarsson et al., 2013] and Zhao et al. [Zhao et al., 2016]. However, most often, the energy needed for production is obtained entirely from an external provider. In this case, a possible way to integrate energy in the lot-sizing problem consists in taking explicitly into account its cost in the objective function: see e.g. [START_REF] Özdamar | A hierarchical planning system for energy intensive production environments[END_REF], Uzel [Uzel, 2004], [START_REF] Heck | Lot-Size Planning with Non-linear Cost Functions Supporting Environmental Sustainability[END_REF] and Tang et al. [Tang et al., 2012]. Note that these approaches all assume that the energy is available in an unlimited quantity. We found a single approach integrating energy-related constraints in a lot-sizing problem. Masmoudi [START_REF] Masmoudi | Lot-sizing in a multi-stage flow line production system with energy consideration[END_REF]] consider a single-item multi-stage multi-resource lotsizing problem and introduce constraints setting a strict limit on the maximum instantaneous power demand used in each period. They propose two heuristic solution approaches for their problem and compare their performance to the one of the mathematical programming solver CPLEX. In contrast, we focus on providing structural results of the optimal solutions and on developing exact polynomial time algorithms.

A simplified version of the energy-LSP problem is very similar to the singleitem lot sizing problem with capacity acquisition and subcontracting studied by [START_REF] Atamtürk | Capacity acquisition, subcontracting, and lot sizing[END_REF]. The authors propose an exact algorithm running in O(T 3 ) time for the case with nonspeculative linear production costs, nonspeculative linear subcontracting costs and concave capacity acquisition costs. We show (see Section 4 and Appendix) that the same problem can be solved with a reduced time complexity of O(T log T ) in case of infinite subcontracting costs. Note that Li and Meissner [START_REF] Li | Capacitated dynamic lot sizing with capacity acquisition[END_REF] also study an integrated lot sizing and capacity acquisition problem in which subcontracting is not allowed but, in contrast to our work, they assume timeindependent fixed-charge production costs and focus on developing a heuristic polynomial time algorithm rather than an exact polynomial time algorithm. In the literature one can find many papers dealing with capacity planning in supply chains, but most of them focus on strategic decisions such as integrated facility location and capacity planning. For more details on the mathematical models proposed in capacity planning in manufacturing, see the review from Martìnez-Costa et al. [Martínez-Costa et al., 2014]. To the best of our knowledge, in the literature, such a capacity acquisition problem has not been yet studied under energy constraints.

Problem formulation

The energy-LSP consists in planning the production of a single-item, assuming a deterministic demand (d t in each period t) over a finite horizon of T periods. The production takes place on parallel and identical machines, each one with a maximum capacity U . The energy available in period t, E t , is used both for the production of units and for the start-up of the machines. We denote p t the amount of energy needed to produce one unit of product in period t and w t the amount of energy needed to start-up one machine in period t. The cost to start m machines in a period is f (m). The start-up cost functions f () are assumed to be concave and stationary. The concavity of f () allows us to capture in our model economies of scale that can be encountered in industry. We also assume that f (m) can be evaluated in constant time O(1) for any value of m. The production and inventory holding costs are linear: we denote c t the unit production cost and h t the unit holding cost for storing one unit between periods t and t + 1. Notice that we have no classical setup cost (paid in each period with a positive production amount), but only start-up costs to be incurred when one or several machines are switched on in compliance with small bucket multiresource CSLP we consider. We summarize below the different parameters of the problem:

c t unit
x t ≤ U m t ∀t ∈ {1..T } (2) p t x t + w t m + t ≤ E t ∀t ∈ {1..T } (3) m + t ≥ m t -m t-1 ∀t ∈ {1..T } (4) s t ≥ 0, x t ≥ 0, m t ∈ Z + , m + t ∈ Z + ∈ Z + ∀t ∈ {1..T } (5)
The objective function aims at minimizing the overall cost of machine startup, production and inventory. Constraints (1) are the balancing constraints linking inventory, ordering quantities and demand over T periods. Constraints (2) ensure that, in each period, the production quantity stays below the total production capacity which depends on the number of machines on. Constraints (3) guarantee that the energy consumed by both production and starting-up activities is limited by the available energy in each period. Constraints (4) link the number of machines started-up at the beginning of period t to the number of machines on in periods t and t -1, thus ensuring that the start-up costs f (m + t ) are correctly incurred in the objective function. The remaining constraints (5) define the feasibility domain of each decision variable.

In the following of the paper we assume that production costs follow Wagner-Whitin cost structure, also called non-speculative motives. It implies that for any period t of the time horizon, we have c t + h t ≥ c t+1 . Notice that, without loss of generality, we can consider that the holding costs are null, by substituting unit production costs c t with ct = c t + T u=t h u , see for instance [START_REF] Pochet | Production Planning by Mixed Integer Programming[END_REF]. The non-speculative production cost assumption implies that ct ≥ ct+1 holds for all periods, that is, modified unit production costs are non-increasing over time. Finally, recall that we focus on the case of stationary and concave start-up cost function f (). In the next section, we relax the energy constraints (3) of the problem to focus on what we call the pure capacity acquisition problem. We develop a very efficient algorithm for this case in Section 4, and we extend it in Section 5 to two cases with only one non-null energy parameter at a time (either start-up or production consumption).

Pure capacity acquisition problem

We begin with the study of a simplified version of the energy-LSP in which all energy consumption parameters are assumed to be null, that is, w t = p t = 0 for all periods t, and propose an algorithm in O(T 2 ) for this special case. This algorithm will be generalized to deal with positive energy parameters in Section 5. In Appendix, we also show that it can be implemented to run in time complexity O(T log T ).

Preliminaries

When all energy consumption parameters are null, energy-LSP reduces to a pure lot-sizing problem with capacity acquisition. Namely, the machines on do not incur any cost in a given period if they do not produce at that period. Thus, it is dominant to start-up all the necessary machines at the beginning of the planning horizon and to keep them on till the end. As a consequence, one simply has to decide the optimal value of m + 1 and set m t = m + 1 , ∀t, m + t = 0, ∀t = 2...T . The main decision thus consists in deciding how many machines to start-up at the beginning of the planning horizon, i.e. in determining the capacity to be acquired at the beginning of the planning horizon.

The energy-LSP problem with w = p = 0 is a special case of the problem investigated by [START_REF] Atamtürk | Capacity acquisition, subcontracting, and lot sizing[END_REF] in which there is no subcontracting option. Atamtürk and Hochbaum [START_REF] Atamtürk | Capacity acquisition, subcontracting, and lot sizing[END_REF] consider a continuous capacity model in which capacity can be acquired at any non-negative level. They show that, under the assumptions of non-speculative linear production and subcontracting costs and of concave capacity acquisition costs, the problem is polynomially solvable in time O(T 3 ). In this section, we show that there exists an efficient algorithm in O(T 2 ) to solve the special case with no subcontracting option under the same assumptions. The time complexity of our algorithm can even be improved into an O(T log T ) algorithm, using more sophisticated data structures (see Appendix).

We thus consider in what follows energy-LSP with w = p = 0. The problem is a capacitated lot-sizing problem: in each period, the amount produced cannot exceed the production capacity of the system. However, the production capacity C = U m + 1 is not a given parameter but is itself a decision variable of the problem. It means that we have to simultaneously determine the capacity C to be acquired (i.e. the number of machines m + 1 to start at the beginning of the planning horizon) and the production plan which should be feasible with respect to this installed capacity C. For the sake of simplicity, we consider in the description of the proposed algorithm a continuous capacity C and concave acquisition cost f (C). As explained at the end of Subsection 4.2, our results can be applied without any difficulty to the discrete capacity problem in which the value of C would be restricted to the set of integer multiples of U . Note that the problem may become infeasible if the production capacity C chosen is too small. We denote by C min the smallest capacity value for which the problem is feasible. Clearly, there exists a feasible planning for any given capacity C ≥ C min . Though the value of C min can be easily determined, the algorithm we propose does not need to compute it explicitely.

We now provide some insights of the ideas underlying the proposed algorithm.

• Since there are no setup costs and production costs are assumed to be nonspeculative, one should seek to plan production using a lot-for-lot strategy in order to reduce the inventory holding costs. However, as the production level is limited by the installed capacity C (which makes part of the decision), the demands exceeding C should be anticipated and smoothed to the previous periods. The optimal production plans have a particular structure: in all the periods of a subplan, except the first one, production occurs at full capacity, that is, x t = C. For a given value of the installed capacity C, determining the optimal production plan can be done in linear time O(T ), see [START_REF] Atamtürk | Capacity acquisition, subcontracting, and lot sizing[END_REF].

• Searching for the optimal solution can thus be reduced to the search of the optimal value of the installed capacity C. In their approach, Atamtürk and Hochbaum studied the polyhedral structure of the problem, and established that the number of extreme points can be bounded by O(T 2 ). Since each extreme point can be evaluated in linear time, it leads to an O(T 3 ) algorithm. In our approach, we show that among the O(T 2 ) extreme points, only a small number needs to be considered to find an optimal solution. Specifically, we focus on the function OP T (C) representing the cost of an optimal solution for an installed capacity C, and show that OP T (C) is a piecewise concave function comprising at most T breakpoints. Hence, the number of values of C to be considered for finding the minimum value of OP T (C) and the corresponding optimal solution is bounded by T , leading to an O(T 2 ) time algorithm.

Structure of an optimal solution

Clearly, a solution is entirely specified by its vector x = (x 1 , . . . , x T ) of production. For the ease of the presentation, we also use the associated vector s of inventory levels, where

s t ≡ (x 1 + . . . + x t ) -(d 1 + . . . + d t )
represents the amount of products in stock at the end of period t. We focus on solutions with a particular structure. For a given capacity C, we call a period t a full production period if x t = C. A period t with a null entering stock (s t-1 = 0) is called in the literature a regeneration period. We introduce the following definition :

Definition 1 We say that a solution x is no-slack if and only if s t-1 (C-x t ) = 0 holds for each period t = 1, . . . , T . That is, each period of the planning horizon is either a full production period or a regeneration period.

In addition, we say that solution x is dominant if it is feasible and s T = 0, that is, its final inventory is null.

Such a no-slack solution x defines in particular a subset B = {t|x t = C} of periods, corresponding to the periods where the capacity constraint is saturated. We will establish that conversely, each subset B defines a no-slack solution, and that only a small number of subsets need to be considered to find an optimal solution.

For a set B of periods, let us define A = {t | t / ∈ B and (t+1) ∈ B}. Set A is thus the set of periods immediately preceding a sequence of saturated periods. If a period t belongs to A , by definition period (t + 1) belongs to set B. Hence we can define v(t) as the largest index such that t + 1, • • • , v(t) belongs to B. We call set {t, • • • , v(t)} the block associated with a period t ∈ A . We denote by D i,j ≡ j t=i d t the demands to satisfy over the periods i upto j. For a given set B of periods, we define the associated no-slack solution x(B, C) as follows:

x t (B, C) =    C if t ∈ B D t,v(t) -(v(t) -t)C if t ∈ A d t otherwise
We illustrate the definition of a no-slack solution on the following example. We consider a planning horizon of 6 periods with demand d = (10,12,6,20,8,14) to satisfy. The no-slack solution x(B, C) associated with capacity C = 13.5 and set B = {4, 6} is represented in Figure 1a. The full production periods (set B) are represented in red and periods in set A are represented in grey. The network flow representation of this solution is given in Figure 1b, where the saturated arcs are represented by bold red arrows. In accordance with the definition of x(B, C), there is no entering stock in periods {1, 2, 3, 5}. It is easy to see that this solution is feasible and has a null ending stock level : It is also a dominant solution. Notice that solution x(B, C) may be infeasible, since the production in a period t / ∈ B may exceed capacity C. It may also be feasible but not dominant if it produces too much, for instance if C is a large value and set B = {1, . . . , T }. However, if x(B, C) is dominant, that is, is feasible and produces exactly the demand, we have the following property:

Property 1 If solution x(B, C
) is a dominant no-slack solution, then it is an optimal solution for the problem with capacity C.

Proof. Recall that ct represents the modified unit production cost after the addition of the linear holding cost. Since unit production costs ct are nonincreasing, a solution producing "at latest" is clearly optimal. For short, let us denote simply by x the solution x(B, C) and let X t = x 1 + • • • + x t-1 denote its cumulative production up to period t -1. It is sufficient to prove that for any feasible solution z, the cumulative production Z t is greater than or equal to X t for all periods. For a period t that does not belong to set B, this is immediate, since s t-1 = 0 implies that X t = D 1,t-1 . For a period t ∈ B, let v be the last period of its block. By construction, if v < T , we have s v = 0 since (v + 1) does not belong to set B. Notice that s v = 0 also holds if v = T , since we assume that x is dominant. Hence, we have

X t + (v -t + 1)C = D 1,v . Flow conservation imposes for any feasible solution that D 1,v ≤ Z t +z t +• • •+z v ≤ Z t +(v -t+1)C. The result follows.
Recall that we denote by OPT(C) the cost of an optimal solution relatively to an installed capacity C. For a set B of periods, Property1 shows that OPT(C) is realized by solution x(B, C) if it is dominant. Next property will allow us to establish that OPT(C) is a piecewise concave function:

OPT(C ) -OPT(C) = f (C ) -f (C) + t∈A ct (x t -x t ) + t∈B ct (C -C) = f (C ) -f (C) + t∈A ct (v(t) -t)(C -C ) + v(t) u=t+1 cu (C -C) = f (C ) -f (C) -(C -C) t∈A v(t) u=t+1 (c t -cu )
The term t∈A v(t) u=t+1 (c t -cu ) is fixed for a given set B, which means that the modified production costs vary linearly with C over [a, b]. Writing the expression OPT(C ) -OPT(C) for the particular value C = a, we obtain :

OP T (C ) = OP T (a) -f (a) + f (C ) -(C -a) t∈A v(t) u=t+1 (c t -cu )
The quantity OP T (a) -f (a) being fixed, OP T is the sum of a concave function and of a linear function and thus is itself a concave function. The result follows.

Next section describes how the capacity breakpoints of OP T (C) can be found very efficiently.

Description of the algorithm

Our algorithm constructs a series of capacities k) such that the no-slack solution x(B (i) , C (i) ) is dominant for each capacity value in [C (i+1) , C (i) ]. The algorithm stops when the first period is a full production period, that is, 1 ∈ B (k) . We claim that the capacity C (k) is then equal to the smallest capacity value, C min , such that the problem is feasible. Indeed, a necessary and sufficient condition for the problem with a capacity C to be feasible is that D 1,t ≤ tC for all periods t = 1, • • • , T . If a subset B contains the first period and solution x(B, C) is feasible, then we have D 1,v = vC for v = v(1). This implies that C = C min . As a consequence, any capacity C admitting a feasible plan belongs to (at least) one of the intervals [C (i+1) , C (i) ]. We will demonstrate that at most T different capacities C (i) need to be considered, that is, k is bounded by T .

C (1) ≥ C (2) ≥ • • • ≥ C (k) and a series of subsets B (1) ⊂ B (2) ⊂ • • • ⊂ B (
Let us define capacity C (1) = max t {d t } and B (1) = {t|d t = C (1) }. Clearly the no-slack solution associated with B (1) is dominant, and thus realizes OPT(C (1) ). Consider now that we have determined a capacity C (i) and a subset B (i) such that x(B (i) , C (i) ) is dominant, and that 1 / ∈ B (i) (otherwise the algorithm stops). With immediate notation, we denote by A (i) the set of A periods associated with B (i) . We want to determine the smallest value C ≤ C (i) for which x(B (i) , C) is dominant. Consider a capacity C ≤ C (i) . Using its definition, solution x(B (i) , C) is dominant if and only if the two following conditions hold:

(i) ∀t / ∈ A (i) ∪ B (i) , C ≥ d t (ii) ∀t ∈ A (i) , C ≥ D t,v(t) /|{t, . . . , v(t)}|
We define

C (i+1) = max max t / ∈A (i) ∪B (i) d t , max t∈A (i) D t,v(t) |{t, . . . , v(t)}| (6)
Notice that by construction, x(B (i) , C (i+1) ) is dominant. However, in this solution, at least one period is both a full production period and a regeneration period. That is, there exists a full production period not belonging to B (i) . We define:

B (i+1) = {t | x(B (i) , C (i+1) ) = C (i+1) }
By construction we have i+1) . As a consequence, set B (i) contains at least i periods, and thus the algorithm stops after at most T steps. Since each step can be performed in time O(T ), determining the C (i) 's and the B (i) 's takes O(T 2 ) time.

B (i) ⊂ B (
To output the optimal solution, we use Property 2. Since OPT(C) is concave over each interval [C (i+1) , C (i) ], the optimal value is reached at one of the capacities {C (1) , • • • , C (k) }. Notice that determining the optimal cost OPT(C (i) ) for a given capacity can be performed in linear time O(T ). Hence, the algorithm runs in time complexity O(T 2 ). We can state the following theorem:

Theorem 1 An optimal solution OPT for the pure capacity acquisition problem can be found in time O(T 2 ).

In case of a discrete capacity problem, where C is restricted to be a multiple of U , we consider for each index i = 1, • • • , k the capacities C (i) /U U and C (i) /U U . It represents at most 2k different values. Again, due to Property 2, we can assert that the discrete capacity minimizing OPT(C) belongs to this set. Evaluating all these 2k values to find the minimal value can be performed in time complexity O(T 2 ).

In the next section we illustrate the idea of the algorithm on a simple example.

Illustrative example

We consider the same instance as in Figures 1a and1b : we have 6 periods with demand d = (10,12,6,20,8,14) to satisfy. In the Figures 2 to 6 below, we have represented the optimal policy related with the current capacity C. We used different colors to distinguish periods belonging to A or B: A red rectangle corresponds to a period in B (saturated periods), a grey rectangle corresponds to a period in A . Other periods are represented by white rectangles. Notice that for these periods the production is equal to the demand. Finally, the red circles indicate the periods that are saturated but not yet in the current set B. Note that those periods will enter set B in the next step.

First step (i = 1)

The algorithm starts with an empty set B and finds the first capacity C (1) to consider, which corresponds to the highest demand in period 4, here C (1) = max t {d t } = 20. By construction, period 4 is saturated (red circle in Figure 2) and it enters set B. Thus, we have B (1) = {4}. Clearly x(B (1) , C (1) ) is dominant (hence optimal for C (1) ). Second step (i = 2)

We search now the lowest capacity C (2) ≤ C (1) for which solution x(B (1) , C (2) ) is dominant. We compute C (2) using Equation (6). Note that A (1) = {3}. We obtain:

C (2) = max max t / ∈{3}∪{4} d t , D 3,4 4 -3 + 1 = max 14, 26 2 = 14
This value corresponds to the demand of period 6, hence, the new capacity C (2) is obtained by hitting the second highest demand. See Figure 3. The new

B (2) is thus B (2) = B (1) ∪ {6} = {4, 6}.
Third step (i = 3) Again, using Equation ( 6), we find C (3) = max 12, max 26 2 , 22 2 = 13. This time, a period of set A (2) becomes saturated, here, period 3 (red circle in Figure 4). Hence,

B (3) = B (2) ∪ {3} = {3, 4, 6}.
Fourth step (i = 4) C (4) = max 12, max 38 3 , 22 2 = 12.66. Period 2 in set A (3) becomes saturated (red circle in Figure 5). Thus, B (4) = B (3) ∪ {2} = {2, 3, 4, 6}. In the next section, we generalize this algorithm to energy-LSP when only one of the energy consumption parameter is null. We also show that it can be implemented to run in time complexity O(T log T ).

Polynomial time algorithms for energy-LSP

In this section we extend the algorithm presented for the pure capacity acquisition problem and give two polynomial time algorithms, for the cases when one of the consumption energy parameters of our model, either start-ups or unit production, is time-varying and the other is negligible. In both cases, we show that the problem is equivalent to a capacity acquisition problem with additional production restrictions. That is, in addition to the nominal capacity C to be acquired, in each period t we also have a restriction on the production quantity, due to the limited energy availability. We will call this limitation F t . Hence, in each period t, the production is limited by the minimum between C and F t . We show that with null setup costs and under non-speculative motive, we can find the nominal capacity C to install and the optimal production plan in time O(T 2 ), using the same approach as the one presented in the previous section for the pure capacity acquisition problem.

5.1 Case with w = 0, p t > 0, E t > 0: Production energy consumption

We consider that the unit energy consumption p t and the available energy E t may vary from one period to another and that the energy needed to start a machine is negligible, that is, w t = 0 for all periods t. Hence, whatever the installed capacity C, i.e. whatever the number of machines started up in period 1, the production quantity x t in period t is limited by E t /p t . We thus have to decide on a nominal capacity C for the system, that is on the number of machines to start at period 1, and on a production plan where the production quantities comply with the restrictions imposed by both the installed capacity C and the energy availability E t /p t . In case the acquisition cost function f is concave in C, the resulting problem can be seen as a capacity acquisition problem with additional production restrictions F t ≡ Et pt .

5.2 Case with p = 0, w t > 0, E t > 0: Start-up energy consumption

We now assume that the start-up of a machine is the operation requiring the greatest amount of energy (think for instance of furnaces), and that unit production energy consumption is negligible compared to it. More precisely, we consider that the start-up energy consumption w t and the available amount of energy E t are time-varying, and that the unit production energy consumption is null for all periods t. This problem seems to significantly differ from problems described in Sections 4 and 5.1 since it may not be possible to install the nominal capacity C, i.e. to start all the machines needed, in the first period of the planning horizon due to the energy constraints. However, under the assumption of a linear acquisition cost f , that is, if f (C) = f C, the problem can be seen as a special case of the problem described in 5.1. For the ease of the discussion, we focus on the discrete capacity case, that is, when we have to choose the number m + t of machines to start in each period. Since the production energy consumption is null and the acquisition cost f is linear and stationary, it is dominant to start as many machines as possible in each period, till the desired nominal capacity of the system is obtained. In a period t, we can start at most E t /w t machines. We define:

N t = t u=1 E u /w u
Clearly, in any solution, at most N t machines can be available for production in period t. Hence x t ≤ U N t is a dominant inequality for the problem, that is, production in period t is limited by F t = U N t . Again, we have a capacity acquisition problem with additional production restrictions F t .

Extension of the pure capacity acquisition algorithm

We consider the capacity acquisition problem with additional production capacities F t . These additional capacities are imposed by the energy parameters, as explained in 5.1 and 5.2. As a consequence, in each period t, the production is limited by C t ≡ min{C, F t }. Notice that this production capacity is timevarying. We show how the algorithm described in Section 4 can be adapted to deal with these additional constraints. As in the previous section, we assume null setup costs, non-speculative motive, and a stationary and concave acquisition cost function f (). Recall that for start-up energy consumption parameters, the acquisition function is restricted to be linear, see Section 5.2.

Again, we start by considering the case of a continuous nominal capacity: C can be any real positive value. To deal with the time-dependent capacities C t = min{C, F t }, we extend the definition of the no-slack solution x(B, C) given in the previous section as follows:

x t (B, C) =    C t if t ∈ B D t,v(t) - v(t) u=t+1 C u if t ∈ A d t otherwise.
It is easy to check that Property 1 still holds with this definition for our case of energy-LSP. To extend Property 2 to Property 3 below, we need the notion of F -interval. For a given nominal capacity C, we define F (C) = {t|F t < C} as the set of periods whose capacity is limited by F t . Let us denote by l) the l ≤ T distinct values appearing in set F . We also add the dummy value F (0) = 0. We call an F -interval a real interval of the form [F (j) , F (j+1) ]. Notice that set F (C) does not change inside an F -interval. We have the following property:

F (1) ≤ • • • ≤ F (
OPT(C ) -OPT(C) = f (C ) -f (C) + t∈A ct (x t -x t ) + t∈B ct (C t -C t ) = f (C ) -f (C) + t∈A ct v(t) u=t+1 (C u -C u ) + v(t) u=t+1 cu (C u -C u ) = f (C ) -f (C) + t∈A v(t) u=t+1 (c t -cu )(C u -C u ) = f (C ) -f (C) + t∈A v(t) u=t+1 (c t -cu )(min{C, F u } -min{C , F u })
Now, since [a, b] is included in an F -interval, we have the same set of periods where the capacity is limited by F u for C and C . More precisely, C u ≡ min{C, F u } is equal to F u if period u belongs to F (b), and to C otherwise. Hence, we can write that:

OPT(C ) -OPT(C) = f (C ) -f (C) -(C -C) t∈A u∈{t,...,v(t)}\F (b) (c t -cu )
Hence, OP T is the sum of a concave function and of a term which varies linearly with C. It is thus a concave function with respect to the capacities over the interval [a, b].

In the same way as in Section 4, our algorithm constructs a series of capacities

C (1) > C (2) > • • • > C (k) and a series of subsets B (1) ⊆ B (2) ⊆ • • • ⊆ B (k)
such that the no-slack solution associated with B (i) is dominant for any capacity in [C (i+1) , C (i) ]. However, contrary to the algorithm for the pure acquisition problem, we can not assert that B (i) is a strict subset of B (i+1) . But we will prove the invariant that, either B (i) ⊂ B (i+1) , or C (i+1) does not belong to the same F -interval as C (i) . It will allow us to establish that k ≤ 2T , that is, at most 2T different capacities C (i) need to be considered.

To avoid straightforward cases, we assume that a feasible solution exists if C ≥ max t {F t }. We can compute an initial no-slack optimal solution x (0) for capacities (F t ) t=1,•••,T in time O(T ) using a simple backward algorithm to shift the demands. If S t denotes the amount of demand which has to be shifted backward (i.e. whose production should be anticipated) due to the limited production capacity, then for all period t = 1, . . . , T , we have x (0) t = min(S t+1 + d t ; F t ) and S t = S t+1 + d t -x t , initializing the backward computation with S T +1 = 0. Let us define capacity C (1) = max t {x (0) t } and B (1) = {t|x (0) t = min{F t , C (1) }}. Hence x (0) is the no-slack solution associated with set B (1 ), and by construction it is optimal for capacity C = C (1) and its value is equal to OPT(C (1) ).

Consider now that we have determined a capacity C (i) and a subset B (i) such that x(B (i) , C (i) ) is dominant, and that 1 / ∈ B (i) . Using the same approach as in Section 4, we want to determine the smallest value C ≤ C (i) for which x(B (i) , C) is dominant. For a capacity C ≤ C (i) , solution x(B (i) , C) is dominant if and only if the two following conditions hold:

(i) ∀t / ∈ A (i) ∪ B (i) min{C, F t } ≥ d t (ii) ∀t ∈ A (i) min{C, F t } ≥ D t,v(t) - v(t) u=t+1 min{C, F u }
To determine C, we strengthen our requirement by imposing that C must be in the same F -interval as C (i) . Let us define F (i) = max t {F t |F t < C (i) }. By convention, we let F (i) = 0 if all the additional capacities F t are greater than or equal to C (i) . Clearly, interval [ F (i) , C (i) ] is included in an F -interval. Inside this interval, solution x(B (i) , C) remains dominant if and only if:

(C1) ∀t / ∈ A (i) ∪ B (i) C ≥ d t (C2) ∀t ∈ A (i) v(t) u=t min{C, F u } ≥ D t,v(t) (C3) C ≥ F (i)
Namely, note that, as x(B (i) , C (i) ) is a feasible solution, we have: i) and t ∈ F (C (i) ). By construction, C is restricted to stay within the same F -interval as C (i) so that F (C) = F (C (i) ). We thus have d t ≤ F t ≤ C for all periods t s.t. t / ∈ A (i) ∪ B (i) and t ∈ F (C). Hence, whatever the value of C ≥ F (i) , condition (i) holds for the periods t / ∈ A (i) ∪ B (i) s.t. min(C, F t ) = F t . It can thus be simplified into condition (C1).

d t ≤ F t ≤ C (i) for all periods t s.t. t / ∈ A (i) ∪ B (
Moreover, observe that min{C,

F u } is equal to F u if u ∈ F (C (i)
) and to C otherwise. Hence condition (C2) can be rewritten for a period t ∈ A (i) as: i) for all the periods u of a block {t, . . . , v(t)}, certainly solution x(B (i) , C) remains dominant for any capacity C ≥ F (i) . Also, the right hand side of the previous inequality must be negative in this case. With the convention that the ratio of a negative number divided by 0 is equal to -∞, we define:

u∈{t,...,v(t)}\F (C (i) ) C ≥ D t,v(t) - u∈{t,...,v(t)}∩F (C (i) ) F u Notice that if F u < C (
C (i+1) = max max t / ∈A ∪B d t , max t∈A D t,v(t) -u∈{t,...,v(t)}:Ft<C (i) F u |{u ∈ {t, . . . , v(t)} : C (i) ≤ F t }| , max t {F t | F t < C (i) } (7) By construction, x(B (i) , C (i+1)
) is dominant. We define:

B (i+1) = {t | x t (B i , C (i+1) ) = min(F t , C (i+1) )}
That is, B (i+1) is constituted of the periods producing at full capacity in solution x(B (i) , C (i+1) ). By construction we have B i ⊆ B i+1 . We have two possible situations, depending on which inequality (C1), (C2) or (C3) is saturated by

C (i+1) : • Either C (i+1) = F (i) . It implies that F (i+1) < F (i) .
• Otherwise, in solution x(B (i) , C (i+1) ), at least one period is both a full production period and a regeneration period. That is, there exists a full production period not belonging to B (i) . It implies that

|B (i+1) | > |B (i) |.
Since a set B can contain at most T periods and there are at most T distinct values F t , the algorithm stops after at most 2T steps. Since each step can be performed in time O(T ), determining the C (i) 's and the B (i) 's takes O(T 2 ) time.

To output the optimal solution, we use Property 3. By construction, each capacity interval [C (i+1) , C (i) ] is included in an F -interval, and x(B (i) , C) is dominant for any capacity C ∈ [C (i+1) , C (i) ]. As a consequence, for a continuous nominal capacity, the optimal value is reached at one of the capacities {C (1) , • • • , C (k) }. Notice that determining the optimal cost OPT(C (i) ) for a given capacity can be performed in linear time O(T ). Hence, the algorithm runs in time complexity O(T 2 ). In Appendix, we show that this algorithm can be implemented to run in time complexity O(T log T ) using efficient data structures. We have the following theorem:

Theorem 2 An optimal solution can be found in time O(T log T ) for energy-LSP with either null start-up energy consumption and concave start-up costs f () or null unit production energy consumption and linear start-up costs f (). The algorithm requires at most O(T ) evaluations of function f ().

In case of a discrete capacity problem, where C is restricted to be a multiple of U , we consider for each index i = 1, • • • , k the capacities C (i) /U U and C (i) /U U . It represents at most 4T different values. Again, due to Property 3, we can assert that the discrete capacity minimizing OPT(C) belongs to this set. Evaluating all these 2k values to find the minimal value can be performed in time complexity O(T ) as all the breakpoints of OPT(C) are already computed. Hence, the complexity of the algorithm is not modified in the discrete capacity case.

Extension to more general start-up cost functions

In this section, we show that the algorithms presented in this paper can be easily adapted, by slightly increasing their time complexity, to deal with the case where the start-up (capacity acquisition) cost function f is convex, or more generally, is a piecewise function, with a mix of convex and concave pieces. We restrict this discussion to the case of discrete capacities, that is, the admissible values of C belong to set {0, U, 2U, . . . , M U }, where M is an input and represents the maximal number of available machines. Several real applications can motivate the use of convex start-up cost function. For instance, when the existing capacity in a shop floor is not enough to cover the excess demand, it may be necessary to call for external or remote resources. Those resources are typically more costly due to the supplementary costs of transportation, material handling, and the spot market costs. Another possible application corresponds to capacity reservation contracts, where a certain capacity is reserved a priori at the supplier to obtain a more advantageous price. Once this level is exceeded, the cost to acquire additional capacities become more expensive.

Consider first the case where the start-up cost function f is a convex function of C. In the previous sections, f was assumed to be concave, and, as a consequence, the optimal cost OPT(C) was a piecewise concave function of the installed capacity C. Looking at the proof of Property 3, we know that on each interval [C (i+1) , C (i) ], the optimum cost OPT(C) is of the form f (C)+(αC +β), that is, we add an affine term to function f . Since the sum of two convex functions is a convex function, we have now that OPT(C) is a piecewise convex function, whose breakpoints are the capacities C (i) determined by our algorithm. The only difference is that we can not assert that the optimum of OPT(C) is realized at one of its breakpoints. However, since a convex function is in particular unimodal, we can use a Ternary Search algorithm (see [?]) to find the minimum value of OPT(C) on each interval [C (i+1) , C (i) ]. Ternary Search algorithm is a divide & conquer algorithm, similar to a Binary Search for monotonic functions, except that the search domain is divided into three pieces, and at each step the domain is reduced by at least one third. Note that the running time of a Ternary Search algorithm when restricted to integer values is in O(log l+1) on an interval containing l integers. Hence, OPT(C) can be determined on each interval [C (i+1) , C (i) ] in O(log M ) operations. As a consequence, the overall time complexity of the algorithm is in O(T (log T + log M )), where the first term O(T log T ) corresponds to the determination of the breakpoints C (i) , and the second term O(T log M ) corresponds to finding the optimum value OPT. Hence the time complexity is only slightly increased. However, notice that the number of evaluations of function f is now in O(T log M ), while the number of evaluations is linear for concave functions, see Theorem 2. Now, let us consider a more general case where f is a piecewise function of C, where each piece is either concave or convex. Let R be the number of pieces and {b 0 , b 1 , ..., b R } be the set of corresponding breakpoints. We assume in the following that the list b 1 ≤ ... ≤ b R of its breakpoints are given in the inputs in the description of f . We also assume that on each interval, f can be evaluated, possibly through an oracle, in constant time O(1). As already discussed in Section 5, OP T (C) is the sum of the start-up costs given by the piecewise function f plus the production and inventory costs, which are a piecewise linear function of the installed capacity. This can be exploited in order to find the optimal solution in polynomial time thanks to the following algorithm:

• Step 1. Use the algorithm described in Section 5 to identify the set of capacities C (1) , . . . , C (k) where the block structure of the production plan changes, that is, the set of breakpoints of the production/inventory holding cost function, together with the values OPT(C (1) ), . . ., OPT(C (k) ). This step can be performed in time O(T log T ), and is independent of the startup cost function f .

•

Step 2. Merge all the breakpoints {C (1) , ...C (k) } and {b 1 , ..., b R } to obtain an ordered list a 0 , ..., a r , ..., a l . Since we assume that both the C (i) 's and the b j 's are indexed in non-decreasing order, this can be done in linear time O(T + R) (i) , C (i+1) ] for some index i, we can evaluate OPT(C) in constant time using the expression :

OPT(C) = f (C)+OP T (C (i) )-f (C (i) )+ OPT(C (i+1) ) -OP T (C (i) ) C (i+1) -C (i) (C-C (i) )
Hence the optimum on each interval can be determined in time complexity O(1) if f is concave on the interval, and in time complexity at most O(log M ) if f is convex, using a Ternary Search algorithm.

• Step 4. Output the best value found in the course of the algorithm.

The overall time complexity of the algorithm is dominated by Steps 1 and 3 and is in O(T log T + l log M ), with l the number of breakpoints to consider, which is bounded by (T + R), and by M , since intervals that do not contain an integer can be ignored. Hence, if

M ≥ T + R, the time complexity is in O((T + R) log M ). Conversely, if M ≤ T + R, the time complexity is in O(T log T + R log M ).
Notice that if f is piecewise concave, that is, without convex pieces, then the complexity of the algorithm is reduced to O(T log T +R).

Conclusion and perspectives

In this paper we studied a single-item lot sizing problem integrated with energy consumption constraints. We showed that, without energy parameters, the associated problem becomes a pure capacity acquisition problem, as stated in Atamtürk and Hochbaum [START_REF] Atamtürk | Capacity acquisition, subcontracting, and lot sizing[END_REF], but without the subcontracting option. We developed an original solution approach and proposed an O(T log T ) time algorithm to solve this problem. The same algorithm can be easily extended to some special cases of energy-LSP, with only one non null energy parameter, without altering its time complexity. This algorithm allows to tackle with quite general acquisition cost functions as well, such as concave, convex, or piecewise functions with a mix of concave and convex pieces, by slightly modifying its computational complexity.

We note that the same integrated energy & lot sizing problem can also be seen in the context of a carbon cap to not exceed in each period, and the start-up & production activities that emit carbon units. Another possible application is a certain limited manpower to not exceed in each period, together with a capacity limitation of parallel resources. Notice that the common property of those systems is the existence of a non-storable and periodical capacitated resources.

Among the possible research perspectives of this work, one can investigate whether our approach can be successfully extended to the capacity acquisition problem with subcontracting option or to more general cost structure, keeping a low resolution time complexity. Regarding the energy-LSP problem, one open question is the existence and the design of a polynomial time algorithm for the case with both positive start-up and production energy consumption. Another extension of the model would be to consider that a started machine consumes energy not only while actively producing but also while being idle.

[ Shoreh et al., 2016]

Appendix: Improvement of the time complexity

We show in this appendix that the algorithm described in Section 4 and 5 can be implemented to run in time complexity O(T log T ), as stated in Theorem 2.

We focus on how determining OPT(C). Recall that OPT(C) is a (continuous) piecewise concave function with k ≤ 2T breakpoints, C (1) ≥ C (2) ≥ . . . ≥ C (k) . As a consequence, the optimum is realized at one of the breakpoints, and thus it suffices to evaluate OPT(C (i) ) for i = 1, . . . , k, where k is the number of step of the algorithm. We will establish that the breakpoints C (i) 's can be determined in overall complexity O(T log T ) using efficient data structures. Notice that evaluating an optimal solution for a given capacity C can be performed in linear time. However, as we may have up to 2T values to evaluate in the course of the algorithm, it is still a too large computational effort. Instead, we compute the optimal cost OPT(C (i+1) ) incrementally from the value of OPT(C (i) ). As established in Property 3, the optimum value OPT(C) is concave on each interval [C (i+1) , C (i) ]. More precisely, let us denote by α (i) the quantity:

α (i) = t∈A (i) u∈{t,...,v(t)}:C (i) ≤Fu (c t -cu )
Then, following the proof of Property 3, we have:

OPT(C (i+1) ) = OPT(C (i) ) + f (C (i+1) ) -f (C (i) ) + α (i) (C (i+1) -C (i) )
As a consequence, value OPT(C (i+1) ) can be computed in constant time, if OPT(C (i) ) and the slope α (i) are known.

We now show how these different values can be updated at each step of the algorithm. For a given set B, recall that a block is a set of consecutive periods {t, . . . , v} such that t ∈ A and v = v(t). This corresponds to a subplan of the associated no-slack solution. We maintain in an array all the blocks of the current set B. Initially, there is one block for each single period. In addition, each period refers to the block where it belongs to, simply by using an index. For a current capacity C and a current set B, each block B contains the following data:

• Its first and last periods, t(B) and v(B), respectively.

• D(B), the total demand of the periods of the block, that is, t∈B d t .

• F (B), defined as u∈B:Fu<C F u if set {u ∈ B : F u < C} is not empty, and 0 otherwise.

• Its pseudo-cost c(B), defined as u∈B:C≤Fu cu if set {u ∈ B : C ≤ F u } is not empty, and 0 otherwise.

• Its pseudo-size n(B), defined as |{u ∈ B : C ≤ F u }|, that is, the number of periods of the block with an additional capacity F t being greater than or equal to the nominal capacity C. By convention, |∅| = 0.

We need these different quantities in order to be able to update efficiently the slope α (i) of the current no-slack solution. Consider a block B and let t = t(B) be its first period. If period t ∈ A (i) , then B = {t, . . . , v(t)} and we have:

u∈{t,...,v(t)}:C (i) ≤Fu (c t -cu ) = n(B)c t -c(B)
If period t does not belong to set A (i) , then B = {t}, that is, t is a lot for lot period where demand d t is produced. In this case, we clearly have n(B)c tc(B) = 0. Hence, we can rewrite the slope α (i) as:

α (i) = B:B a block of B (i) n(B)c t(B) -c(B)
Now, consider a step of the algorithm. Suppose that we have found the current breakpoint C (i) . Finding the next breakpoint consists in:

(i) Determining C (i+1) according to Equation (7) (ii) Updating the blocks of the new set B (i+1) (iii) Evaluating OPT(C (i+1) )

As step (iii) can be computed in constant time, we focus on the two first steps.

For the first step (i), we need to find the largest value of C saturating one of the constraints (C1), (C2) or (C3), see Equation ( 7). For constraints (C1) and (C3), we can sort via a precomputation phase all the F t and d t values in non-increasing order. This precomputation phase can be performed in time O(T log T ). During the execution of the algorithm, we maintain an index of the current position in each of the arrays. Clearly, at each step, we can find the next F t value lower than C (i) in constant time. Recall that each period stores the block it belongs to, and hence can be determined in constant time if it belongs to the current set A (i) or B (i) (by comparing its index with the values t(B) and v(B) of its block). From the current position in the sorted array of d t values, we scan the array until we encounter a period not in the current set A (i) or B (i) . Notice that it may require up to O(T ) operations. Nevertheless, since a period entering A ∪ B does not leave this set in the following steps of the algorithm, each period is scaned only once. Hence the overall time complexity is in O(T ). Finally, if C (i+1) saturates constraint (C2), we need to determine the largest value of: e t ≡ D t,v(t) -u∈{t,...,v(t)}:Ft<C (i) F u |{u ∈ {t, . . . , v(t)} : C (i) ≤ F t }| over all the periods t ∈ A (i) . Notice that e t is also equal, with our notations, to:

e t = D(B) -F (B) n(B)
with B = {t, . . . , v(t)} being the block beginning at period t. We note that for a given block B, value e t does not change if capacity C remains in the same F -interval. Hence, at each step, we need to retrieve the largest value e t and to update (delete and re-insert) at most one e t , corresponding to the block updated in the current set B. In order to maintain the e t values of the blocks in a data structure H, we can use for instance a heap or a balanced binary tree such that these operations can be performed in logarithmic time in the number of elements in the structure. Notice that the number of periods of A is at most T , and thus each operation can be performed in time complexity O(log T ). In conclusion, step (i) of the algorithm requires O(T log T ) operations in total.

For the second step (ii), we need to update the blocks of the current set B (i+1) . We distinguish 3 cases depending which constraint (C1), (C2) or (C3) has been saturated by C (i+1) : Case 1. Constraint (C1) is saturated. Let t be a period not in A (i) ∪B (i) such that i+1) . In this case period t belongs to the current solution B (i+1) . This corresponds to the merge of the block B = {t} with its preceding block B = {t(B ), . . . , t -1}. Block B can be found in constant time, since period (t -1) keeps track of its index. Merging B and B into a new block B = {t(B ), . . . , v(B)} can be done in constant time, since each attribute (F (B ), n(B ), . . . ) can be easily deduced from the attributes of blocks B and B . In order to update the new belonging of each period of B to this block, we chose to index B with the index of B , such that only period t needs to be updated.

d t = C (
Case 2. Constraint (C2) is saturated. Let t be a period in A (i) such that e t = C (i+1) . We are in a situation similar to the previous case when constraint (C1) was saturated: this corresponds to the merge of the block B = {t, . . . , v(t)} with its preceding block B = {t(B ), . . . , t -1}. We have the same operations to perform. However, updating the index of its block for each period of B may now be the most costly part, since both B and B may be large sets of periods. For that, we choose for the index of B the index of the largest block between B and B , such that we need to reindex only the periods of the smallest block. Of course, we may still have up to O(T ) periods to reindex. However, we claim that the overall complexity for the algorithm is bounded by O(T log T ). To prove that, let g(m) be the total number of reindexations incurred by the algorithm from the beginning to obtain a block of size m. That is, for each period of the block, we count for the number of times it has been reindexed by the algorithm. We prove by induction that g(m) ≤ m log m. For m = 1, this is true since no merging is performed for a block of a single period. Case 3. Constraint (C3) is saturated. Let t be a period such that C (i+1) = F t . In this case, the structure of the blocks is not modified, that is, B (i+1) = B (i) . However, the slope α (i+1) is modified. More precisely, let B be the block where period t belongs. We update the attributes of block B by setting F (B) ← F (B) -F t , c(B) ← c(B) + ct and n(B) = n(B) + 1. The slope α (i+1) is then equal to α (i+1) = α (i) + (c t(B) -ct ). In this case, we update the e value associated with block B by removing it from H and by inserting its new value.

The time complexity of the algorithm is thus dominated by the precomputation phase, the updating of the e values in the data structure H, and the reindexing of the periods to keep track of their blocks, each one requiring O(T log T ) operations. Hence the algorithm can be implemented in overall time complexity O(T log T ).

  Figure 1: (a) The no-slack solution x(B, C) associated with capacity C = 13.5 and set B = {4, 6} and (b) Its network flow representation.

Property 2

 2 Let B be a set of periods and [a, b] a real interval. If for all capacities C in [a, b] the solution x(B, C) is dominant, then OPT(C) is concave with respect to C over [a, b]. As a consequence, its minimum over [a, b] is reached at one of the extremities of the interval. Proof. Consider C and C two capacities of interval [a, b]. Let x = x(B, C) and x = x(B, C ) be the solutions associated to C and C , respectively. Notice that x and x differ only by their production in periods of set A ∪ B. Precisely, we have:

Figure 2 :

 2 Figure 2: First step (i=1): Search of C (1) and B (1)

Figure 3 :

 3 Figure 3: Second step (i=2): Search of C (2) and B (2)

Figure 4 :

 4 Figure 4: Third step (i=3): Search of C (3) and B (3)

Figure 5 :

 5 Figure 5: Fourth step (i=4): Search of C (4) and B (4)

Figure 6 :

 6 Figure 6: Fifth step (i=5): Search of C (5) and B(5)

Property 3

 3 Let [a, b] be an interval included in an F -interval. If there exists a set B of periods such that x(B, C) is dominant for all capacities C in [a, b], then OPT(C) is concave with C over [a, b]. As a consequence, the minimum of OPT(C) over [a, b] is reached at one of the extremities of the interval. Proof. Consider C and C two capacities of interval [a, b]. Let x = x(B, C) and x = x(B, C ) be the solutions associated to C and C , respectively. Notice that x and x differ only by their production in periods of sets A and B. Precisely we have:

  Consider now a block of m periods, requiring g(m) reindexation, obtained by merging two blocks A and B. Without loss of generality, assume that |A| ≤ |B|. We have:g(n) ≤ |A| + g(|A|) + g(|B|) ≤ |A| + |A| log |A| + |B| log |B| ≤ |A| + |A|(log m -1) + |B| log m ≤ (|A| + |B|) log m = m log mThe first inequality is due to the fact that the algorithm reindexes only the smallest block, A, when merging A and B. The induction hypothesis implies the second inequality. The third inequality simply uses that |A| ≤ |B| and hence |A| ≤ m/2. This allows us to establish that at most T log T reindexations are performed by the algorithm.In both cases, the values e eventually associated with blocks B and B (if they are not reduced to a single period) are removed from data structure H, and the new value e associated with B is inserted, which can be done in time O(log T ). Finally, we update the slope by setting α (i+1) = α (i) + (n(B )c t(B ) -c(B )) -(n(B )c t(B ) -c(B )) -(n(B)c t(B) -c(B)).

•

  Step 3. On each interval [a j , a j+1 ], OPT(C) is either convex or concave (depending if f itself is either convex or concave). Since each interval [a j , a j+1 ] is included by construction in an interval [C

  Shoreh, M. H., Siano, P., Shafie-khah, M., Loia, V., and

	Catal ão, J. P. (2016). A survey of industrial applications of demand re-
	sponse. Electric Power Systems Research, 141:31 -49.
	[Tang et al., 2012] Tang, L., Che, P., and Liu, J. (2012). A stochastic pro-
	duction planning problem with nonlinear cost. Computers & Operations Re-
	search, 39(9):1977 -1987.
	[Uzel, 2004] Uzel, E. (2004). A mathematical modeling approach to energy
	cost saving in a manufacturing plant. Master's thesis, Izmir Institute of
	Technology.
	[Wagner and Whitin, 1958] Wagner, H. M. and Whitin, T. M. (1958). Dynamic
	version of the economic lot size model. Management Science, 5(1):89-96.
	[Waldemarsson et al., 2013] Waldemarsson, M., Lidestam, H., and Rudberg, M.
	(2013). Including energy in supply chain planning at a pulp company. Applied
	Energy, 112:1056 -1065.

[Zhao et al., 2016] Zhao, H., Ierapetritou, M. G., and Rong, G. (2016)

. Production planning optimization of an ethylene plant considering process operation and energy utilization. Computers & Chemical Engineering, 87:1 -12.

Acknowledgments.

This project was supported by Conseil Régional de Lorraine (2015)(2016). We are grateful to the anonymous referees, whose comments helped to improve the presentation of this paper.


Mouna Kchaou Boujelben

Celine Gicquel

A stochastic model for the design of electric vehicle charging infrastructure

Keywords: Flow refueling location model, electric vehicle charging station network design, stochastic range

We study the problem of designing electric vehicle (EV) charging infrastructure. The main decision that we consider is the location of the charging stations in a way that EV drivers can drive along the road network without running out of charge. We take into account the uncertainty of the driving range, which is the maximum distance that a fully charged EV can travel before its battery runs empty. We thus propose a stochastic model and maximize the expected coverage of the recharging demand. We first formulate a new mixed-integer linear program (MILP) for the stochastic problem and compare it with a previously published one. We then develop a Tabu search heuristic method to solve large-size instances of the problem. We carry out our numerical experiments using randomly generated road networks and we show the performance of the new formulation as well as the quality of the solutions provided by the Tabu search heuristic approach.

Introduction

The problem of optimizing the design of an EV charging infrastructure is mostly related to a well-known problem in the Operations Research field, which is the facility location problem. Basically, facility location models aim at finding the best possible locations for facilities in order to optimize a performance criterion such as cost minimization or demand coverage maximization. Reviews on the academic literature on facility location can be found e.g. in [START_REF] Melo | Facility location and supply chain management -A review[END_REF]. In this paper, we study a facility location model for the design of an EV charging infrastructure. The aim of this model is to determine the best locations for EV charging stations in a way to maximize the coverage of the charging demand. Different approaches have been used in the literature to model this type of demand. The first approach assumes that the demand is located at a set of fixed points in space, i.e. is expressed at the nodes of the underlying network. An example of such models can be found in [START_REF] Lisbon | Optimal locations of charging stations for electric vehicles in a neighborhood[END_REF]. These models implicitly consider that drivers will carry out a special-purpose round trip from their home or workplace to the charging station in order to recharge the battery. However, in some cases, drivers do not specifically travel to the station to recharge their vehicles but rather recharge it while on their way to another destination. For EVs, this occurs in particular during longdistance trips exceeding the vehicle range. In this case, demand should be represented as a set of origin-destination paths rather than a set of nodes. Such a representation was first proposed by Hodgson [START_REF] Hodgson | A flow capturing location-allocation model[END_REF] within his flow capturing location model. The main issue in this model is that it does not take into account the limited range of EVs, i.e. the maximum distance that a fully charged EV can travel before its battery runs empty. This was the motivation behind the development of the flow refueling location model by Kuby and Lim in [START_REF] Kuby | The flow-refueling location problem for alternative-fuel vehicles[END_REF]. In this model, a flow on a given path is considered as refueled if and only if there is an adequate number of stations on the path, allowing an EV driver to refuel his journey from the origin to the destination. Different extensions of this model have been proposed in the literature in order to study more realistic assumptions such as e.g. the capacity of charging stations [START_REF] Upchurch | A Model for Location of Capacitated Alternative-Fuel Stations[END_REF] or the deviations from the shortest paths [START_REF]The deviation-flow refueling location model for optimizing a network of refueling stations[END_REF]. However, little attention has been given to studying the uncertainties on the input data of the optimization problem, in particular the driving range uncertainty. To the best of our knowledge, the uncertainty of the driving range has been introduced in only two recent works. The first work [START_REF] De Vries | Incorporating driving range variability in network design for refueling facilities[END_REF] uses a probabilistic approach to model the range uncertainty and the second work [START_REF] Capar | An arc-cover path-cover formulation and strategic analysis of alternative-fuel station locations[END_REF] introduces two stochastic models, one maximizing the expected EV flow covered and the other one is a chance constrained model. In this paper, we present a new formulation for the problem of locating EV charging stations under uncertainties on the driving range, with the objective of maximizing the expected flow coverage.

Problem description and modeling

We consider a road network G(N,A), where N denotes a set of nodes and A denotes a set of arcs linking these nodes. The charging demand to be satisfied by the stations is modeled as a set of flows denoted by Q. Each flow qϵQ is described by its origin O q , its destination D q and the number of EVs traveling along it f q . Drivers belonging to flow q are assumed to follow the shortest path between O q and D q . All EVs have a limited range R, which is subject to uncertainty due to e.g. the traffic conditions, the weather or the age of the battery. A flow q is said to be covered if there is an adequate number of stations on this flow that allows EVs to travel from O q to D q without running out of charge. In other words, the flow q is covered when the distance between each pair k,l of consecutive stations on the flow does not exceed the range R, otherwise, it is not covered. In order to model the problem that determines the best locations for a predetermined number of stations while maximizing the expected flow covered on the network, we use three types of decision variables. Variable x j is a binary variable equal to 1 if a station is opened at node j, 0 otherwise. Variable z q ∈ [0,1] is a continuous variable equal to the probability of covering flow q. Finally, variable w kl q is a binary variable equal to 1 if a vehicle traveling along flow q is recharged at a station located at node k to cover the vehicle trip up to node l and 0 otherwise. The objective function of our MILP model can be written as follows: maximize ∑ 𝑓 𝑧 ∈ , i.e. maximizing the sum of EV flows on the network, weighted by their probability of coverage. We include in the model different sets of constraints that link variables z q and variables w kl q , define the relationship between variables w and variables x and limit the number of charging stations that must be opened to a predetermined number p representing the limited investment budget. We assume that for a given realization of the random conditions w, the value of the range R(w) is the same for the whole network and that R(w) is randomly distributed following the cumulative density function 𝐺: 𝑅 → [0,1]. By using these assumptions and the fact that G is a nondecreasing function, we can write the range constraints as follows: 𝑧 ≤ 1 -∑ 𝐺 𝑡 𝑞 (𝑘, 𝑙) ∈ 𝑤 ∀𝑞 ∈ 𝑄, 𝑙 ∈ 𝑁 \𝑂 , where 𝑡 (𝑘, 𝑙) is the total length of arcs visited when traveling from node k to node l on flow q and 𝑁 is the set of nodes situated along flow q before node l on a trip from Oq to Dq.

Tabu search heuristic approach

For large size instances of the stochastic EV charging station location problem, the computation time required to get an optimal solution might become prohibitively long. Therefore, we propose a Tabu search procedure in order to obtain good quality solutions in short computation times. The algorithm starts by building an initial solution and setting the best feasible solution to the initial solution. Then, two different steps are alternated in the iterations of the algorithm. Step1 consists of selecting a station to be opened among the N-p closed stations. The station should not be Tabu (not recently closed) and should lead to the highest expected coverage among all possible openings. Step2 consists of selecting a station to be closed among the p + 1 opened stations. This station should not be Tabu (not recently opened) and should lead to the highest expected coverage among all possible closings. The Tabu search procedure stops when the number of iterations without improvement of the best objective reaches a maximum limit.

Main conclusions from the numerical experiments

In our numerical experiments, we compare the performance of the new formulation proposed in this paper with the one of the previously published formulation in [START_REF] De Vries | Incorporating driving range variability in network design for refueling facilities[END_REF] for the expected flow refueling location model. We also analyze the quality of the solutions provided by the Tabu search heuristic approach. The road networks that we use in the tests were randomly generated following a procedure similar to the one proposed in [START_REF] Capar | An arc-cover path-cover formulation and strategic analysis of alternative-fuel station locations[END_REF]. We consider two different instances sizes: (N=100, Q=1225) and (N=200, Q=4950) and different values for the number of stations p. The random range R is represented using a Gamma distribution, with a shape parameter of 50 and a scale parameter of 5. For each instance size and each value of p, we generate 5 instances and consider the average output value. We employed the C++ language to implement the model and the commercial solver CPLEX 12.6.2 to solve it.

All tests were carried out on a PC with Intel i5-3210M Core 2 Duo (2.50 GHz) with 8GB of RAM. The CPU time of CPLEX solver was limited to 10 hours. Examples of results are shown in Table 1. The results show that our formulation performs better than the existing one as it leads to significantly decreasing the CPU time. Moreover, when using the existing formulation, CPLEX failed at finding the optimal solution within the time limit for several instances. In those cases, the optimality gap went up to 13%. The last column of Table 1 reports the gap between the solution found by the Tabu Search heuristic and the optimal solution. It shows that the heuristic performs very well as it provides good quality solutions in short computation times.

Conclusion

In this paper, we studied the problem of EV charging infrastructure planning under range uncertainty. We proposed a new MILP formulation for the problem, explicitly including the driving range. Then we developed a Tabu search heuristic approach to solve large-size instances. Our numerical experiments showed the performance of the new formulation as well as the quality of the solutions provided by the heuristic method. However, in our model, we assumed that for any realization of the random conditions, the value of the driving range is the same for all flows. An interesting direction for further research would thus be to relax this assumption and to study the more realistic case where the driving range realization is different on each cycle segment of a flow.

Table 1 .

 1 Average performance of the three methods (5 replications).

	Instance	New formulation	Existing	TABU CPU	TABU Gap
		CPU (s)	formulation CPU (s)	(s)	(%)
	N100Q1225p10	165	583	9	0.2
	N100Q1225p25 609	11008	24	0.1
	N200Q4950p10 8212	32775	58	0.8
	N200Q4950p25 26385	36000	184	0.6