Time-Varrying Source Doesn't Like Small Standard Deviation

I'm trying to simulate the flow of items into our process. It has a low level constant flow of items throughout the day except for a 3 to 4 hour period where it peaks like a gaussian distribution. I want to be able to control the peak amplitude and std deviation and have this model throw Flow Items out into our other modelled sections.

So I have two Sources, one that sends out a continuous, slow output throughout the day, which I've temporarilly disabled by putting its Inter-Arrivaltime to a grossly high number. This lets me concentrate on the other "Peak" Source. I finally wrote the equation for the Inter-Arrivaltime by hand because I couldn't find anything else to work for me. I used the standard A * e ^ [ (t - t_pk)^2 / std_dev^2 ] equation and used the "time( )" function for the input (t). I created some Labels on the Source to specify the amplitude (tubes/hr), center point (hr of the day that the peak happens), and the std_dev of the center point. Some of the constants in the equation are computed at Reset so that it doesn't try to calculate them every update....it just uses the Label variables, does some simple calculations and save the results in global variables...then the InterArrivaltime uses the global variables during the simulation.

Anyway, I have the peak centered at the 12th hour (43200 seconds). I also have a chart that looks at the last few 5 minutes segments to get an average "items per hour" that are arriving and plots those. The graph looks right so long as I specify a std dev of 3 hours or more. If I go to 2.9 or 2.8 for the std dev, it puts out zero Flow Items until about 1/4th to 1/2 of the way into the simulation and then it starts putting them out as expected. If I specify a std dev of 2.5, it never puts anything out. It's like the equation for the Source is "clipping" or something????

Here's the output (items/hour) when the std dev is 3 hours:



And here it is when the std dev is 2.8 hours:


I've looked over the math several times now and don't see anything wrong. Can anyone tell me what's going on? (If you run it, you might want to put a Stop Time of 86400 seconds and crank up the Run Speed).
Thanks,
David White
[email protected]

Put some print statements into your interarrival time code to test your custom distribution and you can see what is going on:

pt("INTERARRIVALTIME\n");
pt("g_PeakAmpl ");pf(g_PeakAmpl);pr();
pt("time() ");pf(time());pr();
pt("(time()/3600.0) ");pf((time()/3600.0));pr();
pt("(time()/3600.0) - g_PeakCenter ");pf((time()/3600.0) - g_PeakCenter);pr();
pt("sqr((time()/3600.0) - g_PeakCenter) ");pf(sqr((time()/3600.0) - g_PeakCenter));pr();
pt("pow(2.71828, sqr((time()/3600.0) - g_PeakCenter) / g_Peak_2SigmaSq) ");pf(pow(2.71828, sqr((time()/3600.0) - g_PeakCenter) / g_Peak_2SigmaSq));pr();
pt("( (g_PeakAmpl * pow(2.71828, sqr((time()/3600.0) - g_PeakCenter) / g_Peak_2SigmaSq) ) ) ");pf(( (g_PeakAmpl * pow(2.71828, sqr((time()/3600.0) - g_PeakCenter) / g_Peak_2SigmaSq) ) ));pr();



INTERARRIVALTIME
g_PeakAmpl 3.600000
time() 0.000000
(time()/3600.0) 0.000000
(time()/3600.0) - g_PeakCenter -12.000000
sqr((time()/3600.0) - g_PeakCenter) 144.000000
pow(2.71828, sqr((time()/3600.0) - g_PeakCenter) / g_Peak_2SigmaSq) 9736.795048
( (g_PeakAmpl * pow(2.71828, sqr((time()/3600.0) - g_PeakCenter) / g_Peak_2SigmaSq) ) ) 35052.462172
INTERARRIVALTIME
g_PeakAmpl 3.600000
time() 35052.462172
(time()/3600.0) 9.736795
(time()/3600.0) - g_PeakCenter -2.263205
sqr((time()/3600.0) - g_PeakCenter) 5.122097
pow(2.71828, sqr((time()/3600.0) - g_PeakCenter) / g_Peak_2SigmaSq) 1.386336
( (g_PeakAmpl * pow(2.71828, sqr((time()/3600.0) - g_PeakCenter) / g_Peak_2SigmaSq) ) ) 4.990809

The first time it calls it, when the time() is 0, you subtract 12 (g_PeakCenter), then square that, giving you 144, which is significantly higher than 15.68 (g_Peak_2SigmaSq) so when you divide it, you get a high number instead of the expected number between 0 and 1. So when you take 2 to that huge number power, you get a huge number.

Your equation returns exponentially larger values when time() == 0 as the PeakStdDev label is set lower.

3.00 returns 10731.391 for the first item to arrive
2.8 returns 35052.462
2.6 returns 152027.128
2.5 returns 362553.053

The Source isn't "clipping" anything. It fires the interarrival code at time 0 for the first item to arrive and your code is returning ridiculously large values for the first item to arrive.

I see what you're saying, but at time 0, I *want* a huge number for the inter-arrival time because I don't want any Flow Items from that Source around time 0....or around time 24hr (86400 seconds) for that matter. I only want Flow Items around the 12 hour mark. And when time is near that (time ~ 43,000), the exponent will start getting closer and closer to 0, and the exponential-function will start getting closer and closer to 1, and then the Flow Items will start showing up closer and closer to the rate specified by PeakAmpl. I hope I explained that well.

In case I didn't, explain to me why having the std dev at 3 (g_Peak_2SigmaSq = 18) gives a good/right gaussian shape, but having the std dev of 2.8 (g_Peak_2SigmaSq = 15.68) doesn't? In fact, not only is it not what is expected, but it's not even symmetric, which it should be because it's an even function. There's something else going wrong.

<<< hope I wasn't editting this reply while you were reading it >>>

Thanks,
David White

Put this in before your return statement if you want it to start around around a certain time period. Your equation doesn't work the way you want it to at time 0.

if(time()==0) return 10000;

David,

Here is a model that has some code on the source that allows you to specify different inter arrival times according to the time of day. You could take this source and adjust the code so that it looked for a distribution as opposed to a static time value. This might help you get around the problems you are having with the two sources and getting the distribution to behave the way you want.

Good Luck,
Brandon

Thanks for both of your comments. Brandon, I found that method when I did a search originally....I like that ability and will keep it in mind for the future, but I didn't think it gave me the resolution and ease. If i wanted to adjust the shape of the bump, I'd have to edit the table every time. But I just want to specify the height (items/hour) and the width (std dev).

Phil, that's interesting...I do see a difference, but I still don't understand why. I actually made it a little different that what you suggested. I now return "min(5000, << exponential eq. >>)". That gets it looking much better. It still has a weird artifact in it though as seen below (this has the std_dev set to 0.5):



The left side still has a sharp rising...not like the right side where it asymtopically drops down to 0. Overall I think this will work. The weird effects are down in the noise...especially when I turn on the other source. So thanks for your help...but I'd still like to know why it's not symmetric. The function is returning the same values for times to the left and right of the 12h mark....so the rate of deliveries should show up the same but they don't.

Anyway, thanks again,
David White

Hi,

What I would do when I get in this type of problem is to record the time the inter arrival triggers
and the value of your formula, export this information to Excel.
In Excel I would let Excel calculate the formula and compare it with the result from Flexsim.
If they are the same Flexsim is not probably doing anything wrong, it is something with your formula.

Lars-Olof

The problem is the discrete nature of the simulation, and your formula is a very continuous formula. Your formula will calculate the correct interarrival time (1/arrival rate) for any time. The problem is that when the standard deviation is 2.8, the peak is at t=12 hours (over 4 standard deviations away), you get a very low arrival rate and a very high interarrival time. This tells the source to create an item at time 35052.462172. It does not recalculate the interarrival time until this first item is created, at which time the interarrival time is much lower, so it looks discontinuous. You have this same problem with a standard deviation of 3, but the first interarrival time is far less, so the discontinuity isn't as big of a jump so you don't notice it.

A solution (maybe not the best): force it to create an item and recalculate the interarrival time in a more continuous manner, at max 1000 time units, like this:

return /**/        min(1000,(  (g_PeakAmpl * pow(2.71828, sqr((time()/3600.0) - g_PeakCenter) / g_Peak_2SigmaSq) )     ))/**/;

Here's a different, possibly better solution. Calculate the interarrival time now and at time()+(interarrival time) and return the average of the two. When there is going to be a big jump in the interarrival times, it will decrease the interarrival time a little to recalculate it several times before the jump. If there isn't a jump, this won't affect it much.

 double interArrivalTimeAtCurrentTime=g_PeakAmpl * pow(2.71828, sqr((time()/3600.0) - g_PeakCenter) / g_Peak_2SigmaSq);
double interArrivalTimeAtNextTime=g_PeakAmpl * pow(2.71828, sqr(((time()+interArrivalTimeAtCurrentTime)/3600.0) - g_PeakCenter) / g_Peak_2SigmaSq);
return (interArrivalTimeAtCurrentTime+interArrivalTimeAtNextTime)/2;

That makes sense..... and on the other side of the hump, the inter-arrival time is getting larger each time, so it's getting re-computed faster than they arrive, therefore that side looks smooth and right.


Thanks Alex.
Dave