Paradoxical as it may seem, but still quite often in the enterprise there are tasks that are different from building another personal account, another monitoring or another document flow. If you think a little, and do not grab at once to code or look for specialized software, you can write a compact, very elegant and fast solution using the Monte Carlo method.

Enterprise tasks are compact enough to iterate over and do not require 100 decimal places. We are not launching rockets or reactors and building a scientific theory of everything.

Let's consider further on the example of one of the non-standard tasks.

It is a continuation of a series of previous publications .

Formulation of the problem

The challenge has its roots in IoT for agriculture. Namely, the placement of sensors in a potato field with a circular irrigation pattern. Let's ask Google for some pictures: "The answer to the mystery of round fields: everything is more interesting than you think!" ... It is necessary to provide the necessary characteristics of the mesh network coverage, taking into account the permissible distances between neighbors. At the same time, it is necessary to optimize the budget and issue gps coordinates for the placement of sensors and form the shortest bypass scheme.

Decision

We connect packages

Packages

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

Decomposition

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

:

(- );
.

s = a t^{2} / 2 = (F / m) t^{2} / 2

$s = at^2/2 = (F/m)t^2/2$

Δ s = F Δ t

$\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

Bypass route

, . GPS, GPS . .

, . « », . tidyverse

$10^3$-$10^4$ . . , , C++. , .

data.table

. Base R .
, . .
data.frame

. . tidyverse

.
, , .
Monte Carlo is a very good approach. A quick initial application can give a useful result, as well as look at the solution to the problem and, possibly, find some simplifications and analytics.
Feel free to use analogy methods. They can make it possible to construct a simplified model of the original problem, which is computationally much simpler than the original one and can be easily transferred to Monte Carlo.

Previous publication - "Children, Russian and R" .

R, Monte Carlo and Enterprise Problems, Part 2

Formulation of the problem

Decision

We connect packages

Decomposition

-

More articles: