Can I use an objective function that uses multiple functions in binary programming in gekko?

Question

I am trying to minimize the L2 norm of the peak flows in a drainage system. I want to use gekko since it can handle integer-binary programming.

In my program, my objective is to minimize the peak flows. I created a function that

1. Inputs the vector in the inp file of the drainage model, runs, and saves it as a new file.
2. Extract the peak flows in the new file.
3. Returns the L2-norm of the peak flows.

Because of the complexity, I created functions for each step in my objective function. Thus, my obj function has functions. I tested my obj function with different inputs, and it runs smoothly.

The input is a vector with 8 elements. Each element must be 1 or 0. The only constraint is that the sum of the elements is 4.

Unfortunately, when I use gekko, it keeps having the error:

Exception: @error: Model Expression
 *** Error in syntax of function string: Invalid element: <functionfuncat0x12fb5
 b2e0>
 
Position: 1                   
 <functionfuncat0x12fb5b2e0>
 ?

I even tried having integer = False and just added a converter in my objective function. The error is still the same. I guess that the error came from my objective function. Can anyone help me? Your response would be well appreciated.

from gekko import GEKKO

m = GEKKO(remote=False)

x = m.Array(m.Var,8,lb=0,ub=1,integer=True)

x0,x1,x2,x3,x4,x5,x6,x7 = x

def Func(x):
#Inputs x in the inp file, determines the peak flows with the new x, and saves a new file.
    Mod_File('Example1.inp',x) 
#Creates the file name of the modified file. 
    Mod_File_name = Mod_File_Name('Example1.inp') 
#Extracts the Peak flows in the mod file and computes the L2 norm.
    L2_Norm = Norm_Peak_Flows(Mod_File_name ,Dict_Links('Example1.inp')) 

    return L2_Norm
    
x0.value = 0
x1.value = 1
x2.value = 0
x3.value = 1
x4.value = 0
x5.value = 1
x6.value = 0  
x7.value = 1

m.Minimize(Func)
m.Equation(sum(x)==4)
m.options.SOLVER=1
m.solve()
print(x)

Answer 1

Gekko can't use black-box models like Func(x) inside the optimization model because it needs symbolic expressions to construct gradients and solve the optimization problem. Before running, it compiles the symbolic expressions into byte-code with operator overloading to get exact gradients of the objective and constraints.

It is possible to build a surrogate model of the objective and then optimize that instead. This can be a linear regression model.

ŷ = β₀ + β₁x₀ + β₂x₁ + ... + β₇x₇

or any other nonlinear model such as a Gaussian Process Regression, Neural Network, etc. Gekko supports model imports from scikit-learn, gpflow, keras / tensorflow, etc.

ŷ = f(x₀,x₁, ... ,x₇)

Here is an example of how to generate training data from the Func() function.

import numpy as np
from sklearn.linear_model import LinearRegression
from itertools import product

# Generate binary combinations (x in {0,1}) – limit size for feasibility
n_samples = 100  # adjust depending on how expensive Func is
X_data = []
y_data = []

# Generate random binary samples
np.random.seed(1)
for _ in range(n_samples):
    x_sample = np.random.randint(0, 2, 8).tolist()
    y_sample = Func(x_sample)
    X_data.append(x_sample)
    y_data.append(y_sample)

X_data = np.array(X_data)
y_data = np.array(y_data)

Next, fit the surrogate model:

from sklearn.linear_model import LinearRegression

model = LinearRegression()
model.fit(X_data, y_data)

# Coefficients and intercept
coef = model.coef_     # β₁ to β₈
intercept = model.intercept_  # β₀

Then build the surrogate model in Gekko and optimize.

from gekko import GEKKO

m = GEKKO(remote=False)

# Create binary decision variables x0 to x7
x = m.Array(m.Var, 8, lb=0, ub=1, integer=True)

# Linear regression expression: yp = β₀ + β₁x₀ + ... + β₇x₇
yp = m.Var()

# Apply the linear equation
m.Equation(yp == intercept + sum(coef[i]*x[i] for i in range(8)))

# Constraint: total number of "on" variables is 4
m.Equation(sum(x) == 4)

# Objective: minimize the surrogate output
m.Minimize(yp)

# Solve
m.options.SOLVER = 1
m.solve(disp=True)

print('Optimized x:', [xi.value[0] for xi in x])

from sklearn.metrics import r2_score

y_pred = model.predict(X_data)
print("R² of surrogate model:", r2_score(y_data, y_pred))

If you only have 8 decision variables that are each 0 or 1, it is also easy just to evaluate all 2^8 = 256 combinations and select the lowest objective value. Because we're also constraining the sum of x to be 4, there are only 70 valid combinations. Even up to 20 variables, there are only ~1M options and the summation constraint would reduce that further. Mixed integer programming or heuristics are better if there are many more decision variables.

from itertools import combinations
import numpy as np

# Store best values
best_x = None
best_obj = float('inf')

def Func(x):
    Mod_File('Example1.inp',x) 
    Mod_File_name = Mod_File_Name('Example1.inp') 
    L2_Norm = Norm_Peak_Flows(Mod_File_name ,Dict_Links('Example1.inp')) 
    return L2_Norm

# Generate all binary vectors of length 8 with exactly 4 ones
for ones_indices in combinations(range(8), 4):
    x = [0]*8
    for idx in ones_indices:
        x[idx] = 1
    obj = Func(x)
    if obj < best_obj:
        best_obj = obj
        best_x = x

print("Best solution:", best_x)
print("Best objective value:", best_obj)