# Changeset 710:8b0df68370a4 in lemon

Ignore:
Timestamp:
05/12/09 12:06:40 (12 years ago)
Branch:
default
Phase:
public
Message:

Fix the GEQ/LEQ handling in NetworkSimplex? + improve doc (#291)

• Fix the optimality conditions for the GEQ/LEQ form.
• Fix the initialization of the algortihm. It ensures correct solutions and it is much faster for the inequality forms.
• Fix the pivot rules to search all the arcs that have to be allowed to get in the basis.
• Better block size for the Block Search pivot rule.
• Improve documentation of the problem and move it to a separate page.
Files:
1 added
3 edited

Unmodified
Added
Removed
• ## doc/Makefile.am

 r634 doc/mainpage.dox \ doc/migration.dox \ doc/min_cost_flow.dox \ doc/named-param.dox \ doc/namespaces.dox \

 r707 /** @defgroup min_cost_flow Minimum Cost Flow Algorithms @defgroup min_cost_flow_algs Minimum Cost Flow Algorithms @ingroup algs This group contains the algorithms for finding minimum cost flows and circulations. The \e minimum \e cost \e flow \e problem is to find a feasible flow of minimum total cost from a set of supply nodes to a set of demand nodes in a network with capacity constraints (lower and upper bounds) and arc costs. Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$, \f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and upper bounds for the flow values on the arcs, for which \f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$, \f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the signed supply values of the nodes. If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$ supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with \f$-sup(u)\f$ demand. A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution of the following optimization problem. \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq sup(u) \quad \forall u\in V \f] \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or negative in order to have a feasible solution (since the sum of the expressions on the left-hand side of the inequalities is zero). It means that the total demand must be greater or equal to the total supply and all the supplies have to be carried out from the supply nodes, but there could be demands that are not satisfied. If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand constraints have to be satisfied with equality, i.e. all demands have to be satisfied and all supplies have to be used. If you need the opposite inequalities in the supply/demand constraints (i.e. the total demand is less than the total supply and all the demands have to be satisfied while there could be supplies that are not used), then you could easily transform the problem to the above form by reversing the direction of the arcs and taking the negative of the supply values (e.g. using \ref ReverseDigraph and \ref NegMap adaptors). However \ref NetworkSimplex algorithm also supports this form directly for the sake of convenience. A feasible solution for this problem can be found using \ref Circulation. Note that the above formulation is actually more general than the usual definition of the minimum cost flow problem, in which strict equalities are required in the supply/demand contraints, i.e. \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) = sup(u) \quad \forall u\in V. \f] However if the sum of the supply values is zero, then these two problems are equivalent. So if you need the equality form, you have to ensure this additional contraint for the algorithms. The dual solution of the minimum cost flow problem is represented by node potentials \f$\pi: V\rightarrow\mathbf{Z}\f$. An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$ node potentials the following \e complementary \e slackness optimality conditions hold. - For all \f$uv\in A\f$ arcs: - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$; - if \f$lower(uv) • ## lemon/network_simplex.h  r690 #define LEMON_NETWORK_SIMPLEX_H /// \ingroup min_cost_flow /// \ingroup min_cost_flow_algs /// /// \file namespace lemon { /// \addtogroup min_cost_flow /// \addtogroup min_cost_flow_algs /// @{ /// constraints of the \ref min_cost_flow "minimum cost flow problem". /// /// The default supply type is \c GEQ, since this form is supported /// by other minimum cost flow algorithms and the \ref Circulation /// algorithm, as well. /// The \c LEQ problem type can be selected using the \ref supplyType() /// function. /// /// Note that the equality form is a special case of both supply types. /// The default supply type is \c GEQ, the \c LEQ type can be /// selected using \ref supplyType(). /// The equality form is a special case of both supply types. enum SupplyType { /// This option means that there are "greater or equal" /// supply/demand constraints in the definition, i.e. the exact /// formulation of the problem is the following. /** \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq sup(u) \quad \forall u\in V \f] \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] */ /// It means that the total demand must be greater or equal to the /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$must be zero or /// negative) and all the supplies have to be carried out from /// the supply nodes, but there could be demands that are not /// satisfied. /// supply/demand constraints in the definition of the problem. GEQ, /// It is just an alias for the \c GEQ option. CARRY_SUPPLIES = GEQ, /// This option means that there are "less or equal" /// supply/demand constraints in the definition, i.e. the exact /// formulation of the problem is the following. /** \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq sup(u) \quad \forall u\in V \f] \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] */ /// It means that the total demand must be less or equal to the /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f\$ must be zero or /// positive) and all the demands have to be satisfied, but there /// could be supplies that are not carried out from the supply /// nodes. LEQ, /// It is just an alias for the \c LEQ option. SATISFY_DEMANDS = LEQ /// supply/demand constraints in the definition of the problem. LEQ }; int _node_num; int _arc_num; int _all_arc_num; int _search_arc_num; // Parameters of the problem const CostVector &_pi; int &_in_arc; int _arc_num; int _search_arc_num; // Pivot rule data _source(ns._source), _target(ns._target), _cost(ns._cost), _state(ns._state), _pi(ns._pi), _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), _next_arc(0) {} bool findEnteringArc() { Cost c; for (int e = _next_arc; e < _arc_num; ++e) { for (int e = _next_arc; e < _search_arc_num; ++e) { c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); if (c < 0) { const CostVector &_pi; int &_in_arc; int _arc_num; int _search_arc_num; public: _source(ns._source), _target(ns._target), _cost(ns._cost), _state(ns._state), _pi(ns._pi), _in_arc(ns.in_arc), _arc_num(ns._arc_num) _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num) {} bool findEnteringArc() { Cost c, min = 0; for (int e = 0; e < _arc_num; ++e) { for (int e = 0; e < _search_arc_num; ++e) { c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); if (c < min) { const CostVector &_pi; int &_in_arc; int _arc_num; int _search_arc_num; // Pivot rule data _source(ns._source), _target(ns._target), _cost(ns._cost), _state(ns._state), _pi(ns._pi), _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), _next_arc(0) { // The main parameters of the pivot rule const double BLOCK_SIZE_FACTOR = 2.0; const double BLOCK_SIZE_FACTOR = 0.5; const int MIN_BLOCK_SIZE = 10; _block_size = std::max( int(BLOCK_SIZE_FACTOR * std::sqrt(double(_arc_num))), std::sqrt(double(_search_arc_num))), MIN_BLOCK_SIZE ); } int cnt = _block_size; int e, min_arc = _next_arc; for (e = _next_arc; e < _arc_num; ++e) { for (e = _next_arc; e < _search_arc_num; ++e) { c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); if (c < min) { const CostVector &_pi; int &_in_arc; int _arc_num; int _search_arc_num; // Pivot rule data _source(ns._source), _target(ns._target), _cost(ns._cost), _state(ns._state), _pi(ns._pi), _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), _next_arc(0) { // The main parameters of the pivot rule _list_length = std::max( int(LIST_LENGTH_FACTOR * std::sqrt(double(_arc_num))), std::sqrt(double(_search_arc_num))), MIN_LIST_LENGTH ); _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length), min = 0; _curr_length = 0; for (e = _next_arc; e < _arc_num; ++e) { for (e = _next_arc; e < _search_arc_num; ++e) { c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); if (c < 0) { const CostVector &_pi; int &_in_arc; int _arc_num; int _search_arc_num; // Pivot rule data _source(ns._source), _target(ns._target), _cost(ns._cost), _state(ns._state), _pi(ns._pi), _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost) _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost) { // The main parameters of the pivot rule _block_size = std::max( int(BLOCK_SIZE_FACTOR * std::sqrt(double(_arc_num))), std::sqrt(double(_search_arc_num))), MIN_BLOCK_SIZE ); _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size), int limit = _head_length; for (int e = _next_arc; e < _arc_num; ++e) { for (int e = _next_arc; e < _search_arc_num; ++e) { _cand_cost[e] = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); _arc_num = countArcs(_graph); int all_node_num = _node_num + 1; int all_arc_num = _arc_num + _node_num; _source.resize(all_arc_num); _target.resize(all_arc_num); _lower.resize(all_arc_num); _upper.resize(all_arc_num); _cap.resize(all_arc_num); _cost.resize(all_arc_num); int max_arc_num = _arc_num + 2 * _node_num; _source.resize(max_arc_num); _target.resize(max_arc_num); _lower.resize(_arc_num); _upper.resize(_arc_num); _cap.resize(max_arc_num); _cost.resize(max_arc_num); _supply.resize(all_node_num); _flow.resize(all_arc_num); _flow.resize(max_arc_num); _pi.resize(all_node_num); _succ_num.resize(all_node_num); _last_succ.resize(all_node_num); _state.resize(all_arc_num); _state.resize(max_arc_num); // Copy the graph (store the arcs in a mixed order) Cost ART_COST; if (std::numeric_limits::is_exact) { ART_COST = std::numeric_limits::max() / 4 + 1; ART_COST = std::numeric_limits::max() / 2 + 1; } else { ART_COST = std::numeric_limits::min(); _last_succ[_root] = _root - 1; _supply[_root] = -_sum_supply; _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST; _pi[_root] = 0; // Add artificial arcs and initialize the spanning tree data structure for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { _parent[u] = _root; _pred[u] = e; _thread[u] = u + 1; _rev_thread[u + 1] = u; _succ_num[u] = 1; _last_succ[u] = u; _cost[e] = ART_COST; _cap[e] = INF; _state[e] = STATE_TREE; if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) { _flow[e] = _supply[u]; _forward[u] = true; _pi[u] = -ART_COST + _pi[_root]; } else { _flow[e] = -_supply[u]; _forward[u] = false; _pi[u] = ART_COST + _pi[_root]; } if (_sum_supply == 0) { // EQ supply constraints _search_arc_num = _arc_num; _all_arc_num = _arc_num + _node_num; for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { _parent[u] = _root; _pred[u] = e; _thread[u] = u + 1; _rev_thread[u + 1] = u; _succ_num[u] = 1; _last_succ[u] = u; _cap[e] = INF; _state[e] = STATE_TREE; if (_supply[u] >= 0) { _forward[u] = true; _pi[u] = 0; _source[e] = u; _target[e] = _root; _flow[e] = _supply[u]; _cost[e] = 0; } else { _forward[u] = false; _pi[u] = ART_COST; _source[e] = _root; _target[e] = u; _flow[e] = -_supply[u]; _cost[e] = ART_COST; } } } else if (_sum_supply > 0) { // LEQ supply constraints _search_arc_num = _arc_num + _node_num; int f = _arc_num + _node_num; for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { _parent[u] = _root; _thread[u] = u + 1; _rev_thread[u + 1] = u; _succ_num[u] = 1; _last_succ[u] = u; if (_supply[u] >= 0) { _forward[u] = true; _pi[u] = 0; _pred[u] = e; _source[e] = u; _target[e] = _root; _cap[e] = INF; _flow[e] = _supply[u]; _cost[e] = 0; _state[e] = STATE_TREE; } else { _forward[u] = false; _pi[u] = ART_COST; _pred[u] = f; _source[f] = _root; _target[f] = u; _cap[f] = INF; _flow[f] = -_supply[u]; _cost[f] = ART_COST; _state[f] = STATE_TREE; _source[e] = u; _target[e] = _root; _cap[e] = INF; _flow[e] = 0; _cost[e] = 0; _state[e] = STATE_LOWER; ++f; } } _all_arc_num = f; } else { // GEQ supply constraints _search_arc_num = _arc_num + _node_num; int f = _arc_num + _node_num; for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { _parent[u] = _root; _thread[u] = u + 1; _rev_thread[u + 1] = u; _succ_num[u] = 1; _last_succ[u] = u; if (_supply[u] <= 0) { _forward[u] = false; _pi[u] = 0; _pred[u] = e; _source[e] = _root; _target[e] = u; _cap[e] = INF; _flow[e] = -_supply[u]; _cost[e] = 0; _state[e] = STATE_TREE; } else { _forward[u] = true; _pi[u] = -ART_COST; _pred[u] = f; _source[f] = u; _target[f] = _root; _cap[f] = INF; _flow[f] = _supply[u]; _state[f] = STATE_TREE; _cost[f] = ART_COST; _source[e] = _root; _target[e] = u; _cap[e] = INF; _flow[e] = 0; _cost[e] = 0; _state[e] = STATE_LOWER; ++f; } } _all_arc_num = f; } // Check feasibility if (_sum_supply < 0) { for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE; } } else if (_sum_supply > 0) { for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE; } } else { for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { if (_flow[e] != 0) return INFEASIBLE; } for (int e = _search_arc_num; e != _all_arc_num; ++e) { if (_flow[e] != 0) return INFEASIBLE; } } } // Shift potentials to meet the requirements of the GEQ/LEQ type // optimality conditions if (_sum_supply == 0) { if (_stype == GEQ) { Cost max_pot = std::numeric_limits::min(); for (int i = 0; i != _node_num; ++i) { if (_pi[i] > max_pot) max_pot = _pi[i]; } if (max_pot > 0) { for (int i = 0; i != _node_num; ++i) _pi[i] -= max_pot; } } else { Cost min_pot = std::numeric_limits::max(); for (int i = 0; i != _node_num; ++i) { if (_pi[i] < min_pot) min_pot = _pi[i]; } if (min_pot < 0) { for (int i = 0; i != _node_num; ++i) _pi[i] -= min_pot; } } } return OPTIMAL;
Note: See TracChangeset for help on using the changeset viewer.