subastas-derechos-firmes/simsdt/model.py

# -*- coding: utf-8 -*-
"""
Modulo principal para ejecutar el modelo de subasta
"""

import logging
from logging.handlers import RotatingFileHandler
from os import path
import queue

import pandas as pd
import pyomo.environ as pe
import pyutilib.subprocess.GlobalData
from numpy import array, zeros
from openpyxl import load_workbook
from pyomo.environ import SolverFactory
from pyomo.kernel import value

from simsdt.mct.makeMCT import linmct, readmct, set_dir_flujo
from simsdt.ofertas.readBids import *
from simsdt.red.create import branchnames, setbranch, setbus
from simsdt.red.makeBdc import makeBdc
from simsdt.red.makePTDF import makePTDF
from simsdt.red.read import excel2net
from simsdt.utils.arr2dict import arr2dict
from simsdt.utils.idx_brch import RT_MAX, RT_MIN, R

from qhandler import QueueHandler

pyutilib.subprocess.GlobalData.DEFINE_SIGNAL_HANDLERS_DEFAULT = False


# log_queue = queue.Queue()
# qh = QueueHandler(log_queue)
# formatter = logging.Formatter(
#     '%(asctime)s: %(message)s', datefmt='%Y-%m-%dT%H:%M:%S')
# qh.setFormatter(formatter)

fh = RotatingFileHandler(
    'log/simsdt_model_monitor.log', maxBytes=10240, backupCount=10)
fmt = logging.Formatter(
    '%(asctime)s - %(name)s - %(levelname)s - %(message)s', datefmt='%Y-%m-%dT%H:%M:%S')
fh.setFormatter(fmt)
# ch = logging.StreamHandler()
# ch.setFormatter(fmt)

# logging.basicConfig(handlers=[fh, ch, qh], level=logging.INFO)
model_logger = logging.getLogger('simsdt.model')
# model_logger.addHandler(qh)
model_logger.addHandler(fh)


class ModeloSubasta:

    def __init__(self, file, solver_path):
        self.file = file
        self.solver_path = solver_path

    def setmodel(self):
        model_logger.info('Modelo de Subasta de Derechos Firmes')

        # ============================================================================
        # Importación de datos desde archivo Excel
        # ============================================================================

        model_logger.info("Inicio del Script de Subasta")
        model_logger.info("Leyendo información de red")

        net = excel2net(self.file)

        model_logger.info(
            "Leyendo información de ofertas y derechos firmes existentes")

        ex = readexistentes(self.file)
        nex = ex.shape[0]
        of = readofertas(self.file)
        nof = of.shape[0]

        # ============================================================================
        # Construcción de vectores y matrices para los SETS y PARAMETERS
        # ============================================================================

        model_logger.info("Creando información de buses y líneas")
        # Datos de Red
        bus = setbus(net)
        branch = setbranch(net, bus)

        nb = bus.shape[0]
        nbr = branch.shape[0]

        model_logger.debug(
            "Se han cargado {0} buses y {1} líneas".format(nb, nbr))
        model_logger.info("Leyendo información de Resistencia de líneas")

        brnames = branchnames(bus, branch)
        r = branch[:, R]

        # Datos de límites superior e inferior asociados a las lineas
        model_logger.info(
            "Obteniendo datos de límites superior e inferior de líneas...")
        bu = branch[:, RT_MAX]/100
        bl = branch[:, RT_MIN]/100

        # Lineas para límites de MCT
        lin_mct = linmct(brnames)
        mct = readmct(self.file)/100
        dirf = set_dir_flujo()

        model_logger.info("Cálculando la matriz H")
        # Cálculo de la matriz H
        H = makePTDF(bus, branch)

        model_logger.info("Cálculando la matriz de Incidencia")
        # Matriz de Incidencia
        _, _, A = makeBdc(bus, branch)
        inc = A.toarray()**2

        model_logger.info("Creando Vectores de Derechos Firmes existentes")
        # Datos de los Derechos Firmes Existentes
        vite = setVITE(bus, ex)/100
        vrte = setVRTE(bus, ex)/100
        vitex = setVITEX(bus, ex)/100
        te = setTE(vite, vrte)

        model_logger.info("Creando vectores de Ofertas de Derechos firmes")
        # Datos de las ofertas de Derechos Firmes
        c = of.oferta.values
        cper = of.cper.values
        perk = of.per.values/100
        vit = setVIT(bus, of)/100
        vrt = setVRT(bus, of)/100
        vitx = setVITX(bus, of)/100
        t = setT(vit, vrt)

        # ============================================================================
        # Modelo de optimización subasta de derechos Firmes
        # ============================================================================

        model_logger.info("Inicio del problema de Optimización")
        # Inicio del modelo de optimización

        model = pe.ConcreteModel()

        # SETS

        model_logger.info("Creando Sets del problema")
        model.c = pe.Set(initialize=range(0, nbr))   # Número de circuitos
        model.n = pe.Set(initialize=range(0, nb))    # Número de nodos
        model.o = pe.Set(initialize=range(0, nex))   # Número de DF existentes
        model.k = pe.Set(initialize=range(0, nof))   # Número de Ofertas DF
        model.i = pe.Set(initialize=lin_mct.keys())  # Interconexiones
        # Sentidos de interconexiones
        model.s = pe.Set(initialize=['sn', 'ns'])

        # PARAMETERS

        model_logger.info("Creando Parametros del problema")
        # Ofertas derechos firmes
        model_logger.info("Parametros de Ofertas de Derechos Firmes")
        model.C = pe.Param(model.k, initialize=dict(enumerate(c)))
        model.Cper = pe.Param(model.k, initialize=dict(enumerate(cper)))
        model.Perk = pe.Param(model.k, initialize=dict(enumerate(perk)))

        # Parametros de la red

        model_logger.info("Parametros de Red")
        model.Bu = pe.Param(model.c, initialize=dict(enumerate(bu)))
        model.Bl = pe.Param(model.c, initialize=dict(enumerate(bl)))
        model.R = pe.Param(model.c, initialize=dict(enumerate(r)))
        model.Inc = pe.Param(model.c, model.n, initialize=arr2dict(inc))

        # elemento H[c,n]

        model_logger.info("Parametros de Matriz H")
        model.He = pe.Param(model.c, model.n, initialize=arr2dict(H))

        model_logger.info("Parametros de Derechos Firmes Existentes")
        # Vector de inyecciones df existentes
        model.VITE = pe.Param(model.n, model.o, initialize=arr2dict(vite))
        # Vector de retiros df existentes
        model.VRTE = pe.Param(model.n, model.o, initialize=arr2dict(vrte))
        # Vector de perididas df existentes
        model.VITEX = pe.Param(model.n, model.o, initialize=arr2dict(vitex))
        # DF existentes
        model.TE = pe.Param(model.n, model.o, initialize=arr2dict(te))

        model_logger.info("Parametros de Inyecciones y Retiros de Ofertas")
        # Vector inyecciones ofertas de compra df
        model.VIT = pe.Param(model.n, model.k, initialize=arr2dict(vit))
        # Vector retiros ofertas de compra df
        model.VRT = pe.Param(model.n, model.k, initialize=arr2dict(vrt))
        # Vector inyecciones perdidas oferta de compra df
        model.VITX = pe.Param(model.n, model.k, initialize=arr2dict(vitx))
        # Ofertas de compora DF
        model.T = pe.Param(model.n, model.k, initialize=arr2dict(t))

        model_logger.info("Parametros de Limites de flujos en las líneas")

        # Lim. Superior flujos asociados a linea c
        def bfu_param(model, c):
            return model.Bu[c] - sum(max(0, sum(model.He[c, n]*model.TE[n, o] for n in model.n)) for o in model.o)

        model.BFu = pe.Param(model.c, initialize=bfu_param)

        # LIm. Inferior flujos asociados a linea c
        def bfl_param(model, c):
            return model.Bl[c] - sum(max(0, -sum(model.He[c, n]*model.TE[n, o] for n in model.n)) for o in model.o)

        model.BFl = pe.Param(model.c, initialize=bfl_param)

        # Máximas Capacidades de Trasferencia
        model.Mct = pe.Param(model.i, model.s, initialize=mct.mct.to_dict())
        model.DirF = pe.Param(model.i, initialize=dirf)

        # VARIABLES

        model_logger.info("Creando Variables del problema")
        # Variables primales

        model.alpha_k = pe.Var(
            model.k, within=pe.NonNegativeReals, bounds=(0, 1))
        model.gamma_k = pe.Var(model.k, within=pe.NonNegativeReals)

        # Flujos

        model.F = pe.Var(model.c)
        model.PL = pe.Var(model.c)
        model.Per = pe.Var(model.n)
        model.F_k = pe.Var(model.c, model.k)
        model.Fp_k = pe.Var(model.c, model.k, within=pe.NonNegativeReals)
        model.Fp = pe.Var(model.c, within=pe.NonNegativeReals)
        model.Fn_k = pe.Var(model.c, model.k, within=pe.NonNegativeReals)
        model.Fn = pe.Var(model.c, within=pe.NonNegativeReals)

        # EQUATIONS

        model_logger.info("Creando Ecuaciones del problema")
        # Ecuación Objetivo
        # =============================================================================

        model_logger.info("Ecuación de Función Objetivo Max Z")

        def eq_Z(model):
            return sum(model.alpha_k[k]*model.C[k] - model.gamma_k[k]*model.Cper[k] for k in model.k)

        model.Z = pe.Objective(rule=eq_Z, sense=pe.maximize)

        # Restricciones
        # =============================================================================

        model_logger.info("Restricciones del Modelo de Optimización")
        model_logger.info("Límites de Gamma")

        def eq_gamma_k_rule(model, k):
            return model.gamma_k[k] <= model.alpha_k[k]

        model.eq_gamma_k = pe.Constraint(model.k, rule=eq_gamma_k_rule)

        # Factibilidad de derechos firmes

        model_logger.info("Factibilidad de derechos Firmes")

        def eq_F_k_rule(model, c, k):
            return model.F_k[c, k] == model.alpha_k[k]*sum(model.He[c, n]*model.T[n, k] for n in model.n)

        model.eq_F_k = pe.Constraint(model.c, model.k, rule=eq_F_k_rule)

        # Flujos positivos

        model_logger.info("Flujos positivos")

        def eq_Fp_k_rule(model, c, k):
            return model.Fp_k[c, k] >= model.F_k[c, k]

        model.eq_Fp_k = pe.Constraint(model.c, model.k, rule=eq_Fp_k_rule)

        def eq_Fp_rule(model, c):
            return model.Fp[c] == sum(model.Fp_k[c, k] for k in model.k)

        model.eq_Fp = pe.Constraint(model.c, rule=eq_Fp_rule)

        def eq_Fpl_rule(model, c):
            return sum(model.Fp_k[c, k] for k in model.k) <= model.BFu[c]

        model.eq_Fpl = pe.Constraint(model.c, rule=eq_Fpl_rule)

        # Flujos negativos

        model_logger.info("Flujos negativos")

        def eq_Fn_k_rule(model, c, k):
            return model.Fn_k[c, k] >= -model.F_k[c, k]

        model.eq_Fn_k = pe.Constraint(model.c, model.k, rule=eq_Fn_k_rule)

        def eq_Fn_rule(model, c):
            return model.Fn[c] == sum(model.Fn_k[c, k] for k in model.k)

        model.eq_Fn = pe.Constraint(model.c, rule=eq_Fn_rule)

        def eq_Fnl_rule(model, c):
            return sum(model.Fn_k[c, k] for k in model.k) <= model.BFl[c]

        model.eq_Fnl = pe.Constraint(model.c, rule=eq_Fnl_rule)

        # Suficiencia Financiera

        model_logger.info("Suficiencia Financiera")

        def eq_sf1_rule(model, c):
            return model.F[c] <= model.Bu[c]

        model.eq_sf1 = pe.Constraint(model.c, rule=eq_sf1_rule)

        def eq_sf2_rule(model, c):
            return model.F[c] >= -model.Bl[c]

        model.eq_sf2 = pe.Constraint(model.c, rule=eq_sf2_rule)

        # Balance de Pérdidas

        model_logger.info("Balance de Pérdidas")

        def eq_PL_rule(model, c):
            return model.PL[c] == model.R[c]*(model.F[c]*model.F[c])

        model.eq_PL = pe.Constraint(model.c, rule=eq_PL_rule)

        def eq_Per_rule(model, n):
            return model.Per[n] == sum(model.Inc[c, n]*(model.PL[c]/2) for c in model.c)

        model.eq_Per = pe.Constraint(model.n, rule=eq_Per_rule)

        # Balance de Flujos

        model_logger.info("Balance de Flujos")

        def eq_Balance_rule(model, c):
            return model.F[c] == sum(model.He[c, n]*(
                sum(model.alpha_k[k]*model.T[n, k] for k in model.k) +
                sum(model.TE[n, o] for o in model.o) +
                sum(model.gamma_k[k]*model.VITX[n, k] for k in model.k) +
                sum(model.VITEX[n, o]for o in model.o) -
                model.Per[n]) for n in model.n)

        model.eq_Balance = pe.Constraint(model.c, rule=eq_Balance_rule)

        # Perdidas en estado base

        model_logger.info("Perdidas en estado Base")

        def eq_perdidas_base_rule(model):
            return sum(model.PL[c] for c in model.c) == sum(
                sum(model.gamma_k[k]*model.VITX[n, k] for k in model.k) +
                sum(model.VITEX[n, o] for o in model.o) for n in model.n)

        model.eq_perdidas_base = pe.Constraint(rule=eq_perdidas_base_rule)

        # Máximas capacidades de trasferencia entre áreas de control

        model_logger.info(
            "Máxima Capacidad de Trasnferencia entre áreas de control")

        def eq_mct_ns_rule(model, i):
            return sum(model.DirF[i]*model.F[c] for c in lin_mct[i]) <= \
                model.Mct[i, 'ns']

        model.eq_mct_ns = pe.Constraint(model.i, rule=eq_mct_ns_rule)

        def eq_mct_sn_rule(model, i):
            return sum(model.DirF[i]*model.F[c] for c in lin_mct[i]) >= \
                -model.Mct[i, 'sn']

        model.eq_mct_sn = pe.Constraint(model.i, rule=eq_mct_sn_rule)

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

        model_logger.info("Construcción del modelo terminada.")

        filename = path.basename(self.file).split('.')
        filepath = path.dirname(self.file)

        try:
            model_logger.info("Enviando modelo algebraico al solver")
            opt = SolverFactory(
                'ipopt',
                executable=self.solver_path
            )

            model_logger.info("Buscando una solución optima al problema")
            result = opt.solve(
                model, tee=True, logfile="log/solver_{}.log".format(filename[0]))

            model.solutions.store_to(result)
            result.write(filename='results.json', format='json')
            model_logger.info("Solucion encontrada")

        except Exception as e:
            model_logger.info("No se ejecutó el problema.")
            model_logger.error("Error: {}".format(e))
            return

        # Cálculo de Precios Nodales
        # =============================================================================
        model_logger.info("Calculando Precios Nodales")

        Sigma = zeros(nbr)
        for c in model.c:
            Sigma[c] = model.dual[model.eq_Balance[c]]

        Lambda = model.dual[model.eq_perdidas_base]

        PON = (H.T @ Sigma) + Lambda

        # Cálculo de Pagos por Derechos Firmes
        # =============================================================================
        model_logger.info("Calculando Pagos Asignados por Derechos Firmes")

        alpha = array(list(model.alpha_k.get_values().values()))
        psi = array(list(model.gamma_k.get_values().values()))

        PDF_k = zeros(nof)

        for i in range(0, nof):
            PDF_k[i] = -PON.T @ (alpha[i]*t[:, i] + psi[i]*vitx[:, i])

        # Construcción de array flujos en interconexiones
        # =============================================================================

        intercon = {}

        for i in model.i:
            intercon[i] = [pe.value(model.eq_mct_sn[i].lower),
                           pe.value(model.eq_mct_ns[i].body),
                           pe.value(model.eq_mct_ns[i].upper)]

        flujos_inter = pd.DataFrame.from_dict(intercon)
        flujos_inter = flujos_inter * 100

        # Construcción de array para grabar
        # =============================================================================
        filename = path.basename(self.file).split('.')
        filepath = path.dirname(self.file)

        model_logger.info(
            "Escribiendo resultados en carpeta {0}".format(filepath))

        # Resultados de Asignación

        resultado = of.copy()

        resultado['alpha'] = alpha
        resultado['gamma'] = psi
        resultado['MW_asign'] = resultado['MWo'] * resultado['alpha']
        resultado['per_asign'] = resultado['per'] * resultado['gamma']
        resultado['monto_asign'] = PDF_k.round(2)
        resultado['precio_mw'] = resultado['monto_asign'] / \
            (resultado['MW_asign'] * 24 * 31)  # Se toman meses de 31 días

        # Flujos de Potencia

        flujos = net.copy()
        flujos = flujos.drop(['x', 'r', 'max', 'min'], axis=1)
        f = array(list(model.F.get_values().values()))
        flujos['linea'] = brnames
        flujos['flujo'] = f
        flujos['perdidas'] = array(list(model.PL.get_values().values()))
        flujos['sigma'] = Sigma

        # Precios por Nodo

        pon = bus.copy()
        pon['PON'] = PON

        # Grabar en excel
        # ================================================================================

        # new_filename = filename[0]+'_results.'+filename[1]
        #
        # new_file = path.join(filepath, new_filename)
        #
        # book = load_workbook(new_file)
        # writer = ExcelWriter(new_file, engine='openpyxl')
        # writer.book = book
        # writer.sheets = dict((ws.title, ws) for ws in book.worksheets)
        #
        # resultado.to_excel(writer, sheet_name='asignacion', index=False, header=True)
        #
        # flujos.to_excel(writer, sheet_name='flujos', index=False, header=True)
        #
        # pon.to_excel(writer, sheet_name='pon', index=False, header=True)
        #
        # writer.save()
        # writer.close()

        # Grabar en csv
        # =================================================================================
        csv_path = filepath
        model_logger.info("Escribiendo archivo de asignaciones")
        resultado.to_csv(
            path.join(csv_path, filename[0] + '_asignacion.csv'), sep=',', index=False)

        model_logger.info("Escribiendo archivo de flujos de potencia")
        flujos.to_csv(
            path.join(csv_path, filename[0] + '_flujos.csv'), sep=',', index=False)

        model_logger.info("Escribiendo archivo de flujos de precios")
        pon.to_csv(
            path.join(csv_path, filename[0] + '_pon.csv'), sep=',', index=False)

        model_logger.info(
            "Escribiendo archivo de flujos de potencia entre areas de control")
        flujos_inter.to_csv(
            path.join(csv_path, filename[0] + '_flujos_inter.csv'), sep=',', index=False)

        model_logger.info("Script Finalizado ")