There are two major implications for the model itself:
- We can use integers directly, i.e. \(x_{i,j} \in \{1,\dots,u_{i,j}\}\) instead of \(x_{i,j,k} \in \{0,1\}\).
- We have access to constructs like Distinct and != (unequal).
However most of the code is devoted to data manipulation as we will see below. Although the approaches are different (different languages, different solver technologies), the actual code is remarkably similar.
Data entry: GAMS
First we need to get the data into GAMS. Unfortunately we need to enter two matrices: one for the blocks and one for the initial values.
| set i 'rows' /i1*i7/ j 'columns' /j1*j9/ k 'values' /1*9/ b 'blocks' /b1*b15/ ; table blockno(i,j) j1 j2 j3 j4 j5 j6 j7 j8 j9 i1 1 1 2 2 2 3 3 4 6 i2 1 1 2 2 3 3 4 4 6 i3 1 7 8 8 3 4 4 5 6 i4 7 7 7 8 8 9 9 6 6 i5 7 10 10 8 12 9 9 14 15 i6 11 11 10 10 12 12 9 15 15 i7 11 11 10 12 12 13 13 15 15 ; table initval(i,j) j1 j2 j3 j4 j5 j6 j7 j8 j9 i1 1 1 i2 2 3 i3 4 4 1 4 i4 1 i5 4 5 4 5 1 i6 3 i7 2 5 ; |
The set k is a bit of an over-estimate: we actually need only values 1 through 5 for this problem (the maximum size of a block is 5).
Data entry: Python
Entering the same boards in Python can look like:
from z3 import *
blockno = [ [ 1, 1, 2, 2, 2, 3, 3, 4, 6], [ 1, 1, 2, 2, 3, 3, 4, 4, 6], [ 1, 7, 8, 8, 3, 4, 4, 5, 6], [ 7, 7, 7, 8, 8, 9, 9, 6, 6], [ 7,10,10, 8,12, 9, 9,14,15], [11,11,10,10,12,12, 9,15,15], [11,11,10,12,12,13,13,15,15] ] initval = [ [ 1, 0, 0, 1, 0, 0, 0, 0, 0], [ 0, 2, 0, 0, 0, 0, 0, 3, 0], [ 4, 0, 4, 0, 0, 0, 0, 1, 4], [ 0, 0, 0, 1, 0, 0, 0, 0, 0], [ 4, 0, 5, 0, 4, 0, 5, 1, 0], [ 0, 0, 0, 0, 0, 0, 0, 3, 0], [ 0, 0, 0, 0, 2, 0, 0, 0, 5] ] m = len(blockno) # number of rows n = len(blockno[0]) # number of columns R = range(m) C = range(n) bmax = max([blockno[i][j] for i in R for j in C]) B = range(bmax) |
The order is a bit reversed here: first enter the boards with the block numbers and the initial values and then calculate the dimensions.
Derived data: GAMS
We need to calculate a number of things before we can actually write the equations.
alias (i,i2),(j,j2); sets bmap(b,i,j) 'maps between b <-> (i,j)' bk(b,k) 'allowed (b,k) combinations' ijk(i,j,k) 'allowed values' n(i,j,i2,j2) 'neighbors' ; parameters len(b) 'number of cells in block' ; bmap(b,i,j) = blockno(i,j)=ord(b); len(b) = sum(bmap(b,i,j),1); bk(b,k) = ord(k)<=len(b); ijk(i,j,k) = sum((bmap(b,i,j),bk(b,k)),1); n(i,j,i,j+1) = blockno(i,j)<>blockno(i,j+1); n(i,j,i+1,j-1) = blockno(i,j)<>blockno(i+1,j-1); n(i,j,i+1,j) = blockno(i,j)<>blockno(i+1,j); n(i,j,i+1,j+1) = blockno(i,j)<>blockno(i+1,j+1); |
Note that the set n (indicating neighbors we need to inspect for different values) only looks “forward” from cell \((i,j)\). This is to prevent double-counting: we don’t want to compare neighboring cells twice in the equations.
Derived data: Python
The Python/Z3 model needs to calculate similar derived data:
# derived data
maxb = [0]*bmax # number of cells in each block bij = [[] for b in B] # list of cells in each block nij = [ [ [] for j in C ] for i in R ] # neighbors for i in R: for j in C: bn = blockno[i][j]-1 maxb[bn] += 1 bij[bn] += [(i,j)] if j+1 nij[i][j] += [(i,j+1)] if i+1 nij[i][j] += [(i+1,j-1)] if i+1 nij[i][j] += [(i+1,j)] if i+1 nij[i][j] += [(i+1,j+1)] # upperbound on X[i,j] maxij = [ [ maxb[blockno[i][j]-1] for j in C ] for i in R ] |
I use here nested lists. bij has a list of cells for each block. nij contains a list of neighboring cells for each cell.
Model constraints: GAMS
The decision variables and linear constraints are described in [1]. The GAMS representation is very close the mathematical formulation.
variable dummy 'objective variable';
binary variable x(i,j,k); * fix initial values x.fx(i,j,k)$(initval(i,j)=ord(k)) = 1; equation oneval(i,j) 'pick one value per cell' unique(b,k) 'unique in a block' neighbors(i,j,i2,j2,k) 'neighboring cells must be different' edummy 'dummy objective' ; oneval(i,j).. sum(ijk(i,j,k), x(i,j,k)) =e= 1; unique(bk(b,k)).. sum(bmap(b,i,j), x(i,j,k)) =e= 1; neighbors(i,j,i2,j2,k)$(n(i,j,i2,j2) and ijk(i,j,k) and ijk(i2,j2,k)).. x(i,j,k)+x(i2,j2,k) =l= 1; edummy.. dummy =e= 0; |
Notice we added a dummy objective to the model.
Model constraints: Z3
# Model
X = [ [ Int("x%d.%d" % (i+1, j+1)) for j in C ] for i in R ] Xlo = And([X[i][j] >= 1 for i in R for j in C]) Xup = And([X[i][j] <= maxij[i][j] for i in R for j in C]) Uniq = And([ Distinct([X[i][j] for (i,j) in bij[b]]) for b in B]) Nbor = And([ And([X[i][j] != X[i2][j2] for (i2,j2) in nij[i][j] ]) for i in R for j in C if len(nij[i][j])>0 ]) Xfix = And([X[i][j] == initval[i][j] for i in R for j in C if initval[i][j]>0]) Cons = [Xlo,Xup,Uniq,Nbor,Xfix] |
Most of this is boilerplate. We use Distinct to make the contents of the cells unique within a block. The != operator is used to make sure neighbors of different blocks do not have the same value. Both these constructs are not so easy to implement in a MIP model, but there we were rescued by using binary variables instead.
Solving and printing results: GAMS
We only look for a feasible solution (any feasible integer solution will be globally optimal immediately). This means we don’t have to worry about setting the OPTCR option in GAMS. The default gap of 10% is irrelevant here as the solver will stop after finding the first integer solution.
model m /all/;
solve m minimizing dummy using mip; parameter v(i,j) 'value of each cell'; v(i,j) = sum(ijk(i,j,k), x.l(i,j,k)*ord(k)); option v:0; display v; |
To report a solution we recover the integer cell values. The results look like:
---- 81 PARAMETER v value of each cell j1 j2 j3 j4 j5 j6 j7 j8 j9 i1 1 5 4 1 5 1 5 2 1 |
Solving and printing results: Python/Z3
# solve and
print
s = Solver() s.add(Cons) if s.check() == sat: m = s.model() sol = [ [ m.evaluate(X[i][j]) for j in C] for i in R] print(sol) |
The solution looks like:
[[1, 5, 4, 1, 5, 1, 5, 2, 1], [3, 2, 3, 2, 3, 2, 4, 3, 5], [4, 1, 4, 5, 4, 1, 5, 1, 4], |
Check for uniqueness: GAMS
To check if the solution is unique we add a cut and resolve. If this is infeasible we know there was only one solution.
*
* test if solution is unique * equation cut 'forbid current solution'; cut.. sum(ijk(i,j,k), x.l(i,j,k)*x(i,j,k)) =l= card(i)*card(j)-1; model m2 /all/; solve m2 minimizing dummy using mip; abort$(m2.modelstat=1 or m2.modelstat=8) "Unexpected solution found",x.l; |
We abort if the model status is not "optimal" (1) or "integer solution found" (8). (It should be enough to check on optimal only, we err on the safe side here).
Check for uniqueness: Z3
# Check for
uniqueness
Cut = Or([X[i][j] != sol[i][j] for i in R for j in C]) s.add(Cut) if s.check() == sat: print("There is an alternative solution.") else: print("Solution is unique.") |
No comments:
Post a Comment