Una introduccion al metodo del escalado af´in para programacion lineal

´ AL METODO ´ UNA INTRODUCCION DEL ESCALADO ´ LINEAL AF´IN PARA PROGRAMACION Jordi Castro Statistics & Operations Research Dept. UPC DATE 2/2000 DR 2000/02

Una introducci´ on al m´ etodo del escalado af´ın para programaci´ on lineal Jordi Castro Dept. d’Estad´ıstica i Investigaci´o Operativa Universitat Polit`ecnica de Catalunya Pau Gargallo 5, 08028 Barcelona Abstract: El m´etodo del escalado af´ın fue el primero de los algoritmos denominados de punto interior para la soluci´on de problemas de programaci´on lineal. A pesar de que puede ser introducido mediante sencillos razonamientos geom´etricos, los aspectos m´as te´oricos de dicho algoritmo (como la demostraci´on de su convergencia) suelen ser m´as sutiles que los de m´etodos alternativos. Este hecho ha motivado que la mayor´ıa de textos bien presenten el m´etodo de forma muy superficial, bien reduzcan su exposici´on a la demostraci´on de una lista de teoremas y proposiciones. En este trabajo, de car´acter principalmente docente, se ha pretendido compensar ambas carencias, de forma que, manteniendo el rigor, la presentaci´on sea adecuada tanto para estudios de Ingenier´ıa como de Matem´aticas. Se incluye un apartado con una serie de problemas resueltos sobre el tema, entre los que aparece una implementaci´on del algoritmo en el lenguaje Matlab. Esta implementaci´on se utiliza para solucionar un problema de ingenier´ıa de pocas variables, el cual se propone como ejercicio pr´actico a realizar. Tambi´en se presenta una comparativa computacional entre una implementaci´on eficiente del m´etodo del escalado af´ın y una del algoritmo del s´ımplex en la soluci´on de un total de 90 problemas est´andar de programaci´on lineal. Los resultados obtenidos ratifican la eficiencia del algoritmo del escalado af´ın. An introduction to the affine scaling algorithm for linear programming. Keywords: Algoritmo de escalado af´ın, M´etodos de punto interior, Programaci´on lineal. Clasificaci´ on AMS: 90C05, 90C06.

1. Introducci´ on Los m´etodos de punto interior han sido utilizados durante los u ´ltimos a˜ nos con gran ´exito para solucionar problemas de programaci´on lineal. Este ´exito ha provocado que actualmente se consideren una alternativa m´as efectiva que el m´etodo del s´ımplex para una gran variedad de problemas lineales. En este documento vamos a desarrollar uno de estos m´etodos: el del escalado af´ın. Este m´etodo fue el primero de los denominados de punto interior, y data del 1967 (Dikin (1967)), aunque no fue conocido en Occidente hasta dos d´ecadas m´as tarde. En la actualidad existen ya diversos textos que describen dicho algoritmo. La mayor´ıa de ellos, sin embargo, no realizan una presentaci´on que balancee los aspectos te´oricos con los m´as intuititivos y geom´etricos del algoritmo. As´ı, por ejemplo, mientras las obras de Arbel (1993) y Vanderbei (1996) realizan una presentaci´on de car´acter muy introductorio, en Tsuchiya (1996) la exposici´ on se reduce principalmente a un conjunto de teoremas y proposiciones. Algunos textos, como Bertsimas y Tsitsiklis (1997), s´ı realizan el esfuerzo de presentar ambos aspectos, aunque, en este caso, la presentaci´on se centra en la versi´on del algoritmo del escalado af´ın denominada de paso corto, la cual ha mostrado ser la menos eficiente. A diferencia de estos textos, en este documento nos hemos centrado en la variante m´as eficiente —de paso largo—, presentando

1

una demostraci´on de convergencia de la misma, sin olvidar los aspectos m´as geom´etricos del algoritmo. Iniciaremos la exposici´on presentando la notaci´on que se utilizar´a (secci´on 2) e introduciendo el problema lineal que pretendemos solucionar as´ı como los requisitos necesarios para el seguimiento de la exposici´on (secci´on 3). En la secci´on 4 se dar´a una visi´on global de las diferentes t´ecnicas de punto interior para programaci´on lineal. Acto seguido pasaremos a exponer en la secci´on 5 el m´etodo del escalado af´ın. Se realiza una presentaci´on constructiva del m´etodo, introduciendo en cada apartado sus diferentes etapas algor´ıtmicas. Se incluyen diversos ejemplos para ilustrar cada una de estas etapas. La secci´on 5 finaliza presentando los principales resultados sobre la convergencia del algoritmo, y una comparativa de su rendimiento computacional respecto el algoritmo del s´ımplex. En la secci´on 6 se listan una serie de ejercicios propuestos, con sus respectivas soluciones. Durante la presentaci´on del m´etodo se hacen diversas referencias a estos ejercicios. En la secci´on 7 se sugiere un trabajo pr´actico donde debe solucionarse un problema de optimizaci´on real de peque˜ na dimensi´on, utilizando una implementaci´on del algoritmo del escalado af´ın en el lenguaje Matlab. Se adjunta la soluci´on de dicho trabajo pr´actico. Las secciones 8 y 9 listan determinados aspectos del algoritmo no tratados en el documento y la bibliograf´ıa utilizada, respectivamente. 2. Notaci´ on Se seguir´a el siguiente convenio a lo largo del documento: Los t´erminos con supra´ındice (como xk ) har´an referencia a los diferentes puntos obtenidos por el algoritmo. Los t´erminos que tengan un ∗ como supra´ındice (como x∗ ) denotar´an puntos ´optimos o puntos soluci´on. Las expresiones con sub´ındices (como xi ) denotar´an componentes de vectores. Las letras griegas representar´an escalares. Las letras may´ usculas (A, X, etc.) ser´an matrices. Las letras min´ usculas podr´an representar tanto vectores (x) como dimensiones de vectores (n, m), y se distinguir´a seg´ un el contexto. La letra k se reserva para indicar el n´ umero de iteraci´on del algoritmo (xk denotar´a, por tanto, la soluci´on que se tiene en la iteraci´o k-´esima). La letra e denotar´a el vector todo de unos e = (1 , . . . , 1)T . Se seguir´a la convenci´on de que los vectores ser´an vectores columna, y sus traspuestos vectores fila. 3. Presentaci´ on del problema y requisitos necesarios El problema lineal que queremos solucionar (y que denotaremos por (P)) estar´a en forma est´andar: min cT x x

sujeto a

Ax = b

(P)

x≥0 , donde x ∈ IRn es nuestro vector de variables, A ∈ IRm×n es la matriz de restricciones del problema (donde m < n), b ∈ IRm el vector de t´erminos independientes, y c ∈ IRn el vector de costes lineales. El conjunto de puntos {x | Ax = b, x ≥ 0} se denomina regi´on factible (o pol´ıtopo factible). Un punto que pertenece a este conjunto se denomina punto factible. Consideraremos las siguientes hip´otesis acerca de (P). Hip´ otesis 1. La matriz A es de rango completo, esto es, rango(A)= m. Hip´ otesis 2. El problema (P) no es degenerado. Esto significa que todo punto x factible para (P) tiene como m´ınimo m componentes diferentes de 0.

2

La interpretaci´on de no-degeneraci´on coincide con la del m´etodo del s´ımplex, donde se denomina degenerada a una soluci´on b´asica x = [xTB xTN ]T (xB ∈ IRm , xN ∈ IRn−m , xN = 0, xB ≥ 0) que tiene alguna componente de xB igual a 0. Esta hip´otesis facilita la presentaci´on y justificaci´on de determinadas etapas del m´etodo, as´ı como su convergencia. Entre los requisitos necesarios para desarrollar el m´etodo del escalado af´ın hallamos algunos conceptos b´asicos de teor´ıa de la dualidad. En primer lugar debemos saber que el problema (P) tiene asociado un problema dual (D), que tiene la siguiente formulaci´on: max y,z

sujeto a

bT y AT y + z = c

(D)

z ≥ 0, donde y ∈ IRm es el vector de variables duales, y z ∈ IRn es el vector denominado de holguras duales. Al igual que para el problema primal, consideraremos que el dual no es degenerado (esta hip´otesis ser´a u ´nicamente utilizada para demostrar la convergencia del m´etodo). Hip´ otesis 3. El problema (D) no es degenerado. Esto significa que todo punto (y, z) factible para (D) tiene como m´aximo m componentes de z iguales a 0. Las Hip´otesis 2 y 3 garantizan que los problemas primal y dual, si son factibles, tienen un u ´nico punto ´optimo. Las relaciones entre los problemas (P) y (D) nos vienen dadas por los Teoremas de la Dualidad D´ebil, Dualidad Fuerte, y de la Holgura Complementaria, que tambi´en ser´an utilizados durante la exposici´on: Teorema 1. (Teorema de la Dualidad D´ebil) Si x es factible para el problema primal (P), e y es factible para el problema dual (D), entonces cT x ≥ bT y.

Teorema 2. (Teorema de la Dualidad Fuerte) Si el problema primal (P) tiene una soluci´on ´optima x∗ , entonces el problema dual tambi´en tiene una soluci´on ´optima y ∗ de forma que cT x∗ = bT y ∗ .

Teorema 3. (Teorema de la Holgura Complementaria) Si x es factible para el problema primal (P), y el par (y, z) es factible para el problema dual (D), entonces son ´optimos de sus respectivos problemas si y solamente si xi zi = 0

para todo i = 1, . . . , n.

La distancia cT x − bT y entre las funciones objetivo primal y dual se denominar´a gap dual. En el ´optimo de ambos problemas, el gap dual es 0. En determinados puntos del desarrollo se utilizar´an tambi´en algunos conceptos de programaci´on no lineal, como las condiciones de optimalidad de problemas sujetos a restricciones de igualdad (m´etodo de los multiplicadores de Lagrange). Estos conceptos se usar´an para justificar formalmente determinados pasos del algoritmo. En particular, usaremos el siguiente resultado:

3

Teorema 4. (Condiciones necesarias y suficientes de optimalidad para un problema con restricciones de igualdad) Dado el problema min f (x) x

sujeto a h(x) = 0, donde f : IRn → IR y g : IRn → IRm son continuas hasta sus segundas derivadas, y considerando la funci´on Lagrangiana definida como L(x, y) = f (x) − y T h(x), donde y ∈ IRm recibe el nombre de vector de multiplicadores de Lagrange, entonces es condici´on necesaria para que x∗ sea m´ınimo del problema que, si los gradientes de las restricciones en x∗ son independientes, exista un vector de multiplicadores y ∗ tal que ½ ∇f (x∗ ) − ∇h(x∗ )T y ∗ = 0 ∇(x,y) L(x∗ , y ∗ ) = 0 ⇐⇒ h(x∗ ) = 0 . Y es condici´on suficiente para que sea m´ınimo que, adem´as de la condici´on necesaria anterior, se verifique v T [∇2xx L(x∗ , y ∗ )]v > 0

para todo v 6= 0 tal que ∇x h(x∗ )v = 0 .

4. M´ etodos de punto interior para programaci´ on lineal Los m´etodos de punto interior constituyen una familia de t´ecnicas no-s´ımplex para programaci´on lineal. De hecho, puede decirse que se basan en la aplicaci´on a problemas lineales de m´etodos que cl´asicamente hab´ıan sido considerados de programaci´on no lineal. Por ello la mayor´ıa de obras recientes de programaci´on no lineal incluyen este tipo de t´ecnicas (Bertsekas (1995)). Se denominan m´etodos de punto interior precisamente porque los puntos generados por estos algoritmos se hallan en el interior de la regi´on factible. Esta es una clara diferencia respecto el m´etodo del s´ımplex, el cual avanza por la frontera de dicha regi´on movi´endose de un punto extremo a otro. La Fig. 1 muestra gr´aficamente ambos comportamientos. Puede comprobarse como el m´etodo de punto interior requiere que el punto inicial x0 sea tambi´en un punto interior a la regi´on factible. x

x

x*

2

x*

2

x0

x0 a)

x

b)

1

Figura 1. Trayectoria seguida para alcanzar el punto ´optimo x∗ desde el punto inicial x0 por a) el m´etodo del s´ımplex b) un m´etodo de punto interior

4

x

1

Otra de las diferencias fundamentales respecto el m´etodo del s´ımplex es que existen algoritmos de punto interior polin´omicos. Esto significa que obtienen el punto ´optimo en un n´ umero de iteraciones que es funci´on polin´omica del tama˜ no del problema. El m´etodo del s´ımplex, en el peor de los casos, tiene un coste exponencial. En la actualidad los m´etodos de punto √ interior m´as eficientes tienen una complejidad de O( nL), donde n es el n´ umero de variables del problema y L una medida del tama˜ no del problema (a saber, el n´ umero de bits necesarios para representar los datos del problema). Podemos clasificar los m´etodos de punto interior en cuatro grandes categor´ıas: M´ etodos de escalado af´ın. Son los m´as sencillos, y transforman el problema original mediante un escalado af´ın de las variables. En la secci´on 5 se presentan con todo detalle. M´ etodos basados en transformaciones proyectivas. Al igual que los anteriores transforman el problema original, aunque de forma m´as elaborada, utilizando lo que se denomina una transformaci´on proyectiva. Estas transformaciones son del tipo: x e=

X −1 x eT X −1 x

,

donde x ∈ IRn es el punto actual, x e representa la imagen de x en el nuevo espacio de variables, X es una matriz diagonal que tiene como t´erminos diagonales las componentes de x (X = diag(x1 , . . . , xn )), y eT es el vector n-dimensionald de unos (eT = (11 , . . . , 1n )). El ahora ya famoso primer algoritmo de punto interior, desarrollado por Karmarkar en 1984, pertenece a esta categor´ıa. Los denominados m´ etodos path-following (se hace dif´ıcil traducir este t´ermino, por lo que preferimos mantener el original en ingl´es). Estos m´etodos utilizan una t´ecnica para problemas no lineales, denominada barrera logar´ıtmica. En vez de solucionar (P) directamente, solucionan una secuencia de problemas de la forma min cT x − µ x

sujeto a

n X

ln(xj )

j=1

(Pµ )

Ax = b .

Puede probarse que, a medida que µ → 0 la secuencia de soluciones de (Pµ ) se acerca a la soluci´on del problema (P) original. El conjunto de soluciones de (Pµ ) para diferentes µ proporciona el denominado camino central (central path). Los m´etodos path-following siguen una trayectoria cercana a este camino central en su b´ usqueda del ´optimo. El algoritmo primal-dual de punto interior pertenece a esta categor´ıa. Este algoritmo ha mostrado ser, en general, el m´etodo de punto interior m´as eficiente. Hoy en d´ıa es considerado como la mejor opci´on para solucionar problemas lineales de gran dimensi´on. Los m´ etodos de reducci´ on de potencial, que intentan reducir una determinada funci´on de m´erito o de potencial. Nosotros nos centraremos en la primera categor´ıa de m´etodos: los de escalado af´ın. 5. M´ etodo del escalado af´ın Este tipo de m´etodos son los m´as sencillos de todos los de punto interior. A pesar de ello, muestran un excelente rendimiento en la soluci´on de problemas de gran dimensi´on. El m´etodo del escalado af´ın fue sugerido originalmente por el matem´atico ruso Dikin (Dikin (1967)) en 1967. Este trabajo, sin embargo, pas´o desapercibido hasta que un par de d´ecadas m´as tarde el m´etodo presentado fue redescubierto en Occidente de forma independiente por diferentes 5

investigadores como una simplificaci´on del algoritmo de las transformaciones proyectivas de Karmarkar (Barnes (1986), Vanderbei et al. (1986)). Las principales caracter´ısticas del m´etodo del escalado af´ın son: Es de comprensi´on sencilla. Puede ser f´acilmene implementado. Por lo general, proporciona muy buenos resultados. Puede verse como una aplicaci´on de una variante de un m´etodo de programaci´on no lineal —denominado del gradiente proyectado— a la soluci´on de un problema lineal. 5.1. Algoritmo b´ asico Nuestro objetivo es solucionar el problema lineal formulado en (P) minx cT x sujeto a Ax = b, x ≥ 0 (recordemos que contamos con n variables y m restricciones). Para desarrollar el algoritmo del escalado af´ın consideraremos que disponemos de un punto interior y factible x0 ∈ IRn . Que sea factible significa que satisface las restricciones del problema. Que sea interior quiere decir que no contacta con las paredes de la regi´on factible, esto es, que x > 0. Por tanto, nuestro punto interior y factible x0 inicial garantiza: Ax0 = b

y

x0 > 0.

EJEMPLO 1

Dado el problema lineal siguiente (a la izquierda se muestra el problema original y a la derecha su forma est´ andar tras haber introducido las variables x3 , x4 y x5 ) min

min

− 3x1 − 2x2

− 3x1 − 2x2

sujeto a

sujeto a 4x1 3x1 x1

−2x2 +4x2 + x2

x1 ≥ 0

≤ ≥ ≤

5 1 2

4x1 3x1 x1

=⇒

−2x2 +4x2 + x2

+x3 − x4 + x5

= = =

5 1 2

xi ≥ 0 para toda i = 1, . . . , 5

x2 ≥ 0

puede comprobarse como el punto x0 = (1/2 1/2 4 5/2 1)T es interior (todas las componentes son mayores que 0) y factible (satisface las restricciones del problema en forma est´ andar): 4 · 1/2 −2 · 1/2 +4 = 5 3 · 1/2 +4 · 1/2 − 5/2 = 1 1/2 + 1/2 +1 = 2 La Fig. 2 muestra la situaci´ on del punto x0 dentro de la regi´ on factible, resaltada en oscuro. La frontera de la regi´ on factible viene delimitada por las tres restricciones (de desigualdad). Se consideran u ´nicamente las dos primeras componentes de x0 , puesto que usamos la formulaci´ on original del problema, no la forma est´ andar, para permitir su representaci´ on en el plano.

El algoritmo b´asico seguido por el m´etodo del escalado af´ın consiste en, comenzando por x , generar una secuencia de puntos {xk } que nos lleven hasta el punto ´optimo x∗ del problema. As´ı, partiendo de x0 avanzaremos hasta x1 , de ´este hasta x2 , y as´ı sucesivamente, hasta que se verifique alguna condici´on que nos indique que el punto actual xk es soluci´on de nuestro problema. Cada nuevo punto se obtendr´a del anterior, siguiendo el siguiente procedimiento iterativo xk+1 = xk + αk ∆xk , αk ≥ 0. (1) 0

6

x2

2 rest. 1 rest. 2 rest. 3

1.75 1.5 1.25 1 0.75 x0 0.5 0.25 0

0

0.25

0.5

0.75

1

1.25

1.5

1.75

2 x1

Figura 2. Situaci´on del punto x0 dentro de la regi´on factible delimitada por las tres restricciones del ejemplo 1: rest. 1 ≡ 4x1 − 2x2 ≤ 5, rest. 2 ≡ 3x1 + 4x2 ≥ 1, rest. 3 ≡ x1 + x2 ≤ 2.

Aqu´ı ∆xk ∈ IRn representa nuestro vector o direcci´on de movimiento, mientras que el escalar αk , denominado longitud de paso, es positivo y nos indica cuanto nos alejamos de xk a lo largo de la direcci´on calculada. Este esquema es com´ un a la mayor´ıa de m´etodos de Optimizaci´on. Por tanto debemos centrarnos en las dos siguientes cuestiones: i) ¿C´omo obtener la direcci´on de movimiento ∆xk ? y ii) ¿C´omo calcular la longitud de paso αk ? La primera cuesti´on va a ocupar la mayor parte del resto de secciones. El c´alculo de la longitud de paso es sencillo una vez se ha obtenido la direcci´on de movimiento, tal y como veremos. 5.2. C´ alculo de la direcci´ on de movimiento La direcci´on de movimiento ∆x calculada a cada iteraci´on del algoritmo no puede ser una direcci´on cualquiera. Concretamente, ∆x debe verificar dos condiciones para ser considerada una “buena” direcci´on: 1) Debe preservar la factibilidad del nuevo punto (direcci´on factible). 2) Debe mejorar el valor de la funci´on objetivo (mejorar significa disminuir en este caso, puesto que el problema (P) es de minimizaci´on). Estas condiciones son tratadas en los dos siguientes apartados. (N´ otese que en el p´arrafo anterior hemos escrito ∆x y no ∆xk , como en la secci´on previa. Para simplificar la notaci´on preferimos eliminar el supra´ındice k . Por tanto ∆x indicar´a la direcci´on de movimiento en el punto actual xk . Seguiremos este acuerdo para la mayor´ıa de t´erminos que aparecer´an de ahora en adelante —entre ellos la longitud de paso, que se denotar´a por α.) 5.2.1. Obtenci´ on de una direcci´ on factible Supongamos que nos hallamos en el punto factible e interior xk . En el siguiente punto xk+1 queremos que se mantenga la factibilidad de nuestra soluci´on. De momento nos preocuparemos u ´nicamente por las restricciones de igualdad Ax = b; las restricciones x ≥ 0 pueden ser satisfechas de forma simple, y ser´an consideradas m´as adelante. Por tanto, el nuevo punto xk+1 debe garantizar Axk+1 = b.

7

Sabemos por (1) que el punto xk+1 se obtiene como xk+1 = xk + α∆x, y que el punto xk verifica Axk = b por ser factible. Entonces, la condici´on de factibilidad de xk+1 se reduce a  b = Axk+1     k = A(x + α∆x)  ⇒ A∆x = 0. = Axk + αA∆x      = b + αA∆x Conclu´ımos que, para que xk+1 sea factible, debe verificarse que la direcci´on de movimiento debe pertenecer al espacio nulo de A definido como N (A) = {v | Av = 0, v ∈ IRn }. Sin embargo, dado un vector cualquiera d ∈ IRn , ´este siempre puede ser proyectado sobre N (A), si A tiene rango completo, premultiplic´andolo por la siguiente matriz de proyecci´on ortogonal P = In − AT (AAT )−1 A,

(2)

donde In denota la matriz identidad de dimensi´on n. En el ejercicio 1 se plantea la deducci´on de este resultado. Esta matriz de proyecci´on ortogonal verifica las siguientes propiedades AP = 0 P = PT

(simetr´ıa)

2

P =P

(3)

(idempotencia).

(ver el ejercicio 2 para una comprobaci´on de las mismas). EJEMPLO 2

Consideremos el siguiente problema lineal 1 2 x1 − x2 3 3 x1 + x2 + x3 = 1

min

−

sujeto a

x1 ≥ 0, x2 ≥ 0, x3 ≥ 0. En este problema m = 1, n = 3 y A = ( 1 1 1 ), y el ´ optimo viene dado claramente por x∗ = ( 0 1 0 )T . El punto actual factible es xk = ( 1/3 1/3 1/3 )T . La matriz de proyecci´ on ortogonal sobre el espacio nulo de A se calcula como:

à T

T −1

P = In − A (AA )

à =

A=

!

1 1

à − 1/3

1

1 1 1 1

1 1 1

Ã

2 −1 −1

1 1 1

−

!

−1 2 −1

1 1 1

à = 1/3

−1/2

à !!−1

à !Ã

1

Dado un vector cualquiera d = ( 0 espacio nulo de A P d = 1/3

!

1

(1 2 −1 −1

1

−1 2 −1

1 1 1

1) −1 −1 2

!

(1

1

1)

.

0 )T , podemos calcular su proyecci´ on sobre el

−1 −1 2

!Ã

8

0 −1/2 0

!

à =

1/6 −1/3 1/6

! .

Podemos comprobar que la direcci´ on P d as´ı obtenida es factible observando que el punto xk + P d satisface la restricci´ on de igualdad del problema:

à k

x + Pd =

1/3 1/3 1/3

!

à +

1/6 −1/3 1/6

!

à =

1/2 0 1/2

! 1/2 + 0 + 1/2 = 1.

La Fig. 3 a) ilustra gr´ aficamente este ejemplo. x3

x3

Pd d

xk

xk -c x*

x* x2

x2 -Pc

x1

x1 a)

b)

Figura 3. Representaci´on gr´afica de los ejemplos 2 y 3. Se detalla la regi´on factible {x1 + x2 + x3 = 1, xi ≥ 0}, y sobre ´esta el punto actual xk , y las curvas de nivel de la funci´ on objetivo. a) Proyecci´ on ortogonal del vector d = ( 0 −1/2 0 )T . El nuevo punto xk + P d es ( 1/2 0 1/2 )T y pertenece a la regi´ on factible. b) Proyecci´ on ortogonal del vector −c = ( 1/3 2/3 0 )T . El nuevo punto xk − P c es ( 1/3 2/3 0 )T y mejora el valor de la funci´ on objetivo.

Por tanto, dada una direcci´on de movimiento cualquiera d ∈ IRn , siempre se verificar´a que su proyecci´on ortogonal sobre el espacio nulo de A, P d, es una direcci´on factible para nuestro problema. Esto es, que AP d = 0 (se deduce inmediatamente de la primera propiedad de (3)). Sin embargo, cualquier direcci´on no siempre nos permitir´a mejorar nuestra funci´on de coste. As´ı, el punto xk + P d obtenido en el ejemplo 2 tiene un coste superior al del punto original xk : cT (xk + P d) = −1/3 · 1/2 − 2/3 · 0 = −1/6 > cT xk = −1/3 · 1/3 − 2/3 · 1/3 = −1/3. Tambi´en se observa este hecho en la Fig. 3 a): nos hemos movido hacia curvas de nivel de coste superior al del punto xk . Nuestro objetivo es ahora escoger, de entre todas las direcciones posibles, una que nos garantice una mejora de la funci´on objetivo. 5.2.2. Mejora de la funci´ on objetivo En un problema de minimizaci´on como el que estamos considerando, en el nuevo punto xk+1 debemos tener un coste inferior al del punto xk . Esto significa que cT xk+1 ≤ cT xk . Utilizando que xk+1 = xk + α∆x, y teniendo en cuenta que α es siempre positivo, podemos escribir la condici´on anterior como cT xk+1 = cT (xk + α∆x) = cT xk + αcT ∆x ≤ cT xk Esta condici´on se denomina condici´on de descenso.

9

=⇒

cT ∆x ≤ 0.

Consideremos la siguiente direcci´on ∆x ∆x = −P c,

(4)

que no es sino la proyecci´on de −c (menos el gradiente de nuestra funci´on lineal) sobre el espacio nulo de A. Por el apartado anterior sabemos que esta direcci´on es factible. Y utilizando la simetr´ıa e idempotencia de P puede demostrarse que tambi´en es de descenso: cT ∆x = −cT P c = −cT P 2 c = −cT P T P c = −||P c||2 ≤ 0. A pesar de aparecer un ≤ en la u ´ltima desigualdad anterior, en la pr´actica siempre ser´a de 0, por lo que la expresi´on (7) nos proporciona el m´ınimo de nuestro problema. Esto justifica el por qu´e hemos considerado un signo positivo para el t´ermino λ en la funci´on Lagrangiana: si hubi´eramos utilizado un signo menos, la soluci´on obtenida hubiera sido la direcci´on contraria a (7), verific´andose entonces que ∇∆x∆x = −2λ, lo que nos indicar´ıa que esta direcci´on contraria corresponde a un m´aximo de nuestro problema en vez de a un m´ınimo. Dejando a parte el escalar ρ/||P c||, la direcci´on obtenida coincide con la presentada en (4). ´ Esta es la direcci´on que nos garantiza el m´aximo descenso a partir de un punto manteniendo la factibilidad del siguiente. Recibe el nombre de direcci´on de (menos) el gradiente proyectado (omitiremos el “menos” al referirnos a esta direcci´on de aqu´ı en adelante), y coincide con la utilizada en el m´etodo del gradiente proyectado para Programaci´on No Lineal. La Fig. 4 muestra la proyecci´on del gradiente sobre la regi´on factible de nuestro problema. 5.3. Escalado El gradiente proyectado nos proporciona una excelente direcci´on, en el sentido que entre todas las direcciones factibles es la que tiene la m´axima disminuci´on de cT x por unidad de longitud de paso. Sin embargo, y seg´ un la expresi´on (1) que rige el proceso iterativo del m´etodo, si nos hallamos muy cerca de las paredes de la regi´on factible definidas por las restricciones x ≥ 0, esta longitud de paso α que podremos realizar ser´a peque˜ na (de otro modo, violar´ıamos la no-negatividad de alguna variable). No es suficiente, por tanto, con tener una buena direcci´on de descenso; conviene tambi´en estar en un punto centrado respecto de la regi´on factible (esto es, alejado de las paredes delimitadas por x ≥ 0), lo que nos permitir´a tener α À 0. Este hecho se observa claramente en la Fig. 5, donde el movimiento desde el punto xa es peque˜ no comparado con el realizado desde el punto xb . El m´etodo del escalado af´ın solventa la dificultad de no tener un punto centrado mediante el siguiente procedimiento de 3 etapas, que se repiten a cada iteraci´on del algoritmo: 1. Escalar el problema de forma que el punto actual se encuentre alejado de las paredes delimitadas por x ≥ 0. El punto actual xk queda transformado en x ek , en el nuevo espacio de variables x e.

12

x* -Pc xa xb

Figura 5. Movimiento realizado en la direcci´on del gradiente proyectado −P c desde un punto “centrado” xb y desde un punto “no-centrado” xa .

2. Calcular la direcci´on del gradiente proyectado en el nuevo problema escalado. En este nuevo espacio de variables el gradiente proyectado se denota por ∆e x. 3. Trasladar la direcci´on ∆e x obtenida al problema original, deshaciendo el escalado previo. A partir de ∆e x obtenemos por tanto la direcci´on ∆x en el espacio original de variables. Ahora debemos determinar qu´e tipo de escalado debe utilizarse. Cualquier escalado no sirve. As´ı, por ejemplo, un escalado que multiplicase todas las componentes xki i = 1, . . . , n por un mismo valor no ser´ıa v´alido para nuestros prop´ositos, porque se mantendr´ıa la distancia relativa de todas las variables respecto las paredes definidas por x ≥ 0. La estrategia seguida por el m´etodo del escalado af´ın consiste en transformar cada come. En este ponente xki i = 1, . . . , n del espacio de variables x en un 1 en el nuevo espacio x nuevo espacio todas las componentes del punto x ek se encuentran alejadas por igual de 0. Para conseguir esto debemos dividir cada componente xki por ella misma. El cambio de variables que se realiza es por tanto: xi x ei = k i = 1, . . . , n, xi donde xk es el punto actual de iteraci´on. Este cambio se deshace de forma trivial: xi = x ei xki

i = 1, . . . , n.

Si denotamos por X k la matriz diagonal que tiene por t´erminos diagonales las componentes del punto actual xki ,  xk  1 k x2   , Xk =  ..   . xkn podemos escribir las operaciones anteriores de cambio de variable en forma matricial: x e = (X k )−1 x x = Xkx e.

(8)

Observemos que este escalado est´a bien definido porque siempre nos movemos por puntos interiores de la regi´on factible, con lo que xki > 0, siendo entonces X k no singular y definida positiva.

13

Veamos ahora en los tres apartados siguientes como se realizan a cada iteraci´on del algoritmo las 3 etapas antes expuestas. Supondremos que nos hallamos en un punto xk factible e interior determinado. 1. Escalado del problema original: x → x e El problema primal (P) que estamos solucionando se escala utilizando el cambio de variables x = Xkx e de la siguiente manera: min cT x x

sujeto a

(P) −→

Ax = b

min cT X k x e e x sujeto a AX k x e=b

x≥0

−→

min e cT x e e x ex = b sujeto a Ae

Xkx e≥0

e (P)

x e≥0 ,

donde, en el u ´ltimo paso, se ha utilizado la notaci´on e c = Xkc e = AX k . A

(9)

Tambi´en conviene notar que en el u ´ltimo paso las restricciones X k x e ≥ 0 se han transformado en x e ≥ 0. Ambos grupos de restricciones son equivalentes gracias a que las componentes del punto actual de iteraci´on xk son positivas. El punto xk que ten´ıamos en (P) se transforma en el punto x ek = e = (1, . . . , 1)T que tiene e La etapa 1 de la Fig. 6 muestra este proceso de escalado. Se una posici´on centrada en (P). observa como el escalado modifica las curvas de nivel de la funci´on objetivo a lo largo de la regi´on factible. Problema original

(x)

Problema escalado

( ~x )

~ x3

1)

~ = (X k ) -1 x x

x3

e

x* xk

∆x

2) ∆ ~x

x2

~* x ~x 2

x~1

x1

3) ∆ x = (X k ) ∆ ~x

Figura 6. Representaci´on gr´afica de las tres etapas del proceso de escalado realizado a cada iteraci´ on del m´etodo del escalado af´ın. 1) Escalado del problema original, pasando del espacio de variables x al nuevo espacio e x. El punto xk queda transformado en el punto e = (1 , . . . , 1)T en el nuevo espacio. eec en el 2) C´ alculo de la direcci´ on de movimiento ∆e x = −P nuevo espacio de variables. 3) Retorno al problema original, obteniendo la direcci´ on de movimiento ∆x a partir de ∆e x aplicando el escalado inverso al utilizado en la etapa 1.

14

2. C´ alculo de ∆e x En el nuevo espacio de variables x e la direcci´on de descenso ∆e x ser´a la del gradiente proyectado, calculada seg´ un la expresi´on (4). Sin embargo ahora utilizaremos la nueva matriz de e y el nuevo vector de costes e eye restricciones A c. Utilizando la definici´on de A c en (9), y teniendo k en cuenta que X es una matriz diagonal, el c´alculo de ∆e x es como sigue: ∆e x = −Pee c eT (A ec eA eT )−1 A)e = −(In − A = −(In − (AX k )T (AX k (AX k )T )−1 AX k )X k c = −X k (In − AT (A(X k )2 AT )−1 A(X k )2 )c.

(10)

La etapa 2 de la Fig. 6 muestra la direcci´on ∆e x en el problema escalado. 3. Retorno al problema original: x e→x Si nos movi´esemos en el problema escalado lo har´ıamos de la siguiente forma (ignoremos por el momento la longitud de paso): x ek+1 = x ek + ∆e x. Si multiplicamos por X k la expresi´on anterior deshacemos el cambio de variables realizado, seg´ un (8): Xkx ek+1 = X k x ek + X k ∆e x xk+1 = xk + X k ∆e x. Por tanto, el movimiento ∆x realizado en el espacio original de variables es igual a X k ∆e x. Utilizando la expresi´on (10) de ∆e x obtenida en la etapa 2 tenemos: ∆x = X k ∆e x = −(X k )2 (In − AT (A(X k )2 AT )−1 A(X k )2 )c = −D(c − AT (ADAT )−1 ADc)

[D = (X k )2 ]

= −D(c − AT y)

[y = (ADAT )−1 ADc]

= −Dz

[z = c − AT y],

(11)

donde D, y y z se han definido como como D = (X k )2 y = (ADAT )−1 ADc

(12)

T

z = c − A y. La etapa 3 de la Fig. 6 muestra la direcci´on ∆x obtenida en el problema original. En el ejercicio 4 se pide comprobar que esta direcci´on del gradiente proyectado con escalado es factible y de descenso. N´ otese que para obtener la direcci´on ∆x debemos factorizar la matriz ADAT = A(X k )2 AT . Para ello debe verificarse que AX k sea una matriz de rango completo. Sin embargo, sabemos que A es de rango completo (Hip´otesis 1) y que X k tiene como m´ınimo m componentes > 0 por ser el problema no-degenerado (Hip´otesis 2), concluyendo que A(X k )2 AT es no singular. La siguiente proposici´on, que no demostramos, recoge este resultado.

15

Proposici´ on 1. Si (P) no es degenerado y la matriz de restricciones A tiene rango completo, entonces A(X k )2 AT es no singular. Adem´as, los t´erminos y y z definidos en (12) son una funci´on continua de xk . A cada iteraci´on del algoritmo las tres etapas anteriores se realizan impl´ıcitamente, y u ´nicamente es necesario calcular la direcci´on ∆x seg´ un la expresi´on (11). Esta ser´a la direcci´on a utilizar en el proceso iterativo del m´etodo presentado en (1). El siguiente ejemplo ilustra todo lo expuesto hasta el momento. EJEMPLO 4

Consideremos el problema siguiente (que ya fue utilizado en el ejemplo 1): min

− 3x1 − 2x2

min

sujeto a

− 3x1 − 2x2

sujeto a 4x1 3x1 x1

−2x2 +4x2 + x2

x1 ≥ 0

≤ ≥ ≤

5 1 2

4x1 3x1 x1

=⇒

x2 ≥ 0

−2x2 +4x2 + x2

+x3 − x4 + x5

= = =

5 1 2

xi ≥ 0 para toda i = 1, . . . , 5

y el punto actual xk = ( 1.25 0.1 0.2 3.15 0.65 )T . Este punto es factible e interior —ver Fig. 7 a)—, y tiene un coste de −3 · 1.25 − 2 · 0.1 = −3.95. En este problema tenemos que

à A=

4 3 1

−2 4 1

1 0 0

0 −1 0

0 0 1

!

 −3   −2   c=  0  0 0



 1.25 

Xk =  

0.1

 . 

0.2 3.15 0.65

Utilizando la matriz diagonal X k anterior, y aplicando la ecuaci´ on (9), podemos construir nuestro problema escalado: min

− 3.75 e x1 − 0.2 e x2

sujeto a 5e x1 3.75 e x1 1.25 e x1

−0.2 e x2 +0.4 e x2 +0.1 e x2

+0.2 e x3 −3.15 e x4 +0.65 e x5

= = =

5 1 2

e xi ≥ 0 para toda i = 1, . . . , 5. Utilizando las expresiones anteriores de A, c y X k , la direcci´ on de movimiento en este nuevo espacio de variables se calcula seg´ un (10)

 0.0183  0.319 ∆e x = −X k (In − AT (A(X k )2 AT )−1 A(X k )2 )c =   −0.140

0.0623 −0.0844

  . 

La regi´ on factible del problema escalado y la direcci´ on de movimiento ∆e x obtenida se representan en la Fig. 7 b), considerando las variables e x3 , e x4 y e x5 como holguras. N´ otese la diferencia de escala entre el eje de ordenadas y el de abcisas, debido al valor tan cercano a 0 que tiene la componente xk2 = 0.1.

16

x2

~ x2

2 rest. 1 rest. 2 rest. 3

1.75

30 rest. 1 rest. 2 rest. 3

25

1.5 20 1.25 1

15

0.75

10

x*

0.5 0.25 0

xk

x k+1

∆~ x

5

~ xk

xn 0

0

0.25

0.5

0.75

1

1.25

1.5

1.75

0

2 x1

b) Problema escalado

Representaci´ on gr´ afica del ejemplo 4. a) Problema original. Situaci´ on de los puntos xk , k+1 n x y x dentro de la regi´ on factible (en oscuro). En gris se resaltan las curvas de nivel de la funci´ on objetivo −3x1 − 2x2 . b) Problema escalado. Se muestra el punto e xk = T (1, . . . , 1) , la direcci´ on ∆e x, y las curvas de nivel de la funci´ on objetivo del problema escalado −3.75e x1 − 0.2e x2 .

Como ya contamos con ∆e x, la direcci´ on de movimiento en el espacio de variables original la obtenemos de la siguiente forma:

 0.0229   0.0318   ∆x = X k ∆e x=  −0.0280  . 0.196 −0.0548

Podr´ıamos haber calculado directamente ∆x a traves de (11). En esta direcci´ on podemos avanzar, sin violar la no-negatividad de las variables, hasta el punto

xk+1

2 ~ x1

a) Problema original

Figura 7.

1

 1.414    0.328 −7  =  5 · 10  . 4.55 0.258

Veremos m´ as adelante en la subsecci´ on 5.4 como se determina en la pr´ actica este nuevo punto. El coste de la funci´ on objetivo en xk+1 es menor que en el punto xk : cT xk+1 = −3 · 1.414 − 2 · 0.328 = −4.8984 < cT xk = −3.95. Puede comprobarse que si no hubi´eramos realizado ning´ un tipo de escalado en el problema original, hubi´eramos obtenido la siguiente direcci´ on de movimiento (del gradiente proyectado):

 0.0984   0.0685   ∆xn = −P c =   −0.257  . 0.569 −0.167

17

Y avanzando por esta direcci´ on, hasta donde nos permitan los l´ımites x ≥ 0 de las variables, hubi´eramos obtenido el punto  1.32   0.153  −7  xn =   5 · 10  , 3.59 0.520 que tiene un coste superior al del punto xk+1 cT xn = −3 · 1.32 − 2 · 0.153 = −4.2870 > cT xk+1 = −4.8984. Por tanto, el escalado ha sido una estrategia efectiva, puesto que nos ha permitido reducir m´ as nuestra funci´ on objetivo. En la Fig. 7 a) se muestran los puntos xk+1 y xn obtenidos.

Del ejemplo anterior debemos destacar dos cosas. La primera es que el gradiente proyectado sin escalado nos ha proporcionado una direcci´on de movimiento peor que con escalado, en el sentido que hemos obtenido un coste superior. Sin embargo tiene un inconveniente todav´ıa m´as grave. Si volvi´eramos a iterar con el gradiente proyectado, obtendr´ıamos la misma direcci´on de movimiento, puesto que ´esta es siempre constante e igual a −P c. La consecuencia directa de esto es que no podremos avanzar m´as en esta direcci´on una vez hallamos llegado a alguna pared definida por x ≥ 0. As´ı, pues, en el punto xn obtenido en el ejemplo 4 se habr´ıan acabado todas nuestras opciones de mejorar la funci´on objetivo. El escalado, por contra, evita esta situaci´on, ajustando la direcci´on de movimiento con la topolog´ıa de la regi´on factible. El segundo aspecto a destacar del ejemplo anterior es que no se ha explicitado como calcular el nuevo punto xk+1 una vez se tiene ∆x, esto es, qu´e longitud de paso α utilizar. Sin embargo lo que s´ı hemos observado es que se ha avanzado tanto como se ha podido en la direcci´on ∆x, hasta llegar a anular la componente de alguna variable. Esta decisi´on parece razonable, ya que la disminuci´on de coste que obtenemos al pasar de xk a xk+1 es de αcT ∆x (se deduce directamente premultiplicando a ambos lados de (1) por cT ), y cuanto mayor sea α, m´as disminuiremos nuestro coste. Utilizando este hecho, la siguiente subsecci´on nos muestra c´omo obtener la longitud de paso. 5.4. Determinaci´ on de la longitud de paso Nuestro objetivo ser´a movernos desde xk a lo largo de ∆x tanto como podamos mientras no violemos la no-negatividad de las variables. Esto es, debe garantizarse que xk+1 ≥ 0

=⇒

xk + α∆x ≥ 0.

Como α > 0, las componentes ∆xi ≥ 0 no pueden llevar la variable xi a 0. Por tanto u ´nicamente debemos preocuparnos por los t´erminos ∆xi < 0. El paso m´aximo, denotado por α que podemos realizar viene dado entonces por ¾ ½ xki para toda i tal que ∆xi ≤ 0 . (13) α = min − ∆xi Sin embargo, si defini´esemos α de esta manera la componente i-´esima de xk+1 , donde i es el ´ındice de la componente que ha proporcionado el m´ınimo −xki /∆xi con ∆xi < 0, se anular´ıa. Esto har´ıa que xk+1 no fuese interior, y por tanto la matriz X k+1 de escalado ser´ıa singular. La estrategia a seguir para evitar este inconveniente consiste en reducir el paso m´aximo calculado

18

seg´ un la expresi´on anterior multiplic´andolo por ρ ∈ (0, 1) (en la pr´actica se utilizan valores entre 0.95 y 0.9995). Entonces α queda definido como ½ ¾ xk α = ρ · α = ρ · min − i para toda i tal que ∆xi ≤ 0 ρ ∈ [0.95, 0.9995]. (14) ∆xi El ejercicio 5 presenta una estrategia de comprobaci´on de la ilimitaci´on de un problema que est´a estrechamente relacionada con el c´alculo de α. EJEMPLO 5

En el ejemplo 4 se iter´ o a partir del punto xk = ( 1.25 0.1 0.2 3.15 0.65 )T obteni´endose una direcci´ on ∆x = ( 0.0229 0.0318 −0.0280 0.196 −0.0548 )T . Utilizando un valor de ρ = 0.95, la longitud de paso α se calcular´ıa como:

½ α = ρ · min

−

xk3 xk ,− 5 ∆x3 ∆x5

¾

n = 0.95 · min −

0.2 0.65 ,− −0.0280 −0.0548

o

= 0.95 · min {7.143, 11.861} = 0.95 · 7.143 = 6.7885. El nuevo punto ser´ıa por tanto

 1.25 

xk+1

 0.0229   1.406   0.1   0.0318   0.316       = xk + α∆x =   0.2  + 6.7885  −0.0280  =  0.0101  . 3.15 0.65

0.196 −0.0548

4.483 0.278

de 0 N´ otese como en el nuevo punto se mantiene m´ as alejada la tercera componente xk+1 3 que en el ejemplo 4. Esto se debe a la utilizaci´ on del t´ermino ρ de reducci´ on de la longitud de paso.

5.5. Una digresi´ on: el escalado af´ın de paso corto Se ha observado como el m´etodo del escalado af´ın presentado intenta avanzar tanto como puede en la direcci´on de movimiento calculada. Por este motivo dicho m´etodo se conoce como escalado af´ın de paso largo. Es la variante m´as eficiente y fue descubierta hacia mediados de los ochenta (Vanderbei et al. (1986)). Recibe este nombre para diferenciarse del denominado escalado af´ın de paso corto, que, cron´ologicamente, apareci´o en primer lugar. Esta variante de paso corto es, de hecho, el m´etodo del escalado af´ın original, presentado por Dikin hacia finales de los 60 (Dikin (1967)). Algunos textos actuales utilizan esta versi´on para introducir las t´ecnicas del escalado af´ın (por ejemplo, Bertsimas y Tsitsiklis (1997)). En esta subsecci´on presentaremos brevemente dicho m´etodo. En el resto de la presentaci´on, sin embargo, continuaremos tratando u ´nicamente la versi´on de paso largo. El m´etodo del escalado af´ın de paso corto, al igual que el de paso largo introducido con anterioridad, utiliza el cambio de variables x e = (X k )−1 x presentado en (8), obteniendo entonces el siguiente problema en el espacio de variables x e min e cT x e e x ex = b sujeto a Ae

(15)

x e≥0 , e se han definido en (9). Este problema ya fue denotado con antedonde los t´erminos e c y A e rioridad como (P). Sabemos que con este escalado el punto xk en el problema original queda

19

transformado en x ek = eT = (1 , . . . , 1)T en el nuevo espacio. Consideremos ahora la siguiente esfera e x, ρ) = {e ex = b}. B(e x | ||e x − e|| ≤ ρ, Ae (16) ex = b. Dado que e es factible para formada por la intersecci´on de ||e x − e|| ≤ ρ y el subespacio Ae e x, ρ) estar´a incluido dentro de la regi´on factible {Ae ex = b, x (15), tenemos que B(e e ≥ 0} de (15) e si ρ ≤ 1. El conjunto de puntos B(e x, ρ) puede ser transformado al espacio de variables original x deshaciendo el escalado mediante (8) y (9): B(x, ρ) = {x | ||(X k )−1 x − e|| ≤ ρ, AX k (X k )−1 x = b} = {x | ||(X k )−1 x − e|| ≤ ρ, Ax = b}.

(17)

Lo que en el espacio x e era una esfera se ha convertido en un elipsoide en el espacio original x, e x, ρ = 1) estaba incluido en la gracias a la matriz de escalado (X k )−1 . Adem´as, al igual que B(e regi´ on factible de (15), se verifica que el elipsoide B(x, ρ = 1) tambi´en est´a incluido en la regi´on factible de (P). Esto lo garantiza el siguiente resultado: Proposici´ on 2. El conjunto de puntos definidos por el elipsoide B(x, ρ) son factibles para el problema (P) si ρ ≤ 1. Esto es B(x, ρ ≤ 1) ⊆ {x | Ax = b, x ≥ 0}. Si ρ < 1, adem´as de factibles son puntos interiores. Demostraci´on. Por definici´on de B(x, ρ) sabemos que todos sus puntos x verifican que Ax = b. Debemos ver que adem´as x ≥ 0. Sabemos que si x ∈ B(x, ρ) entonces ||(X k )−1 x − e|| ≤ ρ. Utilizando que para cualquier vector v, vi ≤ ||v||, tenemos |

xi − 1| ≤ ||(X k )−1 x − e|| ≤ ρ ≤ 1. xki

Multiplicando la relaci´on anterior por xki podemos escribirla como |xi − xki | ≤ ρxki ≤ xki , obteniendo el resultado deseado: xi ≥ 0. Si ρ < 1 las relaciones anteriores se verifican como estrictas desigualdades, teniendo que xi > 0 y por tanto interior. La Fig. 8 muestra gr´aficamente la relaci´on entre la esfera (16) y el elipsoide (17) en sus respectivos espacios para ρ = 1. Se observa como el elipsoide es factible para (P), tal y como se ha demostrado en la Proposici´on 2. Como los puntos del elipsoide (17) son factibles para (P), el m´etodo del escalado af´ın de paso corto itera sustituyendo las restricciones x ≥ 0 por las de pertenencia a dicho elipsoide. En el problema escalado, en vez del elipsoide considera una esfera, planteando el siguiente problema de optimizaci´on para obtener ∆e x. min e cT (e xk + ∆e x) ∆e x e xk + ∆e sujeto a A(e x) = b

(18)

||∆e x||2 = ρ2 . e x, ρ), Notemos que la regi´on factible de este problema coincide con la frontera de la esfera B(e k dado que x e = e. Si consideramos ρ ≤ 1 dicha regi´on est´a adem´as incluida en la regi´on 20

Problema original

(x)

Problema escalado

( ~x )

x3 ~ x3

e

xk

x~* ~x 2

x* x2 x~1

x1

e (e Figura 8. Representaci´on de la esfera B x, ρ = 1) = {e x | ||e x − e|| ≤ e ρ, Ae x = b} en el problema escalado y el elipsoide B(x, ρ = 1) = {x | ||(X k )−1 x − e|| ≤ ρ, Ax = b} en el problema original. Se considera un problema con una u ´nica restricci´ on y tres variables. Se observa como los puntos del elipsoide son factibles para el problema original.

factible de (15), como hemos argumentado anteriormente. Podr´ıamos haber escrito la segunda restricci´on como ||∆e x||2 ≤ ρ2 , teniendo entonces exactamente como regi´on factible de (18) e x, ρ). Sin embargo el ´optimo del problema se hallar´a en la frontera, dado que la la esfera B(e funci´on objetivo es lineal, por lo que preferimos formular directamente la segunda restricci´on en forma de igualdad. En este problema estamos buscando la direcci´on ∆e x que nos proporciona la m´axima disminuci´on de nuestra funci´on objetivo en un radio de longitud ρ alrededor de x ek = e. Puede verse que es equivalente al problema (5) que planteamos en el apartado 5.2.3. Y sabemos por (7) que su soluci´on viene dada por ∆e x = −ρ

Pee c . e ||P e c||

El problema (18) podr´ıa haber sido planteado en el espacio de variables x directamente, deshaciendo el escalado, obteniendo entonces: min cT (xk + ∆x) ∆x

(19)

sujeto a A(xk + ∆x) = b k −1

||(X )

2

2

∆x|| = ρ .

La regi´on factible corresponde ahora a la frontera del elipsoide B(x, ρ), y sabemos por la Proposici´ on 2 que este elipsoide est´a incluido en la regi´on factible del problema original (P) si ρ ≤ 1, y que es interior a ´esta si ρ < 1 . La soluci´on de (19) se obtiene a partir de la de (18), deshaciendo el escalado: X k Pee c . ∆x = X k ∆e x = −ρ e ||P e c|| El nuevo punto ser´a entonces xk+1 = xk − ρ

21

X k Pee c . e ||P e c||

x* x2

x0 x

1

Figura 9. Secuencia de puntos generados por el algoritmo del escalado af´ın de paso corto. Se muestran los diferentes elipsoides considerados a cada iteraci´ on, y como el nuevo punto se halla sobre la frontera de los mismos. Con l´ınea discontinua se muestra, en dos casos, el movimiento que realizar´ıa la variante de paso largo, avanzando hasta encontrar alguna de las paredes de la regi´ on factible.

Para garantizar que sea interior, u ´nicamente es necesario utilizar un ρ < 1. En (10) vimos como desarrollar el t´ermino X k Pee c. Observemos que en el algoritmo del escalado af´ın de paso corto no se realiza el c´alculo de α, puesto que la longitud de paso viene determinada por ρ/||Pee c|| y nos conduce hasta la frontera del elipsoide. El t´ermino ρ juega un papel similar al mismo t´ermino del algoritmo de paso largo: nos acorta el paso realizado, aqu´ı alej´andonos de la frontera del elipsoide, y en el algoritmo de paso largo de alguna pared definida por x ≥ 0. La Fig. 9 muestra la secuencia de puntos que generar´ıa el m´etodo para un problema concreto. Los elipsoides se adaptan, segun el escalado, a la forma de la regi´on factible. Naturalmente, si en vez de quedarnos en la frontera del elipsoide continuamos avanzando hasta alcanzar alguna pared definida por x ≥ 0, tal y como se realiza en la variante de paso largo, disminuiremos nuestra funci´on objetivo. 5.6. Detecci´ on del ´ optimo (finalizaci´ on del algoritmo) Sabemos como calcular ∆x y α, pudiendo generar la sucesi´on de puntos {xk } a trav´es del proceso iterativo (1). Ahora debemos establecer un criterio para determinar cuando el punto actual xk est´a lo suficientemente pr´oximo a x∗ . Para ello utilizaremos el Teorema 2 (de la Dualidad Fuerte). Si somos capaces de obtener un estimador de la soluci´on dual y asociada con xk , entonces podremos detener nuestro proceso iterativo cuando el gap dual est´e por debajo de una determinada tolerancia ε. Nuestra condici´on para detener el algoritmo ser´a entonces: |cT xk − bT y| |gap dualk | = ≤ ε. 1 + |cT xk | 1 + |cT xk |

(20)

El denominador en la expresi´on anterior se incluye para considerar el gap dual relativo al valor de la funci´on objetivo, y se le suma 1 para evitar problemas de precisi´on cuando el coste ´optimo sea cercano a 0. En la pr´actica se utilizan valores ε ∈ [10−6 , 10−8 ]. Ahora debemos obtener un estimador razonable de y. Y para ello utilizaremos el Teorema 3 (de la Holgura Complementaria). Sabemos que las holguras duales z verifican z = c − AT y, z ≥

22

0, y que en el ´optimo de los problema primal y dual xi zi = 0, i = 1, . . . , n. Dejando a parte la restricci´on z ≥ 0, podemos estimar y mediante el valor proporcionado por el siguiente problema 1 ||X k z||2 2 sujeto a z = c − AT y. min z,y

Introduciendo la funci´on Lagrangiana L(z, y, u) =

1 ||X k z||2 − uT (z − c + AT y), 2

donde los multiplicadores u ∈ IRn , e igualando a 0 sus derivadas respecto z, y y u, obtenemos las condiciones que debe satisfacer el punto soluci´on: (X k )2 z − u = 0 Au = 0 z − c + AT y = 0. Multiplicando la primera ecuaci´on por A, y utilizando que Au = 0, obtenemos A(X k )2 z = Au = 0. Multiplicando la tercera ecuaci´on por A(X k )2 y utilizando el resultado anterior obtenemos ) 0 = A(X k )2 z − A(X k )2 c + A(X k )2 AT y ⇒ y = (A(X k )2 AT )−1 A(X k )2 c. (21) = −A(X k )2 c + A(X k )2 AT y Esta ser´a nuestra estimaci´on de y. An´alogamente, el vector z = c − AT y corresponder´ıa a la estimaci´on de las holguras duales. Haciendo uso de la notaci´on D = (X k )2 introducida en apartados anteriores, la estimaci´on y puede escribirse como y = (ADAT )−1 ADc. Observemos que coincide con el vector y definido en (12), que interviene en el c´alculo de la direcci´on de movimiento ∆x. Por tanto no es necesario hacer ninguna operaci´on adicional para obtener dicho estimador. Tambi´en es interesante notar que la expresi´on (21) coincide con el vector de multiplicadores (6) de los problemas de optimizaci´on (5), (18) y (19), que nos proporcionaban la direcci´on de movimiento a cada xk . Por tanto, nuestra secuencia de estimaciones duales no es m´as que la secuencia de multiplicadores de las ecuaciones Ax = b de estos problemas. En el ejercicio 6 se pide comprobar que esta estimaci´on de las variables duales equivale a las del m´etodo del s´ımplex a medida que nos acercamos a un punto extremo de la regi´on factible (en este caso, sin embargo, nos acercamos por puntos interiores, en vez de por la frontera como en el s´ımplex). Dado que consideramos que el problema no es degenerado, la Proposici´on 1 nos garantiza que A(X k )2 AT es no singular, con lo cual el c´alculo de los estimadores duales y est´a bien definido. Si, por contra, el problema fuese degenerado, podemos hallarnos con una matriz (A(X k )2 AT ) casi singular, y si el punto soluci´on x∗ es degenerado, esta situaci´on ir´a a peor a medida que iteremos (en la subsecci´on 5.9 se presenta muy brevemente detalles sobre la convergencia del m´etodo y se comentan los problemas derivados por la degeneraci´on). En este caso las estimaciones de y obtenidas por (21) no ser´an buenas (pueden no converger al y ∗ de (D)), y nuestro criterio de parada (20) poco acertado. Un criterio alternativo que puede 23

usarse en estas situaciones consiste en parar cuando la mejora relativa en la funci´on objetivo sea peque˜ na, esto es |cT xk − cT xk+1 | ≤ ε. (22) 1 + |cT xk | 5.7. Obtenci´ on de una soluci´ on inicial factible Sabemos como iterar y como finalizar el m´etodo. Sin embargo hemos supuesto que part´ıamos de un punto inicial x0 factible e interior. Hallar tal punto para (P) no es trivial. A continuaci´on presentamos un m´etodo para solventar esta dificultad, denominado m´etodo Big-M (equivalente a la t´ecnica del mismo nombre para el algoritmo del s´ımplex). Consideremos un punto x0 > 0 cualquiera (por ejemplo, x0 = (1, . . . , 1)T ). Este punto escogido al azar no verificar´a probablemente Ax0 = b. Podemos calcular su vector de infactibilidades como r = b − Ax0 r ∈ IRm . Ahora, utilizando un valor M ∈ IR, donde M À 0, definamos el siguiente problema alternativo a (P) · ¸ x min [cT M ] x xn+1 · ¸ x sujeto a [A r] =b (23) xn+1 · ¸ x ≥0 . xn+1 Este nuevo problema tiene una variable m´as que (P), y gracias a ello podemos obtener un punto factible e interior de forma sencilla. Este punto ser´a  x0  1

 x02   .   . .  .   0 xn 1

(24)

Por ser x0 interior, (24) tambi´en lo es. Y utilizando la definici´on del vector r comprobamos su factibilidad para (23): · 0¸ x [A r] = Ax0 + r = Ax0 + (b − Ax0 ) = b. 1 Como M es un valor grande, es de esperar que en el punto ´optimo de (23) se verifique = 0. Cuando esto no sea as´ı concluiremos que el problema original (P) es infactible, puesto que no hemos podido eliminar la componente asociada con el vector de infactibilidades r. Los principales inconvenientes de este m´etodo son: x∗n+1

Debemos obtener un valor M adecuado, ni demasiado grande, ni demasiado peque˜ no. Si es muy grande podemos introducir inestabilidad num´erica en los c´alculos. Si es muy peque˜ no, penalizaremos poco la variable asociada al vector r, con lo que podemos obtener soluciones con x∗n+1 6= 0 sin ser el problema infactible. En la literatura se han propuesto t´ecnicas de adaptaci´on de M a medida que avanza el algoritmo. La columna r introducida en la matriz de restricciones de (23) por lo general ser´a densa. Esto provoca que la matriz ADAT sea tambi´en densa. La factorizaci´on de ADAT a cada

24

iteraci´on del algoritmo es el paso m´as costoso a realizar, desde un punto de vista computacional. El hecho que esta matriz sea densa incrementa considerablemente el tiempo de ejecuci´on, especialmente en la resoluci´on de problemas de gran dimensi´on. Entre sus ventajas destacamos que proporciona el ´optimo solucionando un u ´nico problema. En esto se diferencia de m´etodos alternativos que realizan la optimizaci´on en dos fases. 5.8. Algoritmo del escalado af´ın Estamos ahora en disposici´on de presentar todas las etapas del algoritmo del escalado af´ın ´ primal. Estas se presentan en el Algoritmo 1, el cual no requiere pr´acticamente comentarios. En la zona de inicializaciones se amplia el problema seg´ un el m´etodo Big-M. El punto inicial de iteraci´on utilizado —correspondiente a (24)— tiene todas las componentes a 1 (l´ınea 4 del algoritmo). En las inicializaciones tambi´en se obtiene la estimaci´on dual de x0 (l´ınea 5). En la zona iterativa se comprueba en primer lugar la condici´on de finalizaci´on del algoritmo (l´ınea 7), para pasar a calcular el nuevo punto (l´ıneas 8–12). A continuaci´on se estiman las variables duales para el nuevo punto (l´ıneas 14–16), y se vuelve a iterar. Al finalizar el proceso iterativo se comprueba si el problema es infactible o si hemos hallado su ´optimo, seg´ un el valor de la componente de x asociada con la columna de infactibilidades. En el ejercicio 7 se pide programar el Algoritmo 1 en un lenguaje de c´alculo matricial como Matlab. Veamos un ejemplo de aplicaci´on del Algoritmo 1 a la resoluci´on de un problema lineal de peque˜ na dimensi´on. EJEMPLO 6

Solucionemos el problema lineal min

− 3x1 − 2x2

sujeto a 4x1 3x1 x1

−2x2 +4x2 + x2

+x3

= = =

− x4 + x5

5 1 2

xi ≥ 0 para toda i = 1, . . . , 5 aplicando el algoritmo del escalado af´ın primal. Este problema ya fue utilizado en los ejemplos 1 y 4, y sabemos que x0 = (1/2 1/2 4 5/2 1)T es un punto interior factible. Usaremos este punto evitando tener que aplicar el m´etodo Big-M. Utilizaremos un valor de ρ = 0.95. La primera iteraci´ on del algoritmo ser´ıa como sigue: La funci´ on objetivo en el punto actual es cT x0 = −2.5. Calculamos y solucionando (ADAT )y = ADc, donde D = (X k )2 , X k = diag(0.5, 0.5, 4, 2.5, 1).

Ã

21 1 0.5

1 12.5 1.75

0.5 1.75 1.5

!

à y=

−2 −4.25 −1.25

!

à =⇒

y=

−0.070718 −0.264115 −0.501627

!

Comprobamos que el punto actual no es ´ optimo. gap dual =

|cT x0 − bT y| = | − 2.5 − (−1.620)|/3.5 = 0.25114. 1 + |cT x0 |

Calculamos z = c − AT y:

 −3 

 −1.576842 

 −1.423158 

 −2   −1.416651   −0.583349       z = c − AT y =   0  −  −0.070718  =  0.070718  . 0 0.264115 −0.264115 0 −0.501627 0.501627 25

.

Algoritmo 1. Algoritmo del escalado af´ın

1 2 3 4 5 6

INICIALIZACIONES Calcular infactibilidades: r = b − Aen , donde el = (11 , . . . , 1l )T Ampliar A con la columna adicional r: A ← [A r], n ← n + 1 Ampliar vector de costes: c ← [cT M ]T , M ∈ IR grande x0 = en+1 ; x0 es interior (x0 > 0) y factible (Ax0 = b) D = I, y = (ADAT )−1 ADc k = 0, ε = 10−6 , ρ ∈ [0.95, 0.9995] PROCESO ITERATIVO

7 8 9 10 11 12 13 14 15 16 17

|cT xk − bT y| > ε hacer 1 + |cT xk | Calcular z: z = c − AT y Calcular ∆x: ∆x = −Dz si (∆x ≥ 0) entonces FIN: ilimitado ½ Problema ¾ xki Calcular α: α = ρ · min − ∀i ∆xi ≤ 0 ∆xi Calcular xk+1 : xk+1 = xk + α∆x k ←k+1 Calcular X k : X k = diag(xk1 , , . . . , xkn ) Calcular D: D = (X k )2 Calcular y: (ADAT )y = ADc fin mientras mientras

´ FINALIZACION k 18 si xn+1 6= 0 entonces 19 FIN: Problema infactible 20 sino ´ 21 FIN: Optimo hallado: x∗ ← (xk1 , . . . , xkn )T 22 fin si Calculamos ∆x = −Dz:

 0.35579   0.14584   ∆x = −Dz =   −1.13148  . 1.65072 −0.50163

Determinamos α: α = ρ · min{−x0i /∆xi ∀i tal que ∆xi < 0} = 0.95 · min{4/1.13148, 1.0/0.50163} = 0.95 · 1.9935. Finalmente obtenemos el nuevo punto x1 = x0 + α∆x:

 0.5 

 0.35579   1.173803   0.5   0.14584   0.776190       x1 =   4  + 0.95 · 1.9935  −1.13148  =  1.857169  . 2.5 1

1.65072 −0.50163

26

5.626170 0.050007

La funci´ on objetivo en este nuevo punto es cT x1 = −5.0738 < cT x0 = −2.5. El vector de estimadores duales y en el nuevo punto x1 viene dado por y = ( −0.1361242

0.0017763

−2.4025868 )T .

El gap dual contin´ ua siendo elevado |cT x1 − bT y| = 0.067538, 1 + |cT x1 | por lo que deber´ıamos seguir iterando. La secuencia de puntos obtenida, junto con el gap dual, es la siguiente k

xk1

xk2

xk3

xk4

xk5

gap dual

0 1 2 3 4 5 6

0.5 1.173808 1.475381 1.487639 1.499064 1.4995 1.5000

0.5 0.776192 0.497191 0.510988 0.499914 0.50042 0.50000

4 1.857154 0.092858 0.071418 0.003570 0.002723 0.000103

2.5 5.626192 5.414905 5.506874 5.496851 5.500300 5.500000

1 0.050000 0.027429 0.001371 0.001020 0.000051 0.000038

0.25116 0.06753 0.01181 0.00235 0.00045 8.82 · 10−5 1.72 · 10−5 ,

y su representaci´ on gr´ afica se muestra en la Fig. 10. En ´esta, las tres u ´ltimas variables se toman como holguras, considerando las restricciones del problema como desigualdades. x2

2 rest. 1 rest. 2 rest. 3

1.75 1.5 1.25 1 x1 0.75 x3

x0

0.5

2

x

x6 5 x4 x

0.25 0

0

0.25

0.5

0.75

1

1.25

1.5

1.75

2 x1

Figura 10. Representaci´on gr´afica de la secuencia de puntos {xk } obtenidos al solucionar el problema del ejemplo 6.

5.9. Convergencia del algoritmo El an´alisis de convergencia del algoritmo del escalado af´ın es m´as complejo que el de otros m´etodos de punto interior m´as “elaborados” y eficientes, contrariamente a lo que se podr´ıa suponer. Aqu´ı u ´nicamente demostraremos su convergencia para el caso particular que estamos considerando, donde los problemas primal y dual son no degenerados (Hip´otesis 2 y 3). En las Proposiciones 3 y 4 se presenta dicha demostraci´on de convergencia. El resto de resultados u ´nicamente ser´an presentados, indicando donde pueden hallarse sus respectivas demostraciones.

27

En primer lugar vamos a mostrar que si la secuencia de puntos xk , y k y z k converge, lo hace a los puntos ´optimos de (P) y (D) (los t´erminos y k y z k son los estimadores de las variables duales definidos en (12), a˜ nadiendo ahora por conveniencia el supra´ındice k de la iteraci´on). Proposici´ on 3. Dado un problema no ilimitado (P) y su dual (D), si la secuencia de puntos xk , y k y z k generada por el m´etodo del escalado af´ın converge hacia x∗ , y ∗ y z ∗ , entonces x∗ es ´optimo para (P), y el par (y ∗ , z ∗ ) es ´optimo para (D). Demostraci´on. Consideremos dos puntos xk y xk+1 de la secuencia. Como ambos puntos son factibles tenemos que A(xk − xk+1 ) = b − b = 0. Utilizando este hecho y las definiciones de z k en (12), ∆x en (11), y del proceso iterativo (1), la disminuci´on de coste que se obtiene al pasar de xk a x+1 puede escribirse como ∆coste = cT (xk − xk+1 ) = cT (xk − xk+1 ) − (y k )T A(xk − xk+1 ) = (c − AT y k )T (xk − xk+1 ) = (z k )T (xk − xk+1 ) = −(z k )T ρα∆x = (z k )T ρα(X k )2 z k

(25)

= ρα||X k z k ||2 . Por la definici´on de α en (13), y de ∆x en (11), tenemos ½ ¾ ¾ ½ xk 1 k α = min − i para toda i tal que ∆xi ≤ 0 = min ≥ 0 para toda i tal que z i ∆xi xki zik 1 © k k ª. = (26) max xi zi para toda i tal que zik ≥ 0 Utilizando que para cualquier vector v, vi ≤ ||v||, podemos acotar inferiormente α: α=

max

©

xki zik

1 1 ª≥ . k k ||X z k || para toda i tal que zi ≥ 0

Sustituyendo esta u ´ltima expresi´on en (25) se tiene ∆coste = ρα||X k z k ||2 ≥ ρ

||X k z k ||2 = ρ||X k z k ||. ||X k z k ||

Sabemos que la sucesi´on de costes cT xk es decreciente y que est´a acotada inferiormente, por ser (P) ilimitado. Por tanto es una sucesi´on convergente, verific´andose que lim ∆coste = lim cT (xk − xk+1 ) = 0.

k→∞

k→∞

Combinando las dos expresiones anteriores se tiene que lim ρ||X k z k || = 0 ⇒ lim xki zik = 0 para todo i = 1, . . . , n.

k→∞

k→∞

(27)

Como xk converge hacia x∗ y z k hacia z ∗ tenemos que x∗i zi∗ = 0 para todo i = 1, . . . , n.

(28)

Ademas, por ser xk una secuencia de puntos factibles, Ax∗ = b.

28

(29)

Comprobemos ahora que z ∗ ≥ 0. Supongamos que existe un zi∗ < 0. Por tanto existe un ´ındice K tal que para todo k > K se verifica zik < 0. Escribiendo el proceso iterativo que genera la secuencia de puntos xk para la componente i-´esima como xk+1 = xik − α(xki )2 zik , i y considerando un k > K, conclu´ımos que xk+1 > xki > 0 (ya que α > 0, xki > 0, y hemos i supuesto que zik < 0). Por tanto x∗i zi∗ < 0 contradiciendo (28). Conclu´ımos entonces que z ∗ ≥ 0.

(30)

Como z ∗ = c − AT y ∗ (por (12)), por las condiciones (28), (29) y (30), y el Teorema 3 de la Holgura Complementaria, se concluye directamente que x∗ y el par (y ∗ , z ∗ ) son los ´optimos de (P) y (D) respectivamente. El resultado anterior supone que las secuencias xk , y k y z k convergen. La siguiente Proposici´on demuestra que, si se verifican las Hip´otesis 2 y 3 de no-degeneraci´on de los problemas primal y dual, las secuencias anteriores son, en efecto, convergentes. Proposici´ on 4. Si los problemas (P) y (D) no son degenerados, y (P) no es ilimitado, entonces para todo ρ < 1 la secuencias xk , y k generadas por el m´etodo del escalado af´ın son convergentes. Demostraci´on. La demostraci´on consta de tres partes. En la primera se muestra que la secuencia xk est´a acotada. En la segunda, que cualquier subsecuencia de xk converge hacia una soluci´on b´asica factible. En la tercera se demuestra que la soluci´on b´asica factible a la cual converge cualquier subsecuencia es u ´nica, garantiz´andose entonces que xk (y a su vez tambi´en y k ) es convergente. En primer lugar mostraremos que la secuencia xk est´a acotada. Dado que (P) no es ilimitado y que, al igual que (D), es no degenerado, tiene una u ´nica soluci´on ´optima. Sea x∗ dicha soluci´on. Consideremos ahora el siguiente problema: max x

sujeto a

n X

xi

i=1

Ax = b cT x ≤ cT x∗ x ≥ 0.

El problema anterior tiene una soluci´on ´optima finita, puesto que x∗ es el u ´nico punto factible. Esto implica que el siguiente problema tambi´en tiene una soluci´on ´optima finita: max x

sujeto a

n X

xi

i=1

Ax = b cT x ≤ cT x0 x ≥ 0,

donde x0 es el primer punto generado por el algoritmo del escalado af´ın (la implicaci´on anterior puede comprobarse intuitivamente utilizando propiedades de los pol´ıtopos convexos y observando que en el segundo problema u ´nicamente hemos desplazado el hiperplano de la restricci´on cT x ≤ cT x∗ , de forma que si el primer problema no era ilimitado, el segundo tampoco puede 29

serlo; una demostraci´on formal de esta afirmaci´on puede hallarse en el Teorema 4.14 del texto Bertsimas y Tsitsiklis (1997)). Por tanto el conjunto de puntos {x|Ax = b, cT x ≤ cT x0 , x ≥ 0} Pn est´a acotado (de lo contrario el m´aximo de i=1 xi no ser´ıa finito). Dado que la secuencia de puntos xk generada por el algoritmo del escalado af´ın pertenece al conjunto anterior (ya que todos ellos son factibles y verifican cT xk ≤ cT x0 ), concluimos que es una secuencia acotada. Por ser xk una secuencia acotada, y utilizando un resultado b´asico de an´alisis matem´atico, podemos concluir que existe una subsecuencia xkj que converge a un punto l´ımite x. Vamos a comprobar que este punto l´ımite es una soluci´on b´asica factible. Consideremos el conjunto de ´ındices B = {i|xi > 0}. Como (P) es no degenerado, la cardinalidad de B es mayor que m, el n´ umero de restricciones del problema. El Teorema fundamental de la Programaci´on Lineal nos garantiza que si existe una soluci´on factible, existe tambi´en una soluci´on b´asica factible, y las componentes de la soluci´on b´asica factible son un subconjunto de las componentes positivas de la soluci´on factible. Por tanto podemos escoger un conjunto de m ´ındices B1 ,. . . ,Bm ∈ B, tal que las columnas AB1 ,. . . ,ABm son linealmente independientes y la soluci´on b´asica asociada x b es factible, y adem´as no degenerada. Sea B = [AB1 | . . . |ABm ] la base asociada y cB = (cB1 , . . . , cBm ). Usando la definici´on de z k , z k = c − AT y k , y la ecuaci´on (27) de k la Proposici´on 3, tenemos que limj→∞ xi j (ci − ATi y kj ) = 0 para toda componente i. Las k componentes xBji , i = 1 . . . m, convergen hacia xBi y son positivas porque Bi ∈ B, i = 1 . . . m. Por tanto cB − By kj debe converger hacia 0 para garantizar (27). Puesto que B es no singular, se tiene que y kj converge hacia B −1 cB . Por tanto, para toda componente i, ci −ATi y kj converge hacia ci −ATi B −1 cB , valor que coincide con el coste reducido de la variable i-´esima en la soluci´on b´asica x b. La no-degeneraci´on del problema dual (D) nos garantiza que ci − ATi B −1 cB 6= 0 para toda variable i no b´asica (i 6∈ {B1 , . . . , Bm }). Por tanto, para las variables no b´asicas i 6∈ {B1 , . . . , Bm } se tiene que ci − ATi y kj converge a un valor diferente de 0, por lo que, usando una vez m´as (27), tenemos que xi = 0 para i 6∈ {B1 , . . . , Bm }. Por tanto concluimos que el conjunto de ´ındices de las variables con valor estrictamente positivo es B = {B1 , . . . , Bm }, y que x = x b, garantizando que el punto l´ımite x es una soluci´on b´asica factible. Hemos comprobado que xk tiene como m´ınimo un punto l´ımite, y que cada punto l´ımite es una soluci´on b´asica factible. Ahora comprobaremos que ´este es el u ´nico punto l´ımite que existe. Sea δ > 0 un valor real tal que cualquier variable b´asica en cualquier soluci´on b´asica factible es mayor o igual que δ. Este valor δ siempre existir´a ya que, por la no-degeneraci´on de (P), las variables b´asicas son estrictamente positivas. Sea ² = δ/3. Dado que cualquier punto l´ımite de xk es una soluci´on b´asica factible, existe un valor K tal que para todo k ≥ K, xk se halla a una distancia inferior a ² de una determinada soluci´on b´asica x e (esto es, ||xk − x e|| ≤ ², k ≥ K). Supongamos ahora que existen dos soluciones b´asicas diferentes x y x e tales que para k ≥ K, ||xk − x|| ≤ ²

||xk+1 − x e|| ≤ ².

y

(31)

Consideremos que la variable i es no b´asica en x (xi = 0) y b´asica en x e (e xi ≥ δ). Sustituyendo xi = 0 en la primera expresi´on de (31) (y usando que vi ≤ |vi | ≤ ||v|| para cualquier vector v) tenemos que xki ≤ ². (32) An´alogamente, sustituyendo x ei ≥ δ en la segunda expresi´on de (31), obtenemos que |xk+1 −x ei | ≤ ² ⇒ −² ≤ xk+1 −x ei ≤ ² ⇒ xk+1 ≥x ei − ² ≥ δ − ² = 2². i i i

(33)

Utilizando la ecuaci´on de c´alculo del nuevo iterado (1), la definici´on de ∆x en (11), y la expresi´on de la longitud de paso m´aximo (26), obtenemos ¡ xk+1 = xki + ρα(−(xki )2 zik ) = xki 1 − ρ i

max

©

xki zik

30

¢ 1 ª xki zik . k para toda i tal que zi ≥ 0

(34)

Usando un argumento similar al utilizado al final de la Proposici´on 3 se puede afirmar que zik ≥ 0 (de otra forma no se garantizar´ıa que xi = 0). A partir de (34), sabiendo por las hip´ otesis que ρ < 1, y utilizando (32) y (33) obtenemos finalmente xk+1 ≤ xki (1 + ρ) < 2xki ≤ 2² ≤ xk+1 , i i llegando a una contradicci´on. Por tanto s´olo existe un u ´nico punto l´ımite x y la secuencia xk converge a ´el. Tal y como se ha mostrado anteriormente, y k converge tambi´en a B −1 cB , la soluci´ on dual asociada a x. Sin las hip´otesis de no-degeneraci´on, tambi´en puede demostrarse la convergencia del algoritmo, aunque en este caso nos vemos obligados a reducir m´as la longitud de paso a trav´es de la disminuci´on de ρ. El siguiente resultado puede hallarse demostrado en Tsuchiya (1996): Proposici´ on 5. Si ρ ≤ 2/3, y (P) no es ilimitado, la secuencias xk , y k generadas por el m´etodo del escalado af´ın son convergentes (independientemente de la degeneraci´on o no de los problemas primal y dual). A continuaci´on apuntamos un resultado sobre el orden de convergencia del m´etodo. Dado que avanzamos en la direcci´on de un gradiente proyectado (con un escalado) no debemos esperar un orden de convergencia superior al lineal. El siguiente resultado, de nuevo demostrado en Tsuchiya (1996), ratifica esta afirmaci´on Proposici´ on 6. Si ρ ≤ 1, y (P) no es ilimitado, la secuencia cT xk generada por el m´etodo del escalado af´ın √ converge linealmente hacia cT x∗ , con una velocidad de convergencia no inferior a a 1−ρ/(2 n). x*

Figura 11. Trayectoria seguida por el m´etodo del escalado af´ın para alcanzar el punto ´ optimo x∗ desde diferentes puntos iniciales, considerando valores 0 < ρ ≈ 0.

Acabamos este apartado indicando que todav´ıa no est´a resuelto el problema de demostrar si este m´etodo de punto interior, al igual que otros, tiene un coste polin´omico (esto es, si obtiene una soluci´on en un n´ umero de iteraciones que est´a acotado por un polinomio funci´on del n´ umero de variables n). Existe la sospecha, sin embargo, de que esto no es as´ı. Tal argumento se basa en el hecho que, para valores de ρ infinitesimamente peque˜ nos la trayectoria seguida por el algoritmo del escalado af´ın consistir´ıa en una curva que une el punto inicial x0 con el punto soluci´ on x∗ . La Fig. 11 muestra el aspecto t´ıpico que tendr´ıan algunas de estas curvas partiendo desde puntos iniciales distintos aunque pr´oximos. La principal caracter´ıstica que tienen estas 31

trayectorias es que, cuando se acercan a una faceta de la regi´on factible se mantienen tangentes a ellas. Esto significa que cuando una variable xi tiende hacia 0, su componente en la direcci´on de movimiento tambi´en se hace 0: lim ∆xi = 0. xi →0

Puede comprobarse utilizando la definici´on de ∆x en (11). En estas situaciones la trayectoria realizada por el m´etodo del escalado af´ın es similar a la del algoritmo del s´ımplex, avanzando a trav´es de la frontera de la regi´on factible de v´ertice en v´ertice, tal y como muestra la Fig. 11. Esto hace suponer que, al igual que el s´ımplex, el escalado af´ın puede no tener una complejidad polin´omica. 5.10. Rendimiento computacional del algoritmo Se ha comparado el rendimiento computacional de una implementaci´on del escalado af´ın con el de una implementaci´on del algoritmo del s´ımplex, en la soluci´on de una bater´ıa est´andar de 90 problemas de programaci´on lineal, conocida como Netlib (Gay (1985)). La implementaci´on del algoritmo del s´ımplex utilizada ha sido la del paquete Minos (Murtagh y Saunders (1983)). Como implementaci´on del algoritmo del escalado af´ın se ha tomado la del sistema Lpabo (Park et al. (1997)). El sistema Lpabo se distribuye en c´odigo fuente, est´a totalmente escrito en C, y permite al usuario escoger entre diferentes algoritmos de punto interior, lo cual hace de ´el un c´odigo especialmente apropiado para tareas docentes. La Tabla 1 muestra los resultados obtenidos en la soluci´on de los 90 problemas utilizados. Para cada problema se indica su nombre (columna “Problema”), y el n´ umero de restricciones, variables y elementos diferentes de cero de la matriz A (columnas “m”, “n”, y “ne” respectivamente). Para los dos c´odigos utilizados, Lpabo y Minos, se muestra el n´ umero de iteraciones realizado (columnas “N. iter”) y el tiempo de CPU en segundos requerido para la soluci´on del problema (columna “CPU”). Todas las ejecuciones han sido realizadas sobre un ordenador personal con un procesador Pentium a 133MHz, utilizando el sistema operativo Linux. La Tabla 2 resume la informaci´on presentada en la Tabla 1. Para cada c´odigo se indica el n´ umero total de problemas que pudieron ser solucionados, y se realiza una peque˜ na estad´ıstica descriptiva donde se muestra la media, desviaci´on est´andar, valor m´aximo y valor m´ınimo del n´ umero de iteraciones y tiempo de CPU de los dos c´odigos. Para obtener estos estad´ısticos u ´nicamente se han utilizado los 85 problemas que pudieron ser solucionados por los dos c´odigos. Puede comprobarse como Lpabo no proporcion´o la soluci´on ´optima en 4 problemas (debido fundamentalmente, al uso del m´etodo Big-M), mientras que Minos tuvo un comportamiento m´ as robusto y u ´nicamente en un caso no obtuvo el punto ´optimo. Tambi´en se observa como Lpabo tuvo un mejor comportamiento que Minos en lo que a tiempo de ejecuci´on se refiere. La comparaci´on del n´ umero de iteraciones invertido por cada c´odigo destaca una de las propiedades de los m´etodos de punto interior: el n´ umero de iteraciones es muy inferior al requerido por el algoritmo del s´ımplex. Finalmente, en las Figuras 12 y 13 se muestra la evoluci´on del n´ umero de iteraciones de los dos algoritmos y de su cociente de tiempos de CPU (tiempo de CPU de Minos dividido por el de Lpabo), seg´ un el n´ umero de variables de cada problema, respectivamente. En la Fig. 12 se observa claramente una de las propiedades del algoritmo del escalado af´ın, y en general de los m´etodos de punto interior: el n´ umero de iteraciones realizado no aumenta con la dimensi´on del problema, y normalmente es un valor relativamente bajo (en pocas ocasiones se superan las 100 iteraciones). En la Fig. 13 se observa como Lpabo fue, en la mayor´ıa de problemas, m´as eficiente que Minos (en la figura se separan las iteraciones favorables a Lpabo y Minos mediante la l´ınea horizontal situada en el valor de abscisas 1, donde ambos c´odigos tienen id´entico rendimiento).

32

Se observa tambi´en una ligera tendencia positiva, a favor de Lpabo, en los problemas de mayor dimensi´on. Tabla 1. Rendimiento de Lpabo y Minos Problema

m

25fv47 80bau3b adlittle afiro agg agg2 agg3 bandm beaconfd blend bnl1 bnl2 boeing1 boeing2 bore3d brandy capri cycle czprob d2q06c d6cube degen2 degen3 dfl001 e226 etamacro fffff800 finnis fit1d fit2d forplan ganges gfrd-pnc greenbea greenbeb

822 2263 57 28 489 517 517 306 174 75 644 2325 351 167 234 221 272 1904 930 2172 416 445 1504 6072 224 401 525 498 25 26 162 1310 617 2393 2393

Dimensiones n ne 1571 9799 97 32 163 302 302 472 262 83 1175 3489 384 143 315 249 353 2857 3523 5167 6184 534 1818 12230 282 688 854 614 1026 10500 421 1681 1092 5405 5405

11127 29063 465 88 2541 4515 4531 2659 3476 521 6129 16124 3865 1339 1525 2150 1786 21322 14173 35674 43888 4449 26230 41873 2767 2489 6235 2714 14430 138018 4916 7021 3467 31499 31499

Lpabo N. Iter. CPU

Minos N. Iter. CPU

41 59 25 19 33 38 78 32 27 24 51 52 40 31 47 26 37 53

7901 12230 145 15 167 247 310 505 118 101 1513 7164 692 194 175 354 288 2818 2001 58928 106649 1331 8040

57 31 39 158 34 37 67 45 33 32 45 36 49 147 70

33

10.6 23.3 0.1 0.0 1.8 4.7 9.1 0.7 0.2 0.1 4.3 44.3 1.5 0.4 0.3 0.5 0.8 16.5 — — 52.3 3.6 61.5 19668.3 1.1 3.6 6.0 2.0 4.5 41.5 2.4 10.9 1.4 71.6 34.6

91.0 267.6 0.4 0.1 2.0 2.9 3.1 3.2 1.4 0.4 12.6 151.9 3.7 0.9 1.3 1.9 1.8 59.2 25.3 1078.4 1288.9 8.9 194.7 —

618 849 290 609 3125 36532 347 782 719 31008 16309

2.7 4.7 4.1 4.2 8.3 336.2 2.5 12.9 7.0 902.4 516.1

Tabla 1. (cont.) Rendimiento de Lpabo y Minos Problema grow15 grow22 grow7 israel kb2 lotfi maros modszk1 nesm perold pilot pilot.ja pilot.we pilot4 pilot87 pilotnov recipe sc105 sc205 sc50a sc50b scagr25 scagr7 scfxm1 scfxm2 scfxm3 scorpion scrs8 scsd1 scsd6 scsd8 sctap1 sctap2 sctap3 seba

Dimensiones m n ne 301 441 141 175 44 154 847 688 663 626 1442 941 723 411 2031 976 92 106 206 51 51 472 130 331 661 991 389 491 78 148 398 301 1091 1481 516

645 946 301 142 41 308 1443 1620 2923 1376 3652 1988 2789 1000 4883 2172 180 103 203 48 48 500 140 457 914 1371 358 1169 760 1350 2750 480 1880 2480 1028

5665 8318 2633 2358 291 1086 10006 4158 13988 6026 43220 14706 9218 5145 73804 13129 752 281 552 131 119 2029 553 2612 5229 7846 1708 4029 3148 5666 11334 2052 8124 10734 4874

Lpabo N. Iter. CPU

Minos N. Iter. CPU

22 24 22 76 30 36 51

555 895 230 328 62 284 2524 1269 3319 4807 18709 8531 6357 1753 28105 2830 29 44 94 35 15 407 125 482 890 1369 153 967 549 1315 3470 304 812 1123 397

1.9 3.1 0.9 11.2 0.2 0.5 7.0 —

38 107 64 350 55 79 58 42 19 20 20 21 18 28 22 38 42 43 25 43 28 27 25 32 34 35 34

34

12.8 34.5 444.8 220.2 9.6 14.9 1422.7 31.2 0.2 0.2 0.3 0.1 0.1 0.8 0.2 1.5 3.4 5.7 0.3 2.5 0.8 1.4 2.5 0.9 4.9 6.4 47.6

6.1 12.0 2.0 1.7 0.2 1.2 28.4 12.2 31.2 43.7 794.1 107.2 62.1 13.2 1876.3 39.8 0.5 0.4 0.8 0.2 0.2 3.3 0.6 2.8 7.8 15.7 1.7 7.0 2.5 6.8 31.9 2.0 10.9 18.1 4.2

Tabla 1. (cont.) Rendimiento de Lpabo y Minos Problema share1b share2b shell ship04l ship04s ship08l ship08s ship12l ship12s sierra stair standata standgub standmps stocfor1 stocfor2 tuff vtp.base wood1p woodw

Dimensiones m n ne 118 97 537 403 403 779 779 1152 1152 1228 357 360 362 468 118 2158 334 199 245 1099

225 79 1775 2118 1458 4283 2387 5427 2763 2036 467 1075 1184 1075 111 2031 587 203 2594 8405

1182 730 4900 8450 5810 17085 9501 21597 10941 9252 3857 3038 3147 3686 474 9492 4523 914 70216 37478

Lpabo N. Iter. CPU

Minos N. Iter. CPU

156 24 39 34 32 40 35 37 34 82 27 55 55 75 27 44 27 20

285 132 302 388 209 617 331 1212 561 1378 726 171 116 460 113 3399 771 172 1081 5375

1.7 0.3 1.9 2.0 1.3 4.1 1.9 4.9 2.0 10.6 6.1 2.1 2.0 3.4 0.2 8.8 2.4 0.2 —

47

19.5

1.0 0.5 4.8 5.9 4.0 13.5 7.6 26.6 11.9 18.2 6.4 2.6 2.6 4.0 0.6 68.3 4.6 0.9 39.3 118.4

Tabla 2. Resumen de los resultados obtenidos con Lpabo y Minos Lpabo problemas solucionados media desviaci´on est´andar m´aximo m´ınimo

86

89

N. Iter

CPU

46.3 41.1 350 18

32.8 161.9 1422.7 0.04

35

Minos N. Iter

CPU

4126.3 86.3 13009.0 279.5 106649 1876.3 15 0.12

N. iteraciones (escala log)

100000

Minos Lpabo

10000

1000

100

10

100

1000 n (escala log)

10000

Figura 12. N´umero de iteraciones requerido por Lpabo y Minos en la soluci´ on de los problemas Netlib, seg´ un el n´ umero de variables de cada problema. Ambos ejes est´ an en escala logar´ıtmica. 25 cociente tiempos de CPU Minos y Lpabo

Minos vs Lpabo 20

15

10

5

0 100

1000 n (escala log)

10000

Figura 13. Cociente (tiempo de CPU de Minos)/(tiempo de CPU de Lpabo), seg´ un el n´ umero de variables de cada problema. El eje horizantal est´ a en escala logar´ıtmica.

6. Ejercicios propuestos 1. Deducid la expresi´on de la matriz de proyecci´on ortogonal P = In − AT (AAT )−1 A sobre el espacio nulo de A ∈ IRm×n , que fue definida en (2). Para ello utilizad que un vector cualquiera d ∈ IRn puede descomponerse en d = u + v, donde u pertenece al espacio nulo de A, y v al espacio de rango, tal y como muestra la Fig. 14. ´ SOLUCION

Como u y v pertenecen al espacio nulo y de rango de A respectivamente tenemos que AT y = v,

Au = 0

donde y ∈ IRm . Podemos sustituir el valor anterior de v en d = u + v, obteniendo u = d − AT y.

36

d

v

u

Figura 14. Descomposici´on del vector d como suma de los vectores ortogonales u y v. Multiplicando por A en ambos lados de la expresi´ on anterior, considerando que A es de rango completo, y utilizando que Au = 0, obtenemos y: 0 = Au = Ad − AAT y

⇒

y = (AAT )−1 d.

Este valor de y lo sustitu´ımos ahora en la expresi´ on u = d − AT y: u = d − AT y = d − AT (AAT )−1 d = (In − AT (AAT )−1 A)d = P d. Por tanto P nos proporciona la proyecci´ on ortogonal u de d sobre el subespacio nulo de A. Puede realizarse una deducci´ on equivalente, observando la Fig. 14 y planteando el problema de obtenci´ on del vector v de norma m´ınima de forma que u pertenezca al espacio nulo de A, esto es: 1 1 min ||v||2 = ||d − u||2 2 2 sujeto a Au = 0. Planteando la funci´ on Lagrangiana de este problema L(u, y) =

1 (d − u)T (d − u) − y T Au, 2

y deriv´ andola respecto u e y se obtienen las condiciones AT y = d − u Au = 0 , que coinciden con las obtenidas por el primer m´etodo.

2. Comprobad que la matriz de proyecci´on ortogonal P = In − AT (AAT )−1 A definida en (2) verifica las tres propiedades de (3): AP = 0, es sim´etrica, y es idempotente. ´ SOLUCION

Para ver que AP = 0 basta sustituir la definici´ on de P : AP = A(In − AT (AAT )−1 A) = A − (AAT )(AAT )−1 A = A − A = 0. La simetr´ıa se comprueba aplicando propiedades b´ asicas de c´ alculo matricial: P T = (In − AT (AAT )−1 A)T = InT − (AT (AAT )−1 A)T = In − (AT ((AAT )−1 )T (AT )T = In − (AT ((AAT )T )−1 A) = In − AT (AAT )−1 A = P. La idempotencia se observa simplemente elevando al cuadrado P : P 2 = (In − AT (AAT )−1 A) · (In − AT (AAT )−1 A) = In − AT (AAT )−1 A − AT (AAT )−1 A + AT (AAT )−1 (AAT )(AAT )−1 A = In − AT (AAT )−1 A = P.

37

3. Supongamos que la direcci´on de movimiento ∆x = −P c introducida en (4) verifica que ||P c||2 = 0. Comprobad que en estas circunstancias cualquier punto factible es ´optimo. ´ SOLUCION

Que P c sea igual a 0 significa que la proyecci´ on de c sobre el espacio nulo de A ´es 0. O lo que es lo mismo, que c pertenece al espacio de rango o imagen de A, por lo que existe un vector y ∈ IRm (considerando que A ∈ IRm×n ) tal que cT = y T A (esto es, el vector de costes c ´es combinaci´ on lineal de las filas de la matriz de restricciones A). Consideremos ahora un punto factible cualquiera. Postmultiplicando la ecuaci´ on anterior por este punto obtenemos cT x = y T Ax = y T b. En la expresi´ on anterior hemos utilizado la factibilidad de x, Ax = b. Por tanto, la funci´ on objetivo es constante en la regi´ on factible, concluy´endose que cualquier punto factible es o ´ptimo. Puede darse una justificaci´ on alternativa observando que y es un punto factible para el problema dual (D) max y T b sujeto a AT y +z = cT , z ≥ 0, con z = 0. Dado un x factible primal cualquiera, hemos visto que y T b = cT x. Por el Teorema 2 de la Dualidad Fuerte podemos concluir que y y x son ´ optimos de sus respectivos problemas dual y primal. Como que no hemos impuesto ninguna condici´ on especial sobre x (´ unicamente que sea factible), se tiene que cualquier x factible es ´ optimo. Podemos observar esta situaci´ on en el problema min sujeto a

3x1 + 4x2 x1 + x2 + x3 = 1 2x1 + 3x2 − x3 = 1 x1 ≥ 0

x2 ≥ 0

x3 ≥ 0,

donde y = ( 1 1 )T . En este caso se puede comprobar que cualquier punto factible verifica que 3x1 + 4x2 = 2.

4. Comprobad que la direcci´on de movimiento ∆x = −D(In − AT (ADAT )−1 AD)c definida en (11), donde D = (X k )2 , preserva la factibilidad del nuevo punto y es de descenso, esto es 1) A∆x = 0. 2) cT ∆x ≤ 0. ´ SOLUCION

1) Directamente multiplicando comprobamos la factibilidad: A∆x = −AD(In − AT (ADAT )−1 AD)c = −(ADIn − (ADAT )(ADAT )−1 AD)c = −(AD − AD)c = 0.

eec, donde P e es la matriz de proyecci´on ortogonal sobre 2) Sabiendo que ∆x = −X k P e e e es idempotente y sim´etrica (ver el A, que e c y A se definen como en (9), y que P

38

ejercicio 2), entonces:

eec cT ∆x = −cT X k P e Xkc = −cT X k P eT P e Xkc = −cT X k P e X k c||2 ≤ 0. = −||P

5. Mostrad que si en una iteraci´on del m´etodo del escalado af´ın se tiene una direcci´on ∆x ≥ 0 donde ∆x 6= 0, entonces el problema es ilimitado (esto es, podemos disminuir el valor de la funci´on objetivo tanto como queramos). ´ SOLUCION

Sabemos que, si ∆x 6= 0, la direcci´ on ∆x calculada es siempre de descenso, verific´ andose que cT ∆x < 0 (ver ejercicio 4). El nuevo punto se obtiene como xk+1 = xk + α∆x, y la mejora en nuestra funci´ on objetivo respecto el punto xk viene dada por αcT ∆x. Ahora bien, como ∆x ≥ 0 no podemos llevar ninguna componente de xk > 0 a 0, por lo que el valor de α (que es la longitud de paso que podemos realizar manteniendo la no-negatividad de xk ) ser´ıa infinito, y la reducci´ on de coste que obtendr´ıamos αcT ∆x tambi´en, concluyendo que el problema es ilimitado.

6. Considerad una partici´on de las columnas de A = [B N ] en una parte b´asica B ∈ IRm×m y otra no b´asica N ∈ IRm×(n−m) de forma an´aloga a como se realiza en el m´etodo del s´ımplex. Utilizando una partici´on equivalente para xk = [xTB xTN ]T y c = [cTB cTN ]T , y suponiendo que no existe degeneraci´on (esto es, que xB > 0), mostrad que el vector de estimadores duales y = (A(X k )2 AT )−1 A(X k )2 c definido en (21) verifica lim y = (B T )−1 cB ,

xN →0

siendo entonces equivalente con el vector de las variables duales obtenido en el m´etodo del s´ımplex. ´ SOLUCION

En primer lugar expresemos y seg´ un el particionamiento de matrices y vectores del enunciado: y = (A(X k )2 AT )−1 A(X k )2 c

µ

=

(B

µ

N)

2 XB

¶µ

2 XN

BT NT

¶¶−1

µ (B

N)

2 XB

¶µ 2 XN

cB cN

¶

2 2 2 2 = (BXB B T + N XN N T )−1 (BXB cB + N XN cN ). 2 A partir de la expresi´ on anterior, y utilizando que BXB B T es no singular por ser xB > 0, se obtiene el resultado deseado: 2 2 2 −1 −1 2 lim y = (BXB B T )−1 (BXB cB ) = (B T )−1 (XB ) B BXB cB = (B T )−1 cB .

xN →0

7. Programad el algoritmo del escalado af´ın tal y como se presenta en el Algoritmo 1 de la subsecci´on 5.8. Utilizad alg´ un lenguaje de alto nivel que permita operaciones de c´alculo matricial, como Matlab.

39

´ SOLUCION

La codificaci´ on del algoritmo en el lenguaje de Matlab es como sigue: [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56]

function xopt= eap(A,b,c,info) % uso: eap(A,b,c,info) % A es una matriz (m x n) % b es un vector columna de dimension m % c es un vector columna de dimension n % info indica nivel de informacion: % si info>0 se muestra informacion a cada iteracion % si info=0 se muestra la solucion obtenida unicamente % Halla la solucion del problema de programacion lineal % min c’x % suj.a Ax=b % x>=0 % aplicando el algoritmo del escalado afin (primal) % comprobamos parametro de informacion if (nargin==3) info=0; end%if % comprobamos que las dimensiones son correctas [m,n]= size(A); [mc,nc]= size(c); [mb,nb]= size(b); % comprobamos si c y b son vectores columna if ((nc ~= 1) | (nb ~= 1)) error(’Error: c y b deben ser vectores columna’); end%if % comprobamos si las dimensiones de b y c concuerdan con la matriz A if ((mc ~= n) | (mb ~= m)) error(’Error: dimensiones incompatibles para c, b y A’); end%if % comprobamos rango de A if (rank(A) ~= m) error(’Error: la matriz A no tiene rango completo’); end%if % parametros optgap= 1.0e-6; myeps= 1.0e-12; rho= 0.95; M= 100*(max(abs(c))); % calculamos punto factible inicial, ampliamos el problema e iteramos x= ones(n,1); A= [A b-A*x]; x= [x;1]; c= [c;M]; D= diag(x)^2; y= inv(A*D*A’)*A*D*c; i=0; dualgap= abs(c’*x-b’*y)/(1.0+abs(c’*x)); while ( dualgap > optgap )

40

[57] [58] [59] [60] [61] [62] [63] [64] [65] [66] [67] [68] [69] [70] [71] [72] [73] [74] [75] [76] [77] [78] [79] [80] [81] [82] [83] [84]

z= c-A’*y; dx= -D*z; if all(dx>myeps) x= inf*ones(n,1); fprintf(’Problema ilimitado\n’); break; end%if alph= min(-x(dx1.0e-4) fprintf(’Problema infactible\n’); xopt=nan*ones(1,n); else xopt=x(1:n); end%if %endfunction

Puede comprobarse que se ha utilizado un valor de ε = 10−6 para detectar la finalizaci´ on del algoritmo (l´ınea 44), de ρ = 0.95 (l´ınea 46), y una constante M = 100 · max{|ci | i = 1, . . . , n} como coste de la columna de infactibilidades del m´etodo Big-M (l´ınea 47).

8. Solucionad el problema lineal min 3x1 + x2 sujeto a 2x1 + x2 ≥ 2 3x1 + 4x2 ≤ 12 x1 ≥ 0

x2 ≥ 0.

aplicando el algoritmo del escalado af´ın. ´ SOLUCION

El punto soluci´ on es x∗ = (0 2)T . Se dejan los c´ alculos al lector.

7. Trabajo pr´ actico: resoluci´ on de un problema de ingenier´ıa mediante el escalado af´ın En esta secci´on sugerimos un ejercicio de car´acter pr´actico. Consiste en solucionar un problema de ingenier´ıa mediante una implementaci´on del Algoritmo 1 similar a la presentada en el ejercicio 7. El problema consiste en la optimizaci´on de un sistema turbo-generador, y constituye una simplificaci´on del presentado en Edgar y Himmelblau (1988). Las subsecciones 7.1 y 7.2 presentan el enunciado del ejercicio. La subsecci´on 7.3 muestra la soluci´on del mismo.

41

(Nota: la notaci´on utilizada en la presentaci´on del ejercicio para designar las variables y constantes del problema difiere de la del resto del documento. As´ı, por ejemplo, las variables se denotar´ an en min´ uscula y las constantes en may´ uscula; los t´erminos con sub´ındice no corresponder´an a componentes de vectores). 7.1. Presentaci´ on del problema: optimizaci´ on de un sistema turbo-generador Las t´ecnicas de Optimizaci´on son a menudo usadas para el dise˜ no y determinaci´on del r´egimen de funcionamiento de sistemas de vapor y turbo-generaci´on en la industria qu´ımica. Este tipo de sistemas funcionan con vapor de agua, el cual es utilizado para generar electricidad mediante el uso de turbinas. En este caso concreto intentaremos determinar los valores de las corrientes de vapor y flujos de potencia que minimizan los costes de un sistema de estas caracter´ısticas. La Fig. 15 muestra el esquema del sistema de vapor y turbo-generaci´on el´ectrica que ser´a considerado. El agua se calienta en la caldera, usando un determinado combustible. El vapor generado se incorpora a una corriente de vapor de alta presi´on de va Kg/h . Con este vapor de alta presi´on se alimentan dos turbinas (denotadas por turb. 1 y turb. 2 en el esquema). Los flujos de entrada a cada turbina son de i1 e i2 Kg/h, y con ´estos se genera una potencia el´ectrica de p1 y p2 Kw, respectivamente. En el caso que no pueda satisfacerse una determinada demanda de potencia el´ectrica mediante las dos turbinas ser´a necesario adquirir pe Kw de un proveedor externo. El vapor de salida de las turbinas, o1 y o2 , se incorpora a una corriente de vapor de media presi´on de vm Kg/h. Tambi´en es posible traspasar directamente vam Kg/h de vapor de la corriente de alta a la de media presi´on a trav´es de la v´alvulaam . combustible

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

caldera

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

agua

va (vapor a alta presi´on)

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

i1 v´alvulaam

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

vam ... .. ...... ....

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

turb. 1 ...... o1

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

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

i2

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

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

p1 ............ .........

turb. 2 ....

c o2

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

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

potencia generada ............ ...........

p2 pe (potencia externa)

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

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

vm (vapor a media presi´on) .... ....... ...........

Figura 15. Esquema del sistema turbo-generador. La turbina 1 es de doble extracci´on. Esto significa que realiza una doble transformaci´on de la energ´ıa mec´anica en energ´ıa el´ectrica, primero pasando el vapor de alta a media presi´on, y despu´es condensado este vapor de media presi´on, obteni´endose la corriente de c Kg/h de agua mostrada en la Fig. 15. La turbina 2, por su parte, es de extracci´on simple, y u ´nicamente tiene una etapa de generaci´on de energ´ıa el´ectrica (pasando de alta a media presi´on). Nuestro objetivo es satisfacer, con el coste m´ınimo, una determinada demanda de potencia el´ectrica generada y de vapor a media presi´on. En este problema el coste se debe u ´nicamente a la obtenci´on de la corriente de alta presi´on en la caldera y a la adquisici´on de potencia el´ectrica externa.

42

7.2. Modelizaci´ on del problema Una vez detallado el funcionamiento del proceso, pasamos a realizar la modelizaci´on de ´ la funci´on objetivo y de las diferentes restricciones que intervienen. Estas se han agrupado 6 categor´ıas: restricciones asociadas con la turbina 1 y con la turbina 2, balance de flujos del problema, balances de potencia del sistema, y demanda de energ´ıa y vapor a media presi´on. Funci´ on objetivo La funci´on objetivo a minimizar es en este caso lineal, y u ´nicamente debe contemplar la suma de los costes de producci´on de la corriente de vapor a alta presi´on y de adquisici´on de potencia externa. Supondremos que el coste de obtenci´on de un Kg de vapor a alta presi´on ya ha sido calculado previamente, teniendo en cuenta el precio del combustible y la eficiencia de la caldera. El coste, en pesetas por hora de funcionamiento del sistema, viene dado entonces por funci´on objetivo = Cva va + Cpe pe pts/h, donde el t´ermino Cva (pts/Kg) es el coste de generaci´on de un Kg de vapor a alta presi´on, y Cpe (pts/Kwh) es el coste de adquisici´on de un Kwh de energ´ıa. Restricciones asociadas con la turbina 1 Las restricciones que modelizan el funcionamento de la turbina 1 son: L´ımite sobre el flujo m´aximo de entrada i1 : i1 ≤ I1 . L´ımite sobre el condensado producido c: c ≤ C. L´ımite sobre la potencia p1 generada por la turbina: p1 ≤ P1 . Restricciones asociadas con la turbina 2 Las restricciones que determinan el comportamiento de la segunda turbina son: L´ımite sobre el flujo m´aximo de entrada i2 : i2 ≤ I2 . L´ımite sobre la potencia p2 generada por la turbina: p2 ≤ P2 . Restricciones de balance de flujos Las ecuaciones de balance de flujos nos permiten modelizar las relaciones que hay entre los diferentes flujos y corrientes de vapor del sistema, seg´ un el esquema de la Fig. 15. Las unidades de todos los flujos son Kg/h. Estas restricciones son: Balance de flujos general sobre el vapor de entrada a alta presi´on, el de salida a media presi´on, y el condensado obtenido: va = vm + c.

43

Balance de flujos en la turbina 1. El vapor de entrada debe ser igual a todo el vapor y condensado de salida: i1 = o1 + c. Balance de flujos en la turbina 2. De nuevo, el vapor de entrada debe ser igual a todo el vapor de salida: i2 = o2 . A la vista de esta ecuaci´on podemos eliminar una de las dos variables. Eliminaremos o2 y mantendremos i2 . Balance de flujos en la corriente de alta presi´on. Lo que entra en esta corriente (va ) debe ser igual a lo que sale: va = i1 + i2 + vam . Balance de flujos en la corriente de media presi´on. Como antes, todo lo que entra en esta corriente debe ser igual a todo lo que sale. vm = vam + o1 + o2 . Puede observarse que esta ecuaci´on es combinaci´on lineal de las anteriores, por lo que ser´a eliminada (garantizando que la matriz de restricciones tenga rango completo). Restricciones de balance de potencias Las ecuaciones de balance de potencias nos permitir´an ligar las potencias p1 y p2 generadas con los diferentes flujos de vapor de entrada y salida de cada turbina. Para obtener la potencia de una determinada corriente de vapor (expresada en Kg/h), s´olo debemos multiplicar su valor por una constante que nos proporciona la energ´ıa por cada Kg de vapor. Estas constantes se denotar´an por Ea , Em y Ec (en Kwh/Kg) para la corriente de alta presi´on, de media presi´on, y el condensado, respectivamente. Las ecuaciones de balance de potencias para cada turbina ser´an: Balance de potencias para la turbina 1: Ea i1 = Em o1 + Ec c + p1 . Balance de potencias para la turbina 2: Ea i2 = Em o2 + p2 . Teniendo en cuenta que i2 = o2 , tal y como hemos visto antes, podemos formular esta restricci´on como (Ea − Em )i2 = p2 . Restricciones de demanda del sistema Finalmente s´olo nos queda indicar cuanto vapor y potencia debe producir el sistema que estamos considerando. • Demanda de vapor generado a media presi´on. Este valor vm no debe ser inferior a su demanda Dvm : vm ≥ Dvm . • Demanda de potencia. La suma de la potencia generada y la externa adquirida no debe ser inferior a la potencia demandada Dp : p1 + p2 + pe ≥ Dp . 44

En la modelizaci´on anterior hemos obtenido un total de 10 variables y 12 restricciones. De ´estas, 5 eran de igualdad, 5 de ≤ y 2 de ≥. Todas las variables est´an acotadas inferiormente por 0, dado que representan potencias o flujos de vapor. A˜ nadiendo variables de holgura hj j = 1, . . . , 5 para las restricciones de ≤ y de exceso fj j = 1, 2 para las de ≥, el problema finalmente obtenido consta de 17 variables y 12 restriciones, y puede ser formulado en forma est´andar como sigue: min

Cva va + Cpe pe

sujeto a i1 + h1 = I1 c + h2 = C p1 + h3 = P1 i2 + h4 = I2 p2 + h5 = P2 va − vm − c = 0 i1 − o1 − c = 0 va − i1 − i2 − vam = 0 Ea i1 − Em o1 − Ec c − p1 = 0 (Ea − Em )i2 − p2 = 0 vm − f1 = Dvm p1 + p2 + pe − f2 = Dp

      

Turbina 1

) Turbina 2        )

Balance de flujos

(35)

Balance de potencias ) Demanda del sistema

va , vm , vam , i1 , i2 , o1 , c, p1 , p2 , pe ≥ 0 hj ≥ 0 j = 1, . . . , 5

f1 , f2 ≥ 0

El ejercicio consiste en solucionar el problema (35) mediante una implementaci´on del algoritmo del escalado af´ın similar a la del ejercicio 7. Para ello se utilizar´an los datos particulares presentados en la Tabla 3. Tabla 3. Datos del problema Coeficiente

Valor

Unidad

Cva Cpe I1 I2 C P1 P2 Dvm Dp Ea Em Ec

1 5 192000 150000 40000 12000 6000 15000 2000 0.8779 0.8185 0.1246

pts/Kg pts/Kwh Kg/h Kg/h Kg/h Kw Kw Kg/h Kw Kwh/Kg Kwh/Kg Kwh/Kg

45

7.3. Resoluci´ on del problema El problema (35) se ha solucionado utilizando la implementaci´on del Algoritmo 1 realizada en el ejercicio 7. Para mejorar el condicionamiento num´erico de los c´alculos (especialmente a la hora de factorizar la matriz ADAT ) se han escalado los t´erminos independientes I1 , C, P1 , I2 , I2 , Dvm y Dp , dividi´endolos por 1000 (por tanto los flujos de vapor son de 1000·Kg/h, las potencias 1000·Kw, y la funci´on objetivo est´a en 1000·pts/h durante el proceso de optimizaci´on). Hay que hacer notar que sin este escalado previo de los datos el programa abortaba por problemas de inestabilidad num´erica. Una implementaci´on m´as robusta, y realizada en un lenguaje de m´as bajo nivel, evitar´ıa este inconveniente. El algoritmo del escalado af´ın ha realizado 14 iteraciones para encontrar el punto ´optimo. La soluci´on proporcionada se muestra en la Tabla 4. Los datos se presentan ya en las unidades originales (Kw, Kg/h) tras haber multiplicado por 1000 el ´optimo obtenido. Se observa como hay 12 valores significativamente diferentes de 0, mientras que 5 de ellos son pr´acticamente nulos (si hubi´eramos utilizado el algoritmo del s´ımplex ser´ıan 0 y corresponder´ıan a las variables no b´asicas). El coste ´optimo es de 66475.99 pts/h. En la Fig. 16 se representa gr´aficamente la soluci´ on obtenida. Pueden observarse algunas discrepancias en los balances de flujo, motivadas por el hecho que ninguna variable, por ser interior, llega a tener el valor de 0. Tabla 4. Soluci´on ´optima proporcionada por el m´etodo del escalado af´ın primal Variable Valor ´optimo

Holgura

Valor ´optimo

va vm vam i1 i2

30929.92 15000.01 0.01 15929.92 14999.99

h1 h2 h3 h4 h5

176070.07 0.00 24070.09 135000.01 5109.21

o1

0.01

c

15929.91

Exceso

Valor ´optimo

p1 p2

12000.00 890.79

f1 f2

0.00 0.01

pe

7109.21

caldera

combustible

agua

30929.92 Kg/h

turb. 1

potencia generada

14999.99 Kg/h

15929.92 Kg/h

turb. 2

12000.0 Kw

890.79 Kw

7109.21 Kw

15929.91 Kg/h

0.01 Kg/h 0.01 Kg/h

14999.99 Kg/h

15000.01 Kg/h

Figura 16. Representaci´on de la soluci´on ´optima obtenida. 46

cT x - bT y

25 20 15 10 5 0 -5 -10 -15 -20 -25 2

4

6

8

10

12

14

k

Figura 17. Evoluci´on del gap dual durante la optimizaci´on del sistema turbogenerador por el m´etodo del escalado af´ın.

Finalmente, en la Fig. 17 se muestra la evoluci´on del gap dual cT x − bT y durante las 14 iteraciones realizadas por el algoritmo. Se observa como tiende hacia 0 a medida que nos acercamos al punto ´optimo. Tambi´en puede comprobarse como el punto (v´ertice) soluci´on ha sido identificado ya en la iteraci´on n´ umero 8. En el resto de iteraciones el algoritmo u ´nicamente realiza movimientos de aproximaci´on a este punto. 8. Aspectos del algoritmo no tratados Existen diversos aspectos del m´etodo que no han sido tratados y que podr´ıan incluirse en una presentaci´on m´as extensa. Entre estos destacamos los siguientes: La versi´on del m´etodo del escalado af´ın presentado soluciona el problema (P) en forma primal. Por este motivo se denomina normalmente m´etodo del escalado af´ın primal. Existe una versi´on, denominada escalado af´ın dual, que soluciona el problema dual. En este caso se itera con las variables duales y y z, y se realizan estimaciones de las variables primales x. La principal ventaja del m´etodo dual es que puede hallarse un punto inicial de forma directa para aquellos problemas donde las variables primales est´en acotadas superior e inferiormente, sin utilizar la t´ecnica Big-M. El principal inconveniente es que, al iterar en el espacio dual, proporciona u ´nicamente estimaciones de las variables x originales del problema. Para una descripci´on de este m´etodo puede consultarse el cap´ıtulo 6 de Arbel (1993). A la hora de realizar una implementaci´on eficiente existen diversos detalles a tener en cuenta. A cada iteraci´on debe factorizarse el sistema A(X k )2 AT , que es sim´etrico, definido positivo y tiene el mismo patr´on de esparsidad para toda X k . Esto permite que sea posible aplicar t´ecnicas espec´ıficas para sistemas sim´etricos y definidos positivos de gran dimensi´on. Para una presentaci´on de este tipo de t´ecnicas se podr´ıa utilizar el texto cl´asico de George y Liu (1981). El hecho de contar con una columna densa en A tambi´en debe ser tenido en cuenta, puesto que puede degradar completamente la esparsidad de A(X k )2 AT . Esto ocurre, por ejemplo, al a˜ nadir la nueva columna de infactibilidades en el m´etodo BigM. Entonces pueden aplicarse tanto t´ecnicas de descomposici´on del sistema en bloques (descomposici´on de Schur), como m´etodos iterativos de soluci´on de sistemas sim´etricos y definidos positivos (como el del gradiente conjugado). En el trabajo de Andersen et al. (1996) pueden obtenerse m´as detalles sobre este tema. El algoritmo del escalado af´ın (dual) ha sido utilizado con ´exito en la optimizaci´on de problemas de flujos en redes, solucionando el sistema ADAT de forma iterativa (mediante 47

un gradiente conjugado). Para una descripci´on de esta aplicaci´on concreta del escalado af´ın puede consultarse Resende y Veiga (1993). M´etodos alternativos de punto interior tambi´en han sido aplicados con ´exito a problemas de flujos en redes denominados multiart´ıculo —diferentes art´ıculos comparten la misma red—, solucionando el sistema de ecuaciones ADAT mediante un gradiente conjugado precondicionado. Pueden hallarse m´as detalles en Castro (1999). Los resultados de convergencia del algoritmo presentados en la subsecci´on 5.9 podr´ıan ampliarse. En particular, podr´ıa realizarse un an´alisis m´as profundo de las implicaciones de la degeneraci´on del problema y como ´esta afecta a la convergencia del m´etodo. El trabajo de Tsuchiya (1996) podr´ıa servir de base para una presentaci´on de estas caracter´ısticas. 9. Referencias [1] E.D. Andersen, J. Gondzio, C. M´esz´aros y X. Xu (1996), “Implementation of interior point methods for large scale linear programming”, Interior Point Methods of Mathematical Programming, ed. T. Terlaky, Kluwer Academic Publishers, The Netherlands. [2] A. Arbel (1993), Exploring Interior-Point Linear Programming: Algorithms and Software, The MIT Press, Cambridge, USA. [3] E.R. Barnes (1986), “A variation on Karmarkar’s algorithm for solving linear programming problems”, Mathematical Programming, 36, pp. 174–182. [4] D. Bertsimas y J.N. Tsitsiklis (1997), Introduction to Linear Optimization, Athena Scientific, Belmont MA. [5] D.P. Bertsekas (1995), Nonlinear Programming, Athena Scientific, Belmont, USA. [6] J. Castro (1999), “A specialized interior point algorithm for multicommodity flows”, SIAM Journal on Optimization (pendiente de publicaci´on). [7] I.I. Dikin (1967), “Iterative solution of problems of linear and quadratic programming”, Soviet Mathematics Doklady 8(3), pp. 674–675. [8] T.F. Edgar y D.M. Himmelblau (1988), Optimization of Chemical Processes, McGraw-Hill, New York. [9] D.M. Gay (1985), “Electronic mail distribution of linear programming test problems”. Mathematical Programming Society COAL Newsletter, 13 10–12. [10] J.A. George y J.W.H. Liu (1981), Computer Solution of Large Sparse Positive Definite Systems, Prentice-Hall, Englewood Cliffs, USA [11] B.A. Murtagh y M.A. Saunders. (1983), “MINOS 5.0. User’s guide”, Dept. of Operations Research, Stanford University, USA. [12] S.G. Nash y A. Sofer (1996), Linear and Nonlinear Programming, McGraw-Hill, Singapore. [13] S. Park et al. (1997), “Lpabo 5.2 User’s Manual”, Department of Industrial Engineering, Seoul National University. [14] M. Resende y G. Veiga (1993), “An implementation of the dual affine scaling algorithm for minimum cost flow on bipartite uncapacitated networks”, SIAM Journal of Optimization, 3, pp. 516–537. [15] T. Tsuchiya (1996), “Affine Scaling Algorithm” en Interior Point Methods of Mathematical Programming, ed. T. Terlaky, Kluwer Academic Publishers, The Netherlands.

48

[16] R.J. Vanderbei (1996), Linear Programming: Foundations and Extensions, Kluwer Academic Publishers, Boston, USA. [17] R.J. Vanderbei, M.S. Meketon y B.A. Freedman (1986), “A modification of Karmarkar’s linear programming algorithm”, Algorithmica, 1, pp. 395–407.

49