Scientia Et Technica ISSN: 0122-1701 [email protected] Universidad Tecnológica de Pereira Colombia

GARCÉS RUIZ, ALEJANDRO; TORO O., ELIANA MIRLEDY; GALVIS M., JUAN CARLOS MÉTODO DE PUNTOS INTERIORES APLICADO AL PROBLEMA DE TRANSPORTES Scientia Et Technica, vol. XI, núm. 27, abril, 2005, pp. 49-54 Universidad Tecnológica de Pereira Pereira, Colombia

Disponible en: http://www.redalyc.org/articulo.oa?id=84911698011

Cómo citar el artículo Número completo Más información del artículo Página de la revista en redalyc.org

Sistema de Información Científica Red de Revistas Científicas de América Latina, el Caribe, España y Portugal Proyecto académico sin fines de lucro, desarrollado bajo la iniciativa de acceso abierto

49

Scientia et Technica Año XI, No 27, Abril 2005. UTP. ISSN 0122-1701

MÉTODO DE PUNTOS INTERIORES APLICADO AL PROBLEMA DE TRANSPORTES RESUMEN El problema de transportes es un problema de programación lineal que tradicionalmente ha sido resuelto usando el algoritmo de transportes el cual aprovecha las características topológicas del problema. En los últimos años ha tomado vigencia una nueva forma de solucionar problemas de programación lineal: los métodos de puntos interiores. En este artículo se muestra la aplicación del método de puntos interiores al problema de transportes bajo una nueva formulación matemática que adapta el problema al método aprovechando las características de uno y otro. Adicionalmente se muestra una aplicación con base en datos estadísticos del transporte de carga en Colombia en el que compara las dos técnicas. PALABRAS CLAVES: Problema de transportes, Método de puntos interiores, programación lineal. ABSTRACT Transportation problem has traditionally been resolved as linear programming where topological characteristics are taking into account. A new called Interior Points technique has recently showed an alternative way to solve these kinds of problems. This paper shows how interior point technique can be applied to transportation problem. An adaptation of the mathematical model is proposed in order to apply this technique. Results comparing linear programming and interior point are also shown applied to Colombian load transport. KEYWORDS: programming. 1.

Transportation

problem,

interior

point

INTRODUCCIÓN

El problema de transporte es una estructura especial de programación lineal (PL) que puede resolverse por el método simplex [1]. Basado en las características especiales del problema se desarrolló un algoritmo llamado de transportes [2] que hace que el método simplex resulte ineficiente. Recientemente, ha cobrado vigencia la aplicación de métodos de puntos interiores para la solución de todo tipo de problemas de programación lineal y programación no lineal. Convencionalmente, se ha llevado problemas a una formulación estándar para aplicar una metodología determinada [3]. No obstante, el presente artículo plantea un nuevo enfoque en la solución de problemas de PL: adaptar el método al problema particular. Esto implica la reformulación del método de puntos interiores para explotar las características del problema de transportes. En la literatura sobre el tema[4] no hay una adaptación de la teoría de puntos interiores al problema de transportes, en este artículo se propone un desarrollo matemático al respecto. Inicialmente, se explica el problema de transportes y su importancia en aplicaciones prácticas. Posteriormente, se da una breve reseña de los métodos de puntos Fecha de Recepción: 31 Enero de 2005 Fecha de Aceptación: 16 Marzo de 2005

ALEJANDRO GARCÉS RUIZ. Profesor Catedrático Maestría en Ingeniería Eléctrica. Área de Planeamiento. [email protected] ELIANA MIRLEDY TORO O Profesor Catedrático Maestría en Ingeniería Eléctrica. Área de Investigación Operativa. [email protected] JUAN CARLOS GALVIS M. Profesor Catedrático Maestría en Ingeniería Eléctrica. Área de Planeamiento. [email protected] Grupo de Planeamiento en Sistemas eléctricos. Área de Investigación de Operaciones.

method, lineal

interiores y sus principales características en la solución de problemas lineales. Finalmente, se muestra la formulación matemática propuesta, sus ventajas y el algoritmo implementado en la solución de un problema de gran tamaño. Igualmente, se hace una comparación del algoritmo de transporte y el algoritmo de punto interior desarrollado. 2.

EL PROBLEMA DE TRANSPORTES.

El algoritmo de transporte sólo puede aplicarse cuando el sistema esta balanceado, es decir la oferta debe ser igual a la demanda, adicionalmente se debe tener una solución inicial, el método simplex la consigue a través de las variables de holgura o a través de las variables artificiales o empleando un método de dos fases, el algoritmo de transporte se inicializa por medio de métodos auxiliares como son: el método de la esquina noroeste, el método del costo mínimo o el método de aproximación de Voguel, siendo este último el más adecuado debido a que este nos garantiza un punto cercano al óptimo e incluso en algunos casos encuentra el óptimo. Se supone que Ns orígenes (S) tienen que surtir a Nd centros de consumo (D) con un cierto producto (Fig 1). La capacidad de oferta del origen i es Si i(i=1,...Ns) y la demanda en el centro de consumo j es Dj (j=1,...Nd). Se

Scientia et Technica Año XI, No 27, Abril 2005. UTP

50 supone que Tij es el costo de enviar una unidad del producto del origen i al centro de consumo j. El problema se reduce a determinar cuántas unidades del producto deben enviarse del origen i al centro de consumo j, tal que minimicen los costos totales de distribución, se satisfaga la demanda del centro de consumo j y no se exceda la capacidad de oferta del origen i. Tij Si

Dj Tim

. . .

. . .

Tnj Tnm

Sn

Dm

Sea Xij esta variable de decisión. Entonces la formulación del problema lineal es Ns Nd

∑∑ T i =1 j =1

ij

⋅ X ij

sa Nd

∑X j =1

ij

= Si , i = {1,.., Ns}∀i ≠ r

ij

= D j , j = {1,..., Nd }

Ns

∑X i =1

(1)

X ij ≥ 0 Ns: Número de nodos fuente. Nd: Número de nodos de demanda. Si: Cantidad a enviarse desde el nodo i. Dj: Cantidad a recibir en el nodo j. r : Nodo de referencia. Esta formulación se denomina estructura de transporte. Una condición necesaria y suficiente para que la estructura de transporte tenga solución es que la oferta total sea igual a la demanda total es decir:

∑

i =1

Si =

Nd

∑

j =1

Dj

MÉTODOS DE PUNTOS INTERIORES.

Los métodos de puntos interiores surgen como una alternativa eficiente a la solución de todo tipo de problemas en programación lineal, sus principales características son: - Los puntos no son necesariamente factibles. - No requiere de una fase de inicialización. - El camino de solución solo debe ser interior a las restricciones de desigualdad. - Se puede establecer la exactitud deseada. Esta última característica es especialmente importante ya que es posible controlar el tiempo de cálculo; el método puede encontrar una solución cercana a la óptima en pocas iteraciones.

en donde:

Ns

- El rango de A es Ns+Nd-1. Esto puede probarse fácilmente mostrando que la suma de los primeros renglones es igual a la suma de los últimos Nd renglones, y que cualquier submatriz cuadrada de A de orden Ns+Nd-1 es no singular. - La matriz A es unimodular, es decir, que cualquier submatriz cuadrada de A de orden Ns+Nd-1 tiene un determinante que es igual a 0, 1 ó -1 3.

Figura 1. El problema de transportes.

min

La estructura de la matriz de coeficientes tecnológicos de A de Ns+Nd renglones y Ns·Nd columnas permite el desarrollo del algoritmo de transporte y tiene las siguientes propiedades:

(2)

en donde Si representa la cantidad de producto generado por cada una de las fuentes (S), Dj representa la cantidad de producto solicitado por cada una de las demandas (D) y Tij es el costo de transportar el producto desde el nodo generador Si hasta el nodo de demanda Dj.

Los métodos de puntos interiores se basan en la aplicación de la barrera logarítmica, según iasco McCormick [5] la solución de un problema de programación lineal estándar como se plantea en la ecuación 4 puede ser encontrada mediante la solución sucesiva de sub-problemas usando una función de penalidad logarítmica como se muestra en la ecuación 5.

min cT ⋅ X sa A⋅ X = b

(4)

X ≥0 La secuencia de soluciones a los sub-problemas (Ec 5) tiende a la solución del problema primal (Ec 4) si los valores de µ(k) son elegidos tal que µ → 0 si cada punto encontrado es un punto interior. n

El problema de transporte puede escribirse en forma condensada como: Min Z = cX sujeto a AX = d X ≥0

max bT ⋅ Y + µ ( k ) ⋅ ∑ ln( Zi ) i =1

sa A ⋅Y + Z = c Zi > 0 T

(3)

(5)

Scientia et Technica Año XI, No 27, Abril 2005. U.T.P

51

Cabe anotar que Y representa las variables duales del problema (Ec 4) y Z las variables de holgura, así mismo, la ecuación 5 es una formulación dual a la cual se le ha agregado un factor logarítmico de penalidad, este factor introduce una no-linealidad al problema por lo cual debe ser resuelto por técnicas no-lineales (Método de Newton). La metodología de solución consiste en resolver la ecuación 5 para distintos valores de µ en donde µ debe ser disminuido con algún criterio que garantice la convergencia bajo un tiempo de cálculo apropiado. La función Lagrangiana para la ecuación 5 es la siguiente:

(

i =1

− X ⋅ A ⋅Y + Z − C T

Como se expuso anteriormente, el planteamiento tradicional de los problemas solucionados con puntos interiores consiste en formular un sistema estándar como el mostrado en la ecuación 3 y solucionarlo mediante la aplicación sucesiva del sistema 8. No obstante, en este artículo se muestra un planteamiento matemático adaptado al problema de transportes. 4.

n

L = bT ⋅ Y + µ (k ) ⋅ ∑ ln(Z i ) T

necesario chequear la factibilidad primal, la factibilidad dual y la cercanía entre las soluciones primal y dual. - El sistema puede ser solucionado de manera más eficiente si se consideran términos de segundo orden (Métodos de alta orden).

FORMULACIÓN MATEMÁTICA MODELO DE TRANSPORTES.

(6)

)

El modelo de transportes en su presentación primal es mostrado en la ecuación 1 El modelo dual puede ser planteado de la siguiente forma:

Las condiciones de Karush-Kuhn-Tucker (KKT) son:

∂L = 0 ∴ AT ⋅ Y + Z = C ∂X ∂L = 0 ∴ A⋅ X = b ∂Y ∂L = 0 ∴ X ⋅ Z = µ ⋅ eˆ ∂Z

AT 0

(7)

Rp = b − A ⋅ X

i =1

j =1

(9)

Yi + Y j ≤ Tij en donde: Yi: Variables duales correspondientes al nodo fuente i. Yj: Variables duales correspondiente al nodo de demanda j. Adicionando las variables de holgura se obtiene: Ns

Nd

i =1

j =1

max ∑ Si ⋅ Yi + ∑ D j ⋅ Y j (10)

sa Yi + Y j + Z ij = Tij

0 ⎤ ⎡ ∆X ⎤ ⎡ R p ⎤ I ⎥⎥ ⋅ ⎢⎢ ∆Y ⎥⎥ = ⎢⎢ Rd ⎥⎥ X ⎥⎦ ⎢⎣ ∆Z ⎥⎦ ⎢⎣ Rc ⎥⎦

Donde

Nd

sa

Para resolver el sistema (7) se procede a encontrar un sistema de newton linealizado:

0

Ns

max ∑ Si ⋅ Yi + ∑ D j ⋅ Y j

Un punto interior a un espacio de soluciones es aquel que cumpla con las restricciones de desigualdad aunque no cumpla con las restricciones de igualdad, por tal motivo deben ser respetadas las condiciones de no-negatividad ( X > 0 ) y ( Z > 0 ) para garantizar que el siguiente punto sea interior.

⎡A ⎢0 ⎢ ⎢⎣ Z

DEL

Usando el concepto de barrera logarítmica puede ser planteado la ecuación lagrangiana para la aplicación del método primal dual: (8)

Rc = µ ⋅ eˆ − X ⋅ Z ⋅ eˆ

Ns

Nd

Ns Nd

i =1

j =1

i =1 j =1

L = ∑ S i ⋅ Yi + ∑ D j ⋅ Y j +µ ⋅ ∑∑ ln( Z ij ) − ∑ ∑ X ij ⋅ (Y j + Yi + Z ij − Tij ) Ns Nd

(11)

i =1 j =1

Rd = c − A ⋅ Y − Z T

Este sistema tiene las siguientes características: - La matriz mostrada en (8), aunque dispersa, es de gran tamaño por lo que su inversa representa el mayor esfuerzo computacional del método. - El algoritmo de solución se acerca tanto por la formulación primal como por la dual, por tal razón es

Finalmente, derivando respecto a cada una de las variables se obtienen las condiciones de primer orden de Karush Kuhn Tucker (KKT): Condición Dual:

Scientia et Technica Año XI, No 27, Abril 2005. UTP

52

∂L = Yi + Y j + Z ij − Tij = 0 ∂X ij

Nd

∑F

(12)

Condición Primal

Ns

∂L = Si − ∑ X ij = 0 ∂Yi j =1

∑F

Nd

∂L = D j − ∑ X ij = 0 ∂Y j i =1 Ns

µ ∂L = − X ij = 0 ∂Z ij Z ij

(13)

ij

= RPSi

i =1

+ Wij ⋅ (∆Yi + ∆Y j ) = RPDj

(24)

Los sistemas de ecuaciones (23) y (24) pueden ser reescritos de la siguiente forma

(14)

Ns ⎛ Ns ⎞ ∆ Y ⋅ W + ∆ Y ⋅ W = R − Fij ⎜ ⎟ ∑ ∑ ij ⎠ PDj ∑ i ij j i =1 i =1 ⎝ i =1

(25)

Ns

(26)

Finalmente, una representación matricial de (25) y (26) es la siguiente:

[Wi ]⋅ [∆Yi ] + [Wij ]⋅ [∆Y j ] = [Gi ]

(15)

(27)

[W ] ⋅ [∆Y ] + [W ]⋅ [∆Y ] = [G ] ij

(16)

i

j

j

j

(28)

En donde:

NS

∑ ∆X

(23)

T

Nd

j =1

+ Wij ⋅ (∆Yi + ∆Y j ) = RPSi

Nd ⎛ Nd ⎞ Nd ∆Yi ⋅ ⎜⎜ ∑ Wij ⎟⎟ + ∑ ∆Y j ⋅ Wij = RPSi − ∑ Fij j =1 ⎝ j =1 ⎠ j =1

Este sistema no lineal de ecuaciones puede ser solucionado por medio de un algoritmo numérico convencional siempre que se tengan en cuenta las condiciones de no negatividad. Aplicando NewtonRaphson se tiene:

∆Yi + ∆Y j + ∆Z ij = RDij

ij

i =1

Condición de Complementariedad

∑ ∆X

ij

j =1

ij

= RPDj

Z ij ⋅ ∆X ij + X ij ⋅ ∆Z ij = RCij

[Wi] es una matriz diagonal de Ns x Ns (17)

⎛ Nd

⎞

⎝

⎠

[Wi ] = diag ⎜⎜ ∑Wij ⎟⎟

(18)

j =1

(29)

[Wj] es una matriz diagonal de Nd x Nd

donde:

RDij = Tij − Yi − Y j − Z ij

[W ] = diag ⎛⎜ ∑W Ns

j

Nd

RPSi = S i − ∑ X ij = 0 j =1

Ns

(19)

⎝ i =1

⎞ ⎟ ⎠

(30)

Y los vectores columna:

RPDj = D j − ∑ X ij = 0

Nd

i =1

[Gi ] = RPSi − ∑ Fij

(31)

[G ] = R

(32)

j =1

RCij = µ − X ij ⋅ Z ij De la ecuación (15) y la ecuación (18) se obtiene X como función de Y:

∆X ij = Fij + Wij ⋅ (∆Yi + ∆Y j )

ij

j

(20)

Ns

PDj

− ∑ Fij i =1

De las ecuaciones (27) y (28) se tiene

donde

Fij = Wij =

RCij − X ij ⋅ RDij Z ij X ij Z ij

Reemplazando en (16) y (17) se obtiene:

( [W ]− [W ]

T

j

(21)

[ ] )⋅ ∆Y

⋅ [Wi ] ⋅ Wij −1

[G ] − [W ] ⋅ [W ] ⋅ [G ] −1

T

j

(22)

ij

ij

i

j

= (33)

i

Las inversas de las matrices [Wi] y [Wj] tienen un bajo esfuerzo computacional al ser matrices diagonales. Una vez encontrado el vector [∆Yj] se puede hallar fácilmente el vector [∆Yi]:

Scientia et Technica Año XI, No 27, Abril 2005. U.T.P

[ ] [ ])

∆Yi = [Wi ] ⋅ ([Gi ] − Wij ⋅ ∆Y j −1

53 (34)

Una vez calculadas las variaciones en las variables duales (∆Y) pueden ser calculadas las de las variables primales (∆X) y las de las variables de holgura de la formulación dual (∆Z) usando las ecuaciones (20) y (15) respectivamente. 5. ACTUALIZACIÓN DE LAS VARIABLES La aplicación del método de Newton-Raphson permite la solución del problema planteado para un valor de µ . Sin embargo, para encontrar una solución óptima factible es necesario determinar un paso máximo que asegure que el nuevo punto sea un punto interior. Como se puede observar en las ecuaciones (1) y (5) el problema de transportes exige que X > 0 , Z > 0 para garantizar que un punto sea un punto interior, por tal razón el paso escogido en la actualización de las variables debe garantizar la no-negatividad de estas variables; para ello se escoge un valor de αp y αd así:

⎧⎪ min ⎧⎪ − Xij ⎫⎪⎫⎪ α p = min ⎨1, ⎨ ⎬⎬ ⎪⎩ ∆X ij < 0 ⎪⎩ ∆X ij ⎪⎭⎪⎭

(35)

⎧⎪ min ⎧⎪ − Zij ⎫⎪⎫⎪ α d = min ⎨1, ⎨ ⎬⎬ ⎪⎩ ∆Z ij < 0 ⎪⎩ ∆Z ij ⎪⎭⎪⎭

(36)

Las variables son entonces actualizadas de acuerdo con las ecuaciones 37 a 40:

Yi

( k +1)

= X ij

( k +1)

= Yi

( k +1)

= Yj

Yj

Z ij

( k +1)

(k )

(k )

= Zij

(k )

+ γ ⋅ α p ⋅ ∆X ij

+ γ ⋅ α d ⋅ ∆Yi

(k )

+ γ ⋅ α d ⋅ ∆Y j

(k )

(k )

( 0)

(37)

=

σ Tij

( 41)

En donde σ es un parámetro de centralidad mayor a uno, para este caso se usó 1x108. - Las variables duales (Yi,Yj) se inicializan en cero mientras que las variables de holgura son inicializadas de acuerdo a los costos de Transportes ( Zo = T ) 7. CRITERIO DE CONVERGENCIA Como el método hace una búsqueda en puntos interiores de la región primal y dual, la convergencia debe garantizar optimalidad primal y dual, para el caso del problema de transportes las ecuaciones son: •

Convergencia Primal:

Cp = •

•

Norm2 {R PS , RPD } 1 + Norm2 {X }

(42)

Convergencia Dual:

Norm2 {RD } (43) 1 + Norm2 {Yi, Yj} + Norm2 {Z }

Criterio de Optimalidad

Co =

P

1 + Norm 2 {Pd }

(44)

En donde P indica la diferencia entre la función objetivo del problema primal y dual: Nd ⎛ Ns Nd ⎞ ⎛ Ns ⎞ P = ⎜⎜ ∑∑ Tij ⋅ X ij ⎟⎟ − ⎜⎜ ∑ Si ⋅ Yi + ∑ D j ⋅ Y j ⎟⎟ (45) j =1 ⎝ i =1 j =1 ⎠ ⎝ i =1 ⎠

Se deben cumplir los tres criterios para finalizar el proceso iterativo.

(38)

(k )

+ γ ⋅ α d ⋅ ∆Z ij

X ij

Cd =

Con el fin de acelerar la convergencia, αp y αd tienen un rango de valores entre 0 y 1. Se ha observado que valores superiores a 1 hacen más lenta la convergencia del método de puntos interiores.

X ij

Las variables primales son inicializadas de tal forma que las rutas de menor costo tengan mayor cantidad de producto así:

(k )

(39) (40)

en donde γ es un parámetro que hace al método más conservador evitando que problemas de redondeo o de cercanía a las restricciones de igualdad puedan ocasionar divergencia. Para éste trabajo se usó un valor de 0.995 6. INICIALIZACIÓN DE LAS VARIABLES. Una de las mayores ventajas de los métodos de puntos interiores es que no requiere de un punto inicial factible, solo basta asegurar que sea un punto interior. En el caso particular del problema de transportes se tiene: - Variables primales Xo

8. REDUCCIÓN BARRERA.

DEL

PARÁMETRO

DE

Un paso importante para asegurar la convergencia del método consiste en encontrar una función adecuada para la reducción del parámetro de barrera µ. En el caso del problema de transportes se opta por una reducción acorde con las condiciones de complementariedad.

⎛ Ns Nd ⎞ ⎜ ∑∑ Zij ⋅ X ij ⎟ i =1 j =1 ⎠ P=⎝ Ns ⋅ Nd

(46)

µ ( k +1) = φ ⋅ P

(47)

Scientia et Technica Año XI, No 27, Abril 2005. UTP

54 En donde φ asume el valor de 0.1 si disminuye el valor de P y 2 en caso contrario. 9. EJEMPLO DE APLICACIÓN. Se tomó como base el informe de distribución de carga en toneladas por origen y destino elaborado por el ministerio de transporte colombiano, estadísticas de la subdirección operativa, y los costos por tonelada transportada entre los distintos puntos de origen y destino de una compañía que opera en Colombia. Se hizo el modelamiento a una estructura de transporte. Número de nodos fuente: 37 Número de nodos demanda 37 Intervalo de los valores de [320x103~15.9x109 ] generación de producto (s) kg Intervalo de los valores de [ 0~13.1x109 ] kg demanda (D) Intervalo de los Costos de [ 99 ~ 20x103 ] $/kg Transportes (T) Tabla 1. Características del ejemplo de prueba.

Método Punto Interior 20 1.133x1010

Nº Iteraciones Función Objetivo Tabla 2. Resultados del método.

Algoritmo Transportes 210 1.133x1010

Las configuraciones encontradas pueden ser distintas en caso de presentarse óptimos alternativos ya que mientras el algoritmo de transportes exige que la solución sea un vértice el método de puntos interiores puede encontrar soluciones sobre la línea de isoutilidad sin ser necesariamente un vértice. El algoritmo de transporte inicia desde un punto cercano al óptimo mientras que el método de puntos interiores inicia desde un punto interior cualquiera. 10. CONCLUSIONES. Se presentó una nueva metodología para resolver el problema de transportes aplicando el método de puntos interiores.

Los resultados aplicando el algoritmo de transportes y el método de puntos interiores con una tolerancia de 1x10-6 son mostrados en la tabla 2 El algoritmo fue implementado usando Matlab 6.5 en un computador con procesador Pentium 4 de 2.66GHz. En la figura se muestran los tres criterios de convergencia del método de puntos interiores.

La metodología general del método fue particularizada al problema de transportes con el fin de estudiar su comportamiento. Esta metodología puede ser mejorada aplicando métodos de puntos interiores de alta orden los cuales pueden lograr mayor velocidad de convergencia.

11. BIBLIOGRAFÍA. 10

10

10

10

10

10

10

10

15

[1] BAZARAA M, SHERALI H, SHETTY C. Linear programming and network flows. John Wiley and Sons. 1983.

Primal Dual Optimalidad

10

5

[2] PRAWDA W, Juan. Métodos y modelos de investigación de operaciones. Volumen 1. Editorial Limusa. México D.F. 1986.

0

[3] RIDER, Marcos Julio. Método de puntos interiores para optimización en sistemas eléctricos. En: Seminario de optimización en sistemas de potencia. Pereira, Noviembre de 2004.

-5

-10

-15

-20

0

5

10

15

20

25

[4] R.J. Vanderbei. Linear programing: Foundations and Extensions. Kluwer Academic Publishers. Boston 2001. 2nd Edition. In English. Available in internet: http:www.princeton.edu/~rvdb/LPbook

Figura 2. Convergencia del método.

[5] S.J. Wright, Primal dual Interior point methods. SIAM Philadelphia. PA. 1997.

Se puede observar que criterio dual se cumple en la primera iteración mientras el criterio de convergencia primal se cumple en la iteración 15. El último en llegar a la convergencia es el criterio de optimalidad el cual se cumple en la iteración 20.

[6] TORO OCAMPO, Eliana; ESCOBAR, Antonio; CARREÑO, Edgar. Optimización de sistemas lineales usando métodos de punto interior. En. Scientia et Técnica. Año X Nº 24. Mayo de 2004.