下面详细介绍了我无法解决的优化问题的简化版本。
目标是最大限度地降低通过卡车运输水的组织的成本函数，并使用该公式生成最小化成本的卡车交付计划。
该组织全年为约10,000个家用水箱提供水。
这些坦克的最大容量为300加仑，最小期望限制为100加仑 - 也就是说，坦克应该在低于100之前加满300。
例如，如果第2周的油箱为115加仑，并且估计在第3周使用20加仑，则需要在第3周重新加油。
费用包括：
(200 * 10 + 1000 * 1)。如果201交付，成本显着跳跃到4,010美元(201 * 10 + 1000 * 2)。不同家庭和数周的用水量不同。高峰用水量在夏季。如果我们在达到100加仑最低限度之前盲目地遵循规则重新填充，那么如果交付分散到夏季的“肩膀”，卡车的峰值数量可能会高于所需的数量。
我创建了每个家庭每周每周用水量的估算值。此外，我像家庭一样分组，以减少优化问题的规模（约10万户家庭减少到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)
}
决策变量
一年中的每个星期每个小组交付与否（二进制）
约束
没有部分卡车：卡车数量必须是整数
卡车容量：卡车交付/周<= 200
坦克不能低于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包创建的模型的开头（虽然使用另一个包可以）：
我对如何使用上面的数据进行设置感到困惑。三个明显的问题：
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的任何值。
ceiling函数是一个困难的非线性函数（不可微分，不连续），应该不惜一切代价避免。但是，它可以使用常规整数变量很容易地建模。对于非负变量x >= 0，我们可以制定
y = ceiling(x)
如
x <= y <= x+1
y integer
这是完全线性的，在OMPR（或任何其他LP / MIP工具）中实现是微不足道的。
细节笔记。这个公式将允许模型在y=x假定整数值的特殊情况下选择y=x+1或x。如果你想对这种情况挑剔，你可以这样做：
x+0.0001 <= y <= x+1
y integer
我不担心这个。
有了天花板功能，对于爬山优化器来说，这似乎是一个难题。我认为遗传算法更适合。每周每个房子的交付或不交付矩阵都是一个很好的基因组。
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)
}