Expected improvement optimization with GEKKO

Viewed 46

I am trying to optimize the expected improvement function for Bayesian optimization applications. For this, I am using the scikit-learn Gaussian Process model embedded into the GEKKO optimization suite. When solving the optimization model the following error is shown:

@error: Model Expression *** Error in syntax of function string: Missing operator

Position: 128
((0.5)((1+(((2/pi))(atan(((((2)((((v1-i320))/(((2)(sqrt(2))))))))((1+(((((v1-i320))/(((2)(sqrt(2))))))^(4))))))))-(0.0)))=0)

The code is below

import numpy as np
import pandas as pd
from gekko import GEKKO 
from gekko.ML import Gekko_GPR
from gekko.ML import CustomMinMaxGekkoScaler 
import sklearn.gaussian_process as gpr

# Training data 
x_train = np.array([0.6, 0.9, 0.3, 0.45, 1.05, 0.75, 0.15, 
                    0.225, 0.825, 1.125]).reshape(-1,1)
y_train = np.array([-0.809016994, 0.809016994, -0.309016994, -0.951056516, 
                    0.951056516, -1.83772E-16, 0.587785252, 0.156434465, 
                    0.4539905, 0.707106781]).reshape(-1,1)

# Additional information
lb = [0.0] # lower bound
ub = [1.2] # upper bound
n_dim = len(lb) # number of dimension
n_train = x_train.shape[0] # size of the training set

# Function to fit the Gaussian process
def gp_fit(data_s, gp_reg):    
    d_array = data_s.to_numpy()
    x_tr = d_array[:,1].reshape(-1,1)
    y_tr = d_array[:,-1].reshape(-1,1)
    gp_model = gp_reg.fit(x_tr, y_tr)
    return gp_model # it delivers the gp model object

# gekko scaler definition
data = pd.DataFrame(np.hstack((x_train, y_train)), columns=['x', 'y'])
features = ['x']
label = ['y']
scaler = CustomMinMaxGekkoScaler(data,features,label)
data_s = scaler.scaledData() # data scaled

# kernel and gp regressor definition
bounds_m   = (1e-4, 3000) # bounds for the hyperparameters
kernel_main = gpr.kernels.Matern(length_scale=np.ones(n_dim), 
                            length_scale_bounds=bounds_m, 
                            nu=2.5)
constant_kernel = gpr.kernels.ConstantKernel(1.0, constant_value_bounds=bounds_m)
white_kernel = gpr.kernels.WhiteKernel(1.0, noise_level_bounds=(1.13e-07, 1.83e-02))
K_cov = constant_kernel*kernel_main + white_kernel 
gp_regressor = gpr.GaussianProcessRegressor(kernel=K_cov, alpha=1e-8,
                                            optimizer='fmin_l_bfgs_b',
                                            n_restarts_optimizer=50,
                                            random_state=20)
# gp_model creation
gp_model = gp_fit(data_s, gp_regressor) # trainig the model with the data scaled

# gekko model definition and solution 
m = GEKKO(remote=False)  # model definition
x = m.Var(0.4, lb=0, ub=1) # definition of variables scaled
y, std = Gekko_GPR(gp_model, m).predict(x, return_std=True) # gp prediction with std    
# constants
epsilon = m.Const(0.01, 'epsilon')
best_y = m.Const(1.0, 'best_y')
pi_m = m.Const(np.pi, 'pi')
# equations
Z = (y - best_y - epsilon)/std == 0.0 
pdf = 1/(std*m.sqrt(2*pi_m))*m.exp(-0.5*((x-y)/std)**2) == 0.0
erf = 2/pi_m*m.atan(2*((x-y)/(2*m.sqrt(2)))*(1+((x-y)/(2*m.sqrt(2)))**4)) == 0.0
cdf = 0.5*(1+erf) == 0
m.Equations([Z, pdf, erf, cdf])
# objective function
ei = Z*std*cdf + std*pdf
m.Maximize(ei)
m.options.IMODE = 3 # steady state optimization
m.solve(disp=True)

1 Answers

I was able to fix your error, but I am unable to get it fully working. Here is what I suggest:

for your objective function and cdf function, you are using Gekko equations from variables like erf. I suggest reformatting some of that with gekko Intermediate values, like so:

# equations
tZ = m.Intermediate((y - best_y - epsilon)/std)
Z = tZ == 0.0 
tpdf = m.Intermediate(1/(std*m.sqrt(2*pi_m))*m.exp(-0.5*((x-y)/std)**2))
pdf = tpdf == 0.0
terf = m.Intermediate(2/pi_m*m.atan(2*((x-y)/(2*m.sqrt(2)))*(1+((x-y)/(2*m.sqrt(2)))**4)))
erf = terf == 0.0
tcdf = m.Intermediate(0.5*(1+terf))
cdf = tcdf == 0.0
m.Equations([Z, pdf, erf, cdf])
# objective function
ei = tZ*std*tcdf + std*tpdf

Changing this causes Gekko to throw a "TOO_FEW_DEGREES_OF_FREEDOM" Error, as you are trying to solve 4 equations with 1 variable. I suggest making these equations a soft constraint (trying to minimize them rather than set them to 0) or add additional variables to the problem statement.

Related