# -*- 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 ")