方框約束可以設置爲目標函數本身的懲罰約束,如下所示。
lambda
乘數的值表示目標函數和約束之間的相對重要性。 對於lambda
的更高值,求解器將優先考慮目標函數的約束條件。
在這裏,我們解決對於x1^2 + X 2^2 = 25×2 < X1爲約束, 的整數解的簡單方程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)
}))
浮點表示的期望的輸出:
由於浮點表示,大多數情況下的結果不會完全等於您的預期整數輸出。 根據您的域名必須定義一個可接受的閾值,低於該閾值我們將數字視爲整數。
詳情參閱?.Machine
和 此R floating point precision
輸出收斂和驗證:
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
對於較低水平的λ的求解器忽略約束和與拉姆達增加求解器分配更大的權重到 約束相對於目標函數,因此約束得到滿足。
對於您的特定問題,您需要修改lambda,NP,itermax,F,CR
的值。
這工作完美。 – Patty