Sunday, December 17, 2017

A production scheduling 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:


  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).


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


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:

&\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

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.


  1. Similar to cutting stock problem but not quite,
  2. Multiplication of a continuous and a binary variable,
  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.