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)