Facility Location (FLP)

Principles learned

  • Create a generic model that uses data

  • Define the actual set of decision variables

  • Use non-linearities: ternary operator, minimizing a sum of “min”

Problem

../_images/facilitylocation.svg

The facility location problem is defined as follows: given a set of N locations and a transportation cost W between each pair of locations, select a subset S of p locations, which become facilities, in order to minimize transportation costs. The transportation cost for a location is defined as the distance to the closest facility, that is, the closest location in S. Thus, it aims to provide an optimal placement of facilities to minimize transportation costs. This problem is also known as the p-median problem.

Download the example


Data

The data files pmed1.in, pmed2.in, …, pmed40.in provided come from the OR-LIB. They are the 40 test problems from Table 2 of J.E.Beasley “A note on solving large p-median problems” European Journal of Operational Research 21 (1985) 270-273.

Each instance file contains the number of locations n, the number of edges in the original instance, the number of facilities to select p and the distance matrix computed from the original file using Floyd’s algorithm.

Program

The general idea behind the model is that once one knows the facilities that are in the subset S, one can deduce the distance between each location and the closest facility in S. There is no need to actually know which facility is the closest in the model, just the closest distance. It can still be obtained through post-processing if one needs to know this information to print out the results.

Thus, the only decision variables are the booleans x[i], equal to 1 if the ith location is in the subset S. The expressions costs[i][j] are created to store the transportation cost between location i and j. This cost is:

  • edge weight between i and j if j is in the subset S, i.e is a facility

  • infinity otherwise (represented by 2 times the maximum edge weight)

The transportation cost between each location i and the closest facility in S is represented by cost[i]: it is the minimum value among costs[i][j] for all j. The objective is then to minimize the sum of these transportation costs.

Execution:
hexaly facility_location.hxm inFileName=instances/pmed1.in [hxTimeLimit=] [solFileName=]
use io;

/* Read instance data */
function input() {
    local usage = "Usage: hexaly facility_location.hxm "
            + "inFileName=inputFile [solFileName=outputFile] [hxTimeLimit=timeLimit]";

    if (inFileName == nil) throw usage;

    local inFile = io.openRead(inFileName);
    N = inFile.readInt();
    inFile.readInt(); // Skip number of edges
    p = inFile.readInt();

    w[i in 0...N][j in 0...N] = inFile.readInt();
    wmax = max[i in 0...N][j in 0...N] (w[i][j]);
}

/* Declare the optimization model */
function model() {
    // One variable for each location: 1 if facility, 0 otherwise
    x[i in 0...N] <- bool();

    // No more than p locations are selected to be facilities
    constraint sum[i in 0...N] (x[i]) <= p;

    // Costs between location i and j is w[i][j] if j is a facility or 2 * wmax if not
    costs[i in 0...N][j in 0...N] <- x[j] ? w[i][j] : 2 * wmax;

    // Cost between location i and the closest facility
    cost[i in 0...N] <- min[j in 0...N] (costs[i][j]);

    // Minimize the total cost
    totalCost <- sum[i in 0...N] (cost[i]);
    minimize totalCost;
}

/* Parametrize the solver */
function param() {
    if (hxTimeLimit == nil) hxTimeLimit = 20;
}

/* Write the solution in a file with the following format: 
 * - value of the objective
 * - indices of the facilities (between 0 and N-1) */
function output() {
    if (solFileName == nil) return; 
    local solFile = io.openWrite(solFileName);
    solFile.println(totalCost.value);
    for [i in 0...N : x[i].value == 1]
        solFile.print(i, " ");
    solFile.println("");
}
Execution (Windows)
set PYTHONPATH=%HX_HOME%\bin\python
python facility_location.py instances\pmed1.in
Execution (Linux)
export PYTHONPATH=/opt/hexaly_13_0/bin/python
python facility_location.py instances/pmed1.in
import hexaly.optimizer
import sys

def read_integers(filename):
    with open(filename) as f:
        return [int(elem) for elem in f.read().split()]

#
# Read instance data
#
def read_instance(instance_file):
    file_it = iter(read_integers(instance_file))
    # Number of locations
    N = next(file_it)
    next(file_it) # Skip number of edges
    # Size of the subset S of facilities
    p = next(file_it)

    # w: Weight matrix of the shortest path between locations
    # wmax: Maximum distance between two locations
    wmax = 0
    w = [None] * N
    for i in range(N):
        w[i] = [None] * N
        for j in range(N):
            w[i][j] = next(file_it)
            if w[i][j] > wmax:
                wmax = w[i][j]

    return N, p, wmax, w

def main(instance_file, output_file, time_limit):
    N, p, wmax, w = read_instance(instance_file)

    with hexaly.optimizer.HexalyOptimizer() as optimizer:
        #
        # Declare the optimization model
        #
        m = optimizer.model

        # One variable for each location: 1 if facility, 0 otherwise
        x = [m.bool() for i in range(N)]

        # No more than p locations are selected to be facilities
        opened_locations = m.sum(x[i] for i in range(N))
        m.constraint(opened_locations <= p)

        # Costs between location i and j is w[i][j] if j is a facility or 2 * wmax if not
        costs = [None] * N
        for i in range(N):
            costs[i] = [None] * N
            for j in range(N):
                costs[i][j] = m.iif(x[j], w[i][j], 2 * wmax)

        # Cost between location i and the closest facility
        cost = [None] * N
        for i in range(N):
            cost[i] = m.min(costs[i][j] for j in range(N))

        # Minimize the total cost
        total_cost = m.sum(cost[i] for i in range(N))
        m.minimize(total_cost)

        m.close()

        # Parameterize the optimizer
        optimizer.param.time_limit = time_limit

        optimizer.solve()

        #
        # Write the solution in a file with the following format:
        # - value of the objective
        # - indices of the facilities (between 0 and N-1)
        #
        if output_file is not None:
            with open(output_file, 'w') as f:
                f.write("%d\n" % total_cost.value)
                for i in range(N):
                    if x[i].value == 1:
                        f.write("%d " % i)
                f.write("\n")

if __name__ == '__main__':
    if len(sys.argv) < 2:
        print("Usage: python facility_location.py instance_file [output_file] [time_limit]")
        sys.exit(1)

    instance_file = sys.argv[1]
    output_file = sys.argv[2] if len(sys.argv) >= 3 else None
    time_limit = int(sys.argv[3]) if len(sys.argv) >= 4 else 20
    main(instance_file, output_file, time_limit)
Compilation / Execution (Windows)
cl /EHsc facility_location.cpp -I%HX_HOME%\include /link %HX_HOME%\bin\hexaly130.lib
facility_location instances\pmed1.in
Compilation / Execution (Linux)
g++ facility_location.cpp -I/opt/hexaly_13_0/include -lhexaly130 -lpthread -o facility_location
./facility_location instances/pmed1.in
#include "optimizer/hexalyoptimizer.h"
#include <fstream>
#include <iostream>
#include <sstream>
#include <vector>

using namespace hexaly;
using namespace std;

class Facilitylocation {
public:
    // Number of locations
    int N;

    // Size of the subset S of facilities
    int p;

    // Weight matrix of the shortest path between locations
    vector<vector<int>> w;

    // Maximum distance between two locations
    int wmax;

    // Hexaly Optimizer
    HexalyOptimizer optimizer;

    // Decisions variables
    vector<HxExpression> x;

    // Objective
    HxExpression totalCost;

    // Read instance data
    void readInstance(const string& fileName) {
        ifstream infile;
        infile.exceptions(ifstream::failbit | ifstream::badbit);
        infile.open(fileName.c_str());

        int tmp;
        infile >> N;
        infile >> tmp; // Skip number of edges
        infile >> p;

        w.resize(N);
        wmax = 0;
        for (int i = 0; i < N; ++i) {
            w[i].resize(N);
            for (int j = 0; j < N; ++j) {
                infile >> w[i][j];
                if (w[i][j] > wmax) {
                    wmax = w[i][j];
                }
            }
        }
    }

    // Declare the optimization model
    void solve(int limit) {
        HxModel m = optimizer.getModel();

        // One variable for each location: 1 if facility, 0 otherwise
        x.resize(N);
        for (int i = 0; i < N; ++i)
            x[i] = m.boolVar();

        // No more than p locations are selected to be facilities
        HxExpression openedLocations = m.sum(x.begin(), x.end());
        m.constraint(openedLocations <= p);

        // Costs between location i and j is w[i][j] if j is a facility or 2 * wmax if not
        vector<vector<HxExpression>> costs(N);
        for (int i = 0; i < N; ++i) {
            costs[i].resize(N);
            for (int j = 0; j < N; ++j) {
                costs[i][j] = m.iif(x[j], w[i][j], 2 * wmax);
            }
        }

        // Cost between location i and the closest facility
        vector<HxExpression> cost(N);
        for (int i = 0; i < N; ++i) {
            cost[i] = m.min(costs[i].begin(), costs[i].end());
        }

        // Minimize the total cost
        totalCost = m.sum(cost.begin(), cost.end());
        m.minimize(totalCost);
        m.close();

        // Parameterize the optimizer
        optimizer.getParam().setTimeLimit(limit);

        optimizer.solve();
    }

    /* Write the solution in a file with the following format: 
     * - value of the objective
     * - indices of the facilities (between 0 and N-1) */
    void writeSolution(const string& fileName) {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());

        outfile << totalCost.getValue() << endl;
        for (int i = 0; i < N; ++i) {
            if (x[i].getValue() == 1)
                outfile << i << " ";
        }
        outfile << endl;
    }
};

int main(int argc, char** argv) {
    if (argc < 2) {
        cerr << "Usage: facility_location inputFile [outputFile] [timeLimit] " << endl;
        return 1;
    }

    const char* instanceFile = argv[1];
    const char* solFile = argc > 2 ? argv[2] : NULL;
    const char* strTimeLimit = argc > 3 ? argv[3] : "20";

    try {
        Facilitylocation model;
        model.readInstance(instanceFile);
        model.solve(atoi(strTimeLimit));
        if (solFile != NULL)
            model.writeSolution(solFile);
        return 0;
    } catch (const exception& e) {
        cerr << "An error occurred: " << e.what() << endl;
        return 1;
    }
}
Compilation / Execution (Windows)
copy %HX_HOME%\bin\Hexaly.NET.dll .
csc Facilitylocation.cs /reference:Hexaly.NET.dll
Facilitylocation instances\pmed1.in
using System;
using System.IO;
using System.Collections.Generic;
using Hexaly.Optimizer;

public class FacilityLocation : IDisposable
{
    // Number of locations
    int N;

    // Size of the subset S of facilities
    int p;

    // Weight matrix of the shortest path beween locations
    int[][] w;

    // Maximum distance between two locations
    int wmax;

    // Hexaly Optimizer
    HexalyOptimizer optimizer;

    // Decision variables
    HxExpression[] x;

    // Objective
    HxExpression totalCost;

    // List of facilities
    List<int> solution;

    public FacilityLocation()
    {
        optimizer = new HexalyOptimizer();
    }

    // Read instance data
    public void ReadInstance(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            var tokens = input.ReadLine().Split(' ');
            N = int.Parse(tokens[0]);
            p = int.Parse(tokens[2]);

            w = new int[N][];
            wmax = 0;
            for (int i = 0; i < N; ++i)
            {
                tokens = input.ReadLine().Split(' ');
                w[i] = new int[N];
                for (int j = 0; j < N; ++j)
                {
                    w[i][j] = int.Parse(tokens[j]);
                    if (w[i][j] > wmax)
                        wmax = w[i][j];
                }
            }
        }
    }

    public void Dispose()
    {
        if (optimizer != null)
            optimizer.Dispose();
    }

    // Declare the optimization model
    public void Solve(int limit)
    {
        optimizer = new HexalyOptimizer();
        HxModel model = optimizer.GetModel();
        x = new HxExpression[N];

        // One variable for each location: 1 if facility, 0 otherwise
        for (int i = 0; i < N; ++i)
            x[i] = model.Bool();

        // No more than p locations are selected to be facilities
        HxExpression openedLocations = model.Sum(x);
        model.Constraint(openedLocations <= p);

        // Costs between location i and j is w[i][j] if j is a facility or 2 * wmax if not
        HxExpression[][] costs = new HxExpression[N][];
        for (int i = 0; i < N; ++i)
        {
            costs[i] = new HxExpression[N];
            for (int j = 0; j < N; ++j)
                costs[i][j] = model.If(x[j], w[i][j], 2 * wmax);
        }

        // Cost between location i and the closest facility
        HxExpression[] cost = new HxExpression[N];
        for (int i = 0; i < N; ++i)
            cost[i] = model.Min(costs[i]);

        // Minimize the total cost
        totalCost = model.Sum(cost);
        model.Minimize(totalCost);

        model.Close();

        // Parameterize the optimizer
        optimizer.GetParam().SetTimeLimit(limit);
        optimizer.Solve();

        solution = new List<int>();
        for (int i = 0; i < N; ++i)
            if (x[i].GetValue() == 1)
                solution.Add(i);
    }

    /* Write the solution in a file with the following format:
     * - value of the objective
     * - indices of the facilities (between 0 and N-1) */
    public void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine(totalCost.GetValue());
            for (int i = 0; i < N; ++i) {
                if (x[i].GetValue() == 1)
                    output.Write(i + " ");
            }
            output.WriteLine();
        }
    }

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: Facilitylocation inputFile [outputFile] [timeLimit]");
            Environment.Exit(1);
        }
        string instanceFile = args[0];
        string outputFile = args.Length > 1 ? args[1] : null;
        string strTimeLimit = args.Length > 2 ? args[2] : "20";

        using (FacilityLocation model = new FacilityLocation())
        {
            model.ReadInstance(instanceFile);
            model.Solve(int.Parse(strTimeLimit));
            if (outputFile != null)
                model.WriteSolution(outputFile);
        }
    }
}
Compilation / Execution (Windows)
javac Facilitylocation.java -cp %HX_HOME%\bin\hexaly.jar
java -cp %HX_HOME%\bin\hexaly.jar;. Facilitylocation instances\pmed1.in
Compilation / Execution (Linux)
javac Facilitylocation.java -cp /opt/hexaly_13_0/bin/hexaly.jar
java -cp /opt/hexaly_13_0/bin/hexaly.jar:. Facilitylocation instances/pmed1.in
import java.util.*;
import java.io.*;
import com.hexaly.optimizer.*;

public class Facilitylocation {
    // Number of locations
    private int N;

    // Size of the subset S of facilities
    private int p;

    // Weight matrix of the shortest path between locations
    private int[][] w;

    // Maximum distance between two locations
    private int wmax;

    // Hexaly Optimizer.
    private final HexalyOptimizer optimizer;

    // Decisions variables
    private HxExpression[] x;

    // Objective
    private HxExpression totalCost;

    private Facilitylocation(HexalyOptimizer optimizer) {
        this.optimizer = optimizer;
    }

    // Read instance data
    private void readInstance(String fileName) throws IOException {
        try (Scanner input = new Scanner(new File(fileName))) {
            N = input.nextInt();
            input.nextInt(); // Skip number of edges
            p = input.nextInt();

            w = new int[N][N];
            wmax = 0;
            for (int i = 0; i < N; ++i) {
                for (int j = 0; j < N; ++j) {
                    w[i][j] = input.nextInt();
                    if (w[i][j] > wmax)
                        wmax = w[i][j];
                }
            }
        }
    }

    private void solve(int limit) {
        // Declare the optimization model
        HxModel model = optimizer.getModel();

        // One variable for each location: 1 if facility, 0 otherwise
        x = new HxExpression[N];
        for (int i = 0; i < N; ++i) {
            x[i] = model.boolVar();
        }

        // No more than p locations are selected to be facilities
        HxExpression openedLocations = model.sum(x);
        model.constraint(model.leq(openedLocations, p));

        // Costs between location i and j is w[i][j] if j is a facility or 2 * wmax if not
        HxExpression[][] costs = new HxExpression[N][N];
        for (int i = 0; i < N; ++i) {
            for (int j = 0; j < N; ++j) {
                costs[i][j] = model.iif(x[j], w[i][j], 2 * wmax);
            }
        }

        // Cost between location i and the closest facility
        HxExpression[] cost = new HxExpression[N];
        for (int i = 0; i < N; ++i) {
            cost[i] = model.min(costs[i]);
        }

        // Minimize the total cost
        totalCost = model.sum(cost);
        model.minimize(totalCost);

        model.close();

        // Parameterize the optimizer
        optimizer.getParam().setTimeLimit(limit);
        optimizer.solve();
    }

    /* Write the solution in a file with the following format:
     * - value of the objective
     * - indices of the facilities (between 0 and N-1) */
    private void writeSolution(String fileName) throws IOException {
        try (PrintWriter output = new PrintWriter(fileName)) {
            output.println(totalCost.getValue());
            for (int i = 0; i < N; ++i) {
                if (x[i].getValue() == 1)
                    output.print(i + " ");
            }
        }
    }

    public static void main(String[] args) {
        if (args.length < 1) {
            System.err.println("Usage: java Facilitylocation inputFile [outputFile] [timeLimit]");
            System.exit(1);
        }

        String instanceFile = args[0];
        String outputFile = args.length > 1 ? args[1] : null;
        String strTimeLimit = args.length > 2 ? args[2] : "20";

        try (HexalyOptimizer optimizer = new HexalyOptimizer()) {
            Facilitylocation model = new Facilitylocation(optimizer);
            model.readInstance(instanceFile);
            model.solve(Integer.parseInt(strTimeLimit));
            if (outputFile != null) {
                model.writeSolution(outputFile);
            }
        } catch (Exception ex) {
            System.err.println(ex);
            ex.printStackTrace();
            System.exit(1);
        }
    }
}