⚠

EXPERIMENTAL: These features are experimental and should not be used in production systems.

Mathematical Optimization: The Facility Location Problem

Imagine that you are the logistics manager at a large retail company. Your company is planning to open a new distribution center to improve delivery service to its customers. You have identified three potential sites for the new distribution center, labeled as "F1", "F2", and "F3".

You also have information about three major customer locations in the cities "C1", "C2", and "C3". Each of these customer locations has a demand level - this could be the amount of goods they require in a certain period, or the frequency of deliveries they need.

However, opening a distribution center is a costly affair. You can only afford to open one center at this time. You also have information about the distance from each potential distribution center to each customer location. It is reasonable to assume that the transportation cost will be directly proportional to this distance. The question then becomes: which location should you choose for the new distribution center to minimize the overall transportation cost?

This scenerio is called the facility location problem, where the goal is to select a single facility location that minimizes the total weighted transportation cost to all customer locations. The weights here represent the demand at each customer location.

The facility location problem is best modeled and solved as a mixed-integer linear programming problem (MIP), and so, we can use Rel’s mathematical optimization library to solve it.

To begin, we will first need to import the mathopt:Solver:Solve DSL, which allows us to define the objective function and constraints using a set of mathematical operators.

// model
 
// Import the operators from mathopt:Solver:Solve DSL.
// These operators will be used to construct the model constraints and the model objective function.
with mathopt:Solver:Solve use sum, foreach, +, -, *, ≼, ≽, /, =

Next suppose that our facilities, customers, distances, and demands are specified by the following data:

// model
 
// Define the data.
module data
    def facilities   = "F1" ; "F2" ; "F3"
    def customers    = "C1" ; "C2" ; "C3"
    def distances = {
        ("F1", "C1", 10.0) ; ("F1", "C2", 20.0) ; ("F1", "C3", 30.0) ;
        ("F2", "C1", 20.0) ; ("F2", "C2", 10.0) ; ("F2", "C3", 20.0) ;
        ("F3", "C1", 30.0) ; ("F3", "C2", 20.0) ; ("F3", "C3", 10.0)
    }
    def demands      = { ("C1", 5.0) ; ("C2", 6.0) ; ("C3", 4.0) }
end
 
// Bring the data into scope.
with data use facilities, customers, distances, demands

Now that we have imported the mathopt:Solver:Solve operators and defined our data, we next construct a optimization model consisting of variables, constraints, and an objective function. We begin by first defining our model variables:

// model
 
// Define the model variables.
def model:variables = rel:mathopt:variables[
        {:y, "binary", (facilities)};
        {:x, "continuous", (facilities, customers)};
        {:x, "lower_bound", 0.0};
        {:x, "upper_bound", 1.0}
    ]
// Bring the model variables into scope.
with model:variables use x, y

In this case, y[f] is a binary variable that indicates whether facility f is selected (1) or not (0), and x[f, c] is a continuous variable that represents the fraction of customer c’s demand that is served by facility f.

Next we define our model objective function. In this case, the objective is to minimize the total weighted transportation cost, which is the sum of the distances from the chosen facility to the customer locations, each multiplied by the proportion of demand it serves and the demand itself.

// model
 
// Define the sense of the optimization and the objective function.
def model:minimize = sum[distances[f, c] * demands[c] * x[f, c] for f in facilities, c in customers]

Next we define the model constraints. In this case, we have three constraints:

Facility Selection Constraint: This constraint ensures that only one facility location is selected. In other words, the sum of the binary variables y[f] for all facilities f should be equal to 1.

// model
 
// Define the facility selection constraint.
// This ensures that only one facility location is selected.
def model:constraints:facility_selection = sum[y[f] for f in facilities] = 1

Demand Constraints: These constraints ensure that each customer’s demand is fully satisfied. This means the sum of the continuous variables x[f, c] for all facilities f should be equal to 1 for each customer c.

// model
 
// Define the demand constraints.
// This ensures that each customer's demand is fully satisfied.
def model:constraints:demand = foreach[c in customers : sum[x[f, c] for f in facilities] = 1]

Linking Constraints: These constraints ensure that if a facility is not selected, it cannot serve any customers. In other words, the continuous variable x[f, c] should be less than or equal to the binary variable y[f] for each facility f and customer c.

// model
 
// Define the linking constraints.
// This ensures that if a facility is not selected, it cannot serve any customers.
def model:constraints:linking = foreach[f in facilities, c in customers : x[f, c] ≼ y[f]]

This formulation is a classical mixed-integer linear programming (MIP) problem, which can be efficiently solved using MIP solvers. Rel’s mathematical optimization library provides a convenient way to define and solve such problems in a high-level, declarative manner. The following line of code invokes the rel:mathopt:optimize relation which optimizes the model that we’ve defined.

// model
 
// Optimize the model.
def result = rel:mathopt:optimize[model, {}]

Once the model is solved, we can extract the results. The following line of code uses the rel:mathopt:extract_result relation to retrieve the solution of the optimization problem, including the optimal values of the decision variables and the optimal objective value.

// read query
 
// Extract the result of the optimization process.
def output = rel:mathopt:extract_result[result, model:variables]

Full Code Implementation

The full code implementation of the facility location problem is shown below.

// model
 
// Import the mathopt library DSL operators.
with mathopt:Solver:Solve use sum, foreach, +, -, *, ≼, ≽, /, =
 
// Define the data.
module data
    def facilities   = "F1" ; "F2" ; "F3"
    def customers    = "C1" ; "C2" ; "C3"
    def distances = {
        ("F1", "C1", 10.0) ; ("F1", "C2", 20.0) ; ("F1", "C3", 30.0) ;
        ("F2", "C1", 20.0) ; ("F2", "C2", 10.0) ; ("F2", "C3", 20.0) ;
        ("F3", "C1", 30.0) ; ("F3", "C2", 20.0) ; ("F3", "C3", 10.0)
    }
    def demands      = { ("C1", 5.0) ; ("C2", 6.0) ; ("C3", 4.0) }
end
// Bring the data into scope.
with data use facilities, customers, distances, demands
 
// Define the model variables.
def model:variables = rel:mathopt:variables[
        {:y, "binary", (facilities)};
        {:x, "continuous", (facilities, customers)};
        {:x, "lower_bound", 0.0};
        {:x, "upper_bound", 1.0}
    ]
// Bring the model variables into scope.
with model:variables use x, y
 
// Define the sense of the optimization problem and the objective function.
def model:minimize = sum[distances[f, c] * demands[c] * x[f, c] for f in facilities, c in customers]
 
// Define the facility selection constraint.
// This ensures that only one facility location is selected.
def model:constraints:facility_selection = sum[y[f] for f in facilities] = 1
 
// Define the demand constraints.
// This ensures that each customer's demand is fully satisfied.
def model:constraints:demand = foreach[c in customers : sum[x[f, c] for f in facilities] = 1]
 
// Define the linking constraints.
// This ensures that if a facility is not selected, it cannot serve any customers.
def model:constraints:linking = foreach[f in facilities, c in customers : x[f, c] ≼ y[f]]
 
// Optimize the model.
def result = rel:mathopt:optimize[model, {}]
 
// Extract the result of the optimization process.
def output = rel:mathopt:extract_result[result, model:variables]

Alternatively, you can wrap the model in a module that keeps the imported DSL operators local to the module. This can be beneficial when you want to use Rel base operators with the same name as the DSL operators.

// model
 
// Define the data.
module data
    def facilities   = "F1" ; "F2" ; "F3"
    def customers    = "C1" ; "C2" ; "C3"
    def distances = {
        ("F1", "C1", 10.0) ; ("F1", "C2", 20.0) ; ("F1", "C3", 30.0) ;
        ("F2", "C1", 20.0) ; ("F2", "C2", 10.0) ; ("F2", "C3", 20.0) ;
        ("F3", "C1", 30.0) ; ("F3", "C2", 20.0) ; ("F3", "C3", 10.0)
    }
    def demands      = { ("C1", 5.0) ; ("C2", 6.0) ; ("C3", 4.0) }
end
 
// Define the optimization model in a module.
module model
    // Import the mathopt library DSL operators.
    with mathopt:Solver:Solve use sum, foreach, +, -, *, ≼, ≽, /, =
    // Bring the data into scope.
    with data use facilities, customers, distances, demands
 
    // Define the model variables.
    def variables = rel:mathopt:variables[
            {:y, "binary", (facilities)};
            {:x, "continuous", (facilities, customers)};
            {:x, "lower_bound", 0.0};
            {:x, "upper_bound", 1.0}
        ]
    // Bring the model variables into scope.
    with variables use x, y
 
    // Define the sense of the optimization problem and the objective function.
    def minimize = sum[distances[f, c] * demands[c] * x[f, c] for f in facilities, c in customers]
 
    // Define the facility selection constraint.
    // This ensures that only one facility location is selected.
    def constraints:facility_selection = sum[y[f] for f in facilities] = 1
 
    // Define the demand constraints.
    // This ensures that each customer's demand is fully satisfied.
    def constraints:demand = foreach[c in customers : sum[x[f, c] for f in facilities] = 1]
 
    // Define the linking constraints.
    // This ensures that if a facility is not selected, it cannot serve any customers.
    def constraints:linking = foreach[f in facilities, c in customers : x[f, c] ≼ y[f]]
end
 
// Optimize the model.
def result = rel:mathopt:optimize[model, {}]
 
// Extract the result of the optimization process.
def output = rel:mathopt:extract_result[result, model:variables]

References

Operations Research An Introduction 10th Edition, Hamdy A. Taha.
Numerical Optimization 2nd Edition, Jorge Nocedal and Stephen J. Wright.
Facility Location with JuMP (opens in a new tab)
Facility Location Problem (opens in a new tab)

Was this doc helpful?