I read an interesting paper by Zhang, Beltran-Royo and Ma on the excellent Optimization Online site this morning. As an exercise I wrote a Solver Foundation program to reproduce some of the results. In previous posts I showed how to write a C# program to solve Quadratic Assignment Problems (QAP). I wrote everything from scratch, including the Gilmore-Lawler bound computation. In this post (and the ones that follow) I am going to write everything using Solver Foundation Services in an effort to show how much easier it is!

If you are coming at this from the QAP angle, I recommend reading my previous QAP posts above. If you are coming at this from the Solver Foundation angle, look at page 4 of the paper (equations 2.12 – 2.15 in particular) and look at how easy it is to write down this rather complicated model using Solver Foundation Services. If you want a short Solver Foundation, check this out, or check out our new and improved documentation.

The paper talks about ways to use mixed integer linear programs (such as those built into Solver Foundation) to solve linearizations of QAP. One popular linearization is the one proposed by Kaufman and Broeckx. This formulation introduces n^2 continuous variables in addition to the standard integer variables x_ij (which represent the optimal assignment of facilities to locations). If you look at equations 2.12 – 2.15 on page 4, you will see that the goal function involves auxiliary real variables z. 2.13 is the most important constraint. 2.14 just enforces the non-negativity of the z variables, and we also need to make sure that the standard constraints on x hold: the row and column sums are equal to 1.

The code to compute the Kaufman-Broeckx linearization is below. I hardcoded a sample problem in A and B. You can replace it to read a standard problem from QAPLIB. Notice that the goals and constraints are basically transcriptions of the equations! The equations are complicated, but the way they map to the SFS APIs is straightforward. The other issue is how to use data binding to establish the values of the “q” and “a” parameters. This is done with the help of an auxiliary Cell class, and helper methods that return the “q” and “a” items. The SetBinding calls are pretty straightforward once you’ve got those.

If I were better at LINQ I could have simplified some of the code in GetCells and GetSums. But I’m not, so I didn’t. Viewing the output shows that the goal function is 50 – which appears to be correct (that’s the optimal value for the nug05 problem). Looking at the x values you can determine the optimal permutation according to the linearization.

using System;

using System.Collections.Generic;

using System.Linq;

using System.Text;

using Microsoft.SolverFoundation.Services;

namespace QapLinearization {

class Program {

// Here is the data for nug05.dat from QAPLIB.

private double[][] A = new double[][] {

new double[] { 0, 5, 2, 4, 1 },

new double[] { 5, 0, 3, 0, 2 },

new double[] { 2, 3, 0, 0, 0 },

new double[] { 4, 0, 0, 0, 5 },

new double[] { 1, 2, 0, 5, 0 }

};private double[][] B = new double[][] {

new double[] { 0, 1, 1, 2, 3 },

new double[] { 1, 0, 2, 1, 2 },

new double[] { 1, 2, 0, 1, 2 },

new double[] { 2, 1, 1, 0, 1 },

new double[] { 3, 2, 2, 1, 0 }

};

static void Main(string[] args) {

new Program();Program p =

p.KaufmanBroeckx();

}

public void KaufmanBroeckx() {

SolverContext context = SolverContext.GetContext();

Model model = context.CreateModel();

new Set(Domain.Integer, "N");Set N =

new Parameter(Domain.Real, "q", N, N, N, N);Parameter q =

"Value", "I", "J", "K", "L");q.SetBinding(GetCells(),

new Parameter(Domain.Real, "a", N, N);Parameter a =

"Value", "I", "J");a.SetBinding(GetSums(),

model.AddParameters(q, a);

new Decision(Domain.IntegerRange(0, 1), "x", N, N);Decision x =

new Decision(Domain.RealNonnegative, "z", N, N);Decision z =

model.AddDecisions(x, z);

"g", GoalKind.Minimize,model.AddGoal(

Model.Sum(Model.ForEach(N, i =>

Model.Sum(Model.ForEach(N, j => z[i, j]))

))

);

"c",model.AddConstraint(

Model.ForEach(N, i =>

Model.ForEach(N, j =>

z[i, j] >= Model.Sum(Model.ForEach(N, k =>

Model.Sum(Model.ForEach(N, l => q[i, j, k, l] * x[k, l] - a[i, j] * (1 - x[i, j])))

))

))

);

"rowsum",model.AddConstraint(

Model.ForEach(N, i =>

Model.Sum(Model.ForEach(N, j => x[i, j])) == 1

));

"colsum",model.AddConstraint(

Model.ForEach(N, j =>

Model.Sum(Model.ForEach(N, i => x[i, j])) == 1

));

var solution = context.Solve();

Console.WriteLine(solution.GetReport());

}

// This is a helper class for data binding.

private class Cell {

public int I { get; set; }

public int J { get; set; }

public int K { get; set; }

public int L { get; set; }

public double Value { get; set; }

}

// This returns the q_ijkl as described on page 2 of the paper.

private IEnumerable<Cell> GetCells() {

for (int i = 0; i < A.Length; i++) {

for (int j = 0; j < A[i].Length; j++) {

for (int k = 0; k < B.Length; k++) {

for (int l = 0; l < B[k].Length; l++) {

yield return new Cell { I = i, J = j, K = k, L = l, Value = A[i][k] * B[j][l] };

}

}

}

}

}

// This returns the a_ij as described on page 4 of the paper.

private IEnumerable<Cell> GetSums() {

for (int i = 0; i < A.Length; i++) {

for (int j = 0; j < A[i].Length; j++) {

double sum = 0;

for (int k = 0; k < B.Length; k++) {

for (int l = 0; l < B[k].Length; l++) {

if (B[k][j] != 0) {

sum += A[i][k] * B[j][l];

}

}

}yield return new Cell { I = i, J = j, Value = sum };

}

}

}

}

}