R中的优化：具有二进制调度变量的成本函数？答案

【问题标题】：Optimization in R: cost function with binary scheduling variables?R中的优化：具有二进制调度变量的成本函数？
【发布时间】：2019-04-30 01:04:39
【问题描述】：

下面详细介绍了我无法解决的优化问题的简化版本。

目标是最小化通过卡车运送水的组织的成本函数，并使用该等式生成最小化成本的卡车交付时间表。

该组织全年向约 10,000 个家庭水箱供水。

油箱的最大容量为 300 加仑，所需的最小限制为 100 加仑——也就是说，油箱在低于 100 之前应加满 300。

例如，如果油箱在第 2 周为 115 加仑，预计在第 3 周使用 20 加仑，则需要在第 3 周重新填充。

费用包括：

每次送货费 10 美元
卡车的每周成本。一辆卡车每周的成本是 1,000 美元。因此，如果一周内交付 200 次，则成本为 3,000 美元 (200 * 10 + 1000 * 1)。如果交付 201 次，则成本将大幅上升至 4,010 美元 (201 * 10 + 1000 * 2)。

用水量因家庭和周而异。夏季用水高峰。如果我们盲目地遵循规则在达到 100 加仑的最低限度之前重新加注，那么如果交货分散到夏天的“肩膀”，那么卡车的高峰数量可能会超过所需的数量。

我为每个家庭创建了每周用水量的估算值。此外，我对类似家庭进行了分组以减少优化问题的规模（约 10k 家庭减少到 8 组）。

重申目标：此优化器的输出应该是：对于每个家庭组，在一年中的 52 周中的每一周是否交付。

简化数据（即 8 组 12 周）：

df.usage <-  structure(list(reduction.group = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 
                                                1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 
                                                3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 
                                                5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 
                                                7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 
                                                8, 8, 8), week = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 
                                                                   2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 
                                                                   10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 
                                                                   5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 
                                                                   12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 
                                                                   7, 8, 9, 10, 11, 12), water_usage = c(46, 50, 42, 47, 43, 39, 
                                                                                                         38, 32, 42, 36, 42, 30, 46, 50, 42, 47, 43, 39, 38, 32, 42, 36, 
                                                                                                         42, 30, 46, 50, 43, 47, 43, 39, 38, 32, 42, 36, 42, 30, 46, 50, 
                                                                                                         43, 47, 43, 39, 38, 32, 42, 36, 42, 30, 29, 32, 27, 30, 27, 25, 
                                                                                                         24, 20, 26, 23, 27, 19, 29, 32, 27, 30, 27, 25, 24, 20, 26, 23, 
                                                                                                         27, 19, 29, 32, 27, 30, 28, 25, 25, 21, 27, 23, 27, 19, 29, 32, 
                                                                                                         27, 30, 28, 25, 25, 21, 27, 23, 27, 20), tank.level.start = c(115, 
                                                                                                                                                                        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 165, NA, NA, NA, 
                                                                                                                                                                        NA, NA, NA, NA, NA, NA, NA, NA, 200, NA, NA, NA, NA, NA, NA, 
                                                                                                                                                                        NA, NA, NA, NA, NA, 215, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                                                                                                                                        NA, NA, 225, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 230, 
                                                                                                                                                                        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 235, NA, NA, NA, 
                                                                                                                                                                        NA, NA, NA, NA, NA, NA, NA, NA, 240, NA, NA, NA, NA, NA, NA, 
                                                                                                                                                                        NA, NA, NA, NA, NA)), row.names = c(NA, 96L), class = "data.frame")

油箱等级补充规则

以下是一组嵌套循环，用于使用“重新填充”逻辑确定罐液位随时间的变化：

library(dplyr)

reduction.groups <- unique(df.usage$reduction.group)
df.after.refill.logic <- list()

for (i in reduction.groups) {

  temp <- df.usage %>% filter(reduction.group == i)
  temp$refilled <- 0
  temp$level <- temp$tank.level.start

  n <- nrow(temp)

  if (n > 1) for (j in 2:n) {
    temp$level[j] <- ( temp$level[j-1] - temp$water_usage[j] )
    if(temp$level[j] < 100) {
      temp$level[j] <- 300
      temp$refilled[j] <- 1
    }
  }
  df.after.refill.logic <- bind_rows(df.after.refill.logic, temp)
}

决策变量

一年中的每个星期是否向每个组交付或不交付（二进制）

约束

没有部分卡车：卡车的数量必须是整数

卡车容量：卡车交付/周

坦克不能低于 100 加仑：level >= 100

交付必须是二进制的

常量

1600 # truck_weekly_costs
10 # cost_per_delivery
200 # weekly_delivery_capacity_per_truck

成本函数示例

weekly_cost_function <- function(i){
  cost <- (ceiling(sum(i)/200)) * 1600 + (sum(i) * 10)
  cost
}

**example cost for one week with i = 199 deliveries:**
weekly_cost_function(i = 199)
[1] 3590

尝试使用 OMPR 对问题进行建模

以下是使用 OMPR 包创建的模型的开头（尽管使用其他包也可以）：

我对如何使用上述数据进行设置感到困惑。三个明显的问题：

如何在 OMPR 代码中包含 示例成本函数 中表达的上限逻辑？
下面的模型没有将数据合并到上面的数据框中 (df.usage)。目标是优化器根据四个变量（reduction.group、week、water_usage、tank_level_start）以及常量生成“refill”和“level”变量的值。
我在上面的“确定油箱液位”循环中编写的加注逻辑未包含在内。是否应该将其添加为约束？如果有，怎么做？

num_groups <- length(unique(df.usage$reduction.group))
num_weeks <- length(unique(df.usage$week))

MIPModel() %>%
  add_variable(x[i,w],                         # create decision variable: deliver or not by...
               i = 1:num_groups,               # group,
               w = 1:num_weeks,                # in week.
               type = "integer",               # Integers only
               lb = 0, ub = 1) %>%             # between 0 and 1, inclusive 
  set_objective(sum_expr( x[i,w]/200 * 1600 + x[i,w] * 10,
                          i = 1:num_groups, 
                          w = 1:num_weeks),
                sense = "min") %>%
  # add constraint to achieve ceiling(x[i,w]/200), or should this be in the set_objective call?
  add_constraint(???) %>%
  solve_model(with_ROI("glpk"))

期望的输出

这是head() 输出的示例：


 reduction.group   week   water.usage  refill   level
               1      1            46       0     115
               1      2            50       1     300
               1      3            42       0     258
               1      4            47       0     211
               1      5            43       0     168
               1      6            39       0     129

重要的是，refill 值将最小化成本函数并使level 保持在 100 以上。

【问题讨论】：

标签： r optimization linear-programming nonlinear-optimization lpsolve

【解决方案1】：

ceiling 函数是一个困难的非线性函数（不可微分，不连续），应该不惜一切代价避免。然而，它可以很容易地用一般整数变量建模。对于非负变量x >= 0我们可以制定

y = ceiling(x)

作为

x <= y <= x+1
y integer

这是完全线性的，在 OMPR（或任何其他 LP/MIP 工具）中实现起来很简单。

详细说明。此公式将允许模型在 x 假定整数值的特殊情况下选择 y=x 或 y=x+1。如果你想对这个案例挑剔，你可以这样做：

x+0.0001 <= y <= x+1
y integer

我不会担心这个。

【讨论】：

我看到了很多问题，其中大部分是合乎逻辑的。首先写下数学模型。只有在您对数学模型的正确性有信心后才开始编码。

【解决方案2】：

使用上限函数，这对于爬山优化器来说似乎是一个难题。我认为遗传算法更合适。每周每所房子的交付或不交付矩阵构成了一个不错的基因组。

library(dplyr)

# Original given sample input data.
df.usage <-  structure(list(reduction.group = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 
                                                1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 
                                                3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 
                                                5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 
                                                7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 
                                                8, 8, 8), week = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 
                                                                   2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 
                                                                   10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 
                                                                   5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 
                                                                   12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 
                                                                   7, 8, 9, 10, 11, 12), water_usage = c(46, 50, 42, 47, 43, 39, 
                                                                                                         38, 32, 42, 36, 42, 30, 46, 50, 42, 47, 43, 39, 38, 32, 42, 36, 
                                                                                                         42, 30, 46, 50, 43, 47, 43, 39, 38, 32, 42, 36, 42, 30, 46, 50, 
                                                                                                         43, 47, 43, 39, 38, 32, 42, 36, 42, 30, 29, 32, 27, 30, 27, 25, 
                                                                                                         24, 20, 26, 23, 27, 19, 29, 32, 27, 30, 27, 25, 24, 20, 26, 23, 
                                                                                                         27, 19, 29, 32, 27, 30, 28, 25, 25, 21, 27, 23, 27, 19, 29, 32, 
                                                                                                         27, 30, 28, 25, 25, 21, 27, 23, 27, 20), tank.level.start = c(115, 
                                                                                                                                                                       NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 165, NA, NA, NA, 
                                                                                                                                                                       NA, NA, NA, NA, NA, NA, NA, NA, 200, NA, NA, NA, NA, NA, NA, 
                                                                                                                                                                       NA, NA, NA, NA, NA, 215, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                                                                                                                                       NA, NA, 225, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 230, 
                                                                                                                                                                       NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 235, NA, NA, NA, 
                                                                                                                                                                       NA, NA, NA, NA, NA, NA, NA, NA, 240, NA, NA, NA, NA, NA, NA, 
                                                                                                                                                                       NA, NA, NA, NA, NA)), row.names = c(NA, 96L), class = "data.frame")

# Orginal given delivery cost function.
weekly_cost_function <- function(i){
  cost <- (ceiling(sum(i)/200)) * 1600 + (sum(i) * 10)
  cost
}

# Calculate the list of houses (reduction.groups) and number of delivery weeks (weeks).
reduction.groups <- unique(df.usage$reduction.group)
temp             <- df.usage %>% filter(reduction.group == 1)
weeks            <- nrow(temp)

# The genome consists of a matrix representing deliver-or-not to each house each week.
create_random_delivery_schedule <- function(number_of_houses, number_of_weeks, prob = NULL) {
  matrix(sample(c(0, 1), number_of_houses * number_of_weeks, replace = TRUE, prob = prob), number_of_houses)
}

# Generate a population of random genes.
population_size <- 100
schedules <- replicate(population_size, create_random_delivery_schedule(length(reduction.groups), weeks), simplify = FALSE)

# Calculate fitness of an individual.
fitness <- function(schedule) {

  # Fitness is related to delivery cost.
  delivery_cost <- sum(apply(schedule, 2, weekly_cost_function))

  # If the schedule allows a tank level to drop below 100, apply a fitness penalty.
  # Don't make the fitness penalty too large.
  # If the fitness penalty is large enough to be catastrophic (essentially zero children)
  # then solutions that are close to optimal will also be likely to generate children
  # who fall off the catastropy cliff so there will be a selective pressure away from
  # close to optimal solutions.
  # However, if your optimizer generates a lot of infeasible solutions raise the penalty.
  for (i in reduction.groups) {

    temp <- df.usage %>% filter(reduction.group == i)
    temp$level <- temp$tank.level.start

    if (weeks > 1) for (j in 2:weeks) {
      if (1 == schedule[i,j]) {
        temp$level[j] <- 300
      } else {
        temp$level[j] <- ( temp$level[j-1] - temp$water_usage[j] )

        if (100 > temp$level[j]) {
          # Fitness penalty.
          delivery_cost <- delivery_cost + 10 * (100 - temp$level[j])
        }
      }
    }
  }

  # Return one over delivery cost so that lower cost is higher fitness.
  1 / delivery_cost
}

# Generate a new schedule by combining two parents chosen randomly weighted by fitness.
make_baby <- function(population_fitness) {

  # Choose some parents.
  parents <- sample(length(schedules), 2, prob = population_fitness)

  # Get DNA from mommy.
  baby <- schedules[[parents[1]]]

  # Figure out what part of the DNA to get from daddy.
  house_range <- sort(sample(length(reduction.groups), 2))
  week_range  <- sort(sample(weeks, 2))

  # Get DNA from daddy.
  baby[house_range[1]:house_range[2],week_range[1]:week_range[2]] <- schedules[[parents[2]]][house_range[1]:house_range[2],week_range[1]:week_range[2]]

  # Mutate, 1% chance of flipping each bit.
  changes <- create_random_delivery_schedule(length(reduction.groups), weeks, c(0.99, 0.01))
  baby <- apply(xor(baby, changes), c(1, 2), as.integer)
}

lowest_cost <<- Inf

# Loop creating and evaluating generations.
for (ii in 1:100) {
  population_fitness <- lapply(schedules, fitness)
  lowest_cost_this_generation <- 1 / max(unlist(population_fitness))
  print(sprintf("lowest cost = %f", lowest_cost_this_generation))

  if (lowest_cost_this_generation < lowest_cost) {
    lowest_cost <<- lowest_cost_this_generation
    best_baby   <<- schedules[[which.max(unlist(population_fitness))]]
  }

  schedules <<- replicate(population_size, make_baby(population_fitness), simplify = FALSE)
}

【讨论】：

我只是简要回顾了这个答案，但它看起来很棒。我现在就给悬赏，快到期了，以后有问题再来！谢谢鲍勃。