R/BSASoptim.R
BSASoptim.RdYou could find more information about BSAS in https://arxiv.org/abs/1807.10470.
BSASoptim(fn, init = NULL, lower = c(-6, 0), upper = c(-1, 2), k = 5, constr = NULL, d0 = 0.001, d1 = 3, eta_d = 0.95, l0 = 0, l1 = 0, eta_l = 0.95, step = 0.8, eta_step = 0.95, n = 200, seed = NULL, trace = T, steptol = 0.01, p_min = 0.2, p_step = 0.2, n_flag = 2, pen = 1e+05)
| fn | objective function; function need to be optimized |
|---|---|
| init | default = NULL, it will generate randomly; Of course, you can specify it. |
| lower | lower of parameters to be estimated; Default = c(-6,0) because of the test on Michalewicz function of which thelower is c(-6,0); By the way, you should set one of init or lower parameter at least to make the code know the dimensionality of your problem. |
| upper | upper of parameters; Default = c(-1,2). |
| k | a positive integer.k is the number of beetles for exploring in every iteration. |
| constr | constraint function. For example, you can formulate \(x<=10\) as \(constr = function(x) return(x - 10)\). |
| d0 | a constant to gurantee that sensing length of antennae d doesn't equal to zero. More specifically, $$d^t = \eta_d * d^{t-1} + d_0$$where attenuation coefficient \(\eta_d\) belongs to \([0,1]\) |
| d1 | initial value of antenae length. You can specify it according to your problem scale |
| eta_d | attenuation coefficient of sensing length of antennae |
| l0 | position jitter factor constant.Default = 0. |
| l1 | initial position jitter factor.Default = 0.$$x = x - step * dir * sign(fn(left) - fn(right)) + l *random(npars)$$ |
| eta_l | attenuation coefficient of jitter factor.$$l^t = \eta_l * l^{t-1} + l_0$$ |
| step | initial step-size of beetle |
| eta_step | attenuation coefficient of step-size.$$step^t = \eta_step * step^{t-1}$$ |
| n | iterations times |
| seed | random seed; default = NULL ; The results of BAS depend on random init value and random directions.
Therefore, if you set a random seed, for example, |
| trace | default = T; trace the process of BAS iteration. |
| steptol | default = 0.01; Iteration will stop if step-size in current moment is less than steptol. |
| p_min | a constant belongs to [0,1]. If random generated value is lager than p_min and there are better positions in k beetles than current position, the next position of beetle will be the best position in k beetles. $$x_{best} = argmin(fn(x_i^t))$$ where \(i belongs [1,2,...,k]\) If random generated value is smaller than(<=) p_min and there are better positions in k beetles than current position, the next position of the beetle could be the random position within better positions. |
| p_step | a constant belongs to [0,1].If no better position in k beetles than current position, there is still a samll probability that step-size doesn't update. If you set a little k, you could set p_step slightly large. |
| n_flag | an positive integer; default = 2; If step-size doesn't update for successive n_flag times because p_step is larger than random generated value, the step-size will be forced updating. If you set a large p_step, set a small n_flag is suggested. |
| pen | penalty conefficient usually predefined as a large enough value, default 1e5 |
A list including best beetle position (parameters) and corresponding objective function value.
X. Y. Jiang, and S. Li, BAS: beetle antennae search algorithm for optimization problems, arXiv:1710.10724v1.
#======== examples start ======================= # >>>>>> example without constraint: Michalewicz function <<<<<< library(rBAS) mich <- function(x){ y1 <- -sin(x[1])*(sin((x[1]^2)/pi))^20 y2 <- -sin(x[2])*(sin((2*x[2]^2)/pi))^20 return(y1+y2) } result <- BSASoptim(fn = mich, lower = c(-6,0), upper = c(-1,2), seed = 1, n = 100,k=5,step = 0.6, trace = FALSE) result$par#> [1] -4.970202 1.578791result$value#> [1] -1.963534# >>>> example with constraint: Mixed integer nonlinear programming <<<< pressure_Vessel <- list( obj = function(x){ x1 <- floor(x[1])*0.0625 x2 <- floor(x[2])*0.0625 x3 <- x[3] x4 <- x[4] result <- 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 ) } ) result <- BSASoptim(fn = pressure_Vessel$obj, k = 10, lower =c( 1, 1, 10, 10), upper = c(100, 100, 200, 200), constr = pressure_Vessel$con, n = 200, step = 100, d1 = 4, pen = 1e6, steptol = 1e-6, n_flag = 2, seed = 2,trace = FALSE) result$par#> [1] 14.077163 7.087842 45.335198 140.284669result$value#> [1] 6090.567# >>>> example with constraint: Himmelblau function <<<< himmelblau <- list( obj = function(x){ x1 <- x[1] x3 <- x[3] x5 <- x[5] result <- 5.3578547*x3^2 + 0.8356891*x1*x5 + 37.29329*x[1] - 40792.141 }, con = function(x){ x1 <- x[1] x2 <- x[2] x3 <- x[3] x4 <- x[4] x5 <- x[5] g1 <- 85.334407 + 0.0056858*x2*x5 + 0.00026*x1*x4 - 0.0022053*x3*x5 g2 <- 80.51249 + 0.0071317*x2*x5 + 0.0029955*x1*x2 + 0.0021813*x3^2 g3 <- 9.300961 + 0.0047026*x3*x5 + 0.0012547*x1*x3 + 0.0019085*x3*x4 c( -g1, g1-92, 90-g2, g2 - 110, 20 - g3, g3 - 25 ) } ) result <- BSASoptim(fn = himmelblau$obj, k = 5, lower =c(78,33,27,27,27), upper = c(102,45,45,45,45), constr = himmelblau$con, n = 200, step = 100, d1 = 10, pen = 1e6, steptol = 1e-6, n_flag = 2, seed = 11,trace = FALSE) result$par # 78.01565 33.00000 27.07409 45.00000 44.95878#> [1] 78.01565 33.00000 27.07409 45.00000 44.95878result$value # -31024.17#> [1] -31024.17#======== examples end =======================