The box constraints could be setup as penalty constraints to the objective function itself as below.
The value of lambda multiplier dictates the relative importance between objective function and constraints.
With higher values of lambda the solver will prioritize constraints over the objective function.
Here we solve a simple equation for integer solutions of x1^2 + x2^2 = 25 with x2 < x1 as constraint,
with expected output of c(4,3)
library("DEoptim")
#run DEoptim and set a seed first for replicability
#note the results are sensitive to seed values and parameters (lambda,NP,itermax,F,CR)
set.seed(1234)
#create a vector/grid of lambda values
lambdaVec = sapply(0:12,function(x) 10^x)
lambdaVec
#[1] 1e+00 1e+01 1e+02 1e+03 1e+04 1e+05 1e+06 1e+07 1e+08 1e+09 1e+10
#For each value of lambda compute the output of optimization function and combine the results
optimSummary = do.call(rbind,lapply(lambdaVec, function(lambdaParam) {
fnCustom = function(x,lambda=lambdaParam) {
x1 <- x[1]
x2 <- x[2]
#integer param constraints
firstIntPenalty <- (x1-round(x1))^2
secondIntPenalty <- (x2-round(x2))^2
# x2 < x1, note the sign is opposite of desired constraint
inEqualityConstraint <- sum(x2>x1)
100 * (x1^2 + x2^2 - 25)^2 + lambda * ( firstIntPenalty + secondIntPenalty + inEqualityConstraint )
}
lower <- c(0,0)
upper <- c(5,5)
#you will have to tinker with values of NP,F and CR and monitor convergence, see ?DEoptim last paragraph
outDEoptim <- DEoptim(fnCustom, lower, upper, DEoptim.control(NP = 80, itermax = 1000, F = 1.2, CR = 0.7,trace=TRUE))
#output data.frame of optimization result
optimRes <- data.frame(lambda = lambdaParam ,param1 = outDEoptim$optim$bestmem[1],param2 = outDEoptim$optim$bestmem[2])
rownames(optimRes) <- NULL
return(optimRes)
} ))
Floating Point Representation:
Due to floating point represenation the results in most cases will not to be exactly equal to your expected integer output.
Depending on your domain will have to define an acceptable threshold below which we consider the numbers as integer.
For more details refer to ?.Machine and
this R floating point precision
Output Convergence and verification:
threshold = 1e-6
expectedOut = c(4,3)
#optimization summary
optimSummary
# lambda param1 param2
#1 1e+00 4.999996 0.0002930537
#2 1e+01 4.000000 3.0000000000
#3 1e+02 4.000000 3.0000000000
#4 1e+03 4.000000 3.0000000000
#5 1e+04 4.000000 3.0000000000
#6 1e+05 4.000000 3.0000000000
#7 1e+06 4.000000 3.0000000000
#8 1e+07 4.000000 3.0000000000
#9 1e+08 4.000000 2.9999999962
#10 1e+09 3.999999 2.9999998843
#11 1e+10 4.000000 2.9999999569
#12 1e+11 4.000000 3.0000000140
#13 1e+12 4.000000 3.0000000194
#Verify output
#1)With constraintBreach1 and constraintBreach2 we check if difference between output and expected result
#has breached threshold
#2)With constraintBreach3 we check if x1 > x2 condition is violated
#3)Columns with TRUE observations indicate breach of respective constraints at particular lambda values
verifyDF = data.frame(lambdaVec,constraintBreach1 = abs(optimSummary$param1-expectedOut[1]) > threshold
, constraintBreach2 = abs(optimSummary$param2-expectedOut[2]) > threshold
, constraintBreach3 = optimSummary$param1 < optimSummary$param1)
verifyDF
# lambdaVec constraintBreach1 constraintBreach2 constraintBreach3
#1 1e+00 TRUE TRUE FALSE
#2 1e+01 FALSE FALSE FALSE
#3 1e+02 FALSE FALSE FALSE
#4 1e+03 FALSE FALSE FALSE
#5 1e+04 FALSE FALSE FALSE
#6 1e+05 FALSE FALSE FALSE
#7 1e+06 FALSE FALSE FALSE
#8 1e+07 FALSE FALSE FALSE
#9 1e+08 FALSE FALSE FALSE
#10 1e+09 FALSE FALSE FALSE
#11 1e+10 FALSE FALSE FALSE
#12 1e+11 FALSE FALSE FALSE
#13 1e+12 FALSE FALSE FALSE
For lower levels of lambda the solver ignores constraints and with increasing lambda the solver assigns more weight to
constraints relative to objective function and hence the constraints are satisfied.
You will have to tinker with values of lambda,NP,itermax,F,CR for your particular problem.
