EfficientGlobalOptimization examples

The EGO algorithm [jones1998] is an adaptative optimization method based on kriging.

An initial design of experiment is used to build a first metamodel. At each iteration a new point that maximizes a criterion is chosen as optimizer candidate.

The criterion uses a tradeoff between the metamodel value and the conditional variance.

Then the new point is evaluated using the original model and the metamodel is relearnt on the extended design of experiments.

from openturns.usecases import branin_function
from openturns.usecases import ackley_function
import openturns as ot
import openturns.viewer as viewer
from matplotlib import pylab as plt

ot.RandomGenerator.SetSeed(0)
ot.ResourceMap.SetAsString("KrigingAlgorithm-LinearAlgebra", "LAPACK")
ot.Log.Show(ot.Log.NONE)

Ackley test-case

We first apply the EGO algorithm on the Ackley function.

Define the problem

The Ackley model is defined in the usecases module in a data class AckleyModel :

am = ackley_function.AckleyModel()

We get the Ackley function :

model = am.model

We specify the domain of the model :

dim = am.dim
lowerbound = [-4.0] * dim
upperbound = [4.0] * dim

We know that the global minimum is at the center of the domain. It is stored in the AckleyModel data class.

xexact = am.x0

The minimum value attained fexact is :

fexact = model(xexact)
fexact

class=Point name=Unnamed dimension=1 values=[4.44089e-16]

graph = model.draw(lowerbound, upperbound, [100] * dim)
graph.setTitle("Ackley function")
view = viewer.View(graph)

We see that the Ackley function has many local minimas. The global minimum, however, is unique and located at the center of the domain.

Create the initial kriging

Before using the EGO algorithm, we must create an initial kriging. In order to do this, we must create a design of experiment which fills the space. In this situation, the LHSExperiment is a good place to start (but other design of experiments may allow one to better fill the space). We use a uniform distribution in order to create a LHS design. The length of the first LHS is set to ten times the problem dimension as recommended in [jones1998].

listUniformDistributions = [
    ot.Uniform(lowerbound[i], upperbound[i]) for i in range(dim)
]
distribution = ot.JointDistribution(listUniformDistributions)
sampleSize = 10 * dim
experiment = ot.LHSExperiment(distribution, sampleSize)
inputSample = experiment.generate()
outputSample = model(inputSample)

graph = ot.Graph(
    "Initial LHS design of experiment - n=%d" % (sampleSize), "$x_0$", "$x_1$", True
)
cloud = ot.Cloud(inputSample)
graph.add(cloud)
view = viewer.View(graph)

We now create the kriging metamodel. We selected the MaternModel covariance model with a constant basis as recommended in [leriche2021].

covarianceModel = ot.MaternModel([1.0] * dim, [0.5], 2.5)
basis = ot.ConstantBasisFactory(dim).build()
kriging = ot.KrigingAlgorithm(inputSample, outputSample, covarianceModel, basis)
kriging.run()

Create the optimization problem

We finally create the OptimizationProblem and solve it with class:~openturns.EfficientGlobalOptimization.

problem = ot.OptimizationProblem()
problem.setObjective(model)
bounds = ot.Interval(lowerbound, upperbound)
problem.setBounds(bounds)

In order to show the various options, we configure them all, even if most could be left to default settings in this case.

The most important method is setMaximumCallsNumber which limits the number of iterations that the algorithm can perform. In the Ackley example, we choose to perform 30 iterations of the algorithm.

algo = ot.EfficientGlobalOptimization(problem, kriging.getResult())
algo.setMaximumCallsNumber(30)
algo.run()
result = algo.getResult()

result.getIterationNumber()

30

result.getOptimalPoint()

class=Point name=Unnamed dimension=2 values=[-0.0062633,-0.0052574]

result.getOptimalValue()

class=Point name=Unnamed dimension=1 values=[0.0249089]

fexact

class=Point name=Unnamed dimension=1 values=[4.44089e-16]

Compared to the minimum function value, we see that the EGO algorithm provides solution which is accurate. Indeed, the optimal point is in the neighbourhood of the exact solution, and this is quite an impressive success given the limited amount of function evaluations: only 20 evaluations for the initial DOE and 30 iterations of the EGO algorithm, for a total equal to 50 function evaluations.

graph = result.drawOptimalValueHistory()
optimum_curve = ot.Curve(ot.Sample([[0, fexact[0]], [29, fexact[0]]]))
graph.add(optimum_curve)
view = viewer.View(graph)

inputHistory = result.getInputSample()

graph = model.draw(lowerbound, upperbound, [100] * dim)
graph.setTitle("Ackley function")
cloud = ot.Cloud(inputSample, "initial")
cloud.setPointStyle("bullet")
cloud.setColor("black")
graph.add(cloud)
cloud = ot.Cloud(inputHistory, "solution")
cloud.setPointStyle("diamond")
cloud.setColor("forestgreen")
graph.add(cloud)
view = viewer.View(graph)

We see that the initial (black) points are dispersed in the whole domain, while the solution points are much closer to the solution.

Branin test-case

We now take a look at the Branin-Hoo function.

Define the problem

The Branin model is defined in the usecases module in a data class BraninModel :

bm = branin_function.BraninModel()

We load the dimension,

dim = bm.dim

the domain boundaries,

lowerbound = bm.lowerbound
upperbound = bm.upperbound

and we load the model function and its noise :

objectiveFunction = bm.model
noise = bm.noiseModel

We build a sample out of the three minima :

xexact = ot.Sample([bm.xexact1, bm.xexact2, bm.xexact3])

The minimum value attained fexact is :

fexact = objectiveFunction(xexact)
fexact

	y0
0	-1.04741
1	-1.04741
2	-1.04741

graph = objectiveFunction.draw(lowerbound, upperbound, [100] * dim)
graph.setTitle("Branin function")
view = viewer.View(graph)

The Branin function has three local minimas.

Create the initial kriging

distribution = ot.JointDistribution([ot.Uniform(0.0, 1.0)] * dim)
sampleSize = 10 * dim
experiment = ot.LHSExperiment(distribution, sampleSize)
inputSample = experiment.generate()
outputSample = objectiveFunction(inputSample)
noiseSample = noise(inputSample)

graph = ot.Graph(
    "Initial LHS design of experiment - n=%d" % (sampleSize), "$x_0$", "$x_1$", True
)
cloud = ot.Cloud(inputSample)
graph.add(cloud)
view = viewer.View(graph)

covarianceModel = ot.MaternModel([1.0] * dim, [0.5], 2.5)
basis = ot.ConstantBasisFactory(dim).build()
kriging = ot.KrigingAlgorithm(inputSample, outputSample, covarianceModel, basis)

kriging.setNoise([x[0] for x in noiseSample])
kriging.run()

Create and solve the problem

We define the problem :

problem = ot.OptimizationProblem()
problem.setObjective(objectiveFunction)
bounds = ot.Interval(lowerbound, upperbound)
problem.setBounds(bounds)

We configure the algorithm, with the model noise:

algo = ot.EfficientGlobalOptimization(problem, kriging.getResult(), noise)
algo.setMaximumCallsNumber(30)

We run the algorithm and get the result:

algo.run()
result = algo.getResult()

result.getIterationNumber()

30

result.getOptimalPoint()

class=Point name=Unnamed dimension=2 values=[0.528177,0.167468]

result.getOptimalValue()

class=Point name=Unnamed dimension=1 values=[-1.04293]

fexact

	y0
0	-1.04741
1	-1.04741
2	-1.04741

inputHistory = result.getInputSample()

graph = objectiveFunction.draw(lowerbound, upperbound, [100] * dim)
graph.setTitle("Branin function")
cloud = ot.Cloud(inputSample, "initial")
cloud.setPointStyle("bullet")
cloud.setColor("black")
graph.add(cloud)
cloud = ot.Cloud(inputHistory, "solution")
cloud.setPointStyle("diamond")
cloud.setColor("forestgreen")
graph.add(cloud)
view = viewer.View(graph)

We see that the EGO algorithm reached the different optima locations.

graph = result.drawOptimalValueHistory()
optimum_curve = ot.Curve(ot.Sample([[0, fexact[0][0]], [29, fexact[0][0]]]))
graph.add(optimum_curve)
view = viewer.View(graph, axes_kw={"xticks": range(0, result.getIterationNumber(), 5)})
plt.show()

Reset default settings

ot.ResourceMap.Reload()