always-cautious/ortools-clone
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Export math_opt from google3

Browse Source
This commit is contained in:
Corentin Le Molgat
2023-07-13 20:22:58 +02:00
parent 944fd974a5
commit a29d1eed56
29 changed files with 415 additions and 2707 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
ortools/math_opt/BUILD.bazel
View File

@@ -110,6 +110,7 @@ proto_library(
deps = [
":solution_proto",
"//ortools/gscip:gscip_proto",
"//ortools/pdlp:solve_log_proto",
"@com_google_protobuf//:duration_proto",
],
)

4
ortools/math_opt/cpp/BUILD.bazel
View File

@@ -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",

1
ortools/math_opt/cpp/solution.h
View File

@@ -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"

42
ortools/math_opt/cpp/solve_result.cc
View File

@@ -20,13 +20,12 @@
#include <utility>
#include <vector>
#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<TerminationReason> 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> Termination::FromProto(
const TerminationProto& termination_proto) {
const std::optional<TerminationReason> reason =
@@ -327,6 +357,9 @@ absl::StatusOr<SolveResultProto> 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> 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;
}

71
ortools/math_opt/cpp/solve_result.h
View File

@@ -17,6 +17,7 @@
#ifndef OR_TOOLS_MATH_OPT_CPP_SOLVE_RESULT_H_
#define OR_TOOLS_MATH_OPT_CPP_SOLVE_RESULT_H_
#include <initializer_list>
#include <optional>
#include <ostream>
#include <string>
@@ -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<TerminationReason> 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.
//

49
ortools/math_opt/result.proto
View File

@@ -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.

48
ortools/math_opt/samples/BUILD.bazel
View File

@@ -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",

5
ortools/math_opt/samples/advanced_linear_programming.cc → ortools/math_opt/samples/advanced_linear_programming_mo.cc
View File

@@ -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;

88
ortools/math_opt/samples/area_socp_mo.cc Normal file
View File

@@ -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 <cmath>
#include <iostream>
#include <limits>
#include <ostream>
#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<double>::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;
}

72
ortools/math_opt/samples/basic_example.cc
View File

@@ -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 <iostream>
#include <limits>
#include <ostream>
#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;
}

17
ortools/math_opt/samples/basic_example_mo.cc
View File

@@ -15,6 +15,7 @@
#include <iostream>
#include <limits>
#include <ostream>
#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

380
ortools/math_opt/samples/cocktail_hour.cc
View File

@@ -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 <iostream>
#include <limits>
#include <ostream>
#include <string>
#include <vector>
#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<std::string>, existing_ingredients, {},
"Ingredients you already have (ignored in analysis mode).");
ABSL_FLAG(std::vector<std::string>, unavailable_ingredients, {},
"Ingredients you cannot get (ignored in analysis mode).");
ABSL_FLAG(std::vector<std::string>, required_cocktails, {},
"Cocktails you must be able to make (ignored in analysis mode).");
ABSL_FLAG(std::vector<std::string>, 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<std::string> ingredients;
};
std::vector<Cocktail> 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<std::string> ingredients;
std::vector<Cocktail> cocktails;
};
absl::StatusOr<Menu> SolveForMenu(
const int max_new_ingredients, const bool enable_solver_output,
const absl::flat_hash_set<std::string>& existing_ingredients,
const absl::flat_hash_set<std::string>& unavailable_ingredients,
const absl::flat_hash_set<std::string>& required_cocktails,
const absl::flat_hash_set<std::string>& blocked_cocktails) {
const std::vector<Cocktail> all_cocktails = AllCocktails();
math_opt::Model model("Cocktail hour");
absl::flat_hash_map<std::string, math_opt::Variable> 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<std::string, math_opt::Variable> 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<std::string> SetFromVec(
const std::vector<std::string>& 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> 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<Cocktail>& cocktails,
absl::string_view title = "Cocktail Hour") {
std::vector<std::string> 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<std::string>({"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;
}

14
ortools/math_opt/samples/cocktail_hour_mo.cc
View File

@@ -33,6 +33,9 @@
// will create a PDF in /tmp.
#include <iostream>
#include <limits>
#include <ostream>
#include <string>
#include <vector>
#include "absl/container/flat_hash_set.h"
#include "absl/status/status.h"
@@ -252,14 +255,7 @@ absl::StatusOr<Menu> 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<Cocktail>& cocktails,
const std::string& title = "Cocktail Hour") {
absl::string_view title = "Cocktail Hour") {
std::vector<std::string> lines;
lines.push_back("\\documentclass{article}");
lines.push_back("\\usepackage{fullpage}");

277
ortools/math_opt/samples/cutting_stock.cc
View File

@@ -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 <cmath>
#include <iostream>
#include <limits>
#include <ostream>
#include <utility>
#include <vector>
#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<double>::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<int> piece_sizes;
std::vector<int> 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<int> pieces;
std::vector<int> 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<Configuration> configurations;
std::vector<int> 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<std::pair<Configuration, double>> BestConfiguration(
const std::vector<double>& piece_prices,
const std::vector<int>& 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<math_opt::Variable> 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<int>(
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<CuttingStockSolution> SolveCuttingStock(
const CuttingStockInstance& instance) {
math_opt::Model model("cutting_stock");
model.set_minimize();
const int n = instance.piece_sizes.size();
std::vector<math_opt::LinearConstraint> 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<std::pair<Configuration, math_opt::Variable>> 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<double> 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<int>(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;
}

124
ortools/math_opt/samples/cutting_stock_mo.cc
View File

@@ -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 <cmath>
#include <iostream>
#include <limits>
#include <ostream>
#include <utility>
#include <vector>
@@ -86,11 +94,11 @@ namespace {
namespace math_opt = operations_research::math_opt;
constexpr double kInf = std::numeric_limits<double>::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<int> piece_sizes;
std::vector<int> piece_lengths;
std::vector<int> 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<std::pair<Configuration, double>> BestConfiguration(
const std::vector<double>& piece_prices,
const std::vector<int>& piece_sizes, const int board_size) {
const std::vector<int>& 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<math_opt::Variable> 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<int>(
@@ -158,7 +161,7 @@ absl::StatusOr<CuttingStockSolution> 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<int>(instance.piece_lengths.size());
std::vector<math_opt::LinearConstraint> demand_met;
for (int i = 0; i < n; ++i) {
const int d = instance.piece_demands[i];
@@ -169,17 +172,17 @@ absl::StatusOr<CuttingStockSolution> 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<CuttingStockSolution> 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<CuttingStockSolution> 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<CuttingStockSolution> 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;
}

683
ortools/math_opt/samples/facility_lp_benders.cc
View File

@@ -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 <algorithm>
#include <iostream>
#include <limits>
#include <memory>
#include <ostream>
#include <utility>
#include <vector>
#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<double>::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<int, int>;
// 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<Edge>& edges() const { return edges_; }
const std::vector<Edge>& edges_incident_to_facility(
const int facility) const {
return facility_edge_incidence_[facility];
}
const std::vector<Edge>& 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<Edge> edges_;
absl::flat_hash_map<Edge, double> edge_costs_;
std::vector<std::vector<Edge>> facility_edge_incidence_;
std::vector<std::vector<Edge>> 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<std::vector<Edge>>(num_facilities);
location_edge_incidence_ =
std::vector<std::vector<Edge>>(num_locations, std::vector<Edge>());
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<int>(0, num_locations - 1)(bitgen);
const std::pair<int, int> 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<int>(0, num_facilities - 1)(bitgen);
const std::pair<int, int> 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<math_opt::Variable> z;
for (int j = 0; j < num_facilities; j++) {
z.push_back(model.AddContinuousVariable(0.0, kInf));
}
// Flow variables
absl::flat_hash_map<Edge, math_opt::Variable> 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<math_opt::Variable> 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<double> 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<std::unique_ptr<SecondStageSolver>> New(
FacilityLocationInstance instance, math_opt::SolverType solver_type);
absl::StatusOr<std::pair<double, Cut>> Solve(
const std::vector<double>& z_values, double w_value,
double fist_stage_objective);
private:
SecondStageSolver(FacilityLocationInstance instance,
math_opt::SolveParameters second_stage_params);
absl::StatusOr<Cut> OptimalityCut(
const math_opt::SolveResult& second_stage_result);
absl::StatusOr<Cut> 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<Edge, math_opt::Variable> x_;
std::vector<math_opt::LinearConstraint> supply_constraints_;
std::vector<math_opt::LinearConstraint> demand_constraints_;
std::unique_ptr<math_opt::IncrementalSolver> solver_;
};
absl::StatusOr<math_opt::SolveParameters> 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<std::unique_ptr<SecondStageSolver>> 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<SecondStageSolver> second_stage_solver =
absl::WrapUnique<SecondStageSolver>(
new SecondStageSolver(std::move(instance), parameters));
ASSIGN_OR_RETURN(std::unique_ptr<math_opt::IncrementalSolver> 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<std::pair<double, Cut>> SecondStageSolver::Solve(
const std::vector<double>& 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<Cut> 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<Cut> 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<math_opt::IncrementalSolver> first_stage_solver,
math_opt::IncrementalSolver::New(&first_stage.model, solver_type));
// Setup second stage solver.
ASSIGN_OR_RETURN(std::unique_ptr<SecondStageSolver> second_stage_solver,
SecondStageSolver::New(instance, solver_type));
// Start Benders
int iteration = 0;
double best_upper_bound = kInf;
std::vector<double> 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;
}

44
ortools/math_opt/samples/facility_lp_benders_mo.cc
View File

@@ -49,6 +49,8 @@
#include <algorithm>
#include <iostream>
#include <limits>
#include <memory>
#include <ostream>
#include <utility>
#include <vector>
@@ -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<std::vector<std::pair<int, int>>>(num_facilities);
location_edge_incidence_ = std::vector<std::vector<std::pair<int, int>>>(
num_locations, std::vector<std::pair<int, int>>());
facility_edge_incidence_ = std::vector<std::vector<Edge>>(num_facilities);
location_edge_incidence_ =
std::vector<std::vector<Edge>>(num_locations, std::vector<Edge>());
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<int, int> 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<math_opt::Variable> 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<math_opt::SolveParameters> 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<Cut> 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<Cut> 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]);
}

9
ortools/math_opt/samples/graph_coloring.cc → ortools/math_opt/samples/graph_coloring_mo.cc
View File

@@ -110,14 +110,7 @@ absl::StatusOr<GraphColoringSolution> 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;

85
ortools/math_opt/samples/integer_programming.cc
View File

@@ -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 <cmath>
#include <iostream>
#include <limits>
#include <ostream>
#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<double>::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;
}

25
ortools/math_opt/samples/integer_programming_mo.cc
View File

@@ -13,8 +13,10 @@
// Simple integer programming example
#include <cmath>
#include <iostream>
#include <limits>
#include <ostream>
#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

493
ortools/math_opt/samples/lagrangian_relaxation.cc
View File

@@ -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 <math.h>
#include <algorithm>
#include <iostream>
#include <limits>
#include <memory>
#include <ostream>
#include <string>
#include <vector>
#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<Arc> arcs;
int source;
int sink;
};
struct FlowModel {
FlowModel() : model(std::make_unique<math_opt::Model>("LagrangianProblem")) {}
std::unique_ptr<math_opt::Model> model;
math_opt::LinearExpression cost;
math_opt::LinearExpression resource_1;
math_opt::LinearExpression resource_2;
std::vector<math_opt::Variable> 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<math_opt::LinearExpression> out_flow(graph.num_nodes);
std::vector<math_opt::LinearExpression> 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<Arc> 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<FlowModel> 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<math_opt::LinearExpression> out_flow(graph.num_nodes);
std::vector<math_opt::LinearExpression> 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<double> mu;
std::vector<math_opt::LinearExpression> 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<double>::infinity();
double lower_bound = -std::numeric_limits<double>::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<double>& 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<double> 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;
}

75
ortools/math_opt/samples/lagrangian_relaxation_mo.cc
View File

@@ -82,6 +82,7 @@
#include <iostream>
#include <limits>
#include <memory>
#include <ostream>
#include <string>
#include <vector>
@@ -239,25 +240,19 @@ absl::StatusOr<FlowModel> 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<double>& vars_val = result.variable_values();
bool feasible = true;

93
ortools/math_opt/samples/linear_programming.cc
View File

@@ -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 <iostream>
#include <limits>
#include <ostream>
#include <string>
#include <vector>
#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<double>::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<math_opt::Variable> x;
for (int j = 0; j < 3; j++) {
x.push_back(model.AddContinuousVariable(0.0, kInf, absl::StrCat("x", j)));
}
// Constraints
std::vector<math_opt::LinearConstraint> 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;
}

8
ortools/math_opt/samples/linear_programming_mo.cc
View File

@@ -15,6 +15,7 @@
#include <iostream>
#include <limits>
#include <ostream>
#include <string>
#include <vector>
@@ -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();

6
ortools/math_opt/samples/linear_regression.cc → ortools/math_opt/samples/linear_regression_mo.cc
View File

@@ -163,11 +163,7 @@ absl::StatusOr<LinearModel> 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)};
}

0
ortools/math_opt/samples/time_indexed_scheduling.cc → ortools/math_opt/samples/time_indexed_scheduling_mo.cc
View File

359
ortools/math_opt/samples/tsp.cc
View File

@@ -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 <cmath>
#include <iostream>
#include <optional>
#include <ostream>
#include <string>
#include <utility>
#include <vector>
#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<int>;
// 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<std::vector<math_opt::Variable>> variables_;
};
// Produces a random TSP problem where cities have random locations that are
// I.I.D Uniform [0, 1].
std::vector<std::pair<double, double>> RandomCities(int num_cities) {
absl::BitGen rand;
std::vector<std::pair<double, double>> cities;
for (int i = 0; i < num_cities; ++i) {
cities.push_back({absl::Uniform<double>(rand, 0.0, 1.0),
absl::Uniform<double>(rand, 0.0, 1.0)});
}
return cities;
}
std::vector<std::pair<double, double>> 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<std::vector<double>> DistanceMatrix(
const std::vector<std::pair<double, double>>& cities) {
const int num_cities = cities.size();
std::vector<std::vector<double>> distance_matrix(
num_cities, std::vector<double>(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<std::vector<bool>> EdgeValues(
const EdgeVariables& edge_vars,
const math_opt::VariableMap<double>& var_values) {
const int n = edge_vars.num_cities();
std::vector<std::vector<bool>> edge_values(n, std::vector<bool>(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<Cycle> FindCycles(
const std::vector<std::vector<bool>>& 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<Cycle> result;
std::vector<bool> visited(n, false);
for (int i = 0; i < n; ++i) {
if (visited[i]) {
continue;
}
std::vector<int> cycle;
std::optional<int> 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<int> cycle_as_set(cycle.begin(), cycle.end());
std::vector<int> 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<Cycle> SolveTsp(
const std::vector<std::pair<double, double>>& cities) {
const int n = cities.size();
const std::vector<std::vector<double>> 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<Cycle> 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<Cycle> 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<std::pair<double, double>>& cities,
const Cycle& cycle) {
constexpr int image_px = 1000;
constexpr int r = 5;
constexpr int image_plus_border = image_px + 2 * r;
std::vector<std::string> svg_lines;
svg_lines.push_back(absl::StrCat("<svg width=\"", image_plus_border,
"\" height=\"", image_plus_border, "\">"));
std::vector<std::string> polygon_coords;
for (const int city : cycle) {
const int x =
static_cast<int>(std::round(cities[city].first * image_px)) + r;
const int y =
static_cast<int>(std::round(cities[city].second * image_px)) + r;
svg_lines.push_back(absl::StrCat("<circle cx=\"", x, "\" cy=\"", y,
"\" r=\"", r, "\" fill=\"blue\" />"));
polygon_coords.push_back(absl::StrCat(x, ",", y));
}
std::string polygon_coords_string = absl::StrJoin(polygon_coords, " ");
svg_lines.push_back(
absl::StrCat("<polygon fill=\"none\" stroke=\"blue\" points=\"",
polygon_coords_string, "\" />"));
svg_lines.push_back("</svg>");
return absl::StrJoin(svg_lines, "\n");
}
void RealMain() {
std::vector<std::pair<double, double>> cities;
if (absl::GetFlag(FLAGS_test_instance)) {
cities = TestCities();
} else {
cities = RandomCities(absl::GetFlag(FLAGS_num_cities));
}
absl::StatusOr<Cycle> 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;
}

48
ortools/math_opt/samples/tsp_mo.cc
View File

@@ -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 <cmath>
#include <iostream>
#include <optional>
#include <ostream>
#include <string>
#include <utility>
#include <vector>
#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<Cycle> 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<int>& nodes, const EdgeVariables& edge_vars) {
const int n = edge_vars.num_cities();
const absl::flat_hash_set<int> cycle_as_set(cycle.begin(), cycle.end());
std::vector<int> not_in_cycle;
const absl::flat_hash_set<int> node_set(nodes.begin(), nodes.end());
std::vector<int> 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<Cycle> 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<Cycle> cycles =
FindCycles(EdgeValues(edge_vars, result.variable_values()));

1
ortools/math_opt/solvers/BUILD.bazel
View File

@@ -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",