Yet Another Math Programming Consultant

Monday, October 7, 2019

Scipy linear programming: a large but easy LP

Scipy.optimize.linprog [1] recently added a sparse interior point solver [2]. In theory we should be able to solve some larger problems with this solver. However the input format is matrix based. This makes it difficult to express LP models without much tedious programming. Of course if the LP model is very structured things are a bit easier. In [3] the question came up if we can solve some reasonable sized transportation problems with this solver. I claimed this new interior point solver should be able to tackle reasonably large transportation problems. As transportation problems translate into large but easy LPs (very sparse, network structure) this is a good example to try out this solver. It should not require too much programming.

An LP model for the transportation problem can look like:

Transportation Model
\[ \begin{align} \min \> & \sum_{i,j} \color{darkblue}c_{i,j} \color{darkred} x_{i,j} \\ & \sum_j \color{darkred} x_{i,j} \le \color{darkblue}s_i &&\forall i\\ & \sum_i \color{darkred} x_{i,j} \ge \color{darkblue}d_j &&\forall j\\ & \color{darkred}x_{i,j}\ge 0\end{align} \]

Here $i$ indicate the supply nodes and $j$ the demand nodes. The problem is feasible if total demand does not exceed total supply (i.e. $\sum_i s_i \ge \sum_j d_j$).

Even if the transportation problem is dense (that is each supply node can serve all demand nodes or in other words each link $ i \rightarrow j$ exists), the LP matrix is sparse. There are 2 nonzeros per column.

LP Matrix

The documentation mentions we can pass on the LP matrix as a sparse matrix. Here are some estimates of the difference in memory usage:

	100x100	500x500	1000x1000
Source Nodes	100	500	1,000
Destination Nodes	100	500	1,000
LP Variables	10,000	250,000	1,000,000
LP Constraints	200	500	2,000
LP Nonzero Elements	20,000	500,000	2,000,000
Dense Memory Usage (MB)	15	1,907	15,258
Sparse Memory Usage (MB)	0.3	7.6	30.5

For the $1000\times 1000$ case we see that a sparse storage scheme will be about 500 times as efficient.

Sparsity is very important both inside the solver as at the modeling level. Larger but sparser models are often to be preferred above smaller but denser models. Exploiting sparsity will not only save memory but also increase performance: by skipping zero elements we prevent doing all kind of useless work (like multiplying by zero or adding zeros).

Solving a 1000x1000 transportation problem: Implementation

The package scipy.sparse [4] is used to form a sparse matrix. We use three parallel arrays to populate the sparse matrix: one integer array with the row numbers, one integer array with the column numbers and one floating point array with the values. All these arrays have the same length: the number of nonzeros in the LP matrix.
Scipy.optimize.linprog does not allow for $\ge$ constraints. So our model becomes: \[\begin{align} \min &\sum_{i,j} c_{i,j} x_{i,j}\\ & \sum_j x_{i,j} \le s_i &&\forall i \\ & \sum_i -x_{i,j} \le -d_j &&\forall j\\ & x_{i,j}\ge 0\end{align}\]

When I run this, I see:

Primal Feasibility  Dual Feasibility    Duality Gap         Step             Path Parameter      Objective
1.0                 1.0                 1.0                 -                1.0                 4999334.387281
0.01096610265509    0.01096610265504    0.01096610265504    1.0              0.01096610265523    3423127.924532
0.007470719084731   0.007470719084695   0.007470719084695   0.3369198212982  0.007470719084826   1045138.710249
0.007375696439705   0.007375696439669   0.007375696439669   0.01405171378191 0.007375696439798   946062.4541516
0.006900523710037   0.006900523710004   0.006900523710004   0.07151611989327 0.006900523710125   631457.8940984
0.003392688227185   0.003392688227169   0.003392688227169   0.5542765654086  0.003392688227229   106030.5627759
0.002716216726218   0.002716216726205   0.002716216726205   0.2210823772546  0.002716216726252   77660.93708537
0.00151605426328    0.001516054263272   0.001516054263272   0.4706161702772  0.001516054263299   39012.6976106
0.001238382883199   0.001238382883193   0.001238382883193   0.2007381529847  0.001238382883215   31262.77924434
0.0006888763719364  0.000688876371933   0.0006888763719331  0.4711955496918  0.0006888763719452  16884.5788155
0.0004045311601541  0.0004045311601521  0.0004045311601522  0.4504577243574  0.0004045311601593  9812.570668161
0.0003278435563858  0.0003278435563842  0.0003278435563842  0.2062071599936  0.00032784355639    7943.50442653
0.0001938174872602  0.0001938174872593  0.0001938174872593  0.4304958950603  0.0001938174872627  4718.01892459
0.0001272127336263  0.0001272127336257  0.0001272127336257  0.371775562858   0.000127212733628   3126.320160308
7.325610966318e-05  7.325610966282e-05  7.325610966283e-05  0.4526986333113  7.325610966411e-05  1837.061691682
6.047737643405e-05  6.047737643373e-05  6.047737643375e-05  0.1896942778068  6.047737643482e-05  1530.292617672
3.301112106729e-05  3.301112106712e-05  3.301112106713e-05  0.4758440911431  3.301112106771e-05  870.6399411648
2.231615463384e-05  2.231615463375e-05  2.231615463374e-05  0.3562669388094  2.231615463413e-05  613.0954966036
1.300693055479e-05  1.300693055474e-05  1.300693055474e-05  0.4437694284722  1.300693055496e-05  388.3160007487
7.533045251385e-06  7.533045251368e-06  7.533045251357e-06  0.4485635094836  7.533045251489e-06  255.9636413848
3.799832196644e-06  3.799832196622e-06  3.799832196633e-06  0.5264643380152  3.7998321967e-06    165.5742065953
2.01284588862e-06   2.012845888624e-06  2.012845888615e-06  0.5006416028336  2.01284588865e-06   122.2520897954
1.143491145379e-06  1.143491145387e-06  1.143491145377e-06  0.4712206751612  1.143491145397e-06  101.1678772704
5.277850584407e-07  5.277850584393e-07  5.277850584402e-07  0.5711139027487  5.277850584494e-07  86.20125171613
3.125695105059e-07  3.125695105195e-07  3.125695105058e-07  0.4315945026363  3.125695105113e-07  80.96090171621
1.118500099738e-07  1.118500099884e-07  1.118500099743e-07  0.6743118743189  1.118500099763e-07  76.06812425522
4.412565084911e-08  4.412565086951e-08  4.412565085053e-08  0.6297257004579  4.412565085131e-08  74.41374033755
6.833044779903e-09  6.833044770544e-09  6.833044776856e-09  0.8682333453577  6.833044776965e-09  73.50145019804
3.3755500974e-10    3.375549807043e-10  3.375549865145e-10  0.9528386773998  3.375549866004e-10  73.34206256371
1.066148223577e-13  1.065916625724e-13  1.066069704785e-13  0.9998776987355  1.066069928771e-13  73.3337765897
7.763476236577e-18  3.543282811637e-17  5.469419174887e-18  0.9999500035089  5.330350298476e-18  73.3337739491
Optimization terminated successfully.
         Current function value: 73.333774
         Iterations: 30
Filename: transport.py

Line #    Mem usage    Increment   Line Contents
================================================
    59     70.6 MiB     70.6 MiB   @profile
    60                             def run():
    61                                # dimensions
    62     70.6 MiB      0.0 MiB      M = 1000  # sources
    63     70.6 MiB      0.0 MiB      N = 1000 # destinations
    64     78.3 MiB      7.7 MiB      data = GenerateData(M,N)
    65    108.9 MiB     30.5 MiB      lpdata = FormLPData(data)
    66    122.6 MiB     13.7 MiB      res = opt.linprog(c=np.reshape(data['c'],M*N),A_ub=lpdata['A'],b_ub=lpdata['rhs'],options={'sparse':True, 'disp':True})

This proves we can actually solve a $1000 \times 1000$ transportation problem (leading to an LP with a million variables) using standard Python tools.

Notes

It is noted that a dense transportation problem (with all links $i \rightarrow j$ allowed) produces a sparse LP. It is also possible the transportation problem is sparse: only some links are allowed. Sparse transportation problems are a little bit more difficult to set up: the LP matrix is less structured, so we need more advanced data structures (probably a dict to establish a mapping from each existing link $i \rightarrow j$ to a column number). A good modeling tool may help here.

References

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html
https://docs.scipy.org/doc/scipy/reference/optimize.linprog-interior-point.html
Maximum number of decision variables in scipy linear programming module in Python, https://stackoverflow.com/questions/57579147/maximum-number-of-decision-variables-in-scipy-linear-programming-module-in-pytho
https://docs.scipy.org/doc/scipy/reference/sparse.html

Wednesday, October 2, 2019

Running a MIP solver on Raspberry Pi

Raspberry Pi

The Raspberry Pi [1] is a small single board computer. It comes with an ARM based CPU (64 bit, quad core). You can buy it for $35 (no case included). The 4GB RAM version retailes for $55. Raspberry Pi runs some form of Linux. It is mainly used for educational purposes.

SCIP

SCIP [2] is a solver for MIP (and related) models. It is only easily available for academics, using a somewhat non-standard license. As a result, I don't see it used much outside academic circles. So it can not really be called open source.

SCIP on Raspberry Pi

In [3] SCIP is used on the Rapberry Pi with 4GB of RAM. They call it an example of "Edge Computing": bring the algorithm to where it is needed [4] (opposed to moving the data to say a server).

On average SCIP is 3 to 5 times slower on an (standard or overclocked) Rasberry Pi than on a MacBook Pro laptop.

Of course the small amount of RAM means we can only solve relatively small problems. (These days what we call a small MIP problem is actually not so small).

References

Saturday, September 28, 2019

Duplicate constraints in Pyomo model

Introduction

Pyomo [1] is a popular Python based modeling tool. In [2] a question is posed about a situation where a certain constraint takes more than 8 hours to generate. As we shall see, the reason is that extra indices are used.

A simple example

The constraint \[y_i = \sum_j x_{i,j} \>\>\>\forall i,j\] is really malformed. The extra $\forall j$ is problematic. What does this mean? There are two obvious interpretations:

One could say, this is just wrong. Below we will see that GAMS and AMPL follow this approach.
We can also interpret this differently. Assume the inner $j$ is scoped (i.e. local). Then we could read this as: repeat the constraint $y_i = \sum_j x_{i,j}$, $n$ times. Here $n=|J|$ is the cardinality of set $J$. This is what Pyomo is doing.

LaTeX: no checks

Journal papers are often typeset using LaTeX. Not much checking there, and as a result authors are free to write equations like [3]:

Here we have the same thing. We have a summation over $i$ as well as a $\forall i$.

GAMS: syntax error

The GAMS fragment corresponding to this example, shows GAMS will object to this construct:

11 equation e(i,j);

12 e(i,j).. y(i) =e= sum(j, x(i,j));

**** $125

**** 125 Set is under control already

**** 1 ERROR(S) 0 WARNING(S)

AMPL: syntax error

AMPL also gives a syntax error. The message could be improved a bit:

D:\etc>ampl x.mod

d:x.mod, line 8 (offset 155):

syntax error

context: e{i in I, j in J}: Y[i] = sum{j >>> in <<< J} X[i,j];

Pyomo: create duplicate constraints

The Pyomo equivalent can look like:

def eqRule(m,i,j):

return m.Y[i] == sum(m.X[i,j] for j in m.J);

model.Eq = Constraint(model.I,model.J,rule=eqRule)

This fragment is a bit more difficult to read, largely due to syntactic clutter. But in any case: Python and Pyomo accepts this constraint as written. To see what is generated, we can use

model.Eq.pprint()

This will show something like:

Eq : Size=6, Index=Eq_index, Active=True

Key : Lower : Body : Upper : Active

('i1', 'j1') : 0.0 : Y[i1] - (X[i1,j1] + X[i1,j2] + X[i1,j3]) : 0.0 : True

('i1', 'j2') : 0.0 : Y[i1] - (X[i1,j1] + X[i1,j2] + X[i1,j3]) : 0.0 : True

('i1', 'j3') : 0.0 : Y[i1] - (X[i1,j1] + X[i1,j2] + X[i1,j3]) : 0.0 : True

('i2', 'j1') : 0.0 : Y[i2] - (X[i2,j1] + X[i2,j2] + X[i2,j3]) : 0.0 : True

('i2', 'j2') : 0.0 : Y[i2] - (X[i2,j1] + X[i2,j2] + X[i2,j3]) : 0.0 : True

('i2', 'j3') : 0.0 : Y[i2] - (X[i2,j1] + X[i2,j2] + X[i2,j3]) : 0.0 : True

We see for each $i$ we have three duplicate constraints. Instead of 6 constraints, we just should have 2. The way to fix this is to remove the function argument $j$ from eqRule:

def eqRule(m,i):

return m.Y[i] == sum(m.X[i,j] for j in m.J);

model.Eq = Constraint(model.I,rule=eqRule)

After this, model.Eq.pprint() produces

Eq : Size=2, Index=I, Active=True

Key : Lower : Body : Upper : Active

i1 : 0.0 : Y[i1] - (X[i1,j3] + X[i1,j2] + X[i1,j1]) : 0.0 : True

i2 : 0.0 : Y[i2] - (X[i2,j3] + X[i2,j2] + X[i2,j1]) : 0.0 : True

This looks much better.

The original problem

The constraint in the original question was:

def period_capacity_dept(m, e, j, t, dp):

return sum(a[e, j, dp, t]*m.y[e,j,t] for (e,j) in model.EJ)<= K[dp,t] + m.R[t,dp]

model.period_capacity_dept = Constraint(E, J, T, DP, rule=period_capacity_dept)

Using the knowledge of the previous paragraph we know this should really be:

def period_capacity_dept(m, t, dp):

return sum(a[e, j, dp, t]*m.y[e,j,t] for (e,j) in model.EJ)<= K[dp,t] + m.R[t,dp]

model.period_capacity_dept = Constraint(T, DP, rule=period_capacity_dept)

Pyomo mixes mathematical notation with programming. I think that is one of the reasons this bug is more difficult to see. In normal programming, adding an argument to a function has an obvious meaning. However in this case, adding e,j means in effect: $\forall e,j$. If $e$ and $j$ belong to large sets, we can easily create a large number of duplicates.

Notes: Debugging Pyomo models

When debugging Pyomo models there are a few tricks to know about:

Have a small data set available. Debugging with large data sets is very difficult and plain unpleasant.
Use model.symbol.pprint()
Create an LP file of the model. When using the pyomo script, use the convert option (don't forget the --symbolic-solver-labels flag). From Python code, use model.write(filename, io_options={'symbolic_solver_labels': True}).

References

http://www.pyomo.org
Pyomo: Simple inequality constraint takes unreasonable time to build, https://stackoverflow.com/questions/58034244/pyomo-simple-inequality-constraint-takes-unreasonable-time-to-build
Checking math, https://yetanothermathprogrammingconsultant.blogspot.com/2012/03/checking-math.html

Wednesday, September 18, 2019

Demo problem with constraint on standard deviation

In [1] a hypothetical demo problem is shown. I don't think it is a real problem, but rather contrived as an example. Nevertheless, there are things to say about it.

The problem is:

Original Problem
\[\begin{align}\min\>&\sum_i \color{darkred}x_i\\ & \mathbf{sd}(\color{darkred}x) \lt \color{darkblue}\alpha\\ & \color{darkred}x_i \in \{0,1\}\end{align}\]

Notes:

Here sd is the standard deviation.
We assume $x$ has $n$ components.
Of course, $\lt$ is problematic in optimization. So the equation should become a $\le$ constraint.
The standard formula for the standard deviation is: \[ \sqrt{\frac{\sum_i (x_i-\bar{x})^2}{n-1}}\] where $\bar{x}$ is the average of $x$.
This is an easy problem. Just choose $x_i=0$.
When we use $\max \sum_i x_i$ things are equally simple. In that case choose $x_i = 1$.
There is symmetry: $\mathbf{sd}(x) = \mathbf{sd}(1-x)$.
A more interesting problem is to have $\mathbf{sd}(x)\ge\alpha$.

Updated problem

A slightly different problem, and somewhat reformulated is:

MIQCP problem
\[\begin{align}\min\>&\bar{\color{darkred}x} \\ & \bar{\color{darkred}x}= \frac{\sum_i \color{darkred}x_i}{\color{darkblue}n}\\ & \frac{\sum_i (\color{darkred}x_i - \bar{\color{darkred}x})^2}{\color{darkblue}n-1} \ge \color{darkblue}\alpha^2 \\ & \color{darkred}x_i \in \{0,1\}\end{align}\]

First, we replaced $\lt$ by $\ge$ to make the problem more interesting. Furthermore I got rid of the square root. This removes a possible problem with being non-differentiable at zero. The remaining problem is a non-convex quadratically constrained (MIQCP=Mixed Integer Quadratically Constrained Problem). The non-convexity implies we want a global solver.

This model solves easily with solvers like Baron or Couenne.

Integer variable

When we look at the problem a bit more, we see we are not really interested in which $x_i$'s are zero or one. Rather, we need only to worry about the number. Let $k = \sum_i x_i$, Obviously $\bar{x}=k/n$. But more interestingly: \[\mathbf{sd}(x) = \sqrt{\frac{k (1-\bar{x})^2+(n-k)(0-\bar{x})^2}{n-1}}\] The integer variable $k$ is restricted to $k=0,2,\dots,n$.

Thus we can write:

MINLP problem
\[\begin{align}\min\>&\color{darkred}k \\ & \frac{\color{darkred}k (1-\color{darkred}k/\color{darkblue}n)^2+(\color{darkblue}n-\color{darkred}k) (\color{darkred}k/\color{darkblue}n)^2}{\color{darkblue}n-1} \ge \color{darkblue}\alpha^2 \\ & \color{darkred}k = 0,1,\dots,\color{darkblue}n \end{align}\]

The constraint can be simplified into \[\frac{k-k^2/n}{n-1}\ge \alpha^2\] This is now so simple we can do this by enumerating $k=0,\dots,n$, check the constraint, and pick the best.

Because of the form of the standard deviation curve (note the symmetry), we can specialize the enumeration loop and restrict the loop to $k=1,\dots,\lfloor n/2 \rfloor$. Pick the first $k$ that does not violate the constraint (and when found exit the loop). For very large $n$ we can use something like a bisection to speed things up even further.

So this example optimization problem does not really need to use optimization at all.

References

Constrained optimisation with function in the constraint and binary variable, https://stackoverflow.com/questions/57850149/constrained-optimisation-with-function-in-the-constraint-and-binary-variable
Another problem that minimizes the standard deviation, https://yetanothermathprogrammingconsultant.blogspot.com/2017/09/minimizing-standard-deviation.html

Wednesday, August 28, 2019

Octeract

I don't know what or who this is.

This seems to be a parallel deterministic solver for non-convex MINLPs. Some things I noticed:

!np.easy : cute (but nonsense of course: some problems just remain difficult).
"The first massively parallel Deterministic Global Optimization solver for general non-convex MINLPs."
Symbolic manipulation: like some other global solvers they need the symbolic form of the problem so they can reason about this. I.e. no black-box problems.
Support for AMPL and Pyomo
"Octeract Engine has implementations of algorithms that guarantee convergence to a global optimum in finite time. The correctness of the majority of the calculations are ensured through a combination of interval arithmetic and infinite precision arithmetic."
It looks like the benchmarks [3] are against itself (so always a winner).
"Massively parallel" makes me think about thousands of cores. The benchmarks do not seem to report anything like this (mostly 1,4 and 8 cores).
I don't see any names on the web site. The About Company section is unusually vague. I received some more links with additional background info [4,5].

Some of the competing solvers are Baron, Couenne, and Antigone.

References

https://octeract.com/
Manual: https://octeract.com/wp-content/uploads/2019/08/user_manual.pdf
Benchmarks: https://octeract.com/benchmarks/
Presentation: https://www.youtube.com/watch?v=ayZ3iHmQUak
PhD thesis: https://spiral.imperial.ac.uk/bitstream/10044/1/45359/1/Kazazakis-N-2017-PhD-Thesis.pdf

Monday, August 26, 2019

Python-MIP

This is another modeling tool for Python.

There are quite a few modeling tools available for Python: Pyomo, PuLP, and most commercial LP/MIP solvers come with some Python modeling layer.

This is what caught my eye when reading about Python-MIP:

The name is rather unimaginative.
Looks like the authors are from Brazil.
Supported solvers are CBC and Gurobi.
The Python-MIP is compatible with the just-in-time compiler PyPy, which can lead to substantial performance improvements.
It is claimed that with PyPy jit, python-MIP can be 25 times as fast than the Gurobi modeling tool.
There are some interesting facilities supported by Python-MIP:

Cuts can be provided using a call-back mechanism
Support for MIPSTART (initial integer solution)
Solution pool

Interesting question

Related to Python-MIP, in [3] an interesting question came up. The value of an integer variable is often slightly non-integer. E.g. something like 0.0000011625. This is the result of the integer feasibility tolerance that a solver applies.

Background: during the MIP search quite a lot of integer variables are automatically (or by accident) integer valued. To exploit this, and not wasting time branching on these variables, a MIP solver wants to consider variables that are close to being integer as integer. Basically, skipping these variables when branching.

In the discussion [3] the remark is made:

Is rounding integer solutions automatically a good idea? No MIP solver or modeling system is doing this AFAIK The exceptions are modeling tools mainly targeting Constraint Programming solvers. Here are some arguments against this rounding strategy.

Rounding integer solutions can lead to larger infeasibilities. With some aggressive presolve/scaling these infeasibilities can sometimes be large after postsolve/unscaling. Also: some equations may have long summations of binary variables. This would accumulate a lot of rounding errors. And then there are these big-M constraints....

It also means that using the steps:

solve
fix solution (or fix integers) to optimal solution
solve

may lead to "feasible" for the first solve, but "infeasible" for the second solve. E.g. when we want duals for the fixed LP we use this "fix integers" step.

Safer would be to tighten the integer feasibility tolerance. Cplex even allows epint=0 (epint is Cplex's integer feasibility tolerance). Of course tightening the integer feasibility tolerance will likely lead to longer solution times.

Indeed, modeling systems and solvers currently handle this by offloading the problem to the user. The solver is probably the right place to deal with this. Not sure developers are eager to work on this. Taking responsibility by simple-minded rounding may be asking for more problems than it solves.

On the other hand, these slightly fractional values certainly cause confusion. Especially for beginners in MIP modeling. This is the main argument in favor of automatic rounding.

So the question remains: is rounding integer variables a good idea?

References

Wednesday, August 21, 2019

R/Python + C++

In some recent projects, I was working on using algorithms implemented in C++ from R and Python. Basically the idea is: Python and R are great languages for scripting but they are slow as molasses. So, it may make sense to develop the time consuming algorithms in C++ while driving the algorithm from R or Python.

R and C++: Rcpp.

The standard way to build interfaces between R and C++ code is to use Rcpp [1,2].

It is possible to directly interface R with low-level C code, This will require a lot of code and knowledge of R internals. Rcpp automates a lot of this. E.g. Rcpp will take care of translating an R vector into a C++ vector.

Rcpp supports both small fragments of C++ code passed on as an R string to more coarse-grained file based approach [3]. For Windows, you need to download the GNU compilers [4].

If you are new to both Rcpp and to building your own R packages [5] things may be a bit overwhelming.

Rstudio can help a lot. It supports a lot of very useful tasks:

Syntax coloring for C++ code.
Building projects.
Git version control.
Documentation tools (rmarkdown and bookdown). My documentation is also a good test: it executes almost all of the code when building the document.

Editing C++ code in RStudio

Basically I never have to leave RStudio.

I have added an alternative driver file for my C++ code so I can debug it in Visual Studio. I used it only a few times: most of the time I just used RStudio.

Python and C++: pybind11

pybind11 [6] is in many respects similar to Rcpp, although it requires a little bit more programming to bridge the gap between Python and C++.

In the beginning of the above Youtube video [7], the presenter compares pybind11 with some of the alternatives:

SWIG: author of SWIG says: don't use it
ctypes: call c functions but not c++
CFFI: call c functions
Boost.Python: support for older C++ standards, but not much maintained
pybind11: modern

As with Rcpp, calling the compiler is done through running a build or setup script. For Rcpp I used the GNU compilers, while pybind11/pip install supports the Visual Studio C++ compiler. This also means that if you have little experience with both pybind11 and creating packages, the learning curve may be steep.

References

http://www.rcpp.org
Dirk Eddelbuettel, Seamless R and C++ Integration with Rcpp, Springer, 2013
Chapter: Rewriting R code in C++, https://adv-r.hadley.nz/rcpp.html
Hadley Wickham, R packages, O'Reilly, 2015
https://cran.r-project.org/bin/windows/Rtools/
https://pybind11.readthedocs.io/en/master/
Robert Smallshire, Integrate Python and C++ with pybind11, https://www.youtube.com/watch?v=YReJ3pSnNDo

Thursday, August 15, 2019

Scheduling teachers

From [1]:

I am trying to find a solution to the following problem:

There are 6 teachers

Each teacher can only work for 8 hours a day

Each teacher must have a 30 minute break ideally near the middle of his shift

The following table shows the number of teacher needed in the room at certain times:

7am - 8am 1 teacher

8am - 9am 2 teacher

10am - 11am 5 teacher

11am - 5pm 6 teacher

5pm - 6pm 2 teacher

What is a good way to solve this (ideally in python and Google OR-Tools) ?

Thank you

Initial analysis

From the demand data we see that we need all 6 teachers from 11am-5pm working without lunch break. This is 6 hours. So there is no way to allow a lunch break near the middle of the shift.

We have more problems. The picture of the demand data indicates we have no demand between 9am and 10am. That does not look right. Bad data is fact of life. Optimization models are especially sensitive to bad data.

I believe that looking critically at your data is essential for successful optimization applications. You can learn a lot from just a bit of staring.

Alternative problem 1

We can ask a few different questions. If we only allow shifts of the form: 4 hours work, 0.5 hour lunch, 4 hours work (i.e. lunch perfectly in the middle of the shift), how many teachers do we need to meet the demand. To model this, we can assume time periods of half an hour. The number of possible shifts is small:

With an enumeration of the shifts, we can model this as a covering problem:

Covering Model
\[\begin{align} \min\> &\color{darkblue}z= \sum_s \color{darkred} x_s \\ & \sum_{s\|\color{darkblue}{\mathit cover}(s,t)} \color{darkred} x_s \ge \color{darkblue}{\mathit demand}_{t} && \forall t \\ &\color{darkred} x_i \in \{0,1,2,\dots\}\end{align} \]

Here $\mathit{cover}(s,t)=\text{True}$ if shift $s$ covers time period $t$. If we would interpret $\mathit cover(s,t)$ as a binary (data) matrix \[\mathit cover_{s,t} = \begin{cases} 1 & \text{if shift $s$ covers time period $t$}\\ 0 & \text{otherwise}\end{cases}\] we can also write:

Covering Model (alternative interpretation)
\[\begin{align} \min\> &\color{darkblue}z= \sum_s \color{darkred} x_s \\ & \sum_s \color{darkblue}{\mathit cover}_{s,t} \cdot \color{darkred} x_s \ge \color{darkblue}{\mathit demand}_{t} && \forall t \\ &\color{darkred} x_i \in \{0,1,2,\dots\}\end{align} \]

Our new assumptions are:

we only allow shifts with a lunch break in the middle
we add to our demand data: 4 teachers needed between 9am and 10am

After solving this easy MIP model we see:

----     55 VARIABLE x.L  number of shifts needed

shift1 2.000,    shift4 2.000,    shift5 2.000,    shift6 2.000


----     55 VARIABLE z.L                   =        8.000  total number of shifts

I.e. we need 8 teachers to handle this workload.

The picture shows we are overshooting demand in quite a few periods. Note: this picture illustrates the left-hand side (orange bars) and right-hand side (blue line) of the demand equation in the Covering Model.

Alternative problem 2

Lets make things a bit more complicated. We allow now the following rules for a single shift:

The work period before the lunch bread is between 3 and 5 hours (or between 6 and 10 time periods)
There is a lunch break of 0.5 or 1 hour.
After lunch there is another work period of between 3 and 5 hours. The total number of working hours is 8.

When we enumerate the shifts, we see:

We now have 55 different shifts to consider.

The results look like:

----     72 VARIABLE x.L  number of shifts needed

shift1  1.000,    shift22 1.000,    shift26 1.000,    shift39 1.000,    shift49 1.000,    shift52 1.000
shift55 1.000


----     72 VARIABLE z.L                   =        7.000  total number of shifts

We see the number of teachers needed is 7. We are closer to the demand curve:

We see that if we add more flexibility we can do a bit better. Achieving 6 teachers is almost impossible. We would need to introduce shifts like: work 2 hours, lunch, work 6 hours. The teachers union would object.

It is noted there are quite a few other methods to solve models like this where we don't need to enumerate all possible shifts [2,3]. For larger problems it may not be feasible to employ the shift enumeration scheme we used here.

References

Constraint programming scheduling problem, https://stackoverflow.com/questions/57473085/constraint-programming-scheduling-problem
Employee Scheduling II: Column Generation, http://yetanothermathprogrammingconsultant.blogspot.com/2017/01/employee-scheduling-ii-column-generation.html
Employee Scheduling IV: Direct Optimization, http://yetanothermathprogrammingconsultant.blogspot.com/2017/03/employee-scheduling-iv-direct.html

Wednesday, August 7, 2019

Finding the central point in a point cloud

Problem

The problem is easy to describe: given the coordinates of points $p_i$, find the central point $x$ that minimizes the sum of the distances from $x$ to $p_i$ [1]. This point is sometimes called the geometric median [2].

Data

We just use some random data on the unit square $[0,1]\times[0,1]$:

----     13 PARAMETER p  points (randomly generated)

                  x           y

point1   0.17174713  0.84326671
point2   0.55037536  0.30113790
point3   0.29221212  0.22405287
point4   0.34983050  0.85627035
point5   0.06711372  0.50021067
point6   0.99811763  0.57873338
point7   0.99113304  0.76225047
point8   0.13069248  0.63971876
point9   0.15951786  0.25008053
point10  0.66892861  0.43535638
point11  0.35970027  0.35144137
point12  0.13149159  0.15010179
point13  0.58911365  0.83089281
point14  0.23081574  0.66573446
point15  0.77585761  0.30365848
point16  0.11049229  0.50238487
point17  0.16017276  0.87246231
point18  0.26511455  0.28581432
point19  0.59395592  0.72271907
point20  0.62824868  0.46379787
point21  0.41330699  0.11769536
point22  0.31421227  0.04655151
point23  0.33855027  0.18209959
point24  0.64572713  0.56074555
point25  0.76996172  0.29780586
point26  0.66110626  0.75582167
point27  0.62744750  0.28386420
point28  0.08642462  0.10251467
point29  0.64125115  0.54530950
point30  0.03152485  0.79236064
point31  0.07276700  0.17566105
point32  0.52563261  0.75020767
point33  0.17812371  0.03414099
point34  0.58513117  0.62122998
point35  0.38936190  0.35871415
point36  0.24303462  0.24642154
point37  0.13050280  0.93344972
point38  0.37993791  0.78340046
point39  0.30003426  0.12548322
point40  0.74887411  0.06923246
point41  0.20201556  0.00506586
point42  0.26961305  0.49985148
point43  0.15128587  0.17416945
point44  0.33063773  0.31690605
point45  0.32208696  0.96397664
point46  0.99360221  0.36990306
point47  0.37288857  0.77197833
point48  0.39668414  0.91309632
point49  0.11957773  0.73547889
point50  0.05541847  0.57629980
mean     0.39042709  0.47299042

The last entry labeled mean is just the average of the $x$- and $y$-coordinates.

NLP Model

A straightforward non-linear programming model can look like:

Unconstrained NLP Model
\[\min \sum_i \sqrt{\sum_c (\color{darkred}x_c-\color{darkblue}p_{i,c})^2 }\]

We use $c = \{x,y\}$, i.e. we have $x$ and $y$-coordinates. Note that we use $x$ in two different contexts: element of set $c$, being the $x$-coordinate, and the decision variable $x_c$.

We can use the mean as a very good starting point to help the NLP solver. I.e. \[x_c := \frac{\displaystyle\sum_{i=1}^n p_{i,c}}{n}\]

The picture below shows why the mean is such a good starting point:

Optimal center point is close to mean point

The numeric values are here:

----     45 PARAMETER results  x(center) vs mean

               x           y     sumdist

mean  0.39042709  0.47299042 18.02414861
x     0.37048298  0.43857572 17.96891985

The sumdist column shows the objective values for these two points.

This is an easy NLP problem. Most NLP solvers just need a few iterations. With a system like GAMS or AMPL we get exact gradients automatically. That is much preferable to finite differences which seems the prevalent method people use in an R or Python environment.

It was pointed out to me that the problem is not smooth: the gradient of $\sqrt{x}$ is not defined at zero, This can only happen when our point $x$ exactly coincides with a data point (i.e. all coordinates are exactly the same). We can protect against this by using: \[\min \sum_i \sqrt{\varepsilon+\sum_c (x_c-p_{i,c})^2 }\] for a very small constant $\varepsilon\gt 0$.

Cone programming I

The above problem can also be written as a cone programming problem. This will allow us to use a different class of solvers to work on this problem. Here we use CVXPY [3] to express the model. The Python code can look like:

import cvxpy as cp  

x = cp.Variable(2)   # center point

obj = cp.Minimize( cp.sum( [ cp.norm(x-p[i,:]) for i in range(N) ] ) )
prob = cp.Problem(obj)

objval = prob.solve(solver=cp.SCS, verbose=True)

This is very much a straight translation of our unconstrained NLP model. Although we only declared two $x$ variables, behind the scenes the model is blown up to a rather large one. We can see from the log:

----------------------------------------------------------------------------
    SCS v2.1.1 - Splitting Conic Solver
    (c) Brendan O'Donoghue, Stanford University, 2012
----------------------------------------------------------------------------
Lin-sys: sparse-direct, nnz in A = 150
eps = 1.00e-04, alpha = 1.50, max_iters = 5000, normalize = 1, scale = 1.00
acceleration_lookback = 0, rho_x = 1.00e-03
Variables n = 52, constraints m = 150
Cones:  soc vars: 150, soc blks: 50

The generated SOCP (second order cone programming) model is larger, but also very sparse. The solver has no problem solving it very quickly. With SOCP solvers we usually don't worry about an initial point like we used in the NLP model.

Cone programming II

If we don't use a modeling tool that can do these transformations automatically, we can use a DIY approach. Second-order cone constraints can be stated as: \[||A_i x+b_i||_2 \le c_i^Tx + d_i \>\>\forall i\] This would imply that we can write our model as:

SOCP Model Attempt
\[\begin{align} \min & \sum_i \color{darkred} d_i \\ & \color{darkred} d_i^2 \ge \sum_c (\color{darkred}x_c-\color{darkblue}p_{i,c})^2 && \forall i \\ &\color{darkred} d_i \ge 0, \color{darkred} x_c \text{ free}\end{align} \]

Unfortunately this will yield a message like: CPLEX Error 5002: 'e(point1)' is not convex or Constraint 'e(point1)'(0) is not convex. Q should be positive semidefinite in a constraint with finite upper bound. We can repair this as follows:

Repaired SOCP Model
\[\begin{align} \min & \sum_i \color{darkred} d_i \\ & \color{darkred} d_i^2 \ge \sum_c \color{darkred}y_{i,c}^2 && \forall i \\ & \color{darkred} y_{i,c} = \color{darkred} x_c -\color{darkblue}p_{i,c} && \forall i,c \\ &\color{darkred} d_i \ge 0, \color{darkred} x_c \text{ free}, \color{darkred} y_{i,c} \text{ free}\end{align} \]

This now solves quickly. We can now understand that CVXPY did quite a few steps before passing the model to the solver. As argued in [4], it is much better if the modeling system takes care of these reformulations. Some of them are not immediate obvious, and hand-crafted reformulations can be error-prone.

The NLP model has just 2 variables $x_c$. This SOCP model has 152 variables. Using $n=50$ data points, we added $3n$ variables. When looking at the number of constraints we see a similar thing. The NLP model has no constraints, but this SOCP model has $3n=150$ constraints ($n$ of them cones).

Conclusion

The min sum distance problem has a simple NLP formulation which can be improved by using a good initial point. It can also be formulated as a SOCP problem. Using high-level modeling tools this is not difficult. Without automatic reformulations things become a bit less obvious.

References

The point that minimizes the sum of euclidean distances to a set of n points, https://stackoverflow.com/questions/57277247/the-point-that-minimizes-the-sum-of-euclidean-distances-to-a-set-of-n-points
Geometric median, https://en.wikipedia.org/wiki/Geometric_median
https://www.cvxpy.org/
Victor Zverovich, Robert Fourer, Automatic Reformulation of Second-Order Cone Programming Problems, https://ampl.com/MEETINGS/TALKS/2015_01_Richmond_2E.2.pdf