Симуляция Монте-Карло в R

Я пытаюсь выполнить следующее упражнение, но на самом деле понятия не имею, как начать это делать. Я нашел код в своей книге, который выглядит так, но это совершенно другое упражнение, и я не знаю, как связать их друг с другом. Как я могу начать имитировать прибытие и как узнать, когда они закончатся? Я знаю, как их хранить и вычислять a, b, c, d в соответствии с этим. Но я не знаю, как мне на самом деле нужно имитировать симуляцию Монте-Карло. Может ли кто-нибудь помочь мне начать? Я знаю, что это не место, где ваши вопросы получают ответы для вас, а только решаются вместо этого. Но проблема в том, что я не знаю с чего начать.

Служба поддержки ИТ-службы представляет собой систему очередей с пятью помощниками, принимающими звонки от клиентов. Звонки происходят в соответствии с пуассоновским процессом со средней скоростью одного звонка каждые 45 секунд. Время обслуживания для 1-го, 2-го, 3-го, 4-го и 5-го помощников - все экспоненциальные случайные величины с параметрами λ1 = 0,1, λ2 = 0,2, λ3 = 0,3, λ4 = 0,4 и λ5 = 0,5 мин − 1 соответственно ( у помощника справочной службы λk = k / 10 мин − 1). Помимо клиентов, которым оказывается помощь, до десяти других клиентов могут быть приостановлены. В тех случаях, когда эта емкость достигнута, новые абоненты получают сигнал занятости. Используйте методы Монте-Карло для оценки следующих характеристик производительности,

(а) доля клиентов, которые получают сигнал занятости;

(б) ожидаемое время ответа;

(в) среднее время ожидания;

(d) часть клиентов, обслуживаемых каждым помощником службы поддержки;

РЕДАКТИРОВАТЬ: что я до сих пор (не очень):

pa = 1/45sec-1

jobs = rep(1,5); onHold = rep(1,10);

jobsIndex = 0;

onHoldIndex = 0;

u = runif(1)
for (i in 1:1000) {

    if(u  <= pa){ # new arrival

        if(jobsIndex < 5) # assistant is free, #give job to assistant

            jobsIndex++;

        else #add to onHold array

            onHoldIndex++;
    }
}

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

Лучший способ построить такое моделирование - это дизайн сверху вниз.

На самом высоком уровне код должен выглядеть примерно так

initialize(...)
while (process(get.next.event())) {}

(Этот и все последующие примеры являются исполняемым R кодом, а не просто псевдокодом.) Цикл представляет собой симуляцию, управляемую событиями : get.next.event()находит любое интересующее «событие» и передает его описание process, которое что-то с ним делает (включая запись любого информация об этом). Он возвращается TRUEдо тех пор, пока все идет хорошо; после определения ошибки или окончания симуляции возвращается FALSE, завершая цикл.

Если мы представим физическую реализацию этой очереди, например, люди, ожидающие получения свидетельства о браке в Нью-Йорке или водительского удостоверения или билета на поезд почти в любом месте, мы представим два вида агентов: клиенты и «помощники» (или серверы). , Клиенты объявляют себя, показывая; Помощники сообщают о своей доступности, включив свет или знак или открыв окно. Это два вида событий для обработки.

Идеальная среда для такого моделирования - это действительно объектно-ориентированная среда, в которой объекты являются изменяемыми : они могут изменять состояние, чтобы реагировать независимо на окружающие их вещи. Rабсолютно ужасен для этого (даже Fortran был бы лучше!). Тем не менее, мы можем использовать его, если позаботимся. Хитрость заключается в том, чтобы поддерживать всю информацию в общем наборе структур данных, которые могут быть доступны (и изменены) многими отдельными, взаимодействующими процедурами. Я приму соглашение об использовании имен переменных IN ALL CAPS для таких данных.

Следующий уровень нисходящего дизайна - это код process. Он отвечает на один дескриптор события e:

process <- function(e) {
  if (is.null(e)) return(FALSE)
  if (e$type == "Customer") {
    i <- find.assistant(e$time)
    if (is.null(i)) put.on.hold(e$x, e$time) else serve(i, e$x, e$time)
  } else {
    release.hold(e$time)
  }
  return(TRUE)
}

Он должен отвечать на нулевое событие, когда get.next.eventнет событий для отчета. В противном случае processреализуются «бизнес-правила» системы. Это практически пишет себя из описания в вопросе. Как это работает, нужно немного комментировать, за исключением указания на то, что в конечном итоге нам нужно будет кодировать подпрограммы put.on.holdи release.hold(реализация очереди, удерживающей клиента) и serve(реализация взаимодействия клиент-помощник).

Что такое "событие"? Он должен содержать информацию о том, кто действует, какие действия они предпринимают и когда они происходят. Поэтому мой код использует список, содержащий эти три вида информации. Однако get.next.eventнужно только проверять время. Он отвечает только за поддержание очереди событий, в которых

1. Любое событие может быть помещено в очередь, когда оно получено и

2. Самое раннее событие в очереди может быть легко извлечено и передано вызывающей стороне.

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

get.next.event <- function() {
  if (length(EVENTS$time) <= 0) new.customer()               # Wait for a customer$
  if (length(EVENTS$time) <= 0) return(NULL)                 # Nothing's going on!$
  if (min(EVENTS$time) > next.customer.time()) new.customer()# See text
  i <- which.min(EVENTS$time)
  e <- EVENTS[i, ]; EVENTS <<- EVENTS[-i, ]
  return (e)
}

Есть много способов, которыми это могло быть закодировано. Окончательная версия, показанная здесь, отражает выбор, который я сделал при кодировании, как processреагировать на событие «Помощник» и как new.customerработает: get.next.eventпросто выводит клиента из очереди удержания, затем бездействует и ждет другого события. Иногда будет необходимо искать нового клиента двумя способами: во-первых, чтобы увидеть, ждет ли кто-нибудь за дверью (как это было бы), и во-вторых, вошел ли кто-то, когда мы не смотрели.

Понятно new.customerи next.customer.timeважны рутины , поэтому позаботимся о них дальше.

new.customer <- function() {  
  if (CUSTOMER.COUNT < dim(CUSTOMERS)[2]) {
    CUSTOMER.COUNT <<- CUSTOMER.COUNT + 1
    insert.event(CUSTOMER.COUNT, "Customer", 
                 CUSTOMERS["Arrived", CUSTOMER.COUNT])
  }
  return(CUSTOMER.COUNT)
}
next.customer.time <- function() {
  if (CUSTOMER.COUNT < dim(CUSTOMERS)[2]) {
    x <- CUSTOMERS["Arrived", CUSTOMER.COUNT]
  } else {x <- Inf}
  return(x) # Time when the next customer will arrive
}

CUSTOMERS2D-массив с данными по каждому клиенту в столбцах. Он имеет четыре строки (выступающие в качестве полей), которые описывают клиентов и записывают их опыт во время моделирования : «Прибыл», «Обслуживаем», «Продолжительность» и «Ассистент» (положительный числовой идентификатор ассистента, если таковой имеется, который обслуживал их, а в остальном -1для занятых сигналов). При очень гибком моделировании эти столбцы будут генерироваться динамически, но из-за того, как им Rнравится работать, удобно с самого начала генерировать всех клиентов в единой большой матрице, причем время их прибытия уже генерируется случайным образом. next.customer.timeможете заглянуть в следующий столбец этой матрицы, чтобы увидеть, кто будет следующим. Глобальная переменнаяCUSTOMER.COUNTуказывает на последнего клиента, чтобы прибыть. С помощью этого указателя клиенты управляются очень просто, продвигаясь к нему, чтобы получить нового клиента, и не обращая внимания на него (не продвигаясь), чтобы посмотреть на следующего клиента.

serve реализует бизнес-правила в симуляции.

serve <- function(i, x, time.now) {
  #
  # Serve customer `x` with assistant `i`.
  #
  a <- ASSISTANTS[i, ]
  r <- rexp(1, a$rate)                       # Simulate the duration of service
  r <- round(r, 2)                           # (Make simple numbers)
  ASSISTANTS[i, ]$available <<- time.now + r # Update availability
  #
  # Log this successful service event for later analysis.
  #
  CUSTOMERS["Assistant", x] <<- i
  CUSTOMERS["Served", x] <<- time.now
  CUSTOMERS["Duration", x] <<- r
  #
  # Queue the moment the assistant becomes free, so they can check for
  # any customers on hold.
  #
  insert.event(i, "Assistant", time.now + r)
  if (VERBOSE) cat(time.now, ": Assistant", i, "is now serving customer", 
                   x, "until", time.now + r, "\n")
  return (TRUE)
}

Это просто. ASSISTANTSэто информационный фрейм с двумя полями: capabilities(с указанием их скорости обслуживания) и available, который обозначает следующий раз, когда помощник будет свободен. Клиент обслуживается путем генерации случайной продолжительности обслуживания в соответствии с возможностями помощника, обновления времени, когда помощник становится доступным в следующий раз, и регистрации интервала обслуживания в CUSTOMERSструктуре данных. VERBOSEФлаг удобен для тестирования и отладки: если верно, испускает поток английских фраз , описывающих ключевые моменты обработки.

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

find.assistant <- function(time.now) {
  j <- which(ASSISTANTS$available <= time.now)
  #if (length(j) > 0) {
  #  i <- j[ceiling(runif(1) * length(j))]
  #} else i <- NULL                                    # Random selection
  #if (length(j) > 0) i <- j[1] else i <- NULL         # Pick first assistant
  #if (length(j) > 0) i <- j[length(j)] else i <- NULL # Pick last assistant
  if (length(j) > 0) {
    i <- j[which.min(ASSISTANTS[j, ]$available)]
  } else i <- NULL                                     # Pick most-rested assistant
  return (i)
}

Остальная часть симуляции на самом деле является обычным упражнением для убеждения Rвнедрить стандартные структуры данных, главным образом циклический буфер для очереди на удержании. Поскольку вы не хотите разыгрывать глобалы, я собрал все это в одну процедуру sim. Его аргументы описывают проблему: количество имитируемых клиентов ( n.events), коэффициент прибытия клиентов, возможности помощников и размер очереди удержания (которую можно установить равной нулю, чтобы полностью исключить очередь).

r <- sim(n.events=250, arrival.rate=60/45, capabilities=1:5/10, hold.queue.size=10)

Возвращает список структур данных, поддерживаемых во время моделирования; Наибольший интерес представляет CUSTOMERSмассив. Rпозволяет довольно легко представить необходимую информацию в этом массиве интересным способом. Вот один вывод, показывающий последние клиентов в более длительной симуляции из клиентов.$50$$250$

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

При запуске verbose=TRUE, текстовый вывод симуляции выглядит так:

...
160.71 : Customer 211 put on hold at position 1 
161.88 : Customer 212 put on hold at position 2 
161.91 : Assistant 3 is now serving customer 213 until 163.24 
161.91 : Customer 211 put on hold at position 2 
162.68 : Assistant 4 is now serving customer 212 until 164.79 
162.71 : Assistant 5 is now serving customer 211 until 162.9 
163.51 : Assistant 5 is now serving customer 214 until 164.05 
...

(Числа слева - это время, когда было отправлено каждое сообщение.) Вы можете сопоставить эти описания с частями графика «Клиенты», лежащими между и .$160$$165$

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

(Разве все эти графики не стали бы прекрасной панелью мониторинга в реальном времени для тех, кто управляет этой очередью услуг!)

Интересно сравнить графики и статистику, которые вы получаете при изменении параметров, передаваемых sim. Что происходит, когда клиенты приходят слишком быстро для обработки? Что происходит, когда очередь удержания уменьшается или удаляется? Что меняется, когда помощники выбираются по-разному? Как количество и возможности помощников влияют на качество обслуживания клиентов? Каковы критические моменты, когда некоторые клиенты начинают отворачиваться или начинают задерживаться на долгое время?

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

( Обработка на этом сайте испортит отступ в любых строках, содержащих символ , но читаемый отступ должен быть восстановлен, когда код вставлен в текстовый файл.)$\TeX$$\$$

sim <- function(n.events, verbose=FALSE, ...) {
  #
  # Simulate service for `n.events` customers.
  #
  # Variables global to this simulation (but local to the function):
  #
  VERBOSE <- verbose         # When TRUE, issues informative message
  ASSISTANTS <- list()       # List of assistant data structures
  CUSTOMERS <- numeric(0)    # Array of customers that arrived
  CUSTOMER.COUNT <- 0        # Number of customers processed
  EVENTS <- list()           # Dynamic event queue   
  HOLD <- list()             # Customer on-hold queue
  #............................................................................#
  #
  # Start.
  #
  initialize <- function(arrival.rate, capabilities, hold.queue.size) {
    #
    # Create common data structures.
    #
    ASSISTANTS <<- data.frame(rate=capabilities,     # Service rate
                              available=0            # Next available time
    )
    CUSTOMERS <<- matrix(NA, nrow=4, ncol=n.events, 
                         dimnames=list(c("Arrived",  # Time arrived
                                         "Served",   # Time served
                                         "Duration", # Duration of service
                                         "Assistant" # Assistant id
                         )))
    EVENTS <<- data.frame(x=integer(0),              # Assistant or customer id
                          type=character(0),         # Assistant or customer
                          time=numeric(0)            # Start of event
    )
    HOLD <<- list(first=1,                           # Index of first in queue
                  last=1,                            # Next available slot
                  customers=rep(NA, hold.queue.size+1))
    #
    # Generate all customer arrival times in advance.
    #
    CUSTOMERS["Arrived", ] <<- cumsum(round(rexp(n.events, arrival.rate), 2))
    CUSTOMER.COUNT <<- 0
    if (VERBOSE) cat("Started.\n")
    return(TRUE)
  }
  #............................................................................#
  #
  # Dispatching.
  #
  # Argument `e` represents an event, consisting of an assistant/customer 
  # identifier `x`, an event type `type`, and its time of occurrence `time`.
  #
  # Depending on the event, a customer is either served or an attempt is made
  # to put them on hold.
  #
  # Returns TRUE until no more events occur.
  #
  process <- function(e) {
    if (is.null(e)) return(FALSE)
    if (e$type == "Customer") {
      i <- find.assistant(e$time)
      if (is.null(i)) put.on.hold(e$x, e$time) else serve(i, e$x, e$time)
    } else {
      release.hold(e$time)
    }
    return(TRUE)
  }#$
  #............................................................................#
  #
  # Event queuing.
  #
  get.next.event <- function() {
    if (length(EVENTS$time) <= 0) new.customer()
    if (length(EVENTS$time) <= 0) return(NULL)
    if (min(EVENTS$time) > next.customer.time()) new.customer()
    i <- which.min(EVENTS$time)
    e <- EVENTS[i, ]; EVENTS <<- EVENTS[-i, ]
    return (e)
  }
  insert.event <- function(x, type, time.occurs) {
    EVENTS <<- rbind(EVENTS, data.frame(x=x, type=type, time=time.occurs))
    return (NULL)
  }
  # 
  # Customer arrivals (called by `get.next.event`).
  #
  # Updates the customers pointer `CUSTOMER.COUNT` and returns the customer
  # it newly points to.
  #
  new.customer <- function() {  
    if (CUSTOMER.COUNT < dim(CUSTOMERS)[2]) {
      CUSTOMER.COUNT <<- CUSTOMER.COUNT + 1
      insert.event(CUSTOMER.COUNT, "Customer", 
                   CUSTOMERS["Arrived", CUSTOMER.COUNT])
    }
    return(CUSTOMER.COUNT)
  }
  next.customer.time <- function() {
    if (CUSTOMER.COUNT < dim(CUSTOMERS)[2]) {
      x <- CUSTOMERS["Arrived", CUSTOMER.COUNT]
    } else {x <- Inf}
    return(x) # Time when the next customer will arrive
  }
  #............................................................................#
  #
  # Service.
  #
  find.assistant <- function(time.now) {
    #
    # Select among available assistants.
    #
    j <- which(ASSISTANTS$available <= time.now) 
    #if (length(j) > 0) {
    #  i <- j[ceiling(runif(1) * length(j))]
    #} else i <- NULL                                    # Random selection
    #if (length(j) > 0) i <- j[1] else i <- NULL         # Pick first assistant
    #if (length(j) > 0) i <- j[length(j)] else i <- NULL # Pick last assistant
    if (length(j) > 0) {
      i <- j[which.min(ASSISTANTS[j, ]$available)]
    } else i <- NULL # Pick most-rested assistant
    return (i)
  }#$
  serve <- function(i, x, time.now) {
    #
    # Serve customer `x` with assistant `i`.
    #
    a <- ASSISTANTS[i, ]
    r <- rexp(1, a$rate)                       # Simulate the duration of service
    r <- round(r, 2)                           # (Make simple numbers)
    ASSISTANTS[i, ]$available <<- time.now + r # Update availability
    #
    # Log this successful service event for later analysis.
    #
    CUSTOMERS["Assistant", x] <<- i
    CUSTOMERS["Served", x] <<- time.now
    CUSTOMERS["Duration", x] <<- r
    #
    # Queue the moment the assistant becomes free, so they can check for
    # any customers on hold.
    #
    insert.event(i, "Assistant", time.now + r)
    if (VERBOSE) cat(time.now, ": Assistant", i, "is now serving customer", 
                     x, "until", time.now + r, "\n")
    return (TRUE)
  }
  #............................................................................#
  #
  # The on-hold queue.
  #
  # This is a cicular buffer implemented by an array and two pointers,
  # one to its head and the other to the next available slot.
  #
  put.on.hold <- function(x, time.now) {
    #
    # Try to put customer `x` on hold.
    #
    if (length(HOLD$customers) < 1 || 
          (HOLD$first - HOLD$last %% length(HOLD$customers) == 1)) {
      # Hold queue is full, alas.  Log this occurrence for later analysis.
      CUSTOMERS["Assistant", x] <<- -1 # Busy signal
      if (VERBOSE) cat(time.now, ": Customer", x, "got a busy signal.\n")
      return(FALSE)
    }
    #
    # Add the customer to the hold queue.
    #
    HOLD$customers[HOLD$last] <<- x
    HOLD$last <<- HOLD$last %% length(HOLD$customers) + 1
    if (VERBOSE) cat(time.now, ": Customer", x, "put on hold at position", 
                 (HOLD$last - HOLD$first - 1) %% length(HOLD$customers) + 1, "\n")
    return (TRUE)
  }
  release.hold <- function(time.now) {
    #
    # Pick up the next customer from the hold queue and place them into
    # the event queue.
    #
    if (HOLD$first != HOLD$last) {
      x <- HOLD$customers[HOLD$first]   # Take the first customer
      HOLD$customers[HOLD$first] <<- NA # Update the hold queue
      HOLD$first <<- HOLD$first %% length(HOLD$customers) + 1
      insert.event(x, "Customer", time.now)
    }
  }$
  #............................................................................#
  #
  # Summaries.
  #
  # The CUSTOMERS array contains full information about the customer experiences:
  # when they arrived, when they were served, how long the service took, and
  # which assistant served them.
  #
  summarize <- function() return (list(c=CUSTOMERS, a=ASSISTANTS, e=EVENTS,
                                       h=HOLD))
  #............................................................................#
  #
  # The main event loop.
  #
  initialize(...)
  while (process(get.next.event())) {}
  #
  # Return the results.
  #
  return (summarize())
}
#------------------------------------------------------------------------------#
#
# Specify and run a simulation.
#
set.seed(17)
n.skip <- 200  # Number of initial events to skip in subsequent summaries
system.time({
  r <- sim(n.events=50+n.skip, verbose=TRUE, 
           arrival.rate=60/45, capabilities=1:5/10, hold.queue.size=10)
})
#------------------------------------------------------------------------------#
#
# Post processing.
#
# Skip the initial phase before equilibrium.
#
results <- r$c
ids <- (n.skip+1):(dim(results)[2])
arrived <- results["Arrived", ]
served <- results["Served", ]
duration <- results["Duration", ]
assistant <- results["Assistant", ]
assistant[is.na(assistant)] <- 0   # Was on hold forever
ended <- served + duration
#
# A detailed plot of customer experiences.
#
n.events <- length(ids)
n.assistants <- max(assistant, na.rm=TRUE) 
colors <- rainbow(n.assistants + 2)
assistant.color <- colors[assistant + 2]
x.max <- max(results["Served", ids] + results["Duration", ids], na.rm=TRUE)
x.min <- max(min(results["Arrived", ids], na.rm=TRUE) - 2, 0)
#
# Lay out the graphics.
#
layout(matrix(c(1,1,2,2), 2, 2, byrow=TRUE), heights=c(2,1))
#
# Set up the customers plot.
#
plot(c(x.min, x.max), range(ids), type="n",
     xlab="Time", ylab="Customer Id", main="Customers")
#
# Place points at customer arrival times.
#
points(arrived[ids], ids, pch=21, bg=assistant.color[ids], col="#00000070")
#
# Show wait times on hold.
#
invisible(sapply(ids, function(i) {
  if (!is.na(served[i])) lines(x=c(arrived[i], served[i]), y=c(i,i))
}))
#
# More clearly show customers getting a busy signal.
#
ids.not.served <- ids[is.na(served[ids])]
ids.served <- ids[!is.na(served[ids])]
points(arrived[ids.not.served], ids.not.served, pch=4, cex=1.2)
#
# Show times of service, colored by assistant id.
#
invisible(sapply(ids.served, function(i) {
  lines(x=c(served[i], ended[i]), y=c(i,i), col=assistant.color[i], lty=assistant[i])
}))
#
# Plot the histories of the assistants.
#
plot(c(x.min, x.max), c(1, n.assistants)+c(-1,1)/2, type="n", bty="n",
     xlab="", ylab="Assistant Id", main="Assistants")
abline(h=1:n.assistants, col="#808080", lwd=1)
invisible(sapply(1:(dim(results)[2]), function(i) {
  a <- assistant[i]
  if (a > 0) {
    lines(x=c(served[i], ended[i]), y=c(a, a), lwd=3, col=colors[a+2])
    points(x=c(served[i], ended[i]), y=c(a, a), pch="|", col=colors[a+2])
  }
}))
#
# Plot the customer waiting statistics.
#
par(mfrow=c(1,1))
i <- is.na(served)
plot(served - arrived, xlab="Customer Id", ylab="Minutes",
     main="Service Wait Durations")
lines(served - arrived, col="Gray")
points(which(i), rep(0, sum(i)), pch=16, col="Red")
#
# Summary statistics.
#
mean(!is.na(served)) # Proportion of customers served
table(assistant)