Complementaria 10: SDP en Python

Complementaria 10: SDP en Python

#

En este tutorial nos concentramos en la implementación del modúlo de SDP de jmarkov en Python, a través de dos ejemplos. Este módulo de jmarkov permite encontrar la política óptima de un proceso de decisión en el tiempo estocástico (SDP) mediante programación dinámica.

Problema 1: Máquinas

Se resolverá el primer problema del archivo Complementaria 11 (Q).pdf que se encuentra en Bloque Neón.

Vamos a modelar el problema de la máquina mediante un modelo de programación dinámica estocástica. El problema está definido para un horizonte de tiempo de 3 épocas (semanas), así, el conjunto de las épocas es: \( \text{Épocas} = \{1,2,3\} \)

Para poder determinar las posibles decisiones en una época, se necesita conocer el estado de la máquina al inicio de cada semana. Por ende, se define una única variable: \( X_{t} = \text{Estado de la máquina al inicio de la semana }t \) Esta variable puede tomar los siguientes estados: \( S_{X} = \{ \text{Excelente}, \text{Bueno}, \text{Promedio}, \text{Malo}\} \)

Todas las semanas se puede tomar la decisión de reemplazar o no reemplazar, así que las decisiones posibles no dependerán de la época. Sin embargo, no es posible tomar la decisión de reemplazar si la máquina se encuentra en el estado Excelente, por lo que sí dependerán del estado. \( A_{t}(i) = \{\text{Reemplazar},\text{No Reemplazar}\} \forall t \in E, i \in S_{X}|i\neq \text{Excelente} \) \( A_{t}(i) = \{\text{No Reemplazar}\} \forall t \in E, i = \text{Excelente} \)

Los retornos inmediatos dependen del estado que tenga la máquina al inicio de la semana y están dados por: \( r_{t}(i,a) = \text{Ganancia semanal}-\text{Costo asociado a la decisión} \) Dado que si se decide reemplazar la máquina, ésta queda en perfectas condiciones, los ingresos que recibirá esa semana corresponden a una máquina en excelente estado ($100). Los retornos se pueden entender como:

\( r_{t}(i,a)= \begin{array}{|c|c|c|} \hline & \text{Reemplazar}& \text{No Reemplazar} \\ \hline \text{Excelente} & -1000 & 100 \\ \text{Bueno} & -100 & 80 \\ \text{Promedio} & -100 & 50 \\ \text{Malo} & -100 & 10 \\ \hline \end{array} \)

El costo de reemplazar cuando la máquina está en excelentes condiciones toma un valor negativo de penalización, ya que ésta es una decisión infactible. De este modo, nos aseguramos que esta decisión nunca sea tomada.

Las probabilidades de transición están dadas por el enunciado, como se ve a continuación:

\( P_{i \to j}^{(a = No Reemplazar)}= \begin{array}{|c|cccc|} \hline & \text{E}& \text{B} & \text{P} & \text{M} \\ \hline \text{E} & 0.7 & 0.3 & 0 & 0 \\ \text{B} & 0 & 0.7 & 0.3 & 0 \\ \text{P} & 0 & 0 & 0.7 & 0.3 \\ \text{M} & 0 & 0 & 0 & 1 \\ \hline \end{array} \)

\( P_{i \to j}^{(a = Reemplazar)}= \begin{array}{|c|cccc|} \hline & \text{E}& \text{B} & \text{P} & \text{M} \\ \hline \text{E} & 0 & 0 & 0 & 0 \\ \text{B} & 0.7 & 0.3 & 0. & 0 \\ \text{P} & 0.7 & 0.3 & 0 & 0 \\ \text{M} & 0.7 & 0.3 & 0 & 0 \\ \hline \end{array} \)

Notemos que la matriz de Reemplazar cuenta con una condición particular: las entradas de primera fila son todas iguales a 0, en vez de sumar a 1; esto se debe a que la decisión de reemplazar no es fáctible para el estado excelente.

Para empezar con la implementación, llamamos las librerías necesarias.

import numpy as np
from jmarkov.sdp.dtsdp import dtsdp

Crearemos el espacio de las épocas en Python como un arreglo de numpy, de la siguiente manera:

E = np.array([i for i in range(1,4)])

Creamos el espacio de estados.

S = np.array(['Excelente','Bueno','Promedio','Malo']) # Estado de la máquina al inicio de la semana t

Creamos el espacio de acciones.

A = np.array(['Reemplazar','No Reemplazar'])

Creamos los retornos inmediatos.

R = np.zeros((len(E),len(S),len(A)))

# Recorremos sobre las épocas
for t in range(len(E)):
    # Recorremos sobre los estados:
    for s_index, i in enumerate(S):
        # Recorremos sobre las decisiones:
        for a_index, a in enumerate(A):
            if i=='Excelente' and a=='Reemplazar':
                R[t,s_index,a_index] = -1000
            elif i=='Excelente' and a=='No Reemplazar':
                R[t,s_index,a_index] = 100
            elif(i=="Bueno" and a=="Reemplazar"):
                R[t,s_index,a_index]=-100
            elif(i=="Bueno" and a=="No Reemplazar"):
                R[t,s_index,a_index]=80
            elif(i=="Promedio" and a=="Reemplazar"):
                R[t,s_index,a_index]=-100
            elif(i=="Promedio" and a=="No Reemplazar"):
                R[t,s_index,a_index]=50
            elif(i=="Malo" and a=="Reemplazar"):
                R[t,s_index,a_index]=-100
            elif(i=="Malo" and a=="No Reemplazar"):
                R[t,s_index,a_index]=10
            

Creamos las matrices de transición. Tendremos una matriz para cada época y para cada decisión.

matNoReemplazar = np.array([[0.7,0.3,0,0],
                          [0,0.7,0.3,0],
                          [0,0,0.7,0.3],
                          [0,0,0,1]])

matReemplazar = np.array([[1,0,0,0],
                          [0.7,0.3,0,0],
                          [0.7,0.3,0,0],
                          [0.7,0.3,0,0]])

probs = {}
for t in E:  # Iterar sobre cada época
    decisiones_dict = {}
    for posA, a in enumerate(A):
        if a == "Reemplazar":
            decisiones_dict[a] = matReemplazar
        elif a == "No Reemplazar":
            decisiones_dict[a] = matNoReemplazar
    probs[t] = decisiones_dict

probs
{1: {'Reemplazar': array([[1. , 0. , 0. , 0. ],
         [0.7, 0.3, 0. , 0. ],
         [0.7, 0.3, 0. , 0. ],
         [0.7, 0.3, 0. , 0. ]]),
  'No Reemplazar': array([[0.7, 0.3, 0. , 0. ],
         [0. , 0.7, 0.3, 0. ],
         [0. , 0. , 0.7, 0.3],
         [0. , 0. , 0. , 1. ]])},
 2: {'Reemplazar': array([[1. , 0. , 0. , 0. ],
         [0.7, 0.3, 0. , 0. ],
         [0.7, 0.3, 0. , 0. ],
         [0.7, 0.3, 0. , 0. ]]),
  'No Reemplazar': array([[0.7, 0.3, 0. , 0. ],
         [0. , 0.7, 0.3, 0. ],
         [0. , 0. , 0.7, 0.3],
         [0. , 0. , 0. , 1. ]])},
 3: {'Reemplazar': array([[1. , 0. , 0. , 0. ],
         [0.7, 0.3, 0. , 0. ],
         [0.7, 0.3, 0. , 0. ],
         [0.7, 0.3, 0. , 0. ]]),
  'No Reemplazar': array([[0.7, 0.3, 0. , 0. ],
         [0. , 0.7, 0.3, 0. ],
         [0. , 0. , 0.7, 0.3],
         [0. , 0. , 0. , 1. ]])}}

Finalmente, creamos el problema como un SDP.

sdpMaquinas = dtsdp(E,S, A, probs, R, 0.9)

Ahora, solucionamos el SDP. El método solve de la librería recibe el sentido del problema, es decir, si estamos minimizando o maximizando.

sdpMaquinas.solve(minimize=False)
(array([[255.151, 184.6  , 100.   ],
        [193.391, 143.9  ,  80.   ],
        [108.176,  84.2  ,  50.   ],
        [ 55.151,  19.   ,  10.   ]]),
 array([['N', 'N', 'N'],
        ['N', 'N', 'N'],
        ['N', 'N', 'N'],
        ['R', 'N', 'N']], dtype='<U1'))

Notamos que nos devuelve dos resultados. El primero se refiere al valor que toma cada una de las funciones de valor, para cada época y estado. Y el segundo, corresponde a la política óptima. En este caso, únicamente se debe reemplazar cuando se está en la primera época y la máquina se encuentra en un mal estado.

💡 Reto: Animate a implementar el segundo ejercicio propuesto!

Universidad de los Andes | Vigilada Mineducación. Reconocimiento como Universidad: Decreto 1297 del 30 de mayo de 1964. Reconocimiento personería jurídica: Resolución 28 del 23 de febrero de 1949 Minjusticia. Departamento de Ingeniería Industrial Carrera 1 Este No. 19 A 40 Bogotá, Colombia Tel. (57.1) 3324320 | (57.1) 3394949 Ext. 2880 /2881 http://industrial.uniandes.edu.co