# Simulating continuous time semi-Markov state machine and changing transition probability on the fly

The problem that I'm trying to solve (well, I think that I almost did, but need a review from someone more experienced) is about changing probability of the transition for semi-Markov state machine during simulation.

Let's assume that M=<S,T> defines machine with S states and T probabilistic transitions between states. Elements of T are triplets (s_out, s_in, lambda). So the probability of moving from state s_out to s_in is exponential and defined as exp(-lambda*t).

To simulate it I get random value from [0;1), transform it to exponential and save as next event. What happens is that after saving it but before invoking another event may occur (so it is not a pure Markov state machine) and change the lambda. So the event time should be recalculated.

The problem is that it cannot be just recalculated as such changes may occur quite frequently constantly shifting the event time to the future.

What I think is correct is to adjust recorded event time on ratio time_to_event_from_now * lambda_old / lambda_new. Does it sounds right to you? Maybe the problem description recalls the paper/book where this problem solved?

• I'm not sure I follow you. Why can't the event time just be recalculated? Why have you rejected that -- what undesirable consequences does it have? What is the desired behavior? It's not clear what you want to have happen, when lambda changes. You ask about an alternative formula, but I don't know how to determine whether that is a reasonable solution, because you don't tell us what properties you want the solution to have. How will you evaluate possible solutions? – D.W. Sep 29 '15 at 20:22
• If event is recalculated on every transition change that occurs every day (model time) and the transition time is month ahead, then it will have significant impact on the result. After every recalculation the new event will be one month ahead however system had been working one day already. The desired behaviour should be without this effect. – Alex Netkachov Sep 29 '15 at 21:34
• I'm not clear on why this outcome is wrong/undesirable or why it has any effect (isn't the exponential distribution memoryless? so if you recalculate at any time without changing lambda there's no effect), but can you clarify what the desired behavior is? Do you want the change to lambda to have retroactive effect, and if so, retroactive back to when? You've told us that just recalculating is not what you want, so you've told us what you don't want to happen, but can you describe what you do want to happen? Without that information I'm not seeing how this question can be answered. – D.W. Sep 29 '15 at 21:57
• Let's say that I have engine with two modes of operation: normal and turbo. Probability of failure is different for this modes - for example, turbo is 15 days, normal is 4 months. The modes are changed frequently and 10% engine works in turbo node and 90% in normal. In this case the probability can be calculated analytically, say it is X. Can we represent the system as state machine with two states Ok and Fail and transition between them and another state machine with states normal and turbo that changes this transition probability from ok to fail? – Alex Netkachov Sep 29 '15 at 22:08
• The desired behaviour of these two machines should be similar to analytical (e.g. MTBF obtained by simulation should match MTBF obtained analytically). If we just recalculate transition every time mode is changed and it is changed quite often, the transition will never invoke. – Alex Netkachov Sep 29 '15 at 22:13

As it happens, recalculating the event time is exactly the right thing to do. This follows from the fact that the exponential distribution is memoryless; or from the fact that exponential distribution is the time until the next event in a Poisson process. I realize this might sound miraculous and unbelievable, so let me explain why, in detail.

Let's suppose that at time $t_0$ you have lambda equal to $\lambda_0$; then at time $t_1$, lambda changes to the value $\lambda_1$. In other words, at every point in the time interval $[t_0,t_1)$, the value of lambda is $\lambda_0$; at every point in the time interval $[t_1,t_2)$, the value of lambda is $\lambda_1$.

You are simulating a Poisson process with parameter lambda. During the time interval $[t_0,t_1)$, the Poisson process has parameter $\lambda_0$. During time interval $[t_1,t_2)$, the Poisson process has parameter $\lambda_1$. You want to simulate this process. In particular, you want to simulate the time of the first event.

The brute-force way to simulate this process would be to divide time up into very small increments, say of length $\epsilon$, so that time is broken up into small intervals $[t_0,t_0+\epsilon)$, $[t_0+\epsilon,t_0+2\epsilon)$, $[t_0+2\epsilon,t_0+3\epsilon)$, ..., $[t_2-\epsilon,t_2)$. In each interval, you have a value for lambda: its value is $\lambda_0$ for the first $(t_1-t_0)/\epsilon$ small intervals and $\lambda_1$ for the next $(t_2-t_1)/\epsilon$ small intervals. In each small interval, you draw from a Poisson distribution with parameter $\lambda \epsilon$ (i.e., $\lambda_0 \epsilon$ or $\lambda_1 \epsilon$, depending upon whether the small interval falls before or after $t_1$). Because $\lambda_i \epsilon$ will be a very small value, this is essentially equivalent to flipping a biased coin for each small interval; with probability $\lambda_i \epsilon$, you say that the event has happened in that small interval, and with probability $1-\lambda_i \epsilon$ nothing happens in that small interval.

Thus, this gives you a brute-force strategy to simulate your process. This brute-force strategy will give you the right distribution of events. However, it has the disadvantage that it is computationally expensive: it requires $(t_2-t_0)/\epsilon$ coin flips.

Fortunately, there is an equivalent way to do it that yields the same distribution of events, but is much more computationally efficient. Here is the optimized strategy for how to simulate the time of the first event:

• Randomly draw a value $d$ from the exponential distribution with parameter $\lambda_0$. If $d<t_1-t_0$, then output that the first event happens at time $t_0+d$.

• Otherwise (if $d \ge t_1-t_0$), randomly draw a value $d'$ from the exponential distribution with parameter $\lambda_1$ (independent of $d$). Now, if $d'<t_2-t_1$, then output that the first event happens at time $t_1+d'$; otherwise, output that no event happens during $[t_0,t_2)$.

This optimized strategy lets you simulate the time of the first event. (If you want to simulate the time of all events in the interval $[t_0,t_2)$, you can simulate the time of the first event, if any, and then iterate to figure out the time of the next event if there was a first event, repeating until done.) It is efficient, as you can see.

Moreover, this optimized strategy gives the correct distribution. It is equivalent to the brute-force method described above: it gives the same distribution on the time of the first event (as $\epsilon \to 0$). The equivalence follows from the fact that, if you have a Poisson process with parameter $\lambda$, then the time of the first event is exponentially distributed with parameter $\lambda$. (See Wikipedia.)

Finally, you can note that my optimized strategy is exactly equivalent to recalculating whenever lambda changes. It follows that recalculating whenever lambda changes does give you the right answer. Amazing, huh?

Don't expect this to work for all distributions. This is a special property of the exponential distribution.

• @AlexAtNet, simply recalculating should work fine. It's not true that if there are many switches between $\lambda_0$ and $\lambda_1$ then the event will almost never happen, with the recalculating strategy. I suspect a bug or error in your code. If there's a switch in lambda you need to throw away the old delay and compute a new delay from the new exponential distribution (and the new delay is relative to now, i.e., when the switch happened). – D.W. Sep 30 '15 at 3:22
• play.golang.org/p/0f4APJW7qs - I now understand your solution and it seems to really work fine. Amazing. I need some time to understand that is wrong with the code of my simulation. Thanks! – Alex Netkachov Sep 30 '15 at 10:28
• Just to let you know that I did a lot of tests and it seems that discretisation, recalculation and scaling work equally well :) play.golang.org/p/vGfVJxFA7k Thanks again! – Alex Netkachov Oct 1 '15 at 1:33