Wrong convergence of a complex MINLP

Question

Good morning, I have a problem with my optimization problem. runncfspop is a function that simulates a cascade failure in a power grid where the variables x are the new lines to be added to an existing grid. It returns the cascade risk improvement compared to a reference grid. The function is quite complex, because it runs in matlab and uses a trained deep neural network to solve it. My problem is that once it evaluates the objective function, which is performed correctly, it gives me a value error as the value of this function is not integer. If I round this value using a "int(round(" of the objective function, gekko solves the problem and gives a solution. The problem is that this solution makes sense contraint-wise, but does not update the objective function. Namely, it finds a solution that has maximum 8 lines built (upper_bound_lines), and a cost that is correctly constrained, but it does not perform an evaluation of the objective function, whose value always stays the same as the one of the first evaluation, and has nothing to do with the one of the final solution. I can assure you that the loading of the files and the directories at the beginning of the code work correctly. I cannot upload all the other files involved as they are too many, but I can tell you that also the first evaluation with the initial variable values is correct. My questions therefore are, 1) why does it give me error if the objective function is not integer?. 2) Is it possible to optimize any function with gekko, or a deep neural network is something too complicated ?

from gekko import GEKKO
import numpy as np
import scipy.io
from gekko import GEKKO
import numpy as np
import matlab.engine
import torch
import utils
from runCFSsurrpop_gekko import runcfspop
import pandas as pd
import os


# load cost and branch characteristics
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
branch, cost, eng, model, ccdf_ref = utils.data_init_ccdf_ref(device)
root_dir = "C:\\Users\Asus\\Desktop\\semester_project\\TEP_Metaheuristics"
root_dir_mat = "C:\\Users\Asus\\Desktop\\semester_project\\TEP_Metaheuristics\\CascadesTEPforpy"
system = scipy.io.loadmat(root_dir_mat + "\\13_Power_system\\IEEE118\\sys.mat")
branch = system["sys"]["branch"][0][0]
cost = np.multiply(branch[:, 21], branch[:, 22])
system = scipy.io.loadmat(root_dir_mat + "\\13_Power_system\\IEEE118\\sys.mat")
branch = system["sys"]["branch"][0][0]

# Define model
m = GEKKO() 
m.options.SOLVER = 1  # APOPT is an MINLP solver
m.options.DIAGLEVEL = 0
m.options.IMODE = 3  # Steady state optimization
m.solver_options = ['minlp_maximum_iterations 5000', \
                    # 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 100', \
                    # 1 = depth first, 2 = breadth first
                    'minlp_branch_method 1', \
                    # maximum deviation from whole number
                    'minlp_integer_tol 0.00001', \
                    # covergence tolerance
                    'minlp_gap_tol 0.0001']

# Define variables
x = m.Array(m.Var, 186, value=1, lb=0, ub=1, integer=True)

# Define paramters
lower_bound_lines = m.Param(value=1)
upper_bound_lines = m.Param(value=8)
lower_bound_cost = m.Param(value= 0)
upper_bound_cost = m.Param(value= 0.2e+8)

# Define equations
m.Equation(m.sum(x) <= upper_bound_lines)
m.Equation(m.sum(x) >= lower_bound_lines)

m.Equation(m.sum(cost*x) >= lower_bound_cost)
m.Equation(m.sum(cost*x) <= upper_bound_cost)

# Define objective
m.Obj((runcfspop(x, eng, branch, device, model, ccdf_ref)))


# Solve
m.solve(disp=True)

print('')
print('Results')
print(x)


# print contraint
print('Objective: ' + str(m.options.objfcnval))

I tried do solve the issue with the non integer value of the function by approximate it with it's closest integer. This allowed the solver to work, but it does not work correctly.

Thank you so much, Federico.

Answer 1

The function runcfspop needs to be written in terms of Gekko variables and equations. A call to an external black-box function is not supported in Gekko because of the requirements for operator overloading for automatic differentiation. Gekko compiles the model into byte-code with compilation of the model to give solvers sparse gradient information. Other solver options such as scipy.optimize.minimize() support arbitrary black-box functions. Another option is to create a surrogoate model of the function evaluation with a 1D cspline, 2D bspline, or higher dimensional machine learned models.

Here is an example with scipy.optimize.minimize():

import numpy as np
from scipy.optimize import minimize

def objective(x):
    return x[0]*x[3]*(x[0]+x[1]+x[2])+x[2]

def constraint1(x):
    return x[0]*x[1]*x[2]*x[3]-25.0

def constraint2(x):
    sum_eq = 40.0
    for i in range(4):
        sum_eq = sum_eq - x[i]**2
    return sum_eq

# initial guesses
n = 4
x0 = np.zeros(n)
x0[0] = 1.0
x0[1] = 5.0
x0[2] = 5.0
x0[3] = 1.0

# show initial objective
print('Initial Objective: ' + str(objective(x0)))

# optimize
b = (1.0,5.0)
bnds = (b, b, b, b)
con1 = {'type': 'ineq', 'fun': constraint1} 
con2 = {'type': 'eq', 'fun': constraint2}
cons = ([con1,con2])
solution = minimize(objective,x0,method='SLSQP',\
                    bounds=bnds,constraints=cons)
x = solution.x

# show final objective
print('Final Objective: ' + str(objective(x)))

# print solution
print('Solution')
print('x1 = ' + str(x[0]))
print('x2 = ' + str(x[1]))
print('x3 = ' + str(x[2]))
print('x4 = ' + str(x[3]))

And the same example with Python Gekko as a comparison:

from gekko import GEKKO    
import numpy as np
m = GEKKO()
x = m.Array(m.Var,4,value=1,lb=1,ub=5)
x1,x2,x3,x4 = x
# change initial values
x2.value = 5; x3.value = 5
m.Equation(x1*x2*x3*x4>=25)
m.Equation(x1**2+x2**2+x3**2+x4**2==40)
m.Minimize(x1*x4*(x1+x2+x3)+x3)
m.solve()
print('x: ', x)
print('Objective: ',m.options.OBJFCNVAL)