Modelo_de_Predespacho/Predespacho.py

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
import xlsxwriter

def setmodel(file):
        
    print("{:=^100}".format(""))
    print("{:^100}".format(" Modelo de Predespacho"))
    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
    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)

    #Calculo de parametros
    if NOR>0 & NCFF>0:
        fens=3*max(max(ex_oor['precio_oor1']),max(ex_cnfff['precio_i1']))#Valor de los retiro firmes no suministrados
    elif NOR>0 or NCFF==0:
        fens=3*max(ex_oor['precio_oor1'])
    elif NOR==0 or NCFF>0:
        fens=3*max(ex_cnfff['precio_i1'])
    ffenes=0.5*fens#Valor de la energía no flexible de los contratos físicos flexibles no suministrados 

    # 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)))
    
    #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
    
    #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) -  # └
        fens*sum(model.pf_cortada[CF] for CF in model.CF) -                               # ┌ FOENS
        ffenes*sum(model.pff_cortada[CFF] for CFF in model.CFF))                          # └
    model.OBJ= Objective(rule=objfunc, sense=maximize)

    print("Restricciones del Modelo de Optimización")

    #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) - 
        sum(model.pf_cortada[CF]*model.bin_cfr[i,CF] for CF in model.CF) -
        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)== 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)

    #Referencia angular
    #def ref_angular_restriccion(model,i):
    #    return model.ref_angular[0]==0
    #model.ref_angular_constraint =Constraint(model.i, rule=ref_angular_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('cplex')
    result = opt.solve(model)
    model.solutions.store_to(result)
    #print(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 = DataFrame()
    f = array(list(model.rtmw_c.get_values().values()))
    flujos['linea'] = brnames
    flujos['flujo'] = f
    #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['Inyeccion'] = result_inyeccion
    pon['Retiro'] = result_retiro
    pon['Ex-antes'] = Sigma*-1
    #print(pon)
    
    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_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=DataFrame()
    result_foo['Pr Bloque 1']=result_pr[:,0]
    result_foo['Pr Bloque 2']=result_pr[:,1]
    result_foo['Pr Bloque 3']=result_pr[:,2]
    result_foo['Pr Bloque 4']=result_pr[:,3]
    result_foo['Pr Bloque 5']=result_pr[:,4]
    result_foo['Pi Bloque 1']=result_pi[:,0]
    result_foo['Pi Bloque 2']=result_pi[:,1]
    result_foo['Pi Bloque 3']=result_pi[:,2]
    result_foo['Pi Bloque 4']=result_pi[:,3]
    result_foo['Pi Bloque 5']=result_pi[:,4]
    #print(result_foo)
    
    result_foff=DataFrame()
    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['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,'flujos',index=False)
    pon.to_excel(writer,'pon',index=False)
    result_foo.to_excel(writer,'result_foo',index=False)
    result_foff.to_excel(writer,'result_foff',index=False)
    foo_ret_iny.to_excel(writer,'result_foo_ret_iny',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(""))
    
    return 0