Thursday, October 19, 2017

Pareto fronts in Z3: a large data set

In [1] a small multi-objective MIP is solved with Z3. Here we try a larger problem, an engineering design problem. It has three objectives and 2029 non-dominated solutions. The solution in the F space (the objectives) looks like:

image

image

To make life easier for Z3, I converted the objectives to integer values. The good news is that Z3 finds the same number of solutions. The bad news: for this problem it is not very fast:

image

We see it is getting slower and slower:

image

The first 500 solutions are found in a rate of about a second per solution. Later on this becomes 10 seconds per solution.

References
  1. http://yetanothermathprogrammingconsultant.blogspot.com/2017/10/multi-objective-mixed-integer.html

Wednesday, October 11, 2017

Multi-objective Mixed Integer Programming and Z3


Vilfredo Federico Damaso Pareto (1848-1923) [link]

The tool \(\nu Z\) [1,2] is an optimization extension to the SMT solver Z3 [3]. It can actually generating Pareto optimal solutions for MIP problems, using an algorithm based on [6]. In [4,5] a small example model is presented:

\[\bbox[lightcyan,10px,border:3px solid darkblue]{ \begin{align}\max\>&x_1 – 2x_2\\
\max\>&-x_1+3x_2\\
&x_1 \le 2x_2\\
&x_1,x_2 \in \{0,1,2\}\end{align}}\]

There are 5 different Pareto optimal solutions. Let’s see how we can implement this in Z3 using Python:

image

This looks promising. I have some large instances lying around. I am curious how this solver will behave on those.

A real surprise

We have to be very careful when using floating point coefficients. If we divide the first objective by 2, we should get the same Pareto front. However what we see is:

image

image

What happened here? No, error messages, no warnings, but obviously the results are wrong. We can see what happened by printing:

image

Our 0.5 has become 0! (I am yelling now: keep your hands off my data!).

After some experimentation, I came up with the following workaround:

image

Now we get:

image

The results for z1 indicate the solver actually reports rational values. Indeed the printed model has:

image

(A bit funny. I would expect 0.5 or (/1 2) but not this mixture.)

References
  1. Nikolaj Bjørner, Anh-Dung Phan, and Lars Fleckenstein, \(\nu Z\) - An Optimizing SMT Solver, https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/nbjorner-nuz.pdf
  2. Nikolaj Bjørner and Anh-Dung Phan, \(\nu Z\) - Maximal Satisfaction with Z3, in Proceedings of International Symposium on Symbolic Computation in Software Science (SCSS), 2014, https://easychair.org/publications/open/xbn
  3. The Z3 Theorem Prover, https://github.com/Z3Prover/z3 
  4. Generating all non-dominated solutions in a multi-objective integer programming model, http://yetanothermathprogrammingconsultant.blogspot.com/2010/02/generating-all-non-dominated-solutions.html
  5. Sylva, Crema: A method for finding the set of non-dominated vectors for multiple objective integer linear programs, EJOR 158 (2004), 46-55.
  6. D. Rayside, H.-C. Estler, and D. Jackson, The Guided Improvement Algorithm for Exact, General-Purpose, Many-Objective Combinatorial Optimization, Technical Report MIT-CSAIL-TR-2009-033, MIT, 2009, https://dspace.mit.edu/bitstream/handle/1721.1/46322/MIT-CSAIL-TR-2009-033.pdf

Saturday, October 7, 2017

Modern Distributed Optimization

Summary
Matt Adereth talks about the Black-box optimization techniques, what’s actually going on inside of these black-boxes and discusses an idea of how they can be used to solve problems today. He deep dives into a few of the most popular ones, such as Distributed Nelder-Mead and Bayesian Optimization, and discusses their trade-offs.

The video of this presentation can be found in [1]. The title is a bit misleading: the talk is really only about two algorithms for unconstrained derivative-free optimization (a small subset of the total optimization field).

For reference, see the excellent review paper [2].

References
  1. Modern Distributed Optimization, https://www.infoq.com/presentations/black-box-optimization
  2. Rios, L.M. & Sahinidis, N.V., Derivative-free optimization: a review of algorithms and comparison of software implementations, J Glob Optim (2013) 56: 1247, https://doi.org/10.1007/s10898-012-9951-y

Friday, October 6, 2017

von Koch and MiniZinc

image

A MiniZinc [1] version of the graceful labeling model in [3] was sent to me [4]. This looks actually quite nice, so therefore I reproduce it here:

include "globals.mzn";
 
%
% Definition of data
%
 
% Number of nodes in the graph
int: nodes;
set of int: Nodes = 1..nodes;
 
% Connections
int: connections;
set of int: Connections = 1..connections;
array[Connections,1..2] of int: ConnectionEdges;
 
%
% Instance from http://yetanothermathprogrammingconsultant.blogspot.se/2017/09/von-koch-and-z3.html
% The following data could be located in a minizinc data file instead.
%
nodes = 7;
connections = 6;
ConnectionEdges =
         [|  1,  2,
          |  1,  4,
          |  1,  7,
          |  2,  3,
          |  2,  5,
          |  5,  6,
          |];
 
%
% Model proper
%
 
array[Nodes] of var 1..nodes: NodeLabel;
array[Connections] of var 1..connections: ConnectionLabel;
 
constraint alldifferent(NodeLabel) :: domain;
 
constraint alldifferent(ConnectionLabel) :: domain;
 
constraint forall(connection in Connections) (
    ConnectionLabel[connection] = abs(NodeLabel[ConnectionEdges[connection, 1]] - NodeLabel[ConnectionEdges[connection, 2]])
);

 
solve satisfy;
 
output [
  "NodeLabels: " ++ show(NodeLabel) ++
  " ConnectionLabels: " ++ show(ConnectionLabel) ++
  "\n"
  ];

MiniZinc comes with an IDE so that makes it easy to edit and run the model:

image

This is solved with the Gecode [2] Constraint Programming solver.

Generating all solutions is just a question of changing a setting:  

image

Indeed we get all solutions:

image

References
  1. http://www.minizinc.org/
  2. http://www.gecode.org/
  3. von Koch and Z3, http://yetanothermathprogrammingconsultant.blogspot.com/2017/09/von-koch-and-z3.html
  4. Mikael Zayenz Lagerkvist, private communication

Tuesday, October 3, 2017

Modeling environments for the 21st century?

image

Is this a prediction for the period through December 31, 2100?

Anyway, some interesting points are brought forward. A comparison is made between GAMS, AMPL, AIMMS, solver specific C++ code, Pyomo and JuMP:

image

Note: my GAMS models are typically >90% data manipulation and <10% model equations (does X means this is not possible or allowed?).

References
  1. http://egon.cheme.cmu.edu/ewo/docs/EWO_Seminar_03_10_2017.pdf

Saturday, September 30, 2017

Numberlink Models

In the Numberlink puzzle [1] we need to connect end-points by paths over a board with cells such that:

  1. Each cell is part of a path
  2. Paths do not cross

E.g. the puzzle on the left has a solution as depicted on the right [1]:

 

The picture on the left provides the end-points. The solution in the picture on the right can be interpreted as:

image

The main variable we use is:

\[x_{p,k} = \begin{cases}1 & \text{ if cell $p$ has value $k$}\\
                                       0 & \text{ otherwise}\end{cases}\]

The main thought behind solving this model as a math programming model is to look at the neighbors of a cell. The neighbors are the cells immediately on the left, on the right, below or above. If a cell is an end-point of a path then it has exactly one neighbor with the same value. If a cell is not an end-point it is on a path between end-points and has two neighboring cells with the same value.

image

This leads to the following high-level MIQCP (Mixed Integer Quadratically Constrained Model):

\[\bbox[lightcyan,10px,border:3px solid darkblue] {\begin{align}
                      &\sum_k x_{p,k} = 1\\
                      &\sum_{q|N(p,q)} \sum_k x_{p,k}\cdot x_{q,k} =\begin{cases} 1&\text{if cell $p$ is an end-point}\\ 2&\text{otherwise}\end{cases}\\
&x_{p,k} = c_{p,k} \>\>\> \text{if cell $p$ is an end-point}\\&x_{p,k}\in \{0,1\}\end{align}}\]

Here \(N(p,q)\) is true if \(p\) and \(q\) are neighbors and \(c_{p,k}\) indicates if cell \(p\) is an end-point with value \(k\). The model has no objective: we are just looking for a feasible integer solution. Unfortunately, this is a non-convex model. The global solver Baron can solve this during its initial local search heuristic phase:

  Doing local search
  Solving bounding LP
  Starting multi-start local search
  Done with local search
===========================================================================
    Iteration    Open nodes         Time (s)    Lower bound      Upper bound
  *         1             0            23.09      0.00000          0.00000   
            1             0            23.11      0.00000          0.00000   

  Cleaning up

                         *** Normal completion ***           

  Wall clock time:                    23.73
  Total CPU time used:                23.12

  Total no. of BaR iterations:       1
  Best solution found at node:       1
  Max. no. of nodes in memory:       1

  All done

Of course we can try to linearize this model so it becomes a MIP model (and thus convex). The product of two binary variables \(y=x_1\cdot x_2\) can be linearized by:

\[\begin{align}
&y \le x_1\\
&y\le x_2\\
&y \ge x_1+x_2-1\\
&y \in \{0,1\}\end{align}\]

We can relax \(y\) to a continuous variable between 0 and 1. In our model we can exploit some symmetry. This reduces the number of variables and equations:

\[\begin{align}
&y_{p,q,k} \le x_{p,k}&\forall p<q\\
&y_{p,q,k} \le x_{p,k}&\forall p<q\\
&y_{p,q,k} \ge x_{p,k}+x_{q,k}-1&\forall p<q\\
&y_{p,q,k} \in \{0,1\}&\forall p<q\end{align}\]

The non-linear counting equation can now be written as:

\[\sum_{q|N(p,q)}  \sum_k \left( y_{p,q,k}|p<q + y_{q,p,k}|q<p \right)= b_{p}\]

Here \(b_p\) is the right-hand side of the counting constraint. This solves very fast. Cplex shows:

Tried aggregator 2 times.
MIP Presolve eliminated 6178 rows and 1840 columns.
MIP Presolve modified 40 coefficients.
Aggregator did 394 substitutions.
Reduced MIP has 11167 rows, 3892 columns, and 26610 nonzeros.
Reduced MIP has 3892 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.19 sec. (42.60 ticks)
Found incumbent of value 0.000000 after 0.36 sec. (79.41 ticks)

Root node processing (before b&c):
   Real time             =    0.36 sec. (79.70 ticks)
Sequential b&c:
   Real time             =    0.00 sec. (0.00 ticks)
                           ------------
Total (root+branch&cut) =    0.36 sec. (79.70 ticks)
MIP status(101): integer optimal solution
Cplex Time: 0.36sec (det. 79.72 ticks)

Cplex solves this in the root node.

Unique solution

A well-formed puzzle has exactly one solution. We can check this by adding a constraint that forbids the solution we just found. After adding this cut the model should be infeasible. Let \(\alpha_{p,k} = x^*_{p,k}\) be our feasible solution. Then our cut can look like:

\[\sum_{p,k} \alpha_{p,k} x_{p,k} \le |P|-1\]

where \(|P|\) is the number of cells. Indeed, when we add this constraint, the model becomes infeasible. So, we can conclude this is a proper puzzle.

Alternative formulation

In the comments Rob Pratt suggests a formulation that works a better in practice. The counting constraint is replaced by:

\[\begin{align}
&\sum_{q|N(p,q)} x_{q,k}=1 &&\text{if cell $p$ is an end-point with value $k$}\\
&2x_{p,k}\le \sum_{q|N(p,q)} x_{q,k} \le 2x_{p,k}+M_p(1-x_{p,k})&&\text{if cell $p$ is not an end-point}
\end{align}\]

where \(M_p\) can be replaced by \(M_p=|N(p,q)|=\sum_{q|N(p,q)} 1\) (the number of neighbors of cell \(p\)). This is at most 4 so not really a big big-M.

This formulation seems to work better. We can solve with this something like:

image

image

This problem has also exactly one solution.

Historical note

In [3] an early version of this puzzle is mentioned, appearing in

image

image

Here the goal is to draw paths from each house to its gate without the paths overlapping [4].

Large problems

For larger problems MIP solvers are not the best tool. From notes in [2] it looks like SAT solvers are most suitable for this type problem. In a CP/SAT framework we could use directly an integer variable \(x_p \in \{1,\dots,K\}\), and write the counting constraints as something like:

\[\sum_{q|N(p,q)} 1|(x_p=x_q) = \begin{cases} 1&\text{if cell $p$ is an end-point}\\ 2&\text{otherwise}\end{cases}\]
References
  1. Numberlink, https://en.wikipedia.org/wiki/Numberlink
  2. Hakan Kjellerstrand, Numberlink puzzle in Picat, https://github.com/hakank/hakank/blob/master/picat/numberlink.pi (shows a formulation in the Picat language)
  3. Aaron Adcock, Erik D. Demaine, Martin L. Demaine, Michael P. O'Brien, Felix Reidl, Fernando Sanchez Villaamil, and Blair D. Sullivan, Zig-Zag Numberlink is NP-Complete, 2014, https://arxiv.org/pdf/1410.5845.pdf
  4. Facsimiles of The Brooklyn Daily Eagle,  http://bklyn.newspapers.com/image/50475607/, https://bklyn.newspapers.com/image/50475838/
  5. Neng-Fa Zhou, Masato Tsuru, Eitaku Nobuyama, A Comparison of CP, IP, and SAT Solvers through a common interface, 2012, http://www.sci.brooklyn.cuny.edu/~zhou/papers/tai12.pdf

Wednesday, September 27, 2017

von Koch and Z3


Replica of the Zuse Z3 [5].

In [1] we showed a MIP formulation to label a tree graph in a “graceful” way [2]. Here we use the SMT (Satisfiability Modulo Theories [3]) solver Z3 from Microsoft [4]. Z3 is a library and for this example we use the Python bindings. (Z3 is also the name of an early computer, see the picture above).

The problem can be stated as follows. Using a tree graph with \(n\) nodes and \(n-1\) arcs, assign unique integer labels \(1,\dots,n\) to the nodes and \(1,\dots,n-1\) to the arcs such that for an arc \((i,j)\) we have

\[\mathit{ArcLabel}(i,j) = |\mathit{NodeLabel}(i)-\mathit{NodeLabel}(j)|\]

We use again the small dataset:

enter image description here

On the right is one possible solution.

So our model needs to handle the following conditions:

  1. \(1 \le \mathit{NodeLabel(i)} \le n\) and all-different
  2. \(1 \le \mathit{ArcLabel(i,j)} \le n-1\) and all-different
  3. \(\mathit{ArcLabel}(i,j) = |\mathit{NodeLabel}(i)-\mathit{NodeLabel}(j)|\)

The Z3 model is fairly high-level, and we can model these conditions quite naturally.

from z3 import *

nodes = ['a','b','c','d','e','f','g']

network = [('a','b'),('a','d'),('a','g'),
          ('b','c'),('b','e'),
          ('e','f')]

numnodes = len(nodes)
numarcs = len(network)

# nodes
X = [ Int('%s' % nodes[i]) for i in range(numnodes) ]

X_AllDifferent = Distinct( [ X[i] for i in range(numnodes) ] )
X_LB = And( [ X[i] >= 1 for i in range(numnodes) ] ) 
X_UB = And( [ X[i] <= numnodes for i in range(numnodes) ] ) 
X_cons = [ X_AllDifferent,X_LB,X_UB ]

# arcs
A = [ Int('%s.%s' % network[i]) for i in range(numarcs) ]

A_AllDifferent = Distinct( [ A[i] for i in range(numarcs) ] )
A_LB = And( [ A[i] >= 1 for i in range(numarcs) ] ) 
A_UB = And( [ A[i] <= numarcs for i in range(numarcs) ] ) 
A_cons = [ A_AllDifferent,A_LB,A_UB ]

# get node numbers of arcs
node1 = [nodes.index(network[i][0]) for i in range(numarcs)]
node2 = [nodes.index(network[i][1]) for i in range(numarcs)]

# no built-in absolute value
def Abs(x):
    return If(x >= 0,x,-x)

# difference
XA_cons = And([ A[i] == Abs(X[node1[i]]-X[node2[i]]) for i in range(numarcs) ])

all_cons = X_cons + A_cons + [ XA_cons ] 

solve(all_cons)

The main issue is that Z3 does not have a built-in Abs function. So, we implement our Abs using an If.

The output looks like:

[f = 5,
 b = 2,
 a = 7,
 c = 6,
 d = 1,
 e = 3,
 g = 4,
 e.f = 2,
 b.e = 1,
 b.c = 4,
 a.g = 3,
 a.d = 6,
 a.b = 5]
Enumerating all solutions

Like a MIP solver Z3 returns just one solution. We can enumerate all solutions by adding a constraint that forbids the current solution. For the above solution such a constraint can look like Or[f!=5,b!=2,a!=7,c!=6,d!=1,e!=3,g!=4]. A complete algorithm is:

s = Solver()
s.add(all_cons)
k = 0
res = []
while s.check() == sat:
    m = s.model()
    k = k + 1
    sol = [m.evaluate(X[i]) for i in range(numnodes)]
    res = res + [(k,sol)]
    forbid = Or([X[i] != sol[i] for i in range(numnodes)])
    s.add(forbid)     
print(res)

Indeed we also get our 52 solutions (only the node labels are reproduced here):

[(1, [4, 1, 6, 5, 7, 3, 2]), (2, [1, 5, 2, 6, 3, 4, 7]), (3, [1, 5, 2, 7, 3, 4, 6]), 
(4, [2, 7, 5, 3, 1, 4, 6]), (5, [1, 7, 4, 2, 3, 5, 6]), (6, [2, 6, 4, 5, 1, 7, 3]),
(7, [7, 3, 6, 1, 5, 4, 2]), (8, [7, 3, 6, 2, 5, 4, 1]), (9, [7, 1, 4, 6, 5, 3, 2]),
(10, [6, 1, 3, 5, 7, 4, 2]), (11, [7, 1, 3, 4, 5, 6, 2]), (12, [7, 1, 3, 2, 5, 6, 4]),
(13, [7, 1, 4, 2, 5, 3, 6]), (14, [7, 2, 4, 1, 5, 6, 3]), (15, [7, 2, 4, 3, 5, 6, 1]),
(16, [6, 2, 4, 5, 7, 1, 3]), (17, [6, 2, 4, 3, 7, 1, 5]), (18, [5, 1, 4, 6, 7, 2, 3]),
(19, [4, 1, 5, 6, 7, 2, 3]), (20, [3, 1, 5, 6, 7, 2, 4]), (21, [3, 1, 5, 4, 7, 2, 6]),
(22, [4, 1, 5, 3, 7, 2, 6]), (23, [2, 1, 6, 5, 7, 3, 4]), (24, [2, 1, 6, 4, 7, 3, 5]),
(25, [6, 1, 3, 2, 7, 4, 5]), (26, [4, 1, 6, 2, 7, 3, 5]), (27, [5, 1, 4, 3, 7, 2, 6]),
(28, [5, 7, 3, 2, 1, 6, 4]), (29, [5, 2, 6, 3, 7, 1, 4]), (30, [2, 6, 4, 3, 1, 7, 5]),
(31, [3, 7, 4, 2, 1, 6, 5]), (32, [1, 6, 4, 7, 3, 2, 5]), (33, [6, 7, 2, 3, 1, 5, 4]),
(34, [1, 7, 5, 6, 3, 2, 4]), (35, [3, 6, 2, 5, 1, 7, 4]), (36, [1, 6, 2, 7, 5, 3, 4]),
(37, [1, 6, 4, 5, 3, 2, 7]), (38, [1, 7, 5, 4, 3, 2, 6]), (39, [1, 6, 2, 4, 5, 3, 7]),
(40, [5, 2, 6, 4, 7, 1, 3]), (41, [3, 7, 4, 5, 1, 6, 2]), (42, [7, 2, 6, 1, 3, 5, 4]),
(43, [7, 2, 6, 4, 3, 5, 1]), (44, [3, 6, 2, 4, 1, 7, 5]), (45, [4, 7, 3, 2, 1, 6, 5]),
(46, [2, 7, 5, 6, 1, 4, 3]), (47, [1, 7, 4, 6, 3, 5, 2]), (48, [4, 7, 2, 3, 1, 5, 6]),
(49, [4, 7, 3, 5, 1, 6, 2]), (50, [5, 7, 3, 4, 1, 6, 2]), (51, [4, 7, 2, 6, 1, 5, 3]),
(52, [6, 7, 2, 4, 1, 5, 3])]

The last model was infeasible and caused the loop to end. We can in fact print this model and see how all the cuts where added:

[Distinct(a, b, c, d, e, f, g),
 And(a >= 1, b >= 1, c >= 1, d >= 1, e >= 1, f >= 1, g >= 1),
 And(a <= 7, b <= 7, c <= 7, d <= 7, e <= 7, f <= 7, g <= 7),
 Distinct(a.b, a.d, a.g, b.c, b.e, e.f),
 And(a.b >= 1,
     a.d >= 1,
     a.g >= 1,
     b.c >= 1,
     b.e >= 1,
     e.f >= 1),
 And(a.b <= 6,
     a.d <= 6,
     a.g <= 6,
     b.c <= 6,
     b.e <= 6,
     e.f <= 6),
 And(a.b == If(a - b >= 0, a - b, -(a - b)),
     a.d == If(a - d >= 0, a - d, -(a - d)),
     a.g == If(a - g >= 0, a - g, -(a - g)),
     b.c == If(b - c >= 0, b - c, -(b - c)),
     b.e == If(b - e >= 0, b - e, -(b - e)),
     e.f == If(e - f >= 0, e - f, -(e - f))),
 Or(a != 4, b != 1, c != 6, d != 5, e != 7, f != 3, g != 2),
 Or(a != 1, b != 5, c != 2, d != 6, e != 3, f != 4, g != 7),
 Or(a != 1, b != 5, c != 2, d != 7, e != 3, f != 4, g != 6),
 Or(a != 2, b != 7, c != 5, d != 3, e != 1, f != 4, g != 6),
 Or(a != 1, b != 7, c != 4, d != 2, e != 3, f != 5, g != 6),
 Or(a != 2, b != 6, c != 4, d != 5, e != 1, f != 7, g != 3),
 Or(a != 7, b != 3, c != 6, d != 1, e != 5, f != 4, g != 2),
 Or(a != 7, b != 3, c != 6, d != 2, e != 5, f != 4, g != 1),
 Or(a != 7, b != 1, c != 4, d != 6, e != 5, f != 3, g != 2),
 Or(a != 6, b != 1, c != 3, d != 5, e != 7, f != 4, g != 2),
 Or(a != 7, b != 1, c != 3, d != 4, e != 5, f != 6, g != 2),
 Or(a != 7, b != 1, c != 3, d != 2, e != 5, f != 6, g != 4),
 Or(a != 7, b != 1, c != 4, d != 2, e != 5, f != 3, g != 6),
 Or(a != 7, b != 2, c != 4, d != 1, e != 5, f != 6, g != 3),
 Or(a != 7, b != 2, c != 4, d != 3, e != 5, f != 6, g != 1),
 Or(a != 6, b != 2, c != 4, d != 5, e != 7, f != 1, g != 3),
 Or(a != 6, b != 2, c != 4, d != 3, e != 7, f != 1, g != 5),
 Or(a != 5, b != 1, c != 4, d != 6, e != 7, f != 2, g != 3),
 Or(a != 4, b != 1, c != 5, d != 6, e != 7, f != 2, g != 3),
 Or(a != 3, b != 1, c != 5, d != 6, e != 7, f != 2, g != 4),
 Or(a != 3, b != 1, c != 5, d != 4, e != 7, f != 2, g != 6),
 Or(a != 4, b != 1, c != 5, d != 3, e != 7, f != 2, g != 6),
 Or(a != 2, b != 1, c != 6, d != 5, e != 7, f != 3, g != 4),
 Or(a != 2, b != 1, c != 6, d != 4, e != 7, f != 3, g != 5),
 Or(a != 6, b != 1, c != 3, d != 2, e != 7, f != 4, g != 5),
 Or(a != 4, b != 1, c != 6, d != 2, e != 7, f != 3, g != 5),
 Or(a != 5, b != 1, c != 4, d != 3, e != 7, f != 2, g != 6),
 Or(a != 5, b != 7, c != 3, d != 2, e != 1, f != 6, g != 4),
 Or(a != 5, b != 2, c != 6, d != 3, e != 7, f != 1, g != 4),
 Or(a != 2, b != 6, c != 4, d != 3, e != 1, f != 7, g != 5),
 Or(a != 3, b != 7, c != 4, d != 2, e != 1, f != 6, g != 5),
 Or(a != 1, b != 6, c != 4, d != 7, e != 3, f != 2, g != 5),
 Or(a != 6, b != 7, c != 2, d != 3, e != 1, f != 5, g != 4),
 Or(a != 1, b != 7, c != 5, d != 6, e != 3, f != 2, g != 4),
 Or(a != 3, b != 6, c != 2, d != 5, e != 1, f != 7, g != 4),
 Or(a != 1, b != 6, c != 2, d != 7, e != 5, f != 3, g != 4),
 Or(a != 1, b != 6, c != 4, d != 5, e != 3, f != 2, g != 7),
 Or(a != 1, b != 7, c != 5, d != 4, e != 3, f != 2, g != 6),
 Or(a != 1, b != 6, c != 2, d != 4, e != 5, f != 3, g != 7),
 Or(a != 5, b != 2, c != 6, d != 4, e != 7, f != 1, g != 3),
 Or(a != 3, b != 7, c != 4, d != 5, e != 1, f != 6, g != 2),
 Or(a != 7, b != 2, c != 6, d != 1, e != 3, f != 5, g != 4),
 Or(a != 7, b != 2, c != 6, d != 4, e != 3, f != 5, g != 1),
 Or(a != 3, b != 6, c != 2, d != 4, e != 1, f != 7, g != 5),
 Or(a != 4, b != 7, c != 3, d != 2, e != 1, f != 6, g != 5),
 Or(a != 2, b != 7, c != 5, d != 6, e != 1, f != 4, g != 3),
 Or(a != 1, b != 7, c != 4, d != 6, e != 3, f != 5, g != 2),
 Or(a != 4, b != 7, c != 2, d != 3, e != 1, f != 5, g != 6),
 Or(a != 4, b != 7, c != 3, d != 5, e != 1, f != 6, g != 2),
 Or(a != 5, b != 7, c != 3, d != 4, e != 1, f != 6, g != 2),
 Or(a != 4, b != 7, c != 2, d != 6, e != 1, f != 5, g != 3),
 Or(a != 6, b != 7, c != 2, d != 4, e != 1, f != 5, g != 3)]
References
  1. The von Koch conjecture, http://yetanothermathprogrammingconsultant.blogspot.com/2017/09/the-von-koch-conjecture.html
  2. Graceful Labeling, https://en.wikipedia.org/wiki/Graceful_labeling
  3. Satisfiability modulo theories, https://en.wikipedia.org/wiki/Satisfiability_modulo_theories
  4. https://github.com/Z3Prover/z3
  5. Z3 (computer), https://en.wikipedia.org/wiki/Z3_(computer)

Thursday, September 21, 2017

The von Koch conjecture

Helge von Koch.jpg
Helge von Koch (1870-1924)

Consider a tree graph \(G=(V,E)\) with \(n\) nodes and \(n-1\) arcs. We can label the nodes with unique numbers \(1,\dots,n\) and the arcs with unique numbers \(1,\dots,n-1\) such that each arc has a label equal to the (absolute value) of difference of the labels of the attached nodes. This is also called “Graceful Labeling” [1].

Here is a picture from [2]:

enter image description here

We see that:

  • Each node has a unique label between 1 and 7 (the black numbers),
  • Each arc has a unique label between 1 and 6 (the red numbers),
  • Each arc has a label value equal to the difference in the labels of the associated nodes.

In [3] an interesting MIP formulation is proposed to find such a graceful labeling. There is no objective, i.e. we just look for a feasible solution. The uniqueness of node and arc labels is established by using an assignment block using binary variables (i.e. constraints we recognize from an assignment problem). A formulation can look like:

\[ \bbox[lightcyan,10px,border:3px solid darkblue]{ \begin{align}
\text{Node labels}\\&\sum_i x_{i,v} = 1\>\>\forall v&&\text{assignments of node labels}\\
&\sum_v x_{i,v} = 1 \>\>\forall i\\
&v^x_i  = \sum_v v \cdot x_{i,v}&&\text{value of node labels}\\
\hline
\text{Arc labels}\\&\sum_{(i,j)\in E} y_{i,j,e} = 1\>\forall e&&\text{assignment of arc labels}\\
&\sum_e  y_{i,j,e} = 1 \>\>\forall (i,j)\in E\\
&v^y_{i,j}  = \sum_e e \cdot y_{i,j,e}&&\text{value of arc labels}\\
\hline
\text{Relate arc/node}\\&v^y_{i,j} = | v^x_i - v^x_j | \>\>\forall (i,j)\in E&&\text{absolute value}&\\
\hline
\text{Variables}\\&x_{i,v} \in \{0,1\}&&\text{assign label $v$ to node $i$} \\
&y_{i,j,e} \in \{0,1\}&&\text{assign label $e$ to arc $(i,j)$}\\
&v^x_i \in \{1,\dots,n\}&&\text{value of node labels}\\
&v^y_{i,j} \in \{1,\dots,n-1\}&&\text{value of arc labels}
\end{align}} \]

It makes sense to form the set \(E\) only with arcs \(i<j\) to prevent duplication (if we have an arc \(a \rightarrow b\) we don’t want to check also \(b \rightarrow a\)).

The absolute value equation still needs to be linearized, which we detail next.

Linearizing the absolute value

A standard formulation for \(y=|x|\) uses “\(y= -x\) or \(y=x\)”, i.e.

\[\begin{align} & –x – \delta M \le y \le –x + \delta M\\
& x –(1-\delta) M \le y \le x + (1-\delta) M\\&\delta \in \{0,1\} \end{align} \]

We can slightly improve, by using:

\[\begin{align} & –x \le y \le –x + \delta M\\
& x \le y \le x + (1-\delta) M\\&\delta \in \{0,1\} \end{align} \]

For our model we can get good values for the big-M's:

\[\begin{align}
&v^y_{i,j} \ge v^x_i - v^x_j\\
&v^y_{i,j} \ge v^x_j - v^x_i\\
&v^y_{i,j} \le v^x_j - v^x_i  + 2 n \delta_{i,j} \\
&v^y_{i,j} \le v^x_j - v^x_j  + 2 n (1-\delta_{i,j}) \\
& \delta_{i,j} \in \{0,1\}
\end{align}\]

This is the formulation presented in [3]. The solution is not unique. I see:

----     74 VARIABLE vx.L  node labels

a 7,    b 1,    c 4,    d 2,    e 5,    f 3,    g 6


----     74 VARIABLE vy.L  arc labels

            b           c           d           e           f           g

a           6                       5                                   1
b                       3                       4
e                                                           2

This is different from the solution in the picture.

Alternatives?

It seems tempting to try using an alternative approach to the absolute value constraint. May be something like:

\[\begin{align}\min\>&\sum_{(i,j)\in E} v^y_{i,j}\\
&v^y_{i,j} \ge v^x_i - v^x_j\\
&v^y_{i,j} \ge v^x_j - v^x_i\end{align}\]

These formulations do not work (in this case the objective is actually constant and does not really drive down the individual \(v^y\)’s).

Enumerate solutions

We already established that the solution is not unique. How many solutions are there for this small example problem?

We can add a simple cut:

\[\sum_{i,v} \alpha_{i,v} x_{i,v} \le n-1\]

where \(\alpha_{i,v}=x^*_{i,v}\)  is a previously found solution. This will forbid this solution to be found again. If we solve again we will either find a new, different solution, or we see the model has become infeasible. In the latter case we can conclude there are no more solutions to be found. 

If we repeat this, we find a surprisingly large number of solutions (I expected to see a dozen or so). Below are the node and arc labels we collected in the solve loop:

----    124 PARAMETER vxall  collected node labels

              a           b           c           d           e           f           g

k1            7           1           4           2           5           3           6
k2            5           7           3           4           1           6           2
k3            1           6           4           5           3           2           7
k4            3           6           2           4           1           7           5
k5            7           3           6           2           5           4           1
k6            1           7           5           4           3           2           6
k7            1           5           2           7           3           4           6
k8            6           2           4           5           7           1           3
k9            7           2           4           3           5           6           1
k10           2           7           5           6           1           4           3
k11           4           7           3           2           1           6           5
k12           3           7           4           2           1           6           5
k13           2           1           6           5           7           3           4
k14           4           7           3           5           1           6           2
k15           7           3           6           1           5           4           2
k16           6           1           3           2           7           4           5
k17           7           2           6           4           3           5           1
k18           5           1           4           3           7           2           6
k19           5           1           4           6           7           2           3
k20           6           7           2           3           1           5           4
k21           1           6           2           4           5           3           7
k22           7           2           6           1           3           5           4
k23           1           6           4           7           3           2           5
k24           1           6           2           7           5           3           4
k25           4           7           2           3           1           5           6
k26           7           1           4           6           5           3           2
k27           7           2           4           1           5           6           3
k28           5           7           3           2           1           6           4
k29           7           1           3           4           5           6           2
k30           1           7           5           6           3           2           4
k31           4           1           5           3           7           2           6
k32           3           7           4           5           1           6           2
k33           1           5           2           6           3           4           7
k34           3           6           2           5           1           7           4
k35           4           1           6           2           7           3           5
k36           2           7           5           3           1           4           6
k37           4           1           5           6           7           2           3
k38           1           7           4           6           3           5           2
k39           6           1           3           5           7           4           2
k40           6           7           2           4           1           5           3
k41           4           7           2           6           1           5           3
k42           6           2           4           3           7           1           5
k43           2           1           6           4           7           3           5
k44           4           1           6           5           7           3           2
k45           7           1           3           2           5           6           4
k46           5           2           6           3           7           1           4
k47           3           1           5           6           7           2           4
k48           5           2           6           4           7           1           3
k49           2           6           4           3           1           7           5
k50           3           1           5           4           7           2           6
k51           1           7           4           2           3           5           6
k52           2           6           4           5           1           7           3


----    124 PARAMETER vyall  collected arc labels

            a.b         a.d         a.g         b.c         b.e         e.f

k1            6           5           1           3           4           2
k2            2           1           3           4           6           5
k3            5           4           6           2           3           1
k4            3           1           2           4           5           6
k5            4           5           6           3           2           1
k6            6           3           5           2           4           1
k7            4           6           5           3           2           1
k8            4           1           3           2           5           6
k9            5           4           6           2           3           1
k10           5           4           1           2           6           3
k11           3           2           1           4           6           5
k12           4           1           2           3           6           5
k13           1           3           2           5           6           4
k14           3           1           2           4           6           5
k15           4           6           5           3           2           1
k16           5           4           1           2           6           3
k17           5           3           6           4           1           2
k18           4           2           1           3           6           5
k19           4           1           2           3           6           5
k20           1           3           2           5           6           4
k21           5           3           6           4           1           2
k22           5           6           3           4           1           2
k23           5           6           4           2           3           1
k24           5           6           3           4           1           2
k25           3           1           2           5           6           4
k26           6           1           5           3           4           2
k27           5           6           4           2           3           1
k28           2           3           1           4           6           5
k29           6           3           5           2           4           1
k30           6           5           3           2           4           1
k31           3           1           2           4           6           5
k32           4           2           1           3           6           5
k33           4           5           6           3           2           1
k34           3           2           1           4           5           6
k35           3           2           1           5           6           4
k36           5           1           4           2           6           3
k37           3           2           1           4           6           5
k38           6           5           1           3           4           2
k39           5           1           4           2           6           3
k40           1           2           3           5           6           4
k41           3           2           1           5           6           4
k42           4           3           1           2           5           6
k43           1           2           3           5           6           4
k44           3           1           2           5           6           4
k45           6           5           3           2           4           1
k46           3           2           1           4           5           6
k47           2           3           1           4           6           5
k48           3           1           2           4           5           6
k49           4           1           3           2           5           6
k50           2           1           3           4           6           5
k51           6           1           5           3           4           2
k52           4           3           1           2           5           6

Note that there are \(7!=5040\) ways to order seven different labels, so from that perspective this is indeed a small subset of feasible node labels.

Instead of doing this solve loop ourselves we can also use the Solution Pool technology available in some solvers such as Cplex and Gurobi. This automates this loop (and does it way smarter and more efficient).

References
  1. Graceful Labeling, https://en.wikipedia.org/wiki/Graceful_labeling
  2. The von Koch conjecture, https://codegolf.stackexchange.com/questions/119263/the-von-koch-conjecture
  3. Mike Appleby, https://github.com/appleby/graceful-tree
  4. Helge von Koch, https://en.wikipedia.org/wiki/Helge_von_Koch
  5. Barbara M. Smith, Jean-Francois Puget, Constraint Models for graceful graphs, Constraints (2010), vol. 15, pp.64-92, https://link.springer.com/article/10.1007/s10601-009-9071-6. This paper discusses some interesting alternative models when attacking this as a Constraint Programming problem.