Linear Programming

Problem definition.

Linear programming (LP) is not itself an approximation problem, as it can be solved exactly in polynomial time. The method however is used in many approximation algorithms (especially in all Iterative Rounding algorithms).

Linear programs can be expressed in canonical form:

\begin{eqnarray*} \mbox{maximize:} & \mathbf{c}^{T} \mathbf{x} \\ \mbox{subject to:} & A\mathbf{x} \leq \mathbf{b} \\ \mbox{and:} & \mathbf{x} \geq \mathbf{0} \\ \end{eqnarray*}

where \(\mathbf{x}\) represents the vector of variables to be determined, \(\mathbf{b}\) and \(\mathbf{c}\) are vectors of given coefficients and \(A\) is a given matrix of coefficients.

Solution

We provide an interface for creating, modifying and (re)solving LP instances using the primal and dual simplex methods. As the underlying implementation we use the GLPK (GNU Linear Programming Kit) library. Although the simplex method is not guaranteed to work in polynomial time, it usually works equally fast or faster then polynomial algorithms.

Example

#include "paal/lp/glp.hpp"

#include <iostream>

using LP = paal::lp::glp;

void print_solution(paal::lp::problem_type status, const LP &lp_instance) {
    std::cout << "-----------------------------------------" << std::endl;
    if (status == paal::lp::OPTIMAL) {
        std::cout << "Optimal solution cost: " << lp_instance.get_obj_value()
                  << std::endl;
    } else {
        std::cout << "Optimal solution not found" << std::endl;
    }
    std::cout << lp_instance << std::endl;
    // TODO example should show how to read a valuation of a variable, one
    // has to search documentation quite extensively to find that
}

int main() {
    // sample problem
    LP lp_instance;

    lp_instance.set_optimization_type(paal::lp::MAXIMIZE);
    // add_column(column_cost, lower_bound, upper_bound, symbolic_name)
    auto X = lp_instance.add_column(500, 0, paal::lp::lp_traits::PLUS_INF, "x");
    auto Y = lp_instance.add_column(300, 0, paal::lp::lp_traits::PLUS_INF, "y");

    // TODO example should should show how to create an expression with
    // 0 or 1 variable. This is the most standard use case, we use that
    // every time we create an expression in the loop
    auto expr = X + Y;
    std::cout << expr << std::endl;
    std::cout << (expr >= 7) << std::endl;
    lp_instance.add_row(expr >= 7);
    auto row = lp_instance.add_row(expr <= 10);
    lp_instance.add_row(15 <= 200 * X + 100 * Y <= 1200);

    // solve the LP
    auto status = lp_instance.solve_simplex();
    print_solution(status, lp_instance);

    // add new row
    expr += Y;
    lp_instance.add_row(expr == 12);
    // modify row bound
    lp_instance.set_row_upper_bound(row, 20);
    // modify column bound
    lp_instance.set_col_lower_bound(X, -50);
    // modify row expression
    lp_instance.set_row_expression(row, X + Y * 1.01);

    // resolve the LP
    status = lp_instance.resolve_simplex(paal::lp::DUAL);
    print_solution(status, lp_instance);

    // delete row
    lp_instance.delete_row(row);

    // resolve the LP
    status = lp_instance.resolve_simplex();
    print_solution(status, lp_instance);

    return 0;
}

complete usage example can be found in file lp_example.cpp

Row generation

Some of the LP formulations used in different algorithms have an exponential number of rows. To solve such instances we use a technique called row generation, which uses the concept of a "separation oracle". The separation oracle is an algorithm, which given the values of LP variables checks if they form a feasible solution to the LP instance and if not, it returns some violated LP row.

Using a separation oracle we can create a simple method of solving LPs with exponential number of rows:

 solve_lp_row_generation()
 {
     init(LP)
     solve(LP)
     while (not SeparationOracle.feasible_solution(LP))
     {
         LP <- SeparationOracle.add_violated_row(LP)
         resolve(LP)
     }
 }

In other words we start with a subset (maybe empty) of rows and solve the LP, and then while the obtained solution is infeasible add a violated row to the LP and resolve it.

Row generation interface

The row generation LP (re)solving is implemented by two classes: paal::lp::row_generation_solve_lp and paal::lp::row_generation_resolve_lp, both of which implement the following interface:

 row_generation_solve_lp {
     paal::lp::problem_type operator(Problem & problem, LP & lp);
 }

Here LP is the linear programming instance type.
Problem is the problem instance type: a class representing the problem described by the LP.
Problem must implement the following interface:

 ProblemArchetype {
     Oracle & get_oracle();
 }

Oracle is the Separation Oracle concept class. It is responsible for checking if the current LP solution is a feasible solution to the given Problem and if not, adding some violated constraint to the LP. It must implement the following interface:

 OracleArchetype {
    bool feasible_solution(const Problem & problem, const LP & lp);

    void add_violated_constraint(const Problem & problem, LP & lp);
 }

Separation oracle strategies

Typical separation oracles iterate over a set of the constraints that could be violated and choose one violated row to be added to the LP. The library provides a number of typical separation oracle strategies using this model:

1. paal::lp::most_violated_separation_oracle
2. paal::lp::first_violated_separation_oracle
3. paal::lp::random_violated_separation_oracle

To use one of those oracles the user must provide a ViolationsChecker concept class specific to the given problem. This class is responsible for:

1. providing the oracle with a set of constraints (candidates) that could be violated,
2. checking whether a given constraint (candidate) is violated, and if so by how much,
3. adding a given violated constraint to the LP.

ViolationsChecker must implement the following interface:

 ViolationsCheckerArchetype {
     CandidateIteratorRange get_violation_vandidates(const Problem & problem, const LP & lp);

     boost::optional<double> check_violation(Candidate candidate, const Problem & problem);

     void add_violated_constraint(Candidate candidate, const Problem & problem, LP & lp);
 }

CandidateIteratorRange is an iterator range of possible constraints (candidates) that could be violated. Candidate is a single constraint (candidate) that could be violated.

References

http://www.gnu.org/software/glpk/ - GLPK library (used in PAAl for LP solving)