Inventory Management
Find the best order policy to replenish stock levels, i.e. minimize the total cost incurred.
The (Q, R) inventory model is as follows:
- If the inventory falls below threshold R, place an order of quantity Q.
- Lead times: when an order is placed it takes some time to be delivered
- When inventory becomes zero, additional demand will cause a backlog. Backlogged demand will be fulfilled when replenishments arrive but at a cost.
Cost Contribution
- Fixed ordering cost (Fixed cost)
- Holding cost related to inventory (Holding cost)
- Penalties related to backlogs (Shortage cost)
Choose quantity Q and ordering threshold R to strike a good balance among these three costs!
- To save on the shortage cost, we want large R
- To save on the holding cost, we want small Q and small R
- To save on the fixed cost, we want large Q
There are different approaches to this problem. Here we look at a Mixed Integer Programming approach.
Model
A very accessible model description can be found here: Optimal (Q,R) Inventory Policy as a MIP:
$$ \begin{aligned} \min > & z = orderCost \cdot \sum_t order_t+holdCost \cdot \sum_t inv_t+backlogCost \cdot\sum_t back_t && (1) \\ & inv_t-back_t = inv_{t-1}-back_{t-1} - demand_t + repl_t + initInv_t && (2) \\ & inv_t \le invCap \cdot \delta_t \\ & back_t \le maxBackLogged \cdot (1-\delta_t) \\ & inv_t \le R + invCap \cdot (1-low_t) && (3) \\ & inv_t \ge R + 1 - (invCap+1) \cdot low_t \\ & order_t \le 1-low_{t-1} && (4) \\ & order_t \le low_t \\ & order_t \ge low_t -low_{t-1} \\ & repl_t \le maxQ \cdot order_{t-leadTime} && (5) \\ & repl_t \le Q \\ & repl_t \ge Q - maxQ \cdot (1-order_{t-leadTime}) \\ & Q \in [0,maxQ] \\ & R \in [0,invCap] \\ & inv_t \in [0,invCap] \\ & back_t \in [0,maxBackLogged] \\ & \delta_t \in {0,1} \\ & low_t \in {0,1} \\ & order_t \in {0,1} \\ & repl_t \in [0,maxQ] \end{aligned} $$
Implementation
Using Pyomo as modelling framework results in the following Python code for the constraints.
The numbers in the code comments match the constraint numbering in the MIP model above.
# 2a inventory balance equation:
def inv_balance_c(model, t):
if t == 1:
return (
model.inv[t] - model.backlog[t]
== -model.demand[t] + model.repl[t] + model.initInv[t]
)
return (
model.inv[t] - model.backlog[t]
== model.inv[t - 1]
- model.backlog[t - 1]
- model.demand[t]
+ model.repl[t]
+ model.initInv[t]
)
model.inv_balance_c = Constraint(model.H, rule=inv_balance_c)
# 2b only one of these variables can be non-zero
model.inv_backlog_excl_c1 = Constraint(
model.H, rule=lambda model, t: model.inv[t] <= self.invCap * model.delta[t]
)
model.inv_backlog_excl_c2 = Constraint(
model.H,
rule=lambda model, t: model.backlog[t]
<= self.maxBackLogged * (1 - model.delta[t]),
)
# 3. detect low inventory
model.low_inv_c1 = Constraint(
model.H,
rule=lambda model, t: model.inv[t]
<= model.r + self.invCap * (1 - model.low[t]),
)
model.low_inv_c2 = Constraint(
model.H,
rule=lambda model, t: model.inv[t]
>= model.r + 1 - (self.invCap + 1) * model.low[t],
)
# 4. reorder event detection
def order_event_c1(model, t):
if t == 1:
return Constraint.Skip
return model.order[t] <= 1 - model.low[t - 1]
model.order_event_c1 = Constraint(model.H, rule=order_event_c1)
def order_event_c2(model, t):
if t == 1:
return Constraint.Skip
return model.order[t] >= model.low[t] - model.low[t - 1]
model.order_event_c2 = Constraint(model.H, rule=order_event_c2)
model.order_event_c3 = Constraint(
model.H, rule=lambda model, t: model.order[t] <= model.low[t]
)
# 5. inventory replenishment takes place after an order was placed and when the lead time passed.
def repl_c1(model, t):
if t <= self.leadTime:
return model.repl[t] == 0
return model.repl[t] <= self.maxQ * model.order[t - self.leadTime]
model.repl_c1 = Constraint(model.H, rule=repl_c1)
def repl_c2(model, t):
if t <= self.leadTime:
return model.repl[t] == 0
return model.repl[t] >= model.q - self.maxQ * (
1 - model.order[t - self.leadTime]
)
model.repl_c2 = Constraint(model.H, rule=repl_c2)
model.repl_c3 = Constraint(
model.H, rule=lambda model, t: model.repl[t] <= model.q
)
model.q_c = Constraint(rule=lambda model: model.q <= self.maxQ)
model.r_c = Constraint(rule=lambda model: model.r <= self.invCap)
model.inv_c = Constraint(
model.H, rule=lambda model, t: model.inv[t] <= self.invCap
)
model.backlog_c = Constraint(
model.H, rule=lambda model, t: model.backlog[t] <= self.maxBackLogged
)
model.repl_c = Constraint(
model.H, rule=lambda model, t: model.repl[t] <= self.maxQ
)
Translating the MIP Model into Pyomo code is straight-forward. It really is a 1:1 translation of equations and constraints into Pyomo expressions.
As with any MIP modelling language special attention is needed for border conditions, i.e. when indices run out of their definition ranges. However, here you can apply all the facilities Python as a full programming language provides.
Results
CBC solved the model in 3min 27s. With a commercial solver solution time comes down to roundabout 30s.
- Total Cost: 48990
- Number of constraints : 1400
- Number of variables : 602
Q: 512.0 (quantity to order per order event)
R: 344 (threshold when order event triggers)
The results match the results in 1, so the Pyomo approach in combination with CBC seem to be valid. Once the model is in the form of a MIP, it is actually easy to solve it via Pyomo. To model the problem is the real challenge, of course.
The real learning is here: You can actually solve decent MIP problems with Open Source software.