Implementing a Python Solver with Customised Constraints

user_9148880 0
Viewed 2592

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:

  1. 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).
  2. 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.

linear-programming scipy python
1 Answers

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.

user_20383830
Related