Implementing a Python Solver with Customised Constraints

Question

I have two variables who are related to each other and I want to find an optimal solution, which in this case is the minimum of their sum. For now, let's call them X and Y, and along with pre-defined constants, they add up to a set of "variables" s1 and s2 (that later feed the constraints):

105896649.59 + X = s1
    -6738.82 + Y = s2

While searching the SciPy docs, I have come across the linear programming solution, where I have the minimizing function (in this case, X + Y) and a set of inequality and equality constraints to which my variables are bound. In my case, they are as follow:

X >= 0, Y >= 0
s1 >= 1, s2 >= 1
s2 / (s1 + s2) = 0.0001%

For this specific case, the code was easily implementable:

from scipy.optimize import linprog

lstConst = [105896649.59, -6738.82]

# function to minimise: X + Y
c= [1, 1]

# left-hand side of the equation for s2 / (s1 + s2) = 0.0001%
# i.e., -0.000001 * X + 0.999999 * Y
Aeq = [[-0.000001, 0.999999]]

# right-hand side of the equation
beq = [0.000001 * (lstConst[0] + lstConst[1]) - lstConst[1]]

# ensures minimum can't be a negative number
minX = max(1, max(1 -lstConst[0], 0))
minY = max(1, max(1 -lstConst[1], 0))

X_bounds = (minX, None)

Y_bounds = (minY, None)

res = linprog(c, A_eq=Aeq, b_eq=beq, bounds=[X_bounds, Y_bounds])

So we have the values for X and Y to minimize the function on the x parameter:

In [1]: res.x
Out[1]: array([1.00000000e+00, 6.84471676e+03])

I would like to build upon this approach:

In fact, there is another set of restrictions: s1 and s2 must also be integers (note that X and Y have no problem with being floats).
Instead of defining a single value for the ratio between s1 and s2, I would supply a list of different possible ratios.

In essence, I would like to find the minimum values for the X + Y function given several different ratios between s1 and s2. This could be achievable by either iterating over a list to define Aeq and beq on each iteration or defining additional restrictions (if possible).

However, I'm clueless as to the integer restriction and how to make the linear programming algorithm take it into account.

If anyone has an alternative suggestion that uses a library/optimizer other than SciPy and linprog, that's also welcome.

FYI there is another typo in your code around max(1, max(1 -lstConst[0], 0)) which is equivalent to just max(1 -lstConst[0], 1) but the comments say this is suppose to only ensure the number is non negative. — spinkus
– spinkus, Commented Jun 11, 2019 at 23:07
Didn't realise that there was no need for being more verbose, I thought the sign of lstConst would affect the boundary if I didn't set it that way, thanks — manoelpqueiroz
– manoelpqueiroz, Commented Jun 12, 2019 at 0:36
In your implementation you do ~need it to code the inequalities, it's just that it should be minX = max(1 -lstConst[0], 0) not max(1, max(1 -lstConst[0], 0)) AFAICT. I used linprog inequality constraints for these. — spinkus
– spinkus, Commented Jun 12, 2019 at 1:35

spinkus · Accepted Answer · 2019-06-13 23:17:08Z

3

Firstly, restating the problem:

minimize x + y, subject to:

    k1 + x = s1
    k2 + y = s2
    x >= 0
    y >= 0
    s1 >= 1
    s2 >= 1
    s2 / (s1 + s2) = k3

Where:

    k1 = 105896649.59
    k2 = -6738.82
    k3 = 0.000001

Note, you don't need the s1 and s2 vars to code the problem in linprog. Without the s1 and s2 auxiliary vars the problem is:

minimize x + y, subject to:

  x >= 0
  y >= 0
  x + k1 >= 1,
  y + k2 >= 1,
  (1-k3)y - k3x = (k1 + k2)k3 - k2

Which is a bit easier to read and code in linprog:

import numpy as np
from scipy.optimize import linprog
k1, k2, k3 = 105896649.59, -6738.82, 0.000001
A_ub = -np.eye(2)
b_ub = [k1-1, k2-1]
A_eq = [[-k3, (1-k3)]]
b_eq = (k1 + k2)*k3 -k2
res = linprog([1,1], A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=[[0,None], [0, None]])
print(res)

This gives [0., 6844.71675549] where as you had x=1, because you've actually set the lower bounds on x and y to 1 (I think that is an typo ...) but it doesn't matter in the context of the question being asked:

ONTO THE QUESTION:

... I'm clueless as to the integer restriction and how to make the linear programming algorithm take it into account.

If anyone has an alternative suggestion that uses a library/optimizer other than SciPy and linprog, that's also welcome.

What your asking for is mixed integer linear programming (MILP). MILP and linear programming (LP), are typically solved with different algorithms, and a MILP problem is typically harder to solve exactly. SciPy Optimize doesn't support MILP. There are a number of open source tools that do such as OrTools and PySCIPOpt which is a Python wrapper over SCIP.

EXAMPLE IN PySCIPOpt:

PySCIPOpt is nice coz it has a constraint programming type API. In PySCIPOpt your problem is pretty easy to state in a readable form. Reintroducing the auxiliary vars, we can pretty much type in the constraints word for word:

from pyscipopt import Model

k1, k2, k3 = 105896649.59, -6738.82, 0.000001
model = Model()
x = model.addVar(vtype="CONTINUOUS", name="x", lb=0)
y = model.addVar(vtype="CONTINUOUS", name="y", lb=0)
s1 = model.addVar(vtype="CONTINUOUS", name="s1", lb=None, ub=None)
s2 = model.addVar(vtype="CONTINUOUS" name="s2", lb=None, ub=None)
o = model.addVar(vtype="CONTINUOUS", name="Objective Value", lb=0, ub=None)
model.addCons(k1 + x == s1)
model.addCons(k2 + y == s2)
model.addCons(s1 >= 1)
model.addCons(s2 >= 1)
model.addCons(s2/(s1+s2) == k3)
model.addCons(x + y == o)
model.setObjective(o, "minimize")
model.optimize()
print('x + y = o -> (%.4f + %.4f = %.4f)' % (model.getVal(x), model.getVal(y), model.getVal(o)))

Gives the same answer as linprog since it is just a linear program. However, since SCIP supports MILP we can introduce integer variables. To handle your case #1 just change s1 and s2 to integers:

...
s1 = model.addVar(vtype="INTEGER", name="s1", lb=None, ub=None)
s2 = model.addVar(vtype="INTEGER", name="s2", lb=None, ub=None)

Gives:

...
SCIP Status        : problem is solved [optimal solution found]
Solving Time (sec) : 0.00
Solving Nodes      : 1
Primal Bound       : +1.10089229999989e+05 (1 solutions)
Dual Bound         : +1.10089229999989e+05
Gap                : 0.00 %
x + y = o -> (103244.4100 + 6844.8200 = 110089.2300)

Which is quite a different solution ... but this is why MILP is not LP.

From the example above, and by reading the docs, you should be able figure out how to code your #2 case - basically something like 1/k3 becomes another integer variable in your model.

edited Jun 13, 2019 at 23:17

answered Jun 9, 2019 at 23:55

spinkus

8,7704 gold badges45 silver badges74 bronze badges

Sign up to request clarification or add additional context in comments.

5 Comments

manoelpqueiroz Over a year ago

Had a bit of trouble installing pyscipopt, but finally managed to make it running. Does the complexity of MILP make it run relatively slower? I ran for 15 minutes and it still had a 1.5% gap between the primal and dual bounds

spinkus Over a year ago

The exact code?? In general MILP is NP-hard (see wikipedia article linked ni answer) and yes harder than LP. But in this case you've only got two variables it should be easy to solve. It took 0.00 seconds on my system. I updated answer to show more of output. It could also be because of the availability of a solver that I have installed that you don't have. I'm not sure.

manoelpqueiroz Over a year ago

Yes, I actually reproduced your code to check it out for performance. Could it be because I'm using a specific distribution? In my case, since I have Anaconda, I had to install pyscipopt via pip, since there is no channel in Anaconda offering it for Windows

spinkus Over a year ago

I am on linux python 3.5 with pyscipopt 2.1.4. I also tried python 2.7 with pyscipopt 2.1.5 with the same result. I also have SCIP installed on this system, and I thought that it was a prerequisite, but maybe it works with degraded performance without it. Not sure. I would suggest reading the install instructions again github.com/SCIP-Interfaces/PySCIPOpt/blob/master/INSTALL.md. and then if you can't figure it out ask another SO question. Cheers.

manoelpqueiroz Over a year ago

I thought it most likely could be Windows itself. I'll be moving over to Linux eventually, but nevertheless, my question has been answered, so I'll be closing it. Thanks again!

Collectives™ on Stack Overflow

Implementing a Python Solver with Customised Constraints

1 Answer 1

5 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

5 Comments

Your Answer

Sign up or log in

Post as a guest

Related