DDP

#Definir los componentes del modelo

E = c(1,2,3)                 #Epocas (1:Bloqueador, 2:Agua, 3:Enlatados)
S = seq(0, 3000, by = 500)   #Espacio de estados
A = 0:6                      #Decisiones
W = c(500,1000,1500)         #Pesos
B = c(8,15,20)               #Beneficios


#Definir funcion de costo inmediato

costo = function(t, i, a){
  if (a*W[t] <= i) {
    resultado = a*B[t] #Si hay espacio disponible para meter 'a objetos (decisión factible)
  } else {
    resultado = -1e6 #Si no hay espacio disponible para meter 'a objetos (decisión infactible)
  }
  return(resultado)
}


#Matriz con los valores de las ecuaciones de Bellman
matF <- matrix(0, nrow = length(S), ncol = length(E))
rownames(matF) = S

#Matriz con las decisiones optimas
matDec = matrix('-', nrow = length(S), ncol = length(E))
rownames(matDec) = S

#----------------------------Recursion hacia atras. Ultima epoca--------------------------#


#Definir los valores de las ecuaciones de Bellman en la ultima epoca para cada estado

for (i in S) {
  valor = 0
  for (a in A) { #Iterar sobre las decisiones
    
    #Calcular el costo inmediato de tomar la decision 'a' en el estado 'i'
    valor = costo(length(E), i, a)
    
    #Debemos revisar si la decision 'a' es mejor a las otras decisiones
    if (valor > matF[as.character(i),length(E)]) {
      
      #Actualizar el valor de la ecuación de Bellman
      matF[as.character(i),length(E)] = valor
      
      #Actualizar la matriz de decisiones optimas
      matDec[as.character(i),length(E)] = a
      
    }
    
  }
}


#Recursion hacia atras. Desde la penultima epoca hasta la primera epoca.

for (t in (length(E)-1):1) { #Recorriendo las epocas
  for (i in S) {             #Recorriendo los estados
    
    valor = -1e6
    for (a in A) {           # Recorrer las decisiones
      
      #Solo considerar las decisiones factibles
      if (a*W[t] <= i) {
        
        #Calcular la suma del costo inmediato y el valor de la ecuación de Bellman en t+1
        valor = costo(t, i, a) + matF[as.character(i-a*W[t]),t+1]
        
        #Evaluar si la decision 'a' es mejor que las anteriores decisiones
        if (valor >= matF[as.character(i),t]) {
          
          #Actualizar el valor de la ecuación de Bellman
          matF[as.character(i),t] = valor
          
          #Actualizar la matriz de las decisiones optimas
          matDec[as.character(i),t] = a
        }
      }
    }
  }
}
print(matDec)
print(matF)

#Mapa calor
library(reshape2)
library(ggplot2)

# Melt the matrix
decisionesOptimasLista <- melt(matDec)

# Rename columns
colnames(decisionesOptimasLista) <- c("Estado", "Epoca", "Decision")


# Crear el mapa de calor
ggplot(data = decisionesOptimasLista, aes(x = Epoca, y = Estado, fill = Decision)) +
  geom_tile() +
  scale_fill_brewer(palette = "Reds")