Solving a Global Optimization problem using Differential Evolutionary Algorithm using R

Question

I need to determine the global optimum results of this objective function. I define the problem by minimizing the squared difference as represented in function $f(q_1,q_2,\alpha_1,\alpha_2)$

The objective function: $$f(q_1,q_2,\alpha_1,\alpha_2)= \frac{1}{5}\sum\limits_{i = 1}^{5} (b_{i} - b_{mod_{i}})^2$$

With bounds:

lower = c(1000, 1000, 0.8, 1.4) upper = c(1e10, 1e35, 1.5, 15.0)

where, $$b_{mod,i} = q_{1} \lambda_i^{-\alpha_1} + q_{2} \lambda_i^{-\alpha_2}$$

and $\lambda_i = c(375, 470, 528, 625, 880)$ are the wavelengths of 5 channels and $b_i = c(228,124,98,67,44)$ represents the light absorption measurement corresponding to each wavelength.

Here $b_i$ is a measured value from a device and $b_{mod,i}$ is the modelled value achieved by the optimization process.

Approach 1: Excel Solver

In the Excel solver, I used the non-linear solver function with the inbuilt "Evolutionary" algorithm. The optimized values $q_1,\alpha_1,q_2,\alpha_2$ seems reasonable.

Approach 2: in R using DEoptim

I found that the Evolutionary algorithm in Excel is based on the Differential Evolutionary algorithm and in R there is a package called DEoptim (Link). Also, I need to run the optimization process for a large data set where each timestamp of data represents $b$. However, DEoptim results were found to be very close to the upper bound and often are very sensitive to the bounds provided.

Queries

I am trying to use this technique in R to solve a physical problem.

• Is DEoptim a good choice for the defined problem?
• Are there any other optimization techniques to solve the present problem?

Answer 1

(This is too long for a comment.)

What I meant with my comment is this: given two candidate solutions to an optimization model, it is easy to figure out which is better. Insert each solution into the objective function; the solution with the lower objective function value is better (as you minimize).

As pointed out by Enrique Gabriel Baquela in the comments, DE and similar methods require certain parameters to be set. Default parameters are rarely good choices, though in the case of DE, they often work remarkably robustly. But anyway, that leaves the number of iterations to be set. (See http://enricoschumann.net/R/remarks.htm for a discussion.)

In your model, the two terms can be very similar. If the alphas are around 1.4, they can be identical. That means there is an identification problem; the objective function will be flat. DE can handle such cases, but typically requires a large population to do so efficiently (much larger than the 5-10 times the number of parameters, as is often a default choice).

What is more, the parameter values can be quite extreme. For instance, if $q_2$ is at $10^{30}$, say, the first term will vanish under standard computer arithmetic, because of rounding error.

If you fixed the alpha parameters, then the model could be solved with Least Squares (if I understand it correctly). That suggests an alternative solution strategy: run a grid search over $\alpha_1$ and $\alpha_2$, and for each vector of values compute the $q$ parameters. Here is a sketch in R; for the grid search, I use the implementation in package NMOF (which I maintain):

lambda <- c(375, 470, 528, 625, 880)
b <- c(228, 124, 98, 67, 44)

par <- lower <- c(1000, 1000, 0.8, 1.4)
par <- upper <- c(1e10, 1e35, 1.5, 15.0)

## original objective function
fun <- function(par, b, lambda) {
    res <- par[1]*lambda^(-par[3]) + par[2]*lambda^(-par[4])  - b
    sum(res*res)
}

## LS objective function
fun.ls <- function(par, b, lambda) {
    QR <- qr(cbind(lambda^(-par[1]), lambda^(-par[2])))
    coef <- qr.solve(QR, b)

    if (any(coef < lower[1:2] |  coef > upper[1:2]))
        ans <- 1e6  ## mark infeasible solutions with a large value
    else {
        ans <- qr.resid(QR, b)
        ans <- sum(ans*ans)
    }
    attr(ans, "coefficients") <- coef
    ans
}

## run a grid search
par <- c(0.8, 4.010101)
library("NMOF")
gs <- gridSearch(fun.ls,
                 lambda = lambda,
                 b = b,
                 lower = c(0.8, 1.4),
                 upper = c(1.5, 15),
                 n = 100)
## $minfun
## [1] 12.68621
## 
## $minlevels
## [1] 0.800000 4.010101

fun.ls(gs$minlevels, b, lambda)
## [1] 12.68621
## attr(,"coefficients")
## [1] 8.598193e+03 3.212813e+12

par <- c(8598.19, 3212813450654, 0.8, 4.010101)
fun(par, b, lambda)
## [1] 12.68621

data.frame(b,
           fit = par[1]*lambda^(-par[3]) + par[2]*lambda^(-par[4]))
##     b       fit
## 1 228 228.04507
## 2 124 124.49555
## 3  98  95.85731
## 4  67  69.58422
## 5  44  42.91865