From b8de690cb6b700e566c6caba55ad4a3fe4e2dd4c Mon Sep 17 00:00:00 2001 From: codeMGL Date: Fri, 8 May 2026 20:14:49 +0200 Subject: [PATCH 01/11] Matriz Recompensas v1 --- .gitignore | 4 + ControlModule.py | 223 +++++++++++++++++++++++++++++++++++++++++++---- main.py | 13 +-- 3 files changed, 219 insertions(+), 21 deletions(-) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..b2d0b1f --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +__pycache__/ControlModule.cpython-312.pyc +__pycache__/ControlModule.cpython-312.pyc +__pycache__/ControlModule.cpython-312.pyc +__pycache__/Reactor.cpython-312.pyc diff --git a/ControlModule.py b/ControlModule.py index ebe085c..2ef535a 100644 --- a/ControlModule.py +++ b/ControlModule.py @@ -8,37 +8,227 @@ def __init__(self): """Dummy constructor to use the Python Class as a namespace""" pass - @staticmethod # Hace que la función no necesite un argumento + @staticmethod def generate_P(probs) -> np.ndarray: """Function that generates the probabilities (transition) matrix""" ### TO BE COMPLETED BY THE STUDENTS ### # 3x3 o 3x3x3 o 100x100x3 - matriz_P = np.zeros((3, 10, 10), dtype=np.float64) + matrix_P = np.zeros((3, 10, 10), dtype=np.float64) # cambiar a 100x100 + probs_decrease = probs[0] probs_maintain = probs[1] - print(probs_maintain) + probs_increase = probs[2] + + # ---------------- DECREASE ---------------- + for i in range(10): + for j in range(10): + if i == j: + matrix_P[0][i][j] = probs_decrease[2] + elif i - 1 == j: + matrix_P[0][i][j] = probs_decrease[1] + elif i - 2 == j: + matrix_P[0][i][j] = probs_decrease[0] + # ---------------- MAINTAIN ---------------- for i in range(10): for j in range(10): if i == j: - matriz_P[1][i][j] = probs_maintain[1] + matrix_P[1][i][j] = probs_maintain[1] elif i + 1 == j: - matriz_P[1][i][j] = probs_maintain[2] - elif i-1 == j: - matriz_P[1][i][j] = probs_maintain[0] - + matrix_P[1][i][j] = probs_maintain[2] + elif i - 1 == j: + matrix_P[1][i][j] = probs_maintain[0] + # ---------------- INCREASE ---------------- + for i in range(10): + for j in range(10): + if i == j: + matrix_P[2][i][j] = probs_increase[0] + elif i + 1 == j: + matrix_P[2][i][j] = probs_increase[1] + elif i + 2 == j: + matrix_P[2][i][j] = probs_increase[2] - print("Probabilidades:\n", matriz_P.shape) - print(matriz_P) - return matriz_P - ... + print("Probabilidades:\n") + print(matrix_P) + return matrix_P @staticmethod - def generate_R() -> np.ndarray: + def generate_R(estado_actual, demand, n_states: np.int32 = 100) -> np.ndarray: """Function that generates the rewards (costs) matrix""" - ### TO BE COMPLETED BY THE STUDENTS ### - ... + demand = np.floor(demand[0] * 100) + # demand = 2 + """ + Coste_inicial = abs(estado_actual - demand) + matriz_calculo_costes = np.zeros((3, 3)) + + print(matriz_calculo_costes) # debería printearse una matriz 3x3 llena de ceros + + # Se mantiene constanteeeeeee + matriz_cambios_estado = np.array([[-2, -1, 0], [-1, 0, 1], [0, 1, 2]]) # 3 x 3 + matriz_P = np.zeros((3, 100, 100), dtype=np.float64) # cambiar a 100x100 + # Calculas la matriz que guarda los costes a partir de cada cambio de estado a partir de las posibles acciones + for estado_inicial in range(3): + for state in range(100): + for j in range(3): + coste_final_asociado = matriz_cambios_estado[estado_inicial][ + j + ] # 10 x 3 + + next_state = state + coste_final_asociado + next_state = np.clip(next_state, 0, 100 - 1) + matriz_calculo_costes[estado_inicial][j] = coste_final_asociado + matriz_P[estado_inicial][state][next_state] = 1 + print(matriz_P) + + # -------------------------------Territorio de Angel------------------- + # for i in range(3): + # for j in range(3): + # coste_accion_asociado = matriz_cambios_estado[i][j]# 10 x 3 + # next_state = state + coste_final_asociado + # matriz_calculo_costes[i][j] = coste_final_asociado + # print (matriz_P) + + print(matriz_calculo_costes) + """ + + ######################################################################## + matrix_R = np.zeros((3, 10, 10)) + # cambiar a 100x100 + # quitar DEBUGGING + # son recompensas, poner costes negativos + + # Se calcula la matriz de distancias entre el estado actual s y el estado futuro s' + # Aunque no se pueden alcanzar algunos estados (ej, pasar de s a s + 40), + # como en la matriz de probabilidades esa probabilidad es nula, el coste es indiferente + # Entonces, calculamos la distancia entre 'estado_inicial' y 'estado_final' y la duplicamos + # si la acción elegida se aleja del objetivo 'demand' + + # ---------------- DECREASE ---------------- + for estado_inicial in range(10): + for estado_final in range(10): + + delta_t = demand - estado_final + + # Comprobamos si la acción nos aleja del estado final (demanda) + # Si la distancia hasta la demanda es mayor en el estado actual vs. en el estado final, + # nos hemos alejado. Hay que multiplicar x2 la distancia + if delta_t > 0: + # print( + # "Nos alejamos (decrease)", + # estado_inicial, + # estado_final, + # "d_t", + # delta_t, + # ) + delta_t = -abs(delta_t * 2) + else: + delta_t = abs(delta_t) # Para DEBUGING + + # Rellenamos la matriz para estado_final x estado_inicial + matrix_R[0][estado_final][estado_inicial] = delta_t + + # ---------------- MAINTAIN ---------------- + for estado_inicial in range(10): + for estado_final in range(10): + + delta_t = demand - estado_final + + # Comprobamos si la acción nos aleja del estado final (demanda) + # Si la distancia hasta la demanda es mayor en el estado actual vs. en el estado final, + # nos hemos alejado. Hay que multiplicar x2 la distancia + if delta_t != 0: + # print( + # "Nos alejamos (maintain)", + # estado_inicial, + # estado_final, + # "d_t", + # delta_t, + # ) + delta_t = -abs(delta_t * 2) + else: + delta_t = abs(delta_t) # Para DEBUGING + + matrix_R[1][estado_final][estado_inicial] = delta_t + + # ---------------- INCREASE ---------------- + for estado_inicial in range(10): + for estado_final in range(10): + + delta_t = demand - estado_final + + # Comprobamos si la acción nos aleja del estado final (demanda) + # Si la distancia hasta la demanda es mayor en el estado actual vs. en el estado final, + # nos hemos alejado. Hay que multiplicar x2 la distancia + if delta_t < 0: + # print( + # "Nos alejamos (increase)", + # estado_inicial, + # estado_final, + # "d_t", + # delta_t, + # ) + delta_t = -abs(delta_t * 2) + else: + delta_t = abs(delta_t) # Para DEBUGING + + matrix_R[2][estado_final][estado_inicial] = delta_t + + print() + print("DEMANDA ACTUAL: ", demand) + print(matrix_R) + + # Creamos unos tests + print() + for estado_inicial in range(1, 3): + for estado_final in range(estado_inicial - 1, estado_inicial + 2): + coste0 = matrix_R[0][estado_final][estado_inicial] + coste1 = matrix_R[1][estado_final][estado_inicial] + coste2 = matrix_R[2][estado_final][estado_inicial] + print( + f"ACCIÓN: decrease. Para ir desde {estado_inicial} hasta {estado_final} --> Coste: {coste0}" + ) + print( + f"ACCIÓN: maintain. Para ir desde {estado_inicial} hasta {estado_final} --> Coste: {coste1}" + ) + print( + f"ACCIÓN: increase. Para ir desde {estado_inicial} hasta {estado_final} --> Coste: {coste2}" + ) + """ + [[[ *4. *4. *4. *4. *4. *4. *4. *4. *4. *4.] + [ *2. *2. *2. *2. *2. *2. *2. *2. *2. *2.] + [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] + [ 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] + [ 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.] + [ 3. 3. 3. 3. 3. 3. 3. 3. 3. 3.] + [ 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.] + [ 5. 5. 5. 5. 5. 5. 5. 5. 5. 5.] + [ 6. 6. 6. 6. 6. 6. 6. 6. 6. 6.] + [ 7. 7. 7. 7. 7. 7. 7. 7. 7. 7.]] + + [[ *4. *4. *4. *4. *4. *4. *4. *4. *4. *4.] + [ *2. *2. *2. *2. *2. *2. *2. *2. *2. *2.] + [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] + [ *2. *2. *2. *2. *2. *2. *2. *2. *2. *2.] + [ *4. *4. *4. *4. *4. *4. *4. *4. *4. *4.] + [ *6. *6. *6. *6. *6. *6. *6. *6. *6. *6.] + [ *8. *8. *8. *8. *8. *8. *8. *8. *8. *8.] + [*10. *10. *10. *10. *10. *10. *10. *10. *10. *10.] + [*12. *12. *12. *12. *12. *12. *12. *12. *12. *12.] + [*14. *14. *14. *14. *14. *14. *14. *14. *14. *14.]] + + [[ 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.] + [ 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] + [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] + [ *2. *2. *2. *2. *2. *2. *2. *2. *2. *2.] + [ *4. *4. *4. *4. *4. *4. *4. *4. *4. *4.] + [ *6. *6. *6. *6. *6. *6. *6. *6. *6. *6.] + [ *8. *8. *8. *8. *8. *8. *8. *8. *8. *8.] + [*10. *10. *10. *10. *10. *10. *10. *10. *10. *10.] + [*12. *12. *12. *12. *12. *12. *12. *12. *12. *12.] + [*14. *14. *14. *14. *14. *14. *14. *14. *14. *14.]]] + """ + return matrix_R @staticmethod def control_iteration() -> np.int32: @@ -60,3 +250,6 @@ def control_loop( ### DUMMY BEHAVIOUR TO PREVENT CRASHING (MUST BE DELETED AFTER THE FULL IMPLEMENTATION) ### return np.zeros_like(a=demand, dtype=np.float64) ### ### + + +#:) diff --git a/main.py b/main.py index 94ed703..46481dc 100644 --- a/main.py +++ b/main.py @@ -52,19 +52,20 @@ def main() -> None: probs = np.array([reactor.probabilities['decrease'], reactor.probabilities['maintain'], reactor.probabilities['increase']], dtype=np.float64) - + matriz_P = ControlModule.generate_P(probs) - + # Make a radar-plot with the reactor probabilities - plot_reactor_as_radar(probs=probs) - + # plot_reactor_as_radar(probs=probs) + # Generate a random power demand demand = generate_demand(n_samples=512) + matriz_R = ControlModule.generate_R(0, [0.02]) + # Define the number of MDP's states, actions and the discount factor (gamma) n_states = 100 n_actions = 3 - # Get the response time-series (answer to the demand time-series) response = ControlModule.control_loop(demand=demand, @@ -72,7 +73,7 @@ def main() -> None: n_states=n_states, n_actions=n_actions, gamma=gamma) - + # Plot the original power demand plot_demand(demand=demand) From d2e84fe3efa8a396d5d69b1301452954299118cc Mon Sep 17 00:00:00 2001 From: codeMGL Date: Fri, 8 May 2026 20:29:01 +0200 Subject: [PATCH 02/11] Cambio de variable para matrix_P --- .gitignore | 5 +-- ControlModule.py | 86 +++++++++++++++++++++++++++++++++--------------- 2 files changed, 61 insertions(+), 30 deletions(-) diff --git a/.gitignore b/.gitignore index b2d0b1f..c18dd8d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1 @@ -__pycache__/ControlModule.cpython-312.pyc -__pycache__/ControlModule.cpython-312.pyc -__pycache__/ControlModule.cpython-312.pyc -__pycache__/Reactor.cpython-312.pyc +__pycache__/ diff --git a/ControlModule.py b/ControlModule.py index 2ef535a..9979471 100644 --- a/ControlModule.py +++ b/ControlModule.py @@ -18,39 +18,73 @@ def generate_P(probs) -> np.ndarray: probs_decrease = probs[0] probs_maintain = probs[1] probs_increase = probs[2] - # ---------------- DECREASE ---------------- - for i in range(10): - for j in range(10): - if i == j: - matrix_P[0][i][j] = probs_decrease[2] - elif i - 1 == j: - matrix_P[0][i][j] = probs_decrease[1] - elif i - 2 == j: - matrix_P[0][i][j] = probs_decrease[0] + for estado_inicial in range(10): + for estado_final in range(10): + if estado_inicial == estado_final: + matrix_P[0][estado_inicial][estado_final] = probs_decrease[2] + elif estado_inicial - 1 == estado_final: + matrix_P[0][estado_inicial][estado_final] = probs_decrease[1] + elif estado_inicial - 2 == estado_final: + matrix_P[0][estado_inicial][estado_final] = probs_decrease[0] # ---------------- MAINTAIN ---------------- - for i in range(10): - for j in range(10): - if i == j: - matrix_P[1][i][j] = probs_maintain[1] - elif i + 1 == j: - matrix_P[1][i][j] = probs_maintain[2] - elif i - 1 == j: - matrix_P[1][i][j] = probs_maintain[0] + for estado_inicial in range(10): + for estado_final in range(10): + if estado_inicial == estado_final: + matrix_P[1][estado_inicial][estado_final] = probs_maintain[1] + elif estado_inicial + 1 == estado_final: + matrix_P[1][estado_inicial][estado_final] = probs_maintain[2] + elif estado_inicial - 1 == estado_final: + matrix_P[1][estado_inicial][estado_final] = probs_maintain[0] # ---------------- INCREASE ---------------- - for i in range(10): - for j in range(10): - if i == j: - matrix_P[2][i][j] = probs_increase[0] - elif i + 1 == j: - matrix_P[2][i][j] = probs_increase[1] - elif i + 2 == j: - matrix_P[2][i][j] = probs_increase[2] + for estado_inicial in range(10): + for estado_final in range(10): + if estado_inicial == estado_final: + matrix_P[2][estado_inicial][estado_final] = probs_increase[0] + elif estado_inicial + 1 == estado_final: + matrix_P[2][estado_inicial][estado_final] = probs_increase[1] + elif estado_inicial + 2 == estado_final: + matrix_P[2][estado_inicial][estado_final] = probs_increase[2] print("Probabilidades:\n") print(matrix_P) + """ + [[[0.8 0. 0. 0. 0. 0. 0. 0. 0. 0. ] + [0.025 0.8 0. 0. 0. 0. 0. 0. 0. 0. ] + [0.175 0.025 0.8 0. 0. 0. 0. 0. 0. 0. ] + [0. 0.175 0.025 0.8 0. 0. 0. 0. 0. 0. ] + [0. 0. 0.175 0.025 0.8 0. 0. 0. 0. 0. ] + [0. 0. 0. 0.175 0.025 0.8 0. 0. 0. 0. ] + [0. 0. 0. 0. 0.175 0.025 0.8 0. 0. 0. ] + [0. 0. 0. 0. 0. 0.175 0.025 0.8 0. 0. ] + [0. 0. 0. 0. 0. 0. 0.175 0.025 0.8 0. ] + [0. 0. 0. 0. 0. 0. 0. 0.175 0.025 0.8 ]] + + [[0.6 0.35 0. 0. 0. 0. 0. 0. 0. 0. ] + [0.05 0.6 0.35 0. 0. 0. 0. 0. 0. 0. ] + [0. 0.05 0.6 0.35 0. 0. 0. 0. 0. 0. ] + [0. 0. 0.05 0.6 0.35 0. 0. 0. 0. 0. ] + [0. 0. 0. 0.05 0.6 0.35 0. 0. 0. 0. ] + [0. 0. 0. 0. 0.05 0.6 0.35 0. 0. 0. ] + [0. 0. 0. 0. 0. 0.05 0.6 0.35 0. 0. ] + [0. 0. 0. 0. 0. 0. 0.05 0.6 0.35 0. ] + [0. 0. 0. 0. 0. 0. 0. 0.05 0.6 0.35 ] + [0. 0. 0. 0. 0. 0. 0. 0. 0.05 0.6 ]] + + [[0. 0.2 0.8 0. 0. 0. 0. 0. 0. 0. ] + [0. 0. 0.2 0.8 0. 0. 0. 0. 0. 0. ] + [0. 0. 0. 0.2 0.8 0. 0. 0. 0. 0. ] + [0. 0. 0. 0. 0.2 0.8 0. 0. 0. 0. ] + [0. 0. 0. 0. 0. 0.2 0.8 0. 0. 0. ] + [0. 0. 0. 0. 0. 0. 0.2 0.8 0. 0. ] + [0. 0. 0. 0. 0. 0. 0. 0.2 0.8 0. ] + [0. 0. 0. 0. 0. 0. 0. 0. 0.2 0.8 ] + [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.2 ] + [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]]] + """ + return matrix_P @staticmethod @@ -97,7 +131,7 @@ def generate_R(estado_actual, demand, n_states: np.int32 = 100) -> np.ndarray: # cambiar a 100x100 # quitar DEBUGGING # son recompensas, poner costes negativos - + # Se calcula la matriz de distancias entre el estado actual s y el estado futuro s' # Aunque no se pueden alcanzar algunos estados (ej, pasar de s a s + 40), # como en la matriz de probabilidades esa probabilidad es nula, el coste es indiferente From d22d5fe47b65e122dade5f59d746b7af64d7de27 Mon Sep 17 00:00:00 2001 From: codeMGL Date: Fri, 15 May 2026 00:17:26 +0200 Subject: [PATCH 03/11] Inicio control_interation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit TO DO - Comprobar que la matriz R tiene recompensas y no costes - Comprobar si la demanda es [0, 1] o [0, 100] - Hacer tests con diferentes demandas y matriz 10x10 - Revisar codigo base para ver args de funciones y orden de ejecución - Añadir método para subir/bajar de estado, una vez se ha elegido la acción --- .gitignore | 1 + ControlModule.py | 167 +++++++++++++++++++++++++---------------------- main.py | 87 +++++++++++++++--------- 3 files changed, 146 insertions(+), 109 deletions(-) diff --git a/.gitignore b/.gitignore index c18dd8d..5e43c01 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ __pycache__/ +/__pycache__ diff --git a/ControlModule.py b/ControlModule.py index 9979471..aa5a7a5 100644 --- a/ControlModule.py +++ b/ControlModule.py @@ -9,18 +9,19 @@ def __init__(self): pass @staticmethod - def generate_P(probs) -> np.ndarray: + def generate_P(probs, n_states) -> np.ndarray: """Function that generates the probabilities (transition) matrix""" ### TO BE COMPLETED BY THE STUDENTS ### - # 3x3 o 3x3x3 o 100x100x3 - matrix_P = np.zeros((3, 10, 10), dtype=np.float64) # cambiar a 100x100 + matrix_P = np.zeros( + (3, n_states, n_states), dtype=np.float64 + ) # cambiar a 100x100 probs_decrease = probs[0] probs_maintain = probs[1] probs_increase = probs[2] # ---------------- DECREASE ---------------- - for estado_inicial in range(10): - for estado_final in range(10): + for estado_inicial in range(n_states): + for estado_final in range(n_states): if estado_inicial == estado_final: matrix_P[0][estado_inicial][estado_final] = probs_decrease[2] elif estado_inicial - 1 == estado_final: @@ -29,8 +30,8 @@ def generate_P(probs) -> np.ndarray: matrix_P[0][estado_inicial][estado_final] = probs_decrease[0] # ---------------- MAINTAIN ---------------- - for estado_inicial in range(10): - for estado_final in range(10): + for estado_inicial in range(n_states): + for estado_final in range(n_states): if estado_inicial == estado_final: matrix_P[1][estado_inicial][estado_final] = probs_maintain[1] elif estado_inicial + 1 == estado_final: @@ -39,8 +40,8 @@ def generate_P(probs) -> np.ndarray: matrix_P[1][estado_inicial][estado_final] = probs_maintain[0] # ---------------- INCREASE ---------------- - for estado_inicial in range(10): - for estado_final in range(10): + for estado_inicial in range(n_states): + for estado_final in range(n_states): if estado_inicial == estado_final: matrix_P[2][estado_inicial][estado_final] = probs_increase[0] elif estado_inicial + 1 == estado_final: @@ -88,46 +89,12 @@ def generate_P(probs) -> np.ndarray: return matrix_P @staticmethod - def generate_R(estado_actual, demand, n_states: np.int32 = 100) -> np.ndarray: + def generate_R(demand_t: np.float64, n_states: np.int32 = 100) -> np.ndarray: """Function that generates the rewards (costs) matrix""" - demand = np.floor(demand[0] * 100) - # demand = 2 - """ - Coste_inicial = abs(estado_actual - demand) - matriz_calculo_costes = np.zeros((3, 3)) - - print(matriz_calculo_costes) # debería printearse una matriz 3x3 llena de ceros - - # Se mantiene constanteeeeeee - matriz_cambios_estado = np.array([[-2, -1, 0], [-1, 0, 1], [0, 1, 2]]) # 3 x 3 - matriz_P = np.zeros((3, 100, 100), dtype=np.float64) # cambiar a 100x100 - # Calculas la matriz que guarda los costes a partir de cada cambio de estado a partir de las posibles acciones - for estado_inicial in range(3): - for state in range(100): - for j in range(3): - coste_final_asociado = matriz_cambios_estado[estado_inicial][ - j - ] # 10 x 3 - - next_state = state + coste_final_asociado - next_state = np.clip(next_state, 0, 100 - 1) - matriz_calculo_costes[estado_inicial][j] = coste_final_asociado - matriz_P[estado_inicial][state][next_state] = 1 - print(matriz_P) - - # -------------------------------Territorio de Angel------------------- - # for i in range(3): - # for j in range(3): - # coste_accion_asociado = matriz_cambios_estado[i][j]# 10 x 3 - # next_state = state + coste_final_asociado - # matriz_calculo_costes[i][j] = coste_final_asociado - # print (matriz_P) - - print(matriz_calculo_costes) - """ + demand = demand_t + print("Demanda:", demand) - ######################################################################## - matrix_R = np.zeros((3, 10, 10)) + matrix_R = np.zeros((3, n_states, n_states)) # cambiar a 100x100 # quitar DEBUGGING # son recompensas, poner costes negativos @@ -139,8 +106,8 @@ def generate_R(estado_actual, demand, n_states: np.int32 = 100) -> np.ndarray: # si la acción elegida se aleja del objetivo 'demand' # ---------------- DECREASE ---------------- - for estado_inicial in range(10): - for estado_final in range(10): + for estado_inicial in range(n_states): + for estado_final in range(n_states): delta_t = demand - estado_final @@ -156,15 +123,15 @@ def generate_R(estado_actual, demand, n_states: np.int32 = 100) -> np.ndarray: # delta_t, # ) delta_t = -abs(delta_t * 2) - else: - delta_t = abs(delta_t) # Para DEBUGING + # else: + # delta_t = abs(delta_t) # Para DEBUGING # Rellenamos la matriz para estado_final x estado_inicial - matrix_R[0][estado_final][estado_inicial] = delta_t + matrix_R[0][estado_final][estado_inicial] = -delta_t # ---------------- MAINTAIN ---------------- - for estado_inicial in range(10): - for estado_final in range(10): + for estado_inicial in range(n_states): + for estado_final in range(n_states): delta_t = demand - estado_final @@ -180,14 +147,14 @@ def generate_R(estado_actual, demand, n_states: np.int32 = 100) -> np.ndarray: # delta_t, # ) delta_t = -abs(delta_t * 2) - else: - delta_t = abs(delta_t) # Para DEBUGING + # else: + # delta_t = abs(delta_t) # Para DEBUGING - matrix_R[1][estado_final][estado_inicial] = delta_t + matrix_R[1][estado_final][estado_inicial] = -delta_t # ---------------- INCREASE ---------------- - for estado_inicial in range(10): - for estado_final in range(10): + for estado_inicial in range(n_states): + for estado_final in range(n_states): delta_t = demand - estado_final @@ -203,31 +170,31 @@ def generate_R(estado_actual, demand, n_states: np.int32 = 100) -> np.ndarray: # delta_t, # ) delta_t = -abs(delta_t * 2) - else: - delta_t = abs(delta_t) # Para DEBUGING + # else: + # delta_t = abs(delta_t) # Para DEBUGING - matrix_R[2][estado_final][estado_inicial] = delta_t + matrix_R[2][estado_final][estado_inicial] = -delta_t print() print("DEMANDA ACTUAL: ", demand) print(matrix_R) # Creamos unos tests - print() - for estado_inicial in range(1, 3): - for estado_final in range(estado_inicial - 1, estado_inicial + 2): - coste0 = matrix_R[0][estado_final][estado_inicial] - coste1 = matrix_R[1][estado_final][estado_inicial] - coste2 = matrix_R[2][estado_final][estado_inicial] - print( - f"ACCIÓN: decrease. Para ir desde {estado_inicial} hasta {estado_final} --> Coste: {coste0}" - ) - print( - f"ACCIÓN: maintain. Para ir desde {estado_inicial} hasta {estado_final} --> Coste: {coste1}" - ) - print( - f"ACCIÓN: increase. Para ir desde {estado_inicial} hasta {estado_final} --> Coste: {coste2}" - ) + # print() + # for estado_inicial in range(1, 3): + # for estado_final in range(estado_inicial - 1, estado_inicial + 2): + # coste0 = matrix_R[0][estado_final][estado_inicial] + # coste1 = matrix_R[1][estado_final][estado_inicial] + # coste2 = matrix_R[2][estado_final][estado_inicial] + # print( + # f"ACCIÓN: decrease. Para ir desde {estado_inicial} hasta {estado_final} --> Coste: {coste0}" + # ) + # print( + # f"ACCIÓN: maintain. Para ir desde {estado_inicial} hasta {estado_final} --> Coste: {coste1}" + # ) + # print( + # f"ACCIÓN: increase. Para ir desde {estado_inicial} hasta {estado_final} --> Coste: {coste2}" + # ) """ [[[ *4. *4. *4. *4. *4. *4. *4. *4. *4. *4.] [ *2. *2. *2. *2. *2. *2. *2. *2. *2. *2.] @@ -265,9 +232,32 @@ def generate_R(estado_actual, demand, n_states: np.int32 = 100) -> np.ndarray: return matrix_R @staticmethod - def control_iteration() -> np.int32: + def control_iteration(P, R, estado_actual, gamma, max_iter=1000) -> np.int32: """Function that computes one control-iteration""" - ### TO BE COMPLETED BY THE STUDENTS ### + # print("P y R", P.shape, R.shape) + # print() + # print(P[0][:3][:3]) + # print("P\n", P) + # print(R[0][:3][:3]) + # print() + + # -- Comprobacion provisional de si la matriz P es estocastica + # La última fila tiene todo ceros, rellenamos con uno en la posicion final + P[2][-1][-1] = 1 + + for a in range(3): + P[a] = P[a] / P[a].sum(axis=1, keepdims=True) # renormaliza + + # print() + # print("P normalizada\n", P) + # print(check_stochastic(P)) + # print() + + pi = mdptoolbox.mdp.PolicyIteration(P, R, gamma, max_iter=max_iter) + pi.setVerbose() + pi.run() + print("Policy: ", pi.policy) + return pi.policy[estado_actual] ... @staticmethod @@ -287,3 +277,22 @@ def control_loop( #:) + + +def check_stochastic(P, tol=1e-6): + P = np.array(P) + assert P.ndim == 3, f"P debe ser (A, S, S), tiene forma {P.shape}" + A, S, _ = P.shape + ok = True + for a in range(A): + row_sums = P[a].sum(axis=1) + bad = np.where(np.abs(row_sums - 1.0) > tol)[0] + if len(bad) > 0: + ok = False + for s in bad: + print(f" Acción {a}, Estado {s}: suma = {row_sums[s]:.8f}") + if ok: + print("✅ Matriz estocástica: todas las filas suman 1.") + else: + print("❌ Matriz NO estocástica.") + return ok diff --git a/main.py b/main.py index 46481dc..75870c6 100644 --- a/main.py +++ b/main.py @@ -8,39 +8,49 @@ from Metrics import * from Plotter import * -def get_args() -> tuple[Reactor, np.float64, np.float64]: # 3 salidas (?) + +def get_args() -> tuple[Reactor, np.float64, np.float64]: # 3 salidas (?) # Define the parser object parser = argparse.ArgumentParser() # Define the expected arguments to parse and their data types - parser.add_argument("--input-reactor", "-i", type=str, help="Path of the reactor's JSON file") - parser.add_argument("--gamma", "-g", type=float, help="Discount factor used in the MDP") - parser.add_argument("--random-seed", "-r", type=int, help="Pseudo-random number generator seed") + parser.add_argument( + "--input-reactor", "-i", type=str, help="Path of the reactor's JSON file" + ) + parser.add_argument( + "--gamma", "-g", type=float, help="Discount factor used in the MDP" + ) + parser.add_argument( + "--random-seed", "-r", type=int, help="Pseudo-random number generator seed" + ) # Parse the arguments args = parser.parse_args() - + # Some verbose to check the correct parsing of the input arguments print(f"Loading reactor from file: {args.input_reactor}") print(f"Using gamma (discount factor): {args.gamma}") print(f"Using {args.random_seed} as random seed") # Build the Reactor object by reading the reactor's JSON file - with open(args.input_reactor, 'r', encoding='utf-8') as file: + with open(args.input_reactor, "r", encoding="utf-8") as file: json_data = json.load(fp=file) - reactor = Reactor(model=json_data['model'], - effective_section=float(json_data['effective_section']), - neutron_flux=float(json_data['neutron_flux']), - core_volume=float(json_data['core_volume']), - fision_energy=float(json_data['fision_energy']), - probabilities=dict(json_data['probabilities'])) - + reactor = Reactor( + model=json_data["model"], + effective_section=float(json_data["effective_section"]), + neutron_flux=float(json_data["neutron_flux"]), + core_volume=float(json_data["core_volume"]), + fision_energy=float(json_data["fision_energy"]), + probabilities=dict(json_data["probabilities"]), + ) + # Some verbose of the reactor loaded print(reactor) # Overloaded in the __str__ method of Reactor's class - + # Return the Reactor object, gamma and the random seed return reactor, args.gamma, args.random_seed + def main() -> None: # Parse the main arguments reactor, gamma, random_seed = get_args() @@ -49,30 +59,46 @@ def main() -> None: np.random.seed(random_seed) # Get the probabilities from the reactor's dynamics - probs = np.array([reactor.probabilities['decrease'], - reactor.probabilities['maintain'], - reactor.probabilities['increase']], dtype=np.float64) + probs = np.array( + [ + reactor.probabilities["decrease"], + reactor.probabilities["maintain"], + reactor.probabilities["increase"], + ], + dtype=np.float64, + ) - matriz_P = ControlModule.generate_P(probs) + matriz_P = ControlModule.generate_P(probs, 100) # Make a radar-plot with the reactor probabilities # plot_reactor_as_radar(probs=probs) # Generate a random power demand - demand = generate_demand(n_samples=512) + # demand = generate_demand(n_samples=512) + # print("Demand:\n", demand) # Array con 512 demands + + estado_actual = 5 + # demand = 0.02 + # for i in range(0, 100, 20): + i = 90 + if True: + demand = i / 100 + matriz_R = ControlModule.generate_R(demand) - matriz_R = ControlModule.generate_R(0, [0.02]) + # Ejecutamos para un solo estado + response = ControlModule.control_iteration(matriz_P, matriz_R, estado_actual, gamma) + print("Response 'unitaria':", response) # Define the number of MDP's states, actions and the discount factor (gamma) - n_states = 100 + n_states = 100 n_actions = 3 # Get the response time-series (answer to the demand time-series) - response = ControlModule.control_loop(demand=demand, - probs=probs, - n_states=n_states, - n_actions=n_actions, - gamma=gamma) + response = ControlModule.control_loop( + demand=demand, probs=probs, n_states=n_states, n_actions=n_actions, gamma=gamma + ) + + # print("Response:\n", response) # Array con 512 respuestas # Plot the original power demand plot_demand(demand=demand) @@ -87,9 +113,9 @@ def main() -> None: plot_correlation(demand=demand, response=response) # Print the four regression metrics for the current demand-response data - _MAE = MAE(y_true=demand, y_pred=response) - _MSE = MSE(y_true=demand, y_pred=response) - _R2 = R2(y_true=demand, y_pred=response) + _MAE = MAE(y_true=demand, y_pred=response) + _MSE = MSE(y_true=demand, y_pred=response) + _R2 = R2(y_true=demand, y_pred=response) _Corr = Corr(y_true=demand, y_pred=response) print(f"MAE={_MAE:.6f}") print(f"MSE={_MSE:.6f}") @@ -102,5 +128,6 @@ def main() -> None: # Plot the R2 and the Corr in a bar-plot plot_r2_and_pearson(R2=_R2, Pearson=_Corr) -if __name__ == '__main__': + +if __name__ == "__main__": main() From 93ae8d2085eea3849db3c9b57e1123151ab3f613 Mon Sep 17 00:00:00 2001 From: Angel Date: Fri, 15 May 2026 16:05:56 +0200 Subject: [PATCH 04/11] Loop_control, Matrix_P cambiada Ahora la matriz_P no se genera bien --- ControlModule.py | 153 +++++++++++++++------- __pycache__/ControlModule.cpython-313.pyc | Bin 2086 -> 8842 bytes main.py | 23 ++-- 3 files changed, 119 insertions(+), 57 deletions(-) diff --git a/ControlModule.py b/ControlModule.py index aa5a7a5..41812c4 100644 --- a/ControlModule.py +++ b/ControlModule.py @@ -9,17 +9,35 @@ def __init__(self): pass @staticmethod - def generate_P(probs, n_states) -> np.ndarray: + def generate_P(probs: np.ndarray, n_states: np.int32 = 100) -> np.ndarray: """Function that generates the probabilities (transition) matrix""" - ### TO BE COMPLETED BY THE STUDENTS ### - matrix_P = np.zeros( - (3, n_states, n_states), dtype=np.float64 - ) # cambiar a 100x100 + matrix_P = np.zeros((3, n_states, n_states), dtype=np.float64) # cambiar a 100x100 + #saca las probabilidades del Json probs_decrease = probs[0] probs_maintain = probs[1] probs_increase = probs[2] + + #Las probabilidades en los bordes se pierden, así que las filas no suman 1. (Y creo que normalizar es un apaño que no debería estar bien) + #Ej: COn decrease en estado 1, sumas las probabilidad de bajar a 0, pero no a la de bajar a -1 porque nunca se cumple en los IFs + + #Tantas iteraciones creo que son ineficientes, y mi pc es una patata, así que creo que hay que aprovechar mejor el hecho que + #sabemos que solo hay que modificar 3 IJs, no hace falta comprobar los 100, el 97% de las comprobaciones son inncesarias + # ---------------- DECREASE ---------------- + """ + + DECREASE ANGEL + for estado_inicial in range(n_states): + dest_0 = max(0, estado_inicial - 2) + dest_1 = max(0, estado_inicial - 1) + dest_2 = estado_inicial + + matrix_P[0][estado_inicial][dest_0] += probs_decrease[0] + matrix_P[0][estado_inicial][dest_1] += probs_decrease[1] + matrix_P[0][estado_inicial][dest_2] += probs_decrease[2] + + """ for estado_inicial in range(n_states): for estado_final in range(n_states): if estado_inicial == estado_final: @@ -49,6 +67,33 @@ def generate_P(probs, n_states) -> np.ndarray: elif estado_inicial + 2 == estado_final: matrix_P[2][estado_inicial][estado_final] = probs_increase[2] + # ---------------- CORRECCIÓN DE BORDES ---------------- + # Si una acción intenta llevar el reactor por debajo del estado 0 + # o por encima del estado n_states - 1, esa probabilidad no se pierde: + # se acumula en el borde correspondiente. + + # DECREASE: + # En estado 0, los movimientos -2 y -1 se quedan en 0. + matrix_P[0][0][0] += probs_decrease[0] + probs_decrease[1] + + # En estado 1, el movimiento -2 se queda en 0. + matrix_P[0][1][0] += probs_decrease[0] + + # MAINTAIN: + # En estado 0, el movimiento -1 se queda en 0. + matrix_P[1][0][0] += probs_maintain[0] + + # En el último estado, el movimiento +1 se queda en el último estado. + matrix_P[1][n_states - 1][n_states - 1] += probs_maintain[2] + + # INCREASE: + # En el penúltimo estado, el movimiento +2 se queda en el último estado. + matrix_P[2][n_states - 2][n_states - 1] += probs_increase[2] + + # En el último estado, los movimientos +1 y +2 se quedan en el último estado. + matrix_P[2][n_states - 1][n_states - 1] += probs_increase[1] + probs_increase[2] + + print("Probabilidades:\n") print(matrix_P) """ @@ -91,11 +136,10 @@ def generate_P(probs, n_states) -> np.ndarray: @staticmethod def generate_R(demand_t: np.float64, n_states: np.int32 = 100) -> np.ndarray: """Function that generates the rewards (costs) matrix""" - demand = demand_t + demand = np.float64(demand_t) print("Demanda:", demand) - matrix_R = np.zeros((3, n_states, n_states)) - # cambiar a 100x100 + matrix_R = np.zeros((3, n_states, n_states), dtype=np.float64) # (3x)100x100 # quitar DEBUGGING # son recompensas, poner costes negativos @@ -105,16 +149,19 @@ def generate_R(demand_t: np.float64, n_states: np.int32 = 100) -> np.ndarray: # Entonces, calculamos la distancia entre 'estado_inicial' y 'estado_final' y la duplicamos # si la acción elegida se aleja del objetivo 'demand' + + #al final hacer otra iteración exterior, y probar con np.where # ---------------- DECREASE ---------------- for estado_inicial in range(n_states): for estado_final in range(n_states): - delta_t = demand - estado_final + delta_t = demand - (estado_final/ n_states) + coste = abs(delta_t) # Comprobamos si la acción nos aleja del estado final (demanda) # Si la distancia hasta la demanda es mayor en el estado actual vs. en el estado final, # nos hemos alejado. Hay que multiplicar x2 la distancia - if delta_t > 0: + if demand > (estado_inicial/ n_states): # print( # "Nos alejamos (decrease)", # estado_inicial, @@ -122,46 +169,34 @@ def generate_R(demand_t: np.float64, n_states: np.int32 = 100) -> np.ndarray: # "d_t", # delta_t, # ) - delta_t = -abs(delta_t * 2) + coste *= 2 # else: # delta_t = abs(delta_t) # Para DEBUGING - # Rellenamos la matriz para estado_final x estado_inicial - matrix_R[0][estado_final][estado_inicial] = -delta_t + matrix_R[0][estado_inicial][estado_final] = -coste # ---------------- MAINTAIN ---------------- for estado_inicial in range(n_states): for estado_final in range(n_states): - delta_t = demand - estado_final + delta_t = demand - (estado_final/ n_states) - # Comprobamos si la acción nos aleja del estado final (demanda) - # Si la distancia hasta la demanda es mayor en el estado actual vs. en el estado final, - # nos hemos alejado. Hay que multiplicar x2 la distancia - if delta_t != 0: - # print( - # "Nos alejamos (maintain)", - # estado_inicial, - # estado_final, - # "d_t", - # delta_t, - # ) - delta_t = -abs(delta_t * 2) - # else: - # delta_t = abs(delta_t) # Para DEBUGING + coste = abs(delta_t) - matrix_R[1][estado_final][estado_inicial] = -delta_t + + matrix_R[1][estado_inicial][estado_final] = -coste # ---------------- INCREASE ---------------- for estado_inicial in range(n_states): for estado_final in range(n_states): - delta_t = demand - estado_final + delta_t = demand - (estado_final/ n_states) + coste = abs(delta_t) # Comprobamos si la acción nos aleja del estado final (demanda) # Si la distancia hasta la demanda es mayor en el estado actual vs. en el estado final, # nos hemos alejado. Hay que multiplicar x2 la distancia - if delta_t < 0: + if demand < (estado_inicial/ n_states): # print( # "Nos alejamos (increase)", # estado_inicial, @@ -169,11 +204,11 @@ def generate_R(demand_t: np.float64, n_states: np.int32 = 100) -> np.ndarray: # "d_t", # delta_t, # ) - delta_t = -abs(delta_t * 2) + coste *= 2 # else: # delta_t = abs(delta_t) # Para DEBUGING - matrix_R[2][estado_final][estado_inicial] = -delta_t + matrix_R[2][estado_inicial][estado_final] = -coste print() print("DEMANDA ACTUAL: ", demand) @@ -234,6 +269,10 @@ def generate_R(demand_t: np.float64, n_states: np.int32 = 100) -> np.ndarray: @staticmethod def control_iteration(P, R, estado_actual, gamma, max_iter=1000) -> np.int32: """Function that computes one control-iteration""" + + # ECUACION BELLMAN: V(s) = max_a [ sum_s'( P[s,s',a] * (R[s,s',a] + gamma * V[s']) ) ] + + # print("P y R", P.shape, R.shape) # print() # print(P[0][:3][:3]) @@ -246,40 +285,57 @@ def control_iteration(P, R, estado_actual, gamma, max_iter=1000) -> np.int32: P[2][-1][-1] = 1 for a in range(3): - P[a] = P[a] / P[a].sum(axis=1, keepdims=True) # renormaliza + P[a] = P[a] / P[a].sum(axis=1, keepdims=True) # renormaliza # Cambiar si hay tiempo, en vez de distribuir las probabilidades uniformemente, truncarlas a los bordes # print() # print("P normalizada\n", P) # print(check_stochastic(P)) # print() + #politica optima con la libreria .... pi = mdptoolbox.mdp.PolicyIteration(P, R, gamma, max_iter=max_iter) pi.setVerbose() pi.run() print("Policy: ", pi.policy) return pi.policy[estado_actual] - ... + + @staticmethod - def control_loop( - demand: np.ndarray, - probs: np.ndarray, - n_states: np.int32, - n_actions: np.int32, - gamma: np.float64, - ) -> np.ndarray: + def control_loop(demand: np.ndarray, + probs: np.ndarray, + n_states: np.int32, + n_actions: np.int32, + gamma: np.float64,) -> np.ndarray: + """Function that computes all the required iterations (control-loop) to satisfy the power demand""" - ### TO BE COMPLETED BY THE STUDENTS ### + ControlModule._probs = probs + ControlModule._n_states = n_states + ControlModule._n_actions = n_actions + + respuesta = np.zeros_like(a=demand, dtype=np.float64) # Almacena las acciones para cada demanda + current_state = np.int32(0) # Estado inicial, se puede modificar si se desea empezar en otro estado + + action_deltas = [np.array([-2, -1, 0], dtype=np.int32), np.array([-1, 0, 1], dtype=np.int32),np.array([ 0, 1, 2], dtype=np.int32)] + ControlModule._P = ControlModule.generate_P(probs, n_states) + + for t in range(demand.shape[0]): # La demanda cambia en el tiempo - ### DUMMY BEHAVIOUR TO PREVENT CRASHING (MUST BE DELETED AFTER THE FULL IMPLEMENTATION) ### - return np.zeros_like(a=demand, dtype=np.float64) - ### ### + ControlModule._demand_t = np.float64(demand[t]) + ControlModule._current_state = current_state + ControlModule._R = ControlModule.generate_R(demand_t = ControlModule._demand_t, n_states) # esta demanda concreta. + action = ControlModule.control_iteration(P= ControlModule._P, R= ControlModule._R, current_state = ControlModule._current_state, gamma=gamma) -#:) + state_increment = np.random.choice( a=action_deltas[action], p=probs[action]) + current_state = current_state + state_increment + current_state = np.int32(np.clip(a=current_state, a_min=0, a_max=n_states - 1)) + respuesta[t] = current_state / 100.0 + return respuesta def check_stochastic(P, tol=1e-6): + P = np.array(P) assert P.ndim == 3, f"P debe ser (A, S, S), tiene forma {P.shape}" A, S, _ = P.shape @@ -296,3 +352,6 @@ def check_stochastic(P, tol=1e-6): else: print("❌ Matriz NO estocástica.") return ok + + + diff --git a/__pycache__/ControlModule.cpython-313.pyc b/__pycache__/ControlModule.cpython-313.pyc index 87888329d161bbae779494d33dad907291b3e813..43629b4449631ce9d92e6dd72d39eb26fd96968e 100644 GIT binary patch literal 8842 zcmeG>ZEzdMb$5ra0||g%A_bBXd88zZkVx7TDVu~oB1kkvMItTXj0m=9p-13Eg2o4V z2iPL(H0n&+0WGDXJaI%iX~}px4gII3rs-Io$xM`{Q`?=+6e!h@bKwt9qv=$gc7_Vs z`d`xb_HYLf0wt@HpWOlW_IBTU`}Xa=kA1gcx0g^*ep_Z4EB8>;ukb}H#%yEb5H#MW zc#5abQW8y`jKnypW{PL(DBf_GQ6VJ*Z{m%0)EF%pM`)M%A$AW2)j_qkXJ0n{l17?C zgr3CkG~fa)$v9-<%X!0)+3UzhQ`@`|`pb*-n~KACnRzpCV;rLoXOSc;p>F~B9YyFo zJAJl*ekste;E|sqe0DLswE*4*@Fm5$e`YC^7N3>3ZL`|#q-Nvo$fYhbJMWMjBW12i z#eOOtOGxp^Ks=O;h!4#m!XWIUUa)O|{{|kN_bE67%v5)8G8&!ag7KJ~kdnbfT;dXO zE-8y#VpQbLPbNm=G451EkY!GQ9~Tp%qC74H#lep#7p)j&F>*yl!kp3@bPxGtQIdx| zZ*>fG4E4r?$*34h$U~>YiT>p9P&72&CW!)0(iTewBcdP;ZDF~6d{SWpfp9FG2n3Xx ztz%lh1tybCP+gxQY@G~`Upg!Jjn7xu>SX{ zOQaJkZ$4DE7nFzgS9BTjpfrWk@hDKR9wx^4K*^Aym;qWw(lSEJL|Uc>>JqX{ijO{= z-HfEY~x*Q(I zP*;^F2YppWeN~TR{q9#~(q~im?CGllvz*=C$1$|M(&eD9%B0KUajb$>X{SKtj5i!( zjxl6LZNCjF2W>D*8#n>kATxNjbUeb}<9G$IStz8O2b-0CF<2k1o?F68wypD)-RlB8 zb*(f{PH$1XtX#IiBj;lqJarc|`~q2AvO)(byjjy?j>`4=m|SlINNz9XvP>Q8OU8nU zFxb$KD~Gw9TkU#oq!b`!1=ez=)^z@~R~{^Zb4-JiHKRedRT z-AvUr&Q8oZQ;yZzhJwa=rOj&FtfLzDF7*FdJTOBiWrh5P4(dFuwo(}Xc+XHRbsiqem8K}ng%^m;G zg(i!2j}WQIKn-LSdE0W)kn4zo`l z*0K9%CvKg~uuTgsix+;_y39Vmbo^dZrn5K0_QFe!%v0c=O7Gquq`z@7bK0Nj^k-N< zklq9oN7~VJpKV@a?WxflW$R3n)wbS8Q{~Q7nPek+7ZRgFEEGs=iBHfdk@#d5RN(3H ziI39kO^msOxTNqj*da6zP{2B=4XC{s{S*|uAQZ$nSb*e7e9NzC&? zHWoRVjYV|Zyd}cwd;)h(UzNKo52_c>=V4El2Z+QmF4{&4C!WB=vr7o1QdB}JU-{e> zN=QE;+Qk7l*86zNmQ$nm`wZF{<0$ZH2JmU@)S)N&G)a70kV4?if^j(^YmRJ+?IkB( za7z`i{SqvLlB%)711I5SZ>qfa!^7cWGX!5yv%dukedcB8KYoh7+c(3`R@^F`9n6#;Ts*g2-m~G-%;|Mh zXt-^Wcvg#V(MDW?F}3hU?G(eUYZuXb5=Ht+?L|rjg@^c(H}^1Tx%?(#aT>B(uHOQS zSyZoXhjIk{w7-kD^cZ7ipV_P1zub6EA7NZez~qLh-1wz3t(#{7wnz#iXZeed`-qf4 zi>*_wTeUbCkB%oXFpbAV%m!c{p)H)iXcYodqM<}#AO02!QQ}u!jA9h7hT+=%rYMew z!ciF#0uWMO4IthWdp;fs2PZ*~ty!IwlToDvk^qT#JTe@=su-ZHl&N56w5gO5SxkIG zl!oIFh8m<~jHte19w*dX7KJ{q&|bx!4LBj2kc9ZtG$KT!f&@m7QfNV8#=|lmk}WZ4 zMRr+)J+dsNace007=Q!Lko+zXx<)|=YZ&#*Y#wodQzV2y$}pm)@@c&4(H9w)0bE52Q&7Ai=L%@Y5T#n z{cMU^WlL__r|mb&uwkAy-(Xd1(md~)bKPh6u9j8YoSdGV-JL1hpJw-y2ts-e7ix_K zVEM2;vUCLAx2OmnM4}@90rmTokIsbx5nVBSOkS_<`h9sxkBF`umoTU&#S5;yuryOX zDh~pqoPxsB5Js5tV2WfAF(A=Y6@rD_j75%mGxa=NGB411P(>*6#vUdQVmmsf-w_@2 zSE2K~A!f4D0uZ(YKokLi|A6ngl>=&ALT{oSWr|`sL2sb-a42(HctvgcBh&0 zvuI&8th{-3ljJk!_3E&9iXEOZ;PP(B5)>b!ovswlwGK%K@!7 z_XLeH^9&?ltQb1;B|Um5hK__UC9cU&T+^>r*W@#SVYf}av5jFDA`y~xki>~(SQ0~= zZrw?=rds?+JU;Hi46qCaQob^&<|*QDixQ^>ST06_JRv3FsU^Xi%kc0w=n`In2e zdF_O&MnadYm;-7K*c#B2z$MVs3?xR9N(sq~1R~)#MG{piCP;5|9s-h*Buq+82*U)< ztM1Aqj|$_WVuhrLdRbE(fh}xGxN89Km5yVhEa2VZO&|vicf5F1F$YKE;h?A(gOTvK ztCZO8ToOh)0m#IDD^>|o>q&HM(kY~9SD6O}Abgb-I-!`=9;J+=xgq-k=@^n?$G9Tv z?C_x`U2K{^!(CPCws6;;Kh}Y6?%Pm-n=7ZdeK+j)*PpF026 z`46r4FaO@k|B|Zk3NBSd7zrujgsy)x7>}&Nc zEA^e3`p%`7m+O0x6DRuJw{E|4>zxJ9t#AFJVt>l>;JKC`jeR(_RCA|x`MD!E2JYMT ztRC#fXUl`e)|JMWGL0|YH7_@w!LD6v)y*r_U76~xpV(Kr`!n7B|MI5@{T;U^K-31)gyuqaU9yD9sIDX3_p(G2N?J`QW6M6$!8vD z+X8_J2!gUBWr4tzup}oU;g}eU2Lcj04#kEiCLD~4@Ifa;?p2B<7E;gsSpXo1U%hJK zu_i}aLW`+X2Z&aNg8?A{Ny%Z15`ii*sIrQo2GK!MVV*($4OD+j{hDD+cDS{eD<3hx zPn(>x&GQH54y;qyUSd9G@33n7-l0$3KXY%A{tW|wuM_a~DQrJtPBMV%oqv7q^>qr{ zONT!0zSFI??+t%C_Or1~(w{{|koA5eL2WR#e8&6}`#HNwK{E>gExcObG`mf;8f9p% z*6*Ho&$-tP(Eb88d30`c-2`n&N^#BeSLd#-TkwsgU~3W1s=l?h|CReA_XfU^88(ot zDjKv*TXH<8rx;Rg{3mGJGv#<7!15xwg2L1I#U}(Nb|(Z+0E6-`C=RNV(QXkU+?9&K zB!R2(Ju%4DLy)H5Vsbfc2qUJ@srgABq}V-%hkShxy(P3M8ffIX;i7I4vOi{Yd0Hsq zqut-mz3~=wv=Gi<9)T1et3mP>%pj5PBOnu=jhTEb@C~U0aO`@R#an%52RKxp73xye zeeH<2;pb}Q;1d8`yspHhsKC9dx;o%C42FOH*D-F& z+RD95Qt{Aqb8y2FxMSRup}qS`vI8Lh^#^~%4dAzgDGvNfJoxh;%3uWqH@Nr^e9wjN zy4;m8)`%J7I@*)Bue|R zZHS2K{SDC2+bf1dJfhGZh2|AHpfK?_RoBi+@wWrOf~**Zg^)tavbs#orK$d`JUA)_ z-vnzE4~`1B4ru~tdeM3P8#Jy_58wvHzRS*3Tra)mSuM5Qd~y24nG4sCq>X*)x}$e5 z{Pf}{7w?`)JNuT6eW@d>jQtjRDPmbJ^)9$nBnsdJ>9{`(CQ~Rit za@72fgJMfaGQ@MAJqcFF=D2od(?Al+2k*13tEKj9XNX~O(IemaufJLN^6!_9!$%)e zEEyf2RQKf&2tqZ+D4PN+t*}Q`r*T}8VH5;M$-ALiH_$ZwE2{FBRK=I(GbXzB5ry@C M`%5Zl2bu7H0p9`;5&!@I literal 2086 zcma)6&1)M+6rWvb<&`4Yv6}?fP3yJOLj-k7nxw?01X`Omq{IX*q2P)x0+gq7l0_mZ-i73!(-`ka3i^xd^>Fsacym|A!d$n4H z!20p;f1H+^mK< zE2p*oruUhJnwpswY5SHepaxlNI)0&9YK$Mhl;f+mUz{P|7O(Y>r29%)XGpf(mpwC5 zZ)r&O6v&Q^w!wE#;Fp0vLbqQzfqx44(=Xd{qVH9|I%H$bugUSv@%n_REeo!th;NHl z8nb7RIub%`j~3zr{Ck+hT>^QrA1Nz(B`!IhS2sD<>j@KVJ*y~Vsc9la&Ew;WTY zJPJRL(}bxu4VkxkNa~g;Di&`j)a;p+<)!9b#iVL}djIn6%gq%LrU~O(HNS}T&2+t) zwAu@jQG8{Ar(w*fY!2CcsogPd5JWuEL14}epAGmWcub*lJLHL5?ExUBB$-YnA1bg> z?%w!3<)MziC%r|rx5+q@R5O)%+m>RTu19gCBj8-rl5z!D@Vb{!Eu(F3RgS@HIjJHl zNy z63H7`w=is7D3W%H!4#ZfMPV&2L>eQE5g#4A$-ZE4Bxq)bWB8K5xgje}sl^h?TgJis zDdX}$Y0Rl9^8jiMo>s=$q)9?UteU|>oOQ!Z0K4St0iQJblk129V4mDPE-{#b$3%Q< z(SDymM}2cbP-q-yRVLZJG?J|49XXt3FppL&M7xe1s32x)qm$RAc)+BWd+F6}V^QO@ zbf?Yg1!MEJsXSy-s34BMV`gk47F2(HweCn9D)Mb4Gf3V+avn(5J07mQe5->1ZM+L05mlW!r*iMYP< zk<7|r9=R8gdp5}OAB6$c(6sAVk)Stg+@qeAn8A!-8K~cYd{3U+wo`-VEKPOof>YhI z4=VeW9s$@Z44B=@0KM*59&GP#kG@@Y7WWzlYx`?m0`SlC?;YIQztt@Oelk0U*SZei zZi&oY_{II%?T+EjCGXB=a^-yY!KKGbhf6)&W#;v4KzRbtxq8t!JWbl2JbV?Hy;_KP tU2NyR568o}pX3R;AwPu6I6Bo?Al-sxSx?E-Uu5#%;+(bkoB+ul`42su>Qw*$ diff --git a/main.py b/main.py index 75870c6..026f352 100644 --- a/main.py +++ b/main.py @@ -74,20 +74,23 @@ def main() -> None: # plot_reactor_as_radar(probs=probs) # Generate a random power demand - # demand = generate_demand(n_samples=512) - # print("Demand:\n", demand) # Array con 512 demands + # Generate a random power demand + demand = generate_demand(n_samples=512) + # print("Demand:\n", demand) estado_actual = 5 - # demand = 0.02 - # for i in range(0, 100, 20): i = 90 if True: - demand = i / 100 - matriz_R = ControlModule.generate_R(demand) - - # Ejecutamos para un solo estado - response = ControlModule.control_iteration(matriz_P, matriz_R, estado_actual, gamma) - print("Response 'unitaria':", response) + demand_test = i / 100 + matriz_R = ControlModule.generate_R(demand_test) + + response_test = ControlModule.control_iteration( + matriz_P, + matriz_R, + estado_actual, + gamma + ) + print("Response 'unitaria':", response_test) # Define the number of MDP's states, actions and the discount factor (gamma) n_states = 100 From 0c660e18cdafb21cca96732957f1288f2677d36e Mon Sep 17 00:00:00 2001 From: Angel Date: Fri, 15 May 2026 16:29:54 +0200 Subject: [PATCH 05/11] =?UTF-8?q?P=5Fmax,=20K=20y=20P=20en=20funci=C3=B3n?= =?UTF-8?q?=20de=20P=5Fmax?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Reactor.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/Reactor.py b/Reactor.py index 2fcfb3a..9a65604 100644 --- a/Reactor.py +++ b/Reactor.py @@ -31,20 +31,19 @@ def __str__(self) -> str: def compute_max_power(self) -> np.float64: """ Computes the maximum power of a reactor based on its physical features """ - ### TO BE COMPLETED BY THE STUDENTS ### - ... + return np.float64(self.effective_section * self.neutron_flux * self.core_volume * self.fision_energy) def compute_k(self) -> np.float64: """ Computes the value of the k-constant """ - ### TO BE COMPLETED BY THE STUDENTS ### - ... + self.P_max = self.compute_max_power() + return - np.log(10** -6/self.P_max) def compute_power(self, control_bars_insertion: np.float64) -> np.float64: """ Computes the power delivered (%) by the reactor based on the % of control-bars inserted """ - ### TO BE COMPLETED BY THE STUDENTS ### - ... + k = self.compute_k() + return np.float64(self.P_max * (np.e**(-k * control_bars_insertion))) def compute_control_bars_insertion(self, power: np.float64) -> np.float64: """ Computes the % of controls-bars inserted based on the % of power delivered by the reactor """ ### TO BE COMPLETED BY THE STUDENTS ### - ... + From 93894bb246c9ce1b52ad6e4fb0cdf6c22992bc40 Mon Sep 17 00:00:00 2001 From: Angel Date: Fri, 15 May 2026 17:11:42 +0200 Subject: [PATCH 06/11] Update Metrics.py --- Metrics.py | 28 ++++++++-------------------- 1 file changed, 8 insertions(+), 20 deletions(-) diff --git a/Metrics.py b/Metrics.py index 2515656..a4c25ea 100644 --- a/Metrics.py +++ b/Metrics.py @@ -3,32 +3,20 @@ def MAE(y_true: np.ndarray, y_pred: np.ndarray) -> np.float64: """ Implementation of the Mean Absolute Error (MAE) """ - ### TO BE COMPLETED BY THE STUDENTS ### - - ### DUMMY BEHAVIOUR TO PREVENT CRASHING (MUST BE DELETED AFTER THE FULL IMPLEMENTATION) ### - return 0.0 - ### ### + return np.float64((np.sum(abs(y_pred - y_true)))/len(y_pred)) + # return 0.0 DUMMY def MSE(y_true: np.ndarray, y_pred: np.ndarray) -> np.float64: """ Implementation of the Mean Squared Error (MSE) """ - ### TO BE COMPLETED BY THE STUDENTS ### - - ### DUMMY BEHAVIOUR TO PREVENT CRASHING (MUST BE DELETED AFTER THE FULL IMPLEMENTATION) ### - return 0.0 - ### ### + return np.float64(np.sum(abs((y_pred - y_true)**2))/len(y_pred)) + # return 0.0 DUMMY def R2(y_true: np.ndarray, y_pred: np.ndarray) -> np.float64: """ Implementation of the R2 metric """ - ### TO BE COMPLETED BY THE STUDENTS ### - - ### DUMMY BEHAVIOUR TO PREVENT CRASHING (MUST BE DELETED AFTER THE FULL IMPLEMENTATION) ### - return 0.0 - ### ### + return np.float64(1-(np.sum((y_true - y_pred)**2)/np.sum((y_true - np.mean(y_true))**2))) + # return 0.0 DUMMY def Corr(y_true: np.ndarray, y_pred: np.ndarray) -> np.float64: """ Implementation of the Pearson's Correlation Coefficient """ - ### TO BE COMPLETED BY THE STUDENTS ### - - ### DUMMY BEHAVIOUR TO PREVENT CRASHING (MUST BE DELETED AFTER THE FULL IMPLEMENTATION) ### - return 0.0 - ### ### + return np.float64(np.sqrt((np.sum((y_true - np.mean(y_true))*(y_pred - np.mean(y_pred)))/len(y_true))/((np.sum((y_true - np.mean(y_true))**2))/len(y_true)* (np.sum((y_pred - np.mean(y_pred))**2)/len(y_pred))))) + # return 0.0 DUMMY \ No newline at end of file From 0121b03cc5d60836bc11a654e01b7acccedb9fce Mon Sep 17 00:00:00 2001 From: codeMGL Date: Fri, 15 May 2026 17:17:34 +0200 Subject: [PATCH 07/11] quitado demand_t= --- ControlModule.py | 4 ++-- __pycache__/ControlModule.cpython-312.pyc | Bin 2812 -> 8984 bytes __pycache__/Reactor.cpython-312.pyc | Bin 3065 -> 3065 bytes 3 files changed, 2 insertions(+), 2 deletions(-) diff --git a/ControlModule.py b/ControlModule.py index 41812c4..535b54d 100644 --- a/ControlModule.py +++ b/ControlModule.py @@ -136,7 +136,7 @@ def generate_P(probs: np.ndarray, n_states: np.int32 = 100) -> np.ndarray: @staticmethod def generate_R(demand_t: np.float64, n_states: np.int32 = 100) -> np.ndarray: """Function that generates the rewards (costs) matrix""" - demand = np.float64(demand_t) + demand = np.float64(demand_t) # Debug print("Demanda:", demand) matrix_R = np.zeros((3, n_states, n_states), dtype=np.float64) # (3x)100x100 @@ -324,7 +324,7 @@ def control_loop(demand: np.ndarray, ControlModule._demand_t = np.float64(demand[t]) ControlModule._current_state = current_state - ControlModule._R = ControlModule.generate_R(demand_t = ControlModule._demand_t, n_states) # esta demanda concreta. + ControlModule._R = ControlModule.generate_R(ControlModule._demand_t, n_states) # esta demanda concreta. action = ControlModule.control_iteration(P= ControlModule._P, R= ControlModule._R, current_state = ControlModule._current_state, gamma=gamma) state_increment = np.random.choice( a=action_deltas[action], p=probs[action]) diff --git a/__pycache__/ControlModule.cpython-312.pyc b/__pycache__/ControlModule.cpython-312.pyc index 5d3bcb1dfc4efbcbc1022030cd887f2a5079af54..cfc526d4f316d0cdc7c797ad73561597a13e7191 100644 GIT binary patch literal 8984 zcmeHMYfKzjcCM;^RYP|J25cUF6n>0p8?%1Rn8A;A+Zf{=gMsmkZGzotcNKKE`vFzO z477@bjU1_&A#$TcGRQjuZxSsdTMAY-8jrNnGMXqelSB$_h>R-y7-jt<{;^W~5lxgI zo1Am2tGa>4-jPQ6m1|Rd?>*<0%h_EX6Z7Dc*FC(Kx$IyoER4q`S;RwA=a=E8F}UN&7l8g%=Yv_9&3v zx)`3mNp+b7GwdpE2d#t?WrsZCyK^IHu z7eL=dJj#(rz8;2Wb9f7Z*PdVE?>+N^{3G&(>&$Y!m2&Won`UFB979fbvFbb>i^^gw z+#2Ht!h#<*Yeh0a=>?vS+Z1dx6K`l5h(v}se=I7=;(%X{iJTnc1|)%#0|M7RBnM(q z?sV8ENt_St2gIm9b6Jpyvcex(&BL&qw<+0J=};~g zVWr^uDe$kyBxlML71FP$QA1Jp8T~pAGYVYJN}vs!AZI2yGvq8JXR%W<3Qm$j-^+~M zL|wfIEv)O%nt*0gtcvYMEk%h%3S4>YO<}Z=6e}s&jzQhnCdJB|;{Y%@+(e4J-1We1 z$#P5d*b(TR$e)+L*-VYSZ~h_XVAj?KM$x2OI|rPhK4Y^nGrD!$CO+HMO$ni%O7>McrJ)CyT+srwBiv(yE1*hx!B&Vdp6#!Bfn|i`(_~ypI>8 z#sarRWukqmB@O^AskWYQ%qPEfShYaEy#iuhFeWHy`@d!CyCi-ndxWoM6*ds%qOr|71nNwN$xvRVKZ` zVRyW&qPFdvZn+bCPTB08bfRsU+ch2d%cAF$)!9iW9ZRJ()7_6s4=kN*`?^1MvLo3t z(J}tcv@2EiqgltJvg1qTJEnc9^4&AX9+e+?J=4PaOwXJ+V%mgaye#{FsgC`gYnA3G zdr7inv@q!zbtU^3*vg4(GoI9z`q`~>bn3vV)Rt2z_7qgZVaU&zV(XF2zM5>bZkgyG zKc8Z2rfX+9?$#}^Kbk!;SCcx_lwzBpknQF*e6Q_q?vLi*zmz(AId$lAioFb7!%?@p z=3O<5?Dl2WnFx#(rI{MLBi&3>#TAJn5q$ua<%NhZ%6sJ%77AK5#6tN90O55OibCo5 zG|b4@k|eNO44VYdENkPS^|LDR%o-afk0o^9a5;ueIp~;c!<+y+#k>x!9%vSYQ7j4_ z2S7lZ1MD4y29gv5d%@CS6tl*411ulWqS!#2fD-`}TOnmwKN(afR?mL0OHP7aVqp0S z+v>n8;X<@nhIAJyc$E`hq^u!5*Nt9FSDwbv0HHC<~$9G>)R@B`uuEdYG-4r>S)sTq@sG_>>Wqa0{B!`1qgXL zwdKu6rEjH+sq!jNX+KGx0m_Nf^uy*PJ5e%TIMJ3W-aB*tQSsZ)+}LBoMHeqKkHyVBaU(To=Ky}OpYd-JdreV?DI=1daY;DTHR=;t=Vebr| z56~A770?mFD4IR1eGx^nD9M}kvk#t`h7gpZS*V%c5aRQy#Dn#E->TmwL^%bx2|fWr z$Yz~Alm~-D26E|Kq&2OW6!U%Cy-eF-t6~iqn3?6w-fN)i0iP_H%F<=LgCy}R@Zx1$ zm2~-ODnDHz>jr3GkJqh41Ai>iH-JWPEGl4N0;3uGf-+jhVEmGsPA>k_|A9on;G&yR z&A!1Pxa>Csp^pznB#2Fb`8w#uc5${n77qG{01xb$yDtqy)B*?~$gx=XdTdZNL0&D= z*v{!gsRfcCzb}Z_V_@={#DOS*5Y^g8dUM-Ux?QC^RcFSyg*e3k*t(WpUj(8bb+8^P z?Ngb)poFga3VXLCGc6(>nO1z;O0#zjgA-1iGzlGzP+PF zpj;KVu8m%M?A)7j?w#?>^v~{^tDSf5U2wjWV3yc|Th3AESP^Edqt-E21F7v(?n(C| zyK||iw<(3r#zed@ zM=dE#POV`@C^_*HuW@s-+B$^B(5j?RGP)QfVr}dLMY{}^v z?Dv9y`*S*H+&sPAH|yiRwlB*M*q8P9@2l4|Y|HLK2zBkM1oobx^Rt{@V;5vpv?M81 z2-WYZq+w?kBa`%1I03_K&Whv-C<7X@h_-sWoRJO8Qmi?(#_DJ5Lfck}0>vC^)mm!a znHsp4(-tqls+2=5TSbCm4wZW)b@N*43RJh}oVcB9c|*I`)+si{j%x`!Nq+CM`l%c) zfJk8hMtoMBc4IBGedX*{ZMIg=K~#-_IFTJK_4@_)GPhRWJJdvcA<=Obg+(2gq{qx= z(08C29Pz8`Xa-+6OhP20&_56q1)ek112G#ll@X4``rH`JmH=U;o*^x$68lgPInB~@ zGhzdbD3;&{LB?%*`a9GqZoPt7nQl`~1=^WVVc48R6S>XeVQlgyd{n1@Koo_jti_HA z{$&PI={};ppgoSOT72FJ1W)kCH`og`kgl9`my4*TRjXHvLfgGY7`gzmID;j>X9@hI?;S*msso(2X`8;anpz)EPc!gydiODZj^ zR;@-YB5`#H4?(;r0&69qF_O`c_@hOI9J&(3os3uISCB~CkGk<^@B!UeOmVx$oQv${ zrR{t0ciYmAI{f9HK*XqZymh+ym)k$v^I*?n+3|#ZxuR-{on)s=?i41RPb#aYDkdv_ zUiZuEp9LNS7At@J_Zz;d_@d$;>lUy5;p1yPscSupp5Dbk=y9Mg73h1^Km3xSnwb{n zJN#*5o@qZ&8Kf^TX^N_BXA;ik>e|QEhf>vtX3soqSg3A7St?Lne{fqFSEfDVf4oq# z`vp~KKalV|d876hp}V2kjSuMul?!hi9cx{5>{!~{fWNg*w$(k}b|kg!$giA_8!n_8 zF3fjaS=iQvw9U)q+aH%7PL&@-<)5!UZoQgny}B4vK)OEqJW7W@Al60Z3Nv4R z7@j&db}@HGM6 zdJ{}5T|pIB&a)N7JZN*%PkGpDHy!*BfC9G@A87DR$?H`Myxs_Tw1T|D>+J_)N~4IV zY*h!InV>%+z>^eD+&a}3<+X#}3=HJ7Yn~pSV{(SY_b`WcvSys1B|vT$`Cz^DJNVqB zo->Ta3BHhZ!!ssMTPh~DPwko9lcq2~%Y4Q@V72_*!LJ&>XnaoUUzmVBP1w_en19Bc zVxVc~RM%uzn!^0-!Ot2VG-&y`>tBVw2t6nDnWo5+ZZVUlTP?L;GylSV!#<}Vn_++) zIwLsD7E7hxGh~;lw@fuoHl|IG{}whdFd0Z&AP?a#Zu`{W&vyxkG)z z=pM$OfM=DDryC?1W`<3}@M@52uYq<3DKUYkHHbP3tOOEa$>-Afj5#+vjJZ$sS*_6%wmk`y*eLr!?1xy-W=(PBYvusVkAvSDZO zcEx&!;tLczd<%6=%{%UxKQW?QG&Kfh^jI7(s{q{|j_!4#l;S;x1|o|CP05_~BCPsYRMNNIt;vSexMq#CYOq5~s09t(<=kfc4 zU;b^Bi`(nCGvpBfvW*=0TR!eM7dO>6^b8yT&VT>epK`7EiV^2PVZ{7j-j(pdw-NL; z4-d`o;LP;|@rx~^+=2RmYHan%nBJ=G!it9V?gA~;a~9ZCa}<0)@f>7{>J#^2rRBqb zAPVAk%$vf3=FLGfRNbmcj)hg)qtfrHv{z+fH#8N^im?y9AcCZtuKReEmLzSOx~|d` zRk1%H_-}y1kNE>WoQD{IK0DxCN&f|z5$XvzP%LF<|~>P&CN@$qFcSAy+1uW;u+~r)GXPZBd3;~81yR8x~4;B-mwF@vXpU;4j46AB=x65tu2Rj4o{2H)~35I<(+AJlm0S9Z8s$ z3tb86;|9RyV@o9)6Q`ES(PebD(Mfu_xO}2u`oQ##=_@n+Gyd7ihYfR{&yUVY^NsJ! zZ#@66#jVdw$eAug8J@W)wt$3BJd5lpKplr`r8OP1HUd6X4!Oy{vUo!WgM#3y-On2$ZIF|vubO2&3cN;S#$ zJ?|8ENhU!uzPENg(SR-ieUF>;w9eG#|L`D6KFASsGr8=j*0UG=<(4z-M25x`CKTuN;Lf`T}~UV+$k!J{H9 zc-L8FdfS1WuGKU>87-ajA{0$W-LXwU=n)@v!W>tEzBmZHIo_kTuvOw8*6!z!*A!9p zSVSE`SH6W4LId6i@M0yQ6uSluNm@(L9aKY!t9`>=>R1dqXTJyx2P$vx2@C=~zdPKEdz#TgyTdGHl*^IDBQ@Ba>lNJC_!J_UJ z^LUenn~^qYItc~-zOm>+UDdHWn9Rfj-p&W24<+mdiAmE@y_25?VxAbb7l@$mDh@<< zO3mYpNXgEG8%RvIuZ#sPjG$O5NW(ASjDVYFwI|a>Bv2Y8BcGHpzq`I~rN93`;=CAJ`T#%k0e5X9K+P%>7OiC}<>(b#r+q#y>-Gfz>Xc}A? zde+%f?!Ws(6~)toOWel(jxu-W3p zTvqR78%r$BEMF=Q-Tl6S#LF}&stBO_RTpY%dD`;s{ed6PJ~+EB^-}|<7F#!34=l72 znivg!b&GW!ZKB%9>VI1e5PU-=I%HpST+fXXKT6_oJ_IcJuZ4QNVW`KC!oe+Ck~dLB zLQSdvXX7XJHfoNQFVHe|Y)oe@9ngJ-M{5_#b8qX~YI4f9oPVUg9^Mf6m^vv_K%ETY zI<{0>C*5=-*^2AYRtcJH1x709bp!FLYFVmAXIr#rpjpJf0%Qffea@B8f6O^9>TBp; z-oMfyDynI+1}1o>U=Xa!+b*+2hT$?DC!e8@lM6v}CyQarJ2x;NJG{BNE3Qseg& z=6HLBKOm;cf`W>;tTW;ctcH}ERzdj`YNN2*>pCu From 04868bd72b48ebf2c57a959e602d31644ceb1da4 Mon Sep 17 00:00:00 2001 From: Arleking558 Date: Fri, 15 May 2026 18:18:23 +0200 Subject: [PATCH 08/11] Errores resueltos --- ControlModule.py | 40 ++++++++++++++++++++-------------------- Metrics.py | 9 ++++++--- 2 files changed, 26 insertions(+), 23 deletions(-) diff --git a/ControlModule.py b/ControlModule.py index 535b54d..a3a5277 100644 --- a/ControlModule.py +++ b/ControlModule.py @@ -94,8 +94,7 @@ def generate_P(probs: np.ndarray, n_states: np.int32 = 100) -> np.ndarray: matrix_P[2][n_states - 1][n_states - 1] += probs_increase[1] + probs_increase[2] - print("Probabilidades:\n") - print(matrix_P) + """ [[[0.8 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [0.025 0.8 0. 0. 0. 0. 0. 0. 0. 0. ] @@ -137,7 +136,6 @@ def generate_P(probs: np.ndarray, n_states: np.int32 = 100) -> np.ndarray: def generate_R(demand_t: np.float64, n_states: np.int32 = 100) -> np.ndarray: """Function that generates the rewards (costs) matrix""" demand = np.float64(demand_t) # Debug - print("Demanda:", demand) matrix_R = np.zeros((3, n_states, n_states), dtype=np.float64) # (3x)100x100 # quitar DEBUGGING @@ -210,9 +208,7 @@ def generate_R(demand_t: np.float64, n_states: np.int32 = 100) -> np.ndarray: matrix_R[2][estado_inicial][estado_final] = -coste - print() - print("DEMANDA ACTUAL: ", demand) - print(matrix_R) + # Creamos unos tests # print() @@ -282,22 +278,28 @@ def control_iteration(P, R, estado_actual, gamma, max_iter=1000) -> np.int32: # -- Comprobacion provisional de si la matriz P es estocastica # La última fila tiene todo ceros, rellenamos con uno en la posicion final - P[2][-1][-1] = 1 - - for a in range(3): - P[a] = P[a] / P[a].sum(axis=1, keepdims=True) # renormaliza # Cambiar si hay tiempo, en vez de distribuir las probabilidades uniformemente, truncarlas a los bordes - # print() # print("P normalizada\n", P) # print(check_stochastic(P)) # print() - + """ #politica optima con la libreria .... pi = mdptoolbox.mdp.PolicyIteration(P, R, gamma, max_iter=max_iter) pi.setVerbose() pi.run() print("Policy: ", pi.policy) - return pi.policy[estado_actual] + return pi.policy[estado_actual]""" + + mdp = mdptoolbox.mdp.ValueIteration( + transitions=P, + reward=R, + discount=gamma, + max_iter=max_iter + ) + + mdp.run() + + return np.int32(mdp.policy[estado_actual]) @@ -309,9 +311,7 @@ def control_loop(demand: np.ndarray, gamma: np.float64,) -> np.ndarray: """Function that computes all the required iterations (control-loop) to satisfy the power demand""" - ControlModule._probs = probs - ControlModule._n_states = n_states - ControlModule._n_actions = n_actions + respuesta = np.zeros_like(a=demand, dtype=np.float64) # Almacena las acciones para cada demanda current_state = np.int32(0) # Estado inicial, se puede modificar si se desea empezar en otro estado @@ -325,15 +325,15 @@ def control_loop(demand: np.ndarray, ControlModule._current_state = current_state ControlModule._R = ControlModule.generate_R(ControlModule._demand_t, n_states) # esta demanda concreta. - action = ControlModule.control_iteration(P= ControlModule._P, R= ControlModule._R, current_state = ControlModule._current_state, gamma=gamma) + action = ControlModule.control_iteration(P= ControlModule._P, R= ControlModule._R, estado_actual = ControlModule._current_state, gamma=gamma) state_increment = np.random.choice( a=action_deltas[action], p=probs[action]) current_state = current_state + state_increment current_state = np.int32(np.clip(a=current_state, a_min=0, a_max=n_states - 1)) - respuesta[t] = current_state / 100.0 + respuesta[t] = current_state / n_states return respuesta - +""" def check_stochastic(P, tol=1e-6): P = np.array(P) @@ -353,5 +353,5 @@ def check_stochastic(P, tol=1e-6): print("❌ Matriz NO estocástica.") return ok - +""" diff --git a/Metrics.py b/Metrics.py index a4c25ea..e31c215 100644 --- a/Metrics.py +++ b/Metrics.py @@ -17,6 +17,9 @@ def R2(y_true: np.ndarray, y_pred: np.ndarray) -> np.float64: # return 0.0 DUMMY def Corr(y_true: np.ndarray, y_pred: np.ndarray) -> np.float64: - """ Implementation of the Pearson's Correlation Coefficient """ - return np.float64(np.sqrt((np.sum((y_true - np.mean(y_true))*(y_pred - np.mean(y_pred)))/len(y_true))/((np.sum((y_true - np.mean(y_true))**2))/len(y_true)* (np.sum((y_pred - np.mean(y_pred))**2)/len(y_pred))))) - # return 0.0 DUMMY \ No newline at end of file + """Implementation of the Pearson's Correlation Coefficient""" + + if np.std(y_true) == 0 or np.std(y_pred) == 0: + return np.float64(0.0) + + return np.float64(np.corrcoef(y_true, y_pred)[0, 1]) \ No newline at end of file From 66f5452b11180ce15a5df80bdaa32a6bc4597ed9 Mon Sep 17 00:00:00 2001 From: Arleking558 Date: Fri, 15 May 2026 19:34:27 +0200 Subject: [PATCH 09/11] Codigo Limpio --- ControlModule.py | 280 ++++++---------------------------------- Reactor.py | 6 +- main.py | 34 ++--- mdp-reactor-nuclear.zip | Bin 0 -> 50845 bytes 4 files changed, 50 insertions(+), 270 deletions(-) create mode 100644 mdp-reactor-nuclear.zip diff --git a/ControlModule.py b/ControlModule.py index a3a5277..931d5df 100644 --- a/ControlModule.py +++ b/ControlModule.py @@ -11,33 +11,13 @@ def __init__(self): @staticmethod def generate_P(probs: np.ndarray, n_states: np.int32 = 100) -> np.ndarray: """Function that generates the probabilities (transition) matrix""" - matrix_P = np.zeros((3, n_states, n_states), dtype=np.float64) # cambiar a 100x100 + matrix_P = np.zeros((3, n_states, n_states), dtype=np.float64) - #saca las probabilidades del Json probs_decrease = probs[0] probs_maintain = probs[1] probs_increase = probs[2] - #Las probabilidades en los bordes se pierden, así que las filas no suman 1. (Y creo que normalizar es un apaño que no debería estar bien) - #Ej: COn decrease en estado 1, sumas las probabilidad de bajar a 0, pero no a la de bajar a -1 porque nunca se cumple en los IFs - - #Tantas iteraciones creo que son ineficientes, y mi pc es una patata, así que creo que hay que aprovechar mejor el hecho que - #sabemos que solo hay que modificar 3 IJs, no hace falta comprobar los 100, el 97% de las comprobaciones son inncesarias - # ---------------- DECREASE ---------------- - """ - - DECREASE ANGEL - for estado_inicial in range(n_states): - dest_0 = max(0, estado_inicial - 2) - dest_1 = max(0, estado_inicial - 1) - dest_2 = estado_inicial - - matrix_P[0][estado_inicial][dest_0] += probs_decrease[0] - matrix_P[0][estado_inicial][dest_1] += probs_decrease[1] - matrix_P[0][estado_inicial][dest_2] += probs_decrease[2] - - """ for estado_inicial in range(n_states): for estado_final in range(n_states): if estado_inicial == estado_final: @@ -68,290 +48,104 @@ def generate_P(probs: np.ndarray, n_states: np.int32 = 100) -> np.ndarray: matrix_P[2][estado_inicial][estado_final] = probs_increase[2] # ---------------- CORRECCIÓN DE BORDES ---------------- - # Si una acción intenta llevar el reactor por debajo del estado 0 - # o por encima del estado n_states - 1, esa probabilidad no se pierde: - # se acumula en el borde correspondiente. - - # DECREASE: - # En estado 0, los movimientos -2 y -1 se quedan en 0. matrix_P[0][0][0] += probs_decrease[0] + probs_decrease[1] - - # En estado 1, el movimiento -2 se queda en 0. matrix_P[0][1][0] += probs_decrease[0] - # MAINTAIN: - # En estado 0, el movimiento -1 se queda en 0. matrix_P[1][0][0] += probs_maintain[0] - - # En el último estado, el movimiento +1 se queda en el último estado. matrix_P[1][n_states - 1][n_states - 1] += probs_maintain[2] - # INCREASE: - # En el penúltimo estado, el movimiento +2 se queda en el último estado. matrix_P[2][n_states - 2][n_states - 1] += probs_increase[2] - - # En el último estado, los movimientos +1 y +2 se quedan en el último estado. matrix_P[2][n_states - 1][n_states - 1] += probs_increase[1] + probs_increase[2] - - - """ - [[[0.8 0. 0. 0. 0. 0. 0. 0. 0. 0. ] - [0.025 0.8 0. 0. 0. 0. 0. 0. 0. 0. ] - [0.175 0.025 0.8 0. 0. 0. 0. 0. 0. 0. ] - [0. 0.175 0.025 0.8 0. 0. 0. 0. 0. 0. ] - [0. 0. 0.175 0.025 0.8 0. 0. 0. 0. 0. ] - [0. 0. 0. 0.175 0.025 0.8 0. 0. 0. 0. ] - [0. 0. 0. 0. 0.175 0.025 0.8 0. 0. 0. ] - [0. 0. 0. 0. 0. 0.175 0.025 0.8 0. 0. ] - [0. 0. 0. 0. 0. 0. 0.175 0.025 0.8 0. ] - [0. 0. 0. 0. 0. 0. 0. 0.175 0.025 0.8 ]] - - [[0.6 0.35 0. 0. 0. 0. 0. 0. 0. 0. ] - [0.05 0.6 0.35 0. 0. 0. 0. 0. 0. 0. ] - [0. 0.05 0.6 0.35 0. 0. 0. 0. 0. 0. ] - [0. 0. 0.05 0.6 0.35 0. 0. 0. 0. 0. ] - [0. 0. 0. 0.05 0.6 0.35 0. 0. 0. 0. ] - [0. 0. 0. 0. 0.05 0.6 0.35 0. 0. 0. ] - [0. 0. 0. 0. 0. 0.05 0.6 0.35 0. 0. ] - [0. 0. 0. 0. 0. 0. 0.05 0.6 0.35 0. ] - [0. 0. 0. 0. 0. 0. 0. 0.05 0.6 0.35 ] - [0. 0. 0. 0. 0. 0. 0. 0. 0.05 0.6 ]] - - [[0. 0.2 0.8 0. 0. 0. 0. 0. 0. 0. ] - [0. 0. 0.2 0.8 0. 0. 0. 0. 0. 0. ] - [0. 0. 0. 0.2 0.8 0. 0. 0. 0. 0. ] - [0. 0. 0. 0. 0.2 0.8 0. 0. 0. 0. ] - [0. 0. 0. 0. 0. 0.2 0.8 0. 0. 0. ] - [0. 0. 0. 0. 0. 0. 0.2 0.8 0. 0. ] - [0. 0. 0. 0. 0. 0. 0. 0.2 0.8 0. ] - [0. 0. 0. 0. 0. 0. 0. 0. 0.2 0.8 ] - [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.2 ] - [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]]] - """ - return matrix_P @staticmethod def generate_R(demand_t: np.float64, n_states: np.int32 = 100) -> np.ndarray: """Function that generates the rewards (costs) matrix""" - demand = np.float64(demand_t) # Debug - - matrix_R = np.zeros((3, n_states, n_states), dtype=np.float64) # (3x)100x100 - # quitar DEBUGGING - # son recompensas, poner costes negativos - - # Se calcula la matriz de distancias entre el estado actual s y el estado futuro s' - # Aunque no se pueden alcanzar algunos estados (ej, pasar de s a s + 40), - # como en la matriz de probabilidades esa probabilidad es nula, el coste es indiferente - # Entonces, calculamos la distancia entre 'estado_inicial' y 'estado_final' y la duplicamos - # si la acción elegida se aleja del objetivo 'demand' + demand = np.float64(demand_t) + matrix_R = np.zeros((3, n_states, n_states), dtype=np.float64) - - #al final hacer otra iteración exterior, y probar con np.where # ---------------- DECREASE ---------------- for estado_inicial in range(n_states): for estado_final in range(n_states): - - delta_t = demand - (estado_final/ n_states) - + delta_t = demand - (estado_final / n_states) coste = abs(delta_t) - # Comprobamos si la acción nos aleja del estado final (demanda) - # Si la distancia hasta la demanda es mayor en el estado actual vs. en el estado final, - # nos hemos alejado. Hay que multiplicar x2 la distancia - if demand > (estado_inicial/ n_states): - # print( - # "Nos alejamos (decrease)", - # estado_inicial, - # estado_final, - # "d_t", - # delta_t, - # ) + + if demand > (estado_inicial / n_states): coste *= 2 - # else: - # delta_t = abs(delta_t) # Para DEBUGING - # Rellenamos la matriz para estado_final x estado_inicial + matrix_R[0][estado_inicial][estado_final] = -coste # ---------------- MAINTAIN ---------------- for estado_inicial in range(n_states): for estado_final in range(n_states): - - delta_t = demand - (estado_final/ n_states) - + delta_t = demand - (estado_final / n_states) coste = abs(delta_t) - matrix_R[1][estado_inicial][estado_final] = -coste # ---------------- INCREASE ---------------- for estado_inicial in range(n_states): for estado_final in range(n_states): - - delta_t = demand - (estado_final/ n_states) - + delta_t = demand - (estado_final / n_states) coste = abs(delta_t) - # Comprobamos si la acción nos aleja del estado final (demanda) - # Si la distancia hasta la demanda es mayor en el estado actual vs. en el estado final, - # nos hemos alejado. Hay que multiplicar x2 la distancia - if demand < (estado_inicial/ n_states): - # print( - # "Nos alejamos (increase)", - # estado_inicial, - # estado_final, - # "d_t", - # delta_t, - # ) + + if demand < (estado_inicial / n_states): coste *= 2 - # else: - # delta_t = abs(delta_t) # Para DEBUGING matrix_R[2][estado_inicial][estado_final] = -coste - - - # Creamos unos tests - # print() - # for estado_inicial in range(1, 3): - # for estado_final in range(estado_inicial - 1, estado_inicial + 2): - # coste0 = matrix_R[0][estado_final][estado_inicial] - # coste1 = matrix_R[1][estado_final][estado_inicial] - # coste2 = matrix_R[2][estado_final][estado_inicial] - # print( - # f"ACCIÓN: decrease. Para ir desde {estado_inicial} hasta {estado_final} --> Coste: {coste0}" - # ) - # print( - # f"ACCIÓN: maintain. Para ir desde {estado_inicial} hasta {estado_final} --> Coste: {coste1}" - # ) - # print( - # f"ACCIÓN: increase. Para ir desde {estado_inicial} hasta {estado_final} --> Coste: {coste2}" - # ) - """ - [[[ *4. *4. *4. *4. *4. *4. *4. *4. *4. *4.] - [ *2. *2. *2. *2. *2. *2. *2. *2. *2. *2.] - [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] - [ 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] - [ 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.] - [ 3. 3. 3. 3. 3. 3. 3. 3. 3. 3.] - [ 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.] - [ 5. 5. 5. 5. 5. 5. 5. 5. 5. 5.] - [ 6. 6. 6. 6. 6. 6. 6. 6. 6. 6.] - [ 7. 7. 7. 7. 7. 7. 7. 7. 7. 7.]] - - [[ *4. *4. *4. *4. *4. *4. *4. *4. *4. *4.] - [ *2. *2. *2. *2. *2. *2. *2. *2. *2. *2.] - [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] - [ *2. *2. *2. *2. *2. *2. *2. *2. *2. *2.] - [ *4. *4. *4. *4. *4. *4. *4. *4. *4. *4.] - [ *6. *6. *6. *6. *6. *6. *6. *6. *6. *6.] - [ *8. *8. *8. *8. *8. *8. *8. *8. *8. *8.] - [*10. *10. *10. *10. *10. *10. *10. *10. *10. *10.] - [*12. *12. *12. *12. *12. *12. *12. *12. *12. *12.] - [*14. *14. *14. *14. *14. *14. *14. *14. *14. *14.]] - - [[ 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.] - [ 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] - [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] - [ *2. *2. *2. *2. *2. *2. *2. *2. *2. *2.] - [ *4. *4. *4. *4. *4. *4. *4. *4. *4. *4.] - [ *6. *6. *6. *6. *6. *6. *6. *6. *6. *6.] - [ *8. *8. *8. *8. *8. *8. *8. *8. *8. *8.] - [*10. *10. *10. *10. *10. *10. *10. *10. *10. *10.] - [*12. *12. *12. *12. *12. *12. *12. *12. *12. *12.] - [*14. *14. *14. *14. *14. *14. *14. *14. *14. *14.]]] - """ return matrix_R @staticmethod def control_iteration(P, R, estado_actual, gamma, max_iter=1000) -> np.int32: """Function that computes one control-iteration""" - - # ECUACION BELLMAN: V(s) = max_a [ sum_s'( P[s,s',a] * (R[s,s',a] + gamma * V[s']) ) ] - - - # print("P y R", P.shape, R.shape) - # print() - # print(P[0][:3][:3]) - # print("P\n", P) - # print(R[0][:3][:3]) - # print() - - # -- Comprobacion provisional de si la matriz P es estocastica - # La última fila tiene todo ceros, rellenamos con uno en la posicion final - # print() - # print("P normalizada\n", P) - # print(check_stochastic(P)) - # print() - """ - #politica optima con la libreria .... - pi = mdptoolbox.mdp.PolicyIteration(P, R, gamma, max_iter=max_iter) - pi.setVerbose() - pi.run() - print("Policy: ", pi.policy) - return pi.policy[estado_actual]""" - mdp = mdptoolbox.mdp.ValueIteration( transitions=P, reward=R, discount=gamma, - max_iter=max_iter + max_iter=max_iter, ) mdp.run() return np.int32(mdp.policy[estado_actual]) - - @staticmethod - def control_loop(demand: np.ndarray, - probs: np.ndarray, - n_states: np.int32, - n_actions: np.int32, - gamma: np.float64,) -> np.ndarray: - + def control_loop( + demand: np.ndarray, + probs: np.ndarray, + n_states: np.int32, + n_actions: np.int32, + gamma: np.float64, + ) -> np.ndarray: """Function that computes all the required iterations (control-loop) to satisfy the power demand""" + respuesta = np.zeros_like(a=demand, dtype=np.float64) + current_state = np.int32(0) + action_deltas = [ + np.array([-2, -1, 0], dtype=np.int32), + np.array([-1, 0, 1], dtype=np.int32), + np.array([0, 1, 2], dtype=np.int32), + ] - respuesta = np.zeros_like(a=demand, dtype=np.float64) # Almacena las acciones para cada demanda - current_state = np.int32(0) # Estado inicial, se puede modificar si se desea empezar en otro estado - - action_deltas = [np.array([-2, -1, 0], dtype=np.int32), np.array([-1, 0, 1], dtype=np.int32),np.array([ 0, 1, 2], dtype=np.int32)] ControlModule._P = ControlModule.generate_P(probs, n_states) - for t in range(demand.shape[0]): # La demanda cambia en el tiempo - + for t in range(demand.shape[0]): ControlModule._demand_t = np.float64(demand[t]) - ControlModule._current_state = current_state + ControlModule._current_state = current_state + ControlModule._R = ControlModule.generate_R(ControlModule._demand_t, n_states) - ControlModule._R = ControlModule.generate_R(ControlModule._demand_t, n_states) # esta demanda concreta. - action = ControlModule.control_iteration(P= ControlModule._P, R= ControlModule._R, estado_actual = ControlModule._current_state, gamma=gamma) + action = ControlModule.control_iteration( + P=ControlModule._P, + R=ControlModule._R, + estado_actual=ControlModule._current_state, + gamma=gamma, + ) - state_increment = np.random.choice( a=action_deltas[action], p=probs[action]) + state_increment = np.random.choice(a=action_deltas[action], p=probs[action]) current_state = current_state + state_increment current_state = np.int32(np.clip(a=current_state, a_min=0, a_max=n_states - 1)) respuesta[t] = current_state / n_states - return respuesta - -""" -def check_stochastic(P, tol=1e-6): - - P = np.array(P) - assert P.ndim == 3, f"P debe ser (A, S, S), tiene forma {P.shape}" - A, S, _ = P.shape - ok = True - for a in range(A): - row_sums = P[a].sum(axis=1) - bad = np.where(np.abs(row_sums - 1.0) > tol)[0] - if len(bad) > 0: - ok = False - for s in bad: - print(f" Acción {a}, Estado {s}: suma = {row_sums[s]:.8f}") - if ok: - print("✅ Matriz estocástica: todas las filas suman 1.") - else: - print("❌ Matriz NO estocástica.") - return ok - -""" + return respuesta diff --git a/Reactor.py b/Reactor.py index 9a65604..691c072 100644 --- a/Reactor.py +++ b/Reactor.py @@ -45,5 +45,7 @@ def compute_power(self, control_bars_insertion: np.float64) -> np.float64: def compute_control_bars_insertion(self, power: np.float64) -> np.float64: """ Computes the % of controls-bars inserted based on the % of power delivered by the reactor """ - ### TO BE COMPLETED BY THE STUDENTS ### - + min_power = 10**-6 / self.max_power + power = np.clip(a=power, a_min=min_power, a_max=1.0) + insertion = -np.log(power) / self.k + return np.float64(np.clip(a=insertion, a_min=0.0, a_max=1.0)) diff --git a/main.py b/main.py index 026f352..66815b1 100644 --- a/main.py +++ b/main.py @@ -9,7 +9,7 @@ from Plotter import * -def get_args() -> tuple[Reactor, np.float64, np.float64]: # 3 salidas (?) +def get_args() -> tuple[Reactor, np.float64, np.float64]: # Define the parser object parser = argparse.ArgumentParser() @@ -45,7 +45,7 @@ def get_args() -> tuple[Reactor, np.float64, np.float64]: # 3 salidas (?) ) # Some verbose of the reactor loaded - print(reactor) # Overloaded in the __str__ method of Reactor's class + print(reactor) # Return the Reactor object, gamma and the random seed return reactor, args.gamma, args.random_seed @@ -68,41 +68,25 @@ def main() -> None: dtype=np.float64, ) - matriz_P = ControlModule.generate_P(probs, 100) - # Make a radar-plot with the reactor probabilities - # plot_reactor_as_radar(probs=probs) + plot_reactor_as_radar(probs=probs) - # Generate a random power demand # Generate a random power demand demand = generate_demand(n_samples=512) - # print("Demand:\n", demand) - - estado_actual = 5 - i = 90 - if True: - demand_test = i / 100 - matriz_R = ControlModule.generate_R(demand_test) - - response_test = ControlModule.control_iteration( - matriz_P, - matriz_R, - estado_actual, - gamma - ) - print("Response 'unitaria':", response_test) - # Define the number of MDP's states, actions and the discount factor (gamma) + # Define the number of MDP's states and actions n_states = 100 n_actions = 3 # Get the response time-series (answer to the demand time-series) response = ControlModule.control_loop( - demand=demand, probs=probs, n_states=n_states, n_actions=n_actions, gamma=gamma + demand=demand, + probs=probs, + n_states=n_states, + n_actions=n_actions, + gamma=gamma, ) - # print("Response:\n", response) # Array con 512 respuestas - # Plot the original power demand plot_demand(demand=demand) diff --git a/mdp-reactor-nuclear.zip b/mdp-reactor-nuclear.zip new file mode 100644 index 0000000000000000000000000000000000000000..0a0e399aebc2fe6fe35c2e95f5cddefbc151e48f GIT binary patch literal 50845 zcmb5UL$GK;x2?Ht+qU=Fwr$(CZQHhO+qP}nXIpjtTNQC{-PU`nk!xkdh-~Csjrq+n zh!j$rbZeoonMc{ za)YPehvo>$Dik#H3#%*Ey%7u;g9|P{dO4Sq_sEo z&e~2Avl7WvJt`=f4AaKibv>RWUj|UA`)tQpPo8+`^77SheVfA{ebo(8N^y}``*BDV zjrd^^P23w8d#^vqUTK0F&O|;V>Wi}jXOnGjvoGFD z#$ZqEoou(B-H%wTsxC^813zy*+65;6bgEc7LpjI(voO@v!D*o|F(76zObQ@ z6)f&BJ);4bBq)X81V65l9ceT{&Y+NEVhQhBpBSTT_v=+p_h2ns#wWfUJ(Td%76aev z_}++qr`#Hdxmh{Gb9HK!dMz*PVC@ve-lksg-k#gW+#bc=mcGI7WSz;_*?mL(r+WWy zjdPOf{(m)w`ByzXJ$nx$10!=2J-z=?XM*TDZM*7=^K5tk01RRP0LuTR6tc5*cC@pW zu`_nDHlZ`J_i#41v!!KX{6qKev>OdAYiu>oT+CURFtr_%OGAreDVoj3)wbbCZ;ho? z6Az<_#KYlMsw`Y5Q6iIFW+q7z!(~r&U4KdCbU|qKK(_jBn+7QB{6C!{IQju9@rFs( z(8|rq>FJk*=N|>W>i)eblppms9A@l7lO!h~;#+C4CrrQDzrSJ3u?gJ2IPX!#cPEg6SPaJ=@peHTb3h*_2494j<#ft!O?$yG zlv!b|r~1+wXE;o>1U5*yHQJ*I1I-+*ZTeRg+n4XQovaePwGwpEu3~c9^vGk-Hpy{r z6;~XHkWxp)>$wKwKtb_yS;XQXWeFqDK-0~MR#zF(pk(TM#nlr>81rxtH8=)iC4i_u zL0jDfU&IElBM8VPXuyBg8PJi6DbA2Rc}KU~B_##z#RF67e{%sYB3=gENCc9>gtE`0 zU|b;^sjpw*{$Ys4aXW>HM&TSQzZdY5bIQOt}|I zZ1OjQY8Bik_OEj`26?{mw1MfZ!2}~e2yapxdp1*BW-43H^=l%nMm=0(drhDfQzG@* zn*!)pc{Huca@^*8+;)GQ_Svr@W}`@+RqpB;*4%_jxOFY4cTY^mOapO_kpX}+fT0N= zc@qF6C?XIXkRF8Kb_zcO0XV;G6OE$*B|yxMb$HT9{~Iy{AN~>5$P}`IK=Awn$m#=# zWW%EFPTrZjfMr`JWeS9t3KbZDy5fw`00eoT1AhsB#}LqPH@~_oQz0ySRJ@I{Xw~Re z6(=Hrr^f&|fEv=@I|Oj;e!p)ngwpQJCYw^5R=ZpcbKTm#-zuLz+uvrNmO`f4VvFil2Uq7= zRJT$~EM26<6W6G35!%&uo8BwcV@%W5zef+bEUus3zw}bM2F*qIa_dT{GEYJ9)Gwv6 z&)yuF+^cElnKxM;R=XM1Ujv2T1}6#1E2|!K8ZAZ#_d}IIb*!3{wttrat)#Q`nRJzZ z{Hk4@L< z#T(Q4d)F`1&YGAlQ8k&wI1GO;ky;c^7|OXFs#j3!5@{~To`fKvpY}XEKDVgS`S5^lLuXf!@mOXa*G;@T<&%WbY=Dm zbv18ZbZ@o=ZoVe1|8ABefQp6Ca8iv{&d1@n%i`JTwkbw?k2J@=O8dX2ulL5d%c1O^ z3F9&3(i2T2%DOU^0@pay>d*whs6Zt_rT0sYI0}2=PIJd>tY!S4&9_v@dQ8}Mp38xWo!^s*8+bTo9N}-W@{DCLz&kaod@EQ(2;o%O~ zy<7Ou4dhap5P-S`$Bg|z{qw&rkaQW z?@#!yBoKb*7k50e$B*+W(2Kr^0RirX^aoF!-c1`7tYjG$I}s;wC;d~K>0g~{Zd&$S zsa!=eeg}1!0s&zmbEu9X5H2Z}kg-*L*h95j$N6V~_AF_Pzgh~`sQ`hzm0u}aDZ-rp zkJ$o{Kj=99l=GC8Gxcga9@8TXjy|b`(|Ils2w!a`JDy!8(b2b(bfNd;PkS{d$U{{YdP_Y04&B5~|%OkHc+HY%`>KFU44PEp}QN z9iQX1B&RE_JDMnZ8$Am>mun6CZX`10sOJ#lNjhBs5RZr)6PCD5$4#zN-JdN*QWG|^ zQZiFEaw%6DCa!r`#J0g7jX`^6{wmUk{VA{U!n*=bXS>4v4fg@7Yn+hr36??#Qu5Y& zEH3h?2x3*5I|YV}g%iOv^H+OfyaGo_DrX+E?0vEeoGBnHFARc#Ifsmm@@B7M-8<^3 z!rrmQb347&dW*$u%$gV||}CIY)ZG!2Mqxn5vgO&A&d^nf0T3AETV82wCS zYxHHNz)dKe$X#E zxEl)nd)4aATW~i%ag44(Dg~%k*C#8w`Xk5YSx+rcu(O~bb^(kRJXbr@wKKOQ`E7); zMJQm7+CX_^68E6#19)A^!IDiJRiZwDQ$_XX+%0WFi$i?35m%ZlMVqX`rsPo4zQ6$2x+X` z#vp$GpjNmmzl!CdoWX()o{z-9TLC*hq<2Y@AkMhAXqiW?wTHHPHnnR;1j;(69l?Ut zLQ#yWRPd|s_Q9i2Ms5c~pm#{Rd$BAQk?^6T@_T$3D7bhs*~~P%gqiNqw-$kg_fsJ$ zTIBH|-9u7cK!X+9Vhg@q(kbxs!Mg^ka23b?!i4YgGtV9`_NsO4JtU4|uDd)NB zp#VKgwbBTgYAn_!5oAxO+ZY4}ck;wz`-Ki#B0w~5mEl%N`&drNzS&CK*aXP^|?O^Y)R^XTSq}JK6|mHV4CbQJ(W?LaOE1YPVKgzJy^~3VRSGx zs(%CY=`}rW7F^t_e6B9nvcB2+Zm&O!>x>T;bF6T=phk_cmiG@Lv_aVIEez0L^fUX_KT@!Euk(Ujh>CzPAVT;HE=! zTbj$THq+he7AAel`n-QQ?!DK=(y%qxGHXZY$*gpt&yroik94hFV%xgSxO}Vp($DU3 z^f7cCzdsH>hD7PX*4pB$SNs;ayH(bD96THEQ5=BS`&c-=^EE9$)(y>LKwrOk}E+mv&1SZwavcNqi$er9Q<38;%By-!}k>%&NfZF(xRzWZkpPvy#j5|6)n! z|IOx;0mQ`VSf_0pI~!dxeH+n#h8>?R^l54Sk@+0EncH7;su|%&s`>D*?&4{@aLA1G zE88L-e}3;&`fxQ?T}}9YB6r@X#wTiLi{WF@m%B+HeTm2 z6+yUtcn8Y+wMmcB&|FVC*9IAOMaynue}TYeYH5XyxdBc6{bDkB0@QoT9tPuxz=i5Z@eJMQsWEbyzbdlWsXDvn98+ym>AK_q4G4ek2gL;l3Uti`wgs0D#euAB7aQ4Ab*PgRuC+v&2~T+K$%jm|>>PZcqLWZv zo_Q8lxqO1|t~B(%cAZwaQJN@jZ&EU_~uf^4KN+y5tbw9uqRpN2S znajT0r(N*YD6}atw93|SO5MGed;b;ee&3z;YWfL@?6{vS-H9qaF8!vS1&;EOI2&lU zm^qD>!fyq%=FO(6ip7^x86Yky0+-vjDp6@#eeZeFI-{qHrh7xP!fU4W;N@v1p>6AH z(e^pblj_RL%s1Ht+qz$lU{rPXIn*$T1zv98s>|^;5L8>~Qq4$Bo^iNg20~-w{W9zO z?RvNKZE|0hsZe&OLCKjnty}>x^5d2KC3jLM4k|i?(5%v@IdH*YZ##UutAF$xnCva+ zqUb)Q1Y1JC)10y=cl@>c?@`Ww5dTH#mO}a;;zttwpCdlge<8l6)jz~{{V&9q>oK`B zupoIS+)Omoa3orYW3drQT`M)T8e~DHG!;7*FBXKW$S^Rh*rlrGbh)V3~C+g#jzuV@B z^P2O?e#-H>AcqS>1)j5$oJKCuh4~~^y^rsBO)AQNqZ>8A@dZq>|4T+P7{X$NqY&j# zY7|UIk-;8GY}Ej%O*H1f_6>mqE)-(J8O{0&6=oBg|3VTtM2c)JMI=U1o6Hn5FC=Td zHI(45s~Iedns!&pfu$}ayFccLG7-4R2K+v`C-FJ0bhE|(X@aBDonsr=#7$}k{XgoLSu3_qNKNjh&9iQYa zj`0QIbc96--N1;+I&JAgI2W#mM&dF2sajU&R7{U3IRdJXAT}$K60m$O&Yf+Y27>3h zNLMPz7?38~n!cToy_l}gKqZ=@FJHAg$#iZypCKJ^$oT?l0ze)m*eH0@GwbbJ|te=LcRy^u(IT6$`cXhw2X> z{4rb+QB-u~bX?>Yl&HvwS5dB$va&B;Iqn+N7AqAMANq2Ws7ZC1LN19c5{hd5ZpRLg zmxY^tx^-ei5a_2ranuzLsTA5x(6+{E+|`weI{kc61Sx`O};2}ZuWbBI6`e#~g0B#OKy zY6SWT1&V8vgKk3M8s32YGh(fRfe96*ID<4#pbvLo%f+HpsBbzvbmQ1%DYKxo<#ICt zsqqXqn$-Pd)D9NtCN8>ubaZ57RszAe6oLk`=pC86gi|2^77R2YLcYnsTkI1*=oeFO zQQ-hOe?Ba^zBru7d%EeIPn?JdFBy|Pl4BeQxp)&6F}P1h>*%vgEUc(&jDLbp&;V7{ z-QVstCO40SiDkF%w{tn}Jon|uALgWzhRxBA2@kc+4%MZ5Wl}BMwskhkhEbm=V82{( z_}ZpVbYee^KUR8G=9(4f+FU$u8^xXwW!r$7|%~jV+-}doEB*jA0TBq|Xs1@}zENtm74-&Z&^=Nm(vr6 zKDz_F^jAEem3v+bL!A~CSdKdNdX2Vr5SGunEV=2c$4NkC6Do}BrJk>&>pAt z2ug4Hg$nZ&25RnDn~1rgrUcxi_8@Z1s=6L!FFl(b0V$2&7FeAh!`lJi(n~{;m522z zFv-2;CTUtLK6ouJZ`JPwRkpX=uQyF7;j$?$Qh+XZ4Jnc6O(p6ov=CK#c|Mr)aycV+ z&vSpAOxkRmf-El%N}{MsdgQj>_5F@JKJ(qot4Q`#pVjp3z#&rEX7peF$hsm*=WT^& zeCkhMQ-@~|LT!0ZixyvtuAvZ|28gSsdlh=ajz#atpuK7HJaV8Y{B;};w((+**K>5H z2lS!K|5Os1jT_$~x{NpBqV)&w5dZ zcxaKsMvTyJ7h16A(sYS14OS+Qul2-xb#|oHsyllH2~xwtFv|?crDg=3>+ZQ^(`4Wy ztN@mbKl2hviSY3UR*tZtZBwY;;yvL^f)x9Ed?L_0Im_(WQMWTlU#?PDke@do?E3;t*TU}n7B@mtAIJ-b^nmY{NLjTM#Uq1msgek5u2BYQBYV`%@V*bHioxmXh zECJMU_9$`^kNF~K<0$JB(`(TE^R$N@t>basLypaY+`1gIRv|t)F8Xv~CxJy8P2IRfa@S>n&u0 z$mrtcaNDU;Sw!2D+HIe$4x#H>tmmx`0^uxSsR0_?bvmZwU0G;pLW&B3{vm%JTUhC4 zz8=4*+h?Gk0z+@U8ARW}rvBmhrk6}s9$2JKMR zj7Ukr^cgXZTEn3lJE7RME#%R-O-*;NTZ; zR<+o2&|t}&_&O;S&I0gNF9Q;%V8uj|@JL^utSh0Rq)@7l;x~P6szP)>8|}QPHQBBz zQY1iO`99!Qkf*f%)qPX{YEcWkPE5$rLOUpms^#+Xa%Wsn8A&8QQWgOu#O~rDR8q0G zVt`F1IIB#i_7Ulecsf3koWPXC{dqLz1#`%9(nBoVNRmQA0b-#B^F2Ccr}fGt5z z*r+3lT&cipX#xGD*+Pwl82G}!Dy#M)MJ|!5fiUxRM)Rm$HpthkfOEF{nA$2et?a7; zy5zwfc16=@$9C60S+U^K-vRQ@wUlTzIR?wC%k1KI<<<3vtA5#FjH_KfPZ{YkM6(4L1WTds;o5rGb2RZ$ab3 zmL%`=`*x36*chvFFbH~sYXlGMzmbP^2n6OU+C+=7&z-9)opdd=R_`Y%oxJ7|(1*8L zZo3`dvxToOZN9p$3^zXR`Ojg<__8!=do>yy_FtuU0X-%+_VC_qyw=uQHr^68-WPa1 zi|%-dmupp-YE_*pmACj@s$O8rxdDEjMk{OuA7rXN4lC&3sySrd^OZ}(5D)5pCY7B5 zA}e+X7e8UJ*MiQxotyIotj;uTRbq2KuGddV8KJS#Lsl1RtCu^L-w(}pJ-r47`#}MH z7#{~I@O3mkPra{A{ABl=XT+C3YigHaOT8_ft0hL3I|+Z^vC;AG1}3$d=&GB^hT5I| z!#1Y!vHeg&;rzvZhd^{m@5d0|fTbBAUR)Vs(U6vJ2KGtB10Ayy_3de>CM0R7pbUvH zlB><>`*D$BM?gsvA+C1y(S|;H`vADCC-GbMTSc`r0}+*~puk*q2TNb2i_9d|mxV7U#R|)FxK9hyU$-ScZ&Uo~snHat#onj9tDu z7w%`uYsShTKLhEH#jV+YsarW}j;XwM`9kM@%AC3SG1B}FMp>LD(XY(oH4Dw5mhdW?K($z??b^BpGy7{{Ae(y0B{0aRR6iTcQ(x1Hp z{0PEtK`8?HiYQts*dIMj0j8VOZ609Q)fA=gJGj7R$tpblr=q%ZWoiO)nNVHX=p#P5ZT zh|9*pNKcfUQ2x1oKX{~xe8WOXY)-urS^Q@!?+5WJh>d88rawja+>9?ifcyy$;K4xG z)M$TC9VJy>$(V~v(MA4h&K65EsLop9+Ft8fDkF-4Ii869St4k(CVDx-{^=&10CLgS z`8W#l;~lpAnIJGR+@f7=hyxbhQut5?cN3HH!4)B$z^YV%xB8&bT4ZQD(yXK)9Kn+b zVkY8#9>)UdmIy-kH+I$>egHDkBQ1lr%b@F4r)k5khyCgB`UHE7$=I*-Ww>XGuj6qd z#p{fz_6N0*&&PX4!_msy^b3AY=e0Czl_vJeY!&yUDaQ{<{OFwDt?bNo&Tzw8uC$PJ zrkYF2PGmQxn|{aRs;BLgNDB4oa5DGBB&oB^YGk9wTf?O{nYCss+u~hxj^;N9*~(<) zf~l|Eapv%m0k)svUQS8vvx&bX5Pr+XRPP}F6{2>Z19M~{3| zeZeIB70?Mx7`N-`-(4Vg%JY8tTWB?L)aHo=jKswGo%Q{?K1}(k^PZG4OG=Lt09OXXW9CFZ;a*95>i)cY?3aqPl`;Uwr3i?p^~j= z_W+g6BUedQRe%;$N*}RGso^C!b9dUua76t?3`i3pT&4smB%F=I1Pi(Z0ca(abcyCP z|8jKFV#dL9eFxP|94XgB6HcRPFGM;^1m{RIQ^!MbUrRusmY_rX>2PPAmOxyc%khM7 zXB;7MC;)`4*saxOUtIj8jX*BhwNFN;dv$%#&cg6cxUxb1T3{HP3S7QI3^9P~9dYSy zVZ#Z(M&dy0aYDy`4r5oH-b`gJ{ur+pWRi}DpM)t&E;iES8cHwVL8^dzV%lsC7Om?P zbOSHx0HT2%G~{yyedEX9{qg|-s&4_x9tC5hf04Pyx!)%^dOV~uoD4*$dft#Jr-mVf z!bA-5l_>BKP1MLpLXkW^E+6lL1U(~;%z&>i9EYPej!(nJ-fkC!g)mG{LFF2HTMz~b z^ko-?q}=irDNc~Z9;)Z2C-9l4HlglmjNXVh!!mfNquxiv7!^M&)Qq*k+hc4cE?BiY zXuHsAR_)umSD4yClM$U3@MVsTT})dh3Gv>kYBPHlfITPfgvNplQf9?B^h)v}X3s{I zS8I(HfmjbLX-S`-E@_EB$zQAHd#^7&&OQrz!52)vWHB2FB{aS5o!`LTZd(5C5GRxep+cGl_#eGYA!OUQUf9qVsR39=k9@XB* ze%_dD>3Foj@1@r{n{4SJ3XOhQ-AKDqAd0k`LxT2MAKh(JG*cvREnrM%tz8N$soQI* zr@0$id*JmBk^O@1i~AD18|oeYtY_JSUu`YH(JV^L1HMt%r3VMaOT<-m)!BEK6{gfU z;iZsi`W}#%hZ{`Q)G+p-3=iI+OOF}ZpZ7`^1Uo)}IAe_fwKVOUjtDGH2NwjQ=G!Wc zub%>)m5J&=Kz)jbsBtmjmaLToP6s#(^ANFRfw{=Crd#<2uB@F(z3r!&0R#3vzWq*N_i-wL;0wh$jN}TRT1XC)Po@Ps7B1N-?HT)F+^J+wF9#n(UtIie~m|p^$VHSO5 z2JI%(?{ql(Opu=Ph;nNBJMJG zkGTKkN7%^kovOEwDv5T#B7@?+=&bk^b_a3V>G|S9Z&uC2Lq-68=X`%=2CoJX)GU2ODWluE%4ysbquqxlBQs>nZ`GljQ>_`zZ`)rn z^gG>U<3AZ(6-=*QZ=UI0@aH&Wfps%`P4(X1GO2RagUPc5u61t(k}**b6iwI9J$8qXRF2`VyhT) zxz91=PTUgM4dk2a>CF%2uD^L7-S-i^V{DgB=?jo7MQqh63nI6qsxbZ9x9B7qi%ae5 zk%w{2!X_M4R0>~6)_ovPuSykPG8TY2lSjG^OoOci(>VKrf2PKvBK|sAof_X&a%_BJ zLH(7b$w-2Zs!QupVcld@s;3pDMd4^PGiUj&9p~_)c9K~gZEGE#dWlI#&ps&{+ii%6 zw)JrB^EPcg-w@`tPac4I{2EtS_q~UB5%`4Q(X^*E`>V|DcQ%v6^uDL=SoXNN zbma@)5?R}wT$}jY+XG1-iX^mU83DM4EVP6!4%4-=#I%sRx#D2+jZ$Mp-zJDnQr4Sh zM~aQE@WA*({$Zq9W4W3bik|rl`w9H4UiAg3)fzT!9(eR*z&mIbm}Td96dTXGXf&I= z9EV;dKws15QB80Z7!mMy{m+(<1z@!usHOjE33$04gjgX3FKDjQdrSB_YO(?aL8*7hgKJGw0h`53>tzg@^7U>s^J6djIZDQvOp7BVi(R zb!)p!JvO}-nq$;%j=w<#3WnR9ra9~gZnqf|R_TSVuI zlm{;Lrl`lXMm3Q@s{K??u<&uPrQ^k8fof-tgp7nBnLG;xL{g~C7IH5y{WczTr-(c! zVN9d&re&8t*4H}cq$Pq#Qc`kgnK&j1q|B0P!j_goEEM%oV|2CsUOAgH`V?*cezdG0 zBRbhx0za(_jIQXE`@w6b9w})V*#+o60m>1Mw=ZfO{;}eWY#9%nq}h&=bxI?DbGpOSXuf*hVlj1_GP3_xDuyFlZP||* zH##C=5g%h88M5G~MBNe1M-_=fZ>6_5>?yNLjig!M`f7Av#HwRf3q?&tS(gWKjhjfNktamzI|deL?f5$NV<_!6@C=y0|* zi|umvx`)@z!mj0c{u+YM;dk~Wv}3*3cJ={Q`Uz&LML!Lvn~fd)h54i7c|r(q;|~4< zu_JpBwF%1#)pHvHI(=s?`nnl%+1YW~v$W#H*F4AKCd1|?lf%oT z7Ur- zUhM-&wzsctUt8(hS?TMj)_t(vW5&zDprXR3w(Hhq z=e*~AB74=7Tan9@_jlHoi%ZFYZk0n%7n|NL8n@5W(E*vFT83S~_X4FRCqA6wP_Cl^ z9I#osVSmqHxzpy{iRt5Vh9NS=VDPedpfkpRP9QT~_z}aDfI`xdV%({_pP^7B$E}$& z@P~mM9tzNfqsmHIC-0DhEvOMcgiTpZ|84o|Opq|t!3Zltfbcnc)Ii!Z zU^-QJ*g!W3`V604cJpA4xbhyc_bb3H!XFm0h46u>AaKe|TmnG_IZTGC{Fg#Pus-!( zp9%;wX7f!>`;4^`F}^Cvn;8dl)qTB*5<0;pGaU>1ukTC<_(%YWMvA;UIRwlQFL+2^ zV4#>Xc!U66X)3*X6l2onH{X9_&0J)i_Au-_ugAYkx#K^#k@|myiG@vU3~Y_XOl(aY z4V>*9|970YMpZ`|djmC>E4Lx0fiRu_P2&96E@&ej%{46_#O8+}23x;PR&XgehUo$3 z1;J~ygGn~pgi7%SSgS=do(i;C(?FM%fQFXahDb%S_=QkGZBjlyG1k32Sye@silTz* zxrE1U6wK`1LfM68hR4Z!Zm#F`^OmFU)mZy|(K^4)uF$M?zA4|$cKMY9^t0?5Rh0r2Cu12^yd0~#@JL)~4R{8#N|y7keWzh-if^tGusX-92w>Go70Csqk-wgSCq;=_>>3X3gL=ypaG*yYP^yAA0tY$1WjW<>0 z!Fd!__v&M>L7O0ktHuB$uVwu`E}D<4v@38DZ5z%;%`(iY^Z5&nXwwQMs!HP5IBszPG>AzuP@tyFFY= zyekj+=MNQOcso5_13s?hMjNszDsVZtww+hxc#~q8qQ#)_uWnEcGDGJhvzZO04O&O| z(XJ+(r7yA%OW@>OXKp|G>%87q=ps|B>CM;0+q(s|Dvr|agWC(~7<%+po~x9DAIVR-YL`t_0{QrZeTP)Vhcz9*<6-`RC~N_`ViXT3)N;5C$prFx&{3Nb?}!^<|ibsph;(3O%v4|*90 zze*}Q!Cu?nA`!F70Jc$XzavTw$~s&M9J(1EW(}KZHb0#vohC2Whh1Z0a4>Di_qrha zP)g|JcC%@0SGT7sKo2FAnUG`cz|xhw3G5~Ok zzp$5CXe`FNi_9GvES_BRrWwzbg8lUMTOQRR?twj)c1Lz# zH7kuh#i^*qr@Q7t;}CqpFSVhQoR6y(WNHU5LLo-^gB(gc*?G-CjQLq3e@fz(h=|r) zDjVb?Es0J)P_oa;u#`|}1=%&@M$#R##1ScR4nFE3x|4-~RVSVtOkt&1PT4I0BZ>VF z4Ix#IY{RuNpr{<`MDpbmsoKz!QuG^jcItj$;cJOq<^$l7SPrOn1$2v<`b~~olDgs2 zM-~q+44w*Ov&LB@tBZpibP)&?hKfj<-36zb?9mMfHtqxqRn?=Q9@`3(Lq_sjK-)2Y ztO%Zd&Htc59Mj+*{O+%}C@Yqzm8Q@@&KiUyhbRNA%YWlPV#_#o#*@!Do8 z<*uJ;b)tLG@!fUs?m789An+}JAIcxex9Rr!9bfm%F3d`AoQrWf0VFH2ktqr75YxVH zif;We{|h1=z}Ig!QtNdlI%lGJ+a}XGQ#gU#xA9tjU@&wz8(UCI8t1m7KhTt?-08|3 zooe|ijf2aVd6Dt)V7+m@(s6eIE=Rwbi}}L;F-xU*xrtcevcjok($JS2H&z7Pkuszk zoIYdgU)Q??3v@WI!x>{-3-F$0U%Jxyn4KEwHB-Up!P|4_IrLAT}=zGrb019Jkl6$8I28wXRZoYnV&873 z+kM0R!>9h0_x_c$fO&6%Vlmsxw^%#~T9CwdO3(d8tJK;zNUTfGn?lRojE`f7D?^vE zS>;*)7T52Nzy7{&LnDjGolhY3gF0nQUiL99Uq!3)mT+}Aroh`fJ_~9vHI_tZ`1Mds z)DxRtGh|b+VKny`Pr#?C37TBR^WM=ITUF>qJy1)+PbzFNKYE zBQhT{O(MsY(5oK?V#L1TKhh z8*n*#FYLG-t;Z(ica+I=w8m4A*O9APxJ$2P^ygVF5{uvL9M#HQruo$J7x=$B7_f*O z?Pq@{HIjcF%m&>5O$YN|X{j1jJ8A4QOusvCH`Y;B{U79fbV;C62Q_g*esu}S{Ca=| zZk5~gBMZ;6)N#0tYV1?F(R2sjS5L;!KGr= ztw-48($DsuyByfrw~J{e!xbuAraD`TiX`#c61HJqzOxz467n#CaEJ~7E+uS68(4~`CSk_@9MS zzemQhi-D=n9<~^X{Z%KSlAm3m4d+$po-=p|n$4Cu(QgL=h=6m>G&2{5hQX`uao>xW zi&HTuZtS87lY#7KC<;2NP8Oiv+2dsI}9 zZ9?R&pdg}eNTi_x>JlaL6MYZPUovo^Bm5$^J77YA=r>xRW#j;9z~n3mw*WLO;)#}( zgv=4u(pacjV$4wM+Fiy}3YiuE6aB=PBl;P>Z$RLS6#i`mWhlMmyvZOyE)4pa!ourg zp6|9;2O*!)?cL;4?eiFkBRL8xp+dHYxN>m3Y|Nc~U-%iSjOkfoOhI;U{}_oNR=I~r z0{9sI<7p^c#~K59%SB>t8PXz--YODco(`_Z)UWd4u#X;4k|5`u*W31@ySl-{7hFD@ z7$?Sywc)yKhWniHxH1ry#+H*=r0WDT1AqG`zxq=3Q%>yu;xXg}^Z}aOE6QC{q3!QM z&g@9T^gdOetEmaI*25_6X_|-^+CFEG6)@cE*SgcSwz|a|Azirw2M*ORS&o=}OhWiq za2$Jd-Rw3^g74!eh}d5+fc>om=m8K9q+G^Bj>fQb0y#nFTxv!H1_lg_{vvdPHq3+7 z`DadwW#doT%xiLUBp~M112bvb82!PFvjYmg13`<+?h>Xm$N>KJy90Zq{`I*7?uu~5 zVA!aoB_LsY0~we|3td#iHx)vZlwtU!a}dWaijntatODe2TZLu&O;(5k7_W^~XB}uO zUMOlR@!-fKiCA%nQiC{^GNZ?800T&0Zs^27Fy0a+hk!JKsEC&LqKH&DS{&jln_IbS?bkqaeV_=-b;o?9lY*3LPb=tyE&NqpPEvK z4TrG@t~5~-yMw(B`?d4OIwvw>(7m$Wc>rkk%n2AJx0D zL`I<0GM5na01e?|_&(V_^oeUlRNcXn9OMTTx%UHMrtakf!c2?wGX~se#k(MZ%R8r- zIn(!t42;WflanloDo%xJu);81oeh{(Gt;0R*%V45AOw~`>w&-5Lf_k>>xCZ-srB~r z)T`XISu}6%k=Mdau?3bZmoi>P@)dM7w{A;*;x+SrXj&1y)_V!<7&mGk=IX#(>-zFU z`p|rRzbvBEsQp|F&R@RQK?Hp6yGo zboV6pI2#s!&q~26B;M|V4!8E?cIctdj6<+QEXIsSOB z4l&2>i_b0wRb~+eb&J}gh#kpV8uo2HD&Yh%NJ%#DEoK%yh=lN@=ABye2s33O#zrlE zcf9bw)rC%2MpjM|Gi9z&!rW;uHE-QD z;9RlL`-b@>~}KOu~0@(q#uf=vN5=ny9J&x7A-Y1g(>fRFS3SH z?-A~&3OlyJXHM6{E6^K^i$2=OpLDo9T#QIiQOYYrN7Gn%QYYbGy0;fj)6sZvf0jn84y4T)gmu}hSd z8grOxo5NyIO;c}u96GXNBJT#O&7%n^y+iwwLJJ`10lng>NVYTG>spewK+iw&zW4W9 zu%&9@#(CjCV@chNer~Z9ZeTBh4$w#^xtMGvtw9kzw!GYF0jDAU^WM31`QS@06;v_|4p;;Uq0}Nnw~ZC z1m>Q`-DYe<&VGDfaDBj+lXxsKEFtU#D8xhp@q~ky{t_3PX#?Af<4dibD{O3@Vhg2qmYauH%19_f0uh_=%EdZav>2Q zbEe)9Q~6A$4ffsXP5BMO7R=p2PPvY0s|KER>DO+V`jaq7K2rJ>`wUB79^n;exkTsP z849Y159~4*0h~JmmGv!vOo}(R4)bg^gRR#DyArdl9;?KkPS=py2Y2{PLpw3lk@8nX_tQ3%>JJY(+4Uc@IygP%IPz)RJ%hRvAh<3NtFc zUlJYca~NC&k+RT4bG5OM-PEBj1unXfSBf|%+0uX0$mJTPIg?bZKWkMLybf z;!;8|ML+82OQeKAY%k>`5sLh;RfWY4E*x}YwRs>}=ukZlZJbTPEs3gLF_r%vCLLYB zRK`=r|HIfj1z8$5>z-xXw!7@=vTfV8ZQHhOv&&t!?OJ79)8E+_`^=mZF*{=2z3YFM zk$m!(;5xmMyH48(pztP0sy0Ax@KyeD4hOA_2k zM!}ou)_@xQD#-9Z{x!~(QsC`}uFviYD%i^s;5K8$8JBZ;Eh!OIIsJtf(h^k^q6&0L zVKnzZa;Zj5343ucRanPpnZW>|R5DD|g(@0oX3AelG*0RkM}dySaKo?P`yk zF4{to0+4RnVyR6svAu+mKCsfNYa(bi(O0koQ` zkV#R0wY0Ul;)3;P6^$)mFcFa0%lc!}$^{QB87@=W4-SX^aQc`%()_%6tM5#Dznh6E z)R%uA&c}-pcH&-l2`A4ta)q==Y}Sz zujYT)w{OqK$tUci{p>qFSQ~P!&u{r@y3CBdXQ$(&ZnrnXQ!{|1qrmGH<(&V?amYX) zHw9m>J*py33x^+0l2nZTr3fH%7qJ8Dx(kq9E+7;@WJf24G*Z{FLtpmTV?(XwV}!^{ zyjR}U`C5MmSNii_clXEhaitOWZ;a#N=LokkT|gL=0O%eN)wX;5oBTZFS#AKNZ+K$p zvk|OqiO{1}z(~T$8GhLeQZJ54`s83?K`oZdK%wcJR;jjtVbj5CAGxTu?6nw*G0)`3{cCKVp5ZtCi`7a`}1(Nv_cck_Z*$)5OLROiR zaMB1*7m65kM#GvwxiaK+iN=zG-%#&dbwK_co)C0X?%Ln)EC==ewYE9V`K~{s*~V2! zDXXkL7nSeQ+wYy{&dU$AZ@+(*uV=meSG{e?{+OQHuz;|ExI7*}BL3=jzP|Q18s@L8 z&3DdmCfeUH9i+wB;s5hxXzD$)9jD?wng5D2V?w`q%FH>hZy+__)X=ZkJ{RrrZs^xgH8ni--Up{OrvheR_|-{T%}gFu_rc*Mk`V3MxJZr-C>%j0E@D%Qu##EB`Dan#Z|8g+M;_cKc&{-ShqG2TJRFunwKpQ<9?JybeKw& z!!4ekNqE5_PMSmEI%yw(DMe99ffwmDcT9pc{jm}pN2{e*>YT|O@=TROLI8vAswBm- zM|$6Bm$AyjUz5Z9aeAdCR;;3YP5ve8jJ`D>x#;X42{kRF!{CT=Q72%?0nYZ{{P+$j z&vyakih3q9uybYfUU;`~v-HQ(#7Tz8*kSDtaUjx-z}nedC{KHxqBSG;a?y(H%D_Iv zrt5+%Ea-odKe4FgV8RJ?6_WBi*j5UCP&N7TwWL1 z<=3>cbOrcTTS$j2kLfW_elP--*C*1`(I4HKx-@~QSrPO2EQ|o0=nWgeyAXH_R8jiP zf|FQJRM{mVQ=Mz8^Qc9C2X##0jx|yf`Gy)aIb6c|GG=bq4^0om*uOg2DFoyH4euI( zpqk)8bmuc3a`%F1)Z8?hC*RM(gGy_5iSmPOPFIYV*S6RJ6DJ;EDCEdwjz&#c_mH75 z$>fO8vFIXlf}{F6*W6*8`3g~oym%l40#BQsVSq;Co$ABX{?Si`iV$~Pu6?B_0tj{v z(m&4}74e9qfNtqxPj&jS9+?(ZTgrzlyL}pq>~a7h10Ockcsm3d+ZC{=#59mKslIo{ zN-9ToG>wbUg7Ds##cZl-h0{SAA8Ah9$rRQ?SNn+3dvMGHMsj~A(tX0OTr6GtGry_f1y%j}wmtqZ!*ERMe^yZ45*=W^s~ zs+Z+onVRdWhieDd8f#o^Pk&8os@LuU*vxdfhZYJRWy|*A{XqWP82{5;$k*L@7utO$e>DOuf ziaT_nFstKH@y%;yz#RQT0tJs4kwp=G3z^G-MW+7|d0M*siMlDBw$LVu1dB-u>gXjL zc4~-X90|*rKkC(tkA{-gi0NY5USQ}`{{^p)hcZ6$n?`AyS zECqriC}~4Qv?U~@HW)|b{=s>f{1qcn2Z8|bu7%o#>E|)WCh%{hz(wHlhfbu?c`e!H z{0}gn-gmFk&Z6CNv$1vGuRwGjI1bnHD>_i7ser*I+$?q`0sAB%+QmO-uElVufv_>r z$#LiyI|ZEaC(*`rZ~ENrWoz?#p%C6##vC}#pgO=#JX2qez@PDE#TKqq%&fEZ$qJm4 zaIcTpTIiQ#eC1J5wJ_+ZTLwia1FXdSyJQ(eN=BP(VxR$KUqF@sNX`NnCwKR4y|Vl6#e%25?FFCaJjWSvcR|W- z$-d{-pME&cyv^L?Ip^2bI)Z$Pd_KGTnFDvgMiG?FT|XJ$0YiZTBY@(=$dDqfi6Pam zVh0X+p)YF=gt>Ara}5MOSyib+H-LNOaFZ=|QVpYE+;7~*lfUg2zXcObJ^I24z`fu@ z+V@K<1QdoRV;c*MlQN=T+9c-wYISd-Gl<SW7`n4WRuYOrEjk^dr(QzLnwushY{iL)S0Jo}*R?y<9wOi_oNrds?(}4s#U4l#0S4-7H_2fw! zg`Xx8<|cn($wMmr6-bEYARqe)lIK$koy549)FfV$E|pQ!QHo`~qojr@0x4pQhplON zMNLI|+kwogWnxpT`YyZLY{h+=BOYVH`D_HBba5qDyl&0f*TB_x_5JKK^A)v!Z_dSE zMQQrAs~9K2Pm=m9pdy?mT3LD_>Y7%;7oi+*fJQY)tq`JI096ho@)SROcJYH&&bWB} zn!b0=16FtEu1vy*LD0k&w=M)<;&83xGYVQUPGf!w>{Etu#O3Ljv)6n*o6BR(vu)hr z{l~mr|7#}I84;(hke1_7ghjvX8f7|V zWfxZiTN9HIzLaCDX26!>FKY((xI9^|{Z&^no!Vrb##n<8~CY zgJ4qsIF1S-{-8l=m@5f8x9#}}Z#$abrhosF?91?7`K|cH9;PZJY&V!GTZQ=l`YZ(-dlwhe|3fQurS%UzxSE-VorkSdBC%~M-B_AP z%3X%v*|_SeyxLlLejt&kDxSK%G++1o09(-*k%5aPi7pMeVo>y|C?K;hL<_B9fes0p zL;;)%iJ_t2s}Yt}RKFrCELISPb`IQ4_ffE7GriAc)<29R(^v}OU8&%~;U5bvpT+$9 zg4qA~#hdH>^VR*QoQ}hJ%f2&8oMuVh{vPiIIGk(hVhcTQS!mc5{Em?VU?kE*vd|`d zT)*P1iYS(mFRIZKS=W&tooF>DwI*xHo1eTQY)d5=X74=}qKXU2X z50=Mbx5%!b#c$_>P*zD{%TfwYdq0py$;e4?ME4R>6*7g@NLC}3Le)lAl=VP{bCCZQ zPD($Qe9$xxd)1W7C8$f522)c3{kjaZ#CLdS^}UY+yNnpX5-5eGoJhem(NxQsa4xyr`>wM5@$dzbxiMLy8r+gnCPOiR_W*x`v7P78k9!0wOB^ljTg&eyT zlt>i#MJ8jp$m%o~VfBFg7L5v5wi+d95yPVSg^Xxj#Cx!k(SYE@pGQzFgq*xCk=>b$ zC@M0vY#fzXYcpCt*TAPZ!ej=NE-sYh~Ro!A8;VhCDRORet3+enWiSLJ@ zgvwB~u8Sa$vH7=zIvNuzi_eAEeu@h26d-Z0(05u@OV)}`I<7cEOJr3o zM$W9`E4P_f)9oThb=D-GzCu#1C{*$h&NXlI%h9s_B#3ui3Y@_>SBS2jw5b^<6U_N& z=t*PjeN_}P%NK!_iuwS^Zj8N-6f&yL3TEr5gV1wz43sTa1AA3Fc4#^d4{@2*h-xi~ zQ*|gJ)F`B)Yp^J7NK69Qe+zPjl!~SMQKaDtuEsoaO(&-MP~HHSb?X=%M+pxMRE6E` z;|eG|6l7Ui!$N@xUAl*_Jkx4U^s~QJZp4rcZIRecZWu~UvdB^`XQe19$Im=r%j`BPwiw$DMX$qxF$Ixwfm6TS2QpBKB z7l~dF>J$`kl|V(Cj%jcJDFF`?a2qF}5@3dw-~t|pl`?I%zO&ax;nh%PqP>G?Z8enE z37JU9;7vO&lC?8o?>t}}Q$KUySUc?HSv%yH(h*+Fy}fccx!8eO`q|{u;?k?Ubz;r6KDQ%hfIT1^GciLhkC5*`AEDg zY{celdVe@^6pWe>$01de`_FLNuJu!ohqb)Y%qkhQmk_cEUX*rL zQAL_9>Uup?w*$GmN%q1Byv-K3CdL+m)^5D5PMX_E4C~f*{aM#R%mdMG$h#STv%x#3 z;Sr9TyVJD~WOYB_&=$|{EK}K|5-0j)|I0Mytok@`g$q(gT{fB?yc@;*7wvD0Fy`X3 zENDk_n>2U~W;zKG^h_#^#@H^KjOZ>JNL>tmi#tqOzhMNvflWB&=<7miW1ce6X-5z{ za=^R3_j;|f36ewb6gE#_yl7a)RcWLRPuZ`KE~qg~+IlAbA-L_>0eY?e5+T0R+e^>XD1dtG7ET+xV$#Yov}+~c{@tZ-AP z81|TJ5Xh|;UOsZf#tjTWn{ptbm5t?kC$d4722)q#CM^D7fm_LxB9x zP(FE_32RKmf<&Ldku#TgF&*g2G^AD9^g6Lh%80WpL*{lPPZ)7S4Jd2k#VQsR=Pc9P zsa_2WwMExS*kTSqq7F+^?XC!LKZnxm}>oi zv@E9AXA3zM2b}vjCVGw{?uS??N4{E&B?YNT*i0^1yi2TS3BQDnfxlHomjAiXztwWN zGU~-2<$cmSVMouS^3n53|F@Z5Y0sa2(jP#X|6`y(EnAOhAszqkJhcT%`#%q3yS-;O zetS2P@XB(v4CSMM7JxLOxtz)SVsUoGFE$yvPlT`0$2F-3Dfzlqchsk{QGnTJtox{4 zd~M5*&uv~NKUYRQp%#wRQ~8)NnXkW%CGD!qxu3HYZ9IhG-si%b_R!fry6F=k`;Khx zY%zm++;d%DPvOywZb)SN1gdfM7jE(k{Sd(5xWd|Xe}^?jFJK;48gn6>nJyFbtk`o& zRGgsjS#~K;Z!CnF>??OU63?pq3{7A_v}<8^tT@j6OE)yx8M80^Ay!yyO%!K%)y+<` zF9)FVo2z^n0yeeto9HDzcZ!aSlE-IZ*^pAH^pWB$=;$o`1Ihf|+U#RbL;9p0sS}xZOM-)0?|v+`tmr%kX-Tz5{kG+IM#U+5s}~0)KabI`y8L86BHk6mG~Z~M>(x#RIM0Vi3vus&(9=WnkMe09}xq5oX(1;aJE6NqQ^ zF}2(CWp|2*$KV(JXN3Qn@3F$$qS+g;$l;Vu|8qaB^D-BF-R_(hCo@~m`_ZNo;MAk? zu0+tQcen#sd=3W~S@Nwed{i8dnN21xA7<0}B7XjLs((5a>&wlK<`r@L`%y~j>EhyU z8FC@?gGATIq{7^&b1b{a_d&@<0Baby4R+V!c4gpnm4O)wPHRMcK`)+z?qadJTDVYG zFU;HBcjD4D_QXYUVO_odX(-9)(jDsFiYZRVInj8IPYG@mrS@?t z(`_L*-5jRS!bqO0Grdx5@s_bDO9WQ)0O^rpb;}LwxCpJpv)fC>O3YJOFRQ0|A8YGR zl0OJ>xUoMwa6HvH>sFLl1A59MOowNxyhot-r9*~}871k&sGd`t_C zt%p;%+%DoNz5fHTS2x+_Cd$T;CmL-O;E88LnJ#7*x0PHmgVV5+VPS7ALiC3|CSM&3 zb8k~WoCHP`E-HyAi4$sfAQDyP8OB`JcR}C{ri?ZgAARBVx<{zFa-@i-Yq#kmIb;n5 zvU;OhsO6Z}JJ? zxXlHv-5{7;?3DCQbkmki6b;p_!kXppLGaj zpPWbL=jw6)d*}7<-7-u>{WiWXf$0$LbiRw<)_B&-JYk_x`0dYs_n7`uuZcpS93%d# z*ItSKC-vHYapW|uWwclE|Gd1v2M$2k>`3k<0YyrV2Pd_dSdvWq+Akf(RKl2B8|ll~ zd`)$(C+_&z6Umh-s&xD%8)}83Qmv9flPZ-wQ6f26B5D0+9gJ>Xl}6R7G}hi`y6&p00$+i{&yt~1sM%8+T+};PC0O$bDqxWnp}9cWIFS%q0czaqMY+Ory8a? zXL1%fr1ip39D2zlY0737?m)9xPT5(k0#9Frn=ryQ$u?Vt>#D)(h?Ywuxz`#?3c^BR z_z|p2W!jbOvbig|SiJsSI9jH1nW$VhuUrCGpmIsT8cDPSis^L*?ORaL?no`OHUo|X%rfalTtzgOt_h73Q5FXkg76f<({F~A4ckTRW~eoWV-I8 z(4)BA1Feq5y8Xgsp=G4!UscfZ$bv;yDr6%JRW?`*YHTFLfFaxa79RxWd@f#T*3kE?lzvdzh6I&n6s~KP^?cw1U%v+l z(z(?~Hn~7VITBJ=d`vY5?snGD`zW# zMwS!{O+Z5)buQ=Li}*fXnqvv#CovPjMnB}k6edGj4yph`j}wLr5YbFyl)apC#&REc z8+*zBrTUzleOcL~RO;IIJ`nmsPrFK)3y)|UJK4TcgnN0>#k|@u&VG{2(Z(e1XuX0pdChb^xNS=slLr^W zNN1CF1TX?(5)#H?JxWgUNeVRGS&}p)|T zyBbUE7SjKI11SP!>8R>a${JV!7>Zl=B`_Rd@2uEO2@)=ex<9z}WiVKpJ2$j?Y4oYM zqi?|)AfC>vVVeXKB)J2WY$u9qHPTtUPTn>qP=zEOVTW@$DZIfNy6Flibnp}fv}4K? zCItPldM9)=@HT>64b)RG8fxuJsA#0C@fUayYvQ(kD_7X07eG}xx@jqiqyw)9SOo%} zBJ&ameKsthJVnvy@Dqa8ko#GEJW82kSmoZqXJb{xaE-}4RN^Ox>X=iZ!&KSP&`fP} zPpJJR6lU5gDgK1`Qs-MG1hLTY*5H7Lgm5h*B?NII+Yu+As=!1B!^z|VISHXQ&{c(i zVWFxDgV>xhn?6=c#aB`IG(-_qk7{dFzN;MKnKVJqUcI+G zH8=m^>Rw6a`7RbYw)S`Wd3z_Ur@T@(+RpI$JlC}>wCr@5Z>6QMj-$N#YrO5DD>nIn zU;EAd<+vDS`2|RS!leZ`7X1}Rd9uCe#kEP#>+%QQuE!aoVx>M6;Q+Tpy`=`F*dtRr zRX-!jT|F~;wLz~k#t`s~&0K6n{!Qng@rApo=?>g_Aq}gQ0NgsB-H^!G&Qt*NB7%FH zMdvYn37GUq$t=~fnfvXXdem-?uG*+r*9={F>SRDi-$u)erJ7y%GzeULE0t8XZq>%~ncrE5>0Q97uma_g1jl z`^8igRWD~q<8Ca5>0TPfBdm~mAp<6qs0Edi%mZ8&*SU%L6q74mu(VsnOQDp(`kGiH zil8QnjuQ$tyrh6Rq@Sn*&VSi#yDy>(XjQc36$3!FB%4Bk06#?hC!oPD9KS-1i&>I& zDRv!@NKA*csTo@n!m}7uBSA33jfn=5B)&tl7!zPhgQh~nsA@^5rx!?QeXTiOEmPeB ziT@WXzb$Ea#KO@WjQBP}*+iFCZS4EDB>=)e7igu0`(j87Cp+BF0y_i@s1Ffy|Kr?f z$!+}?%r`!B7+hG)?mNeKTNI0PIbs1=E?FL{UNIg9PtimO`SKoSutsS3XHpH4cY&6Y zIB|82`aP9ZAT13l83l!#plLbTAN^LL3pyV%4wF`oa+yIR%~pY-9DSIAD*S3gF-!@Y zyr2vRohwX7IUNs8x+DW-)Jh!QLgqobb-0{lThWPe{8BeQ;xW_^@j0?<12nQdz(U+!b9T|)I;&nkx9bIU*`T_EAJ8uD{o^f=`ByZE%3u$ZA?pV z-)9|2_TAjiV`}Lw%e*bif4vw>n@Aghw`^&3EOY@8FXZbR8i;DM-b+$x_>)hOO<`o4 zNVxsc4XIm%9UKihcbj#tv9Z|JUgmqx0!5FPpyD<}<8Fe&H&;7vd8UTP z$F13oJ7k-Amev|vO%!DGdVwd^+|$4Ct%ojFje6II)Zu=yT%gXeL)iwmLmtY)&V&x# z9cJ6m*?A7Dk*!C)P6e3*kL+qWQWp5Z6nwRR9ENfJ8CT=C&S|1e;&OfGGK}Q1P2Mqm zMo*jb)+Wgd>nGjBI19z$Nyhi$u9QYM5EXL~PEagYgtY+7lOhen8I96{eW|8UK0<7Y zpYFK?!|tT$i6f|S410X{JwA)(p%G~q=tnR&|0B<=jq0Yj5(2%FY=51Fv$N~ks1isk zwYWPBZMs&TG|KlGsY^aujh>feyi`nG4Fm1f3uE5N7Eb>H_wB$?Q5c(*g)xXVLAGQP z62b5(iXJip?T~Q(!JbCL8W%QEV{}-`^>dw&C&9m_w(X6!eu+bDH@a(7CMF3c1Ladh$g0GsHRCc1?8kl$umQv*XH+98xFTJp;3Axz3T%Je#E@PE2}h{d@Zya zOFn<9X27(+hfl`54Bnz^9rOHjRv#XEH}S?#ta+)en-TS75lrBrFeu zmTO~Z*IE%zlWb>H*c`k8Kfj|4MdM9DmeA2SL2rs1k~%y9EG2tUPpt(mpF2^~=W)wq z;%|(a*`O`LB-8T;NKp`3ywzpc#%-16Lp+=Zn9~uAXd90ShRMi&3jlUtk@6k7fnBG+ znb=VPiE>YnIbeh#D7zd8K+b%-rr4=KP#^om>z4A@+K9Jbu-Dz|^g~U?I!$D;6S8l( zm4kh49)D&AnmOLKsg5;x#;z%~Tl4PJ7l)AfT5Z_}+`Gq#>eXNEbH9cYl!G5u5-ACZ zjN2w+u-_~8T5Zs?a&gqHg5`lMX@BILR<8#e*&%^7{Y30RwU#d?Vbs1-)3HdPXvjit z+S-oX2V=!dOmPcTI%x#0UX>P&ObQXSIduf_DtLXyrpf#%W#?08B5=u+UA@r?P7Gw9 zYg#510t{U3G@U{|2 z1btoOzsJ`_R~=umwlL-w!=s4@pvV6{=7iGgarP|rto6w5A4px8pQ_C!9(WBYkF`w# z-MuUA8WCFRIr}I+c*0Zj8T>$FbA|lxD6s#Cz`dXw!Ht61-x_jJ)=h3Z+4#D9=jrzAXgytOMIN#< zzafK@bR#V*6;!E^Xh!`~t@@>E^UxPhgRvok-5IzJw%J=hD*&A zx$bVZ^S0%g|Br%T>Gb;&?DOsW8PTvWf{2fl*)QeU2PIT828&#yDLz5+vqUb-qzfli z(f1_TY+XmTq0A|$bTZN)(8h3Dq0L}@LaWebs3~?junB*X%Z5JXp~1!+U%}0YD}6fL z2{`RRqbbb+olmS8%Ox=?^_VN47K9KQR>Ibc^^>L0l3-}j1+{o-z8Y*l+_XD_gz*MSm5P1aW*%`=z#epEr5jliY zQ%HHO*dZ4`uNeiAOx}W-R*bhkx&(Qo@Ur_rT*+#rQ#2s;zr{D^g_73&-W>#Pk0B1=>(t7f@jHL2Fy)b>_j zqgWPzBoV_iz;a$`!Wp@br3rO6eLDu-GklpGB9H;wm0@Z@bSNZ%Qsm>}{*r=~lLF2( z1!cI!Ln8MOFdy%R={JaRLc}o&iSC8$kvgko>K_@ox5$`n_IBruzgn;j8|BFmnR79mGsopYwUQ zzTWl{*XwMP==1VgRqQt3>%=Na=v~tGPd8RD$;y?dpI?C6GPMj``_{T&9h~LS`ZL>t z&*^$N4?k999DCZ?NFmM3vjr*BI^2YvF`rfpjE=ik8i5sdvl$hDRFDqrPOrY($?I_% zAS^$}l%j*fenesuFn(2p_=~?nNB98$*732J3ghO@v9wER0*b=IN}JK}XS83WO>qpq z*oF^gJS!hMhIST6Y^5G7#%2~=j0-8osrX7hSh3lU5t(Luk>n6_^sgrnwLQKNFzcDc z({G2&BrU8~BG-=j5J4us+ZsJJnC2gPHT3e1J+gc;a)q!*cOvbAhTMQ&z*h0$dV!&& zBt)8r@*(RPg_E(f0#`x-q8T`P>pj4mq4#0INVciNZdvhHMq{1?NmJ%yTa6tYZkaYV5W05068~0 zHJ!`d=QY|ruZ53+DK6~vxzfMmhrJSZhMyv{MPW^E(&^aAXmt4bAEVpyTWcT5vik5D z{O{T?Z)Sc@dwEwQim%&geMK^ zoFOgOyzDUfJvZ$z7>+h`+hTH?xH?T)f!c#uEz;Hm+sj*q_U)WPO{pBr>gWxLz1gUZrK+X@9FHZjSS&(wa)P z(!)mk>6Y4%Kr>a08Ja+|E|zzrtmcGTCu1Rngv3SSxI-bjim)OgJWx*l@&ydiaTdY6 zvp^Tq!=dj8s*vM9)w=cgH&|fei6e+ZTFxSHPC!n0zoDXZv+#zHd<}nVaaW9B3z_45=W_T6K=WGz6zd988@XfPuSj-bo$-fmw@MzKTq%JaD zOxyraZWuc_9YPJ-?nPI1xfVC_`z{7|;h=IldjuVV2+X$aUs3R60D6@<@6H#kS(0`a zk7J`64nhB4x@zqz-Gd|Xh=m6%l?SD$zlPkm8+WEoq{~VxyY*Dua8h&oe7xEW>pbjc zZIpxRg2)j8O&>1nmrFVJ2vdE1!=q=+01DD&<#hS!EHi!CIdu7CH}?AHBlZj&uZATx z)^hI#dNzwPjaKz&27MaOq2SoCA#=>NW~5GB)g^0@YY`EI=}4ySg5JIaofEE{YB8)7`o1L7 zS%PEF1a$8(TCj-)u2&W9wpEEWpptJl%h3km-RJ7jIYyXwh;y#o#F-uYFJs;42K=hs zCGsiG`SsXfk*Buz(1p)z&y^F)MF2$_A#bvo&h*T)l+1*pNFG7eARipf1MEJ~RH!Kt zL?#mgQ{yn`a$-(O2T(lk80gruR)a3dMD9u8InT8z3pe{UL?M>|G4(BtR9t~uX|7lm z5FNx{v|@BC8Atuz=r@AJ^gN|r0!^3WY7z7Jdx4Z5!F>?K-wF@gR0H%Tsq%256s9}X z)041}$vtQDE8KCKpQ1Qm%HKwYV&Qr{k}%W~Ehn^9WNL9~y?!P&&}8@n-$h{@?ke=0 zJMOgWteB4C<~j7fM^qCL+TY`f{;@t9+qVT>AAYVcl$1XXGXiBsESK}?nQd+x*wn#CBus zEUIHLeC-~V{@%=x%LjywS8vlS(kj!8yO4GmS97ZR;NPRVMTM{wb$T3HfmWObpQvZN#~#{YsNWTKtDMIL&F z%cr;;)h~s2eP{y8!1o;mN(P|I>`Mc%dMIfykp0_+_Zn zWSj|$j?XyJKpSCpzJ+_>!19f*#XA3G)zw~z>;n`Ar(Tn@tnZH|mwa+3pd{}`Ah-J( z_BU1*RKs_Ak4&UsOcZ7LB>a4$h|v?LQs&i0Ux?9C=D&=2$7HAS0eKu*D?8!*X zYr!%V=5BA?G8ZzzcI>;AR_tkVW>)faVWx0O;Y94GGDWNd6V~nm#ub#zfgc4!#e6I& z6-WrtuRBZ`JnT`vUuHfr5FW&nVAndjbype^kMa<2>HYDn&8bK&wV4u0pD|SLDHY!- zYj5DTl+J(F5@C}kj51i68P7-x?H|BFr}rmOwkf`0|9cAepPNqaJOJ3`e>*L2{}-jV zlBuEb|AnxYswQiLDvq|JTYoh@Dwh+oa9hlXKo*OPFswj@VJiyS5{V;hnW$yu+PaQE zL~)_(?uH+cD#hy%QKfcX8>x_TghuQYfXe%PIw#D~R*_;MUyRk&EX;F))yKznM}*PmZ$-?+`IKv#@s6p=H9|4O zwKKuqZ;nZdSN=2RV%9vQJ9&@bQc8ohS`aN$M}cETesbnbKDLJ7yj?-UEkVM;LLCwK z-TnG-RbWV`oDB$TL3FLd0s5d`+kaYvHBq*OjlM%Zz1YD%x7I`dn~M9EnMM}VlMw;;e2y8g$RGp- zhgq&gU`FzNNttz97#;bR1GVogkZ=9XHfFx@DOROj&Y^(Og%S{19sL~<__-1iSs#6< z0L%)x8O?;v$MtHhQU~yQ>#wHiRp5^fC4L+>_(=-UEGmr=tKo^ zi3o&85|X(<0ucZ%r*2dx1DY-E@$sh8Ku>7o;!49~?IHqqE);i9=+o9UrIt!~&&^G} z0dLlfUNd83lV>KC*)vF6nQe1Iu0v^}p-ZTRr;JPpJgXZYx7nJ;s-1!GyopQRnwDH} zYF@99pJu4Lx=BjSUkfyJ{7EuzCTBX9-)`T7$2NL&WRd2w?<+W79gG)HNhR1>bdseR z^-pZ96c?aaeI`{*Wk0Q+-`0+#d0D6anD6WdlD}^@BM6I~Z^V2uWuq{ra-2S zEu3V@0F>2Ojl$ z8&#eb?cyuca+d%HqBR4K3-r;~3#gfEE^`;cGo){J$T0l8aKR%F}kPrX6nD7pD z{|<$RQWTzpCY_tU6&KVK2NDK^19_D?I+`RFId|QfxF->l6j)xUm7`4+lGcM-C#rw$ zBPgTZn0|j>H#teU{ecljOQ_4!<=}8^#l_A#YY&y^^wJ;I#phAS>`}*adxA)qZk$dt zNWS`bbvu#wLCT46F@hegxZGMObKc3jRMLyD za+2%pU2g9d#DP@3Sgj0>@A1_Da>^j+j^9afntl@ zT$A)-d#W3?V9UkS7Q~CKj;xmAhkKhIrJMum-T^B_!FBhgOuk)gWnzT26qN0di#dU1 zfOau-`?#uE_die0DG=F1HnH4&mdx`MFZZ`rK3x&GNeG9srqm#d;T63k4Dfk+&ORN2e zJIo8slLITeD3M}nN>&!C&txIeRFT6%mz=yB1^sSkitoY?RcU@Q-r2qu?aKZD7NW!6 z^PxBei!&Iel8E(}(QcggbJC~V$C!{Tw>4Symw_gm*}*tV`Cs0c?)qfSTbi7G?oO~uh@-+&j<*d z#@G;=>dmJn3fz_M>q~lQz9$t&iE2HXKrbLUDh?;9cw2_v4Hy&)dBsIx=VqxA9d4Ej z*Dpv*aiiyH@7<}MmlN+T~jC#*B>@q=GYXsRFQfWtfGpB_#4! z2&w`8V`t#Gb`i=%XNK~1v-N(U2kjKSO&A`whMhrdK5VDsc6+@x_^<51fs<{}YaH@S{Y zEHsb(*QGlq`KVdyWb4d(zdmK;3JhpzvKEVv9-?e2e??gw)E9Ci>o(F^0F3mS3AuEe zfYBXC<}QN{L>{UliHeGTIan@7{y(j~Wk8(Gk~RzkcXvr}cL?qh+$}(GcXxsX4G@C6 zySqbhcXxMpC*P1~bC#U!?sTh9oVxMDsn#1tHbD2)$NyMmm&cIfq^hzceufQz_VI0(FXgTjR5XdfjAut@-cYe^0 zLgjFI^aAC0avE<6VAAWP^Po}#uyVdfj313Cqvm0aXB5!6RE?8pSn{DFtq37)S|m^6 zmBGsT6<8vD8(>GXksmkc0?i(!E}4O~HZGZT>Y(~4OEnJ*0#tu$V46hbFj9X7WBp|4 zwG=LvTQOBk#ty>FEh&H?wE z?`<1*iI@BuzVN<+EFOM?Wf*1x)w{{Fts2C^3UR@_NR%zR+|4{qE;bwM_ipMJTen|W zw!K25+)&g*9*AZriDvqP5IjGT_gE>`n)H@Mx}07|_e5)G36GD9@_RM9u47 zx2b}CjUEX*X;BE%DlSve9LO$SdrI%?AJvafoquY=9Z{T>1a$DN8ai}D9z|BgFeyUM zEJ%shfkTM8kzuxG^Tke)O=SCTOxtLI#$o7Dh-C|vzC}tJAf@RS#y|z%s7NWrFn?1a zR^&ClSwz`~k~%kDH(LZ1SOiO+_;fyxjqnW216z-Vm(>pF-)%Gcw$Q{X(ZZ_+DCXYk zSN*`{^5w=fKO-0Us9h*J0A=2BjEIq>u!|`qBsC~DDs%+F!aS%ocY{8Io}gi;Cv@c{ z4wSK}`GJ_ZT!H^-aQ5(dA_6PEVuy2=++#hS zWn!V@nv&l-N@=Dlw(c9U>4#4Wy9(TE2tqaEobgA4X3#-O@KOhRLxUn6t~VP-Ofp6V z&QB8eH-n`;@yps{q2Ldxmz=;$w5Td;ak-{JGdiX9YOV0OrjN|!E&I79_asvf{#e;W zpB0P=-JV~Cf+k0nYM2C=L%{;h7W#zx_qNjC6NG=x6jIeRUzkSQE)^Q6AmMcMOrj=P znhnzBbkwI)^D8pL;8a`ezM52|Wt_pKQ3yf~o}|sEM!)_R0DKW#L{kwky8yJ{CdLFC z_ds72euO{G-CqV~B!tgu|6?=&m=Si+ZIL@RWAH=81{|gcUqI{*?x#=77al9kiZXTz zUhO@pCHGdJ$N1#5t&?mZ0;!y#*U(dy5a>DlkjwtKtqBi=oFBSsp>E# zTp!4dWYrPMK}{x3zBaEhsxh-?u&9 zosJ=gjC|%EdvdmSsp=H`AsUn%a`(l(4?zXw>U)*5+nQrOa0y7%<@v@B+q6kY6@Kyk zdGq#Snmobr6AI{!G~&7|Gl9^cvQJP_fgGA92)mNmPX)uu^Bb<&<0vaLF15*N-HCY# zdU?oGPbq^j^O7=#z6V5;`sd>$f;}X|ZIeMF;#HRUx-O)`^2#>U$b6%`$UlBOK~@2? z_yB*Lr4b{F99p{@CcoJ%Iuk^jgG;I{92PL?>^JB{A?c9fj_6D&Jo9BkCcDL9tydx6Nt5kSk)x8}q3Upl2ql(r)pA#~#si-@vmzXENp%Wcy|V%cK1*!Pl7pMXG>H$j~%vByfHu zSt{76caSbLli)dcIL-xg+aOA-3#@I_xNmcl8R{W+)~WeA+>MgR94g9t!yfhxVW#HT8i9>q94g;9RN&Y+rhRgZ${+4> zhTm?z*w#eT$xaYEc(2KI0G*c}*wrUcL7ijE5h093C4-*r5JJqZEu+Us?kdmH)$|j) z*F>^vkMH?MU}3(lnt~o1zocQ#{V@3XkX~3*@Y_-b1}98}57ZJR6xSiS)}Is~JIFGy zpf(d}b*|9xuDX(M(<0mWcs#-+fFSoQ<{;E_su^$u*Gpq1G=Sd3sx8Q+dnB7|P}u${4E{ z*qS)lG!8MJ2j}Y;%UumJW*Ow_7^-aq)tl4NIIF5I(OTDbQ|5AaEzZtG+^2x#t|r-DL=Jz#APHk? zL{DvtZ09I^=6}3KXJmHy$WOYQvid|;XB@lNuN*OoON zU7Pciu&R!iG*%sZ&MJWpZ{kiznb<9EiF&YCH5R%k@`k^qU2~tk%e&6Jn8Cz-->se;p_sEIQk?Xy7-$iqST2SgrW2a-3R*+KQ$^#(VY zjUf!~AlQYeLOUbnQl5ztSM(3`x+OLFtF$B-p`7HgN`P{-n%s zPNA-AlJNZ9bZHxWci<8oEN1b_>SbOH4LC0d2#7E^M&s9u|M~&&set~=fwFrw{l)_P zne;O3Q5E2ACSV{S#Mgf&Pfu%VXJx7MRSL$357GaG|0n;o`9MjIPEW~LwpnimD-KKp z&H==i+Z7xYCY(5|xSjag6-zE^`FPH#&E1HDJzLV2_o@8i_ms=N5UhZ5VZLR@onASw_ZWLkiP!u%hIej4hLId2T-rJ{f@U-wO;c z8m3qxI@SO*qy0-mhF=;Qg#EXMTn80}YIEg@La`<9LJ0JzcP)zwQP)3Pv33Q68*ZO+ zT`+dIZfVztQl>j1gt@V&$asFQOJPFEqa^>j>s{+e6@bYK_^ zWpcPC6^9w1j5u?UTU_OgE;e}N%cCGRV&L1aG>;`pDO(xnB77n4acwv;vc|n1aon;& z!}C0XzSCc|RFk5NYs+9)i{Q?QWopUwwl+26Q7301ZJmMrspC4DR#hWFN0z^IWc=$> ze98x~%jU_$fxqo{bFFgW}t=38EcerDH2SLfh4i zKX${WlfN}{hOjACfJRm0dlW_-Wfc#uBVTyWh!EFUCP^8{saAvbmg8IUo7Ntm8C48l zpsUBXjzMo&^l34Y+q(>xk8khU4DZED~uoKy~4f@u*ADP zI5xXrtH_KOxymjvA>x?8aOa*NTB-u+@65XA8&B@SzVLTvF*{7_CPSH68%0Q7w~}r) zgpX&B%6y*px;dO~={ndO*eNJmQ`7XRd=R4% zlcJR)7b_;0`!FylMx*>;P(qIU1C96qUBA+)jY996@gVkir67{{JC`PnkJ{NaL43_k z8n2FXec6*CBm<^K17P3+o-cjJ3zmtIrIoEgNH|@dq_pJlhtTkuovm>I)6^FD)j!t* zGJ@ikw7OvdXn+g6|MVS#R+jd*R_0Px`VQs>fXtsq>O)}*y+|!Mr<|yUh_EJhsB(I4 z?*;Vr9g7?cb;7}1X86(|1rVmrFXQyh=AyVicJwUT4zRFq0$(F;#uTm)%X}}`wu+oQ znN2h=tLRF<-znVR-`90Odq5oat^cYUczXkbbVeb0mlU9Y@9I<7=p-bML+f*g9Oq-5 z5(_rrLY@7Y!s>WA{d3Alln)bDcf@7tz+U^#Qlk?PuL6&6pO7D6iX?T{Q+%@W=+5fu zY9AqF%xa%(*I0Xq0cg43mvuxf8%*jADy;T+5m9(dKk8M}VXiPum?ZXT6de1XoM@G* z&cHn-6FvbmT$@NCclz2xL@xOaY)Q&&26CM-@O~)!U6Jxg-DG&X-V2rs){k|Qqxw$= z!n<7{cOj*!Nhq=jb{TE7U+8pENmOk*1xWJB;ebPkZ1hf3QKi5~iFA2SKTr;EP|Q*u zn?vNB*$%tST5e0A^SWxk!S(6!k!T~-?!1EYWU?9Ek!1)6hA!7O>(A+(iF!B)$)LvrOcE)ON zeb4(0xk6Wm-y=E75|sw%z}hUluORqz^E?w3cF1X+?{0*{ zQt6r7QaBsO*6nLKtsM$O*26mIsY{d8g^}i&R&;h-_my0MR^i2Mz#@|%yVmc|f@$z1 zz)j+bY+P77+;QRD+|BfV8?A*kDB|nP`-pFfpwdNdvMkhV1fpV=7p)KF?NuKOo3K{E z7p+6mx20J+EC17g&A)cBQD7UyGt>;>Y7P4HfwKRz#yMX>R1uw>7oRj3qt}zn>Lbp$ zLcs;wU>pIH#UiE8 z@GS0D&6X8Omsvi#TsQd{ul7Dtj$=l@f~(Da6T_7rP&xr_J-`lym& z1cl(I*eBVq(s=jbJf8=uwEUy4t|9S0oNLcjPaeCqQ&bJ8%D>OB(D>0G%J`IVC43Cd zQEP{9 zsa;pe*@5TVH7xG0QgUi2rp$tAP)~hHhrl)1hGb<@^y)qq5S2*g7;!AxjUCKBM=%P& zg=V(nX>U#`inP%VmFYRXEqvyX%{k*jwHtpF0WaX?Aen2$LC*q$qRV7}1%0u@sY*o7 znAXn-q3uMC7O00t6VAQ9Oo%OhA)GqHLS{kw+>?<-`J&^@N?>CKnzT#&hrBce<{Nn9K9>u+^XnyakpZ54`a z>Z`;X$7gx4o(+%`!rzW5$M8h7FoVj!oih%8ZidVj2bO(&oKI)Q&m5wW%|U23aUQKn zW+aA}`7U8WuRiovskWNg9XEYh#MnkQ=!ew7SJ5gf8kF5k9GN``vPc*MZg*u!;^>|T zEP6H{gHKs{;EH@p=yqDppt_RLK&`WTcB!o*JXfR(J2QR$nx8(|CD5^nqaFor1ITVn zX?wf+5}=G>sxidlaqpe!^o1zWnSnRbJ=4B7H@5S4%+Gd3K7VF7G`I$CzJH99(%m&U zIK?bt94-y%nw;d0Kf^DRFvymLQ;Qa4lVFl#UfTwZNI++bG-<*YFa=j7T9?o<9Y#|q z|72uBCg0&jtIcGj(eRz3yz11G7!4XEisB*f2hj zz}G<9O01VE!z3Wj*f_jcsrr(=uGg8^u$bCZmw@fSH9pCF!ah1L*dM~UVMam9!N35g zv7S(_){T3Yc?k6e?Sxf{kH0*lP>~amCAO$V=Q8A2@c4miQl3U~Bgqwu>NtmRfJ5yHe<9tF^2MWLDs=1&5HEmO(%FkODC`}_8$x8hC#;nM!wfI z4M&yAA&rbCx5VmN!pK$-**o<94%wLR&r|6;LA)|Eq=jonLP=;Wk!hWIzD#`nuA2K~ zAjzHY;Tdx{aE-P}#_x75Ehx`}3cqr2kZgVI375+-W*6mbcCJ^CMFMa_ewB8t1bYg| z4F}Eamae@t2-Skt;4m$kgAZsuIV~HSAaEY^e3OxlLu9X%NY<|ykPPs-qbhN95&y_|aGGhcT1{9nC@>#oe@DyQO z9uW*U(N-Q&KF*pK_cL+oK(_2mZ^ zIKh`NVoNAbQUz4IpJAaI(@pDIay8+^eyHs^ctBIPqCua;mFit*pt@mug?09f3(mma z)EgujQF{>h@9E(p!i1q959EV}%agzh3A2T?PlwRM+>c}p6j3cX%==EoFe~RklxoB>V@;U$&ZW&>-esWLL8w?+;t8Yo61-0ixpkNUywD!p*zj_~^7YgYue zle*S04{nyPK+9&LUZjVYP4MR-$9?xrp+Y9Y4T>@Virm382xl zh9C~HB4%@htnnAGHEtqh0>kx*@a>anLnIhgYdz(h(y^I!1Xh~8Xn$9*<<83X)K#NN za0+FwlMAxA=imNA&lnE4kTRH3b)G?p5`9Ypl!%lr)PW{Dt=z!I!6llFcAczwHvpER zXqBVm8r^qxN64h|%dlyObNlc$zs44BQ|mp@mqQ-xJG-wO*r@vT-t!!zX?ig=b#*Ol z6-|;*BagegWFAJN1grI`5~cMIqns$s1^)e6^MVvnDn!dkhd41=Nkn3Ta096h%0U>T zgmQAl6ox8Et<&4p4ZNQlG*&9uX=fA}n~mhDG-q~4LBC1vL7>%*g2ZrU@?PzR3sXL` zhBgbfVl$PYthn?&@M-sKpFn&e4O zwSu=zE`5<3JeuSXt7*5vkl+ZCSI@2|nIYjqIWdN9kovSTeGBGZ@kV~`vO)GsSzS{g z)&6rGBTo`z;yFa2u?3V?2O+lP$;F(>sMSidWFUo{kSvIhz%gb`kH=7=&@weY@<4u}Jj_a5p^c&p;XMyb24 zTxwUV8QV{FYCbxHCr*>)N~V zIVJ%;Ku})W>uRgedn@0MzQR|xaQd2*o(j)F6F$zhO|+9%0oEqVlp%Z)bEw!t*MrRB z0QBc6YGz*E82&?7$mBeeRXkBfb0+;dLe)HZ@KMYL8Ev~|AT3@fVS;YdMPl3j#!kdx zf6%~fyU!&Tn2r1wf=1`o6ihbpp}wQYr31-$<$8F74lpLZdogwi&WSjFI55*`R~qa@ z+jZRgVeD(xenjdVl2j1n1+f|)iP6X=Nn#gUYL%@BL8Kj%05_SO>C2P5JRcrv8R%*u zyN_w2E#S$^8R^f+)(Y=f8unp&hSoz_;JN!ak1{Sm}@NCM|JH)$GD!Gt#ta z?0ya@bkUx)&1_f`v&eK>ur*uTx4SdgsGQU@Hf5Wa!FiTw3F#4%+Wo9#=kI)I;*!*# z1kbd($N4;a>VExdS^H&OD_hvK4+IA=RV-dQ(r5g)KfPsn z74h{;C1!MC(%>|iP;p|Wt!re|_(zeq%A!AnX3e__|9m*UB~+;J&4xlAXlmS94$l?c ztSwQQM;`i#<+sg}?=oUAWI$5K&Wn;l8Bnree?aq{*1j!~#G-ki0v zTDDl4U#a)4yjY=_i#bQpeIY5i0NW5iyfb(3Xu zwFX(wt@;JHHZ>T5WXs!u)IhxZm=$nAm{Yr{s^GbumD~|4Iw%*R6y%PNa)Qj~sTey( zeyeL}R(_kE%Qov789@MAarD)twPF~*dh=*w`o?$(9*I^`YBJw?s;ie*R@o1WKqdUI zUOT-+%=z&KoyTsL@&;7OGtB*4DL%Y&wiNojv)R`zhG4isMb?nn<*Sc)Ru7EoOm3oz z8?aDK;c$0P+e2pKfiXt-wr70FoYM}2d-KAYfYf{o(e%;SStOyQNo`&W0vty(DhH*L zIqx_M)Uk)naur0A1;Vgr{1K|5i&WRtdq*v)S;-FgBRA{p*)zMNoRO80DUuDe;~yZW z91vUN@W=@2mMJ@Bvl=$0+#r-Gir-C0YxHWuHcCutD3rDZHT6Bq7Kto>!&&dgW3m(D z$=r^P%<+CCQJ3ucX&bEZ@xv#X%+^4}F!%7hozt~U%wjsWw~tN)&7D_T;A6VY=WH+n zZr_gZ@vOrZ6^!xMGk=)+jYSbWfMa!ob_~ixDjND)#PI6FRP}&u%HptP=9$HfbC;%p zA2jVrY^vkl>dFonTgd}9EAkU+w*nm#>wOO`7z{R~-be={Bw;@)j$?SFK;o-{NL)KI z2u8Sj&5@f-z^8@*%kP&{E28HTD8iuGQzsyXli=w5zz&8hUmzHe?(G?GR#vY;D=xDR zkE+fiZAaNxytQvVog*fUd;cld0gD_a%KeRe3ehm4xK5RvIEUqlVyJOq`D8&1mYrE6njEp>#KSDa$p{*2Y?MY-LxJ30>4Cx}~u!t-t3&}tyvvJX%THS6ba zk|XbRl1U*(YBoGsL-)3yl%#W0Re1a=5#p|E79}{%c1O-M3G9rP{lMFoQ0rY5cm2u3 z&urL2`fA4d&PaZ^wp2sDgJU*X^6s!P#0Kx2CjJiY_}mI=W9}ln+eS~Vn7odm4Ba~Q zP7gSomj`MD?oI(=@AevV(ue*O>YS$gk`KOqb}=fjE?ggU;?WP!43s60c$=mD@kCkb|Y&eFmll-G-nIeKT%qyChAdB=K32LEDTYv#mESo%Lv2 z7ZsGf9}!&K>85{MEb>nExGov7nGAE1*~O?HN-8=+7uLmIjRZ4dZRXfC5eu{Yxf<+- z3l@QN{3?rRm{9Ei@*Ja-YYA0d+bpMZA+=DhxV`QI`)SzuY?>ks4h=HB~u;2Ycm++Zj!kO1h4>lF#3)(&1?--B>@ z)f~2>VgoToahv*dW*Cf&6`YZg*RFRuAo&2{L(QTbU@K|s%|F{8ztgK&yHrOGnRd`4 z`kjEcc6`UB{PM__bBXxjT^9SLHaQe`l)R5_X?|NwQndT&f=bOH>$h+6;_ly##s=L5 zu-onAC>i;jix}HMk*VcK$v0rDPYf(63p0xNjlX4|+fpwhuV#%O`K#(@Qwats3#Cz8w|&O;p}6Yvvg-%nHegyjbpL8#3d z_TFvLP(UGqi@{*9#6gMPW^_==3#!gD=|U~6B9h$+yv0TC~x`88Wm4MF3SBEC}Tb15w z!h5=HgMJLutmvi*3-E2rI_0>Wcn|vVs?%X&2^>9(?=9Tav4~z71ob8xo3yv*nVqu2E4j7_O zKVxe2_aDu(ZGQ~b38#`ehSU%rJ__(+5Lu44Kys;$uCifW3~!C7vI=`|Vq3*B`)SDd z)8$}@p;lWLZDhZf?uu=RDhSkoYC;+Fw&aa;-A>LOEtQw)hsuPSOu)t}WDB;z7mtcZ z$++ddk9~N9n4zJHvd6Ts&5_Gu?bF&fbSO2MCIXKus1D)dew4x?z_6bKMi>nymZ4U2 zS57+6#Et3I!9826McPjBa7=SX)!rsp$y0(mHpU8CiI}zYp%LBY-f@TY1|aglg2Sy? z6JhS9u82qA&?>CvSAOr~R>&@XPE$FlO+gW!fTO}A*rOSmb#KD-+qy;s)nq^dGqS|0 z$~wUfxNOZ>u_zW6Yc*8^WAuNQ1V8@(x}a_UA*4c%>z$U;@JHt6=m*nqT~ZYjh+~N^ zds2lV4iH3vE_V}E|LC;P@HvSOx2{8|5<~n&&oY7VW%W};r+)pb;V|+N4f4w|rN-6% zYGTf6X~ieCikpYIUF-h|iwayU26wEN;9B2P+_& z?gR&h9^RlbiQO!5aR#!iJhFR@V^!%SzvkS9N?2_tTixc*{a!tb5fT2QEK3y~Eu$7( zTzF*|W-U{dz1OXsqhc2p3p1oQ2@bl%J0;Y^+C_TeC5dQyG%5U;Bu>j+oBXL;b8k|0 zS^j)05@YuXCHobMl68#fwDX(g@B7&qs=CAqw%Ms#HjbQ{Hf^C6iMB_91r6#7O(4l1 zQz7P2eX6NW^bW-*k}>FdD=Cu7F|#e=XmhQ_Tue`K+PhF`ozqwd4Sl9Q9-D=U#rqKQ?ch?5C zA2r(#4jngXtnRIrw-C?Gq43k5Do5{9CMIch!S2C9YB*<;AU@UX8$e#d6KcWU^^;~; z&JR8yF6(`j=sq_jAmL)|%>SxE5#2l17-49rR+B#BLTx^$qHRXBK~v2qD4%5yIJgag zfjQsCkaR;N#Eb+#wycxRKGbFer81OsmM25oZqZ-Xe&`DuayH?@xG`5k6X*A$#K~NZ zdZb+Y_s*S(tn&!&uExS+gwCDif;83IL@JFAI!=yk2Hq`#)K<*-$i~hMyQ>eE@1932 zOuaqtk7Y9u-9@_%N_gi)TRCg9M8wK9WGFJXrRnH2Kq zvQEu!xtJ9>n-@8|#*mdloYHA9-<{w0$e?ic$6QZY1Xx+b5{%@IZlZkJ)q;ma2m#(< zsfjo~8u+wN4oGZ4x6`)1d)x-s>qLnWY8QxMQl%YW4}xZC$G2&VQ^(hE3R=d7`IDV~VI?jTZvh8d>NHg6A7Xw=m~TU-b1C{B(GK1>fKN!{~Qa&p0Lp%XZHjU=nw*)USD zsFdre-kE5*GDtWl>cHQfH6QJ#z~n<^(3@sksgB%fbxkz?=LwzL7=<>eGoiQ7XPhZG z!LJVpN&+ysmu;UH!2n}w4}f{}(!Kv?Z(F-0E2&#e(4%^mPX4r>gr>RrCYo?bir8tW-FNv=*R%RTn?xJ(rjta?alu zaI^+NaEJIlEIR~uAU=FIFC)7I_#*x9p@{r9?-FStf))C(uG+~xIQ*}$dY1i#O}>Wn z3lb{YA}9HtmcXVDZ&Q*|b}PL; zQD2DVe9-P3RwIP*i`g7(DuzJ+ym=A+m2uliyRP zhB~ynoMBQHvKxe#0Zxe-J#U9YvF5(1vI147x$lGWz;xW<5x5;msMlzKZ3U{yzJw8^ zLc)kJ;}k|Zp^22%^d*Chq(Nwh}5P;s74f1K*#=J16TzF@AGlW zRcC#l1JIqw(MQ+0SaBq4=$PP3XwRP@IyZieiX7gydyY9t74n~M>RJYb-fSTs&+pH7 z^x&^FZc87}i!aPPvc*v3#=ChuTJYof$u2<_OO3b-vFdMqSZDShR2z5$OgHWM#jn9s z$xR6Q1AeW~Ex8XOkx~#vKt>Z~?o{QSbdM|N5uPIldbMh_YeNx z+b7FhT#96!SJ>dmZO=c|uC zl{DZ%pI3PPj4X6LX*X_>`rcR6X6jpf<_tTaUZ?tF>*#7Jx69TB3`0xj&J*}%fKDYO zajUso5uPEv>w*D-hUf-YIy4{~B$jU~n;@c7;9}4gq*uXaJ3uuan0si1bYi{;9<~6| zGJO%&GL~ZIGmd%mC1Rc=SF+>*h0-wZe>3%7)_)~2A!nEYH_6NM7gJA|Ur0)r)1-vN%3_|eh zp8J~y2nj?CfOvT@{8IqYI@;+01kgbPz5VqYU?5B&E5OdN+4v$b+B+*j|II~U{GU@W zFZF*U141)^Ac0;k`&+RdN+`g0t^fq?q5MnoiyXkIVE>YAuH#^-XZ(ws^%~^2&>t-C zj{^Pz`G>{$5)uNC^sfi(zd&9te-i+$rDg4+r=w?Vpr!SP3Zjb0$1hTUmy-$#Klf(- zRN-Zn$f{^>Hwgg<2osRRhvHw--vXq=15)0-mi{NTNzdBF-q^~LhKYd@u&AW>S~0~c zI6N-^boYG1ukwF`zE(IUmjplu6aG8&8DPok{{{M2DCBEsBUSi6pyz`A*17>e+iH}2 zU;Ygi20e)InV{VVk5HMAmH-uM@Q*1Y!c6EEfO`s+6HlG?Hp?vBR=5(vl^ zpgGka23M#36Z+pbLH-Yuql6q~djdcO--zeFl)vk*!CxXxue|J3hymbpfPI!f!SUYz zzk>f|bh`$nUrmc5J*me70A8B>x6%L3oBsPgd}(UI0QTH#qyO)k|7CWL!%+U0#`f>! zYt8@8DEb=w=)<4jzbC-|KbZT^%#p7_T_ye#=>KYF|Nei>{M+oWHT^rA;A_wYss9A} zmwBCi%ny73<~9G+{j2<4e+~N5u>MR3_8K%?{!h@~)7<~Vyw3c}zZ%qUO_yo^ssB6Z zpJ&>?2IbZMPoRI9*HG~v5NWLBU#<9e{Wa*{2S2_B4Y&ADpnsVc*Ej4}>;3ogHRww( z{pYzMuR*)q{sjH)X8s@M&3*W@ll-4F{W~oEHRyu(e**o>yvZ*Am^bO$Kj!^if34}? z0g$gjw*&qY=wIgjJ^0lNM*cAv{sa2Y5c_waz;Do(A(koPU!J?&ALA=jY0w`of{Xg} z|L7w88C>vN^h>Ar@wMpx3N?5=p!`eX%KEQ1@~`?Y_vXKykUz=&uMxPuzq)w&f2jWd z6GHw_4fu<1@V_DaNv(g4@HX-_!aqd&*9e`(|5)WeBK*N=e~kbY^%~(HhWr170J2v5 zONZZxuiMp&tM(_O{58T60O3!5`D^*Vy)@l%%U`abtj^!A;L9EQCyV^Ge6{G;obuPQ z83LgH@K%1;Utjy*B=gs@V*s-M^vkSWUL$m~*N?q8BlMuznSp@vw>Ga4Uaa#c|Kc^m zQpW$p#(0hKTV@ipGuGdJ>_6@F#X^7XN4~~j$otps)Q}?de`{m+!p8vkZ z{$9a;Em*Y#2>XWJGu| z$nW}VqyN1^@LG0R^}iZDz`1@cyVH#R*JL3`GU(TJ|0(;Q9_{Z5)oa;HbpYAFr&q5< h{}%O@*i&L- Date: Fri, 15 May 2026 22:26:51 +0200 Subject: [PATCH 10/11] Normalizacion y reactores --- .vscode/launch.json | 33 +++++++++ ControlModule.py | 159 ++++++++++++++++++++++++++-------------- Reactor.py | 70 +++++++++++------- Reactors/Reactor_1.json | 24 ++++++ Reactors/Reactor_2.json | 24 ++++++ main.py | 2 +- 6 files changed, 230 insertions(+), 82 deletions(-) create mode 100644 Reactors/Reactor_1.json create mode 100644 Reactors/Reactor_2.json diff --git a/.vscode/launch.json b/.vscode/launch.json index 411a02c..095ed2d 100644 --- a/.vscode/launch.json +++ b/.vscode/launch.json @@ -1,6 +1,7 @@ { "version": "0.2.0", "configurations": [ + { "name": "Test_00", // Name of the current test "type": "debugpy", @@ -55,6 +56,38 @@ "0.9", // Discount factor employed in the MDP "--random-seed", "42"] // Random seed for the experiment + }, + + { + "name": "Custom_Test_1", // Name of the current test + "type": "debugpy", + "request": "launch", + "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) + "console": "integratedTerminal", + "args": [ + "--input-reactor", + "${workspaceFolder}/Reactors/Reactor_1.json", // Path to the reactor's JSON definition + "--gamma", + "0.9", // Discount factor employed in the MDP + "--random-seed", + "42" + ] // Random seed for the experiment + }, + + { + "name": "Custom_Test_2", // Name of the current test + "type": "debugpy", + "request": "launch", + "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) + "console": "integratedTerminal", + "args": [ + "--input-reactor", + "${workspaceFolder}/Reactors/Reactor_2.json", // Path to the reactor's JSON definition + "--gamma", + "0.9", // Discount factor employed in the MDP + "--random-seed", + "42" + ] // Random seed for the experiment } ] } \ No newline at end of file diff --git a/ControlModule.py b/ControlModule.py index 931d5df..ddf866f 100644 --- a/ControlModule.py +++ b/ControlModule.py @@ -9,8 +9,12 @@ def __init__(self): pass @staticmethod - def generate_P(probs: np.ndarray, n_states: np.int32 = 100) -> np.ndarray: + def generate_P( + probs: np.ndarray, n_states: np.int32 = 100, n_actions: np.int32 = 3 + ) -> np.ndarray: """Function that generates the probabilities (transition) matrix""" + + # Initialization of the Probability matrix matrix_P = np.zeros((3, n_states, n_states), dtype=np.float64) probs_decrease = probs[0] @@ -18,60 +22,70 @@ def generate_P(probs: np.ndarray, n_states: np.int32 = 100) -> np.ndarray: probs_increase = probs[2] # ---------------- DECREASE ---------------- - for estado_inicial in range(n_states): - for estado_final in range(n_states): - if estado_inicial == estado_final: - matrix_P[0][estado_inicial][estado_final] = probs_decrease[2] - elif estado_inicial - 1 == estado_final: - matrix_P[0][estado_inicial][estado_final] = probs_decrease[1] - elif estado_inicial - 2 == estado_final: - matrix_P[0][estado_inicial][estado_final] = probs_decrease[0] + for initial_state in range(n_states): + for final_state in range(n_states): + if initial_state == final_state: + matrix_P[0][initial_state][final_state] = probs_decrease[2] + elif initial_state - 1 == final_state: + matrix_P[0][initial_state][final_state] = probs_decrease[1] + elif initial_state - 2 == final_state: + matrix_P[0][initial_state][final_state] = probs_decrease[0] # ---------------- MAINTAIN ---------------- - for estado_inicial in range(n_states): - for estado_final in range(n_states): - if estado_inicial == estado_final: - matrix_P[1][estado_inicial][estado_final] = probs_maintain[1] - elif estado_inicial + 1 == estado_final: - matrix_P[1][estado_inicial][estado_final] = probs_maintain[2] - elif estado_inicial - 1 == estado_final: - matrix_P[1][estado_inicial][estado_final] = probs_maintain[0] + for initial_state in range(n_states): + for final_state in range(n_states): + if initial_state == final_state: + matrix_P[1][initial_state][final_state] = probs_maintain[1] + elif initial_state + 1 == final_state: + matrix_P[1][initial_state][final_state] = probs_maintain[2] + elif initial_state - 1 == final_state: + matrix_P[1][initial_state][final_state] = probs_maintain[0] # ---------------- INCREASE ---------------- - for estado_inicial in range(n_states): - for estado_final in range(n_states): - if estado_inicial == estado_final: - matrix_P[2][estado_inicial][estado_final] = probs_increase[0] - elif estado_inicial + 1 == estado_final: - matrix_P[2][estado_inicial][estado_final] = probs_increase[1] - elif estado_inicial + 2 == estado_final: - matrix_P[2][estado_inicial][estado_final] = probs_increase[2] - - # ---------------- CORRECCIÓN DE BORDES ---------------- - matrix_P[0][0][0] += probs_decrease[0] + probs_decrease[1] - matrix_P[0][1][0] += probs_decrease[0] - - matrix_P[1][0][0] += probs_maintain[0] - matrix_P[1][n_states - 1][n_states - 1] += probs_maintain[2] - - matrix_P[2][n_states - 2][n_states - 1] += probs_increase[2] - matrix_P[2][n_states - 1][n_states - 1] += probs_increase[1] + probs_increase[2] - + for initial_state in range(n_states): + for final_state in range(n_states): + if initial_state == final_state: + matrix_P[2][initial_state][final_state] = probs_increase[0] + elif initial_state + 1 == final_state: + matrix_P[2][initial_state][final_state] = probs_increase[1] + elif initial_state + 2 == final_state: + matrix_P[2][initial_state][final_state] = probs_increase[2] + + # ---------------- EDGE CORRECTION ---------------- + # We count how many non-zero elements are to normalize the whole line + # Returns the index of the position where the prob is not null + for action in range(n_actions): + # Transposing the matrix + for initial_state in range(n_states): + non_zero_elts = np.where(matrix_P[action, initial_state, :] != 0)[0] + # We divide all elements by the sum of the total probabilities + if len(non_zero_elts) > 0: + alpha_val = 0 + for idx in non_zero_elts: + alpha_val += matrix_P[action][initial_state][idx] + + matrix_P[action][initial_state] *= 1 / alpha_val + # Adding an average probability if all elts are null + if len(non_zero_elts) == 0: + for end_state in range(n_states): + matrix_P[action][initial_state][end_state] = 1 / n_states + + check_stochastic(matrix_P) return matrix_P @staticmethod def generate_R(demand_t: np.float64, n_states: np.int32 = 100) -> np.ndarray: """Function that generates the rewards (costs) matrix""" - demand = np.float64(demand_t) + # Initialization of the Rewards matrix matrix_R = np.zeros((3, n_states, n_states), dtype=np.float64) # ---------------- DECREASE ---------------- for estado_inicial in range(n_states): for estado_final in range(n_states): - delta_t = demand - (estado_final / n_states) + delta_t = demand_t - (estado_final / n_states) coste = abs(delta_t) - if demand > (estado_inicial / n_states): + if demand_t > (estado_inicial / n_states): coste *= 2 matrix_R[0][estado_inicial][estado_final] = -coste @@ -79,7 +93,7 @@ def generate_R(demand_t: np.float64, n_states: np.int32 = 100) -> np.ndarray: # ---------------- MAINTAIN ---------------- for estado_inicial in range(n_states): for estado_final in range(n_states): - delta_t = demand - (estado_final / n_states) + delta_t = demand_t - (estado_final / n_states) coste = abs(delta_t) matrix_R[1][estado_inicial][estado_final] = -coste @@ -87,10 +101,10 @@ def generate_R(demand_t: np.float64, n_states: np.int32 = 100) -> np.ndarray: # ---------------- INCREASE ---------------- for estado_inicial in range(n_states): for estado_final in range(n_states): - delta_t = demand - (estado_final / n_states) + delta_t = demand_t - (estado_final / n_states) coste = abs(delta_t) - if demand < (estado_inicial / n_states): + if demand_t < (estado_inicial / n_states): coste *= 2 matrix_R[2][estado_inicial][estado_final] = -coste @@ -107,6 +121,7 @@ def control_iteration(P, R, estado_actual, gamma, max_iter=1000) -> np.int32: max_iter=max_iter, ) + # Running the algorithm to converge to a single best action mdp.run() return np.int32(mdp.policy[estado_actual]) @@ -120,32 +135,70 @@ def control_loop( gamma: np.float64, ) -> np.ndarray: """Function that computes all the required iterations (control-loop) to satisfy the power demand""" - respuesta = np.zeros_like(a=demand, dtype=np.float64) + # Array of zeros of size 'demand' + response = np.zeros_like(a=demand, dtype=np.float64) current_state = np.int32(0) + # Matrix with the positions each action can apply action_deltas = [ np.array([-2, -1, 0], dtype=np.int32), np.array([-1, 0, 1], dtype=np.int32), np.array([0, 1, 2], dtype=np.int32), ] - ControlModule._P = ControlModule.generate_P(probs, n_states) + P = ControlModule.generate_P(probs, n_states, n_actions) for t in range(demand.shape[0]): - ControlModule._demand_t = np.float64(demand[t]) - ControlModule._current_state = current_state - ControlModule._R = ControlModule.generate_R(ControlModule._demand_t, n_states) + # Generation of the matrix R + R = ControlModule.generate_R(demand[t], n_states) + + # Calculation of the best action action = ControlModule.control_iteration( - P=ControlModule._P, - R=ControlModule._R, - estado_actual=ControlModule._current_state, + P=P, + R=R, + estado_actual=current_state, gamma=gamma, ) + # Normalizing probs array + sum = np.sum(probs[action]) + if sum == 0: + for i in range(len(probs[action])): + probs[action][i] = 1 / sum + else: + probs[action] /= sum + + # Execution of the action (taking account uncertainty) state_increment = np.random.choice(a=action_deltas[action], p=probs[action]) current_state = current_state + state_increment - current_state = np.int32(np.clip(a=current_state, a_min=0, a_max=n_states - 1)) - respuesta[t] = current_state / n_states - return respuesta + # The current state must be within the borders + current_state = max(0, min(n_states - 1, current_state)) + + # We add the next state to the array to return it + response[t] = current_state / n_states + + return response + + +import numpy as np + + +def check_stochastic(P, tol=1e-6): + P = np.array(P) + assert P.ndim == 3, f"P debe ser (A, S, S), tiene forma {P.shape}" + A, S, _ = P.shape + ok = True + for a in range(A): + row_sums = P[a].sum(axis=1) + bad = np.where(np.abs(row_sums - 1.0) > tol)[0] + if len(bad) > 0: + ok = False + for s in bad: + print(f" Acción {a}, Estado {s}: suma = {row_sums[s]:.8f}") + if ok: + print("✅ Matriz estocástica: todas las filas suman 1.") + else: + print("❌ Matriz NO estocástica.") + return ok diff --git a/Reactor.py b/Reactor.py index 691c072..874505a 100644 --- a/Reactor.py +++ b/Reactor.py @@ -1,51 +1,65 @@ # Import required dependencies import numpy as np + class Reactor: - def __init__(self, - model: str, - effective_section: np.float64, - neutron_flux: np.float64, - core_volume: np.float64, - fision_energy: np.float64, - probabilities: dict): - """ Constructor of the Reactor class """ - self.model = model + def __init__( + self, + model: str, + effective_section: np.float64, + neutron_flux: np.float64, + core_volume: np.float64, + fision_energy: np.float64, + probabilities: dict, + ): + """Constructor of the Reactor class""" + self.model = model self.effective_section = effective_section - self.neutron_flux = neutron_flux - self.core_volume = core_volume - self.fision_energy = fision_energy - self.probabilities = probabilities - self.max_power = self.compute_max_power() - self.k = self.compute_k() + self.neutron_flux = neutron_flux + self.core_volume = core_volume + self.fision_energy = fision_energy + self.probabilities = probabilities + self.max_power = self.compute_max_power() + self.k = self.compute_k() def __str__(self) -> str: - """ Overloading of the native __str__ function to print the class instances """ - _str = f"Model: {self.model}\n" + """Overloading of the native __str__ function to print the class instances""" + _str = f"Model: {self.model}\n" _str += f"Effective section: {self.effective_section} cm^-1\n" _str += f"Neutron flux: {self.neutron_flux} neutrons / (cm^2 · s)\n" _str += f"Core volume: {self.core_volume} cm^3\n" _str += f"Fision energy: {self.fision_energy} J\n" - _str += f"Probabilities: {self.probabilities}" + _str += f"Probabilities: {self.probabilities}\n" return _str def compute_max_power(self) -> np.float64: - """ Computes the maximum power of a reactor based on its physical features """ - return np.float64(self.effective_section * self.neutron_flux * self.core_volume * self.fision_energy) - + """Computes the maximum power of a reactor based on its physical features""" + # Returns the max power of the reactor (Equation 1) + return np.float64( + self.effective_section + * self.neutron_flux + * self.core_volume + * self.fision_energy + ) + def compute_k(self) -> np.float64: - """ Computes the value of the k-constant """ + """Computes the value of the k-constant""" self.P_max = self.compute_max_power() - return - np.log(10** -6/self.P_max) - + # Applies the formula for 'k' and returns it + return -np.log(10**-6 / self.P_max) + def compute_power(self, control_bars_insertion: np.float64) -> np.float64: - """ Computes the power delivered (%) by the reactor based on the % of control-bars inserted """ + """Computes the power delivered (%) by the reactor based on the % of control-bars inserted""" k = self.compute_k() - return np.float64(self.P_max * (np.e**(-k * control_bars_insertion))) - + # Applies the function to compute the power and returns it (Equation 2) + return np.float64(self.P_max * (np.e ** (-k * control_bars_insertion))) + def compute_control_bars_insertion(self, power: np.float64) -> np.float64: - """ Computes the % of controls-bars inserted based on the % of power delivered by the reactor """ + """Computes the % of controls-bars inserted based on the % of power delivered by the reactor""" min_power = 10**-6 / self.max_power + # np.clip() forces 'power' to stay within the limits [min_power, 1.0] power = np.clip(a=power, a_min=min_power, a_max=1.0) + + # Returns the percentage of bars to be inserted (B, inverse Equation 2) insertion = -np.log(power) / self.k return np.float64(np.clip(a=insertion, a_min=0.0, a_max=1.0)) diff --git a/Reactors/Reactor_1.json b/Reactors/Reactor_1.json new file mode 100644 index 0000000..810dbea --- /dev/null +++ b/Reactors/Reactor_1.json @@ -0,0 +1,24 @@ +{ + "model": "UP-REACTOR", + "effective_section": 17.6, + "neutron_flux": 5e6, + "core_volume": 9.42e6, + "fision_energy": 3.2e-12, + "probabilities": { + "decrease": [ + 0.34, + 0.33, + 0.33 + ], + "maintain": [ + 0.05, + 0.6, + 0.35 + ], + "increase": [ + 0.05, + 0.9, + 0.05 + ] + } +} \ No newline at end of file diff --git a/Reactors/Reactor_2.json b/Reactors/Reactor_2.json new file mode 100644 index 0000000..3c62d63 --- /dev/null +++ b/Reactors/Reactor_2.json @@ -0,0 +1,24 @@ +{ + "model": "BAD-REACTOR", + "effective_section": 5, + "neutron_flux": 5e13, + "core_volume": 9.42e6, + "fision_energy": 3.2e-11, + "probabilities": { + "decrease": [ + 0, + 0.33, + 0 + ], + "maintain": [ + 0.05, + 0.6, + 0.35 + ], + "increase": [ + 0.05, + 0.9, + 0.05 + ] + } +} \ No newline at end of file diff --git a/main.py b/main.py index 66815b1..52fe7b2 100644 --- a/main.py +++ b/main.py @@ -69,7 +69,7 @@ def main() -> None: ) # Make a radar-plot with the reactor probabilities - plot_reactor_as_radar(probs=probs) + # plot_reactor_as_radar(probs=probs) # Generate a random power demand demand = generate_demand(n_samples=512) From 9ff4ba7beeb2ecf5cb4348a0c506824a91d7bb1d Mon Sep 17 00:00:00 2001 From: codeMGL Date: Fri, 15 May 2026 23:35:32 +0200 Subject: [PATCH 11/11] Reactores --- .vscode/launch.json | 233 ++++++++++++++++++++++++++-------------- ControlModule.py | 25 +---- Reactors/Reactor_3.json | 24 +++++ Reactors/Reactor_4.json | 24 +++++ Reactors/Reactor_5.json | 24 +++++ Reactors/Reactor_6.json | 24 +++++ main.py | 4 +- 7 files changed, 250 insertions(+), 108 deletions(-) create mode 100644 Reactors/Reactor_3.json create mode 100644 Reactors/Reactor_4.json create mode 100644 Reactors/Reactor_5.json create mode 100644 Reactors/Reactor_6.json diff --git a/.vscode/launch.json b/.vscode/launch.json index 095ed2d..8c52285 100644 --- a/.vscode/launch.json +++ b/.vscode/launch.json @@ -1,93 +1,162 @@ { "version": "0.2.0", "configurations": [ + { + "name": "Test_00", // Name of the current test + "type": "debugpy", + "request": "launch", + "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) + "console": "integratedTerminal", + "args": [ + "--input-reactor", + "${workspaceFolder}/Reactors/R0.json", // Path to the reactor's JSON definition + "--gamma", + "0.9", // Discount factor employed in the MDP + "--random-seed", + "42" + ] // Random seed for the experiment + }, - { - "name": "Test_00", // Name of the current test - "type": "debugpy", - "request": "launch", - "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) - "console": "integratedTerminal", - "args": ["--input-reactor", - "${workspaceFolder}/Reactors/R0.json", // Path to the reactor's JSON definition - "--gamma", - "0.9", // Discount factor employed in the MDP - "--random-seed", - "42"] // Random seed for the experiment - }, + { + "name": "Test_01", // Name of the current test + "type": "debugpy", + "request": "launch", + "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) + "console": "integratedTerminal", + "args": [ + "--input-reactor", + "${workspaceFolder}/Reactors/R1.json", // Path to the reactor's JSON definition + "--gamma", + "0.9", // Discount factor employed in the MDP + "--random-seed", + "42" + ] // Random seed for the experiment + }, - { - "name": "Test_01", // Name of the current test - "type": "debugpy", - "request": "launch", - "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) - "console": "integratedTerminal", - "args": ["--input-reactor", - "${workspaceFolder}/Reactors/R1.json", // Path to the reactor's JSON definition - "--gamma", - "0.9", // Discount factor employed in the MDP - "--random-seed", - "42"] // Random seed for the experiment - }, + { + "name": "Test_02", // Name of the current test + "type": "debugpy", + "request": "launch", + "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) + "console": "integratedTerminal", + "args": [ + "--input-reactor", + "${workspaceFolder}/Reactors/R2.json", // Path to the reactor's JSON definition + "--gamma", + "0.9", // Discount factor employed in the MDP + "--random-seed", + "42" + ] // Random seed for the experiment + }, - { - "name": "Test_02", // Name of the current test - "type": "debugpy", - "request": "launch", - "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) - "console": "integratedTerminal", - "args": ["--input-reactor", - "${workspaceFolder}/Reactors/R2.json", // Path to the reactor's JSON definition - "--gamma", - "0.9", // Discount factor employed in the MDP - "--random-seed", - "42"] // Random seed for the experiment - }, + { + "name": "Test_03", // Name of the current test + "type": "debugpy", + "request": "launch", + "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) + "console": "integratedTerminal", + "args": [ + "--input-reactor", + "${workspaceFolder}/Reactors/R3.json", // Path to the reactor's JSON definition + "--gamma", + "0.9", // Discount factor employed in the MDP + "--random-seed", + "42" + ] // Random seed for the experiment + }, - { - "name": "Test_03", // Name of the current test - "type": "debugpy", - "request": "launch", - "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) - "console": "integratedTerminal", - "args": ["--input-reactor", - "${workspaceFolder}/Reactors/R3.json", // Path to the reactor's JSON definition - "--gamma", - "0.9", // Discount factor employed in the MDP - "--random-seed", - "42"] // Random seed for the experiment - }, + { + "name": "Custom_Test_1", // Name of the current test + "type": "debugpy", + "request": "launch", + "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) + "console": "integratedTerminal", + "args": [ + "--input-reactor", + "${workspaceFolder}/Reactors/Reactor_1.json", // Path to the reactor's JSON definition + "--gamma", + "0.9", // Discount factor employed in the MDP + "--random-seed", + "42" + ] // Random seed for the experiment + }, - { - "name": "Custom_Test_1", // Name of the current test - "type": "debugpy", - "request": "launch", - "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) - "console": "integratedTerminal", - "args": [ - "--input-reactor", - "${workspaceFolder}/Reactors/Reactor_1.json", // Path to the reactor's JSON definition - "--gamma", - "0.9", // Discount factor employed in the MDP - "--random-seed", - "42" - ] // Random seed for the experiment - }, + { + "name": "Custom_Test_2", // Name of the current test + "type": "debugpy", + "request": "launch", + "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) + "console": "integratedTerminal", + "args": [ + "--input-reactor", + "${workspaceFolder}/Reactors/Reactor_2.json", // Path to the reactor's JSON definition + "--gamma", + "0.9", // Discount factor employed in the MDP + "--random-seed", + "42" + ] // Random seed for the experiment + }, - { - "name": "Custom_Test_2", // Name of the current test - "type": "debugpy", - "request": "launch", - "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) - "console": "integratedTerminal", - "args": [ - "--input-reactor", - "${workspaceFolder}/Reactors/Reactor_2.json", // Path to the reactor's JSON definition - "--gamma", - "0.9", // Discount factor employed in the MDP - "--random-seed", - "42" - ] // Random seed for the experiment - } + { + "name": "Custom_Test_3", // Name of the current test + "type": "debugpy", + "request": "launch", + "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) + "console": "integratedTerminal", + "args": [ + "--input-reactor", + "${workspaceFolder}/Reactors/Reactor_3.json", // Path to the reactor's JSON definition + "--gamma", + "0.9", // Discount factor employed in the MDP + "--random-seed", + "42" + ] // Random seed for the experiment + }, + + { + "name": "Custom_Test_4", // Name of the current test + "type": "debugpy", + "request": "launch", + "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) + "console": "integratedTerminal", + "args": [ + "--input-reactor", + "${workspaceFolder}/Reactors/Reactor_4.json", // Path to the reactor's JSON definition + "--gamma", + "0.9", // Discount factor employed in the MDP + "--random-seed", + "42" + ] // Random seed for the experiment + }, + { + "name": "Custom_Test_5", // Name of the current test + "type": "debugpy", + "request": "launch", + "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) + "console": "integratedTerminal", + "args": [ + "--input-reactor", + "${workspaceFolder}/Reactors/Reactor_5.json", // Path to the reactor's JSON definition + "--gamma", + "0.9", // Discount factor employed in the MDP + "--random-seed", + "42" + ] // Random seed for the experiment + }, + { + "name": "Custom_Test_6", // Name of the current test + "type": "debugpy", + "request": "launch", + "program": "${workspaceFolder}/main.py", // Main Python file (do not modify) + "console": "integratedTerminal", + "args": [ + "--input-reactor", + "${workspaceFolder}/Reactors/Reactor_6.json", // Path to the reactor's JSON definition + "--gamma", + "0.9", // Discount factor employed in the MDP + "--random-seed", + "42" + ] // Random seed for the experiment + } ] } \ No newline at end of file diff --git a/ControlModule.py b/ControlModule.py index ddf866f..2fe1e76 100644 --- a/ControlModule.py +++ b/ControlModule.py @@ -70,7 +70,6 @@ def generate_P( for end_state in range(n_states): matrix_P[action][initial_state][end_state] = 1 / n_states - check_stochastic(matrix_P) return matrix_P @staticmethod @@ -78,7 +77,7 @@ def generate_R(demand_t: np.float64, n_states: np.int32 = 100) -> np.ndarray: """Function that generates the rewards (costs) matrix""" # Initialization of the Rewards matrix matrix_R = np.zeros((3, n_states, n_states), dtype=np.float64) - + # Method explained at the memory of the project # ---------------- DECREASE ---------------- for estado_inicial in range(n_states): for estado_final in range(n_states): @@ -180,25 +179,3 @@ def control_loop( response[t] = current_state / n_states return response - - -import numpy as np - - -def check_stochastic(P, tol=1e-6): - P = np.array(P) - assert P.ndim == 3, f"P debe ser (A, S, S), tiene forma {P.shape}" - A, S, _ = P.shape - ok = True - for a in range(A): - row_sums = P[a].sum(axis=1) - bad = np.where(np.abs(row_sums - 1.0) > tol)[0] - if len(bad) > 0: - ok = False - for s in bad: - print(f" Acción {a}, Estado {s}: suma = {row_sums[s]:.8f}") - if ok: - print("✅ Matriz estocástica: todas las filas suman 1.") - else: - print("❌ Matriz NO estocástica.") - return ok diff --git a/Reactors/Reactor_3.json b/Reactors/Reactor_3.json new file mode 100644 index 0000000..ad91be1 --- /dev/null +++ b/Reactors/Reactor_3.json @@ -0,0 +1,24 @@ +{ + "model": "AVERAGE-REACTOR", + "effective_section": 5, + "neutron_flux": 5e13, + "core_volume": 9.42e6, + "fision_energy": 3.2e-11, + "probabilities": { + "decrease": [ + 1, + 1, + 1 + ], + "maintain": [ + 1, + 1, + 1 + ], + "increase": [ + 1, + 1, + 1 + ] + } +} \ No newline at end of file diff --git a/Reactors/Reactor_4.json b/Reactors/Reactor_4.json new file mode 100644 index 0000000..a23f36d --- /dev/null +++ b/Reactors/Reactor_4.json @@ -0,0 +1,24 @@ +{ + "model": "ALWAYS-MAINTAIN", + "effective_section": 5, + "neutron_flux": 5e13, + "core_volume": 9.42e6, + "fision_energy": 3.2e-11, + "probabilities": { + "decrease": [ + 0.1, + 0.1, + 0.8 + ], + "maintain": [ + 0.05, + 0.9, + 0.05 + ], + "increase": [ + 0.7, + 0.15, + 0.1 + ] + } +} \ No newline at end of file diff --git a/Reactors/Reactor_5.json b/Reactors/Reactor_5.json new file mode 100644 index 0000000..62d9711 --- /dev/null +++ b/Reactors/Reactor_5.json @@ -0,0 +1,24 @@ +{ + "model": "IDEAL-REACTOR", + "effective_section": 5, + "neutron_flux": 5e13, + "core_volume": 9.42e6, + "fision_energy": 3.2e-11, + "probabilities": { + "decrease": [ + 0, + 1, + 0 + ], + "maintain": [ + 0, + 1, + 0 + ], + "increase": [ + 0, + 1, + 0 + ] + } +} \ No newline at end of file diff --git a/Reactors/Reactor_6.json b/Reactors/Reactor_6.json new file mode 100644 index 0000000..5aa6ce8 --- /dev/null +++ b/Reactors/Reactor_6.json @@ -0,0 +1,24 @@ +{ + "model": "USUALLY-DOWN", + "effective_section": 1000, + "neutron_flux": 5e23, + "core_volume": 9.42e16, + "fision_energy": 3.2e-5, + "probabilities": { + "decrease": [ + 0.4, + 0.6, + 0 + ], + "maintain": [ + 0.4, + 0.6, + 0 + ], + "increase": [ + 0.6, + 0.4, + 0 + ] + } +} \ No newline at end of file diff --git a/main.py b/main.py index 52fe7b2..d3eb90b 100644 --- a/main.py +++ b/main.py @@ -67,9 +67,9 @@ def main() -> None: ], dtype=np.float64, ) - + # Make a radar-plot with the reactor probabilities - # plot_reactor_as_radar(probs=probs) + plot_reactor_as_radar(probs=probs) # Generate a random power demand demand = generate_demand(n_samples=512)