A Markov Chain is a weighted digraph representing a discrete-time system that can be in any number of discrete states.  The nodes of the digraph represent the states, and the directed edge weight between two states a and b represents the probability (called the transition probability from a to b) that the system will move to state b in the next time period, given that it is currently in state a.  The sum of the transition probabilities out of any node is, by definition, 1.  The set of probabilities is stored in a transition matrix P, where entry (i, j) is the transition probability from state i to state j. Clearly, the sum of each row of P is 1.  

A well-known theorem of Markov chains states that the probability of the system being in state j after k time periods, given that the system begins in state i, is the (i, j) entry of $P^k$.

A common question arising in Markov-chain models is, what is the long-term probability that the system will be in each state?  The vector containing these long-term probabilities, denoted $\mathrm{Pi}$, is called the steady-state vector of the Markov chain.

This Maple application creates a procedure for answering this question.  As a case study, we'll analyze a two-server computer network whose servers have known probabilities of going down or being fixed in any given hour. The goal is to compute the long-term probability that at least one server is working.

We create a Maple procedure called steadyStateVector that takes as input the transition matrix of a Markov chain and returns the steady state vector, which contains the long-term probabilities of the system being in each state.  The input transition matrix may be in symbolic or numeric form.  

The procedure steadyStateVector implements the following algorithm:  Given an n x n transition matrix P, let I be the n x n identity matrix and Q = P - I.  Let e be the n-vector of all 1's, and b be the (n+1)-vector with a 1 in position n+1 and 0 elsewhere.  To compute the steady state vector, solve the following linear system for $\mathrm{Pi}$, the steady-state vector of the Markov chain:

${\mathrm{(Q\; |\; e)}}^T\mathrm{Pi}\=b$

Appending e to Q, and a final 1 to the end of the zero-vector on the right-hand side ensures that the solution vector $\mathrm{Pi}$ has components summing to 1.

Here is the steadyStateVector procedure.  The input is a transition matrix P, and the output is the steady-state vector $\mathrm{pi}$ reflecting the long-term probability of the system being in each state.  Comments are preceded by the # symbol.

The problem is to estimate the long-term probability that at least one server in a two-server computer network is working during any given hour.  We'll model this problem as a Markov chain as follows:  

Assume the network can be in one of three states:

Let:

The transition matrix is then the following:

$P≔\left[\begin{array}{ccc}1-{\mathrm{\mu }}_1& {\mathrm{\mu }}_1& 0\\ {\mathrm{\lambda }}_1& 1-{\mathrm{\lambda }}_1-{\mathrm{\mu }}_2& {\mathrm{\mu }}_2\\ 0& {\mathrm{\lambda }}_2& 1-{\mathrm{\lambda }}_2\end{array}\right]$

Note that the steadyStateVector procedure computes $\mathrm{pi}$ symbolically.  Numerical values for $\mathrm{lambda}$ and $\mathrm{mu}$ are not required.

$\mathrm{\π}≔\left[\begin{array}{c}\frac{{\mathrm{\lambda }}_1{\mathrm{\lambda }}_2}{{\mathrm{\lambda }}_1{\mathrm{\lambda }}_2+{\mathrm{\lambda }}_2{\mathrm{\mu }}_1+{\mathrm{\mu }}_2{\mathrm{\mu }}_1}\\ \frac{{\mathrm{\lambda }}_2{\mathrm{\mu }}_1}{{\mathrm{\lambda }}_1{\mathrm{\lambda }}_2+{\mathrm{\lambda }}_2{\mathrm{\mu }}_1+{\mathrm{\mu }}_2{\mathrm{\mu }}_1}\\ \frac{{\mathrm{\mu }}_2{\mathrm{\mu }}_1}{{\mathrm{\lambda }}_1{\mathrm{\lambda }}_2+{\mathrm{\lambda }}_2{\mathrm{\mu }}_1+{\mathrm{\mu }}_2{\mathrm{\mu }}_1}\end{array}\right]$

$\mathrm{pi}$ is currently a symbolic vector. Let's supply some example values:

${\mathrm{\lambda }}_1≔\frac{1}{400}$

${\mathrm{\lambda }}_2≔\frac{1}{800}$

${\mathrm{\mu }}_1≔\frac{1}{4}$

${\mathrm{\mu }}_2≔\frac{1}{8}$

Ask Maple for the value of $\mathrm{pi}$, and the updated vector is given, reflecting the inputs above.  By inspection, we see the vector $\mathrm{pi}$ is stochastic.

$\mathrm{\π}≔\left[\begin{array}{c}\frac{1}{10101}\\ \frac{100}{10101}\\ \frac{10000}{10101}\end{array}\right]$

The long-term probability that at least one server is operable in any given hour is the sum of the last two components of $\mathrm{pi}$.

$\mathrm{probWorking}≔\frac{10100}{10101}$

If we want a 10-digit floating-point approximation to this probability, use the evalf command.

$0.9999009999$

Erase the current values of lambda and mu.

$\mathrm{\lambda }≔\mathrm{\lambda }$

$\mathrm{\mu }≔\mathrm{\mu }$

Let's generalize this Markov chain, allowing for the possibility of both servers going down at once or both being repaired at once.

$P≔\left[\begin{array}{ccc}1-{\mathrm{\mu }}_1-{\mathrm{\mu }}_3& {\mathrm{\mu }}_1& {\mathrm{\mu }}_3\\ {\mathrm{\lambda }}_1& 1-{\mathrm{\lambda }}_1-{\mathrm{\mu }}_2& {\mathrm{\mu }}_2\\ {\mathrm{\lambda }}_3& {\mathrm{\lambda }}_2& 1-{\mathrm{\lambda }}_2-{\mathrm{\lambda }}_3\end{array}\right]$

$\mathrm{\π}≔\left[\begin{array}{c}\frac{{\mathrm{\lambda }}_2{\mathrm{\lambda }}_1+{\mathrm{\lambda }}_3{\mathrm{\lambda }}_1+{\mathrm{\mu }}_2{\mathrm{\lambda }}_3}{{\mathrm{\lambda }}_2{\mathrm{\lambda }}_1+{\mathrm{\lambda }}_3{\mathrm{\lambda }}_1+{\mathrm{\lambda }}_1{\mathrm{\mu }}_3+{\mathrm{\mu }}_1{\mathrm{\lambda }}_2+{\mathrm{\mu }}_3{\mathrm{\lambda }}_2+{\mathrm{\mu }}_1{\mathrm{\lambda }}_3+{\mathrm{\mu }}_2{\mathrm{\lambda }}_3+{\mathrm{\mu }}_2{\mathrm{\mu }}_1+{\mathrm{\mu }}_2{\mathrm{\mu }}_3}\\ \frac{{\mathrm{\mu }}_1{\mathrm{\lambda }}_2+{\mathrm{\mu }}_3{\mathrm{\lambda }}_2+{\mathrm{\mu }}_1{\mathrm{\lambda }}_3}{{\mathrm{\lambda }}_2{\mathrm{\lambda }}_1+{\mathrm{\lambda }}_3{\mathrm{\lambda }}_1+{\mathrm{\lambda }}_1{\mathrm{\mu }}_3+{\mathrm{\mu }}_1{\mathrm{\lambda }}_2+{\mathrm{\mu }}_3{\mathrm{\lambda }}_2+{\mathrm{\mu }}_1{\mathrm{\lambda }}_3+{\mathrm{\mu }}_2{\mathrm{\lambda }}_3+{\mathrm{\mu }}_2{\mathrm{\mu }}_1+{\mathrm{\mu }}_2{\mathrm{\mu }}_3}\\ \frac{{\mathrm{\lambda }}_1{\mathrm{\mu }}_3+{\mathrm{\mu }}_2{\mathrm{\mu }}_1+{\mathrm{\mu }}_2{\mathrm{\mu }}_3}{{\mathrm{\lambda }}_2{\mathrm{\lambda }}_1+{\mathrm{\lambda }}_3{\mathrm{\lambda }}_1+{\mathrm{\lambda }}_1{\mathrm{\mu }}_3+{\mathrm{\mu }}_1{\mathrm{\lambda }}_2+{\mathrm{\mu }}_3{\mathrm{\lambda }}_2+{\mathrm{\mu }}_1{\mathrm{\lambda }}_3+{\mathrm{\mu }}_2{\mathrm{\lambda }}_3+{\mathrm{\mu }}_2{\mathrm{\mu }}_1+{\mathrm{\mu }}_2{\mathrm{\mu }}_3}\end{array}\right]$

What about a purely numeric matrix?

$\mathrm{P2}≔\left[\begin{array}{ccc}\frac{1}{4}& \frac{1}{2}& \frac{1}{4}\\ \frac{1}{3}& \frac{1}{3}& \frac{1}{3}\\ \frac{1}{4}& \frac{1}{2}& \frac{1}{4}\end{array}\right]$

$\mathrm{\π}≔\left[\begin{array}{c}\frac{2}{7}\\ \frac{3}{7}\\ \frac{2}{7}\end{array}\right]$

$\mathrm{P3}≔\left[\begin{array}{ccc}0.25& 0.45& 0.3\\ 0.13& 0.33& 0.64\\ 0.2& 0.6& 0.2\end{array}\right]$