Оптимизация в R: функция затрат с бинарными переменными планирования? - PullRequest
Новая
       35

Оптимизация в R: функция затрат с бинарными переменными планирования?

4 голосов
/ 30 апреля 2019

Ниже приведена упрощенная версия проблемы оптимизации, решение которой у меня возникло.

Цель - минимизировать функцию затрат для организации, которая доставляет воду через грузовики и используетэто уравнение для составления графика доставки грузовиков, который сводит к минимуму расходы.

Организация поставляет воду в ~ 10 000 бытовых резервуаров в течение года.

Максимальная вместимость танков составляет 300 галлонов, а минимальный желаемый предел - 100 галлонов, т. Е. Танки должны быть заполнены до 300, прежде чем они опустятся ниже 100.

Например,если на 2-й неделе объем бака составляет 115 галлонов, а на неделе 3 - 20 галлонов, его необходимо пополнить на третьей неделе.

Расходы включают в себя:

  1. Плата за доставку в размере 10

  2. Недельная стоимость грузовых автомобилей.Еженедельная стоимость грузовика составляет 1000 долларов.Таким образом, если за одну неделю осуществляется 200 доставок, стоимость составляет 3000 (200 * 10 + 1000 * 1). Если будет произведено 201 поставка, стоимость значительно возрастет до 4 010 (201 * 10 + 1000 * 2).

Расход воды варьируется в зависимости от домохозяйстваи через недели.Пиковое использование воды приходится на лето.Если мы слепо следуем правилу пополнения запасов непосредственно перед достижением минимального предела в 100 галлонов, то, вероятно, максимальное количество грузовых автомобилей будет выше, чем необходимо, если поставки будут распределены по «плечам» лета.

Я создал оценки еженедельного использования воды для каждой недели для каждого домохозяйства.Кроме того, я сгруппировал как домашние хозяйства, чтобы уменьшить размер проблемы оптимизации (~ 10 000 домашних хозяйств до 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 </p>

Танки не могут быть ниже 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 p.ackage (хотя использование другого пакета будет в порядке):

Я не уверен, как настроить это, используя данные выше.Три очевидные проблемы:

  1. Как я могу включить предельную логику, выраженную в Пример функции стоимости в коде OMPR?
  2. Модель, приведенная ниже, не включаетданные в приведенном выше кадре (df.usage).Цель состоит в том, чтобы оптимизатор генерировал значения для переменных «заправка» и «уровень» на основе четырех переменных (extension.group, week, water_usage, tank_level_start) вместе с константами.
  3. Логика пополненияЯ написал в цикле «определение уровня бака» выше, не включены.Должно ли это быть добавлено как ограничение?Если да, то как? 
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.

Ответы [ 2 ]

4 голосов
/ 30 апреля 2019

Функция 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

Я бы не беспокоился об этом.

2 голосов
/ 08 мая 2019

С функцией потолка это кажется сложной задачей для оптимизатора подъема на гору.Я думаю, что генетические алгоритмы лучше подходят.Матрица доставки или нет для каждого дома каждую неделю дает хороший геном. 

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)
}
...