JDEoptim.RdA bespoke implementation of the ‘jDE’ variant by Brest et al. (2006) doi:10.1109/TEVC.2006.872133 .
JDEoptim(lower, upper, fn,
constr = NULL, meq = 0, eps = 1e-05,
NP = 10*length(lower), Fl = 0.1, Fu = 1,
tau_F = 0.1, tau_CR = 0.1, tau_pF = 0.1,
jitter_factor = 0.001,
tol = 1e-15, maxiter = 2000*length(lower), fnscale = 1,
compare_to = c("median", "max"),
add_to_init_pop = NULL,
trace = FALSE, triter = 1,
details = FALSE, ...)numeric vectors of lower and upper
bounds for the parameters to be optimized over. Must be finite
(is.finite) as they bound the hyper-rectangle
of the initial random population.
(nonlinear) objective function to be
minimized. It takes as first argument the vector of
parameters over which minimization is to take place. It must return
the value of the function at that point.
an optional function for specifying the
left-hand side of nonlinear constraints under which
we want to minimize fn.
Nonlinear equalities should be given first and defined to equal zero
(\(h_j(X) = 0\)), followed by nonlinear inequalities defined as
lesser than zero (\(g_i(X) \le 0\)).
This function takes the vector of parameters as its first argument
and returns a real vector with the length of the total number of
constraints. It defaults to NULL, meaning that
bound-constrained minimization is used.
an optional positive integer specifying that the first
meq constraints are treated as equality constraints,
and all the remaining as inequality constraints. Defaults to
0 (inequality constraints only).
maximal admissible constraint violation for equality constraints.
An optional real vector of small positive tolerance values with length
meq used in the transformation of equalities into inequalities
of the form \(|h_j(X)| - \epsilon \le 0\). A scalar value is expanded
to apply to all equality constraints. Default is 1e-5.
an optional positive integer giving the number of candidate
solutions in the randomly distributed initial population. Defaults to
10*length(lower).
an optional scalar which represents the minimum value that the
scaling factor F could take. Default is 0.1,
which is almost always satisfactory.
an optional scalar which represents the maximum value that
the scaling factor F could take. Default is 1,
which is almost always satisfactory.
an optional scalar which represents the probability that
the scaling factor F is updated. Defaults to 0.1,
which is almost always satisfactory.
an optional constant value which represents the probability
that the crossover probability CR is updated. Defaults to
0.1, which is almost always satisfactory.
an optional scalar which represents the probability that
the mutation probability \(p_F\) in the mutation strategy
DE/rand/1/either-or is updated. Defaults to 0.1.
an optional tuning constant for jitter.
If NULL only dither is used. Defaults to 0.001.
an optional positive scalar giving the tolerance for the
stopping criterion. Default is 1e-15.
an optional positive integer specifying the maximum
number of iterations that may be performed before the algorithm is halted.
Defaults to 2000*length(lower).
an optional positive scalar specifying the typical
magnitude of fn. It is used only in the stopping criterion.
Defaults to 1. See ‘Details’.
an optional character string controlling which function
should be applied to the fn values of the candidate solutions
in a generation to be compared with the so-far best one when evaluating
the stopping criterion. If “median” the median
function is used; else, if “max” the max function
is used. It defaults to “median”. See ‘Details’.
an optional real vector of length length(lower)
or matrix with length(lower) rows specifying
initial values of the parameters to be optimized which are appended to
the randomly generated initial population. It defaults to NULL.
an optional logical value indicating if a trace of the
iteration progress should be printed. Default is FALSE.
an optional positive integer giving the frequency of tracing
(every triter iterations) when trace = TRUE. Default is
triter = 1, in which case
iteration : < value of stopping test > ( value of best solution ) best solution { index of violated constraints }
is printed at each iteration.
an optional logical value. If TRUE the output
will contain the parameters in the final population and their
respective fn values. Defaults to FALSE.
optional additional arguments passed to fn
and constr.
The setting of the control parameters of canonical
Differential Evolution (DE) is crucial for the algorithm's performance.
Unfortunately, when the generally recommended values for these parameters
(see, e.g., Storn and Price, 1997) are unsuitable for use,
their determination is often difficult and time consuming.
The jDE algorithm proposed in Brest et al. (2006) employs a simple
self-adaptive scheme to perform the automatic setting of
control parameters scale factor F and crossover rate CR.
This implementation differs from the original description, most notably
in the use of the DE/rand/1/either-or mutation strategy
(Price et al., 2005), combination of jitter with dither
(Storn, 2008), and the random initialization of F and CR.
The mutation operator brings an additional control parameter, the
mutation probability \(p_F\), which is self-adapted in the
same manner as CR.
As done by jDE and its variants (Brest et al., 2021) each worse parent in the current population is immediately replaced (asynchronous update) by its newly generated better or equal offspring (Babu and Angira, 2006) instead of updating the current population with all the new solutions at the same time as in classical DE (synchronous update).
As the algorithm subsamples via sample()
which from R version 3.6.0 depends on RNGkind(*,
sample.kind), exact reproducibility of results from R versions
3.5.3 and earlier requires setting RNGversion("3.5.0").
In any case, do use set.seed() additionally for
reproducibility!
Constraint handling is done using the approach described in Zhang and Rangaiah (2012), but with a different reduction updating scheme for the constraint relaxation value (\(\mu\)). Instead of doing it once for every generation or iteration, the reduction is triggered for two cases when the constraints only contain inequalities. Firstly, every time a feasible solution is selected for replacement in the next generation by a new feasible trial candidate solution with a better objective function value. Secondly, whenever a current infeasible solution gets replaced by a feasible one. If the constraints include equalities, then the reduction is not triggered in this last case. This constitutes an original feature of the implementation.
The performance of any constraint handling technique for metaheuristics is severely impaired by a small feasible region. Therefore, equality constraints are particularly difficult to handle due to the tiny feasible region they define. So, instead of explicitly including all equality constraints in the formulation of the optimization problem, it might prove advantageous to eliminate some of them. This is done by expressing one variable \(x_k\) in terms of the remaining others for an equality constraint \(h_j(X) = 0\) where \(X = [x_1,\ldots,x_k,\ldots,x_d]\) is the vector of solutions, thereby obtaining a relationship as \(x_k = R_{k,j}([x_1,\ldots,x_{k-1},x_{k+1},\ldots,x_d])\). In this way both the variable \(x_k\) and the equality constraint \(h_j(X) = 0\) can be removed altogether from the original optimization formulation, since the value of \(x_k\) can be calculated during the search process by the relationship \(R_{k,j}\). Notice, however, that two additional inequalities $$l_k \le R_{k,j}([x_1,\ldots,x_{k-1},x_{k+1},\ldots,x_d]) \le u_k,$$ where the values \(l_k\) and \(u_k\) are the lower and upper bounds of \(x_k\), respectively, must be provided in order to obtain an equivalent formulation of the problem. For guidance and examples on applying this approach see Wu et al. (2015).
Bound constraints are enforced by the midpoint base approach (see, e.g., Biedrzycki et al., 2019).
Any DE variant is easily extended to deal with mixed integer
nonlinear programming problems using a small variation of the technique
presented by Lampinen and Zelinka (1999). Integer values are obtained by
means of the floor() function only in the evaluation
of the objective function and constraints, whereas DE itself still uses
continuous variables. Additionally, each upper bound of the integer
variables should be added by 1.
Notice that the final solution needs to be converted with floor()
to obtain its integer elements.
The algorithm is stopped if
$$\frac{\mathrm{compare\_to}\{[\mathrm{fn}(X_1),\ldots,\mathrm{fn}(X_\mathrm{npop})]\} - \mathrm{fn}(X_\mathrm{best})}{\mathrm{fnscale}} \le \mathrm{tol},$$
where the “best” individual \(X_\mathrm{best}\) is the
feasible solution with the lowest objective function value in the
population and the total number of elements in the population,
npop, is NP+NCOL(add_to_init_pop).
For compare_to = "max" this is the Diff criterion
studied by Zielinski and Laur (2008) among several other alternatives,
which was found to yield the best results.
A list with the following components:
The best set of parameters found.
The value of fn corresponding to par.
Number of iterations taken by the algorithm.
An integer code. 0 indicates successful completion.
1 indicates that the iteration limit maxiter
has been reached.
and if details = TRUE:
Matrix of dimension (length(lower), npop), with columns
corresponding to the parameter vectors remaining in the population.
The values of fn associated with poppar,
vector of length npop.
It is possible to perform a warm start, i.e., starting from the
previous run and resume optimization, using NP = 0 and the
component poppar for the add_to_init_pop argument.
Babu, B. V. and Angira, R. (2006) Modified differential evolution (MDE) for optimization of non-linear chemical processes. Computers and Chemical Engineering 30, 989–1002. doi:10.1016/j.compchemeng.2005.12.020 .
Biedrzycki, R., Arabas, J. and Jagodzinski, D. (2019) Bound constraints handling in differential evolution: An experimental study. Swarm and Evolutionary Computation 50, 100453. doi:10.1016/j.swevo.2018.10.004 .
Brest, J., Greiner, S., Boskovic, B., Mernik, M. and Zumer, V. (2006) Self-adapting control parameters in differential evolution: A comparative study on numerical benchmark problems. IEEE Transactions on Evolutionary Computation 10, 646–657. doi:10.1109/TEVC.2006.872133 .
Brest, J., Maucec, M. S. and Boskovic, B. (2021) Self-adaptive differential evolution algorithm with population size reduction for single objective bound-constrained optimization: Algorithm j21; in 2021 IEEE Congress on Evolutionary Computation (CEC). IEEE, pp. 817–824. doi:10.1109/CEC45853.2021.9504782 .
Lampinen, J. and Zelinka, I. (1999). Mechanical engineering design optimization by differential evolution; in Corne, D., Dorigo, M. and Glover, F., Eds., New Ideas in Optimization. McGraw-Hill, pp. 127–146.
Price, K. V., Storn, R. M. and Lampinen, J. A. (2005) Differential evolution: A practical approach to global optimization. Springer, Berlin, Heidelberg, pp. 117–118. doi:10.1007/3-540-31306-0_2 .
Storn, R. (2008) Differential evolution research — Trends and open questions; in Chakraborty, U. K., Ed., Advances in differential evolution. SCI 143, Springer, Berlin, Heidelberg, pp. 11–12. doi:10.1007/978-3-540-68830-3_1 .
Storn, R. and Price, K. (1997) Differential evolution - A simple and efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization 11, 341–359. doi:10.1023/A:1008202821328 .
Wu, G., Pedrycz, W., Suganthan, P. N. and Mallipeddi, R. (2015) A variable reduction strategy for evolutionary algorithms handling equality constraints. Applied Soft Computing 37, 774–786. doi:10.1016/j.asoc.2015.09.007 .
Zhang, H. and Rangaiah, G. P. (2012) An efficient constraint handling method with integrated differential evolution for numerical and engineering optimization. Computers and Chemical Engineering 37, 74–88. doi:10.1016/j.compchemeng.2011.09.018 .
Zielinski, K. and Laur, R. (2008) Stopping criteria for differential evolution in constrained single-objective optimization; in Chakraborty, U. K., Ed., Advances in differential evolution. SCI 143, Springer, Berlin, Heidelberg, pp. 111–138. doi:10.1007/978-3-540-68830-3_4 .
Function DEoptim() in the DEoptim package
has many more options than JDEoptim(), but does not allow constraints
in the same flexible manner.
# \donttest{
# NOTE: Examples were excluded from testing
# to reduce package check time.
# Use a preset seed so test values are reproducible.
set.seed(1234)
# Bound-constrained optimization
# Griewank function
#
# -600 <= xi <= 600, i = {1, 2, ..., n}
# The function has a global minimum located at
# x* = (0, 0, ..., 0) with f(x*) = 0. Number of local minima
# for arbitrary n is unknown, but in the two dimensional case
# there are some 500 local minima.
#
# Source:
# Ali, M. Montaz, Khompatraporn, Charoenchai, and
# Zabinsky, Zelda B. (2005).
# A numerical evaluation of several stochastic algorithms
# on selected continuous global optimization test problems.
# Journal of Global Optimization 31, 635-672.
# https://doi.org/10.1007/s10898-004-9972-2
griewank <- function(x) {
1 + crossprod(x)/4000 - prod( cos(x/sqrt(seq_along(x))) )
}
JDEoptim(rep(-600, 10), rep(600, 10), griewank,
tol = 1e-7, trace = TRUE, triter = 50)
#> 50 : < 0.3263563 > ( 1.122444 ) -6.729448 -10.02108 -11.12805 -4.840487
#> 2.996272 0.6737868 -17.28767 -1.521765 -1.649234 6.854039
#> 100 : < 0.712248 > ( 0.2196246 ) 0.1282567 4.383256 -0.2021722 -0.1593782
#> 6.069386 0.2334851 -0.5846669 0.2491573 0.2358098 1.248203
#> 150 : < 0.3593402 > ( 0.1623947 ) -3.2292 -4.379734 0.07603338 0.2827644
#> -0.4572112 -0.7278611 0.1326614 0.161805 -1.176268 -0.2898208
#> 200 : < 0.2351282 > ( 0.04224012 ) -2.972997 0.02048603 0.01127948 0.1631893
#> 6.995231 0.1845278 0.08342469 0.241416 -0.07937427 -0.2440344
#> 250 : < 0.1098993 > ( 0.02653812 ) 3.146412 -4.340857 0.0282541 -0.2611893
#> -0.06540064 0.2789499 0.007484945 0.03182706 0.1376961 0.05464311
#> 300 : < 0.07233788 > ( 0.0082436 ) -0.01934878 0.03627711 0.0241738 0.228187
#> 0.092197 0.04312658 -0.01267154 -0.01943716 0.01010444 0.03959774
#> 350 : < 0.0375359 > ( 0.001272347 ) -0.01340587 0.004539416 -0.02032376
#> 0.006505448 0.03871329 0.01098626 0.007237393 -0.003564126 0.1295924
#> -0.00536604
#> 400 : < 0.01793625 > ( 0.0003173216 ) -0.01340587 0.004539416 -0.004244852
#> 0.006505448 0.03871329 0.01098626 0.007237393 0.001893758 0.02939837
#> -0.00536604
#> 450 : < 0.0004802076 > ( 7.09431e-05 ) 0.00359858 -0.008282913 -0.01293652
#> -0.005235058 0.009704593 0.004982956 0.00563215 -9.934005e-05 -0.006192586
#> 8.710251e-05
#> 500 : < 1.166464e-06 > ( 1.216925e-07 ) -5.232293e-05 -5.227759e-05
#> 0.0002199203 -0.0006180909 -0.0003857748 -0.0002014436 0.0001286359
#> -0.0001717095 -5.556345e-05 0.0009163679
#> $par
#> [1] 5.843883e-05 3.517630e-05 -8.215326e-05 -5.042123e-05 -1.595923e-04
#> [6] -1.260476e-04 7.304026e-05 -2.265423e-04 -5.970416e-05 -3.394538e-04
#>
#> $value
#> [1] 1.693632e-08
#>
#> $iter
#> [1] 518
#>
#> $convergence
#> [1] 0
#>
# Nonlinear constrained optimization
# 0 <= x1 <= 34, 0 <= x2 <= 17, 100 <= x3 <= 300
# The global optimum is
# (x1, x2, x3; f) = (0, 16.666667, 100; 189.311627).
#
# Source:
# Westerberg, Arthur W., and Shah, Jigar V. (1978).
# Assuring a global optimum by the use of an upper bound
# on the lower (dual) bound.
# Computers and Chemical Engineering 2, 83-92.
# https://doi.org/10.1016/0098-1354(78)80012-X
fcn <-
list(obj = function(x) {
35*x[1]^0.6 + 35*x[2]^0.6
},
eq = 2,
con = function(x) {
x1 <- x[1]; x3 <- x[3]
c(600*x1 - 50*x3 - x1*x3 + 5000,
600*x[2] + 50*x3 - 15000)
})
JDEoptim(c(0, 0, 100), c(34, 17, 300),
fn = fcn$obj, constr = fcn$con, meq = fcn$eq,
tol = 1e-7, trace = TRUE, triter = 50)
#> 50 : < 35.81246 > ( 190.0442 ) 0.03673432 16.07114 105.4244 { 1 } { 2 }
#> 100 : < 0.2789816 > ( 189.4313 ) 0.000181572 16.65499 100.0515 { 1 } { 2 }
#> 150 : < 0.1677861 > ( 189.4071 ) 5.266726e-05 16.66677 100.0006 { 1 } { 2 }
#> 200 : < 0.1441106 > ( 189.328 ) 2.801869e-06 16.66668 100.0001 { 1 } { 2 }
#> 250 : < 0.04505788 > ( 189.3337 ) 4.648397e-06 16.66667 100 { 1 } { 2 }
#> 300 : < 0.02396127 > ( 189.3138 ) 9.490506e-08 16.66667 100 { 1 } { 2 }
#> 350 : < 0.004383328 > ( 189.3132 ) 5.463639e-08 16.66667 100 { 1 } { 2 }
#> 400 : < 6.234711e-06 > ( 189.3116 ) 9.364951e-13 16.66667 100 { }
#> $par
#> [1] 1.129659e-15 1.666667e+01 1.000000e+02
#>
#> $value
#> [1] 189.3116
#>
#> $iter
#> [1] 418
#>
#> $convergence
#> [1] 0
#>
# Designing a pressure vessel
# Case A: all variables are treated as continuous
#
# 1.1 <= x1 <= 12.5*, 0.6 <= x2 <= 12.5*,
# 0.0 <= x3 <= 240.0*, 0.0 <= x4 <= 240.0
# Roughly guessed*
# The global optimum is (x1, x2, x3, x4; f) =
# (1.100000, 0.600000, 56.99482, 51.00125; 7019.031).
#
# Source:
# Lampinen, Jouni, and Zelinka, Ivan (1999).
# Mechanical engineering design optimization
# by differential evolution.
# In: David Corne, Marco Dorigo and Fred Glover (Editors),
# New Ideas in Optimization, McGraw-Hill, pp 127-146
pressure_vessel_A <-
list(obj = function(x) {
x1 <- x[1]; x2 <- x[2]; x3 <- x[3]; x4 <- x[4]
0.6224*x1*x3*x4 + 1.7781*x2*x3^2 +
3.1611*x1^2*x4 + 19.84*x1^2*x3
},
con = function(x) {
x1 <- x[1]; x2 <- x[2]; x3 <- x[3]; x4 <- x[4]
c(0.0193*x3 - x1,
0.00954*x3 - x2,
750.0*1728.0 - pi*x3^2*x4 - 4/3*pi*x3^3)
})
JDEoptim(c( 1.1, 0.6, 0.0, 0.0),
c(12.5, 12.5, 240.0, 240.0),
fn = pressure_vessel_A$obj,
constr = pressure_vessel_A$con,
tol = 1e-7, trace = TRUE, triter = 50)
#> 50 : < 1044.137 > ( 7352.375 ) 1.188077 0.6182044 61.1726 30.70314 { }
#> 100 : < 178.2174 > ( 7036.841 ) 1.101541 0.6001186 56.8518 51.83681 { }
#> 150 : < 40.03357 > ( 7028.614 ) 1.101994 0.6000808 57.07092 50.67729 { }
#> 200 : < 8.571083 > ( 7021.713 ) 1.100067 0.6000588 56.972 51.1448 { }
#> 250 : < 2.485194 > ( 7019.718 ) 1.100059 0.6000004 56.98936 51.03379 { }
#> 300 : < 0.5196705 > ( 7019.193 ) 1.100022 0.6000033 56.99455 51.00328 { }
#> 350 : < 0.07812566 > ( 7019.131 ) 1.100068 0.6000009 56.99826 50.98161 { }
#> 400 : < 0.01559852 > ( 7019.034 ) 1.1 0.6 56.99483 51.00124 { }
#> 450 : < 0.001081023 > ( 7019.032 ) 1.1 0.6 56.99483 51.0012 { }
#> 500 : < 0.0002246616 > ( 7019.031 ) 1.1 0.6 56.99482 51.00123 { }
#> 550 : < 2.388982e-05 > ( 7019.031 ) 1.1 0.6 56.99482 51.00125 { }
#> 600 : < 2.414554e-06 > ( 7019.031 ) 1.1 0.6 56.99482 51.00125 { }
#> 650 : < 5.135362e-07 > ( 7019.031 ) 1.1 0.6 56.99482 51.00125 { }
#> 700 : < 1.097051e-07 > ( 7019.031 ) 1.1 0.6 56.99482 51.00125 { }
#> $par
#> [1] 1.10000 0.60000 56.99482 51.00125
#>
#> $value
#> [1] 7019.031
#>
#> $iter
#> [1] 701
#>
#> $convergence
#> [1] 0
#>
# Mixed integer nonlinear programming
# Designing a pressure vessel
# Case B: solved according to the original problem statements
# steel plate available in thicknesses multiple
# of 0.0625 inch
#
# wall thickness of the
# shell 1.1 [18*0.0625] <= x1 <= 12.5 [200*0.0625]
# heads 0.6 [10*0.0625] <= x2 <= 12.5 [200*0.0625]
# 0.0 <= x3 <= 240.0, 0.0 <= x4 <= 240.0
# The global optimum is (x1, x2, x3, x4; f) =
# (1.125 [18*0.0625], 0.625 [10*0.0625],
# 58.29016, 43.69266; 7197.729).
pressure_vessel_B <-
list(obj = function(x) {
x1 <- floor(x[1])*0.0625
x2 <- floor(x[2])*0.0625
x3 <- x[3]; x4 <- x[4]
0.6224*x1*x3*x4 + 1.7781*x2*x3^2 +
3.1611*x1^2*x4 + 19.84*x1^2*x3
},
con = function(x) {
x1 <- floor(x[1])*0.0625
x2 <- floor(x[2])*0.0625
x3 <- x[3]; x4 <- x[4]
c(0.0193*x3 - x1,
0.00954*x3 - x2,
750.0*1728.0 - pi*x3^2*x4 - 4/3*pi*x3^3)
})
res <- JDEoptim(c( 18, 10, 0.0, 0.0),
c(200+1, 200+1, 240.0, 240.0),
fn = pressure_vessel_B$obj,
constr = pressure_vessel_B$con,
tol = 1e-7, trace = TRUE, triter = 50)
#> 50 : < 483.7403 > ( 7394.642 ) 18.29861 10.4097 56.8351 54.28991 { }
#> 100 : < 106.9966 > ( 7234.642 ) 18.98648 10.69023 57.78511 46.62088 { }
#> 150 : < 37.9038 > ( 7208.274 ) 18.18732 10.13377 58.14807 44.51665 { }
#> 200 : < 6.567386 > ( 7199.978 ) 18.35591 10.93516 58.28217 43.77585 { }
#> 250 : < 0.9720159 > ( 7198.407 ) 18.16936 10.31089 58.28979 43.7093 { }
#> 300 : < 0.2751255 > ( 7197.775 ) 18.25372 10.45631 58.28951 43.69634 { }
#> 350 : < 0.057912 > ( 7197.74 ) 18.19731 10.44728 58.29004 43.69339 { }
#> 400 : < 0.009539295 > ( 7197.735 ) 18.32325 10.37728 58.29014 43.69286 { }
#> 450 : < 0.001366935 > ( 7197.729 ) 18.27652 10.6077 58.29015 43.69266 { }
#> 500 : < 0.0001244909 > ( 7197.729 ) 18.56162 10.39549 58.29015 43.69266 { }
#> 550 : < 4.673028e-05 > ( 7197.729 ) 18.08508 10.39535 58.29016 43.69266 { }
#> 600 : < 7.282225e-06 > ( 7197.729 ) 18.46083 10.80006 58.29016 43.69266 { }
#> 650 : < 1.267306e-06 > ( 7197.729 ) 18.91071 10.14678 58.29016 43.69266 { }
#> 700 : < 2.536872e-07 > ( 7197.729 ) 18.27851 10.28549 58.29016 43.69266 { }
res
#> $par
#> [1] 18.05619 10.19147 58.29016 43.69266
#>
#> $value
#> [1] 7197.729
#>
#> $iter
#> [1] 724
#>
#> $convergence
#> [1] 0
#>
# Now convert to integer x1 and x2
c(floor(res$par[1:2]), res$par[3:4])
#> [1] 18.00000 10.00000 58.29016 43.69266
# }