How would you go about modelling this as a Markov chain? There is two machines which break down at different rates, µ$_A$ for machine A and µ$_B$ for machine B. When they break down, a machine can be fixed by one of two repairmen. Assume that two repairmen, X and Y, have different abilities, and they repair in exponential random times with parameters λ$_X$ and λ$_Y$ , respectively. (But the rate doesn’t depend on which machine they repair.)
How would you go about modelling this as a Markov chain?.
 A: Assume that if both machines are working when one fails, then the repairmen have equal probability of being chosen to fix the machine. The state space is
$$
S=\{(1,1), (0_X,1), (0_Y,1), (1,0_X), (1,0_Y), (0_X,0_Y), (0_Y,0_X) \}.
$$
The transition rates are given by
$$
q_{(i,j),(i',j')} = \begin{cases}
\frac{\mu_A}2,& (i,j)=(1,1)\text{ and } (i',j') \in \{(0_X,1),(0_Y,1)\}\\
\frac{\mu_B}2,& (i,j)=(1,1)\text{ and } (i',j') \in \{(1,0_X),(1,0_Y)\}\\
\lambda_X,& (i,j) \in \{(0_X,1),(1,0_X)\} \text{ and } (i',j') = (1,1)\\
\lambda_Y,& (i,j) \in \{(0_Y,1),(1,0_Y)\} \text{ and } (i',j') = (1,1)\\
\mu_A,& (i,j) = (1,0_X)\text{ and } (i',j') = (0_Y,0_X)\\
\mu_A,& (i,j) = (1,0_Y)\text{ and } (i',j') = (0_X,0_Y)\\
\mu_B,& (i,j) = (0_X,1)\text{ and } (i',j') = (0_X,0_Y)\\
\mu_B,& (i,j) = (0_Y,1)\text{ and } (i',j') = (0_Y,0_X)\\
\lambda_X,& (i,j) = (0_X,0_Y)\text{ and } (i',j') = (1,0_Y)\\
\lambda_X,& (i,j) = (0_Y,0_X)\text{ and } (i',j') = (0_Y,1)\\
\lambda_Y,& (i,j) = (0_X,0_Y)\text{ and } (i',j') = (0_X,1)\\
\lambda_Y,& (i,j) = (0_Y,0_X)\text{ and } (i',j') = (1,0_X)\\
0,& \text{otherwise.}
\end{cases}
$$
Let $Z(t)$ be the state of the system at time $t$, then $\{Z(t):t\geqslant 0\}$ is a continuous-time Markov chain with generator matrix
$$
Q = \small\left(
\begin{array}{ccccccc}
 -\left(\mu _A+\mu _B\right) & \frac{\mu _A}{2} & \frac{\mu _B}{2} & \frac{\mu _A}{2} & \frac{\mu _B}{2} & 0 & 0 \\
 \lambda _X & -\left(\mu _B+\lambda _X\right) & 0 & 0 & 0 & \mu _B & 0 \\
 \lambda _Y & 0 & -\left(\mu _B+\lambda _Y\right) & 0 & 0 & 0 & \mu _B \\
 \lambda _X & 0 & 0 & -\left(\mu _A+\lambda _X\right) & 0 & 0 & \mu _A \\
 \lambda _Y & 0 & 0 & 0 & -\left(\mu _B+\lambda _Y\right) & \mu _B & 0 \\
 0 & \lambda _Y & 0 & 0 & \lambda _X & -\left(\lambda _X+\lambda _Y\right) & 0 \\
 0 & 0 & \lambda _X & \lambda _Y & 0 & 0 & -\left(\lambda _X+\lambda _Y\right) \\
\end{array}
\right).
$$
The process has a unique stationary distribution $\pi$ which satisfies
$$
\pi_{(i,j)} = \lim_{t\to\infty} \mathbb P(Z_t = (i,j))
$$
(independent of the distribution of $Z_0$). We can find $\pi$ by computing the matrix exponential $e^{Qt}$ (which is the unique solution to the Kolmogorov backward equation $P'(t)=QP(t)$, $P'(0)=Q$) and taking any of the rows of $\lim_{t\to\infty} e^{Qt}$. More practically, $\pi$ satisfies the system of linear equations $\pi Q=0$. Note that $Q$ is singular (i.e. $\det Q=0$) as its rows all sum to zero, so we must replace one of the equations with $\sum_{(i,j)\in S} \pi_{(i,j)}=1$. However, due to the size of this matrix and the number of parameters, the closed form solution is a bit unwieldy. For example, I found that
$$
\pi_{(1,1)} = \tiny\frac{2 \lambda _X \lambda _Y \left(2 \mu _B+\lambda _X+\lambda _Y\right) \left(\mu _A+\mu _B+\lambda _X+\lambda _Y\right)}{\lambda _X^2 \left(\mu _B \left(3 \mu _A+10 \lambda _Y\right)+\left(\mu _A+2 \lambda _Y\right){}^2+6 \mu _B^2\right)+\lambda _X \left(\mu _B \left(7 \mu _A \mu _B+4 \mu _A^2+5 \mu _B^2\right)+\lambda _Y^2 \left(6 \mu _A+8 \mu _B\right)+\lambda _Y \left(\mu _A+3 \mu _B\right) \left(3 \mu _A+4 \mu _B\right)+2 \lambda _Y^3\right)+\left(\mu _B \left(3 \mu _A+4 \lambda _Y\right)+2 \lambda _Y \left(\mu _A+\lambda _Y\right)+\mu _B^2\right) \left(\mu _B \left(\mu _A+\mu _B\right)+\mu _A \lambda _Y\right)+2 \lambda _X^3 \left(\mu _B+\lambda _Y\right)}
$$
(the denominator is broken into two lines to prevent page stretching).
