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.
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).
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 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)
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.
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)
Thanks, Rob Falck _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
On Mon, Oct 21, 2013 at 3:59 AM, Rob Falck <[hidden email]> wrote:
Hi Rob, I assume you have seen https://github.com/scipy/scipy/pull/218 then? 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.
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.
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?
This looks pretty good from a quick browse (disclaimer: I haven't tried to understand your algorithm in detail).
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)?
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. Cheers, Ralf _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
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 <[hidden email]> wrote:
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?
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.
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".
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.
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.
- Rob Falck _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
This post has NOT been accepted by the mailing list yet.
In reply to this post by robfalck
Hi,
I really need help understand how scipy.optimise.linprog works. I tried replicating the example provided in the mailing list as shown below but I keep getting a different result. I am new to Python and do not seem to be able to see where the problem is coming from: example in mailing list: >>>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 My own script and result: c= [3,2] A_ub = [[2,1], [1,1], [1,0]] b_ub = [10,8,4] res= linprog(c=c,A_ub=A_ub,b_ub=b_ub, options={"disp":True}) print(res) Optimization terminated successfully. Current function value: -0.000000 Iterations: 0 status: 0 slack: array([ 10., 8., 4.]) success: True fun: -0.0 x: array([ 0., 0.]) message: 'Optimization terminated successfully.' nit: 0 Thanks Bayo |
Free forum by Nabble | Edit this page |