Employee Scheduling II: column generation

In (1) and (2) an Employee Scheduling model is considered where all possible shifts are enumerated before we solve the actual model. Although such models models tend to solve quickly (especially considering their size), generating all possible shifts may not be an option for larger problems. One possible way to attack this problem is using a traditional column generation approach. This method will generate attractive candidate solutions and add them to the problem.

Our original problem was:

\[\begin{align}\min\>&\sum_s c_s x_s\\ &\sum_{s|E(e,s)} x_s \le 1 \>\forall e\\ &\sum_{s|S(s,t)} x_s \ge R_t \>\forall t\\ &x_s \in\{0,1\}\end{align}\]

where

\[x_s = \begin{cases}1 & \text{if shift $s$ is selected}\\  0 & \text{otherwise} \end{cases}\]

Instead of generating all possible shifts in advance, we start with a small initial set \(s \in S\). We want this initial set to form a feasible solution for the problem. To make this easier we rewrite the problem as:

\[\begin{align} \min\> & \sum_s c_s x_s + \sum_t c_t s_t\\ &\sum_{s|E(e,s)} x_s \le 1\>\forall e\\ &\sum_{s|S(s,t)} x_s + s_t \ge R_t\>\forall t \\ &x_s \in\{0,1\}\\ & s_t \ge 0\end{align}\]

Here \(s_t\) is a slack with a penalty \(c_t\). We also can consider this as extra, flexible, hourly personnel to hire at a hourly rate of \(c_t\) to cover staffing shortages. Sometimes we call this scheme: making the model “elastic”.

The column generation algorithm iterates between two different problems. First we have a master problem that is essentially a relaxed version of the elastic model we just discussed. We replace the condition \(x_s \in \{0,1\}\) by \(x_s \in [0,1]\). Using this we now have an LP (and thus we can generate duals). Second, there is a sub problem that calculates a candidate column (or candidate shift). 

  

The sub problem uses a particular objective: it minimizes the reduced cost of the new LP column \(A_j\) wrt to the master problem. I.e. the objective is basically

\[\min \sigma_j = c_j – \pi^T A_j\]

where \(\pi\) are the LP duals of the master problem.

This scheme is called column generation as we increase the size of the number of variables \(x_s\) during the algorithm:

Note that the columns \(s_t\) stay fixed.

We stop iterating as soon as we can no longer find a new column with \(\sigma_j < 0\), indicating a “profitable” column. We now have collection of shifts and an LP solution. I.e. we have fractional variable values. To make this a usable schedule we need to solve one more MIP: introduce again the condition \(x_s \in\{0,1\}\). The complete algorithm can look like:

1. Generate an initial (small) set \(S\) of feasible shifts
2. Solve the relaxed master problem (LP)
3. Solve the sub problem
4. If \(\sigma_j < 0\) add shift to the set \(S\) and go to step 2
5. Solve master problem as MIP

Initial Set

For the problem in (1) we generate 100 random feasible shifts. The first 10 shifts are as follows:

Note that we have here two shifts for employee White. This is no problem: the master problem will select the most suitable one to include in the schedule. All these shifts are feasible: they obey the minimum and maximum hours and will not use hours an employee is not available. If the shifts are not sufficient to cover all staffing requirements, no problem: the LP model will start hiring the expensive hourly “emergency” workers. Indeed for these 100 random shifts we need some extra hourly help:

----    214 VARIABLE short.L  staffing shortage in time t t5  1.000,    t6  1.000,    t7  1.000,    t9  2.000,    t13 1.000,    t15 1.000,    t16 1.000,    t17 1.000 t18 1.000,    t19 1.000,    t23 1.000

Master Problem

The master problem expressed in GAMS looks like:

set    s0 'shifts' /s1*s1000/    s(s0) 'dynamic subset will grow during algorithm'    esmap(e,s0) 'mapping between employee and shift'    shifts(s0,t) 'data structire to hold shifts' ; parameters    a(t,s0) 'matrix growing in dimension s'    c(s0)   'cost vector growing in dimension s' ; scalar    shortcost 'cost to hire extra hourly staff' /200/ ; binary variable x(s0) 'main variable: shift is selected'; positive variable short(t) 'staffing shortage in time t'; variable cost; equations   masterObj        'calculate total cost'   oneShift(e)      'each employee can do at most one shift'   coverPeriod(t)   'cover each period' ; masterObj..    cost =e= sum(s, c(s)*x(s)) + sum(t,shortcost*short(t)); oneShift(e)..    sum(esmap(e,s), x(s)) =l= 1; coverPeriod(t)..    sum(shifts(s,t), x(s)) + short(t) =g= requirements(t); model Master /masterObj,oneShift,coverPeriod/;

Sub problem

The sub problem is somewhat more difficult. We need to find a feasible shift and minimize the reduced cost \(\sigma_j\).

set sh 'feasible shift hours' /sh2*sh8/; parameter shlen(sh); shlen(sh) = ord(sh)+1; binary variables    w(t) 'work at time t'    start(t) 'start of shift'    emp(e) 'select shift from employee e'    shft(sh) 'select feasible shift length'    se(sh,e) 'select (sh,e) combination' ; variables    shiftlen  'length of shift'    wagecost    redcost   'reduced cost' ; equations    startShift(t)    'start of shift: w(t) goes from 0 to 1'    oneShiftStart    'exactly one start of shift'    calcLen          'calculate length of shift'    oneEmployee      'pick one employee'    minLength        'observe minimum shift length'    maxLength        'observe maximum shift length'    feasibleShift(t) 'make sure shift is feasible'    calcShift1       'select shift length'    calcShift2       'id'    secomb(sh,e)     'shift/employee combination'    calcwage         'calculate cost of shift'    subObj           'objective' ; startShift(t)..  start(t) =g= w(t) - w(t--1); oneShiftStart..  sum(t,start(t)) =e= 1; calcLen..        shiftLen =e= sum(t, w(t)); oneEmployee..    sum(e, emp(e)) =e= 1; minLength..      shiftLen =g= sum(e, minlen(e)*emp(e)); maxLength..      shiftLen =l= sum(e, maxlen(e)*emp(e)); feasibleShift(t).. w(t) =l= sum(canwork(e,t),emp(e)); calcShift1..     sum(sh,shft(sh)) =e= 1; calcShift2..     sum(sh,shlen(sh)*shft(sh)) =e= shiftlen; secomb(sh,e)..   se(sh,e) =g= shft(sh) + emp(e) - 1; calcwage..       wageCost =e= sum((sh,e),shlen(sh)*wage(e)*se(sh,e)); subObj..   redcost =e= wageCost - sum(e,emp(e)*oneShift.m(e)) - sum(t, w(t)*coverPeriod.m(t)); model sub /subObj,startShift,oneShiftStart,calcLen,oneEmployee,minLength,            maxLength,feasibleShift,calcShift1,calcShift2,secomb,calcWage/;

 

Some of the modeling issues with the sub problem are explained in more detail in (3).

Column Generation Algorithm

The CG algorithm itself can look like:

sets   iter /iter1*iter100/ ; scalars    continue /1/ ; parameter    masterObjVal(iter) ; loop(iter$continue,    solve master using rmip minimizing cost;    masterObjVal(iter) = cost.L;    solve sub using mip minimizing redcost;    if (redcost.l < -0.001,       esmap(e,nexts)$(emp.l(e)>0.5) = yes;       shifts(nexts,t)$(w.l(t)>0.5) = yes;       c(nexts) = wagecost.L;       s(nexts) = yes;       nexts(s0) = nexts(s0-1);    else      continue = 0;    ); ); solve master using mip minimizing cost;

Results

We need to do 36 cycles before we can conclude there are no more interesting columns. The objective of the relaxed master problem decreases to an LP objective representing a cost of $4670. Of course the LP objective is non-increasing: adding a column (i.e. adding a shift) can only help.

After these 36 iterations we have collected 135 shifts (remember we starting in the first iteration with our 100 randomly generated initial shifts). This compares quite favorably to the 1295 shifts we obtained by enumerating all possible shifts. When we reintroduce the integer restrictions \(x_s\in \{0,1\}\) the objective value deteriorates to $4675. This is just a little bit worse than the global optimal objective of $4670. Indeed this method does not guarantee to find absolute best integer solution, however often it finds very good solutions.

The final schedule looks like:

We don’t need any expensive hourly “emergency” staff, i.e. \(s_t = 0\). Again we have no slack in the system: we precisely meet the staffing requirements every hour:

 

References

1. Loren Shure, Generating an Optimal Employee Work Schedule Using Integer Linear Programming, January 6, 2016, http://blogs.mathworks.com/loren/2016/01/06/generating-an-optimal-employee-work-schedule-using-integer-linear-programming/ describes the problem with a Matlab implementation.
2. Employee Scheduling I: Matlab vs GAMS, http://yetanothermathprogrammingconsultant.blogspot.com/2017/01/employee-scheduling-i-matlab-vs-gams.html has a GAMS version of the model in (1).
3. Employee Scheduling III: generating all shifts using the Cplex solution pool, http://yetanothermathprogrammingconsultant.blogspot.com/2017/01/employee-scheduling-iii-generating-all.html
4. Erwin Kalvelagen, Column Generation with GAMS, http://www.amsterdamoptimization.com/pdf/colgen.pdf has some background on Column Generation in GAMS.
5. Employee Scheduling IV: direct optimization, http://yetanothermathprogrammingconsultant.blogspot.com/2017/03/employee-scheduling-iv-direct.html

Employee Scheduling I: Matlab vs GAMS

In (1) a small employee scheduling problem is modeled as a simple MIP model using MATLAB. Basically from the data on employee availability, we need to form all possible shifts. Even for small data sets that can become a large number. Having all these possible shifts we introduce a binary variable:  

\[x_s = \begin{cases}1 & \text{if shift $s$ is selected}\\                               0 & \text{otherwise} \end{cases}\]

We have some requirements with respect to staffing levels, and an employee can only do one shift in the 24 hours we are look at. The standard model for this is:

\[\begin{align}\min\>&\sum_s c_s x_s\\                             &\sum_{s|E(e,s)} x_s \le 1\>\forall e\\                             &\sum_{s|S(s,t)} x_s \ge R_t \> \forall t\\                             &x_s \in\{0,1\}\end{align}\]

The first constraint says each employee can perform at most one shift. The second constraint makes sure that we have enough shifts selected such that we meet staffing requirements. The objective minimizes the cost associated with the selected shifts. This is a very simple model, so actually the main issues are with data manipulation: we need to form a collection of all possible shifts.

Data

We reproduce here the data from (1):

Data entry using Excel is always very convenient. Most users know and have access to Excel. Notice that the last column  (Availability) is a string, which needs to be parsed inside the application.  Matlab has sscanf, GAMS has no similar facilities. If you want to directly read this into GAMS you could use 4 columns, say Start1, End1, Start2, End2 with integer values.

The staffing requirements are:

Matlab code

The Matlab code from (1) and (2) combines generating the shifts and producing the LP matrix. This makes it more difficult to understand the model: it is hidden inside the code. I quickly scanned the code, but decided not to translate this into GAMS but rather start from scratch. Translating code line-by-line is often not as useful.

The blog post (1) displays the non-zero pattern of the LP matrix. In my opinion, in a GAMS environment such a picture is not very useful.

GAMS code: generating shifts

The GAMS code I wrote separates generation of shifts from the model equations. As it turns out the first task is not completely trivial.

The first thing I did was to generate a set canwork(e,t) containing whether a time slot is allowed by an employee:

This was produced from the rules:

canwork(e,t)$(data(e,'start1')=0 and data(e,'end1')=0) = yes; canwork(e,t)$(ord(t)>data(e,'start1') and ord(t)<=data(e,'end1')) = yes; canwork(e,t)$(ord(t)>data(e,'start2') and ord(t)<=data(e,'end2')) = yes;

The next thing I did was to generate a parameter maxstart(e,t) indicating how long the maximum shift can be starting in period t. This looks something like:

It is not simple to compute. Note that near the end of the day we still can start a long shift: shifts wrap from t24 to t1. Here is my code:

set back(t,tt) 'used to loop backwards over time'; back(t,tt)$(ord(tt)+ord(t)=card(t)+1) = yes; * when looping over this set tt runs backwards parameters    count0(e,t)    'max length of remaining part of shift in t1'    maxstart(e,t)  'max length of shift starting in t' ; * * First calculate maximum length of remaining part of shift in t1 * This is used when wrapping from t24 to t1. * Notes: *   We scan backwards. *   count0(e,t) is indexed by t, but we only use t=t1. * count0(e,t) = 0; loop(back(t,tt),    count0(e,'t1')$canwork(e,tt) = min(count0(e,'t1')+1,maxlen(e));    count0(e,'t1')$(not canwork(e,tt)) = 0; ); * * Now calculate maximum length of shift starting in t * maxstart(e,t) = 0; loop(back(t,tt),    maxstart(e,tt)$canwork(e,tt) = min(maxstart(e,tt+1)+1+count0(e,t),maxlen(e)); ); maxstart(e,tt)$(maxstart(e,tt)<minlen(e)) = 0;

From this we can compute the following:

1. The set shifts(s,t) with a shift id s and hour t.
2. A mapping set esmap(e,s), that maps between employee e and shift s.
3. A subset s of used shifts.

The data structure shifts looks like:

We can see here just a part of the shifts. There are 1295 of them. Near the bottom we see a few shifts wrapping from t24 to t1.

Shifts s1 through s24 belongs to employee Smith; Johnson owns shifts s25-s96. This correspondence is stored in the mapping set esmap(e,s). The code to generate this is:

sets   s0 'shifts'  /s1*s10000/   s(s0) 'used shifts'   k(s0) /s1/ ; scalar i; loop((e,t)$maxstart(e,t),    for(i = minlen(e) to maxstart(e,t),       shifts(k,tt)$(ord(tt)>=ord(t) and ord(tt)<ord(t)+i) = yes;       shifts(k,tt)$(ord(tt)>=1 and ord(tt)<ord(t)-card(t)+i) = yes;       esmap(e,k) = yes;       s(k) = yes;       k(s0) = k(s0-1);    ); );

Boy this was not straightforward, and the code is not very intuitive. But we managed. Probably using a standard programming language this shift generation step may be easier.

A very different technique is to use an optimization model in combination with the Cplex solution pool technology to generate all possible shifts. In (4) this approach is applied to this problem.

GAMS code: MIP model

Finally we need to implement optimization model. This is actually very easy now:

parameter shiftlen(s0) 'number of hours in a shift'; shiftlen(s) = sum(shifts(s,t),1); binary variable x(s0); variable cost; equations    oneShift(e)     'each employee can do at most one shift'    coverPeriod(t)  'cover each period'    calcCost        'calculate total cost' ; oneShift(e)..    sum(esmap(e,s), x(s)) =l= 1; coverPeriod(t)..    sum(shifts(s,t), x(s)) =g= requirements(t); calcCost..   cost =e= sum(esmap(e,s),  wage(e)*shiftlen(s)*x(s)); model m /all/; option optcr=0; solve m minimizing cost using mip;

This model follows very closely the mathematical model presented at the beginning of this post. This model solves very fast. Cplex can find the optimal solution (proving optimality) without doing any branching. The problem with this formulation is that it can generate a ton of possible shifts, but an advantage is that it often solves very fast. The results are:

This is different from what is shown in (1) but with the same objective: the cost is $4670. As in (1) there is no slack: we precisely match the staffing requirements:

References

1. Loren Shure, Generating an Optimal Employee Work Schedule Using Integer Linear Programming, January 6, 2016, http://blogs.mathworks.com/loren/2016/01/06/generating-an-optimal-employee-work-schedule-using-integer-linear-programming/
2. http://blogs.mathworks.com/images/loren/2016/makeMILPMatrices.m contains the matrix generation code.
3. Employee Scheduling II : Column Generation, http://yetanothermathprogrammingconsultant.blogspot.com/2017/01/employee-scheduling-ii-column-generation.html uses a Column Generation approach to solve this problem.
4. Employee Scheduling III : generating all shifts using the Cplex solution pool, http://yetanothermathprogrammingconsultant.blogspot.com/2017/01/employee-scheduling-iii-generating-all.html uses the Cplex Solution Pool technology to generate all possible shifts. This is an alternative approach to the method shown in this post.
5. Employee Scheduling IV: direct optimization, http://yetanothermathprogrammingconsultant.blogspot.com/2017/03/employee-scheduling-iv-direct.html

MIP Starts sometimes disappoint

Modern MIP solvers allow to use a good integer solution to “warm start” the branch and bound process. This looks like a very good approach to cut down on total solution times, and it sometimes is. But in many cases it turns out, finding good solutions and passing them on to the MIP solver is not worth the while.

Some reasons can be:

• Solvers on their own are very good in finding good solutions quickly: they contain batteries of heuristics to find good integer solutions and improve them.
• Time you spend on calculating a good initial solution can be used by the solver. If you spent n seconds to construct an initial solution, that is time taken away from the solver.
• The solver spends most time in the tail: getting the gap down to zero (or close).

This is probably best illustrated as follows:

The blue line represents the currently best integer solution. I.e. we are minimizing here. When we provide a good integer solution we basically eat away some time from the left of the picture. But that is exactly the part where MIP solvers are very effective: the blue line decreases there very fast (the gradient with respect to time is large). In the picture above if we would have provided a starting solution with an objective value of 1000, we would essentially cut out the shaded area. The long tail at the right is where we spend most cycles.

The red line is the lower bound: the best possible objective. When we skip the shaded box we also miss out on some improvement in this bound.

It is clear that in this case, if we cannot find good solutions very cheaply ourselves, we may as well leave this to the solver.

Of course there are cases where starting with good solutions may make good sense, e.g.:

• We can obtain good solutions quickly.
• We may even have a feasible operating point by nature of the problem (e.g. when optimizing a certain process we may observe the physical state and use that operating point as a starting point).
• Good solutions found outside the solver may improve reliability of the overall process: if the MIP solver fails completely we still have a good solution.

Update

As suggested in the comments stopping on a gap percentage may be more effective in reducing the solution time:



Note that this is a depiction of quite a typical behavior of MIP lower and upper bounds: almost all non-trivial models show this type of narrowing of the gap. Often things are even more pronounced than in this smaller example. Of course stopping on gap does not likely give you the optimal solution, but in practical cases a gap of say a few percent or so may be quite acceptable. Setting a gap tolerance is a very standard approach when solving large MIP models.

New Book: R for Data Science; Visualize, Model, Transform, Tidy and Import Data

Recommended introductory text. It discusses some packages I was not really familiar with (like tibble).

Contents

Part I. Explore

1. Data Visualization with ggplot2
2. Workflow: Basics
3. Data Transformation with dplyr
4. Workflow: Scripts
5. Exploratory Data Analysis
6. Workflow: Projects

Part II. Wrangle

7. Tibbles with tibble
8. Data Import with readr
9. Tidy with tidyr
10. Relational Data with dplyr
11. Strings with stringr
12. Factors with forcats
13. Dates and Times with lubridate

Part III. Program

14. Pipes with magrittr
15. Functions
16. Vectors
17. Iteration with purrr

Part IV. Model

18. Model Basics with modelr
19. Model Building
20. Many Models with purrr and broom

Part V. Communicate

21. R Markdown
22. Graphics for Communication with ggplot2
23. R Markdown Formats
24. R Markdown Workflow

References

1. Hadley Wickham, Garrett Grolemund, R for Data Science: Import, Tidy, Transform, Visualize, and Model Data, O'Reilly Media; 1 edition (January 5, 2017), 522 pages, price: about $0.08 per page (glossy, with many color pictures) 
2. Website: http://r4ds.had.co.nz/

Solving LP by LU???

Yes, according to this paper [link]:

Abstract This paper presents a new approach for the solution of Linear Programming Problems with the help of LU Factorization Method of matrices. This method is based on the fact that a square matrix can be factorized into the product of unit lower triangular matrix and upper triangular matrix. In this method, we get direct solution without iteration. We also show that this method is better than simplex method.

Basically the paper claims that we can solve an LP

\[\begin{align}\max\>& c^Tx\\ &Ax \le b \\ & x\ge 0\end{align}\]

as follows:

1. Make sure \(A\) is square (by adding rows or column if needed using vacuous constraints or extra slacks).
2. Solve \(Ax=b\) using LU decomposition
3. Done

Wow. This is baloney of course.  But here is someone actually believing this stuff.

Update

Of course we can use this nonsensical method and apply it to a different problem:

Abstract. In this paper the linear fractional programming problem is converted to a linear programming problem and for solving it, LU Factorization method is introduced. Finally, to illustrate the proposed method, it is verified by means of the numerical example and the results conclude that this method gives more reliable solution or the equal solution compared to the existing method.

My goodness.

References

1. S.M.Chincole, A.P.Bhadane, LU Factorization Method To Solve Linear Programming Problem, International Journal of Emerging Technology and Advanced Engineering ,Volume 4, Issue 4, 176-180, http://www.ijetae.com/files/Volume4Issue4/IJETAE_0414_31.pdf. Warning: this is nonsensical. Note that the example has all constraints binding in the optimal solution.
2. R. Sophia Porchelvi, A. Anna Sheela, Solving Linear Fractional Programming Problem Using LU Factorization Method, Asian Journal of Research in Social Sciences and Humanities Vol. 6, No. 10, October 2016, pp. 1929-1938, https://www.aijsh.com/shop/articlepdf/2016/10/1475416103154.pdf. Warning squared: astonishingly this is a paper based on (1).

Solving Takuzu puzzles as a MIP using xor

The Takuzu puzzle (1) can be solved using a MIP model. Interestingly we can use an xor condition in the model.

Assuming an  \(n \times n\) board (with \(n\) even), we need to fill the cells \((i,j)\), with values 0 or 1, with the following restrictions:

1. The number of ones in each row and in each column is equal to \(\displaystyle\frac{n}{2}\).
2. Then we have also: The number of zeroes in each row and in each column is equal to \(\displaystyle\frac{n}{2}\).
3. There can be no duplicate rows and no duplicate columns.
4. In a row or column there cannot be more than two consecutive zeroes or ones.

A small example is given in (1):

Of course our central variable will be \(x_{i,j}\in\{0,1\}\). Condition 1 can be easily modeled using two constraints:

\[\begin{align}   &\sum_j x_{i,j} = \frac{n}{2}\>\> \forall i\\   &\sum_i x_{i,j} = \frac{n}{2}\>\> \forall j \end{align}\]

Condition 2 we don’t need to model explicitly: this follows from the above constraints.

Condition 3 (no duplicate rows or columns) is not very easy. One way to handle this, is to calculate a value for each row and column as follows:

\[\begin{align} &v_i = \sum_j 2^{j-1} x_{i,j}\\ &w_j = \sum_i 2^{i-1} x_{i,j} \end{align}\]

and then make sure there are no duplicate \(v_i\)’s or duplicate \(w_j\)’s. This can be done using some all-different formulation (3). A big-M approach to implement \(v_i \le v_{i’}-1 \textbf{ or } v_i \ge v_{i’}+1, \forall i\ne i’\) can look like:

\[\begin{align} &v_i \le v_{i’}-1 + M \delta_{i,i’}\\ &v_i \ge v_{i’}+1 – M (1-\delta_{i,i’})\\ &\delta_{i,i’} \in \{0,1\} \end{align}\]

However I want to attack this in a different way, using an xor condition (4). The constraints involved can be conceptualized as:

\[\begin{align} &\sum_j |x_{i,j} – x_{i’,j}| \ge 1 \>\> \forall i\ne i’\\ &\sum_i |x_{i,j} – x_{i,j’}| \ge 1 \>\> \forall j\ne j’ \end{align}\]

The value \(z_{i,j,i’,j’}=|x_{i,j}-x_{i’,j’}|\) can be interpreted as \(z_{i,j,i’,j’}=x_{i,j} \textbf{ xor } x_{i’,j’}\): zero if \(x_{i,j}\) and \(x_{i’,j’}\) are identical and one if they are different. This in turn can be written as a set of inequalities (4):

\[\begin{align} &z_{i,j,i’,j’}\le x_{i,j} + x_{i’,j’}\\ &z_{i,j,i’,j’}\ge x_{i,j} - x_{i’,j’}\\ &z_{i,j,i’,j’}\ge x_{i’,j’} - x_{i,j}\\ &z_{i,j,i’,j’}\le 2 - x_{i,j} - x_{i’,j’}\\ &z_{i,j,i’,j’} \in \{0,1\} \end{align}\]

After this we just need:

\[\begin{align} &\sum_j z_{i,j,i’,j} \ge 1 \>\> \forall i\ne i’\\ &\sum_i z_{i,j,i,j’} \ge 1 \>\> \forall j\ne j’ \end{align}\]

Due to the nature of these conditions on \(z\) we can limit the xor constraints to

\[\begin{align} &z_{i,j,i’,j’}\le x_{i,j} + x_{i’,j’}\\ &z_{i,j,i’,j’}\le 2 - x_{i,j} - x_{i’,j’}\\ &z_{i,j,i’,j’} \in \{0,1\} \end{align}\]

Note that we do not need to compute all \(z\) for all combinations \((i,j)\), \((i’,j’)\) (e.g. we never compare \(x_{2,3}\) with \(x_{3,2}\)). Furthermore we can relax \(z\) to be continuous between zero and one. Also note that if we compared row \(i\) and \(i’\) then we don’t need to compare row \(i’\) and row \(i\). Similarly for the columns. So we can do:

\[\begin{align} &\sum_j z_{i,j,i’,j} \ge 1 \>\> \forall i< i’\\ &\sum_i z_{i,j,i,j’} \ge 1 \>\> \forall j< j’ \end{align}\]

Finally we need to handle restriction 4. This can be implemented as:

\[\begin{align} &1\le x_{i,j}+x_{i,j+1}+x_{i,j+2}\le 2\\ &1\le x_{i,j}+x_{i+1,j}+x_{i+2,j}\le 2 \end{align}\]

Usually we need to write this as four inequalities (unless we can specify ranged equations). If the summation was longer I would have introduced extra variables with bounds one and two to reduce the number of constraints.

Note that the lower bound of one can also be inferred from \((1-x_{i,j})+(1-x_{i,j+1})+(1-x_{i,j+2})\le 2 \Rightarrow x_{i,j}+x_{i+1,j}+x_{i+2,j} \ge 1\).

Complete model

The variables are:

The dummy variable is related to a dummy objective: we look for a feasible integer solution only.

The different \(z\) values we need to compare is encoded in a set compare:

Finally the model equations are:

The result looks like:

----     77 VARIABLE x.L              i1          i2          i3          i4 i1                       1           1 i2           1                                   1 i3                                   1           1 i4           1           1

A large problem

The site (2) contains some “very hard” 14 x 14 puzzles. An example is:

It is interesting to see how this model behaves on this problem. It turns out the problem can be solved completely in preprocessing. E.g. Cplex and Gurobi report that zero Simplex iterations were needed and solve this problem blazingly fast:

Tried aggregator 3 times. MIP Presolve eliminated 5491 rows and 2301 columns. MIP Presolve modified 18 coefficients. Aggregator did 396 substitutions. Reduced MIP has 92 rows, 48 columns, and 298 nonzeros. Reduced MIP has 48 binaries, 0 generals, 0 SOSs, and 0 indicators. Presolve time = 0.03 sec. (7.10 ticks) Found incumbent of value 0.000000 after 0.05 sec. (8.15 ticks) Root node processing (before b&c):   Real time             =    0.05 sec. (8.24 ticks) Sequential b&c:   Real time             =    0.00 sec. (0.00 ticks)                           ------------ Total (root+branch&cut) =    0.05 sec. (8.24 ticks) MIP status(101): integer optimal solution Cplex Time: 0.05sec (det. 8.25 ticks)

Verifying Uniqueness

A well-formed puzzle has only one solution. To verify whether a solution is unique we can use a similar strategy as in (6):

1. First solve the model to obtain a solution \(x^*_{i,j}\).
2. Formulate a constraint to forbid this solution.
3. Add this constraint to the model and resolve. The solver should not find a new solution but rather report that this model is infeasible.

Here the cut can be formulated as:

\[\sum_{i,j} x^*_{i,j} x_{i,j} \le \frac{n^2}{2}-1\]

Data Entry

Entering the problem data:

is an interesting issue. In GAMS “zero” and “does not exist” is identical, so we can not just use:

table init(i,j)
    i1 i2 i3 i4
 i1     1     0
 i2        0
 i3     0
 i4  1  1     0
;

I used the following, which is not very clean:

table init(i,j)
    i1 i2 i3 i4
 i1     1     2
 i2        2
 i3     2
 i4  1  1     2
;
x.fx(i,j)$(init(i,j)=1) = 1;
x.fx(i,j)$(init(i,j)=2) = 0;

Not sure if this is better:

parameter init(i,j) / #i.#j NA /;
$onmulti
table init(i,j)
    i1 i2 i3 i4
 i1     1     0
 i2        0
 i3     0
 i4  1  1     0
;
variable x(i,j);
x.fx(i,j)$(init(i,j)<>NA) = init(i,j);

Note that AMPL has a "default" facility, so we can use in a data file:

param init default -1 : i1 i2 i3 i4 :=
                   i1    .  1  .  0
                   i2    .  .  0  .
                   i3    .  0  .  .
                   i4    1  1  .  0
;

It makes sense to use Excel for data entry. The initial configuration in Excel can look like:

To read this into GAMS and retain the zeros we can use the following code:

parameter init(i,j);

$call 'gdxxrw i=takuzu.xlsx sq=n par=init rng=b3 rdim=1 cdim=1'

$oneps
$gdxin takuzu.gdx
$loaddc init

display init;

The Squeeze=no option in the call to GDXXRW will tell it to retain the zeros. The $OnEps command tells GAMS to read the zero as an EPS special value. EPS is very useful: it is physically stored in the sparse data structures, but as soon as you do a numerical operation on it, it becomes zero. This will show:

----     12 PARAMETER init              i1          i2          i3          i4 i1                   1.000                     EPS i2                                 EPS i3                     EPS i4       1.000       1.000                     EPS

This will now allow us to fix variables as:

binary variable x(i,j);
x.fx(i,j)$init(i,j) = init(i,j);

This model turned out to have quite a few interesting angles.

References

1. Takuzu, https://en.wikipedia.org/wiki/Takuzu
2. Online Takuzu puzzles, http://www.binarypuzzle.com/
3. All-different and Mixed Integer Programming, http://yetanothermathprogrammingconsultant.blogspot.com/2016/05/all-different-and-mixed-integer.html
4. xor as Inequalities, http://yetanothermathprogrammingconsultant.blogspot.com/2016/02/xor-as-linear-inequalities.html
5. Similar puzzles are Sudoku and KenKen: http://yetanothermathprogrammingconsultant.blogspot.com/2016/10/mip-modeling-from-sudoku-to-kenken.html
6. Unique Solutions in KenKen, http://yetanothermathprogrammingconsultant.blogspot.com/2016/12/unique-solutions-in-kenken.html