The robust multi-plant capacitated lot-sizing problem

Abstract

In this paper, we study the robust multi-plant capacitated lot-sizing problem with uncertain demands, processing and setup times. This problem consists of a production system with more than one production plant, in which each plant can produce items to meet its demand or transfer items to other plants. The objective is to determine a minimum-cost production and transfer plan considering the compromise between production, inventory, and transfer costs. Using a static robust optimization approach, we propose two different robust mixed-integer programming formulations for the problem. The first formulation applies the standard duality technique to the constraints involving uncertain parameters while the second applies the duality technique only to the time constraints and introduces new parameters, accumulating the worst-case demand realizations, to the inventory balance constraints. This second formulation has the advantage of resulting from a more intuitive and straightforward approach. We perform extensive computational experiments to compare the performance of the formulations and to assess the effect of different budgets of uncertainty on the solutions. Moreover, we observe that demand, processing and setup times have different impacts when taking uncertainty into account.

Appendices

Appendix 1: Compact representation of constraints (10) and (11)

This appendix describes the steps to obtain a compact form for constraints (10) and (11) when considering demand uncertainty, leading to constraints (16)–(18). These steps are described below.

• Step 1) Robust counterpart of the constraints: The robust counterpart of constraints (10) and (11) can be stated as follows:

  $$\begin{aligned}&I_{ikt}^+ \ge \sum _{\tau = 1}^{t} x_{ik\tau } + \sum _{\tau = 1}^{t} \sum _{\begin{array}{c} \ell \in {\mathcal {M}}: \\ \ell \ne k \end{array}} w_{i\ell k\tau } - \sum _{\tau = 1}^{t}\sum _{\begin{array}{c} \ell \in {\mathcal {M}}: \\ \ell \ne k \end{array}} w_{ik\ell \tau } - \sum _{\tau = 1}^{t}\left( d_{ik\tau }+\hat{d}_{ik\tau } \xi _{ik\tau }\right)&\nonumber \\&i \in {\mathcal {N}}, k \in {\mathcal {M}}, t \in {\mathcal {T}}, \forall {\varvec{\xi }} \in {{\mathcal {U}}}^{d}, \end{aligned}$$

  (33)
  $$\begin{aligned}&I_{ikt}^- \ge - \sum _{\tau = 1}^{t} x_{ik\tau } - \sum _{\tau = 1}^{t} \sum _{\begin{array}{c} \ell \in {\mathcal {M}}: \\ \ell \ne k \end{array}} w_{i\ell k\tau } + \sum _{\tau = 1}^{t}\sum _{\begin{array}{c} \ell \in {\mathcal {M}}: \\ \ell \ne k \end{array}} w_{ik\ell \tau } + \sum _{\tau = 1}^{t}\left( d_{ik\tau }+\hat{d}_{ik\tau } \xi _{ik\tau }\right)&\nonumber \\&i \in {\mathcal {N}}, k \in {\mathcal {M}}, t \in {\mathcal {T}}, \forall {\varvec{\xi }} \in {{\mathcal {U}}}^{d}. \end{aligned}$$

  (34)

  These equations impose that constraints (10) and (11) have to be satisfied for every possible realization of ${\varvec{\xi }}\in {\mathcal {U}}^d$. This representation has the disadvantage of resulting in a robust formulation that has infinite constraints (one for each realization of ${\varvec{\xi }}$) and, therefore, it cannot be solved using general-purpose optimization software. For this reason, we apply the following steps to present a compact representation for them.

• Step 2) Worst-case reformulation: Note that the worst-case scenario in constraints (33) and (34) occurs when the accumulated demands take their minimum and maximum values, respectively. The model, therefore, minimizes the inventory and shortage costs using the worst-case scenario for the demands. This can be written as

  $$\begin{aligned}&I_{ikt}^+ \ge \sum _{\tau = 1}^{t} x_{ik\tau } + \sum _{\tau = 1}^{t} \sum _{\begin{array}{c} \ell \in {\mathcal {M}}: \\ \ell \ne k \end{array}} w_{i\ell k\tau } - \sum _{\tau = 1}^{t}\sum _{\begin{array}{c} \ell \in {\mathcal {M}}: \\ \ell \ne k \end{array}} w_{ik\ell \tau } - \sum _{\tau = 1}^{t}d_{ik\tau } - \min _{{\varvec{\xi }} \in {{\mathcal {U}}}^{d}} \left\{ \sum _{\tau = 1}^{t}\hat{d}_{ik\tau } \xi _{ik\tau }\right\}&\nonumber \\&i \in {\mathcal {N}}, k \in {\mathcal {M}}, t \in {\mathcal {T}}, \end{aligned}$$

  (35)
  $$\begin{aligned}&I_{ikt}^- \ge - \sum _{\tau = 1}^{t} x_{ik\tau } - \sum _{\tau = 1}^{t} \sum _{\begin{array}{c} \ell \in {\mathcal {M}}: \\ \ell \ne k \end{array}} w_{i\ell k\tau } + \sum _{\tau = 1}^{t}\sum _{\begin{array}{c} \ell \in {\mathcal {M}}: \\ \ell \ne k \end{array}} w_{ik\ell \tau } + \sum _{\tau = 1}^{t}d_{ik\tau }+\max _{{\varvec{\xi }} \in {{\mathcal {U}}}^{d}} \left\{ \sum _{\tau = 1}^{t}\hat{d}_{ik\tau } \xi _{ik\tau }\right\}&\nonumber \\&i \in {\mathcal {N}}, k \in {\mathcal {M}}, t \in {\mathcal {T}}. \end{aligned}$$

  (36)

  Notice that in (35) the minimization problem can be redefined as a maximization problem, as follows:

  $$\begin{aligned}&I_{ikt}^+ \ge \sum _{\tau = 1}^{t} x_{ik\tau } + \sum _{\tau = 1}^{t} \sum _{\begin{array}{c} \ell \in {\mathcal {M}}: \\ \ell \ne k \end{array}} w_{i\ell k\tau } - \sum _{\tau = 1}^{t}\sum _{\begin{array}{c} \ell \in {\mathcal {M}}: \\ \ell \ne k \end{array}} w_{ik\ell \tau } - \sum _{\tau = 1}^{t}d_{ik\tau } + \max _{{\varvec{\xi }} \in {{\mathcal {U}}}^{d}} \left\{ \sum _{\tau = 1}^{t}\hat{d}_{ik\tau } \xi _{ik\tau }\right\} \nonumber \\&i \in {\mathcal {N}}, k \in {\mathcal {M}}, t \in {\mathcal {T}}. \end{aligned}$$

  (37)

• Step 3) Dualizing the primal protection function: The internal maximization problem appearing in constraints (37) and (36) is often called primal protection function, and can be written in an equivalent form as the following optimization problem:

  $$\begin{aligned}&{\mathbb {B}}_{ikt}^d=\max _{{\varvec{\xi }}} \left\{ \sum _{\tau =1}^t \hat{d}_{ik\tau } \xi _{ik\tau }: -\mathbbm {1} \le {\varvec{\xi }} \le \mathbbm {1}, \, \sum _{\tau =1}^t |\xi _{ik\tau }| \le \Gamma _{ikt} \right\} . \end{aligned}$$

  (38)

  For a given value for the uncertainty budget $\Gamma _{ikt}$, this problem is always feasible and bounded. By strong duality, it follows that the associated dual problem is also feasible and bounded, and their corresponding optimal values are equal. We then determine the dual problem associated with the maximization problem in (38). With this purpose, we rewrite it to eliminate the absolute value function that appears in one of its constraints. This can be done since in the optimal solution of (38) the values of $\xi _{ik\tau }$ will never take negative values given that the $\hat{d}_{ik\tau }$ values are non-negative. The problem can, therefore, be rewritten as

  $$\begin{aligned} {\mathbb {B}}_{ikt}^d=\max _{{\varvec{\xi }} \ge 0} \left\{ \sum _{\tau =1}^t \hat{d}_{ik\tau } \xi _{ik\tau }: \xi _{ik\tau } \le 1, \, \tau =1,\ldots ,t, \sum _{\tau =1}^t \xi _{ik\tau } \le \Gamma _{ikt} \right\} . \end{aligned}$$

  (39)

  The dual problem associated with (39) is as follows:

  $$\begin{aligned} {\mathbb {B}}_{ikt}^d = \min _{{\varvec{\lambda }} \ge 0, {\varvec{\mu }} \ge 0} \left\{ \sum _{\tau =1}^t \mu _{ik\tau t}+\Gamma _{ik\tau }\lambda _{ik\tau }: \mu _{ik\tau t}+\lambda _{ik\tau }\ge \hat{d}_{ik\tau }, \, \tau =1,\ldots , t\right\} . \end{aligned}$$

  (40)

  where $\mu _{ik\tau t}$ and $\lambda _{ik\tau }$ are the dual variables associated with the constraints of the maximization problem (38). We can now replace the maximization problem that appears in constraints (37) and (36) with the minimization problem (40), as shown as follows:

  $$\begin{aligned} I_{ikt}^+ \ge&\sum _{\tau = 1}^{t} x_{ik\tau } + \sum _{\tau = 1}^{t} \sum _{\begin{array}{c} \ell \in {\mathcal {M}}: \\ \ell \ne k \end{array}} w_{i\ell k\tau } - \sum _{\tau = 1}^{t}\sum _{\begin{array}{c} \ell \in {\mathcal {M}}: \\ \ell \ne k \end{array}} w_{ik\ell \tau } - \sum _{\tau = 1}^{t}d_{ik\tau } \nonumber \\&+ \min _{{\varvec{\lambda }} \ge 0, {\varvec{\mu }} \ge 0} \left\{ \sum _{\tau =1}^t \mu _{ik\tau t}+\Gamma _{ik\tau }\lambda _{ik\tau }: \mu _{ik\tau t}+\lambda _{ik\tau }\ge \hat{d}_{ik\tau }, \, \tau =1,\ldots , t \right\} \nonumber \\&i \in {\mathcal {N}}, k \in {\mathcal {M}}, t \in {\mathcal {T}}, \end{aligned}$$

  (41)
  $$\begin{aligned} I_{ikt}^- \ge&- \sum _{\tau = 1}^{t} x_{ik\tau } - \sum _{\tau = 1}^{t} \sum _{\begin{array}{c} \ell \in {\mathcal {M}}: \\ \ell \ne k \end{array}} w_{i\ell k\tau } + \sum _{\tau = 1}^{t}\sum _{\begin{array}{c} \ell \in {\mathcal {M}}: \\ \ell \ne k \end{array}} w_{ik\ell \tau } + \sum _{\tau = 1}^{t}d_{ik\tau }&\nonumber \\&+ \min _{{\varvec{\lambda }} \ge 0, {\varvec{\mu }} \ge 0} \left\{ \sum _{\tau =1}^t \mu _{ik\tau t}+\Gamma _{ik\tau }\lambda _{ik\tau }: \mu _{ik\tau t}+\lambda _{ik\tau }\ge \hat{d}_{ik\tau }, \, \tau =1,\ldots , t \right\}&\nonumber \\&i \in {\mathcal {N}}, k \in {\mathcal {M}}, t \in {\mathcal {T}}. \end{aligned}$$

  (42)

  The minimization that appears in constraints (41) and (42) can be omitted since if they are satisfied for some value of the ${\varvec{\lambda }}$ and ${\varvec{\mu }}$ variables, they will also be satisfied for their minimum values. This leads to the compact version shown in constraints (16)–(18). It is worth mentioning that when we eliminate the minimization term from (41) and (42), the expression $\mu _{ik\tau t}+\lambda _{ik\tau }\ge \hat{d}_{ik\tau }, \, \tau =1,\ldots , t$, will appear twice in the robust reformulation. However, since this would result in redundant constraints, we only consider them once in the reformulation as constraints (18).

Appendix 2: Compact representation of constraints (3)

This appendix describes the application of the three steps described in A to obtain a compact representation for constraints (3) when considering uncertainties in the processing and setup times.

• Step 1) Robust counterpart of the constraints: The robust counterpart of constraints (3) can be written as follows:

  $$\begin{aligned} \sum _{i \in {\mathcal {N}}} \left( b_{ikt}+\hat{b}_{ikt}\eta _{ikt}^b\right) x_{ikt}+\sum _{i \in {\mathcal {N}}} \left( f_{ikt}+\hat{f}_{ikt}\eta _{ikt}^f\right) y_{ikt} \le Q_{kt}, \, k, \in {\mathcal {M}}, t \in {\mathcal {T}}, ({\varvec{\eta }}^b,{\varvec{\eta }}^f) \in {{\mathcal {U}}}^{(b,f)}. \end{aligned}$$

  (43)

• Step 2) Worst-case reformulation: For a given solution $({\varvec{x}}^v,{\varvec{y}}^v)$, the worst-case scenario for the processing and setup times occurs when the left-hand side of (43) reaches its maximum value. Therefore, a worst-case reformulation for constraints (43) is shown below:

  $$\begin{aligned} \sum _{i \in {\mathcal {N}}} b_{ikt}x_{ikt}^v&+ \sum _{i \in {\mathcal {N}}} f_{ikt}y_{ikt}^v \nonumber \\&+ \max _{({\varvec{\eta }}^b,{\varvec{\eta }}^f) \in {\mathcal {U}}^{(b,f)}} \left\{ \sum _{i \in {\mathcal {N}}} \left( \hat{b}_{ikt}x_{ikt}^v\eta _{ikt}^b+\hat{f}_{ikt}y_{ikt}^v\eta _{ikt}^f\right) \right\} \le Q_{kt} \ k \in {\mathcal {M}}, t \in {\mathcal {T}}. \end{aligned}$$

  (44)

• Step 3) Dualizing the primal protection function: the primal protection function ${\mathbb {B}}_{kt}^{(b,f)}$ corresponding to the internal maximization problem in (44) is equivalent to the following optimization problem:

  $$\begin{aligned} {\mathbb {B}}_{kt}^{(b,f)}=\max _{{\varvec{\eta }}^b \ge 0,{\varvec{\eta }}^f \ge 0 } \left\{ \sum _{i \in {\mathcal {N}}} \left( \hat{b}_{ikt}x_{ikt}^v\eta _{ikt}^b+\hat{f}_{ikt}y_{ikt}^v\eta _{ikt}^f\right) : \eta ^b_{ikt} \le 1, \, i \in {\mathcal {N}}, \, \eta ^f_{ikt} \le 1, \, i \in {\mathcal {N}}, \, \sum _{i \in {\mathcal {N}}} \left( \eta ^b_{ikt}+\eta ^f_{ikt}\right) \le \Delta _{kt} \right\} , \end{aligned}$$

  (45)

  whose associated dual problem is

  $$\begin{aligned} {\mathbb {B}}_{kt}^{(b,f)}=\min _{{\varvec{\theta }} \ge 0, {\varvec{\beta }}^ \ge 0, {\varvec{\pi }} \ge 0 } \left\{ \sum _{i \in {\mathcal {N}}} \left( \theta _{ikt} + \beta _{ikt}\right) + \Delta _{kt}\pi _{kt}: \theta _{ikt}+\pi _{kt} \ge \hat{b}_{ikt} x_{ikt}^v, \, i \in {\mathcal {N}}, \, \beta _{ikt}+\pi _{kt} \ge \hat{f}_{ikt} y_{ikt}^v, i \in {\mathcal {N}} \right\} . \end{aligned}$$

  (46)

Finally, the compact representation of constraints (43) (which is given by constraints (19)–(21)) is obtained when the maximization problem in (44) is replaced by the minimization problem in (46), omitting the minimization function in the resulting constraint.