UNB/ CS/ David Bremner/ teaching/ cs6375/ Linear Programming Examples

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
example-5-pyramid.lp
Posted Fri 18 Mar 2016 08:12:01 AM ADT
example-5.1.lp
Posted Wed 02 Mar 2016 05:13:16 PM AST
example-5.2.lp
Posted Thu 03 Mar 2016 09:25:55 PM AST
example-5.2a.lp
Posted Thu 03 Mar 2016 09:22:18 PM AST
example-5.4a.lp
Posted Fri 18 Mar 2016 08:12:01 AM ADT
example-5.4b.lp
Posted Thu 03 Mar 2016 09:42:33 PM AST
Posted Wed 02 Mar 2016 12:00:00 AM AST Tags:
License: by-nc-sa-2.5

A simple example of LP duality, from the lectures.

Posted Thu 15 Nov 2012 12:00:00 AM AST Tags:
License: by-nc-sa-2.5

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

Posted Mon 01 Oct 2012 12:00:00 AM ADT Tags:
License: by-nc-sa-2.5

There are several examples in the glpsol examples collection

Posted Fri 28 Sep 2012 12:00:00 AM ADT Tags:
License: GPL3 or later
  • 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
  • This can be visualized as dot and pdf.
Posted Fri 21 Sep 2012 12:00:00 AM ADT Tags:
License: by-nc-sa-2.5
Posted Thu 20 Sep 2012 12:00:00 AM ADT Tags:
License: by-nc-sa-2.5
Posted Sun 16 Sep 2012 12:00:00 AM ADT Tags:
License: by-nc-sa-2.5
Posted Wed 12 Sep 2012 12:00:00 AM ADT Tags:
License: by-nc-sa-2.5

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.

Posted Fri 07 Sep 2012 12:00:00 AM ADT Tags:
License: by-nc-sa-2.5

Here are some simple examples from the first lecture, as GMPL.

Posted Thu 06 Sep 2012 12:00:00 AM ADT Tags:
License: by-nc-sa-2.5

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.

Posted Sat 12 Dec 2009 09:11:00 PM AST Tags:
License: by-nc-sa-2.5

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);
      }
    }
Tags:

Posted Wed 03 Dec 2008 08:53:00 PM AST
License: by-nc-sa-2.5
Examples from glpsol

There are several examples distributed with glpsol

Posted Sat 11 Oct 2008 08:52:00 PM ADT Tags:
License: GPL3 or later
Matching Students to Topics

In CS3997 this year, I have decided to have "debates" instead of presentations. This means that I need to make sure that each student has has exactly one topic, and every topic is assigned to an even number of students.

Another constraint is that I asked the students to list in order their top 3 preferences. I wanted to maximize (within reason) "student happiness", so I give weight 4 for their first choice, 2 for the second, and 1 for their third.

Finally, the students are numbered 1...26 in the order they sent me their preferences. I decided to enforce "first-come first-serve" in the objective function, so the happiness of student 1 has more weight than student 26. How much more is a bit of a subjective choice.

If you don't want to look at the solution yet, the students preferences are available separately. 'happy[i,j]' measures how happy student i is being assigned topic j

Almost the real solution is available. In actuallity, I first solved the problem for the first 18 students (so they didn't have to wait), and use the following

 printf { i in students, j in topics : x[i,j]=1 } "s.t. fix_%d_%d: x[%d,%d]=1;\n",i,j,i,j; 

to print out some constraints, which I then cut and past into the model, and resolved a week or so later when I had all of the data.

Posted Thu 25 Sep 2008 11:54:00 AM ADT Tags:
License: by-nc-sa-2.5
Cutting Rolls of paper

Here is the Paper roll cutting example from section 2.7.

Posted Wed 17 Sep 2008 12:00:00 AM ADT Tags:
License: by-nc-sa-2.5
Mmm, icecream

Here is the icecream production planner we discussed in class.

Posted Wed 10 Sep 2008 12:00:00 AM ADT Tags:
License: by-nc-sa-2.5
The Diet Example
Posted Mon 08 Sep 2008 12:00:00 AM ADT Tags:
License: GPL3 or later