## Sunday, December 17, 2017

### A production scheduling problem

#### Problem

A post in [1] describes an interesting problem:

I have read about cutting stock problem, but this is a bit different.
We have orders in different variants, and a maximum sized machine to produce it.
Variants would be X,S,XL,L and ordered quantities are 100,40,40,80.
Say, the machine width is 6. This means, we can put 6 different variants
together and produce it.
We can put 2 X,1 S,1 XL,2 L, this means if we produce it 50 times,
output is :
X = 100 (0 waste)
S = 50 (10 waste)
XL = 50 (10 waste)
L = 100 (20 waste)

Total of 40 waste in 300 produced.
Another aproach to reduce waste would be creating 2 different variation. We
can put 4 X, 2 S and produce it 25 times, with 10 waste and make another
setup and put 2 XL,4 L and produce it 20 times with no waste. With total
10 waste we handled this production in 2 setups.
Since setup has a price, we would prefer first setup, or depending on
quantities, we may choose the second one.
I have read about cutting stock and it looks similar to this one, but
ability to divide quantities between different setups, this has more
potential to optimize and therefore more complex.

#### High level model

We can formulate this as a mixed-integer quadratic model. The data looks like:

The optimization model is as follows:

Notes:

1. The variables $$run_r$$ indicate if we execute a production run or not. Within a production run we use the same pattern.
2. We have $$run_r=0 \Rightarrow pattern_{v,r}=0, runlen_r=0$$. We could drop one of these constraints but not both.
3. Equation $$edemand$$ has a nasty quadratic term, It causes the model to be non-linear and non-convex. Of course we can solve this with a global solver.
4. The $$order$$ constraint makes sure all active runs (ie. $$run_r=1$$) are before the inactive ones ($$run_r=0$$). This makes the solution look nicer and will reduce some symmetry.
5. The funny $$\forall r-1$$ means we generate the $$order$$ constraint for all $$r$$ except the first one.
6. The model type is actually MIQCP (Mixed Integer Quadratically Constrained Program). There are commercial solvers for convex problems of this type, but for this non-convex problem I used a global MINLP solver (Baron, Couenne, Antigone).

#### Results

The input data can be summarized as:

----     62 PARAMETER costwaste  cost of waste

X  1.000,    S  2.000,    XL 3.000,    L  4.000

----     62 PARAMETER costsetup            =      100.000  fixed cost of a run

----     62 PARAMETER demand

X  100,    S   40,    XL  40,    L   80


The solution looks like:

----     62 VARIABLE runlen.L  length of a run

run1 20,    run2 40

----     62 VARIABLE pattern.L  production pattern of a run

run1        run2

X            1           2
S                        1
XL           2
L            2           1

----     62 VARIABLE waste.L  waste per variant

( ALL          0. )

----     62 VARIABLE wastecost.L           =        0.000  cost of waste
VARIABLE setupcost.L           =      200.000  cost related to setups
VARIABLE cost.L                =      200.000  total cost


We see we use two production runs. We are under-utilizing the machine (a pattern of 5 and 4 items). It is noted that this solution is not unique. A different solver will likely give a different solution.

#### Shorter  run lengths

I believe we are missing some objective in this model. Suppose we have a demand of X=60 only (zero demand for the other variants). We can organize in a single run in at least the following ways:

• Pattern X=1. This gives a run length of 60. Total cost = 100 (cost of one setup).
• Pattern X=6. This gives a run length of 10. Total cost = 100 (cost of one setup).
I would guess the second configuration is better. Probably we should add some cost for the total sum of the run lengths.

If we add a unit cost of 1 for operating the machine, the total operating cost can be expressed as $$\sum_r runlen_r$$. If we add this to the total cost we get the following solution:

----     68 VARIABLE runlen.L  length of a run

run1  4,    run2 40

----     68 VARIABLE pattern.L  production pattern of a run

run1        run2

X            5           2
S                        1
XL                       1
L                        2

----     68 VARIABLE waste.L  waste per variant

( ALL          0. )

----     68 VARIABLE wastecost.L           =        0.000  cost of waste
VARIABLE setupcost.L           =      200.000  cost related to setups
VARIABLE opcost.L              =       44.000  operating cost
VARIABLE cost.L                =      244.000  total cost


A visualization of these two production schedules can look like:

 Production Schedules

#### Linearization

If we want to use a high-performance MIP solver, we need to linearize equation $$edemand$$. I don't know how to linearize an integer variable times an integer variable. However, if we can replace an integer variable by a series of binary variables, then we are in more favorable territory.

The machine can handle up 6 variants at the time. This means $$pattern_{v,r}\le 6$$. This makes $$pattern$$ a good candidate for a binary expansion:

\begin{align}&pattern_{v,r} = \sum_k k \cdot \delta_{v,r,k}\\&\sum_k \delta_{v,r,k} \le 1\\& \delta_{v,r,k} \in \{0,1\}\end{align}

where $$k$$ runs from 1 to $$cap=6$$. Now we can write the $$edemand$$ equation as:

$\sum_{r,k} k \cdot ( runlen_r \cdot \delta_{v,r,k}) = demand_v + waste_v$

The multiplication of a binary variable and a continuous variable can be written as linear inequalities [2]. This results in:

\begin{align} & \sum_{r,k} k \cdot w_{v,r,k} = demand_v + waste_v\\ & w_{v,r,k} \le \delta_{v,r,k} \cdot maxrunlen\\ & w_{v,r,k} \le runlen_r\\ & w_{v,r,k} \ge runlen_r - maxrunlen \cdot (1-\delta_{v,r,k}) \\ &w_{v,r,k}\ge 0 \end{align}

We can also go the other route: expand $$runlen_r$$. This would lead to:

\begin{align} &\sum_{r,\ell} \ell \cdot v_{v,r,\ell} = demand_v + waste_v\\ &\sum_{\ell} \ell \cdot \rho_{r,\ell} = runlen_r\\ &\sum_{\ell} \rho_{r,\ell} \le 1\\ & v_{v,r,\ell} \le cap \cdot \rho_{r,\ell}\\ & v_{v,r,\ell} \le pattern_{v,r}\\ & v_{v,r,\ell} \ge pattern_{v,r} - cap \cdot (1-\rho_{r,\ell})\\ & \rho_{r,\ell} \in \{0,1\}\\ & v_{v,r,\ell} \ge 0 \end{align}

With these lengthy reformulations we end up with a linear mixed-integer program. This means we can use standard MIP solvers instead of a global nonlinear solver. In addition, MIP formulations often perform better than their non-linear counterparts. The first MIP formulation seems to work much better than the second.

----    136 PARAMETER report

obj        time       nodes

minlp.baron     244.000       2.120      58.000
mip1.cplex      244.000       1.218     673.000
mip2.cplex      244.000      12.063    3382.000


The linearization is somewhat convoluted, so having available a smaller nonlinear model can help in constructing and debugging these models. As is the case in database design (logical vs. physical design) and programming (premature optimization), starting with a simpler model that focuses on concepts is often helpful in developing optimization models.

#### References

1. Similar to cutting stock problem but not quite, https://math.stackexchange.com/questions/2557776/similar-to-cutting-stock-problem-but-not-quite
2. Multiplication of a continuous and a binary variable, http://yetanothermathprogrammingconsultant.blogspot.com/2008/05/multiplication-of-continuous-and-binary.html
3. Gislaine Mara Melega, Silvio Alexandre de Araujo and Raf Jans, Comparison of MIP Models for the Integrated Lot-Sizing and One-Dimensional Cutting Stock Problem, Pesquisa Operacional (2016) 36(1): 167-196.