[SciPy-User] Linear Programming via Simplex Algorithm

Tue Oct 22 21:35:44 EDT 2013

Hello Ralf,

First, my apologies that this went to SciPy Users, it was intended for the
Developers list and I'm redirecting this message there.

On Tue, Oct 22, 2013 at 2:08 PM, Ralf Gommers <ralf.gommers at gmail.com>wrote:
>
>
> On Mon, Oct 21, 2013 at 3:59 AM, Rob Falck <robfalck at gmail.com> wrote:
>
>> Hello,
>>
>> I've spent some time recently polishing a simplex-based linear
>> programming function.  I've seen traffic from time to time about including
>> such functionality in scipy.optimize but it always seems to have been
>> closed without inclusion.
>>
>
> Hi Rob, I assume you have seen https://github.com/scipy/scipy/pull/218then? Looks like it's functionality that we'd like to see in scipy.optimize
> but needs some more effort than anyone has been willing to spend so far.
>

I've seen pull request 218 and parts of it are very good.  The simplex
method itself is lean, but it doesn't support equality or lower-bound
constraints.  It also differs in that it returns a new class LPResult.  In
my approach I had tried to stay with the general scipy.optimize.Result
setup and use only those fields which are applicable.  Do you have any
guidance on which way would be preferable?

> My intent was to develop a linear-programming routine that, given a
>> non-standard linear programming problem, generates a canonical tableau and
>> solves it using the simplex algorithm. By nonstandard I mean that variables
>> may have negative values, and three types of constraints are supported
>> (lower-bound, upper-bound, and equality).
>>
>> I've named a top-level routine "linprog".  I suspect in the future that
>> scipy may include other algorithms besides the simplex to solve linear
>> programming problems, in which case linprog would serve as the main
>> function (similarly to the way minimize serves as the interfact to all nlp
>> routines).
>>
>
> A generic `linprog` interface sounds like a good idea. It looks like the
> API and the current implementation are fairly specific to your lpsimplex
> algorithm however. Compare how little minimize() does to the 400 LoC in
> linprog(). If you want to make linprog() generic maybe it would help if you
> would consider how you could integrate the algorithm of
> https://github.com/scipy/scipy/pull/218 into it without changing the API
> besides adding a "method" keyword.
>

I'd be happy to incorporate both, although theres a good bit of overlap
between the two contributions. I'd hate for that author's work to go
unused, so perhaps I can replace the core simplex solving routine of my
code with that one.  It's leaner and honestly I like a few things about it
better than mine.  That way my top-level simplex routine would rework a
problem with lower-bound and equality constraints into standard form, to be
solved by that underlying simplex routine.

Also, my implementation of linprog is long because it contains a great deal
of error checking.  I will rework a top-level linprog interface which
leaves most of the error-checking to the underlying linear programming
routines.

> For example, consider a relatively simple problem:
>>
>> Maximize:  f = 3*x[0] + 2*x[1]
>> Subject to:  2*x[0] + 1*x[1] <= 10
>>                   1*x[0] + 1*x[1] <= 8
>>                   1*x[0] + 0*x[1] <= 4
>>
>> where: 0 <= x[0] <= inf
>>           0 <= x[1] <= inf
>>
>> The inputs are such that the objective is specified by array c where f =
>> dot(c,x).
>> We have three upper-bound constraints, which we define using dot(A_ub,x)
>> <= b_ub
>>
>> >>>c = [3,2]
>> >>>b_ub = [10,8,4]
>> >>>A_ub = [[2,1],
>> >>>            [1,1],
>> >>>            [1,0]]
>> >>>res=linprog(c=c,A_ub=A_ub,b_ub=b_ub,objtype='max',disp=True)
>> Optimization terminated successfully.
>>          Current function value: 18.000000
>>          Iterations: 3
>>
>> I've written a suite of unit tests and have documented all functions.
>>  It's somewhat non-standard but I've also included two callback functions.
>>  linprog_verbose_callback prints the simplex tableau, pivot element, and
>> basic variables at each iteration of the simplex method.
>> linprog_terse_callback simply prints the value of the solution vector x at
>> each iteration.
>>
>
> Do the callbacks need to be separate public functions? My impression is
> no. And what about lpsimplex()? If linprog() is the generic interface maybe
> lpsimplex can be private?
>
>

I felt like a verbose callback that displays the simplex tableau at each
iteration would be useful for people who want it.  I can remove it and
include it as an example callback in documentation or something like that.

lpsimplex will be made private (or more likely replaced by the
implementation from pull request 218)

The only top-level function will be linprog, which like minimize, will
support a "method" argument.  For now, the only valid method will be
"simplex".

>  The code is currently in the linprog branch at
>> https://github.com/robfalck/scipy/tree/linprog/scipy
>> (relevant files are scipy/optimize/linprog.py,
>> scipy/optimize/__init__.py, and scipy/optimize/tests/test_linprog.py)
>>
>
> This looks pretty good from a quick browse (disclaimer: I haven't tried to
> understand your algorithm in detail).
>
>
>>  I welcome constructive criticism of the algorithm.  It's been designed
>> to detect degenerate cases due to unbounded problems, and it has a
>> relatively basic way to avoid cycling in the simplex solution.
>>
>
> Have you benchmarked your algorithm against other implementations? And/or
> checked the efficiency (typical should be 2N to 3N iterations, with N the
> size of the constraint matrix)?
>

I've checked my algorithm against some textbook examples and a few randomly
generated cases from
http://demonstrations.wolfram.com/TwoPhaseSimplexMethod/

Basically there are two important things I had to get right.  First, was
the linprog routine (to become _linprog_simplex) creating a correct tableau
given the problem at hand?  Secondly, given that the tableau was correct
was the algorithm choosing the correct pivots to get to the solution as
quickly as possible.  In both cases, from all my testing thus far, the
answer has been yes.  I've checked a basic cycling case and was able to
successfully break out of the cycling to converge to the correct solution,
although I wouldn't say my method of avoiding cycling is advanced.  If
someone has a relatively complex example to test this against, I'd be happy
to try it.

I doubt this implementaiton will be fast with dozens of variables and
constraints, but I don't have a good feeling for where the maximum
appropriate problem size should be.

> If there's interest in including this I'd be happy to proceed with
>> submitting a PR.  I still have a bit of cleanup to perform but it's
>> relatively close to being ready (pep8 compliance, etc)
>>
>
> Sounds good to me!
>
> In terms of completing the cleanup, I think fixing line length to 80 chars
> and breaking up linprog() into smaller logical units would be good to
> before submitting imho. Also, linprog.py should be renamed to _linprog.py.
>

I'll will work on this and submit a pull request when I'm satisfied with
the state of things.  I will keep you posted.

> Cheers,
> Ralf
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>

-- 
- Rob Falck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20131022/bc09a073/attachment.html>