Modelo_de_Predespacho/Predespacho_24horas/predespacho_24h.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
from mct.makeMCT import linmct, readmct, set_dir_flujo
from os import path

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

    # Lineas para límites de MCT
    br_t = brnames['Total'].to_numpy()
    lin_mct = linmct(br_t)
    dirf = set_dir_flujo()
    mct = readmct(file)

    #Informacion de los despachos nacionales
    dg_t =read_D_G(file)

    #Ofertas de todos los periodos
    ex_cnfff_t = readofertas_cnfffs(file)
    ex_ooi_t = readofertas_oois(file)
    ex_oor_t = readofertas_oors(file)
    ex_cf_t = readofertas_cfs(file)

    #Dataframe para almacenar la informacion de cada predespacho
    flujos_t=DataFrame()
    iep_t=DataFrame()
    pon_t = DataFrame()
    result_foo_i_t = DataFrame()
    result_foo_r_t = DataFrame()
    result_fof_t = DataFrame()
    result_foff_t = DataFrame()
    foo_ret_iny_t = DataFrame()
    #El range del for se modifica dependiendo los periodos que se quieren correr.
    for PERIODO in range(0,24):

        print("Leyendo información de los despacho nacionales. Periodo: " + str(PERIODO))
        dg=dg_t[dg_t.periodo.isin([PERIODO])]#Dependiendo de los valores de la columna periodo del df son los valores que toma 
        dg.reset_index(drop=True, inplace=True)#Resetea valores de ordenamiento del DF
        dg_n=set_dgnacional(bus,dg)#Ordena los valores que tiene segun los nodos de rtr

        print("Leyendo información de ofertas. Periodo: " + str(PERIODO))   
        #Funcion leer contratos no firmes fisicos flexibles
        ex_cnfff=ex_cnfff_t[ex_cnfff_t.periodo.isin([PERIODO])]#Dependiendo de los valores de la columna periodo del df son los valores que toma 
        ex_cnfff.reset_index(drop=True, inplace=True)#Resetea valores de ordenamiento del DF
        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_n['Demanda'])
        gen=setvariable(dg_n['Generacion'])

        #Funcion leer parametros oferta de oportunidad inyeccion
        ex_ooi = ex_ooi_t[ex_ooi_t.periodo.isin([PERIODO])]#Dependiendo de los valores de la columna periodo del df son los valores que toma 
        ex_ooi.reset_index(drop=True, inplace=True)#Resetea valores de ordenamiento del DF
        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=ex_oor_t[ex_oor_t.periodo.isin([PERIODO])]
        ex_oor.reset_index(drop=True, inplace=True)#Resetea valores de ordenamiento del DF
        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=ex_cf_t[ex_cf_t.periodo.isin([PERIODO])]#Dependiendo de los valores de la columna periodo del df son los valores que toma 
        ex_cf.reset_index(drop=True, inplace=True)#Resetea valores de ordenamiento del DF
        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
        model.inter = Set(initialize=lin_mct.keys())  # Interconexiones
        model.sen = Set(initialize=['sn', 'ns'])    # Sentidos de interconexiones
        model.per=Set(initialize=range(0, 24))#numero de periodos

        #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 PERIODOria de inyeecion 
        model.pff_ret_fisico=Var(model.CFF, domain=NonNegativeReals)#Componente fisica de energia PERIODOria 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

        # Máximas Capacidades de Trasferencia
        model.Mct = Param(model.inter, model.sen, model.per, initialize=mct.mct.to_dict())#Maxima transferenica(interconexion,sentido,periodo)
        model.DirF = Param(model.inter, initialize=dirf)

        print("Ecuación de Función Objetivo Max. Periodo: " + str(PERIODO))

        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. Periodo: "+str(PERIODO))

        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] - 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_cortada_constraint=Constraint(model.CF, rule=pf_cortada_restriccion)

        print('Restricciones de transmision. Periodo: ' + str(PERIODO))
        #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)

        # Máximas capacidades de trasferencia entre áreas de control
        def eq_mct_ns_rule(model, inter):
            if inter=='ES':#Maxima exportacion de El Salvador
                return (-(model.DirF[inter]*model.rtmw_c[lin_mct[inter][0]]+model.DirF[inter]*model.rtmw_c[lin_mct[inter][1]]) +
                (model.DirF[inter]*model.rtmw_c[lin_mct[inter][2]]+model.DirF[inter]*model.rtmw_c[lin_mct[inter][3]])) <= model.Mct[inter, 'ns', PERIODO] + 0.001
            elif inter=='HON':#Maxima importacion de Honduras
                return (-(model.DirF[inter]*model.rtmw_c[lin_mct[inter][0]]+model.DirF[inter]*model.rtmw_c[lin_mct[inter][1]] +
                model.DirF[inter]*model.rtmw_c[lin_mct[inter][2]]) + (model.DirF[inter]*model.rtmw_c[lin_mct[inter][3]]+
                model.DirF[inter]*model.rtmw_c[lin_mct[inter][4]])) <= model.Mct[inter, 'sn', PERIODO] + 0.001
            else:
                return sum(model.DirF[inter]*model.rtmw_c[c] for c in lin_mct[inter]) <= model.Mct[inter, 'ns', PERIODO] + 0.001
        model.eq_mct_ns = Constraint(model.inter, rule=eq_mct_ns_rule)

        def eq_mct_sn_rule(model, inter):
            if inter=='ES':#Maxima importacion de El Salvador
                return (-(model.DirF[inter]*model.rtmw_c[lin_mct[inter][0]]+model.DirF[inter]*model.rtmw_c[lin_mct[inter][1]]) +
                (model.DirF[inter]*model.rtmw_c[lin_mct[inter][2]]+model.DirF[inter]*model.rtmw_c[lin_mct[inter][3]])) >= -model.Mct[inter, 'sn', PERIODO] + 0.00001
            elif inter=='HON':#Maxima importacion de Honduras
                return (-(model.DirF[inter]*model.rtmw_c[lin_mct[inter][0]]+model.DirF[inter]*model.rtmw_c[lin_mct[inter][1]] +
                model.DirF[inter]*model.rtmw_c[lin_mct[inter][2]]) + (model.DirF[inter]*model.rtmw_c[lin_mct[inter][3]]+
                model.DirF[inter]*model.rtmw_c[lin_mct[inter][4]])) >= -model.Mct[inter, 'ns', PERIODO] + 0.00001
            else:
                return sum(model.DirF[inter]*model.rtmw_c[c] for c in lin_mct[inter]) >= -model.Mct[inter, 'sn', PERIODO] + 0.00001
        model.eq_mct_sn = Constraint(model.inter, rule=eq_mct_sn_rule)

        model.dual = Suffix(direction=Suffix.IMPORT)

        print("Construcción del modelo terminada. Periodo: "+str(PERIODO))

        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. Periodo:" + str(PERIODO))
        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()
        flujos['Periodo']=set_periodo(nbr,PERIODO)
        flujos['BUS I'] = brnames['BUS I']
        flujos['BUS J'] = brnames['BUS J']
        flujos['CKT'] = brnames['CKT']
        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
        flujos_t=concat([flujos_t,flujos])#Guarda los resultados con los anteriores

        pon = DataFrame()
        result_inyeccion = array(list(model.inyeccion.get_values().values()))
        result_retiro = array(list(model.retiro.get_values().values()))
        pon['Periodo']=set_periodo(nb,PERIODO)
        pon['nodo']= bus
        pon['Precio Exante'] = Sigma*-1
        pon_t=concat([pon_t,pon])#Guarda los resultados con los anteriores
        #print(pon)
        
        iep=DataFrame()
        iep['Periodo']=set_periodo(nb,PERIODO)
        iep['nodo']= bus
        iep['Inyeccion'] = result_inyeccion
        iep['Retiro'] = result_retiro
        iep_t=concat([iep_t,iep])#Guarda los resultados con los anteriores 

        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['Periodo']=ex_oor['periodo']
        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_r_t=concat([result_foo_r_t,result_foo_r])#Guarda los resultados con los anteriores

        result_foo_i=DataFrame()
        result_foo_i['Periodo']=ex_ooi['periodo']
        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]
        result_foo_i_t=concat([result_foo_i_t,result_foo_i])#Guarda los resultados con los anteriores
        
        result_fof=DataFrame()
        result_fof['Periodo']=ex_cf['periodo']
        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_fof_t=concat([result_fof_t,result_fof])#Guarda los resultados con los anteriores

        result_foff=DataFrame()
        result_foff['Periodo']=ex_cnfff['periodo']
        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]
        result_foff_t=concat([result_foff_t,result_foff])#Guarda los resultados con los anteriores

        foo_ret_iny=DataFrame()
        foo_ret_iny['Periodo']=ex_cnfff['periodo']
        foo_ret_iny['N°']=ex_cnfff['N°']
        foo_ret_iny['Inyeccion']=result_pff_iny
        foo_ret_iny['Retiro']=result_pff_ret
        foo_ret_iny_t=concat([foo_ret_iny_t,foo_ret_iny])#Guarda los resultados con los anteriores

    print("Escribiendo resultados en carpeta.")
    filename = path.basename(file).split('.')
    filepath = path.dirname(file)
    new_filename = filename[0]+'_results.'+filename[1]
    new_file = path.join(filepath, new_filename)
    writer=ExcelWriter(new_file)
    flujos_t.to_excel(writer,'IEP-RTR',index=False)
    iep_t.to_excel(writer,'IEP-TOTAL',index=False)
    pon_t.to_excel(writer,'PEXANTES',index=False)
    result_foo_i_t.to_excel(writer,'TOP-I',index=False)
    result_foo_r_t.to_excel(writer,'TOP-R',index=False)
    result_fof_t.to_excel(writer,'TCP-CF',index=False)
    result_foff_t.to_excel(writer,'TCP-CNFFF-1',index=False)
    foo_ret_iny_t.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(""))
        
    return 0