Constrained Exploration#
This section implements constrained exploration process where the next infill point is obtained by maximizing the uncertainty in model prediction while taking into account the constraints. The optimization problem statement is written as
where \(\hat{\sigma}(x)\) is the uncertainty in model prediction, \(\hat{g}_j(x)\) is the constraint model and \(J\) is the number of constraints. Below code imports required packages, defines modified branin function, constraint, and creates plotting data:
# Imports
import numpy as np
from smt.surrogate_models import KRG
from smt.sampling_methods import LHS, FullFactorial
import matplotlib.pyplot as plt
from pymoo.core.problem import Problem
from pymoo.algorithms.soo.nonconvex.de import DE
from pymoo.optimize import minimize
from pymoo.config import Config
Config.warnings['not_compiled'] = False
def modified_branin(x):
dim = x.ndim
if dim == 1:
x = x.reshape(1, -1)
x1 = x[:,0]
x2 = x[:,1]
b = 5.1 / (4*np.pi**2)
c = 5 / np.pi
t = 1 / (8*np.pi)
y = (x2 - b*x1**2 + c*x1 - 6)**2 + 10*(1-t)*np.cos(x1) + 10 + 5*x1
if dim == 1:
y = y.reshape(-1)
return y
def constraint(x):
dim = x.ndim
if dim == 1:
x = x.reshape(1, -1)
x1 = x[:,0]
x2 = x[:,1]
g = -x1*x2 + 30
if dim == 1:
g = g.reshape(-1)
return g
# Bounds
lb = np.array([-5, 0])
ub = np.array([10, 15])
# Plotting data
sampler = FullFactorial(xlimits=np.array( [[lb[0], ub[0]], [lb[1], ub[1]]] ))
num_plot = 400
xplot = sampler(num_plot)
yplot = modified_branin(xplot)
gplot = constraint(xplot)
Differential evolution (DE) from pymoo is used for minimizing the surrogate model. Below code defines problem class and initializes DE. Note how problem class uses the objective prediction uncertainty and constraint surrogate model.
# Problem class
class ConstrainedExploration(Problem):
def __init__(self, sm_func, sm_const):
super().__init__(n_var=2, n_obj=1, n_ieq_constr=1, xl=lb, xu=ub)
self.sm_func = sm_func # store the function surrogate model
self.sm_const = sm_const # store the constraint surrogate model
def _evaluate(self, x, out, *args, **kwargs):
out["F"] = - np.sqrt(self.sm_func.predict_variances(x)) # Standard deviation
out["G"] = self.sm_const.predict_values(x)
# Optimization algorithm
algorithm = DE(pop_size=100, CR=0.8, dither="vector")
Below block of code creates 10 training points and performs sequential sampling using constrained exploration. The maximum number of iterations is set to 30 and a convergence criterion is defined based on the value of maximum uncertainty.
sampler = LHS(xlimits=np.array( [[lb[0], ub[0]], [lb[1], ub[1]]] ) , criterion='ese')
# Training data
num_train = 10
xtrain = sampler(num_train)
ytrain = modified_branin(xtrain)
gtrain = constraint(xtrain)
# Variables
itr = 0
max_itr = 30
tol = 0.1
max_unc = [1]
bounds = [(lb[0], ub[0]), (lb[1], ub[1])]
corr = 'squar_exp'
fs = 12
# Sequential sampling Loop
while itr < max_itr and tol < max_unc[-1]:
print("\nIteration {}".format(itr + 1))
# Initializing the kriging model
sm_func = KRG(theta0=[1e-2], corr=corr, theta_bounds=[1e-6, 1e2], print_global=False)
sm_const = KRG(theta0=[1e-2], corr=corr, theta_bounds=[1e-6, 1e2], print_global=False)
# Setting the training values
sm_func.set_training_values(xtrain, ytrain)
sm_const.set_training_values(xtrain, gtrain)
# Creating surrogate model
sm_func.train()
sm_const.train()
# Find the minimum of surrogate model
result = minimize(ConstrainedExploration(sm_func, sm_const), algorithm, verbose=False)
# Computing true function value at infill point
y_infill = modified_branin(result.X.reshape(1,-1))
if itr == 0:
max_unc = [-result.F]
else:
max_unc.append(-result.F)
print("Max uncertainty: {}".format(max_unc[-1]))
# Appending the the new point to the current data set
xtrain = np.vstack(( xtrain, result.X.reshape(1,-1) ))
ytrain = np.append( ytrain, y_infill )
gtrain = np.append( gtrain, constraint(result.X) )
itr = itr + 1 # Increasing the iteration number
Iteration 1
Max uncertainty: [23.36932637]
Iteration 2
Max uncertainty: [19.30934963]
Iteration 3
Max uncertainty: [9.85329274]
Iteration 4
Max uncertainty: [13.43157144]
Iteration 5
Max uncertainty: [11.11581009]
Iteration 6
Max uncertainty: [2.85847076]
Iteration 7
Max uncertainty: [3.95954151]
Iteration 8
Max uncertainty: [4.30972981]
Iteration 9
Max uncertainty: [2.96315553]
Iteration 10
Max uncertainty: [2.49902351]
Iteration 11
Max uncertainty: [0.88637202]
Iteration 12
Max uncertainty: [0.54919759]
Iteration 13
Max uncertainty: [0.40213456]
Iteration 14
Max uncertainty: [0.36601961]
Iteration 15
Max uncertainty: [0.17449173]
Iteration 16
Max uncertainty: [0.10732988]
Iteration 17
Max uncertainty: [0.04599643]
Below code plots the convergence of the exploration process.
####################################### Plotting convergence history
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(np.arange(itr) + 1, max_unc, c="red", marker=".")
ax.plot(np.arange(itr) + 1, [tol]*(itr), c="black", linestyle="--", label="Tolerance")
ax.set_xlabel("Iterations", fontsize=14)
ax.set_ylabel(r"$\sigma _{max}$", fontsize=14)
ax.set_xlim(left=1, right=itr)
ax.grid()
ax.set_yscale("log")
ax.set_title("Constrained Exploration".format(itr), fontsize=15)
Text(0.5, 1.0, 'Constrained Exploration')
Note that process terminated before reaching maximum number of iterations since the \(\hat{\sigma}(x)\) is below given tolerance. This essentially means that model is more or less globally accurate and will be confident about its prediction. Next block minimizes the surrogate model prediction taking into account constraint model using differential evolution.
# Problem class
class SurrogatePrediction(Problem):
def __init__(self, sm_func, sm_const):
super().__init__(n_var=2, n_obj=1, n_constr=1, xl=lb, xu=ub)
self.sm_func = sm_func
self.sm_const = sm_const
def _evaluate(self, x, out, *args, **kwargs):
out["F"] = self.sm_func.predict_values(x)
out["G"] = self.sm_const.predict_values(x)
result = minimize(SurrogatePrediction(sm_func, sm_const), algorithm, verbose=False)
print("Minimum of surrogate:\n")
print("x* = {}".format(result.X))
print("f(x*) = {}".format(result.F))
print("g(x*) = {}".format(result.G))
Minimum of surrogate:
x* = [9.13142199 3.28535951]
f(x*) = [47.53360888]
g(x*) = [-1.96431561e-07]
The obtained optimum is quite close to the true global minimum which shows that a globally accurate model of objective and constraint functions is created using exploration process. But the number of infill points needed for accurate model is high which makes the process computationally expensive. The optimum can be obtained using less number of infill points by using a more efficient infill criterion.
Below block of code plots the infill points and minimum of the objective surrogate.
####################################### Plotting initial samples and infills
# Reshaping into grid
reshape_size = int(np.sqrt(num_plot))
X = xplot[:,0].reshape(reshape_size, reshape_size)
Y = xplot[:,1].reshape(reshape_size, reshape_size)
Z = yplot.reshape(reshape_size, reshape_size)
G = gplot.reshape(reshape_size, reshape_size)
# Level
levels = np.linspace(-17, -5, 5)
levels = np.concatenate((levels, np.linspace(0, 30, 7)))
levels = np.concatenate((levels, np.linspace(35, 60, 5)))
levels = np.concatenate((levels, np.linspace(70, 300, 12)))
fig, ax = plt.subplots(figsize=(8,6))
# Plot function
CS=ax.contour(X, Y, Z, levels=levels, colors='k', linestyles='solid', alpha=0.5, zorder=-10)
ax.clabel(CS, inline=1, fontsize=8)
# Plot constraint
ax.contour(X, Y, G, levels=[0], colors='r', linestyles='dashed')
ax.contourf(X, Y, G, levels=np.linspace(0,G.max()), colors="red", alpha=0.3, antialiased = True)
ax.annotate('$g_1$', xy =(2.0, 11.0), fontsize=14, color='b')
# Plot minimum
ax.scatter(xtrain[0:num_train,0], xtrain[0:num_train,1], c="red", label='Initial samples')
ax.scatter(xtrain[num_train:,0], xtrain[num_train:,1], c="black", label='Infills')
ax.plot(9.143, 3.281, 'g*', markersize=15, label="Constrained minimum")
ax.plot(result.X[0], result.X[1], 'bo', label="Minimum of surrogate")
# Asthetics
ax.set_xlabel("$x_1$", fontsize=14)
ax.set_ylabel("$x_2$", fontsize=14)
ax.set_title("Constrained Modified Branin function", fontsize=15)
ax.legend()
<matplotlib.legend.Legend at 0x7f8bab08b560>
The infills are added all over the feasible region and not close to minimum value which results in a globally accurate model. As mentioned earlier, exploration process is usually used to build a globally accurate model but is useful to escape a local minimum when combined with exploitation which will enable efficient optimization.
NOTE: Due to randomness in differential evolution, results may vary slightly between runs. So, it is recommend to run the code multiple times to see average behavior.
Final result:
Parameter |
True minimum |
Minimum of surrogate |
|---|---|---|
\(x_1^*\) |
9.143 |
9.131 |
\(x_2^*\) |
3.281 |
3.285 |
\(f(x_1^*, x_2^*)\) |
47.560 |
47.534 |
\(g_1(x_1^*,x_2^*)\) |
0.0 |
-1.964 \(\cdot 10^{-7}\) |