On this page you will find examples of using Linear Programming, some discussed in class.
Examples from chapter 5, in CPLEX lp format.
To read into polymake, use e.g.
$P=lp2poly<Rational>('example-5.1.lp');
To solve with glpsol, use
glpsol --lp example-5.1.lp
A simple example of LP duality, from the lectures.
The ACM programming contest problem Color a tree turns out to be equivalent to scheduling with tree precedence constraints.
There is a fast solution to this first worked out (although not analyzed) by Horn in 1972.
I wanted to check my solution, so I modelled this as an integer program. Unlike Horn's algorithm, this takes no advantage of the special tree structure of the constraints. This is both good and bad: more general constraint DAGs can be modelled, but it is much slower to solve the integer program. Here some example data files, all trees
There are several examples in the glpsol examples collection
graph.mod shows the fanciest input format
color.mod is more interesting from a modelling point of view
vertex-cover.mod uses the graph from graph.mod for a vertex cover example.
independent-set.mod uses the graph from graph.mod for an independent-set example. In this example the difference between the LP relaxation and the integer optimal is not too bad, compared to the worst case.
- Here is the bipartite matching example we discussed in class.
This file can be used to find all of the vertices (solutions). Run with
lrs < bipartite.ine
It turns out there are only two vertices, both integer.
1 1 0 0 1 1 1 0 1 0 1 1 0 1 0 1
- Solving with a uniform objective, and
choosing
--interior
with glpsol, we get the non-integer optimal solution
x[1,A].val = 0.5 x[1,B].val = 0.5 x[2,A].val = 0.5 x[2,D].val = 0.5 x[3,C].val = 1.0 x[4,B].val = 0.5 x[4,D].val = 0.5
- Here is the largest disk example we discussed in class.
Here is the line fitting example we discussed in class.
Here is the simple linear separation example. The quadratic version can be visualized using racket (pdf)
We can also fix other parameters to obtain the model from the book, and another one. The solutions are visualized using racket (pdf)
Here is the simple flow example we discussed in class. LP format
Here is the a slightly fancier example from the glpk source. The original example was a directed flow; this has been brutally symmetrized (to make it consistent with the undirected flows in our book) by adding reverse arcs. A less inefficient symmetrization simply changes the capacity constraints. There is also a dot drawing of the network and the residual graph.
All of the examples from the AMPL Book can be found on ampl.com
To run these examples with glpsol, you will need a command line like
glpsol -m prod.mod -d prod.dat
note that solve
and display
statements can be added to model file
in GMPL. There is also printf
in GMPL for more control over
output. See Chapter 4 in the
GNU MathProg Reference
for more information. See also triangle for a simple example of using display.
Here are some simple examples from the first lecture, as GMPL.
Recently I was asked how to read mps (old school linear programming input) files. I couldn't think of a completely off the shelf way to do, so I write a simple c program to use the glpk library.
Of course in general you would want to do something other than print it out again.
Recently I suggested to some students that they could use the Gnu Linear Programming Toolkit from C++. Shortly afterwards I thought I had better verify that I had not just sent people on a hopeless mission. To test things out, I decided to try using GLPK as part of an ongoing project with Lars Schewe
The basic idea of this example is to use glpk to solve an integer program with row generation.
The main hurdle (assuming you want to actually write object oriented c++) is how to make the glpk callback work in an object oriented way. Luckily glpk provides a pointer "info" that can be passed to the solver, and which is passed back to the callback routine. This can be used to keep track of what object is involved.
Here is the class header
#ifndef GLPSOL_HH
#define GLPSOL_HH
#include "LP.hh"
#include "Vektor.hh"
#include "glpk.h"
#include "combinat.hh"
namespace mpc {
class GLPSol : public LP {
private:
glp_iocp parm;
static Vektor<double> get_primal_sol(glp_prob *prob);
static void callback(glp_tree *tree, void *info);
static int output_handler(void *info, const char *s);
protected:
glp_prob *root;
public:
GLPSol(int columns);
~GLPSol() {};
virtual void rowgen(const Vektor<double> &candidate) {};
bool solve();
bool add(const LinearConstraint &cnst);
};
}
#endif
The class LP
is just an abstract base class (like an interface for
java-heads) defining the add
method. The method rowgen
is virtual
because it is intended to be overridden by a subclass if row
generation is actually required. By default it does nothing.
Notice that the callback method here is static; that means it is
essentially a C
function with a funny name. This will be the
function that glpk calls when it wants from help.
#include <assert.h>
#include "GLPSol.hh"
#include "debug.hh"
namespace mpc{
GLPSol::GLPSol(int columns) {
// redirect logging to my handler
glp_term_hook(output_handler,NULL);
// make an LP problem
root=glp_create_prob();
glp_add_cols(root,columns);
// all of my variables are binary, my objective function is always the same
// your milage may vary
for (int j=1; j<=columns; j++){
glp_set_obj_coef(root,j,1.0);
glp_set_col_kind(root,j,GLP_BV);
}
glp_init_iocp(&parm);
// here is the interesting bit; we pass the address of the current object
// into glpk along with the callback function
parm.cb_func=GLPSol::callback;
parm.cb_info=this;
}
int GLPSol::output_handler(void *info, const char *s){
DEBUG(1) << s;
return 1;
}
Vektor<double> GLPSol::get_primal_sol(glp_prob *prob){
Vektor<double> sol;
assert(prob);
for (int i=1; i<=glp_get_num_cols(prob); i++){
sol[i]=glp_get_col_prim(prob,i);
}
return sol;
}
// the callback function just figures out what object called glpk and forwards
// the call. I happen to decode the solution into a more convenient form, but
// you can do what you like
void GLPSol::callback(glp_tree *tree, void *info){
GLPSol *obj=(GLPSol *)info;
assert(obj);
switch(glp_ios_reason(tree)){
case GLP_IROWGEN:
obj->rowgen(get_primal_sol(glp_ios_get_prob(tree)));
break;
default:
break;
}
}
bool GLPSol::solve(void) {
int ret=glp_simplex(root,NULL);
if (ret==0)
ret=glp_intopt(root,&parm);
if (ret==0)
return (glp_mip_status(root)==GLP_OPT);
else
return false;
}
bool GLPSol::add(const LinearConstraint&cnst){
int next_row=glp_add_rows(root,1);
// for mysterious reasons, glpk wants to index from 1
int indices[cnst.size()+1];
double coeff[cnst.size()+1];
DEBUG(3) << "adding " << cnst << std::endl;
int j=1;
for (LinearConstraint::const_iterator p=cnst.begin();
p!=cnst.end(); p++){
indices[j]=p->first;
coeff[j]=(double)p->second;
j++;
}
int gtype=0;
switch(cnst.type()){
case LIN_LEQ:
gtype=GLP_UP;
break;
case LIN_GEQ:
gtype=GLP_LO;
break;
default:
gtype=GLP_FX;
}
glp_set_row_bnds(root,next_row,gtype,
(double)cnst.rhs(),(double)cnst.rhs());
glp_set_mat_row(root,
next_row,
cnst.size(),
indices,
coeff);
return true;
}
}
All this is a big waste of effort unless we actually do some row generation. I'm not especially proud of the crude rounding I do here, but it shows how to do it, and it does, eventually solve problems.
#include "OMGLPSol.hh"
#include "DualGraph.hh"
#include "CutIterator.hh"
#include "IntSet.hh"
namespace mpc{
void OMGLPSol::rowgen(const Vektor<double>&candidate){
if (diameter<=0){
DEBUG(1) << "no path constraints to generate" << std::endl;
return;
}
DEBUG(3) << "Generating paths for " << candidate << std::endl;
// this looks like a crude hack, which it is, but motivated by the
// following: the boundary complex is determined only by the signs
// of the bases, which we here represent as 0 for - and 1 for +
Chirotope chi(*this);
for (Vektor<double>::const_iterator p=candidate.begin();
p!=candidate.end(); p++){
if (p->second > 0.5) {
chi[p->first]=SIGN_POS;
} else {
chi[p->first]=SIGN_NEG;
}
}
BoundaryComplex bc(chi);
DEBUG(3) << chi;
DualGraph dg(bc);
CutIterator pathins(*this,candidate);
int paths_found=
dg.all_paths(pathins,
IntSet::lex_set(elements(),rank()-1,source_facet),
IntSet::lex_set(elements(),rank()-1,sink_facet),
diameter-1);
DEBUG(1) << "row generation found " << paths_found << " realized paths\n";
DEBUG(1) << "effective cuts: " << pathins.effective() << std::endl;
}
void OMGLPSol::get_solution(Chirotope &chi) {
int nv=glp_get_num_cols(root);
for(int i=1;i<=nv;++i) {
int val=glp_mip_col_val(root,i);
chi[i]=(val==0 ? SIGN_NEG : SIGN_POS);
}
}
}
So ignore the problem specific way I generate constraints, the key
remaining piece of code is CutIterator
which filters the generated
constraints to make sure they actually cut off the candidate
solution. This is crucial, because row generation must not add
constraints in the case that it cannot improve the solution, because
glpk assumes that if the user is generating cuts, the solver doesn't
have to.
#ifndef PATH_CONSTRAINT_ITERATOR_HH
#define PATH_CONSTRAINT_ITERATOR_HH
#include "PathConstraint.hh"
#include "CNF.hh"
namespace mpc {
class CutIterator : public std::iterator<std::output_iterator_tag,
void,
void,
void,
void>{
private:
LP& _list;
Vektor<double> _sol;
std::size_t _pcount;
std::size_t _ccount;
public:
CutIterator (LP& list, const Vektor<double>& sol) : _list(list),_sol(sol), _pcount(0), _ccount(0) {}
CutIterator& operator=(const Path& p) {
PathConstraint pc(p);
_ccount+=pc.appendTo(_list,&_sol);
_pcount++;
if (_pcount %10000==0) {
DEBUG(1) << _pcount << " paths generated" << std::endl;
}
return *this;
}
CutIterator& operator*() {return *this;}
CutIterator& operator++() {return *this;}
CutIterator& operator++(int) {return *this;}
int effective() { return _ccount; };
};
}
#endif
Oh heck, another level of detail; the actual filtering actually
happens in the appendTo
method the PathConstraint class. This is
just computing the dot product of two vectors. I would leave it as an
exercise to the readier, but remember some fuzz is neccesary to to
these kinds of comparisons with floating point numbers. Eventually,
the decision is made by the following feasible
method of the
LinearConstraint
class.
bool feasible(const
Vektor<double> & x){ double sum=0; for (const_iterator
p=begin();p!=end(); p++){ sum+= p->second*x.at(p->first); }
switch (type()){
case LIN_LEQ:
return (sum <= _rhs+epsilon);
case LIN_GEQ:
return (sum >= _rhs-epsilon);
default:
return (sum <= _rhs+epsilon) &&
(sum >= _rhs-epsilon);
}
}