我想使用库mco中的函数mco来解决一个多目标问题并找到Pareto,但是我无法正确地设置约束。
目标函数如下。问题的背景是项目选择,即我有五个项目由x1、x2、.可以选择x5,并且只能选择一些。例如,如果选择了项目1,那么x1=1 (如果没有选择x1=0 ),对于所有项目( xn的值是离散的,无论是1还是0),这都是正确的。另一个限制因素是,选定项目的总预算应少于100项。运行nsga2函数后,Solution中的参数似乎不正确,因为参数不是1或0。我的约束有错吗?如何找到x1到x5的最佳值?谢谢!
# objective functions to minimize
ObjFun <- function (x){
f1 <- -0.02*x[1] + 0.01*x[2] + 0.02*x[3] + -0.01*x[4] + 0.02*x[5]
f2 <- 0.17*x[1] + -0.08*x[2] + 0.10*x[3] + 0.09*x[4] + 0.07*x[5]
c(f1, f2) }
# The constraints
Constr <- function(x){
100 >= 20*x[1] + 30*x[2] + 20*x[3] + 33*x[4] + 60*x[5] # Total budget >= total project costs
x[1:5] == 1
x[1:5] == 0 }
library(mco)
Solution <- nsga2(ObjFun, 5, 2, lower.bounds=c(0,0,0,0,0), upper.bounds=c(1,1,1,1,1), constraints = Constr)
# plot(Solution)
Solution$par发布于 2016-08-09 23:23:15
由于x[i]只能是1或0,所以您正在处理一个组合优化问题,其中必须对其进行优化的空间是离散的:
通常,数值优化程序是在连续空间(R^n的子集)上构造的。然而,在您的例子中,离散空间是小的,问题本身是一个简单的蛮力方法,在所有32个可能的点上对ObjFunc进行评估。帕累托边界也是离散的。
## objective functions to minimize
ObjFun <- function (x){
f1 <- -0.02*x[1] + 0.01*x[2] + 0.02*x[3] + -0.01*x[4] + 0.02*x[5]
f2 <- 0.17*x[1] + -0.08*x[2] + 0.10*x[3] + 0.09*x[4] + 0.07*x[5]
c(f1=f1, f2=f2)
}
## space of all 32 feasible solutions
space <- expand.grid(data.frame(matrix(0:1, nrow=2, ncol=5)))
## brute force evaluation of ObjFun on all the 32 feasible solutions
val <- sapply(data.frame(t(space)), ObjFun)
tmp <- sol <- cbind(space, t(val))
## returns indices of all rows which are Pareto dominated
## by the i-th row
which.are.dominated <- function(i, tmp){
s1 <- tmp$f1[i]
s2 <- tmp$f2[i]
with(tmp,
which( (s1 <= f1) &
(s2 <= f2) &
( (s1 < f1) |
(s2 < f2) )
))
}
## For each feasible solution i, remove all feasible solutions which are Pareto dominated by feasible solutions i
i <- 1
repeat{
remove <- which.are.dominated(i, tmp)
if(length(remove)>0) tmp <- tmp[-remove, ]
if(i>=nrow(tmp)) break
i <- i+1
}
with(sol, plot(f1, f2))
points(tmp$f1, tmp$f2, pch=20, col=2)
legend("topright", col=2, pch=20, "Pareto frontier")参考文献:
从几年前开始使用R.我可能第一次使用repeat语句.
编辑:一种非暴力的方法是使用nsga2 :D,当我设置它时,在n维立方体0,1^n中,n是项目的数目,搜索x的解决方案;该算法产生许多解决方案(在我的示例中为200),然后您可以使用round将其“离散化”到0或1。对于更多的项目,要获得更精确的接近帕累托边界,你必须使用更多的世代(例如600)。在最后的地块中,如果考虑了12个以上的项目,则只绘制一个成本样本。
##n.projects <- 12
n.projects <- 50
if(n.projects>25) generations=600
set.seed(1)
vecf1 <- rnorm(n.projects)
vecf2 <- rnorm(n.projects)
vcost <- rnorm(n.projects)
n.solutions <- 200
library(mco)
ObjFun <- function (x){
f1 <- sum(vecf1*x)
f2 <- sum(vecf2*x)
c(f1=f1, f2=f2)
}
Constr <- function(x){
c(100 - sum(vcost*x)) # Total budget >= total project costs
}
Solution <- nsga2(ObjFun, n.projects, 2,
lower.bounds=rep(0,n.projects), upper.bounds=rep(1,n.projects),
popsize=n.solutions, constraints = Constr, cdim=1,
generations=generations)
selected.project.combinations <- unique(round(Solution$par))
selected.project.combinations.costs <- sapply(data.frame(t(selected.project.combinations)), ObjFun)
## final plotting of results
max.n.proj.plot <- 12
if(n.projects <= max.n.proj.plot){
xsamp <- expand.grid(data.frame(matrix(0:1, nrow=2, ncol=n.projects)))
}else{
xsamp <- matrix(sample(0:1, n.projects*2^max.n.proj.plot, replace=TRUE), ncol=n.projects)
}
fsamp <- sapply(data.frame(t(xsamp)), ObjFun)
par(mfrow=c(1,2))
plot(Solution)
points(fsamp[1, ], fsamp[2, ])
points(t(selected.project.combinations.costs), col=3, pch=20)
legend("bottomleft", bty="n", pch=c(20,1), col=c(3,1),
c("Costs of optimal\nproject combinations",
"Costs of discarded\nproject combinations"),
y.intersp=1.8
)
plot(t(fsamp), xlim=range(Solution$value[ ,1], fsamp[1, ]),
ylim=range(Solution$value[ ,2], fsamp[2, ]))
points(Solution$value, col=2, pch=".")
points(t(selected.project.combinations.costs), col=3, pch=20)https://stackoverflow.com/questions/38860287
复制相似问题