Constraining a Mixed-Integer Non-Linear optimization problem by L0-norm/number of non-zero elements in GEKKO

Question

I want to minimize the cost for achieving a given set of non-negative integer values b that are computed linearly from two sets of non-negative integer variables x, y in GEKKO. If have stated my problem in a way, that the b's are constraints on the x & y. My cost function is non-linear: quadratic using a conditional/minimum. In addition to these standard constraints, I have a constraints which requires the number of non-zero elements in x is at least as large as the largest element in x (e.g. that the L0-norm is equal to the LInifity-norm).

My difficulty is now twofold, since I'm fairly new to optimization and a novice in GEKKO.

1. I saw that GEKKO supports numpy arrays which would make the problem statement rather concise, but I struggle to get it to work - resulting in a lot of list comprehensions instead of vectorized operations.
2. I managed to define the L0-norm constraint and GEKKO actually runs with it, but it fails to find a solution. I recognize that L0 problems are really hard (e.g. combinatoric), but somehow the solution is fairly easy to find by "hand". I feel just that I'm doing something wrong.

I would appreciate any help! Here is what I have done so far:

from gekko import GEKKO
import numpy as np

# Setup gekko (taken from their MINLP tutorial with more iterations).
m = GEKKO()
m.options.SOLVER = 1  # APOPT is an MINLP solver
m.solver_options = ('minlp_maximum_iterations 500',
                    'minlp_max_iter_with_int_sol 10',
                    'minlp_as_nlp 0',
                    'nlp_maximum_iterations 50',
                    'minlp_branch_method 1',
                    'minlp_integer_tol 0.05',
                    'minlp_gap_tol 0.01')

# Define variables as arrays.
y = m.Array(m.Var, (7), value=1, lb=1, ub=12, integer=True)
x = m.Array(m.Var, (18), value=0, lb=0, ub=8, integer=True)

# Example of the user-defined target values b as constraints (actually input args).
m.Equation(x[2] + y[1] == 7)
m.Equation(x[12] + y[2] == 5)

# This is the L0 constraint.
# I thought using m.Array would make this a nice definition like 
#   m.Equation(np.count_nonzero(x) >= np.max(x))
# but it didn't, since these numpy functions are not supported in GEKKO.
# So I defined the following , which feels wrong:
m.Equation(m.sum([int(x_i.value > 0) for x_i in x]) >= max(x_i.value for x_i in x))

# Finally, the objective function (intermediates for readability).
k = np.array([m.min2(y_i, 3) for y_i in y])
x_cost = m.Intermediate(m.sum(x * (x + 1)))
y_cost = m.Intermediate(m.sum((k - 1) * (k + 2) + 2.5 * (y - k) * (y + k - 3)))
m.Obj(x_cost + y_cost)

# Solve.
m.solve(disp=True, debug=True)

Answer 1

Gekko uses gradient-based solvers so the equations shouldn't change iteration-to-iteration. There is a way to get the count of non-zero variables that is compatible with gradient based solvers. Here is an alternative form to give you count and max for your x vector. This uses if3, max3, and min3 as found in the documentation on Model Building Functions with logical conditions.

from gekko import GEKKO
import numpy as np
# Setup gekko (taken from their MINLP tutorial with more iterations).
m = GEKKO(remote=False)

# Define variables as arrays.
y = m.Array(m.Var, (7), value=1, lb=1, ub=12, integer=True)
x = m.Array(m.Var, (18), value=0, lb=0, ub=8, integer=True)

# Example of the user-defined target values b as constraints (actually input args).
m.Equation(x[2] + y[1] == 7)
m.Equation(x[12] + y[2] == 5)

# This is the L0 constraint.
#   m.Equation(np.count_nonzero(x) >= np.max(x))
eps = 0.05 # decision point for a "zero" value
count = m.sum([m.if3(x_i-eps,0,1) for x_i in x])
max_x = 0
for x_i in x:
    max_x = m.Intermediate(m.max3(max_x,x_i))
m.Equation(count >= max_x)

# Finally, the objective function (intermediates for readability).
k = np.array([m.min3(y_i, 3) for y_i in y])
x_cost = m.Intermediate(m.sum(x * (x + 1)))
y_cost = m.Intermediate(m.sum((k - 1) * (k + 2) + 2.5 * (y - k) * (y + k - 3)))
m.Minimize(x_cost + y_cost)

# Solve.
m.options.SOLVER = 3  # Initialize with IPOPT
m.solve(disp=True)

m.options.SOLVER = 1  # APOPT is an MINLP solver
m.solver_options = ('minlp_maximum_iterations 500',
                    'minlp_max_iter_with_int_sol 10',
                    'minlp_as_nlp 0',
                    'nlp_maximum_iterations 50',
                    'minlp_branch_method 1',
                    'minlp_integer_tol 0.05',
                    'minlp_gap_tol 0.01')
m.solve(disp=True)

I used IPOPT to give a non-integer solution for initialization and then APOPT to find the optimal Integer solution. This approach produces a successful solution with x

>>> x
array([[0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0],
       [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0]],
      dtype=object)

and y

>>> y
array([[1.0], [7.0], [5.0], [1.0], [1.0], [1.0], [1.0]], dtype=object)

You probably aren't intending a zero solution for x so you may need to add constraints such as m.Equation(count>=n) or change eps=0.05 to find a non-zero value or push the solver away from zero at a local minimum. With m.Equation(count>=n) and eps=0.5 it finds a better solution:

n=3, obj=63
n=5, obj=52

Here is the solution of x and y when obj=52 (best solution found).

>>> x
array([[0.0], [1.0], [4.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0],
       [0.0], [1.0], [0.0], [2.0], [0.0], [0.0], [0.0], [1.0], [0.0]],
      dtype=object)

>>> y
array([[1.0], [3.0], [3.0], [1.0], [1.0], [1.0], [1.0]], dtype=object)

You may need to adjust eps that is used in if3 to adjust when a value is counted as non-zero or adjust minlp_integer_tol 0.05 that is a solver option to determine the integer tolerance. Here is the final script with the additional inequality constraint that gives the best solution.

from gekko import GEKKO
import numpy as np
# Setup gekko (taken from their MINLP tutorial with more iterations).
m = GEKKO(remote=False)

# Define variables as arrays.
y = m.Array(m.Var, (7), value=1, lb=1, ub=12, integer=True)
x = m.Array(m.Var, (18), value=0, lb=0, ub=8, integer=True)

# Example of the user-defined target values b as constraints (actually input args).
m.Equation(x[2] + y[1] == 7)
m.Equation(x[12] + y[2] == 5)

# This is the L0 constraint.
#   m.Equation(np.count_nonzero(x) >= np.max(x))
eps = 0.5 # threshold for a "zero" value
count = m.sum([m.if3(x_i-eps,0,1) for x_i in x])
max_x = 0
for x_i in x:
    max_x = m.Intermediate(m.max3(max_x,x_i))
m.Equation(count >= max_x)

# Finally, the objective function (intermediates for readability).
k = np.array([m.min3(y_i, 3) for y_i in y])
x_cost = m.Intermediate(m.sum(x * (x + 1)))
y_cost = m.Intermediate(m.sum((k - 1) * (k + 2) + 2.5 * (y - k) * (y + k - 3)))
m.Minimize(x_cost + y_cost)

m.Equation(count>=5)

# Solve.
m.options.SOLVER = 3  # Initialize with IPOPT
m.solve(disp=True)

m.options.SOLVER = 1  # APOPT is an MINLP solver
m.solver_options = ('minlp_maximum_iterations 500',
                    'minlp_max_iter_with_int_sol 10',
                    'minlp_as_nlp 0',
                    'nlp_maximum_iterations 50',
                    'minlp_branch_method 1',
                    'minlp_integer_tol 0.05',
                    'minlp_gap_tol 0.01')
m.solve(disp=True)

You may be able to adjust n or some solver options to get a better solution. I hope this helps point you in the right direction.