Issue with signalling

[Sam Schaffel]

Hello,

I am attempting to simulate an environment that is similar to the example in the State-Dependent Service Rates section of in vignette. It this case, I want the state of those being served by resource 1 to be affected by the number in entities being served by resource 1 and resource 2. Below I have two trajectories, the first being a basic seize and release after a timeout of the two resources and the second is meant to recreate the first trajectory using signalling. I expect both trj_1 and trj_2 to produce the same plot. I test each trajectory by changing the trajectory  argument in the add_generator function.

library(simmer)

library(simmer.plot)

tmax <- 100

#set Timeouts

TIME_OUT <- 5  

TIME_OUT_2  <- 20

INTAKE_TIMES <- at(seq(from = 0, to = tmax, by = 0.1))

#Basic seize/release trajectory with constant timeout

trj_1 <- trajectory() %>%

  seize("resrc")%>%

  timeout(TIME_OUT)%>%

  release("resrc")%>%

  seize("resrc2")%>%

  timeout(TIME_OUT_2)%>%

  release("resrc2")

#Basic Trajectory with signalling and "do nothing" handler

handlr <- trajectory()%>%

  set_attribute(c("start", "delay"), function() {

    left <- sum(get_attribute(envs, c("start", "delay"))) - now(envs)

    delay <- left

    c(now(envs), delay)

  })%>%

  timeout_from_attribute("delay")

trj_2 <- trajectory()%>%

  seize("resrc")%>%

  set_attribute(c("start","delay"), function()

    c(now(envs), TIME_OUT))%>%

  trap("sgnl", handlr)%>%

  send("sgnl")%>%

  untrap("sgnl") %>%

  release("resrc")%>%

  send("sgnl")%>%

  

  seize("resrc2")%>%  #seize resource 2

  send("sgnl")%>%     #update the delay of entities seizing resource 1          

  timeout(TIME_OUT_2)%>%                

  release("resrc2")%>%

  send("sgnl")        

envs <- simmer() 

envs%>%

  add_resource("resrc", Inf, Inf)%>%

  add_resource("resrc2", Inf, Inf)%>%

  

  add_generator("entity", trj_2, distribution = INTAKE_TIMES)%>%

  

  run(until = tmax)%>%

  

  get_mon_resources()%>%

  plot(metric = "usage")

Here is the plot for trajectory 1.

and for trajectory 2.

We can see that trajectory 1 hits the desired steady state of 50 for resource 1 while trajectory 2 does not.

What I have done so far:

1. I tried replicating the problem with only using resource 1. 

Here I was not able to recreate the same problem. I tried sending the signal multiple times after untrapping and releasing resource 1. However the trj_1 and trj_2 plots remained identical with only 1 resource.

2. Removing lines of code to isolate the problem. 

In trajectory 2, if I remove any one of the send("sgnl") commands after releasing resource 1 or the timeout(TIME_OUT_2) command. It results in a plot identical to trajectory 1. I also discovered here that the duration of TIME_OUT_2 affects the resource 1 usage curve.

3. I tried changing the INTAKE_TIMES data. Here I discovered that if the increment of the arrival sequence is less than or equal to 0.1 then the two plots will differ. But if the increment is greater than 0.1 then the plots for trajectory 1 and trajectory 2 will be identical.

4. I tried logging each entity as they went through the handler. Here I discovered an interesting pattern. The first 10 entities sent each time send("sgnl") is called appear to be entities from earlier in the simulation which should have been untrapped and unsubscribed from the signal. Every 5th entity from 274 to 319 seems to have not unsubscribed.

  

99.9: entity274: in the handler

99.9: entity279: in the handler

99.9: entity284: in the handler

99.9: entity289: in the handler

99.9: entity294: in the handler

99.9: entity299: in the handler

99.9: entity304: in the handler

99.9: entity309: in the handler

99.9: entity314: in the handler

99.9: entity319: in the handler

99.9: entity949: in the handler

99.9: entity950: in the handler

99.9: entity951: in the handler

99.9: entity952: in the handler

...

99.9: entity996: in the handler

99.9: entity997: in the handler

[Iñaki Ucar]

Hi, thanks for a detailed assessment, comments below:

On Wed, 22 Jul 2020 at 23:46, Sam Schaffel <sam.sch...@gmail.com> wrote:
>
> In trajectory 2, if I remove any one of the send("sgnl") commands after releasing resource 1 or the timeout(TIME_OUT_2) command. It results in a plot identical to trajectory 1.

This is the main issue: there are two "simultaneous" (because there is
no timeout between them) send() activities. This typically leads to
undefined behaviour, because it's impossible to order the events
properly anymore.

More specifically, the second consecutive send() is triggering a
secondary handler when the arrivals didn't reach the
timeout_from_attribute yet, and thus they end up executing *two*
timeouts instead of one.

TL;DR, the send() activity between resources is superfluous, and the
solution is to remove it, as you found.

--
Iñaki Úcar

[Iñaki Ucar]

Actually, on second thought, if the number of arrivals in resource 1
is going to affect the timeout too, it's not superfluous anymore,
because as soon as resource 2 is full, the arrival will end up in its
queue without executing any send()...

So the most general solution is to insert timeout(0) after the first
send() and before seizing the second resource. The timeout() activity
has the lowest priority, so even a null timeout effectively splits up
a simultaneous event chain.

--
Iñaki Úcar

[Sam Schaffel]

Thank you for the response Iñaki,

This solution worked perfectly for the example I gave. Yet somehow makes the problem worse in my original program, perhaps the problem is something else entirely. Here I think I narrowed the problem down to when I had stochastic effects to the rates and initial conditions to the resources. Interestingly enough, in this version of my code, the problem seems to be in reverse. The signalling trajectory matches the basic trajectory only when all 4 send() cmds are including, and otherwise the graph deviates. Adding the timeout(0) in this case changes the behavior dramatically.

library(simmer)

library(simmer.plot)

x0 <- c(25, 100)

xtarget <- c(50, 200) 

c <- 0.05 

a <- c * xtarget[2]

b <- a / xtarget[1] 

tmax <- 100

TIME_OUT <- rexp(1, b)

TIME_OUT_2 <- function()  rexp(1, c)

INTAKE_TIMES <- function() rexp(1, a)

#List of test trajectories

#Basic seize/release trajectory with contant timeout

trj_1.2 <- trajectory()%>%

  seize("resrc2")%>%

  timeout(TIME_OUT_2)%>%

  release("resrc2")

trj_1 <- trajectory() %>%

  seize("resrc")%>%

  timeout(function() rexp(1, b))%>%

  release("resrc")%>%

  join(trj_1.2)

#Basic Trajectory with signalling and "do nothing" handler

handlr <- trajectory()%>%

  set_attribute(c("start", "delay"), function() {

    left <- sum(get_attribute(envs, c("start", "delay"))) - now(envs)

    delay <- left

    c(now(envs), delay)

  })%>%

  timeout_from_attribute("delay")

trj_2.2 <- trajectory()%>%

  seize("resrc2")%>%

  send("sgnl")%>%               

  timeout(TIME_OUT_2)%>%              

  release("resrc2")%>%

  send("sgnl")

trj_2 <- trajectory()%>%

  seize("resrc")%>%

  set_attribute(c("start","delay","multiplier"), function()

    c(now(envs), rexp(1, b), 1))%>%

  trap("sgnl", handlr)%>%

  send("sgnl")%>%

  untrap("sgnl") %>%

  release("resrc")%>%

  send("sgnl")%>%   #comment out this line for picture 3

  #timeout(0)%>%    #uncomment this line for picture 4

  join(trj_2.2)

envs <- simmer() 

envs%>%

  add_resource("resrc", Inf, Inf)%>%

  add_resource("resrc2", Inf, Inf)%>%

  

  add_generator("initial resrc", trj_1.2, at(numeric(x0[2]))) %>%

  add_generator("initial resrc2", trj_1, at(numeric(x0[1]))) %>% 

  add_generator("entity", trj_1, distribution = INTAKE_TIMES)%>%

  run(until = tmax)

 

 

resrcs <- get_mon_resources(envs)

plot(resrcs, metric = "usage")

This code results in the following graphs:

Trajectory 1.

Trajectory 2 (as in the code):

Trajectory two (with one of the send() cmds removed):

Trajectory 2 (with timeout(0)):

I am confused why these 3 implementations of trajectory 2 produce such different graphs, since I expect the handler and timeout(0) to essentially do nothing.

Finally, the main difference between that code and my program is the handler. The objective of the handler is to keep the overall rate of entities leaving resource 1 to be "constant" when it reaches a saturation point. In other words, when the number of entities served by resource 1 is below 1/3 * the number of entities served by resource 2 I want the behaviour as above. But when the number of entities served by resource 1 reaches or exceeds that saturation point, I want the overall rate of entities releasing resource 1 to no longer depend on the number of entities served by resource 1. Instead I want the over rate to be 1/3* the current number of entities served by resource 2. I also want each entity in resource 1 to receive equal attention from the resource. i.e if a new entity enters the trajectory when it is above the saturation point, all of the entities served by resource 1 must have their timeout extended by the same factor.

I apologize if the is confusing.

The new handler in below:

handlr2 <- trajectory() %>%

  set_attribute(c("start", "delay", "multiplier"), function() {

    cap <- 1/3*get_server_count(envs, "resrc2")

    multiplier <- get_attribute(envs, "multiplier")

    num_mentees <- get_server_count(envs, "resrc")

    left <- sum(get_attribute(envs, c("start", "delay"))) - now(envs)

    new_multiplier <- cap/num_mentees

    #return new values

    if (num_mentees <= cap)

      delay <- left

    else

      delay <- left*multiplier / new_multiplier

    c(now(envs), delay, new_multiplier)

  })%>%

  timeout_from_attribute("delay")

I've also modeled this simulation using just a Monte Carlo method. But I wanted to be able to track the individual entities in the system which is why I'm reimplementing using simmer. However the plots for the Simmer and Monte Carlo methods do not overlap. Below is a graph showing the mean of the Simmer and Monte Carlo programs after being run 100 times. 

Is there a better way to achieve what I want?

Thanks,

Sam

[Iñaki Ucar]

On Mon, 27 Jul 2020 at 16:45, Sam Schaffel <sam.sch...@gmail.com> wrote:
>
> Thank you for the response Iñaki,
> This solution worked perfectly for the example I gave. Yet somehow makes the problem worse in my original program, perhaps the problem is something else entirely. Here I think I narrowed the problem down to when I had stochastic effects to the rates and initial conditions to the resources. Interestingly enough, in this version of my code, the problem seems to be in reverse. The signalling trajectory matches the basic trajectory only when all 4 send() cmds are including, and otherwise the graph deviates. Adding the timeout(0) in this case changes the behavior dramatically.

Of course, silly me... timeout() has the lowest priority, but there is
a timeout_from_attribute() in the handler, which *also* has the lowest
priority. Therefore, we have the correct ordering sometimes, sometimes
not. So, this may sound silly, but add *another* timeout(0), like
this:

trj_2 <- trajectory()%>%
seize("resrc")%>%
set_attribute(c("start","delay","multiplier"), function()
c(now(envs), rexp(1, b), 1))%>%
trap("sgnl", handlr)%>%
send("sgnl")%>%
untrap("sgnl") %>%
release("resrc")%>%
send("sgnl")%>%

timeout(0)%>%
timeout(0)%>%
join(trj_2.2)

And now you should get the same* result, because now the second
timeout(0) is executed necessarily after any timeout_from_attribute()
in the handlers. In the future, we should try to rework this because I
don't like this workaround.

*You won't get exactly the same result unless you load dataframes with
all your arrival and service times and then use add_dataframe instead
of add_generator. In this way, you ensure that all values are
generated exactly with the same random state.

> Finally, the main difference between that code and my program is the handler. The objective of the handler is to keep the overall rate of entities leaving resource 1 to be "constant" when it reaches a saturation point. In other words, when the number of entities served by resource 1 is below 1/3 * the number of entities served by resource 2 I want the behaviour as above. But when the number of entities served by resource 1 reaches or exceeds that saturation point, I want the overall rate of entities releasing resource 1 to no longer depend on the number of entities served by resource 1. Instead I want the over rate to be 1/3* the current number of entities served by resource 2. I also want each entity in resource 1 to receive equal attention from the resource. i.e if a new entity enters the trajectory when it is above the saturation point, all of the entities served by resource 1 must have their timeout extended by the same factor.
> I apologize if the is confusing.
> The new handler in below:
>
> handlr2 <- trajectory() %>%
> set_attribute(c("start", "delay", "multiplier"), function() {
> cap <- 1/3*get_server_count(envs, "resrc2")
> multiplier <- get_attribute(envs, "multiplier")
> num_mentees <- get_server_count(envs, "resrc")
> left <- sum(get_attribute(envs, c("start", "delay"))) - now(envs)
> new_multiplier <- cap/num_mentees
> #return new values
> if (num_mentees <= cap)
> delay <- left
> else
> delay <- left*multiplier / new_multiplier
> c(now(envs), delay, new_multiplier)
> })%>%
> timeout_from_attribute("delay")
>
> I've also modeled this simulation using just a Monte Carlo method. But I wanted to be able to track the individual entities in the system which is why I'm reimplementing using simmer. However the plots for the Simmer and Monte Carlo methods do not overlap. Below is a graph showing the mean of the Simmer and Monte Carlo programs after being run 100 times.
>

> Is there a better way to achieve what I want?

Please, try the workaround above with this handler.

--
Iñaki Úcar

[Sam Schaffel]

Thanks again for such a quick response,

I tried the workaround you suggested. The result I got was that 

trj_2 <- trajectory()%>%

  seize("resrc")%>%

  set_attribute(c("start","delay","multiplier"), function()

    c(now(envs), rexp(1, b), 1))%>%

  trap("sgnl", handlr)%>%

  send("sgnl")%>%

  untrap("sgnl") %>%

  release("resrc")%>%

  send("sgnl")%>%

  timeout(0)%>%

  timeout(0)%>%

  

  join(trj_2.2)

is equivalent to 

trj_2 <- trajectory()%>%

  seize("resrc")%>%

  set_attribute(c("start","delay","multiplier"), function()

    c(now(envs), rexp(1, b), 1))%>%

  trap("sgnl", handlr)%>%

  send("sgnl")%>%

  untrap("sgnl") %>%

  release("resrc")%>%

  # send("sgnl")%>%

  # timeout(0)%>%

  # timeout(0)%>%

  

  join(trj_2.2)

However now both of which result in resource 1 having a steady state somewhere between 65-70 (as opposed to 50 from trajectory 1).

The following gives a plot with steady state 50, however I don't actually want to send the signal when an entity is in between resources, since in my simulation when an entity is done with resource 1 it immediately goes to resource 2 with no time delay. And I wouldn't want to update the delay for entities served by resource 1 based on a state that *shouldn't exist*.

trj_2 <- trajectory()%>%

  seize("resrc")%>%

  set_attribute(c("start","delay","multiplier"), function()

    c(now(envs), rexp(1, b), 1))%>%

  trap("sgnl", handlr)%>%

  send("sgnl")%>%

  untrap("sgnl") %>%

  release("resrc")%>%

  send("sgnl")%>%

  # timeout(0)%>%

  # timeout(0)%>%

  

  join(trj_2.2)

I tried the workaround with handlr2. Again I find that the usage curve for resource 1 is higher than I expected/wanted.

Thanks again for dedicating your time to helping me,

Sam

[Iñaki Ucar]

They are as long as resource 2 has enough capacity. If you are going
to simulate with infinite capacity, that's fine. If not, then I would
suggest using send > timeout(0) > timeout(0) instead.

> However now both of which result in resource 1 having a steady state somewhere between 65-70 (as opposed to 50 from trajectory 1).

Maybe a single consecutive run. But run many replications and you'll
see that on average it's completely equivalent to trajectory 1.

> I tried the workaround with handlr2. Again I find that the usage curve for resource 1 is higher than I expected/wanted.

You mean compared to the MC method? Are you taking a single simmer
simulation or the average of many replications?

--
Iñaki Úcar

[Sam Schaffel]

Hello,

Maybe a single consecutive run. But run many replications and you'll
see that on average it's completely equivalent to trajectory 1.

I tried 35 replications and got the following result:

Here all 25 seem to reach a steady state of above 60.

You mean compared to the MC method? Are you taking a single simmer 
simulation or the average of many replications? 

Yes compared to the MC method. The original plot I sent comparing Simmer to MC was the mean over 100 iterations. When testing the workaround I plotted the mean, median and quartiles against the MC for 3 iterations (because the run time is quite long with all the signalling). Here the difference was too vast for me to try again with 100 iterations (this takes about an hour for me).

Thanks,

Sam

[Iñaki Ucar]

On Mon, 27 Jul 2020 at 23:31, Sam Schaffel <sam.sch...@gmail.com> wrote:
>
> Hello,
>
>> Maybe a single consecutive run. But run many replications and you'll
>> see that on average it's completely equivalent to trajectory 1.
>
> I tried 35 replications and got the following result:

> Here all 25 seem to reach a steady state of above 60.

That is not what I get. Code and plot attached.

--
Iñaki Úcar

[Sam Schaffel]

Thanks Iñaki, the code you sent seems to work well.

I do have one more question. The only difference between your code and the one I have is where the initial entities at time 0 are sent.

I have

    add_generator("initial resrc", trj_2.2, at(numeric(x0[2]))) %>%

    add_generator("initial resrc2", trj_2, at(numeric(x0[1]))) %>%

    add_generator("entity", trj_2, distribution = INTAKE_TIMES)%>%

while you have

    add_generator("initial resrc", trj_1.2, at(numeric(x0[2]))) %>%

    add_generator("initial resrc2", trj_1, at(numeric(x0[1]))) %>%

    add_generator("entity", trj_2, distribution = INTAKE_TIMES)%>%

This means that the initial entities do not update the state of the system when they change resources, and the initial entities in resource 1 never get sent to the handler to have their delay updated.. However in my version of the code the plot is completely different. I assume this is because multiple signals are being sent simultaneously at time 0. Is there any way around that?

Thanks,

Sam

[Iñaki Ucar]

On Tue, 28 Jul 2020 at 20:43, Sam Schaffel <sam.sch...@gmail.com> wrote:
>
> Thanks Iñaki, the code you sent seems to work well.
> I do have one more question. The only difference between your code and the one I have is where the initial entities at time 0 are sent.
> I have
> add_generator("initial resrc", trj_2.2, at(numeric(x0[2]))) %>%
> add_generator("initial resrc2", trj_2, at(numeric(x0[1]))) %>%
> add_generator("entity", trj_2, distribution = INTAKE_TIMES)%>%
> while you have
> add_generator("initial resrc", trj_1.2, at(numeric(x0[2]))) %>%
> add_generator("initial resrc2", trj_1, at(numeric(x0[1]))) %>%
> add_generator("entity", trj_2, distribution = INTAKE_TIMES)%>%
>
> This means that the initial entities do not update the state of the system when they change resources, and the initial entities in resource 1 never get sent to the handler to have their delay updated.. However in my version of the code the plot is completely different. I assume this is because multiple signals are being sent simultaneously at time 0. Is there any way around that?

Oh, I see. I thought you wanted the initial arrivals to have fixed
times. Then, yes, the issue is all the stuff going on at the same time
at t=0. The workaround is to randomize those arrivals too.

--
Iñaki Úcar

[Sam Schaffel]

Just to clarify, I can't have multiple arrivals at time 0 subscribe to the signal and have their delay updated?

[Iñaki Ucar]

I don't think so. Arrivals execute their handler multiple times, thus
ending up executing several timeouts (and that's why you see a higher
occupation over time).

> --
> You received this message because you are subscribed to the Google Groups "simmer-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to simmer-devel...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/simmer-devel/6a910358-9f2b-4313-9cc4-120d4a693678o%40googlegroups.com.



--
Iñaki Úcar

[Sam Schaffel]

I see, thanks for all the help!

Sam

> To unsubscribe from this group and stop receiving emails from it, send an email to simmer...@googlegroups.com.

[Sam Schaffel]

Hi Iñaki,

I have to admit I'm still quite a bit confused what is going on behind the scenes here, and I'm hoping to understand a little better. I've narrowed the problem I'm having even further. Here I have just one resource and I only include the 25 initial conditions at time 0. All of the entities should leave the system at time 5. However as seen from the graph, most of the entities do not get untrapped at t = 5. By tracking the delays, I can see that after t = 5, the delay attribute is calculated to be negative (because it's calculating when it should have left the system from the current time). Then when timeout_from_attribute() is called, the delay is coerced to positive, and so these entities never leave the system. 

Uncommenting the add_dataframe() command  sends more signals so it is easier to see how these delays are being calculated. Looking at the attrbs dataframe i see the delay is actually because calculating correctly up until t = 5 for these initial entities (delay + current time = 5). 

This behaviour is fixed by uncommenting either the first send() command in trj or uncommenting the delay <- max(0,delay) line in the handler.

My question is why does adding another send() command at t = 0 affect whether or not these entities get untrapped 5 time units later? And why do the entities stay stuck in the system without that send() command?

library(simmer)

library(simmer.plot)

tmax <- 100

handlr <- trajectory()%>%

  set_attribute(c("start", "delay"), function() {

    left <- sum(get_attribute(envs, c("start", "delay"))) - now(envs)

    delay <- left

    #delay <- max(0,delay)

    c(now(envs), delay)

  })%>%

  timeout_from_attribute("delay")

trj <- trajectory()%>%

  #send("sgnl") %>%

  seize("resrc")%>%

  set_attribute(c("start","delay") , function()

    c(now(envs), 5))%>%

  trap("sgnl", handlr)%>%

  send("sgnl")%>%

  untrap("sgnl") %>%

  release("resrc")

envs <- simmer() 

data <- data.frame(time = replicate(1000, rexp(1, 10)))

init <- replicate(25, 0)

envs%>%

  

  add_resource("resrc", Inf, Inf)%>%

  add_generator("initial entities", trj, at(init), mon = 2) %>%

  

  #add_dataframe("signallers", trj, data)%>%

  run(until = tmax)

resrcs <- get_mon_resources(envs)

plot(resrcs, metric = "usage")

attrbs <- get_mon_attributes(envs)

plot(attrbs)

Thanks again,

Sam

[Iñaki Ucar]

As we discussed in previous examples, simultaneous arrivals are not
going to work in general with this pattern, because there are so many
things happening at the same time, and event order is not guaranteed
anymore. In other words, signaling with simultaneous events causes
undefined behaviour. But this is a general problem of DES, and the
only way to properly fix it is to avoid simultaneous arrivals
altogether.

And if you really really want to know what's going on, you can set
`envs <- simmer(verbose=TRUE)`.

> To unsubscribe from this group and stop receiving emails from it, send an email to simmer-devel...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/simmer-devel/e58f52ff-4812-4662-8dff-5dc4121c4900o%40googlegroups.com.



--
Iñaki Úcar