Solves the resource constrained shortest path problem with idol's branch-and-price algorithm on its Dantzig-Wolfe reformulation using the Boost graph Dijkstra algorithm for the subproblem.
Resource-constrained shortest path problem (RSPP)
Given a graph \(G = (V, A)\) where \(V\) is a set of nodes and \(A\) is a set of arcs, the RSPP can be modeled as
\[\begin{align}
\min_x \quad & \sum_{(i,j)\in A} c_{ij} x_{ij} \\
\text{s.t.}\quad
& \sum_{j:(i,j)\in A} x_{ij} - \sum_{j:(j,i)\in A} x_{ji}
=
\begin{cases}
1 & i = s, \\
-1 & i = t, \\
0 & \text{otherwise},
\end{cases}
\quad \text{for all } i \in V, \\
& \sum_{(i,j)\in A} r_{ij} x_{ij} \le C, \\
& x_{ij} \in \{0,1\}, \quad \text{for all } (i,j)\in A.
\end{align}
\]
Here, \(c_{ij} > 0\) denotes the travel cost of arc \((i,j)\in A\), \(r_{ij} > 0\) its capacity consumption while \(C\) denotes the maximum resource consumption of a path. The binary variable \(x_{ij}\) equals \(1\) if and only if arc \((i,j)\in A\) is part of the path. The goal is to minimize the overall cost of a path.
Branch-and-Price algorithm
Let \(\Omega\) be the set of all feasible paths connecting source \(s\) to sink \(t\) in graph \(G\) that satisfy the flow conservation constraints. For each path \(p \in \Omega\), we define
- \(c_p = \sum_{(i,j) \in p} c_{ij}\) the total cost of path \(p\),
- \(r_p = \sum_{(i,j) \in p} r_{ij}\) the total resource consumption of path \(p\).
We introduce a binary decision variable \(\lambda_p\) for each path \(p \in \Omega\), where \(\lambda_p = 1\) if path \(p\) is selected, and \(0\) otherwise. The relationship with the original arc variables is given by \(x_{ij} = \sum_{p \in \Omega} \delta_{ij}^p \lambda_p\), where \(\delta_{ij}^p = 1\) if arc \((i,j)\) belongs to path \(p\), and \(0\) otherwise.
The Master Problem
By replacing the original arc variables with path variables, the master problem is formulated as:
\[\begin{align}
\min_{\lambda} \quad & \sum_{p \in \Omega} c_p \lambda_p \\
\text{s.t.} \quad
& \sum_{p \in \Omega} r_p \lambda_p \le C, \quad && (\mu \le 0) \\
& \sum_{p \in \Omega} \lambda_p = 1, \quad && (\pi) \\
& \lambda_p \in \{0,1\} \quad \forall p \in \Omega,
\end{align}
\]
in which \(\mu\) is the dual variable associated with the resource constraint and \(\pi\) is the dual variable associated with the convexity constraint.
Pricing Subproblem
Since the set \(\Omega\) is exponentially large, a restricted master problem (RMP) is solved over a subset of paths \(\bar{\Omega} \subset \Omega\). A path with a negative reduced cost indicates that it can improve the current solution.
The reduced cost \(\bar{c}_p\) of a path \(p\) is expressed using the optimal dual variables \(\mu\) and \(\pi\) from the RMP as
\[\bar{c}_p = c_p - \mu r_p - \pi.
\]
Expressing \(c_p\) and \(r_p\) back in terms of the original arcs, we obtain
\[\bar{c}_p = \sum_{(i,j) \in p} (c_{ij} - \mu r_{ij}) - \pi.
\]
The pricing subproblem aims to find a path \(p \in \Omega\) that minimizes this reduced cost. Dropping the constant \(\pi\), the optimization subproblem becomes
\[\min_{p \in \Omega} \sum_{(i,j) \in p} (c_{ij} - \mu r_{ij}).
\]
Therefore, the pricing subproblem reduces to a standard shortest path problem on graph \(G\) from \(s\) to \(t\) with modified arc weights. In this example, the subproblem will be solved using the Dijkstra's algorithm from Boost Graph library.
Graph definition
We consider a simple instance defined on a directed graph \(G = (V, A)\).
digraph G {
rankdir=LR;
node [shape=circle];
s -> A [label="c=2, r=4"];
s -> B [label="c=3, r=2"];
A -> t [label="c=2, r=3"];
B -> t [label="c=2, r=2"];
A -> B [label="c=1, r=2"];
B -> A [label="c=1, r=2"];
}
The node set is given by
\[V = \{s, A, B, t\}
\]
and the arc set is
\[A = \{(s,A), (s,B), (A,t), (B,t), (A,B), (B,A)\}.
\]
Each arc \((i,j) \in A\) is associated with a cost \(c_{ij}\) and a resource consumption \(r_{ij}\):
| Arc \((i,j)\) | Cost \(c_{ij}\) | Resource \(r_{ij}\) |
| \((s,A)\) | 2 | 4 |
| \((s,B)\) | 3 | 2 |
| \((A,t)\) | 2 | 3 |
| \((B,t)\) | 2 | 2 |
| \((A,B)\) | 1 | 2 |
| \((B,A)\) | 1 | 2 |
The resource capacity is \(C = 6\).
Implementation in idol
#include <iostream>
#include "idol/mixed-integer/modeling/variables/Var.h"
#include "idol/mixed-integer/modeling/models/Model.h"
#include "idol/mixed-integer/modeling/objects/Env.h"
#include "idol/mixed-integer/modeling/expressions/operations/operators.h"
#include "idol/mixed-integer/optimizers/wrappers/GLPK/GLPK.h"
#include "idol/mixed-integer/optimizers/dantzig-wolfe/DantzigWolfeDecomposition.h"
#include "idol/mixed-integer/optimizers/branch-and-bound/BranchAndBound.h"
#include "idol/mixed-integer/optimizers/branch-and-bound/node-selection-rules/factories/BestBound.h"
#include "idol/mixed-integer/optimizers/branch-and-bound/branching-rules/factories/MostInfeasible.h"
#include "idol/mixed-integer/optimizers/dantzig-wolfe/infeasibility-strategies/FarkasPricing.h"
#include "idol/general/optimizers/LambdaOptimizer/LambdaOptimizer.h"
#include "idol/general/optimizers/LambdaOptimizer/LambdaContext.h"
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/dijkstra_shortest_paths.hpp>
using namespace idol;
int main() {
using Graph = boost::adjacency_list<
boost::vecS,
boost::vecS,
boost::directedS,
boost::no_property,
boost::property<boost::edge_weight_t, double>
>;
using Vertex = boost::graph_traits<Graph>::vertex_descriptor;
using EdgeKey = std::pair<Vertex, Vertex>;
Graph g;
const auto& s = add_vertex(g);
const auto& A = add_vertex(g);
const auto& B = add_vertex(g);
const auto& t = add_vertex(g);
const std::vector<char> vertex_names = {'s', 'A', 'B', 't'};
add_edge(s, A, 2.0, g);
add_edge(s, B, 3.0, g);
add_edge(A, t, 2.0, g);
add_edge(B, t, 2.0, g);
add_edge(A, B, 1.0, g);
add_edge(B, A, 1.0, g);
const auto& x = model.add_vars(
Dim<1>(boost::num_edges(g)), 0, 1, idol::VarType::Binary, 0.,
"x");
std::map<EdgeKey, Var> edge_to_var;
auto [it_edges, it_edges_end] = boost::edges(g);
for (unsigned int i = 0; it_edges != it_edges_end; ++it_edges, ++i) {
edge_to_var.emplace(EdgeKey{boost::source(*it_edges, g), boost::target(*it_edges, g)}, x[i]);
}
std::map<Var, EdgeKey> var_to_edge;
for (const auto& [edge, var] : edge_to_var) {
var_to_edge.emplace(var, edge);
}
for(const auto& v: boost::make_iterator_range(boost::vertices(g))){
for (auto e : boost::make_iterator_range(boost::edges(g))) {
if (boost::target(e,g) == v) {
expr -= edge_to_var.at({boost::source(e,g), boost::target(e,g)});
}
if (boost::source(e,g) == v) {
expr += edge_to_var.at({boost::source(e,g), boost::target(e,g)});
}
}
if(v == s){
model.add_ctr(expr == 1, "flow_out_" + std::to_string(v));
}
else if(v == t){
model.add_ctr(expr == -1, "flow_in_" + std::to_string(v));
}
else{
model.add_ctr(expr == 0, "flow_null_" + std::to_string(v));
}
}
const std::vector<double> resources{4., 2., 3., 2., 2., 2.};
const double capacity = 6.;
model.add_ctr(idol_Sum(k,
Range(boost::num_edges(g)), resources[k] * x[k]) <= capacity,
"resource");
auto weight_map = get(boost::edge_weight, g);
for (const auto& e : boost::make_iterator_range(edges(g))) {
obj_expr += weight_map[e] * edge_to_var.at(EdgeKey{boost::source(e, g), boost::target(e, g)});
}
model.set_obj_expr(obj_expr);
const Annotation decomposition(env,
"decomposition", MasterId);
for(const auto& ctr : model.ctrs()){
if(ctr.name().rfind("flow_", 0) == 0){
ctr.set(decomposition, 0);
}
}
const auto lambda = [&g, &s, &t, &edge_to_var](
LambdaContext& t_ctx){
const auto& obj_func = t_ctx.get_model().get_obj_expr().affine().linear();
auto weight_map = get(boost::edge_weight, g);
for (const auto& edge : boost::make_iterator_range(boost::edges(g))) {
const auto& var = edge_to_var.at(EdgeKey{boost::source(edge, g), boost::target(edge, g)});
boost::put(boost::edge_weight, g, edge, obj_func.get(var));
}
std::vector<double> dist(num_vertices(g));
std::vector<Vertex> pred(num_vertices(g));
dijkstra_shortest_paths(
g,
s,
boost::weight_map(weight_map)
.predecessor_map(&pred[0])
.distance_map(&dist[0])
);
std::vector<Vertex> path;
for (Vertex v = t; v != s; v = pred[v]) {
path.push_back(v);
}
path.push_back(s);
std::ranges::reverse(path);
for (unsigned int i = 0; i+1 < path.size(); ++i) {
const auto& [e, found] = boost::edge(path[i], path[i+1], g);
assert(found);
auto const& var = edge_to_var.at(EdgeKey{boost::source(e, g), boost::target(e, g)});
t_ctx.set_var_primal(var, 1.);
}
t_ctx.set_status(Optimal);
const auto best_obj = dist[t] + t_ctx.get_model().get_obj_expr().affine().constant();
t_ctx.set_best_obj(best_obj);
t_ctx.set_best_bound(best_obj);
};
column_generation.with_master_optimizer(GLPK::ContinuousRelaxation());
column_generation.with_default_sub_problem_spec(subproblem_specifications);
column_generation.with_hard_branching(false);
column_generation.with_logs(true);
branch_and_bound.with_node_selection_rule(
BestBound());
branch_and_bound.with_logs(true);
const auto branch_and_price = branch_and_bound + column_generation;
model.use(branch_and_price);
model.optimize();
const auto status = model.get_status();
const auto reason = model.get_reason();
std::cout << "Solution status: " << status << std::endl;
std::cout << "Reason: " << reason << std::endl;
if (status == Optimal || status == Feasible) {
const auto n_solutions = model.get_n_solutions();
std::cout << "Number of solutions: " << n_solutions << std::endl;
for (unsigned int i = 0 ; i < n_solutions ; ++i) {
model.set_solution_index(i);
std::cout << "Solution " << i << std::endl;
const auto& primals = save_primal(model) ;
for (const auto& [var, value] : primals) {
const auto& edgeKey = var_to_edge[var];
std::cout << vertex_names[edgeKey.first] << " -> " << vertex_names[edgeKey.second] << std::endl;
}
}
}
return 0;
}
1.16.1