A Stochastic Multi-item Lot-sizing Problem

with Bounded Number of Setups

Etienne de Saint Germain

1,2

, Vincent Leclère

and Frédéric Meunier

Argon Consulting, 122 Rue Édouard Vaillant, 92300 Levallois-Perret, France

CERMICS, Cité Descartes, 6-8 Avenue Blaise Pascal, 77455 Champs-sur-Marne, France

Keywords:

Lot-sizing, Stochastic Optimization, Sample Average Approximation, Simulation.

Abstract:

Within a partnership with a consulting company, we address a production problem modeled as a stochastic

multi-item lot-sizing problem with bounded numbers of setups per period and without setup cost. While this

formulation seems to be rather non-standard in the lot-sizing landscape, it is motivated by concrete missions

of the company. Since the deterministic version of the problem is NP-hard and its full stochastic version

clearly intractable, we turn to approximate methods and propose a repeated two-stage stochastic programming

approach to solve it. Using simulations on real-world instances, we show that our method gives better results

than current heuristics used in industry. Moreover, our method provides lower bounds proving the quality of

the approach. Since the computational times are small and the method easy to use, our contribution constitutes

a promising response to the original industrial problem.

1 INTRODUCTION

Fixing the production level for the forthcoming week

is a basic decision to be taken when managing a pro-

duction line. Usually, a demand has to be satisﬁed

at due dates but the limited capacity of the line pre-

vents last minute production. On the other hand, too

early productions may lead to unnecessary high in-

ventory costs. The challenge of this kind of pro-

blems, known as lot-sizing problems in the operatio-

nal research community, consists in ﬁnding a trade-

off between demand satisfaction and holding costs.

This is a well-studied topic, with many variations (de-

terministic/stochastic, single/multi item, etc.). Re-

cent surveys have been proposed: see (Gicquel et al.,

2008, Quadt and Kuhn, 2008) for the determinis-

tic version and (Mula et al., 2006, Aloulou et al.,

2014, Díaz-Madroñero et al., 2014) for the stochastic

version. When several references can be produced on

a same line – the so-called multi-item lot-sizing pro-

blem –, the capacity is often all the more reduced as

the number of distinct references produced over the

current week is high. Indeed, changing a reference in

production stops the line for a moment. This additi-

onal capacity reduction is usually modeled by setup

costs contributing to the total cost.

The present work introduces a stochastic multi-

item lot-sizing problem met by the authors within

a partnership with a consulting company. A non-

standard feature of the problem is that the capacity

reduction due to reference changes is not modeled by

setup costs but instead by an explicit upper bound on

the total number of references that can be produced

over a week. According to the consulting company,

many clients aim at minimizing mainly their inven-

tory costs while keeping the number of distinct refe-

rences produced over each week below some thres-

hold. This is essentially because, contrary to inven-

tory costs, setup costs are hard to quantify and a max-

imal number of possible changes per week is easy

to estimate. To the best of authors’ knowledge, the

problem addressed in the present work is original and

such a bound on the number of distinct references pro-

duced over a week has not been considered by acade-

mics yet, with the notable exception of (Rubaszewski

et al., 2011) but, contrary to our problem, their bound

is an overall bound for the whole horizon and they

still consider setup costs.

We consider two versions of the problem. In a ﬁrst

one, backorder costs are present and the objective is

to minimize the overall inventory costs (sum of hol-

ding and backorder costs). In a second version, there

are only holding costs and there is a service level con-

straint to be satisﬁed over the whole horizon. We pro-

pose for these two versions a method that can be ea-

sily used and maintained in practice. The efﬁciency

de Saint Germain E., Leclère V. and Meunier F.

Stochastic

Multi-item

Lot-sizing

Problem

with

Bounded

Number

Setups.

DOI:

10.5220/0006622501060114

Proceedings

the

7th

International

Conference

Operations

Research

and

Enterprise

Systems

(ICORES

2018),

pages

106-114

ISBN:

978-989-758-285-1

of the method is proved via extensive numerical ex-

periments on real industrial data. In particular, we

compare our results with current heuristics used in in-

dustry and with lower bounds.

2 PROBLEM FORMULATION

AND MODEL

The assembly line produces a set R of references over

T weeks. The number of distinct references produced

over a week cannot exceed N. There is also an upper

bound on the total week production (summed over all

references). We normalize all quantities so that this

upper bound is equal to 1.

The production must satisfy a random demand.

The demand of reference r over week t is a random

parameter d

, whose realization is known at the end

of week t. When production of a reference r is not

used to satisfy the demand, it can be stored but incurs

a unit holding cost h

> 0 per week. When a demand

is not satisﬁed by the production of the current week

or by inventory, it can be satisﬁed later but incurs a

unit backorder cost γ per week for some coefﬁcient

γ > 0. Note that there is no setup cost, as discussed

in the introduction. For each reference r, there is an

initial inventory s

∈ R

At the beginning of each week, before the demand

of each reference is revealed, the production of the

week has to be ﬁxed. The objective is to minimize

the total expected inventory cost (holding cost plus

backorder cost) over the whole horizon of T weeks.

Regarding randomness, we assume that for any r

and t, realizations of (d

t+1

,...,d

) have ﬁnite ex-

pectation and can be efﬁciently sampled, knowing a

realization of d

t−1

Before turning to the modeling of the problem, we

make two remarks. The ﬁrst one is about its position

in the literaure. Our problem (in its deterministic ver-

sion) is a variation of the Capacitated Lot-Sizing Pro-

blem (see (Karimi et al., 2003) for a review of mo-

dels and algorithms). The single differences with the

usual basic version of that problem is that there is no

setup cost and that there is this upper bound on the

total number of distinct references that can be produ-

ced over a week. The second remark is about the so-

mehow non-conventional assumption that the demand

is revealed at the end of the week. There are indeed

industrial cases where this occurs, like in computer

manufacturing. However, the methods developed in

this paper can be easily adapted to the case where the

demand is revealed at the beginning of the week.

The production problem at the beginning of week

t can be modeled as the following stochastic program

(S):

min E

∑

r∈R

+ γb

)

s.t.

− b

∈ Jt,TK, r ∈ R

= s

−1

+ q

− d

∈ Jt,TK, r ∈ R

∑

r∈R

≤ 1 t

∈ Jt,TK

≤ x

∈ Jt,TK, r ∈ R

∑

r∈R

≤ N t

∈ Jt,TK

∈

{

0,1

}

∈ Jt,TK, r ∈ R

, b

≥ 0 t

∈ Jt,TK, r ∈ R

σ(q

) ⊂ σ



,...,d

−1

)

∈R



∈ Jt,TK, r ∈ R ,

where the variable

(resp. b

) models the inventory

(resp. the backorder) of reference r at the end of week

, the variable q

models the quantity of reference r

produced over week t

, and the variable x

takes the

value 1 if the reference r is produced over week t

and 0 otherwise. All these variables are random. The

last constraint of the program, written as a measurabi-

lity constraint, means that the values of the variables

can only depend on the values taken by the de-

mand before time t

(the planner does not know the

future). The reader would have probably expected si-

milar constraints for the other variables, but it is easy

to see that adding them changes neither the optimal

value, nor the production decisions. Note that in this

program the inventory s

t−1

is not a variable and has a

known deterministic value at the beginning of week t.

A feasible solution (q

)

≥t,r∈R

of (S) provides a

deterministic production (q

)

r∈R

for the current week

3 FOCUS ON THE

DETERMINISTIC VERSION

The deterministic version of (S) is obtained by remo-

ving the measurability constraint and by considering

all d

’s deterministic, denoted then simply by d

. The

variables are then also denoted with non-bold charac-

ters. The main message of this section is that this ver-

sion is already difﬁcult.

The deterministic version is NP-hard in the strong

sense for any ﬁxed N ≥ 3, since there is a straight-

forward reduction from 3-PARTITION. Given a set

of 3n positive integers {a

,...,a

}, this latter pro-

blem consists in deciding whether this set admits a

partition into n triples of same sum. It is a notorious

strongly NP-complete problem (Garey and Johnson,

1979). The problem 3-PARTITION can be reduced to

the case N = 3 where each index in J1,3nK is a dis-

tinct reference r whose demand d

is equal to 0 for all

period t, except for the last period T = n where it is

equal to

∑

. (The cases N ≥ 4 are dealt with

via a similar reduction.) The complexity status of the

deterministic version when N = 1 or N = 2 seems to

be a challenging open question.

It is also worth noting that the optimal value of

the continuous relaxation is independent of N. This

can be seen by considering the mixed integer program

obtained by removing the constraints

∑

r∈R

≤ N:

any feasible solution of that program provides a fe-

asible solution of the original program of same va-

lue, simply by redeﬁning x

= q

for all t

and all r if

necessary. In an attempt of improving the quality of

the continuous relaxation, one may consider the ex-

tended formulation with the binary variables y

for

each p ∈





and each t ∈ J1,T K in place of the x

’s

(indicating the references produced on period t) but it

can be shown with a nontrivial proof that it does not

improve the quality of the continuous relaxation. Si-

milarly, the extended formulation with the binary va-

riables z

for each τ ⊆ J1,TK (indicating the periods

of production of the reference r) does not improve

the quality of the relaxation (this time, with an easy

proof).

4 METHOD

As explained in Section 3, the deterministic version of

(S) is difﬁcult. Therefore, we cannot expect a quick

algorithm solving exactly the problem, and this holds

especially for the full stochastic version. Moreover,

one of the requests of the partner was to have an easy

to understand method, which can be used and main-

tained in practice, with short computation times. We

propose a two-stage approximation consisting in re-

placing the measurability constraint by



σ(q

) ⊂ σ(∅) r ∈ R

σ(q

) ⊂ σ



,...,d



∈R



≥ t + 1, r ∈ R ,

which provides a relaxation of the initial program: the

production decisions for the current week t can still

not depend on the future, but now the subsequent pro-

duction decisions depend on the future demand. We

denote this relaxation by (2SA).

This approximation is a two-stage approximation

as we distinguish between two levels of information

over the uncertainty: production decisions for the ﬁrst

week are the ﬁrst stage variables, while all other de-

cisions are second stage variables. Three-stages or

more generally multistage approximation would give

better approximations of (S) but increases exponenti-

ally the number of variables. We chose for practicabi-

lity reasons to stick to the two-stage approximation.

The (2SA) relaxation is then solved by a classical

sample average approximation, see (Kleywegt et al.,

2002) for a presentation of the method. We build a

set Ω of m scenarios sampled uniformly at random.

Each of these scenarios is a possible realization of

t+1

,...,d

) for each r. The parameter m is ﬁxed

prior to the resolution.

We get the following mixed integer program

(2SA-m), solved by any standard MIP solver.

min

∑

ω∈Ω

∑

r∈R



˜s

,ω

+ γb

,ω



s.t.

,ω

= ˜s

,ω

− b

,ω

∈ Jt,TK, r ∈ R , ω ∈ Ω

,ω

= s

−1,ω

+ q

,ω

− d

,ω

∈ Jt,TK, r ∈ R , ω ∈ Ω

∑

r∈R

,ω

≤ 1 t

∈ Jt,TK, ω ∈ Ω

,ω

≤ x

,ω

∈ Jt,TK, r ∈ R , ω ∈ Ω

∑

r∈R

,ω

≤ N t

∈ Jt,TK, ω ∈ Ω

t,ω

= x

r ∈ R , ω ∈ Ω

t,ω

= q

r ∈ R , ω ∈ Ω

, x

,ω

∈

{

0,1

}

∈ Jt,TK, r ∈ R , ω ∈ Ω

, q

,ω

, ˜s

,ω

, b

,ω

≥ 0 t

∈ Jt,TK, r ∈ R , ω ∈ Ω

At week t, the production is then set to be the solution

)

r∈R

found by the solver.

The validity of this method for solving (2SA) is

supported by the following proposition.

Proposition 1. The following three properties hold

when m goes to inﬁnity:

(i) The value of (2SA-m) converges almost surely to

the optimal value of (2SA).

(ii) For every m, we consider the values

( ˆq

t,m

, ˆx

t,m

)

r∈R

of the decision variables for

week t of an optimal solution of (2SA-m). Any

limit point of these values is an optimal solution

of (2SA).

(iii) Let ε > δ > 0. Assume that the random demand

)

≥t,r∈R

is such that

∃C,K, ∀u ∈ R, E[e

ukdk

] ≤ Ce

. (1)

Denote by Q

(resp. Q

) the set of all possible va-

lues of ( ˆq

t,m

)

r∈R

in a δ-optimal solution of (2SA-

m) (resp. in an ε-optimal solution of (2SA)). Then

for every α ∈ (0,1), we have P(Q

⊆ Q

) > 1 −α

for m large enough.

If the random demand d is bounded or Gaussian

then it satisﬁes (1). (Some companies work with pos-

sibly negative demands. Assuming a normal distribu-

tion is hence not irrelevant.)

The proof of part (iii) relies on the following

technical lemma.

Lemma 2. Consider g(d) = inf

y∈Y

G(y,d), whereY is

non-empty and where the function G is non-negative

and κ-Lipschitz with respect to d. If the random va-

riable d satisﬁes (1), then g(d) also satisﬁes (1) with

= max{1,e

G(y

,0)

C} and K

= Kκ + G(y

,0), for

∈ Y .

Proof. For u ≥ 0 we take y

∈ Y yielding g(d) ≤

G(y

,0) + κkdk. We then have E[e

ug(d)

] ≤

E[e

uG(y

,0)

κukdk

] ≤ C

. For u ≤ 0, by non-

negativity of G we have g(d) ≥ 0, hence E[e

ug(d)

] ≤

Proof of Proposition 1. Let Q

be the (bounded) set

of feasible values for the ﬁrst-stage variables q

)

r∈R

in (2SA). Denote by F(q

,d) the minimal

cost of (2SA) that can be reached when the ﬁrst stage-

variables are ﬁxed to q

∈ Q

and the realization of

the demand is d = (d

)

≥t,r∈R

. We introduce the

map f : q

7→ E[F(q

,d)], which associates to a gi-

ven choice of q

∈ Q

the expected minimal cost, and

similarly the map

, which associates to a q

∈ Q

the minimal cost of (2SA-m) when the ﬁrst-stage va-

riables are set to this q

The map f is continuous and Q

is compact. Thus

f is bounded, and, by (Shapiro et al., 2009, Theorem

7.48), we have that (

))

m∈Z

converges to f (q

)

uniformly on Q

. Then (i) and (ii) are direct conse-

quences of (Shapiro et al., 2009, Theorem 5.3).

By Lemma 2, there exist K and C such that for

any q

∈ Q

, F(q

,d) satisfy (1). Consequently there

exists σ > 0 such that for all q

∈ Q

, the random

variable [F(q

,d) − f (q

)] − [F(q

,d) − f (q

)] is σ-

subgaussian, (see e.g., (Vershynin, 2010)). Further-

more, for any demand d, the map F(·, d) is Lipschitz-

continuous on Q

. Then, according to (Shapiro et al.,

2009, Theorem 5.18), for every α ∈ (0,1), there exists

M ∈ Z

such that for m ≥ M, we have P(Q

⊆ Q

) >

1 − α.

5 NUMERICAL EXPERIMENTS

5.1 Instances

C++11 has been chosen for the implementations and

Gurobi 6.5.1 was used to solve the model on a PC

with Intel(R) Core(TM) i7-3770 CPU @ 3.40GHz

and 8Go RAM.

The instances used are realistic and have been pro-

vided by a client of the partner. The client gave actu-

ally the ﬁgures of seven assembly lines but we give

the results for only two of them: L2 and L6. The line

L2 experiences overcapacity: the ratio expected de-

mand over total capacity is smaller than 1 (actually

equal to 0.7). The line L6 experiences undercapacity:

the ratio expected demand over total capacity is lar-

ger than 1 (actually equal to 1.1). The horizon T is

the typical one used in practice by this client, namely

T = 13 weeks (a quarter). The demand is obtained

via d

= (

+ e

)

, where the

’s are historical data

and where e

is a generalized autoregressive process

so that e

= 0.25e

t−1

+ 0.75ε

, where ε

is a Gaussian

white noise process with zero mean and standard de-

viation equal to v

where v is the “volatility” and is

chosen in

{

0.2,0.5

}

. The initial inventory is set to

(

). The other parameters are pro-

vided in Table 1. In particular, for each value of v,

we have considered three possible values of the unit

backorder cost γ, which have been determined follo-

wing a procedure described in Section 6. At the mo-

ment, we do not discuss further these values and take

them as part of the input, as required by the problem

formulation. The parameter C is the capacity of the

line before normalization. (Recall that problem and

the model have been formulated in Section 2 after nor-

malization.) In the column

, we indicate the range

of the holding costs before normalization. We obtain

the h

’s by dividing these costs by C.

The number m of scenarios used to solve (2SA-m)

is ﬁxed to 20, determined by preliminary experiments

showing that it is a good trade-off between accuracy

and tractability. The time limit of the solver has been

set to 90 seconds.

5.2 Other Heuristics

Our method is compared with three other heuristics.

The ﬁrst heuristic is the deterministic version

of (S), where the random demand is replaced by its

expectation.

The second one, the lot-size heuristic, consists in

determining before the ﬁrst week once and for all a

value `

∗

for each reference r ∈ R . At time t, if the in-

ventory of reference r is below a precomputed safety

level, the quantity q

is chosen so that the inventory of

reference r exceeds the safety level of exactly `

∗

. In

case of capacity issues, the production is postponed

and thus backorder costs appear. In addition, if some

capacity issues are easily anticipated, the production

of a reference r can be activated even if the inventory

Table 1: Instance characteristics.

Instances Instance characteristics

max(

) C N

v γ

L2_v20_13 21 4992 10562 7 35–61 0.2 13

L2_v20_81 81

L2_v20_203 203

L2_v50_48 0.5 48

L2_v50_154 154

L2_v50_341 341

L6_v20_3 22 8640 13299 8 16–23 0.2 3

L6_v20_19 19

L6_v20_55 55

L6_v50_11 0.5 11

L6_v50_42 42

L6_v50_98 98

is not below the safety level.

The third one, the cover-size heuristic is almost

the same, but instead of precomputing a ﬁxed quantity

for each reference, a duration T

∗

is ﬁxed before the

ﬁrst week. When the inventory of reference r is below

the safety level, the quantity q

is computed so that the

inventory of reference r exceeds the safety level of the

expected demand for the next T

∗

weeks.

The values `

∗

and T

∗

are determined as follows.

∗

)

r∈R

is actually chosen to be the optimal solution

of the following convex program, which somehow

considers the problem at a “macroscopic” level. (Si-

milar convex programs in the same context have been

considered in the literature; see (Ziegler, 1982) for ex-

ample.)

min

∑

r∈R

s.t.

∑

r∈R

≤ N

> 0 r ∈ R ,

where

= E



∑

t=1



The parameter `

∗

of the lot-size heuristic is then

set to

∗

The safety levels have been provided by the part-

ner and are those used in practice.

The cover-size heuristic adapts the production to

the realization of the demand, contrary to the lot-

size heuristic. According to our partner, it makes the

cover-size heuristic more suitable for situations with

low short term volatility of demand or for overcapa-

citated lines, while the lot-size heuristic is expected

to behave better with high short term volatility of de-

mand or for undercapacitated lines.

Notice that the backorder costs are not taken into

account at all for determining the values of the pa-

rameters `

∗

and T

∗

. But playing with safety levels

allows to prevent too large backorder costs. However,

in real life it is usually the other way round: the com-

pany does not associate costs to backorder and aims

at keeping the total amount of unsatisﬁed demand be-

low some predetermined level. We come back to this

point later in Section 6.

5.3 Results

The results are provided in Table 2. All quantities

are in Me and are given with a conﬁdence interval at

95%.

The column LB provides the lower bound obtai-

ned by the optimal value at time t = 1 of program

(2SA-m) (with m = 1000 and a time limit of 24 hours

for the solver).

The column 2SA-m is the cost of the method pro-

posed in Section 4. (We remind the reader that we

propose m = 20 in this case.) The three next columns

provide the results for the three heuristics described

in Section 5.2. For these four columns, between 27

and 30 runs have been used for each instance.

Complementary numerical results of these expe-

riments are provided in Tables 3 and 4, discussed in

Section 6 where other criteria are considered. Table 3

can already by of interest since it deals with holding

costs obtained for the L2 instances.

5.4 Comments

Our method clearly outperforms lot-size and cover-

size heuristics and is better than the deterministic ap-

proximation for all but one instance. By running our

method instead of a usual heuristic at the beginning

of each week, the inventory costs can be reduced of-

ten by more than 50%. For the instance L2_v20_13,

the inventory costs have been divided by more than 6

(which corresponds to several Me).

Table 2: Results - Inventory costs (in Me).

Instances LB 2SA-m Det. Cover-size Lot-size

L2_v20_13 0.53 0.89 ± 0.03 1.17 ± 0.10 6.95 ± 0.17 7.79 ± 0.14

L2_v20_81 0.94 2.29 ± 0.06 2.36 ± 0.07 8.12 ± 0.19 9.65 ± 0.14

L2_v20_203 1.00 3.05 ± 0.07 3.25 ± 0.08 9.35 ± 0.29 10.99 ± 0.19

L2_v50_48 0.97 2.73 ± 0.11 3.06 ± 0.21 8.03 ± 0.26 8.37 ± 0.21

L2_v50_154 1.36 4.54 ± 0.20 5.06 ± 0.33 10.83 ± 0.53 11.20 ± 0.38

L2_v50_341 1.51 5.91 ± 0.25 7.90 ± 0.66 15.17 ± 1.21 14.65 ± 0.77

L6_v20_3 0.54 0.61 ± 0.01 0.70 ± 0.02 1.71 ± 0.08 1.74 ± 0.08

L6_v20_19 1.41 1.81 ± 0.06 1.86 ± 0.06 3.51 ± 0.12 3.20 ± 0.08

L6_v20_55 2.67 3.57 ± 0.24 3.71 ± 0.30 7.49 ± 0.39 6.24 ± 0.34

L6_v50_11 1.33 2.00 ± 0.11 2.14 ± 0.12 3.42 ± 0.15 3.03 ± 0.13

L6_v50_42 2.99 4.45 ± 0.53 4.48 ± 0.51 7.99 ± 0.62 6.57 ± 0.60

L6_v50_98 6.13 8.29 ± 1.23 7.94 ± 1.04 16.34 ± 1.61 12.96 ± 1.45

We may note that the cover-size heuristic behaves

better than the lot-size heuristic on the L2 instances,

while the converse holds on the L6 instances. It is

in line with our partner’s opinion regarding their be-

haviors with respect to the capacity of the line (see

Section 5.2). Regarding their behaviors with respect

to the volatility, it does not seem to be possible to

draw any concrete conclusion.

The 2SA-m algorithm requires 90 seconds to out-

put a solution, while lot-size and cover-size heuris-

tics take less than a second and the deterministic ap-

proximation less than 10 seconds. The 2SA-m algo-

rithm is thus slower, but note that 90 seconds to be

run only once at the beginning of each week remains

very short. Moreover, even with improved computers,

the lot-size and cover-size heuristics and the determi-

nistic approximation will not change their output (in

all our experiments, Gurobi always found the optimal

solution of the deterministic approximation). This is

not the case for 2SA-m, which means that it would

beneﬁt from improved computational capacities.

Otherwise, complementary ﬁgures not provided

here show that the deterministic approximation and

2SA-m use almost every available setups (between

97% and 100%) whereas lot-size and cover-size heu-

ristics cannot take full advantage of this ﬂexibility

(only 84% to 98% of setups are used).

Our method dramatically reduces the inventory

costs, while remaining a quite simple approach. It

provides without any doubt a positive answer to the

industrial request. Moreover, as explained in the next

section, our method can be used to address the pro-

blem of minimizing the holding costs, while keeping

the backorder at a reasonable level.

6 CONTROLLING THE FILL

RATE SERVICE LEVEL

6.1 The Problem

In practice, except when they are enshrined through

contracts with the client, backorder costs can be hard

to estimate. Thus, we consider here an alternative ver-

sion of the problem, which often meets the objecti-

ves in industry. Everything remains the same, except

that there is no backorder cost anymore, which means

that we remove the corresponding part in the objective

function of (S) and that the following constraint is ad-

ded:

∑

r∈R

∑

t=1

min(d

t−1

)

∑

r∈R

≥ β, (2)

where β ∈ [0,1] is the desired ﬁll rate service level de-

ﬁned by the company, and where d

∑

t=1

. This

constraint implies that, in expectation, the proportion

of the production delivered on-time is at least β. In-

deed, at week t, we deliver on-time the minimum bet-

ween the demand and the sum of the real inventory at

the end of week t − 1 and of the production of week t.

First, note that this constraint does not make a dif-

ference between a one-week delay and a two-week

delay. Second, this constraint is a global constraint to

satisfy over the whole horizon: it does not depend on

the current week. This is a matter of modeling point

of view since we were not able to get a precise formu-

lation from our partner and from its clients. A draw-

back of this formulation with respect to the one given

in Section 2 is that this additional constraint cannot

always be satisﬁed. In the next section, we propose a

way to address this new constraint in an approximate

way, while being able to provide solutions in any case.

6.2 Surrogate Backorder Costs

We address the alternative formulation by deﬁning a

surrogate backorder coefﬁcient γ before the ﬁrst week,

with the idea to use the algorithm of Section 4 and to

heuristically entice it to choose solutions satisfying

constraint (2).

γ := max

r∈R

P[d

≤ q

(β)]

P[d

> q

(β)]

(β) := min



q ∈ R





min(d

,q)



≥ β



(3a)

(3b)

(3c)

where d

∑

t=1

is the demand of reference r ag-

gregated over time. Since d

is integrable, q

(β) is

well-deﬁned (we set

= β so that references with no

demand would not impact the constraint). Computing

an approximate value of q

(β) at an arbitrary preci-

sion can easily be performed by binary search.

To justify this choice, consider the original pro-

blem of Section 2 with only one reference and for a

horizon of one week. Assuming no initial inventory,

it takes then the form of the famous newsvendor pro-

blem (see e.g., (Shapiro et al., 2009, Chapter 1))

min

q≥0



(q − d

)

+ γ

− q)



, (4)

where γ

is a unit backorder cost speciﬁc to reference

r. The next proposition means that with the right

choice for γ

, the alternative version with a desired

ﬁll rate service level β (which takes the form of (3b)

since h

> 0) is equivalent to the original one (4). Of

course, it holds only for the case with one reference

and a horizon of one week. Many references lead to

several possible values γ

for the surrogate backorder

cost. The choice in (3a) of taking the maximum of

them for our actual γ aims at facing the additional un-

certainty induced by the limited number of setups per

week and the length of the horizon.

Proposition 3. Deﬁne γ

as in (3b). Then q

(β) is the

smallest optimal solution to (4).

Proof. The aggregated production problem (4)

is known to have the optimal solution q

r∗

−1

(γ

/(γ

)), where F

−1

is the left-inverse of the

cumulative distribution function of d

, i.e., F

−1

(κ) =

inf{q | P(d

≤ q) ≥ κ}. Since we have set γ

P[d

≤q

(β)]

P[d

(β)]

, we have q

r∗

= inf{q | P(d

≤ q) ≥ P(d

≤

(β)}, which implies that q

r∗

≤ q

(β). Now if

this inequality were strict, then it would mean that

P(d

∈ (q

r∗

(β))) = 0, which contradicts the mi-

nimality assumption in the deﬁnition of q

(β) (Equa-

tion (3c)).

Note that this formulation does not take into ac-

count the capacity constraint. It is therefore probably

better suited to overcapacited production lines.

Remark. If instead of controlling the ﬁll rate service

level, we want to control the cycle service level, deﬁ-

ned as the probability of satisfying the whole demand,

then we can choose

γ =

1 − β

max

r∈R

. (5)

Indeed, in this case, the optimal solution q

r∗

of (4)

satisﬁes P(q

r∗

≥ d

) = β. Interestingly, Equation (5)

does not depend on the distribution of the demand,

which contrasts with the ﬁll rate service level.

6.3 Numerical Results

We discuss brieﬂy the results of the experiments des-

cribed in Section 5 in the new context of this section.

Tables 3 and 4 provide the numerical results: the ﬁrst

table shows the holding costs obtained when using the

method proposed in this work; the second table shows

the average (over the 30 runs) of the ﬁll rate service

level measured ex post, which we deﬁne by the reali-

zation of

∑

r∈R

∑

t=1

min(d

)

∑

r∈R

We plot on Figure 1 the values of the holding costs

obtained for each of the 30 runs for the instance L2

with β = 95% and v = 0.5. On Figure 2, we plot the

average of these values for three desired ﬁll rate ser-

vice levels (85%, 95%, and 98%).

Again, the inventory costs (which are in this ver-

sion also the holding costs) are dramatically better

than the ones obtained by the heuristics used by the

clients of our partner. These two heuristic are howe-

ver able to provide almost always very good ﬁll rate

service level but it comes at the price of very high

holding costs

When β ∈ {95%, 98%}, the deterministic approx-

imation reaches better holding costs than 2SA-m but

worse ﬁll rate service levels. This is especially true

when the volatility is high (v = 0.5). Further experi-

ments would probably be necessary to understand bet-

ter how the deterministic approximation and 2SA-m

compare but since they both use the surrogate backor-

der cost, their behavior regarding the ﬁll rate service

level shows that the way we compute it is relevant.

For high desired ﬁll rate service levels (≥ 95%),

the method 2SA-m behaves well: the ﬁll rate service

level measured in practice is close to the desired level.

Since constraint (2) is an approximation of what is

sought in practice, these results can be considered as

Table 3: Results - Holding costs (in Me).

Instances Results

name β γ 2SA-m Det. Cover-size Lot-size

L2_v20 85% 13 0.44 ± 0.04 0.82 ± 0.12 6.85 ± 0.18 7.70 ± 0.15

95% 81 1.58 ± 0.07 1.56 ± 0.11 7.49 ± 0.20 9.14 ± 0.17

98% 203 2.34 ± 0.10 1.97 ± 0.10 7.79 ± 0.21 9.73 ± 0.19

L2_v50 85% 48 1.55 ± 0.12 1.78 ± 0.18 6.96 ± 0.30 7.58 ± 0.24

95% 154 2.94 ± 0.19 2.55 ± 0.19 7.47 ± 0.32 8.80 ± 0.28

98% 341 3.87 ± 0.17 2.68 ± 0.19 7.75 ± 0.34 9.30 ± 0.34

Table 4: Results - Fill rate service level (in %).

Instances Results

name β γ 2SA-m Det. Cover-size Lot-size

L2_v20 85% 13 76.3 ± 0.7 80.3 ± 0.8 96.2 ± 0.4 94.7 ± 0.6

95% 81 92.5 ± 0.5 92.0 ± 0.7 96.2 ± 0.4 95.1 ± 0.6

98% 203 96.8 ± 0.4 94.9 ± 0.5 96.2 ± 0.4 95.1 ± 0.5

L2_v50 85% 48 80.5 ± 1.1 81.2 ± 1.6 92.0 ± 1.1 90.7 ± 1.0

95% 154 90.9 ± 0.9 87.5 ± 1.2 92.1 ± 1.1 91.3 ± 0.9

98% 341 94.6 ± 0.6 88.2 ± 1.2 92.1 ± 1.1 91.3 ± 0.9

80 85 90 95 100

Fill rate service level (in %)

Holding costs (in Me)

2SA-m

Det.

Cover-size Lot-size

Figure 1: Results for the 30 runs of L2 with v = 0.5 and

β = 95%.

positive and the method of Section 4 combined with

our rule of thumb for determining γ forms an effective

way to compute the production levels in real-world

settings.

7 CONCLUSION

In this paper, we focus on a multi-item lot-sizing pro-

blem with a constraint not yet addressed by academic

75 80 85 90 95 100

85%

95%

98%

85%

95%

98%

85%

95%

98%

85%

95%

98%

Fill rate service level (in %)

Holding costs (in Me)

2SA-m

Det.

Cover-size Lot-size

Figure 2: Average results for L2 for various values of desi-

red service level β and v = 0.5.

works, while being often met in practice: an explicit

upper bound on the weekly number of setups. We mo-

del this problem as a stochastic program that includes

this constraint and we propose a repeated two-stage

approximation to solve it. Our method proves its ef-

ﬁciency on real-word instances and outperforms the

heuristics currently used.

The costs in the previous model originate from

holding and backorder. According to the consulting

partner we work with, companies approach their pro-

duction decisions by considering only holding costs,

while trying to keep the backorder below some thres-

hold. Building on the method proposed for the previ-

ous model, we propose a rule of thumb for determi-

ning surrogate backorder costs that bias the solutions

towards the desired service levels. This rule relies on

the famous newsvendor problem, which is the special

case of our problem when there is one reference and

one week. There exist generalizations of the newsven-

dor problem with more references or more weeks; it

would be worth checking whether they could improve

the way this surrogate cost is computed.

REFERENCES

Aloulou, M. A., Dolgui, A., and Kovalyov, M. Y. (2014).

A bibliography of non-deterministic lot-sizing mo-

dels. International Journal of Production Research,

52(8):2293–2310.

Díaz-Madroñero, M., Mula, J., and Peidro, D. (2014). A

review of discrete-time optimization models for tacti-

cal production planning. International Journal of Pro-

duction Research, 52(17):5171–5205.

Garey, M. and Johnson, D. (1979). Computers and Intrac-

tability: A Guide to the Theory of NP-Completeness.

Series of books in the mathematical sciences. W. H.

Freeman.

Gicquel, C., Minoux, M., and Dallery, Y. (2008). Capaci-

tated lot sizing models: a literature review. working

paper or preprint.

Karimi, B., Ghomi, S. F., and Wilson, J. (2003). The ca-

pacitated lot sizing problem: a review of models and

algorithms. Omega, 31(5):365 – 378.

Kleywegt, A. J., Shapiro, A., and de Mello, T. H. (2002).

The sample average approximation method for sto-

chastic discrete optimization. SIAM Journal on Op-

timization, 12(2):479–502.

Mula, J., Poler, R., García-Sabater, J. P., and Lario, F.-C.

(2006). Models for production planning under uncer-

tainty: A review. International Journal of Production

Economics, 103(1):271–285.

Quadt, D. and Kuhn, H. (2008). Capacitated lot-sizing with

extensions: a review. 4OR: A Quarterly Journal of

Operations Research, 6(1):61–83.

Rubaszewski, J., Yalaoui, A., Amodeo, L., Yalaoui, N., and

Mahdi, H. (2011). A capacitated lot-sizing problem

with limited number of changes. In International Con-

ference on Industrial Engineering and Systems Mana-

gement, IESM, Metz, France, May 25-27.

Shapiro, A., Dentcheva, D., and Ruszczy

nski, A. (2009).

Lectures on stochastic programming: modeling and

theory. SIAM.

Vershynin, R. (2010). Introduction to the non-asymptotic

analysis of random matrices. arXiv preprint

arXiv:1011.3027.

Ziegler, H. (1982). Solving certain singly constrained con-

vex optimization problems in production planning.

Operations Research Letters, 1(6):246 – 252.