Discrete optimization in Mathematica 4.2

Daniel Lichtblau

Abstract

Acknowledgements and remarks

Much of the NMinimize code on which this notebook rests was developed by Sergei Shebalov (UIUC Math department, WRI intern, summers 1999, 2000) and Brett Champion (WRI); they also did some of the development on the specific examples we present.

There is a discussion of the technology behind NMinimize in Brett Champion's talk for this conference. The primary focus of that talk is to describe the general technology and provide examples of continuous optimization problems, the notebook for that talk also contains some examples of discrete optimization.

All timings indicated below are for a 1.4 GHz machine running a development Mathematica 4.2 kernel under the Linux operating system.

Initialization

Needs["NumericalMath`NMinimize`"];
Off[General::spell,  General::spell1];

SetOptions[NMinimize, AccuracyGoal -> Infinity, 
        PrecisionGoal -> Infinity, Compiled -> False];

A simple example

We start this talk by searching for Pythagorian triples (integer triples {x, y, z} such that (x^2) + (y^2) == (z^2)). Our objective function to minimize will be just the square of the difference of lefthand side and right hand side in the desired equality above. Note that this is to some extent arbitrary, and minimizing the absolute value rather than the square would suffice just as well for our purpose. We will specify explicitly that the variables are integer valued.

NMinimize[{(x^2 + y^2 - z^2)^2,  {x, y, z} ∈ Integers, 
        5 ≤ x ≤ 25, 5 ≤ y ≤ 25, 
        5 ≤ z ≤ 25, x ≤ y} ,  {x, y, 
        z}]

We see that NMinimize is able to pick an appropriate method by default. Indeed, it uses DifferentialEvolution (an evolutionary method) when given discrete variables. We will say more about use of this method later. Now we show how to obtain different solutions by specifying that the random seed used by the DifferentialEvolution method change for each run.

Timing[
      Table[NMinimize[{(x^2 + y^2 - z^2)^2,  {x, y, z} ∈ 
              Integers, 5 ≤ x ≤ 25, 
            5 ≤ y ≤ 25, 
            5 ≤ z ≤ 25, x ≤ y} ,  {x, y, 
            z}, Method -> {DifferentialEvolution, 
              RandomSeed -> Random[Integer, 1000]}],  {20}]]

Another simple example

We now tackle a coin problem. We are given 143267 coins in pennies, nickels, dimes, and quarters, of total value $12563.29, and we are to determine how many coins might be of each type. There are several ways one might set up such a problem in NMinimize. We will try to minimize the sum of squares of differences between actual values and desired values of the two linear expressions implied by the information above. For our search space we will impose obvious range constraints on the various coin types. Again we will want to alter the seeding of the random number generator so we specify the method with this option added.

Timing[
      Table[{min, sol} = 
          NMinimize[{(((pennies + 5 nickels + 10 dimes + 25 quarters - 1256329))^2) + (((pennies + nickels + dimes + quarters - 143267))^2), {pennies, nickels, dimes, quarters} ∈ Integers, 
              0 ≤ pennies ≤ 1256329, 
              0 ≤ nickels ≤ (1256329/5), 
              0 ≤ dimes ≤ (1256329/10), 
              0 ≤ quarters ≤ (1256329/25)}, {pennies, 
              nickels, dimes, quarters}, MaxIterations -> 400, 
            Method -> {DifferentialEvolution, 
                RandomSeed -> Random[Integer, 1000]}], {10}]]

We obtained several solutions and "near" solutions. As with the previous example, this type of problem is one of constraint satisfaction. An advantage in working with these is that we know from the proposed solution whether it is valid; those are exactly the cases for which we get a value of zero, with all constraints satisfied.

A set partitioning problem

This is an old example from computational folklore. We are to partition the integers from 1 to 100 into two sets of 50, such that the sums of the square roots in each set are as close to equal as possible. There are various ways to set this up as a problem for NMinimize, and we illustrate two such below.

method 1

The first illustrates a simple way to choose 50 elements from the set of 100, and is an approach we often use for adapting permutation problems for NMinimize. We work with 100 real values from 0 to 1 and their ordering (from the Mathematica Ordering function) determines which is to be regarded as "first", which as "second", and so on. Note that this is different from previous examples in that we now work with continuous variables even though the problem involves a discrete set.

splitRange[vec_] := 
      With[{newvec = Ordering[vec], 
          halflen = Floor[Length[vec]/2]}, {Take[newvec, halflen], 
          Drop[newvec, halflen]}]

Once we have a way to associate a pair of subsets to a given set of 100 values in the range from 0 to 1, we form our objective function. As in many above examples this will simply be a square of a difference.

obfun[vec : {__Real}] := 
      With[{vals = splitRange[vec]}, 
        Abs[(Apply[Plus, Sqrt[N[First[vals]]]] - 
              Apply[Plus, Sqrt[N[Last[vals]]]])]]

We now put these components together into a function that provides our set partition.

getHalfSet[n_, opts___Rule] := 
      Module[{vars, xx, ranges, nmin, vals}, 
vars = 
          Array[xx, n]; 
ranges = 
          Map[{#, 0, 1} &, vars]; 
{nmin, vals} = 
          NMinimize[obfun[vars], ranges,  opts]; 
{nmin, 
          Map[Sort, splitRange[vars /. vals]]}]

Timing[{min, subsets} = getHalfSet[100]]

While this is not terrible one might hope for a far better result. It turns out that the default method used in this case is NelderMead, because there are no discrete variables in the problem specification. While better automatic determination of the most appropriate method remains an open issue, we can certainly force use of the DifferentialEvolution method. We will also disable the automated PostProcess attempt as that also does not always determine situations in which it would be of little use, and it can be costly in time. Finally, we set this to run many iterations with alot of search points.

SeedRandom[11111]
Timing[{min, {s1, s2}} = 
        getHalfSet[100, 
          Method -> {DifferentialEvolution, 
              CrossProbability ->  .9, SearchPoints -> 100, 
              PostProcess -> False},  
          MaxIterations -> 10000]]

method 2

Another approach to this problem is as follows. We take the full set, Range[100], and pick 100 corresponding random 0-1 values. An element in the range is put into one or the other subset according to the value of the bit corresponding to that element. We also require a constraint that the size of each subset is half the total size.

obfun2[vec : {__Integer}] := 
      Exp[10*Abs[(2*vec - 1) . Sqrt[N[Range[Length[vec]]]]]]

\(getHalfSet2[n_, opts___] := 
      Module[
{vars, x, nmin, vals, ranges, s1}, 
        vars = Array[x, n]; 
ranges = 
          Map[\((0 ≤ # ≤ 1)\) &, vars]; {nmin, vals} = 
          NMinimize[{obfun2[vars], 
              Join[ranges, {Element[vars, Integers], 
                  Apply[Plus, vars] == n/2}]}, vars, 
            opts]; 
s1\  = 
          Select[Inner[Times, Range[n], \((vars /. vals)\), 
              List], # ≠ 0 &]; 
{nmin, {s1, 
            Complement[Range[n], s1]}}
]\)

Timing[{min, {s1, s2}} = 
        getHalfSet2[100, MaxIterations -> 3000]]

We check that we have a valid partition (both sets of equal size) and also see how close are the sums of square roots of elements in each set.

Length[s1] == Length[s2]
Abs[Apply[Plus, Sqrt[N[s1]]] - Apply[Plus, Sqrt[N[s2]]]]

While promising in principle we see that this method did not fare so well in practice.

A subset covering problem

The problem below was posed in the news group comp.soft-sys.math.mathematica. We are given a set of sets, each containing integers between 1 and 64, and whose union is exactly Range[64]. We want to find a set of 12 subsets that cover the entire range of integers 1 through 64.

subsets = {{1, 2, 4, 8, 16, 32, 64}, {2, 1, 3, 7, 15, 31, 63}, {3, 4, 
            2, 6, 14, 30, 62}, {4, 3, 1, 5, 13, 29, 61}, {5, 6, 8, 4, 12, 28, 
            60}, {6, 5, 7, 3, 11, 27, 59}, {7, 8, 6, 2, 10, 26, 58}, {8, 7, 
            5, 1, 9, 25, 57}, {9, 10, 12, 16, 8, 24, 56}, {10, 9, 11, 15, 7, 
            23, 55}, {11, 12, 10, 14, 6, 22, 54}, {12, 11, 9, 13, 5, 21, 
            53}, {13, 14, 16, 12, 4, 20, 52}, {14, 13, 15, 11, 3, 19, 
            51}, {15, 16, 14, 10, 2, 18, 50}, {16, 15, 13, 9, 1, 17, 
            49}, {17, 18, 20, 24, 32, 16, 48}, {18, 17, 19, 23, 31, 15, 
            47}, {19, 20, 18, 22, 30, 14, 46}, {20, 19, 17, 21, 29, 13, 
            45}, {21, 22, 24, 20, 28, 12, 44}, {22, 21, 23, 19, 27, 11, 
            43}, {23, 24, 22, 18, 26, 10, 42}, {24, 23, 21, 17, 25, 9, 
            41}, {25, 26, 28, 32, 24, 8, 40}, {26, 25, 27, 31, 23, 7, 
            39}, {27, 28, 26, 30, 22, 6, 38}, {28, 27, 25, 29, 21, 5, 
            37}, {29, 30, 32, 28, 20, 4, 36}, {30, 29, 31, 27, 19, 3, 
            35}, {31, 32, 30, 26, 18, 2, 34}, {32, 31, 29, 25, 17, 1, 
            33}, {33, 34, 36, 40, 48, 64, 32}, {34, 33, 35, 39, 47, 63, 
            31}, {35, 36, 34, 38, 46, 62, 30}, {36, 35, 33, 37, 45, 61, 
            29}, {37, 38, 40, 36, 44, 60, 28}, {38, 37, 39, 35, 43, 59, 
            27}, {39, 40, 38, 34, 42, 58, 26}, {40, 39, 37, 33, 41, 57, 
            25}, {41, 42, 44, 48, 40, 56, 24}, {42, 41, 43, 47, 39, 55, 
            23}, {43, 44, 42, 46, 38, 54, 22}, {44, 43, 41, 45, 37, 53, 
            21}, {45, 46, 48, 44, 36, 52, 20}, {46, 45, 47, 43, 35, 51, 
            19}, {47, 48, 46, 42, 34, 50, 18}, {48, 47, 45, 41, 33, 49, 
            17}, {49, 50, 52, 56, 64, 48, 16}, {50, 49, 51, 55, 63, 47, 
            15}, {51, 52, 50, 54, 62, 46, 14}, {52, 51, 49, 53, 61, 45, 
            13}, {53, 54, 56, 52, 60, 44, 12}, {54, 53, 55, 51, 59, 43, 
            11}, {55, 56, 54, 50, 58, 42, 10}, {56, 55, 53, 49, 57, 41, 
            9}, {57, 58, 60, 64, 56, 40, 8}, {58, 57, 59, 63, 55, 39, 
            7}, {59, 60, 58, 62, 54, 38, 6}, {60, 59, 57, 61, 53, 37, 
            5}, {61, 62, 64, 60, 52, 36, 4}, {62, 61, 63, 59, 51, 35, 
            3}, {63, 64, 62, 58, 50, 34, 2}, {64, 63, 61, 57, 49, 33, 
            1}};
Union[Flatten[subsets]] == Range[64]

method 1

We set up our objective function as follows. We represent a set of 12 subsets of this "universe" set by a set of 12 integers in the range from 1 to the number of subsets (which in this example is also 64). This set is allowed to contain repetitions. Our objective function to minimize will be based on how many elements from 1 through 64 are "covered". Specifically it will be 2 raised to the #(elements not covered) power. The code below does this. As it is a problem that may require many iterations and a relatively large search space, we allow for nondefault option values for these parameters.

f[n : {___Integer}, set_, mx_Integer] := 
      2^Length[Complement[Range[mx], 
            Union[Flatten[set[[n]]]]]]

spanningSets[set_, nsets_, iter_, sp_, cp_] := 
      Module[{vars, rnges, max = Length[set], nmin, vals}, 
        vars = Array[xx, nsets]; 
rnges = 
          Map[(1 ≤ # ≤ max) &, 
            vars]; 
{nmin, vals} = 
          NMinimize[{f[vars, set, max], 
              Append[rnges, Element[vars, Integers]]}, vars, 
            MaxIterations -> iter, 
            Method -> {DifferentialEvolution, SearchPoints -> sp, 
                CrossProbability -> cp}]; 
vals = 
          Union[vars /. vals]; 
{nmin, vals}]

Timing[{min, sets} = 
        spanningSets[subsets, 12, 700, 200,  .94]]

Length[Union[Flatten[subsets[[sets]]]]]

This example has an unfortunate drawback. It is requires a nondefault value for an obscure option, CrossProbability, in the DifferentialEvolution method of NMinimize. It is not entirely trivial to find useful values for such options, or even to know which such options to alter from default values. We hope to ameliorate this need in future work by providing an adaptive evolution that modifies its own parameters as it proceeds.

method 2

This problem may also be tackled in other ways. One is to write a simple "greedy" algorithm. Experience indicates that such an approach will fail for this particular set. Another, as was pointed out in the news group, is to cast it as a knapsack problem. We show details below. First we transform our set of subsets into a  "bit vector" representation; each subset is represented by a positional list of zeros and ones.

densevec[spvec_, 
        len_]:=Module[{vec = Table[0, {len}]}, 
Do[
          vec[[spvec[[j]]]] = 1, {j, 
            Length[spvec]}]; 
vec
]

mat = Map[densevec[#, 64] &, subsets];

We now work with 0-1 variables and minimize their sum, subject to the constraint that an appropriate inner product be strictly positive.

spanningSets2[set_, iter_, sp_, seed_, cp_:  .5] := 
      Module[{vars, rnges, max = Length[set], nmin, vals}, 
        vars = Array[xx, max]; 
rnges = 
          Map[(0 ≤ # ≤ 1) &, 
            vars]; 
{nmin, vals} = 
          NMinimize[{Apply[Plus, vars], 
              Join[rnges, {Element[vars, Integers]}, 
                Thread[vars . set ≥ Table[1, {max}]]]}, vars, 
            MaxIterations -> iter, 
            Method -> {DifferentialEvolution, 
                CrossProbability -> cp, SearchPoints -> sp, 
                RandomSeed -> seed}]; 
vals = 
          vars /. vals; 
{nmin, vals}]

Timing[{min, sets} = spanningSets2[mat, 5000, 100, 0,  .1]]

Timing[{min, sets} = spanningSets2[mat, 2000, 100, 0,  .1]]

We have a result that uses 14 subsets.

Map[Min[#, 1] &, sets . mat]

As we already know this solution is not optimal, and it is slow to obtain. We conclude that atpresent it appears NMinimize will not do too well with this formulation. All the same it is encouraging to realize that one may readily set this up as a standard knapsack problem; over time NMinimize internals may improve to the point where it can more effectively handle the problem in this way.

An assignment problem

Our next example is a benchmark from the literature of discrete optimization. We are given two matrices. We want a permutation p such that when applied to the rows and columns of the second matrix, multiplied element-wise with corresponding elements of the first, and all elements summed, we get a minimum. The matrices we use are 25x25. This is known as the "NUG25" example and it is an example of a "quadratic assignment problem" (QAP). The optimal result is known and was verified by a large parallel computation. The actual matrices are in a hidden cell below. One should recognize that the methods of handling this problem can, with minor modification, be applied to related ones such as the travelling salesman problem. In general, problems that require a "best permutation" may be amenable to the methods described below.

aa = {{0, 1, 2, 3, 4, 1, 2, 3, 4, 5, 2, 3, 4, 5, 6, 3, 4, 5, 6, 7, 4, 
            5, 6, 7, 8}, {1, 0, 1, 2, 3, 2, 1, 2, 3, 4, 3, 2, 3, 4, 5, 4, 3, 
            4, 5, 6, 5, 4, 5, 6, 7}, {2, 1, 0, 1, 2, 3, 2, 1, 2, 3, 4, 3, 2, 
            3, 4, 5, 4, 3, 4, 5, 6, 5, 4, 5, 6}, {3, 2, 1, 0, 1, 4, 3, 2, 1, 
            2, 5, 4, 3, 2, 3, 6, 5, 4, 3, 4, 7, 6, 5, 4, 5}, {4, 3, 2, 1, 0, 
            5, 4, 3, 2, 1, 6, 5, 4, 3, 2, 7, 6, 5, 4, 3, 8, 7, 6, 5, 4}, {1, 
            2, 3, 4, 5, 0, 1, 2, 3, 4, 1, 2, 3, 4, 5, 2, 3, 4, 5, 6, 3, 4, 5, 
            6, 7}, {2, 1, 2, 3, 4, 1, 0, 1, 2, 3, 2, 1, 2, 3, 4, 3, 2, 3, 4, 
            5, 4, 3, 4, 5, 6}, {3, 2, 1, 2, 3, 2, 1, 0, 1, 2, 3, 2, 1, 2, 3, 
            4, 3, 2, 3, 4, 5, 4, 3, 4, 5}, {4, 3, 2, 1, 2, 3, 2, 1, 0, 1, 4, 
            3, 2, 1, 2, 5, 4, 3, 2, 3, 6, 5, 4, 3, 4}, {5, 4, 3, 2, 1, 4, 3, 
            2, 1, 0, 5, 4, 3, 2, 1, 6, 5, 4, 3, 2, 7, 6, 5, 4, 3}, {2, 3, 4, 
            5, 6, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 1, 2, 3, 4, 5, 2, 3, 4, 5, 
            6}, {3, 2, 3, 4, 5, 2, 1, 2, 3, 4, 1, 0, 1, 2, 3, 2, 1, 2, 3, 4, 
            3, 2, 3, 4, 5}, {4, 3, 2, 3, 4, 3, 2, 1, 2, 3, 2, 1, 0, 1, 2, 3, 
            2, 1, 2, 3, 4, 3, 2, 3, 4}, {5, 4, 3, 2, 3, 4, 3, 2, 1, 2, 3, 2, 
            1, 0, 1, 4, 3, 2, 1, 2, 5, 4, 3, 2, 3}, {6, 5, 4, 3, 2, 5, 4, 3, 
            2, 1, 4, 3, 2, 1, 0, 5, 4, 3, 2, 1, 6, 5, 4, 3, 2}, {3, 4, 5, 6, 
            7, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 1, 2, 3, 4, 
            5}, {4, 3, 4, 5, 6, 3, 2, 3, 4, 5, 2, 1, 2, 3, 4, 1, 0, 1, 2, 3, 
            2, 1, 2, 3, 4}, {5, 4, 3, 4, 5, 4, 3, 2, 3, 4, 3, 2, 1, 2, 3, 2, 
            1, 0, 1, 2, 3, 2, 1, 2, 3}, {6, 5, 4, 3, 4, 5, 4, 3, 2, 3, 4, 3, 
            2, 1, 2, 3, 2, 1, 0, 1, 4, 3, 2, 1, 2}, {7, 6, 5, 4, 3, 6, 5, 4, 
            3, 2, 5, 4, 3, 2, 1, 4, 3, 2, 1, 0, 5, 4, 3, 2, 1}, {4, 5, 6, 7, 
            8, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 0, 1, 2, 3, 
            4}, {5, 4, 5, 6, 7, 4, 3, 4, 5, 6, 3, 2, 3, 4, 5, 2, 1, 2, 3, 4, 
            1, 0, 1, 2, 3}, {6, 5, 4, 5, 6, 5, 4, 3, 4, 5, 4, 3, 2, 3, 4, 3, 
            2, 1, 2, 3, 2, 1, 0, 1, 2}, {7, 6, 5, 4, 5, 6, 5, 4, 3, 4, 5, 4, 
            3, 2, 3, 4, 3, 2, 1, 2, 3, 2, 1, 0, 1}, {8, 7, 6, 5, 4, 7, 6, 5, 
            4, 3, 6, 5, 4, 3, 2, 5, 4, 3, 2, 1, 4, 3, 2, 1, 0}};
bb = {{0, 3, 2, 0, 0, 10, 5, 0, 5, 2, 0, 0, 2, 0, 5, 3, 0, 1, 10, 0, 
            2, 1, 1, 1, 0}, {3, 0, 4, 0, 10, 0, 0, 2, 2, 1, 5, 0, 0, 0, 0, 0, 
            1, 6, 1, 0, 2, 2, 5, 1, 10}, {2, 4, 0, 3, 4, 5, 5, 5, 1, 4, 0, 4, 
            0, 4, 0, 3, 2, 5, 5, 2, 0, 0, 3, 1, 0}, {0, 0, 3, 0, 0, 0, 2, 2, 
            0, 6, 2, 5, 2, 5, 1, 1, 1, 2, 2, 4, 2, 0, 2, 2, 5}, {0, 10, 4, 0, 
            0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 5, 0, 2, 1, 0, 
            2}, {10, 0, 5, 0, 2, 0, 10, 10, 5, 10, 6, 0, 0, 10, 2, 10, 1, 5, 
            5, 2, 5, 0, 2, 0, 1}, {5, 0, 5, 2, 0, 10, 0, 1, 3, 5, 0, 0, 2, 4, 
            5, 10, 6, 0, 5, 5, 5, 0, 5, 5, 0}, {0, 2, 5, 2, 0, 10, 1, 0, 10, 
            2, 5, 2, 0, 3, 0, 0, 0, 4, 0, 5, 0, 5, 2, 2, 5}, {5, 2, 1, 0, 0, 
            5, 3, 10, 0, 5, 6, 0, 1, 5, 5, 5, 2, 3, 5, 0, 2, 10, 10, 1, 
            5}, {2, 1, 4, 6, 0, 10, 5, 2, 5, 0, 0, 1, 2, 1, 0, 0, 0, 0, 6, 6, 
            4, 5, 3, 2, 2}, {0, 5, 0, 2, 0, 6, 0, 5, 6, 0, 0, 2, 0, 4, 2, 1, 
            0, 6, 2, 1, 5, 0, 0, 1, 5}, {0, 0, 4, 5, 0, 0, 0, 2, 0, 1, 2, 0, 
            2, 1, 0, 3, 10, 0, 0, 4, 0, 0, 4, 2, 5}, {2, 0, 0, 2, 0, 0, 2, 0, 
            1, 2, 0, 2, 0, 4, 5, 0, 1, 0, 5, 0, 0, 0, 5, 1, 1}, {0, 0, 4, 5, 
            0, 10, 4, 3, 5, 1, 4, 1, 4, 0, 0, 0, 2, 2, 0, 2, 5, 0, 5, 2, 
            5}, {5, 0, 0, 1, 2, 2, 5, 0, 5, 0, 2, 0, 5, 0, 0, 2, 0, 0, 0, 6, 
            3, 5, 0, 0, 5}, {3, 0, 3, 1, 0, 10, 10, 0, 5, 0, 1, 3, 0, 0, 2, 
            0, 0, 5, 5, 1, 5, 2, 1, 2, 10}, {0, 1, 2, 1, 0, 1, 6, 0, 2, 0, 0, 
            10, 1, 2, 0, 0, 0, 5, 2, 1, 1, 5, 6, 5, 5}, {1, 6, 5, 2, 2, 5, 0, 
            4, 3, 0, 6, 0, 0, 2, 0, 5, 5, 0, 4, 0, 0, 0, 0, 5, 0}, {10, 1, 5, 
            2, 0, 5, 5, 0, 5, 6, 2, 0, 5, 0, 0, 5, 2, 4, 0, 5, 4, 4, 5, 0, 
            2}, {0, 0, 2, 4, 5, 2, 5, 5, 0, 6, 1, 4, 0, 2, 6, 1, 1, 0, 5, 0, 
            4, 4, 1, 0, 2}, {2, 2, 0, 2, 0, 5, 5, 0, 2, 4, 5, 0, 0, 5, 3, 5, 
            1, 0, 4, 4, 0, 1, 0, 10, 1}, {1, 2, 0, 0, 2, 0, 0, 5, 10, 5, 0, 
            0, 0, 0, 5, 2, 5, 0, 4, 4, 1, 0, 0, 0, 0}, {1, 5, 3, 2, 1, 2, 5, 
            2, 10, 3, 0, 4, 5, 5, 0, 1, 6, 0, 5, 1, 0, 0, 0, 0, 0}, {1, 1, 1, 
            2, 0, 0, 5, 2, 1, 2, 1, 2, 1, 2, 0, 2, 5, 5, 0, 0, 10, 0, 0, 0, 
            2}, {0, 10, 0, 5, 2, 1, 0, 5, 5, 2, 5, 5, 1, 5, 5, 10, 5, 0, 2, 
            2, 1, 0, 0, 2, 0}};

method 1

Our first problem is to decide how one can make a set of values into a permutation. One approach, as with the set partition problem, is to take their Ordering. In order that this avoid collisions we thus work with real values.

permuteMatrix[mat_, perm_] := mat[[perm, perm]]

QAP[mat1_, mat2_, cp_, it_, sp_, sc_] := 
      Module[{len = Length[mat1], obfunc, vars, vv, nmin, vals, rnges}, 
        vars = Array[vv, len]; 
rnges = 
          Map[{#, 0, 1} &, vars]; 
obfunc[vec : {__Real}] := 
          Apply[Plus, 
            Flatten[mat1*
                permuteMatrix[mat2, 
                  Ordering[vec]]]]; 
{nmin, vals} = 
          NMinimize[obfunc[vars], rnges, MaxIterations -> it, 
            PostProcess -> False, 
            Method -> {DifferentialEvolution, SearchPoints -> sp, 
                CrossProbability -> cp, 
                ScalingFactor -> sc}]; 
{nmin, 
          Ordering[vars /. vals]}]

Again we face the issue that this problem requires nonstandard values for options to the DifferentialEvolution method, in order to achieve a reasonable result. While this is regretable it is clearly better than having no recourse at all. The specific values we use were found by Brett Champion by running many short trials and noticing which parameter values seemed to work best (it was also his idea to try nondefault values in the first place).

The idea behind having CrossProbability relatively small is that we do not want many "genetic crossovers" in mating a pair of vectors. This in turn is because of the way we define a permutation. In particular it is not just values but relative values across the entire vector that give us the permutation, and thus disrupting more than a few, even when mating a pair of "good" vectors, is likely to give a "random" bad vector.

Below is the optimal permutation and the objective function value we obtain therefrom. We also show as baseline the result from applying no permutation, and also the result of applying several random permutations. This gives some idea of how to gauge the results below.

p = {5, 11, 20, 15, 22, 2, 25, 8, 9, 1, 18, 16, 3, 6, 19, 24, 21, 14, 
          7, 10, 17, 12, 4, 23, 13};
best = Apply[Plus, Flatten[aa*permuteMatrix[bb, p]]]
baseline = Apply[Plus, Flatten[aa*bb]]
randomvals = 
      Table[
perm = 
          Ordering[Table[Random[], {25}]]; 
Apply[Plus, 
          Flatten[aa*
              permuteMatrix[bb, perm]]], 
{10}]

A much more strenuous run over random permutations may be indicative of how hard it is to get good results by randomized search.

SeedRandom[1111];
Timing[randomvals = 
          Table[
perm = 
              Ordering[Table[Random[], {25}]]; 
Apply[Plus, 
              Flatten[aa*
                  permuteMatrix[bb, 
                    perm]]], 
{1000000}];]
Min[randomvals]

We see that the baseline permutation (do nothing) and random permutations tend to be far from optimal, and even a large sampling will get us only about half way from baseline to optimal. A relatively brief run with good values for the algorithm parameters, on the other hand, yields something better still.

SeedRandom[11111];
Timing[{min, perm} = QAP[aa, bb,  .06, 200, 40,  .6]]

We now try a longer run.

SeedRandom[11111];
Timing[{min, perm} = QAP[aa, bb,  .06, 5000, 100,  .6]]

method 2

Once again we have a very different approach to this problem. The idea is to generate a permutation as a "shuffle" of a set of integers. We have for vector a set of integers from 1 to len, the length of the set in question. The range restriction is the only stipulation and in particular it may contain repeats. We associate to it a unique permutation as follows. We initialize a "deck" to contain len zeros . The first element in deck is then set to vector[[1]]. We also have a "marker" set telling us that vector[[1]] is now "used". We iterate over subsequent elements in deck, setting them to the corresponding values in vector provided those values are not yet used. Once done with this iteration we go through the elements that have no values, assigning them in sequence the values that have not yet been assigned.

This notion of associating a list with repeats to a distinct shuffle has a clear drawback insofar as earlier elements are more likely to be assignedcorresponding values in vector. All the same this provides a reasonable way to make a "chromosome" vector containing repeats correspond to a permutation. Moreover one can see that any sensible mating process of two chromosomes will less drastically alter the objective function than would be the case in method 1, as the corresponding permutation now depends far less on overall ordering in the chromosomes. The advantage is thus that this method will be less in need of intricate crossover-probability parameter tuning.

getPerm = 
        Compile[{{vec, _Integer, 1}}, 
          Module[{p1, p2, len = Length[vec], k}, 
            p1 = p2 = Table[0, {len}]; 
Do[
              k = vec[[j]]; 
If[
                p2[[k]] == 0, 
p2[[k]] = 
                  j; 
p1[[j]] = 
                  k;], 
{j, len}]; 
k = 
              1; 
Do[
              If[p1[[j]] == 0, 
                While[p2[[k]] ≠ 
                    0, k++]; 
p1[[j]] = 
                  k; 
p2[[k]] = j], {j, 
                len}]; 
p1]];

QAP2[mat1_, mat2_, cp_, it_, sp_] := 
      Module[{len = Length[mat1], obfunc, vars, vv, nmin, vals, 
          constraints}, 
vars = 
          Array[vv, len]; 
constraints = 
          Prepend[Map[(1 ≤ # ≤ len) &, vars], 
            Element[vars, Integers]]; 
obfunc[
            vec : {__Integer}] := 
          Apply[Plus, 
            Flatten[mat1*
                permuteMatrix[mat2, 
                  getPerm[vec]]]]; 
{nmin, vals} = 
          NMinimize[{obfunc[vars], constraints}, vars, 
            Method -> {DifferentialEvolution, SearchPoints -> sp, 
                CrossProbability -> cp}, PostProcess -> False, 
            MaxIterations -> it]; 
{nmin, 
          getPerm[vars /. vals]}]

Though we suspect that tuning is less important for this method, we anyway show details of a "tuning run". We take several short runs with a given crossover probability parameter value, average them, and see what values seem to work best. We then try a longer run using one such value.

SeedRandom[111111];
Timing[
      vals = Table[{j, 
            Table[{min, perm} = QAP2[aa, bb, j/100, 200, 10]; min, {5}]}, {j, 
            5, 95, 5}]]

v2  =  
      Map[{#[[1]], Apply[Plus, #[[2]]]/5. } &,  vals]

From this we see that large values seem to work at least slightly better than other values. A similar test run over values in the range .90-.98 indicates that .96 might be best for our purposes.

SeedRandom[111111];
{min, perm} = QAP2[aa, bb,  .96, 20000, 200]

This gets us quite close to the global minimum with a scant two dozen lines of code. While it is mildly more complicated than method 1 above, it has one clear advantage. Both intuition and practical experience indicate that it is less in need of careful tuning of relatively obscure algorithm parameters. Also note that much of the code complexity is in associating a unique permutation to a given vector in our "search space". This part is likely to be problematic for any algorithm. The fact that we employ a genetic one makes it all the more so because we then need it to behave well with respect to a mating process, the internals of which are a black box. Hence a method that is relatively stable with respect to tuning parameters is all the more desirable.

Brett Champion has achieved even better results for this problem, using method 1 above. He has in fact found a permutation that give a function value of 3752. But the important point is that we have obtained a quite practical heuristic improvement in the objective function relative to random values, and this idea may apply even in cases where the absolute minimizer is not already known.

what lies beyond

This example comes from a family, the hardest of which is NUG30. The salient features are that the baseline and "typical" random permutations give values around 8000, and the minimizing permutation gives 6124. With some tuning of DifferentialEvolution parameters and reasonably long runs we have obtained values within 100 of that minimum.

Summary

We have seen several examples of discrete optimization problems and how one might tackle them in Mathematica. The underlying method we utilize is a member of the genetic family of optimization algorithms and we have seen that it is reasonably well adapted for these problems. We have discussed various ways to cast our examples so that simple Mathematica code using NMinimize may be applied. Indeed, we noticed that there are often very different ways to approach the same combinatorial optimization problem.

This uses technology presently under development and so it is to be expected that both the internals and interface might evolve over time. That said, it is clear from the examples that there is already some powerful capability for handling nontrivial combinatorial optimization problems. We hope that people find this functionality to be useful and look forward to developing it further.