Kamila M Yunusova

Denis Konstantinov

Alexei D Chepelianskii

Helene Bouchiat

Coupling between Rydberg states and

published or not. The documents may come   L'archive ouverte pluridisciplinaire

The realization of high purity two dimensional electron systems (2DES) has lead to the discovery of fundamental new states in condensed matter physics like integer and fractional quantum Hall effects [1][2][START_REF] Jain | Composite Fermions[END_REF][START_REF] Novoselov | [END_REF][5], as well as to the more recent discovery of two dimensional topological insulators [6]. Electrons on liquid helium were one of the first historical realisations of the 2DES [7][8][9][10]. This system is formed due to the attractive interaction between electrons and their image charge inside liquid helium, it achieves an exceptional purity and still gives the best known electronic mobilities for a 2DES [11]. Electrons on helium enabled the first observation of Wigner crystallization [12,13], edge magnetoplasmons [14] and other exciting many-electron phenomena [15][16][17][18][19]. Considerable efforts were also devoted to study the interaction between electrons on helium and millimeter-wave photons aiming for applications in quantum computing [20,21]. This research direction recently revealed a rich nonequilibrium physics, showing microwave induced oscillations (MIRO) [22][23][24][25][26], zero-resistance states [27][28][29] and incompressible electronic behavior [30][31][32] under excitation by millimeter-wave photons. In the present Letter we show that electrons on helium also allow to realize a model system for an atom interacting with an oscillator (cavity) and to explore its physical properties directly controlling their coupling with a parallel magnetic field. Such systems have been embodied for example in atomic physics [33] and quantum optics as well as with superconducting circuits [34,35]. We demonstrate that, for weakly coupled electrons, the quantum electro-dynamics (QED) Hamiltonian reproduces quantitatively the spectroscopic properties of our system. This opens a doorway to study quantum phenomena in an ensemble of interacting atoms in a cavity systems by tuning the strength of electron-electron interactions.

Before presenting our experimental results, we describe the derivation of the QED Hamiltonian for electrons on helium. The electric field of an electron polarizes the liq-uid helium around it and creates an image charge that attracts it towards the helium surface; a steep electronvolt high energy barrier prevents it from penetrating inside the liquid helium. The interaction with the image charge gives rise to a one-dimensional Coulomb potential which leads to the quantization of the vertical motion and to the formation of a Rydberg series of bound states for a one-dimensional hydrogen-like atom. This series will play the role of the atomic degree of freedom in our QED model. A pressing perpendicular static electric field E ⊥ is also present in the experiments, it allows to shift the Rydberg levels through linear Stark effect [36]. The spectroscopic positions of the Rydberg states is well described by a one-dimensional Schrödinger equation for vertical motion:

H a = - h2 2m ∂ 2 ∂z 2 + V a (z) = α ε α |α α| (1)
where we introduced z the vertical distance of the electrons to the helium surface, the eigenstates for the vertical motion |α and their eigenenergies α . Above the helium surface, for z > 0, the confinement potential V a (z) is the sum of the interaction with the image charge and with the perpendicular electric field

V 0 a (z) = -Λ/z -eE ⊥ z with Λ = e 2 (ε-1)
16π(ε+1) and where ε is liquid helium's dielectric constant. On the energy scale of the bound states (∼ 7 K) we can set V a (z) = ∞ inside liquid helium for z < 0. We introduced a subscript V 0 a (z) to the potential since we will show later that V a (z) is renormalized when an in plane magnetic field is present. For usual pressing electric fields E ⊥ ∼ 2 V mm -1 the main contribution to the confinement potential for the lowest eigenstates comes from the interaction with the image charge.

In addition to their vertical motion, electrons on helium move horizontally as free particles -electrons with their bare electronic mass m. A perpendicular magnetic field applied to 2DES induces the Landau quantization of horizontal motion and the formation of equidistant Landau levels, the Hamiltonian for horizontal motion (up to a constant) then becomes H l = hω c â+ â where ω c = eB z /m is the cyclotron frequency. This term has the same form as the Hamiltonian of a resonant cavity in QED, the Landau level index then plays the role of the number of light quanta in the cavity. With only a perpendicular magnetic field and in the limit of weak electron-electron interaction the Landau levels and Rydberg states are not coupled. A tunable coupling can be introduced by applying an in plane magnetic field [37,38]. Indeed a magnetic field applied in the y direction will tend to turn a vertical velocity towards the x direction due to cyclotron motion along the y axis induced by the parallel field. This coupling has been investigated in double quantum wells in a regime with many occupied Landau levels [39][40][41], in this letter we focus instead on the limit where only the lowest Landau level is occupied. The quantitative form of the interaction induced by the in plane field can be obtained as follows. We write the total Hamiltonian H = (p-eA) 2 2m + V 0 a (z), using the Landau Gauge A = B y ze x + B z xe y where the vector potential doesn't have any component along the z-axis motion. This Hamiltonian can be expanded in the powers of B y : to the lowest order we have Ĥ = Ĥa + Ĥl , the first order in B y introduces an atom-cavity interaction term Ĥc = -eB y z px /m. Writing Ĥc in terms of the Landau level creation and annihilation operators we obtain the following expression:

Ĥ = hω c â+ â + α ε α |α α| + hω y √ 2 â+ + â ẑ B (2)
In this equation we introduced the cyclotron frequency along the in plane field ω y = eB y /m as well as the magnetic length for the perpendicular field B = h/(mω c ), the notation ẑ stands for the matrix elements of the z operator on the vertical eigenstates |α , it plays here the role of the dipole moment operator in quantum electrodynamics. The QED Hamiltonian Eq. ( 2) appears in models where a photon mode/harmonic oscillator (Landau levels in our experiment) is coupled to an atom/qubit provided by Rydberg states. The strength of the interaction, which would be the vacuum Rabi splitting in atomic physics, is directly controllable and proportional to B y allowing in principle couplings of arbitrary strength. This Hamiltonian may seem valid to first order in B y , however the second order diamagnetic term mω 2 y z 2 /2 only renormalizes the vertical confinement potential V a (z) = V 0 a (z)+mω 2 y z 2 /2. Thus Eq. ( 2) remains valid for arbitrary interaction strength keeping in mind that the in-plane magnetic field then not only controls the coupling strength between the atom and Landau levels but also changes the atom energies α and the dipole momentum matrix ẑ, which can still be obtained easily by solving the one dimensional Schrödinger equation Eq. ( 1) in the modified confinement potential. To check if QED Hamiltonian quantitatively describes the energy of the transitions between Rydberg states, we perform Stark effect spectroscopy. Electrons on helium form a static dipole with their image charge and the Stark effect due to the perpendicular electric field E ⊥ leads to a linear displacement of the Rydberg energy levels which can bring these atomic transitions in resonance with the external millimeter-wave excitation. At resonance a change of resistivity occurs due to MIRO allowing to detect the position of the energy levels. Our experimental setup consists of a cavity with Corbino electrodes, a layout is shown in Fig. 1. This cavity is halffilled with liquid helium by condensing helium vapor and monitoring the capacitance between top and bottom electrodes. Electrons are then deposited by thermal emission from a heated filament and are trapped on the surface by the pressing electric field E ⊥ . They form a 2DES which behaves, for in plane transport, as an effective resistance R placed between two contacts with capacitance C. This resistance can then be determined by measuring the admittance of the cell Y between the two inner Corbino contacts from the top electrodes at frequencies comparable with the RC relaxation time (we used 1137 Hz). To extract the MW dependent (MIRO) admittance δY , MW power is modulated at a frequency of 17 Hz and a double demodulation technique is used. Real and imaginary parts of δY give very similar lineshapes. In our measurements the electron gas density was n e 1.5 × 10 7 cm -2 with a total number of 4 × 10 7 electrons trapped in the cloud.

The conversion between Stark shifts and transition energies is obtained from a calibration experiment where we excite the electrons with photons at different energies and measure the electric field at resonance (see Fig. 2). For weak parallel magnetic fields the transition from the ground |g to the first excited Rydberg state |e manifests as a resonance of the microwave induced change in admittance as function of E ⊥ . The resonance position at energy hν 0 = eg shifts linearly with E ⊥ , the slope can be obtained from the Schrödinger equation Eq. ( 1) with small deviations due to uncertainties on geometrical parameters. This slope is almost independent of B y (see Fig. 2.a), indeed for B y ≤ 1 Tesla the coupling term of the QED Hamiltonian

hωy √ 2 e|ẑ|g B < ∼ 10
GHz is small compared to hν 0 140 GHz and does not change the vertical dipole moment significantly. While the slope as function of E ⊥ remains unchanged, an overall energy shift δ is visible. It appears due to the coupling between Rydberg states and Landau levels at finite B y . In the following, we present a careful experimental investigation of the coupling induced energy shift and show that it can be understood quantitatively from the QED Hamiltonian. To study the evolution of the |g → |e transition with B y , we take advantage of the linear dependence of the energy shifts on E ⊥ , which enables us to fix the excitation frequency to f = 139 GHz and change only E ⊥ . We define δ as the detuning induced by the Shark shift due to the deviation of E ⊥ from its resonant value at the excitation frequency f for B y = 0; δ is thus the difference between E ⊥ and its value at resonance 2 Vmm -1 times (minus) the slope measured in Fig. 2.a. All the collected data is then transformed into a map where the change in admittance is plotted as function of B y and of the detuning δ. The maps that we obtained for B z = 1.3, 1.05 are displayed on Fig. 3(a,b). When increasing B y we can resolve two transitions ∆ 0 and ∆ 1 . The energy of the more intense transition called ∆ 0 increases with B y , quadratically at weak B y with a cross-over to a more linear dependence at the highest fields. In addition, to this main transition a weaker transition ∆ 1 splits off from the main transition as B y becomes stronger (as will be shown later ∆ 1 is visible only at highest microwave excitation power which was here fixed to its maximal value of 10 mW). This ∆ 1 transition becomes as visible as the main transition at lower B z (Fig. 3(c,d) for B z = 0.85, 0.73 T) giving two mutually inverse curves with a characteristic "butterfly" pattern. The coupling strength dependence of ∆ 0 is almost the same for all B z in our data set, on the contrary for ∆ 1 the slope of the transition line as function of the coupling strength increases significantly with B z .

The splitting of the Rydberg transition can be understood from the energy level diagram in Fig. 3(e) which shows how the energy levels from the QED Hamiltonian evolve with the coupling strength. For each atomic state |α , the manifold of dressed states consists of a ladder of Landau levels |α, m with MW exciting transitions that conserve the Landau level number m (between states with the same number of photons in the cavity). Without parallel magnetic field the energy of the |g, m → |e, m transition does not dependent on m. The coupling lifts this degeneracy, making transitions associated to different Landau levels spectroscopically distinguishable. In the special case of the m = 0 transition |g, 0 → |e, 0 with energy ∆ 0 the renormalization of the transition energy is due to an interaction with the lowest Landau level. It can be seen as an effective Lamb shift in analogy with atomic physics, circuit QED [START_REF] Milonni | The quantum vacuum: an introduction to quantum electrodynamics[END_REF][START_REF] Wallraff | [END_REF] and physics of electrons coupled to phonons or ripplons [44,45]. A similar renormalization occurs for all the transitions |g, m → |e, m , and simulations are thus needed to identify the observed spectroscopic lines as one of the transitions ∆ m .

The QED Hamiltonian gives a quantitative prediction on the renormalization of the transition energies ∆ m , we emphasize that all the parameters appearing in the model involve E ⊥ , the applied magnetic field, the liquid Helium dielectric constant and fundamental constants, there are thus no fitting parameters. The values of ∆ m can be obtained from the numerical diagonalization of Eq. ( 2), to obtain accurate values we had to include Rydberg states and Landau levels at an energy scale higher than hν 0 from |g, 0 >, in the simulations shown here we used a basis set of 100 Landau levels and 20 Rydberg states. Results of our simulations for transitions ∆ 0,1 are overlaid on top of the experimental data. We see that they reproduce accurately both upper and lower "butterfly wings", including the striking increase of ∆ 1 (B y ) with B z which contrasts with the ∆ 0 transition that is almost B z independent. The transitions ∆ m between states |g, m → |e, m can only be observed if the initial state |g, m is pop-ulated. At the experiment temperature T = 0.3 K only the ground state |g, 0 is populated in equilibrium (the thermal population of |g, 1 at T=0.3 K and B z = 1 T is only 1%). As a consequence, the transitions ∆ m with m ≥ 1 require an external excitation to become visible. In Fig. 4(a) we show that indeed these transitions physically appear only when the MW power is high enough, as opposed to the ∆ 0 transition which is present even at low MW power. Two possible mechanisms to populate the |g, 1 level may be taken into consideration. The first assumes in plane component of the MW electric field populating a |g, 1 level non-resonantly from the initial |g, 0 level. The second one is illustrated with dashed lines on the Fig. 4(b). When the energy of the transitions ∆ 0,1 are sufficiently close (at low B y ) a ∆ 1 energy MW photon can also excite a transition into the |e, 0 state, scattering can then transfer some population into a nearby |g, m level, leading (after relaxation) to a finite population in the |g, 1 state which makes the transition ∆ 1 visible. In Fig. 3 we see that as B y increases the ∆ 1 transition disappears faster than the ∆ 0 transition. This observation can be understood within the population mechanism for |g, 1 shown in Fig. (4b). Indeed at larger B y the ener-gies of the ∆ 0,1 transitions become different and the MW excitation at the ∆ 1 frequency can no longer excite the ∆ 0 transition which populates the |g, 1 level. The state |g, 1 then remains empty leading to the disappearance of the ∆ 1 wings.

In conclusion we have shown that 2DES on liquid helium can form a prototype quantum system of an atom coupled to an oscillator with the coupling directly controllable by the parallel magnetic field. Our spectroscopic results compare very accurately to theoretical predictions with no adjustable parameters. Control of the population transfer between dressed states could enable tunable mm-wave lasers and future experiments at high electron densities may reveal a rich quantum many body physics.

FIG. 1 .

 1 FIG. 1. Schematic diagram of the experimental cell. The top electrodes 1 and 2 are DC grounded and are used for the AC measurements. A positive DC voltage V d is applied to 4 and 5 confining the electrons into the center of the cell and fixing E ⊥ = V d /h. The electrodes 3 and 6 are used as a guard with negative potential. To ensure an homogeneous E ⊥ we fixed V6 -V3 to V d (and V d -V6 = 6V). The admittance Y of the cell is obtained by applying a 10 mV AC voltage at 1137 Hz to the segmented electrode 2 and measuring the induced pickup voltage on 1 with a lock-in amplifier. It depends on the inplane conductivity of the electrons under magnetic field as obtained from Corbino measurements with Ohmic contacts on conventional 2DES. MW power was sent into the cell through a waveguide and 7 is the filament (e -source).

FIG. 2 .

 2 FIG. 2. The panel a shows the shift of the |g → |e transition energy due to Stark effect at Bz = 0.73 T for By = 0 and By = 0.25 T. The slope is the same for both lines but a small energy shift δ is observed. The panel b illustrates the Stark shift as seen from raw δY data at By = 0. The dashed line displays the shift of the resonance with By for f = 141.6 GHz, δ can also be deduced from the value of this Stark shift.

δ 1 FIG. 3 .FIG. 4 .

 134 FIG. 3. Evolution of the |g → |e resonance with By showing that this resonance splits into two branches. The colour maps represent δY as function of By and detuning δ for Bz = 1.3 (a), 1.05 (b), 0.85 (c) and 0.73 T (d). Red and black curves give the QED Hamiltonian predictions for the ∆0 and ∆1 transitions drawn in panel (e) between states |n, m where the first quantum number gives the atomic state and m is Landau level number (number of photons in the cavity). The calculated evolution of individual levels (rescaled for visibility) with By up to 1 T is shown by coloured lines.

We acknowledge discussions with K. Kono, D.L. Shepelyansky who stimulated these experiments as well as support from N. Vernier, V.P. Dvornichenko, Okinawa Institute of Science and Technology (OIST) Graduate University and ANR SPINEX.


Optimizing Make-To-Stock policies through a robust lot-sizing model

Keywords: Lot-sizing, Make-To-Stock, Robust optimization, Mixed-integer linear programming

In this paper we consider a practical lot-sizing problem faced by an industrial company. The company plans the production for a set of products following a Make-To-Order policy. When the productive capacity is not fully used, the remaining capacity is devoted to the production of those products whose orders are typically quite below the established minimum production level. For these products the company follows a Make-To-Stock (MTS) policy since part of the production is to fulfill future estimated orders. This yields a particular lotsizing problem aiming to decide which products should be produced and the corresponding batch sizes. These lot-sizing problems typically face uncertain demands, which we address here through the lens of robust optimization.

First we provide a mixed integer formulation assuming the future demands are deterministic and we tighten the model with valid inequalities. Then, in order to account for uncertainty of the demands, we propose a robust approach where demands are assumed to belong to given intervals and the number of deviations to the nominal estimated value is limited. As the number of products can be large and some instances may not be solved to optimality, we propose two heuristics. Computational tests are conducted on a set of instances generated from real data provided by our industrial partner. The heuristics proposed are fast and provide good quality solutions for the tested instances. Moreover, since they are based on the mathematical model and use simple strategies to reduce the instances size, these heuristics could be extended to solve other multi-item Agostinho Agra ; Michael Poss ; Micael Santos lot-sizing problems where demands are uncertain.

Introduction

In this paper we consider a practical problem occurring in an aluminium extrusion industrial company. The company produces two main families of products: a family of products representing the main production activity of the company where a Make-To-Order (MTO) policy is followed (MTO family), and a family of products whose orders are typically quite below the established minimum production level. For this family, the company follows a Make-To-Stock (MTS) policy (MTS family). The production planning procedure for the MTO family is well established. However for the MTS family, as the orders are below the minimum production level, the company must find a solution between the two extreme cases: wait for new orders of the same product until the minimum production level is attained, or produce at least at the minimum production level of that item to satisfy the pending orders and store the leftovers in inventory. Both alternatives have their pros and cons. The first alternative has the advantage of avoiding stocks. On the other hand, the backlogging of demand orders may lead to intangible losses. Conversely, the second alternative has the advantage of a ready satisfaction of customer needs but generates high holding costs.

Currently, the company gives priority to the MTO family by planing its production first, and when extra production capacity is available, then it solves a lot-sizing problem to decide which products from the MTS family should be produced and defining the corresponding lot-sizes. This particular lot-sizing problem takes into account not only the pending orders of each product but also future ones, as the excess quantity produced will remain in stock until new orders are received. Therefore, it is necessary to estimate those future client orders. The uncertainty related to forecasting such future demands represents a risk for the planners since the inventory costs will depend greatly on such unknown demands. For industries where holding costs are high (as in the case of our industrial partner) it is desirable to derive robust solutions that take into account possible future deviations from the estimated demand values.

Here we address this lot-sizing problem defined for the MTS family of products, using the available production capacity. We consider both the deterministic and the robust cases where demands are assumed to belong to an uncertainty set and we look for the production plan that optimizes the worst-case scenario.

For the production of the MTS family, we produce at most one batch of each product, hence, we allow at most one set-up. Therefore this particular lot-sizing problem is denoted by LS1S (Lot-Sizing with 1 Set-up). The robust problem is denoted by RLS1S.

Multi-product lot-sizing problems have been receiving a great attention, for recent publications, see e.g. [START_REF] Cunha | On multiitem economic lot-sizing with remanufacturing and uncapacitated production[END_REF][START_REF] Macedo | Hybrid manufacturing and remanufacturing lot-sizing problem with stochastic demand, return, and setup costs[END_REF][START_REF] Sifaleras | Variable neighborhood descent heuristic for solving reverse logistics multi-item dynamic lot-sizing problems[END_REF]. Frequently, due to the variety of products and their demand patterns, the companies follow different production polices for the different products. In some cases, different policies can even be considered for the same product (see [START_REF] Zhang | Dynamic pooling of make-to-stock and make-to-order operations[END_REF]) in order to satisfy the different demand streams. The decision between the MTO and the MTS policies was investigated by [START_REF] Zaerpour | Maketo-order or make-to-stock decision by a novel hybrid approach[END_REF] and Altendorfer and Minner (2014). For an overview on comparison of such approaches see [START_REF] Olhager | The impact of manufacturing and supply chain improvement initiatives: A survey comparing make-to-order and make-tostock firms[END_REF]. However, both MTO and MTS producing processes may share common resources forcing the production planners to coordinate the MTO and MTS policies [START_REF] Rafiei | Capacity coordination in hybrid make-tostock/make-to-order production environments[END_REF]. Examples of problems combining MTO-MTS policies can be found in different industries, such as food production systems [START_REF] Soman | Combined make-to-order and make-to-stock in a food production system[END_REF] and steel plants [START_REF] Zhang | Multi-level inventory matching and order planning under the hybrid make-to-order/make-to-stock production environment for steel plants via particle swarm optimization[END_REF].

Several approaches have been proposed, mostly from last decade, regarding the integration of MTS and MTO policies. Beemsterboer et al. (2016) study the benefits of not prioritizing policies within a hybrid planning MTO-MTS approach. In (Beemsterboer et al., 2017a), the authors analyse the benefits of considering flexible lot sizing policies in a hybrid MTO-MTS approach for a two-product system. In (Beemsterboer et al., 2017b), the authors propose four methods of integrating make-to-stock items in the control of a job shop, which they evaluate using discrete event simulation. [START_REF] Kaminsky | Combined make-to-order/make-to-stock supply chains[END_REF] propose heuristics for a multi-item problem where the manufacturer and the supplier have to decide which items to produce to stock and which to produce to order. [START_REF] Kalantari | A decision support system for order acceptance/rejection in hybrid mts/mto production systems[END_REF] present a decision support system for order acceptance/rejection in a hybrid MTO-MTS production environment. [START_REF] Perona | Combining make-to-order and maketo-stock inventory policies: An empirical application to a manufacturing SME[END_REF] develop a decision-making approach to support inventory management decisions in a MTO-MTS environment for small and medium sized enterprises. [START_REF] Renna | Production control policies for a multistage serial system under mto-mts production environment[END_REF] considers a multistage manufacturing serial system, where a production control strategy is performed to release MTO and MTS orders. [START_REF] Rafiei | Novel bi-level hierarchical production planning in hybrid mts/mto production contexts[END_REF] propose a hierarchical production planning approach for a hybrid MTO-MTS system that includes both mid-term and short-term production planning levels. [START_REF] Rafiei | Multi-site production planning in hybrid make-to-stock/make-to-order production environment[END_REF] propose a genetic algorithm for a multi-site production planning of a hybrid MTO-MTS manufacturing system.

The MTS planning carries the risk that the forecasted orders may not materialize. Such risk has been identified before, see [START_REF] Tang | Identifying risk issues and research advancements in supply chain risk management[END_REF]. When it is possible, delaying product differentiation can be an interesting intermediate solution [START_REF] Gupta | Make-to-order, make-to-stock, or delay product differentiation? A common framework for modeling and analysis[END_REF]), but that is not possible in most practical cases as the one faced by our industrial partner. For those cases, handling with uncertainty is of main relevance on MTS environments. To the best of our knowledge only [START_REF] Khakdaman | Tactical production planning in a hybrid make-to-stock-make-to-order environment under supply, process and demand uncertainties: a robust optimisation model[END_REF] applied a robust multi-objective approach based on a set of scenarios to a hybrid MTO-MTS problem where uncertainty is considered in suppliers, processes and customers.

The problem considered in this paper occurs as a subproblem of a hybrid MTO-MTS manufacture system where a hierarchic approach is followed and priority is given to MTO. The problem focuses on solving the MTS planning considering the remaining manufacturing capacity. From its nature, the MTS subproblem considers medium/long-term horizons where demand uncertainty plays a crucial role when defining lot-sizings.

A large number of publications has been devoted to the study of robust lot-sizing problems with demand uncertainty. One of the first papers on the topic is [START_REF] Bertsimas | A robust optimization approach to inventory theory[END_REF], which proposes a simple conservative approximation of the robust constraints and studies the structure of the optimal policies. In parallel to that work, another paper introduced affine decision rules (Ben-Tal et al., 2004), having the advantage of better approximating the robust constraints. The theoretical strength of affine decision rules has been studied in subsequent papers, among which [START_REF] Iancu | Supermodularity and affine policies in dynamic robust optimization[END_REF]. More recent works have sought to solve the robust problem exactly, by using decomposition algorithms and dynamically adding constraints to the problem, see (Agra et al., 2016;[START_REF] Bienstock | Computing robust basestock levels[END_REF][START_REF] Gorissen | Robust counterparts of inequalities containing sums of maxima of linear functions[END_REF]. Robust lot-sizing problems and their variants are also addressed in more general papers dealing with multi-stage robust optimization, see [START_REF] Delage | Robust multi-stage decision making[END_REF] for a survey on these problems. More generally, we refer to [START_REF] Peidro | Quantitative models for supply chain planning under uncertainty: A review[END_REF] for a survey on papers dealing with uncertainty on supply chains.

Although motivated by a practical problem, we aim to incorporate the recent robust optimization techniques into this particular lot-sizing problem in order to close the gap between the robust techniques for classical lot-sizing problems and the robust techniques for MTS problems within hybrid MTO-MTS manufacture systems.

The contributions of this paper are more specifically detailed below. We introduce a mathematical model for the deterministic case where future demands are assumed to be known. Our model is different from the the classical ones (see for instance [START_REF] Pochet | Production Planning by Mixed Integer Programming[END_REF]) mainly because we suppose that each product has at most one set-up. A proof that this particular problem is NP-hard is given. The model is tightened with valid inequalities.

We develop a robust mixed integer model where demands are considered uncertain and belong to intervals. The uncertainty set is further constrained by budget constraints that limit the number of possible periods where a demand can deviate from its nominal value preventing the solutions to be too conservative, obtaining the well-known budgeted uncertainty set introduced in [START_REF] Bertsimas | The price of robustness[END_REF]. We approximate the resulting robust constraints using the conservative approach of [START_REF] Bertsimas | A robust optimization approach to inventory theory[END_REF], rather than the computationally demanding affine decision rules from (Ben-Tal et al., 2004) or exact approaches used in (Agra et al., 2016;[START_REF] Bienstock | Computing robust basestock levels[END_REF].

Since the problem is NP-hard, and we aim to develop approaches that can be used both with commercial and non-commercial (slower but free) solvers, we propose two heuristics. The first heuristic, called Elite Heuristic, is based on a pre-selection of a set of candidate products. The problem is solved for that restricted set of products using a mixed integer linear programming solver based on the strengthened formulation. The heuristic incorporates the practical rules used by the company to choose the products to produce. The second heuristic, denoted as the Tournament Heuristic, runs in several iterations. At each iteration, the set of candidate products is partitioned into smaller subsets and the problem is solved optimally for each subset. Only the selected products of each subset are considered in the next iteration. The process is repeated until a final subset of products is solved or a number of iterations is attained.

To test the deterministic and robust formulations and the matheuristics we use the non-commercial solver Cbc from Coin-OR (2016), which is referred to as one of the fastest solvers among the non-commercial ones [START_REF] Meindl | Analysis of commercial and free and open source solvers for linear optimization problems[END_REF]. The test set was built from the real data provided by our industrial partner.

As the proposed heuristics use simple strategies to reduce the number of items and, consequently, the size of the instances, such heuristics can be easily adapted to other multi-item lot-sizing problems. It suffices to adapt the mathematical model to the particularities of the other problems. We also show, that in order to derive solutions that take into account future demands variations, robust strategies could be embedded into the mathematical model, and therefore into the heuristics, but of course such strategies would need further computational testing in other cases and contexts.

The outline of the paper is as follows. In Section 2 we introduce a mixedinteger formulation to model the practical LS1S problem assuming the demands are deterministic. The formulation is enhanced and a proof of NP-hardness is given. Then, in Section 3, we derive the robust model for the case where demands belong to an uncertainty set. In Section 4 we present the two heuristics.

Computational experiments are reported in Section 5. Final conclusions are given in Section 6.

Formulation

In this section we introduce a mixed integer formulation for the LS1S problem. The formulation is presented in a generic format in order to establish connections to related models and existent literature. Let m denote the number of items considered and n denote the number of time periods of the planning horizon, and define the sets M = {1, . . . , m} and N = {1, . . . , n}. We split the time horizon into two sub horizons N 1 = {1, . . . , n 1 } and N 2 = {n 1 + 1, . . . , n}.

The first horizon is for production planning (where the extra production capacity is available) while the second horizon is considered for the inventory management aspects. The demand of item i ∈ M, in time period t ∈ N 1 , denoted by

d e
it , is assumed to be known and, in our case, corresponds to pending orders.

The demand for item i ∈ M in time period t ∈ N 2 is denoted by d p it and is forecasted. If item i is produced, then the amount produced must be comprised between Q i and Q i . For each item i, parameters p i and q i represent the unit production cost and the fixed production cost, respectively. Such parameters may be negative if we allow them to incorporate, for instance, the selling price.

For each product i ∈ M , and for each time period t ∈ N , parameters h it and g it are assumed to be nonnegative and model the unit inventory cost and the unit backlogging cost, respectively. Finally, S is the maximum inventory capacity.

To formulate the problem, we define the following decision variables: x it is the production of item i in period t; z it indicates whether there is production of item i in period t; y i is the set-up variable which is 1 if x it > 0 for some t ∈ N 1 , and 0 otherwise; s it is the inventory of item i at the end of time period t, and r it is the backlogged demand of item i at the end of period t. The mixed integer programming formulation for LS1S is described below.

min i∈M t∈N1 p i x it + i∈M t∈N h it s it + i∈M t∈N g it r it + i∈M q i y i (1) s.t. x it + s i,t-1 + r it = d e it + s it + r i,t-1 ,∀i ∈ M, t ∈ N 1 , (2)
s i,t-1 + r it = d p it + s it + r i,t-1 , ∀i ∈ M, t ∈ N 2 , ( 3
) i∈M s it ≤ S, ∀t ∈ N, (4)
Q i z it ≤ x it ≤ Q i z it , ∀i ∈ M, t ∈ N 1 , (5)
y i = t∈N1 z it , ∀i ∈ M, (6)
z it ∈ {0, 1}, ∀i ∈ M, t ∈ N 1 , (7)
y i ∈ {0, 1}, ∀i ∈ M, (8)
s it , r it ≥ 0, ∀i ∈ M, t ∈ N, (9)
r i0 = s i0 = 0, ∀i ∈ M, ( 10
) i∈M x it = V t , ∀t ∈ N 1 . (11)
The objective function (1) aims to minimize the sum of the production costs ( i∈M t∈N1 p i x it ), the inventory cost ( i∈M t∈N h it s it ), the backlogging cost ( i∈M t∈N g it r it ), and the fixed production cost ( i∈M q i y i ) over the planning horizon.

Constraints (2) are the inventory balanced constraints written for each item and each production period t ∈ N 1 , while constraints (3) are the inventory balanced constraints for time periods t ∈ N 2 . Constraints (4) impose an upper bound on the stock level. Constraints (5) are the variable lower and upper bound constraints. They impose a lower and an upper bound on the quantity produced of each product at each period and link the corresponding continuous variables to the set-up variables. Equations ( 6) establish the number of setups for each product. Together with (8) they ensure that at most one setup can occur. Constraints ( 7) and ( 8) define the set-up variables as binary.

Constraints (9) ensure non-negativity of the inventory and backlog variables.

Constraint (11) represents the additional constraints related to the available production capacity.

Summing up the equations ( 2) from 1 to n 1 , one obtains

t∈N1 x it + s i0 -r i0 = t∈N1 d e it + s in1 -r in1 , ∀i ∈ M (12)
Denoting t∈N1 d e it by D e i and using s i0 = r i0 = 0, then ( 12) can be written as

t∈N1 x it = D e i + s in1 -r in1 , ∀i ∈ M. (13)
Similarly, summing up equations (3) from n 1 + 1 to ∈ {n 1 + 1, . . . , n}, we obtain

s in1 -r in1 = t=n1+1 d p it + s i -r i , ∀i ∈ M, ∈ {n 1 + 1, . . . , n}, (14)
Using (13) to eliminate s in1 and r in1 , then

s i -r i = t∈N1 x it -D e i - t=n1+1 d p it , ∀i ∈ M, ∈ {n 1 + 1, . . . , n}. ( 15
)
Since s i and r i are nonnegative, we obtain

s it ≥ t∈N1 x it -D e i - t =n1+1 d p i , ∀i ∈ M, t ∈ N 2 , (16)
r it ≥ - t∈N1 x it + D e i + t =n1+1 d p i , ∀i ∈ M, t ∈ N 2 . ( 17
)
As the holding and backlogging costs are assumed nonnegative, we may replace constraints (3) by ( 16) and ( 17). The resulting deterministic model for LS1S, defined by (1), ( 2), ( 4) -( 11), ( 16), (17), will be denoted by DLS1S.

When N 2 = ∅, we obtain a classical lotsizing model with the additional con-straint that at most one setup is allowed for each item. In a MTS environment the stocks may last for a longer time horizon than the planning one. Under deterministic assumptions the second time horizon can be easily dropped since the holding cost in the last time period can be easily computed. However, in a robust setting it is desirable to consider the second time horizon explicitly since it allows to incorporate different seasonal behaviours of demand and allow to account for different magnitudes of deviation to the estimated demands.

Strengthening the formulation

It is well known that the inclusion of valid inequalities can improve the model significantly, see [START_REF] Pochet | Production Planning by Mixed Integer Programming[END_REF] for details. Constraints ( 16) can be strengthened as follows.

s it ≥ t∈N1 x it -D e i y i - t =n1+1 d p i y i , ∀i ∈ M, t ∈ N 2 , (18)
When y i = 1 inequality (18) coincides with ( 16) and when y i = 0 the right-hand side of (18) becomes negative. Thus ( 16) is valid for the set of feasible solutions.

Moreover,

D e i y i + t =n1+1 d p i y i ≤ d e i + t =1 d p i ⇔ t∈N1 x it -D e i y i - t =n1+1 d p i y i ≥ t∈N1 x it -D e i - t =1 d p i
which implies that ( 18) is stronger than ( 16). In fact we can replace ( 16) by ( 18) in the formulation for LS1S.

A family of valid inequalities for the value of backlog variables follows.

r it ≥ D e i (1 -y i ) + t =n1+1 d p i (1 -y i ), ∀i ∈ M, t ∈ N 2 . ( 19
)
If there is no setup for period i, i.e. y i = 0, inequality (19) forces the demand The strengthened deterministic model for LS1S, denoted by SDLS1S is given by ( 1), ( 2), ( 4) -( 11), ( 17), ( 18), ( 19).

D e i + t
In the practical case provided to us by our industrial partner we have the following assumptions: (i) only one production period is considered (n 1 = 1); (ii) only the backlog of the effective demand is penalized (g it = 0, t ∈ N 2 ); and (iii) the minimum production quantity is at least the effective demand

(Q i ≥ d e i ).
The three assumptions imply that constraints (17) can be eliminated. For completeness we give below the resulting strengthened model, after simplification, for the practical case.

min i∈M t∈N1 p i x it + i∈M t∈N h it s it + i∈M g i1 r i1 + i∈M q i y i (20) s.t. s it ≥ t∈N1 x it -D e i y i - t =1 d p i y i , ∀i ∈ M, t ∈ N, (21)
r i1 ≥ D e i (1 -y i ), ∀i ∈ M, (22)
s i1 ≤ S, (23)
Q i y i ≤ t∈N1 x it ≤ Q i y i , ∀i ∈ M, (24)
i∈M t∈N1

x it = V 1 , ∀t ∈ N, (25)
s it , r it ≥ 0, ∀i ∈ M, (26)
y i ∈ {0, 1}, ∀i ∈ M. ( 27
)
Notice that constraints (23) are given just for t = 1 since the stock level will decrease in subsequent periods. We denote model ( 20)-( 27) by PLS1S.

Complexity analysis

The original problem LS1S is NP-hard as it generalizes the classical lot-sizing problem with varying capacities [START_REF] Pochet | Production Planning by Mixed Integer Programming[END_REF]). Here we focus on the particular case of the practical problem PLS1S. We show it is NP-hard for the simplified version with two time periods (implying it is NP-hard for the general case). The proof is done by reducing the partition problem to PLS1S.

Proposition 2.1. The problem defined by PLS1S is NP-hard for the particular case n 1 = 1, n = 2.

Proof: The decision problem, denoted by D-PLS1S, associated with the optimization problem asks whether there is a solution to ( 21)-( 27) whose objective function value given by ( 20) is greater than L.

Next we reduce the partition problem to D-PLS1S. Recall that in the partition problem we are given k positive integers a i , i ∈ K = {1, . . . , k} and wish to determine whether there exists a partition (S, K \ S) of K such that

i∈S a i = i∈K\S a i = i∈K a i /2. For the reduction consider k = m, Q i = Q i = a i and V 1 = i∈K a i /2. Further, for each i ∈ M , we set D e i = min j∈K a j , d p i1 = a i -D e i , p i = q i = g i1 = 0, h it = 0, t ∈ N, L = 0.
As t∈N1 x it = a i y i , any feasible solution to D-PLS1S must satisfy i∈M a i y i = i∈K a i /2. Hence, there is a one to one correspondence between a feasible solution (x * , y * ) of D-PLS1S, with the specified parameters, and a feasible solution of the partition problem, where S = {j ∈ K|y * j = 1}.

Robust model for the demand uncertainty case

Clearly, it is not possible to know precisely the value of d p it for all i and t since these orders have not been made yet. At best, we can rely on historical data to draw a set of plausible values for these demands. A popular approach (see (Bertsimas andSim, 2003, 2004

∈N t |d p i -d p i | dp i ≤ Γ t . (28)
where N t = {n 1 + 1, . . . , t}. Formally, the uncertainty sets obtained with this approach can be written as

D it = d p i : d p i = d p i + δ + i dp i -δ - i dp i , ∈ N t , δ + i , δ - i ∈ [0, 1] , ∈ N t , ∈N t (δ + i + δ - i ) ≤ Γ t ,
for each i ∈ M and t ∈ N , where δ + and δ -are auxiliary vectors that ease the linearization of constraint (28). Parameter Γ t is often denoted as the budget of uncertainty. Taking a small value of Γ t yields a small uncertainty set, while increasing Γ t yields larger and larger uncertainty sets. The two extremes are

Γ t = 0 for which D it is reduced to the singleton {d p i } and Γ t = t for which D it is equal to the box ∈N t [d p i -dp i , d p i + dp i ].
In general, for each t, Γ t is comprised between 0 and t, and Γ t ≤ Γ t+1 . For simplicity we assume Γ t is integer for each t, but all the results derived in this section could be extended to the case where parameters Γ t are fractional.

Next we present a robust model for RLS1S. The model is based on the deterministic formulation DLS1S where constraints ( 16) are replaced by ( 18).

For the remaining models discussed in the previous section (models D-LS1S, SDLS1S and PLS1S) the robust model can be derived in a similar way.

min i∈M t∈N1 p i x it + i∈M t∈N h it s it + i∈M t∈N g it r it + i∈M q i y i (29) s.t. s it ≥ t∈N1 x it -D e i y i - ∈N t d p i y i ,∀i ∈ M, t ∈ N 2 , d p i ∈ D it (30)
r it ≥ - t∈N1 x it + D e i + ∈N t d p i , ∀i ∈ M, t ∈ N 2 , d p i ∈ D it (31)
(2), ( 4) -(11), (17).

Constraints ( 30) and ( 31) are inspired by the approach used in [START_REF] Bertsimas | A robust optimization approach to inventory theory[END_REF]. Their main advantage is to lead to a fairly simple robust counterpart. Namely, we show in the rest of the section that problem (2), ( 4) -( 11), ( 29) -( 31), can be reformulated as a compact MILP that is essentially of the same order of difficulty as the deterministic version LS1S. In contrast, the more advanced methods used in (Agra et al., 2016;Ben-Tal et al., 2004) provide more accurate solutions however at a high computational cost.

Formulation RLS1S contains an infinite number of constraints, yielding a semi-infinite MILP. We show next how to reformulate the problem as a compact MILP by using a well-known technique from robust optimization (e.g. (Ben-Tal and Nemirovski, 1998; [START_REF] Bertsimas | The price of robustness[END_REF][START_REF] Bertsimas | A robust optimization approach to inventory theory[END_REF]).

First, we realize that for (30), positive deviations of d p i are not increasing the cost of the solution since they would only decrease the stock more rapidly; hence, positive deviations can be neglected in (30). Similarly, we can restrict ourselves to negative deviations of d p i in (31) since the positive deviations will never increase the amount of unmet demands. Formally, we can introduce the smaller uncertainty sets

D + it = d p i : d p i = d p i + δ i dp i , ∈ N t , δ i ∈ [0, 1] , ∈ N t , ∈N t δ i ≤ Γ t , D - it = d p i : d p i = d p i -δ i dp i , ∈ N t , δ i ∈ [0, 1] , ∈ N t , ∈N t δ i ≤ Γ t ,
and we replace constraints ( 30) and ( 31) with the equivalent constraints

s it ≥ t∈N1 x it -D e i y i - ∈N t d p i y i , ∀i ∈ M, t ∈ N 2 , d p i ∈ D - it , (32)
r it ≥ - t∈N1 x it + D e i + ∈N t d p i , ∀i ∈ M, t ∈ N 2 , d p i ∈ D + it . (33)
Next we see that the infinite numbers of constraints ( 32) and ( 33) can be substituted by the following non-linear constraints

s it ≥ t∈N1 x it -D e i y i -min d p i ∈D - it ∈N t d p i y i , ∀i ∈ M, t ∈ N 2 , (34)
r it ≥ - t∈N1 x it + D e i + max d p i ∈D + it ∈N t d p i , ∀i ∈ M, t ∈ N 2 . ( 35
)
Next we follow the classical dualization approach introduced in (Ben-Tal and Nemirovski, 1998) and used to lot-sizing problems in [START_REF] Bertsimas | A robust optimization approach to inventory theory[END_REF]. Let i ∈ M and t ∈ N be fixed and let us focus on constraint (34).

The inner minimization problem over variables d p i in (34) can be replaced by a minimization problem over variables δ i . Moving the term y i ∈N t d p i outside of the minimization because it does not involve the variable δ i , and changing the minimization problem by a maximization problem, we obtain

y i ∈N t d p i -max ∈N t δ i dp i y i s.t. ∈N t δ i ≤ Γ t , (36)
0 ≤ δ i ≤ 1, ∀ ∈ N t . ( 37
)
Let us denote the dual variables of constraints ( 36) and ( 37) as z - it and u - it , respectively. Since constraints ( 36) and ( 37) define a bounded and non-empty polytope, we can apply strong linear programming duality to replace the maximization problem by its dual:

y i ∈N t d p i -min Γ t z - it + ∈N t u t- i s.t. z - it + u t- i ≥ dp i y i , ∀ ∈ N t , (38)
z - it , u t- i ≥ 0 ∀ ∈ N t . ( 39
)
Plugging the above minimization problem into the original constraint (34) for the fixed i and t, we obtain

s it ≥ t∈N1 x it -D e i y i -y i ∈N t d p i +       min Γ t z - it + ∈N t u t- i s.t. z - it + u t- i ≥ dp i y i , ∀ ∈ N t z - it , u t- i ≥ 0 ∀ ∈ N t       . ( 40
)
Finally, notice that (40) is feasible if and only if there exist vectors z - it ≥ 0 and u t- i ≥ 0 that satisfy the dual constraints (38) and such that

s it ≥ t∈N1 x it -D e i y i -y i ∈N t d p i + Γ t z - it + ∈N t u t- i . (41)
Therefore, the robust constraint ( 34) is equivalent to the constraints ( 38), (39), and (41).

Introducing dual variables z + it and u + it to handle with constraints (35), we can reformulate RLS1S as the following compact MILP, denoted by RLS1S.

min i∈M t∈N1 p i x it + i∈M t∈N h it s it + i∈M t∈N g it r it + i∈M q i y i s.t. s it ≥ t∈N1 x it -D e i y i -y i ∈N t d p i + Γ t z - it + ∈N t u t- i ,∀i ∈ M, t ∈ N 2 z - it + u t- i ≥ dp i y i , ∀i ∈ M, t ∈ N 2 , ∈ N t r it ≥ - t∈N1 x it + D e i + ∈N t d p i + Γ t z + it + ∈N t u t+ i , ∀i ∈ M, t ∈ N 2 z + it + u t+ i ≥ dp i , ∀i ∈ M, t ∈ N 2 , ∈ N t
(2), ( 4) -( 11), ( 17),

z - it , u t- i , z + it , u t+ i ≥ 0, ∀i ∈ M, t ∈ N 2 , ∈ N t .
The above approach is easy to apply because it amounts to solve a unique problem. Yet, the dimension of the new MILP is larger than the dimension of the original problem.

An alternative approach to the dualization is to compute a priori the minimum and maximum occurring in the right-hand side of inequalities ( 34) and ( 35), respectively. As explained above, these optimization problems can be converted into a maximization problem where the feasible set is defined by ( 36) -( 37). We can observe that the extreme solutions of this set satisfy δ i ∈ {0, 1}.

Thus, the sets of extreme feasible solutions correspond to uniform matroids.

Consequently, these maximization problems can be solved by a greedy algo-rithm that chooses the highest deviations. Hence, the following equalities hold.

min d p i ∈D - it ∈N t d p i y i = ∈Nt d p i - Γt r=1 dp i t r y i max d p i ∈D - it ∈N t d p i = ∈Nt d p i + Γt r=1 dp i t r
where dp

i t r
is the r th largest deviation among the first t deviations dp i1 , . . . , dp it .

Hence, ( 34) and ( 35) are replaced by

s it ≥ t∈N1 x it -D e i y i - ∈Nt d p i - Γt r=1 dp i t r y i , ∀i ∈ M, t ∈ N 2 , ( 42
)
r it ≥ - t∈N1 x it + D e i + ∈Nt d p i + Γt r=1 dp i t r , ∀i ∈ M, t ∈ N 2 . ( 43
)
For the computational results we use this second approach since it is more efficient computationally.

Heuristics

Companies seek for quick approaches to find good solutions to their problems. While compact, formulations SDLS1S and RLS1S can take too much time to be solved to optimality using exact algorithms such as the Branch and Cut implemented in both commercial and open-source optimization software.

The running time is even more relevant in our case since the number of items considered can be quite large, up to few hundreds. Hence, for a practical use, we present two heuristics to solve problems LS1S and RLS1S approximately.

A first one, called Elite Heuristic extends the current practice of the company which is based on the selection of items accordingly to some criteria. Here we select a larger number of items accordingly to the same criteria and solve the models restricted to the selected items. The second heuristic, called Tournament Heuristic, selects iteratively small subsets of items until it reaches a final and small subset.

Elite Heuristic

The first heuristic is to select a specific subset of items and solve the proposed formulation considering this subset of items.

The choice of the subset is as follows. Choose

• m 1 items with the highest known demand values, D e i ;

• m 2 items with the highest values of the known demand plus the forecasted demand for the first time period, D e i + d p i,n1+1 ;

• m 3 items with the highest values of the known demand plus the forecasted demand for the first two time periods, D e i + n1+2 t=n1+1 d p it ;

• m 4 items with the highest values of the known demand plus the forecasted demand for the first three time periods, D e i + n1+3 t=n1+1 d p it ;

Here we consider m 1 = m 2 = m 3 = m 4 .

After selecting and joining these four list of items, the restricted model (the deterministic SDLS1S or the robust RLS1S) is solved. As the number of elite items is small the resulting model can be solved to optimality easily. As example, if we are given a set with 200 items and take m 1 = m 2 = m 3 = m 4 = 5, then five items will be selected from the initial set of 200 items using each one of the four criteria given above. The resulting restricted model with the selected items is solved using a solver.

Notice that the selection criteria do not take into account the costs. If the costs vary from product to product significantly, which is not the case of the instances we consider, then other criteria taking into account the cost structure (production costs, inventory costs, backlogging costs, fixed production cost) should be used. For instance, if the backlogging costs vary significantly, a fifth criteria should be considered which consists in selecting the m 5 items with the highest backlogging cost. Similarly, to the remaining cost components.

Tournament Heuristic

The second heuristic is to decompose the original problem into smaller and easier subproblems which can be solved quickly. The items selected from these subproblems are the input to the subproblems of the next step. The process is repeated until a final and smaller subset of items is obtained. Next we detail this heuristic.

The first step is to divide the set of m items into r 1 subsets with cardinalities m 11 , m 12 , . . . , m 1r1 , respectively, where m 11 + m 12 + . . . + m 1r1 = m. Then solve the restricted problem for each one of the subsets. As the resulting subproblems are simpler they are solved to optimality. The items that are produced in the optimal solution of each one of the subproblems are selected to the next iteration.

Let m 1 denote the number of items selected in the first step. Then split this set into r 2 subsets with cardinalities m 21 , m 22 , . . . , m 2r2 , respectively, where m 21 +m 22 +. . .+m 2r2 = m 1 . Each one of the subproblems is solved to optimality again in order to identify the new selected (produced) items. The process is repeated until a final subset with cardinality less or equal to a threshold is obtained or a maximum number of iterations is attained. The final solution is the optimal solution to the problem restricted to the final set of items. Table 1 outlines Heuristic 2.

In order to exemplify the decomposition process, assume we are given a set with 200 items and consider r 1 = 6. Thus, the set of items is split into six subsets as follows: four subsets with 33 elements (m 11 = m 12 = m 13 = m 14 = 33) and two subsets with 34 elements (m 15 = m 16 = 34). Then, each . . .

m 32 . . . • • • m n-1 m n1 m 3r3 m 2r2 m 1r1
problem corresponding to a different subset is solved to optimality and the items produced in each one of the six problems are selected to the next round.

Notice that the solution to each subproblem is feasible to the original problem. Hence, it is expected that the quality of the solutions obtained will improve in each iteration, culminating in the solution of the final iteration.

The quality of the final solution may depend on the successive partitions of the set of items. The selection of items for each subset can be done randomly and the process can be repeated several times. Here we will not explore such possible improvement.

Numerical experiments

In this section we report the computational tests. The objectives of these numerical experiments are three-fold. First, we assess the difficulty of the deterministic and robust instances, reporting also the gains obtained by using the strengthened models proposed in Section 2.1. Second, we evaluate the objective function values for the robust and deterministic models to test the importance of using robust approaches. Last, we test the efficiency of the two proposed heuristics.

All tests were conducted on a computer Intel(R) Core(TM) i3-3250 CPU, 3.50GHz with 4 cores, using the open-source solver Cbc 2.9 (Coin-or branch and cut) from [START_REF]Computational Infrastructure for Operations Research[END_REF]. The heuristics were implemented in Julia, using the package JuMP [START_REF] Lubin | Computing in operations research using julia[END_REF].

Instances were generated using data provided by our industrial partner for the MTS family of products. Since the number of items considered depends on the pending order quantities (as items with few orders are not considered), the number of items varies weekly and can go up to a few hundred products (the complete MTS family). Also, the production capacity, V 1 , varies weekly as it represents the residual capacity after the production of MTO products have been considered. Based on the information provided by our partner, suggesting to pick-up 80 items and to use a residual capacity around 2000, we generate fifteen instances for the deterministic problem by varying the number of items and the residual capacity. In relation to these two parameters (number of products and capacities), the fifteen instances aim to simulate realistic instances faced by the company currently and in the future. For the number of items we consider five possible values 80,100,150,200,300. Notice that for testing purposes it is not interesting to consider small size and, therefore, easy instances. For the production capacities, V 1 , three values are considered: 1000, 2000, 3000 (since both sides with higher and lower capacities are relevant). The remaining data is taken as follows. The known demands D e i are given by real data and the forecasted demands d p it are obtained from the historical average demand of each item. The following additional parameters are considered (established by our industrial partner, see [START_REF] Santos | Otimização do dimensionamento de lotes na extrusal[END_REF] for details): T is set to 24, Q i = max{250, D e i },

h it = H * (1 + J) t
, where H = 3.52, J = 0.001651, p i = 0.991003, q i = 669.11.

A unit penalty cost of 1 is assumed for each item (g i1 = 1). For the robust settings, d In relation to the production capacity, we can observe that the medium capacity instances seem to be a bit harder than the other ones.

Although not reported in the tables, the average of root gap is 60% for the weak formulation and 5% for the strong formulation. Nevertheless, all unsolved instances have a final integrality gap less than or equal to 3%. For the strong formulation the lower bound for the geometric average running time of the solved instances is 2.94 (it is a lower bound because instances faster than 1 second are set to 1 second).

We also solved these instances using CPLEX 12.7 rather than Cbc 2.9, keeping the time limit of 1800 seconds. Unsurprisingly, the former is much faster than the latter, the weak formulation solving already all but 11 instances (vs 34 for Cbc) while the strong formulation leaves only 3 instances unsolved (vs 17 for Cbc). Fortunately, we will show below that our heuristics, based on Cbc perform very well on our instances, solving nearly all of them to optimality.

Determining the price of robustness

Here we discuss the price for considering robust solutions. As explained in the introduction, estimating the future demands by the historic average values doesn't immunize the solution for possible deviations in the demand values that can lead to higher costs than the estimated ones. For the practical case considered here, the costs that may be underestimated are the inventory costs which are based on the estimated future demands, since for the backlogged demand only the pendent orders are penalized.

In order to report the computational results we define C (i,j) as the cost of the solution obtained for Γ = i assessed when Γ = j. For instance, C (0,2) is the cost obtained for the optimal deterministic solution (obtained for Γ = 0) when facing an uncertainty level of Γ = 2, that is, when we allow the demand values for two time periods to suffer a maximum deviation, either of 20% or 40%, from the estimated nominal values.

In tables 4 and 5 we report for the two deviation levels the value C (i,j) , which is the average of the percentage gaps between the parameters C (i,j) and C (i,i) , defined formally by

C (i,j) = C (i,j) -C (i,i) C (i,i) .
For example, C (0,2) = 1.76 means that the deterministic optimal solutions are, on average, 1.76% more expensive than the optimal solutions obtained for the robust model with Γ = 2 when facing an uncertainty level of Γ = 2. We see from these tables the robust solutions are more expensive, on average, than the deterministic solutions in the deterministic context. For the case where the deviation is allowed to be 40% we can see that protecting a solution for Γ = 3 deviation periods will increase the cost by 3% in relation to the deterministic solution. The deterministic solution is more expensive when Γ ∈ {1, 2} but not when Γ = 3. In the worst case (two deviations), not protecting the deterministic solution will imply an increase of the cost of 3.67%.

We can also observe that robust solutions with Γ = 3 are not interesting.

This justifies why we have not included results for larger values of Γ. Of course, the final choice of which model to use (which value for Γ) depends on the riskaverseness of the decision maker.

Heuristic performances

Finally we test the two proposed heuristics. As discussed above, running a MILP solver for a given time limit acts as a heuristic for those instances that were not solved to optimality. However, from a practical viewpoint, it is more appealing to have a tool which enables the decision maker to derive good solutions very quickly in order to allow him/her to test different parameters before taking a decision.

Table 6 summarizes the results obtained with the Elite heuristic. We can observe that when m increases the optimality gap decreases as expected, being zero for m greater or equal to 8. The running times are very small indicating that this heuristic performs very well for such values of m. This means that selecting at least 32 items and solving the restricted model is fast and gives good quality solutions.

Table 7 reports the results obtained with the Tournament heuristic. This heuristic is run until a final set of m ≤ 50 items is obtained or until a maximum of 6 iterations is attained. In each iteration the original set of items is split into r subsets. The first line gives the value of r. The second line gives the average number of iterations. The third line gives the average optimality gap and the last line gives the average running time in seconds. Although the heuristic is fast and gives solutions that in average have very small optimality gaps, we could not find parameters that lead to optimality in all the tested cases. Overall, the Tournament heuristic was better than the Elite heuristic only for the case where the number of elite elements selected was small ( m = 5), otherwise the Elite heuristic outperformed the Tournament heuristic.

Conclusions

We consider a practical problem faced by a company that plans the production for a set of products following a Make-To-Order policy and uses the remaining production capacity to produce items for which the quantities ordered are small. The problem is concerned with the use of this extra production capacity in a given time period. Namely, deciding which items to produce and the corresponding production level. A penalty is associated with the backlogged demands of the items that are not produced. For the produced items, as the amount produced is in general greater than the pending orders, a holding cost is incurred.

We present a general model and establish the relation between this model and the classical lot-sizing models. We propose several enhancements for the formulations. As the holding costs depend on the value of the estimated future demands, we derive a robust model from the enhanced model where the future demands are considered uncertain and can vary in a given interval centered in the historical average values. The uncertainty set is the classical budget polytope introduced by [START_REF] Bertsimas | The price of robustness[END_REF].

For practical purposes we propose two heuristic schemes that are based on the enhanced models (deterministic and robust). The models and the heuristics are tested using a set of instances generated from real data provided by our industrial partner. The computational experiments show that most instances cannot be solved to optimality within reasonable running time limit. In particular, they show that the Elite heuristic, which selects a small subset of items (elite items), following the company criteria, is quite fast and in general obtains the optimal solution.

It is well-known that robust optimization approaches can generate too conservative solutions. That is, solutions that are good when the worst case scenario occurs, in our case when the orders for future demands are lower than expected, but their quality may not be so good when other demand scenarios are observed.

While this falls beyond the scope of the current work, it could be interesting to compare our robust approach with the optimal solutions of stochastic programming approaches where one also takes the probabilities of the scenarios into account. Similarly, it would be interesting to test the effect of our model and the quality of the heuristics with different cost structures and/or on more complex problems where the set of feasible production plans is restricted by

  varies in {0, 1, 2, 3}. This gives a total using the strengthened formulation instead of the weak one. There is no clear correlation between the difficulty of the instances and the value of Γ parameter.

  For each test we consider m 1 = m 2 = m 3 = m 4 = m. The first line gives the value of m, the second line gives the average of the optimality gaps. The optimality gap is defined as 100 * ((ZH -Z * )/Z * ), where Z * is the optimum value and ZH is the solution cost returned by the heuristic. The third line gives the average number of items selected when combining the four lists, and the last line gives the average running times in seconds.

  =1 d pi to be backlogged. Otherwise, if y i = 1, inequality (19) simply imposes nonnegativity on the backlog variables. Contrary to the previous case,

	(19) cannot replace inequalities (17).

  ;[START_REF] Bertsimas | A robust optimization approach to inventory theory[END_REF])) considers the nominal value d p it and the deviation dp it for each i and t which could, for instance, be the mean value and the variance of the available historical data.

	The approach then supposes that the unknown parameter d p it can take any
	value in the interval [d p it -dp it , d p it + dp

it ] and that, for each item i and each time period t, the number of demands taking an extreme value is bounded by a given parameter Γ t > 0:

Table 1 :

 1 Scheme of the Tournemant Heuristic.

	1 st Step	2 nd Step	3 rd Step	• • •	n th Step
		m 11
		m 12	m 21
			m 22	m 31
	m	. . .	m 1	m 2

Table 4 :

 4 Costs C (i,j) (expressed in %) considering a deviation 0.2.

		0	1	2	3
	0	0	1.19 1.76 0
	1 1.49	0	0.81 0
	2 1.20	0	0	0
	3 1.19	0	0	0

Table 5 :

 5 Costs C (i,j) (expressed in %) considering a deviation 0.4.

		0	1	2	3
	0	0	2.64 3.67 0
	1 2.96	0	0.88 0
	2 2.99 0.33	0	0
	3 2.99 0.33	0	0

Table 6 :

 6 Average statistics for the Elite heuristic.

	Value of m	5	6	7	8	9	10
	Optimality gap	5.6	0.2 0.04	0	0	0
	Number of items selected 8.6	11	12	13.4 15.4	17
	Running time	0.04 0.07 0.09 0.11 0.12 0.13

Table 7 :

 7 Average statistics for the Tournament heuristic.

	Value of r	6	5	4
	Number of iterations	2	2.62 3.33
	Optimality gap	0.35 0.49 0.12
	Running time	0.12 1.81 0.14

of 120 instances for the robust problem. Notice that the deterministic instances (Γ = 0) are considered for both levels of deviations dp i as the assessment of the price of robustness is different for both levels.

Testing formulations and accessing instances difficulty

In the first experiments the instances are solved with a time limit of 1800 seconds. Table 2 gives the number of instances that were not solved within the given time limit using the model PLS1S with the improvements discussed in Section 2.1. Column capacity gives the production capacity divided by 1000, the second column with Γ = 0 is for the deterministic case, columns 3-8 consider the robust case where Γ t = min{t, Γ}. The first three columns assume a maximum deviation of 20% from the nominal value and the next three columns assume a maximum deviation of 40% of that value. The last column gives the total number of instances that were not solved to optimality. Table 3 provides similar information however for the model PLS1S without the improvements (weak model).

We can see that the number of unsolved instances drops from 34 to 17 by additional constraints (e.g. storage capacity, inclusion of set-up times). In particular, the heuristics would need to be adapted to these other contexts and further computational testing would be required.

Aknowledgements

The