1 1.1 Elementos del problema de Optimización Hemos establecido que la optimización consiste en la búsqueda de un grupo de variables manipulables con la finalidad de maximizar o minizar una característica de un sistema. Desde un punto de vista de Ingeniería, la optimización es posible cuando hemos modelado la carcaterística deseable del sistema, también denominada función objetivo. La mayoría de los problemas de optimización tienen una estructura que es similar, lo que permite darles un tratamiento generalizado. Los elementos de un problema de optimización son: •
1º Categoría Una función objetivo a minimizar o maximizar. Esta función objetivo, dentro del problema de optimización, recibe el nombre de modelo económico.
•
2º Categoría Un conjunto de restricciones de igualdad (ecuaciones), que corresponden a identidades matemáticas impuestas por el diseño, el equipo o el modelo matemático del sistema.
•
3º Categoría Un conjunto de restricciones de desiguladad (inecuaciones), que corresponden a limitaciones físicas (o de operación) impuestas por el diseño, el equipo o el modelo matemático del sistema. Durante el curso utilizaremos la siguiente notación:
minimizar: f ( x )
sujeto a : h( x ) = 0 : g( x ) ≥ 0
función objetivo vector de restricciones independientes de igualdad. vector de restricciones independientes de desigualdad.
en esta notación: x = x1 x2 ... xn T
en
general
es
es el vector de variables manipulables en un problema de optimización.
suficiente revisar el problema de minimización, pues má xτ ( x) ≡ mín - τ ( x) . Por n denotaremos el número de variables en optimización, por m1 denotaremos el número de restricciones de igualdad, y por m2 el número de restricciones de desigualdad.
2
Se tienen los siguientes casos:
Caso Interpretación m1 < n número de incógnitas es mayor que el número de ecuaciones que las relaciona. El sistema es subdeterminado, o con infinitas soluciones. De ese conjunto de soluciones, sólo una satisface el problema de optimización. m1 = n número de incógnitas es igual al número de ecuaciones que las relaciona. El sistema es determinado y sólo se requiere resolver el sistema de restricciones de igualdad. Esta solución es, por definición, la solución del problema de optimización. m1 > n número de incógnitas es menor que el número de ecuaciones que las relaciona. El sistema es sobredeterminado, ningún conjunto de soluciones satisface todas las retricciones de igualdad. Lo mismo sucederá para el problema de optimización. Figura 1.1 Gráfica de una función objetivo bidimensional
20 18 16 12 10
3
4
2 1 3
2
1
0
0 -1 -2 -1
-2
2
8 6 4 2 0 4
x
f (x1, x2)
14
-3 -3
-4
-4
x1
Las restricciones operan entre variables de optimización, de modo que ellas aparecen en el plano de nivel x1 , x2 . Las restricciones pueden ser funciones de constantes o simplemente funciones lineales o no lineales.
3
Figura 1.2 Restricciones en el plano de nivel x1-x2 4 g3(x) g4(x) 2 g2(x) Dominio factible
x2
g5(x)
h1(x)
0
g1(x) -2 Curva de soluciones satisfactorias
-4 -4
-3
-2
-1
0
1
2
3
4
x1
Ejemplo Nº 1:
Planteamiento de un problema de optimización:
Se dispone del siguiente equipo existente, para separar un componente valioso:
4 se necesitan determinar las condiciones de operación que maximizarán la recuperación del sistema.
Etapa Nº 1: se formula la función objetivo y el modelo que relaciona las variables.
R=
m1,out V2 y1,2 Fψ1ψ 2 y1,2 ψ1ψ 2 y1,2 = = = m1,in Fz1,1 Fz1,1 z1,1
deseamos obtener la máxima recuperación, es decir, nuestra función objetivo puede escribirse:
minimizar: f ( x ) = −
ψ1ψ 2 y1, 2 z1,1
recordamos ahora que, suponiendo que las corrientes alcanzan, las fracciones vaporizadas pueden escribirse en términos de las constantes de vaporización como: ψ=
z1 z2 − K1 − 1 1 − K2
y las composiciones de equilibrio pueden calcularse como:
xi , j =
(K
i, j
zi , j
)
−1 ψ +1
; yi , j =
(K
zi , j Ki , j
i, j
)
−1 ψ +1
suponiendo ahora solución ideal, tenemos que:
Pi sat Ki = ⇒ Ki = Ki ( P, T ) P de acuerdo con el modelo indicado, tenemos que:
(
f ( x ) = f y1, 2 , z1,1 , T1 , T2 , P1 , P2
)
Etapa Nº 2: Establecer las restricciones de igualdad.
Puesto que el equipo existe, su geometría está completamente definida. Las características de su alimentación estarán impuestas por la planta que precede a la separación:
5
⎫ z1,1 = z1º,1 h1 = z1,1 − z1º,1 = 0 ⎪ P1 = P º − ∆P1 ⎬ ⇒ h2 = P1 − P º + ∆P1 = 0 P2 = P1 − ∆P2 ⎪⎭ h3 = P2 − P1 + ∆P2 = 0 Obviamente tenemos acá ecuaciones y modelos para el cálculo de la pérdida de carga. Etapa Nº 3: Establecer las restricciones de desigualdad.
al separar el producto, podemos necesitarlo con una pureza mínima:
y1,2 ≥ esp. ⇒ g1 = y1,2 − esp. ≥ 0 de igual forma, la flexibilidad del equipo (su diseño estructural, o su capacidad de intercambio de calor), podrá soportar un rango determinado de temperaturas y presiones, lo que define rangos factibles para T y P.
Una vez que se ha definido completamente el problema se procede a seleccionar una técnica de optimización, lo que dice estricta relación con la naturaleza de las funciones involucradas. La función objetivo puede ser continua (o continua por secciones), para lo que nos podría servir una técnica con derivadas. Si es lineal, existen técnicas especiales para tratarla, como la programación lineal. También puede ser discreta, por ejemplo: el diámetro de una cañería o el área de intercambio de calor. En este caso podemos seleccionar un modelo que represente continuamente a la función objetivo, y seleccionar el valor más próximo al diámetro o al área. Existen casos discretos que no pueden resolverse con la técnica anterior, por ejemplo, las piezas de una máquina son discretas. Si se quiere resolver el problema de dar la mayor durabilidad de operación al mínimo costo, es necesario combinar calidades. Existe para este caso una técnica especial de optimización, denominada programación entera. En Ingeniería Química, muchos problemas clásicos pueden resolverse utilizando funciones continuas o continuas por tramos. Dedicaremos una parte abundante del programa al tratamiento de tales funciones.
6
1.2 Teoría básica de Funciones Continuidad Una función unidimensional es continua en un punto xº si:
a. f ( x0 ) existe. b. lím f ( x ) existe, y es único aproximándose direccionalmente a xº x → x0
c. lím f ( x ) = f ( x0 ) x → x0
el concepto de continuidad tiene gran importancia en optimización, pues permite seleccionar una técnica ad-hoc.
Función discontinua (a)
Función continua (b) 70
50
60
40
50 40 f(x)
f(x)
30 20
30 20
10
Punto de discontinuidad
0
10 0
0
2
4
6 x
8
10
12
0
2
4
6
x
El caso (a) podría producirse por un cambio de tecnología, por ejemplo, la eficiencia de dos compresores recomendados para ciertos rangos de flujo. Se observa la presencia de un máximo que podría ser alcanzado con técnicas diferenciales, dependiendo del valor inicial. El caso (b), corresponde a un caso continuo, este tipo de funciones se obtienen comúnmente al aplicar valor absoluto (por ejemplo, en ajuste de parámetros). Aunque la función es comtinua por definición, no es diferenciable. Una técnica de derivadas puede fallar.
7
Funciones uni y multimodales.
Otro aspecto de gran importancia cuando se buscan puntos candidatos al óptimo, dice relación con la monotonía de las funciones. En optimización, la monotonía se denomina modalidad de una función. Una función se define unimodal en un dominio definido cuando presenta un solo punto estacionario (máximo o mínimo) en el dominio de estudio. Por otro lado, se dice que una función es multimodal en un dominio definido cuando presenta varios puntos estacionarios. El caso (a), por ejemplo, es unimodal, pues presenta un solo punto estacionario, que correponde a un máximo. El caso (b) presenta mínimos y máximos, dando origen a una función multimodal. Cuando las funciones son continuas y de clase C2 , la unimodalidad puede ser estudiada analizando la segunda derivada. La función será unimodal cuando el signo de su segunda derivada no cambia. Desde el punto de vista de la optimización, el caso más sencillo ocurre cuando las funciones objetivo son unimodales.
Funciones cóncavas y convexas
Una función es matemáticamente cóncava en una regíon R, si al tomar dos puntos de la región, se cumple que:
f ( θxa + ( 1 − θ ) xb ) ≥ θf ( xa ) + ( 1 − θ ) f ( xb ) ; 0 ≤ θ ≤ 1 desde un punto de vista gráfico, esto signfica que el hiperplano que une dos puntos de la función evaluada se sitúa por debajo de la función. Una función convexa se define invirtiendo la desigualdad:
f ( θxa + ( 1 − θ ) xb ) ≤ θf ( xa ) + ( 1 − θ ) f ( xb ) ; 0 ≤ θ ≤ 1 en funciones objetivo unidimensionales, la situación gráfica es la siguiente:
8
Función cóncava
Función convexa
16
0
14 12
-2 xa
-4 -6 xb
8
f(x)
f(x)
10
-8
6
-10
4
-12
2
-14
0
-16 0
1
2
3
4
5
6
7
xb
xa
0
1
2
x
3
4
5
6
7
x
el concepto de concavidad y convexidad, permite clasificar el punto extremo (máximo o mínimo) como un extremo global. La técnica más poderosa para estudiar la concavidad o convexidad de una función de clase C2 , es el análisis de la segunda derivada.
⎛d 2 f ⎜ 2 ⎝ dx
⎧> 0 función convexa ⎞⎪ ⎟: ⎨< 0 función cóncava ⎠⎪ ⎩= 0 punto de silla (o de inflexión)
las funciones multidimensionales pueden ser aproximadas localmente por expansiones en serie de Taylor de funciones cuadráticas:
f ( x ) = f ( x 0 ) + ∇ T f ( x 0 )( x − x 0 ) +
T 1 x − x 0 ) H( x − x 0 ) ( 2
Donde H es la matriz de segundas derivadas parciales, o Hessiana:
⎡ ∂2f ⎢ ∂x 2 ⎢ 21 ⎢ ∂ f ( ) H x = ⎢ ∂x2∂x1 ⎢ 2M ⎢ ∂ f ⎢⎣ ∂xn∂x1
∂2f ∂x1∂x2 ∂2f ∂x22 M
∂2f ∂xn∂x2
∂2f ⎤ ... ∂x1∂xn ⎥ ⎥ ∂2f ⎥ ... ∂x2∂xn ⎥
O ...
M
⎥
∂2f ⎥ ∂xn2 ⎥⎦
9 recordamos que para funciones diferenciables, se cumple la identidad: ∂2 f ∂2 f ≡ ∂xi∂x j ∂x j∂xi
de donde se deduce que la Hessiana es una matriz simétrica. El estudio de la concavidad o convexidad de una función multidimensional en un punto x0 , puede extenderse de lo que ya sabemos para las ecuaciones cuadráticas. La función es convexa cuando el Hessiano es definido positivo y cóncava cuando es definido negativo. Definiciones:
1. Una matriz es definida positiva ssi
x Mx > 0 ∀ x ≠ 0
2. Una matriz es definida negativa ssi
x Mx < 0 ∀ x ≠ 0
3. Una matriz es indefinida ssi
⎧> 0 para x a ≠ 0 x T M x⎨ ⎩< 0 para x b ≠ 0 ; x b ≠ x a
4. Una matriz es semidefinida positiva ssi
x Mx ≥ 0 ∀ x ≠ 0
5. Una matriz es semidefinida negativa ssi
x Mx ≤ 0 ∀ x ≠ 0
T T
T
T
En base a lo anterior, y haciendo la analogía con funciones cuadráticas a partir de la expansión en serie de Taylor, la definición del Hessiano permite averiguar la característica cóncava y convexa de una función. • • • •
Una función es cóncava si H es semidefinido negativo. Una función es estrictamente cóncava si H es definido negativo. Una función es convexa si H es semidefinido positivo. Una función es estrictamente convexa si H es definido positivo.
Teorema:
Una función es estricamente convexa en un punto x0 , cuando los valores propios de la Hessiana en ese punto son positivos. Demostración:
Puesto que una función es convexa cuando es definida positiva, se tiene: x0 H ( x0 ) x0 = χ > 0 T
10
donde χ es una constante escalar positiva, que también puede ser escrita como: χ = x0 λ x 0 T
puesto que χ es positivo, λ también debe serlo, por aparecer como multiplicador de una norma de vectores. Notamos ahora que: x0 H ( x0 ) x0 = x0 λ x0 ⇒ H ( x0 ) x0 = λ Ix0 T
T
donde I es la matriz identidad de dimensión equivalente a la Hessiana. En esta última ecuación, puesto que x0 es un vector no nulo, se deduce que:
H − λI = 0 Esta última ecuación define el problema de los valores característicos de una matriz. Tomando determinante a ambos lados, se obtiene el polinomio característico de la matriz H, cuyas raíces son precisamente los valores propios. Recordamos que el polinomio característico tiene el mismo orden que su matriz asociada, por tanto el número de valores característicos corresponde a la dimensión de la matriz. Como la Hessiana es simétrica, todos sus valores propios son reales. Corolarios del teorema: • • • •
Una función es convexa cuando todos los valores propios son mayores o iguales que cero (Hessiana semidefinida positiva). Una función es estrictamente cóncava todos los valores propios son negativos. Una función es cóncava cuando todos los valores propios son menores o iguales que cero. Cuando una función combina simultáneamente valores propios negativos y positivos, o bien, cuando todos los valores propios son nulos, su geometría es silla de montar.
11
Gráfica de una función indefinida
Ejemplo : Determine si la siguiente función es cóncava o convexa
f ( x ) = 2 x12 + 2 x1 x2 + 1.5 x22 + 7 x1 + 8 x2 + 24 Tomando las primeras derivadas parciales, se tiene:
∂f = 4 x1 + 2 x2 + 7 ∂x1 ∂f = 2 x1 + 3x2 + 8 ∂x2
12 De las segundas derivadas parciales: ∂2 f =4 ∂x12 ∂2 f =3 ∂x22 ∂2 f ∂2 f = =2 ∂x1∂x2 ∂x2∂x1 de esta forma, la Hessiana asociada a f (x) será:
⎡ 4 2⎤ H =⎢ ⎣ 2 3⎥⎦ es decir, una matriz de constantes para todos los reales. Calculamos ahora los valores propios de H:
2 ⎤ 2 ⎡4 − λ = λ − 7λ + 8 = 0 det H − λ I = det⎢ 3 − λ⎥⎦ ⎣ 2
(
)
los valores propios son: λ1 = 5.5616 λ 2 = 1. 4384 puesto que ambos son positivos, la función es estrictamente convexa en todo ℜ. En otras palabras, esta función tiene un solo mínimo, este mínimo es global y único.
13
Gráfica de una función objetivo bidimensional
140 120
80 60 40 3 2
20
1 3
2
1
0
0 -1 -2 -1
-2
x
0
4
2
f (x1, x2)
100
-3 -3
-4
-4
x1
Una segunda forma que permite determinar la característica de la Hessiana, es el de los determinantes de los menores principales: •
Cuando los elementos de la diagonal de la Hessiana son positivos, y los determinantes de todos los menores principales son positivos, la Hessiana es definida positiva.
Regiones cóncavas y convexas
Cuando se optimiza una función objetivo sujeta a restricciones, existe un dominio factible en donde se encuentra la solución. Se define un dominio convexo como aquel cuyo hiperplano que une dos elementos del dominio, tiene todos sus puntos contenidos en el dominio:
14
Dominio convexo
xb
xa ℜ
Notamos en la figura que:
x = θ x a + ( 1 − θ ) x b ( con 0 ≤ θ ≤ 1) ∈ ℜ En el mismo caso anterior, cuando la región factible de optimización es externa a la elipse, se tendrá: Ejemplo de región cóncava
xa ℜ
xb
en este caso, algunos de los puntos de la línea (hiperplano) que une dos elementos del dominio, no estan contenidos en el dominio. Teoremas (topología matemática):
1. La intersección de dos regiones convexas siempre es convexa. 2. La unión de dos regiones convexas, no necesariamente define una región convexa.
15
ℜa U ℜ b
Las regiones cóncavas y convexas adquieren relativa importancia cuando se combinan restricciones de desigualdad en un mismo problema de optimización.
16 2. Optimización de funciones continuas sin restricciones
El problema de minimización de funciones continuas sin restricciones se reduce a: minimizar: f ( x ) utilizando como dominio factible a todos los reales. En clases anteriores indicamos que las condiciones necesarias y suficientes para la existencia de máximos y mínimos es: a. la función es diferenciable en el mínimo, al menos dos veces. b. ∇f ( x0 ) = 0
c. H ( x ) es definida en el punto candidato.
también hemos indicado que una función minimizable en el sentido de las condiciones expuestas, es aproximable por una expansión de Taylor cuadrática (o de mayor orden) en las vecindades del mínimo:
1 f ( x) = f ( x 0 ) + ∇T f ( x 0 ) ∆ x + ∆ x T H( x 0 ) ∆ x 2 ∆ x = x − x0
(2.1)
La ecuación (2.1) representa exactamente a la función cuando esta es cuadrática, y puede ser utilizada para el cálculo del mínimo. En particular, tenemos:
∇f ( x ) = ∇f ( x 0 ) + H ( x 0 ) ∆ x = 0 ⇒ H ( x 0 ) ∆ x = −∇f ( x 0 )
(2.2)
cuando (2.2) representa exactamente a la función cuadrática, la minimización es directa y no requiere iteraciones. Cuando la función no es cuadrática, la Hessiana y el gradiente cambian con las coordenadas, (2.2) puede ser utilizada en forma iterativa para la localización del mínimo:
( )
( )
H x j ∆ x j +1 = −∇f x j
x j +1 = x j + ∆ x j +1
(2.3)
La ecuación (2.3) es muy eficiente cuando es "facil" obtener un conjunto de derivadas analíticas de la función objetivo, pues en ella se reconoce el algoritmo clásico del método de Newton, cuyas propiedades de convergencia nos son familiares: depende del punto inicial y de la funcionalidad de las derivadas en las variables. En un problema real de optimización, la tarea de obtener derivadas consume bastante tiempo, y no siempre la ecuación (2.3) es una alternativa del todo adecuada. A modo de ejemplo, analicemos la función:
17
f ( x1 , x2 ) = 25 x14 + 18 x22 + 23 x1 x2 + 5 x1 + 3 x2 − 7 para este caso particular, se tiene:
⎡ ∂f / ∂x1 ⎤ ⎡ 5 + 100 x13 + 23 x2 ⎤ ⎡ 0 ⎤ = ∇f ( x ) = ⎢ = ⎣∂f / ∂x2 ⎥⎦ ⎢⎣ 3 + 23 x1 + 36 x2 ⎥⎦ ⎢⎣ 0 ⎥⎦ con la Hessiana respectiva:
⎡ 300 x12 23⎤ H( x) = ⎢ 36 ⎥⎦ ⎣ 23 obviamente ahora la definición de la Hessiana es variable y depende de x1 ; puesto que la función no es cuadrática, la Hessiana tampoco es una matriz de elementos constantes. Figura 2.1: Gráfica de la función
En la figura 2.2 se puede apreciar el desempeño del método de Newton considerando dos puntos de iniciación. El problema converge al mínimo rápidamente al tomar x T0 = ( −0.5 −0.5) , pero una situación diferente se tiene al tomar x T0 = ( 0 1) . En la misma figura es posible apreciar que el método va generando direcciones secuenciales de búsqueda.
18
Figura 2.2 Minimización de una función utilizando la técnica de Newton (Sensitividad al punto de partida)
Una opción interesante al método anterior, en el que la función es analítica, es definir direcciones de búsqueda, y minimizar a lo largo de esas direcciones. Analicemos la ventaja de esta alternativa: Dada s, una dirección de búsqueda, se obtiene la trayectoria de la línea de búsqueda utilizando la ecuación: x = x0 + ts
(2.1)
en donde t es un parámetro de variación. La ecuación (2.1) es la ecuación de la recta en el hiperespacio. Notemos ahora que, al tomar una dirección definida:
19
f ( x ) = f ( x 0 + t s) = f ( t )
(2.2)
la ecuación (2.2) indica que al definir una dirección, y conocido un punto inicial de búsqueda, la función multivariable se transforma en una función de una sola variable. Esta situación es claramente ventajosa, pues minimizar funciones de una dimensión es un problema más sencillo. A modo de ejemplo, considérese la proyección de f (x) en dos direcciones distintas, a partir de x T0 = ( 1 1)
88
sT = ( 1 -1 ) 66
fs(t) 44
22 sT = ( 1 1 )
0
-1
0
1
2
3
t
claramente se aprecia que una de las direcciones sT = ( 1 1) es más conveniente, pues es estrictamente cóncava, en la otra dirección se aprecia un cambio de concavidad, que es justamente la que está provocando problemas en el método de Newton al tomar x T0 = ( 0 1) . El segundo enfoque entrega información direccional de la función, lo que lo hace más adecuado para objetivos generales. Desde el punto de vista de la optimización, la minimización por direcciones definidas es mucho más cómoda. Estrategia de optimización de funciones continuas n-dimensionales 1. 2. 3. 4.
Definir un punto inicial. Definir una dirección. Minimizar a lo largo de la dirección. Repetir los pasos 2. , 3. hasta obtener el mínimo
2.1 Técnicas de minimización de funciones en una dimensión.
20
En la sección anterior, hemos visto que independientemente de la dimensión de la función a optimizar, la minimización en una dimensión es una componente del problema general de optimización. Este capítulo se dedica a la revisión de los métodos de minimización unidireccionales. La técnicas de minimización unidireccional pueden clasificarse en dos familias: •
Métodos directos: son aquellos que minimizan utilizando evaluaciones secuenciales de la función.
•
Métodos indirectos: cuando las funciones son de clase Cn , la búsqueda de puntos estacionarios es equivalente a resolver la nulidad del gradiente. Al ser este el caso, no se requiere directamente evaluar la función, sino más bien, iterar sobre sus derivadas.
La selección de una u otra técnica depende de la diferenciabilidad de la función en el dominio factible, que posibilitará el cálculo y la evaluación de las derivadas. Métodos Directos. a. Métodos de Contracción del Dominio. La base de estos métodos consiste en acotar una región unimodal del dominio, y muestrear la función en un grupo de puntos. Generalmente, estos métodos proceden sobre un máximo de cuatro muestreos, para evitar las evaluaciones excesivas de la función:
Al cumplirse la condición:
( )
( )
( )
f x j < f x j +1 < f x j + 2 con x j < x j +1 <
15 14
f (x1)
13
con j = 1; en el supuesto de que la función es unimodal, es posible eliminar la región ( x1 , x 2 ) , para acotar el mínimo a la región (x2 , x 4 ) y volver a producir una subdivisión de cuatro muestreos en el dominio. Cuando j = 2, es posible eliminar la región ( x 3 , x 4 ) , y reducir el estudio a la región ( x1 , x 3 ) .
f (x)
12 11
f (x2) f (x4)
10 9 f (x3) 8 7 -2
-1
0
1 x
2
3
4
Esta técnica puede ser interpretada como un acotamiento del mínimo, que procederá hasta que: x1 − x 4 ≤ ξ
21 Si bien la subdivisión del dominio unimodal es completamente arbitraria, existen técnicas que permiten seleccionarlo de modo de evaluar la función un mínimo de veces, pues obviamente el número de evaluaciones de la función incide sobre el tiempo requerido para la convergencia. Técnica de la sección áurea
Consideremos un trazo de longitud L, con subdivisiones internas como se indica:
la sección áurea corresponde a la siguiente ley de proporciones: x+y y = y x con x+y= L
(2.3)
cuya solución para raices positivas es:
x=
3− 5 5 −1 L ; y= L 2 2
(2.4)
de esta forma, los puntos de muestreo para una función serán:
xa ; xb =
3− 5 5 −1 ( ( xd − xa ) + xa ; xd xd − xa ) + xa ; xc = 2 2
(2.5)
por cada etapa de reducción de dominio, la sección áurea contrae el dominio de búsqueda al 50( 5 − 1) = 61. 8% . cuando f ( x a ) > f ( x b ) > f ( x c ) se elimina la región x a → x b , revisando los puntos ( xb , xc , xd ) , los que siguen en proporción áurea, bastando introducir un solo nuevo punto. De lo contrario se elimina la
22 región x c → x d , revisando los puntos ( xa , xb , xc ) , los que siguen en proporción áurea, bastando introducir un solo nuevo punto.
Ejemplo Consideremos la minimización de la función f ( x ) = x ; como se sabe, esta función no es diferenciable en cero, donde presenta un mínimo absoluto:
1. se busca una región donde la función sea convexa, de modo de acotar el mínimo. En nuestro caso, el mínimo está acotado por -1 < x < 1. 2. se subdivide el dominio de acuerdo a la regla (2.5)
3− 5 5 −1 ( ( xd − xa ) + xa = 0.2361 xd − xa ) + xa = −.2361; xc = 2 2
xa = −1; xd = 1; xb =
2.0
1.5
y
x
x 1.0 f (x) 0.5
0.0 región eliminada
-0.5 xa
xb
xc
xd
-1.0 -1.5
-1.0
-0.5
0.0
0.5
1.0
1.5
x
xa
xc
xd
como f ( x a ) > f ( x b ) > f ( x c ) no se cumple, entonces se elimina la región x c → x d , y se inicia un nuevo ciclo de eliminación. Los puntos del dominio restante se encuentran en proporción áurea, de modo que para repetir la nomenclatura del ciclo se tiene:
xa1 = xa0
xd1 = xc0
xc1 = xb0
f ( xa1 ) = f ( xa0 )
f ( xd1 ) = f ( xc0 )
f ( xc1 ) = f ( xb0 )
3− 5 1 ( xd − xa1 ) + xa1 2 evaluar → f ( xb1 )
xb1 =
23 si f ( x a ) > f ( x b ) > f ( x c ) se cumple, entonces se elimina la región x a → x b , los puntos vuelven a encontrarse en sección áurea de modo que:
xa1 = xb0
xd1 = xd0
xb1 = xc0
f ( xa1 ) = f ( xb0 )
f ( xd1 ) = f ( xd0 )
f ( xb1 ) = f ( xc0 )
5 −1 1 ( xd − xa1 ) + xa1 2 evaluar → f ( xc1 )
xc1 =
en cualquiera de los dos casos es posible observar que se requiere solo una nueva evaluación por ciclo. xa
i 0 1 2 3 4 5 6 7 8 9 10
-1.0000 -1.0000 → -0.5279 → -0.2361 -0.2361 → -0.1246 → -0.0557 -0.0557 → -0.0294 → -0.0132 -0.0132 →
xb
xc
-0.2361 -0.5279 -0.2361 -0.0557 -0.1246 -0.0557 -0.0132 -0.0294 -0.0132 -0.0031 -0.0069
0.2361 -0.2361 -0.0557 0.0557 -0.0557 -0.0132 0.0132 -0.0132 -0.0031 0.0031 -0.0031
xd
1.0000 → 0.2361 0.2361 0.2361 → 0.0557 0.0557 0.0557 → 0.0132 0.0132 0.0132 → 0.0031
f (xa ) 1.0000 1.0000 → 0.5279 → 0.2361 0.2361 → 0.1246 → 0.0557 0.0557 → 0.0294 → 0.0132 0.0132 →
f (x b ) 0.2361 0.5279 0.2361 0.0557 0.1246 0.0557 0.0132 0.0294 0.0132 0.0031 0.0069
f (xc ) 0.2361 0.2361 0.0557 0.0557 0.0557 0.0132 0.0132 0.0132 0.0031 0.0031 0.0031
f (xd ) 1.0000 → 0.2361 0.2361 0.2361 → 0.0557 0.0557 0.0557 → 0.0132 0.0132 0.0132 → 0.0031
Propiedades de convergencia de la sección áurea
Por cada ciclo iterativo, la sección áurea contrae el dominio al 61.8%, entonces se tiene:
L1 = 0.618 L0 = 0.618( xd0 − xa0 ) L2 = 0.618 L1 = 0.618 2 L0 M Ln = 0.618 n L0
(2.6)
pero Ln = ( xdn − xan ) = ξ , donde ξ es un criterio de convergencia. Tomando en cuenta lo anterior, el número de ciclos para llegar a un error determinado es:
ln ξ − ln( xd0 − xa0 ) n= ln 0.618
(2.7)
para contraer el dominio al 1% se requiere de unas 12 iteraciones, lo que corresponde a 15 evaluaciones de la función. El orden de convergencia del método de la sección áurea es de 3/2, bastante comparable al método de Newton, sumando además la ventaja de que la
24 función no tiene porqué ser diferenciable, y que en algunos casos, ni siquiera se requiere una función continua. Es absolutamente necesario, en todo caso, que la función sea convexa en el dominio inicial, y unimodal (monótona) antes y después del mínimo. b. Métodos de interpolación
Una sucesión de puntos no colineales de una función continua, puede ser aproximada por un polinomio de grado (n-1), donde n corresponde al número de puntos no colineales. En tales casos puede ser de utilidad la interpolación cuadrática a la cúbica. •
Interpolación cuadrática
cuando se conocen 3 puntos de tendencia convexa de una función objetivo uni-dimensional, según se muestra en la figura:
f(x)
f(xa) f(xc) f(xb)
x
es posible ajustar una función cuadrática, donde: ax a2 + bx a + c = f ( x a ) ax 2b + bx b + c = f ( x b ) ax c2 + bx c + c = f ( x c ) cuya solución parcial es:
a=
xa ( f ( xb ) − f ( xc )) + xb ( f ( xc ) − f ( xa )) + xc ( f ( xa ) − f ( xb ))
( xb − xa )( xc − xa )( xc − xb ) xa2 ( f ( xb ) − f ( xc )) + xb2 ( f ( xc ) − f ( xa )) + xc2 ( f ( xa ) − f ( xb )) b= ( xb − xa )( xc − xa )( xc − xb ) una cuadrática como la indicada, tiene punto extremo cuando:
25
1 xa2 ( f ( xb ) − f ( xc )) + xb2 ( f ( xc ) − f ( xa )) + xc2 ( f ( xa ) − f ( xb )) b =− x =− 2a 2 xa ( f ( xb ) − f ( xc )) + xb ( f ( xc ) − f ( xa )) + xc ( f ( xa ) − f ( xb )) *
(2.8)
la ecuación anterior es exacta cuando el ajuste procede sobre la cuadrática, pero puede proceder iterativamente cuando la ecuación original no lo es:
(
x *j +1 = x *j +1 xa , j , xb , j , xc , j
)
(2.9)
Casos :
⎧ si f(x b ) > f ( x * ) y f(x* ) < f ( xc ) nuevo rango: x b x * xc xb < x < xc ⇒ ⎨ * * * ⎩ si f(x b ) < f ( x ) y f(x ) < f ( xc ) nuevo rango: x a xb x *
⎧ si f(x a ) > f ( x * ) y f(x* ) < f ( xb ) nuevo rango: xa x * xb xa < x < xb ⇒ ⎨ * * * ⎩ si f(x a ) > f ( x ) y f(x ) > f ( xb ) nuevo rango: x xb xc *
La velocidad de convergencia del método es variable, depende fundamentalmente del grado cuadrático de la función objetivo. •
Interpolación cúbica
Al igual que en el caso anterior, la base del método es minimizar el ajuste de una función, ahora se trata de un polinomio cúbico:
f ( x) = a 3x3 + a 2 x2 + a1x + a 0
(2.10)
dada la evaluación de 4 puntos de una región convexa del dominio, los coeficientes polinómicos se calculan resolviendo el sistema:
⎡ xa3 xa2 ⎢x3 x2 ⎢ b3 b2 ⎢ xc xc ⎢⎣ xd3 xd2
xa xb xc xd
1⎤⎡ a3 ⎤ ⎡ f a ⎤ 1⎥⎢ a2 ⎥ ⎢ f b ⎥ ⎥⎢ ⎥ = ⎢ ⎥ 1⎥⎢ a1 ⎥ ⎢ f c ⎥ 1⎥⎦⎣ a0 ⎦ ⎣ f d ⎦
(2.11)
el punto estacionario en (2.10) puede obtenerse de:
−2a 2 ± 4a 22 − 12a 3a1 df 2 * = 3a 3x + 2a 2 x + a1 = 0 ⇒ x = dx 6a 3
eligiendo x* de modo que:
(2.12)
26
d2f dx2
= 6a 3x* + 2a 2 > 0
(2.13)
x*
la ecuación (2.13) es representativa de una región convexa para la cúbica de aproximación. Cuando se dispone además de información diferencial, bastan dos puntos para evaluar la cúbica, y :
⎡ f b' + w − z ⎤ ( x − xa ) x = xb − ⎢ ' ⎣ f b − f a' + 2w ⎥⎦ b *
z=
3( f a − f b )
( xb − xa )
w = ( z 2 − f b' f a' )
(2.14)
+ f b' + f a' 1/ 2
obviamente, los puntos a y b son representativos de una región convexa para la que si xa < x b , entonces fa' < 0 y fb' > 0. El requisito de la función es la convexidad y la continuidad. Métodos indirectos
Los métodos indirectos se caracterizan por no evaluar directamente la función a minimizar, sino más bien, su información diferencial (primera y segunda derivada). Todos están basados en el Método de Newton y algunas modificaciones numéricas. Método de Newton
El problema de la minimización es equivalente a la anulación puntual de la primera derivada, con requisitos necesarios y suficientes para que el punto corresponda a un mínimo: f ' ( x) = 0 f '' ( x) ≥ 0
(2.15)
generalmente ocurre que en un modelo no lineal, la primera derivada es también no lineal, de manera que puede ser anulada utilizando el método ordinario de iteración NewtonRaphson. Consideremos una expansión de Taylor de segundo orden, en la vecindad de un punto x k
1 f ( x ) = f ( x k ) + f ' ( x k )( x − x k ) + f '' ( x k )( x − x k )2 2
(2.15)
27
el mínimo debe cumplir la condición necesaria: f ' ( x ) = 0 = f ' ( x k ) + f '' ( x k )( x − x k ) ó x = xk −
f ' (xk ) f '' ( x k )
(2.16)
el mínimo en (2.16), corresponde al mínimo de la aproximación afín de segundo orden. Si f(x) es cuadrática, la solución es exacta; pero si no lo es, se puede utilizar la fórmula recursiva de iteración: x k +1 = x k −
f ' (xk ) f '' ( x k )
(2.17) Método Quasi-Newton
El método Quasi-Newton consiste en aproximar numéricamente las derivadas en (2.17). Se sabe que: f ( x + h) − f ( x − h) h→ 0 2h + − 2 f ( x) + f ( x − h) f x h ( ) f '' ( x ) = lím h→ 0 h2 f ' ( x ) = lím
(2.18)
reemplazando (2.18) en (2.17), con h suficientemente pequeño, se obtiene: x k +1 = x k −
h f ( x k + h) − f ( x k − h) 2 f ( x k + h) − 2 f ( x k ) + f ( x k − h)
(2.19)
notamos que el método requiere tres evaluciones de la función por iteración. Si el costo de tiempo de evalución de la función es elevado (por ejemplo, una función objetivo resultante de un cálculo no lineal), el método es poco recomendable. Otro aspecto que debe tenerse en cuenta, es que la fórmula (2.19) es una aproximación, por tanto, de ella no se obtiene el mínimo exacto, sino que solamente el punto que anula a la primera derivada numérica de la función.
28
Método de la secante (variable métrica unidimensional)
En este caso, la primera derivada de la función a minimizar es conocida y analítica, lo que no se conoce es su segunda derivada. Consideremos la redefinición: g( x) = f ' ( x) (2.20)
g '( x) = f '' ( x) x k +1 = x k −
g( x k ) g'( x k )
la derivada g'(x) puede ser aproximada por una recta que una dos puntos de la curva: g' ( x k ) =
g ( x k ) − g ( x k −1 ) x k − x k −1
(2.21)
de esta forma: x k +1 = x k −
g( x k ) x k − x k −1 g ( x k ) − g ( x k −1 )
(2.21)
Interpretación gráfica del método
g(x)
0
xk
xk+1
xk-1
La partida del método requiere de la definición de dos puntos iniciales, posteriormente la construcción de la secuencia es automática. La ecuación (2.21) se indefine numéricamente en el mínimo exacto, lo que hace recomendable utilizarla como: x k +1 = x k −
g( x k ) x k − x k −1 g ( x k ) − g ( x k −1 ) + ∆
(2.22)
29
puesto que el método de Newton es el que define todas las búsquedas indirectas, es obvio esperar que estos métodos conserven las mismas ventajas y desventajas de la metodología materna.
2.2 Acotado de regiones funcionales convexas
Para todos los métodos anteriores, fundamentalmente los directos, hemos discutido que es necesario acotar un dominio donde la función tiene característica convexa. En general, esto no es sencillo cuando la función no tiene restricciones, y puede ser comparado a bajar un cerro en la oscuridad. Sintomáticamente, la primera información la da la derivada local de un punto de partida, pues si: df > 0 : la función es minimizable en el sentido en que x decrece dx df < 0 : la función es minimizable en el sentido en que x crece dx
(2.23)
posteriormente, es necesario muestrear la función en el sentido de la minimización, hasta encontrar una condición en la que: f ( x k +1 ) > f ( x k ) ∧ f ( x k ) < f ( x k −1 ) cuando x k+1 > x k > x k −1 ó
(2.24)
x k +1 < x k < x k −1
Obviamente, 2.24 no garantiza que la función sea convexa: 40 30 20 10 f(x)
0
-10 -20 -30 xk-1
xk
xk+1
-40 -3
-2
-1
0
1
2
3
x
en este caso, cualquier método podría dar solución sobre un mínimo falso de la función (o relativo). La metodología que veremos aquí, permite solucionar el problema matemático de
30 encontrar un punto estacionario, pero no garantiza si estamos acotando un mínimo relativo o uno absoluto. Generalmente es conveniente perturbar la solución de un problema de minimización, para establecer la calidad del punto encontrado. Una técnica que permite acotar el mínimo, consiste en utilizar el muestreo:
xk +1 = xk + sgn( f '( x0 ))δ × 2k
(2.25)
donde δ es un factor de expansión que lleva el signo de la primera derivada en el punto inicial de búsqueda (o bracketting). Esta técnica puede ser combinada con cualquiera de los métodos que hemos visto, por ejemplo, la sección áurea entre muestreos sucesivos, o información sobre las derivadas superiores de f(x).
31 Capítulo #3 Optimización de funciones N-Dimensionales
En el capítulo anterior, se introdujeron algunas ideas básicas relativas a la optimización de funciones de múltiples dimesiones, lo que en resumen correspondía a utilizar minimizaciones unidimensionales secuenciales a través de múltiples direcciones estrategicamente seleccionadas. Estrategia de optimización de funciones continuas n-dimensionales 1. 2. 3. 4.
Definir un punto inicial. Definir una dirección. Minimizar a lo largo de la dirección. Repetir los pasos 2. , 3. hasta obtener el mínimo
Los diferentes algoritmos de optimización para este tipo de funciones, difieren únicamente en la metodología utilizada para ir generando direcciones secuenciales. Al igual que en el caso de la minimización 1-D, los algoritmos se clasifican como directos e indirectos. 3.1 Métodos directos
Los métodos directos se caracterizan por evaluar directamente una función para obtener información acerca del mínimo. Es necesario, entonces, definir una estrategia para que con el mínimo número de evaluciones de la función, sea posible obtener el punto óptimo.
3.1.a. Métodos de Malla El método de malla consiste en evaluar una función en los vértices de una geometría predefinida. A modo de ejemplo, en dos dimensiones, esta geometría puede ser un cuadrado regular. Se busca el vértice de menor valor para f (x), que constituirá el centroide para la construcción de un nuevo cuadrado regular. La geometría cuadrada regular suele denominarse diseño factorial central, y es muy común para definir una estrategia de experimentación en dos variables. Nótese que por diseño, son necesarias 3n − 1 evaluaciones de la función (n: número de variables de dependecia), situación que en funciones de muchas variables hace prohibitivo el método.
32
x
1
x
0
3.1.a. Método Simplex Es Simplex es una variación del método de malla, en este caso, se utiliza una figura geométrica cerrada de (n+1) vértices y lados iguales, donde n es el número de variables independientes de la función. En el caso de dos dimensiones, esta figura es un triángulo equilátero. la función se evalúa en los (n+1) vétices, y se proyecta como imagen especular del de mayor valor: x2
x1
f(x 3 ) < f(x 1 ) y f(x 2 ) < f(x 1 ) x3
33 En la medida que el método avanza hacia el mínimo, generalmente ocurrirá que la figura geométrica pasará por los mismos puntos. Cuando este es el caso, se procede a la compactación los lados, por ejemplo, bisección, hasta una toleracia aceptable para la presición del mínimo. Uno de los códigos más utilizados en paquetes comerciales de optimización, es el método de Nelder-Mead (1965)1, que toma las ideas expuestas aquí, pero considera métodos más eficientes para la proyección del polígono equilátero, como para su compactación. 3.2 Métodos Direccionales
Los métodos direccionales pueden ser calificados de directos o indirectos de acuerdo a la información que utilicen para avanzar hacia un mínimo. La base de estos métodos consiste en minimizar unidimensionalmente una función, utilizando direcciones predefinidas. Se define la proyección de una función en una dirección s según:
x s = x 0 + qs
( )
(3.2.1)
f s ( x) = f x s = f ( q)
en este caso, q corresponde al parámetro de la función, en tanto que los x s corresponden a la ecuación de una recta en el espacio. Cuando una función N-Dimensional se proyecta en una dirección, se transforma en una función exclusiva del parámetro q, según se observa en (3.2.1), una función unidimensional. El mínimo direccional de f (x) proyectado en s, debe cumplir la condición de mínimo, es decir:
df ( q ) ≡0 dq
d 2 f ( q) ≥0 dq 2
(3.2.2)
Ejemplo
Considere la función:
f ( x1 , x2 ) = 2 x12 + x22 − 3
analice la proyección de la función en la dirección s T = −4 −2 , considerando el punto x 0T = 1 1 , y obtenga también el mínimo direccional. Solución:
1Nelder,
J. and Mead, R., 1965. A Simplex Method for Function Minimization. The Computer Journal,7:308
34 El lugar geométrico de las variables en la dirección, se calcula por:
⎡1⎤ ⎡ −4 ⎤ x = x 0 + qs = ⎢ ⎥ + q⎢ ⎥ ⎣1⎦ ⎣ −2⎦ ó
x1 = 1 − 4 q x2 = 1 − 2q
entonces la proyección de f en s será:
f s ( x ) = 2( 1 − 4q ) + ( 1 − 2q ) − 3 2
2
= 36q 2 − 20q una parábola convexa. Minimizando f en función de q, es posible obtener el mínimo direccional: dfs 5 = 72 q − 20 = 0 ⇒ q = dq 18 y el mínimo direccional de f en s es:
⎡1⎤ 5 ⎡ −4 ⎤ ⎡ −1 / 9 ⎤ x* = ⎢ ⎥ + ⎢ ⎥ = ⎢ ⎣1⎦ 18 ⎣ −2⎦ ⎣ 4 / 9 ⎥⎦
Teorema:
El gradiente de una función n-dimensional evaluado en el mínimo direccional, es ortogonal a la dirección. Demostración:
El mínimo direccional se obtiene al resolver la ecuación: dfs ( q ) =0 dq puesto que: dfs ( q ) df ( x( q )) = dq dq
35
utilizando la regla de la cadena para la derivación, tenemos:
df ( x (q )) ∂f dx1 ∂f dx2 ∂f dxn = + +L+ dq ∂x1 dq ∂x2 dq ∂xn dq ⎡ dx1 ⎤ ⎢ dq ⎥ ⎥ ⎢ ⎡ ∂f ∂f ∂f ⎤⎢ dx2 ⎥ =⎢ L ∂xn ⎥⎦⎢ dq ⎥ ⎣ ∂x1 ∂x2 ⎢ dxM ⎥ ⎢ n⎥ ⎣ dq ⎦ dx = ∇ T f ( x) • dq dado que: x = x 0 + qs
se deduce dx =s dq ó
df ( x (q )) dx = ∇ T f ( x (q )) • = ∇ T f ( x (q )) • s dq dq
en el mínimo direccional, se tiene:
df ( x (q )) = 0 = ∇ T f ( x * (qmín )) • s dq ecuación que indica ortogonalidad. Veamos que esto se verifica para el ejemplo:
f ( x1 , x2 ) = 2 x12 + x22 − 3 el mínimo direccional encontrado a partir de x 0T = 1 1 en la dirección s T = −4 −2 es:
⎡ −1 / 9 ⎤ x* = ⎢ ⎣ 4 / 9 ⎥⎦ por tanto:
36
⎡ 4 x1 ⎤ ⎡ −4 / 9 ⎤ ∇f ( x * ) = ⎢ ⎥ = ⎢ ⎣ 2 x2 ⎦ ⎣ 8 / 9 ⎥⎦ y
⎡ −4 ⎤ ∇ T f ( x * ) • s = [ −4 / 9 8 / 9] • ⎢ ⎥ = 0 ⎣ −2⎦ interpretación geométrica del teorema:
x0
x*
37
x0
38
1.0
Y Axis
0.5
0.0
-0.5
-1.0 -1.0
-0.5
0.0 X Axis
0.5
1.0
39
40
41 3.3 Métodos direccionales directos: Técnica de Powell Hemos demostrado que para funciones exactamente cuadráticas, es posible la construcción de un conjunto conveniente de direcciones, denominadas direcciones conjugadas, sin la participación de derivadas superiores o gradientes. Recordamos en este punto que, cuando se seleccionan dos puntos de partida xa y x b , y una dirección común de minimización s , la dirección generada por xa − x b es exactamente conjugada a s si la función materna es cuadrática.
Figura 3.3.1 Generación de direcciones conjugadas
s
x
x
a
b
x
* a
x
* b
s
Si bien el método puede generar direcciones conjugadas para funciones cuadráticas, es necesario recordar que, puntualmente, una función arbitraria admite ajustes cuadráticos. Por esta razón, una dirección conjugada sera una buena dirección en tanto la función original sea aproximable cuadráticamente en un determinado rango. En el método de Powell se aplica una versión modificada de esta idea, que automáticamente comenzará a generar direcciones conjugadas en la medida que el método avanza. El método parte considerando minimizaciones direccionales en los unitarios cartesianos. En funciones bidimensionales, la construcción de direcciones se aprecia en la figura 3.3.2. En una etapa inicial, se parte con minimizaciones en términos de los unitarios cartesianos, de modo de generar una dirección s03 , que representa una dirección global de la etapa. Minimizando a lo largo de esa dirección, se obtiene el punto inicial de una nueva 1 etapa x0 . En esta nueva etapa, se reemplaza una de las direcciones (generalmente la de
42 j+ 1 mayor progreso) por la dirección global generada en la etapa anterior, de manera que s3 conjugará a s3j .
Figura 3.3.2 Generación de direcciones conjugadas para el método de Powell
s1 3
1 =x 2
x1 0
x12 0
1
s
x10 0
0
s0
3
=x
2
-x
s11 = s01
2
=s
3
x11
0
x02 ( s02 )T = [ 0 1 ]
( s01 )T = [ 1 0 ] x00
x01
Todas las minimizaciones componentes de cada etapa, son unidimensionales, y pueden efecturase por una técnica directa, como la sección áurea, dando origen a un método directo de minimización. Algoritmo del método de Powell Etapa 0 Dada una base ortogonal de "n" vectores, donde "n" corresponde al número de variables en optimización, y dado un punto de partida x00 , se procede a minimizar en todas las direcciones: punto inicial x00
mínimo direccional x10
x10
dirección s = 1 0 L 0 s = 0 1 L 0
M x0n −1
M s = 0 0 L 1
M x0n
0 1 0 2
0 n
la función se evalúa además en el punto: x0n +1 = 2 x0n − x00
x20
43 que se sitúa en la recta de la dirección global. Ir a 2. Etapa k 1. Dado un vector inicial x0k , y un conjunto de "n" direcciones, minimizar en términos de todas las direcciones utilizando una técnica unidimensional directa: punto inicial x0k x1k
dirección s = 1 0 L 0 s = 0 1 L 0
mínimo direccional x1k x2k
M xnk −1
M s = 0 0 L 1
M xnk
k 1 k 2
k n
2. La ilustración gráfica del método de Powell, muestra que la secuencia de direcciones va cambiando por cada etapa de minimización, de manera de generar por cada etapa una nueva dirección conjugada. Si en alguna de las direcciones de una etapa el avance es poco significativo, numéricamente es posible que las direcciones globales se hagan colineales entre dos etapas sucesivas. Cuando este es el caso, el método converge exactamente al mismo mínimo de la etapa anterior, que no necesariamente es el mínimo funcional. Powell sugiere realizar el cálculo de la función en el punto: f ( x nk +1 ) = f ( 2 x nk − x0k ) desde un punto de vista matemático, esta evaluación tiene por objetivo perturbar el mínimo de la etapa, de manera de testear su calidad de mínimo funcional. 3. Se define la variación ∆k según:
{
}
∆k = má x f ( x ik− 1 ) − f ( x ik )
esta operación representa la mayor evolución en el decrecimiento de la función objetivo, asociado a una de las direcciones de la etapa, la que denotaremos por smk . Se definen los siguientes puntos de referencia: f1k = f ( x0k ) f2k = f ( x nk ) f3k = f ( x nk +1 ) = f ( 2 x nk − x0k ) y los siguientes criterios:
λ1 = f3k − f1k
λ2 = ( f 1k − 2 f 2k + f 3k )( f 1k − f 2k − ∆k ) − ∆k ( f 1k − f 2k ) 1 2
2
44 ambos criterios miden la posibilidad de minimizar la función en la dirección global generada por la etapa. 4.a Si λ1 ≥ 0 o λ 2 ≥ 0 , la función no será minimizable en la dirección global asociada a la k-ésima etapa, lo que indica que probablemente la direcciones globales en dos etapas sucesivas se han hecho colineales. Si este es el caso, entonces para la nueva etapa se seleccionan: x0k +1 = x nk skj +1 = skj en otras palabras, se eligen las mismas direcciones de la etapa anterior, para una nueva etapa que se construye sobre el último mínimo direccional anterior. 4.b Si λ1 < 0 y λ 2 < 0 , se procede como sigue: i ). definir la dirección
s*k = x nk − x0k
ii ). se minimiza la función desde xnk utilizando la dirección definida en i ). El punto obtenido en esta minimización será el punto de partida para la nueva etapa, es decir el mínimo direccional será x0k +1 iii ). se reemplaza la dirección smk , representativa del mayor avance de minimimización en la etapa k. La secuencia de direcciones para la etapa k+1 será: s1k +1 s2k +1 L snk +1 = s1k L smk −1 smk +1 L
s*k
de acuerdo al procedimiento descrito aquí, las direcciones globales irán conjugando por etapa. 5. repetir la etapa "k" hasta que: x nk − x0k < ξ
Cuando las direcciones para el método de Powell se obtienen por el procedimiento 4.b, la convergencia del método es cuadrática, pues se están utilizando direcciones conjugadas, que son muy eficientes. En particular, puesto que cada etapa va generando direcciones conjugadas por 4.b, una función cuadrática n-dimesional requiere exactamente de n etapas de Powell. En funciones generales, es posible esperar una convergencia menor que la cuadrática, dependiendo de la definición local de la Hessiana.
45 Brent (1973), propone una técnica similar a la de Powell, pero más eficiente y con mayor carga computacional. En el método de Brent se hace un ajuste cuadrático local de la función, y se utilizan los vectores principales del ajuste. La convergencia de este método es cuadrática, aún en problemas pobremente condicionados.
Ventajas del método de Powell •
La aplicación del método es sencilla, y aconsejable para aquellos problemas en los que es muy difícil obtener información diferencial.
•
La convergencia del método es dominantemente cuadrática en la mayoría de las etapas, debido al uso de direcciones conjugadas. Desventajas
•
La calidad de la técnica de acotado del mínimo direccional, como el método de minimización unidimensional, tienen un gran impacto en la eficiencia global.
•
La minimización unidimensional puede llegar a requerir excesivas evaluaciones de la función objetivo. Si bien, la convergencia del método es cuadrática por etapas, el tiempo requerido en completar una etapa puede anular las ventajas de la convergencia.
•
En funciones objetivo pobremente condicionadas, como aquellas que exhiben Hessianas locales indefinidas, la secuencia de direcciones globales por etapa no es conjugada, originando un avance tediosamente lento del método.
Ejemplo
Una función Test para métodos de optimización es la función de Rosenbrok, definida como:
f ( x ) = 100( x2 − x12 ) + ( 1 − x1 ) 2
2
de apariencia inocente, la función combina puntos de indefinición, y dos mínimos, uno de los cuales es local. El aspecto amenazante de su superficie se oberva en la figura 3.3.3. En este ejemplo se aplicará el método de Powell a la minimización de esta superficie, partiendo del punto xT0 = −1. 2 1 . El mínimo funcional absoluto se encuentra en el punto x*T = 1 1
46
Figura 3.3.3: Superficie de la función f ( x ) = 100( x2 − x12 ) + ( 1 − x1 ) 2
2
En la tabla a continuación, se muestran las iteraciones por etapa del método de Powell. Es interesante notar que la optimización procede sobre un punto de silla, aunque pese a esto, el método puede llegar al mínimo. En la figura 3.3.4 se muestra la trayectoria de minimización, la que procede bordeando la sima de la función. Es importante hacer notar que el proceso completo requiere de aproximadamente 1500 evaluaciones sobre la función objetivo, situación que lo hace poco saludable para funciones objetivos iterables
47
Etapa
x1
x2
1
-1.2000 -0.9950 -0.9950
1.0000 1.0000 0.9900
2
-0.9897 -0.9897 -0.9843
0.9897 0.9794 0.9791
3
-0.7611 -0.7611 -0.7015
0.5397 0.5792 0.4621
4
-0.5028 -0.5383 -0.4656
0.2027 0.2725 0.1776
5
-0.1959 -0.2276 -0.1617
-0.0041 0.0374 -0.0070
6
0.0850 0.0560 0.1164
-0.0284 -0.0089 -0.0141
7
0.3432 0.3177 0.3713
0.0893 0.0915 0.1159
8
0.5739 0.5524 0.5982
0.3081 0.2983 0.3417
9
0.7733 0.7567 0.7929
0.5841 0.5684 0.6184
sT
∆
[1.0000 [0.0000 λ1 =
0.0000] 1.0000] -8.3275
f (x) 24.2000 3.9900 3.9799 λ2 =
[0.0000 [0.9988 λ1 =
1.0000] -0.0488] -0.0426
3.9693 3.9587 3.9481 λ2 =
0.0106 0.0107 0.0000
[0.0000 [0.4529 λ1 =
1.0000] -0.8915] -0.4833
3.2570 3.1013 2.9858 λ2 =
0.1557 0.1155 -0.0174
[0.4529 [0.6082 λ1 =
-0.8915] -0.7938] -0.3740
2.5101 2.3962 2.3014 λ2 =
0.1139 0.0948 -0.0076
[0.6082 [0.8294 λ1 =
-0.7938] -0.5587] -0.2700
1.6102 1.5280 1.4597 λ2 =
0.0823 0.0683 -0.0029
[0.8294 [0.9963 λ1 =
-0.5587] -0.0863] -0.1910
0.9642 0.9055 0.8572 λ2 =
0.0587 0.0483 -0.0010
[0.9963 [0.9099 λ1 =
-0.0863] 0.4148] -0.1231
0.5127 0.4745 0.4434 λ2 =
0.0382 0.0310 -0.0003
[0.9099 [0.7257 λ1 =
0.4148] 0.6880] -0.0697
0.2270 0.2051 0.1876 λ2 =
0.0219 0.0175 -0.0001
[0.7257 [0.5856 λ1 =
0.6880] 0.8106] -0.0303
0.0709 0.0610 0.0535 λ2 =
0.0098 0.0076 0.0000
20.2100 0.0101 -700.7541
Continúa..
48
Etapa
x1
x2
10
0.9321 0.9225 0.9444
0.8626 0.8493 0.8877
11
0.9942 0.9989 1.0000
0.9895 0.9977 1.0000
12
1.0000 1.0000 1.0000
1.0000 1.0000 1.0000
sT
∆
[0.5856 [0.4952 λ1 =
0.8106] 0.8688] -0.0060
f (x) 0.0084 0.0063 0.0048 λ2 =
[0.4952 [0.4396 λ1 =
0.8688] 0.8982] 0.0000
0.0001 0.0000 0.0000 λ2 =
0.0001 0.0000 0.0000
[0.4952 [0.4396 λ1 =
0.8688] 0.8982] 0.0000
0.0000 0.0000 0.0000 λ2 =
0.0000 0.0000 0.0000
0.0022 0.0015 0.0000
49 Figura 3.3.4 : Minimización de la función f ( x ) = 100( x2 − x12 ) + ( 1 − x1 ) Trayectoria del método de Powell (1965) 2
x2
x1
2
50 3.4 Técnicas indirectas de primer orden
En algunos problemas de optimización de funciones continuas, es posible obtener, ya sea analítica o numéricamente, alguna información sobre de carácter diferencial, en particular el gradiente. Como ya se ha explicado, el uso del gradiente negativo como dirección de optimización no es recomendable, pues el avance de esta dirección produce trayectorias escalonadas hacia el mínimo. A pesar de lo anterior, la minimización unidimensional se ve positivamente acelerada por el conocimiento del gradiente al principio de un proceso de minimización. En este punto, es conveniente recordar que las direcciones conjugadas son beneficiosas para llegar al mínimo de las funciones exactamente cuadráticas. De acuerdo con el método del gradiente, la técnica de minimización será:
( )
s k = −∇f x k
(3.4.1)
de modo que las ecuaciones de proyección direccional estén dadas por:
( )
x k + 1 = x k + qs k = x k − q∇f x k
(3.4.2)
la sucesión de direcciones generadas por (3.4.2) es mutuamente ortogonal, lo que genera el avance escalonado al mínimo. Este tipo de avance resulta inconveniente en las vecindades del mínimo en la generalidad de las funciones, pues el grado de avance es muy poco significativo, dando como resultado un número apreciable de evaluaciones de la función. Procesos de bi-ortogonalización (Polak, 1971)1
El objetivo de este preámbulo es mostrar cómo a partir de una sucesión de vectores direccionales mutuamente ortogonales, es también posible la construcción de un grupo de vectores conjugantes de un Hessiano. Recordamos que el método del gradiente permite obtener una sucesión de vectores mutuamente ortogonales, a partir de los cuales queremos asociar un conjunto de direcciones conjugadas. Consideremos una matriz A de grado N y de característica simétrica, y un par de vectores arbitrarios g1 y s1 , de modo que: g1 = s1
(3.4.3)
Polack sugiere construir un conjunto de vectores asociados g, s de acuerdo a las ecuaciones:
1Polak,
E., 1971. Computational Methods in Optimization. New York: Academic Press.
51 gi +1 = gi − τ i Asi si +1 = gi +1 + κ i si
1< i ≤ N
(3.4.4)
de modo que: g i +1 • g i = 0 T
(3.4.5)
si +1 A si = 0 T
de acuerdo a la definición (3.4.5) los vectores g resultan secuencialmente ortogonales, en tanto que los vectores s resultan secuencialmente conjugantes de la matriz A . Haciendo los cálculos correspondientes, se encuentra: T
gi gi +1 = 0 = gi gi − τi gi Asi ⇒ τi = T
T
gi gi
T
T
gi Asi T
s Asi +1 = 0 = s Agi +1 + κ s Asi ⇒ κ i = − T i
T i
T i i
si Agi +1 T
si Asi
(3.4.6)
Polack demuestra que el conjunto de N vectores g, s generados por (3.4.6) y (3.4.4), no solamente son secuencialmente ortogonales y conjugados respectivamente, sino que además son mutuamente ortogonales y conjugados. Esta demostración se puede lograr con mucha paciencia, utilizando inducción matemática. El proceso de bi-ortogonalización es una extensión de método de ortogonalización de Gramm-Smith. Ilustración: considere la matriz simétrica
⎡ 2 1 1⎤ A = ⎢ 1 3 −2 ⎥ ⎢ ⎥ ⎣ 1 −2 4 ⎦ y el vector: g1 = s1 = 1 1 1 T
T
utilice el método de Polack para obtener una secuencia de vectores g ortogonales, y un conjunto de vectores s conjugantes de A .
52 De (3.4.4) y (3.4.6), se deduce: T
gi +1 = gi − τi Asi = gi −
gi g i T
gi Asi
Asi
T
si +1 = gi +1 + κ i si = gi +1 −
si Agi +1 T
si Asi
si
aplicando directamente las ecuaciones, tenemos:
⎡1⎤ 1 1 1]⎢1⎥ [ ⎢ ⎥ ⎡ 4 ⎤ ⎡ −1 / 3 ⎤ ⎡ 2 1 1⎤⎡1⎤ ⎡1⎤ ⎡1⎤ gT g 1 ⎦ ⎣ 1 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ g = g − T1 1 As1 = 1 − 1 3 −2 1 = 1 − 2 = 1 / 3 2 1 ⎥ ⎥⎢ ⎥ ⎢ ⎥ 3 ⎢ ⎥ ⎢ ⎢ ⎥ ⎡ 2 1 1⎤⎡1⎤ ⎢ g As1 1 − 1 2 4 1 1 3 0 ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎣1⎦ [1 1 1]⎢⎢ 1 3 −2⎥⎥⎢⎢1⎥⎥ ⎣ 1 −2 4 ⎦⎣1⎦ ⎡ 2 1 1⎤⎡ −1 / 3⎤ ⎥ ⎥⎢ ⎢ 1 1 1] 1 3 −2 1 / 3 [ ⎥ ⎡1⎤ ⎡ −7 / 27⎤ ⎥⎢ ⎢ ⎡ −1 / 3 ⎤ s1T Ag 0⎦ ⎢ ⎥ ⎢ ⎣ 1 −2 4 ⎦⎣ ⎥ ⎢ 2 s2 = g − T s1 = 1 / 3 − 1 = 11 / 27⎥ 2 ⎥ ⎥ ⎢ ⎡ 2 1 1⎤⎡1⎤ ⎢ ⎥ ⎢ s1 As1 0⎦ ⎣1⎦ ⎣ 2 / 27⎦ ⎣ ⎥ ⎢ ⎥ ⎢ [1 1 1]⎢ 1 3 −2⎥⎢1⎥ ⎣ 1 −2 4 ⎦⎣1⎦ de igual forma, se calcula que: T
g3 = g2 −
g 2 g2 T
g2 As2
As2 = −7 / 23 −7 / 23 14 / 23
T
T
s3 = g3 −
s2 Ag3 T 2
s As2
s2 = −504 / 529 378 / 529 420 / 529
T
es fácil comprobar que los vectores g son mutuamente ortogonales, y los vectores h son mutuamente conjugantes de la matriz A . Supongamos conocida una base de vectores g mutuamente ortogonales. Como veremos a continuación, de acuerdo con el método de Polack, si conocemos la base de vectores ortogonales, podríamos encontrar el grupo de vectores asociados s mutuamente conjugantes.
53 gi +1 = gi − τ i Asi con τi =
gTi gi gTi Asi
si +1 = gi +1 + κ i si con κ i = −
1< i ≤ N
sTi Agi +1 sTi Asi
notemos que:
κi = −
siT Ag
i +1
T i
s Asi
=−
g T Asi i +1 T i
s Asi
=
(
gT g i +1
i +1 T i
)
− g /τi i
s Asi
=
gT g i +1 T i i
i +1
τ s Asi
pero:
τ s Asi = T i i
gT g i T
i
g Asi
s Asi = T i
gT g i
( s − κ s ) As i
i
siT Asi = g T g
i
i i −1
i
i
i
entonces:
si +1 = gi +1 +
gTi +1 gi +1 gTi gi
si
(3.4.7)
s1 = g1 La ecuación (3.4.7) establece que cuando conocemos completamente un conjunto de vectores ortogonales, sin necesidad de conocer la matriz A , podemos obtener un conjunto de direcciones asociadas que son mutuamente conjugantes de tal matriz. Gradiente conjugado (Fletcher and Reeves, 19642 )
Al utilizar un método de minimización direccional, se sabe que la primera dirección atractiva será s0 = −∇f ( x0 ) , esta dirección nos permitirá obtener un mínimo direccional x1 de modo que: s0 • ∇f ( x1 ) = −∇f ( x0 ) • ∇f ( x1 ) = 0
(3.4.8)
es decir, tenemos una secuencia de direcciones mutuamente ortogonales. La fórmula de minimización de Fletcher and Reeves, deducida de la de Polack, en función a que las direcciones negativas del gradiente son mutuamente ortogonales, nos permitirá construir una dirección conjugada, asociada al gradiente:
2
Fletcher, R. and Reeves, C., 1964. Function Minimization by Conjugate Gradients. Comput. J. , 7:149
54 s1 = −∇f ( x1 ) +
∇T f ( x1 )∇f ( x1 ) s0 ∇ T f ( x 0 ) ∇f ( x 0 )
(3.4.9)
que para etapas superiores se generaliza como: ∇T f ( x k +1 )∇f ( x k +1 ) sk +1 = −∇f ( x k +1 ) + sk ∇ T f ( x k ) ∇f ( x k )
(3.4.10)
De acuerdo con las fórmulas, se aprecia que sólo se necesita información de gradiente, pero no de la Hessiana. Si la función es exactamente cuadrática y de orden n, la función será minimizada exactamente en n etapas Fletcher-Reeves, puesto que la Hessiana en estos casos, es una matriz de elementos constantes. Cuando la función no es cuadrática, se requieren más iteraciones. A diferencia del método de Powell, el método de Reeves genera una etapa iterativa con una sola dirección, de modo que cabe esperar su llegada al mínimo con menos evaluaciones de la función objetivo. Esta afirmación, sin embargo, es válida para problemas bien condicionados, puesto que los métodos basados en el gradiente resultan muy sensitivos a las características geométricas de las funciones. Debido a que suponemos que ahora se cuenta con información del gradiente, la minimización direccional puede realizarse de otras maneras:
x *k +1 = x *k + qk* s k
∂x ∂f df =∑ i∑ = sT • ∇f ( x ) = 0 dq q i ∂q i ∂xi
(3.4.11)
*
de acuerdo con esta, bastará encontrar un gradiente ortogonal a s, pero el gradiente ahora es explícito. Ejemplo considrere la función de Rosenbrok en el punto xT0 = −1. 2 1
f ( x ) = 100( x2 − x12 ) + ( 1 − x1 ) 2
2
Desarrolle dos iteraciones del método de Fletcher-Reeves. Vemos que:
⎡ −400( x2 − x12 ) x1 − 2( 1 − x1 ) ⎤ ∇f ( x ) = ⎢ ⎥ ; f ( x 0 ) = 24.2 2 ⎣ 200( x2 − x1 ) ⎦ entonces,
55
⎡ −215.6 ⎤ s0 = −∇f ( x 0 ) = −⎢ ⎣ −88.0 ⎥⎦ generalmente ocurrirá que las direcciones son bastante variables en magnitud, de modo que una buena idea es producir una normalización de las direcciones:
⎡ 0.925848 ⎤ s'0 = ⎢ ⎣ 0.377897⎥⎦ así tenemos que:
⎡ −1.03011⎤ x *1 = ⎢ ; f ( x *1 ) = 4.1281 ⎣ 1.06934 ⎥⎦ a partir de este mínimo direccional, podemos calcular el gradiente:
⎡ −0.67167⎤ ∇f ( x *1 ) = ⎢ ⎣ 1.64475 ⎥⎦ y la dirección conjugada será:
s1 = −∇f ( x1* ) +
∇ f ( x ) ∇f ( x
) ⎡ −0.67167⎤ s = − ⎢⎣ 1.64475 ⎥⎦ + ∇ T f ( x0 ) ∇f ( x0 ) 0
⎡ 0.65912⎤ =⎢ ⎣ −1.64987⎥⎦
T
* 1
* 1
⎡ −0.67167⎤ 1.64475]⎢ ⎣ 1.64475 ⎥⎦ ⎡ −215.6 ⎤ ⎢⎣ −88 ⎥⎦ ⎡ −215.6 ⎤ [ −215.6 −88]⎢⎣ −88⎥⎦
[ −0.67167
56
it
x1
x2
f ( x)
0 1 2 3 4 5 10 15 20 25 30 35 40 45 50 55 60 61 62 63
-1.2000 -1.0301 -0.9295 -0.7850 -0.6795 -0.5941 -0.3012 -0.1028 0.0634 0.2249 0.4017 0.6146 0.8978 1.3645 1.4954 1.0128 1.0001 1.0001 1.0001 1.0000
1.0000 1.0693 0.8282 0.5574 0.3895 0.2713 -0.0162 -0.1089 -0.1234 -0.0820 0.0263 0.2444 0.6841 1.7964 2.2397 1.0267 1.0002 1.0002 1.0001 1.0000
24.2000 4.1281 3.8514 3.5319 3.3421 3.2069 2.8367 2.6440 2.4998 2.3577 2.1830 1.9288 1.4982 0.5619 0.2466 0.0003 0.0000 0.0000 0.0000 0.0000
s 2.16E+02 6.84E-01 9.23E+01 1.88E+02 2.45E+02 2.81E+02 3.40E+02 3.53E+02 3.80E+02 4.37E+02 5.27E+02 6.37E+02 6.91E+02 2.71E+02 -2.16E+00 7.04E-02 -2.70E-04 -5.00E-05 -2.96E-02 3.00E-05
s/ s 8.80E+01 -1.64E+00 -1.73E+02 -2.99E+02 -3.39E+02 -3.41E+02 -2.13E+02 -7.91E+01 4.58E+01 2.02E+02 4.43E+02 8.29E+02 1.34E+03 8.68E+02 -5.70E+00 -7.17E-01 1.00E-05 -1.10E-04 -5.65E-02 -4.00E-05
0.9258 0.3849 0.4709 0.5320 0.5858 0.6358 0.8468 0.9758 0.9928 0.9081 0.7658 0.6089 0.4591 0.2976 -0.3549 0.0977 -0.9993 -0.4138 -0.4641 0.6000
Trayectorias de minimización
gradiente conjugado gradiente directo
2
x2
1
0
-1
0
1
x1
2
0.3779 -0.9229 -0.8822 -0.8467 -0.8104 -0.7719 -0.5319 -0.2187 0.1196 0.4188 0.6431 0.7932 0.8884 0.9547 -0.9349 -0.9952 0.0370 -0.9104 -0.8858 -0.8000
57
3.5 Técnicas indirectas de segundo orden
En ciertas situaciones, existe la posibilidad de evaluar información diferencial de primer y segundo orden para la función objetivo en minimización. Si consideramos una aproximación cuadrática a una función, tendremos:
1 f ( x ) = f ( x 0 ) + ∇ T f ( x 0 ) ∆ x + ∆ x T H∆ x 2
(3.5.1)
luego, tomando el gradiente de la aproximación, obtenemos:
∇f ( x ) = ∇f ( x 0 ) + H∆ x = 0
(3.5.2)
puesto que la dirección ∆x representa una dirección de minimización, ie: x = x0 + qs ⇒ ∆ x = q s
la forma recursiva de obtener una dirección de minimización será:
s k +1 = − H −1 ( x k ) ∇f ( x k )
(3.5.3)
en funciones objetivo exactamente cuadráticas, la ecuación (3.5.3) apunta directamente al centro del sistema coordenado que centra a las cónicas, siendo posible obtener el mínimo en una sola estapa de minimización. En formas no cuadráticas, (3.5.3) permite generar una secuencia de direcciones que otorga a la convergencia una característica tipo Newton, o de segundo orden. Obviamente esta técnica constituye el ideal de la minimización, pero es necesario tener en cuenta que la obtención de la Hessiana analítica no siempre es posible. Al contar con toda la información diferencial posible, existe la alternativa de revisar la característica de la Hessiana, de modo de evitar los puntos de silla asociados a Hessianas indefinidas.
58 3.5 Técnicas indirectas de segundo orden En algunos problemas de optimización, existe la posibilidad de evaluar información diferencial de primer y segundo orden para la función objetivo. Si se considera una aproximación cuadrática a una función, tendremos:
1 f ( x ) = f ( x 0 ) + ∆ x T ∇f ( x 0 ) + ∆ x T H ( x 0 ) ∆ x 2
(3.5.1)
luego, tomando el gradiente de la aproximación, y anulándolo para conseguir una condición necesaria de punto estacionario, obtenemos:
∇f ( x ) = ∇f ( x 0 ) + H∆ x = 0
(3.5.2)
en la ecuación (3.5.2), el vector ∇f ( x 0 ) es un vector de constantes. Si la Hessiana es una matriz de elementos constantes, como es el caso de las formas cuadráticas, entonces el problema tiene solución inmediata, dada por: ∆ x = x − x 0 = − H −1 ∇ f ( x 0 )
(3.5.3)
Ejemplo:
considere la cónica
f ( x1 , x2 ) = 6 x12 − x1 x2 + 8 x22 − 3 en el punto xT0 = 1 1 . Determine el mínimo global del sistema.
⎡ 12 x1 − x2 ⎤ ⎡ 11⎤ = ∇f ( x 0 ) = ⎢ ⎣ − x1 + 16 x2 ⎥⎦ ⎢⎣15 ⎥⎦
⎡ 12 −1⎤ λ1 = 11.76 H (x) = ⎢ ; ⎣ −1 16 ⎥⎦ λ2 = 16.24
luego:
⎡ 12 −1⎤⎡ ∆x1 ⎤ ⎡ 11⎤ ⎢⎣ −1 16 ⎥⎦⎢⎣ ∆x ⎥⎦ = −⎢⎣15 ⎥⎦ 2 como solución del sistema, se obtiene:
⎡ −1⎤ ⎡1⎤ ⎡0⎤ ∆ x = ⎢ ⎥ = x − x0 = x − ⎢ ⎥ ⇒ x = ⎢ ⎥ ⎣ −1⎦ ⎣1⎦ ⎣0⎦ como la Hessiana es estrictamente definida positiva en todo ℜ, el mínimo es absoluto. Notemos que el vector ∆x es una dirección de minimización, a partir de un punto inicial.
59 Este simple problema, nos permite advertir todas las características de un método de segundo orden: • • •
al contar con la Hessiana, es posible hacer un análisis geométrico local de la función. en funciones cuadráticas, el método converge en una sola etapa. recordando que en los métodos direccionales:
x = x0 + qs ⇒ ∆ x = s
=
método cuadrático
− H −1∇f ( x0 ) (q = 1) (3.5.4)
este método no solo entrega la mejor de todas las posibles direcciones, sino que también establece cuánto debe avanzarse por esa dirección. En este sentido es posible, en principio, eliminar las minimizaciones unidimensionales en una dirección. Tratamiento de funciones no cuadráticas
En funciones no cuadráticas, que corresponden a la generalidad de los problemas de optimización, la Hessiana deja de ser una matriz de elementos constantes, sin embargo, el problema de minimización puede resolverse en forma recursiva, utilizando la fómula generatriz (3.5.4): sk = ∆ x k = x k +1 − x k = − H −1 ( x k )∇f ( x k )
(3.5.5)
recordando que esta última fórmula proviene de una aproximación cuadrática, su desempeño será bueno cuando la función es aproximable cuadráticamente en k+1 y k. En rigor, x k+1 no necesariamente es el mínimo direccional en la dirección sk = H −1 ( x k )∇f ( x k ), a menos de que la función original sea cuadrática. Esta situación puede apreciarse en la figura (3.5.1). Obviamente el hecho de no ir minimizando exactamente puede provocar un daño a la convergencia cuadrática, fundamentalmente en aquellos sectores cercanos a mínimos muy dependientes de las variables, donde la anulación del gradiente provoca poco avance; o bien en aquellos casos en que la minimización por aproximación cuadrática entra en regiones indefinidas de la superficie. Las desventajas del método de Newton en funciones no cuadráticas son las siguientes: •
la determinación de la Hessiana y el gradiente en forma anaítica es compleja en la generalidad de los problemas de optimización.
•
la convergencia se ve afectada cuando la aproximación cuadrática local no es adecuada.
•
la determinación de direcciones sucesivas requiere de la inversión de la Hessiana, o bien, la resolución de un sistema lineal, lo que obliga a revisar la condición del problema. el método tendrá problemas de convergencia en regiones indefinidas de la Hessiana, pudiendo conducir a indeterminación o convergencia a puntos de silla.
•
)
0 +qs
f (x
+q t s f(x) =f( x)
0
0
f(q)
0
∇f( x)
+1/ 2
q 2s t
0
H(x )
s
60
0.5
1.0
1.5
q
en todo caso, cabe mencionar que al disponer de información del gradiente, y al verificarse que la función objetivo es continua; hay una gran probabilidad de que un método de segundo orden pueda adecuarse a las dificultades del problema, y constituirse en la mejor alternativa de optimización. Evitado de regiones indefinidas de la superficie
Como ya sabemos, una Hessiana es representativa de un problema convexo (o cóncavo) cuando todos sus valores propios tienen el mismo signo. En el problema de minimización, nos interesa que los valores propios sean positivos. Levenberg y Marquardt proponen el siguiente procedimiento:
61
s k = −( H ( x k ) + βI ) ∇f ( x k ) −1
(3.5.6)
donde β es una constante positiva y no nula cada vez que la Hessiana es indefinida. La forma de evitar la combinación de valores propios, es elegir
β = má x( − λk ,0) + 1
(3.5.7)
k =1,n
cuando al menos un valor propio es negativo, β deja de ser no nulo. Ahora bien, si el valor propio es demasiado negativo, entonces:
s k = −( βI ) ∇f ( x k ) = − β∇f ( x k ) −1
(3.5.8)
es decir, el método tiende a lograr una convergencia escalonada del método del gradiente. Aún así, esta metodología logra que el método pueda avanzar hacia mínimos en regiones indefinidas, esperando una mejor definición de la Hessiana en el avance del problema. Ejemplo:
Minimizar : f ( x ) = 100( x2 − x12 )2 + (1 − x1 )2 partiendo de xT0 = 0. 3 0. 4 Desarrollando los cálculos correspondientes, se tiene:
⎡ −2 + 2 x1 + 400 x13 − 400 x1 x2 ⎤ ⎡ −38.6 ⎤ ∇f ( x ) = ⎢ f ( x ) ∇ = ⎥ 2 0 ⎢⎣ 62.0 ⎥⎦ ⎣ 200 ( x2 − x1 ) ⎦ ⎡ −50 −120 ⎤ ⎡ 2 + 1200 x12 − 400 x2 −400 x1 ⎤ ( ) H x = H (x ) = ⎢ 0 ⎢⎣ −120 200 ⎥⎦ 200 ⎥⎦ ⎣ −400 x1 de acuerdo con la fórmula direccional de segundo orden, podríamos encontrar la primera dirección resolviendo:
⎡ −11475 ⎡ 38.6 ⎤ ⎡ −50 −120 ⎤ . × 10 −4 ⎤ H0 s0 = −∇f ( x 0 ) ⇒ ⎢ s = ⇒ s0 = ⎢ −0.3169 ⎥⎦ ⎣ −120 200 ⎥⎦ 0 ⎢⎣ −62.0 ⎥⎦ ⎣ a continuación, encontramos el próximo vector de:
⎡ 0.2885 ⎤ x 1 = x 0 + s0 = ⎢ ⎣ 0.0833 ⎥⎦ siguiendo este mismo esquema en las siguientes etapas, se obtendrá:
62
it
x1
x2
0 1 2 3 4
0.3000 0.2885 0.9817 0.9819 0.9999
0.4000 0.0831 0.4832 0.9641 0.9999
sT -0.0115 0.6932 0.0002 0.0181 0.0001
f (x) -0.3169 0.4001 0.4809 0.0355 0.0001
10.1000 0.5062 23.0917 0.0003 0.0000
notar que la convergencia tiene un comportamiento oscilatorio e irregular. Haciendo un cálculo de valores propios de la Hessiana asociada al punto de partida, se obtiene: λ1 = 248. 3 λ 2 = −98. 28
de donde se deduce un comportamineto tipo punto de silla en la superficie. Si seguimos el algoritmo de Marquardt, tendremos:
⎡ 49.28 −120 ⎤ H$ k = Hk + má x{− λk ,0} + 1 ∏ ⇒ H$ 0 = ⎢ ⎣ −120 299.28 ⎥⎦
(
)
los valores propios asociados a esta nueva matriz serán positivos. En el método de Marquardt, la generatriz de direcciones es equivalente a la de Newton, pero con una Hessiana forzada convexa:
⎡ 49.28 −120 ⎤ ⎡ 38.6 ⎤ ⎡ 11.8 ⎤ = ⇒ = H0 s0 = −∇f ( x 0 ) ⇒ ⎢ s s ⎢⎣ 4.52⎥⎦ 0 ⎣ −120 299.28 ⎥⎦ 0 ⎢⎣ −62.0 ⎥⎦ esta dirección de partida debe ser tratada con una minimización direccional, pues no corresponde a una dirección original del Newton. Consideramos entonces una minimización de la función en la parametrización:
⎡ 0.760621⎤ x 1 = x 0 + qs0 ⇒ x 1 = ⎢ ⎣ 0.576441⎥⎦ avanzando con el algoritmo, se obtiene: it
x1
x2
0 1 2 3 4
0.3000 0.7606 0.9291 0.9397 0.9987
0.4000 0.5764 0.8349 0.8830 0.9939
sT 11.8000 0.1685 0.0106 0.0589 0.0008
4.5200 0.2584 0.0481 0.1109 0.0050
λ1
λ2
f (x)
248.30 0.85 2.96 0.45 0.68
-98.28 664.82 901.10 908.07 1000.58
10.1000 0.0577 0.0856 0.0036 0.0012
63 Comparando las dos tablas, vemos que el método de Marquardt no reduce las iteraciones, pero si estabiliza la convergencia del método, situación deseable en problemas mal condicionados, con la finalidad de evitar oscilaciones. Las iteraciones 1 y siguientes corresponden a regiones bien definidas de la Hessiana. El cambio observado en los valores propios nos da además una medida de la calidad de la aproximación cuadrática. Claramente el método logra grandes progresos cuando estos valores no cambian mucho, lo que indica que la aproximación cuadrática es buena. Evitado de grandes diferencias en mínimos por aproximación cuadrática y proyección funcional
De acuerdo a su fórmula genérica, el método de Newton adopta la dirección: sk = − H −1 ( x k )∇f ( x k )
(3.5.9)
en rigor, debe minimizarse de forma que: x k +1 = x k + q sk
(3.5.10)
para un método de Newton convencional, q adquiere el valor unitario. De esta forma, es posible evaluar la función en dos puntos de la proyección: q = 0 ⇒ f (x) = f (xk ) q =1
⇒ f ( x ) = f ( x k + sk )
además de esta información, se sabe que: df df = ∇T f ( x ) sk ⇒ = ∇T f ( x0 ) sk dq dq q = 0 con estos tres datos, es posible evaluar una función cuadrática unidimensional en q: f ( q ) = a + bq + cq 2 y minimizar en términos de q:
q
ópt
∇T f ( x k ) sk b =− =− 2a 2 f ( x k + sk ) − f ( x k ) − ∇T f ( x k ) sk
[
]
(3.5.11)
este procedimiento puede repetirse hasta alcanzar una precisión deseable en el mínimo direccional. Obviamente a través de (3.5.11), q será distinto al valor unitario si la función no es exactamente cuadrática. Esta última también puede utilizarse para decidir si el paso de Newton es adecuado (q ópt → 1), pues su construcción es sencilla y depende de información secuencial generada en el proceso de optimización.
64
Estimando la Hessiana, a partir de información diferencial de primer orden: variable métrica n-dimensional (o método de la secante)
Capítulo 4 Programación Lineal (LP) Un número importante de problemas de optimización conduce a funciones objetivo y restricciones lineales. A modo de ejemplo: •
La planificación de producción de una empresa que puede adquirir materias primas de diferente calidad y rendimiento, para dar origen a una variedad de productos de diferente costo, cuyo rendimiento depende de las materias primas, con la finalidad de aumentar el beneficio económico (problema de mezclado).
•
La asignación diaria de actividades a un grupo de empleados, de manera de maximizar su eficiencia productiva.
•
El manejo de la distribución de productos o adquisición de materias primas.
•
etc.
Dentro de los rangos de flexibilidad de operación de una planta química establecida, los balances de materia y energía comúnmente pueden ser linealizados o linealizados por tramos, estableciendo rangos de operación de funcionalidad lineal en las variables de optimización. El problema de optimización de funciones lineales, sometidas a restricciones también lineales, da origen al problema de programación lineal. El término programación lineal no está específicamente referido a los códigos de programación del problema, sino que a que pueden ser tratados con técnicas de resolución y reducción de sistemas lineales ordinarios. De modo de ir conociendo las características y elementos de un problema LP, vamos a considerar en detalle el desarrollo de un problema específico: Una refinería puede adquirir dos calidades de crudo para su operación: crudo #1 y crudo #2. Sus productos de venta son la Gasolina, el Kerosene, el Fuel Oil y otros productos combustibles o asfálticos denominados residuos:
Crudo #1
Gasolina
US$36/bbl
Kerosene
US$24/bbl
Fuel Oil
US$21/bbl
Residuos
US$10/bbl
US$24/bbl Refinería
Crudo #2
US$15/bbl
1
por otro lado, se dispone del siguiente cuadro representativo de los rendimientos para cada crudo, la limitaciones operativas, y los costos de producción: Tabla de datos para el problema
Gasolina Kerosene Fuel Oil Residuos Costo de procesamiento (US$/bbl)
Rendimiento volumétrico (%) Crudo #1 Crudo #2 80 44 5 10 10 36 5 10 0.50
capacidad de producción (bbl/día) 24000 2000 6000
1.00
Se requiere establecer el flujo de materias primas y la distribución de productos, de modo de maximizar las utilidades. Formulación del problema de Optimización. 1. Establecer la nomenclatura para las variables de optimización:
x1 : x2 : x3 : x4 : x5 : x6 :
bbl/día de Crudo #1 bbl/día de Crudo #2 bbl/día de Gasolina. bbl/día de Kerosene. bbl/día de Fuel Oil. bbl/día de residuos.
2. Establecer la función objetivo y su dependencia en las variables de optimización:
⎫ ⎧ costo ⎫ ⎧ costo ⎬−⎨ ⎬ Má xf ( x ) = Má xUtilidad = ventas − ⎨ ⎩ materia prima ⎭ ⎩operativo ⎭ donde: ventas = 36 x 3 + 24 x 4 + 21x5 + 10 x 6 ⎧ ⎫ costo ⎨ ⎬ = 24 x1 + 15x 2 ⎩ materia prima ⎭ ⎧ costo ⎫ ⎨ ⎬ = 0.5x1 + x 2 ⎩operativo ⎭
vemos que la función objetivo depende de 6 variables de optimización.
2
3. Determinar las restricciones de igualdad: del cuadro de rendimientos, se obtiene un balance en volumen para cada crudo: x3 = 0.80 x1 + 0. 44 x2 x4 = 0. 05x1 + 0.10 x2 x5 = 0.10 x1 + 0. 36 x2 x6 = 0. 05x1 + 0.10 x2 haciendo un análisis de grados de libertad, tenemos aquí 4 restricciones de igualdad que relacionan 6 variables. Obviamente, la solución satisfactoria, debe cumplir exactamente todas las relaciones de igualdad. Como habíamos señalado, las restricciones de igualdad permiten reducir las dimensiones del problema de optimización.
6 variables − 4 ecuaciones ⇒ 2 grados de libertad 4. Determinamos ahora las restricciones de desigualdad. En la tabla se aprecia que hay capacidades productivas máximas:
x3 ≤ 24000 x4 ≤ 2000 x5 ≤ 6000 por otro lado, puesto que cada variable representa físicamente a un flujo, estos deben estar restringidos a ser positivos, es decir,
xi ≥ 0 i = 1,6 si además hubieran restricciones en cuanto a la disponibilidad de materias primas, también sería necesario incluirlas como restricciones. Siguiendo este problema, vemos que un problema (LP) tiene una forma general: r
Mín f (x) = ∑ ci xi = c x T
r: número de variables
i =1
sujeto a: r
∑a
x = b j ⇒ Ax = b
ji i
i =1 r
∑ a' i =1
ki
xi ≥ b 'k ⇒ A ' x ≥ b '
j = 1,m m: número de restricciones de igualdad. k = 1,p p: número de restricciones de desigualdad
la dimensión del problema de optimización será: r - m. Siempre es conveniente el reemplazo de las restricciones de igualdad, de modo de reducir las dimensiones del problema de optimización. En un caso general de funciones lineales el objetivo señalado es 3
siempre posible y sencillo, pero no lo es para funciones generales no lineales, como por ejemplo:
Mín: f ( x) = x 1 + x 2 Sa: 0 = x 1 exp( x 2 ) − 4 x 2 donde x2 no puede despejarse en términos de x1 de la restricción de igualdad. En el caso que venimos siguiendo de la refinería, es posible reducir el problema a 2 variables. Análogamente, es posible convertir un problema de maximización en uno equivalente de minimización siguiendo la transformación:
Má x U( x) ⇔ Mín f ( x) = − U( x) al reemplazar las restricciones de igualdad en la función objetivo y en las restricciones de desigualdad, se obtiene:
Máx U( x) ⇔ Mín f ( x) = − U( x) = −8.1x1 − 10.8x2 restringido a: A. Producción máxima de gasolina:
0.80x1 + 0. 44 x2 ≤ 24000 B. Producción máxima de Kerosene
0. 05x1 + 0.10x2 ≤ 2000 C. Producción máxima de fuel oil
0.10x1 + 0. 36x2 ≤ 6000 D. Flujos positivos de crudo, o condición de planta operativa x1 ≥ 0 x2 ≥ 0 y observando las restricciones de igualdad, es posibles constatar que al elegir flujos positivos de crudo, también se obtendrán flujos positivos de todos los productos. Consideremos ahora el problema en un plano de curvas de nivel. Notamos que la función objetivo es lineal:
4
Mín f ( x) = −8.1x1 − 10.8x2 por tanto:
. ⎤ ⎡ −81 ∇f ( x) = ⎢ ⎥ ⎣ −10.8⎦
;
⎡ 0 0⎤ H ( x) = ⎢ ⎥ ⎣ 0 0⎦
claramente, la función por sí misma no tiene solución finita o acotada. El gradiente negativo indica que la función es minimizable cuando x1 y x2 crecen hasta el infinito. El problema LP no puede ser resuelto por técnicas convencionales de optimización cuadrática, y sin la presencia de restricciones, su solución siempre será infinita.
Figura #1 Curvas de nivel de una función objetivo lineal 6e+4
4e+4
x2 (bbl/día crudo #2)
f= f=
2e+4
f=
f=
-2
00
−∇f(x)
00 0
00 00
-1 00
f=
0e+0
-3
0
00 0
0.0
10
00
00
-2e+4 -2e+4
0e+0
2e+4
4e+4
6e+4
x1 (bbl/día crudo #1)
Obviamente, la situación de que las variables de optimización deban cumplir con un conjunto de restricciones, genera un dominio factible de búsqueda de variables. Al incluir las restricciones en el problema anterior, se obtiene la figura #2, en la que se aprecia un dominio factible para solucionar el problema. En la Figura indicada aparecen los siguientes hechos notables: • • • •
la inclusión de un dominio restricto produce un problema de solución finita. la solución aparecerá en algún vértice del dominio factible. la posición de la solución en el dominio factible depende de la pendiente de las curvas de nivel los vértices de la forma geométrica correspondiente al dominio factible corresponden a intersecciones entre las restricciones.
5
Esta última observación es la más importante, pues define la estrategia de búsqueda de soluciones para el problema de la programación lineal: basta analizar el monto funcional en las intersecciones entre restricciones. En el caso específico que se está tratando, la solución del problema ocurre en la intersección de las restricciones A y B, que puede ser obtenida directamente solucionando el correspondiente problema lineal (A ): 0.80x 1 + 0.44 x 2 ≤ 24000 x 1 = 26207 ⇒ f mín = −286764 ( B): 0.05x 1 + 010 . x 2 ≤ 2000 x 2 = 6897 La solución gráfica de un problema de programación lineal es posible cuando el problema depende, globalmente, de dos variables de optimización. Sin embargo, la generalidad de los casos depende de un número considerable de variables, de modo que a continuación presentaremos una una estrategia que permite resolver un problema de cualquier dimensión. Figura #2 Curvas de nivel de una función objetivo lineal 6e+4 A : 0.80 x1 + 0.44 x2 ≤ 24000 B : 0.05 x1 + 0.10 x2 ≤ 2000 C : 0.10 x1 + 0.36 x2 ≤ 6000
x2 (bbl/día crudo #2)
4e+4
(A) (B) 2e+4
dominio factible
(C)
0e+0
-2e+4 -2e+4
0e+0
2e+4
4e+4
6e+4
x1 (bbl/día crudo #1)
4.1 Clasificación de las soluciones en términos de las restricciones •
Problemas de múltiple solución en el LP
Dependiendo de las pendientes de las curvas de nivel de la función objetivo, el problema de programación lineal podría llegar a tener múltiples soluciones. A modo de ejemplo, en el problema:
6
Máx: f (x) = 2 x1 + 0.5x2 Sujeto a: 6x1 + 5x2 ≤ 30 (A) 4 x1 + x2 ≤ 12 (B) x1 , x2 ≥ 0 Notar que si el problema es de minimización, tendría solución única en el vértice ( 0, 0) . Es posible, de este modo, sospechar de problemas de múltiple solución, cuando la función objetivo es una combinación lineal de una de sus restricciones, aunque siempre es necesario desarrollar un análisis previo de la posición de la solución. •
Problemas sin solución acotada.
Estos problemas aparecen cuando las restricciones son insuficientes para producir el acotado de la solución dentro de un dominio factible, por ejemplo: Mín: f ( x) = − x1 − x2 Sujeto a: 3x1 − x2 ≥ 0 (A) x2 ≤ 3 (B) x1 , x2 ≥ 0 Figura #3 Problema de múltiple solución
Figura #4 Problema sin solución acotada
6
6
(B)
5
5 (A)
4
4
3
x2
x2
f decreciente 3 (B) f creciente 2
2
1
1
(A)
0
0 0
1
2
3
4
x1
5
6
0
1
2
3
4
5
6
x1
Gráficamente se observa que x2 permanecerá acotado en 3, pero x1 puede crecer sin restricción hasta el infinito. •
Restricciones incompatibles.
7
Dentro de la multiplicidad de problemas LP, también es posible encontrar aquellos que poseen restricciones incompatibles, no pudiendo generar regiones factibles de solución.
Figura #5 Problema con restricciones incompatibles
x1
x2
..... Las dos últimas categorías señaladas aquí, corresponden a problemas LP mal definidos. Cabe indicar que todo problema LP correctamente definido conduce a dominios factibles de tipo convexo, y a superficies definidas (no estrictamente cóncavas o convexas, sino más bien, cóncavas y convexas). 4.2 Determinación de restricciones de carácter lineal El objetivo de esta sección es facilitar el planteamiento de restricciones de carácter lineal, en problemas productivos relacionados con la Ingeniería Química. Cada vez que se plantea un problema de optimización han de tenerse en cuenta: •
que hay limitaciones productivas que se alcanzan por límites físicos de operación de los equipos, por la capacidad de almacenamiento, o bien, por condiciones restrictivas del mercado.
•
la disponibilidad de materias primas puede estar limitada, debido a problemas de distribución y transporte, almacenamiento, o bien, porque la planta depende de bienes intermediarios producidos por otras plantas.
•
existen restricciones de seguridad u operabilidad, por ejemplo, hay mezclas que son explosivas o contaminantes en determinados rangos de composición, los materiales sufren fatigas térmicas de acuerdo a su calidad. De la misma forma, la resistencia física de un equipo puede estar limitada a ciertos rangos específicos de presión.
8
•
los balances de materia y energía, en general, no son lineales, pero pueden ser considerados lineales dentro de ciertos rangos, donde es posible plantear una proximación afín del problema. Generalmente, el balance de materia puede escribirse en términos de rendimiento, que da origen a expresiones lineales.
•
algunas propiedades físicas de mezclas, pueden escribirse como una contribución ponderada linealmente en composición de todos sus constituyentes. Muchos productos puestos en el mercado requieren de ciertas especificaciones mínimas de calidad. Por ejemplo, los hidrocarburos combustibles forman soluciones aproximadamente ideales, la presión de vapor de una gasolina podría ser escrita como: Nc
P = ∑ xi Pisat v
i =1
4.3 Algoritmo Simplex para problemas LP En 1947, Dantzig diseñó una estrategia para resolver problemas de optimización con funciones objetivo y restricciones de carácter lineal, que consiste simplemente en analizar el valor de la función objetivo en los vértices del dominio factible. Suponemos que el problema de optimización está bien definido, es decir, no hay restricciones incompatibles o soluciones no acotadas, la metodología para ubicar tales inconsistencias será analizada más adelante. De igual forma, la función objetivo y las restricciones de desigualdad, ya contienen las restricciones de igualdad: r
f ( x ) = ∑ ci x i
Mín:
i =1
≥ Sujeto a: ∑ a ji xi ≤ b j ≥ 0 j = 1, p i =1 xi ≥ 0 i = 1, r r
(4.1)
en esta última ecuación, todos los b j deben ser positivos. El primer paso consiste en introducir tantas variables de holgura positivas como restricciones haya (p), observando lo siguiente: Caso 1
r
r
i =1
i =1
∑ a jixi ≥ b j ⇒ ∑ a jixi − sj = b j
(4.2)
si se resta un número positivo en el miembro izquierdo de (7.2), podemos hacer que la suma iguale al valor límite de la restricción. Caso 2
r
r
i =1
i =1
∑ a jixi ≤ b j ⇒ ∑ a jixi + sj = b j
(4.3)
9
si se suma un número positivo en el miembro izquierdo de (7.3), podemos hacer que la suma iguale al valor límite de la restricción. Las variables s j se denominan variables de holgura, y son de carácter artificial. Cuando s j = 0 , estamos exactamente en el contorno de la restricción. Si s j > 0 , estamos en un punto que cumple la restricción, pero este punto no se encuentra en la línea de la restricción, sino que en otra parte del dominio factible. La sustitución previamente indicada “pone” el problema en el cuadrante positivo del sistema coordenado. A modo de ejemplo, consideremos el problema: Mín: Sujeto a:
•
f ( x) = − x1 + x2 2 x1 − x2 ≥ −2 − x1 + 3x2 ≥ −2 − x1 − x2 ≥ −4 x1 , x2 ≥ 0
Etapa 0: Transformar todas las restricciones de desigualdad, de modo que los b j sean positivos.
En el ejemplo, se obtiene: Mín: f ( x) = − x1 + x2 Sujeto a: −2 x1 + x2 ≤ 2 (A) x1 − 3x2 ≤ 2 (B) x1 + x2 ≤ 4 (C) x1 , x2 ≥ 0 •
Problema original
Etapa 1: Introducir las variables de holgura a cada ecuación de restricción, siguiendo la regla (4.2) y (4.3)
En el ejemplo, todas las restricciones de desigualdad p, conducen a relaciones de minorante, de modo que a cada ecuación sumamos una variable de holgura: Mín: Sujeto a:
f ( x ) = − x1 + x2 −2 x1 + x2 + x3 = 2 (A) x1 − 3x2 + x4 = 2 (B) x1 + x2 + x5 = 4 (C) x1 , x2 , x3 , x4 , x5 ≥ 0
Problema aumentado
en este momento, las p restricciones de desigualdad originales, se han transformado en m = p restricciones de igualdad. El problema de optimización depende ahora de 5 variables, 10
relacionadas por 3 ecuaciones (el sistema lineal A, B, C), dando lugar a un problema de dos dimensiones, equivalente al original. Notar que este sistema de 3 ecuaciones en 5 variables, puede resolverse si se especifican dos variables positivas. •
Etapa 2: Encontrar una solución que se encuentre en uno de los vértices del dominio factible.
Si consideramos un análisis de los grados de libertad del nuevo problema de optimización, tenemos un número total de n = r + p variables, donde m = p son nuevas restricciones de igualdad, los grados de libertad del problema, son r. En el caso específico que venimos tratando: n = r+p = 2+3=5 m= p=3 como los grados de libertad son dos, es claro que la especificación de dos de las cinco variables resuelve el problema. En este sentido, tenemos dos variables independientes, que en programación lineal reciben el nombre de variables no básicas, y tres variables dependientes, denominadas variables básicas. Cualquier especificación de variables no básicas capaz de resolver el sistema de restricciones de igualdad, genera una solución básica del problema. Del conjunto de soluciones básicas, existe un subconjunto denominado de soluciones factibles, para las que tanto las variables básicas como no básicas son positivas, encontrándose dentro del dominio de optimización, específicamente en uno de los vértices de la figura geométrica del dominio. Retomando el problema original de optimización, podemos hacer el trazado mostrado en la figura #6. Son soluciones básicas factibles la intersección de (A) con (C), la intersección de (B) con (C), la intersecciones de (A) con el eje x2 y (B) con el eje x1 . Es solución básica, pero no factible, la intersección de (A) y (B). Notemos que una solución factible inicial ocurre cuando x1 = x2 = 0 . En este caso f ( x) = 0 x3 = 2 (A) x4 = 2 (B) x5 = 4 (C)
nótese que las variables de holgura (x3, x4, x5) también son positivas y no nulas por tanto la solución básica y factible no se encuentra en ninguna de las restricciones de igualdad (las variables de holgura son no nulas). Ahora debemos avanzar por los vértices del dominio factible, en busca del óptimo.
11
Figura #6 Función objetivo a minimizar y su dominio factible 6
(A)
x2
4
2
f decreciente
(C)
(B) 0 0
2
4
6
x1
Etapa 3: avanzar a una nueva solución factible, de modo de minimizar la función.
a. La primera solución factible, determinada en la etapa anterior es x1 = x2 = 0 −2 x 1 + x 2 x 1 − 3x 2 x1 + x 2 f ( x) + x 1
+ x3 + x4 + x5 − x2
=2 =2 =4 =0
(A ) ( B) ( C) (*)
en el problema vemos que el aumento de x1 favorece la reducción de la función objetivo, el aumento de x2, en cambio, produce un aumento de la función objetivo. La variable x1 es potencialmente atractiva para el objetivo de minimización. En esta etapa debemos encontrar un nuevo conjunto de variables no básicas (o independientes), al elegir x2 = 0 como variable no básica estamos seguros de que la función debe ser minimizada. La regla de selección establece que debe transformarse en básica aquella variable no básica con el mayor coeficiente positivo en (*), en este caso, x1 b. Cuando x1 crece deja de ser una variable no básica, pues debe crecer conforme a las restricciones, de modo que es necesario reemplazarla otra. Con x2 = 0, vemos que en la restricción (A) x1 puede crecer en forma ilimitada, sin hacer negativa a x3. Por la restricción (B), x1 puede crecer hasta 2, por la restricción (C) x1 puede crecer hasta 4. Como regla general, se incluye como nueva variable no básica a aquella cuya operatoria de coeficientes satisface:
12
⎛ bj ⎞ mín⎜ ⎟ ≥ 0 l: variable transformada en no bá sica en .a (x 1 ) ⎝ a jl ⎠ por este criterio se busca la restricción más desfavorable para el incremento de x1. Analizando las correspondientes, se concluye:
( A ) : 2 / −2 = −1 ( B) : 2 / 1 = 2 ( C) : 4 / 1 = 4 la menor positiva corresponde a la restricción (B). Recordando que x2 = 0, la variable activa con crecimiento de x1 en (B) es x4, y x1 puede crecer con x2 = 0 hasta el valor 2, en el que x4 = 0. De esta forma, las nuevas variables no básicas serán x2 = x4 = 0. Como la variable de holgura de (B) es nula, buscamos una solución en la restricción correspondiente. Luego, el criterio de búsqueda es: a. Escribir la función objetivo despejando el cero como sigue f ( x ) − ∑ ci x i = 0 i
eliminar como variable independiente o no básica a aquella que presente el mayor coeficiente positivo -ci. Tal variable es la que reducirá en mayor magnitud a la función objetivo. b. Introducir como nueva variable no básica a aquella que presente el menor valor positivo de: ⎛ bj ⎞ ⎜ ⎟ j = 1, m ⎝ a jl ⎠ que es representativa de la restricción más desfavorable "o restrictiva" para el crecimiento de la variable mayoritariamente minimizante Etapa 4: Transformar las ecuaciones de acuerdo a la selección de las nuevas variables no básicas.
Como resultado de la etapa 3, hemos seleccionado como variables independientes o no básicas a x2 = x4 = 0.
13
Recordando que el número de grados de libertad del problema de optimización es dos, la selección de dos variables no básicas permite obtener la solución del sistema de restricciones de igualdad. Ahora corresponde escribir las ecuaciones en términos de las nuevas variables no básicas: Variables básicas: x1, x3, x5 Variables no básicas: x2, x4 considerando las relaciones: −2 x 1 + x 2 x 1 − 3x 2 x1 + x 2 f ( x) + x 1
+ x3 + x4 + x5 − x2
=2 =2 =4 =0
(A ) ( B) ( C) (*)
podemos escribir cada una de las variables básicas en términos de variables no básicas. De la restricción más desfavorable se deduce:
( B) : x1 = 2 + 3x2 − x4 puesto que x1 deja de ser no básica, se reemplaza en el resto de las restricciones: (A ): −2( 2 + 3x 2 − x 4 ) + x 2 + x 3 = 2 ⇒ x 3 − 5x 2 + 2 x 4 = 6 x 1 − 3x 2 + x 4 = 2 ( B): (C): ( 2 + 3x 2 − x 4 ) + x 2 + x 5 = 4 ⇒ x 5 + 4x 2 − x 4 = 2 f + ( 2 + 3x 2 − x 4 ) − x 2 = 0 ⇒ f + 2 x 2 − x 4 = −2
notar que dado que las variables básicas son x2 = x4 = 0, hemos logrado que la función objetivo adquiera valor -2. Mirando la línea de la función objetivo, convendría incorporar x2 como una nueva variable básica, dejando a x4 como no básica e incorporando a x5 como la nueva variable no básica. Nótese también que lo único que hemos hecho aquí es eliminar las variables básicas en términos de variables no básicas. El mismo procedimiento puede hacerse en utilizando transformaciones matrizales elementales, lo que introduce la representación tabular del LP: sistema sobredeterminado original −2 x 1 + x 2 x 1 − 3x 2 x1 + x 2 f ( x) + x 1
14
+ x3 + x4 + x5 − x2
=2 =2 =4 =0
(A ) ( B) ( C) (*)
x3 x4 x5
x1
x2
x3
x4
x5
f
b
-2 [1] 1 1
1 -3 1 -1
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
2 2 4 0
$
← fila pivote
columna pivote
la columna pivote se elige determinando la variable con mayor coeficiente positivo en la línea de la función objetivo. La fila pivote se determina buscando la razón b/a positiva más pequeña. Necesitamos eliminar a x1 como variable no básica. Multiplicando la fila (2) por 2 y sumándola a la (1), eliminamos x1 de la primera restricción. Multiplicando la fila (2) por -1 y sumándola a la tercera fila, eliminamos x1 de la tercera restricción. Procediendo también con operaciones elementales, es posible eliminar a x1 de la fila de la función: este es un simple proceso de eliminación de Gauss:
x3 x1 x5
x1
x2
x3
x4
x5
f
b
0 1 0 0
-5 -3 4 2
1 0 0 0
2 1 -1 -1
0 0 1 0
0 0 0 1
6 2 2 -2
esta es la representación matrizal de: ( A ): x3 − 5x2 + 2 x4 = 6 ( B): x1 − 3x2 + x4 = 2 ( C): x5 + 4 x2 − x4 = 2 3x2 − x4 + f = −2 Etapa 5: proceder iterativamente sobre la matriz, hasta que no se logre más mejoramiento en los valores de la función objetivo.
x3 x1 x5
x1
x2
x3
x4
x5
f
b
0 1 0 0
-5 -3 [4] 2
1 0 0 0
2 1 -1 -1
0 0 1 0
0 0 0 1
6 2 2 -2
←
15
↑
indicación: siempre el punto pivote debe quedar como [1] después de las transformaciones elementales:
x3 x1 x2
x1
x2
x3
x4
x5
f
b
0 1 0 0
0 0 1 0
1 0 0 0
0.75 0.25 -0.25 -0.5
1.25 0.75 0.25 -0.5
0 0 0 1
8.5 3.5 0.5 -3
Notar que en la línea de la función objetivo los dos coeficientes alcanzan valores negativos, esto indica que no es posible obtener mejoramiento en la reducción de la función objetivo, por lo que tenemos un mínimo. Quedan como variables no básicas x4 = x5 = 0 . Como x4 es una variable de holgura asociada a la restricción (B) y x5 es una variable de holgura asociada a la restricción (C), la solución está en la intersección de las restricciones (B) y (C), tal como se aprecia en la figura. Además, en el punto solución:
x3 = 8.5 x1 = 3.5 x2 = 0.5 en otras palabras, la restricción (A) representada por la variable de holgura x3 se cumple, pero la solución está en el dominio factible, y no activando la restricción correspondiente.. 4.4 Forma canónica del problema de programación lineal. Según hemos visto, todos los problemas de programación lineal, una vez incluidas las variables de holgura para la aplicación del método Simplex, tienen la forma: Mín: Sujeto a:
f ( x ) = c1x1 + c2 x2 +.... + cn x n a11x1 + a12 x2 +.... + a1n x n = b1 a 21x1 + a 22 x2 +.... + a 2 n x n = b2 M a m1x1 + a m2 x2 +.... + a mn x n = b m xi ≥ 0 ; b j ≥ 0 (i = 1, n ; j = 1, m)
(4.4)
o en una forma más compacta:
Mín: Sujeto a:
f =c x Ax = b T
x≥0; b≥0 16
(4.5)
la dimensión del problema de optimización es n-m, que también es indicativa del número de variables independientes o no básicas que deben especificarse para resolver el sistema de m restricciones de igualdad. Estas variables no básicas toman el valor nulo para generar culaquier solución básica del LP, sometidas a la condición de que el resto de las variables sean positivas. El número de soluciones básicas (no necesariamente factibles) está dado por:
⎛n⎞ n! ψ=⎜ ⎟= ⎝ m⎠ m!( n − m) !
(4.6)
esta relación es representativa del subconjunto de soluciones factibles. Obviamente, al crecer el número de variables, también crecerá el número de soluciones básicas y el número de soluciones factibles. Es posible demostrar que el algoritmo Simplex requiere de entre m a 3 m pruebas de soluciones factibles para llegar al mínimo. 4.5 Obtención de soluciones básicas factibles. Muchas veces ocurre que la obtención de una primera solución factible no es tan evidente, el problema Simplex siempre debe partir de una solución básica factible. Un ejemplo de esta situación es el problema: Mín:
f = x1 + 2 x2
Sujeto a: 3x1 + 4 x2 ≥ 5 x1 + x2 ≤ 4 Siguiendo los pasos del Simplex, es decir, la introducción de variables de holgura, obtendremos: Mín: f = x1 + 2 x2 Sujeto a: 3x1 + 4 x2 − x3 = 5 (A) x1 + x2 + x4 = 4 (B) xi ≥ 0 notamos que la selección de x1 y x2 como variables no básicas (es decir, ambas nulas), viola en la restricción (A) la condición positiva para x3. Esta situación se debe a que la solución no básica indicada no se encuentra en el dominio factible. En estos casos, es necesario agregar m nuevas variables de holgura, generando un procedimiento denominado Fase I - Fase II. En términos simples, en la Fase I se encuentra una solución factible para el problema en términos de la positividad de sus variables básicas, y en una segunda etapa, se parte de la solución factible para obtener soluciones factibles minimizantes. El algoritmo es el siguiente:
17
1. Tomar el problema de optimización original, y transformar todas las restricciones de desigualdad de manera que los b j sean positivos. 2. Agregar las variables de holgura, de manera de transformar las restricciones de desigualdad en restricciones de igualdad de acuerdo a las reglas indicadas. 3. Aumentar el sistema de ecuaciones obtenido en (2), agregando m nuevas variables de holgura: a11x1 + a12 x2 +.. + a1n x n a 21x1 + a 22 x2 +.. + a 2 n xn M a m1x1 + a m2 x2 +.. + a mn x n
+ x n +1 + xn + 2 O
+ xn + m
= b1 = b2 M = bm
las variables básicas serán: xn + i = bi (i = 1, m) las variables no básicas serán: xi = 0 (i = 1, n) Notamos que el problema satisface las restricciones positivas para todas las variables, pues las variables xn+1 a xn+m pueden acomodarse de manera que así sea. 4. Es claro que se recupera el sistema original cuando las variables básicas xn +i = bi (i = 1, m) se anulan. De esta forma, se define la función objetivo correspondiente a la Fase I dada por:
w = xn +1 + xn + 2 +.... + xn+1 escribiendo estas variables no básicas en términos de variables básicas, se obtiene: m ⎛m ⎞ ⎛m ⎞ ⎛m ⎞ w + ⎜ ∑ ai1 ⎟ x1 + ⎜ ∑ ai 2 ⎟ x 2 +....+⎜ ∑ ain ⎟ x n = ∑ bi ⎝ i =1 ⎠ ⎝ i =1 ⎠ ⎝ i =1 ⎠ i =1
en la Fase I se minimiza w utilizando la metodología Simplex. w debe tener mínimo nulo, pues así se recupera el sistema aumentado original. Si w mínimo es no nulo, se tiene un problema mal definido de optimización, es decir, un conjunto de restricciones no compatibles. En el ejemplo, el problema aumentado es: Mín: f = x1 + 2 x2 Sujeto a: 3x1 + 4 x2 − x3 = 5 (A) x1 + x2 + x4 = 4 (B) xi ≥ 0
18
puesto que las variables básicas x1 = x2 = 0 conducen a una solución básica, pero no factible, el problema vuelve a aumentarse: Mín: f = x1 + 2 x2 Sujeto a: 3x1 + 4 x2 − x3 + x5 = 5 (A) x1 + x2 + x4 + x6 = 4 (B) xi ≥ 0 en la Fase 1, se minimiza la función w:
w + 4 x1 + 5x2 − x3 + x4 = 9 para la que tenemos la siguiente tabla inicial para la Fase I
x1
x2
x3
x4
x5
x6
3 1 4
4 1 5
-1 0 -1
0 1 1
1 0 0
0 1 0
x2
x3
x4
x5
x6
1 0 0
-0.25 0.25 0.25
0 1 1
0.25 -0.25 -1.25
0 1 0
iteración 1 x1 0.75 0.25 0.25
$
$
w 0 0 1
b 5 4 9
w 0 0 1
b 1.25 2.75 2.75
←
←
notar que x4 es variable básica, por tanto, pese a presentar el mayor coeficiente para w, no puede ser eliminada como variable no básica.! iteración 2
x1
x2
x3
x4
x5
x6
1 1 0
1 0 0
0 1 0
1 4 0
0 -1 -1
1 4 -1
w 0 0 1
b 4 11 0
en esta iteración, w presenta un mínimo nulo, lo que sugiere seleccionar x1 = x4 = 0 en el problema aumentado:
19
Mín: f = x1 + 2 x2 Sujeto a: 3x1 + 4 x2 − x3 = 5 (A) x1 + x2 + x4 = 4 (B) xi ≥ 0 para las variables no básicas detectadas, se obtiene: x2 = 4 ; x3 = 11. Ambas variables básicas son positivas, por lo que tenemos una solución factible. En la Fase II se minimiza f. f = x1 + 2( 4 − x1 − x 4 ) ⇒ f + x1 + 2 x 4 = 8
la última iteración de la Fase I ya contiene la tabla acomodada en términos de las variables no básicas para la minimización del problema:
x1
x2
x3
x4
1 1 1
1 0 0
0 1 0
1 4 2
f 0 0 1
b 4 11 8
f 0 0 1
b 1.25 2.75 2.50
b 1.67 2.33 1.67
←
↑
Iteración 1:
x1
x2
x3
x4
0.75 0.25 1
1 0 0
-0.25 0.25 -0.5
0 1 0
x1
x2
x3
x4
1 0
1.33 -0.33
-0.33 0.33
0 1
f 0 0
0
-0.67
-0.33
0
1
↑
Iteración 2:
x1 x4
Finalmente la solución es:
x1 = 1. 67 ; x2 = 0 ; x3 = 0 ; x4 = 2. 33 ; f = 1. 67
20
←
Solución gráfica del ejemplo 4
x2
3
2
(B)
1 (A) f(x) 0 0
1
2
3
4
5
x1
4.6 Programación lineal en problemas que pueden incluir soluciones negativas En toda la formulación del problema de la Programación lineal, hemos indicado que las varibles de optimización, como las de holgura, están restrictas a ser positivas. De esta forma el problema pueda ser resuelto en el primer cuadrante del sisema coordenado cartesiano. Sin embargo, existen problemas específicos en los que una variable de optimización podría adquirir valores negativos, por ejemplo, una temperatura en escala relativa de un problema de refrigeración. La variable no restricta puede ser escrita como una combinación lineal de dos variables positivas, técnica que nos permite obtener el problema LP original, con una variable extra. A modo de ejemplo, considérese el problema: Mín: f ( x) = x1 + 4 x2 Sujeto a: x1 + x2 ≤ 3 (A) − x1 + x2 ≤ 1 (B) x2 ≥ 0 en el caso indicado, no se ha puesto restricción para x1 , de modo que puede tener valores positivos o negativos. En este caso:
21
x1 = x3 − x4 x3 ≥ 0 ; x 4 ≥ 0 el problema aumentado será: Mín: Sujeto a:
f ( x ) = x3 − x 4 + 4 x2 x3 − x4 + x2 + x5 = 3 (A) − x3 + x4 + x2 + x6 = 1 (B) x2 , x3 , x4 , x5 , x6 ≥ 0
tenemos un total de 5 variables relacionadas por dos ecuaciones, por tanto, una solución básica factible se obtiene al seleccionar 3 variables no básicas. Por ejemplo:
x2 = x3 = x4 = 0 x2
x3
x4
x5
x6
f
b
1 1 -4
1 -1 -1
-1 1 1
1 0 0
0 1 0
0 0 1
3 1 0
←
↑
Iteración 1
x5 x4
luego:
x2
x3
x4
x5
x6
f
b
2 1 -5
0 -1 0
0 1 0
1 0 0
1 1 -1
0 0 1
4 1 -1
x1 = x3 − x4 = 0 − 1 = −1 x2 = 0 4.7 Sensitividad de la solución del problema de Programación lineal
Una de las grandes ventajas del problema de programación lineal, es que la solución óptima puede ser estudiada con gran facilidad, en la tabla de solución, frente a pequeñas perturbaciones en las restricciones. A modo de ejemplo, consideremos la tabla solución para el problema de la refinería:
22
Mín: f ( x) = −8.1x1 − 10.8x2 Sujeto a: 0.80x1 + 0. 44 x2 ≤ 24000 Producción de Gasolina 0. 05x1 + 0.10x2 ≤ 2000 Producción de Kerosene 0.10x1 + 0. 36x2 ≤ 6000 Producción de FuelOil x1 , x2 > 0 cuyo problema aumentado será: Mín: f ( x) = −8.1x1 − 10.8x2 Sujeto a: 0.80x1 + 0. 44 x2 + x3 = 24000 Producción de Gasolina 0. 05x1 + 0.10x2 + x4 = 2000 Producción de Kerosene 0.10x1 + 0. 36x2 + x5 = 6000 Producción de FuelOil x≥0 la tabla de iteraciones convergente para este problema es:
x5 x1 x2
x1
x2
x3
x4
x5
0 1 0
0 0 1
0.14 1.72 -0.86
-4.21 -7.59 13.79
0
0
-4.66
-87.52
1 0 0
f 0 0 0
b 896.55 26206.90 6896.55
0
1
-286758.5
en términos de las variables no básicas solución, el problema puede escribirse como:
f − 4. 66x3 − 87.52 x4 = −286759 en la que claramente puede observarse que la función depende de dos variables de holgura, y por tanto está en las rectas de dos restricciones (producción de Gasolina y producción de Kerosene). Una pregunta razonable es qué sucedería si el límite productivo de gasolina puede ser aumentado. Mirando la formulación del problema, y recordando que las variables no básicas en la solución valen cero, un aumento en la restricción, es compatible con una disminución de la variable de holgura, por tanto la solución decrece y se hace más negativa. Por otro lado, vemos que:
∂f = 4. 66 ∂x3 ∂f = 87. 52 ∂x4 de donde se deduce que la disminución unitaria de x3 (aumento de la producción de gasolina), disminuye el mínimo en -4.66, y la diminución unitaria de x4 (aumento de la 23
producción de Kerosene) disminuye el mínimo en -87.52 , es decir, el mínimo es sensitivo al aumento de la producción de Kerosene fundamentalmente. Por otro lado, se aprecia que la función solución no es sensitiva a la producción de fuel oil.
24
