1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
|
"""
Optimization of the Rastrigin test function
===========================================
"""
# %%
# The Rastrigin function is defined by:
#
# .. math::
#
# f(\vect{x}) = A + \sum_{i=1}^n \left[x_i^2 - A\cos(2 \pi x_i)\right]
#
# where :math:`A=10` and :math:`\vect{x}\in[-5.12,5.12]^n`.
#
# It has a global minimum at :math:`\vect{x} = \vect{0}` where :math:`f(\vect{x})= - 10`.
#
# This function has many local minima, so optimization algorithms must be run from multiple starting points.
#
# In our example, we consider the bidimensional case, i.e. :math:`n=2`.
#
# **References**:
#
# - Rastrigin, L. A. "Systems of extremal control." Mir, Moscow (1974).
# - Rudolph, G. "Globale Optimierung mit parallelen Evolutionsstrategien". Diplomarbeit. Department of Computer Science, University of Dortmund, July 1990.
#
# %%
# Definition of the problem
# -------------------------
# %%
import openturns as ot
import openturns.viewer as otv
import math as m
def rastriginPy(X):
A = 10.0
delta = [x**2 - A * m.cos(2 * m.pi * x) for x in X]
y = A + sum(delta)
return [y]
dim = 2
rastrigin = ot.PythonFunction(dim, 1, rastriginPy)
print(rastrigin([1.0, 1.0]))
# %%
# Making `rastrigin` into a :class:`~openturns.MemoizeFunction` will make it recall all evaluated points.
# %%
rastrigin = ot.MemoizeFunction(rastrigin)
# %%
# This example is academic and the point achieving the global minimum of the function is known.
# %%
xexact = [0.0] * dim
print(xexact)
# %%
# The optimization bounds must be specified.
# %%
lowerbound = [-4.4] * dim
upperbound = [5.12] * dim
bounds = ot.Interval(lowerbound, upperbound)
# %%
# Plot the iso-values of the objective function
# ---------------------------------------------
# %%
graph = rastrigin.draw(lowerbound, upperbound, [100] * dim)
graph.setTitle("Rastrigin function")
graph.setLegendPosition("upper left")
graph.setLegendCorner([1.0, 1.0])
view = otv.View(graph)
# %%
# We see that the Rastrigin function has several local minima. However, there is only one single global minimum at :math:`\vect{x}^\star=(0, 0)`.
# %%
# Create the problem and set the optimization algorithm
# -----------------------------------------------------
# %%
problem = ot.OptimizationProblem(rastrigin)
# %%
# We use the :class:`~openturns.Cobyla` algorithm and run it from multiple starting points selected by a :class:`~openturns.LowDiscrepancyExperiment`.
# %%
size = 64
distribution = ot.JointDistribution([ot.Uniform(lowerbound[0], upperbound[0])] * dim)
experiment = ot.LowDiscrepancyExperiment(ot.SobolSequence(), distribution, size)
solver = ot.MultiStart(ot.Cobyla(problem), experiment.generate())
# %%
# Visualize the starting points of the optimization algorithm
# -----------------------------------------------------------
# %%
startingPoints = solver.getStartingSample()
graph = rastrigin.draw(lowerbound, upperbound, [100] * dim)
graph.setTitle("Rastrigin function")
cloud = ot.Cloud(startingPoints)
cloud.setPointStyle("bullet")
cloud.setColor("black")
graph.add(cloud)
graph.setLegends([""])
# sphinx_gallery_thumbnail_number = 2
view = otv.View(graph)
# %%
# We see that the starting points are well spread across the input domain of the function.
# %%
# Solve the optimization problem
# ------------------------------
# %%
solver.run()
result = solver.getResult()
xoptim = result.getOptimalPoint()
print(xoptim)
# %%
xexact
# %%
# We can see that the solver found a very accurate approximation of the exact solution.
# %%
# Analyze the optimization process
# --------------------------------
#
# :class:`~openturns.MultiStart` ran an instance of :class:`~openturns.Cobyla` from each starting point.
#
# Let us focus on the instance that found the global minimum. How many times did it evaluate `rastrigin`?
# %%
result.getCallsNumber()
# %%
# Let us view these evaluation points.
# %%
inputSample = result.getInputSample()
graph = rastrigin.draw(lowerbound, upperbound, [100] * dim)
graph.setTitle("Rastrigin function")
cloud = ot.Cloud(inputSample)
cloud.setPointStyle("bullet")
cloud.setColor("black")
graph.add(cloud)
graph.setLegendCorner([1.0, 1.0])
graph.setLegendPosition("upper left")
view = otv.View(graph)
# %%
# How fast did it find the global minimum?
# %%
graph = result.drawOptimalValueHistory()
view = otv.View(graph)
# %%
# Let us now analyze the :class:`~openturns.MultiStart` process as a whole.
#
# Since `rastrigin` is a :class:`~openturns.MemoizeFunction`,
# it has a :meth:`~openturns.MemoizeFunction.getInputHistory` method
# which lets us see all points it was evaluated on since its creation.
# %%
inputSample = rastrigin.getInputHistory()
graph = rastrigin.draw(lowerbound, upperbound, [100] * dim)
graph.setTitle("Rastrigin function")
cloud = ot.Cloud(inputSample)
cloud.setPointStyle("bullet")
cloud.setColor("black")
graph.add(cloud)
graph.setLegendCorner([1.0, 1.0])
graph.setLegendPosition("upper left")
view = otv.View(graph)
# %%
# How many times did all :class:`~openturns.Cobyla` instances combined call `rastrigin`?
print(rastrigin.getInputHistory().getSize())
# %%
otv.View.ShowAll()
|
