From a29d1eed56a91439a548e75b8a4921d51107cd77 Mon Sep 17 00:00:00 2001 From: Corentin Le Molgat Date: Thu, 13 Jul 2023 20:22:58 +0200 Subject: [PATCH] Export math_opt from google3 --- ortools/math_opt/BUILD.bazel | 1 + ortools/math_opt/cpp/BUILD.bazel | 4 +- ortools/math_opt/cpp/solution.h | 1 - ortools/math_opt/cpp/solve_result.cc | 42 +- ortools/math_opt/cpp/solve_result.h | 71 +- ortools/math_opt/result.proto | 49 +- ortools/math_opt/samples/BUILD.bazel | 48 +- ...g.cc => advanced_linear_programming_mo.cc} | 5 +- ortools/math_opt/samples/area_socp_mo.cc | 88 +++ ortools/math_opt/samples/basic_example.cc | 72 -- ortools/math_opt/samples/basic_example_mo.cc | 17 +- ortools/math_opt/samples/cocktail_hour.cc | 380 ---------- ortools/math_opt/samples/cocktail_hour_mo.cc | 14 +- ortools/math_opt/samples/cutting_stock.cc | 277 ------- ortools/math_opt/samples/cutting_stock_mo.cc | 124 ++-- .../math_opt/samples/facility_lp_benders.cc | 683 ------------------ .../samples/facility_lp_benders_mo.cc | 44 +- ...graph_coloring.cc => graph_coloring_mo.cc} | 9 +- .../math_opt/samples/integer_programming.cc | 85 --- .../samples/integer_programming_mo.cc | 25 +- .../math_opt/samples/lagrangian_relaxation.cc | 493 ------------- .../samples/lagrangian_relaxation_mo.cc | 75 +- .../math_opt/samples/linear_programming.cc | 93 --- .../math_opt/samples/linear_programming_mo.cc | 8 +- ..._regression.cc => linear_regression_mo.cc} | 6 +- ...uling.cc => time_indexed_scheduling_mo.cc} | 0 ortools/math_opt/samples/tsp.cc | 359 --------- ortools/math_opt/samples/tsp_mo.cc | 48 +- ortools/math_opt/solvers/BUILD.bazel | 1 + 29 files changed, 415 insertions(+), 2707 deletions(-) rename ortools/math_opt/samples/{advanced_linear_programming.cc => advanced_linear_programming_mo.cc} (95%) create mode 100644 ortools/math_opt/samples/area_socp_mo.cc delete mode 100644 ortools/math_opt/samples/basic_example.cc delete mode 100644 ortools/math_opt/samples/cocktail_hour.cc delete mode 100644 ortools/math_opt/samples/cutting_stock.cc delete mode 100644 ortools/math_opt/samples/facility_lp_benders.cc rename ortools/math_opt/samples/{graph_coloring.cc => graph_coloring_mo.cc} (96%) delete mode 100644 ortools/math_opt/samples/integer_programming.cc delete mode 100644 ortools/math_opt/samples/lagrangian_relaxation.cc delete mode 100644 ortools/math_opt/samples/linear_programming.cc rename ortools/math_opt/samples/{linear_regression.cc => linear_regression_mo.cc} (97%) rename ortools/math_opt/samples/{time_indexed_scheduling.cc => time_indexed_scheduling_mo.cc} (100%) delete mode 100644 ortools/math_opt/samples/tsp.cc diff --git a/ortools/math_opt/BUILD.bazel b/ortools/math_opt/BUILD.bazel index c8fb770d769..3191d3e6936 100644 --- a/ortools/math_opt/BUILD.bazel +++ b/ortools/math_opt/BUILD.bazel @@ -110,6 +110,7 @@ proto_library( deps = [ ":solution_proto", "//ortools/gscip:gscip_proto", + "//ortools/pdlp:solve_log_proto", "@com_google_protobuf//:duration_proto", ], ) diff --git a/ortools/math_opt/cpp/BUILD.bazel b/ortools/math_opt/cpp/BUILD.bazel index 13daabea9f2..b495aeb1377 100644 --- a/ortools/math_opt/cpp/BUILD.bazel +++ b/ortools/math_opt/cpp/BUILD.bazel @@ -170,7 +170,6 @@ cc_library( ":linear_constraint", ":solution", ":variable_and_expressions", - "//ortools/base", "//ortools/base:protoutil", "//ortools/base:status_macros", "//ortools/gscip:gscip_cc_proto", @@ -180,7 +179,6 @@ cc_library( "//ortools/port:proto_utils", "//ortools/util:fp_roundtrip_conv", "//ortools/util:status_macros", - "@com_google_absl//absl/container:flat_hash_map", "@com_google_absl//absl/status:statusor", "@com_google_absl//absl/strings", "@com_google_absl//absl/time", @@ -368,6 +366,7 @@ cc_library( "//ortools/math_opt:parameters_cc_proto", "//ortools/math_opt/solvers:gurobi_cc_proto", "@com_google_absl//absl/status:statusor", + "@com_google_absl//absl/strings", ], ) @@ -385,6 +384,7 @@ cc_library( "//ortools/math_opt:parameters_cc_proto", "//ortools/math_opt/solvers:glpk_cc_proto", "//ortools/math_opt/solvers:gurobi_cc_proto", + "//ortools/math_opt/solvers:highs_cc_proto", "//ortools/port:proto_utils", "//ortools/sat:sat_parameters_cc_proto", "//ortools/util:status_macros", diff --git a/ortools/math_opt/cpp/solution.h b/ortools/math_opt/cpp/solution.h index 18984685733..a3ac7e8f127 100644 --- a/ortools/math_opt/cpp/solution.h +++ b/ortools/math_opt/cpp/solution.h @@ -21,7 +21,6 @@ #include "absl/status/status.h" #include "absl/status/statusor.h" -#include "absl/types/span.h" #include "ortools/math_opt/cpp/basis_status.h" #include "ortools/math_opt/cpp/enums.h" // IWYU pragma: export #include "ortools/math_opt/cpp/linear_constraint.h" diff --git a/ortools/math_opt/cpp/solve_result.cc b/ortools/math_opt/cpp/solve_result.cc index 777a734fac1..48ff99ccf3e 100644 --- a/ortools/math_opt/cpp/solve_result.cc +++ b/ortools/math_opt/cpp/solve_result.cc @@ -20,13 +20,12 @@ #include #include -#include "absl/container/flat_hash_map.h" -#include "absl/status/status.h" #include "absl/status/statusor.h" +#include "absl/strings/str_cat.h" +#include "absl/strings/str_join.h" #include "absl/strings/string_view.h" #include "absl/time/time.h" #include "absl/types/span.h" -#include "ortools/base/logging.h" #include "ortools/base/protoutil.h" #include "ortools/base/status_macros.h" #include "ortools/math_opt/cpp/linear_constraint.h" @@ -170,6 +169,37 @@ bool Termination::limit_reached() const { reason == TerminationReason::kNoSolutionFound; } +absl::Status Termination::ReasonIs(const TerminationReason reason) const { + if (this->reason == reason) return absl::OkStatus(); + return util::InternalErrorBuilder() + << "expected termination reason '" << reason << "' but got " << *this; +} + +absl::Status Termination::ReasonIsAnyOf( + std::initializer_list reasons) const { + for (const TerminationReason reason : reasons) { + if (this->reason == reason) return absl::OkStatus(); + } + return util::InternalErrorBuilder() + << "expected termination reason in {" + << absl::StrJoin(reasons, ", ", + [](std::string* out, TerminationReason reason) { + absl::StrAppend(out, "'"); + absl::StreamFormatter()(out, reason); + absl::StrAppend(out, "'"); + }) + << "} but got " << *this; +} + +absl::Status Termination::IsOptimal() const { + return ReasonIs(TerminationReason::kOptimal); +} + +absl::Status Termination::IsOptimalOrFeasible() const { + return ReasonIsAnyOf( + {TerminationReason::kOptimal, TerminationReason::kFeasible}); +} + absl::StatusOr Termination::FromProto( const TerminationProto& termination_proto) { const std::optional reason = @@ -327,6 +357,9 @@ absl::StatusOr SolveResult::Proto() const { if (gscip_solver_specific_output.ByteSizeLong() > 0) { *result.mutable_gscip_output() = gscip_solver_specific_output; } + if (pdlp_solver_specific_output.ByteSizeLong() > 0) { + *result.mutable_pdlp_output() = pdlp_solver_specific_output; + } return result; } @@ -365,6 +398,9 @@ absl::StatusOr SolveResult::FromProto( case SolveResultProto::kGscipOutput: result.gscip_solver_specific_output = solve_result_proto.gscip_output(); return result; + case SolveResultProto::kPdlpOutput: + result.pdlp_solver_specific_output = solve_result_proto.pdlp_output(); + return result; case SolveResultProto::SOLVER_SPECIFIC_OUTPUT_NOT_SET: return result; } diff --git a/ortools/math_opt/cpp/solve_result.h b/ortools/math_opt/cpp/solve_result.h index cf3ba9b7e30..9222f8a175f 100644 --- a/ortools/math_opt/cpp/solve_result.h +++ b/ortools/math_opt/cpp/solve_result.h @@ -17,6 +17,7 @@ #ifndef OR_TOOLS_MATH_OPT_CPP_SOLVE_RESULT_H_ #define OR_TOOLS_MATH_OPT_CPP_SOLVE_RESULT_H_ +#include #include #include #include @@ -25,8 +26,6 @@ #include "absl/status/statusor.h" #include "absl/time/time.h" -#include "absl/types/span.h" -#include "ortools/base/logging.h" #include "ortools/gscip/gscip.pb.h" #include "ortools/math_opt/cpp/enums.h" // IWYU pragma: export #include "ortools/math_opt/cpp/linear_constraint.h" @@ -103,23 +102,26 @@ struct SolveStats { // contract is clarified. // Solver claims the optimal value is equal or better (smaller for - // minimization and larger for maximization) than best_primal_bound: + // minimization and larger for maximization) than best_primal_bound up to the + // solvers primal feasibility tolerance (see warning below): // * best_primal_bound is trivial (+inf for minimization and -inf - // maximization) when the solver does not claim to have such bound. This - // may happen for some solvers (e.g., PDLP, typically continuous solvers) - // even when returning optimal (solver could terminate with slightly - // infeasible primal solutions). + // maximization) when the solver does not claim to have such bound. // * best_primal_bound can be closer to the optimal value than the objective // of the best primal feasible solution. In particular, best_primal_bound // may be non-trivial even when no primal feasible solutions are returned. - // * best_dual_bound is always better (smaller for minimization and larger - // for maximization) than best_primal_bound. + // Warning: The precise claim is that there exists a primal solution that: + // * is numerically feasible (i.e. feasible up to the solvers tolerance), and + // * has an objective value best_primal_bound. + // This numerically feasible solution could be slightly infeasible, in which + // case best_primal_bound could be strictly better than the optimal value. + // Translating a primal feasibility tolerance to a tolerance on + // best_primal_bound is non-trivial, specially when the feasibility tolerance + // is relatively large (e.g. when solving with PDLP). double best_primal_bound = 0.0; // Solver claims the optimal value is equal or worse (larger for - // minimization and smaller for maximization) than best_dual_bound: - // * best_dual_bound is always better (smaller for minimization and larger - // for maximization) than best_primal_bound. + // minimization and smaller for maximization) than best_dual_bound up to the + // solvers dual feasibility tolerance (see warning below): // * best_dual_bound is trivial (-inf for minimization and +inf // maximization) when the solver does not claim to have such bound. // Similarly to best_primal_bound, this may happen for some solvers even @@ -129,6 +131,28 @@ struct SolveStats { // value than the objective of the best dual feasible solution. For MIP // one of the first non-trivial values for best_dual_bound is often the // optimal value of the LP relaxation of the MIP. + // * best_dual_bound should be better (smaller for minimization and larger + // for maximization) than best_primal_bound up to the solvers tolerances + // (see warning below). + // Warning: + // * For continuous problems, the precise claim is that there exists a + // dual solution that: + // * is numerically feasible (i.e. feasible up to the solvers tolerance), + // and + // * has an objective value best_dual_bound. + // This numerically feasible solution could be slightly infeasible, in + // which case best_dual_bound could be strictly worse than the optimal + // value and best_primal_bound. Similar to the primal case, translating a + // dual feasibility tolerance to a tolerance on best_dual_bound is + // non-trivial, specially when the feasibility tolerance is relatively + // large. However, some solvers provide a corrected version of + // best_dual_bound that can be numerically safer. This corrected version + // can be accessed through the solver's specific output (e.g. for PDLP, + // pdlp_output.convergence_information.corrected_dual_objective). + // * For MIP solvers, best_dual_bound may be associated to a dual solution + // for some continuous relaxation (e.g. LP relaxation), but it is often a + // complex consequence of the solvers execution and is typically more + // imprecise than the bounds reported by LP solvers. double best_dual_bound = 0.0; // Feasibility statuses for primal and dual problems. @@ -285,6 +309,27 @@ struct Termination { // kNoSolutionFound, and limit is not empty). bool limit_reached() const; + // Returns an OkStatus if the reason of this `Termination` is + // `TerminationReason::kOptimal` or `TerminationReason::kFeasible`, or an + // `InternalError` otherwise. + absl::Status IsOptimalOrFeasible() const; + + // Returns an OkStatus if the reason of this `Termination` is + // `TerminationReason::kOptimal`, or an `InternalError` otherwise. + // + // In most use cases, at least for MIPs, `IsOptimalOrFeasible` should be used + // instead. + absl::Status IsOptimal() const; + + // Returns an OkStatus if the reason of this `Termination` is `reason`, or an + // `InternalError` otherwise. + absl::Status ReasonIs(TerminationReason reason) const; + + // Returns an OkStatus if the reason of this `Termination` is in `reasons`, or + // an `InternalError` otherwise. + absl::Status ReasonIsAnyOf( + std::initializer_list reasons) const; + // Sets the reason to kFeasible static Termination Feasible(Limit limit, std::string detail = {}); @@ -340,6 +385,8 @@ struct SolveResult { // Solver specific output from Gscip. Only populated if Gscip is used. GScipOutput gscip_solver_specific_output; + // Solver specific output from Pdlp. Only populated if Pdlp is used. + SolveResultProto::PdlpOutput pdlp_solver_specific_output; // Returns the SolveResult equivalent of solve_result_proto. // diff --git a/ortools/math_opt/result.proto b/ortools/math_opt/result.proto index 7b9643a5c9a..f8f7e2744c1 100644 --- a/ortools/math_opt/result.proto +++ b/ortools/math_opt/result.proto @@ -17,6 +17,7 @@ syntax = "proto3"; package operations_research.math_opt; import "google/protobuf/duration.proto"; +import "ortools/pdlp/solve_log.proto"; import "ortools/gscip/gscip.proto"; import "ortools/math_opt/solution.proto"; @@ -79,21 +80,26 @@ message SolveStatsProto { google.protobuf.Duration solve_time = 1; // Solver claims the optimal value is equal or better (smaller for - // minimization and larger for maximization) than best_primal_bound: + // minimization and larger for maximization) than best_primal_bound up to the + // solvers primal feasibility tolerance (see warning below): // * best_primal_bound is trivial (+inf for minimization and -inf - // maximization) when the solver does not claim to have such bound. This - // may happen for some solvers (e.g. PDLP, typically continuous solvers) - // even when returning optimal (solver could terminate with slightly - // infeasible primal solutions). + // maximization) when the solver does not claim to have such bound. // * best_primal_bound can be closer to the optimal value than the objective // of the best primal feasible solution. In particular, best_primal_bound // may be non-trivial even when no primal feasible solutions are returned. + // Warning: The precise claim is that there exists a primal solution that: + // * is numerically feasible (i.e. feasible up to the solvers tolerance), and + // * has an objective value best_primal_bound. + // This numerically feasible solution could be slightly infeasible, in which + // case best_primal_bound could be strictly better than the optimal value. + // Translating a primal feasibility tolerance to a tolerance on + // best_primal_bound is non-trivial, specially when the feasibility tolerance + // is relatively large (e.g. when solving with PDLP). double best_primal_bound = 2; // Solver claims the optimal value is equal or worse (larger for - // minimization and smaller for maximization) than best_dual_bound: - // * best_dual_bound is always better (smaller for minimization and larger - // for maximization) than best_primal_bound. + // minimization and smaller for maximization) than best_dual_bound up to the + // solvers dual feasibility tolerance (see warning below): // * best_dual_bound is trivial (-inf for minimization and +inf // maximization) when the solver does not claim to have such bound. // Similarly to best_primal_bound, this may happen for some solvers even @@ -103,6 +109,28 @@ message SolveStatsProto { // value than the objective of the best dual feasible solution. For MIP // one of the first non-trivial values for best_dual_bound is often the // optimal value of the LP relaxation of the MIP. + // * best_dual_bound should be better (smaller for minimization and larger + // for maximization) than best_primal_bound up to the solvers tolerances + // (see warning below). + // Warning: + // * For continuous problems, the precise claim is that there exists a + // dual solution that: + // * is numerically feasible (i.e. feasible up to the solvers tolerance), + // and + // * has an objective value best_dual_bound. + // This numerically feasible solution could be slightly infeasible, in + // which case best_dual_bound could be strictly worse than the optimal + // value and best_primal_bound. Similar to the primal case, translating a + // dual feasibility tolerance to a tolerance on best_dual_bound is + // non-trivial, specially when the feasibility tolerance is relatively + // large. However, some solvers provide a corrected version of + // best_dual_bound that can be numerically safer. This corrected version + // can be accessed through the solver's specific output (e.g. for PDLP, + // pdlp_output.convergence_information.corrected_dual_objective). + // * For MIP solvers, best_dual_bound may be associated to a dual solution + // for some continuous relaxation (e.g. LP relaxation), but it is often a + // complex consequence of the solvers execution and is typically more + // imprecise than the bounds reported by LP solvers. double best_dual_bound = 3; // Feasibility statuses for primal and dual problems. @@ -284,8 +312,13 @@ message SolveResultProto { // Statistics on the solve process, e.g. running time, iterations. SolveStatsProto solve_stats = 6; + message PdlpOutput { + pdlp.ConvergenceInformation convergence_information = 1; + } + oneof solver_specific_output { GScipOutput gscip_output = 7; + PdlpOutput pdlp_output = 9; } reserved 1; // Deleted fields. diff --git a/ortools/math_opt/samples/BUILD.bazel b/ortools/math_opt/samples/BUILD.bazel index 9cef4dc98c3..3dcf21eddf2 100644 --- a/ortools/math_opt/samples/BUILD.bazel +++ b/ortools/math_opt/samples/BUILD.bazel @@ -14,8 +14,8 @@ package(default_visibility = ["//visibility:public"]) cc_binary( - name = "basic_example", - srcs = ["basic_example.cc"], + name = "basic_example_mo", + srcs = ["basic_example_mo.cc"], deps = [ "//ortools/base", "//ortools/base:status_macros", @@ -27,8 +27,8 @@ cc_binary( ) cc_binary( - name = "cocktail_hour", - srcs = ["cocktail_hour.cc"], + name = "cocktail_hour_mo", + srcs = ["cocktail_hour_mo.cc"], deps = [ "//ortools/base", "//ortools/base:map_util", @@ -44,8 +44,8 @@ cc_binary( ) cc_binary( - name = "linear_programming", - srcs = ["linear_programming.cc"], + name = "linear_programming_mo", + srcs = ["linear_programming_mo.cc"], deps = [ "//ortools/base", "//ortools/base:status_macros", @@ -58,8 +58,8 @@ cc_binary( ) cc_binary( - name = "integer_programming", - srcs = ["integer_programming.cc"], + name = "integer_programming_mo", + srcs = ["integer_programming_mo.cc"], deps = [ "//ortools/base", "//ortools/base:status_macros", @@ -71,8 +71,8 @@ cc_binary( ) cc_binary( - name = "cutting_stock", - srcs = ["cutting_stock.cc"], + name = "cutting_stock_mo", + srcs = ["cutting_stock_mo.cc"], deps = [ "//ortools/base", "//ortools/base:status_macros", @@ -85,8 +85,8 @@ cc_binary( ) cc_binary( - name = "facility_lp_benders", - srcs = ["facility_lp_benders.cc"], + name = "facility_lp_benders_mo", + srcs = ["facility_lp_benders_mo.cc"], deps = [ "//ortools/base", "//ortools/base:status_macros", @@ -109,8 +109,8 @@ cc_binary( ) cc_binary( - name = "lagrangian_relaxation", - srcs = ["lagrangian_relaxation.cc"], + name = "lagrangian_relaxation_mo", + srcs = ["lagrangian_relaxation_mo.cc"], deps = [ "//ortools/base", "//ortools/base:container_logging", @@ -129,8 +129,8 @@ cc_binary( ) cc_binary( - name = "tsp", - srcs = ["tsp.cc"], + name = "tsp_mo", + srcs = ["tsp_mo.cc"], deps = [ "//ortools/base", "//ortools/base:status_macros", @@ -145,8 +145,8 @@ cc_binary( ) cc_binary( - name = "linear_regression", - srcs = ["linear_regression.cc"], + name = "linear_regression_mo", + srcs = ["linear_regression_mo.cc"], deps = [ "//ortools/base", "//ortools/base:status_macros", @@ -159,8 +159,8 @@ cc_binary( ) cc_binary( - name = "advanced_linear_programming", - srcs = ["advanced_linear_programming.cc"], + name = "advanced_linear_programming_mo", + srcs = ["advanced_linear_programming_mo.cc"], deps = [ "//ortools/base", "//ortools/base:status_macros", @@ -173,8 +173,8 @@ cc_binary( ) cc_binary( - name = "time_indexed_scheduling", - srcs = ["time_indexed_scheduling.cc"], + name = "time_indexed_scheduling_mo", + srcs = ["time_indexed_scheduling_mo.cc"], deps = [ "//ortools/base", "//ortools/base:status_macros", @@ -190,8 +190,8 @@ cc_binary( ) cc_binary( - name = "graph_coloring", - srcs = ["graph_coloring.cc"], + name = "graph_coloring_mo", + srcs = ["graph_coloring_mo.cc"], deps = [ "//ortools/base", "//ortools/base:status_macros", diff --git a/ortools/math_opt/samples/advanced_linear_programming.cc b/ortools/math_opt/samples/advanced_linear_programming_mo.cc similarity index 95% rename from ortools/math_opt/samples/advanced_linear_programming.cc rename to ortools/math_opt/samples/advanced_linear_programming_mo.cc index d28706bc9cd..bf01e644cee 100644 --- a/ortools/math_opt/samples/advanced_linear_programming.cc +++ b/ortools/math_opt/samples/advanced_linear_programming_mo.cc @@ -67,10 +67,7 @@ absl::Status Main() { ASSIGN_OR_RETURN(const math_opt::SolveResult result, Solve(model, math_opt::SolverType::kGlop)); - if (result.termination.reason != math_opt::TerminationReason::kOptimal) { - return util::InternalErrorBuilder() - << "model failed to solve to optimality" << result.termination; - } + RETURN_IF_ERROR(result.termination.IsOptimal()); std::cout << "Problem solved in " << result.solve_time() << std::endl; std::cout << "Objective value: " << result.objective_value() << std::endl; diff --git a/ortools/math_opt/samples/area_socp_mo.cc b/ortools/math_opt/samples/area_socp_mo.cc new file mode 100644 index 00000000000..4c127ea0473 --- /dev/null +++ b/ortools/math_opt/samples/area_socp_mo.cc @@ -0,0 +1,88 @@ +// Copyright 2010-2022 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. + +// Simple SOCP problem showing that minimizing the perimeter of a rectangle +// with fixed area results in equal width and height. +#include +#include +#include +#include + +#include "absl/status/status.h" +#include "ortools/base/init_google.h" +#include "ortools/base/logging.h" +#include "ortools/base/status_builder.h" +#include "ortools/base/status_macros.h" +#include "ortools/math_opt/cpp/math_opt.h" + +ABSL_FLAG(double, area, 9, "Area lower bound."); + +namespace { + +namespace math_opt = ::operations_research::math_opt; + +constexpr double kInf = std::numeric_limits::infinity(); + +// We want to minimize the width plus height of a rectangle with area A. +// +// First we can relax to the area being at least A: +// min width + height +// s.t. width*height >= A (Area) +// width in [0.0, infinity) +// height in [0.0, infinity) +// +// Next we need to reformulate the area constraint as a second order cone +// constraint: +// min width + height +// s.t. ||((width - height)/2, sqrt(A))||_2 <= (width + height)/2 (Area-SOCP) +// width in [0.0, infinity) +// height in [0.0, infinity) +// +// To see how these two problems are equivalent, first note that by squaring +// both sides of constraint (Area-SOCP) we can see that it is equivalent to: +// (width - height)^2/4 + A <= (width + height)^2/4 +// because width + height >= 0. Expanding the two squares and reordering shows +// the equivalence to constraint (Area). +absl::Status Main(const double target_area) { + math_opt::Model model; + const math_opt::Variable width = + model.AddContinuousVariable(0.0, kInf, "width"); + const math_opt::Variable height = + model.AddContinuousVariable(0.0, kInf, "height"); + + model.AddSecondOrderConeConstraint( + {(width - height) / 2, std::sqrt(target_area)}, (width + height) / 2); + model.Minimize(width + height); + ASSIGN_OR_RETURN(const math_opt::SolveResult result, + Solve(model, math_opt::SolverType::kEcos)); + RETURN_IF_ERROR(result.termination.IsOptimalOrFeasible()); + std::cout << "Target area: " << target_area << std::endl; + std::cout << "Area: " + << result.variable_values().at(width) * + result.variable_values().at(height) + << std::endl; + std::cout << "Perimeter = " << result.objective_value() << std::endl; + std::cout << "Width: " << result.variable_values().at(width) << std::endl; + std::cout << "Height: " << result.variable_values().at(height) << std::endl; + return absl::OkStatus(); +} +} // namespace + +int main(int argc, char** argv) { + InitGoogle(argv[0], &argc, &argv, true); + const absl::Status status = Main(absl::GetFlag(FLAGS_area)); + if (!status.ok()) { + LOG(QFATAL) << status; + } + return 0; +} diff --git a/ortools/math_opt/samples/basic_example.cc b/ortools/math_opt/samples/basic_example.cc deleted file mode 100644 index c88de3dd068..00000000000 --- a/ortools/math_opt/samples/basic_example.cc +++ /dev/null @@ -1,72 +0,0 @@ -// Copyright 2010-2022 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. - -// Testing correctness of the code snippets in the comments of math_opt.h. - -#include -#include -#include - -#include "absl/status/status.h" -#include "absl/status/statusor.h" -#include "ortools/base/init_google.h" -#include "ortools/base/logging.h" -#include "ortools/base/status_builder.h" -#include "ortools/base/status_macros.h" -#include "ortools/math_opt/cpp/math_opt.h" - -namespace { - -namespace math_opt = ::operations_research::math_opt; - -// Model the problem: -// max 2.0 * x + y -// s.t. x + y <= 1.5 -// x in {0.0, 1.0} -// y in [0.0, 2.5] -// -absl::Status Main() { - math_opt::Model model("my_model"); - const math_opt::Variable x = model.AddBinaryVariable("x"); - const math_opt::Variable y = model.AddContinuousVariable(0.0, 2.5, "y"); - // We can directly use linear combinations of variables ... - model.AddLinearConstraint(x + y <= 1.5, "c"); - // ... or build them incrementally. - math_opt::LinearExpression objective_expression; - objective_expression += 2 * x; - objective_expression += y; - model.Maximize(objective_expression); - ASSIGN_OR_RETURN(const math_opt::SolveResult result, - Solve(model, math_opt::SolverType::kGscip)); - switch (result.termination.reason) { - case math_opt::TerminationReason::kOptimal: - case math_opt::TerminationReason::kFeasible: - std::cout << "Objective value: " << result.objective_value() << std::endl - << "Value for variable x: " << result.variable_values().at(x) - << std::endl; - return absl::OkStatus(); - default: - return util::InternalErrorBuilder() - << "model failed to solve: " << result.termination; - } -} -} // namespace - -int main(int argc, char** argv) { - InitGoogle(argv[0], &argc, &argv, true); - const absl::Status status = Main(); - if (!status.ok()) { - LOG(QFATAL) << status; - } - return 0; -} diff --git a/ortools/math_opt/samples/basic_example_mo.cc b/ortools/math_opt/samples/basic_example_mo.cc index 3fbe0e4459b..e4ed504748f 100644 --- a/ortools/math_opt/samples/basic_example_mo.cc +++ b/ortools/math_opt/samples/basic_example_mo.cc @@ -15,6 +15,7 @@ #include #include +#include #include "absl/status/status.h" #include "absl/status/statusor.h" @@ -47,17 +48,11 @@ absl::Status Main() { model.Maximize(objective_expression); ASSIGN_OR_RETURN(const math_opt::SolveResult result, Solve(model, math_opt::SolverType::kGscip)); - switch (result.termination.reason) { - case math_opt::TerminationReason::kOptimal: - case math_opt::TerminationReason::kFeasible: - std::cout << "Objective value: " << result.objective_value() << std::endl - << "Value for variable x: " << result.variable_values().at(x) - << std::endl; - return absl::OkStatus(); - default: - return util::InternalErrorBuilder() - << "model failed to solve: " << result.termination; - } + RETURN_IF_ERROR(result.termination.IsOptimalOrFeasible()); + std::cout << "Objective value: " << result.objective_value() << std::endl + << "Value for variable x: " << result.variable_values().at(x) + << std::endl; + return absl::OkStatus(); } } // namespace diff --git a/ortools/math_opt/samples/cocktail_hour.cc b/ortools/math_opt/samples/cocktail_hour.cc deleted file mode 100644 index de07eee34de..00000000000 --- a/ortools/math_opt/samples/cocktail_hour.cc +++ /dev/null @@ -1,380 +0,0 @@ -// Copyright 2010-2022 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. - -// Pick ingredients to buy to make the maximum number of cocktails. -// -// Given a list of cocktails, each of which is made from a list of ingredients, -// and a budget of how many ingredients you can buy, solve a MIP to pick a -// subset of the ingredients so that you can make the largest number of -// cocktails. -// -// This program can be run in three modes: -// text: Outputs the optimal set of ingredients and cocktails that can be -// produced as plain text to standard out. -// latex: Outputs a menu of the cocktails that can be made as LaTeX code to -// standard out. -// analysis: Computes the number of cocktails that can be made as a function -// of the number of ingredients for all values. -// -// In latex mode, the output can be piped directly to pdflatex, e.g. -// blaze run -c opt \ -// ortools/math_opt/examples/cpp/cocktail_hour \ -// -- --num_ingredients 10 --mode latex | pdflatex -output-directory /tmp -// will create a PDF in /tmp. -#include -#include -#include -#include -#include - -#include "absl/container/flat_hash_set.h" -#include "absl/status/status.h" -#include "absl/status/statusor.h" -#include "absl/strings/str_join.h" -#include "absl/strings/str_replace.h" -#include "absl/strings/string_view.h" -#include "ortools/base/init_google.h" -#include "ortools/base/logging.h" -#include "ortools/base/map_util.h" -#include "ortools/base/status_macros.h" -#include "ortools/math_opt/cpp/math_opt.h" -#include "ortools/util/status_macros.h" - -ABSL_FLAG(std::string, mode, "text", - "One of \"text\", \"latex\", or \"analysis\"."); -ABSL_FLAG(int, num_ingredients, 10, - "How many ingredients to buy (ignored in analysis mode)."); -ABSL_FLAG(std::vector, existing_ingredients, {}, - "Ingredients you already have (ignored in analysis mode)."); -ABSL_FLAG(std::vector, unavailable_ingredients, {}, - "Ingredients you cannot get (ignored in analysis mode)."); -ABSL_FLAG(std::vector, required_cocktails, {}, - "Cocktails you must be able to make (ignored in analysis mode)."); -ABSL_FLAG(std::vector, blocked_cocktails, {}, - "Cocktails to exclude from the menu (ignored in analysis mode)."); - -namespace { - -namespace math_opt = ::operations_research::math_opt; - -constexpr absl::string_view kIngredients[] = {"Amaro Nonino", - "All Spice Dram", - "Aperol", - "Bitters", - "Bourbon", - "Brandy", - "Campari", - "Cinnamon", - "Chambord", - "Cherry", - "Cloves", - "Cointreau", - "Coke", - "Cranberry", - "Creme de Cacao", - "Creme de Violette", - "Cucumber", - "Egg", - "Gin", - "Green Chartreuse", - "Heavy Cream", - "Lemon", - "Lillet Blanc", - "Lime", - "Luxardo", - "Mint", - "Orange", - "Orange Flower Water Extract", - "Orgeat", - "Pickle", - "Pineapple Juice", - "Pisco", - "Prosecco", - "Raspberry Vodka", - "Ruby Port", - "Rum", - "Seltzer", - "Simple Syrup", - "Sugar", - "Sweet Vermouth", - "Tequila", - "Tonic Water", - "Vodka"}; - -constexpr std::size_t kIngredientsSize = - sizeof(kIngredients) / sizeof(kIngredients[0]); - -struct Cocktail { - std::string name; - std::vector ingredients; -}; - -std::vector AllCocktails() { - return { - // Aperitifs - {.name = "Prosecco glass", .ingredients = {"Prosecco"}}, - {.name = "Aperol Spritz", .ingredients = {"Prosecco", "Aperol"}}, - {.name = "Chambord Spritz", .ingredients = {"Prosecco", "Chambord"}}, - {.name = "Improved French 75", - .ingredients = {"Prosecco", "Vodka", "Lemon", "Simple Syrup"}}, - // Quick and Simple - {.name = "Gin and Tonic", .ingredients = {"Gin", "Tonic Water", "Lime"}}, - {.name = "Rum and Coke", .ingredients = {"Rum", "Coke"}}, - {.name = "Improved Manhattan", - .ingredients = {"Bourbon", "Sweet Vermouth", "Bitters"}}, - // Vodka - - // Serve with a sugared rim - {.name = "Lemon Drop", - .ingredients = {"Vodka", "Cointreau", "Lemon", "Simple Syrup"}}, - // Shake, then float 2oz Prosecco after pouring - {.name = "Big Crush", - .ingredients = {"Raspberry Vodka", "Cointreau", "Lemon", "Chambord", - "Prosecco"}}, - {.name = "Cosmopolitan", - .ingredients = {"Vodka", "Cranberry", "Cointreau", "Lime"}}, - // A shot, chase with 1/3 of pickle spear - {.name = "Vodka/Pickle", .ingredients = {"Vodka", "Pickle"}}, - - // Gin - {.name = "Last Word", - .ingredients = {"Gin", "Green Chartreuse", "Luxardo", "Lime"}}, - {.name = "Corpse Reviver #2 (Lite)", - .ingredients = {"Gin", "Cointreau", "Lillet Blanc", "Lemon"}}, - {.name = "Negroni", .ingredients = {"Gin", "Sweet Vermouth", "Campari"}}, - // "Float" Creme de Violette (it will sink) - {.name = "Aviation", - .ingredients = {"Gin", "Luxardo", "Lemon", "Creme de Violette"}}, - - // Bourbon - {.name = "Paper Plane", - .ingredients = {"Bourbon", "Aperol", "Amaro Nonino", "Lemon"}}, - {.name = "Derby", - .ingredients = {"Bourbon", "Sweet Vermouth", "Lime", "Cointreau"}}, - // Muddle sugar, water, bitters, and orange peel. Garnish with a Luxardo - // cherry (do not cheap out), spill cherry syrup generously in drink - {.name = "Old Fashioned", - .ingredients = {"Bourbon", "Sugar", "Bitters", "Orange", "Cherry"}}, - {.name = "Boulevardier", - .ingredients = {"Bourbon", "Sweet Vermouth", "Campari"}}, - - // Tequila - {.name = "Margarita", .ingredients = {"Tequila", "Cointreau", "Lime"}}, - // Shake with chopped cucumber and strain. Garnish with cucumber. - {.name = "Midnight Cruiser", - .ingredients = {"Tequila", "Aperol", "Lime", "Pineapple Juice", - "Cucumber", "Simple Syrup"}}, - - {.name = "Tequila shot", .ingredients = {"Tequila"}}, - // Rum - - // Shake with light rum, float a dark rum on top. - {.name = "Pineapple Mai Tai", - .ingredients = {"Rum", "Lime", "Orgeat", "Cointreau", - "Pineapple Juice"}}, - {.name = "Daiquiri", .ingredients = {"Rum", "Lime", "Simple Syrup"}}, - {.name = "Mojito", - .ingredients = {"Rum", "Lime", "Simple Syrup", "Mint", "Seltzer"}}, - // Add bitters generously. Invert half lime to form a cup, fill with - // Green Chartreuse and cloves. Float lime cup on drink and ignite. - {.name = "Kennedy", - .ingredients = {"Rum", "All Spice Dram", "Bitters", "Lime", - "Simple Syrup", "Cloves", "Green Chartreuse"}}, - - // Egg - - {.name = "Pisco Sour", - .ingredients = {"Pisco", "Lime", "Simple Syrup", "Egg", "Bitters"}}, - {.name = "Viana", - .ingredients = {"Ruby Port", "Brandy", "Creme de Cacao", "Sugar", "Egg", - "Cinnamon"}}, - // Add cream last before shaking (and seltzer after shaking). Shake for 10 - // minutes, no less. - {.name = "Ramos gin fizz", - .ingredients = {"Gin", "Seltzer", "Heavy Cream", - "Orange Flower Water Extract", "Egg", "Lemon", "Lime", - "Simple Syrup"}}}; -} - -struct Menu { - std::vector ingredients; - std::vector cocktails; -}; - -absl::StatusOr SolveForMenu( - const int max_new_ingredients, const bool enable_solver_output, - const absl::flat_hash_set& existing_ingredients, - const absl::flat_hash_set& unavailable_ingredients, - const absl::flat_hash_set& required_cocktails, - const absl::flat_hash_set& blocked_cocktails) { - const std::vector all_cocktails = AllCocktails(); - math_opt::Model model("Cocktail hour"); - absl::flat_hash_map ingredient_vars; - for (const absl::string_view ingredient : kIngredients) { - const double lb = existing_ingredients.contains(ingredient) ? 1.0 : 0.0; - const double ub = unavailable_ingredients.contains(ingredient) ? 0.0 : 1.0; - const math_opt::Variable v = model.AddIntegerVariable(lb, ub, ingredient); - gtl::InsertOrDie(&ingredient_vars, std::string(ingredient), v); - } - math_opt::LinearExpression ingredients_used; - for (const auto& [name, ingredient_var] : ingredient_vars) { - ingredients_used += ingredient_var; - } - model.AddLinearConstraint(ingredients_used <= - max_new_ingredients + existing_ingredients.size()); - - absl::flat_hash_map cocktail_vars; - for (const Cocktail& cocktail : all_cocktails) { - const double lb = required_cocktails.contains(cocktail.name) ? 1.0 : 0.0; - const double ub = blocked_cocktails.contains(cocktail.name) ? 0.0 : 1.0; - const math_opt::Variable v = - model.AddIntegerVariable(lb, ub, cocktail.name); - for (const std::string& ingredient : cocktail.ingredients) { - model.AddLinearConstraint(v <= - gtl::FindOrDie(ingredient_vars, ingredient)); - } - gtl::InsertOrDie(&cocktail_vars, cocktail.name, v); - } - math_opt::LinearExpression cocktails_made; - for (const auto& [name, cocktail_var] : cocktail_vars) { - cocktails_made += cocktail_var; - } - model.Maximize(cocktails_made); - const math_opt::SolveArguments args = { - .parameters = {.enable_output = enable_solver_output}}; - ASSIGN_OR_RETURN(const math_opt::SolveResult result, - math_opt::Solve(model, math_opt::SolverType::kGscip, args)); - - switch (result.termination.reason) { - case math_opt::TerminationReason::kOptimal: - case math_opt::TerminationReason::kFeasible: - break; - default: - return util::InternalErrorBuilder() - << "failed to find a solution: " << result.termination; - } - Menu menu; - for (const absl::string_view ingredient : kIngredients) { - if (result.variable_values().at(ingredient_vars.at(ingredient)) > 0.5) { - menu.ingredients.push_back(std::string(ingredient)); - } - } - for (const Cocktail& cocktail : all_cocktails) { - if (result.variable_values().at(cocktail_vars.at(cocktail.name)) > 0.5) { - menu.cocktails.push_back(cocktail); - } - } - return menu; -} - -absl::flat_hash_set SetFromVec( - const std::vector& vec) { - return {vec.begin(), vec.end()}; -} - -absl::Status AnalysisMode() { - std::cout << "Considering " << AllCocktails().size() << " cocktails and " - << kIngredientsSize << " ingredients." << std::endl; - std::cout << "Solving for number of cocktails that can be made as a function " - "of number of ingredients" - << std::endl; - - std::cout << "ingredients | cocktails" << std::endl; - for (int i = 1; i <= kIngredientsSize; ++i) { - const absl::StatusOr menu = SolveForMenu( - i, false, /*existing_ingredients=*/{}, /*unavailable_ingredients=*/{}, - /*required_cocktails=*/{}, /*blocked_cocktails=*/{}); - RETURN_IF_ERROR(menu.status()) - << "Failure when solving for " << i << " ingredients"; - std::cout << i << " | " << menu->cocktails.size() << std::endl; - } - return absl::OkStatus(); -} - -std::string ExportToLaTeX(const std::vector& cocktails, - absl::string_view title = "Cocktail Hour") { - std::vector lines; - lines.push_back("\\documentclass{article}"); - lines.push_back("\\usepackage{fullpage}"); - lines.push_back("\\linespread{2}"); - lines.push_back("\\begin{document}"); - lines.push_back("\\begin{center}"); - lines.push_back(absl::StrCat("\\begin{Huge}", title, "\\end{Huge}")); - lines.push_back(""); - for (const Cocktail& cocktail : cocktails) { - lines.push_back(absl::StrCat(cocktail.name, "---{\\em ", - absl::StrJoin(cocktail.ingredients, ", "), - "}")); - lines.push_back(""); - } - lines.push_back("\\end{center}"); - lines.push_back("\\end{document}"); - - return absl::StrReplaceAll(absl::StrJoin(lines, "\n"), {{"#", "\\#"}}); -} - -absl::Status Main() { - const std::string mode = absl::GetFlag(FLAGS_mode); - CHECK(absl::flat_hash_set({"text", "latex", "analysis"}) - .contains(mode)) - << "Unexpected mode: " << mode; - - // We are in analysis mode. - if (mode == "analysis") { - return AnalysisMode(); - } - - OR_ASSIGN_OR_RETURN3( - Menu menu, - SolveForMenu(absl::GetFlag(FLAGS_num_ingredients), mode == "text", - SetFromVec(absl::GetFlag(FLAGS_existing_ingredients)), - SetFromVec(absl::GetFlag(FLAGS_unavailable_ingredients)), - SetFromVec(absl::GetFlag(FLAGS_required_cocktails)), - SetFromVec(absl::GetFlag(FLAGS_blocked_cocktails))), - _ << "error when solving for optimal set of ingredients"); - - // We are in latex mode. - if (mode == "latex") { - std::cout << ExportToLaTeX(menu.cocktails) << std::endl; - return absl::OkStatus(); - } - - // We are in text mode - std::cout << "Considered " << AllCocktails().size() << " cocktails and " - << kIngredientsSize << " ingredients." << std::endl; - std::cout << "Solution has " << menu.ingredients.size() - << " ingredients to make " << menu.cocktails.size() << " cocktails." - << std::endl - << std::endl; - - std::cout << "Ingredients:" << std::endl; - for (const std::string& ingredient : menu.ingredients) { - std::cout << " " << ingredient << std::endl; - } - std::cout << "Cocktails:" << std::endl; - for (const Cocktail& cocktail : menu.cocktails) { - std::cout << " " << cocktail.name << std::endl; - } - return absl::OkStatus(); -} - -} // namespace - -int main(int argc, char** argv) { - InitGoogle(argv[0], &argc, &argv, true); - const absl::Status status = Main(); - if (!status.ok()) { - LOG(QFATAL) << status; - } - return 0; -} diff --git a/ortools/math_opt/samples/cocktail_hour_mo.cc b/ortools/math_opt/samples/cocktail_hour_mo.cc index 3742dbc6a80..3017f4ccf11 100644 --- a/ortools/math_opt/samples/cocktail_hour_mo.cc +++ b/ortools/math_opt/samples/cocktail_hour_mo.cc @@ -33,6 +33,9 @@ // will create a PDF in /tmp. #include #include +#include +#include +#include #include "absl/container/flat_hash_set.h" #include "absl/status/status.h" @@ -252,14 +255,7 @@ absl::StatusOr SolveForMenu( ASSIGN_OR_RETURN(const math_opt::SolveResult result, math_opt::Solve(model, math_opt::SolverType::kGscip, args)); - switch (result.termination.reason) { - case math_opt::TerminationReason::kOptimal: - case math_opt::TerminationReason::kFeasible: - break; - default: - return util::InternalErrorBuilder() - << "failed to find a solution: " << result.termination; - } + RETURN_IF_ERROR(result.termination.IsOptimalOrFeasible()); Menu menu; for (const absl::string_view ingredient : kIngredients) { if (result.variable_values().at(ingredient_vars.at(ingredient)) > 0.5) { @@ -299,7 +295,7 @@ absl::Status AnalysisMode() { } std::string ExportToLaTeX(const std::vector& cocktails, - const std::string& title = "Cocktail Hour") { + absl::string_view title = "Cocktail Hour") { std::vector lines; lines.push_back("\\documentclass{article}"); lines.push_back("\\usepackage{fullpage}"); diff --git a/ortools/math_opt/samples/cutting_stock.cc b/ortools/math_opt/samples/cutting_stock.cc deleted file mode 100644 index f947e8104b3..00000000000 --- a/ortools/math_opt/samples/cutting_stock.cc +++ /dev/null @@ -1,277 +0,0 @@ -// Copyright 2010-2022 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. - -// The Cutting Stock problem is as follows. You begin with unlimited boards, all -// of the same length. You are also given a list of smaller pieces to cut out, -// each with a length and a demanded quantity. You want to cut out all these -// pieces using as few of your starting boards as possible. -// -// E.g. you begin with boards that are 20 feet long, and you must cut out 3 -// pieces that are 6 feet long and 5 pieces that are 8 feet long. An optimal -// solution is: -// [(6,), (8, 8) (8, 8), (6, 6, 8)] -// (We cut a 6 foot piece from the first board, two 8 foot pieces from -// the second board, and so on.) -// -// This example approximately solves the problem with a column generation -// heuristic. The leader problem is a set cover problem, and the worker is a -// knapsack problem. We alternate between solving the LP relaxation of the -// leader incrementally, and solving the worker to generate new a configuration -// (a column) for the leader. When the worker can no longer find a column -// improving the LP cost, we convert the leader problem to a MIP and solve -// again. We now give precise statements of the leader and worker. -// -// Problem data: -// * l_i: the length of each piece we need to cut out. -// * d_i: how many copies each piece we need. -// * L: the length of our initial boards. -// * q_ci: for configuration c, the quantity of piece i produced. -// -// Leader problem variables: -// * x_c: how many copies of configuration c to produce. -// -// Leader problem formulation: -// min sum_c x_c -// s.t. sum_c q_ci * x_c = d_i for all i -// x_c >= 0, integer for all c. -// -// The worker problem is to generate new configurations for the leader problem -// based on the dual variables of the demand constraints in the LP relaxation. -// Worker problem data: -// * p_i: The "price" of piece i (dual value from leader's demand constraint) -// -// Worker decision variables: -// * y_i: How many copies of piece i should be in the configuration. -// -// Worker formulation -// max sum_i p_i * y_i -// s.t. sum_i l_i * y_i <= L -// y_i >= 0, integer for all i -// -// An optimal solution y* defines a new configuration c with q_ci = y_i* for all -// i. If the solution has objective value <= 1, no further improvement on the LP -// is possible. For additional background and proofs see: -// https://people.orie.cornell.edu/shmoys/or630/notes-06/lec16.pdf -// or any other reference on the "Cutting Stock Problem". -// -// Note: this problem is equivalent to symmetric bin packing: -// https://en.wikipedia.org/wiki/Bin_packing_problem#Formal_statement -// but typically in bin packing it is not assumed that you should exploit having -// multiple items of the same size. -#include -#include -#include -#include -#include -#include - -#include "absl/status/status.h" -#include "absl/status/statusor.h" -#include "ortools/base/init_google.h" -#include "ortools/base/logging.h" -#include "ortools/base/status_builder.h" -#include "ortools/base/status_macros.h" -#include "ortools/math_opt/cpp/math_opt.h" - -namespace { - -namespace math_opt = operations_research::math_opt; -constexpr double kInf = std::numeric_limits::infinity(); - -// piece_sizes and piece_demands must have equal length. -// every piece must have 0 < size <= board_length. -// every piece must have demand > 0. -struct CuttingStockInstance { - std::vector piece_sizes; - std::vector piece_demands; - int board_length; -}; - -// pieces and quantity must have equal size. -// Defined for a related CuttingStockInstance, the total length all pieces -// weighted by their quantity must not exceed board_length. -struct Configuration { - std::vector pieces; - std::vector quantity; -}; - -// configurations and quantity must have equal size. -// objective_value is the sum of the vales in quantity (how many total boards -// are used). -// To be feasible, the demand for each piece type must be met by the produced -// configurations. -struct CuttingStockSolution { - std::vector configurations; - std::vector quantity; - int objective_value = 0; -}; - -// Solves the worker problem. -// -// Solves the problem on finding the configuration (with its objective value) to -// add the to model that will give the greatest improvement in the LP -// relaxation. This is equivalent to a knapsack problem. -absl::StatusOr> BestConfiguration( - const std::vector& piece_prices, - const std::vector& piece_sizes, const int board_size) { - int num_pieces = piece_prices.size(); - CHECK_EQ(piece_sizes.size(), num_pieces); - math_opt::Model model("knapsack"); - std::vector pieces; - for (int i = 0; i < num_pieces; ++i) { - pieces.push_back( - model.AddIntegerVariable(0, kInf, absl::StrCat("item_", i))); - } - model.Maximize(math_opt::InnerProduct(pieces, piece_prices)); - model.AddLinearConstraint(math_opt::InnerProduct(pieces, piece_sizes) <= - board_size); - ASSIGN_OR_RETURN(const math_opt::SolveResult solve_result, - math_opt::Solve(model, math_opt::SolverType::kCpSat)); - if (solve_result.termination.reason != - math_opt::TerminationReason::kOptimal) { - return util::InternalErrorBuilder() - << "Failed to solve knapsack pricing problem to optimality: " - << solve_result.termination; - } - Configuration config; - for (int i = 0; i < num_pieces; ++i) { - const int use = static_cast( - std::round(solve_result.variable_values().at(pieces[i]))); - if (use > 0) { - config.pieces.push_back(i); - config.quantity.push_back(use); - } - } - return std::make_pair(config, solve_result.objective_value()); -} - -// Solves the full cutting stock problem by decomposition. -absl::StatusOr SolveCuttingStock( - const CuttingStockInstance& instance) { - math_opt::Model model("cutting_stock"); - model.set_minimize(); - const int n = instance.piece_sizes.size(); - std::vector demand_met; - for (int i = 0; i < n; ++i) { - const int d = instance.piece_demands[i]; - demand_met.push_back(model.AddLinearConstraint(d, d)); - } - std::vector> configs; - auto add_config = [&](const Configuration& config) { - const math_opt::Variable v = model.AddContinuousVariable(0.0, kInf); - model.set_objective_coefficient(v, 1); - for (int i = 0; i < config.pieces.size(); ++i) { - const int item = config.pieces[i]; - const int use = config.quantity[i]; - if (use >= 1) { - model.set_coefficient(demand_met[item], v, use); - } - } - configs.push_back({config, v}); - }; - - // To ensure the leader problem is always feasible, begin a configuration for - // every item that has a single copy of the item. - for (int i = 0; i < n; ++i) { - add_config(Configuration{.pieces = {i}, .quantity = {1}}); - } - - ASSIGN_OR_RETURN(auto solver, math_opt::IncrementalSolver::New( - &model, math_opt::SolverType::kGlop)); - int pricing_round = 0; - while (true) { - ASSIGN_OR_RETURN(math_opt::SolveResult solve_result, solver->Solve()); - if (solve_result.termination.reason != - math_opt::TerminationReason::kOptimal) { - return util::InternalErrorBuilder() - << "Failed to solve leader LP problem to optimality at " - "iteration " - << pricing_round << " termination: " << solve_result.termination; - } - if (!solve_result.has_dual_feasible_solution()) { - // MathOpt does not require solvers to return a dual solution on optimal, - // but most LP solvers always will, see go/mathopt-solver-contracts for - // details. - return util::InternalErrorBuilder() - << "no dual solution was returned with optimal solution at " - "iteration " - << pricing_round; - } - std::vector prices; - for (const math_opt::LinearConstraint d : demand_met) { - prices.push_back(solve_result.dual_values().at(d)); - } - ASSIGN_OR_RETURN( - (const auto [config, value]), - BestConfiguration(prices, instance.piece_sizes, instance.board_length)); - if (value <= 1 + 1e-3) { - // The LP relaxation is solved, we can stop adding columns. - break; - } - add_config(config); - LOG(INFO) << "round: " << pricing_round - << " lp objective: " << solve_result.objective_value(); - pricing_round++; - } - LOG(INFO) << "Done adding columns, switching to MIP"; - for (const auto& [config, var] : configs) { - model.set_integer(var); - } - ASSIGN_OR_RETURN(const math_opt::SolveResult solve_result, - math_opt::Solve(model, math_opt::SolverType::kCpSat)); - switch (solve_result.termination.reason) { - case math_opt::TerminationReason::kOptimal: - case math_opt::TerminationReason::kFeasible: - break; - default: - return util::InternalErrorBuilder() - << "Failed to solve final cutting stock MIP, termination: " - << solve_result.termination; - } - CuttingStockSolution solution; - for (const auto& [config, var] : configs) { - int use = - static_cast(std::round(solve_result.variable_values().at(var))); - if (use > 0) { - solution.configurations.push_back(config); - solution.quantity.push_back(use); - solution.objective_value += use; - } - } - return solution; -} - -absl::Status RealMain() { - // Data from https://en.wikipedia.org/wiki/Cutting_stock_problem - CuttingStockInstance instance; - instance.board_length = 5600; - instance.piece_sizes = {1380, 1520, 1560, 1710, 1820, 1880, 1930, - 2000, 2050, 2100, 2140, 2150, 2200}; - instance.piece_demands = {22, 25, 12, 14, 18, 18, 20, 10, 12, 14, 16, 18, 20}; - ASSIGN_OR_RETURN(CuttingStockSolution solution, SolveCuttingStock(instance)); - std::cout << "Best known solution uses 73 rolls." << std::endl; - std::cout << "Total rolls used in actual solution found: " - << solution.objective_value << std::endl; - return absl::OkStatus(); -} - -} // namespace - -int main(int argc, char** argv) { - InitGoogle(argv[0], &argc, &argv, true); - const absl::Status status = RealMain(); - if (!status.ok()) { - LOG(QFATAL) << status; - } - return 0; -} diff --git a/ortools/math_opt/samples/cutting_stock_mo.cc b/ortools/math_opt/samples/cutting_stock_mo.cc index a4b4e0ee839..58bbf8a0a27 100644 --- a/ortools/math_opt/samples/cutting_stock_mo.cc +++ b/ortools/math_opt/samples/cutting_stock_mo.cc @@ -24,43 +24,49 @@ // the second board, and so on.) // // This example approximately solves the problem with a column generation -// heuristic. The leader problem is a set cover problem, and the worker is a -// knapsack problem. We alternate between solving the LP relaxation of the -// leader incrementally, and solving the worker to generate new a configuration -// (a column) for the leader. When the worker can no longer find a column -// improving the LP cost, we convert the leader problem to a MIP and solve -// again. We now give precise statements of the leader and worker. +// heuristic. The leader problem is a set cover problem, and the worker is an +// unbounded knapsack problem. We alternate between solving the LP relaxation of +// the leader incrementally, and solving the worker to generate new a +// configuration (a column) for the leader. When the worker can no longer find a +// column improving the LP cost, we convert the leader problem to a MIP and +// solve again. We now give precise statements of the leader and worker. // // Problem data: -// * l_i: the length of each piece we need to cut out. -// * d_i: how many copies each piece we need. +// * P: the set of pieces +// * l_i: the length of each piece we need to cut out, for all i in P. +// * d_i: how many copies of each piece we need, for all i in P. // * L: the length of our initial boards. -// * q_ci: for configuration c, the quantity of piece i produced. +// * C: the set of configurations. A configuration specifies a feasible set of +// pieces to cut from a board (see q_ci below). Note that there are +// exponentially many configurations. +// * q_ci: for configuration c in C, the quantity of piece i in P to cut from a +// board (a nonnegative integer). // // Leader problem variables: -// * x_c: how many copies of configuration c to produce. +// * x_c: how many copies of configuration c in C to produce. // // Leader problem formulation: -// min sum_c x_c -// s.t. sum_c q_ci * x_c = d_i for all i -// x_c >= 0, integer for all c. +// min sum_{c in C} x_c +// s.t. sum_{c in C} q_ci * x_c = d_i, for all i in P +// x_c >= 0, integer for all c in C. // // The worker problem is to generate new configurations for the leader problem // based on the dual variables of the demand constraints in the LP relaxation. // Worker problem data: -// * p_i: The "price" of piece i (dual value from leader's demand constraint) +// * p_i: The "price" of piece i in P (dual value from leader's demand +// constraint) // // Worker decision variables: -// * y_i: How many copies of piece i should be in the configuration. +// * y_i: How many copies of piece i in P should be in the configuration. // // Worker formulation -// max sum_i p_i * y_i -// s.t. sum_i l_i * y_i <= L -// y_i >= 0, integer for all i +// max sum_{i in P} p_i * y_i +// s.t. sum_{i in P} l_i * y_i <= L +// y_i >= 0, integer for all i in P // // An optimal solution y* defines a new configuration c with q_ci = y_i* for all -// i. If the solution has objective value <= 1, no further improvement on the LP -// is possible. For additional background and proofs see: +// i in P. If the solution has objective value <= 1, no further improvement on +// the LP is possible. For additional background and proofs see: // https://people.orie.cornell.edu/shmoys/or630/notes-06/lec16.pdf // or any other reference on the "Cutting Stock Problem". // @@ -68,8 +74,10 @@ // https://en.wikipedia.org/wiki/Bin_packing_problem#Formal_statement // but typically in bin packing it is not assumed that you should exploit having // multiple items of the same size. +#include #include #include +#include #include #include @@ -86,11 +94,11 @@ namespace { namespace math_opt = operations_research::math_opt; constexpr double kInf = std::numeric_limits::infinity(); -// piece_sizes and piece_demands must have equal length. -// every piece must have 0 < size <= board_length. +// piece_lengths and piece_demands must have equal length. +// every piece must have 0 < length <= board_length. // every piece must have demand > 0. struct CuttingStockInstance { - std::vector piece_sizes; + std::vector piece_lengths; std::vector piece_demands; int board_length; }; @@ -117,30 +125,25 @@ struct CuttingStockSolution { // Solves the worker problem. // // Solves the problem on finding the configuration (with its objective value) to -// add the to model that will give the greatest improvement in the LP +// add to the leader model that will give the greatest improvement in the LP // relaxation. This is equivalent to a knapsack problem. absl::StatusOr> BestConfiguration( const std::vector& piece_prices, - const std::vector& piece_sizes, const int board_size) { + const std::vector& piece_lengths, const int board_length) { int num_pieces = piece_prices.size(); - CHECK_EQ(piece_sizes.size(), num_pieces); + CHECK_EQ(piece_lengths.size(), num_pieces); math_opt::Model model("knapsack"); std::vector pieces; for (int i = 0; i < num_pieces; ++i) { pieces.push_back( - model.AddIntegerVariable(0, kInf, absl::StrCat("item_", i))); + model.AddIntegerVariable(0, kInf, absl::StrCat("piece_", i))); } model.Maximize(math_opt::InnerProduct(pieces, piece_prices)); - model.AddLinearConstraint(math_opt::InnerProduct(pieces, piece_sizes) <= - board_size); + model.AddLinearConstraint(math_opt::InnerProduct(pieces, piece_lengths) <= + board_length); ASSIGN_OR_RETURN(const math_opt::SolveResult solve_result, math_opt::Solve(model, math_opt::SolverType::kCpSat)); - if (solve_result.termination.reason != - math_opt::TerminationReason::kOptimal) { - return util::InternalErrorBuilder() - << "Failed to solve knapsack pricing problem to optimality: " - << solve_result.termination; - } + RETURN_IF_ERROR(solve_result.termination.IsOptimal()); Configuration config; for (int i = 0; i < num_pieces; ++i) { const int use = static_cast( @@ -158,7 +161,7 @@ absl::StatusOr SolveCuttingStock( const CuttingStockInstance& instance) { math_opt::Model model("cutting_stock"); model.set_minimize(); - const int n = instance.piece_sizes.size(); + const int n = static_cast(instance.piece_lengths.size()); std::vector demand_met; for (int i = 0; i < n; ++i) { const int d = instance.piece_demands[i]; @@ -169,17 +172,17 @@ absl::StatusOr SolveCuttingStock( const math_opt::Variable v = model.AddContinuousVariable(0.0, kInf); model.set_objective_coefficient(v, 1); for (int i = 0; i < config.pieces.size(); ++i) { - const int item = config.pieces[i]; + const int piece = config.pieces[i]; const int use = config.quantity[i]; if (use >= 1) { - model.set_coefficient(demand_met[item], v, use); + model.set_coefficient(demand_met[piece], v, use); } } configs.push_back({config, v}); }; - // To ensure the leader problem is always feasible, begin a configuration for - // every item that has a single copy of the item. + // To ensure the leader problem is always feasible, begin with one + // configuration per piece, each having a single copy of the piece. for (int i = 0; i < n; ++i) { add_config(Configuration{.pieces = {i}, .quantity = {1}}); } @@ -189,13 +192,8 @@ absl::StatusOr SolveCuttingStock( int pricing_round = 0; while (true) { ASSIGN_OR_RETURN(math_opt::SolveResult solve_result, solver->Solve()); - if (solve_result.termination.reason != - math_opt::TerminationReason::kOptimal) { - return util::InternalErrorBuilder() - << "Failed to solve leader LP problem to optimality at " - "iteration " - << pricing_round << " termination: " << solve_result.termination; - } + RETURN_IF_ERROR(solve_result.termination.IsOptimal()) + << " at iteration " << pricing_round; if (!solve_result.has_dual_feasible_solution()) { // MathOpt does not require solvers to return a dual solution on optimal, // but most LP solvers always will, see go/mathopt-solver-contracts for @@ -209,9 +207,9 @@ absl::StatusOr SolveCuttingStock( for (const math_opt::LinearConstraint d : demand_met) { prices.push_back(solve_result.dual_values().at(d)); } - ASSIGN_OR_RETURN( - (const auto [config, value]), - BestConfiguration(prices, instance.piece_sizes, instance.board_length)); + ASSIGN_OR_RETURN((const auto [config, value]), + BestConfiguration(prices, instance.piece_lengths, + instance.board_length)); if (value <= 1 + 1e-3) { // The LP relaxation is solved, we can stop adding columns. break; @@ -227,15 +225,8 @@ absl::StatusOr SolveCuttingStock( } ASSIGN_OR_RETURN(const math_opt::SolveResult solve_result, math_opt::Solve(model, math_opt::SolverType::kCpSat)); - switch (solve_result.termination.reason) { - case math_opt::TerminationReason::kOptimal: - case math_opt::TerminationReason::kFeasible: - break; - default: - return util::InternalErrorBuilder() - << "Failed to solve final cutting stock MIP, termination: " - << solve_result.termination; - } + RETURN_IF_ERROR(solve_result.termination.IsOptimalOrFeasible()) + << " in final cutting stock MIP"; CuttingStockSolution solution; for (const auto& [config, var] : configs) { int use = @@ -253,12 +244,12 @@ absl::Status RealMain() { // Data from https://en.wikipedia.org/wiki/Cutting_stock_problem CuttingStockInstance instance; instance.board_length = 5600; - instance.piece_sizes = {1380, 1520, 1560, 1710, 1820, 1880, 1930, - 2000, 2050, 2100, 2140, 2150, 2200}; + instance.piece_lengths = {1380, 1520, 1560, 1710, 1820, 1880, 1930, + 2000, 2050, 2100, 2140, 2150, 2200}; instance.piece_demands = {22, 25, 12, 14, 18, 18, 20, 10, 12, 14, 16, 18, 20}; ASSIGN_OR_RETURN(CuttingStockSolution solution, SolveCuttingStock(instance)); - std::cout << "Best known solution uses 73 rolls." << std::endl; - std::cout << "Total rolls used in actual solution found: " + std::cout << "Best known solution uses 73 boards." << std::endl; + std::cout << "Total boards used in actual solution found: " << solution.objective_value << std::endl; return absl::OkStatus(); } @@ -267,10 +258,9 @@ absl::Status RealMain() { int main(int argc, char** argv) { InitGoogle(argv[0], &argc, &argv, true); - absl::Status result = RealMain(); - if (!result.ok()) { - std::cout << result; - return 1; + const absl::Status status = RealMain(); + if (!status.ok()) { + LOG(QFATAL) << status; } return 0; } diff --git a/ortools/math_opt/samples/facility_lp_benders.cc b/ortools/math_opt/samples/facility_lp_benders.cc deleted file mode 100644 index 4e50b8ef64d..00000000000 --- a/ortools/math_opt/samples/facility_lp_benders.cc +++ /dev/null @@ -1,683 +0,0 @@ -// Copyright 2010-2022 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. - -// An advanced benders decomposition example -// -// We consider a network design problem where each location has a demand that -// must be met by its neighboring facilities, and each facility can control -// its total capacity. In this version we also require that locations cannot -// use more that a specified fraction of a facilities capacity. -// -// Problem data: -// * F: set of facilities. -// * L: set of locations. -// * E: subset of {(f,l) : f in F, l in L} that describes the network between -// facilities and locations. -// * d: demand at location (all demands are equal for simplicity). -// * c: cost per unit of capacity at a facility (all facilities are have the -// same cost for simplicity). -// * h: cost per unit transported through an edge. -// * a: fraction of a facility's capacity that can be used by each location. -// -// Decision variables: -// * z_f: capacity at facility f in F. -// * x_(f,l): flow from facility f to location l for all (f,l) in E. -// -// Formulation: -// -// min c * sum(z_f : f in F) + sum(h_e * x_e : e in E) -// s.t. -// x_(f,l) <= a * z_f for all (f,l) in E -// sum(x_(f,l) : l such that (f,l) in E) <= z_f for all f in F -// sum(x_(f,l) : f such that (f,l) in E) >= d for all l in L -// x_e >= 0 for all e in E -// z_f >= 0 for all f in F -// -// Below we solve this problem directly and using a benders decompostion -// approach. - -#include -#include -#include -#include -#include -#include -#include - -#include "absl/container/flat_hash_map.h" -#include "absl/flags/flag.h" -#include "absl/random/random.h" -#include "absl/random/seed_sequences.h" -#include "absl/random/uniform_int_distribution.h" -#include "absl/status/status.h" -#include "absl/status/statusor.h" -#include "absl/strings/str_format.h" -#include "absl/strings/string_view.h" -#include "absl/time/clock.h" -#include "absl/time/time.h" -#include "ortools/base/init_google.h" -#include "ortools/base/logging.h" -#include "ortools/base/status_macros.h" -#include "ortools/math_opt/cpp/math_opt.h" -#include "ortools/util/status_macros.h" - -ABSL_FLAG(int, num_facilities, 3000, "Number of facilities."); -ABSL_FLAG(int, num_locations, 50, "Number of locations."); -ABSL_FLAG(double, edge_probability, 0.99, "Edge probability."); -ABSL_FLAG(double, benders_precission, 1e-9, "Benders target precission."); -ABSL_FLAG(double, location_demand, 1, "Client demands."); -ABSL_FLAG(double, facility_cost, 100, "Facility capacity cost."); -ABSL_FLAG( - double, location_fraction, 0.001, - "Fraction of a facility's capacity that can be used by each location."); -ABSL_FLAG(operations_research::math_opt::SolverType, solver_type, - operations_research::math_opt::SolverType::kGlop, - "the LP solver to use, possible values: glop, gurobi, glpk, pdlp"); - -namespace { - -namespace math_opt = operations_research::math_opt; -constexpr double kInf = std::numeric_limits::infinity(); -constexpr double kZeroTol = 1.0e-3; - -//////////////////////////////////////////////////////////////////////////////// -// Facility location instance representation and generation -//////////////////////////////////////////////////////////////////////////////// - -// First element is a facility and second is a location. -using Edge = std::pair; - -// A simple randomly-generated facility-location network. -class Network { - public: - Network(int num_facilities, int num_locations, double edge_probability); - - int num_facilities() const { return num_facilities_; } - int num_locations() const { return num_locations_; } - - const std::vector& edges() const { return edges_; } - - const std::vector& edges_incident_to_facility( - const int facility) const { - return facility_edge_incidence_[facility]; - } - - const std::vector& edges_incident_to_location( - const int location) const { - return location_edge_incidence_[location]; - } - - double edge_cost(const Edge& edge) const { return edge_costs_.at(edge); } - - private: - int num_facilities_; - int num_locations_; - // No order is assumed for the following lists of edges. - std::vector edges_; - absl::flat_hash_map edge_costs_; - std::vector> facility_edge_incidence_; - std::vector> location_edge_incidence_; -}; - -Network::Network(const int num_facilities, const int num_locations, - const double edge_probability) { - absl::SeedSeq seq({1, 2, 3}); - absl::BitGen bitgen(seq); - - num_facilities_ = num_facilities; - num_locations_ = num_locations; - facility_edge_incidence_ = std::vector>(num_facilities); - location_edge_incidence_ = - std::vector>(num_locations, std::vector()); - for (int facility = 0; facility < num_facilities_; ++facility) { - for (int location = 0; location < num_locations_; ++location) { - if (absl::Bernoulli(bitgen, edge_probability)) { - const Edge edge({facility, location}); - facility_edge_incidence_[facility].push_back(edge); - location_edge_incidence_[location].push_back(edge); - edges_.push_back(edge); - edge_costs_.insert({edge, absl::Uniform(bitgen, 0, 1.0)}); - } - } - } - // Ensure every facility is connected to at least one location and every - // location is connected to at least one facility. - for (int facility = 0; facility < num_facilities; ++facility) { - auto& locations = facility_edge_incidence_[facility]; - if (locations.empty()) { - const int location = - absl::uniform_int_distribution(0, num_locations - 1)(bitgen); - const std::pair edge({facility, location}); - locations.push_back(edge); - location_edge_incidence_[location].push_back(edge); - edges_.push_back(edge); - edge_costs_.insert({edge, absl::Uniform(bitgen, 0, 1.0)}); - } - } - for (int location = 0; location < num_locations; ++location) { - auto& facilities = location_edge_incidence_[location]; - if (facilities.empty()) { - const int facility = - absl::uniform_int_distribution(0, num_facilities - 1)(bitgen); - const std::pair edge({facility, location}); - facilities.push_back(edge); - facility_edge_incidence_[facility].push_back(edge); - edges_.push_back(edge); - edge_costs_.insert({edge, absl::Uniform(bitgen, 0, 1.0)}); - } - } -} - -struct FacilityLocationInstance { - Network network; - double location_demand; - double facility_cost; - double location_fraction; -}; - -//////////////////////////////////////////////////////////////////////////////// -// Direct solve -//////////////////////////////////////////////////////////////////////////////// - -// See file level comment for problem description and formulation. -absl::Status FullProblem(const FacilityLocationInstance& instance, - const math_opt::SolverType solver_type) { - const int num_facilities = instance.network.num_facilities(); - const int num_locations = instance.network.num_locations(); - - math_opt::Model model("Full network design problem"); - - // Capacity variables - std::vector z; - for (int j = 0; j < num_facilities; j++) { - z.push_back(model.AddContinuousVariable(0.0, kInf)); - } - - // Flow variables - absl::flat_hash_map x; - for (const auto& edge : instance.network.edges()) { - const math_opt::Variable x_edge = model.AddContinuousVariable(0.0, kInf); - x.insert({edge, x_edge}); - } - - // Objective Function - math_opt::LinearExpression objective_for_edges; - for (const auto& [edge, x_edge] : x) { - objective_for_edges += instance.network.edge_cost(edge) * x_edge; - } - model.Minimize(objective_for_edges + - instance.facility_cost * math_opt::Sum(z)); - - // Demand constraints - for (int location = 0; location < num_locations; ++location) { - math_opt::LinearExpression incoming_supply; - for (const auto& edge : - instance.network.edges_incident_to_location(location)) { - incoming_supply += x.at(edge); - } - model.AddLinearConstraint(incoming_supply >= instance.location_demand); - } - - // Supply constraints - for (int facility = 0; facility < num_facilities; ++facility) { - math_opt::LinearExpression outgoing_supply; - for (const auto& edge : - instance.network.edges_incident_to_facility(facility)) { - outgoing_supply += x.at(edge); - } - model.AddLinearConstraint(outgoing_supply <= z[facility]); - } - - // Arc constraints - for (int facility = 0; facility < num_facilities; ++facility) { - for (const auto& edge : - instance.network.edges_incident_to_facility(facility)) { - model.AddLinearConstraint(x.at(edge) <= - instance.location_fraction * z[facility]); - } - } - ASSIGN_OR_RETURN(const math_opt::SolveResult result, - math_opt::Solve(model, solver_type)); - if (result.termination.reason != math_opt::TerminationReason::kOptimal) { - return util::InternalErrorBuilder() - << "failed to find an optimal solution: " << result.termination; - } - - std::cout << "Full problem optimal objective: " - << absl::StrFormat("%.9f", result.objective_value()) << std::endl; - return absl::OkStatus(); -} - -//////////////////////////////////////////////////////////////////////////////// -// Benders solver -//////////////////////////////////////////////////////////////////////////////// - -// Setup first stage model: -// -// min c * sum(z_f : f in F) + w -// s.t. -// z_f >= 0 for all f in F -// sum(fcut_f^i z_f) + fcut_const^i <= 0 for i = 1,... -// sum(ocut_f^j z_f) + ocut_const^j <= w for j = 1,... -struct FirstStageProblem { - math_opt::Model model; - std::vector z; - math_opt::Variable w; - - FirstStageProblem(const Network& network, double facility_cost); -}; - -FirstStageProblem::FirstStageProblem(const Network& network, - const double facility_cost) - : model("First stage problem"), w(model.AddContinuousVariable(0.0, kInf)) { - const int num_facilities = network.num_facilities(); - - // Capacity variables - for (int j = 0; j < num_facilities; j++) { - z.push_back(model.AddContinuousVariable(0.0, kInf)); - } - - // First stage objective - model.Minimize(w + facility_cost * Sum(z)); -} - -// Represents a cut if the form: -// -// z_coefficients^T z + constant <= w_coefficient * w -// -// This will be a feasibility cut if w_coefficient = 0.0 and an optimality -// cut if w_coefficient = 1. -struct Cut { - std::vector z_coefficients; - double constant; - double w_coefficient; -}; - -// Solves the second stage model: -// -// min sum(h_e * x_e : e in E) -// s.t. -// x_(f,l) <= a * zz_f for all (f,l) in E -// sum(x_(f,l) : l such that (f,l) in E) <= zz_f for all f in F -// sum(x_(f,l) : f such that (f,l) in E) >= d for all l in L -// x_e >= 0 for all e in E -// -// where zz_f are fixed values for z_f from the first stage model, and generates -// an infeasibility or optimality cut as needed. -class SecondStageSolver { - public: - static absl::StatusOr> New( - FacilityLocationInstance instance, math_opt::SolverType solver_type); - - absl::StatusOr> Solve( - const std::vector& z_values, double w_value, - double fist_stage_objective); - - private: - SecondStageSolver(FacilityLocationInstance instance, - math_opt::SolveParameters second_stage_params); - - absl::StatusOr OptimalityCut( - const math_opt::SolveResult& second_stage_result); - absl::StatusOr FeasibilityCut( - const math_opt::SolveResult& second_stage_result); - - math_opt::Model second_stage_model_; - const Network network_; - const double location_fraction_; - math_opt::SolveParameters second_stage_params_; - - absl::flat_hash_map x_; - std::vector supply_constraints_; - std::vector demand_constraints_; - std::unique_ptr solver_; -}; - -absl::StatusOr EnsureDualRaySolveParameters( - const math_opt::SolverType solver_type) { - math_opt::SolveParameters parameters; - switch (solver_type) { - case math_opt::SolverType::kGurobi: - parameters.gurobi.param_values["InfUnbdInfo"] = "1"; - break; - case math_opt::SolverType::kGlop: - parameters.presolve = math_opt::Emphasis::kOff; - parameters.scaling = math_opt::Emphasis::kOff; - parameters.lp_algorithm = math_opt::LPAlgorithm::kDualSimplex; - break; - case math_opt::SolverType::kGlpk: - parameters.presolve = math_opt::Emphasis::kOff; - parameters.lp_algorithm = math_opt::LPAlgorithm::kDualSimplex; - parameters.glpk.compute_unbound_rays_if_possible = true; - break; - default: - return util::InternalErrorBuilder() - << "unsupported solver: " << solver_type; - } - return parameters; -} - -absl::StatusOr> SecondStageSolver::New( - FacilityLocationInstance instance, const math_opt::SolverType solver_type) { - // Set solver arguments to ensure a dual ray is returned. - ASSIGN_OR_RETURN(math_opt::SolveParameters parameters, - EnsureDualRaySolveParameters(solver_type)); - - std::unique_ptr second_stage_solver = - absl::WrapUnique( - new SecondStageSolver(std::move(instance), parameters)); - ASSIGN_OR_RETURN(std::unique_ptr solver, - math_opt::IncrementalSolver::New( - &second_stage_solver->second_stage_model_, solver_type)); - second_stage_solver->solver_ = std::move(solver); - return std::move(second_stage_solver); -} - -SecondStageSolver::SecondStageSolver( - FacilityLocationInstance instance, - math_opt::SolveParameters second_stage_params) - : second_stage_model_("Second stage model"), - network_(std::move(instance.network)), - location_fraction_(instance.location_fraction), - second_stage_params_(second_stage_params) { - const int num_facilities = network_.num_facilities(); - const int num_locations = network_.num_locations(); - - // Flow variables - for (const auto& edge : network_.edges()) { - const math_opt::Variable x_edge = - second_stage_model_.AddContinuousVariable(0.0, kInf); - x_.insert({edge, x_edge}); - } - - // Objective Function - math_opt::LinearExpression objective_for_edges; - for (const auto& [edge, x_edge] : x_) { - objective_for_edges += network_.edge_cost(edge) * x_edge; - } - second_stage_model_.Minimize(objective_for_edges); - - // Demand constraints - for (int location = 0; location < num_locations; ++location) { - math_opt::LinearExpression incoming_supply; - for (const auto& edge : network_.edges_incident_to_location(location)) { - incoming_supply += x_.at(edge); - } - demand_constraints_.push_back(second_stage_model_.AddLinearConstraint( - incoming_supply >= instance.location_demand)); - } - - // Supply constraints - for (int facility = 0; facility < num_facilities; ++facility) { - math_opt::LinearExpression outgoing_supply; - for (const auto& edge : network_.edges_incident_to_facility(facility)) { - outgoing_supply += x_.at(edge); - } - // Set supply constraint with trivial upper bound to be updated with first - // stage information. - supply_constraints_.push_back( - second_stage_model_.AddLinearConstraint(outgoing_supply <= kInf)); - } -} - -absl::StatusOr> SecondStageSolver::Solve( - const std::vector& z_values, const double w_value, - const double fist_stage_objective) { - const int num_facilities = network_.num_facilities(); - - // Update second stage with first stage solution. - for (int facility = 0; facility < num_facilities; ++facility) { - if (z_values[facility] < -kZeroTol) { - return util::InternalErrorBuilder() - << "negative z_value in first stage: " << z_values[facility] - << " for facility " << facility; - } - // Make sure variable bounds are valid (lb <= ub). - const double capacity_value = std::max(z_values[facility], 0.0); - for (const auto& edge : network_.edges_incident_to_facility(facility)) { - second_stage_model_.set_upper_bound(x_.at(edge), - location_fraction_ * capacity_value); - } - second_stage_model_.set_upper_bound(supply_constraints_[facility], - capacity_value); - } - // Solve and process second stage. - ASSIGN_OR_RETURN(const math_opt::SolveResult second_stage_result, - solver_->Solve(math_opt::SolveArguments{ - .parameters = second_stage_params_})); - switch (second_stage_result.termination.reason) { - case math_opt::TerminationReason::kInfeasible: - // If the second stage problem is infeasible we can construct a - // feasibility cut from a returned dual ray. - { - OR_ASSIGN_OR_RETURN3(const Cut feasibility_cut, - FeasibilityCut(second_stage_result), - _ << "on infeasible for second stage solver"); - return std::make_pair(kInf, feasibility_cut); - } - case math_opt::TerminationReason::kOptimal: - // If the second stage problem is optimal we can construct an optimality - // cut from a returned dual optimal solution. We can also update the upper - // bound. - { - // Upper bound is obtained by switching predicted second stage objective - // value w with the true second stage objective value. - const double upper_bound = fist_stage_objective - w_value + - second_stage_result.objective_value(); - OR_ASSIGN_OR_RETURN3(const Cut optimality_cut, - OptimalityCut(second_stage_result), - _ << "on optimal for second stage solver"); - return std::make_pair(upper_bound, optimality_cut); - } - default: - return util::InternalErrorBuilder() - << "second stage was not solved to optimality or infeasibility: " - << second_stage_result.termination; - } -} - -// If the second stage problem is infeasible we get a dual ray (r, y) such that -// -// sum(r_(f,l)*a*zz_f : (f,l) in E, r_(f,l) < 0) -// + sum(y_f*zz_f : f in F, y_f < 0) -// + sum(y_l*d : l in L, y_l > 0) > 0. -// -// Then we get the feasibility cut (go/math_opt-advanced-dual-use#benders) -// -// sum(fcut_f*z_f) + fcut_const <= 0, -// -// where -// -// fcut_f = sum(r_(f,l)*a : (f,l) in E, r_(f,l) < 0) -// + min{y_f, 0} -// fcut_const = sum*(y_l*d : l in L, y_l > 0) -absl::StatusOr SecondStageSolver::FeasibilityCut( - const math_opt::SolveResult& second_stage_result) { - const int num_facilities = network_.num_facilities(); - Cut result; - - if (!second_stage_result.has_dual_ray()) { - // MathOpt does not require solvers to return a dual ray on infeasible, - // but most LP solvers always will, see go/mathopt-solver-contracts for - // details. - return util::InternalErrorBuilder() - << "no dual ray available for feasibility cut"; - } - - for (int facility = 0; facility < num_facilities; ++facility) { - double coefficient = 0.0; - for (const auto& edge : network_.edges_incident_to_facility(facility)) { - const double reduced_cost = - second_stage_result.ray_reduced_costs().at(x_.at(edge)); - coefficient += location_fraction_ * std::min(reduced_cost, 0.0); - } - const double dual_value = - second_stage_result.ray_dual_values().at(supply_constraints_[facility]); - coefficient += std::min(dual_value, 0.0); - result.z_coefficients.push_back(coefficient); - } - result.constant = 0.0; - for (const auto& constraint : demand_constraints_) { - const double dual_value = - second_stage_result.ray_dual_values().at(constraint); - result.constant += std::max(dual_value, 0.0); - } - result.w_coefficient = 0.0; - return result; -} - -// If the second stage problem is optimal we get a dual solution (r, y) such -// that the optimal objective value is equal to -// -// sum(r_(f,l)*a*zz_f : (f,l) in E, r_(f,l) < 0) -// + sum(y_f*zz_f : f in F, y_f < 0) -// + sum*(y_l*d : l in L, y_l > 0) > 0. -// -// Then we get the optimality cut (go/math_opt-advanced-dual-use#benders) -// -// sum(ocut_f*z_f) + ocut_const <= w, -// -// where -// -// ocut_f = sum(r_(f,l)*a : (f,l) in E, r_(f,l) < 0) -// + min{y_f, 0} -// ocut_const = sum*(y_l*d : l in L, y_l > 0) -absl::StatusOr SecondStageSolver::OptimalityCut( - const math_opt::SolveResult& second_stage_result) { - const int num_facilities = network_.num_facilities(); - Cut result; - - if (!second_stage_result.has_dual_feasible_solution()) { - // MathOpt does not require solvers to return a dual solution on optimal, - // but most LP solvers always will, see go/mathopt-solver-contracts for - // details. - return util::InternalErrorBuilder() - << "no dual solution available for optimality cut"; - } - - for (int facility = 0; facility < num_facilities; ++facility) { - double coefficient = 0.0; - for (const auto& edge : network_.edges_incident_to_facility(facility)) { - const double reduced_cost = - second_stage_result.reduced_costs().at(x_.at(edge)); - coefficient += location_fraction_ * std::min(reduced_cost, 0.0); - } - double dual_value = - second_stage_result.dual_values().at(supply_constraints_[facility]); - coefficient += std::min(dual_value, 0.0); - result.z_coefficients.push_back(coefficient); - } - result.constant = 0.0; - for (const auto& constraint : demand_constraints_) { - double dual_value = second_stage_result.dual_values().at(constraint); - result.constant += std::max(dual_value, 0.0); - } - result.w_coefficient = 1.0; - return result; -} - -absl::Status Benders(const FacilityLocationInstance& instance, - const double target_precission, - const math_opt::SolverType solver_type, - const int maximum_iterations = 30000) { - const int num_facilities = instance.network.num_facilities(); - - // Setup first stage model and solver. - FirstStageProblem first_stage(instance.network, instance.facility_cost); - ASSIGN_OR_RETURN( - const std::unique_ptr first_stage_solver, - math_opt::IncrementalSolver::New(&first_stage.model, solver_type)); - // Setup second stage solver. - ASSIGN_OR_RETURN(std::unique_ptr second_stage_solver, - SecondStageSolver::New(instance, solver_type)); - // Start Benders - int iteration = 0; - double best_upper_bound = kInf; - std::vector z_values; - z_values.resize(num_facilities); - while (true) { - LOG(INFO) << "Iteration: " << iteration; - - // Solve and process first stage. - ASSIGN_OR_RETURN(const math_opt::SolveResult first_stage_result, - first_stage_solver->Solve()); - if (first_stage_result.termination.reason != - math_opt::TerminationReason::kOptimal) { - return util::InternalErrorBuilder() - << "could not solve first stage problem to optimality:" - << first_stage_result.termination; - } - for (int j = 0; j < num_facilities; j++) { - z_values[j] = first_stage_result.variable_values().at(first_stage.z[j]); - } - const double lower_bound = first_stage_result.objective_value(); - LOG(INFO) << "LB = " << lower_bound; - // Solve and process second stage. - ASSIGN_OR_RETURN( - (auto [upper_bound, cut]), - second_stage_solver->Solve( - z_values, first_stage_result.variable_values().at(first_stage.w), - first_stage_result.objective_value())); - math_opt::LinearExpression cut_expression; - for (int j = 0; j < num_facilities; j++) { - cut_expression += cut.z_coefficients[j] * first_stage.z[j]; - } - cut_expression += cut.constant; - first_stage.model.AddLinearConstraint(cut_expression <= - cut.w_coefficient * first_stage.w); - best_upper_bound = std::min(upper_bound, best_upper_bound); - LOG(INFO) << "UB = " << best_upper_bound; - ++iteration; - if (best_upper_bound - lower_bound < target_precission) { - std::cout << "Total iterations = " << iteration << std::endl; - std::cout << "Final LB = " << absl::StrFormat("%.9f", lower_bound) - << std::endl; - std::cout << "Final UB = " << absl::StrFormat("%.9f", best_upper_bound) - << std::endl; - break; - } - if (iteration > maximum_iterations) { - break; - } - } - return absl::OkStatus(); -} -absl::Status Main() { - const FacilityLocationInstance instance{ - .network = Network(absl::GetFlag(FLAGS_num_facilities), - absl::GetFlag(FLAGS_num_locations), - absl::GetFlag(FLAGS_edge_probability)), - .location_demand = absl::GetFlag(FLAGS_location_demand), - .facility_cost = absl::GetFlag(FLAGS_facility_cost), - .location_fraction = absl::GetFlag(FLAGS_location_fraction)}; - absl::Time start = absl::Now(); - RETURN_IF_ERROR(FullProblem(instance, absl::GetFlag(FLAGS_solver_type))) - << "full solve failed"; - std::cout << "Full solve time: " << absl::Now() - start << std::endl; - start = absl::Now(); - RETURN_IF_ERROR(Benders(instance, absl::GetFlag(FLAGS_benders_precission), - absl::GetFlag(FLAGS_solver_type))) - << "Benders solve failed"; - std::cout << "Benders solve time: " << absl::Now() - start << std::endl; - return absl::OkStatus(); -} -} // namespace - -int main(int argc, char** argv) { - InitGoogle(argv[0], &argc, &argv, true); - const absl::Status status = Main(); - if (!status.ok()) { - LOG(QFATAL) << status; - } - return 0; -} diff --git a/ortools/math_opt/samples/facility_lp_benders_mo.cc b/ortools/math_opt/samples/facility_lp_benders_mo.cc index 4e6c5a710c8..e6a7eb441b8 100644 --- a/ortools/math_opt/samples/facility_lp_benders_mo.cc +++ b/ortools/math_opt/samples/facility_lp_benders_mo.cc @@ -49,6 +49,8 @@ #include #include #include +#include +#include #include #include @@ -134,14 +136,13 @@ Network::Network(const int num_facilities, const int num_locations, num_facilities_ = num_facilities; num_locations_ = num_locations; - facility_edge_incidence_ = - std::vector>>(num_facilities); - location_edge_incidence_ = std::vector>>( - num_locations, std::vector>()); + facility_edge_incidence_ = std::vector>(num_facilities); + location_edge_incidence_ = + std::vector>(num_locations, std::vector()); for (int facility = 0; facility < num_facilities_; ++facility) { for (int location = 0; location < num_locations_; ++location) { if (absl::Bernoulli(bitgen, edge_probability)) { - const std::pair edge({facility, location}); + const Edge edge({facility, location}); facility_edge_incidence_[facility].push_back(edge); location_edge_incidence_[location].push_back(edge); edges_.push_back(edge); @@ -219,12 +220,12 @@ absl::Status FullProblem(const FacilityLocationInstance& instance, // Demand constraints for (int location = 0; location < num_locations; ++location) { - math_opt::LinearExpression incomming_supply; + math_opt::LinearExpression incoming_supply; for (const auto& edge : instance.network.edges_incident_to_location(location)) { - incomming_supply += x.at(edge); + incoming_supply += x.at(edge); } - model.AddLinearConstraint(incomming_supply >= instance.location_demand); + model.AddLinearConstraint(incoming_supply >= instance.location_demand); } // Supply constraints @@ -247,10 +248,7 @@ absl::Status FullProblem(const FacilityLocationInstance& instance, } ASSIGN_OR_RETURN(const math_opt::SolveResult result, math_opt::Solve(model, solver_type)); - if (result.termination.reason != math_opt::TerminationReason::kOptimal) { - return util::InternalErrorBuilder() - << "failed to find an optimal solution: " << result.termination; - } + RETURN_IF_ERROR(result.termination.IsOptimal()); std::cout << "Full problem optimal objective: " << absl::StrFormat("%.9f", result.objective_value()) << std::endl; @@ -273,7 +271,7 @@ struct FirstStageProblem { std::vector z; math_opt::Variable w; - FirstStageProblem(const Network& network, const double facility_cost); + FirstStageProblem(const Network& network, double facility_cost); }; FirstStageProblem::FirstStageProblem(const Network& network, @@ -357,6 +355,7 @@ absl::StatusOr EnsureDualRaySolveParameters( case math_opt::SolverType::kGlpk: parameters.presolve = math_opt::Emphasis::kOff; parameters.lp_algorithm = math_opt::LPAlgorithm::kDualSimplex; + parameters.glpk.compute_unbound_rays_if_possible = true; break; default: return util::InternalErrorBuilder() @@ -407,12 +406,12 @@ SecondStageSolver::SecondStageSolver( // Demand constraints for (int location = 0; location < num_locations; ++location) { - math_opt::LinearExpression incomming_supply; + math_opt::LinearExpression incoming_supply; for (const auto& edge : network_.edges_incident_to_location(location)) { - incomming_supply += x_.at(edge); + incoming_supply += x_.at(edge); } demand_constraints_.push_back(second_stage_model_.AddLinearConstraint( - incomming_supply >= instance.location_demand)); + incoming_supply >= instance.location_demand)); } // Supply constraints @@ -505,7 +504,7 @@ absl::StatusOr SecondStageSolver::FeasibilityCut( Cut result; if (!second_stage_result.has_dual_ray()) { - // MathOpt does not require solvers to return a dual solution on optimal, + // MathOpt does not require solvers to return a dual ray on infeasible, // but most LP solvers always will, see go/mathopt-solver-contracts for // details. return util::InternalErrorBuilder() @@ -556,6 +555,9 @@ absl::StatusOr SecondStageSolver::OptimalityCut( Cut result; if (!second_stage_result.has_dual_feasible_solution()) { + // MathOpt does not require solvers to return a dual solution on optimal, + // but most LP solvers always will, see go/mathopt-solver-contracts for + // details. return util::InternalErrorBuilder() << "no dual solution available for optimality cut"; } @@ -606,12 +608,8 @@ absl::Status Benders(const FacilityLocationInstance& instance, // Solve and process first stage. ASSIGN_OR_RETURN(const math_opt::SolveResult first_stage_result, first_stage_solver->Solve()); - if (first_stage_result.termination.reason != - math_opt::TerminationReason::kOptimal) { - return util::InternalErrorBuilder() - << "could not solve first stage problem to optimality:" - << first_stage_result.termination; - } + RETURN_IF_ERROR(first_stage_result.termination.IsOptimal()) + << " in first stage problem"; for (int j = 0; j < num_facilities; j++) { z_values[j] = first_stage_result.variable_values().at(first_stage.z[j]); } diff --git a/ortools/math_opt/samples/graph_coloring.cc b/ortools/math_opt/samples/graph_coloring_mo.cc similarity index 96% rename from ortools/math_opt/samples/graph_coloring.cc rename to ortools/math_opt/samples/graph_coloring_mo.cc index bfd782c1c78..fb5926ebabf 100644 --- a/ortools/math_opt/samples/graph_coloring.cc +++ b/ortools/math_opt/samples/graph_coloring_mo.cc @@ -110,14 +110,7 @@ absl::StatusOr SolveGraphColoring( // solve the model and check the result ASSIGN_OR_RETURN(const math_opt::SolveResult result, Solve(model, math_opt::SolverType::kCpSat)); - switch (result.termination.reason) { - case math_opt::TerminationReason::kOptimal: - case math_opt::TerminationReason::kFeasible: - break; - default: - return util::InternalErrorBuilder() - << "model failed to solve: " << result.termination; - } + RETURN_IF_ERROR(result.termination.IsOptimalOrFeasible()); // build solution from solver output GraphColoringSolution solution; diff --git a/ortools/math_opt/samples/integer_programming.cc b/ortools/math_opt/samples/integer_programming.cc deleted file mode 100644 index 3c2cf6969e0..00000000000 --- a/ortools/math_opt/samples/integer_programming.cc +++ /dev/null @@ -1,85 +0,0 @@ -// Copyright 2010-2022 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. - -// Simple integer programming example - -#include -#include -#include -#include - -#include "absl/status/statusor.h" -#include "absl/time/time.h" -#include "ortools/base/init_google.h" -#include "ortools/base/logging.h" -#include "ortools/base/status_builder.h" -#include "ortools/base/status_macros.h" -#include "ortools/math_opt/cpp/math_opt.h" - -namespace { - -namespace math_opt = ::operations_research::math_opt; - -constexpr double kInf = std::numeric_limits::infinity(); - -// Model and solve the problem: -// max x + 10 * y -// s.t. x + 7 * y <= 17.5 -// x <= 3.5 -// x in {0.0, 1.0, 2.0, ..., -// y in {0.0, 1.0, 2.0, ..., -// -absl::Status Main() { - math_opt::Model model("Integer programming example"); - - // Variables - const math_opt::Variable x = model.AddIntegerVariable(0.0, kInf, "x"); - const math_opt::Variable y = model.AddIntegerVariable(0.0, kInf, "y"); - - // Constraints - model.AddLinearConstraint(x + 7 * y <= 17.5, "c1"); - model.AddLinearConstraint(x <= 3.5, "c2"); - - // Objective - model.Maximize(x + 10 * y); - - ASSIGN_OR_RETURN(const math_opt::SolveResult result, - Solve(model, math_opt::SolverType::kGscip)); - - switch (result.termination.reason) { - case math_opt::TerminationReason::kOptimal: - case math_opt::TerminationReason::kFeasible: - // A feasible solution is always available on termination reason kOptimal, - // and kFeasible, but in the later case the solution may be sub-optimal. - std::cout << "Problem solved in " << result.solve_time() << std::endl; - std::cout << "Objective value: " << result.objective_value() << std::endl; - std::cout << "Variable values: [x=" - << std::round(result.variable_values().at(x)) - << ", y=" << std::round(result.variable_values().at(y)) << "]" - << std::endl; - return absl::OkStatus(); - default: - return util::InternalErrorBuilder() - << "model failed to solve: " << result.termination; - } -} -} // namespace - -int main(int argc, char** argv) { - InitGoogle(argv[0], &argc, &argv, true); - const absl::Status status = Main(); - if (!status.ok()) { - LOG(QFATAL) << status; - } - return 0; -} diff --git a/ortools/math_opt/samples/integer_programming_mo.cc b/ortools/math_opt/samples/integer_programming_mo.cc index ae1c3e798b6..b64ad201991 100644 --- a/ortools/math_opt/samples/integer_programming_mo.cc +++ b/ortools/math_opt/samples/integer_programming_mo.cc @@ -13,8 +13,10 @@ // Simple integer programming example +#include #include #include +#include #include "absl/status/statusor.h" #include "absl/time/time.h" @@ -53,19 +55,16 @@ absl::Status Main() { ASSIGN_OR_RETURN(const math_opt::SolveResult result, Solve(model, math_opt::SolverType::kGscip)); - - switch (result.termination.reason) { - case math_opt::TerminationReason::kOptimal: - case math_opt::TerminationReason::kFeasible: - std::cout << "Problem solved in " << result.solve_time() << std::endl; - std::cout << "Objective value: " << result.objective_value() << std::endl; - std::cout << "Variable values: [x=" << result.variable_values().at(x) - << ", y=" << result.variable_values().at(y) << "]" << std::endl; - return absl::OkStatus(); - default: - return util::InternalErrorBuilder() - << "model failed to solve: " << result.termination; - } + RETURN_IF_ERROR(result.termination.IsOptimalOrFeasible()); + // A feasible solution is always available on termination reason kOptimal, and + // kFeasible, but in the later case the solution may be sub-optimal. + std::cout << "Problem solved in " << result.solve_time() << std::endl; + std::cout << "Objective value: " << result.objective_value() << std::endl; + std::cout << "Variable values: [x=" + << std::round(result.variable_values().at(x)) + << ", y=" << std::round(result.variable_values().at(y)) << "]" + << std::endl; + return absl::OkStatus(); } } // namespace diff --git a/ortools/math_opt/samples/lagrangian_relaxation.cc b/ortools/math_opt/samples/lagrangian_relaxation.cc deleted file mode 100644 index 0638e42f3ce..00000000000 --- a/ortools/math_opt/samples/lagrangian_relaxation.cc +++ /dev/null @@ -1,493 +0,0 @@ -// Copyright 2010-2022 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. - -// Solves a constrained shortest path problem via Lagrangian Relaxation. The -// Lagrangian dual is solved with subgradient ascent. -// -// Problem data: -// * N: set of nodes. -// * A: set of arcs. -// * R: set of resources. -// * c_(i,j): cost of traversing arc (i,j) in A. -// * r_(i,j,k): resource k spent by traversing arc (i,j) in A, for all k in R. -// * b_i: flow balance at node i in N (+1 at the source, -1 at the sink, and 0 -// otherwise). -// * r_max_k: availability of resource k for a path, for all k in R. -// -// Decision variables: -// * x_(i,j): flow through arc (i,j) in A. -// -// Formulation: -// Z = min sum(c_(i,j) * x_(i,j): (i,j) in A) -// s.t. -// sum(x_(i,j): (i,j) in A) - sum(x_(j,i): (j,i) in A) = b_i for all i in N, -// sum(r_(i,j,k) * x_(i,j): (i,j) in A) <= r_max_k for all k in R, -// x_(i,j) in {0,1} for all (i,j) in A. -// -// Upon dualizing a subset of the constraints (here we chose to relax some or -// all of the knapsack constraints), we obtaing a subproblem parameterized by -// dual variables mu (one per dualized constraint). We refer to this as the -// Lagrangian subproblem. Let R+ be the set of knapsack constraints that we -// keep, and R- the set of knapsack constraints that get dualized. The -// Lagrangian subproblem follows: -// -// z(mu) = min sum( -// (c_(i,j) - sum(mu_k * r_(i,j,k): k in R)) * x_(i,j): (i,j) in A) -// + sum(mu_k * r_max_k: k in R-) -// s.t. -// sum(x_(i,j): (i,j) in A) - sum(x_(j,i): (j,i) in A) = b_i for all i in N, -// sum(r_(i,j,k) * x_(i,j): (i,j) in A) <= r_max_k for all k in R+, -// x_(i,j) in {0,1} for all (i,j) in A. -// -// We seek to solve the Lagrangian dual, which is of the form: -// Z_D = max{ z(mu) : mu <=0 }. Concavity of z(mu) allows us to solve the -// Lagrangian dual with the iterates: -// mu_(t+1) = mu_t + step_size_t * grad_mu_t, where -// grad_mu_t = r_max - sum(t_(i,j) * x_(i,j)^t: (i,j) in A) is a subgradient of -// z(mu_t) and x^t is an optimal solution to the problem induced by z(mu_t). -// -// In general we have that Z_D <= Z. For convex problems, Z_D = Z. For MIPs, -// Z_LP <= Z_D <= Z, where Z_LP is the linear relaxation of the original -// problem. -// -// In this particular example, we use two resource constraints. Either -// constraint or both can be dualized via the flags `dualize_resource_1` and -// `dualize_resource_2`. If both constraints are dualized, we have that Z_LP = -// Z_D because the resulting Lagrangian subproblem can be solved as a linear -// program (i.e., the problem becomes a pure shortest path problem upon -// dualizing all the side constraints). When only one of the side constraints is -// dualized, we can have Z_LP <= Z_D because the resulting Lagrangian subproblem -// needs to be solved as an MIP. For the particular data used in this example, -// dualizing only the first resource constraint leads to Z_LP < Z_D, while -// dualizing only the second resource constraint leads to Z_LP = Z_D. In either -// case, solving the Lagrandual dual also provides an upper bound to Z. -// -// Usage: blaze build -c opt -// ortools/math_opt/examples:lagrangian_relaxation -// blaze-bin/ortools/math_opt/examples/lagrangian_relaxation - -#include - -#include -#include -#include -#include -#include -#include -#include - -#include "absl/flags/flag.h" -#include "absl/memory/memory.h" -#include "absl/status/status.h" -#include "absl/status/statusor.h" -#include "absl/strings/str_format.h" -#include "absl/strings/string_view.h" -#include "ortools/base/container_logging.h" -#include "ortools/base/init_google.h" -#include "ortools/base/logging.h" -#include "ortools/base/mathutil.h" -#include "ortools/base/status_builder.h" -#include "ortools/base/status_macros.h" -#include "ortools/math_opt/cpp/math_opt.h" - -ABSL_FLAG(double, step_size, 0.95, - "Stepsize for gradient ascent, determined as step_size^t."); -ABSL_FLAG(int, max_iterations, 1000, - "Max number of iterations for gradient ascent."); -ABSL_FLAG(bool, dualize_resource_1, true, - "If true, the side constraint associated to resource 1 is dualized."); -ABSL_FLAG(bool, dualize_resource_2, false, - "If true, the side constraint associated to resource 2 is dualized."); - -ABSL_FLAG(bool, lagrangian_output, false, - "If true, shows the iteration log of the subgradient ascent " - "procedure use to solve the Lagrangian problem"); - -constexpr double kZeroTol = 1.0e-8; - -namespace { -using ::operations_research::MathUtil; - -namespace math_opt = ::operations_research::math_opt; - -struct Arc { - int i; - int j; - double cost; - double resource_1; - double resource_2; -}; - -struct Graph { - int num_nodes; - std::vector arcs; - int source; - int sink; -}; - -struct FlowModel { - FlowModel() : model(std::make_unique("LagrangianProblem")) {} - std::unique_ptr model; - math_opt::LinearExpression cost; - math_opt::LinearExpression resource_1; - math_opt::LinearExpression resource_2; - std::vector flow_vars; -}; - -// Populates `model` with variables and constraints of a shortest path problem. -FlowModel CreateShortestPathModel(const Graph graph) { - FlowModel flow_model; - math_opt::Model& model = *flow_model.model; - for (const Arc& arc : graph.arcs) { - math_opt::Variable var = model.AddContinuousVariable( - /*lower_bound=*/0, /*upper_bound=*/1, - /*name=*/absl::StrFormat("x_%d_%d", arc.i, arc.j)); - flow_model.cost += arc.cost * var; - flow_model.resource_1 += arc.resource_1 * var; - flow_model.resource_2 += arc.resource_2 * var; - flow_model.flow_vars.push_back(var); - } - - // Flow balance constraints - std::vector out_flow(graph.num_nodes); - std::vector in_flow(graph.num_nodes); - for (int arc_index = 0; arc_index < graph.arcs.size(); ++arc_index) { - out_flow[graph.arcs[arc_index].i] += flow_model.flow_vars[arc_index]; - in_flow[graph.arcs[arc_index].j] += flow_model.flow_vars[arc_index]; - } - for (int node_index = 0; node_index < graph.num_nodes; ++node_index) { - int rhs = node_index == graph.source ? 1 - : node_index == graph.sink ? -1 - : 0; - model.AddLinearConstraint(out_flow[node_index] - in_flow[node_index] == - rhs); - } - - return flow_model; -} - -Graph CreateSampleNetwork() { - std::vector arcs; - arcs.push_back( - {.i = 0, .j = 1, .cost = 12, .resource_1 = 1, .resource_2 = 1}); - arcs.push_back( - {.i = 0, .j = 2, .cost = 3, .resource_1 = 2.5, .resource_2 = 1}); - arcs.push_back( - {.i = 1, .j = 3, .cost = 5, .resource_1 = 1, .resource_2 = 1.5}); - arcs.push_back( - {.i = 1, .j = 4, .cost = 5, .resource_1 = 2.5, .resource_2 = 1}); - arcs.push_back( - {.i = 2, .j = 1, .cost = 7, .resource_1 = 2.5, .resource_2 = 1}); - arcs.push_back( - {.i = 2, .j = 3, .cost = 5, .resource_1 = 7, .resource_2 = 2.5}); - arcs.push_back( - {.i = 2, .j = 4, .cost = 1, .resource_1 = 6.5, .resource_2 = 1}); - arcs.push_back( - {.i = 3, .j = 5, .cost = 6, .resource_1 = 1, .resource_2 = 2.0}); - arcs.push_back( - {.i = 4, .j = 3, .cost = 3, .resource_1 = 1, .resource_2 = 0.5}); - arcs.push_back( - {.i = 4, .j = 5, .cost = 5, .resource_1 = 2.5, .resource_2 = 1}); - const Graph graph = {.num_nodes = 6, .arcs = arcs, .source = 0, .sink = 5}; - - return graph; -} - -// Solves the constrained shortest path as an MIP. -absl::StatusOr SolveMip(const Graph graph, - const double max_resource_1, - const double max_resource_2) { - FlowModel flow_model; - math_opt::Model& model = *flow_model.model; - for (const Arc& arc : graph.arcs) { - math_opt::Variable var = model.AddBinaryVariable( - /*name=*/absl::StrFormat("x_%d_%d", arc.i, arc.j)); - flow_model.cost += arc.cost * var; - flow_model.resource_1 += +arc.resource_1 * var; - flow_model.resource_2 += arc.resource_2 * var; - flow_model.flow_vars.push_back(var); - } - - // Flow balance constraints - std::vector out_flow(graph.num_nodes); - std::vector in_flow(graph.num_nodes); - for (int arc_index = 0; arc_index < graph.arcs.size(); ++arc_index) { - out_flow[graph.arcs[arc_index].i] += flow_model.flow_vars[arc_index]; - in_flow[graph.arcs[arc_index].j] += flow_model.flow_vars[arc_index]; - } - for (int node_index = 0; node_index < graph.num_nodes; ++node_index) { - int rhs = node_index == graph.source ? 1 - : node_index == graph.sink ? -1 - : 0; - model.AddLinearConstraint(out_flow[node_index] - in_flow[node_index] == - rhs); - } - - model.AddLinearConstraint(flow_model.resource_1 <= max_resource_1, - "resource_ctr_1"); - model.AddLinearConstraint(flow_model.resource_2 <= max_resource_2, - "resource_ctr_2"); - model.Minimize(flow_model.cost); - ASSIGN_OR_RETURN(const math_opt::SolveResult result, - Solve(model, math_opt::SolverType::kGscip)); - switch (result.termination.reason) { - case math_opt::TerminationReason::kOptimal: - case math_opt::TerminationReason::kFeasible: - std::cout << "MIP Solution with 2 side constraints" << std::endl; - std::cout << absl::StrFormat("MIP objective value: %6.3f", - result.objective_value()) - << std::endl; - std::cout << "Resource 1: " - << flow_model.resource_1.Evaluate(result.variable_values()) - << std::endl; - std::cout << "Resource 2: " - << flow_model.resource_2.Evaluate(result.variable_values()) - << std::endl; - std::cout << "========================================" << std::endl; - return flow_model; - default: - return util::InternalErrorBuilder() - << "model failed to solve: " << result.termination; - } -} - -// Solves the linear relaxation of a constrained shortest path problem -// formulated as an MIP. -absl::Status SolveLinearRelaxation(FlowModel& flow_model, const Graph& graph, - const double max_resource_1, - const double max_resource_2) { - math_opt::Model& model = *flow_model.model; - ASSIGN_OR_RETURN(const math_opt::SolveResult result, - Solve(model, math_opt::SolverType::kGscip)); - switch (result.termination.reason) { - case math_opt::TerminationReason::kOptimal: - case math_opt::TerminationReason::kFeasible: - std::cout << "LP relaxation with 2 side constraints" << std::endl; - std::cout << absl::StrFormat("LP objective value: %6.3f", - result.objective_value()) - << std::endl; - std::cout << "Resource 1: " - << flow_model.resource_1.Evaluate(result.variable_values()) - << std::endl; - std::cout << "Resource 2: " - << flow_model.resource_2.Evaluate(result.variable_values()) - << std::endl; - std::cout << "========================================" << std::endl; - return absl::OkStatus(); - default: - return util::InternalErrorBuilder() - << "model failed to solve: " << result.termination; - } -} - -absl::Status SolveLagrangianRelaxation(const Graph graph, - const double max_resource_1, - const double max_resource_2) { - // Model, variables, and linear expressions. - FlowModel flow_model = CreateShortestPathModel(graph); - math_opt::Model& model = *flow_model.model; - math_opt::LinearExpression& cost = flow_model.cost; - math_opt::LinearExpression& resource_1 = flow_model.resource_1; - math_opt::LinearExpression& resource_2 = flow_model.resource_2; - - // Dualized constraints and dual variable iterates. - std::vector mu; - std::vector grad_mu; - const bool dualized_resource_1 = absl::GetFlag(FLAGS_dualize_resource_1); - const bool dualized_resource_2 = absl::GetFlag(FLAGS_dualize_resource_2); - const bool lagrangian_output = absl::GetFlag(FLAGS_lagrangian_output); - CHECK(dualized_resource_1 || dualized_resource_2) - << "At least one of the side constraints should be dualized."; - - // Modify the lagrangian problem according to the constraints that are - // dualized. We use a initial dual value different from zero to prioritize - // finding a feasible solution. - const double initial_dual_value = -10; - if (dualized_resource_1 && !dualized_resource_2) { - mu.push_back(initial_dual_value); - grad_mu.push_back(max_resource_1 - resource_1); - model.AddLinearConstraint(resource_2 <= max_resource_2); - for (math_opt::Variable& var : flow_model.flow_vars) { - model.set_integer(var); - } - } else if (!dualized_resource_1 && dualized_resource_2) { - mu.push_back(initial_dual_value); - grad_mu.push_back(max_resource_2 - resource_2); - model.AddLinearConstraint(resource_1 <= max_resource_1); - for (math_opt::Variable& var : flow_model.flow_vars) { - model.set_integer(var); - } - } else { - mu.push_back(initial_dual_value); - mu.push_back(initial_dual_value); - grad_mu.push_back(max_resource_1 - resource_1); - grad_mu.push_back(max_resource_2 - resource_2); - } - - // Gradient ascent setup - bool termination = false; - int iterations = 1; - const double step_size = absl::GetFlag(FLAGS_step_size); - CHECK_GT(step_size, 0) << "step_size must be strictly positive"; - CHECK_LT(step_size, 1) << "step_size must be strictly less than 1"; - const int max_iterations = absl::GetFlag(FLAGS_max_iterations); - CHECK_GT(max_iterations, 0) - << "Number of iterations must be strictly positive."; - - // Upper and lower bounds on the full problem. - double upper_bound = std::numeric_limits::infinity(); - double lower_bound = -std::numeric_limits::infinity(); - double best_solution_resource_1 = 0; - double best_solution_resource_2 = 0; - - if (lagrangian_output) { - std::cout << "Starting gradient ascent..." << std::endl; - std::cout << absl::StrFormat("%4s %6s %6s %9s %10s %10s", "Iter", "LB", - "UB", "Step size", "mu_t", "grad_mu_t") - << std::endl; - } - - while (!termination) { - math_opt::LinearExpression lagrangian_function; - lagrangian_function += cost; - for (int k = 0; k < mu.size(); ++k) { - lagrangian_function += mu[k] * grad_mu[k]; - } - model.Minimize(lagrangian_function); - ASSIGN_OR_RETURN(math_opt::SolveResult result, - Solve(model, math_opt::SolverType::kGscip)); - switch (result.termination.reason) { - case math_opt::TerminationReason::kOptimal: - case math_opt::TerminationReason::kFeasible: - break; - default: - return util::InternalErrorBuilder() - << "failed to minimize lagrangian function: " - << result.termination; - } - - const math_opt::VariableMap& vars_val = result.variable_values(); - bool feasible = true; - - // Iterate update. Takes a step in the direction of the gradient (since the - // Lagrangian dual is a max problem), and projects onto {mu: mu <=0} to - // satisfy the sign of the dual variable. In general, convergence to an - // optimal solution requires diminishing step sizes satisfying: - // * sum(step_size_t: t=1...) = infinity and, - // * sum((step_size_t)^2: t=1...) < infinity - // See details in Prop 3.2.6 Bertsekas 2015, Convex Optimization Algorithms. - // Here we use step_size_t = step_size^t which does NOT satisfy the - // first condition, but is good enough for the purpose of this example. - std::vector grad_mu_vals; - const double step_size_t = MathUtil::IPow(step_size, iterations); - for (int k = 0; k < mu.size(); ++k) { - // Evaluate resource k and evaluate the gradient of z(mu). - const double grad_mu_val = grad_mu[k].Evaluate(vars_val); - grad_mu_vals.push_back(grad_mu_val); - mu[k] = std::min(0.0, mu[k] + step_size_t * grad_mu_val); - if (grad_mu_val < 0) { - feasible = false; - } - } - // Bounds update - const double path_cost = cost.Evaluate(vars_val); - if (feasible && path_cost < upper_bound) { - best_solution_resource_1 = resource_1.Evaluate(vars_val); - best_solution_resource_2 = resource_2.Evaluate(vars_val); - if (lagrangian_output) { - std::cout << "Feasible solution with" - << absl::StrFormat( - "cost=%4.2f, resource_1=%4.2f, and resource_2=%4.2f. ", - path_cost, best_solution_resource_1, - best_solution_resource_2) - << std::endl; - } - upper_bound = path_cost; - } - if (lower_bound < result.objective_value()) { - lower_bound = result.objective_value(); - } - - if (lagrangian_output) { - std::cout << absl::StrFormat("%4d %6.3f %6.3f %9.3f", iterations, - lower_bound, upper_bound, step_size_t) - << " " << gtl::LogContainer(mu) << " " - << gtl::LogContainer(grad_mu_vals) << std::endl; - } - - // Termination criteria - double norm = 0; - for (double value : grad_mu_vals) { - norm += (value * value); - } - norm = sqrt(norm); - if (iterations == max_iterations || lower_bound == upper_bound || - step_size_t * norm < kZeroTol) { - termination = true; - } - iterations++; - } - - std::cout << "Lagrangian relaxation with 2 side constraints" << std::endl; - std::cout << "Constraint for resource 1 dualized: " - << (dualized_resource_1 ? "true" : "false") << std::endl; - std::cout << "Constraint for resource 2 dualized: " - << (dualized_resource_2 ? "true" : "false") << std::endl; - std::cout << absl::StrFormat("Lower bound: %6.3f", lower_bound) << std::endl; - std::cout << absl::StrFormat("Upper bound: %6.3f (Integer solution)", - upper_bound) - << std::endl; - std::cout << "========================================" << std::endl; - return absl::OkStatus(); -} - -void RelaxModel(FlowModel& flow_model) { - for (math_opt::Variable& var : flow_model.flow_vars) { - flow_model.model->set_continuous(var); - flow_model.model->set_lower_bound(var, 0.0); - flow_model.model->set_upper_bound(var, 1.0); - } -} - -absl::Status SolveFullModel(const Graph& graph, double max_resource_1, - double max_resource_2) { - ASSIGN_OR_RETURN(FlowModel flow_model, - SolveMip(graph, max_resource_1, max_resource_2)); - RelaxModel(flow_model); - return SolveLinearRelaxation(flow_model, graph, max_resource_1, - max_resource_2); -} - -absl::Status Main() { - // Problem data - const Graph graph = CreateSampleNetwork(); - const double max_resource_1 = 10; - const double max_resource_2 = 4; - - RETURN_IF_ERROR(SolveFullModel(graph, max_resource_1, max_resource_2)) - << "full solve failed"; - RETURN_IF_ERROR( - SolveLagrangianRelaxation(graph, max_resource_1, max_resource_2)) - << "lagrangian solve failed"; - return absl::OkStatus(); -} -} // namespace - -int main(int argc, char** argv) { - InitGoogle(argv[0], &argc, &argv, true); - const absl::Status status = Main(); - if (!status.ok()) { - LOG(QFATAL) << status; - } - return 0; -} diff --git a/ortools/math_opt/samples/lagrangian_relaxation_mo.cc b/ortools/math_opt/samples/lagrangian_relaxation_mo.cc index f04bd2a8c37..f4c43bd1c71 100644 --- a/ortools/math_opt/samples/lagrangian_relaxation_mo.cc +++ b/ortools/math_opt/samples/lagrangian_relaxation_mo.cc @@ -82,6 +82,7 @@ #include #include #include +#include #include #include @@ -239,25 +240,19 @@ absl::StatusOr SolveMip(const Graph graph, model.Minimize(flow_model.cost); ASSIGN_OR_RETURN(const math_opt::SolveResult result, Solve(model, math_opt::SolverType::kGscip)); - switch (result.termination.reason) { - case math_opt::TerminationReason::kOptimal: - case math_opt::TerminationReason::kFeasible: - std::cout << "MIP Solution with 2 side constraints" << std::endl; - std::cout << absl::StrFormat("MIP objective value: %6.3f", - result.objective_value()) - << std::endl; - std::cout << "Resource 1: " - << flow_model.resource_1.Evaluate(result.variable_values()) - << std::endl; - std::cout << "Resource 2: " - << flow_model.resource_2.Evaluate(result.variable_values()) - << std::endl; - std::cout << "========================================" << std::endl; - return flow_model; - default: - return util::InternalErrorBuilder() - << "model failed to solve: " << result.termination; - } + RETURN_IF_ERROR(result.termination.IsOptimalOrFeasible()); + std::cout << "MIP Solution with 2 side constraints" << std::endl; + std::cout << absl::StrFormat("MIP objective value: %6.3f", + result.objective_value()) + << std::endl; + std::cout << "Resource 1: " + << flow_model.resource_1.Evaluate(result.variable_values()) + << std::endl; + std::cout << "Resource 2: " + << flow_model.resource_2.Evaluate(result.variable_values()) + << std::endl; + std::cout << "========================================" << std::endl; + return flow_model; } // Solves the linear relaxation of a constrained shortest path problem @@ -268,25 +263,19 @@ absl::Status SolveLinearRelaxation(FlowModel& flow_model, const Graph& graph, math_opt::Model& model = *flow_model.model; ASSIGN_OR_RETURN(const math_opt::SolveResult result, Solve(model, math_opt::SolverType::kGscip)); - switch (result.termination.reason) { - case math_opt::TerminationReason::kOptimal: - case math_opt::TerminationReason::kFeasible: - std::cout << "LP relaxation with 2 side constraints" << std::endl; - std::cout << absl::StrFormat("LP objective value: %6.3f", - result.objective_value()) - << std::endl; - std::cout << "Resource 1: " - << flow_model.resource_1.Evaluate(result.variable_values()) - << std::endl; - std::cout << "Resource 2: " - << flow_model.resource_2.Evaluate(result.variable_values()) - << std::endl; - std::cout << "========================================" << std::endl; - return absl::OkStatus(); - default: - return util::InternalErrorBuilder() - << "model failed to solve: " << result.termination; - } + RETURN_IF_ERROR(result.termination.IsOptimalOrFeasible()); + std::cout << "LP relaxation with 2 side constraints" << std::endl; + std::cout << absl::StrFormat("LP objective value: %6.3f", + result.objective_value()) + << std::endl; + std::cout << "Resource 1: " + << flow_model.resource_1.Evaluate(result.variable_values()) + << std::endl; + std::cout << "Resource 2: " + << flow_model.resource_2.Evaluate(result.variable_values()) + << std::endl; + std::cout << "========================================" << std::endl; + return absl::OkStatus(); } absl::Status SolveLagrangianRelaxation(const Graph graph, @@ -365,15 +354,7 @@ absl::Status SolveLagrangianRelaxation(const Graph graph, model.Minimize(lagrangian_function); ASSIGN_OR_RETURN(math_opt::SolveResult result, Solve(model, math_opt::SolverType::kGscip)); - switch (result.termination.reason) { - case math_opt::TerminationReason::kOptimal: - case math_opt::TerminationReason::kFeasible: - break; - default: - return util::InternalErrorBuilder() - << "failed to minimize lagrangian function: " - << result.termination; - } + RETURN_IF_ERROR(result.termination.IsOptimalOrFeasible()); const math_opt::VariableMap& vars_val = result.variable_values(); bool feasible = true; diff --git a/ortools/math_opt/samples/linear_programming.cc b/ortools/math_opt/samples/linear_programming.cc deleted file mode 100644 index aa427ff6feb..00000000000 --- a/ortools/math_opt/samples/linear_programming.cc +++ /dev/null @@ -1,93 +0,0 @@ -// Copyright 2010-2022 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. - -// Simple linear programming example - -#include -#include -#include -#include -#include - -#include "absl/status/statusor.h" -#include "absl/strings/str_cat.h" -#include "absl/strings/str_join.h" -#include "absl/time/time.h" -#include "ortools/base/init_google.h" -#include "ortools/base/logging.h" -#include "ortools/base/status_builder.h" -#include "ortools/base/status_macros.h" -#include "ortools/math_opt/cpp/math_opt.h" - -namespace { - -namespace math_opt = ::operations_research::math_opt; - -constexpr double kInf = std::numeric_limits::infinity(); - -// Model and solve the problem: -// max 10 * x0 + 6 * x1 + 4 * x2 -// s.t. 10 * x0 + 4 * x1 + 5 * x2 <= 600 -// 2 * x0 + 2 * x1 + 6 * x2 <= 300 -// x0 + x1 + x2 <= 100 -// x0 in [0, infinity) -// x1 in [0, infinity) -// x2 in [0, infinity) -// -absl::Status Main() { - math_opt::Model model("Linear programming example"); - - // Variables - std::vector x; - for (int j = 0; j < 3; j++) { - x.push_back(model.AddContinuousVariable(0.0, kInf, absl::StrCat("x", j))); - } - - // Constraints - std::vector constraints; - constraints.push_back( - model.AddLinearConstraint(10 * x[0] + 4 * x[1] + 5 * x[2] <= 600, "c1")); - constraints.push_back( - model.AddLinearConstraint(2 * x[0] + 2 * x[1] + 6 * x[2] <= 300, "c2")); - // sum(x[i]) <= 100 - constraints.push_back(model.AddLinearConstraint(Sum(x) <= 100, "c3")); - - // Objective - model.Maximize(10 * x[0] + 6 * x[1] + 4 * x[2]); - - ASSIGN_OR_RETURN(const math_opt::SolveResult result, - Solve(model, math_opt::SolverType::kGlop)); - if (result.termination.reason != math_opt::TerminationReason::kOptimal) { - return util::InternalErrorBuilder() - << "model failed to solve to optimality" << result.termination; - } - - std::cout << "Problem solved in " << result.solve_time() << std::endl; - std::cout << "Objective value: " << result.objective_value() << std::endl; - - std::cout << "Variable values: [" - << absl::StrJoin(Values(result.variable_values(), x), ", ") << "]" - << std::endl; - - return absl::OkStatus(); -} -} // namespace - -int main(int argc, char** argv) { - InitGoogle(argv[0], &argc, &argv, true); - const absl::Status status = Main(); - if (!status.ok()) { - LOG(QFATAL) << status; - } - return 0; -} diff --git a/ortools/math_opt/samples/linear_programming_mo.cc b/ortools/math_opt/samples/linear_programming_mo.cc index 96ab39845b4..31cbe792742 100644 --- a/ortools/math_opt/samples/linear_programming_mo.cc +++ b/ortools/math_opt/samples/linear_programming_mo.cc @@ -15,6 +15,7 @@ #include #include +#include #include #include @@ -66,16 +67,13 @@ absl::Status Main() { ASSIGN_OR_RETURN(const math_opt::SolveResult result, Solve(model, math_opt::SolverType::kGlop)); - if (result.termination.reason != math_opt::TerminationReason::kOptimal) { - return util::InternalErrorBuilder() - << "model failed to solve to optimality" << result.termination; - } + RETURN_IF_ERROR(result.termination.IsOptimal()); std::cout << "Problem solved in " << result.solve_time() << std::endl; std::cout << "Objective value: " << result.objective_value() << std::endl; std::cout << "Variable values: [" - << absl::StrJoin(result.variable_values().Values(x), ", ") << "]" + << absl::StrJoin(Values(result.variable_values(), x), ", ") << "]" << std::endl; return absl::OkStatus(); diff --git a/ortools/math_opt/samples/linear_regression.cc b/ortools/math_opt/samples/linear_regression_mo.cc similarity index 97% rename from ortools/math_opt/samples/linear_regression.cc rename to ortools/math_opt/samples/linear_regression_mo.cc index 991137fc4ef..e376e3a3bd2 100644 --- a/ortools/math_opt/samples/linear_regression.cc +++ b/ortools/math_opt/samples/linear_regression_mo.cc @@ -163,11 +163,7 @@ absl::StatusOr Train( args.parameters.enable_output = true; ASSIGN_OR_RETURN(const math_opt::SolveResult result, Solve(model, absl::GetFlag(FLAGS_solver_type), args)); - if (result.termination.reason != math_opt::TerminationReason::kOptimal) { - return util::InternalErrorBuilder() - << "Expected termination reason optimal, but termination was: " - << result.termination; - } + RETURN_IF_ERROR(result.termination.IsOptimal()); std::cout << "Training time: " << result.solve_time() << std::endl; return LinearModel{.betas = Values(result.variable_values(), betas)}; } diff --git a/ortools/math_opt/samples/time_indexed_scheduling.cc b/ortools/math_opt/samples/time_indexed_scheduling_mo.cc similarity index 100% rename from ortools/math_opt/samples/time_indexed_scheduling.cc rename to ortools/math_opt/samples/time_indexed_scheduling_mo.cc diff --git a/ortools/math_opt/samples/tsp.cc b/ortools/math_opt/samples/tsp.cc deleted file mode 100644 index 3c6ac26a333..00000000000 --- a/ortools/math_opt/samples/tsp.cc +++ /dev/null @@ -1,359 +0,0 @@ -// Copyright 2010-2022 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. - -// A minimal TSP solver using MathOpt. -// -// In the Euclidean Traveling Salesperson Problem (TSP), you are given a list of -// n cities, each with an (x, y) coordinate, and you must find an order to visit -// the cities in to minimize the (Euclidean) travel distance. -// -// The MIP "cutset" formulation for the problem is as follows: -// * Data: -// n: An integer, the number of cities -// (x_i, y_i): a pair of floats for each i in 1..n, the location of each -// city -// d_ij for all (i, j) pairs of cities, the distance between city i and j. -// * Decision variables: -// x_ij: A binary variable, indicates if the edge connecting i and j is -// used. Note that x_ij == x_ji, because the problem is symmetric. We -// only create variables for i < j, and have x_ji as an alias for -// x_ij. -// * MIP model: -// minimize sum_{i=1}^n sum_{j=1, j < i}^n d_ij * x_ij -// s.t. sum_{j=1, j != i}^n x_ij = 2 for all i = 1..n -// sum_{i in S} sum_{j not in S} x_ij >= 2 for all S subset {1,...,n} -// |S| >= 3, |S| <= n - 3 -// x_ij in {0, 1} -// The first set of constraints are called the degree constraints, and the -// second set of constraints are called the cutset constraints. There are -// exponentially many cutset, so we cannot add them all at the start of the -// solve. Instead, we will use a solver callback to view each integer solution -// and add any violated cutset constraints that exist. -// -// Note that, while there are exponentially many cutset constraints, we can -// quickly identify violated ones by exploiting that the solution is integer -// and the degree constraints are all already in the model and satisfied. As a -// result, the graph n nodes and edges when x_ij = 1 will be a degree two graph, -// so it will be a collection of cycles. If it is a single large cycle, then the -// solution is feasible, and if there multiple cycles, then taking the nodes of -// any cycle as S produces a violated cutset constraint. -// -// Note that this is a minimal TSP solution, more sophisticated MIP methods are -// possible. - -#include -#include -#include -#include -#include -#include -#include - -#include "absl/flags/flag.h" -#include "absl/random/random.h" -#include "absl/random/uniform_real_distribution.h" -#include "absl/strings/str_cat.h" -#include "absl/strings/str_join.h" -#include "ortools/base/helpers.h" -#include "ortools/base/init_google.h" -#include "ortools/base/logging.h" -#include "ortools/base/status_builder.h" -#include "ortools/base/status_macros.h" -#include "ortools/math_opt/cpp/math_opt.h" -#include "ortools/port/proto_utils.h" - -ABSL_FLAG(int, num_cities, 50, "Number of cities in random TSP."); -ABSL_FLAG(std::string, output, "", - "Write a svg of the solution here, or to standard out if empty."); -ABSL_FLAG(bool, test_instance, false, - "Solve the test TSP instead of a random instance."); -ABSL_FLAG(int, threads, 0, - "How many threads to solve with, or solver default if <= 0."); -ABSL_FLAG(bool, solve_logs, false, - "Have the solver print logs to standard out."); - -namespace { - -namespace math_opt = operations_research::math_opt; -using Cycle = std::vector; - -// Creates variables modeling the undirected edges for the TSP. For every (i, j) -// pair in [0,n) * [0, n), a variable is created only for j < i, but querying -// for the variable x_ij with j > i returns x_ji. Querying for x_ii (which does -// not exist) gives a CHECK failure. -// -// The Model object passed in to create EdgeVariables must outlive this. -class EdgeVariables { - public: - EdgeVariables(math_opt::Model& model, const int n) { - variables_.resize(n); - for (int i = 0; i < n; ++i) { - variables_[i].reserve(i); - for (int j = 0; j < i; ++j) { - variables_[i].push_back( - model.AddBinaryVariable(absl::StrCat("e_", i, "_", j))); - } - } - } - - math_opt::Variable get(const int i, const int j) const { - CHECK_NE(i, j); - return i > j ? variables_[i][j] : variables_[j][i]; - } - - int num_cities() const { return variables_.size(); } - - private: - std::vector> variables_; -}; - -// Produces a random TSP problem where cities have random locations that are -// I.I.D Uniform [0, 1]. -std::vector> RandomCities(int num_cities) { - absl::BitGen rand; - std::vector> cities; - for (int i = 0; i < num_cities; ++i) { - cities.push_back({absl::Uniform(rand, 0.0, 1.0), - absl::Uniform(rand, 0.0, 1.0)}); - } - return cities; -} - -std::vector> TestCities() { - return {{0, 0}, {0, 0.1}, {0.1, 0}, {0.1, 0.1}, - {1, 0}, {1, 0.1}, {0.9, 0}, {0.9, 0.1}}; -} - -// Given an n city TSP instance, computes the n by n distance matrix using the -// Euclidean distance. -std::vector> DistanceMatrix( - const std::vector>& cities) { - const int num_cities = cities.size(); - std::vector> distance_matrix( - num_cities, std::vector(num_cities, 0.0)); - for (int i = 0; i < num_cities; ++i) { - for (int j = 0; j < num_cities; ++j) { - if (i != j) { - const double dx = cities[i].first - cities[j].first; - const double dy = cities[i].second - cities[j].second; - distance_matrix[i][j] = std::sqrt(dx * dx + dy * dy); - } - } - } - return distance_matrix; -} - -// Given the EdgeVariables and a var_values containing the value of each edge in -// a solution, returns an n by n boolean matrix of which edges are used (with -// false diagonal elements). It is assumed that var_values are approximately 0-1 -// integer. -std::vector> EdgeValues( - const EdgeVariables& edge_vars, - const math_opt::VariableMap& var_values) { - const int n = edge_vars.num_cities(); - std::vector> edge_values(n, std::vector(n, false)); - for (int i = 0; i < n; ++i) { - for (int j = 0; j < n; ++j) { - if (i != j) { - edge_values[i][j] = var_values.at(edge_vars.get(i, j)) > 0.5; - } - } - } - return edge_values; -} - -// Given an n by n boolean matrix of edge values, returns a cycle decomposition. -// it is assumed that edge values respects the degree constraints (each row has -// only two true entries). Each cycle is represented as a list of cities with -// no repeats. -std::vector FindCycles( - const std::vector>& edge_values) { - // Algorithm: maintain a "visited" bit for each city indicating if we have - // formed a cycle containing this city. Consider the cities in order. When you - // find an unvisited city, start a new cycle beginning at this city. Then, - // build the cycle by finding an unvisited neighbor until no such neighbor - // exists (every city will have two neighbors, but eventually both will be - // visited). To find the "unvisited neighbor", we simply do a linear scan - // over the cities, checking both the adjacency matrix and the visited bit. - // - // Note that for this algorithm, in each cycle, the city with lowest index - // will be first, and the cycles will be sorted by their city of lowest index. - // This is an implementation detail and should not be relied upon. - const int n = edge_values.size(); - std::vector result; - std::vector visited(n, false); - for (int i = 0; i < n; ++i) { - if (visited[i]) { - continue; - } - std::vector cycle; - std::optional next = i; - while (next.has_value()) { - cycle.push_back(*next); - visited[*next] = true; - int current = *next; - next = std::nullopt; - // Scan for an unvisited neighbor. We can start at i+1 since we know that - // everything from i back is visited. - for (int j = i + 1; j < n; ++j) { - if (!visited[j] && edge_values[current][j]) { - next = j; - break; - } - } - } - result.push_back(cycle); - } - return result; -} - -// Given a cycle and an EdgeVariables, returns the cutset constraint for the set -// of nodes in cycle. -math_opt::BoundedLinearExpression CutsetConstraint( - const Cycle& cycle, const EdgeVariables& edge_vars) { - const int n = edge_vars.num_cities(); - const absl::flat_hash_set cycle_as_set(cycle.begin(), cycle.end()); - std::vector not_in_cycle; - for (int i = 0; i < n; ++i) { - if (!cycle_as_set.contains(i)) { - not_in_cycle.push_back(i); - } - } - math_opt::LinearExpression cutset_edges; - for (const int in_cycle : cycle) { - for (const int out_of_cycle : not_in_cycle) { - cutset_edges += edge_vars.get(in_cycle, out_of_cycle); - } - } - return cutset_edges >= 2; -} - -// Solves the TSP by returning the ordering of the cities that minimizes travel -// distance. -absl::StatusOr SolveTsp( - const std::vector>& cities) { - const int n = cities.size(); - const std::vector> distance_matrix = - DistanceMatrix(cities); - CHECK_GE(n, 3); - math_opt::Model model("tsp"); - const EdgeVariables edge_vars(model, n); - math_opt::LinearExpression edge_cost; - for (int i = 0; i < n; ++i) { - for (int j = i + 1; j < n; ++j) { - edge_cost += edge_vars.get(i, j) * distance_matrix[i][j]; - } - } - model.Minimize(edge_cost); - - // Add the degree constraints - for (int i = 0; i < n; ++i) { - math_opt::LinearExpression neighbors; - for (int j = 0; j < n; ++j) { - if (i != j) { - neighbors += edge_vars.get(i, j); - } - } - model.AddLinearConstraint(neighbors == 2, absl::StrCat("n_", i)); - } - math_opt::SolveArguments args; - args.parameters.enable_output = absl::GetFlag(FLAGS_solve_logs); - const int threads = absl::GetFlag(FLAGS_threads); - if (threads > 0) { - args.parameters.threads = threads; - } - args.callback_registration.events.insert( - math_opt::CallbackEvent::kMipSolution); - args.callback_registration.add_lazy_constraints = true; - args.callback = [&edge_vars](const math_opt::CallbackData& cb_data) { - // At event CallbackEvent::kMipSolution, a solution is always present. - CHECK(cb_data.solution.has_value()); - const std::vector cycles = - FindCycles(EdgeValues(edge_vars, *cb_data.solution)); - math_opt::CallbackResult result; - if (cycles.size() > 1) { - for (const Cycle& cycle : cycles) { - result.AddLazyConstraint(CutsetConstraint(cycle, edge_vars)); - } - } - return result; - }; - ASSIGN_OR_RETURN(const math_opt::SolveResult result, - math_opt::Solve(model, math_opt::SolverType::kGurobi, args)); - if (result.termination.reason != math_opt::TerminationReason::kOptimal) { - return util::InternalErrorBuilder() - << "Expected TSP solve terminate with reason optimal, found: " - << result.termination; - } - std::cout << "Route length: " << result.objective_value() << std::endl; - const std::vector cycles = - FindCycles(EdgeValues(edge_vars, result.variable_values())); - CHECK_EQ(cycles.size(), 1); - CHECK_EQ(cycles[0].size(), n); - return cycles[0]; -} - -// Produces an SVG to draw a route for a TSP. -std::string RouteSvg(const std::vector>& cities, - const Cycle& cycle) { - constexpr int image_px = 1000; - constexpr int r = 5; - constexpr int image_plus_border = image_px + 2 * r; - std::vector svg_lines; - svg_lines.push_back(absl::StrCat("")); - std::vector polygon_coords; - for (const int city : cycle) { - const int x = - static_cast(std::round(cities[city].first * image_px)) + r; - const int y = - static_cast(std::round(cities[city].second * image_px)) + r; - svg_lines.push_back(absl::StrCat("")); - polygon_coords.push_back(absl::StrCat(x, ",", y)); - } - std::string polygon_coords_string = absl::StrJoin(polygon_coords, " "); - svg_lines.push_back( - absl::StrCat("")); - svg_lines.push_back(""); - return absl::StrJoin(svg_lines, "\n"); -} - -void RealMain() { - std::vector> cities; - if (absl::GetFlag(FLAGS_test_instance)) { - cities = TestCities(); - } else { - cities = RandomCities(absl::GetFlag(FLAGS_num_cities)); - } - absl::StatusOr solution = SolveTsp(cities); - if (!solution.ok()) { - LOG(QFATAL) << solution.status(); - } - const std::string svg = RouteSvg(cities, *solution); - if (absl::GetFlag(FLAGS_output).empty()) { - std::cout << svg << std::endl; - } else { - QCHECK_OK( - file::SetContents(absl::GetFlag(FLAGS_output), svg, file::Defaults())); - } -} - -} // namespace - -int main(int argc, char** argv) { - InitGoogle(argv[0], &argc, &argv, true); - RealMain(); - return 0; -} diff --git a/ortools/math_opt/samples/tsp_mo.cc b/ortools/math_opt/samples/tsp_mo.cc index 0018ae2db00..d2ca1d9b322 100644 --- a/ortools/math_opt/samples/tsp_mo.cc +++ b/ortools/math_opt/samples/tsp_mo.cc @@ -15,23 +15,25 @@ // // In the Euclidean Traveling Salesperson Problem (TSP), you are given a list of // n cities, each with an (x, y) coordinate, and you must find an order to visit -// the cities in to minimize the (Euclidean) travel distance. +// the cities which minimizes the (Euclidean) travel distance. // // The MIP "cutset" formulation for the problem is as follows: // * Data: // n: An integer, the number of cities -// (x_i, y_i): a pair of floats for each i in 1..n, the location of each -// city -// d_ij for all (i, j) pairs of cities, the distance between city i and j. +// (x_i, y_i): a pair of floats for each i in N={0..n-1}, the location of +// each city +// d_ij for all (i, j) pairs of cities, the distance between city i and j +// (derived from the cities coordinates (x_i, y_i); this function is +// symmetric, i.e. d_ij = d_ji). // * Decision variables: // x_ij: A binary variable, indicates if the edge connecting i and j is // used. Note that x_ij == x_ji, because the problem is symmetric. We // only create variables for i < j, and have x_ji as an alias for // x_ij. // * MIP model: -// minimize sum_{i=1}^n sum_{j=1, j < i}^n d_ij * x_ij -// s.t. sum_{j=1, j != i}^n x_ij = 2 for all i = 1..n -// sum_{i in S} sum_{j not in S} x_ij >= 2 for all S subset {1,...,n} +// minimize sum_{i in N} sum_{j in N, j < i} d_ij * x_ij +// s.t. sum_{j in N, j != i} x_ij = 2 for all i in N +// sum_{i in S} sum_{j not in S} x_ij >= 2 for all S subset N // |S| >= 3, |S| <= n - 3 // x_ij in {0, 1} // The first set of constraints are called the degree constraints, and the @@ -51,15 +53,20 @@ // Note that this is a minimal TSP solution, more sophisticated MIP methods are // possible. +#include #include #include +#include +#include +#include +#include #include "absl/flags/flag.h" #include "absl/random/random.h" #include "absl/random/uniform_real_distribution.h" #include "absl/strings/str_cat.h" #include "absl/strings/str_join.h" -#include "ortools/base/file.h" +#include "ortools/base/helpers.h" #include "ortools/base/init_google.h" #include "ortools/base/logging.h" #include "ortools/base/status_builder.h" @@ -212,22 +219,21 @@ std::vector FindCycles( return result; } -// Given a cycle and an EdgeVariables, returns the cutset constraint for the set -// of nodes in cycle. +// Returns the cutset constraint for the given set of nodes. math_opt::BoundedLinearExpression CutsetConstraint( - const Cycle& cycle, const EdgeVariables& edge_vars) { + const std::vector& nodes, const EdgeVariables& edge_vars) { const int n = edge_vars.num_cities(); - const absl::flat_hash_set cycle_as_set(cycle.begin(), cycle.end()); - std::vector not_in_cycle; + const absl::flat_hash_set node_set(nodes.begin(), nodes.end()); + std::vector not_in_set; for (int i = 0; i < n; ++i) { - if (!cycle_as_set.contains(i)) { - not_in_cycle.push_back(i); + if (!node_set.contains(i)) { + not_in_set.push_back(i); } } math_opt::LinearExpression cutset_edges; - for (const int in_cycle : cycle) { - for (const int out_of_cycle : not_in_cycle) { - cutset_edges += edge_vars.get(in_cycle, out_of_cycle); + for (const int in_set : nodes) { + for (const int out_of_set : not_in_set) { + cutset_edges += edge_vars.get(in_set, out_of_set); } } return cutset_edges >= 2; @@ -285,11 +291,7 @@ absl::StatusOr SolveTsp( }; ASSIGN_OR_RETURN(const math_opt::SolveResult result, math_opt::Solve(model, math_opt::SolverType::kGurobi, args)); - if (result.termination.reason != math_opt::TerminationReason::kOptimal) { - return util::InternalErrorBuilder() - << "Expected TSP solve terminate with reason optimal, found: " - << result.termination; - } + RETURN_IF_ERROR(result.termination.IsOptimal()); std::cout << "Route length: " << result.objective_value() << std::endl; const std::vector cycles = FindCycles(EdgeValues(edge_vars, result.variable_values())); diff --git a/ortools/math_opt/solvers/BUILD.bazel b/ortools/math_opt/solvers/BUILD.bazel index 7b89ad20ede..b7d4380b2a1 100644 --- a/ortools/math_opt/solvers/BUILD.bazel +++ b/ortools/math_opt/solvers/BUILD.bazel @@ -134,6 +134,7 @@ cc_library( "//ortools/math_opt/solvers/gurobi:g_gurobi", "//ortools/math_opt/validators:callback_validator", "//ortools/port:proto_utils", + "@com_google_absl//absl/algorithm:container", "@com_google_absl//absl/container:flat_hash_map", "@com_google_absl//absl/container:flat_hash_set", "@com_google_absl//absl/log",