268 lines
10 KiB
C++
268 lines
10 KiB
C++
// Copyright 2010-2024 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 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:
|
|
// * 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.
|
|
// * 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 in C to produce.
|
|
//
|
|
// Leader problem formulation:
|
|
// 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 in P (dual value from leader's demand
|
|
// constraint)
|
|
//
|
|
// Worker decision variables:
|
|
// * y_i: How many copies of piece i in P should be in the configuration.
|
|
//
|
|
// Worker formulation
|
|
// 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 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".
|
|
//
|
|
// 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/log/check.h"
|
|
#include "absl/status/status.h"
|
|
#include "absl/status/statusor.h"
|
|
#include "absl/strings/str_cat.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_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_lengths;
|
|
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 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_lengths, const int board_length) {
|
|
int num_pieces = piece_prices.size();
|
|
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("piece_", i)));
|
|
}
|
|
model.Maximize(math_opt::InnerProduct(pieces, piece_prices));
|
|
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));
|
|
RETURN_IF_ERROR(solve_result.termination.EnsureIsOptimal());
|
|
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 = 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];
|
|
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 piece = config.pieces[i];
|
|
const int use = config.quantity[i];
|
|
if (use >= 1) {
|
|
model.set_coefficient(demand_met[piece], v, use);
|
|
}
|
|
}
|
|
configs.push_back({config, v});
|
|
};
|
|
|
|
// 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}});
|
|
}
|
|
|
|
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());
|
|
RETURN_IF_ERROR(solve_result.termination.EnsureIsOptimal())
|
|
<< " 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.
|
|
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_lengths,
|
|
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));
|
|
RETURN_IF_ERROR(solve_result.termination.EnsureIsOptimalOrFeasible())
|
|
<< " in final cutting stock MIP";
|
|
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_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 boards." << std::endl;
|
|
std::cout << "Total boards 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;
|
|
}
|