In [1] the question was posed: "Can we formulate a (Q,R) inventory model as an integer programming model?" Of course, I said: Yes. But we need to be a bit careful if we want to make it a linear model. Let's check if I was just bluffing here...
(Q,R) Inventory Policy
The (Q,R) inventory model is as follows: if the inventory falls below R, place an order for Q. In addition, we deal with:
- Lead times: when we placed an order 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.
- Costs are:
- Fixed ordering cost
- Holding cost related to inventory
- Penalties related to backlogs
Data
---- 11 PARAMETER demand random data t1 17, t2 84, t3 55, t4 30, t5 29, t6 22, t7 34, t8 85, t9 6, t10 50 t11 99, t12 57, t13 99, t14 76, t15 13, t16 63, t17 15, t18 25, t19 66, t20 43 t21 35, t22 35, t23 13, t24 15, t25 58, t26 83, t27 23, t28 66, t29 77, t30 30 t31 11, t32 50, t33 16, t34 87, t35 26, t36 28, t37 59, t38 72, t39 62, t40 46 t41 41, t42 11, t43 31, t44 4, t45 33, t46 18, t47 64, t48 56, t49 76, t50 29 t51 66, t52 75, t53 62, t54 28, t55 8, t56 10, t57 64, t58 54, t59 3, t60 79 t61 7, t62 17, t63 52, t64 75, t65 17, t66 3, t67 58, t68 62, t69 38, t70 35 t71 24, t72 24, t73 13, t74 93, t75 37, t76 78, t77 30, t78 12, t79 74, t80 6 t81 20, t83 26, t84 49, t85 15, t86 17, t87 33, t88 31, t89 32, t90 96, t91 99 t92 36, t93 37, t94 77, t95 39, t96 91, t97 11, t98 73, t99 5, t100 57
We assume demand occurs during period \(t\), i.e. it is a "flow" quantity (opposed to "stock variables"). Note that demand for time period 82 is zero (this is not printed here).
The other data for the problem looks like:
---- 28 PARAMETER invCap = 1000.000 inventory capacity PARAMETER maxQ = 1000.000 max order quantity PARAMETER orderCost = 500.000 cost for each order PARAMETER holdCost = 2.000 holding cost (per unit, per time period) PARAMETER leadTime = 10.000 lead time PARAMETER backLogCost = 10.000 backlogged orders penalty (per unit, per time period) PARAMETER maxBackLogged = 1000.000 big-M ---- 28 PARAMETER initInv intial inventory available at t1 t1 500.000
Variables
I need a bunch of decision variables to model this. Let's list them:
| Decision variables | |
|---|---|
| \(\color{darkred}Q \in [0,\color{darkblue}{\mathit{maxQ}}]\) | Reorder quantity. Each time we place an order, this is the size of the order. |
| \(\color{darkred}R \in [0,\color{darkblue}{\mathit{invCap}}]\) | Reorder point. We place an order when the inventory drops below \(\color{darkred}R\). |
| \(\color{darkred}{\mathit{inv}}_t \in [0,\color{darkblue}{\mathit{invCap}}] \) | Inventory at end of period \(t\). |
| \(\color{darkred}{\mathit{back}}_t \in [0,\color{darkblue}{\mathit{maxBackLogged}}]\) | Backlog, unmet demand we fulfill once new inventory arrives, again at end of period \(t\). |
| \(\color{darkred}{\mathit{low}}_t \in \{0,1\}\) | Inventory is low, defined by \(\color{darkred}{\mathit{inv}}_t \le \color{darkred}R\). |
| \(\color{darkred}{\mathit{order}}_t \in \{0,1\}\) | An order is placed. We place an order when \(\color{darkred}{\mathit{inv}}_t\) drops below \(\color{darkred}R\). |
| \(\color{darkred}{\mathit{repl}}_t \in \{0,\color{darkred}Q\}\) | Replenishment. Inventory is replenished by a previously placed order (there is a lead time). This replenishment is also used to handle backlogged demand. \(\color{darkred}{\mathit{repl}}_t\) is either zero or equal to the reorder quantity \(\color{darkred}Q\). |
In the question, it was required that the decision variables are integer-valued. So, I declared variables as integer where needed.
Model equations
1. The objective has three parts: fixed ordering cost, holding cost, and backlog penalities. So we have \[\min \> \color{darkred} z = \color{darkblue}{\mathit{orderCost}}\cdot\sum_t \color{darkred}{\mathit{order}}_t+\color{darkblue}{\mathit{holdCost}}\cdot \sum_t \color{darkred}{\mathit{inv}}_t+\color{darkblue}{\mathit{backlogCost}} \cdot\sum_t \color{darkred}{\mathit{back}}_t\]
2. The inventory balance equation is a bit complex as we allow backlog: \[\color{darkred}{\mathit{inv}}_t-\color{darkred}{\mathit{back}}_t = \color{darkred}{\mathit{inv}}_{t-1}-\color{darkred}{\mathit{back}}_{t-1} - \color{darkblue}{\mathit{demand}}_t + \color{darkred}{\mathit{repl}}_t + \color{darkblue}{\mathit{initInv}}_t \] Instead of just inventory, we update the quantity \(\color{darkred}{\mathit{inv}}_t-\color{darkred}{\mathit{back}}_t \). We require that only one of these variables can be non-zero. That means we need to add the complementarity condition \[\color{darkred}{\mathit{inv}}_t\cdot \color{darkred}{\mathit{back}}_t = 0\] In many cases when we minimize cost this condition is not needed, as the model will automatically have not both variables nonzero. Unfortunately, in this case the model may not always do this, just to get some early reordering, so we need to enforce this explicitly. The complementarity condition can be linearized as: \[\begin{align} & \color{darkred}{\mathit{inv}}_t \le \color{darkblue}{\mathit{invCap}} \cdot \color{darkred}\delta_t \\ & \color{darkred}{\mathit{back}}_t \le \color{darkblue}{\mathit{maxBackLogged}} \cdot (1-\color{darkred}\delta_t) \\ & \color{darkred}\delta_t \in \{0,1\}\end{align} \] In GAMS, when indexing outside the domain, a zero is returned. That means \(\color{darkred}{\mathit{inv}}_{0}=\color{darkred}{\mathit{back}}_{0}=0\). To handle initial inventory, we use the term \(\color{darkblue}{\mathit{initInv}}_t\). This is a sparse vector with only a nonzero value in the first position. This trick makes it easy to have a single constraint for all time periods as opposed to handling the first period differently.
3. To detect low inventory, i.e. when \(\color{darkred}{\mathit{inv}}_t \le \color{darkred}R\), we add a binary variable \(\color{darkred}{\mathit{low}}_t \in \{0,1\}\). We need to make sure: \[\begin{cases}\color{darkred}{\mathit{inv}}_t \le \color{darkred}R \Rightarrow \color{darkred}{\mathit{low}}_t = 1 \\ \color{darkred}{\mathit{inv}}_t \ge \color{darkred}R + 1 \Rightarrow \color{darkred}{\mathit{low}}_t = 0 \end{cases}\] I implemented this as \[\begin{align} & \color{darkred}{\mathit{inv}}_t \le \color{darkred}R + \color{darkblue}{\mathit{invCap}} \cdot (1-\color{darkred}{\mathit{low}}_t ) \\ &\color{darkred}{\mathit{inv}}_t \ge \color{darkred}R + 1 - (\color{darkblue}{\mathit{invCap}}+1)\cdot \color{darkred}{\mathit{low}}_t\end{align} \] This implies \(\color{darkred}R \le \color{darkblue}{\mathit{invCap}} \), which is a proper bound.
4. A reorder event takes place when inventory drops below \(\color{darkred}R\). This means we should trigger a reorder only when \(\color{darkred}{\mathit{low}}_{t-1}=0\) and \(\color{darkred}{\mathit{low}}_{t}=1\). In terms of our model, this condition can be written as \[\color{darkred}{\mathit{order}}_t = (1-\color{darkred}{\mathit{low}}_{t-1}) \cdot \color{darkred}{\mathit{low}}_t\] I implemented a linearization of this as: \[\begin{align} & \color{darkred}{\mathit{order}}_t \le 1-\color{darkred}{\mathit{low}}_{t-1} \\ & \color{darkred}{\mathit{order}}_t \le \color{darkred}{\mathit{low}}_t \\ & \color{darkred}{\mathit{order}}_t \ge \color{darkred}{\mathit{low}}_t -\color{darkred}{\mathit{low}}_{t-1} \end{align}\] Note that this rule implies that after replenishment the inventory must exceed \(\color{darkred}R\) (if not, we cannot reorder).
5. Inventory replenishment takes place after an order was placed and when the lead time passed. We can write this as \[\color{darkred}{\mathit{repl}}_t = \color{darkred}Q \cdot \color{darkred}{\mathit{order}}_{t-\color{darkblue}{\mathit{leadTime}}}\] with \(\color{darkred}{\mathit{repl}}_t \in [0,\color{darkblue}{\mathit{maxQ}}]\). I linearized this as: \[ \begin{align} & \color{darkred}{\mathit{repl}}_t \le \color{darkblue}{\mathit{maxQ}} \cdot \color{darkred}{\mathit{order}}_{t-\color{darkblue}{\mathit{leadTime}}} \\ & \color{darkred}{\mathit{repl}}_t \le \color{darkred}Q \\ & \color{darkred}{\mathit{repl}}_t \ge \color{darkred}Q - \color{darkblue}{\mathit{maxQ}}\cdot (1-\color{darkred}{\mathit{order}}_{t-\color{darkblue}{\mathit{leadTime}}}) \end{align}\] For simplicity, I will assume that no ordering happened just before \(t=1\). In other words, all \(\color{darkred}{\mathit{order}}_t=0\) when \(t \le 0\). That will make a model a little bit messier, but you will need to address this in a real, practical model.
OK, this was a lot of things to digest. Let's summarize the model:
| MIP Model |
|---|
| \[\begin{align} \min \> & \color{darkred} z = \color{darkblue}{\mathit{orderCost}}\cdot\sum_t \color{darkred}{\mathit{order}}_t+\color{darkblue}{\mathit{holdCost}}\cdot \sum_t \color{darkred}{\mathit{inv}}_t+\color{darkblue}{\mathit{backlogCost}} \cdot\sum_t \color{darkred}{\mathit{back}}_t && (1) \\ & \color{darkred}{\mathit{inv}}_t-\color{darkred}{\mathit{back}}_t = \color{darkred}{\mathit{inv}}_{t-1}-\color{darkred}{\mathit{back}}_{t-1} - \color{darkblue}{\mathit{demand}}_t + \color{darkred}{\mathit{repl}}_t + \color{darkblue}{\mathit{initInv}}_t && (2) \\ & \color{darkred}{\mathit{inv}}_t \le \color{darkblue}{\mathit{invCap}} \cdot \color{darkred}\delta_t \\ & \color{darkred}{\mathit{back}}_t \le \color{darkblue}{\mathit{maxBackLogged}} \cdot (1-\color{darkred}\delta_t) \\ & \color{darkred}{\mathit{inv}}_t \le \color{darkred}R + \color{darkblue}{\mathit{invCap}} \cdot (1-\color{darkred}{\mathit{low}}_t ) && (3) \\ &\color{darkred}{\mathit{inv}}_t \ge \color{darkred}R + 1 - (\color{darkblue}{\mathit{invCap}}+1)\cdot \color{darkred}{\mathit{low}}_t \\ & \color{darkred}{\mathit{order}}_t \le 1-\color{darkred}{\mathit{low}}_{t-1} && (4) \\ &\color{darkred}{\mathit{order}}_t \le \color{darkred}{\mathit{low}}_t \\ & \color{darkred}{\mathit{order}}_t \ge \color{darkred}{\mathit{low}}_t -\color{darkred}{\mathit{low}}_{t-1} \\ & \color{darkred}{\mathit{repl}}_t \le \color{darkblue}{\mathit{maxQ}} \cdot \color{darkred}{\mathit{order}}_{t-\color{darkblue}{\mathit{leadTime}}} && (5) \\ & \color{darkred}{\mathit{repl}}_t \le \color{darkred}Q \\ & \color{darkred}{\mathit{repl}}_t \ge \color{darkred}Q - \color{darkblue}{\mathit{maxQ}}\cdot (1-\color{darkred}{\mathit{order}}_{t-\color{darkblue}{\mathit{leadTime}}}) \\ & \color{darkred}Q \in [0,\color{darkblue}{\mathit{maxQ}}] \\ & \color{darkred}R \in [0,\color{darkblue}{\mathit{invCap}}] \\ & \color{darkred}{\mathit{inv}}_t \in [0,\color{darkblue}{\mathit{invCap}}] \\ & \color{darkred}{\mathit{back}}_t \in [0,\color{darkblue}{\mathit{maxBackLogged}}] \\ & \color{darkred}\delta_t \in \{0,1\} \\ & \color{darkred}{\mathit{low}}_t \in \{0,1\} \\ & \color{darkred}{\mathit{order}}_t \in \{0,1\} \\ & \color{darkred}{\mathit{repl}}_t \in [0,\color{darkblue}{\mathit{maxQ}}] \end{align}\] |
In my experiment, all decision variables were binary or integer variables. This is optional: we can use continuous variables for \(\color{darkred}Q\) and \(\color{darkred}R\) and related variables, if this is more appropriate.
Results
---- 163 PARAMETER result demand inv backlog below reorder replenish t1 17 483 t2 84 399 t3 55 344 1 1 t4 30 314 1 t5 29 285 1 t6 22 263 1 t7 34 229 1 t8 85 144 1 t9 6 138 1 t10 50 88 1 t11 99 11 1 t12 57 68 1 t13 99 345 512 t14 76 269 1 1 t15 13 256 1 t16 63 193 1 t17 15 178 1 t18 25 153 1 t19 66 87 1 t20 43 44 1 t21 35 9 1 t22 35 26 1 t23 13 39 1 t24 15 458 512 t25 58 400 t26 83 317 1 1 t27 23 294 1 t28 66 228 1 t29 77 151 1 t30 30 121 1 t31 11 110 1 t32 50 60 1 t33 16 44 1 t34 87 43 1 t35 26 69 1 t36 28 415 512 t37 59 356 t38 72 284 1 1 t39 62 222 1 t40 46 176 1 t41 41 135 1 t42 11 124 1 t43 31 93 1 t44 4 89 1 t45 33 56 1 t46 18 38 1 t47 64 26 1 t48 56 430 512 t49 76 354 t50 29 325 1 1 t51 66 259 1 t52 75 184 1 t53 62 122 1 t54 28 94 1 t55 8 86 1 t56 10 76 1 t57 64 12 1 t58 54 42 1 t59 3 45 1 t60 79 388 512 t61 7 381 t62 17 364 t63 52 312 1 1 t64 75 237 1 t65 17 220 1 t66 3 217 1 t67 58 159 1 t68 62 97 1 t69 38 59 1 t70 35 24 1 t71 24 1 t72 24 24 1 t73 13 475 512 t74 93 382 t75 37 345 t76 78 267 1 1 t77 30 237 1 t78 12 225 1 t79 74 151 1 t80 6 145 1 t81 20 125 1 t82 125 1 t83 26 99 1 t84 49 50 1 t85 15 35 1 t86 17 530 512 t87 33 497 t88 31 466 t89 32 434 t90 96 338 1 1 t91 99 239 1 t92 36 203 1 t93 37 166 1 t94 77 89 1 t95 39 50 1 t96 91 41 1 t97 11 52 1 t98 73 125 1 t99 5 130 1 t100 57 325 1 512 Q 512 R 344 numorder 8 holdcost 37580 ordercost 4000 backlogcost 7410
A picture of this:
This looks reasonable to me.
Note that in our model the reorder point \(\color{darkred}R\) (straight light-blue line) can not be higher than the inventory after replenishment (see point 4 above). The first replenishment is exactly demonstrating this bound. Of course, we can refine the model to relax this requirement.
Conclusion
With some effort, we can form a linear MIP model that determines the optimal values for \((\color{darkred}Q, \color{darkred}R)\) for any (simulated) demand pattern. We needed a few modeling tricks to formulate this problem as a proper linear MIP model, however. This requires a number of big-M constraints. Using indicator constraints and/or quadratic terms (available in some advanced solvers) can make the model simpler. Our model does not make assumptions about special facilities available in the MIP solver, so this model should work with any (good) MIP solver.
References
- Minimize an Integer Programming problem with nested decision variables, https://stackoverflow.com/questions/64652771/minimize-an-integer-programming-problem-with-nested-decision-variables
the reference link seems broken
ReplyDeleteThe parameter `initInv` should not have the index t in the MIP model summary (2).
ReplyDeleteYes it does. It is a sparse vector with one element. That makes it easier to write just one equation for the inventory balance. No special case for the first period. This is somewhat of a standard trick in models dealing with inventory.
DeleteAh yeah, thank you. I could have read this in your description, sorry! I was usually used to something like initInv$(ord(t)=1)
Deletewhat application that you used for coding the model? Cplex or something else.
ReplyDeleteGAMS.
Delete