from __future__ import division from pyomo.environ import * from red.read import excel2net from red.create import * from red.create1 import * from red.makeBdc import makeBdc from utils.arr2dict import arr2dict from pyomo.environ import SolverFactory from pyomo.kernel import value from numpy import zeros, array from pandas import DataFrame, ExcelWriter def setmodel(file): print("{:=^100}".format("")) print("{:^100}".format(" Modelo de Predespacho Regional")) print("{:^100}".format(" Mercados Eléctricos de Centroamérica (c) 2020 ")) print("{:=^100}".format("")) # ============================================================================ # Importación de datos desde archivo Excel # ============================================================================ #Parametros de la linea print("Inicio del Script de Predespacho") print("Leyendo información de red...") net = excel2net(file) bus = setbus(net)#Nodos branch = setbranch(net, bus)#Set lineas bu = branch[:, 5]#potenica max de la linea bl = branch[:, 6]#potenica min de la linea xc = branch[:,3]#Reactancia de la linea rc = branch[:,4]#Resistencia de la linea nb = bus.shape[0] #Numero de nodos nbr = branch.shape[0]#Numero de lineas A = makeBdc(bus, branch)#Matriz incidente brnames = branchnames(bus, branch)#Nombre de las lineas inc = A.toarray()*-1 dg=read_D_G(file) print("Leyendo información de ofertas") #Funcion leer contratos no firmes fisicos flexibles ex_cnfff = readofertas_cnfffs(file) vnfff_ed= setvariable(ex_cnfff['energía_dec'])#Energia declarada con respecto al nodo vnfff_m_i= setvariable_s(ex_cnfff[['magnitud_i1','magnitud_i2','magnitud_i3','magnitud_i4','magnitud_i5']])#Magnitud de energia ofertada -flexibilizacion vnfff_m_r= setvariable_s(ex_cnfff[['magnitud_r1','magnitud_r2','magnitud_r3','magnitud_r4','magnitud_r5']]) vnfff_m_cvt=setvariable_s(ex_cnfff[['magnitud_cvt1','magnitud_cvt2','magnitud_cvt3','magnitud_cvt4','magnitud_cvt5']]) p_cnfffi= setvariable_s(ex_cnfff[['precio_i1','precio_i2','precio_i3','precio_i4','precio_i5']])#Precio de la ofertas de inyeccion p_cnfffr= setvariable_s(ex_cnfff[['precio_r1','precio_r2','precio_r3','precio_r4','precio_r5']])#Precio de la ofertas de retiro p_cnfffcvt= setvariable_s(ex_cnfff[['precio_cvt1','precio_cvt2','precio_cvt3','precio_cvt4','precio_cvt5']])#Precio de la ofertas de inyeccion k_cnfffcvt= setvariable(ex_cnfff['k'])#Precio de la ofertas de inyeccion NCFF = ex_cnfff.shape[0]#Numero de contratos firmes var_bin_cnfffr=MATRIZ_VNFFF_R(bus,ex_cnfff) var_bin_cnfffi=MATRIZ_VNFFF_I(bus,ex_cnfff) dem=setvariable(dg['demanda']) gen=setvariable(dg['generacion']) #Funcion leer parametros oferta de oportunidad inyeccion ex_ooi = readofertas_oois(file) vooi_m=setvariable_s(ex_ooi[['magnitud_ooi1','magnitud_ooi2','magnitud_ooi3','magnitud_ooi4','magnitud_ooi5']])#Ofertas con respecto al nodo p_ooi= setvariable_s(ex_ooi[['precio_ooi1','precio_ooi2','precio_ooi3','precio_ooi4','precio_ooi5']])#Precio de la ofertas NOI = ex_ooi.shape[0]#Numero de ofertas de inyeccion var_bin_ooi=MATRIZ_OOI(bus,ex_ooi) #Funcion leer parametros oferta de oportunidad retiro ex_oor = readofertas_oors(file) voor_m=setvariable_s(ex_oor[['magnitud_oor1','magnitud_oor2','magnitud_oor3','magnitud_oor4','magnitud_oor5']])#Ofertas con respecto al nodo p_oor= setvariable_s(ex_oor[['precio_oor1','precio_oor2','precio_oor3','precio_oor4','precio_oor5']])#Precio de la ofertas - flexibilizacion NOR = ex_oor.shape[0]#Numero de ofertas de retiro var_bin_oor=MATRIZ_OOR(bus,ex_oor) #Funcion leer parametros contratos firmes ex_cf = readofertas_cfs(file) vcf_ed=setvariable(ex_cf['energía_dec'])#Energia declarada con respecto al nodo vcf_pr=setvariable(ex_cf['potencia_req'])#Potencia requerida con respecto al nodo vcf_m=setvariable_s(ex_cf[['magnitu_cf1','magnitu_cf2','magnitu_cf3','magnitu_cf4','magnitu_cf5']])#Magnitud de energia ofertada -flexibilizacion vcf_p= setvariable_s(ex_cf[['precio_cf1','precio_cf2','precio_cf3','precio_cf4','precio_cf5']])#Precio de la ofertas NCF = ex_cf.shape[0]#Numero de contratos firmes var_bin_cfr=MATRIZ_CFR(bus, ex_cf) var_bin_cfi=MATRIZ_CFI(bus, ex_cf) # Inicio del modelo de optimización model=ConcreteModel() #sets model.i=Set(initialize=range(0, nb))#numero de nodos model.c=Set(initialize=range(0, nbr))#Numero de lineas model.OR = Set(initialize=range(0, NOR))#numero de ofertas de oportuniddad retiro model.OI = Set(initialize=range(0, NOI))#numero de ofertas de oportunidad inyeccion model.CFF=Set(initialize=range(0, NCFF))#numero de ofertas de CNFFF model.CF=Set(initialize=range(0, NCF))#numero de ofertas de contratos firmes model.s=Set(initialize=range(0, 5))#Numero de bloques #Parametros #Parametros de la red model.rtmw_min= Param(model.c, initialize=dict(enumerate(bl))) model.rtmw_max= Param(model.c, initialize=dict(enumerate(bu))) model.Inc = Param(model.c, model.i, initialize=arr2dict(inc)) model.Xc = Param(model.c, initialize=dict(enumerate(xc))) model.Rc = Param(model.c, initialize=dict(enumerate(rc))) #Parametros de los predespachos nacionales model.D = Param(model.i, initialize=dict(enumerate(dem))) model.G = Param(model.i, initialize=dict(enumerate(gen))) #Ofertas de oportunidad #Oferta de oportunidad de retiro model.fr= Param(model.OR, model.s, initialize=arr2dict(p_oor))#Oferta bloques 1 model.pr_ofertado = Param(model.OR, model.s, initialize=arr2dict(voor_m))#Magnitud de la oferta MW-h model.bin_pr = Param(model.i, model.OR, initialize=arr2dict(var_bin_oor)) #Oferta de oportunidad de inyeccion model.fi= Param(model.OI, model.s, initialize=arr2dict(p_ooi))#Precio de bloques - Oferta de oportunidad de inyeccion model.pi_ofertado= Param(model.OI, model.s, initialize=arr2dict(vooi_m))#Magnitud de la oferta MW-h model.bin_pi = Param(model.i,model.OI,initialize=arr2dict(var_bin_ooi)) #Contratos firmes model.pf_declarada=Param(model.CF, initialize=dict(enumerate(vcf_ed)))#Energia declarada model.pf_req=Param(model.CF, initialize=dict(enumerate(vcf_pr)))# Potencia requerida - Si no se flexbiliza deberian de ser igual la energia y la potencia #Precio de flexibilidad de contrato model.ffi=Param(model.CF, model.s, initialize=arr2dict(vcf_p))#Precio de bloques - Contrato firme - Oferta de flexibilidad model.pfi_ofertado=Param(model.CF, model.s, initialize=arr2dict(vcf_m))#Magnitud de la oferta - tiene que ser igual a la suma de la energia declarada model.bin_cfi=Param(model.i, model.CF,initialize=arr2dict(var_bin_cfi)) model.bin_cfr=Param(model.i, model.CF,initialize=arr2dict(var_bin_cfr)) #Ofertas de flexibilidad de contratos fisicos flexibles #Ofertas de inyeccion model.pff_declarada=Param(model.CFF, initialize=dict(enumerate(vnfff_ed))) model.pffi_ofertado=Param(model.CFF, model.s, initialize=arr2dict(vnfff_m_i))#Magnitud del bloque model.fffi=Param(model.CFF, model.s, initialize=arr2dict(p_cnfffi))#Precio de inyeccion model.bin_pffi=Param(model.i,model.CFF, initialize=arr2dict(var_bin_cnfffi)) #Oferta de retiro model.pffr_ofertado=Param(model.CFF, model.s, initialize=arr2dict(vnfff_m_r))#Magnitud de bloque de retiro ofertado model.fffr=Param(model.CFF, model.s, initialize=arr2dict(p_cnfffr))#Precio de bloques - Contrato no firme fisico flexible model.bin_pffr=Param(model.i,model.CFF, initialize=arr2dict(var_bin_cnfffr)) #Ofertad de pago maximo por CVT model.k=Param(model.CFF, initialize=dict(enumerate(k_cnfffcvt)))#Indicador de oferta model.pfft_ofertado=Param(model.CFF, model.s, initialize=arr2dict(vnfff_m_cvt))#Magnitud del bloque model.ffft=Param(model.CFF, model.s, initialize=arr2dict(p_cnfffcvt))#Precio de pago maximo CVT #Variabeles #Variable ofertas de oportunidad model.pr= Var(model.OR, model.s, domain=NonNegativeReals)#Parte aceptada de cada bloques de las ofertas de oportunidad de retiro model.pi= Var(model.OI, model.s, domain=NonNegativeReals)#Parte aceptada de cada bloques de las ofertas de oportunidad de inyeccion #Variables CF model.pfi=Var(model.CF, model.s, domain=NonNegativeReals)#Parte aceptada de cada bloques de las ofertas de flexibilidad de contratos firmes #Variables CNFFF model.pffr=Var(model.CFF, model.s, domain=NonNegativeReals)#Parte aceptada de cada bloques de las oferta de flexibilidad de retiro model.pffi=Var(model.CFF, model.s, domain=NonNegativeReals)#Parte aceptada de cada bloques de las ofertas de flexibilidad de inyección de los contratos físicos flexibles model.pfft=Var(model.CFF, model.s, domain=NonNegativeReals)#Parte aceptada de cada bloques de las oferta de pago máximo de CVT de los contratos físicos flexibles model.pff_iny_fisico=Var(model.CFF, domain=NonNegativeReals)#Componente fisica de energia horaria de inyeecion model.pff_ret_fisico=Var(model.CFF, domain=NonNegativeReals)#Componente fisica de energia horaria de retiro #Variables FOENS model.pf_cortada=Var(model.CF, domain=NonNegativeReals)#Energia firme de lo CF model.pf_pre_cortada=Var(model.CF, domain=NonNegativeReals) model.pff_cortada=Var(model.CFF, domain=NonNegativeReals)#Energia firme de los CNFFF model.fens=Var() #Variable problema de optimizacion model.inyeccion= Var(model.i, domain=NonNegativeReals)#Inyeccion por nodo model.retiro= Var(model.i, domain=NonNegativeReals)#Retiro por nodo model.ref_angular= Var(model.i)#Fase del voltaje en el nodo model.rtmw_c= Var(model.c)#Flujo de potencia actica por linea print("Ecuación de Función Objetivo Max") def objfunc(model): return (sum(model.fr[OR,s]*model.pr[OR,s] for OR in model.OR for s in model.s) - # ┌ FOO sum(model.fi[OI,s]*model.pi[OI,s] for OI in model.OI for s in model.s) - # └ sum(model.ffi[CF,s]*model.pfi[CF,s] for CF in model.CF for s in model.s) + # [ FOF sum(model.fffr[CFF,s]*model.pffr[CFF,s] for CFF in model.CFF for s in model.s) - # ┌ sum(model.fffi[CFF,s]*model.pffi[CFF,s] for CFF in model.CFF for s in model.s) + # │ FOFF sum(model.ffft[CFF,s]*model.pfft[CFF,s] for CFF in model.CFF for s in model.s) - # └ 0*model.fens*sum(model.pf_cortada[CF] for CF in model.CF) - # ┌ FOENS 0*model.fens*0.5*sum(model.pff_cortada[CFF] for CFF in model.CFF)) # └ model.OBJ= Objective(rule=objfunc, sense=maximize) print("Restricciones del Modelo de Optimización") def fens_restriccion(model): if NOR==0 & NCF==0 & NCFF>0: return model.fens == 3*max(model.fffr[CFF,s] for CFF in model.CFF for s in model.s) elif NOR==0 & NCF>0 & NCFF==0: return model.fens == 3*max(model.ffi[CF,s] for CF in model.CF for s in model.s) elif NOR==0 & NCF>0 & NCFF>0: return model.fens == 3*max(max(model.ffi[CF,s] for CF in model.CF for s in model.s), max(model.fffr[CFF,s] for CFF in model.CFF for s in model.s)) elif NOR>0 & NCF==0 & NCFF==0: return model.fens == 3*max(model.fr[OR,s] for OR in model.OR for s in model.s) elif NOR>0 & NCF==0 & NCFF>0: return model.fens == 3*max(max(model.fr[OR,s] for OR in model.OR for s in model.s), max(model.fffr[CFF,s] for CFF in model.CFF for s in model.s)) elif NOR>0 & NCF>0 & NCFF==0: return model.fens == 3*max(max(model.fr[OR,s] for OR in model.OR for s in model.s), max(model.ffi[CF,s] for CF in model.CF for s in model.s)) elif NOR>0 & NCF>0 & NCFF>0: return model.fens == 3*max(max(model.fr[OR,s] for OR in model.OR for s in model.s), max(model.fffr[CFF,s] for CFF in model.CFF for s in model.s), max(model.ffi[CF,s] for CF in model.CF for s in model.s)) model.fens_constraint= Constraint(rule=fens_restriccion) #Restrecciones FOO def pi_restriccion(model,OI,s): return ((model.pi[OI,s] <=model.pi_ofertado[OI,s])) model.pi_constraint= Constraint(model.OI, model.s, rule=pi_restriccion) def pr_restriccion(model,OR,s): return ((model.pr[OR,s]<=model.pr_ofertado[OR,s])) model.pr_constraint= Constraint(model.OR, model.s, rule=pr_restriccion) #Restricciones FOF def pfi_restriccion(model,CF,s): return (model.pfi[CF,s]<=model.pfi_ofertado[CF,s]) model.pfi_constraint= Constraint(model.CF, model.s, rule=pfi_restriccion) #Restricciones FOFF def pffr_restriccion(model,CFF,s): return (model.pffr[CFF,s]<=model.pffr_ofertado[CFF,s]) model.pffr_constraint= Constraint(model.CFF, model.s, rule=pffr_restriccion) def pffi_restriccion(model,CFF,s): return (model.pffi[CFF,s]<=model.pffi_ofertado[CFF,s]) model.pffi_constraint= Constraint(model.CFF, model.s, rule=pffi_restriccion) if (model.k[CFF] ==0 for CFF in model.CFF): def pfft_restriccion(model,CFF,s): #if model.k[CFF] ==0: return (model.pfft[CFF,s]<=model.pfft_ofertado[CFF,s]) model.pfft_constraint= Constraint(model.CFF, model.s, rule=pfft_restriccion) #K(cff) vale 0 si hay oferta de pago maximo por CVT def pff_iny_fisico_restriccion(model,CFF): if model.k[CFF] ==0: return (model.pff_iny_fisico[CFF]==sum(model.pfft[CFF,s] for s in model.s)- sum(model.pffr[CFF,s] for s in model.s)) elif model.k[CFF] ==1: return (model.pff_iny_fisico[CFF]==model.pff_declarada[CFF] - model.pff_cortada[CFF]- sum(model.pffr[CFF,s] for s in model.s)) model.pff_iny_fisico_constraint=Constraint(model.CFF, rule=pff_iny_fisico_restriccion) def pff_ret_fisico_restriccion(model,CFF): if model.k[CFF] ==0: return (model.pff_ret_fisico[CFF]==sum(model.pfft[CFF,s] for s in model.s) - sum(model.pffi[CFF,s] for s in model.s)) elif model.k[CFF] ==1: return (model.pff_ret_fisico[CFF]==model.pff_declarada[CFF] - model.pff_cortada[CFF]- sum(model.pffi[CFF,s] for s in model.s)) model.pff_ret_fisico_constraint=Constraint(model.CFF, rule=pff_ret_fisico_restriccion) #Restriccion FOENS def pff_cortada_restriccion(model,CFF): return (model.pff_cortada[CFF]<=model.pff_declarada[CFF]) model.pff_cortada_constraint=Constraint(model.CFF, rule=pff_cortada_restriccion) def pf_cortada_restriccion(model,CF): return (model.pf_cortada[CF]<=model.pf_req[CF] - model.pf_pre_cortada[CF]) model.pf_cortada_constraint=Constraint(model.CF, rule=pf_cortada_restriccion) print('Restricciones de transmision') #Restricciones de transmision def inyec(model,i): return (model.inyeccion[i] == model.G[i] + sum(model.pi[OI,s]*model.bin_pi[i,OI] for OI in model.OI for s in model.s) + sum(model.pfi[CF,s]*model.bin_cfi[i,CF] for CF in model.CF for s in model.s) + sum(model.pff_iny_fisico[CFF]*model.bin_pffi[i,CFF] for CFF in model.CFF)) model.inyec_constraint= Constraint(model.i, rule=inyec) def retiro(model,i): return (model.retiro[i] == model.D[i] + sum(model.pr[OR,s]*model.bin_pr[i,OR] for OR in model.OR for s in model.s) + sum(model.pf_req[CF]*model.bin_cfr[i,CF] for CF in model.CF) + sum(model.pff_ret_fisico[CFF]*model.bin_pffr[i,CFF] for CFF in model.CFF) - 0*sum(model.pf_cortada[CF]*model.bin_cfr[i,CF] for CF in model.CF) - 0*sum(model.pf_pre_cortada[CF]*model.bin_cfr[i,CF] for CF in model.CF)) model.retiro_constraint= Constraint(model.i, rule=retiro) #Los multiplicadore o variable duales de esta restriccion son los precios nodales def balance_inyeccion_retiro(model,i): return (model.inyeccion[i] + sum(model.Inc[c,i]*model.rtmw_c[c] for c in model.c) - # Al dividir entre 100 Rc 0.5*sum(((model.Inc[c,i]*model.rtmw_c[c])**2)*(model.Rc[c]/100) for c in model.c) # toda la ec. se pasa a p.u. == model.retiro[i]) model.balance_inyeccion_retiro_constraint= Constraint(model.i,rule=balance_inyeccion_retiro) def rtmw_min_restriccion(model,c): return (model.rtmw_c[c]>=-model.rtmw_min[c]) model.rtmw_min_constraint= Constraint(model.c, rule=rtmw_min_restriccion) def rtmw_max_restriccion(model,c): return (model.rtmw_c[c]<=model.rtmw_max[c]) model.rtmw_max_constraint= Constraint(model.c, rule=rtmw_max_restriccion) def flujo_potencia_actica(model,c): return ((model.Xc[c])*model.rtmw_c[c]+sum(model.Inc[c,i]*model.ref_angular[i] for i in model.i)== 0) model.flujo_potencia_actica_constraint= Constraint(model.c, rule=flujo_potencia_actica) model.dual = Suffix(direction=Suffix.IMPORT) print("Construcción del modelo terminada.") opt = SolverFactory('ipopt') #opt.options['max_iter']= 10000 result = opt.solve(model) model.solutions.store_to(result) # Cálculo de Precios Nodales # ============================================================================= print("Calculando Precios Nodales") Sigma = zeros(nb) for i in model.i: Sigma[i] = model.dual[model.balance_inyeccion_retiro_constraint[i]] # Construcción de array para grabar # ============================================================================= flujos = brnames.copy() f = array(list(model.rtmw_c.get_values().values())) perdidas=zeros(nbr) for c in model.c: perdidas[c]=(f[c]**2)*rc[c]/100 flujos['Flujo'] = f flujos['Perdidas'] = perdidas #print(flujos) pon = DataFrame() result_inyeccion = array(list(model.inyeccion.get_values().values())) result_retiro = array(list(model.retiro.get_values().values())) pon['nodo']= bus pon['Precio Exante'] = Sigma*-1 #print(pon) iep=DataFrame() iep['nodo']= bus iep['Inyeccion'] = result_inyeccion iep['Retiro'] = result_retiro result_pff_iny=array(list(model.pff_iny_fisico.get_values().values())) result_pff_ret=array(list(model.pff_ret_fisico.get_values().values())) result_pr = setvariable_p(array(list(model.pr.get_values().values())),NOR) result_pi = setvariable_p(array(list(model.pi.get_values().values())),NOI) result_pfi = setvariable_p(array(list(model.pfi.get_values().values())),NCF) result_pffr = setvariable_p(array(list(model.pffr.get_values().values())),NCFF) result_pffi = setvariable_p(array(list(model.pffi.get_values().values())),NCFF) result_pfft = setvariable_p(array(list(model.pfft.get_values().values())),NCFF) result_foo_r=DataFrame() result_foo_r['N°']=ex_oor['N°'] result_foo_r['Pr Bloque 1']=result_pr[:,0] result_foo_r['Pr Bloque 2']=result_pr[:,1] result_foo_r['Pr Bloque 3']=result_pr[:,2] result_foo_r['Pr Bloque 4']=result_pr[:,3] result_foo_r['Pr Bloque 5']=result_pr[:,4] result_foo_i=DataFrame() result_foo_i['N°']=ex_ooi['N°'] result_foo_i['Pi Bloque 1']=result_pi[:,0] result_foo_i['Pi Bloque 2']=result_pi[:,1] result_foo_i['Pi Bloque 3']=result_pi[:,2] result_foo_i['Pi Bloque 4']=result_pi[:,3] result_foo_i['Pi Bloque 5']=result_pi[:,4] #print(result_foo) result_fof=DataFrame() result_fof['N°']=ex_cf['N°'] result_fof['Pfi Bloque 1']=result_pfi[:,0] result_fof['Pfi Bloque 2']=result_pfi[:,1] result_fof['Pfi Bloque 3']=result_pfi[:,2] result_fof['Pfi Bloque 4']=result_pfi[:,3] result_fof['Pfi Bloque 5']=result_pfi[:,4] result_foff=DataFrame() result_foff['N°']=ex_cnfff['N°'] result_foff['Pffr Bloque 1']=result_pffr[:,0] result_foff['Pffr Bloque 2']=result_pffr[:,1] result_foff['Pffr Bloque 3']=result_pffr[:,2] result_foff['Pffr Bloque 4']=result_pffr[:,3] result_foff['Pffr Bloque 5']=result_pffr[:,4] result_foff['Pffi Bloque 1']=result_pffi[:,0] result_foff['Pffi Bloque 2']=result_pffi[:,1] result_foff['Pffi Bloque 3']=result_pffi[:,2] result_foff['Pffi Bloque 4']=result_pffi[:,3] result_foff['Pffi Bloque 5']=result_pffi[:,4] result_foff['Pfft Bloque 1']=result_pfft[:,0] result_foff['Pfft Bloque 2']=result_pfft[:,1] result_foff['Pfft Bloque 3']=result_pfft[:,2] result_foff['Pfft Bloque 4']=result_pfft[:,3] result_foff['Pfft Bloque 5']=result_pfft[:,4] #print(result_foff) foo_ret_iny=DataFrame() foo_ret_iny['N°']=ex_cnfff['N°'] foo_ret_iny['Inyeccion']=result_pff_iny foo_ret_iny['Retiro']=result_pff_ret print("Escribiendo resultados en carpeta") writer=ExcelWriter("Resultados_predespacho.xlsx") flujos.to_excel(writer,'IEP-RTR',index=False) iep.to_excel(writer,'IEP-TOTAL',index=False) pon.to_excel(writer,'PEXANTES',index=False) result_foo_i.to_excel(writer,'TOP-I',index=False) result_foo_r.to_excel(writer,'TOP-R',index=False) result_fof.to_excel(writer,'TCP-CF',index=False) result_foff.to_excel(writer,'TCP-CNFFF-1',index=False) foo_ret_iny.to_excel(writer,'TCP-CNFFF-2',index=False) writer.save() print("{:=^100}".format("")) print("{:^100}".format(" Script Finalizado ")) print("{:^100}".format(" Mercados Eléctricos de Centroamérica (c) 2020 ")) print("{:=^100}".format("")) #print(result) return 0