From 39e525fb24c9c8f1e2b1cf746f13eb6c8b014048 Mon Sep 17 00:00:00 2001 From: Daniel Espinoza Date: Thu, 26 Mar 2020 15:18:46 -0400 Subject: [PATCH] Add uncapacitated facility location example --- examples/cpp/BUILD | 10 + .../cpp/uncapacitated_facility_location.cc | 230 ++++++++++++++++++ makefiles/Makefile.cpp.mk | 1 + 3 files changed, 241 insertions(+) create mode 100644 examples/cpp/uncapacitated_facility_location.cc diff --git a/examples/cpp/BUILD b/examples/cpp/BUILD index fd55fb7bb4..82c1956a3a 100644 --- a/examples/cpp/BUILD +++ b/examples/cpp/BUILD @@ -489,6 +489,16 @@ cc_binary( ], ) +cc_binary( + name = "uncapacitated_facility_location", + srcs = ["uncapacitated_facility_location.cc"], + deps = [ + "//ortools/base", + "//ortools/linear_solver", + ], +) + + cc_binary( name = "variable_intervals_sat", srcs = ["variable_intervals_sat.cc"], diff --git a/examples/cpp/uncapacitated_facility_location.cc b/examples/cpp/uncapacitated_facility_location.cc new file mode 100644 index 0000000000..27e6622733 --- /dev/null +++ b/examples/cpp/uncapacitated_facility_location.cc @@ -0,0 +1,230 @@ +// Copyright 2020 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// +// Uncapacitated Facility Location Problem. +// A description of the problem can be found here: +// https://en.wikipedia.org/wiki/Facility_location_problem. +// The variant which is tackled by this model does not consider capacities +// for facilities. Moreover, all cost are based on euclidean distance factors, +// i.e. the problem we really solve is a Metric Facility Location. For the +// sake of simplicity, facilities and demands are randomly located. Distances +// are assumed to be in meters and times in seconds. + +#include +#include + +#include "google/protobuf/text_format.h" +#include "ortools/base/commandlineflags.h" +#include "ortools/base/random.h" +#include "ortools/base/integral_types.h" +#include "ortools/base/logging.h" +#include "ortools/linear_solver/linear_solver.h" + +DEFINE_int32(verbose, 0, "Verbosity level."); +DEFINE_int32(facilities, 20, "Candidate facilities to consider."); +DEFINE_int32(clients, 100, "Clients to serve."); +DEFINE_double(fix_cost, 5000, "Cost of opening a facility."); + +namespace operations_research { + +typedef struct { + double x{0}; + double y{0}; + } Location; + +typedef struct { + int f{-1}; + int c{-1}; + MPVariable* x{nullptr}; + } Edge; + +static double Distance(const Location& src, const Location& dst) { + return sqrt((src.x-dst.x)*(src.x-dst.x) + (src.y-dst.y)*(src.y-dst.y)); +} + +static void UncapacitatedFacilityLocation(int32 facilities, + int32 clients, double fix_cost, + MPSolver::OptimizationProblemType optimization_problem_type) { + LOG(INFO) << "Starting " << __func__; + // Local Constants + const int32 kXMax = 1000; + const int32 kYMax = 1000; + const double kMaxDistance = 6*sqrt((kXMax*kYMax))/facilities; + int kStrLen{1024}; + // char buffer for names + char name_buffer[kStrLen+1]; + name_buffer[kStrLen] = '\0'; + LOG(INFO) << "Facilities/Clients/Fix cost/MaxDist: " << facilities << "/" + << clients << "/" << fix_cost << "/" << kMaxDistance; + // Setting up facilities and demand points + MTRandom randomizer(/*fixed seed*/20191029); + std::vector facility(facilities); + std::vector client(clients); + for (int i = 0; i < facilities; ++i) { + facility[i].x = randomizer.Uniform(kXMax + 1); + facility[i].y = randomizer.Uniform(kYMax + 1); + } + for (int i = 0; i < clients; ++i) { + client[i].x = randomizer.Uniform(kXMax + 1); + client[i].y = randomizer.Uniform(kYMax + 1); + } + + // Setup uncapacitated facility location model: + // Min sum( c_f * x_f : f in Facilities) + sum(x_{f,c} * x_{f,c} : {f,c} in E) + // s.t. (1) sum(x_{f,c} : f in Facilities) >= 1 forall c in Clients + // (2) x_f - x_{f,c} >= 0 forall {f,c} in E + // (3) x_f in {0,1} forall f in Facilities + // + // We consider E as the pairs {f,c} in Facilities x Clients such that + // Distance(f,c) <= kMaxDistance + MPSolver solver("UncapacitatedFacilityLocation", optimization_problem_type); + const double infinity = solver.infinity(); + MPObjective* objective = solver.MutableObjective(); + objective->SetMinimization(); + + // Add binary facilities variables + std::vector xf{}; + for (int f = 0; f < facilities; ++f) { + snprintf(name_buffer, kStrLen, "x[%d](%g,%g)", f, + facility[f].x, facility[f].y); + MPVariable* x = solver.MakeBoolVar(name_buffer); + xf.push_back(x); + objective->SetCoefficient(x, fix_cost); + } + + // Build edge variables + std::vector edges; + for (int c = 0; c < clients; ++c) { + snprintf(name_buffer, kStrLen, "R-Client[%d](%g,%g)", c, + client[c].x, client[c].y); + MPConstraint* client_constraint = solver.MakeRowConstraint(/* lb */1, + /* ub */infinity, name_buffer); + for (int f = 0; f < facilities; ++f) { + double distance = Distance(facility[f], client[c]); + if (distance > kMaxDistance) continue; + Edge edge{}; + snprintf(name_buffer, kStrLen, "x[%d,%d]", f, c); + edge.x = solver.MakeNumVar(/* lb */0, /*ub */1, name_buffer); + edge.f = f; + edge.c = c; + edges.push_back(edge); + objective->SetCoefficient(edge.x, distance); + // coefficient for constraint (1) + client_constraint->SetCoefficient(edge.x, 1); + // add constraint (2) + snprintf(name_buffer, kStrLen, "R-Edge[%d,%d]", f, c); + MPConstraint* edge_constraint = solver.MakeRowConstraint(/* lb */0, + /* ub */infinity, name_buffer); + edge_constraint->SetCoefficient(edge.x, -1); + edge_constraint->SetCoefficient(xf[f], 1); + } + }// End adding all edge variables + LOG(INFO) << "Number of variables = " << solver.NumVariables(); + LOG(INFO) << "Number of constraints = " << solver.NumConstraints(); + // display on screen LP if small enough + if (clients <= 10 && facilities <= 10) { + std::string lp_string{}; + solver.ExportModelAsLpFormat(/* obfuscate */false, &lp_string); + std::cout << "LP-Model:\n" << lp_string << std::endl; + } + // Set options and solve + solver.SetNumThreads(8); + solver.EnableOutput(); + const MPSolver::ResultStatus result_status = solver.Solve(); + // Check that the problem has an optimal solution. + if (result_status != MPSolver::OPTIMAL) { + LOG(FATAL) << "The problem does not have an optimal solution!"; + } else { + LOG(INFO) << "Optimal objective value = " << objective->Value(); + if (FLAGS_verbose) { + std::vector> solution(facilities); + for (auto& edge : edges) { + if (edge.x->solution_value() < 0.5) continue; + solution[edge.f].push_back(edge.c); + } + std::cout << "\tSolution:\n"; + for (int f = 0; f < facilities; ++f) { + if (solution[f].size() < 1) continue; + assert(xf[f]->solution_value() > 0.5); + snprintf(name_buffer, kStrLen, "\t Facility[%d](%g,%g):", f, + facility[f].x, facility[f].y); + std::cout << name_buffer; + int i = 1; + for (auto c : solution[f]) { + snprintf(name_buffer, kStrLen, " Client[%d](%g,%g)", c, + client[c].x, client[c].y); + if(i++ >= 5) { + std::cout << "\n\t\t"; + i = 1; + } + std::cout << name_buffer; + } + std::cout << "\n"; + } + } + std::cout << "\n"; + LOG(INFO) << ""; + LOG(INFO) << "Advanced usage:"; + LOG(INFO) << "Problem solved in " << solver.DurationSinceConstruction() + << " milliseconds"; + LOG(INFO) << "Problem solved in " << solver.iterations() << " iterations"; + LOG(INFO) << "Problem solved in " << solver.nodes() + << " branch-and-bound nodes"; + } +} + +void RunAllExamples(int32 facilities, int32 clients, double fix_cost) { +#if defined(USE_CBC) + LOG(INFO) << "---- Integer programming example with CBC ----"; + UncapacitatedFacilityLocation(facilities, clients, fix_cost, + MPSolver::CBC_MIXED_INTEGER_PROGRAMMING); +#endif +#if defined(USE_GLPK) + LOG(INFO) << "---- Integer programming example with GLPK ----"; + UncapacitatedFacilityLocation(facilities, clients, fix_cost, + MPSolver::GLPK_MIXED_INTEGER_PROGRAMMING); +#endif +#if defined(USE_SCIP) + LOG(INFO) << "---- Integer programming example with SCIP ----"; + UncapacitatedFacilityLocation(facilities, clients, fix_cost, + MPSolver::SCIP_MIXED_INTEGER_PROGRAMMING); +#endif +#if defined(USE_GUROBI) + LOG(INFO) << "---- Integer programming example with Gurobi ----"; + UncapacitatedFacilityLocation(facilities, clients, fix_cost, + MPSolver::GUROBI_MIXED_INTEGER_PROGRAMMING); +#endif // USE_GUROBI +#if defined(USE_CPLEX) + LOG(INFO) << "---- Integer programming example with CPLEX ----"; + UncapacitatedFacilityLocation(facilities, clients, fix_cost, + MPSolver::CPLEX_MIXED_INTEGER_PROGRAMMING); +#endif // USE_CPLEX +} + +} // namespace operations_research + +int main(int argc, char** argv) { + google::InitGoogleLogging(argv[0]); + gflags::SetUsageMessage( + std::string("This program solve a (randomly generated)\n") + + std::string("Uncapacitated Facility Location Problem. Sample Usage:\n")); + gflags::ParseCommandLineFlags(&argc, &argv, true); + CHECK_LT(0, FLAGS_facilities) << "Specify an instance size greater than 0."; + CHECK_LT(0, FLAGS_clients) << "Specify a non-null client size."; + CHECK_LT(0, FLAGS_fix_cost) << "Specify a non-null client size."; + FLAGS_logtostderr = 1; + operations_research::RunAllExamples(FLAGS_facilities, FLAGS_clients, + FLAGS_fix_cost); + return EXIT_SUCCESS; +} diff --git a/makefiles/Makefile.cpp.mk b/makefiles/Makefile.cpp.mk index decc5d4c3e..ed5dd7a13d 100755 --- a/makefiles/Makefile.cpp.mk +++ b/makefiles/Makefile.cpp.mk @@ -429,6 +429,7 @@ test_cc_cpp: \ rcc_random_tsp \ rcc_slitherlink_sat \ rcc_strawberry_fields_with_column_generation \ + rcc_uncapacitated_facility_location \ rcc_weighted_tardiness_sat $(MAKE) run \ SOURCE=examples/cpp/dimacs_assignment.cc \