Problem in Jacobian matrix calculation in optimization problem using scipy

  python, scipy-optimize-minimize

I’m trying to solve and optimization problem with scipy.optimize.minimize but it seems there is an error in library when it tries to compute the Jacobian matrix.

# Librerias necesarias 
import numpy as np
import matplotlib.pyplot as plt 
from scipy.optimize import minimize  

# Constantes de eficiencia de los conversores del EH
# CHP 
n_W = 0.3 # Electricidad 
n_Q = 0.45 # Termal 

# AB 
n_T = 0.8 # Termal 

# CERG 
COP_CERG = 3 

# WARG 
COP_WARG = 0.7  

# Limites superiores 
v_p_CHP = np.array([0,300,450,0]) 
v_p_AB = np.array([0,0,900,0])
v_p_CERG = np.array([0,0,0,400]) 
v_p_WARG = np.array([0,0,0,400]) 

# Limites inferiores 
v_d_CHP = np.array([0,300/2,450/2,0]) 
v_d_AB = np.array([0,0,0,0])
v_d_CERG = np.array([0,0,0,0]) 
v_d_WARG = np.array([0,0,0,0])  

# Precio del gas 
P_g = 0.02 #kW/h 

# Precio electricidad 
P_e = np.array([0,40,20,35,40,60,110,120,115,75,80,120,90,70,60,40])  
# Demanda de electricidad, calefacción y refrigeración 
D_h = np.array([0,200,320,400,580,550,520,500,470,420,320,200])
D_e = np.array([0,150,210,180,175,180,210,160]) 
D_c = np.array([0,0,25,50,75,125,150,125,75,50,25,0])  

# Matrices de espeficicacion por componente 
M = 4 # Numero de vectores de energia (Electricidad,gas,calefaccion y regrigeracion) 
zero_vect_v = np.zeros((1,4)) 
zero_vect_h = np.zeros((4,1))

# H_CHP  
H_CHP = np.array([[0,0,0,0],
                  [n_W,0,0,0],
                  [n_Q,0,0,0], 
                  [0,0,0,0]]) 
# H_AB 
H_AB = np.array([[0,0,0,0],
                 [0,0,0,0],
                 [n_T,0,0,0], 
                 [0,0,0,0]]) 
# H_WARG  
H_WARG = np.array([[0,0,0,0],
                   [0,0,0,0],
                   [0,0,0,0], 
                   [0,0,COP_WARG,0]])
# H_CERG  
H_CERG = np.array([[0,0,0,0],
                   [0,0,0,0],
                   [0,0,0,0], 
                   [0,COP_CERG,0,0]]) 

# Matrices de interconección 
def T_o_AB(X): 
  return np.diag([0,0,X[5],0])  

def T_o_CHP(X): 
  return np.diag([0,X[2],X[3],0]) 

T_o_CERG = np.diag([0,0,0,1]) 
T_o_WARG = np.diag([0,0,0,1]) 

def T_o_EDS(X): 
  return np.diag([0,X[0],0,0]) 

def T_AB_i(X): 
  return np.diag([1-X[1],0,0,0]) 

def T_CHP_i(X): 
  return np.diag([X[1],0,0,0]) 

def T_CERG_CHP(X): 
  return np.diag([0,1-X[2],0,0]) 

def T_CERG_i(X): 
  return np.diag([0,1,1-X[0],0]) 

def T_WARG_AB(X): 
  return np.diag([0,0,1-X[5],0]) 

def T_WARG_CHP(X): 
  return np.diag([0,0,X[4],0])  

T_zero = np.zeros((M,M)) 
I = np.identity(M)

# Costo de la electricidad por horas
# Dominio del tiempo (rangos de tiempo) 
P_tiempo = np.arange(1,25,1) 
# Definicion del vector de demanda eléctrica 
P_e_vect_1 = np.zeros((P_tiempo.shape[0]))
P_e_vect_1[(P_tiempo>=0) & (P_tiempo<=2)] = P_e[1]  
P_e_vect_1[(P_tiempo>2) & (P_tiempo<=6)] = P_e[2] 
P_e_vect_1[(P_tiempo>6) & (P_tiempo<=7)] = P_e[3]
P_e_vect_1[(P_tiempo>7) & (P_tiempo<=8)] = P_e[4] 
P_e_vect_1[(P_tiempo>8) & (P_tiempo<=9)] = P_e[5]
P_e_vect_1[(P_tiempo>9) & (P_tiempo<=10)] = P_e[6]
P_e_vect_1[(P_tiempo>10) & (P_tiempo<=12)] = P_e[7] 
P_e_vect_1[(P_tiempo>12) & (P_tiempo<=13)] = P_e[8]  

P_e_vect_1[(P_tiempo>13) & (P_tiempo<=15)] = P_e[9]  
P_e_vect_1[(P_tiempo>15) & (P_tiempo<=19)] = P_e[10]  
P_e_vect_1[(P_tiempo>19) & (P_tiempo<=20)] = P_e[11]  
P_e_vect_1[(P_tiempo>20) & (P_tiempo<=21)] = P_e[12]  
P_e_vect_1[(P_tiempo>21) & (P_tiempo<=22)] = P_e[13]  
P_e_vect_1[(P_tiempo>22) & (P_tiempo<=23)] = P_e[14]    
P_e_vect_1[(P_tiempo>23) & (P_tiempo<=24)] = P_e[15]    

P_e_vect_1 = P_e_vect_1/1000

# Demanda 
# Definicion del vector para demanda electricidad
P_e_vect_d_1 = np.zeros((P_tiempo.shape[0]))
P_e_vect_d_1[(P_tiempo>=0) & (P_tiempo<=8)] = D_e[1]  
P_e_vect_d_1[(P_tiempo>8) & (P_tiempo<=12)] = D_e[2] 
P_e_vect_d_1[(P_tiempo>12) & (P_tiempo<=15)] = D_e[3]
P_e_vect_d_1[(P_tiempo>15) & (P_tiempo<=18)] = D_e[4] 
P_e_vect_d_1[(P_tiempo>18) & (P_tiempo<=19)] = D_e[5]
P_e_vect_d_1[(P_tiempo>19) & (P_tiempo<=20)] = D_e[6] 
P_e_vect_d_1[(P_tiempo>20) & (P_tiempo<=24)] = D_e[7]   

# Definicion del vector para demanda refrigeracion
P_c_vect_d_1 = np.zeros((P_tiempo.shape[0])) 
P_c_vect_d_1[(P_tiempo>=0) & (P_tiempo<=9)] = D_c[1] 
P_c_vect_d_1[(P_tiempo>=9) & (P_tiempo<=10)] = D_c[2]
P_c_vect_d_1[(P_tiempo>=10) & (P_tiempo<=12)] = D_c[3]
P_c_vect_d_1[(P_tiempo>=12) & (P_tiempo<=13)] = D_c[4]
P_c_vect_d_1[(P_tiempo>=13) & (P_tiempo<=14)] = D_c[5] 
P_c_vect_d_1[(P_tiempo>=14) & (P_tiempo<=15)] = D_c[6] 
P_c_vect_d_1[(P_tiempo>=15) & (P_tiempo<=16)] = D_c[7] 
P_c_vect_d_1[(P_tiempo>=16) & (P_tiempo<=17)] = D_c[8] 
P_c_vect_d_1[(P_tiempo>=17) & (P_tiempo<=18)] = D_c[9] 
P_c_vect_d_1[(P_tiempo>=18) & (P_tiempo<=19)] = D_c[10]
P_c_vect_d_1[(P_tiempo>=19) & (P_tiempo<=24)] = D_c[11] 

# Definicion del vector para demanda calefaccion 
P_h_vect_d_1 = np.zeros((P_tiempo.shape[0])) 
P_h_vect_d_1[(P_tiempo>=0) & (P_tiempo<=6)] = D_h[1]   
P_h_vect_d_1[(P_tiempo>=6) & (P_tiempo<=7)] = D_h[2] 
P_h_vect_d_1[(P_tiempo>=7) & (P_tiempo<=8)] = D_h[3] 
P_h_vect_d_1[(P_tiempo>=8) & (P_tiempo<=12)] = D_h[4] 
P_h_vect_d_1[(P_tiempo>=12) & (P_tiempo<=13)] = D_h[5] 
P_h_vect_d_1[(P_tiempo>=13) & (P_tiempo<=14)] = D_h[6]  
P_h_vect_d_1[(P_tiempo>=14) & (P_tiempo<=16)] = D_h[7] 
P_h_vect_d_1[(P_tiempo>=16) & (P_tiempo<=17)] = D_h[8] 
P_h_vect_d_1[(P_tiempo>=17) & (P_tiempo<=21)] = D_h[9]  
P_h_vect_d_1[(P_tiempo>=21) & (P_tiempo<=23)] = D_h[10] 
P_h_vect_d_1[(P_tiempo>=23) & (P_tiempo<=24)] = D_h[11] 

# Arreglo de dimensiones 
P_e_vect_d_1 = np.array([P_e_vect_d_1]) 
P_c_vect_d_1 = np.array([P_c_vect_d_1]) 
P_h_vect_d_1 = np.array([P_h_vect_d_1])

# Concatenar los vectores de demanda 
Vout = np.concatenate((np.zeros((1,24)),P_e_vect_d_1,P_h_vect_d_1,P_c_vect_d_1),axis=0)   
Vout = Vout.T 

# Construccion de matriz H 
def H_matrix(X): 
  return T_o_EDS(X)*I + T_o_AB(X)*H_AB*T_AB_i(X)*I + T_o_CHP(X)*H_CHP*T_CHP_i(X)*I + [email protected]_CERG*(T_CERG_i(X)*I + T_CERG_CHP(X)*H_CHP*T_CHP_i(X)*I) + [email protected]_WARG*(T_WARG_AB(X)*H_AB*T_AB_i(X)*I + T_WARG_CHP(X)*H_CHP*T_CHP_i(X)*I) 

save_results = [] # Lista para guardar los resultados  

# Problema de optimizacion 
for t in P_tiempo: 
  # Definicion de la funcion de costo 
  def objective_function(X): 
    return P_e_vect_1[t-1]*X[7] + P_g*X[6]  

  # Restricciones  
  cons = ({'type': 'eq',
        'fun': lambda X: np.dot(H_matrix(X),np.array([[X[6]],[X[7]],[0],[0]])) - Vout[t-1,:].reshape((4,1))}, 
        {'type': 'eq',
        'fun': lambda X: np.dot(H_CHP,np.array([[X[1]*X[6]],[0],[0],[0]])) - np.array([[0],[X[8]],[X[9]],[0]])}, 
        {'type': 'eq',
        'fun': lambda X: np.dot(H_AB,np.array([[(1-X[1])*X[6]],[0],[0],[0]])) - np.array([[0],[0],[X[10]],[0]])},  
        {'type': 'eq',
        'fun': lambda X: np.dot(H_WARG,np.array([[0],[0],[X[4]*X[9]+(1-X[5])*X[10]],[0]])) - np.array([[0],[0],[0],[X[12]]])},
        {'type': 'eq',
        'fun': lambda X: np.dot(H_CERG,np.array([[0],[(1-X[0])*X[7]+(1-X[2])*X[8]],[0],[0]])) - np.array([[0],[0],[0],[X[11]]])},
        {'type': 'ineq',
        'fun': lambda X: -(np.array([[(X[0]-1)*X[7]],[0]]))},
        {'type': 'ineq',
        'fun': lambda X: -([email protected]([[X[1]*X[6]],[0],[0],[0]]) - v_p_CHP.reshape((4,1)))}, 
        {'type': 'ineq',
        'fun': lambda X: -([email protected]([[(1-X[1])*X[6]],[0],[0],[0]]) - v_p_AB.reshape((4,1)))}, 
        {'type': 'ineq',
        'fun': lambda X: -([email protected]([[0],[0],[X[4]*X[9]+(1-X[5])*X[10]],[0]]) - v_p_WARG.reshape((4,1)))}, 
        {'type': 'ineq',
        'fun': lambda X: -([email protected]([[0],[(1-X[0])*X[7]+(1-X[2])*X[8]],[0],[0]]) - v_p_CERG.reshape((4,1)))}, 
        {'type': 'ineq',
        'fun': lambda X: -(v_d_CHP.reshape((4,1)) - [email protected]([[X[1]*X[6]],[0],[0],[0]]))}, 
        {'type': 'ineq',
        'fun': lambda X: ([email protected]([[(1-X[1])*X[6]],[0],[0],[0]]))}, 
        {'type': 'ineq',
        'fun': lambda X: ([email protected]([[0],[0],[X[4]*X[9]+(1-X[5])*X[10]],[0]]))}, 
        {'type': 'ineq',
        'fun': lambda X: ([email protected]([[0],[(1-X[0])*X[7]+(1-X[2])*X[8]],[0],[0]]))})  
  
# Condicion inicial  
  if t == 1:
    x0 = np.array([1,1,1,1,1,1,0,0,0,0,0,0,0]).reshape((13,1))
  else: 
    x0 = save_results[t-2]

  # Bounds 
  bnds = ((0, 1),(0, 1),(0, 1),(0, 1),(0, 1),(0, 1),(0, None),(0, None),(0, None),(0, None),(0, None),(0, None),(0, None)) 
  # Problema de optimizacion
  sol = minimize(objective_function,x0,constraints=cons,bounds=bnds)  
  solution = sol.x
  save_results.append(solution)

It outputs these error ValueError: could not broadcast input array from shape (4,1) into shape (4), I don’t know if someone has experience these error before and could help me with an intuition for what it’s happening, clearly is something about the dimensions of the internal variables while solving the optimization problem but it isn’t clear how should I reformulate the problem to avoid this error, I appreciate the help.

Source: Python Questions

LEAVE A COMMENT