Definiciones para el problema para:
un dia cualquiera: Parte operando A ó B:
falla: a
sigue operando: (1-a)
Parte en reparación A ó B:
se repara: b
sigue en reparación: (1-b)
Estados Xn:
0: No existen partes operativas/ todas fallaron
1: Una parte operativa/ una con falla
2: Dos partes operando/ no existen fallas
El diagrama a plantear con tres estados es:

Considera los eventos con las dos partes, estén operativas o dañadas. Los cambios, transición o pasos se consideran desde un dia para el siguiente día:
estado 0:
sigue en 0, en reparación A y en reparacion B:
(1-b)(1-b) = (1-b)2
pasa a 2, se repara A y se repara B: b*b = b2
pasa a 1, se repara A y B sigue en reparación, ó,
se repara B y A sigue en reparación:
b(1-b) + (1-b)b = 2b(1-b)
estado 1:
sigue en 1, no falla la operativa y
la otra sigue en reparación, ó,
falla la operativa y
se repara la que estaba con falla
(1-a)(1-b) + ab
pasa a 2, no falla la operativa y se repara la otra:
(1-a)b
pasa a 0, falla la operativa y
la otra sigue en reparación: a(1-b)
estado 2:
sigue en 2: no falla ninguna: (1-a)(1-a) = (1-a)2
pasa a 1: falla A y B sigue operando, ó,
falla B y A sigue operando:
a(1-a)+(1-a)a = 2a(1-a)
pasa a 0, falla A y falla B: a*a = a2
que al ponerlo en el digagrama, queda:

en la matriz de transición de estados de un paso P:
P=⎝⎛(1−b)2a(1−b)a22b(1−b)(1−a)(1−b)+ab2a(1−a)b2b(1−a)(1−a)2⎠⎞
que reemplazando los valores para a=0.1 y b=0.7 se verifica que las filas suman 1:
P=⎝⎛0.090.030.010.420.340.180.490.630.81⎠⎞
para la pmf de estado estable o largo plazo Pn:
0.09 Π0 + 0.03 Π1 + 0.01 Π2 = Π0
0.42 Π0 + 0.34 Π1 + 0.18 Π2 = Π1
0.49 Π0 + 0.63 Π1 + 0.81 Π2 = Π2
Π0 + Π1 + Π2 = 1
Resolviendo queda:
Π0 = 0.015625
Π1 = 0.21875
Π2 = 0.765625
La proyección para n partes indica que los exponentes de la matriz se convertirán en n, y otros términos sin exponentes deben aumentar para incluir las otras opciones.
Caso de resolver la matriz con python:
# 3ra Evaluación II Término 2017
# Tema 1
import numpy as np
# INGRESO
a = 0.1
b = 0.7
n = 100
# PROCEDIMIENTO
P = np.array([[(1-b)**2, 2*b*(1-b), b**2],
[a*(1-b), (1-a)*(1-b)+a*b, b*(1-a)],
[a**2, 2*a*(1-a), (1-a)**2]])
Pn = np.linalg.matrix_power(P,n)
# SALIDA
print(P)
print(Pn)
los resultados son
[[ 0.09 0.42 0.49]
[ 0.03 0.34 0.63]
[ 0.01 0.18 0.81]]
[[ 0.015625 0.21875 0.765625]
[ 0.015625 0.21875 0.765625]
[ 0.015625 0.21875 0.765625]]
>>>