Markov models in OpenModelica
/*
BSD-3 licence
Copyright 2022 Felix Felgner
Copyright 2022 Mark Clements
*/
connector Cut
flow Real p_flow;
Real p;
end Cut;
model State
Cut cut;
parameter Real p0 = 0 "initial probability";
Real p(start=p0) "transition probability";
Real l "length of stay";
equation
der(p) = cut.p_flow;
der(l) = p;
cut.p = p;
end State;
model Arc
Cut Start, End;
parameter Real lambda = 1 "transition rate";
Real p_flow;
equation
Start.p_flow + End.p_flow = 0;
Start.p_flow = p_flow;
p_flow = Start.p*lambda;
end Arc;
model Markov1
State s1(p0=1), s2, s3;
Arc arc1(lambda=0.1);
Arc arc2(lambda=0.1);
Arc arc3(lambda=0.05);
equation
// s1 -> s2
connect(arc1.Start,s1.cut);
connect(arc1.End, s2.cut);
// s2 -> s3
connect(arc2.Start,s2.cut);
connect(arc2.End, s3.cut);
// s2 -> s1
connect(arc3.Start,s2.cut);
connect(arc3.End, s1.cut);
end Markov1;
loadFile("simple.mo")
simulate(Markov1, stopTime=10, outputFormat="csv")
simple = read.csv("Markov1_res.csv")
par(mfrow=1:2)
cols = c("s1.p", "s2.p", "s3.p")
matplot(simple$time, simple[,cols], type="l", ylab="Transition probabilities")
legend("topright", legend=cols, lty=1:3, col=1:3, bty="n")
cols = c("s1.l", "s2.l", "s3.l")
matplot(simple$time, simple[,cols], type="l", ylab="Length of stay")
legend("topleft", legend=cols, lty=1:3, col=1:3, bty="n")