# P-median¶

## 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¶

The pmedian problem is defined as follows: given a set I={1...n} of locations and a transportation cost W between each pair of locations. Select a subset S of p locations minimizing the sum of the distances between each location and the closest one in S.

## 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 locations to select p and the distance matrix computed from the original file using Floyd’s algorithm.

## Program¶

In this problem, once you know the locations that are in the subset S, you can deduce which location in S is the closest to each location.

The general idea behind the model is than once one know the locations that are in the subset S, one can deduce the distance between each location and the closest location in S. There is no need to actually know which location is the closest in the model, just the closest distance. It can still be obtained through post-processing if one need 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
• infinity otherwise (represented by 2 times the maximum edge weight)

The transportation cost between each location i and the closest location 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:
localsolver pmedian.lsp inFileName=instances/pmed1.in [lsTimeLimit=] [solFileName=]
```/********** pmedian.lsp **********/

use io;

function input(){
local usage = "Usage: localsolver pmedian.lsp "
+ "inFileName=inputFile [solFileName=outputFile] [lsTimeLimit=timeLimit]";

if (inFileName == nil) throw usage;

wmax = 0;

for [i in 1..N][j in 1..N] {
if(w[i][j] > wmax) wmax = w[i][j];
}
}

/* Declares the optimization model. */
function model(){
// One variable for each location
x[1..N] <- bool();

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

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

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

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

/* Parameterizes the solver. */
function param(){
if(lsTimeLimit == nil) lsTimeLimit = 10;
}

/* Writes the solution in a file following the following format:
* - value of the objective
* - indiced of the selected locations (between 0 and N-1) */
function output() {
if(solFileName == nil) return;
local solFile = io.openWrite(solFileName);
solFile.println(totalCost.value);
for [i in 1..N : x[i].value == 1] {
solFile.print(i-1, " ");
}
solFile.println("");
}
```
Execution (Windows)
set PYTHONPATH=%LS_HOME%\bin\python27\
python pmedian.py instances\pmed1.in
Execution (Linux)
export PYTHONPATH=/opt/localsolver_XXX/bin/python27/
python pmedian.py instances/pmed1.in
```########## pmedian.py ##########

import localsolver
import sys

if len(sys.argv) < 2:
print ("Usage: python pmedian.py inputFile [outputFile] [timeLimit]")
sys.exit(1)

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

with localsolver.LocalSolver() as ls:

#
#

# Number of locations
N = file_it.next()
# Number of edges between locations
E = file_it.next()
# Size of the subset S of locations
p = file_it.next()

# 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] = file_it.next()
if(w[i][j] >  wmax):
wmax = w[i][j];

#
# Declares the optimization model
#
m = ls.model

# Decision variables
x = [m.bool() for i in range(N)]

# No more than p locations are selected
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 selected in S 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 selected location
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()

#
# Parameterizes the solver
#
if len(sys.argv) >= 4: ls.create_phase().time_limit = int(sys.argv[3])
else: ls.create_phase().time_limit = 20

ls.solve()

#
# Writes the solution in a file following the following format:
# - value of the objective
# - indiced of the selected locations (between 0 and N-1) */
#
if len(sys.argv) >= 3:
with open(sys.argv[2], '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")
```
Compilation / Execution (Windows)
cl /EHsc pmedian.cpp -I%LS_HOME%\include /link %LS_HOME%\bin\localsolver.dll.lib
pmedian instances\pmed1.in
Compilation / Execution (Linux)
g++ pmedian.cpp -I/opt/localsolver_XXX/include -llocalsolver -lpthread -o pmedian
./pmedian instances/pmed1.in
```//********* pmedian.cpp *********

#include <iostream>
#include <sstream>
#include <fstream>
#include <vector>
#include "localsolver.h"

using namespace localsolver;
using namespace std;

class Pmedian {
public:
// Number of locations.
lsint N;

// Number of edges between locations.
lsint E;

// Size of the subset S of locations.
lsint p;

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

// Maximum distance between two locations.
lsint wmax;

// LocalSolver.
LocalSolver localsolver;

// Decisions variables
vector<LSExpression> x;

// Objective
LSExpression totalCost;

// Vector of selected locations
vector<int> solution;

void readInstance(const string & fileName) {
ifstream infile;
infile.open(fileName.c_str());

infile >> N;
infile >> E;
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];
}
}
}
}

// Declares the optimization model.
void solve(int limit){
LSModel m = localsolver.getModel();

// Decision variables
x.resize(N);
for(int i = 0; i < N; i++)
x[i] = m.boolVar();

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

// Costs between location i and j is w[i,j] if j is selected in S or 2*wmax if not.
vector<vector<LSExpression> > 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 selected location.
vector<LSExpression> 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();

// Parameterizes the solver.
LSPhase phase = localsolver.createPhase();
phase.setTimeLimit(limit);
localsolver.solve();

solution.clear();
for(int i = 0; i < N; i++){
if(x[i].getValue() == 1){
solution.push_back(i);
}
}
}

// Writes the solution in a file following the following format:
// each line contains the index of a selected location (between 0 and N-1)
void writeSolution(const string& fileName){
ofstream outfile;
outfile.open(fileName.c_str());

outfile << totalCost.getValue() << endl;
for(int i = 0; i < solution.size(); i++)
outfile << solution[i] << " ";
outfile << endl;
}
};

int main(int argc, char** argv){
if(argc < 2){
cerr << "Usage: pmedian 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 {
Pmedian model;
model.solve(atoi(strTimeLimit));
if(solFile != NULL) model.writeSolution(solFile);
return 0;
} catch (const exception& e){
cerr << "Error occurred: " << e.what() << endl;
return 1;
}
}
```
Compilation/Execution (Windows)
copy %LS_HOME%\bin\*net.dll .
csc Pmedian.cs /reference:localsolvernet.dll
Pmedian instances\pmed1.in
```/********** Pmedian.cs **********/

using System;
using System.IO;
using System.Collections.Generic;
using localsolver;

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

// Number of edges between locations.
int E;

// Size of the subset S of locations.
int p;

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

// Maximum distance between two locations.
int wmax;

// LocalSolver.
LocalSolver localsolver;

// Decision variables.
LSExpression[] x;

// Objective.
LSExpression totalCost;

// List of selected locations
List<int> solution;

public Pmedian()
{
localsolver = new LocalSolver();
}

{
{
N = int.Parse(tokens[0]);
E = int.Parse(tokens[1]);
p = int.Parse(tokens[2]);

w = new int[N][];

wmax = 0;

for (int i = 0; i < N; i++)
{
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 (localsolver != null)
localsolver.Dispose();
}

// Declares the optimization model.
public void Solve(int limit)
{
localsolver = new LocalSolver();
LSModel model = localsolver.GetModel();
x = new LSExpression[N];

// Boolean decisions.
for (int i = 0; i < N; i++)
x[i] = model.Bool();

// No more than p locations are selected.
LSExpression openedLocations = model.Sum();
for (int i = 0; i < N; i++)

model.Constraint(openedLocations <= p);

// Costs between location i and j is w[i,j] if j is selected in S or 2*wmax if not.
LSExpression[,] costs = new LSExpression[N, N];
for (int i = 0; i < N; i++)
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 selected location.
LSExpression[] cost = new LSExpression[N];
for (int i = 0; i < N; i++)
{
cost[i] = model.Min();
for (int j = 0; j < N; j++)
{
}
}

// Minimize the total cost.
totalCost = model.Sum();
for (int i = 0; i < N; i++)

model.Minimize(totalCost);

model.Close();

// Parameterizes the solver.
LSPhase phase = localsolver.CreatePhase();
phase.SetTimeLimit(limit);
localsolver.Solve();

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

// Writes the solution in a file following the following format:
// - value of the objective
// - indiced of the selected locations (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 < solution.Count; ++i)
output.Write(solution[i] + " ");
output.WriteLine();
}
}

public static void Main(string[] args)
{
if (args.Length < 1)
{
Console.WriteLine("Usage: Pmedian 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 (Pmedian model = new Pmedian())
{
model.Solve(int.Parse(strTimeLimit));
if (outputFile != null)
model.WriteSolution(outputFile);
}
}
}
```
Compilation / Execution (Windows)
javac Pmedian.java -cp %LS_HOME%\bin\localsolver.jar
java -cp %LS_HOME%\bin\localsolver.jar;. Pmedian instances\pmed1.in
Compilation/Execution (Linux)
javac Pmedian.java -cp /opt/localsolver_XXX/bin/localsolver.jar
java -cp /opt/localsolver_XXX/bin/localsolver.jar:. Pmedian instances/pmed1.in
```/********** Pmedian.java **********/

import java.util.*;
import java.io.*;
import localsolver.*;

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

// Number of edges between locations.
private int E;

// Size of the subset S of locations.
private int p;

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

// Maximum distance between two locations.
private int wmax;

// LocalSolver.
private LocalSolver localsolver;

// Decisions variables
private LSExpression[] x;

// Objective
private LSExpression totalCost;

// List of selected locations
private List<Integer> solution;

private void readInstance(String fileName) throws IOException {
try(Scanner input = new Scanner(new File(fileName))) {
N = input.nextInt();
E = input.nextInt();
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];
}
}
}
}

// Declares the optimization model.
private void solve(int limit) {
localsolver = new LocalSolver();
LSModel model = localsolver.getModel();

// Boolean decisions
x = new LSExpression[N];
for (int i = 0; i < N; i++) {
x[i] = model.boolVar();
}

// No more than p locations are selected
LSExpression openedLocations = model.sum();
for (int i = 0; i < N; i++) {
}

model.constraint(model.leq(openedLocations, p));

// Costs between location i and j is w[i][j] if j is selected in S or 2*wmax if not
LSExpression[][] costs = new LSExpression[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 selected location
LSExpression[] cost = new LSExpression[N];
for (int i = 0; i < N; i++) {
cost[i] = model.min();
for (int j = 0; j < N; j++) {
}
}

// Minimize the total cost
totalCost = model.sum();
for (int i = 0; i < N; i++) {
}
model.minimize(totalCost);

model.close();

// Parameterizes the solver.
LSPhase phase = localsolver.createPhase();
phase.setTimeLimit(limit);
localsolver.solve();

solution = new ArrayList<Integer>();
for (int i = 0; i < N; i++) {
if (x[i].getValue() == 1)
}
}

// Writes the solution in a file following the following format:
// - value of the objective
// - indiced of the selected locations (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 < solution.size(); i++) {
output.print(solution.get(i) + " ");
}
output.println();
}
}

public static void main(String[] args) {
if (args.length < 1) {
System.err.println("Usage: java Pmedian 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] : "10";

try {
Pmedian model = new Pmedian();