universidad de chile facultad de ciencias físicas y matemáticas

Loading...

UNIVERSIDAD DE CHILE FACULTAD DE CIENCIAS FÍSICAS Y MATEMÁTICAS DEPARTAMENTO DE INGENIERÍA INDUSTRIAL

OPTIMIZACIÓN LINEAL ENTERA MIXTA APLICADA A PROBLEMAS DE PLANIFICACIÓN ESTRATÉGICA EN ELECTRICIDAD

TESIS PARA OPTAR AL GRADO DE DOCTOR EN SISTEMAS DE INGENIERÍA

ALEJANDRO ALBERTO ANGULO CÁRDENAS

PROFESOR GUÍA: DANIEL ESPINOZA GONZÁLEZ PROFESOR CO-GUÍA: RODRIGO PALMA BEHNKE

MIEMBROS DE LA COMISIÓN: EDUARDO MORENO ARAYA FERNANDO ORDOÑEZ PIZARRO JOSÉ CORREA HAEUSSLER

Esta investigación fue financiada por becas CONICYT del Gobierno de Chile y el Instituto Milenio de Sistemas Complejos de Ingeniería

SANTIAGO DE CHILE 2015

ii

Resumen En esta tesis se presentan los resultados del trabajo desarrollado por el autor durante el periodo en que fue estudiante de doctorado en el Departamento de Industrias de la Universidad de Chile. El trabajo se centra en la aplicación de técnicas de optimización entera-mixtas a problemas de planificación estratégica del sector eléctrico, donde el problema de corto plazo correspondiente al predespacho de unidades de generación en sistemas térmicos es el tema central en estudio. En lo relativo al modelamiento del problema de predespacho de unidades, se considera el análisis de las distintas formulaciones entera-mixtas disponibles en la literatura junto con una nueva basada en un formulaciones extendidas tipo red. Se investiga su desempeño sobre un conjunto de instancias reales desde el punto de vista de su eficiencia computacional al ser resueltas con softwares comerciales. Lo anterior incluye análisis de tiempos de solución, nodos utilizados e iteraciones de simplex realizadas para distintas tolerancias requeridas. Los experimentos muestran la calidad de la aproximación propuesta, siendo esta completamente competitiva respecto a las ya documentadas. Este resultado era esperable, dada la estructura totalmente unimodular de gran parte de la formulación propuesta, pero para nada justificable debido al tamaño de la misma. Lo anterior muestra que el efecto del preproceso de los softwares comerciales puede ser fundamental en algunas formulaciones. Por otro lado, respecto a la función objetivo del problema de predespacho de unidades, que por lo general se representa como una función cuadrática de la generación, se presenta una nueva manera de linealizar su comportamiento de modo que su inclusión en una formulación entera-mixta lineal tradicional sea eficiente. Esto último debe entenderse a partir de la necesidad que el tamaño de la aproximación no crezca de manera desmedida si el error requerido para la misma decrece. Si bien ya existía la posibilidad de hacer esto mediante la aplicación de la aproximación desarrollada por Ben-Tal y Nemirovsky para conos de segundo orden [2], acá se presenta un método alternativo, con mejores propiedades numéricas, un orden de magnitud mejor en calidad de aproximación, y cuya aplicación a problemas reales de predespacho de unidades genera mejores resultados respecto de las aproximaciones tradicionales. Por último, con el fin de mejorar el desempeño de la formulación entera-mixta presentada, se realiza el análisis poliedral de una de sus subestructuras esperando identificar desigualdades válidas que permitan mejorar su cota dual. Esta subestructura corresponde al knapsack semicontinuo con restricciones adicionales del tipo generalized upper bound. Se demuestra que bajo supuestos simples es posible identificar facetas tipo generalized flow cover en espacios restringidos de dimensión inferior. Luego se llevan estas desigualdades al espacio original utilizando procedimientos de lifting multidimensional independiente de la secuencia [38, 27, 16, 17] y se iii

prueba que con supuestos adicionales también son facetas allí. Experimentos computacionales en instancias derivadas de problemas de UC muestran su eficiencia, donde más de un 50 % del gap integral del nodo raíz se reduce aplicando en promedio solo tres de estos cortes. Además, en este contexto, también se ha implementado un solver ad-hoc para la solución eficiente de las relajaciones lineales de la formulación tipo red, con un speed-up del orden de 4x a 8x respecto a CPLEX barrier optimizer, pero que aún no está documentado.

iv

A mis padres, Yolanda y Hugo

v

Tabla de contenido Tabla de contenido

vi

1. Introducción a los artículos 1.1. Motivación . . . . . . . . . . . . . . . . . . . . 1.1.1. Objetivos . . . . . . . . . . . . . . . . 1.1.2. Estructura de la tesis . . . . . . . . . . 1.2. Introducción a sistemas eléctricos de potencia 1.2.1. Sistema de transmisión . . . . . . . . . 1.2.2. Sistema de distribución . . . . . . . . . 1.2.3. Productores . . . . . . . . . . . . . . . 1.2.4. Consumidores . . . . . . . . . . . . . . 1.2.5. Mercados . . . . . . . . . . . . . . . . 1.3. Planificación estratégica de sistemas eléctricos 1.3.1. Planificación de mediano y largo plazo 1.3.2. Planificación de la operación . . . . . . 1.4. Predespacho de unidades de generación . . . . 1.5. Discusión . . . . . . . . . . . . . . . . . . . .

1 1 2 2 4 4 5 5 5 6 7 7 8 9 12

Bibliografía

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

15

2. Primer artículo: A network flow–based MILP formulation for the thermal unit commitment problem. 18 3. Segundo artículo: A polyhedral–based approach applied to quadratic cost curves in the unit commitment problem. 21 4. Tercer artículo: Sequence independent, simultaneous and multidimensional lifting of generalized flow covers for the semi–continuous knapsack problem with generalized upper bounds constraints. 31 5. Cuarto artículo: Sequence independent lifting for mixed knapsack problems with GUB constraints. 44

vi

Capítulo 1 Introducción a los artículos

1.1.

Motivación

La dinámica que actualmente tienen los sistemas eléctricos de potencia, el impacto económico de su adecuada planificación y control, y la cantidad de actores que se ven involucrados en su funcionamiento, hacen que el estudio desde los puntos de vista técnico, económico y social de su comportamiento sean de vital importancia en las políticas de un país. Dentro de este contexto se encuentra un conjunto de problemas que tienen relación con la adecuada explotación del sistema y que se denominan problemas de planificación estratégica. Como en otros casos, a estos se les asocian distintas escalas de tiempo y según este criterio podrán clasificarse en problemas de largo plazo, que tienen relación con inversiones en infraestructura y uso adecuado del agua; corto plazo, que deciden la operación técnica y económicamente óptima del sistema; y mediano plazo, que realiza la conexión consistente entre la planificación de largo plazo y la operación de corto plazo. Dentro de los problemas de corto plazo, uno de los más importantes es el predespacho de unidades (Unit Commitment, UC) y tiene relación con selección de los generadores que serán utilizados para suplir la demanda y la reserva de un sistema en periodos de tiempo que van desde un día hasta un par de semanas [39]. El impacto que la adecuada solución de este problema tiene sobre los costos del sistema es bien conocido y ha sido documentado en varios reportes [28, 29]. Es posible estimar que el impacto de un punto en la mejora de su solución, es equivalente a un ahorro del orden de 50MUS$ anuales, solo en costos operacionales en sistema medianos como chileno [39]. Para resolver el problema de UC en sistemas reales, durante los último 40 años, por lo general se han utilizado aproximaciones del tipo relajación Lagrangeana [26, 37, 23, 3, 33, 18]. Esto debido principalmente a la alta eficiencia computacional de este tipo de algoritmos. En los último 10 años, debido al incremento en la capacidad de procesamiento, al aumento en la memoria disponible y al desarrollo constante de los algoritmos del tipo Branch&Bound&Cut en software comerciales tipo black-box (CPLEX, Gurobi, XPress), el uso de formulaciones lineales entera-mixtas (Mixed Integer Linear Programming, MILP) y su solución a través 1

de este tipo de herramientas resulta cada vez más atractivo [28, 35]. Lo anterior tiene su origen en que a diferencia de la relajación Lagrangeana, esta última aproximación tiene mayor flexibilidad (permite la incorporación simple de nuevas características en el modelo) y permite la solución a optimalidad del problema si se dispone de tiempo suficiente (mejora continua de cotas primales y duales). Bajo este nuevo paradigma, donde se requieren buenas soluciones del problema de UC en corto tiempo y donde se utilizan formulaciones MILP para su modelamiento (que se resuelven con softwares comerciales), resulta conveniente investigar los distintos aspectos que tradicionalmente permiten mejorar la eficiencia de este tipo de aproximaciones. En este contexto, la evaluación del efecto que tienen las distintas formulaciones (existentes y nuevas) sobre los tiempos de solución del problema resulta fundamental. Además, definida la “mejor formulación”, se debe investigar su estructura de modo de suministrar información adicional al solver comercial que facilite la mejora de las cotas en tiempos menores. Todo esto abre un conjunto de posibilidades de investigación, que son el origen del trabajo aquí presentado.

1.1.1.

Objetivos

Los objetivos definidos para esta tesis tienen relación con la investigación del problema de UC y su solución eficiente en un entorno del tipo Branch&Bound&Cut. Su especificación se muestra a continuación: 1. Formular un modelo eficiente de programación entera-mixta lineal para el problema de UC, que permita la inclusión de nuevas características de las unidades de generación (ciclos combinados por ejemplo), de nuevas restricciones operacionales sistémicas (distintos tipos de reserva) y cuya solución en solvers comerciales sea competitiva respecto a modelos clásicos. 2. Dado que los costos de las unidades se modelan como funciones cuadráticas de la potencia generada, en el contexto de una formulación MILP estos deben ser linealizados. La idea es que la calidad requerida para esta aproximación no incida en gran medida en el tamaño del problema. En este contexto, la investigación de formulaciones extendidas para su representación es de vital importancia. 3. Investigar la solución eficiente de la mejor formulación identificada en un entorno Branch&Bound&Cut como CPLEX. Esto incluye la generación de cortes para mejorar la cota dual y el desarrollo de algoritmos ad-hoc para la solución eficiente de los nodos.

1.1.2.

Estructura de la tesis

Esta tesis se divide en cinco capítulos; en el primero se desarrolla de manera breve la problemática asociada a la planificación estratégica en sistemas eléctricos de potencia, mostrando sus características y particularidades; mientras que en los siguientes cuatro capítulos, debido a que la presente tesis fue desarrollada como un compendio de publicaciones, se incluyen los artículos en inglés que se generaron durante el desarrollo de este trabajo. 2

En el capítulo 2 se presenta la publicación titulada “A network flow–based MILP formulation for the thermal unit commitment problem”, donde se investiga el uso de una formulación extendida basada en redes para modelar el problema de UC y se contrasta con la mejor formulación MILP disponible en la bibliografía. Los resultados muestran que la nueva formulación es competitiva desde el punto de vista de los tiempos computacionales, pero más atractiva desde el punto de vista de sus posibilidades de mejora. Esto último dado que su particular estructura posibilita la incorporación de nuevos cortes y facilita el desarrollo de nuevos algoritmos ad-hoc para la solución de sus relajaciones. Este trabajo se encuentra en carácter de enviado a IEEE Transaction on Power Systems (9-Marzo-2015). Luego, en el capítulo 3 se presenta la publicación titulada “A polyhedral–based approach applied to quadratic cost curves in the unit commitment problem”, donde los costos cuadráticos asociados a la generación de las unidades son convenientemente linealizados con una aproximación inferior eficiente. Este trabajo tiene su origen en la pérdida de precisión de la aproximación para conos de segundo orden de Ben-Tal y Nemirovsky cuando se usa en funciones cuadráticas unidimensionales [2, 15]. Se presenta una nueva formulación extendida que resuelve el problema anterior, que crece solo logarítmicamente respecto de la precisión requerida – ε = (τ −1)12 4ν+1 , donde ν y τ son parámetros que definen el tamaño de la aproximación en forma lineal – y que además es representable mediante números racionales. Experimentos computacionales muestran que en general la nueva aproximación es mejor que la solución directa del problema mediante programación entera-mixta cónica de segundo orden (MixedInteger Second Order Cone Programming, MISCOP) o que aquella obtenida mediante la aplicación de la aproximación de Ben-Tal y Nemirovsky. Este trabajo también se encuentra en carácter de enviado a IEEE Transaction on Power Systems (9-Marzo-2015). Finalmente, en los capítulos 4 y 5 se presenta las publicaciones tituladas “Sequence independent, simultaneous and multidimensional lifting of generalized flow covers for the semi– continuous knapsack problem with generalized upper bounds constraints” y “Sequence independent lifting for mixed knapsack problems with GUB constraints”, que corresponden a la versión original presentada en IPCO-2014 y a su versión extendida publicada en Mathematical Programming Series B. En ellas se analiza una subestructura de la formulación de red para el problema de UC, que corresponde a un “knapsack ” semicontinuo con restricciones adicionales del tipo “generalized upper bound ”. La idea es determinar cortes que eliminen eventuales punto fraccionarios obtenidos de la solución del problema de UC mediante un algoritmo del tipo Branch&Bound&Cut. Para ello, primero se realiza un análisis poliedral del conjunto y se demuestra que la desigualdades del tipo “flow cover ” generalizado, bajo supuestos mínimos, son facetas del poliedro proyectado en un espacio de dimensión menor. Luego, dada la desigualdad anterior, se realiza el lifting simultaneo e independiente de la secuencia de dos variables a la vez y se lleva el resultado a la dimensión del problema original, demostrándose que bajo supuestos adicionales, ellas también son facetas allí. A continuación se presentan procedimientos eficientes para la generación de la desigualdad semilla a partir del punto fraccionario, para el cálculo de las funciones super-aditivas de lifting y para la selección de la(s) desigualdad(es) válida(s) a aplicar; esto último debido a que el procedimiento de lifting desarrollado es multidimensional y podría haber un número exponencial de ellas. Por último, mediante experimentos computacionales aplicados a instancias obtenidas de problema de UC reales, se muestra la calidad de las desigualdades generadas, las cuales cierran un 57,7 % de gap integral del nodo raíz y agregan menos de 3 cortes en promedio.

3

Ambas publicaciones fueron aceptadas, la primera en Integer Programming and Combinatorial Optimization (18-Enero-2014) y la segunda en Mathematical Programming, Series B (20-Febrero-2015).

1.2. Introducción a sistemas eléctricos de potencia Residential electricity costs ||| Fact sheet 1 Con el fin de presentar un lenguaje básico que permita entender cómo funcionan los sisteThe supply chain mas eléctricos de potencia, a continuación se muestran los distintos actores que forman parte de su estructura básica. A diferencia de otros, el mercado eléctrico se desarrolla en tiempo real, debiendo cumplir siempre un conjunto de leyes físicas que acoplan su funcionamiento exact amount charged by retailers to supply electricity to a household varies depending on the retailer tantoThe dinámico como estacionario. Detalles de cómo esto se modela matemáticamente pueden and type of consumer account, but electricity bills generally reflect the costs incurred by participants encontrarse [12], donde autor entrega un breve resumen de and las isleyes constitutivas del across theen electricity industry.elThis fact sheet describes the industry structure the first in a series outlining the costs1 are made and how comparede to la prices paid by consumers such as mundo de lahow electricidad que sonupútiles en they el mundo gestión deother operaciones. Además, en business customers. la Fig. 1.1 se muestra la estructura típica asociada a esta cadena se suministro donde cada uno de sus actores es identificado. The electricity supply chain in New Zealand local gener- ators

direct consumers

Retail

Wholesale

Grid exit points

Grid injection points

Consumers

Distribution

Transmission

Generation

230 volts for residential and business

11, 33 or 66kV

typically 110 or 220kV

11–16kV

Source: 2011 Electricity Authority Electricity in New Zealand.

Figura 1.1: Típica cadena de suministro del sector electricidad (Fuente: Electricity Authority Electricity in Most electricity is produced by generators located well away from where it is eventually used. This is often because of the geographical location of energy sources, for example, rivers used for hydro-generation, geothermal fields, or the location of good wind generation sites. There are also other factors to consider when placing large generation close to consumers and communities.

New Zealand)

Around 90 percent of the total electricity New Zealanders use passes through the high

1.2.1. de that transmisión voltage Sistema transmission system spans the country and is known as the national grid.

90 percent of the total electricity New Zealanders use passes through the high voltage transmission system

The transmission system delivers electricity at high voltage to substations in each area. Large cities may have several substations serving them and in a few cases electricity is Elsupplied sistema de transmisión, también conocido como red, es el elemento de un sistema de directly to large industrial consumers such as the Tiwai Point aluminium smelter. potencia encargado de trasferir energía eléctrica de un lugar a otro. Por lo general, esto se Local distribution systems take the power delivered to each substation and deliver it at lower homes and businesses. The remaining of electricity that pass through the donde se hace voltages desde tolas centrales de generación a 10 laspercent subestaciones de doesn't distribución, desde transmission system is generated by plant that is directly connected to the local distribution system or, in the case suministra el recurso a los consumidores finales, a través de líneas aéreas o cables subterráneos. of some large industrial consumers, by their own on-site generation.

Dado que estas líneas tienen restricciones de capacidad debido a su límite térmico, entonces 1. The term ‘cost’ in this paper includes the costs associated with producing electricity, including any profits. es posible encontrar congestiones en su operación. Hoy en día, en gran parte de los sistemas eléctricos, los negocios de transmisión y generación están separados, siendo generalmente el primero controlado por un operador independiente de sistema y regulado por el estado (debido a su natural carácter monopólico). Este operador mantiene la red y provee de acceso a la energía eléctrica a los distintos actores: generadores, distribuidores y consumidores finales. 4

Además, es este operador el encargado de realizar el balance instantáneo entre demanda, generación y pérdidas del sistema, de modo que este funcione siempre de manera segura.

1.2.2.

Sistema de distribución

Este sector se caracteriza por la existencia de áreas exclusivas de prestación del servicio dentro de un territorio. Debido a la existencia de monopolios geográficos, se hace necesaria la existencia de regulación y mecanismos orientados a incentivar a que las empresas se desarrollen en forma competitiva. Dos estrategias orientadas al logro de este objetivos son la empresa modelo (caso chileno) y la evaluación relativa del desempeño (caso inglés). Los precios de distribución, por tratarse de una actividad regulada, deben permitir cubrir los costos totales de la actividad, que básicamente son de inversión, operación y mantenimiento [5].

1.2.3.

Productores

Como productores de un sistema eléctrico se entienden a todos los entes que inyectan energía en la red. Esto incluye a las unidades de generación térmica e hidráulica, junto con las plantas que generan a partir de energías renovables como las eólicas y las solares. En las plantas térmicas, a menudo se quema algún tipo de combustible (carbón, gas, diesel, nuclear, etc.) para transformarlo en energía eléctrica mediante el uso de turbogeneradores. Las plantas hidráulicas consisten en un conjunto de reservorios, posiblemente conectados en cascada, los cuales suministran un caudal de agua que se utiliza para mover turbinas conectadas a sendos generadores. De este modo la energía potencial del agua se transforma en energía eléctrica que es inyectada a la red. Se debe considerar que dada la capacidad que estos reservorios tienen para almacenar energía, su correcta operación en conjunto con las plantas térmicas requiere de una adecuada coordinación. Esto último incluye un tipo de plantas híbridas donde un reservorio se utiliza en conjunto con una unidad de turbina-generador que también puede operar en modo bomba-motor. De este modo, se impulsa agua al reservorio en horas cuando el costo de la energía es bajo y se libera a la unidad de generación cuando los precios son altos. En general, los productores se presentan en tamaños no uniformes, desde pequeñas plantas de generación diesel o solar a grandes conglomerados que controlan plantas con distintas tecnologías y distribuidas espacialmente sobre la red. Respecto a cómo venden su energía, las formas son variadas también, e incluyen contratos bilaterales, ventas a precio spot o participación en bolsas de energía; dependiendo del procedimiento que utilicen para participar en el mercado, actuarán como tomadores de precio o serán capaces de influir sobre el mismo.

1.2.4.

Consumidores

Los consumidores son los usuarios finales de la energía eléctrica y corresponderán a plantas industriales o a consumos residenciales. La adquisición del bien la hacen mediante contratos bilaterales con los productores, directamente en los mercados de energía o comprando a 5

empresas de distribución o a comercializadores. Por el tamaño de sus consumos, pueden clasificarse como clientes regulados o no regulados, donde los primeros son tomadores de precio y los segundos tienen la posibilidad de negociar este valor con los productores [22].

1.2.5.

Mercados

Como resultado de la interacción de los agentes en el mercado eléctrico, existen varias maneras en que este se organiza. Desde el punto de vista de su operación se distinguen cuatro formas básicas para comprar y vender energía: pool o mancomunidad, bolsa de energía, contratos bilaterales físicos y contratos bilaterales financieros. En el modelo clásico de pool, productores y consumidores renuncian a establecer relaciones comerciales directas entre ellos. Las compras y ventas de energía son determinadas y valorizadas por el operador de mercado en base a una optimización de los costos totales del sistema. El plan de operación resultante es transferido al operador de sistema, quien verifica la factibilidad técnica del mismo y determina los servicios auxiliares requeridos para su ejecución. Otro modelo es una bolsa de energía, que es una entidad que recibe ofertas por la compra y venta de energía y establece la casación entre ellas. Generalmente, es un lugar virtual donde los agentes del mercado se reúnen y puede considerarse como un pool con ciertas características especiales como: estandarización de productos, no hay incidencia en el despacho final real (se requiere de un operador de red), la participación no es obligatoria y la información de los Participants and money flows in the electricity industry agentes es mucho más reservada.

Distribution

Transmission

Pay for spot sales

Generators

Electricity flows

Pay transmission and distribution charges

Spot market Pay retail bill Pay for spot purchases

Money flows

Consumers

Retailers

Payments for contracts and bills

Source: Electricity Authority2

Figura 1.2: Flujos monetarios entre los distintos agentes en un mercado eléctrico (Fuente: Electricity Authority Electricity in New Zealand)

Retailers purchase electricity from generators and charge consumers for its use. The sale and purchase of electricity between generators and retailers is carried out 6 in the wholesale electricity spot market. Many retailers operate their own generation, but because of the way that both consumption and generation change over time, all retailers need to buy electricity through the wholesale market at times. In most areas, retailers also manage

Los contratos bilaterales, ya sean físicos o financieros, son producto de un libre intercambio comercial entre productores y consumidores, ya sea en forma directa o a través de un comercializador. Lo que caracteriza a un contrato bilateral físico es su relación directa con el despacho de la operación resultante. Sin embargo, los contratos bilaterales financieros no afectan al despacho de la operación, ya que ellos tienen por objeto manejar, acorde a una estrategia de mercado, el riesgo del precio futuro de la energía eléctrica. Los mercados reales se forman como combinaciones de alguna de estas modalidades, pudiendo corresponder a una de ellas o a un híbrido que contenga a varias simultáneamente. El esquema mostrado en la Fig. 1.2, que corresponde al caso neozelandés y cuya estructura es muy similar al caso chileno, permite entender el funcionamiento de estos mercados, donde los flujos monetarios identificados tienen directa relación con las distintas formas de tranzar que tienen los agentes [5].

1.3.

Planificación estratégica de sistemas eléctricos

La planificación estratégica de un sistema eléctrico comprende todas las decisiones que deben tomarse relativas a la operación de este, de modo que su correcto funcionamiento esté técnicamente asegurado tanto para el corto como para el largo plazo. Esto debe realizarse considerando que la suma de los costos involucrados, presentes y futuros, debe ser minimizada. Otros criterios para la toma de decisiones también son válidos: mejorar niveles de seguridad, minimizar la ocurrencia de fallas de abastecimiento, minimizar pérdidas, etc. La Fig. 1.3 muestra cómo los sistemas eléctricos son optimizados en el tiempo; en el corto plazo el detalle técnico es relevante, dada la mínima incertidumbre, mas en el mediano-largo plazo la incertidumbre crece y el modelamiento fino de los aspectos técnicos ya no resulta crítico. En el primer caso las decisiones son por lo general operacionales, mientras que en el segundo son de inversión o de valuación de algunos bienes (agua, por ejemplo). La conexión consistente de estos problemas es todavía un problema abierto, habiendo aproximaciones iniciales que, haciendo uso de las nuevas tecnologías, intentan realizar esto coordinando el control en tiempo real con la planificación de más corto plazo (optimal power flow, economic dispatch y unit commitment) [21, 7].

1.3.1.

Planificación de mediano y largo plazo

La planificación de más largo plazo (5-20 años) esta enfocada a sistemas hidrotérmicos donde se determinan las políticas de uso de las reservas de regulación inter-anual de manera tal que se minimice el costo futuro “esperado” de la operación del sistema eléctrico. Por tratarse de un horizonte de tiempo largo, los modelos de optimización usados son de naturaleza estocástica, donde se brinda especial importancia al modelamiento de la incertidumbre de los recursos hídricos. Para sistemas eléctricos térmicos, los horizontes de planificación de largo plazo son más cortos (3-5 años), y allí se definen los grandes contratos de combustibles, las nuevas unidades 7

3

A PRIMER ON OPTIMAL POWER FLOW

Planning Horizon Years

Months

Weeks

Days

Hours

Minutes

asin Incre Decr

ion solut el Re d o gM

easin g

Long-term Generation Scheduling

Unit Commitment

Economic Dispatch

Capacity Expansion Planning

SecurityConstrained Unit Commitment

Classic Optimal Power Flow

SecurityConstrained Unit Commitment

Optimal Reactive Power Flow

Real-Time

Unce rtain

ty

Automatic Generation Control

Real Power Dispatch Reactive Power Dispatch Reactive Power Planning

Voltage Control

Fig. 1.1. Optimization in power system operation via incremental planning. Long-term plan-

Figura 1.3: Optimización de decisions los sistemas eléctricos víamodels. planificación incremental ning procedures make high level based on coarse system Short-term procedures fine [12] tune earlier decisions, using detailed models but a more limited decision space. Bold text indicates planning procedures which incorporate variants of optimal power flow.

generadoras, etc. • linear algebra [16], En 50el mediano plazo la información con que se cuenta es más detallada que para el caso 51 • complex number theory [13, 16], del largo plazo• yanalysis está compuesta por: análisis de la predicción de la demanda de energía 52 of differential equations in the frequency domain [16], and y demanda máxima del sistema, disponibilidad de and las unidades contratos de 53 • linear and nonlinear optimization theory application generadoras, [20, 24]. intercambio de potencia energía entre empresas precio[21] y consumo esperado de 54 Readers who alsoyhave a working knowledge generadoras, of electrical circuits and electric combustibles las centrales coordinación de losofprogramas de mantenimiento de 55 powerdesystems analysis térmicas, [14] will find the development the power flow equations 56 familiar. readers may to expand their understanding consulting a good en pasos las empresas, etc.Other El horizonte dewish estudio es generalmente de un by año, discretizado 57 power systems text as [14] principal or [32]. However, priorla familiarity with power flow iso semanal mensuales o semanales, y elsuch objetivo es realizar programación mensual 58 not strictly required in order to follow the development presented in this primer. de la generación. La coordinación óptima de los programas de mantenimiento también puede 59 The primer begins with a guide to OPF notation in §2, including notational differser incluida aquí. 60 ences between electric power systems engineers and Operations Researchers. Sections 3 and 4 introduce the modeling of electric power systems and the power flow equations; these fundamental topics are omitted in most other introductory OPF texts. 63 upon the previous §5 discusses OPF formulations; this section in1.3.2. Building Planificación de lasections, operación 64 cludes full formulations for several of the decision processes shown in Figure 1.1. §6 65 provides a descriptive guide to two common file formats for exchanging PF and OPF La 66planificación plazo data. Finally,de§7 corto concludes the suele primer.dividirse en planificación semanal, programación 61 62

diaria y programación horaria. En la planificación semanal, se determina la operación más conveniente de los embalses con capacidad de regulación semanal y se define que unidades térmicas participarán del despacho. La modelación del parque generador hidráulico y térmico así como la red de transmisión poseen mayor nivel de detalle que en la programación de largo/mediano plazo. Por lo general, un modelo determinístico se considera suficiente, mas 8

con la penetración de las energías renovables esto es cada vez más discutible. La programación diaria es un reajuste de la programación semanal. En este período el pronóstico de la demanda y de los caudales tiene menor incertidumbre y el sistema - generación hidráulica, térmica y red de transmisión - es modelado con mayor detalle. En el caso de la programación horaria, entendido como un flujo de potencia óptimo, el modelo de los elementos constitutivos del sistema es todavía más fino y como resultado de su desarrollo se determinan las consignas de tensión, potencia activa y reactiva de todas las unidades [5].

1.4.

Predespacho de unidades de generación

El problema de predespacho de unidades de generación consiste en definir el conjunto de centrales que serán utilizadas para suplir la demanda y la reserva de un sistema eléctrico, respetando las restricciones técnicas y minimizando los costos operacionales, en horizontes de planificación que van desde un día hasta un par de semanas [39]. Puede incluir o no la red, y considerar a todas las unidades o solo a parte de ellas. En el caso en que la red no se considera y solo se incluyen las unidades térmicas, el problema se denomina predespacho de unidades térmicas uninodal. Desde el punto de vista de modelación matemática, este es un problema no lineal entero mixto de gran escala (Mixed-Integer Nonlinear Programming, MINLP), cuya complejidad es bien conocida [36] y cuya solución con state of the art softwares, todavía no es posible cuando el número de unidades de generación o el número de períodos considerados crece en demasía [10]. Dado que generalmente los costos de operación de las unidades se representan como funciones cuadráticas dependientes del nivel de generación [6], es posible resolver el problema de UC mediante técnicas más específicas de modelación como la programación cuadrática entera mixta (Mixed-Integer Quadratic Programming, MIQP) [9, 20, 40, 30]. Si bien la solución de este tipo de problemas con softwares comerciales como CPLEX o Gurobi ha mejorado sustancialmente en los últimos años, aún no posibilita la solución de instancias reales. Un camino alternativo al anterior, y que tiene relación con que en la práctica muchas veces las soluciones óptimas no son necesarias, corresponde a la aproximación del problema de UC mediante algún tipo de linealización conveniente. Con esto se posibilita el uso de la programación lineal entera mixta en la solución del problema de UC, incorporando en este proceso todas las ventajas que una tecnología madura tiene [4]. De hecho, como ya se ha comentado, el uso de este tipo de modelos en la operación de sistemas eléctricos reales ya ha generado un impacto económico positivo en ellos [28]. Bajo este contexto de modelación, una formulación típica para el problema de predespacho 9

de unidades térmicas uninodal es la que se muestra a continuación: P P m´ın cg (ptg , utg )

(1.1)

t∈T g∈G

s.t.

P

g∈G

P

ptg = Dt ,

g∈G

∀t ∈ T

pg utg ≥ Dt + Rt ,

ptg , utg ∈ Πg ,

(1.2)

∀t ∈ T

∀t ∈ T , ∀g ∈ G

(1.3) (1.4)

en donde se busca minimizar los costos totales de operación (1.1), sujeto a restricciones de balance entre generación y demanda (1.2), de satisfacción de requerimientos mínimos de reserva (1.3) y de operación de las unidades dentro de sus límites técnicos (1.4). En esta formulación, T son los periodos en estudio, G es el conjunto de todas las unidades de generación, ptg es la potencia inyectada por la unidad g en el periodo t, ugt es una variable binaria que modela el estado de la unidad g en el periodo t, cg (·) es la función de costos de la unidad g convenientemente linealizada, Dt es la demanda del sistema en el periodo t, Rt la reserva mínima requerida por el sistema en el periodo t, pg es la máxima potencia generable por la unidad g y Πg es la región de operación factible de la unidad g que se aproxima linealmente a través de un poliedro. En la Fig. 1.4 se muestra el comportamiento típico de la operación de una central térmica en el tiempo. Allí se observan los procesos de sincronización, arranque, operación, parada y detención de las unidades, los cuales dificultan 1968 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 25, NO. 4, NOVEMBER 2010 la modelación del conjunto Πg , pues inducen relaciones inter-temporales entre las variables.

Fig. 1. Operating states of a thermal unit.

Figura 1.4: Estados operacionales de una unidad térmica [32] straints on the continuous and the binary decision variables associated with the operation of the unit, denoted by the vector Dado variable que en ,elfurther problema detheUC térmico described in following. decision

TABLE I UNIT START-UP MODELING

uninodal se intenta representar de la manera más fiel posible el comportamiento de las unidades (como el mostrado en la Fig. 1.4), por lo A. Unit Operating State Modeling general se asume que su modelación incluye al menos las siguientes características: Fig. 1 shows the different operating states of a thermal unit for hours [31]. After being reserved 1. Condiciones iniciales: Número de horas que la unidad ha estado detenida u operativa. the unit starts-up at hour and remains comuntil it is shut-down at hour mitted 2. Límites de generación: Corresponde. Once a los valores mínimo y máximo entre los cuales la committed, the unit follows four consecutive operation phases, unidad puede ser despachada. 1) synchronization, 2) soak, 3) dispatchable, and 4) desynchro, , , and nization, denoted by binary variables maximum reservation time to cold start of all units) and is fur, respectively (Fig. 1). The first two phases comprise the 10 ther explained in the Appendix. unit start-up phase. During the dispatchable phase, the unit can Once a thermal unit enters a hot, warm, or cold start-up phase, receive dispatch instructions to vary its power output between it should complete the start-up sequence, and enter the dispatchits technical minimum and its nominal power output according

3. Generación mínima de partida/parada: Permite considerar las rampas de subida y bajada cuando la unidad cambia de estado operacional. 4. Rampas de subida y bajada: Cuando una unidad opera por sobre su mínimo técnico, el aumento y la disminución de la generación entre dos periodos consecutivos está limitado. 5. Tiempos mínimos en operación y detenido: Debido a las características termo-mecánicas de las unidades, cuando ellas parten deben estar un número mínimo de horas operativas; por otro lado, cuando se detienen se debe disponer de un tiempo mínimo para realizar su mantenimiento. 6. Unidades siempre operativas: Por su alta eficiencia algunas unidades no participan del despacho, pero su inyección de potencia debe ser considerada. 7. Unidades no operativas: Algunas unidades deberán estar fuera de servicio por algunos periodos debido a un mantenimiento programado a una salida forzada. 8. Unidades de partida rápida: Son unidades que pueden entrar en funcionamiento o ser detenidas en poco tiempo. Si se utilizan, generalmente es en las horas punta del sistema. 9. Restricciones conjuntas sobre unidades: Permite modelar eventuales acoplamientos entre las características operativas de distintas unidades (ciclos combinados o límites en la red, por ejemplo). 10. Costos de partida y parada: Se incluyen los costos asociados a poner en servicio la unidad como los incurridos por el mantenimiento realizado cuando está detenida. 11. Rampas de partida y parada: Incluyen la energía inyectada por las unidades cuando ellas están por debajo de sus mínimos técnicos. Detalles específicos de cómo cada una de estas características se modela pueden encontrarse en [39, 31, 34]. Allí es posible observar un fuerte comportamiento no lineal en algunas de estas características, donde, a modo de ejemplo, en la Fig. 1.5 se muestran las típicas funciones de costo. De la figura de la izquierda se deduce la gran discrepancia que existiría entre la data empírica y una eventual aproximación polinómica para el net heat rate (razón entre los costos de generación y la potencia, donde los primeros son generalmente representados por polinomios cuadráticos). 13

1.6 CHARACTERISTICS OF STEAM UNITS

11

Start-up cost

1.6 CHARACTERISTICS OF STEAM UNITS

1

2

3

Off-line time

Figure 1.8  Approximate net heat rates for a range of simple cycle gas turbine units. Units are fired by natural gas and represent performance at standard conditions of an ambient temperature of 15°C at sea level. (Heat rate data from reference 1 were adjusted by 13% to represent HHVs and auxiliary power needs).

Figura 1.5: Costos de unidades térmicas. De izquierda a derecha se observa: aproximación de la curva net heat rate a partir de data empírica, característica entrada-salida de un conjunto turbina-generador con cuatro válvulas y aproximación de los costos de partida mediante gears. In larger units the generators are driven directly, without any gears. Exhaust función [39]cycle units. In combined cycles, gases are discharged to theexponencial atmosphere in the simple the exhaust gases are used to make steam in a heat-recovery steam generator (HRSG) before being discharged. The early utility applications of simple cycle gas turbines for power generation after World War II through about the 1970s were generally to supply power for peak 11 load periods. They were fairly low-efficiency units that were intended to be available for emergency needs and to insure adequate generation reserves in case of unex1.6  Characteristics of a steam turbine generator with four steam admission valves. pected load peaks or generation outages. Full-load netFigure heat rates were typically

1.5.

Discusión

En las secciones anteriores se ha introducido de manera progresiva el problema de UC térmico, el cual es el centro de estudio del presente trabajo. Se presentaron primero los actores que forman parte de la cadena de suministro del sector electricidad y sus interacciones a nivel de mercado. Luego, de manera breve, se exponen los problemas que comúnmente atañen a la planificación estratégica del sector, tanto en el corto como en el largo plazo, indicando las particularidades que comúnmente presentan respecto a la incertidumbre y la precisión técnica de los modelos. Finalmente, el problema de UC es presentado a través de un modelo entero-mixto que captura gran parte de sus características relevantes. No se considera una formulación más fina, pues el objeto de esta introducción es mostrar la problemática asociada a su modelamiento y no la solución particular que varios autores han dado a ella [14, 24, 26, 37, 8, 23, 3, 18, 1, 11, 6, 13, 10, 20, 25, 30]. De la revisión de estos trabajos se deduce que existe una gran cantidad de desafíos pendientes relativos al a la solución del problema de UC térmico. Esto incluye aspectos tanto de modelamiento como de algoritmos e implementación computacional. En la referencia [19], es posible encontrar un buen resumen de estos desafíos y posibles rutas para poder enfrentarlos. La tabla mostrada en las Fig. 1.6 y 1.7, obtenida de esta fuente, muestra estos desafíos agrupados según su estado de desarrollo y ordenados según su relevancia. Se observa que muchos de ellos tienen que ver con un mejor tratamiento de la incertidumbre y muchos otros con un modelamiento de mercado donde los agentes tienen mayor flexibilidad. Si bien estos desafíos se levantaron hace ya 15 años en un workshop desarrollado para tal efecto, casi todos ellos siguen siendo todavía fuente de fuerte investigación. Particular interés se tiene hoy por los puntos dynamic feasibility, other commodities y uncertainty, variability debido a la fuerte penetración de las energías renovables y el efecto que ellas tienen en la operación segura de los sistemas eléctricos frente a eventuales contingencias (reservas, rampas). Además, el punto time scale – donde se intenta conectar las decisiones secuenciales que se toman en corto plazo (operación) con las de tiempo real (control) – resulta ser cada vez más interesante para los investigadores debido principalmente a que en los último años se ha identificado una modelación-solución eficiente del problema de flujo óptimo a nivel de distribución mediante ADMM (Alternating Direction Method of Multipliers); y a la necesidad de que los sistemas eléctricos a nivel de distribución operen de manera más autónoma e inteligente debido a su evolución natural, donde ahora coexiste la demanda tradicional con una gran cantidad de pequeños medios de generación (smartgrids). En fin, todavía hay mucho por hacer en esta área y la dinámica que los mercados eléctricos tienen ahora y presentarán en el futuro no hacen otra cosa que incentivar un estudio más profundo de cada uno de los aspectos técnicos y económicos que movilizan al sector. Por lo tanto, dadas las especiales características de los sistemas eléctricos de potencia, se abren infinitas puertas para desarrollar investigación relativa al uso de la gestión de operaciones en su explotación, siendo este trabajo solo el inicio de ello y que en el futuro pretende cubrir también la problemática de mediano y largo plazo.

12

Figura 1.6: Desafíos para el problema de UC categorizados y ordenados por importancia (III=requiere investigación fundamental; II=Necesidad de desarrollo e implementación; I=Ya existen softwares para su modelamiento) [19]

13

10

The Next Generation of Unit Commitment Models

Figura 1.7: Continuación tabla

14

Bibliografía [1] J. M. Arroyo and A. J. Conejo. Modeling of start-up and shut-down power trajectories of thermal units. IEEE Trans. on Power Systems, 19(3):1562 – 1568, 2004. [2] A. Ben-Tal and A. Nemirovski. On polyhedral approximations of the second-order cone. Math. Oper. Res., 26(2):193 –Ű 205, 2001. [3] D. P. Bertsekas, G. S. Lauer, N. R. Sandell Jr., and T. A. Posbergh. Optimal short-term scheduling of large-scale power systems. IEEE Trans. on Automatic Control, 28(1):1 – 11, 1983. [4] R. Bixby, Rothberg, and Z. Gu. The latest advances in mixedĄinteger programming solvers, presentation notes from the cirrelt spring school on combinatorial optimization in logistics, montreal, 2010. [5] Carlos Cabrera. Aplicación de algoritmos genéticos al predespacho de unidades térmicas usando flujo óptimo de potencia. PhD thesis, 2008. [6] M. Carrion and J. M. Arroyo. A computationally efficient mixed-integer linear formulation for the thermal unit commitment problem. IEEE Trans. on Power Systems, 21(3):1371 – 1378, 2006. [7] Emiliano Dall’Anese, Sairaj V. Dhople, and Georgios B. Giannakis. Photovoltaic inverter controllers seeking ac optimal power flow solutions, 2014. [8] T. S. Dillon, K. W. Edwin, H. D. Kochs, and R. J. Taud. Integer programming approach to the problem of optimal unit commitment with probabilistic reserve determination. IEEE Trans. on Power Apparatus and Systems, 97(6):2154 – 2166, 1978. [9] A. Frangioni and C. Gentile. A computational comparison of reformulations of the perspective relaxation: Socp vs. cutting planes. Operations Research Letters, 37(3):206 – 210, 2009. [10] A. Frangioni, C. Gentile, and F. Lacalandra. Tighter approximated milp formulations for unit commitment problems. IEEE Trans. on Power Systems, 24(1):105 – 113, 2009. [11] Antonio Frangioni. Solving nonlinear single-unit commitment problems with ramping constraints. Operations Research, 54(4):767 – 775, 2006. [12] Stephen Frank and Steffen Rebennack. A Primer on Optimal Power Flow: Theory,

15

Formulation, and Practical Examples. Working Papers 2012-14, Colorado School of Mines, Division of Economics and Business, October 2012. [13] Y. Fu and M. Shahidehpour. Fast scuc for large-scale power systems. IEEE Trans. on Power Systems, 22(4):2144 – 2151, 2007. [14] L. L. Garver. Power generation scheduling by integer programming - development of theory. AIE Transactions, 2:730 – 735, 1963. [15] François Glineur, Rue De Houdain, and B-Mons. Computational experiments with a linear approximation of second-order cone optimization, 2000. [16] Z. Gu, G. Nemhauser, and M. Savelsbergh. Lifted flow cover inequalities for mixed 0-1 integer programs. Mathematical Programming, 85:439 –Ű 468, 1999. [17] Z. Gu, G. Nemhauser, and M. Savelsbergh. Sequence independent lifting in mixed integer programming. Journal of Combinatorial Optimization, 4:109 –Ű 129, 2000. [18] X. Guan, P. B. Luh, H. Yan, and J. A. Amalfi. An optimization-based method for unit commitment. International Journal of Electrical Power and Energy Systems, 14(1):9 – 17, 1992. [19] B. F. Hobbs, M. H. Rothkopf, R. P. O’Neill, and H. P. Chao. The Next Generation of Electric Power Unit Commitment Models. Kluwer Academic Publishers, 2001. [20] R. A. Jabr. Tight polyhedral approximation for mixed-integer linear programming unit commitment formulations. IEEE Trans. on Generation, Transmission & Distribution, 6(11):1104 – 1111, 2012. [21] M. Kraning, E. Chu, J. Lavaei, and S. Boyd. Dynamic network energy management via proximal message passing, 2014. [22] Trine Kristoffersen. Stochastic programming with applications to power systems. PhD thesis, 2007. [23] G. S. Lauer, N. R. Sandell, D. P. Bertsekas, and T. A. Posbergh. Solution of large-scale optimal unit commitment problems. IEEE Trans. on Power Apparatus and Systems, 101(1):78 – 86, 1982. [24] P. G. Lowery. Generating unit commitment by dynamic programming. IEEE Trans. on Power Apparatus and Systems, 85(5):422 – 426, 1966. [25] G. Morales, J. M. Latorre, and A. Ramos. Tight and compact milp formulation for the thermal unit commitment problem. IEEE Trans. on Power Systems, 28(4):4897 – 4908, 2013. [26] J. Muckstadt and S. Koening. An application of lagrangean relaxation to scheduling in power-generation systems. Operations Research, 25(3):387 – 403, 1977. [27] G. Nemhauser and L. Wolsey. Integer and Combinatorial Optimization. Wiley, 1988. [28] A. L. Ott. Evolution of computing requirements in the pjm market: Past and future. 16

Power and Energy Society General Meeting, 2010. [29] Zhai Qiaozhu and Wu. Several notes on Lagrangian relaxation for unit commitment. In Proccedings of the 29th Chinese Control Conference, pages 1848–1853, Jul 2010. [30] Ran Quan, Jin bao Jian, and Yun dong Mu. Tighter relaxation method for unit commitment based on second-order cone programming and valid inequalities. International Journal of Electrical Power & Energy Systems, 55(0):82 – 90, 2014. [31] S. Rebennack, P.M. Pardalos, M.V.F. Pereira, and N.A. Iliadis. Handbook of Power Systems I. Springer-Verlag Berlin Heidelberg, 2010. [32] Christos K. Simoglou, Pandelis N. Biskas, and Anastasios G. Bakirtzis. Optimal selfscheduling of a thermal producer in short-term electricity markets by milp. IEEE Trans. on Power Systems, 25(4):1965 – 1977, 2010. [33] W. L. Snyder, H. D. Powell, and J. C. Rayburn. Dynamic programming approach to unit commitment. IEEE Trans. on PWRS-2, 1:339 – 350, 1987. [34] Soliman Abdel-Hady Soliman and Abdel-Aal Hassan Mantawy. Modern Optimization Techniques with Applications in Electric Power Systems. Springer-Verlag New York, 2012. [35] EirGrid & SONI. Solver choice in the sem: A comparative study of lagrangian relaxation vs. mixed integer programming, 2010. [36] C.L. Tseng. On Power System Generation Unit Commitment Problems. PhD thesis, 1996. [37] A. Turgeon. Optimal scheduling of thermal generating units. IEEE Trans. on Automatic Control, 23(6):1000 – 1006, 1978. [38] L. A. Wolsey. Valid inequalities and superadditivity for 0/1 integer programs. Mathematics of Operations Research, 2:65 – 77, 1977. [39] A. J. Wood and B. F. Wollemberg. Power Generation, Operation and Control. Wiley, 2nd. ed. edition, 1996. [40] Xiaohui Yuan, Hao Tian, Shuangquan Zhang, Bin Ji, and Yanhong Hou. Second-order cone programming for solving unit commitment strategy of thermal generators. Energy Conversion and Management, 76(0):20 – 25, 2013.

17

Capítulo 2 Primer artículo: A network flow–based MILP formulation for the thermal unit commitment problem.

18

A MANUSCRIPT SUBMITTED TO THE IEEE TRANSACTIONS ON POWER SYSTEMS

1

A Network Flow–Based MILP Formulation for the Thermal Unit Commitment Problem A. Angulo, Member, IEEE, F. Mancilla–David, Member, IEEE, R. Palma, Senior Member, IEEE, and D. Espinoza Abstract—This letter proposes a reformulation for the thermal unit commitment (UC) problem using a network flow–based monolithic mixed–integer linear programming (MILP) scheme. Most current approaches to solve the UC problem require a tight and compact MILP formulation, making the use of network flow impractical since network flows do not feature compactness. We show in this letter that compactness is not necessary because preprocessing techniques available in state–of–the–art commercial solvers reduce the formulation size. Therefore, network flows may be utilized for the UC problem formulation. In contrast with compact formulations, network flows allow for straightforward inclusion of additional/new constraints, making the problem formulation much more flexible. We evaluate the computational performance of the proposed reformulation against the best compact formulation found in the literature. Results indicate that our approach features overall superior performance. We also discuss a pathway to further improve performance under this new paradigm. Index Terms—Thermal unit commitment, mixed integer linear programming, network flow with side constraints.

I. I NTRODUCTION HE aim of the thermal unit commitment (UC) problem is to decide which of the available generating units in a power system will be used to cover the demand and reserve requirements at a minimum cost. The solution of the UC problem has to satisfy the operational constraints of the generating units, with typical planning horizons ranging from one day to a couple of weeks [1]. The state–of–the–art approach to tackle the UC problem is based on monolithic mixed–integer linear programming (MILP) formulations, which are solved solved via a particular branch & bound (BB) scheme. For example, it has been reported by the PJM interconnection that the MILP–BB approach has been able to improve the UC problem’s solution by about 1% when compared against strategies based on Lagrangian relaxations (LR) [2]. Every percent of improvement translates into millions of dollars of savings [3], thus the UC problem continues to be a topic of active research. Most monolithic MILP formulations found in the literature attempt to make the UC problem’s solution space tight and compact. See, for example, [4], [5] and references within. Tight formulations have a better dual bound and reduce the search space in a BB environment, and compact formulations have smaller size, allowing subproblems associated with nodes to be solved more efficiently. While tightness is indeed required, in this letter we question the assumed need for compactness, because in most of the state– of–the–art solvers a MILP problem entered by a user is first preprocessed to improve quality and reduce its size [6], [7].

T

Manuscript received July 22, 2015. A. Angulo and D. Espinoza are with the Department of Industrial Engineering, University of Chile, Santiago, Chile. F. Mancilla–David is with the Department of Electrical Engineering, University of Colorado Denver, Denver, Colorado 80217, USA R. Palma is with the Department of Electrical Engineering, University of Chile, Santiago, Chile.

As a result of relaxing the need for compactness, we propose reformulating the problem via a network flows approach with side constraints, similar to those utilized in LR–based schemes [8]. From a modeling viewpoint, it is well known that network flows provide the best representation of the set of constraints, because the constraints in this case form a polyhedron defined by a totally unimodular matrix (a polyhedron whose vertices are integer numbers) [9]. Furthermore, a network flows approach allows for straightforward inclusion of additional/new constraints, making the formulation much more flexible than its compact counterparts. The proposed network flow–based MILP strategy is evaluated against the best compact MILP approach found in the literature. Using state–of–the–art solvers, we compare execution time and number of nodes and iterations. In the comparison, we consider power systems of different size, and quantify the computational burden of both the preprocessing and solution stages. Results indicate our approach features overall better computational performance. II. P ROPOSED A PPROACH A typical formulation for the thermal UC problem minimizes the total operational cost (1) subject to system load balance (2), reserve requirements (3) and generating units operational constraints (4), P P min cg (ptg , utg ) (1) t∈T g∈G

s.t.

P

g∈G

P

ptg = Dt ,

g∈G

∀t ∈ T ∀t ∈ T

(3)

∀t ∈ T , ∀g ∈ G

(4)

pg utg ≥ Dt + Rt ,

ptg , utg ∈ Πg ,

(2)

In this formulation, T and G are the set of periods and generating units, respectively, ptg is the output power of unit g at time t, ugt is an on–off binary indicator of unit g at time t, cg (·) is the nonlinear cost function of unit g , Dt is the system demand at time t, Rt is the minimum power system reserve at time t, pg is the maximum available power of unit g and Πg is the feasible region of unit g for all periods. The proposed model uses a network to represent Πg , similar to the state diagram reported in [8], where upper/lower bounds constraints, min–up/min–down constraints and start–up/shut–down constraints (ramping and cost) are included. Ramp–up/ramp– down constraints are modeled as side constraints. Other constraints such as: must run/out, generation levels for the first and last hours, user defined constraints over the state of a single unit, or limits on new type of reserves can be included straightforwardly. Fig. 1 illustrates a typical network for a generating unit. The network generated with this approach is replicated in each period, and hence its computer representation is very simple. Furthermore, the variables for the whole network may be grouped in two sets. First, a set of binary variables associated

19

A MANUSCRIPT SUBMITTED TO THE IEEE TRANSACTIONS ON POWER SYSTEMS

to all edges, and second, a set of continuous, binary bounded, nonnegative variables associated to all edges. It is also possible to use more than one continuous variable per edge if other type of resources are constrained (for example, when modeling various type of reserves is required). It is not difficult to show that all previous tight formulations for thermal UC problem found in the literature, where a polytope is present, can be interpreted in terms of network flows [5], [10]. On, ready to shut-down ramp On, wait for minimum on time

2

are used in a monolithic way and solved with state–of–the–art solvers. Unlike compact MILP models, a network flow–based UC problem formulation has the advantage that the inclusion of new features in the modeling is simple and transparent, removing the need for “tricks” in the construction of models. Indeed, utilizing a network flows approach, the user will always have the best polyhedral representation of the problem, and will be able to easily implement the model in an off–the–shelf solver. The possibility to further reduce computational time may be addressed through a dedicated solver to reduce time in the solution of nodes. Such a solver ought to have warm–start capability for an efficient inclusion in a generic B&B scheme. Efforts in this direction have indeed been reported for the open–pit mine scheduling problem [12]. Performance may also be enhanced via user–specific cuts to strengthen the model and reduce its initial duality gap [13].

On, start-up ramp States

TABLE I C OMPARATIVE EVALUATION OF COMPUTATIONAL PERFORMANCE . On, shut-down ramp

# Instances

Time (s)

Nodes

Iterations

Off, wait for minimum off time

BEST–MILP

All >10 s >100 s >1000 s

156 147 123 65

467.8 675.8 1356.7 3698.8

34035 36325 44031 77424

2086768 2228497 2698523 4601620

Off, ready to start-up ramp

NETWORK

All >10 s >100 s >1000 s

156 147 123 65

501.7 650.5 1067.1 3409.3

24969 26118 30454 59054

1290746 1350760 1572014 2884570

Time

Fig. 1. Typical network model of a generating unit. Tmin−on = 3 s, Tmin−off = 2 s,Tstart−up = 2 s, Tshut−down = 1 and Tcold = 4 s.

III. C OMPUTATIONAL S TUDY AND C ONCLUSIONS To evaluate computational performance, our network–based MILP formulation is evaluated against the best compact MILP found in the literature [5]. In the discussion hereafter, the two formulations are referred to as NETWORK and BEST–MILP, respectively. Both NETWORK and BEST–MILP were programmed in the AMPL environment and optimized using CPLEX’s MILP solver [11]. The study evaluates the 156 instances of the thermal UC problem presented in [5]. These instances include between 28 and 1870 generating units, quadratic generation cost modeled as piecewise linear functions with 1, 2 or 5 steps, start–up costs with 1 or 2 levels, and time horizons between one day and one week with a time seed of one hour. All runs use a single thread with an address space and data limit of 2 GB, 8 hours maximum run time and relative tolerance of 0.01% (CPLEX default). The machines run Linux 2.6.18 under x86–64 architecture, with two quad–core IntelrXeonrE5620 processors and with 48 GB of RAM. Results are summarized in Table I. Run times are computed as geometric means and the counting of nodes and iterations are calculated as arithmetic means. In order to analyze complexity within the results, average values were calculated over a subset of problems, defining three categories based on the time required (10, 100 and 1000 seconds) to achieve a specific tolerance. Table I shows two important general results. First, average times of both formulations are similar. Second, NETWORK outperforms BEST–MILP in more complex instances, achieving the required tolerance in less explored nodes/iterations in all cases. The table also shows that the reduction of iterations is greater than the reduction of explored nodes in the NETWORK approach. Therefore, we conclude that NETWORK features overall better performance than that of BEST–MILP, when they

20

R EFERENCES [1] A. J. Wood, B. F. Wollemberg, and G. B. Shebl´e, Power Generation, Operation and Control, 3rd ed. Wiley, 2013. [2] Z. Qiaozhu and Wu, “Several notes on Lagrangian relaxation for unit commitment,” in Proccedings of the 29th Chinese Control Conference, Jul 2010, pp. 1848–1853. [3] A. L. Ott, “Evolution of computing requirements in the PJM market: Past and future,” in Proceedings of the 2010 IEEE Power and Energy Society General Meeting, Jul 2010, pp. 1–4. [4] M. Carrion and J. M. Arroyo, “A computationally efficient mixed–integer linear formulation for the thermal unit commitment problem,” IEEE Trans. Power Systems, vol. 21, no. 3, pp. 1371–1378, Aug 2006. [5] G. Morales–Espa˜na, J. M. Latorre, and A. Ramos, “Tight and compact MILP formulation for the thermal unit commitment problem,” IEEE Trans. Power Systems, vol. 28, no. 4, pp. 4897–4908, Nov 2013. [6] R. E. Bixby and E. Rothberg, “Progress in computational mixed integer programming—a look back from the other side of the tipping point,” Annals of Operations Research, vol. 149, no. 1, pp. 37–41, Feb 2007. [7] G. Gamrath, T. Koch, A. Martin, M. Miltenberger, and D. Weninger, “Progress in presolving for mixed integer programming,” Konrad–Zuse– Zentrum f¨ur Informationstechnik Berlin (ZIB), Tech. Rep. ZR–13–48, Aug 2013. [8] X. Guan, P. B. Luh, H. Yan, and J. A. Amalfi, “An optimization–based method for unit commitment,” International Journal of Electrical Power and Energy Systems, vol. 14, no. 1, pp. 9–17, Feb 1992. [9] G. Nemhauser and L. Wolsey, Integer and Combinatorial Optimization. Wiley, 1988. [10] D. Rajan and S. Takriti, “Minimum up/down polytopes of the unit commitment problem with start–up costs,” IBM, Tech. Rep. RC23628, Aug 2005. [11] IBM CPLEX optimizer v12.4. [Online]. Available: http://www-01.ibm. com/software/commerce/optimization/cplex-optimizer/ [12] D. Bienstock and M. Zuckerberg. (2009) A new LP algorithm for precedence constrained production scheduling. [Online]. Available: http://www.columbia.edu/∼dano/papers/lpprec.pdf [13] A. Angulo, D. Espinoza, and R. Palma, “Sequence independent, simultaneous and multidimensional lifting of generalized flow covers for the semi– continuous knapsack problem with generalized upper bounds constraints.” in Proceedings the 17th Conference on Integer Programming and Combinatorial Optimization, Jun 2014, pp. 64–75.

Capítulo 3 Segundo artículo: A polyhedral–based approach applied to quadratic cost curves in the unit commitment problem.

21

A MANUSCRIPT SUBMITTED TO THE IEEE TRANSACTIONS ON POWER SYSTEMS

1

A Polyhedral–Based Approach Applied to Quadratic Cost Curves in the Unit Commitment Problem Alejandro Angulo, Member, IEEE, F. Mancilla–David, Member, IEEE, Rodrigo Palma, Senior Member, IEEE, and Daniel Espinoza

Abstract—This paper presents a new polyhedral approximation for the quadratic cost curves of thermal units into the unit commitment problem. The approximation is tight, can be represented by rational numbers, and its size grows logarithmically with respect to the required precision. In addition, the scalability is improved compared to the direct resolution of the quadratic problem involved, due to the linear mixed-integer formulation. On the other hand, the particular structure of the unit commitment problem guarantees that the solution of the mixed–integer linear formulation proposed is a ε-optimal of the original problem. This property allows the direct comparison of the proposed approximation with the original formulation. In most cases, analyzed results indicate that the proposed approach has the best performance compared to other formulations where the bounds of the quadratic cost function are known. Index Terms—Thermal unit commitment, quadratic cost function, mixed integer linear programming, polyhedral approach.

N OMENCLATURE G T pgt ugt Πgt Dt Rt cs cg (.) a, b, c p, p ZLB ZU B ZˆF RGAPF AGAPF

T

Set of generating unit. Set of hourly periods. Power output of unit g in period t. On-off indicator of unit g in period t. Feasible region of unit g in period t. System power demand in period t. Minimum system power reserve in period t. Non-served energy cost. Non linear cost function of unit g. Coefficients of the quadratic production cost function of thermal unit. Minimum/maximum capacity of thermal unit. Lower bound for quadratic formulation of the unit commitment problem. Upper bound for quadratic formulation of the unit commitment problem. Optimal value for formulation F of unit commitment problem. Relative gap for formulation F of unit commitment problem. Absolute gap for formulation F of unit commitment problem. I. I NTRODUCTION

HE study to satisfy a electrical power system’s load demand and to meet the forecast reserve requirements while also satisfying the required technical constraints of Manuscript received July 22, 2015.

22

generation units and minimizing the short term operational cost, i.e. from one day to a couple of weeks, is called Unit Commitment (UC) [1]. From a mathematical point of view, the UC problem is a large scale Mixed-Integer Nonlinear Programming (MINLP) problem, which has well-known complexity [2]. The solution to the UC problem using state of the art software is still not possible when the number of generation units, or periods, increases significantly [3]. Given that the operational costs of thermal units is commonly represented as a quadratic function of the generation level [4], the solution of the UC problem with this type of cost requires more specific techniques such as Mixed-Integer Quadratic Programming (MIQP). Although the solution of such problems with commercial software, such as CPLEX or GUROBI, have significantly improved in recent years, it still does not allow the resolution of real instances. Another technique recently incorporated into commercial software is a tool to solve Mixed-Integer problems with Second Order Cone constraints known as Second Order Cone Programming (SOCP). Due to the fact that the UC problem with a quadratic cost function can be conveniently written as an SOCP using perspective formulation, many studies have proposed solutions using the SOCP method [5]–[7]. This method proved to be more efficient than MIQP to solve UC problems [8], but results show that the improvement made by the SOCP method is not enough to obtain the optimal solution of the UC problem in real instances. From an engineering perspective, the optimal solution of the UC is not always needed and faster approximations have more value to engineers nowadays. Therefore, the solution of the UC problem is open to reliable and efficient linear approximation and using Mixed-Integer Linear Programming (MILP) is a valid alternative with which to solve the UC problem. The MILP technique has the benefit of being a mature technology [9] and its application to the operation of electric systems has generated a positive economic impact [10]. In this modeling context, and in order to avoid the use of interior point algorithms (QP, SOCP) for the node solution in a branch-and-bound environment, the cone problem is represented by its linear semi-infinitive equivalent and its solution is approximated by dynamically adding valid inequalities (perspective cuts) in the Mixed-Integer original domain [3], [8], [11]. After applying this technique the problem has been linearized but its implementation still presents clear drawbacks. One drawback is related to the inefficiencies introduced by the use of some control callback, as cutcallback, to

A MANUSCRIPT SUBMITTED TO THE IEEE TRANSACTIONS ON POWER SYSTEMS

add new cuts to the MILP model. For example, because these controls intervene in the branch & cut search, the presence of a control callback in an application will cause CPLEX to turn off some advanced features as dynamic search or deterministic parallelism. A second is related to the disconnection of the error bounds between the original problem and the proposed linear approximation. Perspective cuts are used in formulations where the quadratic costs are approximated by piecewise continuous functions [4], [12]–[14], but choosing arbitrarily the group of nodes where the valid inequalities are applied. Even though these monolithic formulations solve the inefficiency produced by the use of the cutcallback, they don’t handle the lack of connection with the original problem, hence it is not possible to determine the quality of the solution obtained. The aim of this work is to improve the performance of the linear approximations applied to the UC problem focusing on the issues described in previous paragraphs. This work proposes an efficient polyhedral approximation of the quadratic cost functions, using the approach of Ben-Tal and Nemirovski as a starting point as described in their work of higher dimensional polyhedral relaxation of second order cone constraint [15], [16]. The novel linear formulation obtained in this work has the property of increasing its size logarithmically with the precision required. Previous works in this research area [5] used the approach of Ben-Tal and Nemirovski directly, without proposing new approximations. Moreover, the previous study didn’t analyze the connection between the linear approximation and the original quadratic problem, thus the solution of the UC optimization problem is incomplete. The document is organized as follows: In Section II the general formulation of the UC problem used in this study is presented. Section III explains the procedure to linearize the quadratic costs of the generation units obtaining an efficient monolithic MILP approximation of the UC problem. In this section, the quality of the solution is controlled using consistent adjustable parameters to get an ε-optimal solution of the original problem. In Section IV, some practical implementations of the proposed technique are explored and presented together with the computational results that demonstrate its efficiency. Finally, in Section V the conclusions of this study are given. II. UC P ROBLEM F ORMULATION The UC formulation used in this work is the following: ( ) P P min cg (ptg , utg ) (1) t∈T g∈G

s.t.

P

g∈G

P

ptg = Dt ,

g∈G

∀t ∈ T

pg utg ≥ Dt + Rt ,

ptg , utg ∈ Πgt ,

∀t ∈ T

∀t ∈ T , ∀g ∈ G

(2) (3) (4)

where (1) models the objective function of the system operator, which is to minimize the operational costs of the units;

2

(2) represents the constant balance between generation and demand required in power systems; (3) ensures a minimum reserve in order to adequately handle system contingencies; and (4) constrains the operation of the units within their technical limits. In this formulation, Πgt models the feasible operation range of the units, including constraints such as: power limits, rampup and ramp-down limits, start-up and shut-down ramps and costs, minimum up/down times, must run/out, fixed generation levels, etc. Polyhedral formulations of Πgt are common in the literature [4], [13], [14], [17]–[20] and their quality varies depending on how close the linear relaxation is from its convex hull. This work focuses only on the of Πgt part that is related to the power limits constraints, pu ≤ p ≤ pu (indexes are omitted to simplify notation). On the other hand, cg (.) is a non linear function that represents the generation costs of each unit, which includes both the fuel costs and the start-up costs. In literature, these costs are usually modeled as quadratic functions of the generated power and the operational state of the unit [1], cg (p, u) = ap2 + bp + cu con a > 0. Under these assumptions, the UC problem falls in the category of Mixed-Integer Quadratic Problem (MIQP) and by using simple transformations, it can be converted into a mixedinteger Second Order Cone Problem (SOCP) [5]. The SOCP solution is generally more efficient than the MIQP [8], but in both cases the problem maintains its nonlinearity. Because in realistic problems the information is uncertain (or incomplete) and computers have a finite capacity to represent the problem, exact optimal solutions are not sought. Instead, engineers try to find -´optimal solutions that give satisfactory results within reasonable time. Today, the solution of MILP problems is considered a mature technology where it is possible to find -´optimal solutions in a efficient manner. Thus, it is natural to propose an approximation of the quadratic cost function of the UC problem (which is the only part that cannot be modeled in a monolithic MILP) using a polyhedral set in which quality is controlled only by the  parameter. III. Q UADRATIC C OST F UNCTION A PPROXIMATION The set X o is defined as:  (u, p, w) ∈ {0, 1} × R2 : Xo =

ap2 + bp + cu ≤ w pu ≤ p ≤ pu



(5) which, according to preceding paragraphs, is a substructure of the original UC problem associated to a generator and specific periods. Note that when u = 1, i.e. the unit is in operation, the previous set is the epigraph of the quadratic function ap2 + bp + c in the [p, p] interval. If the following linear transformation is considered p = pu + ∆y, where ∆ = p − p, the original set is converted into   (u, y, w) ∈ {0, 1} × R2 : a ¯y 2 + ¯by + c¯u ≤ w XP = 0≤y≤u (6)

23

A MANUSCRIPT SUBMITTED TO THE IEEE TRANSACTIONS ON POWER SYSTEMS

where a ¯ = a∆2 , ¯b = (2apu + b)∆, and c¯ = ap2 u + bp + c. Note that set X P is nonlinear, but it can be conveniently linearized considering that: since u ∈ {0, 1}, then u = u2 and c¯ can be rewritten as c¯ = ap2 + bp + c; and, since 0 ≤ y ≤ u, then y = yu and ¯b can be rewritten as ¯b = (2ap + b)∆. This transformation is similar to that used in [14], but now the continuous variable was scaled to one. The reason for using this type of transformation is to simplify the formulation, because now you have the same number of variables, but one constraint less. The lower limit constraint of the continuous variable is now a lower bound. Modern LP solvers can handle these bounds more efficiently than other constraints, so its inclusion is desirable for speed of calculation. Finally the X and Y sets are defined as:   (y, z) ∈ R2 : y2 ≤ z X= (7) −1 ≤ y ≤ 1 Y =



(u, y, w, z) ∈ {0, 1} × R3 :

a ¯z + ¯by + c¯u ≤ w 0≤y≤u



(8)

Then the X P set can be expressed as shown below X P = Proj {X ∩ Y }

(9)

(u,y,w)

where X represents the nonlinear convex part of X P , and Y the part that can be directly modeled in a traditional MILP formulation. Note that the X set is the epigraph of the y 2 function in the [−1, 1] interval, and this will be the set used for polyhedral approximation. The method proposed in this work outperforms previous studies in the same research area [3], [5], achieving a better quality polyhedral approximation with strong theoretical results. The lower approximation of the epigraph of the convex function f (y) in the [−1, 1] interval is defined as the intersection of the half-spaces as described by the tangent lines to a set of nodes belonging to the former domain {yi }i=0...,2m , where the interval borders are included. The approximation proposed is clearly polyhedral because it is defined based on a finite number of half-spaces. It will be assumed that the best approximation for this set is the distribution of interior nodes that minimizes the error’s infinity norm within the interval [−1, 1]. This problem can be mathematically expressed as follows: min max yi

y,z

s.t.

f (y) − z

z ≥ f (yi ) + f 0 (yi )(y − yi ), i = 0 . . . 2m 1 ≥ y ≥ −1 yi ≥ yi−1 , i = i = 1 . . . 2m y0 = −1 y2m = 1 (10) where the expressions z ≥ f (yi ) + f 0 (yi )(y − yi ) are the half-spaces defined by the tangent lines of the selected nodes. For the problem where f (y) = y 2 , experiments according to Rebennack et al. show that the best distribution of nodes is equispaced [21]. However, a demonstration is not presented

24

3

nor are error boundaries given. The following theorem complements the results of that work. Theorem 1. The approximation of X with 2m + 1 equispaced supporting nodes where the limits of the interval are included, is the best possible approximation and its precision is  = 1 4m2 . Proof. Consider the optimization problem defined in Equation (10), where f (y) = y 2 . min max yi

y,z

s.t.

y2 − z 2yi y − yi2 ≤ z, i = 0 . . . 2m y0 ≤ y ≤ y2m yi ≥ yi−1 , i = 1 . . . 2m y0 = −1, y2m = 1

where the first two constraints plus the inner objective define a simple problem where the optimal solution for a predefined lattice is the midpoint of the biggest interval: ( 2 ) yi − yi−1 ∗ i = argmax 2 i=1...2m Then, by symmetry, the distance between two consecutive points must be the same, the nodes are equispaced, and the 1 optimal value of optimization problem is 4m 2 , as it was defined in Theorem 1. At this point it is convenient to define the set Qm ⊂ R2 :   (y, z) ∈ R2 : 2yi y − yi2 ≤ z, i = 0, . . . , 2m Qm = −1 ≤ y ≤ 1 (11) which resumes the results of theorem 1 and where yi = −1 + 2 mi for i = 0, . . . , 2m. Note that although the Qm approximation of X is the best achievable, in order to reach a quality of 10−4 , 51 nodes are required (which is equivalent to adding 53 variables/restrictions to the model). If you consider that it must be done many times within the problem, then the size of the problem increases considerably and the possibility opens to study new approaches for more efficient solutions. Based on the ideas introduced by Ben-Tal and Nemirovski in [15] for better polyhedral approximation of the Lorentz cone in R3 , it is possible to improve the result presented in Theorem 1 taking the defined polyhedron as a projection of another polyhedron in a higher dimensional space. Defining the set Pν,τ ⊂ R2ν+4 as:   ξ0 ≥ |y|         η0 = z        ξi ≥ |2ξi−1 − 1| , i = 1, . . . , ν  Pν,τ = η = −4ξ + 4η + 1, i = 1, . . . , ν i i−1 i−1      2     η ≥ 2 j−1 ξ − j−1 , j = 1, . . . , τ     ν ν   τ −1 τ −1     −1 ≤ y ≤ 1 (12) where the integers ν ≥ 0 and τ ≥ 2, are constructive parameters of the set. Just like Qm , this set is clearly a polyhedron, because it is defined by a finite number of linear constraints.

A MANUSCRIPT SUBMITTED TO THE IEEE TRANSACTIONS ON POWER SYSTEMS

4

The following theorem shows the particular structure of this set and its relation with the existing best approximation of X.

2

Theorem 2. The projection of the set Pν,τ in the subspace of the (y, z) variables is equal to the set Qm , where m = (τ − 1)2ν .

ηi

1.5

Proof. To proceed with the demonstration, the following transformations must be taken into account    0   0  y y y ≥ |y| M1 : R2 7→ R2 : 7→ : 0 0 z z z =z    0   0  y y y = 2y − 1 M2 : R2 7→ R2 : 7→ : 0 0 z z z = −4y + 4z + 1

and, using the Fourier-Motzkin Elimination Algorithm, the projection in the (y, z) space corresponds to: ) (    2 j−1 , j = 1, . . . , τ z ≥ ±2 τj−1 −1 y − τ −1 P0,τ = −1 ≤ y ≤ 1

which can be conveniently written using a new indexation as follows:     j   z ≥ 2 −1 + y   τ −1    2 j P0,τ = − −1 + τ −1 , j = 0, . . . , 2(τ − 1)       −1 ≤ y ≤ 1 which demonstrates that P0,τ and Qτ −1 are equivalents.

0.5

0 −1.2

−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

0.2

0.4

0.6

0.8

1

1.2

ξi

(a) 2

ηi−1

1.5

1

0.5

0 −1.2

−1

−0.8

−0.6

−0.4

−0.2

0

ξi−1

(b) 2

1.5

z

which can be understood as a polyhedra within the space they operate. It can be seen that M1 maps the space and folds it with respect to y = 0, transforming one point in the (y, z) space into a line in the right semi-plane of the (y 0 , z 0 ) space; for the case where (y, z) ∈ X, the mapping transforms this region into the set {(y 0 , z 0 ) ∈ R2+ } which can be easily demonstrated using the Fourier-Motzkin Elimination Algorithm. M2 is an affine transformation of the (y, z) space that allows taking the set {(y, z) ∈ X : y ≥ 0} to a region {(y 0 , z 0 ) ∈ X}. To show this, use the inverse affine transformation and replace it. Note that this transformation has particularities regarding the mapping process: it maintains the form of X in the interval of interest and allows the boundary points in (x, y) to still be boundary points in (0 y, z 0 ); and in addition has a fixed point in (1, 1) and moves the (0, 0) point to (−1, 1). Fig. 1 shows the operation of this transformation in the set of interest. Based on this definition it is possible to build Pν,τ using the algorithm shown in Fig. 2. In the eleventh step (11) of the Algorithm 2, τ cuts in the (ξν , ην ) space have been added, as indicated in Theorem 1. However, only in this case the ξν ∈ [0, 1] interval is discretized. For the case ν = 0,   ξ0 ≥ |y|         η0 = z     2 P0,τ = j−1 j−1    ην ≥ 2 τ −1 ξν − τ −1 , j = 1, . . . , τ      −1 ≤ y ≤ 1

1

1

0.5

0 −1.2

−1

−0.8

−0.6

−0.4

−0.2

0

y

0.2

0.4

0.6

0.8

1

1.2

(c) Fig. 1. Effects of affine transformations M1 and M2 on set X. (a) Add cuts in the last iteration, (b) apply (M1 ◦ M2 )−1 in intermediate iterations, and finally (c) apply (M1 )−1 to come back to the original variables.

For cases where ν > 0, the same reasoning is valid if you consider that in this procedure the effect of the M2 + M1 transformation is incorporated, which means that the number of nodes (cuts) is duplicated in the left semiplane during each iteration. To prove this, consider the i-th iteration and assume that there are τi equispaced nodes in the [0, 1] interval. At this (i) stage, the polyhedron Pν,τ is:

(i) Pν,τ

  

ξi ≥ |2ξi−1 − 1| ηi = −4ξi−1 + 4ηi−1 + 1 =    2   ηi ≥ 2 τj−1 ξi − τj−1 , i −1 i −1

  

 j = 1, . . . , τi 

and, using the Fourier-Motzkin Elimination Algorithm again,

25

A MANUSCRIPT SUBMITTED TO THE IEEE TRANSACTIONS ON POWER SYSTEMS

5

its projection in space (ξi−1 , ηi−1 ) is: (   2 )  1 1 1 j−1 1 j−1 ξ − η ≥ 2 ± ± , (i) i−1 i−1 2 2 τi −1 2 2 τi −1 Pν,τ = j = 1, . . . , τi

If this is performed recursively from the defined space (ξν , ην ) to the original space (y, z) adding in the end effects of the transformation M1 (step 1 of the Algorithm 2), the result in the end of the process will be a set of (τ − 1)2ν+1 + 1 cuts in equispaced nodes in the interval [−1, 1] plus the |y| ≤ 1 constraint, which by definition corresponds to the Q(τ −1)2ν polyhedron. Require: Integers ν ≥ 0 and τ ≥ 2. i ← 0. j ← 0.  Pν,τ ← (y, z, ξ, η) ∈ R2ν+4 : |y| ≤ 1 . Ensure: Pν,τ . o Tn M1 1: Pν,τ ← Pν,τ (y, z) −−→ (ξ0 , η0 ) . 2: if ν > 0 then 3: loop 4: i ++ o Tn M2 5: Pν,τ ← Pν,τ (ξi−1 , ηi−1 ) −−→ (ξ˜i−1 , η˜i−1 ) . n o T ˜ M1 6: Pν,τ ← Pν,τ (ξi−1 , η˜i−1 ) −−→ (ξi , ηi ) . 7: if i = ν then 8: break 9: loop 10: j ++    2   T j−1 ξ − 11: Pν,τ ← Pν,τ ην ≥ 2 τj−1 ν −1 τ −1 12:

13:

if j = τ then return Set Pν,τ .

Fig. 2. Algorithm for Pν,τ construction.

Note that the construction process of Pν,τ can be interpreted as the refinement of the initial approximation Q2(τ −1) , through the incorporation of new nodes at the midpoints and iterative application of the procedure ν times. Fig. 3 shows the result of the lower approximation proposed in this work. The same result can be obtained from different pairs of (τ, ν) parameters. The following corollary is derived from the application of previous theorems. Corollary 1. The projections of the set Pν,τ in the sub-space of the (y, z) variables is a polyhedral approximation of X with a precision (ν, τ ) = (τ −1)12 4ν+1 . Proof. Directly from Theorem 1 and 2.

It is now possible to implement a more efficient polyhedral approximation compared to the one obtained when applying Theorem 1 directly. For a tolerance near to 10−4 , the parameters ν = 3 and τ = 8 can be selected, which adds the following variables/constraints to the model: 2(ν + 2) = 10 variables, ν + 1 = 4 equalities and 2(ν + 1) + τ = 16 inequalities. This

26

1.2 1

z

0.8 0.6 0.4 0.2 0 −1

−0.8

−0.6

−0.4

−0.2

0

y

0.2

0.4

0.6

0.8

1

Fig. 3. Outer approximation of X for τ = 2 and ν = 1 (or τ = 3 and ν = 0)

implies that the original number of 59 variables/constraints has decreased to 30, which is nearly half the size of the final approximate problem. Furthermore, since the formulation of Pν,τ has equalities, it is possible to eliminate them and reduce the computational burden even more. A more detailed evaluation of this performance can be done if we calculate the ratio between the problem sizes of monolithic formulation with all cuts added and proposed approach. It is shown in Fig. 4. 103

problem size ratio

which can be conveniently indexed as, (   2 )  j−1 j−1 , ηi−1 ≥ 2 τi−1 (i) −1 ξi−1 − τi−1 −1 Pν,τ = j = 1, . . . , τi−1 = 2(τi − 1) + 1

1.4

τ τ τ τ

102

=2 =4 =6 =8

101

100

10−1 −10 10

10−9

10−8

10−7

10−6

10−5



10−4

10−3

10−2

10−1

100

Fig. 4. Problem size ratio between monolithic formulation with all cuts added and proposed approach

It is observed that for high accuracy tolerances, lower than to 10−3 , a smaller problem is obtained; and for low accuracy tolerances, the additional cost in problem size is negligible. Note that although this result is not directly applied to the quadratic cost functions of a specific unit – ap2 + bp + cu – its extent is simple, but dependent on the model parameters. Using the results of Corollary 1, it is easy to show that a valid bound for the absolute error of this cost function approximation o is given by: a ¯ o (ν, τ, a ¯) = (13) (τ − 1)2 4ν+1 Details of how these local error bounds can be used in determining a valid error bound of the entire quadratic problem

A MANUSCRIPT SUBMITTED TO THE IEEE TRANSACTIONS ON POWER SYSTEMS

will be shown in the next section. In order to evaluate the quality of the approach proposed with other extended polyhedral approximations, it will be compared theoretically with the results of Ben-Tal and Nemirovski [15] applied to set X. For this purpose, the set X is rewritten as  2 2    (y, z, u) ∈ R3 : y 2 + z−u ≤ z+u 2 2 (14) X= u=1   −1 ≤ y ≤ 1

The analysis of the Ben-Tal and Nemirovski approximation method over the X set indicates that in  the nodes are located sin θi π 1 ν ¯+1 yi = 1−cos . Two θi , with θi = 2ν¯ i − 2 and i = 1, . . . , 2 elements of inefficiency can be detected in this approximation. First, the set is no longer equispaced, and second, half of the nodes are located outside the interval [−1, 1]. Hence, those nodes are not used in the polyhedral approximation. Fig. 5 presents the effects of observations made. It shows that with the new approach, both the size of final approximation and the error were reduced by one order of magnitude. 0

10

Ben-Tal & Nemirovski Proposed approach

10−1 10−2

the timelimit set by the user is reached (see Table IV [3]). In addition there is no certainty that the tolerance set for the MILP solver is valid on the original quadratic problem. To show this, a problem in our database has taken, “10 5.dat”, and solved using the formulation presented in [3] for a different number of nodes, where the perspective cuts are applied, and where the MILP solver always uses the same relative tolerance of 10−4 . It is possible to observe from the results shown in Table I that the MILP tolerance does not necessarily ensure that the MIQP tolerance and the required number of nodes is not necessarily low when a high accuracy tolerance is required. This work aims to close this gap by defining a priori a relationship between required MIQP tolerance, number of nodes and consistent MILP solver tolerance configuration. As noted above, this will be discussed in the next section. TABLE I M ONOLITHIC PERSPECTIVE CUTS FORMULATION BEHAVIOR FOR AN MILP SOLVER RELATIVE TOLERANCE OF 10−4 nodes

RGAPM ILP

ZLB

ZU B

RGAPM IQP

4 8 16 32 64

9.9940E-05 9.8056E-05 9.5677E-05 9.2610E-05 9.2506E-05

507821.97 508343.98 508443.74 508463.87 508468.10

508668.65 508555.42 508525.66 508518.80 508517.08

1.6645E-03 4.1577E-04 1.6108E-04 1.0802E-04 9.6318E-05

It is important to clarify that the presented polyhedral approximation can be applied to any convex problem where the quadratic terms appear, but it will only ensure the tolerance in the optimum value if the quadratic terms are in the objective function, which is true in the case of the UC problem.



10−3 10−4 10−5 10−6 10−7

6

1

2

3

4

5

ν

6

7

8

9

IV. C OMPUTATIONAL E XPERIENCE

10

Fig. 5. Theoretical error behavior of outer approximations

Additionally, note that the nodes defined in the Ben-Tal and Nemirovski approximation are not rational numbers, so the inequalities induced by these nodes cannot be exactly represented by computational methods. Any inequality expressed in terms of rational numbers can be scaled and rewritten only with integer coefficients. In our case, the maximum value of these coefficients is given by max{4, (τ − 1)2 } and for typical values of τ between 2 − 10, they have exact binary representation in the computer (an IEEE 754 64-bit floating point number – double in C – can exactly represent integers with an absolute value of less than or equal to 253 ). This means that when the tolerance is small, it could lead to the generation of problems in the numerical conditioning of the matrices involved, with concomitant negative effects in the solution of the linear optimization problem. The novel method uses only rational coefficients, therefore it is less susceptible to this type of problem. Regarding the monolithic approach using perspective cuts as introduced in [3], we can see two things: for lower tolerances (0.5%) the addition of very few cuts is required; and for higher tolerances (0.01%), it is no longer clear and usually

In this section, the performance of the proposed methodology will be evaluated against other techniques. To accomplish this task, different formulations of the UC problem were implemented in AMPL and optimized using the state of the art mixed integer optimization solver CPLEX, version 12.4.0.0, to contrast results. CPLEX is an efficient platform to solve linear and quadratic mixed integer problems. The method used in the optimization process was the branch-and-cut, which combines branch-and-bound with the cutting plane method. The efficient formulation presented in [14] was the starting point to model the UC problem. This formulation includes: production cost, startup cost, power system requirements (energy and spinning reserve), minimum up and down times, generation limits and ramping limits. Three reformulations of this problem, where quadratic cost function bounds are known, were evaluated: 1. Second order cone programming (SOCP): Equivalent to the reformulation studied in [3], which includes perspective approximation of the quadratic terms. This is similar to the mathematics presented in Equation (14), although u is not fixed in one, but varies depending on the state of the unit (binary). 2. Ben-Tal & Nemirovski (BTN): Uses the Ben-Tal & Nemirovski polyhedral approximation in the previous cone formulation. Equivalent to the work presented in [5].

27

A MANUSCRIPT SUBMITTED TO THE IEEE TRANSACTIONS ON POWER SYSTEMS

3. Proposed approach (PAP): Use the polyhedral approximation presented in Theorem 2. With the purpose of improving the robustness of the method, perspective cut ideas were included in this formulation. Thus, the Pν,τ set can be strengthened in the next way:   ξ0 ≥ |y|         η0 = z        ξi ≥ |2ξi−1 − u| , i = 1, . . . , ν  Pν,τ = ηi = −4ξi−1 + 4ηi−1 + u, i = 1, . . . , ν     2      j−1    u, j = 1, . . . , τ  ην ≥ 2 τj−1   −1 ξν − τ −1     0≤y≤u (15) where all the constant coefficients are now multiplied by the status of the unit u. Note that when u = 0 and Pν,τ is included in a UC MILP formulation (where z has a nonnegative objective function coefficient), Pν,τ is equal to the zero vector and z = 0. Also, note that the last τ inequalities applied to (ξν , ην ) space correspond to τ equispaced perspective cuts on interval [0, 1]. Then, when ν = 0, we have a similar formulation to the monolithic perspective cuts model presented by Frangioni in [3]. Along with the newly introduced models, formulations with quadratic terms in the objective function were also evaluated, but, as expected [3] [5], the results were highly inefficient and therefore omitted from this study. Monolithic perspective cuts formulation with equispaced node positions are omitted because it is included in the proposed approach as it indicated in the previous paragraph. In this work, 162 daily instances of the UC problem were considered. Regarding the number of generators involved, instances were divided into three groups: 57 cases with 1020 units (small size instances), 69 cases with 50-75-100 units (medium size instances), and 36 cases with 150-200 units (large size instances). Also, to assess the effect of the required tolerance on the performance of the formulations, all instances were executed for six levels of relative tolerance: 0.5%, 0.25%, 0.1%, 0.05%, 0.025% and 0.01% (CPLEX default). All UC cases were generated using a specific code for this purpose available in [22]. Minor code modifications were made and default values for quadratric coefficient limits of thermal units were used. All generated instances ara available on Daniel Espinoza’s home page. All runs where made using a single thread with an address space limit and data limit of 2GB and 24 hours of running time limit. The machines were running Linux 2.6.18 under x86 64 architecture, with two quad-core IntelrXeonrE5620 processors and with 48GB of RAM. In order to allow fair comparison of the polyhedral formulations BTN and PAP with SOCP, local error bounds show in Corollary 1 are now extended and used to define appropriate values of ν, τ and MILP solver relative tolerance RGAPM ILP to ensure the relative tolerance of the original quadratic problem RGAPM IQP . Let ZˆM ILP , yˆtg and zˆtg , best feasible solution of polyhedral formulation obtained from MILP solver where the relative tolerance has been set to RGAPM ILP . From these results the

28

7

lower bound of the original quadratic problem is calculated as follows: ZLB = ZˆM ILP − AGAPM ILP In a similar way, the upper bound of the original quadratic problem is calculated as: ZU B = ZˆM ILP + S P 2 where S = t∈T g∈G a ¯tg (ˆ ytg − zˆtg ). These bounds are correct because, the best bound obtained from the MILP solver is a valid lower bound of the original quadratic problem, both polyhedral formulations are relaxations of this problem, and the feasible solutions in the polyhedral problem are feasible in the original quadratic problem but have to be corrected due to the outer approximation made in the objective function. Based on these properties, it is possible to calculate the equivalent quadratic relative error of the polyhedral approximations as: P

RGAPM IQP =

AGAPM ILP + S ZˆM ILP + S

(16)

which can be conveniently bounded by taking into account that the absolute and relative gaps of the MILP solver are related as AGAPM ILP = RGAPM ILP × ZˆM ILP , and the worst case error for the lower approximation of set X is bounded by 2 result of Corollary 1 as ytg − ztg ≤ (ν, τ ): RGAPM IQP ≤ RGAPM ILP + (ν, τ )

S¯ ˆ ZM ILP

(17)

P P where S¯ = t∈T g∈G a ¯tg . Let α and β parameters defined as: α≥

(ν, τ ) , RGAPM ILP

β≥

S¯ ˆ ZM ILP

(18)

Note that α can be freely defined by the user, but β must be estimated. This can be performed in the preprocessing stage, solving the linear relaxation of the integer problem for a given (ν, τ ), obtaining some valid lower bound for this problem and using this value to replace ZˆM ILP in the previous expression. Using these definitions and the tightest version of the result shown in (17), it is possible to calculate the relative tolerance RGAPM ILP required by the MILP solver (associated with the polyhedral formulations) and the (ν, τ ) parameters required to approximate the quadratic part of the cost functions, in such a way that the relative tolerance of the original quadratic problem RGAPM IQP is satisfied. The equation that relates both relative tolerances and previous defined parameters is: RGAPM ILP =

RGAPM IQP 1 + αβ

(19)

For the PAP model and τ parameter fixed in 2, νPAP was calculated as: & ' 1+αβ log α RGAP M IQP νPAP = −1 (20) log 4

A MANUSCRIPT SUBMITTED TO THE IEEE TRANSACTIONS ON POWER SYSTEMS

For the BTN model, νBTN parameter was determined as:   π log cos−1 1 α RGAPM IQP +1   1+αβ (21) νBTN =  − 1   log 2  

For all our experiments α = 1 and β = 14 . The results of applying this procedure to the set of relative tolerances considered in this study are shown in TABLE II TABLE II PARAMETERS CONFIGURATION FOR POLYHEDRAL APPROACHES RELTOLSOCP

νBTN

νPAP

0.00500 0.00250 0.00100 0.00050 0.00025 0.00010

5 5 6 6 7 7

3 4 5 5 6 6

All CPLEX parameters in linear and quadratic case were fixed at the default settings with the exception of the relative tolerance mipgap and qcpconvergetol which were fixed at 10−8 .

8

TABLE III RUNNING TIME COMPARISON OF DIFFERENT APPROACHES ON 162 INSTANCES OF UC PROBLEM UNDER DIFFERENT STOPPING CRITERIA RGAPM IQP

SOCP

BTN

PAP

Global

# Inst. 162

0.5% 0.25% 0.1% 0.05% 0.025% 0.01%

99.3 177.6 367.7 528.6 651.3 710.2

79.7 176.6 414.4 590.2 1115.3 1279.9

51.9 116.2 175.8 356.5 421.2 722.3

> 10s

130

0.5% 0.25% 0.1% 0.05% 0.025% 0.01%

223.5 446.1 1112.8 1742.5 2228.3 2489.3

167.6 416.5 1084.4 1677.3 3386.3 4014.3

105.3 270.8 450.8 992.1 1212.4 2126.4

> 100s

101

0.5% 0.25% 0.1% 0.05% 0.025% 0.01%

469.7 1117.7 3523.3 6236.8 8456.7 9695.0

305.2 834.1 2543.4 4379.9 9500.1 11661.2

186.1 545.0 1032.3 2471.2 3166.6 6021.4

> 1000s

76

0.5% 0.25% 0.1% 0.05% 0.025% 0.01%

718.7 1933.2 7810.8 16018.9 23107.8 27133.3

399.2 1133.4 3988.7 7775.7 18590.1 23747.5

243.0 755.6 1595.0 4265.4 5812.4 10698.3

> 10000s

38

0.5% 0.25% 0.1% 0.05% 0.025% 0.01%

1449.0 4142.8 21490.5 44929.5 69984.9 80810.6

587.2 1597.6 6315.7 15095.3 49212.4 68181.2

382.1 1184.0 2850.9 9594.3 15015.8 33539.8

A. Comparison of formulations Two metrics were developed to evaluate the performance of the three formulations described in this work; geometric mean running time and mean processing nodes. In both cases, the analysis was made for different stopping criteria (relative tolerance less than 0.5%, 0.1%, 0.05%, 0.025% and 0.01%) and subsets of the instances where at least one of the algorithms required more than 10, 100, 1000 and 10000 seconds to achieve the lowest tolerance (0.01%). These criteria were selected to capture relative difficulty of the problems for the set of formulations under evaluation. TABLE III and TABLE IV show the metrics for the described criteria. Regarding the time to solve the problems, the results show that the PAP formulation outperforms the results of SOCP and BTN. There was only one case where SOCP was better than PAP, however this improvement was insignificant. Taking into account all instances, the speed-up of PAP is 1.55 and 1.87 regrading SOCP and BTN, respectively. This level of improvement is stable when PAP is compared with BTN, independent of the difficulty of the instance. Contrarily, when PAP is compared with SOCP, an increase in the difficulty reduces the speed of SOCP, showing that the level of speedup in this case increases as the difficulty of the instance increases. This behavior can be explained based on three technical aspects. First, CPLEX uses auxiliary variables to add conic constraints to a model and this procedure increases the size of the problem. Second, the performance of interior point algorithms in a branch and bound environment is not good. They do not have the warm-start capabilities, therefore these algorithms require more time to solve the nodes. Third, the size of the branch-and-bound tree for the conic approach is greater than the linear one mainly due to the weaker cut generation engine.

TABLE IV P ROCESSING NODES COMPARISON OF DIFFERENT APPROACHES ON 162 INSTANCES OF UC PROBLEM UNDER DIFFERENT STOPPING CRITERIA SOCP

BTN

PAP

Global

# Inst. 162

0.5% 0.25% 0.1% 0.05% 0.025% 0.01%

5162 12258 34874 63498 94375 118017

1433 2606 7329 17981 34663 51198

1141 2537 6347 15094 26401 55232

> 10s

161

0.5% 0.25% 0.1% 0.05% 0.025% 0.01%

5194 12334 35091 63892 94961 118750

1442 2622 7375 18092 34879 51516

1148 2553 6387 15188 26565 55575

> 100s

146

0.5% 0.25% 0.1% 0.05% 0.025% 0.01%

5723 13596 38691 70450 104710 130943

1576 2875 8113 19931 38441 56787

1256 2801 7027 16728 29273 61259

> 1000s

112

0.5% 0.25% 0.1% 0.05% 0.025% 0.01%

7353 17603 50299 91687 136334 170518

1902 3526 10348 25731 49848 73735

1482 3478 8974 21549 37883 79548

> 10000s

63

0.5% 0.25% 0.1% 0.05% 0.025% 0.01%

11330 27921 83813 153409 218465 277108

2466 4517 14778 41467 83687 125428

1685 4529 13397 34979 63715 124638

29

RGAPM IQP

A MANUSCRIPT SUBMITTED TO THE IEEE TRANSACTIONS ON POWER SYSTEMS

Regarding the effect of the desired tolerance, solution time results show more erratic behavior, but PAP always outperforms SOCP and BTN. For different levels of difficulties, the best speed-up is always obtained for the same tolerance (for both SOCP 0.1% and BTN 0.025% cases). These results could be related with the approximations made in Equation 20 and 21, and how these approximations affect the relative tolerance used to set the MILP solver of the polyhedral approximations. The results of the mean processing nodes metric is similar in both polyhedral approximations (PAP and BTN). Regarding SOCP, the number of explored nodes triples the number of nodes explored by PAP. This behavior is observed independently of the instance difficulty. Additionally, when the tolerance decreases, a relative reduction is observed in nodes explored by SOCP compared to PAP. As discussed above, a reasonable explanation for this observation relates to the way that CPLEX improves the bounds, where only in the case of linear problems, extensive use of cuts are made. All results point out that monolithic polyhedral approximations (PAP and BTN) have an improved scalability compared with the monolithic quadratic formulation (SOCP), regarding to difficulty and size. This justifies their use in the obtainment of ε-approximations of the UC problem. V. C ONCLUSION In this work, a new polyhedral approximation for the UC problem with a quadratic cost function is presented. This novel method outperforms the only similar polyhedral approximation available, which is based on the work of Ben-Tal and Nemirovski, and it is the direct solution of the quadratic problem reformulated as a second order cone problem. The new polyhedral approximation reduces the time to obtain a solution and also reduces the number of nodes processed in the branch and bound tree. Due to the particular structure of the new formulation, the size of the problem grows logarithmically with the precision required, similar to Ben-Tal and Nemirovski’s approximation. Regarding the representation, the novel approach uses only rational numbers which can be represented exactly with modeling tools. Finally, the proposed approximation based on linear mixed-integer programming assures better scalability compared to the quadratic approximation. The approach developed in this work is general and could be applied to other quadratic convex problems including second order cone constraints, which is proposed as future work. Another improvement considered for future work is the development of a platform that automatically generates this type of approximation based on second order problems. This might be possible based on the structures of the approximations proposed. Nevertheless, the connection between the tolerance desired, the quality of the polyhedral approximation required for this purpose, and the effects of quadratic terms in the objective function must be studied. ACKNOWLEDGMENT The authors acknowledge the financial support from Chilean Government Project FONDECYT 1110024 and Millennium

30

9

Nucleus Information and Coordination in Networks ICM/FIC RC13003. R EFERENCES [1] A. J. Wood and B. F. Wollemberg, Power Generation, Operation and Control, 2nd ed. Wiley, 1996. [2] C. Tseng, “On power system generation unit commitment problems,” Ph.D. dissertation, 1996. [3] A. Frangioni, C. Gentile, and F. Lacalandra, “Tighter approximated MILP formulations for unit commitment problems,” IEEE Trans. on Power Systems, vol. 24, no. 1, pp. 105 – 113, 2009. [4] M. Carrion and J. M. Arroyo, “A computationally efficient mixed-integer linear formulation for the thermal unit commitment problem,” IEEE Trans. on Power Systems, vol. 21, no. 3, pp. 1371 – 1378, 2006. [5] R. A. Jabr, “Tight polyhedral approximation for mixed-integer linear programming unit commitment formulations,” IEEE Trans. on Generation, Transmission & Distribution, vol. 6, no. 11, pp. 1104 – 1111, 2012. [6] X. Yuan, H. Tian, S. Zhang, B. Ji, and Y. Hou, “Second-order cone programming for solving unit commitment strategy of thermal generators,” Energy Conversion and Management, vol. 76, no. 0, pp. 20 – 25, 2013. [7] R. Quan, J. bao Jian, and Y. dong Mu, “Tighter relaxation method for unit commitment based on second-order cone programming and valid inequalities,” International Journal of Electrical Power & Energy Systems, vol. 55, no. 0, pp. 82 – 90, 2014. [8] A. Frangioni and C. Gentile, “A computational comparison of reformulations of the perspective relaxation: SOCP vs. cutting planes,” Operations Research Letters, vol. 37, no. 3, pp. 206 – 210, 2009. [9] R. Bixby, Rothberg, and Z. Gu, “The latest advances in mixed integer programming solvers, presentation notes from the CIRRELT spring school on combinatorial optimization in logistics, montreal,” 2010. [10] A. L. Ott, “Evolution of computing requirements in the PJM market: Past and future,” Power and Energy Society General Meeting, 2010. [11] A. Viana and J. P. Pedroso, “A new MILP–based approach for unit commitment in power production planning,” International Journal of Electrical Power & Energy Systems, vol. 44, no. 1, pp. 997 – 1005, 2013. [12] L. Wu, “A tighter piecewise linear approximation of quadratic cost curves for unit commitment problems,” IEEE Trans. on Power Systems, vol. 26, no. 4, pp. 2581 – 2583, 2011. [13] G. Morales, J. M. Latorre, and A. Ramos, “Tight and compact MILP formulation of start-up and shut-down ramping,” IEEE Trans. on Power Systems, vol. 28, no. 2, pp. 1288 – 1296, 2013. [14] ——, “Tight and compact MILP formulation for the thermal unit commitment problem,” IEEE Trans. on Power Systems, vol. 28, no. 4, pp. 4897 – 4908, 2013. [15] A. Ben-Tal and A. Nemirovski, “On polyhedral approximations of the second-order cone,” Math. Oper. Res., vol. 26, no. 2, pp. 193 – 205, 2001. [16] F. Glineur, R. D. Houdain, and B-Mons, “Computational experiments with a linear approximation of second-order cone optimization,” 2000. [17] J. Lee, J. Leung, and F. Margot, “Min-up/min-down polytopes,” Discrete Optimization, vol. 1, pp. 77 – 85, 2004. [18] D. Rajan and S. Takriti, “Minimum up/down polytopes of the unit commitment problem with start-up costs,” IBM Research Report, 2005. [19] K. Hedman, R. P. O’Neill, and S. S. Oren, “Analizing valid inequalities of the generation unit commitment problem,” Power Systems Conference and Exposition, 2009. PSCE ’09., 2009. [20] J. Ostrowski, M. F. Anjos, and A. Vannelli, “Tight mixed integer linear programming formulations for the unit commitment problem,” IEEE Trans. on Power Systems, vol. 27, no. 1, pp. 39 – 46, 2012. [21] S. Rebbenack and J. Kallrath, “Continious piecewise linear δapproximation for MINLP problems. i. minimal breakpoints systems for univariate functions,” 2012. [22] Unit commitment instance generator. [Online]. Available: http: //www.di.unipi.it/∼frangio/

Capítulo 4 Tercer artículo: Sequence independent, simultaneous and multidimensional lifting of generalized flow covers for the semi–continuous knapsack problem with generalized upper bounds constraints.

31

Sequence independent, simultaneous and multidimensional lifting of generalized flow covers for the semi-continuous knapsack problem with generalized upper bounds constraints Alejandro Angulo1 , Daniel Espinoza1 , and Rodrigo Palma2 ? 1 2

Department of Industrial Engineering, Universidad de Chile Department of Electrical Engineering, Universidad de Chile

Abstract. We consider the semi-continuous knapsack problem with generalized upper bound constraints on binary variables. We prove that generalized flow cover inequalities are valid in this setting. We also prove that, under mild assumptions, they are facet-defining inequalities for the full problem. We then focus on simultaneous lifting of pairs of variables. The associated lifting problem naturally induce multidimensional lifting functions, and we prove that a simple relaxation, in a restricted domain, is a superadditive function. We also prove that in many cases this approximation is actually the optimal lifting function. We then analyze the separation problem, which we separate in two phases: first, find a seed inequality, where we evaluate both exact and heuristic methods; secondly, since the lifting is simultaneous, our class of lifted inequalities might contain an exponential number of them. We choose a strategy of maximizing resulting violation. Finally, we test this class of inequalities on instances arising from electricity planning problems. Our test show that the proposed class of inequalities are strong in the sense that adding a few of these inequalities, they close, on average, 57.70% percent of the root integrality gap, and close 97.70% of the relative gap, while adding very few cuts. Keywords: Knapsack problem, sequence independent multidimensional lifting, generalized upper bounds.

1

Introduction

Binary knapsack programs are a common model for choosing between discrete alternatives. If the choice is continuous but limited, we can see the model as the classical single-node capacitated network design model [10], if the choice is semi-continuous, then we must consider mixed-binary knapsack programs; this problem is known as the semi-continuous knapsack problem (SCKP). If, in addition, we have that binary variables are partitioned in sets such that we can choose at most one of them from any given set, we have what we call a semi-continuous knapsack problem with generalized ?

This research was funded by FONDECYT grant 1110024 and Millennium Nucleus Information and Coordination in Networks ICM/FIC P10-024F

32

upper bound constraints (SCKPGUB). This kind of models are common for representing (possibly non-linear) functions (with only one GUB constraint), or when we are looking at the combined non-linear output of several machines (which is the case on production scheduling problems, as is the case in electricity generation, among many others). In this setting, sequential lifting is too limited, in the sense that whenever we have a constraint y ≤ x for x ∈ {0, 1} and y ∈ [0, 1], lifting must be done first on the integer variable and then on the continuous variable, which precludes finding some facet defining inequalities for the complete problem; thus, simultaneous lifting is essential in this setting. In this paper we use generalized flow cover inequalities [12] (GFC), show that they are valid in our setting, and in many cases induce facets on faces of the original problem. We then propose a valid sequence-independent multidimensional lifting scheme, to obtain valid inequalities for SCKPGUB; we show that superadditivity on a restricted set of feasible right hand sides, and show that this condition is enough to obtain sequenceindependent lifting. We also provide sufficient conditions for this lifting to be maximal. The paper is organized as follows: Section 2 cover some of the known facts on semicontinuous knapsack problems, including some known valid inequalities, and some basic results for the semi-continuous knapsack problem with GUB constraints. Section 3 deals with multidimensional lifting for SCKPGUB, specifically on how to obtain valid sequence-independent lifting functions for them. It also proposes simple algorithms for solving the separation problem; both for the seed inequality, and for selecting maximally violated lifted inequalities. Section 4 present experiments designed to show that the proposed heuristic separation provide good results, and that the lifting step is crucial. Instances are randomly generated and are derived from electricity generation problems. Finally, Section 5 presents our conclusions and further questions on this topic.

2 2.1

Definitions and basic polyhedral results The problem

P To simplify P notation we will write a · b(A) to represent (ai bi : i ∈ A), a(A) to represent (ai : i ∈ A) and [n] to represent {1, . . . , n}. We consider the semi continuous knapsack problem with generalized upper bound constraints, given by  M  (x, y) ∈ {{0, 1} × [0, 1]} :   (a · x + m · y) (M ) ≤ b XG =  yk ≤ xk   x(Mg ) ≤ 1

   

∀k ∈ M    ∀g ∈ G

,

(1)

where G is the set of GUB constraints, M = {(gi , ji )}ni=1 , Mg = {(g 0 , j 0 ) ∈ M : g 0 = g}. Note that ak xk + mk yk is a model for a semi-continuous variable with values in {0} ∪ [ak , ak + mk ]. The first constraint is the semi-continuous knapsack constraint, the second constraint ensure semi-continuity, and the last constraint imposes the generalized upper bound condition among disjoint sets of binary variables.

33

2.2

Literature Review

Many special cases of this structure have been studied before. For example, the classical binary knapsack problem was studied by Balas and Jeroslow [3] in a theoretical work where canonical cuts on the unit hypercube were introduced. Based on this work, in 1975, Wolsey [15] and Balas [2] presented facet defining inequalities for the KP using for the first time the notion of cover. In 1978, Balas and Zemel [4] extend previous work, applying lifting procedures to valid inequalities obtained from minimal covers. In 1980, Padberg [11] present (1, k)-configurations as a generalization of minimal covers inequalities. Inclusion of GUB constraint in KP (KPGUB), was studied by Johnson and Padberg [8]; they also shown how to transform a general instance of the problem in one with non-negative coefficients only. In [14], Wolsey defined some valid inequalities for the KPGUB, and proved that they are facet defining under certain conditions. Sherali and Lee [13], apply sequential and simultaneous lifting to valid inequalities for KPGUB deduced from minimal covers. Another special case is when ak = 0 and |Mg | = 1. This case is called single node flow sets (SNFS) and their study has been extended from work of Gu et al. [6] from lifting procedures applied to this set. In 2007, Louveaux and Wolsey [9] give an interesting survey of strong valid inequalities for knapsack and single node flow sets. As can be seen, application of lifting procedures is a fundamental part of cut generation techniques for many specific sets. In 1977, Wolsey [16] presented the first work on this area where the concept of superadditivity was used. In 2000, Gu et al. [7] generalized it and defined sequence independent lifting of general mixed integer programs. In 2004, Atamt¨urk [1] gave similar results. All these works can be seen as one-dimensional lifting, since they consider the perturbation of one constraint only. For multidimensional lifting, applications are scarce, being the work of Zeng and Richard [18, 19, 17] the most interesting. They defined a general framework to derive multidimensional superadditive lifting functions and applied it to the precedence constrained knapsack problem and to the single node flow set. They show that the traditional concept of superadditivity used by Gu et al. [7], can be restricted depending of the problem in which it is applied. We provide a simple proof of this result in the context of SCKPGUB. 2.3

Polyhedral results

Basic results for SCKPGUB. From this point onward we will assume that a ¯ := max{ak : k ∈ M } < b. With this in place, Proposition 1 follows. Proposition 1. 1. XG is full dimensional. 2. Inequality yk ≥ 0 is facet-defining for XG , ∀k ∈ M . 3. If ak + mk ≤ b, inequality yk ≤ xk is facet-defining for XG , ∀k ∈ M . 4. If ag := 0 max {ag0 j } + min {agj } < b, then x(Ng ) ≤ 1 is facet defining for g 6=g,j∈Ng0

j∈Ng

XG , for g ∈ G.

34

Generalized flow cover inequalities for SCKP. Consider the set    (x, y) ∈ {0, 1}n × [0, 1]n+ :  (a · x + m · y) (N ) ≤ b X= .   yj ≤ xj ∀j ∈ N

(2)

Van Roy and Wolsey [12] studied (a generalization of) this polyhedron, and proposed a family of valid inequalities which they called generalized flow cover inequalities (GFC). In our setting, this family of inequalities can be stated as follows: Given X as defined in (2), a pair (C, CU ) with C ⊂ N , CU ⊂ C, satisfying Γ := a(C) + m(CU ) − b > 0 and mCU > 0, is called a generalized o Defining n cover. ξ

ξj = aj for j ∈ C \ CU , ξj = aj + mj for j ∈ Cu , and γj = min 1, Γj ; then m  γ · x(C) + · (y − x) (CU ) ≤ γ(C) − 1, (3) Γ is valid for X. The following theorem give sufficient conditions for (3) to be facet-defining for X. P Theorem 1. Let (C, CU ) be a generalized flow cover, satisfying mj > Γ . j∈CU :ξj >Γ

Then inequality (3), is facet-defining for Xo := X ∩ {xi = 0, ∀i ∈ / C}.

Note that X is a face of XG , where we choose at most one element from every GUB constraint to be active; from this, the previous theorem give simple conditions for having facet-defining inequalities for this face of XG . Moreover, since X can also be seen as a relaxation of XG , we have that (3) also define valid inequalities for XG .

3 Multidimensional lifting for SCKPGUB 3.1

Valid lifting functions

As was noted before, given a generalized flow cover C, CU in XG satisfying |C ∩ {(g, j) : j ∈ Ng }| ≤ 1, ∀g ∈ G, then (3) is valid for XG , and if m(CU+ ) > Γ , where CU+ = {k ∈ CU : ak + mk > Γ }, then (3) it is a facet-defining inequality for XG ∩ {xi = 0, ∀i ∈ / C}. Following Gu et. al [7], we consider the problem of sequentially lifting pairs of variables3 (xk , yk ), k ∈ / C to obtain m  γ · x(C) + (y − x) (CU ) + (α · x + β · y) (C c ) ≤ γ(C) − 1. (4) Γ n−|C|

n−|C|

If we index pairs {ki }i=1 = {(gi , ji )}i=1 = M \ C, and assume that the first i − 1 pair of variables have been lifted, the i-th lifting functions can be written as hki (z, v) = max αki xki + βki yki s.t. aki xki + mki yki = z xki = vgi 0 ≤ yki ≤ xki 3

And then, simultaneously lift them.

xki ∈ {0, 1},

35

(5)

and the functions fki (z, v) are given by fki (z, v) = min γ · (1 − x)(C) −

m

 (y − x) (CU ) − (α · x + β · y) (K i ) − 1

Γ s.t. (a · x + m · y) (C ∪ K i ) ≤ b − z  x Mg ∩ C ∪ K i ≤ 1 − vg , ∀g ∈ G

0 ≤ yk ≤ xk , xk ∈ {0, 1} ∀k ∈ C ∪ K i ,

(6)

where K i = {k1 , . . . , ki−1 }, z ∈ [0, b] and v has the dimension of the right hand sides of XG for GUB constraints and for simple bounds on binary variables, and is defined as vk0 = δk0 ,ki for k ∈ M and vg0 = δg0 ,gi for g ∈ G, where δa,b = 1 if a = b and zero otherwise, i.e. v = (eki , egi ). We are interested in finding αki , βki that ensure that hki (z, uv) ≤ fki (z, uv) for all (z, u) ∈ {(0, 0)} ∪ {([aki , aki + mki ], 1)} and v = {(eki , egi )}. This implies that we are not interested in hki , fki for all possible (z, v) ∈ R × {0, 1}G∪M , but only on the true domain of feasible points of XG . Although the lifting is multidimensional, at every step, there are only two degrees of freedom in the functions, namely z and u ∈ {0, 1}. The analysis for hki is easy, for the case where mki > 0, the optimal value for hki (z, v) is hki (z, v) =



0 ˜ α ˜ + βz

u = 0, z = 0 u = 1, aki ≤ z ≤ aki + mki

a where α ˜ = αki − mkki βki and β˜ = m1k βki . i i For the case where mki = 0, the optimal value of the function is

hki (z, v) =



0 α ˜

u = 0, z = 0 u = 1, z = aki

where α ˜ = αki and β˜ = βki = 0. α ˜ , β˜ are called normalized lifting coefficients. To study fki we start with a simple case in the following propositions: Proposition 2. Let D = {(g, j) ∈ M \ C : ∃(g, j 0 ) ∈ C, ξgj 0 ≥ Γ, agj + mgj ≤ ξgj 0 − Γ }. Then, for all k ∈ D, the maximal lifting coefficients (αk , βk ) are (0, 0). Proof. Since lifting functions are decreasing in the order, assume that the first elements to lift from the seed inequality are from D. Let k be an element in D. It is known that fk (z, v) ≥ 0 and monotone non-decreasing function, and that f (0, 0) = 0. This implies that is enough to find a feasible point for z = ak + mk with objective value equal to zero to prove our result. Let ko ∈ C satisfying ξko > Γ and ak + mk ≤ ξko − Γ , then, setting (x, y) = (1C − eko , 1C − eko ), we have that fk (z, v) ≤ 0 for z ∈ [ak , ak + mk ], v = (ek , egk ).  Proposition 3. If k ∈ M \ (C ∪ D), then for every optimal solution of the problem fk (z, v), it is always possible to find an optimal solution x∗ , y ∗ satisfying yk∗ = 0 for k ∈ CL and x∗k = yk∗ = 0 for k ∈ D .

36

These two propositions allow us to work assuming that D = ∅ and that m(CL ) = 0. Given v ∈ {0, 1}G , and defining Cv = {(g, j) ∈ C : vg = 1}, then, we can re-write the first lifting function f (z1 , v) as   m f1 (z, v) = γ(C) − 1 − max γ · x + · (y − x) (C \ Cv ) Γ s.t. (ξ · x + m · y) (C \ Cv ) ≤ b − z 0 ≤ yk ≤ xk , xk ∈ {0, 1} ∀k ∈ C \ Cv .

(7)

In general, f1 (z, v) is a complex function; however, in many common cases4 f1 (z, v) is equivalent to (or is bounded from below by) the function obtained from the following relaxation5 : s f˜(z, v) = γ(Cv ) − 1+ min x(C + \ Cv ) + Γ s.t. ξ · x(C + \ Cv ) + s ≥ z + Γ − ξ(Cv ) 0 ≤ s, xk ∈ {0, 1}, ∀k ∈ C + \ Cv ,

(8)

where C + = {k ∈ C : ξk > Γ }. In this form it is easy to prove the following results: v Proposition 4. By renaming C + \ Cv = [rv ] while ensuring that ξhv ≥ ξh+1 , defining v v Λh = ξ(Cv ) + ξ ([h − 1]), and defining H(z) = 0 if z ≤ 0 and H(z) = 1 if z > 0, then s∗ X f˜(z, v) = γ(Cv ) − 1 + + (H (z − s∗ − Λvh + Γ ) : h ∈ [rv ]) , (9) Γ rv P where s∗ = (z − Λvrv+1 + Γ )I(z ≥ Λvrv+1 − Γ ) + (z − Λvh ) I(0 ≤ z − Λvh ≤ Γ ), h=1

moreover, the optimal solution for x is given by x∗h = H(z − s∗ − Λvh + Γ ), ∀h ∈ [rv ]. Theorem 2. The function f˜(z, v) is superadditive for (z, v) ∈ [0, +∞) × {0, 1}G . Corollary 1. If, for each pair of variables (xk , yk ), where k = (g, j), we choose the lifting coefficients (αk , βk ) such that hk (z, u) ≤ f˜(z, ueg ) for (z, u) ∈ {([ak , ak + mk ], 1), (0, 0)}, then the lifting process is sequence independent. 3.2

Algorithmic separation

The previous results show that we can use f˜(z, v) to find valid (and in many cases optimal) lifting coefficients for generalized flow cover inequalities (where, for each GUB, at most one binary variable is in the cover) to obtain strong inequalities for XG . In this section we deal with the separation problem of such lifted inequalities. More precisely, given x∗ a fractional solution in the LP relaxation of (1), try to find a violated lifted constraint. We address this problem in two stages; first, we show how to lift a candidate inequality, then we propose an heuristic to identify a candidate seed inequality. Finally, we keep inequalities that are violated. 4 5

A simple sufficient condition is that the mk ≥ Γ for the two smallest ξk coefficients in C This relaxation is obtained by relaxing integrality of x ∈ / C + , discarding yi ≤ xi , aggregating all continuous variables into s, and complementing remaining integer variables.

37

Lifting GUB-constrained flow-cover inequalities Although closed form expressions for f˜(z, v) are possible; it is important to note that given for each pair of lifted variables k ∈ M \ C and its range [ak , ak + mk ], there are several maximal pairs of coefficients αk , βk satisfying hk (z, xk ) ≤ f˜(z, xk eg(k) ). This implies that the number of possible lifted inequalities derived by this method can be exponential; and a proper method to choose which (set of) inequalities to use is crucial, and should depend on the fractional values of the current fractional point x∗ , y ∗ . Fortunately, a complete description of all pairs of maximal lifting coefficients can be obtained using Algorithm 1, whose complexity is O(|C|).In our implementation we choose (α∗ , β ∗ ) ∈ argmax{x∗k α + yk∗ β : (α, β) ∈ H}, where H is the set of maximal lifting coefficients. Algorithm 1 Finding lower envelope of f˜(·, v). Require: Breakpoints of f˜(·, v), B = {zi }m i=1 where zi ≤ zi+1 and |B| ≥ 2. Interval [a, b], actual range for z (we assume z1 ≤ a ≤ b ≤ zn ). Ensure: H = {α ˜ j , β˜j }, pairs of (normalized) maximal lifting coefficients. 1: B[a,b] ← {a} ∪ {zi ∈ B : a < zi < b} ∪ {b} (ordered set) 2: n ← |B[a,b] |, H ← ∅, kl ← 1, kr ←− 2 3: if n = 1 then 4: return {(f˜(a, v), 0)} 5: loop 6: zl ← B[a,b] [kl ], fl ← f˜(zl , v), zr ← B[a,b] [kr ], fr ← f˜(zr , v) −fl ˜ l 7: β˜ ← fzrr −z ,α ˜ ← fl − βz l 8: if kr + 1 > n then ˜ 9: return H ∪ {(α, ˜ β)} 10: z2r ← B[a,b] [kr + 1], f2r ← f˜(z2r , v) ˜ 2r ≤ f2r then 11: if α ˜ + βz ˜ 12: kl ← kr , H ← H ∪ {(α, ˜ β)} 13: kr ← kr + 1

Algorithm 2 Heuristic to find a GFC. Require: Fractional point (x∗ , y ∗ ). Ensure: (C, CU ), generalized flow cover. 1: C ← ∅, CU ← ∅ 2: for g = 1P to G do 3: zg∗ ← k∈Mg (ak x∗k + mk yk∗ ) ¯g from argmin{k ∈ Mg : min{(ak − zg∗ )+ , (zg∗ − ak − mk )+ }} 4: Select k ∗ 5: if zg > 0 then ¯g } 6: C ← C ∪ {k ∗ 7: if |zg − ak¯g | > |zg∗ − ak¯g − mk¯g | then ¯g } 8: CU ← CU ∪ { k P 9: Apply 1-OPT trying to maximize k∈C γk (x∗k − 1) 10: return (C, CU )

Finding generalized flow cover inequalities in proper faces of XG . It is known that finding maximally violated cover inequalities is already N P-hard [6]; and although it

38

is possible to formulate the separation problem of generalized flow cover inequalities as an IP, we propose Algorithm 26 .

4 Numerical experiments 4.1

Instances

To evaluate the performance of the inequalities presented in this paper, we consider a set of 3, 000 random instances inspired by the unit commitment problem. Furthermore we assume that the GUB structure and semi-continuous variables are already identified. In these instances |G| ∈ {5, 10, 20, 40, 80}, and the number of elements in each GUB constraint where randomly chosen as |Mg | ∼ {U[2, 8], U[7, 13], U[17, 23]}. aj ∼ U[10, 150], mj ∼ U[20, 300], ∀j ∈ M . b ∼ U[0.25, 0.95]bmax , where bmax = is the maximum value of the left hand side of the knapsack constraint. The cost coefficients were chosen as cxk ∼ 2, 500ak − U[370, 1000] − U[15, 50]ak , cyk ∼ 2, 500mk − U[15, 50]mk , ∀k ∈ M , to represent typical cost functions in unit commitment instances. To evaluate the effect of cases with mk = 0, half of the instances include GUB constraints where mk = 0 for 40% of the elements in each GUB constraint. 4.2

Quality measures

We use performance profiles (see [5]) on two quality measures: closed root gap (CG) and closed relative gap (CRG) which we define as CG = 100 ×

zLPn − zLPo , zM IP − zLPo

CRG = 100 ×

zLPn , zM IP

where zM IP is the optimal objective value of the mixed integer problem, zLPo is the optimal objective value of the original linear relaxation and zLPn is the optimal objective value of the final LP relaxation (after several rounds of cuts). Note that for all our instances zLPo < zM IP , and that zM IP > 0; thus, we are never dividing by zero. Also, CRG can be seen as what a user will see as the reported gap when using any of the commercial MIP solvers out there (which might be more interesting for practitioners); while CG can be interpreted as the actual improvement in the lower bound due to the given method (which is a proxy on how much we improve the polyhedral representation of the given set for the given objective function). We do not report running times as the separation process is very quick in all instances and the number of calls of the separation heuristic is always less than fourteen. Moreover, we do not evaluate branch and bound performance, because our instances are exactly instances of a single XG problem, and are always very easy to solve; whereas actual unit commitment problems can have anywhere between 24 to 336 such substructures, in addition to of other side constraints. This is why we chose to leave this study for a future work. 6

This heuristic is a simple extension of other heuristics to find maximally violated cover inequalities

39

4.3

The experiments

The general cutting scheme: In each and every case, we apply the cutting scheme described in Algorithm 3. Algorithm 3 General cutting scheme. Require: LP 0 , initial LP relaxation of XG . Ensure: Z, cut generation scheme optimal values. 1: k ← 0, Z ← ∅ 2: loop 3: Solve current relaxation LP k . 4: Obtain optimal value zk∗ and candidate solution (x∗k , yk∗ ). 5: Z ← Z ∪ {zk∗ } 6: if x∗k ∈ {0, 1}M then 7: return Z 8: From (x∗k , yk∗ ) and using Algorithm 2, find base GFC inequality satisfying Γ ≥ 0.1 (it must be not violated). 9: Lift seed inequality, expressed as in (4), while maximizing resulting violation vk . 10: if Γ vk ≥ 0.1 then 11: k ←k+1 12: Add lifted seed inequality to LP k . 13: else 14: return Z

This scheme can be seen as a basic cutting loop at the root node. We will evaluate the following variations of this scheme: IP : exact separation of seed inequality and no lifting. IP+Lift : same as before, but we perform step 9. Heu : construction heuristic for seed inequality, no lifting. Heu+1-opt : same as before, but perform 1-opt optimization of seed inequality. Heu+1-opt+Lift : same as before, but followed by our lifting step. Heu+Lift : construction heuristic for finding seed inequality, and we perform step 9. Effectiveness of the separation heuristic. Figure 1 shows the performance profile for CG and for CRG on 600 instances with five GUB constraints, where we can solve the IP-separation of the base GFC inequality (top), and shows the performance profile for CG and CRG on all 3, 000 instances (bottom)7 . Table 1 has a summary of these results. It is clear that, measured either by CRG or by CG, Heu+1-opt performs very close to the IP separation of the base heuristic on the set of small instances; while maintaining its edge over the basic heuristic on all instances. This, plus the fact that the exact separation is far too costly on running time, justify using the proposed method; however, this is not an exhaustive evaluation, and any practical implementation should deal on this matter in much more detail. 7

In our case, each point (x, y) of the plotted curves mean that for the worst x% of the instances, the given method closes at most y% of absolute root integrality gap (left), and that the method finished with a lower bound of y% of the actual integer optimal solution value or less (right).

40

100 Heu Heu+1-opt Heu+Lift Heu+1-opt+Lift IP IP+Lift

60

80 zLP (n) % zM IP

80

40

20

0

10

20

30

40 50 60 instances %

70

80

90

0

100

Heu Heu+1-opt Heu+Lift Heu+1-opt+Lift

10

20

30

40 50 60 instances %

70

80

90

100

80 zLP (n) % zM IP

80 60

100 ×

zLP (n)−zLP (0) % zM IP −zLP (0)

0

100

100

40

60 LPo Heu Heu+1-opt Heu+Lift Heu+1-opt+Lift

40 20

20 0

LPo Heu Heu+1-opt Heu+Lift Heu+1-opt+Lift IP IP+Lift

40

20 0

100 ×

60

100 ×

100 ×

zLP (n)−zLP (0) % zM IP −zLP (0)

100

0

10

20

30

40 50 60 instances %

70

80

90

0

100

0

10

20

30

40 50 60 instances %

70

80

90

100

Fig. 1. Left: CG performance profile; Right: CRG performance profile; Top: Instances with |G| = 5; Bottom: All instances. Average CG Average CRG Average Ncuts All |G| = 5 All |G| = 5 All |G| = 5 -Lift +Lift -Lift +Lift -Lift +Lift -Lift +Lift -Lift +Lift -Lift +Lift Heu 15.81 29.63 21.70 35.70 95.01 95.53 82.78 84.56 0.31 1.31 0.31 1.16 Heu+1opt 39.47 57.70 53.73 73.46 96.58 97.70 88.38 92.36 1.64 2.08 1.81 2.14 IP – – 55.81 74.13 – – 89.13 93.11 – – 3.17 3.36 Table 1. Summary results of experiments of average CG, CRG and number of cuts, for both, small and all instances for all algorithm variations.

Robustness of the results: Another question to ask, is how robust are the results on the size of the instances. To answer this, we categorize our instances according to the number of GUB constraints (|G|) and on the number of elements in each GUB constraint (|Mg |), and see the average CG and CRG for Heu+1-opt+Lift. Figure 2 represent a graphical representation of how the average CR and CRG vary depending on these two criteria. Although it is expected that CG performance deteriorate as we increase both the number of GUB constraints and the number of elements in each GUB, is surprising that this tendency is reversed for CRG. This might be due to the special cost structure used in these instances; but, if true on a larger scale, it can be beneficial that the final relative integrality gap decreases on larger instances. The effect of lifting: As was noted in Section 3, our seed inequality is already valid for XG ; so a natural question is how much we gain by doing the lifting process. Again, Table 1 is clear on this respect. If we measure CG, the effect ranges between 15% to

41

20% of more closed root gap; and between 0.5% to 4% of extra CRG; which is very impressive. Moreover, if we look at the number of instances where we could not add any cut; in all the variations of our cutting scheme where lifting was performance, at most in two instances we could not find any cut; while for variations without lifting, we could not find cuts for 2,167 instances for Heu, and 174 instances for Heu+1-opt. All this shows a strong impact of lifting. 73.3%

65.6%

|G| = 5

95.6%

93.6%

87.9%

|G| = 10

68.1%

64.4%

61.8%

|G| = 10

98.3%

97.6%

96.4%

|G| = 20

60.5%

54.8%

57.1%

|G| = 20

99.3%

99.2%

98.9%

|G| = 40

52.1%

51.2%

51.5%

|G| = 40

99.7%

99.7%

99.7%

|G| = 80

44.5%

44.2%

35.9%

|G| = 80

99.9%

99.9%

99.9%

|M

g|

=

10 = g|

g|

|M

g|

|M

|M

=

20 =

10 |M

g|

=

5 = g|

|M

20

82.0%

5

|G| = 5

Fig. 2. Left: average CG on categorized instances; Right: average CRG on categorized instances.

Number of added cuts: A common problem with cutting schemes is that they may require too many cutting rounds to achieve the reported quality. Surprisingly enough, in our experiments, we add, on average 2.14 cuts on all instances; for 92.2% of our instances we add up to three cuts; for 99.33% of our instances we add up to six cuts; and at most 14 cuts in the worst case.

5

Final comments

In this paper we study sequence-independent multidimensional lifting of generalized flow cover inequalities to obtain strong inequalities for the so-called semi continuous knapsack problem with GUB constraints. We also prove that, under mild assumptions, the starting inequality is facet-defining on a face of our polyhedron. Also, under simple assumptions, we show that the sequence-independent lifting function is indeed the optimal (maximal) lifting function; which together with the previous result, allow us to obtain high-dimensional facets. Unlike one-dimensional lifting, in our setting, our supper additive lifting function define a large class of valid inequalities. This introduce the problem of selecting the inequality to be added. In our computational study, we choose the inequality to be added by maximizing resulting violation. We use a set of 3, 000 randomly generated instances of different sizes to conduct our experiments. Our experiments show that, although the separation problem is N P-hard, by using simple heuristics and our superadditive lifting, it is possible to close, on average, 57.70% root integrality gap, and 97.70% relative gap. Finally, there are still many open questions: first, can we take advantage of GUBpartitioned binary variables in other classical polyhedral sets to find tighter valid inequalities? Secondly, in our setting, could we extend our analysis for the case where

42

ak might be negative? This is an important question in our application, but is also relevant in other applications. Thirdly, Is it possible to efficiently detect the basic GUB and semi-continuous structure in general problems? probably not, but even if we are given the GUB constraints, can we use the proposed methodology in general problems? Other relevant questions are also how to select the seed inequality; or should we be looking at using (at the same time) several seed inequalities that could better complement each other when we add them to the current LP relaxation? We feel that all these questions are relevant points for the practical use of the proposed inequalities, and we hope to tackle them soon.

References 1. Atamt¨urk, A.: Sequence independent lifting for mixed-integer programming (2004) 2. Balas, E.: Facets of the knapsack polytope. Mathematical Programming 8, 146 – 164 (1975) 3. Balas, E., Jeroslow, R.: Canonical cuts on the unit hypercube. Mathematical Programming 23, 61 – 69 (1972) 4. Balas, E., Zemel, E.: Facets of the knapsack polytope from minimal covers. SIAM J. Appl. Math. 34, 119 – 148 (1978) 5. Dolan, E.D., Mor´e, J.J.: Benchmarking optimization software with performance profiles. Mathematical programming 91(2), 201–213 (2002) 6. Gu, Z., Nemhauser, G., Savelsbergh, M.: Lifted flow cover inequalities for mixed 0-1 integer programs. Mathematical Programming 85, 439 – 468 (1999) 7. Gu, Z., Nemhauser, G., Savelsbergh, M.: Sequence independent lifting in mixed integer programming. Journal of Combinatorial Optimization 4, 109 – 129 (2000) 8. Johnson, E.L., Padberg, M.W.: A note on the knapsack problem with special ordered sets. Operational Research Letters 1, 18 – 22 (1981) 9. Louveaux, Q., Wolsey, L.A.: Lifting, superadditivity, mixed integer rounding and single node flow sets revisited. Annals OR 153, 47 – 77 (2007) 10. Nemhauser, G., Wolsey, L.: Integer and Combinatorial Optimization. Wiley (1988) 11. Padberg, M.: (1,k)-configurations and facets for packing problems. Mathematical Programming 18, 94 – 99 (1980) 12. Roy, T.V., Wolsey, L.: Valid inequalities for mixed 0-1 programs. Discrete Applied Mathematics 14, 199 – 213 (1986) 13. Sherali, H., Lee, Y.: Sequential and simultaneous lifting of minimal cover inequalities for generalized upper bound constrained knapsack polytopes. SIAM J. Disc. Math. 8, 133 – 153 (1995) 14. Wolsey, L.: Valid inequalities for 0-1 knapsack and mips with generalized upper bound constraints. Discrete Applied Mathematics 29, 251 – 261 (1988) 15. Wolsey, L.A.: Facets of linear inequalities in 0-1 variables. Mathematical Programming 8, 165 – 178 (1975) 16. Wolsey, L.A.: Valid inequalities and superadditivity for 0/1 integer programs. Mathematics of Operations Research 2, 65 – 77 (1977) 17. Zeng, B.: Efficient Lifting Methods for Unstructured Mixed Integer Programs with Multiple Constraints. Ph.D. thesis, Purdue University. Industrial Engineering Department (2007) 18. Zeng, B., Richard, J.P.P.: Sequence independent lifting for 0-1 knapsack problems with disjoint cardinality constraints (2006) 19. Zeng, B., Richard, J.P.P.: A framework to derive multidimensional superadditive lifting functions and its applications. Integer Programming and Combinatorial Optimization. Lecture Notes in Computer Science 4513, 210 – 224 (2007)

43

Capítulo 5 Cuarto artículo: Sequence independent lifting for mixed knapsack problems with GUB constraints.

44

Sequence independent lifting for mixed knapsack problems with GUB constraints Alejandro Angulo1 , Daniel Espinoza1 , and Rodrigo Palma2 1 2



Department of Industrial Engineering, Universidad de Chile Department of Electrical Engineering, Universidad de Chile

Abstract. In this paper, we consider the semi-continuous knapsack problem with generalized upper bound constraints on binary variables. We prove that generalized flow cover inequalities are valid in this setting and, under mild assumptions, are facet-defining inequalities for the entire problem. We then focus on simultaneous lifting of pairs of variables. The associated lifting problem naturally induces multidimensional lifting functions, and we prove that a simple relaxation in a restricted domain is a superadditive function. Furthermore, we also prove that this approximation is, under extra requirements, the optimal lifting function. We then analyze the separation problem in two phases. First, finding a seed inequality, and second, select the inequality to be added. In the first step we evaluate both exact and heuristic methods. The second step is necessary because the proposed lifting procedure is simultaneous; from where our class of lifted inequalities might contain an exponential number of these. We choose a strategy of maximizing the resulting violation. Finally, we test this class of inequalities using instances arising from electrical planning problems. Our tests show that the proposed class of inequalities is strong in the sense that the addition of these inequalities closes, on average, 57.70% of the root integrality gap and 97.70% of the relative gap while adding less than three cuts on average. Keywords: Knapsack problem, sequence-independent multidimensional lifting, generalized upper bounds.

1 Introduction Binary knapsack programs are a common model for choosing between discrete alternatives. If the choice is continuous but limited; the resulting model is called a classical single node flow set, as studied in [19, 17]. If the choice is semi-continuous; we must consider mixed-binary knapsack programs. This problem is known as the semi-continuous knapsack problem (SCKP). If binary variables are subjected to independent clique constraints we have what we call a semi-continuous knapsack problem with generalized upper bound constraints (SCKPGUB). This kind of model is common for representing (possibly non-linear) functions (with only one GUB constraint), or for studying the combined non-linear output of several machines. This is the case in production scheduling problems, such as in electricity generation [26, 6, 9, 18], where the cost function of each generating unit can be highly non⋆

This research was funded by FONDECYT grant 1110024 and Millennium Nucleus Information and Coordination in Networks ICM/FIC RC13003

45

linear and even discontinuous. Furthermore, if the original problem has integer variables; these functions are usually approximated by a piece-wise linear function. In this setting, sequential lifting is too limited in the sense that whenever we have a constraint y ≤ x for x ∈ {0, 1} and y ∈ [0, 1], lifting must be carried out first on the integer variable and then on the continuous variable. This precludes finding some facetdefining inequalities for the complete problem; making simultaneous lifting essential in this setting. For a basic reference on simultaneous lifting, see [27, 13], and see [11, 8, 15] for some experiments and results regarding sequential and simultaneous lifting. In this paper, we use generalized flow cover (GFC) inequalities [21]; show that they are valid in our setting and, under mild assumptions, they induce facets or high dimensional faces of the original problem. We then propose a valid sequence-independent multidimensional lifting scheme to obtain valid inequalities for SCKPGUB. We show that the proposed lifting function is superadditive on a restricted set of feasible righthand sides, and show that this condition is sufficient to obtain sequence-independent lifting. Finally, we also provide sufficient conditions for this lifting to be maximal. The paper is organized as follows. Section 2 covers some known facts about semicontinuous knapsack problems; including a class of valid inequalities and basic results for the semi-continuous knapsack problem with GUB constraints. Section 3 deals with multidimensional lifting for SCKPGUB; specifically on how to obtain valid sequenceindependent lifting functions for them. We also propose simple algorithms to solve the separation problem both for the seed inequality and for the selection of maximally violated lifted inequalities. In Section 4 we propose a heuristic separation algorithm that provides good numerical results and shows that the lifting step is crucial. We generate our test-instances from problems in electricity generation where some of the parameters are randomly perturbed. Finally, we present our conclusions and explore further research questions in Section 5.

2

Definitions and basic polyhedral results

In this section we introduce most of the notation used throughout the paper. This includes the precise definition of the polytope with which we will work. We also present some previously known results and state some basic polyhedral results. Finally, we prove that our seed inequality, under mild hypothesis, is also a facet. 2.1

The problem

We consider the semi-continuous knapsack problem with generalized upper bound constraints given by

XG =

 M (x, y) ∈ {{0, 1} × [0, 1]} :  P    (ak xk + mk yk ) ≤ b  k∈M

P yk ≤ x k xk ≤ 1

    

k∈Mg

46

     

∀k ∈ M    ∀g ∈ G  

,

(1)

where G is the set of GUB constraints, for each g P ∈ G, ng ∈ N is the number of elements in the GUB constraint indexed by g, n = ng , M = {(g, j) : g ∈ G, j ∈ g∈G

1, . . . , ng }, and Mg = {(g ′ , j ′ ) ∈ M : g ′ = g}. This implies that each k ∈ M is an ordered pair (g, j); we will use both notations interchangeably. Also, to simplify r notation, when we Phave a vector µ ∈ R (for some r ∈ N) and R ⊆ {1, . . . , r}, we will denote µ(R) := µi : i ∈ R. Note that ak xk +mk yk is a model for a semi-continuous variable with values in {0} ∪ [ak , ak + mk ]. Although ak and mk could in general be negative, we will focus on the case where both quantities are non-negative, i.e., ak , mk ≥ 0, ∀k ∈ M . Note that unlike in the classical knapsack case, this assumption is restrictive, but we choose it nonetheless because it follows what happens in many applications. The first constraint is the semi-continuous knapsack constraint, the second constraint ensures semi-continuity, and the third constraint imposes a generalized upper bound condition among disjoint sets of binary variables. 2.2

Literature review

Several special subsets of this structure have been studied before. For example, the classical binary knapsack problem (KP) was studied by Balas and Jeroslow [4] in a theoretical study where canonical cuts on the unit hypercube were introduced. Based on this work, in 1975, Wolsey [24] and Balas [3] presented facet-defining inequalities for the KP by using the notion of cover for the first time. Hammer, Johnson and Peled [12] also studied facets of regular 0-1 polytopes, which include knapsack problems. This study also characterizes every non-trivial facet with 0-1 coefficients. In 1978, Balas and Zemel [5] extended previous work by applying lifting procedures to valid inequalities obtained from minimal covers. In 1980, Padberg [20] presented (1, k)-configurations as a generalization of minimal cover inequalities. Johnson and Padberg [14] studied the inclusion of GUB constraints in the KP (KPGUB); they also showed how to transform a general instance of the problem into one with only non-negative coefficients. In 1988, Wolsey [23] defined some valid inequalities for the KPGUB and proved that they are facet-defining under certain conditions. Sherali and Lee [22] applied sequential and simultaneous lifting to valid inequalities for KPGUB deduced from minimal covers. Another special case is when ak = 0 and |Mg | = 1. This case is called single-node flow sets (SNFS), and their study has been extended from the work of Gu et al. [10] from lifting procedures applied to this set. In 2007, Louveaux and Wolsey [16] provided a survey of strong valid inequalities for knapsack and single-node flow sets. As can be seen, the application of lifting procedures is a fundamental part of cut generation techniques for many specific sets. In 1977, Wolsey [25] presented the first work in this area and used the concept of superadditivity. In 2000, Gu et al. [11] generalized it and defined sequence-independent lifting of general mixed integer programs. In 2004, Atamt¨urk [2] presented similar results. The above research can be seen as concerned with one-dimensional lifting, since all of these studies consider the perturbation of only one constraint. Applications of multidimensional lifting are scarce, with work by Zeng [28] and Zeng and Richard [29, 30] as the most relevant. They defined a general framework to derive multidimensional and superadditive lifting functions and applied it to the precedence-constrained knapsack

47

problem and the single-node flow set. They showed that the traditional concept of superadditivity used by Gu et al. [11] can be restricted depending on the problem at hand. We provide a simple proof of this result in the context of SCKPGUB. 2.3

Polyhedral results

Basic results for SCKPGUB. We will henceforth assume that a ¯ := max{ak : k ∈ M } < b. With this in place, Proposition 1 follows:

Proposition 1. 1. XG is full-dimensional. 2. Inequality yk ≥ 0 is facet-defining for XG , ∀k ∈ M . 3. If ak + mk ≤ b, the inequality yk ≤ xk is facet-defining for XG , ∀k ∈ M . P 4. If ag := max {ak } + min {ak } < b, then xk ≤ 1 is facet defining for k∈M \Mg

k∈Mg

k∈Mg

XG , for g ∈ G.

Proof. The basic idea of the proof is to find the appropriate number of feasible affinelyindependent points satisfying each inequality at equality. For details see A.1.  Generalized flow cover inequalities for SCKP. Consider the set   n n    (x, y) ∈P{0, 1} × [0, 1] :  (a x + m y ) ≤ b k k k k X= . j∈M     yk ≤ xk ∀k ∈ M

(2)

Van Roy and Wolsey [21] studied (a generalization of) this polyhedron and proposed a family of valid inequalities that they called generalized flow cover (GFC) inequalities. In our setting, we restate this family of inequalities as follows: Given X as defined in (2), we call a pair (C, CU ), with C ⊂ M and CU ⊂ C, satisfying Γ := a(C) + m(CU ) − b > 0 and m(CU ) > 0, a generalized cover. Then the inequality X mj X ξj (yj − xj ) ≤ −1, (3) min{1, }(xj − 1) + Γ Γ j∈CU

j∈C

where ξj = aj for j ∈ C \ CU and ξj = aj + mj for j ∈ CU , is valid for X. Theorem 1 gives sufficient conditions for (3) to be facet-defining for X. P mj > Γ . Theorem 1. Let (C, CU ) be a generalized flow cover satisfying j∈CU :ξj >Γ

Then, inequality (3) is facet-defining for Xo := X ∩ {xi = 0, ∀i ∈ / C}.

Proof. Van Roy and Wolsey [21] proved the validity of inequality (3). We prove that (3) is facet-defining for Xo := X ∩ {xi = 0, i ∈ / C} by constructing a set of 2s affinely independent points in Xo satisfying it at equality. For details see A.2. 

Note that X is a face of XG where we choose at most one element fromPevery GUB mj > constraint to be active. Given this, Theorem 1 ensures that whenever j∈CU :ξj >Γ

Γ , the resulting flow cover inequality is facet-defining for this face of XG . Moreover, since X is also a relaxation of XG , (3) defines valid inequalities for XG .

48

3 Multidimensional lifting for SCKPGUB In this section we deal with the problem of lifting our seed inequality, defined in (3). To achieve this, we first define what is a valid lifting function in a setting that allows simultaneous and multidimensional lifting in general, and apply it to our particular set. We also identify simple conditions under which optimal lifting coefficients are zero; and introduce a superadditive approximation for maximal lifting functions. Since the full separation of the seed inequality is N P-hard, we propose a heuristic algorithm to find seed inequalities, and describe a simple algorithm to apply our sequence-independent lifting function. 3.1

Valid lifting functions

In this section we study the following problem: given a polytope   Ax ≤ b,   x ∈ Rn : 0 ≤ x ≤ u, P = ,   xi ∈ Z, ∀i ∈ I

where I ⊆ {1, . . . , n}; a set N ⊆ {1, . . . , n}, N c = {1, . . . , n} \ N ; and an inequality ax ≤ bo valid for P , satisfying ai = 0, ∀i ∈ N ; we want to find α ∈ Rn satisfying αi = 0, ∀i ∈ / N such that ax + αx ≤ bo is a valid inequality (and hopefully tight) for P , i.e. ax + αx ≤ bo , ∀x ∈ P

Taking advantage of the condition that ai · αi = 0, ∀i ∈ {1, . . . , n}, we can think that x ∈ Rn has two independent components x = (xN c , xN ) = (v, w) and that P = {(v, w) ∈ Rn : A1 v + A2 w ≤ ¯b, v ∈ V, w ∈ W }, where V, W describe the corresponding box-constraints, integrality requirements, and inequalities involving variables indexed by N c and N respectively. Abusing notation, we can re-state our problem as finding α such that av + αw ≤ bo ,

∀A1 v ≤ ¯b − A2 w, v ∈ V, w ∈ W.

Here, we can make a second observation; the interaction between v and w is only through the common inequalities A1 v ≤ ¯b − A2 w; and there, the interaction is only through values of A2 w : w ∈ W for which there exists v ∈ V satisfying A1 v + A2 w ≤ ¯b. Moreover, the problem of ensuring the condition for all feasible points in P is an optimization problem in a different space. Formally, we state this problem as follows: hα (z)

f (z)

}| { z }| { z max αw min bo − ax s.t. A2 w = z ≤ s.t. A1 v ≤ ¯b − z w∈W v∈V

∀z ∈ Z,

where Z = {z : ∃v ∈ V, w ∈ W, z = A2 w, A1 v + z ≤ ¯b}. With these definitions, our problem can be simply stated as finding α such that hα (z) ≤ f (z), ∀z ∈ Z.

49

Note that in the literature, usually Z ⊆ R; however, in our set, Z ⊆ R × {0, 1}G , and the elements in Z will have a continuous component z and a binary vector component v ∈ {0, 1}G . In this sense, our lifting is multidimensional. Also, since at every step N will have two variables, we will be doing simultaneous lifting. To apply these ideas iteratively, in this section, we re-write (3) as X X mj γ j xj + (yj − xj ) ≤ γ(C) − 1, (4) Γ i∈C

j∈CU

ξ

where γj = min{1, Γj }. As noted before, given a generalized flow cover C, CU in XG satisfying |C ∩ {k : k ∈ Mg }| ≤ 1, ∀g ∈ G; (4) is valid for XG and, if m(CU+ ) > Γ , where CU+ = {k ∈ CU : ak + mk > Γ }, (4) is a facet-defining inequality for XG ∩ {xi = 0, ∀i ∈ / C}. Following Gu et al. [11], we consider the problem of sequentially lifting pairs of variables (xk , yk ), k ∈ / C to obtain X X mj X (yj − xj ) + (αk xk + βk yk ) ≤ γ(C) − 1. (5) γ j xj + Γ i∈C

j∈CU

k∈C /

n−|C|

n−|C|

For simplicity we index pairs {ki }i=1 = {(gi , ji )}i=1 = M \ C and assume that the first i − 1 pairs of variables have been lifted. With this, the i-th lifting function reduces to hki (z, v) = max αki xki + βki yki s.t. aki xki + mki yki = z x ki = v g i 0 ≤ yk i ≤ x ki

(6)

xki ∈ {0, 1}, and fki (z, v) reduces to fki (z, v) = min

X i∈C

− s.t.

γj (1 − xj ) −

X

k∈K i

X

X mj (yj − xj ) Γ

j∈CU

(αk xk + βk yk ) − 1

k∈C∪K i

(ak xk + mk yk ) ≤ b − z

X

k∈Mg (C∪K i )

xk ≤ 1 − vg ,

∀g ∈ G

0 ≤ yk ≤ xk , xk ∈ {0, 1} ∀k ∈ C ∪ K i ,

(7)

where K i = {k1 , . . . , ki−1 }, z ∈ [0, b], and v has the dimension of the right-hand sides of XG for GUB constraints; and define vg′ = δg′ ,gi for g ∈ G; where δa,b = 1 if a = b, and zero otherwise, i.e., v = egi . Our objective is to find αki , βki that ensure that hki (z, uv) ≤ fki (z, uv) for all (z, u) ∈ {(0, 0)} ∪ {([aki , aki + mki ], 1)} and

50

v = egi . This implies that we are not interested in hki , fki for all possible (z, v) ∈ R × {0, 1}G∪M , but only in the true domain of feasible points of XG . Although the lifting is multidimensional, there are only two degrees of freedom at each step in the functions, namely z and u ∈ {0, 1}. The analysis of hki is easy because for the case where mki > 0 the optimal value of hki (z, v) is  0 u = 0, z = 0 hki (z, v) = ˜ α ˜ + βz u = 1, aki ≤ z ≤ aki + mki a where α ˜ = αki − mkki βki and β˜ = m1k βki . i i For the case where mki = 0, the optimal value of the function is  0 u = 0, z = 0 hki (z, v) = α ˜ u = 1, z = aki

where α ˜ = αki and β˜ = βki = 0. We call α ˜ and β˜ normalized lifting coefficients. To study fki , we start with a simple case in the following proposition: Proposition 2. Let D = {(g, j) ∈ M \ C : ∃(g, j ′ ) ∈ C, ξgj ′ ≥ Γ, agj + mgj ≤ ξgj ′ − Γ }. Then, for all k ∈ D, the maximal lifting coefficients (αk , βk ) are (0, 0). Proof. Since fki ≤ fki+1 ; we obtain the best possible lifting coefficients for these variables when we lift them first. We will prove that even in this best case these coefficients are zero. So, we assume that the first elements to lift from the seed inequality are from D. Let k be an element in D. It is known that fk (z, v) ≥ 0 is a monotonic, non-decreasing function, and that f (0, 0) = 0. This implies that it is enough to find a feasible point for z = ak + mk with an objective value equal to zero to prove our result. Let ko ∈ C satisfying ξko > Γ and ak + mk ≤ ξko − Γ . Then, setting (x, y) = (1C − eko , 1C − eko ), we obtain fk (z, v) ≤ 0 for z ∈ [ak , ak + mk ], v = egk .  Note that Theorem 1 ensures that there is nothing to be gained from lifting x and/or y in C; whereas Proposition 2 ensures the same for variables in D. The following proposition will allow us to assume that it is enough to consider the case when D = ∅ and where m(CL ) = 0, where CL := C \ CU . Proposition 3. If k ∈ M \ (C ∪ D), then for every optimal solution of the problem fk (z, v), it is always possible to find an optimal solution x∗ , y ∗ satisfying yk∗ = 0 for k ∈ CL and x∗k = yk∗ = 0 for k ∈ D.

Proof. Let y ∗ , x∗ be an optimal solution to fk (z, v). Note that if x∗ , y ∗ is valid for (7), then for any j ∈ M , changing any x∗j , yj∗ to zero maintain feasibility. Thus, we only need to prove that by making this change for j ∈ D and yj∗ , j ∈ CL , the objective function will not deteriorate. For j ∈ D, the optimal lifting coefficients (αj , βj ) = (0, 0). Then, any valid lifting coefficient pairs αj x∗j +βj yj∗ ≤ 0. Replacing x∗j , yj∗ with (0, 0) will then not deteriorate the objective function; thus proving that there exists an optimal solution with these variables set to zero.

51

For j ∈ CL the argument is similar: if k = 1, note that the objective function has a zero coefficient for yj , from where βj ≤ 0. Using the fact that the lifting functions will be decreasing and that the coefficients in the objective function accompanying yj , for all j ∈ CL , are non-positive; we conclude that we can always find an optimal solution with yj∗ = 0 for j ∈ CL .  These two propositions allow us to work with the assumption that D = ∅ and that m(CL ) = 0. This is because Proposition 2 ensures that the best possible lifting coefficients are zero; while Proposition 3 ensures that, if D 6= ∅, there exists an optimal solution for fki with xi , yi = 0 for all i ∈ D. 3.2

A lower bound for lifting functions

Given v ∈ {0, 1}G , and defining Cv = {(g, j) ∈ C : vg = 1}, we can re-write the first lifting function f (z1 , v) as  X  mk (yk − xk ) f1 (z, v) = γ(C) − 1 − max γ k xk + Γ k∈C\Cv X s.t. (ξk xk + mk (yk − xk )) ≤ b − z k∈C\Cv

0 ≤ yk ≤ xk , xk ∈ {0, 1} ∀k ∈ C \ Cv .

(8)

In general, f1 (z, v) is a complex function, and a simpler functional form f˜ is needed satisfying f˜ ≤ f1 . We propose the following relaxation of f1 : first note that b = ξ(C) − Γ ; replace x by 1 − x and y by x − y; and we obtain the following equivalent form for f1 : X  mk  f1 (z, v) = γ(Cv ) − 1 + min γ k xk + yk Γ k∈C\Cv X s.t. (ξk xk + mk yk ) ≥ z + Γ − ξ(Cv ) k∈C\Cv

0 ≤ yk ≤ xk , xk ∈ {0, 1} ∀k ∈ C \ Cv .

P Now we define C + = {k ∈ C : ξk > Γ }, s = (ξk xk + mk yk ) + k∈(C\C + )\Cv P mk yk , discard the inequality yk ≤ xk and the integrality condition of xk for k∈C + \Cv

k ∈ (C \ C + ) \ Cv to obtain the following relaxation: f˜(z, v) = γ(Cv ) − 1+ min

X

xk +

k∈C + \Cv

X

s.t.

k∈C + \C

s Γ

ξk xk + s ≥ z + Γ − ξ(Cv ) v

0 ≤ s, xk ∈ {0, 1}, ∀k ∈ C + \ Cv .

52

(9)

Note that, under mild conditions3 , f1 (z, v) is equivalent to f˜. Now, we will prove that f˜ has a closed form, and then we will prove that it is also superadditive in an appropriated domain. v Proposition 4. By renaming C + \ Cv = {1, . . . , rv } while ensuring that ξhv ≥ ξh+1 , P v v ξi , and defining H(z) = 0 if z ≤ 0 and H(z) = 1 if z > 0, defining Λh = ξ(Cv ) +

we have

i
r

v s∗ X f˜(z, v) = γ(Cv ) − 1 + + H (z − s∗ − Λvh + Γ ) , Γ

(10)

h=1

where

s∗ = (z − Λvrv+1 + Γ )H(z − Λvrv+1 + Γ )+ rv X

h=1

(z − Λvh ) (H(z − Λvh ) − H(Γ + Λvh − z)) .

Moreover, the optimal solution for x is given by x∗h = H(z − s∗ − Λvh + Γ ), ∀h = 1, . . . , rv . Proof. Using the definition of f˜(z, v) given in (9), if we consider z ≥ Λvrv +1 − Γ ; the solution of the continuous relaxation is also integer-feasible (with x∗i = 1 for all i ∈ C + \ Cv ); from where we have that s∗ = z − Λvrv +1 + Γ and z − s∗ − Λvh + Γ = Λvrv +1 − Λvh > 0; thus proving our result for this case. If we now consider z < Λvrv +1 − Γ and restrict ourselves to solutions where s = 0; the resulting problem has an optimal solution given by xh = H(z + Γ − Λvi ),

∀h = 1, . . . , rv ,

whence (10) is directly obtained. To compute f˜(z, v), consider first a simpler optimization problem: n o s q(z) = min g(z, s) := + aH(z − s) ≥ 0. s≥0 b Note that if s ≥ ab, then g(z, s) ≥ g(z, 0) = aH(z), thus proving that we can restrict ourselves to s ∈ [0, ab]. We now compute q(z) by identifying three cases: Case 1. If z ≤ 0, we have s s + aH(z − s) = 0 < + aH(z − s) ⇒ q(z) = 0, s∗ = 0 b b s=0 0 ba, we have s s + aH(z − s) = a < + aH(z − s) b b s=0 0


q(z) = a, s∗ = 0

A sufficient condition is that mk ≥ Γ for the two smallest ξk coefficients in C

53

because z − s > 0. Case 3. If 0 < z ≤ ba, we can select s∗ = z and we have s s + aH(z − s) = a ≥ + aH(z − s) ⇒ b b s=0 s=s∗

q(z) =

z ∗ , s = z. b

We can now write

f˜(z, v) = γ(Cv ) − 1 + min s≥0

(

) rv X s H (z − s − Λvh + Γ ) , + Γ

(11)

h=1

. Note that  ns o v f˜(z, v) ≥ γ(Cv ) − 1 + min min . + H(z − s − Λh + Γ ) s≥0 h=1 Γ rv

(12)

Using the optimal solution and values computed for q(z); we see that the optimal sov lution {s∗h }rh=1 for the lower bound (12) is either the all-zero vector, or it is the case that exactly one component, say h∗ , is non-zero and is in the range ]0, Γ ]. In the first case, setting s = 0 in (11), we attain the lower bound, and thus solve the problem. In the second case, we have 0 ≤ s∗ = z − Λvh∗ ≤ Γ . Since ξh > Γ , for h 6= h∗ , H(z − Λvh + Γ ) = H(z − s∗h∗ − Λvh + Γ ), thus proving that setting s = s∗h in (11) is a feasible solution that attains the lower bound.  Theorem 2. The function f˜(z, v) is superadditive for (z, v) ∈ [0, +∞) × {0, 1}G .

Proof. Note that f˜(z, 0) is exactly the superadditive function defined in [10]. Namely,  −1 z ≤ −Γ    (z − Λ ) /Γ + i − 1 Λ ∀i = 1, . . . , r − 1 i i − Γ ≤ z ≤ Λi , f˜(z, 0) = i−1 Λi ≤ z ≤ Λi+1 − Γ, ∀i = 1, . . . , r    (z − Λr ) /Γ + r − 1 Λr − Γ ≤ z ≤ b

We propose a different proof for this more general case. We start with f˜(z1 , 0) + f˜(z2 , 0) ≤ f˜(z1 + z2 , 0): v Let {xoh }rh=1 , so be the optimal solution of f˜(z1 +z2 , 0), defined as in Proposition 4, v , s1 be the optimal solution of f˜(z1 , 0). By construction, xo ≥ x1 , and and {x1h }rh=1 assume that h∗o , h∗1 ∈ {0, . . . , rv } are the last active elements in xo , x1 respectively. We prove this by analyzing two cases: Case a ξ · x1 + s1 = Γ + z1 : Define x2h = xoh+h∗ for h ≤ h∗o − h∗1 , x2h = 0 for h > h∗o − h∗1 and 0 ≤ s2 = 1 Γ + so − s1 . By construction and Proposition 4, we have f˜(z1 + z2 , 0) = f˜(z1 , 0) + I · x2 + s2 /Γ − 1, where I is the vector of all ones of the appropriate dimension. Thus, it is enough to prove that (x2 , s2 ) is feasible for (9). However, by hypothesis, we have Γ + z2 ≤ ξ(xo − x1 ) + (so − s1 ) + Γ ≤ ξx2 + s2 , thus proving our result. Case b ξ · x1 + s1 > Γ + z1 : In this case, by optimality, s1 = 0 and ξx1 = Λ0h∗ +1 . Define x2h = xoh+h∗ −1 for 1 1 h ≤ 1 + h∗o − h∗1 , x2h = 0 for h ≥ 2 + h∗o − h∗1 and s2 = so . By construction, we

54

have f˜(z1 + z2 , 0) = f˜(z1 , 0) + I · x2 + s2 /Γ − 1. Thus, is enough to prove that (x2 , s2 ) is feasible for (9). However by hypothesis, we have z1 ≥ Λ0h∗ − Γ , whence 1 Γ + z2 ≤ ξxo + so − z1 ≤ ξ(xo − x1 ) + ξh∗1 + s2 ≤ ξx2 + s2 , thus proving our result. This concludes Case 1. The cases f˜(z1 + z2 , ei ) and f˜(z1 + z2 , ei + ej ) are analogous.



Corollary 1. If, for each pair of variables (xk , yk ), where k = (g, j), we choose lifting coefficients (αk , βk ) such that hk (z, u) ≤ f˜(z, ueg ) for (z, u) ∈ {([ak , ak + mk ], 1), (0, 0)}; then the lifting process is sequence-independent. Proof. We only need to prove validity at any intermediate step r, and call Lr the set of variables (not in C) lifted at step r. That is, we need to prove that max s.t.

P

mk Γ (yk

γ k xk +

k∈C P

k∈C∪L P r

k∈Mgr

r  P (αr xr + βr yr ) − xk ) + k∈Lr

(ξj xj + mj (yj − xj )) ≤ b

xk ≤ 1

(13) ∀g ∈ G ∀k ∈ C ∪ Lr ,

0 ≤ yk ≤ xk , xk ∈ {0, 1}

where Mgr = {(g ′ , j) ∈ C ∪ Lr : g ′ = g} is less than or equal to γ(C) − 1. Note that by hypothesis, αk xk + βk yk = h(ξk xk + mk (yk − xk ), xk ) ≤ f˜(ξk xk + mk (yk − xk ), xk eg(k) ), ∀k ∈ Lr . Using this, Equation (13) can be re-written as  P γk xk + mΓk (yk − xk ) + max s.t.

k∈C r P

r k∈L P

f˜ (ξk xk + mk (yk − xk ), xk g(k))

k∈C∪L P r

k∈Mgr

(14)

(ξj xj + mj (yj − xj )) ≤ b

xk ≤ 1

∀g ∈ G

∀k ∈ C ∪ Lr , P (ξk xk + mk (yk − xk )), where g((g ′ , j ′ )) = g ′ for k = (g ′ j ′ ). By defining zg = k∈Lr P xk , and Mg = {(g ′ , j ′ ) ∈ M : g ′ = g}, we bound (14) above by the δg = 0 ≤ yk ≤ xk , xk ∈ {0, 1}

k∈Lr ∩Mg

following expression: P γ k xk + max s.t.

k∈C r P

mk Γ

f˜ (zg , δg eg )

 (yk − xk ) +

g∈G P

P (ξj xj + mj (yj − xj )) + zg ≤ b g∈G k∈C P x k ≤ 1 − δg ∀g ∈ G

k∈Mg ∩C

0 ≤ zg ≤ bδg , δg ∈ {0, 1} 0 ≤ yk ≤ xk , xk ∈ {0, 1}

55

∀g ∈ G ∀k ∈ C.

(15)

Note that by the definition of C, we have set Go = P|C ∩ Mg | ∈ {0, 1}, ∀g ∈ G. We o {g ∈ G : |C ∩ Mg | = 0}, define zo = zg , and identify C with G \ G . With this, g∈Go

Equation (15) is equivalent to P (γg xg + mg (yg − xg )) max g∈G\Go

+f˜(zo , P 0) + f˜(zg , δg eg ) (ξg xg + mg (yg − xg )) + zg ≤ b s.t. zo + g∈G\Go

xg ≤ 1 − δg 0 ≤ zg ≤ bδg , δg ∈ {0, 1} 0 ≤ yg ≤ xg , xg ∈ {0, 1} 0 ≤ zo ≤ b.

∀g ∈ G \ Go ∀g ∈ G \ Go ∀g ∈ G \ Go

(16)

o Using that δ ∈ {0, 1}G\G and the definition of f˜(z, δ); we can bound (16) as

P max γ(δ) − 1 − f˜( zg , δ)+ g∈{o}∪G\Go P ˜ f (zg , δg ) f˜(zo , 0) + g∈G\Go P s.t. zg ≤ b

(17)

g∈{o}∪G\Go

0 ≤ zg ≤ bδg , δg ∈ {0, 1}

∀g ∈ {o} ∪ G \ Go .

Finally, by using the superadditivity of f˜, we bound (17) by n o o max γ(δ) − 1 : δ ∈ {0, 1}G\G = Γ (C) − 1.

which proves our result. 3.3

(18) 

Algorithmic separation

The foregoing results show that we can use f˜(z, v) to find valid lifting coefficients for GFC inequalities for SCKPGUB; and thus obtain strong inequalities for XG . In this subsection we deal with the separation problem of such lifted inequalities. More precisely, given x∗ as a fractional solution in the linear programming (LP) relaxation of (1), we try to find a violated lifted constraint. We address this problem in two stages: we first show how to lift a candidate inequality, and then propose a heuristic to identify a candidate seed inequality. Finally, we only keep strictly violated inequalities.

Lifting GUB-constrained flow cover inequalities It is important to note that for each pair of lifted variables yk , xk , k ∈ M \ C; there are several maximal pairs of coefficients αk , βk satisfying hk (z, xk ) ≤ f˜(z, xk eg(k) ). This implies that the number of possible lifted inequalities derived by this method can be exponential. Fortunately, Algorithm 1 provides a complete description of all pairs of maximal lifting coefficients;

56

Algorithm 1 Finding the lower envelope of f˜(·, v). Require: Breakpoints of f˜(·, v) and B = {zi }m i=1 where zi ≤ zi+1 and |B| ≥ 2. Interval [a, b], actual range for z (we assume z1 ≤ a ≤ b ≤ zn ). Ensure: H = {α ˜ j , β˜j }, pairs of (normalized) maximal lifting coefficients. 1: B[a,b] ← {a} ∪ {zi ∈ B : a < zi < b} ∪ {b} (ordered set) 2: n ← |B[a,b] |, H ← ∅, kl ← 1, kr ←− 2 3: if n = 1 then 4: H ← {(f˜(a, v), 0)} 5: return H 6: else 7: loop 8: zl ← B[a,b] [kl ], fl ← f˜(zl , v), zr ← B[a,b] [kr ], fr ← f˜(zr , v) l ˜ l ,α ˜ ← fl − βz 9: β˜ ← fzrr −f −zl 10: if kr + 1 ≤ n then 11: z2r ← B[a,b] [kr + 1], f2r ← f˜(z2r , v) ˜ 2r ≤ f2r then 12: if α ˜ + βz ˜ 13: kl ← kr , kr ← kr + 1, H ← H ∪ {(α, ˜ β)} 14: else 15: kr ← k r + 1 16: else ˜ 17: H ← H ∪ {(α, ˜ β)} 18: return H

and its complexity is O(|C|). Moreover, we perform this process once, before performing any lifting step. Figure 1 shows an example for f˜(z, v) and its lower envelope (for a given range [a, b]), which also shows how some key variables of the algorithm change from iteration to iteration. It follows that a proper method to choose the (set of) inequalities to be used is crucial and it should depend on the fractional values of the current fractional point x∗ , y ∗ . If we want to maximize the resulting violation of the lifted inequality, a possible criterion is to choose (α∗ , β ∗ ) ∈ argmax{x∗k α + yk∗ β : (α, β) ∈ H}, where H is the set of maximal lifting coefficients. In our implementation, we choose α∗ , β ∗ as defined before, as long as x∗k α∗ + yk∗ β ∗ > 0; otherwise, we take (α∗ , β ∗ ) ∈ argmax{α + β : (α, β) ∈ H}. Finding generalized flow cover inequalities in proper faces of XG . Finding maximally violated cover inequalities is already N P-hard [10]. Although it is possible to formulate the separation problem of generalized flow cover inequalities as an IP; we propose a simple heuristic described in Algorithm 2, this heuristic is a simple extension of other classical heuristics [11] to find maximally violated cover inequalities. In this heuristic, for each GUB constraint g ∈ G, we compute the net contribution of the current fractional solution x∗ , y ∗ to the knapsack constraint a · x∗ + m · y ∗ ≤ b. We call this net contribution zg . If the current GUB constraint is inactive (zg = 0) we discard it from C. Otherwise, we identify the segment (called k¯g ∈ Mg ), whose boundary is closest to zg , and add it to C. If the closest boundary is the upper limit of the segment we also add that segment to CU . This process can be seen as a greedy construction heuristic. The final step is a local search procedure [1] whose objective is to maximize the likelihood that the resulting inequality will be violated. In our imple-

57

58

ple heuristic methods. We also show the effect of the heuristic separation of the seed inequality. For this, we first describe a wide range of instances; then the set of experiments; and also an analysis of closed gap4 with respect to the size of the problem. 4.1

Instances

To evaluate the performance of the inequalities presented in this paper, we consider a set of 3, 000 random instances inspired by the unit commitment problem. We also assume that the GUB structure and semi-continuous variables are already identified. In these instances, |G| ∈ {5, 10, 20, 40, 80} and the number of elements in each GUB constraint were randomly chosen as |Mg | ∼ {U[2, 8], U [7, 13], U [17, 23]}. aj ∼ U[10, 150], mj ∼ U [20, 300], ∀j ∈ M , b ∼ U[0.25, 0.95]bmax , where bmax = is the maximum value of the left-hand side of the knapsack constraint. We chose the cost coefficients as cxk ∼ 2, 500ak − U [370, 1000] − U [15, 50]ak , cyk ∼ 2, 500mk − U[15, 50]mk , ∀k ∈ M ; which represent typical cost functions in unit commitment instances. To evaluate the effect of having mk = 0; half of the instances contain GUB constraints where mk = 0 for 40% of the elements in each GUB constraint. 4.2

Quality measures

We use performance profiles (see [7]) on two quality measures: closed root gap (CG) and closed relative gap (CRG), which we define as CG = 100 ×

zLPn − zLPo , zM IP − zLPo

CRG = 100 ×

zLPn , zM IP

where zM IP is the optimal objective value of the mixed integer problem; zLPo is the optimal objective value of the original linear relaxation; and zLPn is that of the final LP relaxation. Note that for all our instances, zLPo < zM IP , and zM IP > 0. Thus, when computing CG, we never divide by zero. CRG is an aproximation of the reported gap when using any commercial mixed-integer programming (MIP) solver (which might be more relevant for practitioners). CG is the actual improvement in the lower bound due to the given method (which is a proxy for the extent to which we improve the polyhedral representation of the given set for the given objective function). We do not report running times because the separation process is quick in all instances and the number of calls of the separation heuristic is always less than fourteen. We do not evaluate branch and bound performance because our instances are exactly those of a single XG problem and are always easy to solve; whereas unit commitment problems can have between 24 and 336 sub-structures of this sort in addition to other side constraints. This is why we chose to leave this study for a future work. 4.3

The experiments

The general cutting scheme: In each case, we apply the cutting scheme described in Algorithm 3. 4

i.e. the gap between the linear programming optimal value (with cuts and without the cuts) and the integer programming optimal solution value.

59

Algorithm 3 General cutting scheme. Require: LP 0 , initial LP relaxation of XG . Ensure: Z, cut generation scheme optimal values. 1: k ← 0, Z ← ∅ 2: loop 3: Solve current relaxation LP k . 4: Obtain optimal value zk∗ and candidate solution (x∗k , yk∗ ). 5: Z ← Z ∪ {zk∗ } 6: if x∗k ∈ {0, 1}M then 7: return Z 8: From (x∗k , yk∗ ) and using Algorithm 2, find base GFC inequality satisfying Γ ≥ 0.1 (it must be not violated). 9: Lift seed inequality, as expressed in (5), while maximizing resulting violation vk . 10: if Γ vk ≥ 0.1 then 11: k ←k+1 12: Add lifted seed inequality to LP k . 13: else 14: return Z

This scheme can be seen as a basic cutting loop at the root node. We will evaluate the following variations of this scheme: IP : Separation of GFC seed inequality by solving an integer program that maximizes the violation of the resulting GFC inequality (if any) without lifting. IP+Lift : The same as before, except that we lift the resulting inequality as described in (5). Heu : We execute Algorithm 2 without performing step 9 and use the resulting inequality (i.e., no lifting step is carried out). Heu+1-opt : We execute Algorithm 2 and use the resulting inequality (i.e., no lifting step is performed). Heu+1-opt+Lift : We execute Algorithm 3. Heu+Lift : We execute Algorithm 2 without performing step 9 and lift the resulting inequality as described in (5). Effectiveness of the separation heuristic. Figures 2a and 2b shows the performance profiles of CG and for CRG in 600 instances with five GUB constraints where we can solve the IP-separation of the base GFC inequality. Figures 3a and 3b shows the performance profiles of CG and CRG for all 3, 000 instances. 5 Table 1 lists a summary of these results. From these results it is clear that measured by either CRG or CG, Heu+1-opt performs very close to the IP separation of the base heuristic on the set of small instances while maintaining its edge over the basic heuristic for all instances. This, In addition 5

For Figures 2a and 3a each point (x, y) of the plotted curves means that for the worst x% of the instances, the given method closes at most y% of absolute root integrality gap (left). For Figures 2b and 3b each point (x, y) of the plotted curves means that for the worst x% of the instances, the given method concludes with a lower bound of y% or less of the actual integer optimal solution value (right).

60

to the excessive running time cost of an exact separation routine, justifies the use of the proposed method for evaluation purposes; however, any practical implementation should address this matter in much greater detail. Average CG Average CRG Average Ncuts All |G| = 5 All |G| = 5 All |G| = 5 -Lift +Lift -Lift +Lift -Lift +Lift -Lift +Lift -Lift +Lift -Lift +Lift Heu 15.81 29.63 21.70 35.70 95.01 95.53 82.78 84.56 0.31 1.31 0.31 1.16 Heu+1opt 39.47 57.70 53.73 73.46 96.58 97.70 88.38 92.36 1.64 2.08 1.81 2.14 IP – – 55.81 74.13 – – 89.13 93.11 – – 3.17 3.36

Table 1: Summary results of experiments of the average CG, CRG, and number of cuts, for small and all instances for all algorithm variations.

Robustness of the results: The robustness of the results given the size of the instances is an important consideration. For this, we categorize our instances according to the number of GUB constraints (|G|) and the number of elements in each GUB constraint (|Mg |), and calculate the average CG and CRG for Heu+1-opt+Lift. Figure 4 is a graphical representation of the variation in the average CR and CRG values given these two criteria. Although we expeted expected that CG performance deteriorates as we increase the number of GUB constraints and the number of elements in each GUB; it is surprising that this tendency is reversed for CRG. This might be due to the special cost structure used in these instances. However, if this result holds at a larger scale as well, the fact that the final relative integrality gap decreases can be beneficial. The effect of lifting: As noted in Section 3, our seed inequality is already valid for XG . From this, a natural question is how much do we gain by performing the lifting process? Table 1 clearly represents this aspect. If we measure CG, the effect is a 15%20% larger closed root gap, and 0.5%-4% more CRG. Moreover, in all variations of our cutting scheme where lifting was carried out, there were only two instances where we could not find a cut. For variations without lifting, we could not find cuts for 2,167 instances using Heu and 174 instances when using Heu+1-opt. This shows a strong lifting effect. Number of added cuts: A common problem with cutting schemes is that they may require too many cutting rounds to achieve the desired quality. Surprisingly, in our experiments, we added an average of 2.14 cuts to all instances; in 92.2% of instances, we added up to three cuts; for 99.33% of instances, we added up to six cuts. In the worst case, we added 14 cuts.

5 Final comments In this paper, we studied sequence-independent multidimensional lifting of generalized flow cover inequalities to obtain strong inequalities for the so-called semi-continuous

61

100 Heu Heu+1-opt Heu+Lift Heu+1-opt+Lift IP IP+Lift

100 ×

zLP (n)−zLP (0) % zM IP −zLP (0)

80

60

40

20

0

0

10

20

30

40

50 60 instances %

70

80

90

100

(a)

100

100 ×

zLP (n) % zM IP

80

60 Heu Heu+1-opt Heu+Lift Heu+1-opt+Lift IP IP+Lift LPo

40

20

0

0

10

20

30

40

50 60 instances %

70

80

90

100

(b)

Fig. 2: Top: CG performance profile; Bottom: CRG performance profile; for instances with |G| = 5.

62

100 Heu Heu+1-opt Heu+Lift Heu+1-opt+Lift

100 ×

zLP (n)−zLP (0) % zM IP −zLP (0)

80

60

40

20

0

0

10

20

30

40

50 60 instances %

70

80

90

100

(a)

100

100 ×

zLP (n) % zM IP

80

60

40 Heu Heu+1-opt Heu+Lift Heu+1-opt+Lift LPo

20

0

0

10

20

30

40

50 60 instances %

70

80

90

100

(b)

Fig. 3: Top: CG performance profile; Bottom: CRG performance profile; for all instances.

63

73.3%

65.6%

|G| = 5

95.6%

93.6%

87.9%

|G| = 10

68.1%

64.4%

61.8%

|G| = 10

98.3%

97.6%

96.4%

|G| = 20

60.5%

54.8%

57.1%

|G| = 20

99.3%

99.2%

98.9%

|G| = 40

52.1%

51.2%

51.5%

|G| = 40

99.7%

99.7%

99.7%

|G| = 80

44.5%

44.2%

35.9%

|G| = 80

99.9%

99.9%

99.9% 20 = g|

|M

|M

g|

=

5

|M

|M

g|

g|

=

=

20

10 = g|

|M

= g|

|M

10

82.0%

5

|G| = 5

Fig. 4: Left: average CG for categorized instances; Right: average CRG for categorized instances.

knapsack problem with GUB constraints. We also proved that under mild assumptions, the starting inequality is facet-defining on a face of our polyhedron. Moreover, with simple assumptions, we showed that the sequence-independent lifting function is indeed the optimal (maximal) lifting function which, together with the previous result, enabled us to obtain high-dimensional facets. Unlike one-dimensional lifting, in our setting, the superadditive lifting function defines a large class of valid inequalities. This introduces the problem of selecting the inequality to be added. In our study, we chose the inequality to be added by maximizing the resulting violation. We used a set of 3, 000 randomly generated instances of different sizes to conduct our experiments. These experiments show that although the separation problem is N P-hard, using simple heuristics and superadditive lifting function, it is possible to close, on average, 57.70% of the root integrality gap and 97.70% of the relative gap. There remain several open issues, some of them are: – Can we take advantage of GUB-partitioned binary variables in other polyhedral sets to find tight valid inequalities? – In our setting, can we extend our analysis to the case where ak might be negative? – Is it possible to efficiently detect the basic GUB and the semi-continuous structure in general problems? – Even if we are given the GUB constraints; can we use the proposed methodology in general problems? – Other relevant questions relate to the selection of the seed inequality, and whether we should simultaneously use several seed inequalities that can better complement each other when we add them to the current LP relaxation. – Can we prove that some class of seed inequalities are always dominated by others? – More precisely, should we take only a minimal GFC, or is it better to add small elements to the cover? – Are our results too specific for instances derived from the unit commitment problem? We think that all these questions are relevant to the practical use of the proposed inequalities, and we hope to tackle them in the near future.

64

Acknowledgments We thank the reviewers for their thorough review and highly appreciate the comments and suggestions, which significantly contributed to improving the quality of the publication.

References 1. Aarts, E., Lenstra, J.K. (eds.): Local Search in Combinatorial Optimization. John Wiley & Sons, Inc., New York, NY, USA, 1st edn. (1997) 2. Atamt¨urk, A.: Sequence independent lifting for mixed-integer programming (2004) 3. Balas, E.: Facets of the knapsack polytope. Mathematical Programming 8, 146 – 164 (1975) 4. Balas, E., Jeroslow, R.: Canonical cuts on the unit hypercube. Mathematical Programming 23, 61 – 69 (1972) 5. Balas, E., Zemel, E.: Facets of the knapsack polytope from minimal covers. SIAM J. Appl. Math. 34, 119 – 148 (1978) 6. Carrion, M., Arroyo, J.M.: A computationally efficient mixed-integer linear formulation for the thermal unit commitment problem. IEEE Trans. on Power Systems 21, 1371 – 1378 (2006) 7. Dolan, E.D., Mor´e, J.J.: Benchmarking optimization software with performance profiles. Mathematical programming 91(2), 201–213 (2002) 8. Easton, T., Hooker, K.: Simultaneously lifting sets of binary variables into cover inequalities for knapsack polytopes. Discrete Optimization 5(2), 254 – 261 (2008), in Memory of George B. Dantzig 9. Frangioni, A., Gentile, C., Lacalandra, F.: Tighter approximated milp formulations for unit commitment problems. IEEE Trans. on Power Systems 24, 105 – 113 (2009) 10. Gu, Z., Nemhauser, G., Savelsbergh, M.: Lifted flow cover inequalities for mixed 0-1 integer programs. Mathematical Programming 85, 439 – 468 (1999) 11. Gu, Z., Nemhauser, G., Savelsbergh, M.: Sequence independent lifting in mixed integer programming. Journal of Combinatorial Optimization 4, 109 – 129 (2000) 12. Hammer, P., Johnson, E., Peled, U.: Facet of regular 0–1 polytopes. Mathematical Programming 8(1), 179–206 (1975) 13. Hartvigsen, D., Zemel, E.: The complexity of lifted inequalities for the knapsack problem. Discrete applied mathematics 39(2), 113–123 (1992) 14. Johnson, E.L., Padberg, M.W.: A note on the knapsack problem with special ordered sets. Operational Research Letters 1, 18 – 22 (1981) 15. Kaparis, K., Letchford, A.N.: Cover inequalities. Wiley Encyclopedia of Operations Research and Management Science (2010) 16. Louveaux, Q., Wolsey, L.A.: Lifting, superadditivity, mixed integer rounding and single node flow sets revisited. Annals OR 153, 47 – 77 (2007) 17. Nemhauser, G., Wolsey, L.: Integer and Combinatorial Optimization. Wiley (1988) 18. Ott, A.L.: Evolution of computing requirements in the pjm market: Past and future. Power and Energy Society General Meeting (2010) 19. Padberg, M.W., Roy, T.J.v., Wolsey, L.A.: Valid linear inequalities for fixed charge problems. Operations Research 33(4), pp. 842–861 (1985) 20. Padberg, M.: (1,k)-configurations and facets for packing problems. Mathematical Programming 18, 94 – 99 (1980) 21. Roy, T.V., Wolsey, L.: Valid inequalities for mixed 0-1 programs. Discrete Applied Mathematics 14, 199 – 213 (1986) 22. Sherali, H., Lee, Y.: Sequential and simultaneous lifting of minimal cover inequalities for generalized upper bound constrained knapsack polytopes. SIAM J. Disc. Math. 8, 133 – 153 (1995)

65

23. Wolsey, L.: Valid inequalities for 0-1 knapsack and mips with generalized upper bound constraints. Discrete Applied Mathematics 29, 251 – 261 (1988) 24. Wolsey, L.A.: Facets of linear inequalities in 0-1 variables. Mathematical Programming 8, 165 – 178 (1975) 25. Wolsey, L.A.: Valid inequalities and superadditivity for 0/1 integer programs. Mathematics of Operations Research 2, 65 – 77 (1977) 26. Wood, A.J., Wollemberg, B.F.: Power Generation, Operation and Control. Wiley, 2nd. ed. edn. (1996) 27. Zemel, E.: Lifting the facets of zero-one polytopes. Mathematical Programming 15(1), 268– 277 (1978) 28. Zeng, B.: Efficient Lifting Methods for Unstructured Mixed Integer Programs with Multiple Constraints. Ph.D. thesis, Purdue University. Industrial Engineering Department (2007) 29. Zeng, B., Richard, J.P.P.: Sequence independent lifting for 0-1 knapsack problems with disjoint cardinality constraints (2006) 30. Zeng, B., Richard, J.P.P.: A framework to derive multidimensional superadditive lifting functions and its applications. Integer Programming and Combinatorial Optimization. Lecture Notes in Computer Science 4513, 210 – 224 (2007)

A Extended proofs A.1

Proof of Proposition 1

1. ∀j ∈ M , define ej as the index vector of dimension n with a unique one in position j. Let 0n be the zero vector of size n. Let a be the maximum contribution of any binary variable to the knapsack constraint, i.e. a = max aj and 0 < ǫ < b−a. Then, j∈M

∀k ∈ M , construct pk = (ek , 0n ) and q k = (ek , ǫek ). Define s = (0n , 0n ). The points pk , q k for k ∈ M , and s belong to XG because of Equation (1) and a ¯ < b. These points are affinely independent because ∀k ∈ M , pk − s (n points) and q k − s (n points) are linearly independent. Since we have described 2n + 1 affinely independent points in XG ; we have shown that XG is full-dimensional.  2. To prove that yk ≥ 0 is facet-defining, we use the same definitions as before, but eliminate element q k from the set of valid points. Since we have described 2n affinely independent points in XG satisfying these inequalities at equality; we have shown that these inequalities are facet-defining.  3. To prove that yk ≤ xk is facet-defining; we use the same definitions as in Proposition 1, but eliminate element pk from the set of valid points and redefine q j as (ej , ej ) for all j ∈ M . Again, since we have described 2n affinely independent points in XG and these points satisfy our inequality at equality; we have shown that these inequalities are facet-defining.  4. Note that ag is the maximum contribution to the knapsack constraint when we choose any binary variable and the minimum contribution in Mg at the same time. P This allows us to prove that for each g ∈ G, xk ≤ 1 is facet-defining, by k∈Mg

constructing the points p and q , for each k ∈ M as follows:  (ek + eko , 0n ) ∀k ∈ / Mg k p = (ek , 0n ) ∀k ∈ Mg k

k

66

and k

q =



(ek + eko , ǫek ) (ek , ǫek )

∀k ∈ / Mg ∀k ∈ Mg ,

where ko ∈ argmink∈Mg {ak } and 0 < ǫ < b − ag . Since we have described 2n affinely independent points in XG and these points satisfy the inequality at equality; we have shown that these inequalities are facet-defining.  A.2

Proof of Theorem 1

The validity of inequality (3) was shown by Van Roy and Wolsey [21]. We prove that (3) is facet-defining for Xo := X ∩ {xi = 0, i ∈ / C} by constructing a set of 2s affinely independent points in Xo satisfying it at equality, where s = |C|. For this, define CL = C \ CU , Ck+ = {j ∈ Ck : ξj > Γ }, and C¯k+ = Ck \ Ck+ for k = L, U , and recall that Γ = ξ(C) − b > 0 and our hypothesis is m(CU+ ) > Γ . Let to ∈ argminj∈C + {ξj } and U

assume that 01 = ∞. Note that CU+ 6= ∅, and thus to exists. Let ej be the index vector of dimension s with a unique one in position j, forPall j ∈ C. Let 0s beP the zero vector of size s, 1s be the one vector of size s, 1u+ := j∈C + ej , 1u¯+ := j∈C¯ + ej , and U U 1u = 1u+ + 1u¯+ . Then, construct pj as  (1s − ej , 1u ) if j ∈ CL+     (1 − e , 1 − y 1 + ) if j ∈ C¯L+ s j u oj u j p =  (1s − ej , 1u − ej ) if j ∈ CU+    (1s − ej , 1u − ej − yoj 1 + ) if j ∈ C¯ + u

where yoj =

Γ −ξj + . m(CU )

U

Now, construct q j as

 t o if j ∈ CL   p + (0s , δej ) q j = (1s , 1u − y+ǫ 1u+ + εj ej ) if j ∈ CU+   (1s , 1u − y−ǫ 1u+ − εj ej ) if j ∈ C¯U+

where

0 < δ < min{1, min { j∈CL

ξto − Γ }}, mj

0 < ǫ < min {mj } min{Γ, m(CU+ ) − Γ }/m(CU+ ), j∈CU

y±ǫ =

Γ ±ǫ + , m(CU )

and εj =

ǫ mj j

for j ∈ CU . Given the above definitions, it is easy to

prove that points pj and q belong to Xo . Therefore, it is only necessary to evaluate the knapsack constraint of Equation (1). Verifying the validity of the other constraints is straightforward. Let LHSKN be the left-hand side of the knapsack constraint of Equation (1). Case 1. pj with j ∈ CL+ LHSKN = ξ(C) − aj = b + Γ − aj < b

67

Case 2. pj with j ∈ C¯L+ LHSKN = ξ(C) − aj − yoj m(CU+ ) = b + Γ − aj −

Γ − aj + + m(CU ) = b m(CU )

Case 3. pj with j ∈ CU+ LHSKN = ξ(C) − ξj = b + Γ − ξj < b Case 4. pj with j ∈ C¯U+ LHSKN = ξ(C) − ξj − yoj m(CU+ ) = b + Γ − ξj −

Γ − ξj m(CU+ ) = b m(CU+ )

Case 5. q j with j ∈ CL LHSKN = ξ(C) − ξto + δmj = b + Γ − ξto + δmj < b + Γ − ξto + ξto − Γ = b Case 6. q j with j ∈ CU+ LHSKN = ξ(C) − y+ǫ m(CU+ ) + εj mj = b + εj mj − ǫ = b Case 7. q j with j ∈ C¯U+ LHSKN = ξ(C) − y−ǫ m(CU+ ) − εj mj = b − εj mj + ǫ = b By proceeding in the same manner but now evaluating pj and q j in Equation (3), we can show that these points satisfy this inequality at equality. Remember that the leftP P mj ξ hand side of (3) is min{1, Γj } (xj − 1) + Γ (yj − xj ), and its right-hand side is -1.

j∈C

j∈CU

Case 1. pj with j ∈ CL+ Case 3. pj with j ∈ CU+ Case 5. q j with j ∈ CL LHS = − min{1,

ξj } = −1 Γ

Case 2. pj with j ∈ C¯L+ Case 4. pj with j ∈ C¯U+ LHS = − min{1,

ξj Γ }

− yoj

+ m(CU ) Γ

ξ

= − Γj − (1 −

ξj Γ )

= −1

Case 6. q j with j ∈ CU+ LHS = −y+ǫ

m(CU+ ) mj Γ −ǫ+ǫ + εj =− = −1 Γ Γ Γ

68

Case 7. q j with j ∈ C¯U+ LHS = −y−ǫ

m(CU+ ) mj Γ +ǫ−ǫ − εj =− = −1 Γ Γ Γ

Moreover, these points are affinely independent because ∀j ∈ C \ {to }, pj − pto , and q j − pto are linearly independent. To show this, observe the structure of these points:   j   p − pto , j ∈ CL+ eto − ej eto   pj − pto , j ∈ C¯L+   eto − ej eto − yoj 1 + u  j     p − pto , j ∈ C¯ +   eto − ej eto − ej − yoj 1 +  U u  j     p − pto , j ∈ C + \ {to }   eto − ej eto  U    j  = + t   0s q − p o , j ∈ C  δej L       0s  q j − pto , j ∈ C¯ +  δe j L     + j t   eto  q − p o , j ∈ C¯ eto − εj ej − y−ǫ 1u+  U eto eto + εj ej − y+ǫ 1u+ q j − pto , j ∈ CU+ \ {to }

where, if elements xto and yto are placed in the last position in the columns, it is possible to obtain a upper-triangular matrix (using the first 2s − 1 columns) whose diagonal elements are non-zero. Since 2s affinely independent points in Xo have been described and these points satisfy inequality (3) at equality, it has been shown that these inequalities are facet-defining for conv{Xo }. 

69

Loading...

universidad de chile facultad de ciencias físicas y matemáticas

UNIVERSIDAD DE CHILE FACULTAD DE CIENCIAS FÍSICAS Y MATEMÁTICAS DEPARTAMENTO DE INGENIERÍA INDUSTRIAL OPTIMIZACIÓN LINEAL ENTERA MIXTA APLICADA A PRO...

3MB Sizes 1 Downloads 0 Views

Recommend Documents

UNIVERSIDAD AUSTRAL DE CHILE FACULTAD DE CIENCIAS
Oct 11, 2006 - Palabras claves: Pseudalopex griseus, conflictos carnívoros-humanos, actitudes humanas, daños por fauna .

universidad de chile facultad de ciencias físicas y matemáticas
Block analysis is about rearranging rows and columns of adjacency matrices, and finding communities by creating blocks a

universidad de chile facultad de ciencias físicas y matemáticas
Apr 25, 2011 - La solución desarrollada contempló el uso de tecnología y dispositivos móviles para apoyar el proceso

universidad de chile facultad de ciencias físicas y matemáticas
RESUMEN DE LA MEMORIA PARA OPTAR AL. TÍTULO ... puede afectar negativamente otros usos aguas abajo, incluyendo efectos

universidad de chile facultad de ciencias físicas y matemáticas
de la inflación como la caracterización de tensor B-modes vía el cuociente tensor-a-escalar llamado r, está directam

universidad de chile facultad de ciencias físicas y matemáticas
putación verificable, v) Cliente liviano. El objetivo de esta tesis no es ..... ballot, counting blank votes (or markin

universidad de chile facultad de ciencias fisicas y matematicas
RESUMEN. El yacimiento de cobre Mantos Blancos (500 Mt @ 1.0 % Cu), se ubica en la Cordillera de la ... dominado por pla

Pontificia Universidad Católica de Chile Facultad de Ciencias
Me gustaría agradecer a las secretárias del Departamento de Ecología y de la Dirección del Postgrado de la Facultad

universidad austral de chile facultad de ciencias agrarias
9. 1990, Mbata, 2004). La calidad del alimento esta determinada por las características físicas y químicas de las pla

Facultad de Ciencias - Dadun - Universidad de Navarra
Fisiología del Estrés en Plantas (Departamento de Biología Ambiental), Unidad .... 78. 3.4.2. Volumen de la baya …