Tutorial
Moustache problem
\[\begin{array}{rl} (BB) \ \ \ \displaystyle \max_{x,y} & x\\ s.t. & y \in I(x)\\ \end{array}\]
The objective is to maximise $f(x, y) = x$ with $f$ not defined outside of a tight ribbon. At a given $x$, the interval of the admissible $y$ values is denoted $I(x) = [g(x) - \varepsilon(x)$, $g(x) + \varepsilon(x)]$, with
\[g(x) = -(|\cos(x)| + 0.1) \sin(x) + 2 ~~\text{and}~~ \varepsilon(x) = 0.05 + 0.05 \Bigg(1 + \dfrac{1}{1 + |x - 11|}\Bigg)\]
We can rewrite the problem such that it can be solved with NOMAD. We choose as starting point $(0, 2)$ and restrict the domain such that $x \in [0,20]$ and $y \in [0, 4]$.
\[\begin{array}{rrcll} (BB) \ \ \ \displaystyle - \min_{x \in \mathbb{R}^2} && -x_1\\ s.t. & & g(x_1) - \varepsilon(x_1) - x_2 & \leq 0 & \\ & & x_2 - g(x_1) + \varepsilon(x_1) & \leq 0 & \\ & 0 \leq & x_1 & \leq 20 & \\ & 0 \leq & x_2 & \leq 4 & \\ \end{array}\]
using NOMAD
# Objective
function f(x)
return -x[1]
end
# Constraints
function c(x)
g = -(abs(cos(x[1])) + 0.1) * sin(x[1]) + 2
ε = 0.05 + 0.05 * (1 - 1 / (1 + abs(x[1] - 11)))
constraints = [g - ε - x[2]; x[2] - g - ε]
return constraints
end
# Evaluator
function bb(x)
bb_outputs = [f(x); c(x)]
success = true
count_eval = true
return (success, count_eval, bb_outputs)
end
# Define Nomad Problem
p = NomadProblem(2, 3, ["OBJ"; "EB"; "EB"], bb,
lower_bound=[0.0;0.0],
upper_bound=[20.0;4.0])
# Fix some options
p.options.max_bb_eval = 1000 # total number of evaluations
p.options.display_stats = ["BBE", "EVAL", "SOL", "OBJ", "CONS_H"] # some display options
# Solution
result = solve(p, [0.0;2.0])
println("Optimization status: ", result.status)
println("Solution: ", result.x_sol)
println("is feasible: ", result.feasible)BBE EVAL ( SOL ) OBJ CONS_H
1 1 ( 0 2 ) -0 0
10 10 ( 2 1.44 ) -2 0
19 20 ( 3 1.84 ) -3 0
21 22 ( 4 2.5 ) -4 0
30 32 ( 6 2.22 ) -6 0
49 58 ( 6.14 2.12 ) -6.14 0
60 72 ( 6.16 2.22 ) -6.16 0
65 77 ( 6.22 2.14 ) -6.22 0
69 81 ( 6.23 2.07 ) -6.23 0
70 82 ( 6.27 2 ) -6.27 0
71 83 ( 6.32 1.99 ) -6.32 0
72 84 ( 6.4 1.91 ) -6.4 0
73 85 ( 6.45 1.77 ) -6.45 0
75 87 ( 6.58 1.68 ) -6.58 0
76 88 ( 6.74 1.52 ) -6.74 0
79 91 ( 6.79 1.53 ) -6.79 0
80 92 ( 6.96 1.41 ) -6.96 0
85 97 ( 6.99 1.43 ) -6.99 0
86 98 ( 7.12 1.39 ) -7.12 0
91 103 ( 7.14 1.42 ) -7.14 0
92 104 ( 7.23 1.43 ) -7.23 0
95 107 ( 7.27 1.39 ) -7.27 0
97 109 ( 7.38 1.43 ) -7.38 0
101 113 ( 7.39 1.46 ) -7.39 0
102 114 ( 7.45 1.5 ) -7.45 0
105 117 ( 7.48 1.49 ) -7.48 0
107 119 ( 7.55 1.56 ) -7.55 0
108 120 ( 7.64 1.63 ) -7.64 0
111 123 ( 7.67 1.67 ) -7.67 0
112 124 ( 7.77 1.76 ) -7.77 0
113 125 ( 7.9 1.86 ) -7.9 0
117 129 ( 7.93 1.9 ) -7.93 0
130 145 ( 8.14 1.67 ) -8.14 0
135 150 ( 8.35 1.43 ) -8.35 0
150 168 ( 8.37 1.41 ) -8.37 0
151 169 ( 8.39 1.39 ) -8.39 0
161 180 ( 8.77 1.47 ) -8.77 0
163 182 ( 9.35 1.94 ) -9.35 0
171 194 ( 9.93 2.43 ) -9.93 0
178 201 ( 10.07 2.54 ) -10.07 0
BBE EVAL ( SOL ) OBJ CONS_H
183 206 ( 10.14 2.61 ) -10.14 0
189 216 ( 11.14 2.26 ) -11.14 0
196 227 ( 12.14 2.46 ) -12.14 0
205 237 ( 13.14 1.48 ) -13.14 0
208 240 ( 15.14 1.48 ) -15.14 0
209 241 ( 20 1.48 ) -20 0
A termination criterion is reached: No termination (all). Mesh minimum precision reached (Algo)
Best feasible solutions: #36981 ( 20 1.48 ) Evaluation OK f = -20 h = 0 (L2)
#42545 ( 20 1.56 ) Evaluation OK f = -20 h = 0 (L2)
#44321 ( 20 1.47 ) Evaluation OK f = -20 h = 0 (L2)
#44398 ( 20 1.59 ) Evaluation OK f = -20 h = 0 (L2)
#44992 ( 20 1.4825 ) Evaluation OK f = -20 h = 0 (L2)
#45498 ( 20 1.5073 ) Evaluation OK f = -20 h = 0 (L2)
#45499 ( 20 1.5854 ) Evaluation OK f = -20 h = 0 (L2)
#45500 ( 20 1.5527 ) Evaluation OK f = -20 h = 0 (L2)
... A total of 89 feasible solutions were found.
Best infeasible solution: Undefined.
Blackbox evaluations: 400
Total model evaluations: 43620
Cache hits: 97
Total number of evaluations: 497
Optimization status: 1
Solution: [19.9999995, 1.48000001684481]
is feasible: trueLinear equality constraints
It is also possible to define linear equality constraints in a NomadProblem structure. NOMAD handles differently these constraints and solves the original problem on a reduced subspace. For more details we refer the reader to:
We consider the following example:
\[\begin{array}{rl} (BB) \ \ \ \displaystyle \min_{x \in \mathbb{R}^5} & (x_1 -1)^2 + (x_2 - x_3)^2 + (x_4 - x_5)^2\\ s.t. & x_1 + x_2 + x_3 + x_4 + x_5 = 5 \\ & x_3 -2 x_4 - 2 x_5 = -3\\ & -10 \leq x_1, x_2, x_3, x_4, x_5 \leq 10\\ \end{array}\]
This problem can be solved by Nomad the following way:
using NOMAD
# blackbox
function bb(x)
f = (x[1]- 1)^2 * (x[2] - x[3])^2 + (x[4] - x[5])^2
bb_outputs = [f]
success = true
count_eval = true
return (success, count_eval, bb_outputs)
end
# linear equality constraints
A = [1.0 1.0 1.0 1.0 1.0;
0.0 0.0 1.0 -2.0 -2.0]
b = [5.0; -3.0]
# Define blackbox
p = NomadProblem(5, 1, ["OBJ"], # the equality constraints are not counted in the outputs of the blackbox
bb,
lower_bound = -10.0 * ones(5),
upper_bound = 10.0 * ones(5),
A = A, b = b)
# Fix some options
p.options.max_bb_eval = 500
# Define starting points. It must satisfy A x = b.
x0 = [0.57186958424864897665429452899843;
4.9971472653643420613889247761108;
-1.3793445664086618762667058035731;
1.0403394252630473459930726676248;
-0.2300117084673765077695861691609]
# Solution
result = solve(p, x0)
println("status: ", result.status)
println("Solution: ", result.x_sol)
println("Satisfy Ax = b: ", A * result.x_sol ≈ b)
println("And inside bound constraints: ", all(-10.0 .<= result.x_sol .<= 10.0))Computing variable bounds for the reduced problem ...Done
1 9.066529
6 3.831687
17 2.956805
24 1.091294
37 0.499786
67 0.410615
80 0.271021
89 0.098408
96 0.094012
102 0.008064
113 0.001479
149 0.000721
160 0.000703
171 0.000149
173 0.000096
185 0.000065
191 0.000021
195 0.000003
199 0.000001
209 0
A termination criterion is reached: No termination (all). Mesh minimum precision reached (Algo)
Best feasible solutions: #32943 ( 0.087314 -0.0240172 -0.0243684 ) Evaluation OK f = 0 h = 0 (L2)
#38472 ( 0.064314 -0.0165172 -0.0166684 ) Evaluation OK f = 0 h = 0 (L2)
#38474 ( 0.090814 -0.0303172 -0.0305684 ) Evaluation OK f = 0 h = 0 (L2)
#39965 ( 0.083914 -0.0233172 -0.0232684 ) Evaluation OK f = 0 h = 0 (L2)
#39966 ( 0.077714 -0.0239172 -0.0244684 ) Evaluation OK f = 0 h = 0 (L2)
#39967 ( 0.079314 -0.0238172 -0.0241684 ) Evaluation OK f = 0 h = 0 (L2)
#40854 ( 0.082364 -0.0234672 -0.0235684 ) Evaluation OK f = 0 h = 0 (L2)
#41179 ( 0.08735 -0.0239902 -0.0244494 ) Evaluation OK f = 0 h = 0 (L2)
... A total of 104 feasible solutions were found.
Best infeasible solution: Undefined.
Blackbox evaluations: 442
Total model evaluations: 47868
Cache hits: 75
Total number of evaluations: 517
status: 1
Solution: [0.9964434715622237, 0.9453878715871625, 1.0387791045670747, 1.0097703430069789, 1.0096192092765588]
Satisfy Ax = b: true
And inside bound constraints: trueThe reader can take a look at the test folder for more complex examples.
Trade-offs for computational time performance
The default parameters of NOMAD.jl closely follow the default parameters of the NOMAD software. More importantly, NOMAD tries to find the best solution according to the maximum budget of evaluations provided by the user.
However, it happens that the user has a cheap blackbox in terms of computational time and needs a solution in a "short" amount of time. In this case, the user can remove the default quadratic model options. Generally, the computation of a given solution will be faster, (i.e. NOMAD will evaluate more points in a given amount of time) at a potential detriment of the solution quality.
Let illustrate it on the following problem.
using NOMAD
# Objective
function f(x)
return sqrt((x[1]-20)^2 + (x[2]-1)^2)
end
# Constraints
function c(x)
constraints = [sin(x[1]) - 1/10 - x[2], x[2] - sin(x[1])]
return constraints
end
# Evaluator
function bb(x)
bb_outputs = [f(x); c(x)]
success = true
count_eval = true
return (success, count_eval, bb_outputs)
end
p = NomadProblem(2, 3, ["OBJ"; "PB"; "PB"], bb,
lower_bound=[0.0;0.0],
upper_bound=[20.0;4.0])
# Set some options
p.options.display_degree = 2
p.options.max_bb_eval = 1500
# Default
println("This is the default")
@time result_default = solve(p, [0.0;0.0])
# Deactivate quadratic models
p.options.quad_model_search = false # for the search step ..
p.options.direction_type = "ORTHO N+1 NEG" # .. and the computation of the n+1 direction.
# One can also deactivate the sorting of poll directions by quadratic models but it is not
# mandatory as it plays a minimal role in the computational performance.
p.options.eval_queue_sort = "DIR_LAST_SUCCESS"
println("Now with no quadratic models")
@time result_with_no_quad_models = solve(p, [0.0;0.0])(bbo_sol = [16.88804, -0.1, -0.0], feasible = true, status = 1, x_sol = [3.1415925, 4.831242e-8])Note that the deactivation of quadratic models allows the solver to return a solution in a shorter time.
A good rule of thumb is to keep quadratic models if the blackbox possesses smoothness properties even if the derivatives are not available.
For more details about the parameters used in this section, we refer the reader to:
Multiobjective optimization
NOMAD.jl can solve multiobjective optimization problems of the form:
\[\begin{array}{rl} (MBB) \ \ \ \displaystyle \min_{x \in \mathbb{R}^n} & f(x) = \left(f_1(x), f_2(x), \ldots, f_m(x)\right)^\top\\ s.t. & c(x) \leq 0 \\ & lb \leq x \leq ub \\ \end{array}\]
where $m \geq 2$ is the number of objectives.
Unlike single-objective optimization, the set of solutions to such a problem is generally not a singleton. It is composed of several solutions whose objective values represent the best trade-offs that the algorithm may achieve at the end of the optimization.
NOMAD.jl switches to multiobjective mode when a NomadProblem is given more than one objective OBJ as parameter. However, NOMAD.jl only supports up to 5 objectives. Keep in mind that the larger the number of objectives, the more difficult the problem is to solve.
Finally, the multiobjective mode of NOMAD.jl does not support certain options available in single-objective optimization. Specifically, in multiobjective optimization, NOMAD.jl deactivates all search strategies except the speculative search, the NM search and the Quadratic model search. The default eval_sort_type and direction_type options are switched to "DIR_LAST_SUCCESS" and "ORTHO N+1 QUAD", respectively.
We consider the following example, the optimal design of a welded beam.
\[\begin{array}{rl} (WB) \ \ \ \displaystyle \min_{x = (h, l, t, b) \in \mathbb{R}^4} & \left(1.10471 h^2 l + 0.04811 t b (14.0 + l), \dfrac{2.1952}{t^3 b}\right)^\top \\ s.t. & c_1(x) = \tau(x) - 13600 \leq 0 \\ & c_2(x) = \sigma(x) - 30000 \leq 0\\ & c_3(x) = h - b \leq 0\\ & c_4(x) = 6000 - P_c(x) \leq 0\\ & h, b \in [0.125, 5], l, t \in [0.1, 10]\\ \end{array}\]
where $h, l, t, b$ are the four design parameters to optimize.
The following equations describe the physics constraints of the system (stress and buckling terms):
\[\begin{array}{rl} \tau(x) & = \sqrt{(\tau')^2 + (\tau'')^2 + \dfrac{l \tau' \tau''}{\sqrt{0.25 (l^2 + (h + t)^2))}}}, \\ \tau' & = \dfrac{6000}{\sqrt{2} h l}, \\ \tau'' & = \dfrac{6000 (14 + 0.5 l) \sqrt{0.25 (l^2 + (h + t)^2))}}{2 \left(0.707 h l (l^2/12 + 0.25 (h + t)^2)\right)},\\ \sigma(x) & = \dfrac{504000}{t^2 b},\\ P_c(x) & = 64746.022 (1 - 0.0282346 t) t b^3. \end{array}\]
\[f_1\]
corresponds to the construction costs and $f_2$ the end deflection (which corresponds to rigidity) of the beam. Both must be minimized.
using NOMAD
# blackbox
function welded_beam(x)
# Variables
h = x[1]
l = x[2]
t = x[3]
b = x[4]
# Objectives
f1 = 1.10471 * h^2 * l + 0.04811 * t * b * (14.0 + l)
f2 = 2.1952 / (t^3 * b)
# Physics equations
tau_p = 6000.0 / (sqrt(2) * h * l)
tau_pp = 6000 * (14 + 0.5 * l) * sqrt(0.25 * (l^2 + (h + t)^2))
tau_pp /= (2 * (0.707 * h * l * (l^2 / 12.0 + 0.25 * (h + t)^2)))
tau = sqrt(tau_p^2 + tau_pp^2 + l * tau_p * tau_pp / sqrt(0.25 * (l^2 + (h + t)^2)))
sigma = 504000 / (t^2 * b)
Pc = 64746.022 * (1 - 0.0282346 * t) * t * b^3
# Constraints
c1 = tau - 13600
c2 = sigma - 30000
c3 = h - b
c4 = 6000 - Pc
bb_outputs = [f1, f2, c1, c2, c3, c4]
success = true
count_eval = true
return (success, count_eval, bb_outputs)
end
# Define the problem
lb = [0.125, 0.1, 0.1, 0.125]
ub = [5.0, 10.0, 10.0, 5.0]
pb = NomadProblem(4, 6, ["OBJ", "OBJ", "PB", "PB", "PB", "PB"],
welded_beam,
lower_bound=[0.125, 0.1, 0.1, 0.125],
upper_bound=[5.0, 10.0, 10.0, 5.0])
# Set some options
pb.options.display_degree = 2
pb.options.max_bb_eval = 1500
# As for the single-objective case, you could deactivate this
# option to go faster, but the performance may be worse
# pb.options.quad_model_search = false
# Start from a middle box point
result = solve(pb, (lb + ub) / 2)
# The solution set is returned as a matrix of dimensions n x nb_solutions,
# where n is the dimension of the problem.
println("Optimization status: ", result.status)
println("The algorithm has found ", size(result.x_sol, 2), " solutions:")
for (ind, (x, bbo)) in enumerate(zip(eachcol(result.x_sol), eachcol(result.bbo_sol)))
println("Solution ", ind, ": ", x,
"; f(x) = ", bbo[1:2],
"; c(x) = ", bbo[3:end])
end
println("They are ", result.feasible ? "feasible" : "infeasible")
# In multiobjective optimization, it can be interesting to see the
# set of trade-offs in the objective space
using Plots
fig = scatter(result.bbo_sol[1, :], result.bbo_sol[2, :],
xlabel="Cost", ylabel="Deflection",
title="Pareto front approximation")
For more information about the multiobjective algorithm (DMulti-MADS), please refer to the following articles:
The problem is taken from: