Linear Programming in Sage

What is a Linear Program ?

A linear program is the sum of two information :

• A linear function, called the objective, which is to be maximized or minimized ( for example 2 x + y )
• Linear constraints on the variables ( for example, 3 x + y < 2 and 2 x + 3 y < 8 )

What is a Mixed Integer Linear Program ?

It is simply a Linear Program such that some variables are forced to take integer values instead of real values. This difference becomes very important when one learns that solving a Linear Program can be done in polynomial time while solving a general Mixed Integer Linear Program is usually NP-Complete ( = it can take exponential time, according to a widely-spread belief that P is not equal to NP)

Why is Linear Programming so useful ?

Linear Programming is very useful in many Optimization and Graph-Theoretical problems because of its wide range of expression. Most of the time, a natural Linear Program can be easily written to solve a problem whose solution will be quickly computed thanks to the wealth of heuristics already contained in Linear Program Solvers. It is often hard to theoretically find out the execution time of a Linear Program, though they give very interesting results in practice.

For more information, you can consult the Wikipedia page dedicated to Linear Programming : http://en.wikipedia.org/wiki/Linear_programming

How can I solve a linear program using Sage ?

I have been trying to keep the use of MILP in Sage as easy as possible, and to have it look like what we draw on blackboard. If it takes two lines to write it on a blackboard, then it should be the same in Sage.

Sage can solve Linear Programs or Mixed Integer Linear Programs through the class 'MixedIntegerLinearProgram' defined in 'sage.numerical.mip'. To illustrate how it can be used, we will try to solve one easy problem

First, we define our MixedIntegerLinearProgram object as a maximization problem

  sage: p=MixedIntegerLinearProgram( maximization=True ) 

We then define the objective using p[1] and p[2] as variables. Not that we could have picked `p["x1"]` and `p["x2"]`, or even `p[("x",1)]` and `p[("x",2)]` : any hashable Sage object is a valid key, as it would be in an usual Python dictionary !

  sage: p.set_objective( 2*p[1]+p[2] ) 

The two constraints :

  sage: p.add_constraint( 3*p[1]+4*p[2], max = 2.5 ) 
  sage: p.add_constraint( 1.5*p[1]+0.5*p[2], min = 0.5, max = 4 ) 

Our problem is now defined, and we can now ask Sage to solve it.
( be assured that you have installed the Sage package for Coin-OR CBC or GLPK. If you do not, see the Solvers section at the bottom of this page )

  sage: p.solve()
  1.6666666666666665 

The value returned by the ``solve`` method is the optimal value of the objective functions. In order to obtain the optimal values, you have to use the ``get_value``

  sage: print "The optimal values are x_1="+str(p.get_values(p[1]))+"x_2="str(p.get_values(p[2]))

All in all, here is how to solve this system

  sage: p=MixedIntegerLinearProgram( maximization=True )
  sage: p.set_objective( 2*p[1]+p[2] )
  sage: p.add_constraint( 3*p[1]+4*p[2], max=2.5 )
  sage: p.add_constraint( 1.5*p[1]+0.5*p[2], max=4,min=0.5 )
  sage: p.solve()
  1.6666666666666665
  sage: print "The optimal values are x_1 = "+str(p.get_values(p[1]))+", x_2 = "+str(p.get_values(p[2]))
  The optimal values are x_1 = 0.833333333333, x_2 = 0.0

Dictionaries of variables in ``MixedIntegerLinearProgram``

For more complex Linear Programs, it is sometimes useful to associate many parameters to the same object. Instead of using just one dictionary, it is possible to create many of them using the method ``new_variable`` : if several objects `B1, B2, ...` all have both a "value" and a "cost", the two being variables of the same linear program, you can define two dictionaries to store them

  sage: value = p.new_variable()
  sage: cost = p.new_variable()

It is now possible to define constraints mixing the two different kinds of variables

  sage: p.add_constraint( cost[ B1 ] +  3*value[ B2 ], max = 9)

Some famous examples

Vertex Cover in a graph

In the Vertex Cover problem, we are given a graph `G` and we want to find a subset `S` of its vertices of minimal cardinality such that each edge `e` is incident to at least one vertex of `S`. In order to achieve it, we define a binary variable `b_v` for each vertex `v`.

  sage: g=graphs.PetersenGraph() 
  sage: p=MixedIntegerLinearProgram(maximization=False) 
  sage: b=p.new_variable() 
  sage: for (u,v) in g.edges(labels=None): 
  ...:      p.add_constraint(b[u]+b[v],min=1) 
  sage: p.set_objective(sum([b[v] for v in g]))
  sage: p.set_binary(b) 
  sage: p.solve()
  6.0

  sage: b = p.get_values(b)
  sage: print b
  {0: 0.0, 1: 1.0, 2: 0.0, 3: 1.0, 4: 1.0, 5: 1.0, 6: 1.0, 7: 1.0, 8: 0.0, 9: 0.0}

  sage: print [v for v in g if b[v]==1]
  [1, 3, 4, 5, 6, 7]

Maximum matching in a Graph

In the maximum matching problem, we are given a graph `G`, and we are looking for a set of edges `M` of maximum cardinality such that no two edges from `M` are adjacent :

 
  sage: g=graphs.PetersenGraph() 
  sage: p=MixedIntegerLinearProgram() 
  sage: b=p.new_variable() 
  sage: B = lambda x, y : b[Set([x,y])]

  sage: for u in g: 
  ...:    p.add_constraint(sum([B(u,v) for v in g.neighbors(u)]),max=1) 
  
  sage: p.set_objective(sum([B(u,v) for (u,v) in g.edges(labels=None)]))

  sage: p.solve()
  sage: b = p.get_values(b)
  sage: print [(u,v) for (u,v) in g.edges(labels=None) if B(u,v) == 1]
  [(0, 5), (1, 6), (2, 7), (3, 8), (4, 9)]

And the next step is ``p.solve()`` !

Maximum Independent Set

  sage: G = graphs.PetersenGraph()
  sage: LP = MixedIntegerLinearProgram(maximization=True)
  sage: b = LP.new_variable()

  sage: # We define the objective
  sage: LP.set_objective(sum([b[v] for v in G]))

  sage: # For any edge, we define a constraint
  sage: for (u,v) in G.edges(labels=None):
  ...       LP.add_constraint(b[u]+b[v],max=1)

  sage: # The variable b is binary
  sage: LP.set_binary(b)
  
  sage: # The .solve() functions returns the objective value
  sage: LP.solve()
  4.0
  sage: b_sol = LP.get_values(b)
  sage: print b_sol
  {0: 1.0, 1: 0.0, 2: 1.0, 3: 0.0, 4: 0.0, 5: 0.0, 6: 0.0, 7: 0.0, 8: 1.0, 9: 1.0}

Solvers

Sage solves linear programs by calling specific libraries. Two are available for the moment :

• GLPK (website) : A Linear Program solver from GNU
• CBC (website) : Integer Linear Program solver from COIN-OR
• CPLEX (website) : Integer Linear Program solver from IBM ILOG. CPLEX is proprietary but free for researchers and students !

To install them if they are not available on your installation of Sage, type :

  sage: # To install GLPK 
  sage: install_package('glpk') 
  sage: # To install Coin-OR Branch and Cut ( CBC ) 
  sage: install_package('cbc') 

  sage -b

Get help and documentation

• The Sage-support Google group (here)
• Sage's documentation page (here), containing for example :
  • The tutorials (basic explanations on several topics)
  • The reference manual containing help for all of Sage's functions
• Sage's IRC Channel : sage-support@irc.freenode.net
• and .... do not hesitate to send me an email if you have any question ! ;-)