Integer decision variable in non linear programming

Question

I would like to maximize the quotient of two linear functions. I would want my decision variables to be Binary here i.e. they have to be integers and can take values of only 0 and 1.

I wanted to know how can I achieve this? I am looking to use an algorithm like SLSQP and I have looked at scipy but sadly it does not restrict the values of the decision variables to be binary and integer.

Does anyone know of a library with an easy to understand interface which I can use to achieve this? Or if there is any way to achieve this through scipy itself. I have read this question: Restrict scipy.optimize.minimize to integer values

But here out of the three solutions offered, I don't think any of them is efficient. It would be really helpful if any help could be provided.

Answer 1

Since you have no constraints, except that the variables should be binary, the maximization is quite simple. You can just sort the decision variables according to the ratios of the corresponding coefficients in the numerator and the denominator. Assuming that all coefficients are non-negative and there is a bias in the numerator and the denominator (to avoid divison by zero) you can use my implementation below.

import numpy as np

def maximize(numer, denom):
    """ 
    Function that maximizes an expression on the form

    a[0]*x[0] + a[1]*x[1] + ... + a[n-1]*x[n-1]
    -----------------------------------------
    b[0]*x[0] + b[1]*x[1] + ... + b[n-1]*x[n-1]

    where a[i] >= 0, b[i] >= 0, x[i] in [0,1] for 0 < i < n (non-negativity)
    and
    a[0] >= 0, b[0] > 0, x[0] = 1 (no division by zero)
    """

    ratios = numer / denom
    indices, ratios = zip(*sorted(enumerate(ratios), key = lambda x: - x[1]))
    decision = np.zeros_like(numer) 
    decision[0] = 1 # the bias is always enabled
    best_value = np.sum(decision * numer) / np.sum(decision * denom)
    for index, ratio in zip(indices, ratios):
        if index == 0:
            continue
        if ratio > best_value:
            decision[index] = 1 
            best_value = np.sum(decision * numer) / np.sum(decision * denom)
        else:
            # no more ratios can increase the cumulative ratio
            break  
    return decision

and here is a sample usage

if __name__ == "__main__":
    numer = np.array([1, 3, 4, 6])
    denom = np.array([1, 2, 2, 3])
    print("Input: {} / {}".format(",".join([str(x) for x in numer]), ",".join([str(x) for x in denom])))
    decision = maximize(numer, denom)
    print("Decision: {}".format(decision))
    print("Objective: {}".format(np.sum(decision * numer) / np.sum(decision * denom)))

Answer 2

I'm totally doing this off-the-cuff... but here's how I'd do it with mystic.

>>> equations = """
... 3.*x0 + 5.*x1 + 7.*x2 + 9.*x3 = 1.*x0 + 2.*x1 + 3.*x3
... """
>>> bounds = [(0,None)]*4
>>>
>>> def objective(x):
...   return x[0]**2 + 2*x[1] - 2*x[2] - x[3]**2
... 
>>> from mystic.symbolic import generate_penalty, generate_conditions
>>> pf = generate_penalty(generate_conditions(equations))
>>> from mystic.constraints import integers
>>> 
>>> @integers()
... def round(x):
...   return x
... 
>>> from mystic.solvers import diffev2
>>> result = diffev2(objective, x0=bounds, bounds=bounds, penalty=pf, constraints=round, npop=20, gtol=50, disp=True, full_output=True)
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 121
         Function evaluations: 2440
>>> result[0]
array([0., 0., 0., 0.])

Now modify the equations slightly...

>>> equations = """
... 3.*x0 + 5.*x1 + 7.*x2 + 9.*x3 = 5 + 1.*x0 + 2.*x1 + 3.*x3
... """
>>> pf = generate_penalty(generate_conditions(equations))
>>> result = diffev2(objective, x0=bounds, bounds=bounds, penalty=pf, constraints=round, npop=20, gtol=50, disp=True, full_output=True)
Optimization terminated successfully.
         Current function value: 3.000000
         Iterations: 102
         Function evaluations: 2060
>>> result[0]
array([1., 1., 0., 0.])

If you want binary variables instead of integers, then you can either use bounds = [(0,1)]*4 or replace @integers() with @discrete([0.0, 1.0]).

While the above isn't too interesting of a result, there are a few better thought out examples of global optimization with integer programming and generalized constraints on mystic's GitHub: https://github.com/uqfoundation/mystic/blob/master/examples2/integer_programming.py https://github.com/uqfoundation/mystic/blob/master/examples2/olympic.py

Answer 3

There are a few packages available for MINLP solutions in Python including pyomo and gekko. Here is a way to solve an MINLP problem with Python Gekko (a package that I maintain) as a simple example. Install the gekko package that includes the APOPT MINLP solver with pip:

pip install gekko

MINLP Solution

Gekko also solves Mixed Integer Nonlinear Programming (MINLP) problems such as:

from gekko import GEKKO

m = GEKKO(remote=False)
x = m.Array(m.Var,5,lb=0,ub=1,integer=True)

def f(x):
    return ((5+x[0])/(4+x[1])) \
           +(365.54/(3+x[2]))/(375.88/(3+x[3]))\
           +(379.75/(3+x[4]))

m.Minimize(f(x))
m.Equation(sum(x)==2)
m.options.SOLVER=1
m.solve()
print(x)

This gives the solution:

Iter: 1 I: 0 Tm: 0.00 NLPi: 4 Dpth: 0 Lvs: 3 Obj: 9.69E+01 Gap: NaN
--Integer Solution: 9.69E+01 Lowest Leaf: 9.69E+01 Gap: 2.89E-04
Iter: 2 I: 0 Tm: 0.00 NLPi: 1 Dpth: 1 Lvs: 3 Obj: 9.69E+01 Gap: 2.89E-04
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   9.000000001833541E-003 sec
 Objective      :    96.9099912206023     
 Successful solution
 ---------------------------------------------------
 

[[0.0] [1.0] [0.0] [0.0] [1.0]]

The APOPT solver uses a branch and bound solution approach with Nonlinear Programming (NLP) sub-problems to find integer solutions. There are several additional packages listed here: Python Mixed Integer Linear Programming for MILP (and some with MINLP) solvers. The Scipy package will have a Mixed Integer Linear Programming (MILP) solver at the next release, but this doesn't help with your MINLP problem. Gurobi, CPLEX, and Mosel-Xpress are the leaders in MILP/MIQP solvers, but are all commercial solvers. I also recently added an answer to the post you referenced: Restrict scipy.optimize.minimize to integer values in your search for an MINLP solver. If your problem can be reformulated to MILP then this opens up your solution to many other packages.

MILP Solution

Here is a script example that solves a linear programming problem with variables that are restricted to binary values (0 or 1) by specifying upper and lower bounds with integer=True:

from gekko import GEKKO
m = GEKKO()
x,y = m.Array(m.Var,2,integer=True,lb=0,ub=1) 
m.Maximize(y)
m.Equations([-x+y<=1,
             3*x+2*y<=12,
             2*x+3*y<=12])
m.options.SOLVER = 1
m.solve()
print('Objective: ', -m.options.OBJFCNVAL)
print('x: ', x.value[0])
print('y: ', y.value[0])

This generates the solution:

Iter: 1 I: 0 Tm: 0.00 NLPi: 2 Dpth: 0 Lvs: 0 Obj: -1.00E+00 Gap: 0.00E+00
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   1.369999999951688E-002 sec
 Objective      :   -1.00000000000000     
 Successful solution
 ---------------------------------------------------
 
Objective:  1.0
x:  1.0
y:  1.0