A 6 x 6 linprog testcase from glpk: 3 methods x sparse / dense give different results 2019-06-10 Jun 12:45z

0-transp-linprog.md

A 6 x 6 linprog testcase from glpk: 3 methods x sparse / dense give different results

transp-linprog.py below is a 6 x 6 testcase of scipy.optimize.linprog from a GLPK example, transp. In outline:

for method in ["interior-point", "revised simplex", "simplex"]:
    for sparse in [True, False]
        try:
            linprog( ... )

The 6 combinations give rather different results. grep transp-linprog.log:

linprog  method: interior-point  sparse: True 
The algorithm terminated successfully and determined that the problem is infeasible.
final f: 0.649008  x: [0.711 0.57 0.515 0.707 0.557 0.549] 

linprog  method: interior-point  sparse: False 
The algorithm terminated successfully and determined that the problem is unbounded.
final f: 0.649472  x: [0.712 0.57 0.515 0.708 0.557 0.55] 

linprog  method: revised simplex  sparse: True 
ValueError: shapes (1,12) and (1,12) not aligned: 12 (dim 1) != 1 (dim 0)

linprog  method: revised simplex  sparse: False 
final f: 153.675  x: [50 300 0 275 0 275] 

linprog  method: simplex  sparse: True 
ValueError: all the input arrays must have same number of dimensions

linprog  method: simplex  sparse: False 
final f: 153.675  x: [0 300 0 325 0 275] 

The full logfile and transp.* glpk files are under https://gist.github.com/denis-bz transp-linprog .

Transp is the smallest of the GLPK examples/, here made into a standalone testcase. (GLPK has a dozen more plain LPs, and ~ 30 MIPs. In the unlikely event that anyone would like to run them via a hack glp_sparse.py, please let me know.)

cheers
-- denis 10 June 2019

_{Theory and practice are closer in theory than they are in practice.}

transp-linprog.log

================================================================================
from transp-linprog.py Mon Jun 10 13:14:42 2019
versions: numpy 1.16.4  scipy 1.3.0   python 3.7.3 
c A: 
[0.225 0.153 0.162 0.225 0.162 0.126] 
[[0.225 0.153 0.162 0.225 0.162 0.126]
 [1 1 1 0 0 0]
 [0 0 0 1 1 1]
 [-1 0 0 -1 0 0]
 [0 -1 0 0 -1 0]
 [0 0 -1 0 0 -1]] 
b: [1e+10 350 600 -325 -300 -275] 
params: disp True  tol 0.0001 

--------------------------------------------------------------------------------
linprog  method: interior-point  sparse: True 
Primal Feasibility  Dual Feasibility    Duality Gap         Step             Path Parameter      Objective          
1.0                 1.0                 1.0                 -                1.0                 1.053               
0.03970108076208    0.03970100029991    26655.09069371      0.9710884895028  0.0397002956859     0.5294453395671     
2.237111193684e-06  2.237405684e-06     18.77027302619      0.9999436511367  2.236881518392e-06  0.649008009407      
The algorithm terminated successfully and determined that the problem is infeasible.
         Iterations: 2
final f: 0.649008  x: [0.711 0.57 0.515 0.707 0.557 0.549] 

--------------------------------------------------------------------------------
linprog  method: interior-point  sparse: False 
Primal Feasibility  Dual Feasibility    Duality Gap         Step             Path Parameter      Objective          
1.0                 1.0                 1.0                 -                1.0                 1.053               
0.03970102588196    0.03970106469025    67253.82154926      0.9710884895028  0.03970003971184    0.5294751484279     
2.237758684622e-06  2.237995969238e-06  225.9257602917      0.9999436347497  2.240250226874e-06  0.6494724617124     
The algorithm terminated successfully and determined that the problem is unbounded.
         Iterations: 2
final f: 0.649472  x: [0.712 0.57 0.515 0.708 0.557 0.55] 

--------------------------------------------------------------------------------
linprog  method: revised simplex  sparse: True 
transp-linprog.py:64: OptimizeWarning: Unknown solver options: sparse
  method=method, options=options )
Traceback (most recent call last):
  File "transp-linprog.py", line 64, in <module>
    method=method, options=options )
  File "/opt/local/mac/big/py3/lib/python3.7/site-packages/scipy/optimize/_linprog.py", line 532, in linprog
    c, c0=c0, A=A, b=b, x0=x0, callback=callback, _T_o=T_o, **solver_options)
  File "/opt/local/mac/big/py3/lib/python3.7/site-packages/scipy/optimize/_linprog_rs.py", line 543, in _linprog_rs
    mast, pivot, callback, _T_o, disp))
  File "/opt/local/mac/big/py3/lib/python3.7/site-packages/scipy/optimize/_linprog_rs.py", line 51, in _phase_one
    A, b, c, basis, x, status = _generate_auxiliary_problem(A, b, x0, tol)
  File "/opt/local/mac/big/py3/lib/python3.7/site-packages/scipy/optimize/_linprog_rs.py", line 195, in _generate_auxiliary_problem
    cols, rows = _select_singleton_columns(A, r)
  File "/opt/local/mac/big/py3/lib/python3.7/site-packages/scipy/optimize/_linprog_rs.py", line 254, in _select_singleton_columns
    same_sign = A[row_indices, column_indices]*b[row_indices] >= 0
  File "/opt/local/mac/big/py3/lib/python3.7/site-packages/numpy/matrixlib/defmatrix.py", line 220, in __mul__
    return N.dot(self, asmatrix(other))
ValueError: shapes (1,12) and (1,12) not aligned: 12 (dim 1) != 1 (dim 0)

--------------------------------------------------------------------------------
linprog  method: revised simplex  sparse: False 
transp-linprog.py:64: OptimizeWarning: Unknown solver options: sparse
  method=method, options=options )
Phase Iteration Minimum Slack       Constraint Residual Objective          
1     0         -325.0              0.0                 0.0                 
1     1         -300.0              0.0                 73.125              
1     2         -275.0              0.0                 76.95               
1     3         -275.0              0.0                 119.025             
1     4         -225.0              0.0                 127.125             
1     5         0.0                 0.0                 165.6               
2     5         0.0                 0.0                 165.6               
2     6         0.0                 0.0                 155.475             
2     7         0.0                 0.0                 153.675             
Optimization terminated successfully.
         Current function value: 153.675000  
         Iterations: 7
final f: 153.675  x: [50 300 0 275 0 275] 

--------------------------------------------------------------------------------
linprog  method: simplex  sparse: True 
transp-linprog.py:64: OptimizeWarning: Unknown solver options: sparse
  method=method, options=options )
Traceback (most recent call last):
  File "transp-linprog.py", line 64, in <module>
    method=method, options=options )
  File "/opt/local/mac/big/py3/lib/python3.7/site-packages/scipy/optimize/_linprog.py", line 526, in linprog
    c, c0=c0, A=A, b=b, callback=callback, _T_o=T_o, **solver_options)
  File "/opt/local/mac/big/py3/lib/python3.7/site-packages/scipy/optimize/_linprog_simplex.py", line 601, in _linprog_simplex
    row_constraints = np.hstack((A, np.eye(n), b[:, np.newaxis]))
  File "/opt/local/mac/big/py3/lib/python3.7/site-packages/numpy/core/shape_base.py", line 338, in hstack
    return _nx.concatenate(arrs, 0)
ValueError: all the input arrays must have same number of dimensions

--------------------------------------------------------------------------------
linprog  method: simplex  sparse: False 
Optimization terminated successfully.
         Current function value: 153.675000  
         Iterations: 18
final f: 153.675  x: [0 300 0 325 0 275] 

transp-linprog.py

#!/usr/bin/env python
""" linprog testcase from transp.glp: all 3 methods, sparse / dense """

from __future__ import division, print_function
import sys
import time
import traceback  # cf. better-exceptions
import numpy as np
import scipy
from scipy import sparse
from scipy.optimize import linprog  # $scopt13/_linprog.py
    # https://docs.scipy.org/doc/scipy/reference/optimize.linprog-interior-point.html

__version__ = "2019-06-10 june"
__author_email__ = "denis-bz-py t-online.de"

np.set_printoptions( threshold=10, edgeitems=5, linewidth=140,
        formatter = dict( float = lambda x: "%.3g" % x ))  # float arrays %.3g
print( "\n" + 80 * "=" )
print( "from", " ".join( sys.argv ), time.strftime( "%c" ))
print( "versions: numpy %s  scipy %s   python %s " % (
        np.__version__, scipy.__version__ , sys.version.split()[0] ))

#...............................................................................
disp = True
maxiter = 50
tol = 1e-4  # default 1e-8

    # from glpk /examples/transp.glp glp_linprog --
c = np.r_[ 0.225, 0.153, 0.162, 0.225, 0.162, 0.126 ]
lb = np.r_[ 0, 0, 0, 0, 0, 0 ]
ub = np.r_[ 1e10, 1e10, 1e10, 1e10, 1e10, 1e10 ]  # grr inf
b = np.r_[ 1e10, 350, 600, -325, -300, -275 ]

Adense = np.array(
    [[0.225, 0.153, 0.162, 0.225, 0.162, 0.126],
    [1, 1, 1, 0, 0, 0],
    [0, 0, 0, 1, 1, 1],
    [-1, 0, 0, -1, 0, 0],
    [0, -1, 0, 0, -1, 0],
    [0, 0, -1, 0, 0, -1]]
    )

    # to change these params, run this.py  a=1  b=None  'c = expr' ...
    # in sh or ipython --
for arg in sys.argv[1:]:
    exec( arg )


A = [Adense, sparse.csr_matrix( Adense )]
bounds = np.c_[ lb, ub ]
print( "c A: \n%s \n%s \nb: %s " % (c, Adense, b) )
print( "params: disp %s  tol %g " % (disp, tol) )

#...............................................................................
for method in ["interior-point", "revised simplex", "simplex"]:
  for sparse in [True, False]:
    print( "\n" +  80 * "-" )
    print( "linprog  method: %s  sparse: %s " % (method, sparse) )
    options = dict( disp=disp, maxiter=maxiter, tol=tol,
                sparse = sparse )  # rs, simplex OptimizeWarning
    try:
        res = linprog( A_ub=A[sparse], b_ub=b, c=c, bounds=bounds,
                method=method, options=options )
        if not disp:
            print( "linprog:", res.message )
        print( "final f: %g  x: %s " % (res.fun, res.x) )
            # glpk: 50 300 0 275 0 275
    except (IndexError, ValueError):
        traceback.print_exc()

transp.glp

p lp min 6 6 18
n p transp
n z cost
i 1 f
n i 1 cost
i 2 u 350
n i 2 supply[Seattle]
i 3 u 600
n i 3 supply[San-Diego]
i 4 l 325
n i 4 demand[New-York]
i 5 l 300
n i 5 demand[Chicago]
i 6 l 275
n i 6 demand[Topeka]
n j 1 x[Seattle,New-York]
n j 2 x[Seattle,Chicago]
n j 3 x[Seattle,Topeka]
n j 4 x[San-Diego,New-York]
n j 5 x[San-Diego,Chicago]
n j 6 x[San-Diego,Topeka]
a 0 1 0.225
a 0 2 0.153
a 0 3 0.162
a 0 4 0.225
a 0 5 0.162
a 0 6 0.126
a 1 1 0.225
a 1 2 0.153
a 1 3 0.162
a 1 4 0.225
a 1 5 0.162
a 1 6 0.126
a 2 1 1
a 2 2 1
a 2 3 1
a 3 4 1
a 3 5 1
a 3 6 1
a 4 1 1
a 4 4 1
a 5 2 1
a 5 5 1
a 6 3 1
a 6 6 1
e o f

transp.mod

# A TRANSPORTATION PROBLEM
#
# This problem finds a least cost shipping schedule that meets
# requirements at markets and supplies at factories.
#
#  References:
#              Dantzig G B, "Linear Programming and Extensions."
#              Princeton University Press, Princeton, New Jersey, 1963,
#              Chapter 3-3.

set I;
/* canning plants */

set J;
/* markets */

param a{i in I};
/* capacity of plant i in cases */

param b{j in J};
/* demand at market j in cases */

param d{i in I, j in J};
/* distance in thousands of miles */

param f;
/* freight in dollars per case per thousand miles */

param c{i in I, j in J} := f * d[i,j] / 1000;
/* transport cost in thousands of dollars per case */

var x{i in I, j in J} >= 0;
/* shipment quantities in cases */

minimize cost: sum{i in I, j in J} c[i,j] * x[i,j];
/* total transportation costs in thousands of dollars */

s.t. supply{i in I}: sum{j in J} x[i,j] <= a[i];
/* observe supply limit at plant i */

s.t. demand{j in J}: sum{i in I} x[i,j] >= b[j];
/* satisfy demand at market j */

data;

set I := Seattle San-Diego;

set J := New-York Chicago Topeka;

param a := Seattle     350
           San-Diego   600;

param b := New-York    325
           Chicago     300
           Topeka      275;

param d :              New-York   Chicago   Topeka :=
           Seattle     2.5        1.7       1.8
           San-Diego   2.5        1.8       1.4  ;

param f := 90;

end;

transp.sol

c Problem:    transp
c Rows:       6
c Columns:    6
c Non-zeros:  18
c Status:     OPTIMAL
c Objective:  cost = 153.675 (MINimum)
c
s bas 6 6 f f 153.675
i 1 b 153.675 0
i 2 u 350 0
i 3 b 550 0
i 4 l 325 0.225
i 5 l 300 0.153
i 6 l 275 0.126
j 1 b 50 0
j 2 b 300 0
j 3 l 0 0.036
j 4 b 275 0
j 5 l 0 0.00900000000000001
j 6 b 275 0
e o f