10 Linear Programming Interface

Linear programming (LP) is one of the most important general methods of operations research. Countless optimization problems can be formulated and solved using LP techniques. Therefore, developing efficient LP solvers has been of high practical interest for a long time. Nowadays various efficient LP solvers are available, including both open source and commercial software packages. Therefore, LEMON does not implement its own solver, but it features wrapper classes for several known LP packages providing a common high-level interface for all of them.

The advantage of this approach is twofold. First, our C++ interface is more comfortable than the typical native interfaces of the solvers. Second, changing the underlying solver in a certain application using LEMON's LP interface needs no effort. So, for example, one may try her idea using an open source solver, demonstrate its usability for a customer and if it works well, but the performance should be improved, then the customer may decide to purchase and use a better commercial solver.

Currently, the following linear and mixed integer programming packages are supported: GLPK, Clp, Cbc, ILOG CPLEX and SoPlex. However, additional wrapper classes for new solvers can also be implemented quite easily.

10.1 Basic LP Features

The following example demonstrates how simple it is to formalize and solve an LP problem in LEMON. You can find the whole program in lp_demo.cc.

  // Create an instance of the default LP solver class
  // (it will represent an "empty" problem at first)
  Lp lp;

  // Add two columns (variables) to the problem
  Lp::Col x1 = lp.addCol();
  Lp::Col x2 = lp.addCol();

  // Add rows (constraints) to the problem
  lp.addRow(x1 - 5 <= x2);
  lp.addRow(0 <= 2 * x1 + x2 <= 25);
  
  // Set lower and upper bounds for the columns (variables)
  lp.colLowerBound(x1, 0);
  lp.colUpperBound(x2, 10);
  
  // Specify the objective function
  lp.max();
  lp.obj(5 * x1 + 3 * x2);
  
  // Solve the problem using the underlying LP solver
  lp.solve();

  // Print the results
  if (lp.primalType() == Lp::OPTIMAL) {
    std::cout << "Objective function value: " << lp.primal() << std::endl;
    std::cout << "x1 = " << lp.primal(x1) << std::endl;
    std::cout << "x2 = " << lp.primal(x2) << std::endl;
  } else {
    std::cout << "Optimal solution not found." << std::endl;
  }

Lp::Col type represents the columns, i.e. the (primal) variables of the LP problem. The numerical operators can be used to form expressions from columns, which are required for specifying rows (constraints) and the objective function. Due to the suitable operator overloads, a problem can be described in C++ conveniently, directly as it is expressed in mathematics. After specifying all the parameters of the problem, we can solve it using the solve() function. The status of the solution (namely OPTIMAL, UNBOUNDED, INFEASIBLE etc.) can be obtained using primalType() and the results can be obtained using primal().

The above problem has an optimal solution. If you execute this example, it will print out the following results.

Objective function value: 67.5
x1 = 7.5
x2 = 10

However, if you, for example, removed the lines that specify the rows, you would obtain an LP problem, for which the objective function value is unbounded over the feasible solutions. Thus the program would print out this line.

Optimal solution not found.

If you would like to build linear expressions or constraints in more steps, then you can use the classes Lp::Expr and Lp::Constr. For example, this line

lp.addRow(x1 - 5 <= x2);

can also be written as

Lp::Expr e1, e2;
e1 = x1;
e1 -= 5;
e2 = x2;
Lp::Constr c = e1 <= e2;
lp.addRow(c);

These classes make it easy to build more complex expressions and constraints, e.g. using loops. For example, we can sum all the variables (columns) in an expression using the Lp::ColIt iterator, which can be used the same as node and arc iterators for graphs.

Lp::Expr sum;
for (Lp::ColIt col(lp); col != INVALID; ++col) {
  sum += col;
}

After solving the problem, you can query the value of any primal expression, not just the individual columns and the objective function.

Todo:

Check this example after closing ticket #326.

std::cout << lp.primal(sum) << std::endl;

Of course, working with the dual problem is also supported by the LP interface. Lp::Row represents the rows, i.e. the dual variables of the LP problem.

Lp::Row r = lp.addRow(x2 >= 3);

The dual solutions of an LP problem can be obtained using the functions dualType() and dual() (after solving the problem).

std::cout << lp.dual(r) << std::endl;

Lp::DualExpr can be used to build dual expressions from rows and Lp::RowIt can be used to list the rows.

Lp::DualExpr dual_sum;
for (Lp::RowIt row(lp); row != INVALID; ++row) {
  dual_sum += row;
}

The LP solver interface provides several other features, which are documented in the classes LpBase and LpSolver. For detailed documentation, see these classes in the reference manual.

If you would like to use a specific solver instead of the default one, you can use GlpkLp, ClpLp, CplexLp etc. instead of the class Lp.

10.2 Mixed Integer Programming

LEMON also provides a similar high-level interface for mixed integer programming (MIP) solvers. The default MIP solver can be used through the class Mip, while the concrete solvers are implemented in the classes GlpkMip, CbcMip, CplexMip etc.

The following example demonstrates the usage of the MIP solver interface. The whole program can be found in mip_demo.cc.

  // Create an instance of the default MIP solver class
  // (it will represent an "empty" problem at first)
  Mip mip;

  // Add two columns (variables) to the problem
  Mip::Col x1 = mip.addCol();
  Mip::Col x2 = mip.addCol();

  // Add rows (constraints) to the problem
  mip.addRow(x1 - 5 <= x2);
  mip.addRow(0 <= 2 * x1 + x2 <= 25);
  
  // Set lower and upper bounds for the columns (variables)
  mip.colLowerBound(x1, 0);
  mip.colUpperBound(x2, 10);
  
  // Set the type of the columns
  mip.colType(x1, Mip::INTEGER);
  mip.colType(x2, Mip::REAL);
  
  // Specify the objective function
  mip.max();
  mip.obj(5 * x1 + 3 * x2);
  
  // Solve the problem using the underlying MIP solver
  mip.solve();

  // Print the results
  if (mip.type() == Mip::OPTIMAL) {
    std::cout << "Objective function value: " << mip.solValue() << std::endl;
    std::cout << "x1 = " << mip.sol(x1) << std::endl;
    std::cout << "x2 = " << mip.sol(x2) << std::endl;
  } else {
    std::cout << "Optimal solution not found." << std::endl;
  }

Todo:

Check this demo file after closing ticket #326. E.g. could we write mip.sol() instead of mip.solValue()? Or could we write primalValue() for Lp in lp_demo.cc?

In this program, the same problem is built as in the above example, with the exception that the variable x1 is specified to be INTEGER valued. x2 is set to REAL, which is the default option.

// Set the type of the columns
mip.colType(x1, Mip::INTEGER);
mip.colType(x2, Mip::REAL);

Because of this integrality requirement, the results will be different from the LP solution.

Objective function value: 67
x1 = 8
x2 = 9

The documentation of the MIP solver interface can be found in the reference manual at the class MipSolver. The common parts of the LP and MIP interfaces are documented in their common ancestor class LpBase.

10.3 Solving Complex Optimization Problems

The LP and MIP solvers are powerful tools for solving various complex optimization problems, as well. The following example solves a maximum flow problem with linear programming, see the whole program in lp_maxflow_demo.cc. Several other graph optimization problems can also be expressed as linear programs and the interface provided in LEMON facilitates solving them easily (though usually not so efficiently as by a direct combinatorial method, if one exists).

  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
  
  // Create an instance of the default LP solver
  Lp lp;

  // Add a column to the problem for each arc
  typename GR::template ArcMap<Lp::Col> f(g);
  lp.addColSet(f);

  // Capacity constraints
  for (ArcIt a(g); a != INVALID; ++a) {
    lp.colLowerBound(f[a], 0);
    lp.colUpperBound(f[a], capacity[a]);
  }

  // Flow conservation constraints
  for (NodeIt n(g); n != INVALID; ++n) {
    if (n == source || n == target) continue;
    Lp::Expr e;
    for (OutArcIt a(g, n); a != INVALID; ++a) e += f[a];
    for (InArcIt a(g, n); a != INVALID; ++a) e -= f[a];
    lp.addRow(e == 0);
  }

  // Objective function
  Lp::Expr o;
  for (OutArcIt a(g, source); a != INVALID; ++a) o += f[a];
  for (InArcIt a(g, source); a != INVALID; ++a) o -= f[a];
  lp.max();
  lp.obj(o);

  // Solve the LP problem
  lp.solve();

Note

This problem can be solved much more efficiently using common combinatorial optimization methods, such as the preflow push-relabel algorithm.

<< 9 Graph Adaptors | Home | 11 LEMON Graph Format >>