使用GEKKO优化Expected Improvement函数时遇模型表达式语法错误

出现语法错误的核心原因有三点：

• 混淆中间变量与等式约束：直接将==0.0赋值给变量，导致GEKKO解析表达式时出现无效语法（比如erf的表达式中混入了多余的等式判断）。
• EI公式实现错误：PDF和CDF计算时错误使用了x-y，而非贝叶斯优化中Expected Improvement（EI）公式要求的标准化变量Z。
• 括号嵌套不当：erf的近似表达式中多层括号嵌套导致运算符解析异常，触发"Missing operator"错误。

1. 用m.Intermediate()定义所有EI计算所需的中间变量，避免将等式约束混入变量定义。
2. 修正EI公式的变量引用，确保PDF和CDF基于正确的标准化变量Z计算。
3. 简化erf的近似表达式，调整括号结构以符合GEKKO的语法解析规则。
4. 修正训练数据索引与最优值的错误定义，保证模型输入的正确性。

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[:,0].reshape(-1,1)  # 修正：scaledData第一列为x特征，原代码取[:,1]错误
    y_tr = d_array[:,-1].reshape(-1,1)
    gp_model = gp_reg.fit(x_tr, y_tr)
    return gp_model 

# 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) 

# gekko model definition and solution 
m = GEKKO(remote=False)  # model definition
x = m.Var(0.4, lb=0, ub=1) # scaled variable
y, std = Gekko_GPR(gp_model, m).predict(x, return_std=True) # GP prediction with std    

# constants
epsilon = m.Const(0.01, 'epsilon')
# 修正：使用训练数据中y的最大值经缩放后的值作为最优值
best_y_scaled = scaler.scaler_y.transform(np.array([y_train.max()]).reshape(-1,1))[0][0]
best_y = m.Const(best_y_scaled, 'best_y')
pi_m = m.Const(np.pi, 'pi')

# 用Intermediate定义中间变量，避免语法错误
Z = m.Intermediate( (y - best_y - epsilon) / std )
# 标准正态分布PDF
pdf = m.Intermediate( (1/(std*m.sqrt(2*pi_m))) * m.exp(-0.5 * Z**2) )
# 简化erf近似表达式，修正括号结构
erf = m.Intermediate( (2/pi_m) * m.atan( Z * m.sqrt(2) * (1 + ( (Z*m.sqrt(2)/2)**4 )) ) )
# 标准正态分布CDF
cdf = m.Intermediate( 0.5 * (1 + erf) )

# 最大化Expected Improvement
ei = Z * std * cdf + std * pdf
m.Maximize(ei)

m.options.IMODE = 3 # steady state optimization
m.solve(disp=True)

# 输出结果
print(f"最优x（缩放后）: {x.value[0]}")
print(f"最优x（原始尺度）: {scaler.scaler_x.inverse_transform(np.array([x.value]).reshape(-1,1))[0][0]}")
print(f"最大EI值: {ei.value[0]}")

• 训练数据索引：原代码错误取缩放后数据的第二列作为x特征，修正为取第一列[:,0]。
• 最优值定义：将硬编码的1.0替换为训练数据中y最大值的缩放后值，符合贝叶斯优化的逻辑。
• 中间变量处理：所有EI计算的中间变量改用m.Intermediate()定义，彻底避免语法解析错误。
• EI公式修正：PDF和CDF基于标准化变量Z计算，还原了正确的Expected Improvement公式逻辑。

内容的提问来源于stack exchange，提问作者Daniel Rodriguez