False calculation and non-robust behaviour

[forflo]

Hi,

I played around with your Linear.Grammar package and also tried to solve some linear programs with Linear.Simplex.Primal. It seems like either I'm doing something completely wrong or there is a bug in your library.

False calculation

Following code:

import Linear.Grammar
import Linear.Simplex.Primal

calc :: [(String, Rational)]
calc = 
    let light = (12 :: Rational) .*. EVar "x" .+. (16 :: Rational) .*. EVar "y" .=>. ELit (59)
        middle = (1 :: Rational) .*. EVar "x" .+. (7 :: Rational) .*. EVar "y" .=>. ELit (12)
        heavy = (4 :: Rational) .*. EVar "x" .+. (2 :: Rational) .*. EVar "y" .=>. ELit (13)
        consumption = (25 :: Rational) .*. EVar "x" .+. (30 :: Rational) .*. EVar "y" .<=. ELit 205 
        obj = EVar "M" .==. (300 :: Rational) .*. EVar "x" .+. (400 :: Rational) .*. EVar "y" 
    in
        simplexPrimal (standardForm obj) (standardForm <$> [light, middle, heavy, consumption])

produces the output

[("M",17050 % 13),("x",67 % 26),("y",35 % 26),("s0",(-85) % 13),("s1",0 % 1),("s2",0 % 1),("s3",2605 % 26)]

which is incorrect. There is only one solution for this LP namely (x = 2.25, y = 2.0), which I tested using following matlab code.

CoefficientMatrix = [
    -12,    -16
    -1,     -7
    -4,     -2
    25,    30
    -1,      0
     0,     -1
];

ConstraintVector = [-59, -12, -13, 205, 0, 0]';

Objective = [-300, -400];
Solution = [2.25, 2]';

CoefficientMatrix * Solution <= ConstraintVector
;; returns [1,1,1,1,1,1]'

Non-robust behaviour

If I compile following code

import Linear.Grammar
import Linear.Simplex.Primal

calc :: [(String, Rational)]
calc = 
    let light = (12 :: Rational) .*. EVar "x" .+. (16 :: Rational) .*. EVar "y" .=>. ELit (59)
        middle = (1 :: Rational) .*. EVar "x" .+. (7 :: Rational) .*. EVar "y" .=>. ELit (12)
        heavy = (4 :: Rational) .*. EVar "x" .+. (2 :: Rational) .*. EVar "y" .=>. ELit (13)
        consumption = (25 :: Rational) .*. EVar "x" .+. (30 :: Rational) .*. EVar "y" .<=. ELit 205 
       --                                                                                                      --v-- HERE CHANGE
        obj = (300 :: Rational) .*. EVar "x" .+. (400 :: Rational) .*. EVar "y" .==. EVar "M"
    in
        simplexPrimal (standardForm obj) (standardForm <$> [light, middle, heavy, consumption])

simplexPrimal outputs

[("M",0 % 1),("x",0 % 1),("y",0 % 1),("s0",(-59) % 1),("s1",(-12) % 1),("s2",(-13) % 1),("s3",205 % 1)]

which is clearly wront.

Further context

I wrote a little EDSL loosely resembling IBM's OPL (a DSL for optimization problems) where I want to use your simplex algorithm. I'd be very thankful, if you'd fix the problem (if there is an actual bug) Could you also provide a few examples for your Linear.Grammar EDSL, please. How is the syntax for "minimize " and "maximize ...".

Best regards

[athanclark]

Hi there, I've actually tried rebuilding this a second time under this repo, and have tried to more closely follow best principals in solving, however I've actually run into a serious issue regarding simplical optimization itself - cycling. The algorithm I'm using in this repository is actually incorrect, specifically in it's cost ratio. As it turns out, the linear-simplex-haskell repository I've been working on recently cycles infinitely for degenerate cases, and it's pretty unclear how one basic feasible solution grows to degeneracy (there's something about it needing to be in full-row rank, but I'm not good enough with linear algebra to understand this yet).

My next goal for linear optimization is to actually implement an interior-point method, particularly Karmarkar's method (whether I choose the affine or projective variants is uncertain for now). It will be interesting integrated weighting in this method too, but right now I just can't devote all my time to the development on this. Hopefully I can soon :)

I'm sorry that this doesn't work, I hope to have a working implementation of Karmakar's before the new year, but I can't make any promises :\

[forflo]

Ok, thanks for your reply. The solution your incorrect algorithm outputs is close enough for demonstrational purposes. Because of that, and the lack of another working haskell implementation, I'll use it anyway for my implementaiton. However, I strongly suggest you to add some defect statement to your hackage description file and into the documentation of this repository. That would have saved me some 2 hours, but regardless I quite like the interface. So I'm curious about your next version!

Best regards

[athanclark]

Yeah I think I'm just going to deprecate it for now :\ It's pretty unfortunate