Commit inicial

.gitignore

@@ -0,0 +1,142 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+.hypothesis/
+.pytest_cache/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# pyenv
+.python-version
+
+# celery beat schedule file
+celerybeat-schedule
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+
+
+#PyCharm
+### JetBrains template
+# Covers JetBrains IDEs: IntelliJ, RubyMine, PhpStorm, AppCode, PyCharm, CLion, Android Studio and Webstorm
+# Reference: https://intellij-support.jetbrains.com/hc/en-us/articles/206544839
+
+# User-specific stuff:
+.idea/workspace.xml
+.idea/tasks.xml
+.idea/dictionaries
+.idea/vcs.xml
+.idea/jsLibraryMappings.xml
+
+# Sensitive or high-churn files:
+.idea/dataSources.ids
+.idea/dataSources.xml
+.idea/dataSources.local.xml
+.idea/sqlDataSources.xml
+.idea/dynamic.xml
+.idea/uiDesigner.xml
+
+# Gradle:
+.idea/gradle.xml
+.idea/libraries
+
+# Mongo Explorer plugin:
+.idea/mongoSettings.xml
+
+.idea/
+
+# VSCode
+.vscode
+
+# User
+*.xls
+*.csv
+results.json

README.md

@@ -0,0 +1,3 @@
+# Modelo de Subasta de Derechos Firmes
+
+Copyright (c) Mercados Electricos de Centroamerica (MERELEC)

common/log.py

@@ -0,0 +1,3 @@
+#
+# Utility classes for working with the logger
+#

main.py

@@ -0,0 +1,294 @@
+import datetime
+import logging
+import queue
+import signal
+import threading
+import time
+import tkinter as tk
+from logging.handlers import RotatingFileHandler
+from tkinter import (HORIZONTAL, VERTICAL, E, N, S, W, filedialog, messagebox,
+                     ttk)
+from tkinter.scrolledtext import ScrolledText
+
+from PIL import Image, ImageTk
+
+from qhandler import QueueHandler
+from simsdt import runsdt
+
+
+fh = RotatingFileHandler(
+    'log/simsdt_app_monitor.log', maxBytes=1024*50, backupCount=10)
+fmt = logging.Formatter(
+    '%(asctime)s - %(name)s - %(levelname)s - %(message)s', datefmt='%Y-%m-%dT%H:%M:%S')
+fh.setFormatter(fmt)
+fh.setLevel(logging.INFO)
+# logging.basicConfig(format='%(asctime)s - %(name)s - %(levelname)s - %(message)s',
+#                     datefmt='%Y-%m-%dT%H:%M:%S',
+#                     handlers=[RotatingFileHandler(
+#                         'log/simsdt_app_monitor.log', maxBytes=10240, backupCount=10)],
+#                     level=logging.INFO)
+
+main_logger = logging.getLogger("simsdt")
+main_logger.addHandler(fh)
+main_logger.setLevel(logging.DEBUG)
+
+
+class RunsdtThread(threading.Thread):
+    def __init__(self, group=None, target=None, name=None,
+                 args=(), kwargs=None, *, daemon=None):
+        super(RunsdtThread, self).__init__(group=group,
+                                           target=target, name=name, daemon=daemon)
+        self._stop_event = threading.Event()
+        self.args = args
+        self.kwargs = kwargs
+        return
+
+    def run(self):
+        now = datetime.datetime.now()
+        runsdt.main(self.args)
+        # time.sleep(30)
+        main_logger.info('Fin de la función run')
+
+    def stop(self):
+        self._stop_event.set()
+
+
+class Clock(threading.Thread):
+    """Class to display the time every seconds
+
+    Every 5 seconds, the time is displayed using the logging.ERROR level
+    to show that different colors are associated to the log levels
+    """
+
+    def __init__(self):
+        super().__init__()
+        self._stop_event = threading.Event()
+
+    def run(self):
+        # main_logger.setLevel(logging.INFO)
+        main_logger.debug('Clock started')
+        previous = -1
+        while not self._stop_event.is_set():
+            now = datetime.datetime.now()
+            if previous != now.second:
+                previous = now.second
+                if now.second % 5 == 0:
+                    level = logging.ERROR
+                else:
+                    level = logging.INFO
+                main_logger.log(level, now)
+            time.sleep(0.2)
+
+    def stop(self):
+        self._stop_event.set()
+
+
+class ConsoleUi:
+    """Poll messages from a logging queue and display them in a scrolled text widget"""
+
+    def __init__(self, frame):
+        self.frame = frame
+        # Create a ScrolledText wdiget
+        self.scrolled_text = ScrolledText(frame, state='disabled', height=25)
+        self.scrolled_text.grid(row=0, column=0, sticky=(N, S, W, E))
+        # Font config
+        self.scrolled_text.configure(font='TkFixedFont')
+        self.scrolled_text.tag_config('INFO', foreground='black')
+        self.scrolled_text.tag_config('DEBUG', foreground='gray')
+        self.scrolled_text.tag_config('WARNING', foreground='orange')
+        self.scrolled_text.tag_config('ERROR', foreground='red')
+        self.scrolled_text.tag_config(
+            'CRITICAL', foreground='red', underline=1)
+        # Create a logging handler using a queue
+        self.log_queue = queue.Queue()
+        self.queue_handler = QueueHandler(self.log_queue)
+        formatter = logging.Formatter(
+            '%(asctime)s: %(message)s', datefmt='%Y-%m-%dT%H:%M:%S')
+        self.queue_handler.setFormatter(formatter)
+        # self.queue_handler.setLevel(logging.INFO)
+        main_logger.addHandler(self.queue_handler)
+        # Start polling messages from the queue
+        self.frame.after(100, self.poll_log_queue)
+
+    def display(self, record):
+        msg = self.queue_handler.format(record)
+        self.scrolled_text.configure(state='normal')
+        self.scrolled_text.insert(tk.END, msg + '\n', record.levelname)
+        self.scrolled_text.configure(state='disabled')
+        # Autoscroll to the bottom
+        self.scrolled_text.yview(tk.END)
+
+    def poll_log_queue(self):
+        # Check every 100ms if there is a new message in the queue to display
+        while True:
+            try:
+                record = self.log_queue.get(block=False)
+            except queue.Empty:
+                break
+            else:
+                self.display(record)
+        self.frame.after(100, self.poll_log_queue)
+
+
+class FormUi:
+
+    runsdt_thread = None
+
+    def __init__(self, frame, q):
+        self.frame = frame
+        self.q = q
+        image = Image.open('img/play.png')
+        self.photo = ImageTk.PhotoImage(image)
+
+        # Add a text field for file path
+        self.file_path = tk.StringVar()
+        ttk.Label(self.frame, text='Archivo:').grid(column=0, row=0, sticky=W)
+        ttk.Entry(self.frame, textvariable=self.file_path,
+                  width=50).grid(column=1, row=0, sticky=(W, E))
+        self.button = ttk.Button(
+            self.frame, text='...', command=self.open_file_dialog, width=3)
+        self.button.grid(column=2, row=0, sticky=W)
+        self.button_run = ttk.Button(self.frame, image=self.photo, compound="right", text='Ejecutar',
+                                     command=self.run)
+        self.button_run.grid(column=1, row=1, sticky=W)
+
+    def open_file_dialog(self):
+        # Open file dialog
+        file = filedialog.askopenfilename(
+            filetypes=[('Excel', '.xlsx')], title='Seleccione un archivo de subasta...')
+        main_logger.info(f'Archivo seleccionado: {file}')
+        self.file_path.set(file)
+
+    def run(self):
+        if not self.file_path.get():
+            main_logger.error('Debe seleccionar un archivo para ejecutarlo')
+            return
+
+        self.button_run.configure(state='disabled')
+        main_logger.info('Ejecutando')
+        # self.runsdt_thread = threading.Thread(
+        #     target=runsdt.main, args=[self.file_path.get()])
+        # self.runsdt_thread.start()
+        self.runsdt_thread = RunsdtThread(args=(self.file_path.get()))
+        self.runsdt_thread.start()
+
+        self.check_thread()
+
+    def check_thread(self):
+
+        if not self.runsdt_thread.is_alive():
+            self.button_run.configure(state='normal')
+            self.q.put(False)
+        else:
+            self.q.put(True)
+            self.frame.after(100, self.check_thread)
+
+
+class ThirdUi:
+
+    def __init__(self, frame, q):
+        self.q = q
+        self.frame = frame
+        self.pb_max = 100
+        # ttk.Label(self.frame, text='This is just an example of a third frame').grid(
+        #     column=0, row=1, sticky=W)
+        self.progress_bar = ttk.Progressbar(
+            self.frame, mode='determinate')
+        self.progress_bar.pack(expand=True, fill=tk.BOTH, side=tk.TOP)
+
+        self.frame.after(100, self.check_q)
+
+    def check_q(self):
+        # Check every 100ms if there is a new message in the queue to display
+        while True:
+            try:
+                record = self.q.get(block=False)
+            except queue.Empty:
+                break
+            else:
+                if record:
+                    step = 1
+                    self.progress_bar.step(step)
+                    if step >= self.pb_max:
+                        step = 1
+                    else:
+                        step += 1
+                else:
+                    self.progress_bar.config(value=0)
+        self.frame.after(100, self.check_q)
+
+        # max_size = 100
+        # print(self.q.qsize())
+        # while self.q.qsize():
+        #     try:
+        #         active_tread = self.q.get(0)
+        #         print(active_tread)
+        #         main_logger.info(active_tread)
+        #         if active_tread:
+        #             step = 1
+        #             self.progress_bar.step(step)
+        #             if step == max_size:
+        #                 step = 1
+        #             else:
+        #                 step += 1
+        #     except queue.Empty:
+        #         pass
+
+        #     self.frame.after(100, self.check_q)
+
+
+class App:
+
+    def __init__(self, root):
+        self.root = root
+        pbqueue = queue.Queue()
+
+        root.title('SimSDT')
+        root.iconbitmap("app.ico")
+        root.columnconfigure(0, weight=1)
+        root.rowconfigure(0, weight=1)
+        # Create the panes and frames
+        vertical_pane = ttk.PanedWindow(self.root, orient=VERTICAL)
+        vertical_pane.grid(row=0, column=0, sticky="nsew")
+        horizontal_pane = ttk.PanedWindow(vertical_pane, orient=HORIZONTAL)
+        vertical_pane.add(horizontal_pane)
+        form_frame = ttk.Labelframe(horizontal_pane, text="Configuracion")
+        form_frame.columnconfigure(1, weight=1)
+        horizontal_pane.add(form_frame, weight=1)
+        console_frame = ttk.Labelframe(horizontal_pane, text="Console")
+        console_frame.columnconfigure(0, weight=1)
+        console_frame.rowconfigure(0, weight=1)
+        horizontal_pane.add(console_frame, weight=1)
+        third_frame = ttk.Labelframe(vertical_pane, text="Estado")
+        vertical_pane.add(third_frame, weight=1)
+        # Initialize all frames
+        self.form = FormUi(form_frame, pbqueue)
+        self.console = ConsoleUi(console_frame)
+        self.third = ThirdUi(third_frame, pbqueue)
+        self.root.protocol('WM_DELETE_WINDOW', self.quit)
+        self.root.bind('<Control-q>', self.quit)
+        signal.signal(signal.SIGINT, self.quit)
+
+    def quit(self, *args):
+
+        msg = 'Se esta ejecutando un proceso de optimización\n\n¿Desea salir de la aplicación?'
+
+        if self.form.runsdt_thread and self.form.runsdt_thread.is_alive():
+            if messagebox.askokcancel("Salir de SimSDT", msg):
+                main_logger.info('Cerrando aplicacion SimSDT')
+                self.form.runsdt_thread.stop()
+                self.root.destroy()
+        else:
+            main_logger.info('Saliendo de la aplicacion')
+            self.root.destroy()
+
+
+def main():
+    main_logger.info('SimSDT ha iniciado')
+    root = tk.Tk()
+    app = App(root)
+    app.root.mainloop()
+
+
+if __name__ == '__main__':
+    main()

qhandler/__init__.py

@@ -0,0 +1 @@
+from .queuehandler import QueueHandler

qhandler/queuehandler.py

@@ -0,0 +1,19 @@
+import logging
+
+
+class QueueHandler(logging.Handler):
+    """Class to send logging records to a queue
+
+    It can be used from different threads
+    The ConsoleUi class polls this queue to display records in a ScrolledText widget
+    """
+    # Example from Moshe Kaplan: https://gist.github.com/moshekaplan/c425f861de7bbf28ef06
+    # (https://stackoverflow.com/questions/13318742/python-logging-to-tkinter-text-widget) is not thread safe!
+    # See https://stackoverflow.com/questions/43909849/tkinter-python-crashes-on-new-thread-trying-to-log-on-main-thread
+
+    def __init__(self, log_queue):
+        super().__init__()
+        self.log_queue = log_queue
+
+    def emit(self, record):
+        self.log_queue.put(record)

requirements.txt

@@ -0,0 +1,26 @@
+altgraph==0.17
+appdirs==1.4.4
+autopep8==1.5.2
+click==7.1.2
+click-spinner==0.1.10
+et-xmlfile==1.0.1
+future==0.18.2
+jdcal==1.4.1
+nose==1.3.7
+numpy==1.18.4
+openpyxl==3.0.3
+pandas==1.0.3
+pefile==2019.4.18
+Pillow==7.1.2
+ply==3.11
+pycodestyle==2.6.0
+pygubu==0.9.8.6
+PyInstaller==3.6
+Pyomo==5.6.9
+python-dateutil==2.8.1
+pytz==2020.1
+PyUtilib==5.8.0
+pywin32-ctypes==0.2.0
+scipy==1.4.1
+six==1.15.0
+xlrd==1.2.0

simsdt/__init__.py

@@ -0,0 +1,6 @@
+"""
+Modulo de Subastas de Derechos Firmes.
+
+.. moduleauthor:: Oscar A. Leiva
+
+"""

simsdt/mct/makeMCT.py

@@ -0,0 +1,69 @@
+# -*- coding: utf-8 -*-
+# =============================================================================
+#  Copyright (C) 2018 Mercados Electricos de Centroamérica. All rights reserved
+# =============================================================================
+
+"""Crea limites de trasnferencia entre área de Control 
+"""
+from numpy import argwhere
+from pandas import read_excel
+
+# Se definen las líneas que componenn las interconexiones entre las areas de
+# control
+
+__GUAESA = ['1126-29162-1', '1124-29161-1']
+__GUAHON = ['1710-3190-1']
+__ESAHON = ['28181-29181-1', '28181-29182-2']
+__HONNIC = ['3301-4411-1', '3310-4407-1']
+__NICCRC = ['4412-50050-1', '4408-50000-1']
+__CRCPAN = ['6500-56050-1', '6400-58350-1', '6000-56050-1']
+
+__INTER = {'GUAESA': __GUAESA,
+           'GUAHON': __GUAHON,
+           'ESAHON': __ESAHON,
+           'HONNIC': __HONNIC,
+           'NICCRC': __NICCRC,
+           'CRCPAN': __CRCPAN}
+
+
+def linmct(brnames):
+    """Construye un diccionario con el indice que corresponde al nombre del
+    circuito o linea que copone la interconexión entre las áreas de control
+    """
+    lmct = {}
+    
+    for item in __INTER:
+        lmct[item] = []
+        for i in __INTER[item]:
+            l = int(argwhere(brnames == i))
+            lmct[item].append(l)
+    
+    return lmct
+
+
+def readmct(file):
+    """Lee la información de Máximas Capacidades de Trasnferencia del archivo
+    Excel del caso a resolver.
+    """
+    mct = read_excel(file, 'mct', index_col=[0, 1])
+    
+    return mct
+
+
+def set_dir_flujo():
+    """Consruye un diccionario en el cual se define la dirección del circuito
+    que compone la interconexión.
+    
+    En el caso de  la interconexión entre Costa Rica y Panamá, las líneas se
+    encuentran definidas en el sentido Sur-Norte, por lo que su dirección de
+    flujo corresponde a -1
+    """
+    d = {}
+    
+    for key in __INTER.keys():
+        if key == 'CRCPAN':
+            d[key] = -1
+        else:
+            d[key] = 1
+            
+    return d

simsdt/model.py

@@ -0,0 +1,516 @@
+# -*- 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 ")

simsdt/ofertas/readBids.py

@@ -0,0 +1,245 @@
+# -*- coding: utf-8 -*-
+# =============================================================================
+#  Copyright (C) 2018 Mercados Electricos de Centroamérica. All rights reserved
+# =============================================================================
+
+"""Importa la información de inyecciones y retiros de archivo de Excel 
+"""
+
+from pandas import read_excel, concat
+from numpy import zeros
+
+
+def readexistentes(subasta_file):
+    """Lee del archivo Excel la información de Derechos Firmes existente.
+    
+    Retorna un DataFrame con las columnas siguientes:
+        
+        1. cod_agente    : Código del Agente propietario del DF existente.
+        2. tke           : identificador del DF existente
+        3. bus_ie        : Nodo de inyección del DF existente
+        4. bus_je        : Nodo de retiro del DF existente
+        5. MWe           : Potencia asignada al DF existente
+    """
+    exist = read_excel(subasta_file,sheet_name='existentes')
+    exist.columns = ['cod_agente','tke','bus_ie','bus_je','MWe','per']
+    
+    return exist
+
+
+def readofertas(subasta_file):
+    """Lee del archivo Excel la información de Ofertad de Derechos Firmes.
+    
+    Retorna un DataFrame con las columnas siguientes:
+        
+        1. cod_agente    : Código del Agente que oferta.
+        2. tko           : identificador de la oferta
+        3. bus_io        : Nodo de inyección de la oferta
+        4. bus_jo        : Nodo de retiro de la oferta
+        5. MWo           : Potencia ofertada
+        6. oferta        : Oferta en dólares por la compra de DF
+        7. precio        : Precio del MWh ofertado [oferta/(24*365*MWo)]
+        8. per           : Perdidas asignadas a la oferta
+        9. cper          : Oferta en dolares por las pérdidas
+    """
+    ofer = read_excel(subasta_file,sheet_name='oferta')
+    ofer.columns = ['cod_agente','tko','bus_io','bus_jo','MWo','oferta','precio','per','cper']
+    
+    return ofer
+
+
+def setVITE(bus, ex):
+    """Construye el vector VITE de inyecciones de Derechos Firmes existentes
+    
+    Vector cuya dimensión es C{nb}, en donde c{nb} es el número de nodos, y 
+    cuyas componentes son todas nulas, salvo la correspondiente al nodo de 
+    inyección del derecho firme donde tiene el valor correspondiene a la 
+    inyección asociada al derecho
+    
+    La función regresa una matriz de dimensión C{nb x ne}, en donde C{nb} es el
+    número de nodos y C{ne} es el número de derechos firmes existentes. Cada
+    columna de la matriz es el vector VITEe de cada derecho firme e.
+    
+    """
+    
+    _vite = zeros((bus.shape[0],ex.shape[0]))
+    
+    for i in range(0, ex.shape[0]):
+        _vite[bus[bus == ex.iloc[i].bus_ie].index[0],i] = ex.iloc[i].MWe
+        
+    return _vite
+
+
+def setVRTE(bus, ex):
+    """Construye el vector VRTE de retiros de Derechos Firmes existentes
+    
+    Vector cuya dimensión es C{nb}, en donde nb es el número de nodos, y cuyas
+    componentes son todas nulas, salvo la correspondiente al nodo de retiro
+    del derecho firme donde tiene el valor correspondiene a la inyección
+    asociada al derecho
+    
+    La función regresa una matriz de dimensión C{nb x ne}, en donde C{nb} es el
+    número de nodos y C{ne} es el número de derechos firmes existentes. Cada
+    columna de la matriz es el vector VRTEe de cada derecho firme e.
+    
+    """
+    
+    _vrte = zeros((bus.shape[0], ex.shape[0]))
+    
+    for i in range(0, ex.shape[0]):
+        _vrte[bus[bus == ex.iloc[i].bus_je].index[0],i] = ex.iloc[i].MWe
+        
+    return _vrte
+
+
+def setTE(vite,vrte):
+    """Construye el vector TE resta VITE - VRTE
+    
+    La suma de su componentes es nula.
+    
+    La función regresa una matriz de dimensión C{nb x ne}, en donde C{nb} es el
+    número de nodos y C{ne} es el número de derechos firmes existentes. Cada
+    columna de la matriz es el vector TEe de cada derecho firme e.
+    """
+    _te = vite - vrte
+    
+    return _te
+
+
+def setVITEX(bus, ex):
+    """Construye el vector VITEX de perdidas de Derechos Firmes existentes
+    
+    Vector cuya dimensión es C{nb}, en donde c{nb} es el número de nodos, y 
+    cuyas componentes son todas nulas, salvo la correspondiente al nodo de 
+    inyección del derecho firme donde tiene el valor correspondiene a la 
+    inyección asociada al derecho
+    
+    La función regresa una matriz de dimensión C{nb x ne}, en donde C{nb} es el
+    número de nodos y C{ne} es el número de derechos firmes existentes. Cada
+    columna de la matriz es el vector VITEXe de cada derecho firme e.
+    
+    """
+    
+    _vitex = zeros((bus.shape[0],ex.shape[0]))
+    
+    for i in range(0, ex.shape[0]):
+        _vitex[bus[bus == ex.iloc[i].bus_ie].index[0],i] = ex.iloc[i].per
+        
+    return _vitex
+
+
+def setVIT(bus, of):
+    """Construye el vector VIT de inyecciones de Ofertas de Derechos Firmes
+    
+    Vector cuya dimensión es C{nb}, en donde nb es el número de nodos, y cuyas
+    componentes son todas nulas, salvo la correspondiente al nodo de inyeccion
+    del derecho firme donde tiene el valor correspondiene a la inyección
+    ofertada en el derecho k
+    
+    La función regresa una matriz de dimensión C{nb x nk}, en donde C{nb} es el
+    número de nodos y C{nk} es el número de ofertas de derechos firmes. Cada
+    columna de la matriz es el vector VITk de cada oferta k.
+    
+    """
+    
+    _vit = zeros((bus.shape[0], of.shape[0]))
+    
+    for i in range(0, of.shape[0]):
+        _vit[bus[bus == of.iloc[i].bus_io].index[0],i] = of.iloc[i].MWo
+        
+    return _vit
+
+
+def setVRT(bus, of):
+    """Construye el vector VRT de retiros de Ofertas de Derechos Firmes
+    
+    Vector cuya dimensión es C{nb}, en donde nb es el número de nodos, y cuyas
+    componentes son todas nulas, salvo la correspondiente al nodo de retiro
+    del derecho firme donde tiene el valor correspondiene a la inyección
+    ofertada en el derecho k
+    
+    La función regresa una matriz de dimensión C{nb x nk}, en donde C{nb} es el
+    número de nodos y C{nk} es el número de ofertas de derechos firmes. Cada
+    columna de la matriz es el vector VRTk de cada oferta k.
+    
+    """
+   
+    _vrt = zeros((bus.shape[0], of.shape[0]))
+    
+    for i in range(0, of.shape[0]):
+        _vrt[bus[bus == of.iloc[i].bus_jo].index[0],i] = of.iloc[i].MWo
+        
+    return _vrt
+
+
+def setT(vit,vrt):
+    """Construye el vector T resta VIT - VRT
+    
+    La suma de sus componentes es nula.
+    
+    La función regresa una matriz de dimensión C{nb x nk}, en donde C{nb} es el
+    número de nodos y C{nk} es el número de ofertas de derechos firmes. Cada
+    columna de la matriz es el vector Tk de cada oferta k.
+    """
+    _t = vit - vrt
+    
+    return _t
+
+
+def setVITX(bus,of):
+    """Construye el vector VITX de perdidas de Ofertas de Derechos Firmes
+    
+    Vector cuya dimensión es C{nb}, en donde nb es el número de nodos, y cuyas
+    componentes son todas nulas, salvo la correspondiente al nodo de inyección
+    de la oferta de derecho firme donde tiene el valor correspondiene a la 
+    perdida máxima asociada a la oferta k
+    
+    La función regresa una matriz de dimensión C{nb x nk}, en donde C{nb} es el
+    número de nodos y C{nk} es el número de ofertas de derechos firmes. Cada
+    columna de la matriz es el vector VITXk de cada oferta k.
+    
+    """
+    _vitx = zeros((bus.shape[0], of.shape[0]))
+    
+    for i in range(0, of.shape[0]):
+        _vitx[bus[bus == of.iloc[i].bus_io].index[0],i] = of.iloc[i].per
+        
+    return _vitx
+
+
+def setIny(bus, exist, ofer):
+    """Construye el vector de inyecciones para todos los nodos de la red
+    """
+    
+    iny_e = exist.pivot_table(values='MWe', aggfunc='sum', columns='bus_ie').T
+    iny_o = ofer.pivot_table(values = 'MWo', aggfunc='sum', columns='bus_io').T
+    
+    iny=concat([iny_e, iny_o], axis=1).fillna(0)
+    
+    iny['total'] = iny.iloc[:,0]+iny.iloc[:,1]
+    
+    _iny = zeros(bus.shape[0])
+    
+    for i in range(0,iny.shape[0]):
+        _iny[bus[bus == iny.iloc[i].name].index[0]] = iny.iloc[i].total
+    
+    return _iny
+
+
+def setRet(bus, exist, ofer):
+    """ Construye el vector de retiros para todos los nodos de la red
+    """
+    
+    ret_e = exist.pivot_table(values='MWe', aggfunc='sum', columns='bus_je').T
+    ret_o = ofer.pivot_table(values = 'MWo', aggfunc='sum', columns='bus_jo').T
+    
+    ret=concat([ret_e, ret_o], axis=1).fillna(0)
+    
+    ret['total'] = ret.iloc[:,0]+ret.iloc[:,1]
+    
+    _ret = zeros(bus.shape[0])
+    
+    for i in range(0, ret.shape[0]):
+        _ret[bus[bus == ret.iloc[i].name].index[0]] = ret.iloc[i].total
+    
+    return _ret

simsdt/red/create.py

@@ -0,0 +1,46 @@
+# -*- coding: utf-8 -*-
+# =============================================================================
+#  Copyright (C) 2018 Mercados Electricos de Centroamérica. All rights reserved
+# =============================================================================
+
+"""Crea las estructuras con numeración secuencial de buses, circuitos y nombres de circuitos
+"""
+from numpy import array
+from simsdt.utils.idx_brch import BUS_I, BUS_J, CKT
+
+
+def setbus(net):
+    """Lee la información de la red y crea una serie con los códigos de nodos y su numeración correlativa.
+    """
+    b = net['bus_i'].append(net['bus_j']).drop_duplicates().sort_values()
+    b = b.reset_index(drop=True)
+
+    return b
+
+
+def setbranch(net, b):
+    """Utiliza la información de la red, lod códigos de nodos y su numeración correlativa para generar un listado con
+    de los circuitos de acuerdo a la numeración correlativo del nodo.
+    """
+    br = net.copy()
+
+    for i in range(0, net.shape[0]):
+        br.loc[i, 'bus_i'] = b[b == br['bus_i'][i]].index[0]
+        br.loc[i, 'bus_j'] = b[b == br['bus_j'][i]].index[0]
+
+    br = br.sort_values(['bus_i', 'bus_j'])
+
+    return br.values
+
+
+def branchnames(b, br):
+    """Devuelve un array con los nombres de los circuitos.
+    """
+    brnames = []
+    for i in range(0, br.shape[0]):
+        brnames.append("-".join(
+            [str(b.iloc[br[i, BUS_I]]),
+             str(b.iloc[br[i, BUS_J]]),
+             str(br[i, CKT])]))
+    brnames = array(brnames)
+    return brnames

simsdt/red/makeBdc.py

@@ -0,0 +1,69 @@
+# -*- coding: utf-8 -*-
+# =============================================================================
+#  Copyright (C) 2018 Mercados Electricos de Centroamérica. All rights reserved
+# =============================================================================
+
+"""
+Construye las matrices B para el cálculo de flujos DC
+"""
+
+from numpy import ones, r_, zeros
+from scipy.sparse import csr_matrix as sparse
+from simsdt.utils.idx_brch import BUS_I, BUS_J, X
+
+
+def makeBdc(bus, branch):
+    """Construye las matrices B para el cálculo de flujos DC.
+
+    Se construye primer la matriz de incidencia Cft (o matriz A, como se    denomina en la formulación de substas de
+    Derechos Firmes) y se define como:
+
+    El valor será uno (1) para las líneas ij donde su nodo inicial sea igual al nodo i y menos uno (-1) para las lineas
+    ij cuyo nodo final sea igual al nodo j.
+
+    La matriz de incidencias Cft tiene dimenciones C{nbr x nb}, en donde C{nbr}es el numero de líneas (branches) y
+    C{nb} el número de buses.
+
+    Luego se calcula la matriz Bf (ZZ) de dimensiones C{nbr x nb}, en donde C{nbr} es el numero de líneas (branches) y
+    C{nb} el número de buses, igual que la Matriz de Incidencias.
+
+    Contiene los mismos elementos nulos que la matriz Cft y cuyas componentes no nulas son:
+
+        1. Para la componente correspondiente a la fila de cada circuito br,y cuya columna corresponde al nodo de
+        llegada del circuito b
+
+            :math:`Bf[br,b] = 1/Xbr`
+
+        2. Para la componente correspondiente a la fila de cada circuito br, y cuya columna corresponde al nodo de
+        salida del circuito b
+
+            :math:`Bf[br,b] = -1/Xbr`
+
+    La matriz Bbus se calcula como la transpuesta de la matriz de incidencias Cft por la matriz Bf
+
+        :math:`Bbus = Cft.T * Bf`
+
+    .. codeauthor:: Oscar A. Leiva (MERELEC)
+    """
+    # constantes
+    nb = bus.shape[0]
+    nbr = branch.shape[0]
+
+    # construye la matriz de incidencia Cft
+    f = branch[:, BUS_I].astype(float)          # lista de buses de origen
+    t = branch[:, BUS_J].astype(float)          # lista de buses de destino
+    i = r_[range(nbr), range(nbr)]             # Set doble de indices de filas
+
+    # Matriz de incidencia
+    Cft = sparse((r_[ones(nbr), -ones(nbr)], (i, r_[f, t])), (nbr, nb))
+
+    # calculo de las susceptancias en serie
+    b = 1 / branch[:, X].astype(float)          # Susceptancias en serie
+
+    # construye Bf (ZZ en la formulación de Subastas de Derechos Firmes)
+    Bf = sparse((r_[b, -b], (i, r_[f, t])), shape=(nbr, nb))
+
+    # construye la matriz Bbus
+    Bbus = Cft.T * Bf
+
+    return Bbus, Bf, Cft

simsdt/red/makePTDF.py

@@ -0,0 +1,39 @@
+# -*- coding: utf-8 -*-
+# =============================================================================
+#  Copyright (C) 2018 Mercados Electricos de Centroamérica. All rights reserved
+# =============================================================================
+
+"""Construye las matrices PTDF (Matriz H) para el cálculo de flujos DC
+"""
+
+from numpy import zeros
+from numpy.linalg import solve
+from simsdt.red.makeBdc import makeBdc
+
+
+def makePTDF(bus, branch):
+    """Construye la matriz PTDF DC
+
+    Toma como nodo slack el primer nodo.
+
+    La matriz es de dimensiones C{nbr x nb}, en donde C{nbr} es el numero de
+    lineas y C{nb} es el numero de buses o nodos. El nodo slack siempre es el
+    primer nodo.
+
+    @author: Oscar A. Leiva (MERELEC)
+    """
+
+    # constantes
+    nb = bus.shape[0]
+    nbr = branch.shape[0]
+
+    Bbus, Bf, _ = makeBdc(bus, branch)
+
+    Bbus, Bf = Bbus.todense(), Bf.todense()
+
+    H = zeros((nbr, nb))
+
+    H[:, 1:] = solve(Bbus[1:, 1:].T, Bf[:, 1:].T).T
+    #       = Bf[:,1:] * inv(Bbus[1:,1:])
+
+    return H

simsdt/red/read.py

@@ -0,0 +1,51 @@
+# -*- coding: utf-8 -*-
+# =============================================================================
+#  Copyright (C) 2018 Mercados Electricos de Centroamérica. All rights reserved
+# =============================================================================
+
+"""Importa el archivo de red RTR a utilizar para el cálculo 
+"""
+from pandas import read_sql, read_csv, read_excel
+
+# SERVER = "192.168.98.134"
+# DB = "MEIntegrado"
+# conn = pyodbc.connect('DRIVER={SQL Server};SERVER=' + SERVER + ';DATABASE=' + DB)
+# query = "SELECT BUS_I, BUS_J, ID_CKT, X, R, CTOIJ, CTOJI  FROM MATRIZ_RTR WHERE FECHA = '{f}'"
+
+# logger = logging.getLogger('simsdt.red')
+
+
+def sql2net(fecha):
+    """Lee la información de  la red desde una base de datos
+
+    .. warning::
+
+        Función no implementada.
+    """
+    q = query.format(f=fecha)
+    df = read_sql(q, conn)
+    df.columns = ['bus_i', 'bus_j', 'ckt', 'x', 'r', 'max', 'min']
+
+    return df
+
+
+def csv2net(file):
+    """Lee la información de  la red desde un archivo csv
+    """
+    df = read_csv(file)
+    df.columns = ['bus_i', 'bus_j', 'ckt', 'x', 'r', 'max', 'min']
+
+    return df
+
+
+def excel2net(file):
+    """Lee la información de  la red desde un archivo de excel.
+
+    .. note::
+
+        La información debe de estar en una hoja del archivo llamada 'rtr'.
+    """
+
+    df = read_excel(file, sheet_name='rtr')
+    df.columns = ['bus_i', 'bus_j', 'ckt', 'x', 'r', 'max', 'min']
+    return df

simsdt/runsdt.py

@@ -0,0 +1,69 @@
+# -*- coding: utf-8 -*-
+"""
+Modulo principal para ejecutar el modelo de subasta
+"""
+import logging
+from logging.handlers import RotatingFileHandler
+import queue
+import shutil
+import sys
+
+
+from simsdt.model import ModeloSubasta
+from qhandler import QueueHandler
+
+# 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_run_monitor.log', maxBytes=1024*50, 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)
+
+run_logger = logging.getLogger('simsdt.run')
+run_logger.setLevel(logging.INFO)
+run_logger.addHandler(fh)
+# run_logger.addHandler(qh)
+# run_logger.addHandler(ch)
+
+
+def main(file, solver_path='', q=None):
+    """Ejecuta la Simulación de Subastas de Derechos de Transmisión."""
+
+    if not file:
+        run_logger.error('No se incluyo un archivo para la ejecucion')
+        return
+
+    # Check if solver 'ipopt' is on path
+    run_logger.info('Verificando que exista el solver ipopt')
+
+    solver_path = shutil.which('ipopt')
+
+    if(solver_path):
+        run_logger.info('Ubicación del solver: {}'.format(solver_path))
+    else:
+        run_logger.error('El solver no se encuentra en el PATH')
+        return
+
+    run_logger.info("Inicio de la ejecucion del modelo de subasta")
+
+    modelo_subasta = ModeloSubasta(file, solver_path)
+    try:
+        modelo_subasta.setmodel()
+    except Exception as e:
+        run_logger.error(e)
+        return
+
+
+if __name__ == "__main__":
+    if len(sys.argv) < 2:
+        print("Se debe incluir un archivo")
+
+    else:
+        main(sys.argv[1])

simsdt/utils/__init__.py

@@ -0,0 +1,2 @@
+"""Funciones varias para la ejecución del modelo
+"""

simsdt/utils/arr2dict.py

@@ -0,0 +1,16 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jul  5 16:01:56 2018
+
+@author: BI4
+"""
+
+
+def arr2dict(x):
+    """Convierte un array en diccionario."""
+    d = {}
+    for i in range(0, x.shape[0]):
+        for j in range(0, x.shape[1]):
+            d[i, j] = x[i, j]
+            
+    return d

simsdt/utils/idx_brch.py

@@ -0,0 +1,26 @@
+# -*- coding: utf-8 -*-
+# =============================================================================
+#  Copyright (C) 2018 Mercados Electricos de Centroamérica. All rights reserved
+# =============================================================================
+
+"""Define constantes para nombres de las columnas del dataframe de lineas
+
+El indice, nombre y descripción de cada columna es el siguiente:
+
+Columnas 0-5 incluidas en la definición de la red
+    0. C{BUS_I}     Bus de origen
+    1. C{BUS_J}     Bus de destino
+    2. C{X}         Reactancia (p.u)
+    3. C{R}         Resistencia (p.u)
+    4. C{RT_MAX}    Limite superior de transferencia en MW ij
+    5. C{RT_MIN}    Limite inferior de transferencia en MW ji
+"""
+
+# definición de indices
+BUS_I    = 0
+BUS_J    = 1
+CKT      = 2
+X        = 3
+R        = 4
+RT_MAX   = 5
+RT_MIN   = 6

simsdt/utils/map.py

@@ -0,0 +1,63 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Jul 11 09:28:46 2018
+
+@author: BI4
+"""
+
+from mpl_toolkits.basemap import Basemap
+import matplotlib.pyplot as plt
+import matplotlib.patches as mpatches
+
+def plotMCT(MCT):
+    
+    plt.figure(figsize=(24,12))
+    _map = Basemap(projection='merc',
+                  resolution='l',
+                  llcrnrlon=-93, llcrnrlat=7,
+                  urcrnrlon=-75, urcrnrlat=19)
+    
+    _map.drawcoastlines()
+    _map.drawcountries(linewidth=1)
+    _map.drawmapboundary(fill_color='#99ffff')
+    _map.fillcontinents(color='#cc9966',lake_color='#99ffff')
+    
+    
+    _lons = [-90.25,-88.9167,-86.5,-85,-84,-81]
+    _lats = [15.5,13.7,15,13,10,8.2]
+    _inter = {'GUAESA':[0,1,'{0} MW'],
+              'GUAHON':[0,2,'{0} MW'],
+              'ESAHON':[1,2,'{0} MW'],
+              'HONNIC':[2,3,'{0} MW'],
+              'NICCRC':[3,4,'{0} MW'],
+              'CRCPAN':[4,5,'{0} MW']}
+    
+    _x,_y = _map(_lons,_lats)
+    
+    for i in _inter:
+        a = [_x[_inter[i][0]],_x[_inter[i][1]]]
+        b = [_y[_inter[i][0]],_y[_inter[i][1]]]
+        
+        if MCT[i][1] >= MCT[i][2] or MCT[i][1] <= MCT[i][0]:
+            color = 'orangered'
+        elif MCT[i][1] >= 0.8*MCT[i][2] or MCT[i][1] <= 0.8*MCT[i][0]:
+            color = 'yellow'
+        else:
+            color = 'green'
+        
+        plt.plot(a,b,'-',color=color,linewidth=3, label=_inter[i][2])
+    
+        plt.text(sum(a)/2,sum(b)/2,'{:.2f} MW'.format(MCT[i][1]*100), 
+                 fontdict={'fontsize': 12, 'backgroundcolor':'#cc9966', 'fontweight':'bold'},
+                 horizontalalignment='center',
+                 verticalalignment='center',)
+    
+    _map.plot(_x,_y,'o', color='navy',markersize=18)
+    
+    red_patch = mpatches.Patch(color='red', label='Interconexión al 100% de capacidad')
+    yellow_patch = mpatches.Patch(color='yellow', label='Interconexión a mas del 80% de capacidad')
+    green_patch = mpatches.Patch(color='green',label='Interconexión a menos del 80% de capacidad')
+    plt.legend(handles=[red_patch,yellow_patch,green_patch], loc=3, fontsize=16)
+    
+    plt.suptitle('Flujos entre Áreas de Control', fontsize =18, fontweight='bold')
+    plt.show()

simsdt/utils/results.py

@@ -0,0 +1,39 @@
+import json
+
+
+class SolverResult:
+
+    _data = dict()
+    status = None
+    time_elapsed = None
+    condition = None
+    constraints = None
+    objectives = None
+    variables = None
+    solutions = None
+
+    def __init__(self):
+        pass
+
+    def read(self, file):
+
+        with open(file) as f:
+            self._data = json.load(f)
+
+            solver = self._data['Solver'][0]
+            problem = self._data['Problem'][0]
+            solution = self._data['Solution'][0]
+
+            self.status = solver['Status']
+            self.time_elapsed = solver['Time']
+            self.condition = solver['Termination condition']
+            self.status = solver['Status']
+            self.constraints = problem['Number of constraints']
+            self.objectives = problem['Number of objectives']
+            self.variables = problem['Number of variables']
+            self.solutions = solution['number of solutions']
+
+
+sr = SolverResult()
+sr.read('results2.json')
+print(sr.time_elapsed, sr.status, sr.condition)