Ruiwen Liao
email: ruiwen.liao@universite-paris-saclay.fr

Céline Gicquel

Integrated lot-sizing and energy supply planning with onsite generation of intermittent renewable energy

Keywords: Energy-efficient manufacturing, production planning, lot-sizing and scheduling, sequence-dependent changeovers, renewable energy, decentralized energy system

This work considers an industrial production site partially powered by a decentralized energy system based on intermittent renewable energy sources. Our objective is to simultaneous plan the industrial production and the energy supply in this site so as to minimize the total cost. A new way of modelling this combinatorial optimization problem is proposed: it relies on the extension of a multi-product single-resource small-bucket lot-sizing model called the proportional lot-sizing and scheduling problem. This extension involves among others sequence-dependent changeover times overlapping multiple periods and energy-related constraints. Our numerical results show that the resulting mixed-integer linear programming model enables to obtain good-quality production and energy supply plans with a computational effort much smaller than the one required by a previously published large-bucket lot-sizing model.

Introduction

According to a recent report of the International Energy Agency [START_REF] David | Industry[END_REF], industry emitted in 2021 a total of 9.4Gt of CO2, accounting for a quarter of the global CO2 emissions. These emissions should be drastically reduced to a total of 7Gt of CO2 by 2030 in order to achieve carbon neutrality by mid-21st century as set by the COP-21 Paris Agreement. This translates into a huge pressure on industrial companies to lower the carbon footprint of their activities. At the same time, energy cost has become a primary concern for them due to the sharp increase in the price of gas and grid electricity. At a company level, a possible way of simultaneously addressing these two challenges consists in building a decentralized energy system based on renewable sources (e.g. wind, sun) and to use the renewable electricity generated on-site to power, at least partially, the industrial processes. This is already a reality for a growing number of companies. For instance, the British carmaker Bentley Motors equipped its Crewe, UK, plant with more than 30,000 solar panels representing a total energy capacity of 7.7 MW. At peak generation times, up to two-third of the electricity consumed by the factory can be provided by this solar energy system [START_REF] Motors | Largest UK solar car port installed at Bentley factory in Crewe[END_REF][START_REF] Szymkowski | Bentley will add 10000 more solar panels to its plant[END_REF]. Similarly, 30% of the electricity used in the Anheuser-Busch's Budweiser brewery factory, located in California, are supplied by an on-site wind turbine and 6,500 solar panels (Anheuser-Busch 2021).

However, using a decentralized energy system based on renewable sources to power an industrial process poses major difficulties. Namely, due to their intermittence, the amount of energy provided by renewable sources displays strong intraday and seasonal variations. This volatility means that on-site generated renewable electricity alone cannot power an industrial process and that a combination of renewable and grid electricity should be used. Moreover, energy providers usually rely on time-of-use pricing schemes to prompt consumption when electricity supply is high and lower it when electricity supply is low. For an industrial company, the electricity purchasing price thus also displays strong intraday variations. Consequently, the electricity cost of a production site depends on both the amount and the timing of electricity consumption, which themselves depend on production planning decisions. Thus, energy supply and industrial production should be planned simultaneously.

The present work focuses on the short-term industrial production and energy supply planning problem of a factory using on-site generation of renewable electricity. In terms of industrial production planning, a multi-product single-resource lot-sizing problem with sequence-dependent changeover costs and times is considered. Basically, lot-sizing consists in deciding when and how many finished products should be processed on a resource so as to satisfy the customers' demand while minimizing the total changeover and inventory holding costs [START_REF] Drexl | Lot sizing and scheduling -survey and extensions[END_REF]. Lot-sizing arises in production planning whenever the resource can process a single type of finished products at a time and needs to undergo changeover operations (e.g., tool changes) each time the type of finished products processed on it is changed. These changeover operations incur both a changeover cost and a changeover time, whose values depend on the production sequence. As for the inventory holding costs, they represent the cost of keeping finished products in inventory on the time interval between their production and their shipment to the customers. In terms of energy supply planning, three main elements are considered. First, on-site power generation devices produce a time-varying amount of electricity from renewable sources: this electricity is assumed to be free but the generated amount depends mainly on the weather conditions and cannot be controlled. Second, there is an on-site energy storage system with a limited capacity. Finally, this decentralized energy system is connected to the national electric grid and can trade (i.e., buy and sell) electricity with it at a time-varying price. The energy supply planning problem consists in deciding when and how much electricity should be traded with the grid and charged into/discharged from the storage system so as to make sure that energy supply and consumption are balanced at all time. Figure 1 represents the studied integrated energy supply and industrial production system.

Dynamic lot-sizing problems such as the one studied here use a finite planning horizon divided into a discrete set of periods. However, in the present case, defining this discrete set of periods is not straightforward. Namely, the time discretization needed to track the satisfaction of the customers' demand, the one needed to build the industrial production plan and the one needed to manage the balance between energy supply and demand may significantly differ. More precisely, our planning problem involves two exogenous time structures, i.e. two sets of points in time corresponding to externally given events defined by the data of the model [START_REF] Copil | Simultaneous lotsizing and scheduling problems: a classification and review of models[END_REF]. The first time structure is defined by the points in time at which the customers' demand arises. Using these points leads to the definition of rather long periods (typically days or weeks). The second time structure is imposed by the discrete time grid used to track the on-site generation of electricity and the variations in the prices for trading electricity with the energy provider. This time grid usually relies on short periods (typically hours or 15-minutes intervals). In addition to these two exogenous time structures, we may define a third endogenous time structure to plan production. This one corresponds to the points in time at which internal events related to production planning (such as the time at which the processing of a production lot or a changeover between two production lots starts) are captured by decision variables [START_REF] Copil | Simultaneous lotsizing and scheduling problems: a classification and review of models[END_REF].

This modeling difficulty was first tackled by Wichmann, Johannes, and Spengler (2019a) who investigated an energy-oriented lot-sizing problem to plan production in a factory powered by a decentralized power system comprising solar cells and Lithiumion batteries. They used a three-level time structure to build the integrated energy supply and industrial production plan: a fixed time grid based on long macro-periods to track the customers' demand satisfaction, a fixed time grid based on short energyoriented micro-periods to manage the energy supply and consumption, and a flexible time grid based on short production-oriented micro-periods to build the production plan. More precisely, their model can be seen as an Energy-Oriented General Lot-sizing and Scheduling Problem or EOGLSP, i.e., an extension of the General Lot-sizing and Scheduling Problem introduced by [START_REF] Fleischmann | The general lotsizing and scheduling problem[END_REF] to integrate energy supply planning and lot-sizing. However, due to the discrepancy between the energy-oriented and production-oriented time grids, a large number of binary variables and constraints need to be introduced in the mathematical formulation to precisely track in which energy-oriented micro-period the industrial activities planned on the production-oriented time grid actually consume energy. Consequently, the size of the resulting mixed-integer linear program (MILP) formulation increases sharply with the instance size, leading to numerical difficulties during the resolution by a mathematical programming solver.

To overcome these difficulties, we investigate here a new modeling approach and propose to rely on a time structure comprising only two discretization levels. The first time grid is based on long macro-periods to track the customers' demand satisfaction. The second one relies on the short micro-periods imposed by the energy-related input data and is used to both manage the energy supply and plan the industrial production. Planning production using a time structure involving a rather large number of short periods amounts to formulating a 'small bucket' lot-sizing problem, i.e. a lot-sizing problem in which the number of different product types that may be produced in each period is strictly limited. In this work, we propose to limit this number at two, in other words to allow at most one changeover to start in each micro-period. We thus formulate the problem as an Energy-Oriented Proportional Lot-sizing and Scheduling Problem or EOPLSP [START_REF] Drexl | Proportional lotsizing and scheduling[END_REF]. As shown by the results of our numerical experiments, the resulting MILP formulation is much smaller than the one obtained with the EOGLSP model and the computational effort needed to obtain optimal or good-quality feasible solutions for medium-size instances is significantly reduced.

Our contributions are thus threefold. First, we consider the same optimization problem as the one investigated by Wichmann, Johannes, and Spengler (2019a), i.e., we seek to build an integrated industrial production and energy supply plan for a factory equipped with a decentralized energy system using renewable sources. For this problem, we propose a new modeling approach based on a two-level time structure and the use of the fixed time grid imposed by the energy-related input data to plan the industrial production. This approach leads to a drastic reduction in the number of binary variables and constraints involved in the resulting MILP formulation. However, using short micro-periods to determine lot-sizes means that changeover times may exceed the duration of a micro-period and that changeover operations may overlap several periods. Although this problem was investigated by [START_REF] Kaczmarczyk | Explicit modeling of multi-period setup times in proportional lot-sizing and scheduling problem with variable capacity[END_REF] for the PLSP with multi-period sequence-independent changeover times, a direct extension of his work to our setting leads to the formulation of a large-size MILP and to numerical difficulties. Our second contribution thus consists in proposing a new extension of the model recently introduced by [START_REF] Kaczmarczyk | Explicit modeling of multi-period setup times in proportional lot-sizing and scheduling problem with variable capacity[END_REF] in order to accommodate sequence-dependent changeover times and to improve the numerical tractability of the resulting MILP. Finally, we carry out numerical experiments on a large set of small to medium-size instances to compare the proposed EOPLSP model with the EOGLSP one previously published by Wichmann, Johannes, and Spengler (2019a). Our numerical results show that in most cases, our model provides a production plan of the same quality as the one provided by the EOGLSP model but with a significantly reduced computational effort.

The paper is organized as follows. An overview of the literature on energy-efficient lot-sizing is given in Section 2. Section 3 provides a detailed description of the optimization problem under study. Section 4 presents the proposed EOPLSP model based on a two-level time structure. Section 5 discusses the results of the computational experiments conducted on small to medium-size instances. Conclusion and research perspectives are given in Section 6.

Related works

As mentioned by [START_REF] Biel | Systematic literature review of decision support models for energy-efficient production planning[END_REF], the scarcity of resources, rising energy prices, and an increasing awareness of the environmental impact of manufacturing recently triggered research on energy-efficient production planning. Basically, energy-efficient production planning aims at computing production plans considering not only the traditional production-related costs but also energy-related objectives and constraints. The reader is referred to [START_REF] Gahm | Energy-efficient scheduling in manufacturing companies: A review and research framework[END_REF][START_REF] Biel | Systematic literature review of decision support models for energy-efficient production planning[END_REF][START_REF] Gao | A review of energy-efficient scheduling in intelligent production systems[END_REF][START_REF] Bänsch | Energy-aware decision support models in production environments: A systematic literature review[END_REF][START_REF] Terbrack | Energy-oriented production planning in industry: a systematic literature review and classification scheme[END_REF] for comprehensive literature reviews on this field. In what follows, we discuss the works most closely related to dynamic lot-sizing.

A first set of works dealing with energy-aware dynamic lot-sizing consider that energy may only be bought from an external supplier and aim at reducing the corresponding cost. Energy costs should be considered explicitly in lot-sizing when the amount of energy consumption is a non-linear function of the produced quantity and/or when the energy buying price is time-dependent. [START_REF] Özdamar | A hierarchical planning system for energy intensive production environments[END_REF] thus sought to minimize the total inventory holding, production and energy costs in a tile curing industry. As the energy costs are proportional to the number of active curing kilns, the setup and lot-sizing decisions in the production plan have a direct impact on the energy consumption and cost. [START_REF] Heck | Lot-size planning with non-linear cost functions supporting environmental sustainability[END_REF] studied a single-product uncapacitated lot-sizing problem. They introduced in the objective function an additional term, expressed as a non-linear function of the lot size, to take into account the energy costs. Similarly, [START_REF] Tang | A stochastic production planning problem with nonlinear cost[END_REF] considered a single-level single-machine lot-sizing problem in the iron and steel industry. The energy consumed by the heat furnace is a non-linear function of the production quantity so that decisions on the lot sizes directly impact the energy cost. Note that, in all the above-mentioned papers, the energy prices are assumed to be time-invariant. [START_REF] Johannes | Energyoriented production planning with time-dependent energy prices[END_REF] and Wichmann, Johannes, and Spengler (2019b) studied a single-level single-machine lot-sizing problem with time-dependent energy prices. They showed that it is possible to significantly reduce the overall production and energy costs by setting the timing and level of production while taking into account the energy price fluctuations. [START_REF] Perraudat | Stochastic programming approaches for an energy-aware lot-sizing and sequencing problem with incentive[END_REF] recently investigated a stochastic multi-product single-machine lot-sizing problem faced by a company taking part to the demand response program of its energy provider. In such a program, the company receives monetary incentives to reduce its electric load during periods defined by the energy provider.

A second set of works dealing with energy-aware lot-sizing consider that the energy needed to power the industrial process is obtained from both an external supplier (through e.g. the national electric grid) and an on-site power generation system using renewable sources. [START_REF] Rodoplu | Single item lot sizing problem under renewable energy uncertainty[END_REF] addressed a multi-machine multi-level lot-sizing problem with a single end-product and a serial gozinto structure. In their model, the planning horizon is split into a discrete set of rather long planning periods (6 to 8 hours) defined by the points in time at which demand for the endproduct occurs. Within each of these periods, the electricity price is assumed to be constant, the energy supply-demand balance is considered in an aggregate way and the detailed production plan on each resource is built using a continuous-time flow-shop scheduling model. Note that the use of such large time buckets may not be accurate enough to track the variations in the energy price and availability. In contrast, [START_REF] Golpîra | Robust smart energy efficient production planning for a general job-shop manufacturing system under combined demand and supply uncertainty in the presence of grid-connected microgrid[END_REF]) studied lot-sizing for a multi-machine multi-level system powered by a grid-connected micro-grid comprising among others a wind turbine, a combined heat and power generator and an electrical storage system. They used the same time structure based on short one-hour planning periods to track the demand satisfaction and build the industrial and energy production plan. They assumed that the energy consumption depends only on the resource setup status and not on the production quantity and developed a robust optimization approach to deal with the uncertainties of the wind power generation and energy consumption. Finally, Wichmann, Johannes, and Spengler (2019a) studied a single-level single-machine lotsizing problem under intermittent renewable energy and energy storage. The authors proposed an extension of the GLSP to handle this problem [START_REF] Fleischmann | The general lotsizing and scheduling problem[END_REF][START_REF] Meyr | Simultaneous lotsizing and scheduling by combining local search with dual reoptimization[END_REF]. In their model, the planning horizon is discretized following three different time grids. The planning horizon is first divided into macro-periods to track demand fulfillment. Each macro-period is then further subdivided into productionoriented micro-periods of variable length and energy-oriented micro-periods of fixed length. Production-related micro-periods are used to define the production plan as would be done in the GLSP: at most one type of product may be produced per microperiod and the length of each micro-period is not fixed in advance but may be flexibly chosen when building the production plan. Energy-oriented micro-periods are used to track energy generation, trading and consumption. However, the discrepancy between the production and energy time structures makes it possible for a production-oriented micro-period to span over several energy-oriented micro-periods and for an energyoriented micro-period to include several production-oriented micro-periods. Consequently, additional variables and constraints need to be introduced to track in which energy-oriented micro-period the consumption of energy by the production process takes place. They determine, e.g., whether a production-oriented micro-period overlaps a given energy-oriented micro-period, the total duration of this overlapping and the type of activities (manufacturing, changeover, preserving) carried out during the overlapping time. The number of these variables and constraints is thus of order (number of production-oriented micro-periods × number of energy-oriented micro-periods).

This results in the formulation of an MILP which involves a large number of binary variables and big-M type constraints and is consequently computationally intensive to solve with a mathematical programming solver. In what follows, we introduce a new modeling approach for this problem which is intended to be more computationally efficient.

Problem description

We consider an industrial plant producing a set of J finished products to satisfy an external demand and focus on building an industrial production and energy supply plan for this plant over a finite horizon.

The points in time at which the external demand occurs define a first coarse discretization of the planning horizon into a set T = {1, ..., T } of macro-periods. Let l be the fixed length of a macro-period. There is a demand d jt for each product j to be satisfied at the end of macro-period t. The demand cannot be backlogged and should always be met on time. In case some inventory of finished products is kept for later use at the end of macro-period t, an inventory holding cost h jt has to be paid per unit of product j kept in inventory. Let I j,0 be the inventory level of product j at the beginning of macro-period 1.

There is a single production resource. This resource may be in four distinct states: idle, changeover, manufacturing, and preserving. The idle state, represented by a dummy product j = 0, corresponds to the case where the machine is turned on but not ready for production. Let J = {0, 1, 2, . . . , J} denote the set of all products. In order to produce product j, the machine should be in the right configuration, i.e. should be setup for this product. The machine is in the changeover state whenever actions to change its configuration from a product j ′ ∈ J to another product j ∈ J \ {j ′ } are on-going. Re-configuring the machine from j ′ to j incurs a sequence-dependent changeover cost f c j ′ ,j and a sequence-dependent changeover time st j ′ ,j . A changeover cannot be interrupted, i.e., once a changeover process is launched, the reconfiguration actions need to be executed until completely finished. Once the machine is setup for a given product j, it can produce it: this corresponds to the manufacturing state. The production time required to produce one unit of product j is denoted by k j . Finally, we assume that the setup status of the machine for a given product may be conserved even if it is not producing it: this corresponds to the preserving state. The main advantage of using setup status conservation when technically possible is that it allows to restart the manufacturing of the product immediately after an non-production period whereas the other alternative would consist in turning off the machine to the idle state and carrying out a costly and time-consuming changeover from the idle state to the setup status to produce the product. Depending on the production process type, this setup status conservation may or not incur significant preserving costs. Let pc j denote the cost of preserving the machine setup status for product j during one unit of time. We assume w.l.o.g. that the machine is idle at the beginning of the planning horizon.

The machine consumes energy whenever it is in the changeover, manufacturing or preserving state, but not when it is idle. The amount of energy consumed per unit of time during a changeover from product j ′ to product j is given by e F j ′ ,j : the total energy consumed by a changeover from j ′ to j is thus equal to e F j ′ ,j st j ′ ,j . The energy consumed when manufacturing product j is proportional to the produced quantity: let e M j denote the energy needed to produce one unit of product j. Finally, the energy consumption by the machine in the preserving state is proportional to the time during which its setup status is preserved. Let e P j be the energy consumed per unit of time to keep the machine setup for product j without manufacturing it. Note that e M 0,r = e P 0,r = 0 as the machine does not consume energy when in the idle state.

Regarding the energy supply, two sources are available: the on-site power generation devices producing electricity from renewable sources (photovoltaic panels or wind turbines) or the main electric grid. Both the amount of locally generated renewable energy and the prices to trade electricity with the main grid are time-varying. In order to track these variations, we introduce a second exogenously defined time grid. This one relies on a fine discretization of the planning horizon through a set of microperiods sufficiently short to assume that the renewable energy generation and the trading prices with the grid are constant within a micro-period. Thus, each macroperiod t ∈ T is further divided in R micro-periods of length l e such that l = Rl e . Let R = {1, ..., RT } be the set of all micro-periods and R t = {(t -1)R + 1, , ..., tR} be the set of micro-periods belonging to macro-period t ∈ T . The amount of renewable electricity available in micro-period r ∈ R is represented by p R r . This on-site generated electricity is assumed to be free. The unit price to buy (resp. to sell) electricity from (resp. to) the main grid in r ∈ R is denoted by π GP r (resp. π GS r ). Finally, the electricity traded with the main grid must go through an electrical transformer, which induces energy losses. Let η G denote the efficiency of the grid transformer.

The electricity generated on-site or bought from the main grid can be stored locally in a battery. This battery has a storage capacity of c B , a charging efficiency of η C and a discharging efficiency of η D . Moreover, the maximum amount of energy than can be charged into (resp. discharged from) the battery during a micro-period is denoted by m C (resp. m D ). Let P B 0 denote the amount of electricity stored in the battery at the beginning of the planning horizon.

See Table 1 for a summary of the parameters in the problem.

Our aim is to build an integrated energy supply and industrial production plan minimizing the total production and energy procurement cost of the plant. This plan should indicate when and how many finished products should be produced, when and how much energy should be traded with the grid, charged into the battery or discharged from it. Moreover, it should guarantee that the customers' demand is met without any delay at the end of each macro-period and that the balance between energy supply and demand is reached during each micro-period.

Mathematical formulation

This section describes a new modeling approach for this problem. The basic idea of the proposed model is to exploit the energy-oriented fine discretization of the planning horizon to build the production plan, thus avoiding the introduction of a third production-oriented fine discretization as done by Wichmann, Johannes, and Spengler (2019b). The main advantage of this approach is that computing the energy consumed by production-related operations within each micro-period can be done much easier than in the EOGLSP model introduced by Wichmann, Johannes, and Spengler (2019b). However, as discussed in Subsection 4.1, it requires the management of changeover overlapping multiple micro-periods.

Model description

Using an exogenously imposed time structure consisting of a large number of short micro-periods to compute a production plan leads to the definition of a small-bucket lot-sizing problem. Among the small-bucket lot-sizing models investigated in the literature, the Proportional Lot-sizing and Scheduling Problem or PLSP introduced by [START_REF] Drexl | Proportional lotsizing and scheduling[END_REF] is the most flexible one in terms of the number of different product types that may produced in a period. We thus propose in this work to rely on the same modeling assumptions as the ones used in the PLSP. More precisely, we assume that at most one changeover operation may start, i.e., that the machine may be processing for at most two different types of products, within each micro-period.

However, our problem involves sequence-dependent changeover times, the duration of which may be greater than the length of a micro-period. A changeover may thus overlap several micro-periods. Moreover, the number of micro-periods overlapped by a changeover may not be determined in advance but will rather depend on when the changeover starts within the first overlapped period. Extensions of the PLSP with positive changeover times were investigated among others by [START_REF] Drexl | Proportional lotsizing and scheduling[END_REF][START_REF] Tempelmeier | Dynamic multi-machine lotsizing and sequencing with simultaneous scheduling of a common setup resource[END_REF][START_REF] Stadtler | Multi-level single machine lot-sizing and scheduling with zero lead times[END_REF][START_REF] Kaczmarczyk | Modelling set-up times overlapping two periods in the proportional lot-sizing problem with identical parallel machines[END_REF]) but the changeovers were restricted to overlap at most one or two micro-periods. [START_REF] Kaczmarczyk | Explicit modeling of multi-period setup times in proportional lot-sizing and scheduling problem with variable capacity[END_REF] recently introduced an MILP model for the PLSP with changeovers overlapping any number of micro-periods but assumed sequence-independent changeover times. A direct extension of this work to our setting involving sequence-dependent changeover times leads to the formulation of a large-size MILP. Our preliminary computational results showed that this MILP is numerically difficult to solve even for small-size instances. We thus propose in what follows a new extension of the model introduced by [START_REF] Kaczmarczyk | Explicit modeling of multi-period setup times in proportional lot-sizing and scheduling problem with variable capacity[END_REF] to accommodate sequence-dependent changeover times with a reduced number of variables and constraints.

We now discuss the set of binary and continuous decision variables introduced to formulate the problem as an MILP. Let us first focus on the decision variables needed to build the production plan. We start by introducing the continuous decision variables recording the inventory level of each product at the end of each macro-period.

• I j,t : inventory level of product t at the end of macro-period t.

We then introduce a set of binary variables to represent the status of the resource at the end of each micro-period and follow the changeovers between products. The machine is said to be 'processing for product j' if it is either set up for this product (manufacturing it or preserving its setup status for it) or if it is undergoing a changeover from another product j ′ ̸ = j to product j. We thus introduce:

• Y j,r = 1 if the machine is processing for item j at the end of micro-period r, and Y j,r = 0 otherwise.

In the EOPLSP model with multi-period changeovers, the number of periods overlapped by a changeover is not fixed but rather depends on the point in time at which the changeover starts in the first overlapped micro-period, i.e., depends on production planning decisions. Let us denote by qu j ′ ,j = ⌊st j ′ ,j /l e ⌋ the quotient of the modular division of the changeover time from j ′ to j by the duration l e of a micro-period and by re j ′ ,j = st j ′ ,j -qu j ′ ,j • l e the remainder of this division. For the sake of simplicity, we set qu j ′ ,j = ⌊st j ′ ,j /l⌋ -1 and re j ′ ,j = l e in case st j ′ ,j is an integer multiple of l e . A changeover from j ′ to j will overlap qu j ′ ,j + 1 if the time devoted to the changeover in the first overlapped period is greater than re j ′ ,j , and qu j ′ ,j + 2 periods otherwise. In order to distinguish between these two situations, we thus define:

• W 1 j ′ ,j,r = 1 if a changeover from product j ′ ̸ = j to product j begins in microperiod r and overlaps qu j ′ ,j + 1 micro-periods, W 1 j ′ ,j,r = 0 otherwise. • W 2 j ′ ,j,r = 1 if a changeover from product j ′ ̸ = j to product j begins during micro-period r and overlaps qu j ′ ,j + 2 micro-periods, W 2 j ′ ,j,r = 0 otherwise. We then introduce a set of continuous decision variables to determine how the capacity of the machine is used in each micro-period r.

• Q j,r : Production quantity of product j during micro-period r.

• L P j,r : Time spent by the machine preserving its setup status for product j during micro-period r.

• L F j ′ ,j,r : Time devoted to carrying a changeover from product j ′ ̸ = j to product j during micro-period r.

• B j,r : Time reserved for processing product j in micro-period r before a changeover from product j to product j ′ ̸ = j starts, if this happens; B j,r = 0 if no changeover from product j starts during r. • A j,r : Time reserved for processing product j in micro-period r after a changeover from product j ′ ̸ = j to product j starts if this happens; A j,r = 0 if no changeover to product j starts during r.

Note that B j,r (resp. A j,r ) represents the total processing time devoted to j in r before a changeover from (resp. to) j and may include changeover, production and preserving time.

Finally, we define a set of continuous variables to track the energy consumption, trading and storage in each micro-period:

• P U r : Energy consumed by the factory in micro-period r. • P GP r : Electricity purchased from the grid in micro-period r. • P GS r : Electricity sold to the grid in micro-period r. • P C r : Electricity charged into the battery within micro-period r. • P D r : Electricity discharged from the battery within micro-period r. • P B r : Electricity stored in the battery at the end of micro-period r. See Table 2 for a summary of the decision variables introduced in the formulation.

Illustrative example Before going on with the presentation of the MILP formulation, we show on a small illustrative example how the value of the variables involved in the modeling of the multi-period changeovers evolve over time. This small example involves two products (indexed resp. by 1 and 2) and a set of R = 4 micro-periods of length l e = 80. The changeover time from product 2 to product 1 is st 2,1 = 140: we thus have qu 2,1 = 1 and re 2,1 = 60. Assume that the resource is setup for product 2 at the end of micro-period 1 and that a changeover operation from product 2 to product 1 starts during micro-period 2. Depending on when this changeover starts within micro-period 2, the changeover may overlap qu 2,1 + 2 = 3 or qu 2,1 + 1 = 2 micro-periods. These two cases are presented in Figure 2.

• Case 1: see Figure 2(a) and Table 3 (a). In this case, the changeover starts rather late during micro-period 2. Namely, the time devoted to product 1 after the changeover, A 1,2 , is equal to 40 and is thus smaller than re 2,1 = 60. Consequently, the changeover overlaps qu 2,1 + 2 = 3 micro-periods and ends in micro-period 4. We thus have W 1 2,1,2 = 0 and W 2 2,1,2 = 1. B 1,2 = 40 gives the time devoted to processing product 2 in micro-period 2 before the changeover starts whereas the time devoted to carrying out the changeover from 2 to 1 in each micro-period is given by L F 2,1,r . • Case 2: see Figure 2(b) and Table 3(b). In this case, the changeover starts rather early during micro-period 1. Namely, the time devoted to product 1 after the changeover, A 1,2 , is equal to 70 and is thus greater than re 2,1 = 60. Consequently, the changeover overlaps qu 2,1 + 1 = 2 micro-periods and ends in micro-period 3. We thus have W 1 2,1,2 = 1 and W 2 2,1,2 = 0. B 1,2 = 20 gives the time devoted to processing product 2 in micro-period 2 before the changeover starts whereas the time devoted to carrying out the changeover from 2 to 1 in each micro-period is given by L F 2,1,r .

MILP formulation

We now provide the MILP formulation proposed for the integrated energy supply and industrial production planning problem.

Objective function

min r∈R j ′ ∈J j∈J \{j ′ } δ∈{1,2} f c j ′ ,j W δ j ′ ,j,r + t∈T j∈J hc j I j,t + r∈R j∈J pc j L P j,r + r∈R (π GP r • P GP r -π GS r • P GS r ).
(

) 1
The objective is to minimize the sum of the production-related and energy-related costs over the whole planning horizon. The production-related costs comprise the changeover costs, the inventory holding costs and the preserving costs. The energy-related costs are computed as the difference between the cost of buying electricity from the main grid and the revenue obtained by selling extra electricity to the grid.

Inventory balance and demand satisfaction

I j,t-1 + r∈Rt Q j,r = I j,t + d j,t , ∀j ∈ J , t ∈ T . (2)
Equations ( 2) define the inventory balance at the end of each macro-period t. Together with constraints I j,t ≥ 0, they ensure that demand is satisfied on time.

Links between the processing and changeover variables j∈J Y j,r = 1, ∀r ∈ R.

(3)

Y j,r -Y j,r-1 = j ′ ∈J \{j} δ∈{1,2} (W δ j ′ ,j,r -W δ j,j ′ ,r ), ∀j ∈ J , ∀r ∈ R. (4) j∈J j ′ ∈J \{j} δ∈{1,2} W δ j,j ′ ,r ⩽ 1, ∀r ∈ R.
(5)

Y j,r-1 + j ′ ∈J \{j} δ∈{1,2} W δ j ′ ,j,r ⩽ 1, ∀j ∈ J , ∀r ∈ R. (6)
Equalities (3) ensure that the machine is processing for exactly one product at the end of each micro-period, i.e., that it is either set up for product j (manufacturing it or preserving its setup status for it) or undergoing a changeover to this product. Constraints (4) link the processing variables to the changeover variables. Thus, if Y j,r -Y j,r-1 = 0, the processing status of the machine with respect to product j is kept unchanged and no changeover to or from product j may start in micro-period r.

If Y j,r -Y j,r-1 = -1, we have Y j,r-1 = 1 and Y j,r = 0, which means that a changeover from product j to another product starts in r and one of the binary variables W δ j,j ′ ,r should be equal to 1. Finally, if Y j,r -Y j,r-1 = 1, then Y j,r-1 = 0 and Y j,r = 1, which indicates that a changeover to product j starts during micro-period r and one of the binary variables W δ j ′ ,j,r equals one. Inequalities ( 5) enforce that at most one changeover starts during a micro-period, which corresponds to the basic modeling assumption on which the PLSP is based. Constraints (6) are simple valid inequalities. They state that if the machine is already processing for product j at the end of micro-period r -1, then no changeover to this product may start during micro-period r.

Aggregate and disaggregate capacity

j∈J k j Q j,r + j ′ ∈J \{j} L F j ′ ,j,r + L P j,r = l e , ∀r ∈ R. (7)
Constraints ( 7) are aggregate capacity constraints. They state that the total time available on the machine in a micro-period r, l e , is equal to the sum, over all products j ∈ J , of the time spent manufacturing product j, k j Q j,r , the time spent carrying out a changeover to j, j ′ ∈J \{j} L F j ′ ,j,r , and the time spent preserving the setup status for j, L P j,r .

j∈J

(B j,r + A j,r ) = l e • j∈J j ′ ∈J \{j} δ∈{1,2} W δ j,j ′ ,r , ∀r ∈ R. (8) B j,r ⩽ l e j ′ ∈J \{j} δ∈{1,2} W δ j,j ′ ,r , ∀j ∈ J , r ∈ R, (9)
A j,r ⩽ l e j ′ ∈J \{j} δ∈{1,2}

W δ j ′ ,j,r , ∀j ∈ J , r ∈ R. (10)
Constraints ( 8)-( 10) are used to determine the values of variables A j,r and B j,r . Thus, Equalities (8) ensure that the capacity of a micro-period is distributed between the time devoted to the product processed before the changeover and the time devoted to the product processed after the changeover if and only if a changeover starts during micro-period r. Inequalities (9) guarantee that the time reserved for product j before a changeover starts, B j,r , is positive only if a changeover from product j to another product j ′ starts during micro-period r. Similarly, Inequalities (10) guarantee that the time reserved for product j after a changeover starts, A j,r , is positive only if a changeover to product j starts during micro-period r.

k j Q j,r + j ′ ∈J \{j} L F j ′ ,j,r +L P j,r ⩽ l e   Y j,r - j ′ ∈J \{j} δ∈{1,2} W δ j ′ ,j,r   +A j,r +B j,r , ∀j ∈ J , r ∈ R. (11)
Inequalities ( 11) are disaggregate constraints limiting the amount of machine capacity that each product j may use in a given micro-period r. Here, four different situations may occur:

• The machine is processing for product j at the end of micro-period r (Y j,r = 1) and no changeover to j starts during r ( j ′ ∈J \{j} δ∈{1,2} W δ j ′ ,j,r = 0, B j,r = A j,r = 0). The machine was thus already processing for product j at the end of micro-period r -1 and all the capacity available in r, l e , is devoted to j (for manufacturing it, preserving the setup status for it or carrying out a changeover to j started in a previous micro-period).

• The machine is processing for product j at the end of micro-period r (Y j,r = 1) and a changeover to j starts during r ( j ′ ∈J \{j} δ∈{1,2} W δ j ′ ,j,r = 1). We have B j,r = 0 and A j,r ≥ 0. The capacity of the machine devoted to j in r should be less than A j,r .

• The machine is not processing for product j at the end of micro-period r (Y j,r = 0) and a changeover from j starts during r. We have j ′ ∈J \{j} δ∈{1,2} W δ j ′ ,j,r = 0, B j,r ≥ 0 and A j,r = 0. The capacity of the machine devoted to j in r should be less than B j,r .

• The machine is not processing for product j at the end of micro-period r and no changeover to or from j starts during r. We have Y j,r = 0, j ′ ∈J \{j} δ∈{1,2} W δ j ′ ,j,r = 0 and A j,r = B j,r = 0 so that no machine capacity is available for j in r.

Link between the changeover time variables and the number of microperiods overlapped by a changeover

L F j ′ ,j,r ⩾ re j ′ ,j W 1 j ′ ,j,r , ∀(j ′ , j) ∈ J 2 : j ̸ = j ′ , r ∈ R. ( 12
)
L F j ′ ,j,r ⩽ (re j ′ ,j -l e )W 2 j ′ ,j,r + l e , ∀(j ′ , j) ∈ J 2 : j ̸ = j ′ , r ∈ R. (13)
Constraints ( 12)-( 13) are used to enforce the consistency between the number of micro-periods overlapped by a changeover and the time devoted to carrying out this changeover in the period in which it starts. Inequalities (12) thus state that if a changeover from product j ′ ̸ = j to j overlapping qu j ′ ,j + 1 micro-periods starts in micro-period r, then the time reserved for this changeover in r should be greater than the remainder re j ′ ,j . Similarly, Constraints (13) impose that if a changeover from product j ′ ̸ = j to j overlapping qu j ′ ,j + 2 micro-periods starts in micro-period r, then the time devoted to this changeover during r should be less than the remainder re j ′ ,j .

Management of changeovers overlapping multiple periods

qu j ′ ,j τ =0 L F j ′ ,j,r+τ ⩾ st j ′ ,j W 1 j ′ ,j,r , ∀(j ′ , j) ∈ J 2 : j ̸ = j ′ , r ∈ R. ( 14
) qu j ′ ,j +1 τ =0 L F j ′ ,j,r+τ ⩾ st j ′ ,j W 2 j ′ ,j,r , ∀(j ′ , j) ∈ J 2 : j ̸ = j ′ , r ∈ R. ( 15
) qu j ′ ,j τ =0 L F j ′ ,j,r+τ ⩽ (st j ′ ,j -(qu j ′ ,j + 1) * l e ) W 1 j ′ ,j,r + (qu j ′ ,j + 1)l e ∀(j ′ , j) ∈ J 2 : j ̸ = j ′ , r ∈ R. ( 16
)
qu j ′ ,j +1 τ =0 L F j ′ ,j,r+τ ⩽ (st j ′ ,j -(qu j ′ ,j + 1) * l e ) W 2 j ′ ,j,r + (qu j ′ ,j + 1)l e ∀(j ′ , j) ∈ J 2 : j ̸ = j ′ , r ∈ R. ( 17
)
Constraints ( 14)-( 17) make sure that the actual duration of a changeover from j ′ to j equals the prescribed changeover time st j ′ ,j . Inequalities ( 14) thus impose that if a changeover from product j ′ to product j begins in micro-period r and overlaps qu j ′ ,j + 1 micro-periods, then the total time devoted to this changeover over micro-periods r, ..., r + qu j ′ ,j should be greater than st j ′ ,j . Inequalities (15) are similar constraints for the situation where the changeover overlaps qu j ′ ,j + 2 micro-periods. Inequalities (16) (resp. ( 17)) make sure that the cumulative time devoted to a changeover from j ′ ̸ = j to j over micro-periods r, ..., r + qu j ′ ,j (resp. r, ..., r + qu j ′ ,j + 1) does not exceed st j ′ ,j if a changeover from j ′ ̸ = j to j overlapping qu j ′ ,j + 1 (resp. qu j ′ ,j + 2) micro-periods starts during r.

L F j ′ ,j,r ⩾ A j,r -l e (Y j,r -W 2 j ′ ,j,r ), ∀(j ′ , j) ∈ J 2 : qu j ′ ,j = 0, r ∈ R. (18)

L F j ′ ,j,r ⩾ A j,r -l e (Y j,r -W 1 j ′ ,j,r -W 2 j ′ ,j,r ), ∀(j ′ , j) ∈ J 2 : qu j ′ ,j ⩾ 1, r ∈ R. (19) L F j ′ ,j,r ≥ W 2 j ′ ,j,r-1 l e , ∀(j ′ , j) ∈ J 2 : qu j ′ ,j = 1, r ∈ R. (20) L F j ′ ,j,r ≥ qu j ′ ,j τ =1 W 2 j ′ ,j,r-τ + qu j ′ ,j -1 τ =1 W 1 j ′ ,j,r-τ l e , ∀(j ′ , j) ∈ J 2 : qu j ′ ,j ⩾ 2, r ∈ R. ( 21
)
A changeover operation started in r should not be interrupted but rather should continue till its end. This can be enforced by using lower bounds on the time devoted to carrying out this changeover in r and in the subsequent micro-periods. Constraints ( 18)-( 19) thus impose that L F j ′ ,j,r is greater than A j,r , i.e. that all the time available for j during r is devoted to carrying out the changeover (thus forbidding any manufacturing or preserving for j) whenever a changeover from j ′ to j overlapping at least two periods starts during r. Note how these constraints are inactive when r is not the first period of a changeover from j ′ to j. Moreover, recall that Constraints ( 14) and ( 17) above impose L F j ′ ,j,r = st j ′ ,j in the case of a changeover totally included in micro-period r (i.e., in case qu j ′ ,j = 0 and W 1 j ′ ,j,r = 1). Moreover, Constraints ( 20)-( 21) guarantee that micro-period r is fully devoted to carrying out the changeover from j ′ to j, i.e., guarantee that L F j ′ ,j,r is greater than l e , in case r is strictly within the time interval covered by a changeover overlapping at least 3 periods.

Energy consumption

P U r = J j ′ =0 J j=0 e F j ′ ,j • L F j ′ ,j,r + J j=1 e M j • Q j,r + J j=0 e P j • L P j,r , ∀r ∈ R. (22)
Equation ( 22) computes the total energy needed by the factory in each energy-oriented micro-period r as a function of the changeover, manufacturing and preserving activities occurring during r.

Energy supply

P U r + P GS r η G + P C r η C = η G • P GP r + p R r + η D • P D r , ∀r ∈ R. ( 23
)
P B r = P B r-1 + P C r -P D r , ∀r ∈ R. ( 24
) 0 ⩽ P B r ⩽ c B , ∀r ∈ R. (25)
P C r ⩽ m C , ∀r ∈ R. (26)
P D r ⩽ m D , ∀r ∈ R. (27)
Constraints ( 23)-( 27) aim at building a feasible energy supply plan in each energyoriented micro-period. Equalities (23) ensure that the energy supply and demand is balanced within each micro-period. The energy demand consists in the energy consumed by the manufacturing process P U r , the energy sold to the main grid P GS r ηG and the energy charged into the battery P C r ηC . Note how, since part of the energy is lost within the grid transformer, the amount of energy needed to sell P GS r to the main grid is given by P GS r /η G . Similarly, the actual energy needed to charge P C r into the battery is computed by P C r /η C due to the losses during the battery charging process. The energy supply consists in the energy purchased from the grid η G • P GP r , the energy generated from the renewable resources p R r and the energy discharged from the battery η D • P D r . Note that, when we buy P GP r from the main grid, only η G • P GP r energy is really available for the plant. The same applies for the energy discharging process. The energy balance in the battery is defined by Equations ( 24). The amount of energy stored in the battery at the end of micro-period r is equal to the amount stored at the beginning of r plus the amount charged into the battery during r minus the amount discharged from the battery in r. Constraints (25) impose that the amount of energy stored in the battery does not exceed its capacity. Inequalities ( 26) and ( 27) are the battery maximum charging and discharging rate constraints.

Initial and final states of the system I j,0 = I j,0 , ∀j ∈ J , (28) Y 0,0 = 1, (29)

P B 0 = P B 0 . (30)
I j,T ⩾ I j,0 , ∀j ∈ J , (31)

P B T ⩾ P B T . (32)
Constraints ( 28)-( 30) describe the initial state of the system. Constraints (31) ensure that the inventory level of each product at the end of the planning horizon is at least equal to its initial inventory level. Constraint ( 32) is a similar constraint on the final battery state of charge.

Comparison with the EOGLSP model

We now briefly discuss the main differences between the EOPLSP model presented in Subsection 4.2 and the EOGLSP model introduced by Wichmann, Johannes, and Spengler (2019a).

First, thanks to the fact that the EOPLSP model use the energy-oriented fixed time grid to plan the industrial activities, the amount of energy consumed in each energy-oriented micro-period r can be directly computed: see Equalities ( 22). Consequently, there is no need to introduce the |S|R binary variables used in the EOGLSP model, where |S| denotes the total number of micro-periods of flexible length in the production-oriented time grid, to track when the energy consumption actually takes place. This explains the drastic reduction in the MILP size obtained when using the EOPLSP model rather than the EOGLSP model.

Second, the two models rely on different modeling restrictions. Namely, on one hand, the EOPLSP model limits the number of changeovers during an energy-oriented micro-period to one whereas there is no such restriction in the EOGLSP model. The EOGLSP model is thus less restrictive than the EOPLSP model regarding this aspect and may theoretically provide less expensive production plans. However, on the other hand, the EOGLSP model does not allow changeovers overlapping two macro-periods while there is no such restriction in the EOPLSP model. Thus, theoretically speaking, none of the two models has a feasible space strictly included in the feasible space of the other one and is guaranteed to find a better optimal solution. Yet, in practice, the computational results presented in Section 5 show that, on the instances which could be solved to guaranteed optimality using both models, the gap between the optimal values is very small.

Valid inequalities

This subsection discusses valid inequalities that may be used to strengthen the MILP formulation of the EOPLSP model. The objective is to improve the quality of the lower bounds provided by the continuous relaxation of the MILP formulation and to speed up the resolution of the problem by a branch-and-bound type algorithm. [START_REF] Barany | Strong formulations for multi-item capacitated lot sizing[END_REF] proposed a set of (l, S) inequalities strengthening the formulation of the single-product uncapacitated lot-sizing (ULS) problem and enabling to obtain a full description of the convex envelope of its feasible space.

However, the ULS is essentially a large-bucket model in which demand may occur at the end of each period and the setup state of the resource is not conserved between two consecutive periods. Thus, the (l, S) inequalities cannot be directly applied to our small-bucket EOPLSP model. We thus discuss how to extend them.

To gain a first intuition of the proposed extension, let us first consider a given product j ∈ J and a single macro-period t ∈ T . The following simple valid inequality holds:

r∈Rt Q j,r ⩽ d j,t Y j,(t-1)R + r∈Rt j ′ ∈J \{j} δ∈{1,2} W δ j ′ ,j,r + I j,t . (33)
Two situations may namely arise:

• either Y j,(t-1)R + r∈Rt j ′ ∈J \{j} δ∈{1,2} W δ j ′ ,j,r = 0. In this case, the machine is not processing for j at the beginning of macro-period t and no changeover to j occurs during t. Hence, no production for j may occur during t. We thus have Y j,r = 0 and Q j,r = 0 for all r ∈ R t . The inequality reduces to 0 ≤ I j,t and trivially holds.

• or Y j,(t-1)R + r∈Rt j ′ ∈J \{j} δ∈{1,2} W δ j ′ ,j,r ≥ 1. Thanks to the inventory balance constraints (2) and the fact that I j,t-1 ≥ 0, we have r∈Rt Q j,r ≤ d j,t +I j,t . Using the assumption that Y j,(t-1)R

+ r∈Rt j ′ ∈J (W 1 j ′ ,j,r +W 2 j ′ ,j,r ) ≥ 1 gives r∈Rt Q j,r ⩽ d j,t Y j,(t-1)R + r∈Rt j ′ ∈J δ∈{1,2} W δ j ′ ,j,r + I j,t .
Let us now consider a more general case.

Proposition 1. Let j ∈ J be a product, w ∈ T be a macro-period, W = {1, ..., w} and Φ ⊆ W be two subsets of macro-periods. The following inequality holds for the EOPLSP model:

τ ∈Φ r∈Rτ Q j,r ⩽ τ ∈Φ d j,τ w Y j,(τ -1)R + r∈Rτ j ′ ∈J δ∈{1,2} W δ j ′ ,j,r + I j,w , (34)
with d j,τ,w = w t=τ d j,t the cumulative demand for product j over the time interval [τ, w].

Proof. Let us consider a production plan (Y, W, Q) solution of the EOPLSP model.

If Y j,(τ -1)R + r∈Rt j ′ ∈J \{j} δ∈{1,2} W δ j ′ ,j,r = 0 for all r ∈ R τ and τ ∈ Φ, we have Y j,r = 0 for all r ∈ R τ and τ ∈ Φ. No production for j may occur during the macro-periods τ belonging to Φ. We thus have Q j,r = 0 for all r ∈ R τ and τ ∈ Φ. This inequality amounts to I j,w ≥ 0 and is trivially satisfied.

Otherwise, let τ 0 = min{τ ∈ Φ : Y j,(τ -1)R + r∈Rt j ′ ∈J \{j} δ∈{1,2} W δ j ′ ,j,r ⩾ 1} be the first macro-period in Φ in which production for j may occur. We have:

τ ∈Φ r∈Rτ Q j,r ⩽ τ ∈Φ∩{τ0,...,w} r∈Rτ Q j,r (35a)
⩽ w τ =τ0 r∈Rτ Q j,r (35b)
⩽ d j,τ0w + I j,w (35c)
⩽ d j,τ0w Y j,(τ0-1)R + r∈Rτ 0 j ′ ∈J \{j} δ∈{1,2} W δ j ′ ,j,r ) + I j,w (35d)
⩽ τ ∈Φ d j,τ w Y j,(τ -1)R + r∈Rτ j ′ ∈J \{j} δ∈{1,2} W δ j ′ ,j,r + I j,w . (35e)
(35a) follows from the definition of τ 0 and (35b) from the non-negativity of the production variables Q j,r . ( 35c) is obtained by summing up the inventory balance constraints (2) over the time interval [τ 0 , w] and using the fact that I j,τ0-1 ≥ 0. ( 35d) is obtained thanks to the definition of τ 0 as a period in which production for j occurs, i.e., in which Y j,(τ0-1)R + r∈Rτ 0 j ′ ∈J \{j} δ∈{1,2} W δ j ′ ,j,r ≥ 1. (35e) follows from the non-negativity of variables Y and W .

Note that there is a valid inequality for each subset of macro-periods Φ, meaning that the number of valid inequalities (34) is exponential. It is thus not possible to add all of them a priori into the formulation as this would lead to numerical difficulties. We need to solve the associated separation problem which is defined as follows. Given a feasible solution (Y * , W * , Q * ) of the continuous relaxation of the EOPLSP model, either identify an inequality (34) violated by (Y * , W * , Q * ) or prove that no such inequality exists.

In order to more easily solve this separation problem, we reformulate inequalities (34) by using the fact that I j,w = I j,0 + τ ∈W r∈Rτ Q j,r -d j,1w . This gives:

d j,1w -I j,0 ≤ τ ∈W\Φ r∈Rτ Q j,r + τ ∈Φ d j,τ w Y j,(τ -1)R + r∈Rτ j ′ ∈J \{j} δ∈{1,2} W δ j ′ ,j,r (36)
Clearly, for a given product j and macro-period w, the inequality most violated by (Y * , W * , Q * ), if it exists, will be the one with the smallest possible value of the right-hand side of (36).

The following polynomial-time algorithm enables to find, for each product j and macro-period w, the subset Φ minimizing this right-hand side, thus enabling to solve the separation problem defined above.

Note that similar valid inequalities may be obtained for the EOGLSP model introduced by Wichmann, Johannes, and Spengler (2019a) using the same reasoning as the one used above for the EOPLSP model.

Algorithm 1 Separation algorithm for the extended (l, S) inequalities for j = 1, 2, . . . , J do for w = 1, 2, . . . , T do

α j,w := w τ =1 min r∈Rτ Q * j,r , d j,τ w Y * j,(τ -1)R + r∈Rτ j ′ ∈J \{j} δ∈{1,2} W δ, * j ′ ,j,r if α j,w < d j,1w -I j,0 then
The inequality (36) corresponding to the subset defined by Φ =

τ ∈ W : r∈Rτ Q * j,r > d j,τ w Y * j,(τ -1)R + r∈Rτ j ′ ∈J \{j} δ∈{1,2} W δ, * j ′ ,j,r
is violated for product j and macro-period w. else (Y * , W * , Q * ) complies with all valid inequalities (36) defined for product j and macro-period w. end if end for end for

Numerical experiments

Numerical experiments were conducted to assess the performance of the proposed EO-PLSP model and compare it with one of the EOGLSP model investigated by Wichmann, Johannes, and Spengler (2019a).

The procedure used to randomly generate instances is described in Subsection 5.1. The corresponding computational results are reported in Subsection 5.2. Finally, some managerial insights on the cost savings that may be obtained through a decentralized energy system are provided in Subsection 5.3.

Instances

The parameter setting used to randomly generate instances is mostly based on the numerical values provided by Wichmann, Johannes, and Spengler (2019b). Furthermore, as this paper involves energy losses neither in the grid transformer nor during the battery charging/discharging process, we use the values provided in [START_REF] Zhang | Operational optimization of a grid-connected factory with onsite photovoltaic and battery storage systems[END_REF] and [START_REF] Golpîra | Robust smart energy efficient production planning for a general job-shop manufacturing system under combined demand and supply uncertainty in the presence of grid-connected microgrid[END_REF] for these parameters.

Note that the numerical values used in Wichmann, Johannes, and Spengler (2019a) were based on a potential application case of the EOGSLP model in a factory producing automotive steel parts by a warm/hot forming process (Wichmann, Johannes, and Spengler 2019b). In this process, a gas or electrical oven first heats the blanks from the ambient temperature to a high temperature, a forming machine then processes the heated blanks to obtain parts of complex geometric shapes. Steel parts of various size and weight may have to be processed on this resource and changeover operations have to be carried between the processing of lots of different product types. Moreover, in case there is some non-productive time between two production lots for the same product type, it may be preferable to keep the oven at the right temperature, i.e., to preserve its setup status, rather than turning it off and heating it up again from the ambient temperature.

We thus consider a production system involving a single machine producing a set of J products on a finite planning horizon spanning several days. Each day comprises a morning shift (from 6 a.m. to 2 p.m.), an afternoon shift (from 2 p.m. to 10 p.m.), and an evening shift (from 10 p.m. to 6 a.m. of the next day). Demand for the finished products arises at the end of each shift. We thus discretize the planning horizon into T macro-periods, each one corresponding to an eight-hour shift and lasting l = 480 minutes. Macro-periods with an index t ≡ 1 mod 3 correspond to morning shifts, the ones with an index t ≡ 2 mod 3 correspond to afternoon shifts, and the ones with an index t ≡ 0 mod 3 are night shifts. We consider four sets of instances of various sizes: more precisely, we use values of (J, T ) in {(3, 2), (3, 6), (4, 9), (4, 12)}.

Producing one unit of product j ∈ J \ {0} takes k j = 0.05 minutes and consumes e M j = 0.136kWh of electricity. Preserving the machine status for product j during one hour costs pc j = 12 and consumes an energy amount equal to e P j = 5kWh. A changeover from product j ′ ̸ = j to product j costs f c j ′ ,j = 200 and consumes e F j ′ ,j = 10kWh per hour.

Regarding the changeover times, four different patterns are considered:

• Pattern 1: in this pattern, the value of the changeover times are set to the ones used by Wichmann, Johannes, and Spengler (2019b). We thus have st j ′ ,j = 60min for any pair of products j ′ ̸ = j. • Pattern 2: Changeover times remain product and sequence-independent but are set to a larger value than in Pattern 1. We set st j ′ ,j = 0.75 * l e /(J + 1) for any pair of products j ′ ̸ = j. • Pattern 3: Changeover times remain product and sequence-independent but are set to a smaller value than in Pattern 1. We set st j ′ ,j = 0.25 * l e /(J + 1) for any pair of products j ′ ̸ = j. • Pattern 4: Changeover times are product and sequence-dependent. The changeover time from product j ′ ̸ = j to product j is randomly generated following the Normal distribution of mean 0.5l e /(J + 1) and a standard deviation of 0.25l e /(J + 1). st j ′ ,j is set to zero if the randomly generated value is negative.

Thus, we have: st j ′ ,j ∼ max N ( 0.5l e J+1 , 0.25l e J+1 ); 0 . The production capacity of the machine within one macro-period is equal to Cap = l/k j = 480/0.05 = 9600 units of finished product per macro-period. The demand for each item j = 1, ...J is randomly generated as follows. We consider a utilization rate of the resource equal to ρ = 0.5 (resp. ρ = 0.25 and ρ = 0.75) for Patterns 1 and 4 (resp. for Pattern 2 and Pattern 3). The mean value of the total demand to be satisfied in each macro-period t ∈ T is thus set to ρCap units of finished products. The demand d j,t for product j = 1, ..., J is randomly generated according to a Normal distribution of mean ρCap J and standard deviation ρCap 3J . In case the randomly generated value is negative, we replace it by 0. We thus have d j,t ∼ max int(N ( ρCap J , ρCap 3J )); 0 . The initial inventory of product j is randomly generated following the uniform distribution using Īj,0 ∼ randint(0, 2d j,1 ). The unit inventory holding cost of product j is set to hc j = 0.05 per macro-period.

As for the dummy product j = 0, holding and preserving cost are set to zero, i.e., hc j = pc j = 0. The corresponding unit production time, as well as the energy required for the manufacturing and preserving, are set to 0: k 0 = 0, e M 0,r = e P 0,r = 0. Regarding the energy supply and consumption, each macro-period t ∈ T is split evenly into R = 8 energy-oriented micro-periods of l e =60min length. Recall that each day is split into three macro-periods corresponding to a morning, an afternoon and an evening shift and that the production plan is assumed to start with a morning shift. Thus, the micro-period indexed by r = 1 corresponds to a time interval between 6am and 7am, i.e. to the the hour of the day indexed by 6. More generally, the micro-period indexed by r corresponds to the hour of the day indexed by h ≡ (r + 5) mod 24.

The energy price is assumed to display intraday variations but to be otherwise daily periodic. We considered two distinct levels for the energy price:

• Level H corresponding to high energy prices. In this case, we set π

GP r = π ref (r+5)(mod24) with π ref h
the unit price to buy energy from the main grid at the hour of the day indexed by h ∈ {0, .., 23}: see Table 5. The energy selling prices are 0.021 /kWh lower than the buying prices, i.e., π GS r = π GP r -0.021.

• Level L corresponds to low energy prices. In this case, we set the buying and selling prices to be one tenth of the reference prices, i.e., π GP r = 0.1π ref (r+5)(mod24) , and

π GS r = π GP r -0.0021.
The amount of on-site generated electricity depends on the weather and thus varies both within the day and from one day to the next. The value of p R r for micro-period r is randomly generated using a Normal distribution with an expected value and a standard deviation equal to G (r+5)(mod24) : see Table 5. Note how the values of G h are consistent with a power generation by PV panels. We thus have p R r ∼ max N (G (r+5)(mod24) , G (r+5)(mod24) ); 0 . The energy storage system has a capacity of c B = 500kWh. Within each microperiod r, the maximum amount of energy than can be charged into or discharged from the battery is set to m C = m D = 250kWh. The charging and discharging efficiencies, as well as the grid transformer efficiency, are set to η C = η D = η G = 0.95.

In the EOPLSP model, the production-oriented micro-periods are identical to the energy-oriented ones. In the EOGLSP model, each macro-period t may be flexibly divided into a predefined number |S t | of production-oriented micro-periods of variable length: we set |S t | = 7 for each t (refer to Wichmann, Johannes, and Spengler (2019a) for the definition of S t ).

Table 4 provides an overview of this parameters setting.

Numerical results

For each considered instance size, changeover time pattern and energy price level, we randomly generated 10 instances, resulting in a total number of 320 instances. Each instance was solved using either the EOGLSP model or the EOPLSP model (see Section 4) with a cut-and-branch algorithm. This algorithm comprises two steps. We first iteratively add valid inequalities (34) to strengthen the LP relaxation of the MILP formulation using the separation algorithm described in Subsection 4.4. When no more violated valid inequalities can be found, we solve the problem using the strengthened MILP formulation with the mathematical solver CPLEX 20.10. The implementation was done in Python. The computational experiments were carried out on a laptop running under Linux with an Intel(R) Core(TM) i7-1165G7 CPU @ 2.80GHz processor and 32 GB RAM. The time limit was set to 20 minutes.

The results are displayed in Tables 6789. Each table gives the results for the set of 80 instances corresponding to a given changeover time pattern. Within each table, each column provides the average results over the 10 instances of similar size and energy price level. When a feasible solution could not be obtained for the 10 corresponding instances within the computation time limit, we report the average results over the instances for which a feasible solution could be obtained. We provide, for each set of 10 instances:

• # VAR / # BINVAR: the number of variables and binary variables involved in the formulation, • # CONS: the number of constraints in the formulation, • Gap LP 0 : the integrality gap (before adding the extended (l, S) cuts), i.e., the relative difference between the lower bound provided by the linear relaxation of the problem and the value of the best integer feasible solution found by the solver, • # VI: the number of generated extended (l, S) cuts, • Gap LP : the integrality gap after adding the extended (l, S) cuts, • # FEAS: the number of instances for which a feasible solution is found within the time limit, • Z best : the value of the best feasible integer solution found by the solver within the time limit, • Gap M IP : the average optimality gap, i.e., the average gap between the best lower bound and the best upper bound found by the solver before the time limit is reached (note that this gap is equal to 0% in case a guaranteed optimal solution could be found), • σ Gap : the standard deviation of Gap M IP over the 10 considered instances, • Time: the average computation time (in seconds), • σ t : the standard deviation of T ime over the 10 considered instances.

We observe from the results provided in Tables 6-9 that the computational effort needed to solve the problem with an MILP solver is significantly reduced when using the EOPLSP model instead of the EOGLSP model.

Namely, for the 80 instances corresponding to (J, T ) = (3, 2), the average time needed to solve the problem to optimality is significantly reduced from 28s with the EOGLSP model to 2s with the EOPLSP model. Moreover, for the 80 instances corresponding to (J, T ) = (3, 6), we were able to obtain a guaranteed optimal solution within an average computation of 131s when using the EOPLSP model. In contrast, the branch-and-cut algorithm embedded in CPLEX20.10 could not converge to an optimal solution after 20 minutes of computation when using the EOGLSP model. As for the larger instances corresponding to (J, T ) ∈ {(4, 9), (4, 12)}, we could find a feasible solution of acceptable quality for all instances when using the EOPLSP model whereas a feasible solution was found for only 95 out of 180 instances when using the EOGLSP model. Note how the average optimality gap Gap M IP obtained with the EOPLSP model is equal to 5.01% for the 80 instances with (J, T ) = (4, 9) and to 9.26% for the 80 instances with (J, T ) = (4,12).

This improvement in the computational efficiency of the MILP solver is mainly explained by the significant decrease in the size of the MILP formulation obtained when using the EOPLSP model. For instance, for the largest considered instances with (J, T ) = (4, 12), the total number of variables (resp. of binary variables) is drastically reduced from 310,102 (resp. 11,645) in the EOGLSP model to 10,652 (resp. 5,535) in the EOPLSP model. Similarly, the number of constraints involved in the formulation is divided by around 12 when using the EOPLSP model instead of the EOGLSP model. Moreover, the computational efficiency improvement may also come from the fact that the extended (l, S) cuts generated by our cut-and-branch algorithm are slightly more efficient at strengthening the formulation of the EOPLSP model than at strengthening the one of the EOGLSP model: see how the average value of Gap LP is most often smaller for the EOPLSP model than for the EOGLSP model. The use of lower bounds of better quality may help the MILP solver converge more rapidly when using the EOPLSP model.

In terms of solution quality, as explained in Subsection 4.3, none of the EOPLSP or EOGLSP models is theoretically guaranteed to find a planning less expensive than the other model. However, over the 85 instances that could be solved to optimality using both models, we obtain production plans of equal cost for 68 instances and slightly less expensive production plans with the EOPLSP model for 17 instances : see line Z best of Tables 6789. It thus seems that, for the considered instances, the additional flexibility offered by the EOGLSP model to produce more than three products in a microperiod is not useful. This might be explained by the presence of positive changeover times between production lots and a high ratio of the changeover costs over inventory holding costs, which makes it either impossible or sub-optimal to produce using lots of such a small size that three or more of them could fit in a one-hour micro-period. Nevertheless, this observation might not hold for instances featuring zero changeover times, a smaller ratio of the changeover costs over inventory holding costs, and longer energy-oriented micro-periods.

Managerial insights on the cost savings obtained with a decentralized energy system

Installing a decentralized energy system on a production site requires a substantial financial investment and is a strategic long-term decision. Yet, this initial investment may be counterbalanced by cost savings at the operational level through a reduction of the energy supply costs. In this subsection, we try to quantify the energy cost savings that may be obtained at the operational level by installing a decentralized energy system. The resulting managerial insights may be useful for a decision-maker wishing to invest in such a system. To this aim, we consider a reference scenario in which no decentralized energy system is installed on the production site (i.e., there is no on-site generation of electricity: p R r ≡ 0 for all r ∈ R and no battery : c B = 0) and all the energy needed is directly purchased from the main grid. We then define alternative scenarios in which:

• a small or large number of photovoltaic panels are installed, leading to a low/high level of renewable electricity availability. In the low (resp. high) level scenario, p R r is generated as described in Subsection 5.1 using the values of G h (resp. 2G h ) displayed in Table 5.

• a battery of capacity c B ∈ {500, 1000}kWh is installed. This computational study relies on the 80 instances with (J, T ) = (3, 6) introduced in Subsection 5.1. For each instance and each investment scenario, we build a plan using the EOPLSP model presented in Section 4 and record the corresponding total cost Z. We denote by Z ref the cost of the plan obtained in the reference scenario (in which there is no decentralized energy system) and by Z scen the one of the plan obtained in an alternative scenario (in which a battery and/or some PV panels are installed). The relative cost saving obtained for a given instance in a given investment scenario is finally calculated as ∆C scen = (Z ref -Z scen )/Z ref . The average values, over the 80 studied instances, of the total cost of the production plan Z and of ∆C scen for the various investigated scenarios are displayed in Table 10 for the case where electricity selling to the grid is allowed and in Table 11 for the case where it is forbidden.

Results from Table 10 first show that, when energy selling to the main grid is allowed, it is possible to obtain substantial cost savings by investing in a decentralized energy system. The total costs can namely be reduced by up to 43% with a large battery capacity and a high number of PV panels. Moreover, it seems important to simultaneously invest in both PV panels and energy storage capacity as the cost savings obtained by investing in only one of these technologies are significantly smaller (around 20%) than the ones obtained when both technologies are installed. However, a more detailed analysis of these results showed that these cost savings come partly from a speculative behaviour in which the factory buys energy from the grid when its price is low, stores it for a few hours and resells it at a latter point in time when its price is high. This can be seen by looking at the value of the ratio between the total amount of energy not used to supply the industrial process over the total energy consumed by the production process, i.e., α = r=1∈R (ηGP GP r +p R r -P U r ) r=1∈R P U r . Over the three investment scenarios involving a large number of PV panels, α sharply increases from 4.46% when there is no battery (thus indicating that the energy available in the system is mostly used to power the industrial system) to 46.64% (resp. 73.66%) when there is small (resp. large) battery (thus indicating that a large share of the energy available in the system is not used to power the industrial process but rather sold to the main grid.) Nevertheless, results from Table 11 show that, even when energy selling is not allowed, investing in a decentralized energy system may still lead to significant cost savings (nearly 26% in the best case). Moreover, in this case, priority should be placed on installing a large number of PV panels and a battery with a small capacity as using a battery with a large capacity leads to a marginal improvement of the cost savings.

Conclusion and perspectives

We investigated an integrated industrial production and energy supply planning problem faced by a production site partially powered by on-site generated renewable electricity. We proposed to model the industrial production problem as a small-bucket lotsizing problem, more precisely as a PLSP with sequence-dependent changeover times overlapping multiple periods. Our numerical results show that the resulting MILP model enables to obtain good-quality plans with a computational effort much smaller than the one required by a previously published EOGLSP model. Moreover, our results also indicate that substantial cost savings may be obtained at the operational level by investing in a decentralized energy system.

The present work has several limitations. First, we considered a single-machine single-level production system. An interesting research direction could be to extend this work to a more general multi-machine multi-level setting in order to improve its practical relevancy. Second, we assumed that the amount of renewable energy provided by the local system is deterministically known whereas it is in fact difficult to accurately forecast. Thus, it may be worth to study a stochastic extension of this work in which these uncertainties are explicitly considered. Third, we focused on modeling and solving a short-term production planning problem. However, a more strategic problem could be studied to determine whether investing in a decentralized energy system will be profitable for a factory on the long run and to define the energy conversion and storage capacity of this system.

Tables Table 1: Length of a micro-period I j,0

Initial inventory of product j d j,t

Demand for product j at the end of macro-period t k j

Time needed to produce one unit of product j f c j ′ ,j Changeover cost from product j ′ to product j hc j

Unit inventory holding cost for product j pc j

Preserving cost for product j per unit of time st j ′ Changeover time from product j ′ to product j qu j ′ ,j Quotient of the modular division of the changeover time by the length of a micro-period, defined by qu j ′ ,j = ⌊st j ′ ,j /l e ⌋ re j ′ ,j Remainder of the modular division of the changeover time by the length of a micro-period, re j ′ ,j = st j ′ ,j -qu j ′ ,j • l e (if re j ′ ,j = 0, then let qu j ′ ,j = ⌊st j ′ ,j /l e ⌋ -1 and re j ′ ,j = l e ) e M j Energy needed for manufacturing one unit of product j e F j ′ ,j Energy consumed per unit of time during a changeover from product j ′ to product j e P j Energy consumed per unit of time when preserving the machine setup status for product j π GP r Unit price to purchase electricity from the grid in micro-period r π GS r Unit price to sell electricity to the grid in micro-period r p R r Electricity generated by the renewable energy system in micro-period r c B Capacity of the battery m C

Maximum amount of energy than can be charged into the battery during a micro-period m D

Maximum amount of energy that can be discharged from the battery during a micro-period η C

Efficiency of the battery charging process η D

Efficiency of the battery discharging process η G Efficiency of the grid transformer

P B 0
Initial state of charge of the battery Binary variables indicating whether the machine is processing for product j (i.e. manufacturing j, preserving its setup status for j or undergoing a changeover to j) at the end of micro-period r W δ j ′ ,j,r

Binary variables indicating whether a changeover from product j ′ ̸ = j to product j overlapping qu j ′ ,j + δ micro-periods starts during period r Q j,r Production quantity of product j during micro-period r L P j,r

Time devoted to preserving the machine setup status for product j in micro-period r L F j ′ ,j,r Time devoted to carrying out a changeover from product j ′ ̸ = j to product j during micro-period r B j,r Time reserved for processing product j in micro-period r before a changeover from product j to product j ′ starts (if this happens)

A j,r Time reserved for processing product j in micro-period r after a changeover from product j ′ to product j ′ starts (if this happens) 3.0 3.0 3.0 3.0 3.0 3.0 4.8 6.1 6.3 6.0 5.6 4.0 3.7 3.8 4.5 5.1 5.4 5.9 6.4 6.3 5.5 4.5 3.0 3.0 G h 0 0 0 0 0 0 1 4 10 18 25 27 30 30 25 15 5 2 0 0 0 0 0 0 ,395.73 6,405.12 14,145.66 14,320.12 21,426.19 22,378.52 28,029.57

  Electricity stored in the battery at the end of micro-period r Table 3: Decision variables' values for the illustrative example(a) Decision variables' values for Sce-st j ′ ,j = 60min Pattern 2: st j ′ ,j = 0.75 * l e /(J + 1)min Pattern 3: st j ′ ,j = 0.25 * l e /(J + 1)min Pattern 4: st j ′ ,j ∼ max N ( 0.5l e J+1 , 0.25l e J+1 ); 0 min Energy demand for manufacturing e M j = 0.136kWh per unit of product Energy demand for changeover e F j ′ ,j = 10kWh per minute Energy demand for preserving e P j = 5kWh per minute Battery storage capacity c B ∈ {0, 500, 1000}kWh Battery charging maximal rate m C = 250kW Battery discharging maximal rate m D = 250kW Efficiency of the battery charging η C = 0.95 Efficiency of the battery discharging η D = 0.95 Efficiency of the grid transformer η G = 0.95 Initial state of charge of the battery P B 0 = 0kWh

Table 2 :

 2 Decision variablesI j,tInventory of product j at the end of macro-period t Y j,r

Table 5 :

 5 Reference energy price and generated renewable electricity according to the hour of the day

	h	0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
	π ref h

Table 6 :

 6 Results for instances with changeover time pattern 1

	Energy (J, T )	(3, 2)	(3, 6)		(4, 9)	(4, 12)
	price Model	EOPLSP EOGLSP EOPLSP EOGLSP EOPLSP EOGLSP EOPLSP EOGLSP
		# VAR	1,353	7,302	3,753	56,714	7,953	176,656	10,527	310,102
		# BINVAR	676	675	1,828	3,311	4,115	7,232	5,435	11,645
		# CONS	2,698.7	7,268.3	7,939.5	5,9460.5	17,667	183,089	23,515.9	322,782
		Gap LP 0	14.09%	20.96%	15.86%	17.79%	22.31%	26.47%	22.75%	30.37%
	Level H	# VI Gap LP	5.7 11.55%	7.3 15.84%	35.5 9.92%	71.5 9.74%	107 13.09%	202 14.37%	144.4 13.95%	314 17.97%
		# FEAS	10	10	10	10	10	2	10	1
		Z best	6

Table 7 :

 7 Results for instances with changeover time pattern 2

	Energy (J, T )	(3, 2)	(3, 6)	(4, 9)	(4, 12)
	price Model	EOPLSP EOGLSP EOPLSP EOGLSP EOPLSP EOGLSP EOPLSP EOGLSP
		# VAR	1,433	7,302	3,833	56,714	8,078	176,656	10,652	310,102
		# BINVAR	740	675	1,892	3,311	4,215	7,232	5,535	11,645
		# CONS	2,985.8	7,269.7	8,741.3	59,445.9	19,517 183,050.3 25,988.6 322,722.7
		Gap LP 0	24.87%	33.03%	24.09%	28.36%	34.97%	42.72%	34.07%	43.58%
	Level H	# VI Gap LP	5.8 21.87%	8.7 23.99%	37.3 12.16%	56.9 13.07%	107 15.71%	163.3 23.81%	167.6 13.30%	254.7 22.60%
		# FEAS	10	10	10	10	10	9	10	3
		Z best	2,600.05 2,602.68 5,461.52 5,489.81 7,621.79 8,414.75 9,681.79 10,277.47
		Gap M IP	0.00%	0.01%	0.01%	1.83%	5.73%	17.31%	5.68%	17.83%
		σ Gap	0.00%	0.00%	0.00%	2.41%	1.63%	6.19%	1.94%	3.65%
		Time	2.16	74.20	189.28	1,102.43 1,200.86 1,202.16 1,201.06 1,204.23
		σt	0.48	23.05	159.88	205.30	0.13	0.38	0.31	0.36
		# VAR	1,433	7,302	3,833	56,714	8,078	176,656	10,652	310,102
		# BINVAR	740	675	1,892	3,311	4,215	7,232	5,535	11,645
		# CONS	2,988.1	7,269.6	8,759.2	59,462.7 19,567.9 183,109.9 26,051.9 322,820
		Gap LP 0	41.30%	51.34%	48.67%	55.35%	61.45%	68.22%	63.21%	67.91%
	Level L	# VI Gap LP	8.1 28.30%	8.6 31.99%	55.2 25.57%	73.7 29.71%	157.9 28.00%	222.9 37.60%	230.9 27.29%	352 36.43%
		# FEAS	10	10	10	10	10	8	10	2
		Z best	1,159.27 1,159.27 2,496.53 2,515.63 4,066.08 4,584.66 5,240.79 5,899.01
		Gap M IP	0.00%	0.01%	0.01%	3.87%	9.13%	29.09%	17.14%	30.64%
		σ Gap	0.00%	0.00%	0.00%	3.16%	4.52%	5.62%	2.27%	0.96%
		Time	3.17	31.31	167.36	1,202.33 1,201.99 1,202.69 1,201.28 1,204.45
		σt	1.81	10.43	124.08	0.35	0.08	0.21	0.16	0.38

Disclosure statement

The authors report there are no competing interests to declare.

Data availability statement

The data that support the findings of this study are available from the corresponding author, R.L., upon reasonable request.

Notes on contributors

Ruiwen Liao received a master degree in Industrial Engineering in 2021 from Ecole Centrale Pékin, Beihang University, and is currently a Ph.D. candidate in Operations Research in the Laboratory for interdisciplinary research in numerical sciences (LISN) at the Université Paris Saclay (France). Her research focuses on energy-efficient manufacturing, lot-sizing, and combinatorial optimization.

Céline Gicquel received a Ph.D. in Operations Research in 2008 from Ecole Centrale Paris and got an accredition to supervise research from the University Paris Saclay in 2021. She is currently an Associate Professor in Computer Science at the Université Paris Saclay (France) and a member of the Laboratory for interdisciplinary research in numerical sciences (LISN). Her research interests include industrial production planning, facility location, combinatorial optimization and stochastic programming. She is the author of 21 international journal papers.