使用 R 最小化目标函数答案

【问题标题】：Minimise objective function using R使用 R 最小化目标函数
【发布时间】：2016-08-13 23:52:30
【问题描述】：

有人可以帮我解决这个问题，以在 R 中进行多元函数参数优化，我有一个这样的数据集。这只是数据的一个子集，完整数据集的维度是n type * m regions * 12 months。

Month   region  type physics maths allsub 
Jan     r1       1      4     5      9        
Feb     r1       1      3     8      11      
Mar     r1       1      5     4      9
Apr     r1       1      6     7      13
May     r1       1      4     4      8
Jun     r1       1      8     9      17
Jul     r1       1      4     3      7
Aug     r1       1      5     4      9
Sep     r1       1      3     8      11
Oct     r1       1      9     2      11
Nov     r1       1      4     7      11
Dec     r1       1      7     3      10
Jan     r1       2      5     8      13
Feb     r1       2      4     9      13
Mar     r1       2      8     3      11
Apr     r1       2      5     6      11
May     r1       2      6     4      10
Jun     r1       2      7     6      13
Jul     r1       2      3     7      10
Aug     r1       2      4     8      12
Sep     r1       2      4     4      8
Oct     r1       2      8     1      9
Nov     r1       2      2     3      5
Dec     r1       2      1     6      7

...     ...        ..  ...    ...   ....
...     ...        ..  ...    ...   ....

我还有一个数据集，每个区域的物理和数学学生人数最多。我的目标函数是100*(physics) + 65*(maths) >= 0。我想最小化这个功能，我的限制是 1. 物理和数学的总和应始终等于该地区和月份的 allsub。 2. 一个地区每月的物理学生总数应小于该地区物理学生的最大数量。 3. 一个地区每月的数学学生总数应少于该地区的最大数学学生人数。

我正在尝试使用 R。整个想法是在每个地区/类型/月份中找到合适数量的物理和数学学生，以最小化目标函数并满足约束条件。有人可以帮我解决这个问题吗？

EDIT ：根据 cmets 的要求。这是总容量数据集。数据框名称 = totalcap

   Month region physicscap mathscap
1    Jan   r1    9            13
2    Feb   r1    7            17
3    Mar   r1    13           7
4    Apr   r1    11           13
5    May   r1    10           8
6    Jun   r1    15           15
7    Jul   r1    7            10
8    Aug   r1    9            12
9    Sep   r1    7            12   
10   Oct   r1    17           3
11   Nov   r1    6            10
12   Dec   r1    8            9

这是我试过的脚本，

library(dplyr) 
library(MASS)
library(Rsolnp)

Month <- c('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec')
region <- c('r1')
physicscap <- c(5,5,8,6,7,9,5,6,4,10,5,8)
mathscap <- c(5,8,5,8,5,10,5,5,8,5,8,5)

totalcap <- data.frame(Month,region,physicscap,mathscap)
    #Constraints for the optimization.
constraints2 <- function(efforts){
  # constraints are:
  # 1. effort - allsub <= 0 in each region/month
  #     
  efforts$effort_calculated <- efforts$physics + efforts+maths
  reqeff <- summarise(group_by(efforts,region,Month),monthlyeffreg=sum(effort_calculated))
  reqeffallsub <- summarise(group_by(efforts,region,Month),allsubsum=sum(allsub))
  cons1 <- mutate(inner_join(reqeff,reqeffallsub,by=c('region'='region','Month'='Month'))
    ,diff=monthlyeffreg-allsubsum)
  constout <- cons1$diff


  # 2. sum(physics) - total physics available <= 0 in each region/month
  #
  phyreqeff <- summarise(group_by(efforts,region,Month),physicseff=sum(physics))
  cons2 <- mutate(inner_join(totalcap,phyreqeff,by=c('region'='region','Month'='Month')),
                   diff=physicseff-physicscap)
  constout <- c(constout,cons2$diff)


  # 3. sum(maths) - total maths available <= 0 in each region/month
  #
  matreqeff <- summarise(group_by(efforts,region,Month),mathseff=sum(maths))
  cons3 <- mutate(inner_join(totalcap,matreqeff,by=c('region'='region','Month'='Month')),
                   diff=mathseff-mathscap)
  constout <- c(constout,cons3$diff)
  constout
}


#Objective function to minimize the cost function.
objectivefunc <- function(efforts){
  nb_physics <- sum(efforts$physics)
  nb_maths <- sum(efforts$maths)
  objective <- (100*nb_physics + 55*nb_maths - 110)
  objective
}

Out2 <- solnp(pars = efforts,fun=objectivefunc,ineqfun=constraints2,ineqLB = rep(-100000,36), 
              ineqUB = rep(0,36), LB = rep(0,length(u)))

这是我得到的错误，

Error in p0/vscale[(neq + 2):(nc + np + 1)] : 
  non-numeric argument to binary operator

希望这可以解决 cmets 中的问题。我在这里尽力了我的水平，希望有人帮助我解决这个问题。

【问题讨论】：

SO 不是包搜索和代码编写服务。您应该 1) 搜索必要的 R 库，2) 然后发布所有编码函数的 R 代码，安装所需的库，加载它们，并表示所有约束。
这个问题似乎是重复的：stackoverflow.com/questions/5436630/… 尽管我不同意100*(physics) + 65*(maths) >= 0 可能是一个目标函数。这是一个逻辑表达式，只能取两个值之一。
感谢您的 cmets。你是对的，我已经完成了编码，编写了我的约束和目标函数。当我尝试使用 solnp 包运行它时，我收到了这个错误。 p0/vscale[(neq + 2):(nc + np + 1)] 中的错误：二进制运算符的非数字参数。所以我试图了解我是否做错了什么，看看是否有其他选择可以达到这个结果。我在这里给出的只是一个例子，所以不用太担心结构。再次感谢。
@42- 我已经更新了这个问题。希望它能回答你所有的问题。要求任何澄清，我可以进一步更新。
SO 将是一个更好的地方，如果您复制这里的工作，修复所有错误并通过库调用加载 pkg:Rsolnp 重新发布它，那么您将有更好的机会获得答案.我不记得 Rsolnp 是否仍需要外部软件包，但如果需要，您应该警告人们并给他们一个链接。我不是调试这个的人。

标签： r mathematical-optimization

【解决方案1】：

这是lpSolveAPI的一种方法：

dat <- data.frame(
    mon=rep(c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"),2),
    region="r1", 
    type=c(rep("1", 12), rep("2", 12)),
    physicsmin=1,
    mathsmin=1,
    allsub=c(9, 11, 9, 13, 8, 17, 7, 9, 11, 11, 11, 10, 13,13,11,11,10,13,10,12,8,9,5,7),
    stringsAsFactors=FALSE
)
dat
capdat <- data.frame(
    mon=c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"),
    region="r1",
    physicscap=c(9,7,13,11,10,15,7,9,7,17,6,8),
    mathscap=c(13,17,7,13,8,15,10,12,12,3,10,9),
    stringsAsFactors=FALSE
)
capdat

现在要解决每个月/地区组合的优化问题。因此我们将计算包装在一个函数中：

library(lpSolveAPI)
ntypes <- length(unique(dat[,"type"])) # number of types
typemap <- setNames(seq.int(ntypes), unique(dat[,"type"])) # map typename to 1,...,ntypes

solve_one <- function(subdat, capdat) {

    # create object
    lprec <- make.lp(0, ncol=2*ntypes) # for each type, two decision variables

    # By convention, we assume that the first ntypes variables are physics for type 1, ..., ntypes
    # and the second ntypes variables are maths

    # add objective and type
    set.objfn(lprec, obj=c(rep(100, ntypes), rep(65, ntypes))) 
    set.type(lprec, columns=seq.int(2*ntypes), type="integer") # no reals

    # add capacity constraints
    idx <- which(capdat[,"mon"]==subdat[1,"mon"] & capdat[,"region"]==subdat[1,"region"]) # lookup the right cap
    add.constraint(lprec, rep(1, ntypes), type="<=", rhs=capdat[idx,"physicscap"], indices=seq.int(ntypes))
    add.constraint(lprec, rep(1, ntypes), type="<=", rhs=capdat[idx,"mathscap"], indices=seq.int(ntypes+1, 2*ntypes))

    # add allsub equality constraints and minimum constraints
    for (typ in subdat[,"type"]) {
        add.constraint(lprec, c(1,1), type="=", rhs=subdat[typemap[typ], "allsub"], indices=c(typemap[typ], ntypes+typemap[typ]))
        add.constraint(lprec, 1, type=">=", rhs=subdat[typemap[typ],"physicsmin"], indices=typemap[typ])
        add.constraint(lprec, 1, type=">=", rhs=subdat[typemap[typ],"mathsmin"], indices=ntypes+typemap[typ])
    }

    # solution data.frame
    ans <- subdat[, c("mon", "region", "type")]

    # solve       
    if(solve(lprec)==0) {
        sol <- get.variables(lprec)
        for (i in seq.int(nrow(subdat))) {
            ans[i, "physics"] <- sol[typemap[subdat[i,"type"]]]
            ans[i, "maths"] <- sol[typemap[subdat[i,"type"]]+ntypes]
        }
    } else ans[,c("physics", "maths")] <- NA # no solution found

    return(ans)
}

现在我们将该函数应用于每个子数据集，其中包括每个月/区域组合的所有类型。我们在这里使用split/apply/combine 方法：

sp <- split(dat, list(dat[,"mon"], dat[,"region"]))
results <- lapply(sp, solve_one, capdat=capdat)
results <- do.call(rbind, results)
rownames(results) <- NULL
results

代码并不假定每个月/地区组合都存在所有类型（某些类型可能会被省略），但是如果同一月份/地区/类型组合存在多个条目，则解决方案将是错误的。（代码需要为此进行调整）。

【讨论】：

你是否考虑过physicscap和mathscap是区域级别的，而不是每个类型级别的。因此，约束一和二是区域中的 sum(physics) 应小于该区域中的物理上限。数学也是如此。
所以限制是每个月和每个地区，所有类型的物理学生的总和必须低于该地区/月的物理上限？（数学也是如此）
更新了答案，希望对您有所帮助。
我为此添加了一些额外的约束。 dat data.frame 现在还有两个额外的列来设置最小值。
很有趣。祝你有美好的一天。