1 // Copyright 2020 Google LLC
2 // Licensed under the Apache License, Version 2.0 (the "License");
3 // you may not use this file except in compliance with the License.
4 // You may obtain a copy of the License at
5 //
6 // http://www.apache.org/licenses/LICENSE-2.0
7 //
8 // Unless required by applicable law or agreed to in writing, software
9 // distributed under the License is distributed on an "AS IS" BASIS,
10 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11 // See the License for the specific language governing permissions and
12 // limitations under the License.
13
14 //
15 // Uncapacitated Facility Location Problem.
16 // A description of the problem can be found here:
17 // https://en.wikipedia.org/wiki/Facility_location_problem.
18 // The variant which is tackled by this model does not consider capacities
19 // for facilities. Moreover, all cost are based on euclidean distance factors,
20 // i.e. the problem we really solve is a Metric Facility Location. For the
21 // sake of simplicity, facilities and demands are randomly located. Distances
22 // are assumed to be in meters and times in seconds.
23
24 #include <cstdio>
25 #include <vector>
26
27 #include "absl/random/random.h"
28 #include "google/protobuf/text_format.h"
29 #include "ortools/base/commandlineflags.h"
30 #include "ortools/base/integral_types.h"
31 #include "ortools/base/logging.h"
32 #include "ortools/linear_solver/linear_solver.h"
33 #include "ortools/util/random_engine.h"
34
35 ABSL_FLAG(int, verbose, 0, "Verbosity level.");
36 ABSL_FLAG(int, facilities, 20, "Candidate facilities to consider.");
37 ABSL_FLAG(int, clients, 100, "Clients to serve.");
38 ABSL_FLAG(double, fix_cost, 5000, "Cost of opening a facility.");
39
40 namespace operations_research {
41
42 struct Location {
43 double x = 0.0;
44 double y = 0.0;
45 };
46
47 struct Edge {
48 int f = -1;
49 int c = -1;
50 MPVariable* x = nullptr;
51 };
52
Distance(const Location & src,const Location & dst)53 static double Distance(const Location& src, const Location& dst) {
54 return sqrt((src.x - dst.x) * (src.x - dst.x) +
55 (src.y - dst.y) * (src.y - dst.y));
56 }
57
UncapacitatedFacilityLocation(int32_t facilities,int32_t clients,double fix_cost,MPSolver::OptimizationProblemType optimization_problem_type)58 static void UncapacitatedFacilityLocation(
59 int32_t facilities, int32_t clients, double fix_cost,
60 MPSolver::OptimizationProblemType optimization_problem_type) {
61 LOG(INFO) << "Starting " << __func__;
62 // Local Constants
63 const int32_t kXMax = 1000;
64 const int32_t kYMax = 1000;
65 const double kMaxDistance = 6 * sqrt((kXMax * kYMax)) / facilities;
66 const int kStrLen = 1024;
67 // char buffer for names
68 char name_buffer[kStrLen + 1];
69 name_buffer[kStrLen] = '\0';
70 LOG(INFO) << "Facilities/Clients/Fix cost/MaxDist: " << facilities << "/"
71 << clients << "/" << fix_cost << "/" << kMaxDistance;
72 // Setting up facilities and demand points
73 random_engine_t randomizer; // Deterministic random generator.
74 std::vector<Location> facility(facilities);
75 std::vector<Location> client(clients);
76 for (int i = 0; i < facilities; ++i) {
77 facility[i].x = absl::Uniform(randomizer, 0, kXMax + 1);
78 facility[i].y = absl::Uniform(randomizer, 0, kYMax + 1);
79 }
80 for (int i = 0; i < clients; ++i) {
81 client[i].x = absl::Uniform(randomizer, 0, kXMax + 1);
82 client[i].y = absl::Uniform(randomizer, 0, kYMax + 1);
83 }
84
85 // Setup uncapacitated facility location model:
86 // Min sum( c_f * x_f : f in Facilities) + sum(x_{f,c} * x_{f,c} : {f,c} in E)
87 // s.t. (1) sum(x_{f,c} : f in Facilities) >= 1 forall c in Clients
88 // (2) x_f - x_{f,c} >= 0 forall {f,c} in E
89 // (3) x_f in {0,1} forall f in Facilities
90 //
91 // We consider E as the pairs {f,c} in Facilities x Clients such that
92 // Distance(f,c) <= kMaxDistance
93 MPSolver solver("UncapacitatedFacilityLocation", optimization_problem_type);
94 const double infinity = solver.infinity();
95 MPObjective* objective = solver.MutableObjective();
96 objective->SetMinimization();
97
98 // Add binary facilities variables
99 std::vector<MPVariable*> xf{};
100 for (int f = 0; f < facilities; ++f) {
101 snprintf(name_buffer, kStrLen, "x[%d](%g,%g)", f, facility[f].x,
102 facility[f].y);
103 MPVariable* x = solver.MakeBoolVar(name_buffer);
104 xf.push_back(x);
105 objective->SetCoefficient(x, fix_cost);
106 }
107
108 // Build edge variables
109 std::vector<Edge> edges;
110 for (int c = 0; c < clients; ++c) {
111 snprintf(name_buffer, kStrLen, "R-Client[%d](%g,%g)", c, client[c].x,
112 client[c].y);
113 MPConstraint* client_constraint =
114 solver.MakeRowConstraint(/* lb */ 1, /* ub */ infinity, name_buffer);
115 for (int f = 0; f < facilities; ++f) {
116 double distance = Distance(facility[f], client[c]);
117 if (distance > kMaxDistance) continue;
118 Edge edge{};
119 snprintf(name_buffer, kStrLen, "x[%d,%d]", f, c);
120 edge.x = solver.MakeNumVar(/* lb */ 0, /*ub */ 1, name_buffer);
121 edge.f = f;
122 edge.c = c;
123 edges.push_back(edge);
124 objective->SetCoefficient(edge.x, distance);
125 // coefficient for constraint (1)
126 client_constraint->SetCoefficient(edge.x, 1);
127 // add constraint (2)
128 snprintf(name_buffer, kStrLen, "R-Edge[%d,%d]", f, c);
129 MPConstraint* edge_constraint =
130 solver.MakeRowConstraint(/* lb */ 0, /* ub */ infinity, name_buffer);
131 edge_constraint->SetCoefficient(edge.x, -1);
132 edge_constraint->SetCoefficient(xf[f], 1);
133 }
134 } // End adding all edge variables
135 LOG(INFO) << "Number of variables = " << solver.NumVariables();
136 LOG(INFO) << "Number of constraints = " << solver.NumConstraints();
137 // display on screen LP if small enough
138 if (clients <= 10 && facilities <= 10) {
139 std::string lp_string{};
140 solver.ExportModelAsLpFormat(/* obfuscate */ false, &lp_string);
141 std::cout << "LP-Model:\n" << lp_string << std::endl;
142 }
143 // Set options and solve
144 if (optimization_problem_type != MPSolver::SCIP_MIXED_INTEGER_PROGRAMMING) {
145 if (!solver.SetNumThreads(8).ok()) {
146 LOG(INFO) << "Could not set parallelism for " << optimization_problem_type;
147 }
148 }
149 solver.EnableOutput();
150 const MPSolver::ResultStatus result_status = solver.Solve();
151 // Check that the problem has an optimal solution.
152 if (result_status != MPSolver::OPTIMAL) {
153 LOG(FATAL) << "The problem does not have an optimal solution!";
154 } else {
155 LOG(INFO) << "Optimal objective value = " << objective->Value();
156 if (absl::GetFlag(FLAGS_verbose)) {
157 std::vector<std::vector<int> > solution(facilities);
158 for (auto& edge : edges) {
159 if (edge.x->solution_value() < 0.5) continue;
160 solution[edge.f].push_back(edge.c);
161 }
162 std::cout << "\tSolution:\n";
163 for (int f = 0; f < facilities; ++f) {
164 if (solution[f].size() < 1) continue;
165 assert(xf[f]->solution_value() > 0.5);
166 snprintf(name_buffer, kStrLen, "\t Facility[%d](%g,%g):", f,
167 facility[f].x, facility[f].y);
168 std::cout << name_buffer;
169 int i = 1;
170 for (auto c : solution[f]) {
171 snprintf(name_buffer, kStrLen, " Client[%d](%g,%g)", c, client[c].x,
172 client[c].y);
173 if (i++ >= 5) {
174 std::cout << "\n\t\t";
175 i = 1;
176 }
177 std::cout << name_buffer;
178 }
179 std::cout << "\n";
180 }
181 }
182 std::cout << "\n";
183 LOG(INFO) << "";
184 LOG(INFO) << "Advanced usage:";
185 LOG(INFO) << "Problem solved in " << solver.DurationSinceConstruction()
186 << " milliseconds";
187 LOG(INFO) << "Problem solved in " << solver.iterations() << " iterations";
188 LOG(INFO) << "Problem solved in " << solver.nodes()
189 << " branch-and-bound nodes";
190 }
191 }
192
RunAllExamples(int32_t facilities,int32_t clients,double fix_cost)193 void RunAllExamples(int32_t facilities, int32_t clients, double fix_cost) {
194 #if defined(USE_CBC)
195 LOG(INFO) << "---- Integer programming example with CBC ----";
196 UncapacitatedFacilityLocation(facilities, clients, fix_cost,
197 MPSolver::CBC_MIXED_INTEGER_PROGRAMMING);
198 #endif
199 #if defined(USE_GLPK)
200 LOG(INFO) << "---- Integer programming example with GLPK ----";
201 UncapacitatedFacilityLocation(facilities, clients, fix_cost,
202 MPSolver::GLPK_MIXED_INTEGER_PROGRAMMING);
203 #endif
204 #if defined(USE_SCIP)
205 LOG(INFO) << "---- Integer programming example with SCIP ----";
206 UncapacitatedFacilityLocation(facilities, clients, fix_cost,
207 MPSolver::SCIP_MIXED_INTEGER_PROGRAMMING);
208 #endif
209 #if defined(USE_GUROBI)
210 LOG(INFO) << "---- Integer programming example with Gurobi ----";
211 UncapacitatedFacilityLocation(facilities, clients, fix_cost,
212 MPSolver::GUROBI_MIXED_INTEGER_PROGRAMMING);
213 #endif // USE_GUROBI
214 #if defined(USE_CPLEX)
215 LOG(INFO) << "---- Integer programming example with CPLEX ----";
216 UncapacitatedFacilityLocation(facilities, clients, fix_cost,
217 MPSolver::CPLEX_MIXED_INTEGER_PROGRAMMING);
218 #endif // USE_CPLEX
219 LOG(INFO) << "---- Integer programming example with CP-SAT ----";
220 UncapacitatedFacilityLocation(facilities, clients, fix_cost,
221 MPSolver::SAT_INTEGER_PROGRAMMING);
222 }
223
224 } // namespace operations_research
225
main(int argc,char ** argv)226 int main(int argc, char** argv) {
227 google::InitGoogleLogging(argv[0]);
228 absl::SetProgramUsageMessage(
229 std::string("This program solve a (randomly generated)\n") +
230 std::string("Uncapacitated Facility Location Problem. Sample Usage:\n"));
231 absl::ParseCommandLine(argc, argv);
232 CHECK_LT(0, absl::GetFlag(FLAGS_facilities))
233 << "Specify an instance size greater than 0.";
234 CHECK_LT(0, absl::GetFlag(FLAGS_clients))
235 << "Specify a non-null client size.";
236 CHECK_LT(0, absl::GetFlag(FLAGS_fix_cost))
237 << "Specify a non-null client size.";
238 absl::SetFlag(&FLAGS_logtostderr, 1);
239 operations_research::RunAllExamples(absl::GetFlag(FLAGS_facilities),
240 absl::GetFlag(FLAGS_clients),
241 absl::GetFlag(FLAGS_fix_cost));
242 return EXIT_SUCCESS;
243 }
244