# Quantifying a measure of standard deviation in MILP

I am trying to set up a MILP for production scheduling. The specific details I'm not sure are important but in general a plant has M machines running N parts, each part requiring W workers. The model runs over S shifts (usually over 5 day or 10 day modelling periods).

One key feature is it's preferable for the amount of workers to stay relatively constant. One way we have thought to do this is to just put constraints on the minimum and maximum number of workers and adjusting those as tight as we can. The downside to that is that we would need to re-run the model and adjust the parameters.

However, I was wondering if there is a way to somehow quantify either the standard deviation (preferable) or the range (less preferable but still appreciated) of the number of workers per shift such that the model remains a MILP?

A couple of options come to mind. Let $$w_s$$ be a variable representing the number of workers during shift $$s.$$ You can introduce nonnegative variables $$y$$ and $$z$$ to represent the minimum and maximum of $$w_s,$$ along with the constraints $$y \le w_s \le z$$ for all shifts $$s \in S.$$ Now add a penalty term $$\lambda(z-y)$$ to your objective ($$\lambda > 0$$ if your objective is minimization, $$\lambda < 0$$ if maximization). Because you are penalizing the size of the range, the optimal solution will automatically set $$y$$ respectively $$z$$ to the smallest respectively largest values of $$w_s.$$ The virtue of this approach is that it keeps your model linear (assuming it was linear before).

Using the standard deviation would be problematic, but if you are willing to "trade-up" to a quadratic integer program you can penalize the variance of the shift sizes. Introduce a nonnegative variable $$z$$ to represent the mean shift size, together with the constraint $$\frac{1}{\vert S \vert}\sum_{s\in S}w_s = z.$$ Add the penalty term $$\lambda \sum_s \left(w_s - z\right)^2$$ to the objective function. Assuming the constraints were linear before the change, they still are, and many (most?) modern IP solvers can handle minimizing a convex quadratic objective function.

• Excellent, practical options. (+1) Mar 24 at 16:34
• That's very helpful, thank you! Would it be possible to adjust the first formulation to target a specific value V? I was thinking to just add constraints z > V and y < V, or by minimizing z-V and V-y, but I'm interested if there is a smarter way.
– Dano
Mar 24 at 18:49
• @Dano I'm not sure what you have in mind. If you want to penalize the amount that the range exceeds $V,$ you can do it by adding yet another variable $u \ge 0$ along with the constraint $z-y+u \le V.$ You would then penalize $u$ rather than $z-y.$
– prubin
Mar 24 at 20:24
• I don't think standard deviation would be problematic, presuming an MISOCP capable solver (and if applicable, modeling system) were used. For example, Gurobi, Xpress, CPLEX, SCIP, COPT, among others. The standard deviation can be written as a 2-norm. If a standard form having linear objective is needed, epigraph formulation can be used which introduces a new variable t in the objective in lieu of the standard deviation, and add the Second Order Cone constraint: norm(...) $\le$ t, where norm(...) is the standard deviation. Mar 25 at 18:45

Similar to above answer, you can avoid quadratic objective by summing up absolute deviation from the mean with:
$$\vert S \vert w_s - \sum_s w_s \le \vert S \vert z$$
$$\sum_s w_s - \vert S \vert w_s \le \vert S \vert z$$
$$0 \le z$$

There's may not be any need to use $$\vert S \vert$$ on RHS

Then minimize $$z$$

You may choose to have variable $$z$$ for each shift $$s$$, then minimize the sum.

• How would one compute the absolute deviation? For example, I've tried using abs(x) in xpress before and it's reverted to a quadratic solver. Aug 30 at 17:31

Why not try to calculate the overall required operators at each shift or any specific period of time, weekly/monthly, as a predefined parameter, in the best case, for each of $$N$$ parts, and incorporating this in the model directly?

This is a very common in many of the practical scheduling software. Suppose for each part you know how long it takes to complete in each machin. By summing up the required time to complete all of parts on the specific machine and also defining your appreciate shift, the number of required operators can be calculated as:

• $$P_{p,m}$$: processing time of part $$p$$ on machin $$m$$.
• $$D_{p,m}$$: the number of parts should be processed on machin $$m$$.
• $$Opr_{p}$$: the number of required operators for each part.

$$Opr_{p} = Ceil((\sum P_{p,m}D_{p,m}) / |S|)$$

Also, if you would like to incorporate the required operators as a decision variable, what proposed by Prof. Rubin, you can use this number as an upper bound on this variable to improve the solving process.