Блог

Я работаю над оптимизацией на основе суррогатов, где суррогатная модель генерируется методом регрессии гауссовых процессов в scikit-learn. Для части оптимизации я использовал scipy.minimize, но я действительно хочу использовать Pyomo, поскольку он предлагает более надежные решатели и возможность решения задач с целочисленными переменными. Ниже приведен пример кода, в котором я пытаюсь оптимизировать суррогатную модель функции стенда.

 import pyomo.environ as pyo
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn import gaussian_process as gpr
 
# Definition of function to express as surrogate model

def booth_func(x):
    f = (x[0]   2*x[1] - 7)**2   (2*x[0]   x[1] - 5)**2
    return f

# Bounds definition
lb = np.array([-10.0, -10.0])
ub = np.array([10.0, 10.0])

# Definition of training data
n_dim = lb.shape[0]
train_points = 20 
lhs_space  = lhs(n_dim, samples=train_points, criterion='cm') # Sampling via latin hypercube
x_train = lb   (ub-lb)*lhs_space

y_train = np.zeros((train_points, 1))
for i in range(train_points):
    y_train[i] = booth_func(x_train[i])

# Normalization of data
scaler_x = StandardScaler().fit(x_train)
scaler_y = StandardScaler().fit(y_train)
xn_train = scaler_x.transform(x_train)
yn_train = scaler_y.transform(y_train)
lbn = scaler_x.transform(lb.reshape(1, -1))
ubn = scaler_x.transform(ub.reshape(1, -1))


# GP definition
l_initial = np.ones(xn_train.shape[1])
matern = gpr.kernels.Matern(length_scale=l_initial, length_scale_bounds=(1e-5, 1e5), nu=2.5)
kernel = matern

gp = gpr.GaussianProcessRegressor(kernel=kernel, alpha=1e-6,
                                            optimizer="fmin_l_bfgs_b", 
                                            n_restarts_optimizer=20, 
                                            random_state=15)
# Regression via GP
gp_model = gp.fit(xn_train, yn_train)

#Definition of pyomo optimization problem
model = pyo.ConcreteModel()
model.s = pyo.Set(initialize=[1, 2])
model.x = pyo.Var(model.s, within=pyo.Reals, bounds=(lbn, ubn))

def obj_func(m):
    mu = gp_model.predict(m.x, return_std=False) # Let's assume a pure exploitation approach 
    mu = mu[0]
    return mu

model.obj = pyo.Objective(rule=obj_func, sense=pyo.minimize)
status = pyo.SolverFactory('ipopt').solve(model)

 

Есть какие-нибудь советы/помощь в этом?