Optimización de Procesos

Vicente Rico Ramírez

Departamento de
Ingeniería Química
Instituto Tecnológico de
Celaya
 Índice de Contenido
1 Introducción
1.1 La Ingeniería de Procesos
1.2 Modelación y Grados de Libertad
1.3 Representación Matemática Generalizada de un Problema de Optimización
1.4 Tipos de Problemas de Optimización
1.5 Región Factible
1.6 Convexidad

2. Técnicas de Optimización
2.1 Programación Lineal: Método Simplex
2.2 Programación No Lineal
2.2.1 Optimización sin restricciones
2.2.2 Optimización con Restricciones de Igualdad
2.2.3 Optimización con Restricciones de Desigualdad
2.3 La Programación Mixta-Entera en el Diseño de Procesos
2.4 Programación Mixta Entera Lineal: Método de “Branch and Bound”
2.5 Programación Mixta-Entera No Lineal: Método “Outer Approximation”

3. El Ambiente de Modelación GAMS y sus Resolvedores

 Índice de Contenido
4. Aplicaciones en Ingeniería Química

5. Introducción a la Optimización Bajo Incertidumbre

5.1 Tipos de Problemas de Programación Estocástica
5.2 El Método de Descomposición Estocástica

6. Introducción a la Optimización Multiobjetivo

6.1 Métodos de Solución

7. Control Óptimo y Optimización Dinámica

7.1 El Principio del Máximo
7.2 Programación Dinámica
7.3 Programación Dinámica Estocástica
3 Etapas en la Ingeniería de Procesos
”Process Systems Engineering”

• Síntesis (o Diseño)

• Simulación (o Análisis)

• Optimización
Síntesis (o Diseño ) de Procesos

Materia Productos
Prima Bajo
(condiciones Especificación
iniciales)

Determinación de la Estructura del Proceso

para realizar la transformación deseada
Análisis (o Simulación ) de Procesos

Materia Productos
Prima Bajo
(condiciones Especificación
iniciales)

Dadas las condiciones de entrada y la estructura

del proceso, determinar las variables de salida
 Optimización de Procesos

Minimizar Costo
D, Pureza

Alimentación N=?
R=?
P=?

Definir una función objetivo y determinar los

mejores valores de las variables de diseño
Interacción entre Etapas de la Ingeniería de
Procesos
9La optimización requiere de la solución de
problemas de simulación en cada iteración
9La optimización es una herramienta
imprescindible en el diseño de un proceso

Simulación

Diseño Optimización
Introducción: Algunos Conceptos
en la Optimización de Procesos
Previo a la Optimización: Modelación
Representación Matemática de la Fisicoquímica
del proceso:

z Balances de Masa
z Balances de Energía Sistema de
z Relaciones Termodinámicas Ecuaciones
z Ecuaciones de Diseño No Lineales
z Balances de Momentum
z Restricciones Particulares
Análisis de Grados de Libertad

Para un sistema de M Ecuaciones y N Variables,

Variables el
número de grados de libertad, F, está dado por:
F=N-M
Tres casos:
F = 0 El sistema tiene solución UNICA
F > 1 El sistema puede OPTIMIZARSE
F < 0 El sistema está sobre especificado:
MODELO INCORRECTO
 Simulación u Optimización?
Grados de Libertad
F = Número de Variables – Número de Ecuaciones
F=N-M

Simulación ó Análisis Optimización

F=0 F>1
El sistema debe ser consistente Función Objetivo: Maximizar
utilidades, Minimizar costos, etc.

Función Objetivo:
Objetivo
obtención de diseños
óptimos
 Simulación u Optimización?
Simulación Optimización

x1 + x2 = 2 x1 + x2 = 2
x1 = 3 x2 x1 , x2 ≥ 0
x1 , x2 ≥ 0
F = 2 −1 = 1
F = 2−2 = 0 Soluciones posibles:
Solución única min x1 − x2 x1 x2
x1 = 1.5 0 2
Función
x2 = 0.5 objetivo 1 1
1.5 0.5
2 0
Solución óptima M M
 Selección de Variables de Diseño
M = 900 N = 1000

¿ Como seleccionar 100 Variables de Diseño ?

Matriz de Incidencia

x1 x2 x3 x4

f1 = ln( x1 ) − 2 = 0 f1 X
f 2 = x2 − 3 x4 − 5 = 0 f2 X X
f 3 = ( x2 ) − x3 + x4 − 1 = 0
3
f3 X X X
 Trayectorias de Steward

Variable de Diseño: Cualquiera de x2, x3 y x4

x1 x2 x3 x4

f1 X
f2 X X
f3 X X X
 Representación Matemática del
Problema de Optimización
• Variables Discretas y Continuas
• Restricciones (Ecuaciones, Desigualdades)
Lineales y No Lineales

min f ( x, y )
s.t. h( x, y ) = 0
g ( x, y ) ≤ 0
x ∈ R n , y ∈ {0,1}
 El Modelo Matemático
min f ( x, y )
s.t. h( x, y ) = 0
g ( x, y ) ≤ 0
x ∈ R n , y ∈ {0,1}

Límites:
Minimizar Costos 0< Comp <1
Maximizar Utilidades Temperatura,
Presión, etc.

Decisiones discretas:
Balances de Materia, Energía, ¿ Equipo Existe ?
Relaciones de Equilibrio, etc.
 z ¿ Restricción Lineal o No Lineal ?
2x +3 y=1

yx + 3 y = 1 x 2 + ln ( y ) = 1

z ¿ Variable Continua o Discreta ?

y
y=0 y=1
200 K < Temperatura < 500 K
T = 253.75 K Equipo Equipo
T = 493.68 K No Existe Existe
 Clasificación de Técnicas de Optimización

Tipos de Variables

Enteras Continuas Enteras+Continuas

Tipos de Restricciones Tipos de Restricciones

Programación
Entera
Programación Programación
Programación Programación Mixta Entera
Mixta Entera
Lineal No Lineal Lineal No Lineal
Tipos de Problemas de Optimización
(Determinística, Estado Estable)

* Programación Lineal (LP)

* Programación Mixta Entera Lineal (MILP)

* Programación No Lineal (NLP)

* Programación Mixta Entera No lineal (MINLP)

Región Factible y Convexidad

min f ( x)
s.t. h( x) = 0
g ( x) ≤ 0
x ∈ Rn
NLP
 Función Convexa o No Convexa
Función Convexa

f(x) es convexa si para toda x1 y x2 ε R:

f (αx1 + [1 − α ]x2 ) ≤ αf ( x1 ) + (1 − α ) f ( x2 ), ∀α ∈ (0,1)

f(x) f(x)

x1 x2 x1 x2

Convexa No
Convexa
 Función Convexa?
Función f ( x) = x 3 − 6 x 2 + 11x − 6 Intervalo (1,3)

f ( x1 ) = f (1) = 0 f ( x2 ) = f (3) = 0

f (αx1 + [1 − α ]x2 ) ≤ αf ( x1 ) + (1 − α ) f ( x2 )

f (2.5) ≤ 0.25 f (1) + (0.75) f (3)

α = 0.25 No
− 0.375 ≤ 0
Convexa

α = 0.75 f (1.5) ≤ 0.75 f (1) + (0.25) f (3)

0.375 ≤ 0
No se cumple
 Región Factible y Convexa

Región Factible. Definición:

{
FR = x h(x ) = 0, g (x ) ≤ 0, x ∈ R n }
Convexa o No Convexa?
Convexa si para toda x1 y x2 ε FR:

x = αx1 + (1 − α ) x2 ∈ FR, ∀α ∈ (0,1)

x1 x1
x2
x2

Convexa No
Convexa
Técnicas de Optimización
 Técnicas de Optimización

* Métodos Simplex y Puntos Interiores (LP)

* “Branch and Bound” (Ramificación y

Acotamiento) (MILP)

* Estrategia del Conjunto Activo, SQP

(Programación Cuadrática Sucesiva) (NLP)

* “Outer Approximation” y Descomposición de

Benders (MINLP)
Programación Lineal
 Programación Lineal

Forma General
T
min c x Maximize 300 x1 + 200 x2
s.t. A x ≤ b sujeto a 5 x1 + 2 x2 ≤ 180
x≥0 3x1 + 3 x2 ≤ 135
x ∈Rn
x1 ≤ 25

5 2 180
300  x1 
c=  x=  A = 3 3 b = 135
200  x2  1 0  
   25 
 Método Simplex
Si la función objetivo
disminuye en esta dirección

g1=0
g2=0

Este sería el punto óptimo Región

de
búsqueda

El óptimo se encuentra siempre

en un punto “esquina”
 Método Simplex
1) Introducir variables “slack” para convertir desigualdades en
igualdades. El número de grados de libertad no cambia.
2) Utilizar x=0 como punto inicial
3) Construir Matriz de Coeficientes
4) Efectuar pasos de eliminación Gaussiana hasta no obtener
valores negativos en el renglón correspondiente a la función
objetivo

5 x1 + 2 x2 ≤ 180 5 x1 + 2 x2 + s1 = 180
3x1 + 3 x2 ≤ 135 3x1 + 3 x2 + s2 = 135
x1 ≤ 25 x1 + s3 = 25
 Método Simplex
x1 x2 s1 s2 s3 b
Punto inicial
5 2 1 0 0 180 s1
x1 = 0
Matriz de x2=0
3 3 0 1 0 135 s2
Coeficientes s1 = 180
1 0 0 0 1 25 s3
s2=135
-300 -200 0 0 0 0 f
s3 = 25

Eliminación Gaussiana

s3 x2 s1 s2 s3 b

0 2 1 0 -5 55 s1

0 3 0 1 -3 60 s2

1 0 0 0 1 25 x1

0 -200 0 0 300 7500 f

 Pivoteo en el Método
Simplex

El pivote en los pasos de eliminación Gaussiana

se escoge de modo que:

• Columna: el coeficiente de la función

objetivo es el más negativo
• Renglón: la menor proporción entre b y el
coeficiente de la columna seleccionada
 Renglón Pivote
¿Qué ocurre si no se toma la menor proporción?

s2 x2 s1 s2 s3 b

0 -3 1 -5/3 0 -45 s1

1 1 0 1/3 0 45 x1

0 -1 0 -1/3 1 -20 s3

0 300 0 100 0 1350 f

 Método Simplex
Eliminación Gaussiana

s3 s2 s1 s2 s3 b

0 0 1 -2/3 -3 15 s1

0 1 0 1/3 -1 20 x2

1 0 0 0 1 25 x1

0 0 0 200/3 100 11500 f

 Notas en el Método Simplex
1. En cada punto esquina, un número de
variables igual al número de grados de
libertad del problema valen cero
2. Las variables cuyo valor es cero se
denominan no básicas
3. Cada punto de la trayectoria se dice que es
una solución básica
4. El algoritmo concluye cuando ya no hay
coeficientes negativos en el renglón de la
función objetivo
 Método Simplex

Maximize 300 x1 + 200 x2

sujeto a 5 x1 + 2 x2 ≤ 180 x2

3x1 + 3 x2 ≤ 135
x1 ≤ 25 (0,90)

(0,45) (25,27.5)

(25,20)
x1
(0,0)
El óptimo en el punto (25,20) (25,0)
Programación No Lineal
 Condiciones de Optimalidad
Optimización Sin Restricciones
min f ( x )
x ∈ Rn

necesaria x es un punto crítico de f(x)

Condición necesaria:
si se cumple que

 ∂f 
 ∂x1  Se obtiene un sistema
∇f (x ) =  M  = 0 de n Ecuaciones con n
∂f  Variables x
 ∂xn 
 
 Condiciones de Optimalidad
Optimización Sin Restricciones
min f ( x )
x ∈ Rn

 ∂f 
 ∂x1 
∇f (x ) =  M  = 0
∂f 
 ∂xn 
 
Condición suficiente:
suficiente La Matriz Hessiana de la
función objetivo es definida positiva

f (x + ∆ x ) = f (x ) + ∇f (x ) ∆ x + 1 ∆ x H ∆ x
T T
2

 ∂2 f ∂2 f 
 ∂x12 ∂x1 ∂x2  Hessiana para el
H = 2 
∂ f ∂2 f caso de 2 variables
 ∂x ∂ x 2 
 2 1 ∂x1 

1 ∆ xT H ∆ x ≥ 0
2
Optimización Sin Restricciones
min ( x1 ) 2 − 6 x1 + (x2 ) − 2 x2
2

 2 x1 − 6  x1 = 3
∇f ( x ) =   =0
 2 x2 − 2  x2 = 1

2 0 H es positiva definida
H = 
 0 2 
Optimo Global
 Condiciones de Optimalidad
Optimización Con m Restricciones de Igualdad
min f ( x )
s.t. h( x) = 0
x ∈ Rn

Función de Lagrange (función escalar):

m
L ( x, λ ) = f ( x ) + λ h ( x ) = f ( x ) + ∑ λ j h j ( x )
T

j
λ se denominan multiplicadores de
Lagrange y constituyen m variables
adicionales en el problema
 Condición necesaria:Obtener
necesaria: un punto crítico
(estacionario) para la función de Lagrange

∂L( x, λ )
= ∇ f ( x ) + ∑ λ j ∇h j ( x ) = 0 Se obtiene un sistema
∂x j
de n+m Ecuaciones
∂L( x, λ ) con n+m Variables x
= h( x ) = 0 yλ
∂λ

Condición suficiente: La Matriz Hessiana de la

función de Lagrange es definida positiva
 Optimización con
Restricciones de Igualdad
min ( x1 ) 2 − 6 x1 + (x2 )2 − 2 x2
s.t. x1 − x2 − 2 = 0

 2 x1 − 6  1
∇f ( x ) + ∑ λ j ∇h j ( x ) =   + λ1   = 0
j  2 x2 − 2  − 1
h( x ) = x1 − x2 − 2 = 0

2 x1 − 6 + λ1 = 0
2 x2 − 2 − λ1 = 0
x1 − x2 − 2 = 0
 Condiciones de Optimalidad
Optimización Con m Restricciones de Igualdad y r de
Desigualdad
min f ( x)
s.t. h( x) = 0
g ( x) ≤ 0
x ∈ Rn
Función de Lagrange Aumentada (función escalar):
m r
L ( x, λ , µ ) = f ( x ) + λ h ( x ) + µ g ( x ) = f ( x ) + ∑ λ j h j ( x ) + ∑ µ k g k ( x )
T T

j =1 k =1

µ se denominan multiplicadores de Karush-

Kuhn-Tucker y constituyen r variables
adicionales en el problema
Condición necesaria:Obtener
necesaria: un punto crítico
(estacionario) mediante el Teorema de Karush-
Kuhn-Tucker
∂L( x, λ , µ )
= ∇f ( x ) + ∑ λ j ∇ h j ( x ) + ∑ µ k ∇g k ( x ) = 0
∂x j k

∂L(x, λ , µ )
= h( x ) = 0 Se obtiene un sistema
∂λ de n+m+r Ecuaciones
µk (x )⋅ g k (x ) = 0 con n+m+r Variables
x, µ y λ
µk (x ) ≥ 0 g k (x ) ≤ 0

Condición suficiente: La Matriz Hessiana de la

función de Lagrange Aumentada es definida
positiva
 Condición necesaria:

∂L( x, λ , µ )
= ∇f ( x ) + ∑ λ j ∇ h j ( x ) + ∑ µ k ∇g k ( x ) = 0
∂x j k

∂L(x, λ , µ ) µk (x )⋅ g k (x ) = 0
= h( x ) = 0
∂λ
µk (x ) ≥ 0 g k (x ) ≤ 0
 Optimización con
Restricciones de
Desigualdad
min f ( x ) = 1 (x12 + x22 ) − 3 x1 − x2
2
s.t. g1 = − x1 + x2 ≤ 0
g 2 = x1 − 1 x2 − 2 ≤ 0
2 Note:

 x1 − 3 − 1  1 
∇f ( x ) + ∑ µ k ∇hk ( x) =   + µ1   + µ 2   =0
k  x2 − 1 1  − 1 2

x1 − µ1 + µ 2 = 3
x2 + µ1 − 1 µ 2 = 1
2
µ1 (− x1 + x2 ) = 0
( )
µ 2 x1 − 1 x2 − 2 = 0
2
 Programación No Lineal
Estrategia del Conjunto Activo (Active Set Strategy)
Solución al Conjunto de Ecuaciones KKT
∇f ( x ) + ∑ λ j ∇h j ( x ) + ∑ µ k ∇ g k ( x ) = 0
j k

h( x ) = 0
µ k (x ) ⋅ g k (x ) = 0

µk (x ) ≥ 0 g k (x ) ≤ 0

Desigualdades
g k (x ) = 0 g k (x ) < 0

Activa Inactiva
 Programación No Lineal
Estrategia del Conjunto Activo

1) Definir conjunto activo. Inicialmente:

J1 = {k g k = 0}

2) Formular ecuaciones KKT

J1 = ∅ ⇒ g k > 0 µk = 0

∇f ( x ) + ∑ λ j ∇h j ( x ) + ∑ µ k ∇g k ( x ) = 0
j k∈J1

h( x ) = 0

g k ( x ) = 0 k ∈ J1
 Programación No Lineal
Estrategia del Conjunto Activo

3) Si para toda k g k ( x) ≤ 0 y µ k ≥ 0 OK
OK. Se ha obtenido la solución

Si cualquier g k ( x) > 0 y / o µ k < 0

a) Eliminar de J1 la restricción con el µi más negativo

b) Añadir a J1 todas las restricciones violadas g k ( x) > 0

para hacerlas activas

c) Regresar a 2)
 Programación No Lineal
Estrategia del Conjunto Activo: Ejemplo
min f ( x ) = 1 (x12 + x22 ) − 3 x1 − x2
2
s.t. g1 = − x1 + x2 ≤ 0
g 2 = x1 − 1 x2 − 2 ≤ 0
2
g 3 = − x2 ≤ 0
Condiciones de Karush-Kuhn-Tucker (Iteración 1):
J1 = {k g k = 0} ∇f ( x ) + ∑ λ j ∇h j ( x ) + ∑ µ k ∇g k ( x ) = 0
j k∈J1

J1 = ∅ ⇒ g k > 0 µk = 0

Verificando restricciones y µ´s:

 x − 3
∇f ( x ) + ∑ µ k ∇hk ( x) =  1  = 0 g1 = −2 < 0 Activar g2
k  x2 − 1 g 2 = 1 / 2 > 0 !!!
x1 = 3
g 3 = −1 < 0 J 2 = {2}
x2 = 1
µ1 = 0
µ1 = 0
µ2 = 0 µ2 = 0
µ3 = 0 µ3 = 0
 Programación No Lineal
Programación Cuadrática Sucesiva (SQP)
min f ( x )
s.t. h( x ) = 0
g ( x) ≤ 0 Condiciones de optimalidad
x ∈ Rn (Karush-Khun-Tucker)

Sistema de Ecuaciones
No Lineales

Iteración de Newton

Programa cuádratico
Programación Mixta-Entera
 Representación de Procesos en
Términos de Variables Binarias
x1 z1
I

x0
A 10 kmol/hr B

x2 z2
II

 1 if reactor I is selected
y1 =  min C = 7.5 y1 + 6.4 x1 + 5.5 y2 + 6.0 x2
0 if reactor I is not selected
sujeto a 0.8 x1 + 0.67 x2 = 10
 1 if reactor II is selected x1 − 20 y1 ≤ 0 x2 − 20 y2 ≤ 0
y2 = 
0 if reactor II is not selected x1 , x2 ≥ 0 y1 , y2 = 0,1
 Relaciones Lógicas
¬p j
1) NOT
1− y j 7) Teorema de Morgan
pi ∨ p j
2) OR (exclusivo) y + y = 1 ¬( A ∨ B ) ⇔ ¬A ∧ ¬B
¬( A ∧ B ) ⇔ ¬A ∨ ¬B
i j

pi ∨ p j
3) OR (inclusivo)
yi + y j ≥ 1 8) Distribución de
pi ∧ p j AND
4) AND y 1, y 1
j ≥ i ≥
( A ∧ B ) ∨ C ⇔ ( A ∨ C ) ∧ (B ∨ C )
pi → p j ¬pi ∨ p j
5) If-Then
1 − yi + y j ≥ 1 yi ≤ y j
pi ↔ p j
6) Iff-Then
yi = y j
 Representando Alternativas
Big - M Convex Hull
x11 + x12 = x1
x1 − 2 x2 ≤ M (1 − y1 )
x1 − 2 x2 ≥ − M (1 − y1 ) x21 + x22 = x2
x1 − 1 ≤ M (1 − y1 ) x11 ≤ My1
y1 + y2 = 1 x12 ≤ My2
x21 ≤ My1
x1 − 5 x2 ≤ M (1 − y2 )
x22 ≤ My2
x1 − 5 x2 ≥ − M (1 − y2 )
x1 − 1 ≥ − M (1 − y2 ) x11 − 2 x21 = 0
x1 , x2 ≥ 0     x12 − 5 x22 = 0
 y1   y2  x11 ≤ y1
   
= ∨ =
 x1 2 x2   x1 5 x2  x12 ≥ y2
 x1 ≤ 1   x1 ≥ 1  y1 + y2 = 1
 x , x ≥ 0  x , x ≥ 0
 1 2   1 2  x11 , x12 , x21 , x22 ≥ 0
 Programación Mixta-
Entera Lineal
Método de “Branch
“ and Bound” (Ramificación y Acotamiento)

Procedimiento:

1) Resolver el problema como si todas las variables fueran

continuas (problema LP). Es decir, relajar las variables binarias.
A este nivel se le conoce como “Root Node” (Nodo raíz)
• El valor de la función objetivo en el nodo raíz es un límite
inferior a su valor en la solución
• Si la solución del nodo raíz es entera, tal solución es la
solución del problema

2) Realizar una búsqueda de ramificación ordenada mediante la

adición de restricciones enteras
Método de “Branch
“ and Bound” (Ramificación y Acotamiento)

y3=0

0<y2<1 y2=0 y3=1

0<y3<1

y1=0
y2=1
0<y1<1
0<y2<1
0<y3<1 y1=1
Relajado

Para m variables binarias, el número de posibles nodos es 2m+1- 1

3 Reglas en el Método de “Branch
“ and Bound”:

1) La función objetivo en i es un límite inferior a la función

objetivo en k.
2) Si algún nodo resulta en una solución entera, el valor de la
función objetivo en tal nodo es un límite superior de la solución
3) Si i es infactible o ilimitado, entonces k también lo es
Método de “Branch
“ and Bound”: Ejemplo
[1,1,0]
min z = x + y1 + 3 y2 + 2 y3 y2=1
z =9
[1,0.5,0]
s.t. − x + 3 y1 + 2 y2 + y3 ≤ 0 z = 6.5

− 5 y1 − 8 y2 − 3 y3 ≤ −9 infactible
x ≥ 0, y1 , y2 , y3 ∈ {0,1}
y1=1 y2=0
[0,1,1]
y2=1 z =8
[0.2,1,0]
z = 5.8

y3=1
y1=0 [0,0.075,1]
z = 6.75 y =0
m=3 [0,1,0.333] 2

z=6
2 m +1 − 1 = 15 y3=0
infactible

infactible
 Programación Mixta-
Entera No Lineal
Seleccionar un valor inicial para y

Resolver el NLP resultante S(yk)

ZU
Nueva y
Algoritmo Resolver el Problema Maestro MILP
General ZL
SI

ZL<ZU ?

NO

Solución
 Programación Mixta-
Entera No Lineal
Algoritmo “Outer
“ Approximation” (DICOPT+++)

Problema Maestro
Z = min α
T
min c y + f ( x) ( ) ( )( T
sujeto a α ≥ cT y + f x k + ∇f x k x − x k )
g ( x) + B y ≤ 0 ( ) ( )( )T
g x k + ∇g x k x − x k + By ≤ 0  ∀k ∈ T
Ay ≤ a Ay ≤ a 
y ∈ {0,1} x ∈ R n
m
y ∈ {0,1} x ∈ R n
m 

α ∈R 

f(x)
k x k es la solución óptima de S y k ∀
T =
( )
k 
 los posibles valores de y 

x
 Algoritmo “Outer
“ Approximation”

Optimización Semi-Infinita: Problema Relajado

Z = min α
( ) ( )( T
sujeto a α ≥ cT y + f x k + ∇f x k x − x k )
( ) ( )( )T
g x k + ∇g x k x − x k + By ≤ 0  k = 1K K
Ay ≤ a 
y ∈ {0,1} x ∈ R n
m 

α ∈R 

Dada la solución de K subproblemas NLP x k ( )

definidos por y k tal que k = 1K K
 Programación Mixta-
Entera No Lineal
“Outer Approximation”: Ejemplo 1
min z = y1 + 1.5 y2 + 0.5 y3 + x12 + x22
s.t. (x1 − 2 ) − x2 ≤ 0 x1 − (1 − y1 ) ≥ 0
2
y1 + y2 + y3 ≥ 1
x1 − 2 y1 ≥ 0 x2 − y 2 ≥ 0 0 ≤ x1 ≤ 4 y1 , y2 ∈ {0,1}
x1 − x2 − 4(1 − y2 ) ≤ 0 x1 + x2 ≥ 3 y3 0 ≤ x2 ≤ 4

1) NLP y1 = 1 y2 = 1 y3 = 1 x1 = 2 x2 = 2 zU = 11

MILP y1 = 1 y2 = 0 y3 = 0 x1 = 2 x2 = 0 z L = 1

2) NLP y1 = 1 y2 = 0 y3 = 0 x1 = 2 x2 = 0 zU = 5

MILP y1 = 0 y2 = 1 y3 = 0 x1 = 1 x2 = 0 z L = 1.5

3) NLP y1 = 0 y 2 = 1 y3 = 0 x1 = 1 x2 = 1 zU = 3.5

MILP y1 = 0 y2 = 0 y3 = 1 x1 = 2 x2 = 1 z L = 4.5
 Programación Mixta-
Entera No Lineal
“Outer Approximation”: Ejemplo 2, Iteración 1
min z = −2.7 y + x 2
s.t. g1 = − ln (1 + x ) + y ≤ 0
g 2 = − ln ( x − 0.57 ) + y − 1.1 ≤ 0
0≤ x≤2 y ∈ {0,1}
1) Comenzar con y=1 y resolver NLP:
z = 0.2525 x = 1.7183 µ1 = 9.347 µ 2 = 0
( zU )
2) Linealizar el problema MINLP en x= 1.7183 para obtener el problema maestro MILP:

zOA = min α OA
z L < zU
s.t. α OA ≥ −2.7 y + 3.4366 x − 2.9525
y=0
− 0.36787 x + y ≤ 0.36788 Volver a paso 1
zOA = −1.939 con nuevo valor
− 0.87085 x + y ≤ −0.2581
( zL ) y=0
0≤ x y ∈ {0,1}
El Entorno de Modelación
GAMS y sus Resolvedores
Formas de Atacar el Problema de
Modelación
z Corriente Modular Secuencial

* Cada unidad de proceso (módulo) calcula su salida

dadas las entradas
⇓
En cada unidad se resuelve un sistema de ecuaciones no
lineales (subrutina que “convierte” valores de entradas
en valores de salida)
* Las variables de las corrientes de reciclo se suponen
inicialmente y constituyen la guía para proceso iterativo
 Compresor
Divisor de Corriente

Corriente
de
Reciclo
Intercambiador

Separador
Flash
Alimentación
Mezclador Reactor

* Se requieren algoritmos de orden de precedencia y

determinación de las corrientes de reciclo a ser supuestas.
z Corriente de Orientación a las Ecuaciones

* Cada unidad del sistema es representado por un

sistema de ecuaciones (lenguaje declarativo:
declarativo no
restricciones en cuanto a especificación de variables)
* Las ecuaciones son generalmente almacenadas en
librerías
* Para definir el problema, se colectan las ecuaciones
que representan a cada unidad y a las conexiones
entre ellas
Sistema
Resultante

* El sistema resultante se resuelve simultáneamente

(probablemente un problema de miles de ecuaciones)
 Simuladores y Optimizadores
Comerciales
9 Corriente Modular Secuencial: Simuladores de
Procesos ASPEN,
ASPEN PRO/II, HySys
¾ Modelado por configuración

9 Corriente de Orientación a las Ecuaciones: Sistemas

de Modelación gPROMS, SpeedUp, ABACO,
ASCEND
9 GAMS
¾ NLP CONOPT, MINOS, LACELOT, SRQP, LINGO
NLP:
¾ MINLP DICOPT
MINLP:
¾ LP CPLEX, OSL, LINDO, SCICONIC, XA
MILP, LP:
GAMS la Interfase
 Resolvedores en GAMS
Tipos de Problemas

Algoritmos
de Solución
 GAMS: Ejemplo 1
W1
y0 = 0

Q=1000 lb/hr
xF = 0.2
Etapa de Q=1000 lb/hr
Extracción x1

W1
y1

Max Q(xF − x1 ) − λW1

λ = 0.05
Qx F = Qx1 + Wy1
H =1.2
Hx1
y1 =
(H − 1)x1 + 1
 GAMS: Ejemplo 2
zA = 0.6 xDA = 0.8
αAC = 2.3
zB = 0.3
zC = 0.1 αBC =1.3
F = 1000 Kmol/hr
αCC =1.0

q =1.0

N α zj
∑j =1 α
j

−θ
= 1− q
j

N α j x Dj
∑α
j =1 −θ
= 1 + Rmin
j

∑x j =1
D
j =1
Algunas Aplicaciones en
Ingeniería Química
 Optimización de Columnas y
Secuencias de Destilación

Diseño Óptimo de ¿ Cómo separar una mezcla ABC ?

Columnas
A
1)

ABC B
D, Pureza

Alimentación
A
2)
BC
AB
C

Minimizar Costo ABC

B
N=? R=? P=?
C
 Modelo MILP para
Secuencias de Separación
9Supone separaciones del 100% (Sharp splits)
splits
9Supone que las cargas térmicas y costos de inversión son
funciones lineales de la alimentación
9Supone que las cargas de los cambiadores son del mismo
orden de magnitud l∈C top
k ∑x l
F

l∈C ktop
ξ ktop =
∑x
l∈C k
l
F

x iF
k

l ∈ Ck
∑x l
F

l∈C kbot
ξ kbot =
∑x l
F

l ∈ C kbot l∈C k
 Datos de Caso de Estudio
Mezcla Cuaternaria ABCD
F=1000 Kmol/hr: 15% A, 30% B, 35% C, 20% D

Sistema Fijo Variable Carga térmica

αk βk Kk
3 3 6
(10 $/año) (10 $ hr / Kgmol año) (10 KJ / Kgmol)
A/BCD 145 0.42 0.028
AB/CD 52 0.12 0.042

ABC/D 76 0.25 0.054

A/BC 125 0.78 0.024
AB/C 44 0.11 0.039
B/CD 38 0.14 0.040
BC/D 66 0.21 0.047
A/B 112 0.39 0.022
B/C 37 0.08 0.036
C/D 58 0.19 0.044

Costos de utilidades:
Agua de enfriamiento: CW = 1.3 (103$ hr/ 106 KJ año)
Vapor: CH = 34 (103$ hr/ 106 KJ año)
 Superestructura
B/CD
F4

y4 F8
C/D
A/BCD y8

F5

y5
BC/D
F1
y1

Superestructura con ABCD

F2
AB/CD y9
4 Componentes FTOT y2
F9
B/C

F3

y3
A/BC
F6
y6
ABC/D y10
A/B
y7 F10
F7
AB/C
 Modelo Balance Global
Factores de separación FTOT = 1000 = F1 + F2 + F3
A
ξ1A =0.15 ξ6 =0.188

ξ1BCD =0.85 ξ6BC =0.812 Mezclas Intermedias

AB AB
ξ2 =0.45 ξ7 =0.5625
CD C
ξ2 =0.55 ξ7 =0.437
BCD F4 + F5 − 0.85F1 = 0
ξ3ABC =0.8 ξ8C =0.636
ξ3D =0.2 ξ8D =0.364
ξ4B =0.353 ξ9B =0.462 ABC F6 + F7 − 0.8F3 = 0
ξ4CD =0.647 ξ9C =0.538
ξ5BC =0.765 ξ10A =0.333
AB F10 − 0.45 F2 − 0.563F7 = 0
ξ5D =0.235 ξ10B =0.667

BC F9 − 0.765F5 − 0.812 F6 = 0

CD F8 − 0.55 F2 − 0.647 F4 = 0
Flujos (Big-M)

Fk − 1000 y k ≤ 0 Fk ≥ 0 y k = 0,1 k = 1,...,10

Cargas Térmicas Función Objetivo

Qk = K k Fk k = 1,...,10 10 10
min C = ∑ (α k y k + β k Fk ) + (34 + 1.3)∑ Qk
k =1 k =1
 Síntesis de Redes de
Intercambio de Calor
h1
500 ° K

h2 ?
418 ° K

c1
305 ° K
375 ° K

c2

310 ° K 298 ° K
 Modelo MINLP
(SYNHEAT)
Estrategia de Optimización Simultánea

Stage 1 Stage 2

H1-C1 H1-C1

H1
H1-C2 H1-C2

H2-C1 H2-C1 C1

H2
H2-C2 H2-C2

C2
 MODELO
Balance de calor total para cada corriente:
Carga de servicios de enfriamiento y
(TIN i − TOUT i )F i = ∑ ∑q ijk + qcu i , i ∈ HP calentamiento:
k ∈ ST j ∈ CP

(TOUT j − TIN j )F j = ∑ ∑q ijk + qhu j , j ∈ CP (t − TOUT )F = qcu , i ∈ HP

i , NOK +1 i i i
k ∈ ST i∈ HP
(TOUT − t )F = qhu , j ∈ CP
j j ,1 j j
Balance de calor para cada etapa:

(t ik − t i , k + 1 )F i = ∑q ijk , k ∈ ST , i ∈ HP Restricciones lógicas:

j ∈ CP

(t jk − t j , k + 1 )F j = ∑q ijk , k ∈ ST , j ∈ CP qijk − Ωzijk ≤ 0, i ∈ HP, j ∈ CP, k ∈ ST

i∈ HP
qcui − Ωzcui ≤ 0, i ∈ HP
Asignación de las temperaturas de entrada a la
superestructura: qhuj − Ωzhuj ≤ 0, j ∈ CP
zijk , zcui , zhuj = 0,1
t i ,1 = TIN i , i ∈ HP
t j , NOK +1 = TIN j j ∈ CP Cálculo de las diferencias de
temperaturas:
Factibilidad de temperaturas:
t ik ≥ t i , k +1 , k ∈ ST , i ∈ HP dtijk ≤ tik − t jk + Γ(1 − zijk ), k ∈ ST, i ∈ HP, j ∈ CP
t jk ≥ t j , k +1 , k ∈ ST , j ∈ CP dtij,k +1 ≤ ti,k +1 − t j ,k +1 + Γ(1 − zijk ), k ∈ ST, i ∈ HP, j ∈ CP
TOUT i ≤ t i , NOK +1 , i ∈ HP dtcui ≤ ti, NOK+1 − TOUTcu + Γ(1 − zcui ), i ∈ HP
TOUT j ≥ t j ,1 , j ∈ CP dthuj ≤ TOUThu − t j ,1 + Γ(1 − zhu j ), j ∈ CP
 Planeación de la
Producción
La Producción del Petróleo es un Problema
Multiperiódico

Thousands of barrels/day
3200
d1 3100
Demand d3 d4
d2 3000
2900
2800
2700
2600
T1 T2 T3 T4 2500
Jan Mar May Jul Sep Nov Jan
H
Month
 well head
well bore storage Geological Properties:
oil permeability
flow thickness
porosity
well bore etc.
reservoir
oil flow

La extracción produce una disminución de la presión del

pozo con el tiempo
pw q
pup
pflow
pshut
plow

ts t ts t
 Modelos de
Programación Mixta-
Entera

El tiempo de operación H es dividido en NP periodos

de tiempo. Dadas las demandas de producción del
petróleo en cada periodo de tiempo y las constantes de
caracterización de los pozos, determinar:
• Los perfiles de producción y
• Los tiempos operación de cada pozo

en cada periodo de tiempo.

 Modelo MILP
Minimize ∑∑γ (
q T + ∑∑δij yijT + ∑∑αij 1 − yij T
ij ij )
i j i j i j

∑q T ≥ d
i
ij j ∀j ∈ P
 ¬Yij 
   
 Yij   Wij1 Wij 2 
 p f = pin − D  ∨  p f = pin + I  ∨  p f = pup  ∀i ∈W , j ∈ P
 ij ij ij   ij ij ij   ij i 
 in up    
p
 ij + I ij ≤ pi  p
 ij
in
+ I ij > p up
i  
Dij = qij {c1 [ln(T ) + c2 ]} ∀i ∈W , j ∈ P
I ij = qis {c1 [ln(T ) + c2 ]}(1 − yij ) ∀i ∈W , j ∈ P

(
qijmax {c1 [ln(T ) + c2 ]} = pijin − pilow )
∀i ∈W , j ∈ P
qij ≤ qijmax ∀i ∈W , j ∈ P

(
qij ≤ qiup yij + qlow 1 − yij )
∀i ∈W , j ∈ P
qij ≥ qlow ∀i ∈W , j ∈ P

pijin = pijf−1 ∀i ∈W , j ∈P
 Calendarización de
Procesos por Lotes
Gráficas de Gant

Producto A

8 8 8 8
Mezclador
20 20
Reactor 1
20 20
Reactor 2
4 4 4 4
Centrífuga

¾ Determinar el Tamaño y número de unidades paralelas

¾ Determinar la Secuenciación de las unidades
GAMS Instituto Tecnológico de Celaya

Código GAMS
El código de GAMS se puede escribir con cualquier procesador de texto o a
través de la interfase de GAMS. Si se utilizan procesadores especializados como
Word, FrameMaker, PageMaker, etc., asegúrese de guardar el archivo sin
formato (como texto, código ASCII).

Los archivos de GAMS deberán tener la extensión *.gms

Luego de la solución de algún modelo, GAMS crea un archivo de resultados

también en formato de texto y con el mismo nombre que el archivo del código,

pero con extensión *.lst

Como regla general, un modelo de GAMS debe contener las siguientes partes
(se muestra un caso ilustrativo):
1) Título
$TITLE MULTIPRODUCTO

2) Declaración de Conjuntos
SETS
J COMPONENTES /1*3/

3) Declaración de Parámetros
PARAMETERS SA, SB, SC;

4) Declaración de Variables (positivas y generales)

VARIABLES P;
POSITIVE VARIABLES X7, X8, X9, X10, X11, X12;

5) Declaración de Ecuaciones
EQUATIONS RES1, RES2, RES3, INE1, INE2, INE3,OBJ;

6) Ecuaciones del Sistema

RES1.. X11 =E= 0.667*X8 + 0.667 *X9 + 0.5*X10;

Departamento de Ingeniería Química 1

GAMS Instituto Tecnológico de Celaya

Note que el identificador de la ecuación va precedido de dos puntos. En las

ecuaciones el símbolo =E= significa igual, =G= significa mayor que y =L=
significa menor que.

7) Definición de una función objetivo (“Dummy” o verdadera)

OBJ.. P =E= 0.025*X8 + 0.028*X9 + 0.028*X10 - 0.015*X11 -
0.02*X12 - 0.025*X7;

8) Establecimiento de las ecuaciones que componen un modelo en

particular
MODEL PLANTAS /ALL/;

9) Valores de parámetros, estimados iniciales, límites de las variables

SA = 40000;
TETA.L('1')= 1.05;

10) Llamado a la técnica de solución de acuerdo al tipo de problema

SOLVE PLANTAS USING MIP MAXIMIZING P;

Departamento de Ingeniería Química 2

GAMS Instituto Tecnológico de Celaya

Ejemplos Ilustrativos del Uso de GAMS

1. Resolver el problema de programación lineal de la planta multiproducto que
se desarrolló en clase. Recodar que las ecuaciones son:

x11 = 0.667 x8 + 0.667 x9 + 0.5 x10

x12 = 0.333 x8 + 0.333 x9 + 0.167 x10
x7 = 0.333 x10
x11 ≤ 40000
x12 ≤ 30000
x7 ≤ 25000

Mientras que la función objetivo está dada por:

P = 0.025 x8 + 0.028 x9 + 0.028 x10 − 0.015 x11 − 0.02 x12 − 0.025 x7

Departamento de Ingeniería Química 3

GAMS Instituto Tecnológico de Celaya

Código GAMS del Ejemplo 1

$TITLE MULTIPRODUCTO
*
*DEFINICION DE VARIABLES, PARAMETROS Y ECUACIONES
*

VARIABLES P;

POSITIVE VARIABLES X7, X8, X9, X10, X11, X12;

*
* DATOS CONOCIDOS
*
PARAMETERS SA, SB, SC;
*
*ECUACIONES
*
EQUATIONS RES1, RES2, RES3, INE1, INE2, INE3,OBJ;

*
* DEFINICION DE LAS ECUACIONES QUE FORMAN PARTE DEL MODELO
*
RES1.. X11 =E= 0.667*X8 + 0.667 *X9 + 0.5*X10;
RES2.. X12 =E= 0.333*X8 + 0.333*X9 + 0.167 *X10;
RES3.. X7 =E= 0.333*X10;
INE1.. X11 =L= SA;
INE2.. X12 =L= SB;
INE3.. X7 =L= SC;
OBJ.. P =E= 0.025*X8 + 0.028*X9 + 0.028*X10 - 0.015*X11 -
0.02*X12 - 0.025*X7;

MODEL PLANTAS /ALL/;

*
* ASIGNACION DE VALORES A LOS PARAMETROS
*
SA = 40000;
SB = 30000;
SC = 25000;

OPTION LIMROW=0;
OPTION LIMCOL=0;
*
* LLAMADO A LA TECNICA DE SOLUCION
*
SOLVE PLANTAS USING MIP MAXIMIZING P;

Departamento de Ingeniería Química 4

GAMS Instituto Tecnológico de Celaya

Resultados GAMS del Ejemplo 1

4 *DEFINICION DE VARIABLES, PARAMETROS Y ECUACIONES
5 *
6
7 VARIABLES P;
8
9 POSITIVE VARIABLES X7, X8, X9, X10, X11, X12;
10 *
11 * DATOS CONOCIDOS
12 *
13 PARAMETERS SA, SB, SC;
14
15 *
16 *ECUACIONES
17 *
18 EQUATIONS RES1, RES2, RES3, INE1, INE2, INE3,OBJ;
19
20
21 *
22 * DEFINICION DE LAS ECUACIONES QUE FORMAN PARTE DEL MODELO
23 *
24 RES1.. X11 =E= 0.667*X8 + 0.667 *X9 + 0.5*X10;
25 RES2.. X12 =E= 0.333*X8 + 0.333*X9 + 0.167 *X10;
26 RES3.. X7 =E= 0.333*X10;
27 INE1.. X11 =L= SA;
28 INE2.. X12 =L= SB;
29 INE3.. X7 =L= SC;
30 OBJ.. P =E= 0.025*X8 + 0.028*X9 + 0.028*X10 - 0.015*X11 -
0.02*X12 - 0.025*X7;
31
32 MODEL PLANTAS /ALL/;
33
34 *
35 * ASIGNACION DE VALORES A LOS PARAMETROS
36 *
37 SA = 40000;
38 SB = 30000;
39 SC = 25000;
40
41 OPTION LIMROW=0;
42 OPTION LIMCOL=0;
43
44 *
45 * LLAMADO A LA TECNICA DE SOLUCION
46 *
47 SOLVE PLANTAS USING MIP MAXIMIZING P;

MODEL STATISTICS

BLOCKS OF EQUATIONS 7 SINGLE EQUATIONS 7

Departamento de Ingeniería Química 5

GAMS Instituto Tecnológico de Celaya

BLOCKS OF VARIABLES 7 SINGLE VARIABLES 7

NON ZERO ELEMENTS 20

GENERATION TIME = 0.000 SECONDS 1.4 Mb WIN200-121

EXECUTION TIME = 0.000 SECONDS 1.4 Mb WIN200-121

S O L V E S U M M A R Y

MODEL PLANTAS OBJECTIVE P

TYPE MIP DIRECTION MAXIMIZE
SOLVER OSL2 FROM LINE 47

**** SOLVER STATUS 1 NORMAL COMPLETION

**** MODEL STATUS 1 OPTIMAL
**** OBJECTIVE VALUE 705.1354

RESOURCE USAGE, LIMIT 0.070 1000.000

ITERATION COUNT, LIMIT 3 10000

OSL Version 2 Mar 21, 2001 WIN.O2.O2 20.0 [Link] (Jan )

Work space allocated -- 0.09 Mb

LOWER LEVEL UPPER MARGINAL

---- EQU RES1 . . . -0.032

---- EQU RES2 . . . -0.020
---- EQU RES3 . . . -0.026
---- EQU INE1 -INF 40000.000 40000.000 0.017
---- EQU INE2 -INF 13766.923 30000.000 .
---- EQU INE3 -INF 25000.000 25000.000 0.001
---- EQU OBJ . . . 1.000

LOWER LEVEL UPPER MARGINAL

---- VAR P -INF 705.135 +INF .

---- VAR X7 . 25000.000 +INF .
---- VAR X8 . . +INF -0.003
---- VAR X9 . 3691.848 +INF .
---- VAR X10 . 75075.075 +INF .
---- VAR X11 . 40000.000 +INF .
---- VAR X12 . 13766.923 +INF .

**** REPORT SUMMARY : 0 NONOPT

0 INFEASIBLE
0 UNBOUNDED

Departamento de Ingeniería Química 6

GAMS Instituto Tecnológico de Celaya

2. Para el sistema de extracción mostrado en la Figura, utilice el sistema de

modelación GAMS para determinar los valores de las variables W1 y x1 que
maximizan la función:

Q( xF − x1 ) − λW1

donde λ = 0.05. Considere que la relación de equilibrio entre y1 y x1 está

dada por la expresión:

Hx1
y1 =
(H − 1)x1 + 1
Use un valor de H = 1.2. Note también que el balance de masa en el
sistema resulta en la ecuación:

Qx F = Qx1 + Wy1

W1
y0 = 0

Q=1000 lb/hr
xF = 0.2
Etapa de Q=1000 lb/hr
Extracción x1

W1
y1

Figura

Departamento de Ingeniería Química 7

GAMS Instituto Tecnológico de Celaya

Código GAMS del Ejemplo 2

$TITLE EXTRACCION
$OFFSYMXREF
$OFFSYMLIST
*
*DEFINICION DE VARIABLES, PARAMETROS Y ECUACIONES
*
VARIABLES F;

POSITIVE VARIABLES X1, Y1, W1;

PARAMETERS Q, XF, LAMBDA, H;

EQUATIONS MASBAL, EQUILIBRIO, OBJ;

*
*ECUACIONES
*
MASBAL.. Q * XF =E= Q * X1 + W1 * Y1;
EQUILIBRIO.. Y1 =E= (H * X1)/(((H - 1.0) * X1) + 1.0);
OBJ.. F =E= Q * ( XF -X1) - LAMBDA * W1;

*
* DEFINICION DE LAS ECUACIONES QUE FORMAN PARTE DEL MODELO
*
MODEL EXTRACTOR /ALL/;
*
* ASIGNACION DE VALORES A LOS PARAMETROS
*

Q = 1000;
XF = 0.2;
LAMBDA = 0.05;
H = 1.2;
*
* LIMITES Y VALORES INICIALES
*

Y1.L = 0.1;
[Link] = 1.0;
X1.L = 0.1;
[Link] = 0.2;
W1.L = 500;
OPTION LIMROW=0;
OPTION LIMCOL=0;
*
* LLAMADO A LA TECNICA DE SOLUCION
*
SOLVE EXTRACTOR USING NLP MAXIMIZING F;

Departamento de Ingeniería Química 8

GAMS Instituto Tecnológico de Celaya

Resultados GAMS del Ejemplo 2

COMPILATION TIME = 0.000 SECONDS 0.7 Mb WIN194-
116
Model Statistics SOLVE EXTRACTOR USING NLP FROM LINE 55

MODEL STATISTICS

BLOCKS OF EQUATIONS 3 SINGLE EQUATIONS 3

BLOCKS OF VARIABLES 4 SINGLE VARIABLES 4
NON ZERO ELEMENTS 8 NON LINEAR N-Z 3
DERIVATIVE POOL 5 CONSTANT POOL 10
CODE LENGTH 40

GENERATION TIME = 0.110 SECONDS 1.9 Mb WIN194-

116

EXECUTION TIME = 0.110 SECONDS 1.9 Mb WIN194-

116

S O L V E S U M M A R Y

MODEL EXTRACTOR OBJECTIVE F

TYPE NLP DIRECTION MAXIMIZE
SOLVER CONOPT FROM LINE 55

**** SOLVER STATUS 1 NORMAL COMPLETION

**** MODEL STATUS 2 LOCALLY OPTIMAL
**** OBJECTIVE VALUE 58.1881

RESOURCE USAGE, LIMIT 0.391 1000.000

ITERATION COUNT, LIMIT 15 10000
EVALUATION ERRORS 0 0

C O N O P T Wintel version 2.043C-005-039

Copyright (C) ARKI Consulting and Development A/S
Bagsvaerdvej 246 A
DK-2880 Bagsvaerd, Denmark

Using default control program.

** Optimal solution. Reduced gradient less than tolerance.

Departamento de Ingeniería Química 9

GAMS Instituto Tecnológico de Celaya

CONOPT time Total 0.219 seconds

of which: Function evaluations 0.051 = 23.2%
Derivative evaluations 0.000 = 0.0%

Work length = 0.05 Mbytes

Estimate = 0.05 Mbytes
Max used = 0.04 Mbytes

LOWER LEVEL UPPER MARGINAL

---- EQU MASBAL -200.000 -200.000 -200.000 0.463

---- EQU EQUILIBRIO . . . 464.178
---- EQU OBJ 200.000 200.000 200.000 1.000

LOWER LEVEL UPPER MARGINAL

---- VAR F -INF 58.188 +INF .

---- VAR X1 . 0.092 0.200 -2.233E-6
---- VAR Y1 . 0.108 1.000 .
---- VAR W1 . 1002.840 +INF .

**** REPORT SUMMARY : 0 NONOPT

0 INFEASIBLE
0 UNBOUNDED
0 ERRORS

EXECUTION TIME = 0.060 SECONDS 0.7 Mb WIN194-

116

Departamento de Ingeniería Química 10

 GAMS Instituto Tecnológico de Celaya

3. Considere la separación de la mezcla ternaria que se muestra en la figura.

En tal sistema, A es el componente clave ligero (αA,C = 2.3), C es el
componente clave pesado (αC,C = 1.0) y B es el componente intermedio
(αB,C=1.3). Las siguientes ecuaciones permiten la determinación del valor
mínimo de la razón de reflujo y de las composiciones en el destilado de los
componentes B y C a reflujo mínimo.

N α j zj
∑α = 1− q (1)
j =1 j −θ

Ec. De
N α j x Dj
∑
j =1 α j − θ
= 1 + Rmin Underwood
(2)

N (3)
∑x
j =1
D
j =1

Utilice el sistema de modelación GAMS para determinar las dos raíces para θ
en la Ecuación (1), el valor mínimo de la relación de reflujo y los valores de
xDB y xDC. Suponga que q = 1.0.

zA = 0.6 xDA = 0.8

zB = 0.3
zC = 0.1

F = 1000 Kmol/hr

Departamento de Ingeniería Química 11

GAMS Instituto Tecnológico de Celaya

Código GAMS del Ejemplo 3

$TITLE UNDERWOOD

$OFFSYMXREF
$OFFSYMLIST

*
*DEFINICION DE VARIABLES, PARAMETROS Y ECUACIONES
*

SETS
J COMPONENTS /1*3/,
I ROOTS /1*2/;
VARIABLES C;

POSITIVE VARIABLES TETA(I), XD(J), RMIN;

PARAMETERS ALFA(J), Z(J), Q;

*
*ECUACIONES
*

EQ1(I).. SUM(J,((ALFA(J)*Z(J))/(ALFA(J)-TETA(I))))=E= 1.0 - Q;

EQ2(I).. SUM(J,((ALFA(J)*XD(J))/(ALFA(J)-TETA(I))))=E= RMIN + 1.0;
EQ3.. SUM(J,XD(J))=E= 1.0;
OBJ.. C =E= 1.0;

*
* DEFINICION DE LAS ECUACIONES QUE FORMAN PARTE DEL MODELO
*

MODEL UNDEQN /ALL/;

*
* ASIGNACION DE VALORES A LOS PARAMETROS
*

ALFA('1')=2.3;
ALFA('2')=1.3;
ALFA('3')=1.0;

Z('1')=0.6;
Z('2')=0.3;
Z('3')=0.1;

Q = 1.0;

Departamento de Ingeniería Química 12

GAMS Instituto Tecnológico de Celaya

*
* VALORES INICIALES Y LIMITES INFERIOR Y SUPERIOR
*

TETA.L('1')= 1.05;
[Link]('1')= 1.299;
[Link]('1')= 1.001;
TETA.L('2')= 2.1;
[Link]('2')= 2.299;
[Link]('2')= 1.301;
XD.L('2')=0.1;
[Link]('2')=1.0;
XD.L('3')=0.01;
[Link]('3')=1.0;
[Link]('1')=0.8;

OPTION LIMROW=0;
OPTION LIMCOL=0;
*
* LLAMADO A LA TECNICA DE SOLUCION
*

Resultados GAMS del Ejemplo 3

COMPILATION TIME = 0.050 SECONDS 0.7 Mb WIN200-
121
Model Statistics SOLVE UNDEQN USING NLP FROM LINE 72

MODEL STATISTICS

BLOCKS OF EQUATIONS 4 SINGLE EQUATIONS 6

BLOCKS OF VARIABLES 4 SINGLE VARIABLES 7
NON ZERO ELEMENTS 16 NON LINEAR N-Z 10
DERIVATIVE POOL 8 CONSTANT POOL 12
CODE LENGTH 207

GENERATION TIME = 0.050 SECONDS 1.9 Mb WIN200-

121

EXECUTION TIME = 0.110 SECONDS 1.9 Mb WIN200-

121

S O L V E S U M M A R Y

Departamento de Ingeniería Química 13

GAMS Instituto Tecnológico de Celaya

MODEL UNDEQN OBJECTIVE C

TYPE NLP DIRECTION MINIMIZE
SOLVER CONOPT FROM LINE 72

**** SOLVER STATUS 1 NORMAL COMPLETION

**** MODEL STATUS 2 LOCALLY OPTIMAL
**** OBJECTIVE VALUE 1.0000

RESOURCE USAGE, LIMIT 0.488 1000.000

ITERATION COUNT, LIMIT 2 10000
EVALUATION ERRORS 0 0

C O N O P T Windows NT/95/98 version 2.043F-008-043

Copyright (C) ARKI Consulting and Development A/S
Bagsvaerdvej 246 A
DK-2880 Bagsvaerd, Denmark

Using default control program.

** Optimal solution. There are no superbasic variables.

CONOPT time Total 0.160 seconds

of which: Function evaluations 0.000 = 0.0%
Derivative evaluations 0.000 = 0.0%

Work length = 0.05 Mbytes

Estimate = 0.05 Mbytes
Max used = 0.04 Mbytes

---- EQU EQ1

LOWER LEVEL UPPER MARGINAL

1 . . . EPS
2 . . . EPS

---- EQU EQ2

LOWER LEVEL UPPER MARGINAL

1 1.000 1.000 1.000 EPS

2 1.000 1.000 1.000 EPS

Departamento de Ingeniería Química 14

GAMS Instituto Tecnológico de Celaya

LOWER LEVEL UPPER MARGINAL

---- EQU EQ3 1.000 1.000 1.000 EPS

---- EQU OBJ 1.000 1.000 1.000 1.000

LOWER LEVEL UPPER MARGINAL

---- VAR C -INF 1.000 +INF .

---- VAR TETA

LOWER LEVEL UPPER MARGINAL

1 1.001 1.039 1.299 .

2 1.301 1.539 2.299 .

---- VAR XD

LOWER LEVEL UPPER MARGINAL

1 0.800 0.800 0.800 EPS

2 . 0.167 1.000 .
3 . 0.033 1.000 .

LOWER LEVEL UPPER MARGINAL

---- VAR RMIN . 0.450 +INF .

**** REPORT SUMMARY : 0 NONOPT

0 INFEASIBLE
0 UNBOUNDED
0 ERRORS

EXECUTION TIME = 0.000 SECONDS 0.7 Mb

Departamento de Ingeniería Química 15

Introducción a la Optimización
Bajo Incertidumbre
 Reflexión
¾ Compañía petrolera: ¿cuál será el precio y la
demanda del petróleo en 6 meses?
¾ En un proceso continuo
9 Existirá variación en las demandas del producto
9 Calidad de Servicios?
9 Cadena de Suministro de Materias Primas?

El futuro no puede pronosticarse con exactitud

Necesario considerar incertidumbre en algunos

procesos: Procesos Estocásticos
Tipos de Problemas de Optimización
Estocástica

* Programación Lineal Estocástica (SLP)

* Programación Mixta Entera Lineal Estocástica

(SMILP)

* Programación No Lineal Estocástica (SNLP)

* Programación Mixta Entera No lineal Estocástica

(SMINLP)
 Otra Clasificación: Tipos de
Problemas Bajo Incertidumbre

¾ “Wait and see”

see : Esperar ocurrencia de un
evento incierto y entonces optimizar
¾ “Here and now” Optimización inmediata
en base a alguna medida de probabilidad

La mayoría de los algoritmos de

solución utilizan ambas estrategias
Problemas Estocásticos de 2 Etapas
¾ Idea fundamental: Recurso

¾ Recurso en 2 Etapas Primera Etapa

Seleccione la variable de decisión
x

Segunda Etapa
Ocurrencia de un evento incierto

Tomar una acción correctiva (recurso) y

 Un Ejemplo

El problema del vendedor de

periódicos
¾ El vendedor compra x periódicos a un precio c
¾ Entonces vende tantos periódicos como puede a un precio q,
el exceso representa una pérdida
¾ La demanda del periódico cambia día a día (incertidumbre)
¾ Cuando la demanda se conoce, se calculan las ganancias

Cuántos periódicos debe comprar el vendedor para maximizar

sus ganancias ?
El Problema del Vendedor
de Periódicos

Primera Etapa
Seleccionar el número de
periódicos a comprar x

Segunda Etapa
Ocurrencia de un evento incierto (demanda)

Las ganancias se calculan

Se toma una acción correctiva (recurso)

Programación Estocástica Lineal
con Recurso
Representación Matemática Estándar
para Problemas Lineales (SLPwR)
min c T x + Q ( x ) Función de
Primera s. t. Ax = b
Recurso
Etapa
x≥0
donde Q( x ) = Eω [Q( x,ω )] y
Evento
Q( x, ω ) = min q (ω ) y T
incierto
Segunda
s. t. W (ω ) y = h(ω ) − T (ω ) x
Etapa
Matriz de y≥0
Recurso
Clases Especiales de Problemas
Recurso W (ω ) = W
Fijo

Recurso W = ( I ,− I )
Simple y − y = h(ω ) − T (ω ) x
+ −

Recurso
W y = z ∀z , y≥0
Completo
 Reformulación

min cT x + Q( x ) min cT x + θ
s. t. Ax =b s. t. Q( x ) ≤ θ Primera
x≥0 Ax =b Etapa
x≥0

donde Q( x ) = Eω [Q( x,ω )]

Q ( x, ω ) = min q T (ω ) y

s. t. W y = h −T x
y≥0
 Dos Tipos de Cortes en
Algoritmos SLP
Corte de Optimalidad

• Aproximación Lineal de Q(x)

• Basado en el problema dual: proporciona límite inferior

a Q(x)

Corte de Factibilidad

• Asegura que los valores de x (obtenidos en la primera

etapa) no propician infactibilidades en la segunda etapa
 Teorema de la Dualidad (LP)

Multiplicadores
Primo de Lagrange
Dual
min cT y max π T b
s. t. A y = b s. t. π T A ≤ c
y≥0

• Si el dual no es acotado,
acotado el primo es infactible
• Si el dual es infactible, el primo no es acotado
• El valor de la función objetivo del problema dual provee
una cota inferior para la función objetivo del problema
primo.
primo En problemas convexos sus valores son iguales.
 Problema de la Segunda Etapa

Primo Dual
Multiplicadores
min qT y de Lagrange
max π T (h − T x )
s. t. W y = h −T x
s. t. π WT
≤q
y≥0
 Ejemplo Ilustrativo

min − 0.75 x + Eω [Q( x,ω )]

s. t. x ≤ 5
x≥0

Q( x,ω ) = min − y1 + 3 y2 + y3 + y4

s. t. − y1 + y2 − y3 + y4 = ω + 1 2 x
− y1 + y2 + y3 − y4 = 1 + ω + 1 4 x
y1 , y2 , y3 , y4 ≥ 0
 Ejemplo Ilustrativo
c = [− 0.75] x = [x]

A = [1] α →≤ b = [5]
− 1  y1 
3 y 
q=  y =  2
1  y3 
   
1  y4 
Recurso
Fijo
− 1 1 − 1 1 − 1 
 ω 
W =  α → = h(ω ) =  
T = 2
1 
 − 1 1 1 − 1 1 + ω  − 4
Dual del Problema de la Segunda
Etapa
Q( x, ω ) = max π 1 (ω + 1 2 x ) + π 2 (1 + ω + 1 4 x )

s. t. − π 1 − π 2 ≤ − 1

Multiplicadores π1 + π 2 ≤ 3
de Lagrange
− π1 + π 2 ≤ 1

π1 − π 2 ≤ 1
 Corte de Optimalidad:
Aproximación Lineal a Q(x)

Soporte Lineal

10 10

8 8

6 6

4 4

2 2

0 0
0 1 2 3 4 5 0 1 2 3 4 5

10 10

8 8

6 6

4 4

2 2

0 0
0 1 2 3 4 5 0 1 2 3 4 5
 Corte de Optimalidad
• El valor de la función objetivo del problema de la
segunda etapa en cada iteración ν (tomando xν de la
primera etapa) y para el k-ésimo valor de las variables
inciertas, ωk, es:
( ) ( ) (h
Q xν ,ω k = π νk
T
k − Tk xν )
(teorema de la dualidad)

Debido a la convexidad
(dual es Límite inferior)
(
Q x, ω k
) ≥ (π ) (h ν T
k k − Tk x )

• Para una función de probabilidad discreta, teniendo el

valor ωk una probabilidad pk, el valor esperado de la
función objetivo:
( )
ν
Qx =E π [( ) (h − T x )] = ∑ p [ (π ) (h
ν T ν
K

k =1
k
ν T
k k − Tk xν )]
 Corte de Optimalidad
Por lo tanto, debido
a la convexidad

K
Q(x ) ≥ ∑ p k π
k =1
[( ) (h
ν T
k k ] K
− Tk x ) = ∑ p k π
k =1
( )ν T
k
K

k =1
( )T
hk − ∑ p k π νk
T
k x

( )T
K
( )
K
Definiendo e = ∑ pk π ν T
k hk y E = ∑ pk π ν T
k k
k =1 k =1

Se tiene Q( x ) ≥ e − Ex

y dado que θ ≥ Q( x ) entonces Ex + θ ≥ e

θ ≥ e − Ex
 Corte de Factibilidad
• La decisión tomada en la primera etapa xν resulta en un
problema factible en la segunda etapa si existe un vector
finito y tal que las restricciones:

W y = h − T xν
y≥0

se satisfacen. Note: si y es finito, entones qTy es

finito y por lo tanto
Q( x ) < ∞
 Corte de Factibilidad
• Para verificar
factibilidad, resolver el
problema: • Cuyo problema dual
es:

z = min e y + y
T
( + −
) (
max σ T h − T xν )
s. t. W y + y+ − y− = h − T xν s. t. σT W ≤0

y ≥ 0, y + ≥ 0, y − ≥ 0 σ ≤e
Multiplicadores
de Lagrange

• Note que z>0. si z=0 entonces la segunda etapa es factible

 Corte de Factibilidad

• Sin embargo, si z>0 entonces el problema primo:

min qT y
es infactible
s. t. W y = h − T xν
y≥0
• Por el teorema de la dualidad: Si el dual no está acotado,
entonces el primo es infactible. Por lo tanto, el dual:

(
max π T h − T xν ) no estaría acotado
s. t. π W
T
≤q
Note:

(
max π T h − T xν ) no está acotado debido a que (σ ) (h − T x ) ≥ 0
ν T
 Corte de Factibilidad

• Por lo tanto, para asegurar la factibilidad del primo, la restricción:

(σ ) (h − T x ) ≤ 0
ν T
debe añadirse

Feasibility Cut

• Si para algún valor k (evento discreto) el primo es infactible,

entonces se define:
( )T
D = σν
T
k

d = (σ ) h
ν T
k

• Y se incorpora el corte de factibilidad: Dx≥d

 Algoritmos para Recurso
Fijo
• Usa función de probabilidad
discreta para ω
Método “L- • Cálculo exacto del límite
Shaped” inferior de Q(x) (Corte de
Optimalidad)
Optimalidad
• Muestreo de una función de
probabilidad continua para ω
Descomposición
• Estimación del límite inferior
Estocástica (SD)
de Q(x) basado en esperanza
matemática (Corte de
Optimalidad)
Optimalidad
 Método “L-Shaped”

• Supone una función de probabilidad discreta para ωk

• Note la estructura del problema determinístico equivalente:
K
min c x +T
∑p
k =1
k qkT yk

s. t. Ax = b
W yk = hk − Tk x k = 1K K
x ≥ 0, yk ≥ 0 Estructura del Dual
A
T1 W AT T1T T2T … TkT
WT
T2 W
. WT
.
Estructura del primo .
. .
.

Tk … W
WT
 Algoritmo “L-Shaped”

• Paso 0 Haga r = s = ν = 0

• Paso 1 Haga ν = ν+1 y resuelva el problema (Current Problem CP)

min z = cT x + θ
s. t. Ax = b

Corte de Factibilidad Dl x ≥ dl l = 1Kr

Corte de Optimalidad El x + θ ≥ el l = 1K s

xν y θν conforman la solución óptima. Si no hay cortes (iteración

1), haga θν =-∞ y no la considere en el problema
 Algoritmo “L-Shaped”
Paso 2 Para k=1…K (número de realizaciones de un
evento incierto) resuelva el problema:

+ −
z = min eT yk + eT yk
+ −
s. t. W yk + yk − yk = hk − Tk xν
+ −
yk ≥ 0, yk ≥ 0, yk ≥ 0
Si para algún k el valor óptimo es z>0 añada un corte de
factibilidad:
Dr +1 = (σ νk ) Tk
T
Multiplicadores
de Lagrange del Dr +1 x ≥ d r +1
problema ( )
d r +1 = σ ν T
k hk
anterior

Haga r=r+1 y regrese al Paso 1. De otro modo, vaya al

paso 3.
 Algoritmo “L-Shaped”

• Paso 3 Para k=1…K resuelva el problema

min qkT yk

W yk = hk − Tk xν
yk ≥ 0
Y defina:

( )h Es +1 = ∑ pk (π νk ) Tk
K K
Multiplicadores
es +1 = ∑ pk π ν T T

de Lagrange del k k
k =1
k =1
problema
anterior η ν = es +1 − E s +1 xν
si θ ν ≥ ην Pare, xν es la solución óptima
Si no haga s=s+1 , añada el corte de optimalidad
θ = es +1 − E s +1 x Y regrese el paso 1
 Descomposición Estocástica

• Se muestra en cada iteración a partir de una distribución de

probabilidad continua ωk
• Cálculo de límite inferior de Q(x) es aproximado

L-Shaped Descomposición Estocástica

( )
1 ν ν T
eν = ∑ π k hk
( )h
K
ν k =1
es +1 = ∑ pk π ν T
k k

( )
1 ν ν T
k =1
Eν = ∑ π k Tk
ν k =1
( )
K
Es +1 = ∑ pk π νk Tk
T

k =1 Actualización:
ν − 1 ν −1 ν ν − 1 ν −1
eνk = ek Ek = Ek k=1…ν-1
ν ν
 Algoritmo de Descomposición
Estocástica (Higle y Sen)
Recurso completo

• Paso 0 Haga ν = 0, θν =-∞ Suponer x1

• Paso 1 Haga ν = ν+1 y genere una observación de las

variables estocásticas mediante muestreo
 Algoritmo SD

• Paso 2 Determine θν(x) ( v-ésima aproximación lineal a Q(x) )

a) Resuelva el problema de optimización lineal (dual de segunda

etapa): (
max π T hν − Tν xν )
s. t. πT W ≤q
ν
Para obtener πν (v-ésima muestra y al v-ésimo valor del vector x)

Similarmente resuelva el problema v-1 veces:

(
max π T hk − Tk xν )
s. t. πT W ≤q k = 1Kν − 1

Para obtener π νk (k-ésima muestra y al v-ésimo valor de x)

 Algoritmo SD

b) Calcule los coeficientes del corte de optimalidad

e = ∑ (π k ) hk
ν 1 ν ν T
ν
ν k =1 v
v
v
v ( ) (h
1 ν ν
e − E x = ∑ πk
ν k =1
T
k − Tk x )
( )
ν
1
Eνν = ∑
ν k =1
π ν T
k Tk

c) Actualizar los coeficientes de previos cortes

ν − 1 ν −1
eνk = ek
ν
k = 1Kν − 1
ν − 1 ν −1
Ekν = Ek
ν
 Algoritmo SD

• Paso 3 Resuelva el problema de la primera etapa con los cortes

de optimalidad:

min cT x + θν

s. t. Ax = b Q( x )

θν ≥ eνk − Ekν x θν + Ekν x ≥ eνk

k = 1Kν
Para obtener xν+1. Vaya al paso 1

• El algoritmo se detiene si el cambio en la función objetivo es pequeño

 Un Caso
z Suponga que se tiene un SLPwR con 50 restricciones y que
se toman N = 200 muestras del valor de ω
¾ Es necesario resolver 200 problemas de optimización
correspondientes a la primera etapa. El primero posee 50
restricciones, el segundo 51…
51 el último 250 restricciones

¾ Es necesario resolver el problema de la segunda etapa un

número de veces igual a
N

∑i =
i =1
20100

¾ El número de restricciones en la segunda etapa no cambia.

 Implementación del Algoritmo
SD: Técnica de Muestreo HSS

• Monte Carlo puede presentar valores grandes de

varianza
• HSS presenta mejores propiedades de uniformidad

1 1

0.5

Y
Y

0.5

0
0
0 0.5 1
0 0.5 1
X
X

HSS Monte Carlo

 Implementación Computacional
• Integración del entorno de modelación GAMS, el
código de la técnica de muestreo HSS (FORTRAN) y
un programa en C++ como programa maestro

C++ Code Sampling

(FORTRAN)
1) Generation of an
approximation to Q(x):
sampling and multiple
generation and
solution of LP´s GAMS - OSL

2) Addition of
optimality cut and
GAMS - OSL
solution to the 1st
stage problem
Aplicaciones a Ingeniería Química
F uel
B o ile r

6 3 5 p s ig
s tre a m

p ow er

P u rc h a s e d
P re s s u re P1 T u rb in e 2 P2 p ow er
re d u c in g T u rb in e 1
(p o w er) (p o w er)
valve

C o n d e n s ate

1 9 5 p s ig s tre a m

P re s s u re
re d u c in g
valve 6 2 p s ig s trea m

Sistema Turbogenerador
 Aplicaciones a Ingeniería Química

Caso de Estudio PRIMERA ETAPA Segunda Etapa Variables

Renglones Columnas Renglones Columnas Estocásticas
Sistema Turbogenerador 1 2 21 28 4
Planeación de una Refinería 4 5 12 13 8
Planeación de una Planta 1 2 11 15 11
Petroquímica

Caso de Estudio N (MISMO ERROR

PROMEDIO)
HSS MC
Sistema Turbogenerador 60 140
Planeación de una Refinería 190 275
Planeación de una Planta Petroquímica 175 160
 Resultados
Sistema Turbo Generador
1280

1260

1240
Objective

1220

1200

1180

1160

1140
0 20 40 60 80 100
8
Iteration
7
MC HSS
6

% Error
4

0
0 20 40 60 80 100
Iteration

HSS MC
¿Cómo evaluar si el esfuerzo vale la
pena?

z Valor de la Solución Estocástica (VSS) :

Diferencia entre el valor obtenido para la función
objetivo respecto del valor obtenido si se usan
valores promedio para incertidumbres

z Valor de la Información Perfecta (VPI)

Diferencia del resultado con el valor verdadero
luego de la ocurrencia real del evento
 Resultados

Caso de Estudio VSS (%)

Sistema Turbogenerador 0.48
Planeación de una Refinería 6.85
Planeación de una Planta Petroquímica 2.22
 Discusión

z A pesar de la limitación acerca de la linealidad de

las restricciones, existen aplicaciones importantes
en planeación y calendarización de procesos

z Extensión a casos entero y no lineal

¾ BONUS (No lineal)
¾ Desarrollo actual para casos de programación
entera
Programación Estocástica Mixta
Entera Lineal
 Surgen Más Clasificaciones
para Problemas Multi-Etapa
¾ Variables Enteras en la Primera Etapa :
Se utilizan los mismos algoritmos que en
programación lineal estocástica

¾ Variables Enteras en la Segunda Etapa:

Se utilizan variaciones en el método de
“branch and bound” para permitir la
adición iterativa de los cortes de
optimalidad y factibilidad
 Problemas Enteros en la Primera
Etapa: Aplicaciones a Otras Áreas

z Seguridad en Redes de Agua Municipales

ƒ Colocación óptima de sensores para disminuir el
porcentaje de la población en riesgo tras un ataque
químico a la red

z Localización de Estaciones de Desinfección

ƒ Colocación óptima de estaciones para conservar
niveles de cloro bajo especificaciones para aguas
municipales
 Seguridad en Redes de Agua
Municipales
n n β ij
min ∑ ∑ xij + Eω [Q( x, ω )]
i =1 j =1 δ ij

s. t. xij = x ji ∀i = 1K n − 1, i ≤ j
∑(i , j ), j∈E , i≤ j xij ≤ xmax qipj = ωipδ jp
xij ∈ {0,1} ∀i, j ∈ E
n P n
Q( x, ω ) = min ∑ ∑ ∑ q (ω )ipj yipj
i =1 p =1 j =1

s. t. yipi = 1 ∀i = 1K n, p = 1K P
yipk − yipj ≤ xkj ∀(k , j ) ∈ E s.t. f kjp = 1
Localización de Estaciones de Desinfección
nb
min ∑ Wi xi + Eω [Q( x, ω )]
i =1
nb
s. t. ∑ xi ≤ nbmax qik = ωik
i =1
xi ∈ {0,1} nb
1 ni
Q(x, ω ) = min ∑ ∑ q(ω )i yik
k

i =1 ∆Ti k =1
nb ni
s. t. ∑ ∑ α ijkm yik ≤ u j j = 1K nm
i =1 k =1
nb ni m = M K M + nα − 1
− ∑ ∑ α ijkm yik ≤ −l j
i =1 k =1

yik − Yi k xi ≤ 0 i = 1K nb

yik ≥ 0 k = 1K ni
 Programación Estocástica Mixta-
Entera Lineal (Variables Enteras en
la Segunda Etapa)
Check feasibility

Compute Q(xν)
Update z
Check integrality
Check feasibility If θ < Q(xν)
Add Cut
Generate optimality cut
Branch Return to current node
Root Else
Node Check integrality Go to pendant node

z= ∞ Add Cut
θ= − ∞
Otro Tipo de Problemas
Estocásticos
 Chance Constrained Programming
z Hay algunas restricciones para las que sólo existe cierta
probabilidad de que se tengan que satisfacer
z Tales restricciones deben incluir las variables inciertas
dentro de términos lineales

Minimize Z = 4 x1 − x2 Minimize Z = 4 x1 − x2

Sujeto a: Sujeto a:

2 x1 + x2 ≤ 8 2 x1 + x2 ≤ 8
3 x2 ≤ 6
P ( x2 ≤ u ) ≤
7 x1 − x2 ≤ 4
x1 − x2 ≤ 4
x1 , x2 ≥ 0
x1 , x2 ≥ 0
Chance Constrained Programming
Introducción a la Optimización
Multiobjetivo
 Optimización Multiobjetivo (MOP)
¾ Prácticamente en cualquier área y en una variedad de
contextos se presentan problemas con múltiples objetivos
que se contraponen entre sí
¾ A este tema se le conoce también como Optimización
Vectorial y se clasifica en términos del tipo de variables y
restricciones ( MOLP,
MOLP MONLP,
MONLP etc.)

Maximizar Z = (Z1 , Z 2 , Z 3 ,K Z k )

Sujeto a:
h( x) = 0

g( x) ≤ 0
 Un Ejemplo
Un estudiante desea seleccionar la mejor escuela de ingeniería con base a
varios criterios:

Escuelas consideradas
Criterios de Selección

Rangos proporcionados por US News

 Selección de Universidad
Valores normalizados

Análisis de Resultados
 Conjunto Pareto
¾ MIT es mejor que Georgia Tech y que la Universidad de Michigan
en todos los criterios considerados. Sin embargo, Stanford, Cal
Tech, Cornell y Carnegie Mellon son mejores o no que MIT
dependiendo del criterio.
¾ La solución a una problema MOP no es un solo valor, sino un
conjunto de alternativas denominado Conjunto Pareto,
Pareto Conjunto
Preferido o Conjunto No Dominado
¾ Un grupo de 5 escuelas conforman el Conjunto Pareto en el
ejemplo
¾ Conjunto Pareto: Conjunto de alternativas que proporcionan
soluciones potenciales y representan un compromiso entre los
diferentes objetivos
 Otro Ejemplo: Fabricación
de Químicos
Minimize Z1 = 4 x1 − x2 Costo
Minimize Z 2 = −05 x1 + x2 Emisiones

Sujeto a:
Región Factible en
Durabilidad x1 ≥ 1
espacio de Decisión
Almacenamiento 2 x1 + x2 ≤ 8
Disponibilidad x2 ≤ 5
Seguridad x1 − x2 ≤ 4
x1 , x2 ≥ 0
 Otro Ejemplo: Fabricación
de Químicos
Región Factible en
espacio de Objetivos

Valores en Puntos
Frontera BAD constituye el Extremos
Conjunto Pareto
 Métodos de Solución para
MOP
¾ “Métodos Basados en la Preferencia”
Preferencia :
Determinan la solución que mejor satisface la
preferencia de quien toma las decisiones.
Reduce el tiempo y el número de alternativas
pero sufren de subjetividad y falta de información
¾ “Métodos Generadores”
Generadores Determinan el conjunto
Pareto de manera formal

La mejor estrategia es utilizar un método generador

para determinar el Conjunto Pareto y entonces usar
un método basado en la preferencia para seleccionar
la solución óptima final.
 Método Generador: Método de los
Coeficientes de Peso
¾ La idea es asociar cada función objetivo con un
coeficiente de peso y minimizar la suma “pesada” de los
objetivos
¾ El problema se convierte en una serie de problemas de
optimización de una sola función objetivo

k
Optimizar Z mult = ∑wZ
i =1
i i

Sujeto a:
h( x) = 0

g( x) ≤ 0
 Método de los Coeficientes de
Peso: Procedimiento
¾ Encuentre los óptimos individuales para cada objetivo.
Tales puntos representan los extremos del Conjunto No
Dominado. Optimizar Z1
Optimizar Z2
M
Optimizar Zk

¾ Escoja valores no negativos de los pesos y resuelva el

problema: k
Optimizar Z = wZ
mult ∑ i =1
i i

Sujeto a:
h( x) = 0
g( x) ≤ 0

¾ Analice el espacio de la función objetivo y repita

con nuevos pesos de form a que se mueva hacia
la región del conjunto Pareto que se desea
explorar
 Ejemplo Ilustrativo
w1 1
Z2 = − Z1 + Z mult
w2 w2
Minimize Z mult = w1Z1 + w2 Z 2
Sujeto a: Función objetivo “pesada”

x1 ≥ 1
2 x1 + x2 ≤ 8
x2 ≤ 5
x1 − x2 ≤ 4
x1 , x2 ≥ 0
 Método Generador: Método de
Restricciones (Constraint Method)
¾ La idea otra vez es transformar el problema multiobjetivo
a una serie de problemas de un solo objetivo
¾ Se selecciona una función objetivo que se conserva
como tal y el resto se incluye como restricciones de
desigualdad
Minimize Z mult = Z i
Minimize Z = ( Z1 , Z 2 , Z 3 ,K Z k )
Sujeto a:

Z j ≤∈ j ∀j ≠ i

h( x) = 0

g( x) ≤ 0
 Método de Restricciones
Minimize Z1 = 4 x1 − x2
Minimize Z 2 = −05 x1 + x2 Minimize Z1 = 4 x1 − x2

Sujeto a: Sujeto a:

x1 ≥ 1 Z 2 = −0.5 x1 − x2 ≤∈2
2 x1 + x2 ≤ 8 x1 ≥ 1
x2 ≤ 5 2 x1 + x2 ≤ 8
x1 − x2 ≤ 4 x2 ≤ 5
x1 , x2 ≥ 0 x1 − x2 ≤ 4
x1 , x2 ≥ 0
 Método de Restricciones

∈2 = 1
 Método Basado en la Preferencia:
Optimización Mediante Metas (Goal
Programming)
¾ Se define un valor como meta para cada función
objetivo
¾ Se crea una sola función objetivo que minimiza las
desviaciones respecto de las metas definidas

Minimize Z = ( Z1 , Z 2 , Z 3 ,K Z k ) Minimize Z goal = ∑ (δ

i
i
+
+ δ i− )
Z i − G i = δ i+ − δ i− ∀i

Sujeto a:
Se establece una meta Gi h( x) = 0
para cada objetivo Zi
g( x) ≤ 0

δi+ , δi− ≥ 0
 Goal Programming
Minimize Z1 = 4 x1 − x2 G1 = G2 = −5
Minimize Z 2 = −05 x1 + x2
∑ (δ )
2
+
Minimize Z goal = i + δ i−
i =1

Sujeto a:
Sujeto a:
4 x1 − x2 + 5 = δ1+ − δ1−
x1 ≥ 1
− 05 x1 + x2 + 5 = δ 2+ − δ 2−
2 x1 + x2 ≤ 8
x2 ≤ 5 x1 ≥ 1
x1 − x2 ≤ 4 2 x1 + x2 ≤ 8
x1 , x2 ≥ 0 x2 ≤ 5
x1 − x2 ≤ 4
Solución: x1 , x2 ≥ 0

δi+ , δi− ≥ 0
Z 1 = −1 Z 2 = 4
Control Óptimo y Optimización
Dinámica
 Problemas de Control Óptimo
¾ Proceso de solución consiste en encontrar los
perfiles de la variable de control vs tiempo de
modo que se optimice un índice particular de
medida de desempeño del sistema
Maximizar
∫ k (x , θ ) dt + S (T )
T
L=
θ 0

dx
Sujeto a: = f (x , θ ) x(0) = x0
dt

¾ Métodos de Solución Convencionales

• El Principio del Máximo
• Programación Dinámica
• Cálculo de Variación
 Problemas Históricos
¾ La reina Dido planteó el problema isoperimétrico:
isoperimétrico Encuentre el
área mayor que puede ser cubierta con un cordel de longitud fija
(L)
X X 2 y → x1
 dy 
A = ∫ y ( x) dx L = ∫ ds = ∫ (dy ) + (dx )
2 2
=∫ 1 +   dx x→t
o 0  dx 

T
Maximize A = ∫ x1 (t ) dt
o

dx1
=u x1 (0) = 0 x1 (T ) = 0
dt

T
L = ∫ 1 + u 2 dt
0

dx2
= 1+ u2 x2 (0) = 0 x2 (T ) = L
dt
 Problemas Isoperimétrico

T
Maximize A = ∫ x1 (t ) dt
o

dx1
=u x1 (0) = 0 x1 (T ) = 0
dt

dx2
= 1+ u2 x2 (0) = 0 x2 (T ) = L
dt
 Brachistochrone (Tiempo Mas Corto)

Galileo Bernoulli
 Ingeniería Química: Problema de
Destilado Máximo
Maximizar T dD T V
L=∫ dt = ∫ dt
Rt 0 dt 0 Rt + 1

Sujeto a: dxt1 V
=− x01 = Bo = F
dt Rt + 1

dxt2 V ( xt2 − xD(1) )

= x02 = xF(1)
dt Rt + 1 xt1

T V
∫0
xD(1)
Rt + 1
dt
Pureza x =
*
D T V
promedio ∫ 0 Rt + 1
dt
 El Principio del Máximo
¾ La función objetivo se reformula en la forma lineal de
Mayer
¾ Requiere la incorporación de ecuaciones diferenciales
ordinarias adicionales (ecuaciones adjuntas)
adjuntas que
representan la dinámica de las variables adjuntas
(también agregadas al problema)
¾ Se define una función Hamiltoniana (invariante en el
tiempo)
¾ El perfil óptimo se obtiene derivando la función
Hamiltoniana con respecto a la variable de control
¾ El sistema resultante es un problema de valores en la
frontera
 El Principio del Máximo
Forma Lineal
Maximizar Maximizar n
L = ∫ k (x, θ ) dt J = c x (T ) = ∑ ci xi (T )
T T

θ 0 θ i =1

dx dx
=f x(0) = x0 = f x(0) = x0
dt dt

n
Hamiltoniano H = µ f = ∑ µi f i
T

i =1

Ecuaciones y dµ n ∂f j
= −µ T
fx = − ∑ µ µ (T ) = c
Variables Adjuntas dt j =1
j
∂xi

dx
=f x(0) = x0
dt
 Problem a de Destilado
M á ximo: Principio del
M á ximo
Función objetivo es
re-escrita en forma Maximizar
Rt
L=∫
T V
Rt + 1
[ (
1 − λ xD* − xD(1) dt)]
Lagrangiana 0

Sujeto a:
dxt1 V
=− x01 = Bo = F
dt Rt + 1

dxt2 V ( xt2 − xD(1) )

= x02 = xF(1)
dt Rt + 1 xt1
 Problem a de Destilado
M á ximo: Principio del
M á ximo
Para obtener forma x =∫
V
t
3
R +1
[1 − λ (x − x )]dt
t
*
D
(1)
D
lineal de Mayer 0
t

Maximize
x3T
Rt
Sujeto a:
dxt1 V
=− x01 = Bo = F
dt Rt + 1

dxt2 V ( xt2 − xD(1) )

= x02 = xF(1)
dt Rt + 1 xt1

dx t3
dt
=
V
Rt + 1
[ (
1 − λ x D* − x D(1 ) )]
 Problem a de Destilado
M á ximo: Principio del
M á ximo V (x − x )
[1− λ(x )]
2 (1)
V V
Hamiltoniano H = −µ +µ
1
+µ
2 t D 3 *
− xD(1)
t t
Rt +1
t
(Rt +1) x 1
t
t
Rt +1
D

Ecuaciones y Perfil óptimo

Variables Adjuntas  µt2 2 
(
 x1 tx − x (1)
) − µ 1
− λ (x *
− x (1)
) + 1
( )
D t D D
dµt1 2 V xt − x D Rt =   −1
2 (1) t
= µt , µT1 = 0 ∂xD (1)
µ 2
dt ( )
(Rt + 1) xt1 2
∂Rt 
 λ − t1 
xt 
 ∂xD(1) 
V 1 − 2 
dµ t2
 ∂xt  V  ∂xD(1) 
= − µt2
− µ 3
λ 
 
2 
, µ 2
=0
dt (Rt + 1) xt
1 t
(Rt + 1)  ∂xt  T
∂H
=0
∂Rt
dµt3
= 0 , µT3 = 1
dt
 Programación Dinámica
¾ Condición de Optimalidad:
Optimalidad Aplicación del Principio
de Optimalidad de Bellman da como resultado una
ecuación diferencial parcial conocida como Ecuación
Hamilton-Jacobi-Bellman (HJB)

Maximize  ∂L ∂ L dx ti 
0=  + k ( x t ,θ t ) + ∑ 
θt  ∂ t i ∂ x i
t dt 

Maximize  ∂L ∂L 
0=  + k ( x t , θ t ) + ∑ i if
θt  ∂ t i ∂xt 

Maximize
0= [L t + k + L x f ]
θt
 Problem a de Destilado
M á ximo: Program a ción
Ecuación HJB Diná mica
 V ∂L  V  ∂L  V ( xt2 − xD(1) )  
0=
∂L Maximize
∂t
+ 
+
[ 1 − (
λ x *
D − x )]
(1)
D +
∂ 1 
−
+
+ 2  
 ∂xt + 1
Rt  t
R 1 xt  Rt 1  Rt 1 xt  

Perfil óptimo
 ∂L ∂L  xt2 − xD(1)   V   V   ∂xD(1) ∂L 1 ∂xD(1) 
0 = 1 − λ ( xD − xD ) − 1 + 2 
* (1)
 − 2
+  λ ∂R − ∂x 2 x1 ∂R 
 ∂xt ∂xt  xt
1
  ( Rt + 1) R +
  t  1 t t t t 

 2
− (1)
 ∂L
 − 1 − λ (x D* − x D(1) ) + 1
∂L 2  x t x D
∂xt  xt1  ∂xt
Mismo perfil que en el principio
Rt = −1 del máximo si las variables
 ∂L 2  adjuntas son iguales a las
∂x D 
(1)
∂xt  derivadas de la función objetivo
λ −
∂Rt  xt1  (L) con respecto a las variables
  de estado (x)
 Problemas Estocásticos
de Control Óptimo
¾ No es posible despreciar incertidumbres en algunas
aplicaciones prácticas de problemas de control óptimo:
9 En parámetros del modelo
9 En condiciones iniciales

¾ El problema estocástico de control óptimo resultante

puede ser analizado utilizando “Teoría de Opción
Real”:
Real
9 Caracterizando incertidumbres dependientes del tiempo como
Procesos de Ito
9 Usando el Lema de Ito
9 Usando las condiciones de optimalidad de Programación
Dinámica Estocástica
 Procesos de Ito
¾ Las variables estocásticas cambian con el tiempo en una forma
incierta
¾ El denominado proceso Wiener se utiliza como base para modelar
una amplia gama de procesos estocásticos más complicados.
Posee 3 propiedades:
9 Satisface la propiedad de Markov
9 Presenta incrementos independientes
9 Sus cambios en el tiempo se distribuyen normalmente
¾ Un proceso de Ito representa el incremento de una variable
estocástica en el tiempo de acuerdo con:

dx = a( x, t ) dt + b( x, t )dz

a y b son funciones conocidas y dz es el incremento de un

proceso Wiener. Note que E[dz]=0 y E[dz2]=dt
 Procesos de Ito
Algunos parámetros ingenieriles pueden representarse como
procesos de Ito:
¾ Movimiento Browniano ¾ “Mean reverting process”
dx = α dt +σ dz ( )
dx = η xavg − x dt +σ dz

¾ Movimiento Geométrico Browniano

dx = α x dt +σ x dz
 Lema de Ito
¾ Teorema Fundamental del Cálculo Estocástico
¾ Permite derivar e integrar funciones de variables estocásticas que
se comportan como procesos de Ito

dx = a( x, t ) dt + b( x, t )dz

∂F ∂F 1 ∂2F
dF = dt + dx + (dx )2

∂t ∂x 2 ∂x 2

 ∂F ∂F 1 2 ∂2 F  ∂F
dF =  + a( x, t ) + b ( x, t ) 2  dt + b( x, t ) dz
 ∂t ∂x 2 ∂x  ∂x

¾ No se desprecian algunas contribuciones de segundo orden dado

que E[dz2]=dt
Programación Dinámica Estocástica
¾ Se ha desarrollado una extension a las condiciones de
optimalidad de programación dinámica para el caso
estocástico:
Maximize
L = ∫0 k (x t , θ t ) dt
T
θt

Sujeto a: ( )
dx ti = f i x t , θ t dt + σ i dz Procesos de Ito

Condiciones de
Optimalidad:
Optimalidad Maximize  1 
0= k ( xt ,θt ) + E(dL)
θt dt

Maximize  ∂L ∂L σ i2 ∂ 2 L ∂2L 
0=  + k ( xt , t ) + ∑ i fi ( xt , t ) + ∑ + ∑ σ iσ j i j 
θt  ∂t i ∂xt i 2 i 2
(∂xt ) i ≠ j ∂xt ∂xt 
 Principio del Máximo para Problemas
Estocásticos
¾ Con base en las condiciones de optimalidad para programación
dinámica, se pudieron derivar las expresiones correspondientes al
método del principio del máximo

¾ Las variables adjuntas (µ) en el principio del máximo son equivalentes

a las derivadas parciales de la función objetivo con respecto a las
variables de estado (Lx) de programación dinámica

¾ El principal resultado del análisis es la derivación de las ecuaciones

adjuntas

¾ Las derivadas de segundo orden de la función objetivo con especto

a las variables de estado (Lxx) en programación dinámica
estocástica tiene que ser también incluidas y se incorporan en la
formulación a través de las variables adjuntas adicionales, ω
 Principio del Máximo para Problemas
Estocásticos
¾ Se utiliza representación escalar aunque el análisis es válido para
el caso vectorial

H =µ f σ2
H =µ f + ω
2
dx
=f x(0) = x0 dx = f dt + σ dz x ( 0 ) = x0
dt
dµ
dt
= −µ fx µ (T ) = c dµ
dt
= −µ fx −
1
2
σ ( )2
x ω µ (T ) = c

dω
dt
= − 2 ω f x − µ f xx −
1
2
( )
σ 2
xx ω ω (T ) = 0

Determinístico Estocástico
 Versión Estocástica del Problema
de Destilado Máximo
Maximize dD T T V
L=∫ dt = ∫0 dt
Rt dt 0
Rt + 1
Sujeto a:
dxt1 V
T V =− x01 = Bo = F
∫0 Rt + 1 dt *
(1)
x D dt Rt + 1
xDave = T = xD
V
∫0 Rt + 1 dt dx t2 =
V ( x t2 − x D(1) )
dt + x t2 σ 2 dz 2 x 02 = x F(1)
Rt + 1 1
xt
Restricción externa usada Proceso Ito
como criterio de
convergencia
 Volatilidad Relativa como un
Proceso de Ito
3.2
3.15 R igorous
3.1 S im ulation

R e lativ e V o latility
3.05 P ath 1
3
2.95 P ath 2
2.9
2.85 P ath 3
2.8
2.75
0 1 2 3
T im e (H rs)

Ocasiona un comportamiento incierto en las

variables de estado
 Principio del Máximo
 ∂x D(1) 
V 1 − 2 
dµ
= −µ
2
2 V xt − x D (
(1)
) −µ  ∂xt 
− σ 2 2
2 xt ω
dt (Rt + 1) xt1
2
( ) (Rt + 1) xt1
Ecuaciones
 ∂xD(1)  ∂ 2 xD(1)
Adjuntas V 1 − 2  V
dω
= −2 ω  ∂x t 
+µ
∂x 2 2
t 2 ( )
− σ 2 ω − ωµ
V xt2 − xD(1) ( )
dt (Rt + 1) xt
1
(Rt + 1) xt 1
(Rt + 1) xt1
2
( )

( )  2
 ∂σ 2  ( R + 1)
xt1 σ 2 2 xt2 ω   t
Perfil
Rt =
xt1 (
− µ xt2 − xD(1) ) + 
∂ Rt  V

 −1
óptimo ∂xD(1) ∂xD(1)
µ µ
∂Rt ∂Rt
 Perfil Óptimo de la Razón de Reflujo

¾ Se requieren valores de reflujo más grandes debido a la

disminución en el valor de la volatilidad relativa con el tiempo
¾ La desviación respecto al caso determinístico también cambia
con el tiempo debido al efecto de las incertidumbres
 Nuevamente el Problema Isoperimétrico
¾ Considere ahora la versión estocástica del problema
isoperimétrico

Determinístico Estocástico

Maximize
x3 (T )
Movimiento Browniano
u

dx1 dx1 = u dt + σdz x1 (0) = 0 x1 (T ) = 0

=u x1 (0) = 0 x1 (T ) = 0
dt
dx2 σ = 0.5
= 1+ u2 x2 (0) = 0 x2 (T ) = L
dt
dx3
= x1 x3 (0) = 0
dt

Suposición meramente
L = 16
académica
Soluciones al Problema Isoperimétrico

µ1 = −t + c1
µ 2 = c2
µ3 = 1
ω=0
u
µ1 + µ2 = 0
2
1+ u
 Bibliografía
1. Teoría de optimización (determinística) y
Aplicaciones en Ingeniería Química
a) Practical Methods of Optimization; R. Fletcher, 2nd. Ed., Wiley
b) Optimization of Chemical Processes; Edgar, Himmelblau and Larson,
2nd. Ed., McGraw-Hill
c) Nonlinear Programming, Theory and Algorithms; Bazaraa, Sherali and
Shetty, Wiley
d) Systematic Methods for Chemical Process Design; Biegler, Grossmann
and Westerberg, Prentice Hall
e) Linear Programming; Chvatal Vasek, Ed. W. H. Freeman and Co.

2. Programación MultiObjetivo
a) Introduction to Applied Optimization, Diwekar, Kluwer Academic
Publishers
 Bibliografía
3. Programación Estocástica
a) Stochastic Programming, Kall and Wallace, Wiley

4. Control óptimo
a) Batch Distillation, Simulation, Optimal Design and Control; Diwekar, Ed.
Taylor and Francis
b) Optimal Control Theory; Sethi and Thompson, Kluwer Academic
Publishers
c) Investment Under Uncertainty; Dixit and Pindyck, Princeton University
Press