Python GEKKO: optimizing performance of nonlinear optimization

Question

I'm using GEKKO‍ for solving a nonlinear programming problem. My goal is to compare the GEKKO‍ performance with an alternative, and so I want to make sure that I'm getting from GEKKO‍ the best it can offer.

There are n binary variables, each of them is assigned a weight, each of the weights is a number from the interval [0, 1] (i.e., a rational number w satisfying 0<= w <= 1). Each of the constraints is linear. The objective function is non-linear: it is the product of weights of the non-zero variables, and the goal is to maximize the product.

I started with specifying the objective function as

m.Obj(-np.prod([1 - variables[i] + weights[i] * variables[i] for i in range(len(variables))]))

but then I would bump into APM model error: string > 15000 characters. So I switched to auxiliary variables using the if3 function as

aux_variables = [m.if3(variables[i], weights[i], 1) for i in range(len(variables))]
m.Obj(-np.prod(aux_variables))

The only global parameters that I'm setting manually are in the following code.

# initialize model

m = GEKKO(remote=False)


# set global variables

m.options.SOLVER = 1 # APOPT solver
# "APOPT is an MINLP solver"
# "APOPT is also the only solver that handles Mixed Integer problems."

m.options.IMODE = 3 # steady state optimization

m.solver_options = ['minlp_maximum_iterations 500', \
                    # minlp iterations with integer solution
                    'minlp_max_iter_with_int_sol 10', \
                    # treat minlp as nlp
                    'minlp_as_nlp 0', \
                    # nlp sub-problem max iterations
                    'nlp_maximum_iterations 50', \
                    # 1 = depth first, 2 = breadth first
                    'minlp_branch_method 1', \
                    # maximum deviation from whole number
                    'minlp_integer_tol 0.05', \
                    # covergence tolerance
                    'minlp_gap_tol 0.01']

# initialize variables
variables = m.Array(m.Var, (number_of_vars), lb=0, ub=1, integer=True)

# set initial values
for var in variables:
    var.value = 1

Question:

what more, if anything, can I do, in terms of global parameters and formulation of the objective function, to optimize the GEKKO‍'s performance for this particular problem?

At the same time, I would love GEKKO‍ to produce decent results.

John Hedengren John Hedengren · Accepted Answer · 2019-10-15T23:25:52

One way to reformulate your problem for increased speed is to use Intermediate variables.

Original (0.0325 sec with # Var=5)

m.Obj(-np.prod([1 - variables[i] + weights[i] * variables[i] \
      for i in range(len(variables))]))

Modified (0.0156 sec with # Var=5)

ival = [m.Intermediate(1 - variables[i] + weights[i] * variables[i]) \
                       for i in range(len(variables))]
m.Obj(-np.prod(ival))

This should also help you avoid the problem with string length, unless you have number_of_vars that is very large. It seems that the optimal solution will always be variables[i]=1 when weights[i]=1 and variables[i]=0 when weights[i]=0. With np.prod this means that the entire objective function is zero is any one of the product terms are zero. Would it help to set the individual product values equal to 1 instead of using the objective function to find the values? One thing that helps APOPT find a correct solution is to use something like 1.1 in your intermediate declaration instead of 1.0. So that when you are maximizing, it tries to avoid 0.1 values in favor of finding a solution that gives 1.1.

from gekko import GEKKO
import numpy as np
m = GEKKO(remote=False)
number_of_vars = 5
weights = [0,1,0,1,0]
m.options.IMODE = 3
variables = m.Array(m.Var, (number_of_vars), lb=0, ub=1, integer=True)
for var in variables:
    var.value = 1
ival = [m.Intermediate(1.1 - variables[i] + weights[i] * variables[i]) \
                       for i in range(len(variables))]
# objective function
m.Obj(-np.prod(ival))
# integer solution with APOPT
m.options.SOLVER = 1
m.solver_options = ['minlp_maximum_iterations 500', \
                    # minlp iterations with integer solution
                    'minlp_max_iter_with_int_sol 10', \
                    # treat minlp as nlp
                    'minlp_as_nlp 0', \
                    # nlp sub-problem max iterations
                    'nlp_maximum_iterations 50', \
                    # 1 = depth first, 2 = breadth first
                    'minlp_branch_method 1', \
                    # maximum deviation from whole number
                    'minlp_integer_tol 0.05', \
                    # covergence tolerance
                    'minlp_gap_tol 0.01']
m.solve()
print(variables)

It is also much easier for a solver to find a solution to a summation such as m.sum() and it gives the same variables solution as the np.prod() option.

# objective function
m.Obj(-m.sum(ival))

You could add a post-process line to recover product objective function that will be either 0 or 1.

The if3 function isn't a good option for your application because the switching condition is at 0 and slight numerical variations will cause unreliable results. The solver considers 0 to 0.05 and 0.95 to 1 to be integer solutions according to the option minlp_integer_tol=0.05. This is an option that allows integer solutions to be accepted when they are close enough to an integer value. If the variables[i] value is 0.01 then the if3 function will chose the True option when it should select the False option. You could still use the if3 function if you made the switch point between the binary values such as m.if3(variables[i]-0.5, weights[i], 1). However, there are easier ways to solve your problem than using the if3 function.

Python GEKKO: optimizing performance of nonlinear optimization

1 Answers