0

I have a linear system derived from a network of pairs and I would like to optimize the minimum using a solver in python. An example for this system derived from a small network is summarized by the following data

obj_func
{('1', '2'): [1, 1, 1],
 ('2', '3'): [1, 1, 1, 2, 1, 0],
 ('2', '4'): [1, 1, 1, 2, 0, 1],
 ('3', '4'): [0, 1, 1, 1]}

rhs
{('1', '2'): [0.3333487586477922, 0.3333487586477922, 0.3333024827044157, 1.0],
 ('2', '3'): [0.5, 0.5, 0.3333487586477922, 0.3333487586477922, 0.3333024827044157],
 ('2', '4'): [0.49996529223940045, 0.5000347077605996, 0.3333487586477922, 0.3333487586477922, 0.3333024827044157],
 ('3', '4'): [0.49996529223940045, 0.5000347077605996, 0.5, 0.5]}

constraints
defaultdict(<class 'list'>,
  {('1', '2'): [[[1, 0, 0], [0, 1, 0], [0, 0, 1], [1, 1, 1]]],
   ('2', '3'): [[[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]]],
   ('2', '4'): [[[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]]],
   ('3', '4'): [[[1, 1, 0, 0], [0, 0, 1, 1], [1, 0, 1, 0], [0, 1, 0, 1]]]})

Previously, I did the optimization using lpSolve in R but for programming reasons I changed to python3 and I used Scipy.optimize.linprog(). Now I use lp_solve from lpsolve55 but the results are not the same as the one resulted from scipy.optimize.linprog() and LpSolve in R!

Here is the code using scipy.optimize.linprog()

from scipy.optimize import linprog

min_for_pairs = []
for pair in list(network.edges()):
  A = np.reshape(constraints[pair], (-1, len(obj_func[pair]))) 
  res = linprog(obj_func[pair], A_eq=A, b_eq=rhs[pair], method='interior-point', options={'presolve': True})
  min_for_pairs.append(res.fun)

min_for_pairs
[1.0, 0.6666975173156104, 0.666651241372254, 0.5000347083535648]

and this one using lp_solve:

from lpsolve55 import *
from lp_maker import *
from lp_solve import *

min_for_pairs = []
for pair in list(network.edges()):
  A = np.reshape(constraints[pair], (-1, len(obj_func[pair])))
  sense_equality = [0] * len(A)
  lp = lp_maker(obj_func[pair], A , rhs[pair], sense_equality)
  solvestat = lpsolve('solve', lp)
  obj = lpsolve('get_objective', lp)
  min_for_pairs.append(obj)

min_for_pairs
[1.0, 1.3333487586477921, 1.3333487586477921, 1.0]

I would like to know :

1) What is wrong in my code so that I get different results? Is this normal since lp_solve doesn't find the optimum?

2) I changed from scipy (to lpsolve55) because it was too slow when working on huge networks. For example, I used scipy to get the minimum for an objective function derived from a small network of less than 16000 pairs and it took more than 6 hours! Generally, are lp_solve and lp_maker better for huge systems? or do I need to change for another solver?

mattmilten
  • 6,242
  • 3
  • 35
  • 65
Noah16
  • 294
  • 2
  • 13
  • scipy's solvers are not competitive with the usual open-source tools like Clp, GLPK and maybe lp_solve too. In terms of different results, check the objective. That's the only metric the solver cares about. Not always unique optima are available. – sascha Mar 30 '19 at 10:00
  • @sascha Thanks for your comment. I didn't get this part of your comment "In terms of different results, check the objective. " What did you mean? – Noah16 Mar 31 '19 at 11:56
  • You look at the variable values in a solution. But the solver does not care about that (only in regards to constraints). Compare the objectives in both approaches and if they are equal, each solver is correct, without having equal variable solutions. Sometimes there are just multiple solutions. – sascha Mar 31 '19 at 11:58
  • @sascha But this is what I am doing! I didn't look at the variable values, I just compare the value of the objective between solvers. The resulted array ( min_for_pairs) contains the objective for solvers! where each array contains 4 objectives for each solver – Noah16 Apr 01 '19 at 09:04
  • @sascha and because of that I am wondering if the results are correct or not! – Noah16 Apr 01 '19 at 09:05
  • Can you maybe write how you formulate the model? I don't see where `network` is defined. – mattmilten Apr 02 '19 at 12:45
  • Also, you should give `PySCIPOpt` a try. It is generally more stable than `lp_solve` or `scipy.linprog`. https://github.com/SCIP-Interfaces/PySCIPOpt – mattmilten Apr 02 '19 at 12:57

1 Answers1

0

scipy is minimizing your objective.

Whereas the top level linprog module expects a problem of form:

Minimize:

c @ x
Subject to:

A_ub @ x <= b_ub
A_eq @ x == b_eq
 lb <= x <= ub
where lb = 0 and ub = None unless set in bounds.

But lp_maker maximizes the objective.

lp_maker.py
This script is analog to the lp_solve script and also uses the API to create a higher-level function called lp_maker. This function accepts as arguments some matrices and options to create an lp model. Note that this scripts only creates a model and returns a handle. Type help(lp_maker) or just lp_maker() to see its usage:

   LP_MAKER  Makes mixed integer linear programming problems.

   SYNOPSIS: lp_handle = lp_maker(f,a,b,e,vlb,vub,xint,scalemode,setminim)
      make the MILP problem
        max v = f'*x
          a*x <> b
            vlb <= x <= vub
            x(int) are integer

   ARGUMENTS: The first four arguments are required:
            f: n vector of coefficients for a linear objective function.
            a: m by n matrix representing linear constraints.
            b: m vector of right sides for the inequality constraints.
            e: m vector that determines the sense of the inequalities:
                      e(i) < 0  ==> Less Than
                      e(i) = 0  ==> Equals
                      e(i) > 0  ==> Greater Than
          vlb: n vector of non-negative lower bounds. If empty or omitted,
               then the lower bounds are set to zero.
          vub: n vector of upper bounds. May be omitted or empty.
         xint: vector of integer variables. May be omitted or empty.
    scalemode: Autoscale flag. Off when 0 or omitted.
     setminim: Set maximum lp when this flag equals 0 or omitted.

   OUTPUT: lp_handle is an integer handle to the lp created.
kutschkem
  • 7,826
  • 3
  • 21
  • 56