i_shutov 5 апр 2021 в 12:16

R, Монте-Карло и enterprise задачи, часть 2

7 мин

2.1K

Программирование*Data Mining*Алгоритмы*R*

Как ни парадоксально, но пока еще нередко в enterprise встречаются задачи, отличные от построения еще одного личного кабинета, еще одного мониторинга или еще одного документооборота. Если немного подумать, а не хвататься сразу кодировать или искать специализированное ПО, то можно написать компактное, весьма элегантное и быстрое решение, используя метод Монте-Карло .

Задачи в Enterprise достаточны компактны для перебора и не требует точности 100 знаков после запятой. Не ракеты или реакторы запускаем и не научную теорию всего строим.

Рассмотрим далее на примере одной из нестандартных задач.

Является продолжением серии предыдущих публикаций.

Постановка задачи

Задача имеет корни в IoT для сельского хозяйства. А именно, расстановка датчиков на картофельном поле с круговой схемой полива. Попросим у гугла немного картинок: «Разгадка тайны круглых полей: всё интересней, чем ты думаешь!». Нужно обеспечить нужные характеристики покрытия mesh сети, учитывая допустимые расстояния между соседями. При этом надо оптимизировать бюджет и выдать gps координаты на расстановку датчиков и сформировать кратчайшую схему обхода.

Решение

Подключаем пакеты

Пакеты

library(tidyverse)
library(magrittr)
library(scales)
library(data.table)
library(tictoc)
library(stringi)
library(arrangements)
library(ggforce)

Декомпозиця

Сначала попробуем сформировать этапы.

Берем N датчиков.
Ищем оптимальную расстановку.
Если характеристики не достигнуты, увеличиваем N на 1. Повторяем процедуру.
Для найденной расстановки ищем оптимальный маршрут обхода.

План вроде простой. Но как будем решать? Напрашивается метод Монте-Карло.

Функции-хелперы

Оформим функцию для визуального отображения расположения зондов.

Визуализация поля

drawDisk <- function(df) {
  # отрисуем расположение точек и действующие силы
  # если силы не заданы, создадим их по умолчанию равными 0
  for(col in c("force_x", "force_y")){
    if (!(col %in% names(df))) df[, col] <- 0
  }

  ggplot(data = df, aes(x = x, y = y)) +
    ggforce::geom_circle(aes(x0 = 0, y0 = 0, r = 1, colour = "red"), 
                         inherit.aes = FALSE) +
    geom_point(size = 2, fill = "black", shape = 21) +
    geom_text(aes(label = idx), hjust = 0.5, vjust = -1) +
    # рисуем векторное поле
    geom_segment(aes(xend = x + force_x / 10, yend = y + force_y / 10), 
                 colour = "red", 
                 arrow = arrow(length = unit(0.2,"cm")), size = 0.6) + 
    xlim(-1.5, 1.5) +
    ylim(-1.5, 1.5) +
    coord_fixed() +
    labs(x = "Ось X", y = "Ось Y") +
    theme_bw()
}

Генерация первичной расстановки

Конечно же, и для этой задачи можно использовать Монте-Карло. Но только далеко не факт, что на малых выборках подберется оптимальная конфигурация. На больших выборках при неограниченном количестве состояний (а пространство круга содержит бесконечное количество точек) и неудачном коде не хватит никаких вычислительных ресурсов.

Но есть выход — можем вспомнить физику.

Предположим, что мы накидаем в круглую плоскую чашку N одинаково заряженных шариков. Они либо найдут устойчивую конфигурацию с минимум энергии, либо возникнут колебания. Непосредственно в этой задаче с большой вероятностью будет найдено стабильное расположение (из практики). Попробуем провести такую же симуляцию для расстановки зондов. Наделим каждый зонд единичным зарядом с единичной массой, радиус нашего картофельного поля выберем единичным. Для того, чтобы зонды не разлетелись в бесконечность, вокруг нашего поля поставим заряженных неподвижных «охранников». И проведем симуляцию для случайной первичной расстановки.

Генерация первичной расстановки

# Генерим точки-зонды внутри окружности единичного радиуса.
# Считаем, что все частицы единичного заряда, поэтому его опускаем
charges_dt <- tibble(idx = 1:13) %>%
  mutate(angle = runif(n(), min = 0, max = 2*pi), 
         r = runif(n(), min = 0, max = 1),
         x = r * cos(angle), y = r * sin(angle)) %>%
  select(idx, x, y) %>%
  setDT(key = "idx")

# для сходимости задачи генерируем также зафиксированные точки на внешней окружности
keepers_dt <- max(charges_dt$idx) %>% 
  {tibble(idx = (. + 1):(. + 40))} %>%
  mutate(angle = (idx - 1) * (2 * pi / n()),
         x = 1.3 * cos(angle), y = 1.3 * sin(angle)) %>%
  select(idx, x, y) %>%
  setDT(key = "idx")

Посмотрим, как расположились зонды при первичной установке.

Визуализируем

full_dt <- rbindlist(list(charges_dt, keepers_dt))
drawDisk(full_dt)

Поиск оптимального расположения

Здесь задействуем физику, расчитаем малое перемещение частиц за счет действия на них cилы электростатического взаимодействия (закон Кулона).
Для упрощения задачи будем считать:

каждую новую итерацию как движение из статичного состояния (как-бы случайная расстановка зондов);
все зонды обладают единичным зарядом и единичной массой.

$display$

Для малых изменений получаем

$\Delta s = F \Delta t$

Будем итеративно двигать, пока равнодействующие сил на зонды не станут меньше определенного порога.

Поиск оптимального расположения

max_force <- 10

tic("Balancing charges")
# Определяем точность моделирования!
# Будем двигать заряды пока они не застабилизируются 
# Максимальная равнодействуюущая будет близка к 0 
# и точки не стали разлетаться в бесконечность
while (max_force > 0.05 & nrow(charges_dt[x^2 + y^2 > 1.2]) == 0) {
  # общий пул координат частиц на текущую итерацию
  full_dt <- rbindlist(list(charges_dt, keepers_dt), fill = TRUE)

  ff <- function(x0, y0){
    # сила взаимодействия -- 1/r2, заряды единичные;
    # проекция силы на оси, sqrt(r2) -- гипотенуза
    # суперпозиция векторов даст результирующее воздействие на частицу

    dx = full_dt$x - x0
    dy = full_dt$y - y0
    r2 = dx^2 + dy^2
    # na.rm исключает NaN в т.ч.
    list(sum(-dx * r2^(-1.5), na.rm = TRUE),
         sum(-dy * r2^(-1.5), na.rm = TRUE))
  }  

  # проведем расчет сил, итерируем по каждой "неприбитой гвоздями" точке
  charges_dt[, c("force_x", "force_y") := ff(x0 = x, y0 = y), by = idx]
  # определим максимальную силу, действующую на частицу
  max_force <- charges_dt[, max(sqrt(force_x^2 + force_y^2), na.rm = TRUE)]
  force_scale = if_else(max_force > 1, 1 / max_force / 1e2, 1/ max_force / 5e2)

  # проводим передвижение точек
  charges_dt %>%
    .[, `:=`(x = x + force_x * force_scale, 
             y = y + force_y * force_scale)]
}
toc()

full_dt <- rbindlist(list(charges_dt, keepers_dt), fill = TRUE)

Оптимизация маршрута обхода

Для выбора оптимального маршрута опять же используем Монте-Карло. Ряд соображений:

для генерации подмножества обходов используем комбинаторные функции;
для каждого зонда добавим фиктивные точки входа и выхода, которые будут находиться на окружности вдоль радиуса, проходящего через зонд (кратчайший маршрут);
для сокращения временных затрат сделаем однократный расчет матриц расстояний.

Оптимизация маршрута обхода

optimizePath <- function(dt) {
  # попробуем оптимизировать маршрут обхода по предоставленным точкам
  # 1. Выбираем в качестве начальной точки датчик, максимально близко расположенный к краю поля
  dt[, r0 := sqrt(x^2+y^2)] %>%
    setorder(-r0)
  n1 <- dt[1, idx]

  # теперь проводим симуляцию различных вариантов расстановки сенсоров
  # получаем последовательность номеров и убираем n1, его будем принудительно ставить первым
  points_in <- dt[idx != n1, idx]

  # для каждой точки добавим еще ближайшую точку выхода 
  # (перпендикуляр к окружности, которая единичного радиуса)
  augment_tbl <- dt %>%
    mutate_at("idx", `*`, -1) %>%
    mutate(r0 = sqrt(x^2 + y^2)) %>%
    mutate_at(vars(x, y), ~(.x/r0)) %>%
    bind_rows(dt) %>%
    select(idx, x, y)

  # однократно посчитаем матрицу расстояний между зондами
  ll_tbl <- unique(augment_tbl$idx) %>%
    tidyr::expand_grid(l = ., r = .) %>%
    filter(l != r, (l > 0 | r > 0)) %>%
    # построим уникальный идентификатор ребра
    rowwise() %>%
    # mutate(m = list(sort(c(l, r))))
    mutate(edge_id = stri_c(sort(c(l, r)), collapse = "=")) %>%
    ungroup() %>%
    distinct(edge_id, .keep_all = TRUE) %>%
    # подтягиваем координаты левой точки
    left_join(select(augment_tbl, idx, l_x = x, l_y = y), by = c("l" = "idx")) %>%
    # подтягиваем координаты левой точки
    left_join(select(augment_tbl, idx, r_x = x, r_y = y), by = c("r" = "idx")) %>%
    mutate(s = sqrt((l_x - r_x)^2 + (l_y - r_y)^2)) %>%
    arrange(l, r)

  points_seq <- arrangements::permutations(v = points_in, replace = FALSE, 
                                       layout = "column", nsample = 5000)
  # добавляем точку входа в качестве первой и соотв. точку выхода в качестве последней
  routes_lst <- points_seq %>% 
    rbind(-n1, n1, ., -tail(., 1)) %>%
    as.data.frame() %>% as.list()

  # превращаем все пути обхода в последовательности ребер
  routes_dt <- data.table(route_id = seq_along(routes_lst), route = routes_lst) %>%
    .[, .(from = unlist(route)), by = route_id] %>%
    .[, to := shift(from, type = "lead"), by = route_id] %>%
    # выкидываем все терминальные точки
    na.omit() %>%
    # строим нормализованный идентификатор ребра
    .[, edge_id := stri_c(sort(unlist(.BY)), collapse = "="), by = .(from, to)] %>%
    .[, .(route_id, edge_id)] %>%
    # подтянем информацию о длине ребра из справочника
    .[as.data.table(ll_tbl), s := i.s, on = "edge_id"]

  # считаем длину маршрутов, оставляем кратчайший
  best_routes <- routes_dt[, .(len = sum(s)), by = route_id] %>%
    setorder(len) %>%
    head(10) %T>%
    print()

  # сформируем ТП-10 лучших маршрутов
  best_routes %>%
    select(route_id) %>%
    mutate(idx = routes_lst[route_id]) %>%
    tidyr::unnest(idx) %>%
    left_join(augment_tbl) %>%
    tidyr::nest(data = -route_id) %>%
    left_join(best_routes)
}

Получаем табличку подобного рода

    route_id      len
 1:     2070 8.332854
 2:     2167 8.377680
 3:     4067 8.384417
 4:     3614 8.418678
 5:     5000 8.471521
 6:     4495 8.542041
 7:     2233 8.598278
 8:     4430 8.609391
 9:     2915 8.616048
10:     3380 8.695547

И посмотрим результат размещения

Визуализируем

tic("Route optimization")
best_tbl <- optimizePath(charges_dt)
toc()

best_route_tbl <- best_tbl$data[[1]]
full_dt <- rbindlist(list(best_route_tbl, keepers_dt), fill = TRUE)
gp <- drawDisk(full_dt) +
  # добавим маршрут обхода
  geom_path(arrow = arrow(type = "closed"), data = best_route_tbl)

gp

Маршрут обхода

Формирование задания

У нас есть пачка оптимальных расстановок, минимизированных по обходу. Переводим условные координаты в GPS, формируем маршрут обхода и выдаем GPS трекер агроному. Можно ставить зонды.

Полезные трюки

Как всегда, затронем вопросы производительности. Если писать «не думая», то можно отправить машину на счет на часы или дни. tidyverse подход оказывается медленнее в $10^3$-$10^4$ раз. В приведенном выше коде расчеты по блокам занимают доли секунды. Этого достаточно для обычной задачи, но если нужно быстрее, то можно делать вставки на C++. В целом, скоростные характеристики достигаются результат рядом мер и методик.

Для небольших циклических расчетов накладные расходы на разыменовывание переменных даже в data.table могут оказаться значительными. Base R в блоке поиска оптимального расположения дает выйгрыш на порядок.
Если задачу можно распараллелить, то надо применять функциональные подходы. Проще будет сделать последующее распараллеливание.
Для однородных величин работа с матрицами оказывается на несколько порядков быстрее работы с data.frame. Связано это со схемой выделения памяти и адресации к элементам. Про матрицы незаслуженно забывают при погружении в tidyverse.
Все, что можно посчитать однократно и оформить в виде справочной таблицы, должно быть посчитано заранее.
Монте-Карло — очень хороший подход. Быстрое первичное применение может дать полезный результат, а также взглянуть на решение задачи и, возможно, найти какие-то упрощения и аналитики.
Не стесняемся использовать методы аналогии. Они могут позволить построить упрощенную модель исходной задачи, которая вычислительно существенно проще исходной и легко перекладывается на Монте-Карло.

Предыдущая публикация — «Дети, русский язык и R».

Теги:

Хабы: