LIBRO Optimizaciòn Lineal

Optim Optimiza izació ción n Linea Lineall Braulio Gutiérrez Pari Universidad Peruana Unión-Juliaca Facultad de Ingeniería y Arquitectura Ingenieria Enero del 2013 Chullunquiani-Perú

ii

Índice general Introducción

V

1. Preliminares 1.1. El Problema de Programación Lineal . . . . . . . . . . . . . . . . . . . . . . 1.2. El Método Simplex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1. Criterio de Optimalidad . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2. Criterio de la Razón . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3. Finalización . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1. El Algoritmo Simplex . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4. Teoremas de convergencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5. Dualida idad y Condiciones de Optima imalida idad . . . . . . . . . . . . . . . . . . . . 1.5.1. Dualidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5.2. Propiedades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5.3. Condiciones de Optimalidad Optimalidad de Karush-Kuhn-T Karush-Kuhn-Tuck ucker er (KK T ) . . . . 1.6. 1.6. El Métod étodoo de Newto ewton, n, Lagr Lagran ange ge y Penal enaliz izac ació iónn Inte Intern rnaa . . . . . . . . . . . . 1.6.1. Método de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6.2. Método de Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6.3. Método de Barrera Logarímica . . . . . . . . . . . . . . . . . . . . . . 2. Puntos Interior para programación lineal 2.1. .1. Pun Puntos Int Interio eriorr Basa asado en Barre rrera Loga ogarít rítmica ica 2.1.1. Construcción del método . . . . . . . . . 2.1.2. El Algoritmo . . . . . . . . . . . . . . . 2.1.3. Implementación y Experimentos . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

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

1 1 7 8 9 9 10 10 17 17 18 19 20 20 28 32

. . . .

39 43 44 51 52

3. Experimentos Computacionales

63

Conclusiones

74

Apéndice A

76

Apéndice B

79 iii

iv

Bibliografía

ÍNDICE GENERAL GENERAL

80

Introducción Los algoritmos de Puntos interiores surgen, con el trabajo de Karmarkar, como una alternativa de complejidad polinomial al bien establecido método simplex para el caso de programación lineal. En 1987 Kojima-Misuno-Yoshise presentan un algoritmo de Puntos interiores, llamado Primal-Dual, que seguido del trabajo de Mehrotra en 1992, fundamentan las bases de algunos de los algoritmos existentes más e ficientes para programación lineal. El método de Puntos interiores Primal-Dual para programación lineal (1994). Este método es el más eficiente y elegante dentro de la gran variedad de métodos que surgieron después de la publicación formal del primer método de puntos interiores que es, el método de las Elipsoides (1978). La literatura sobre los algoritmos de puntos interiores tipo Primal-Dual es extensa y variada, existen diferentes formulaciones, diversos resultados de convergencia y complejidad. Además de las buenas propiedades teóricas, los algoritmos Primal-Dual han demostrado poseer buen comportamiento práctico y se han utilizado en distintas aplicaciones de manera satisfactoria. Dentro de la optimización, es indiscutible que el Método Simplex (1947) representa un clásico en la resolución de problemas de programación lineal. Más aún, el Simplex ha servido como fuente en la construcción de métodos especializados en la resolución de otros problemas, tales como programación entera. No obstante, cuando se realiza un análisis teórico con respecto a su e ficiencia, el Método Simplex es considerado de tiempo de ejecución en el orden exponencial1 , y por tanto, teóricamente ine ficiente. Por otro lado, los resultados obtenidos en la práctica han mostrado que este método tiene un desempeño razonable y, frecuentemente, apenas requiere pocas iteraciones para resolver problemas de la vida real con un número de variables no muy grande. El motivo para el mal comportamiento del método Simplex frente a pro-blemas particulares, tales como los propuestos por V. Klee y J. Minty, es que se desplaza por los vértices del poliedro definido por la región de factibilidad. A esto se le suman algunas otras desventajas, 1 Es decir, que para problemas de programación lineal en espacios

de los casos, alrededor de 2 iteraciones. n

v

n

dimensionales, puede realizar en el peor

vi

INTRODUCCIÓN

tales como la implementación computacional so fisticada, peligro de ciclaje ante degeneración e ineficiencia ante problemas de gran tamaño. En los últimos treinta años surgieron otros tipos de métodos con una filosofía diferente, desplazarse por el interior de la región factible. Los resultados teóricos fueron sorprendentes, pues mostraron que el problema de programación lineal podía ser resuelto en un tiempo de ejecución polinomial. Dentro de los más representativos están el Método de las Elipsoides de Khachian y el Método de Karmarkar. Por los años noventa surgieron métodos con la misma estrategia, des-plazarse por el interior del poliedro de factibilidad, pero con la diferencia que éstos usaban tanto el primal como el dual para resolver el problema de programación lineal. Estos métodos estaban fundamentados además en el Método de Newton, el Método de Barrera y el Método de Lagrange, estas herramientas matemáticas fusionadas produjeron lo que hoy en día se denomina Método de Puntos Interiores Primal-Dual, y actualmente, es considerada la manera más e ficiente que se conoce para resolver el problema de programación lineal. Existen varias versiones de este método, una de ellas constituirá justamente el motivo de esta exposición.

Capítulo 1 Preliminares 1.1. El Problema de Programación Lineal El problema de programación lineal en la forma estándar corresponde al siguiente modelo matemático. minimizar c t x (1.1) Ax = b x≥0

donde A ∈ Rm n es de rango m, b ∈ Rm y c ∈ Rn . La expresión c t x se llama función objetivo, Ax es el producto de la matriz A con el vector x, Ax = b y x ≥ 0 si y sólo si xi ≥ 0, i = 1, ...n, se denominan las restricciones y las condiciones de no negatividad respectivamente. Los problemas de programación lineal se estudian normalmente en esta forma. Típicamente, n > m. En resumen un problema de programación lineal se dice que está en forma estándar si y sólo si. ×

Es de minimización. Sólo incluye restricciones de igualdad. El vector b es no negativo. Las variables x son no negativas. Existen otras formas de escribir el problema de programación lineal, por ejemplo, la forma canónica.

Minimizar c t x Ax ≥ b x≥0

1

(1.2)

2

CAPÍTULO 1. PRELIMINARES

Es posible transformar el problema (3.1) en uno de la forma (1.1) mediante la introducción de un vector cuyas componentes son denominados variables auxiliares, del siguiente modo: Minimizar c t x (1.3)

Ax − s = b x, s ≥ 0

donde s ∈ Rm , tal que s = Ax − b, estas variables auxiliares se denominan variables de exceso. Debemos resaltar aquí que con la introducción de variables auxiliares, el problema original dado en (3.1) es transformado en otro problema de mayor dimensión, pues fueron añadidas m nuevas variables, que corresponden a las componentes del vector s. En estas condiciones, resolver (1.3) no es lo mismo que resolver (3.1), pero si conseguimos resolver y encontrar una solución [x s]t para el problema en (1.3), entonces x será una solución para (3.1). En los párrafos anteriores asumimos que fueron introducidas exactamente m variables auxiliares de exceso al problema (3.1), pero en realidad no siempre es necesario tal número de variables, pudiéndose añadir variables auxi-liares sólo donde se requiera, es decir, donde las restricciones sean de desigualdad. Por otro lado, para transformar un problema del tipo Minimizar c t x (1.4)

Ax ≤ b x≥0

a la forma estándar, podemos introducir un vector auxiliar de variables de holgura , del siguiente modo: Minimizar c t x (1.5) Ax + h = b x, h ≥ 0

Finalmente, debemos aclarar que la función objetivo en el problema (1.5) es diferente a la función objetivo del problema original, debido a que ahora es de la forma ctx + 0h. Análogamente para el problema dado en (1.3). Cuando formulamos un problema de programación lineal, algunas veces las variables pueden no satisfacer las condiciones de no negatividad, entonces podemos reemplazar por otras variables no negativas mediante las siguientes sustituciones: Si xi es una variable de tipo irrestricta ( xi xi = xi − xi , donde, x i , xi ≥ 0. 0

Si x j

00

0

∈ R),

entonces podemos cambiarla por

00

˜ j = x j − l j , claramente ahora x˜ j ≥ 0. ≥ l j , entonces hacemos x

3

1.1. EL PROBLEMA DE PROGRAMACIÓN LINEAL

Si x j

ˆ j = u j − x j , ≤ u j , entonces hacemos x

observe que xˆ j

≥ 0.

En lo sucesivo y por defecto, asumiremos también que siempre trabajamos con un problema lineal en la forma estándar, salvo se aclare lo contrario.

Conjunto convexo: Un conjunto K es convexo, si la combinación convexa de dos elementos cualesquiera del conjunto está en K . Es decir, si dados x 1 , x2 ∈ K , para todo λ ∈ [0, 1] , se verifica que [x1 , x2 ] = λx1 + (1 − λ)x2 ⊂ K

Figura 1.1: Conjunto convexo y no convexo

Región Factible: Para el problema de programación lineal en la forma estándar, la región factible es el conjunto P = {x ∈ Rn : Ax = b, x ≥ 0} . La definición es análoga para problemas en otros formatos. En programación lineal, la región factible es un poliedro convexo Solución Factible: Es un punto situado en la región factible. Es decir, x es solución factible si x ∈ P. Poliedro: Región definida por la intersección de un conjunto finito de semiespacios. Un poliedro es un conjunto convexo. Politopo: Poliedro no vacío y acotado

Figura 1.2: Polítopo

4


Hiperplano: Un hiperplano en Rn es un conjunto de la forma

©

h = x ∈ Rn : pt x = k

ª

donde p ∈ Rn , p 6 = ¯0, es un vector columna y k es una constante escalar. Un hiperplano es la extensión de lo que significa una recta en R2 o un plano en R3 . Debemos notar que el vector p es normal al hiperplano h

Restricción Activa: Dada la desigualdad pt x ≥ k, donde p ∈ Rn y k es un escalar. Si x ∈ Rn es tal que satisface pt x = k , entonces se dice que x hace activa la desigualdad pt x ≥ k. En programación lineal, estas desigualdades constituyen las restricciones, en esas condiciones, se dice entonces que x hace activa dicha restricción. Para el caso de una restricción de igualdad pt x = k, es claro que cualquier x¯ que satisfaga dicha restricción hará activa ésta. 0

0

0

0

Semiespacio: Un hiperplano divide el espacio Rn en dos subregiones, llamadas semiespacios. Así, un semiespacio es el conjunto

©

ª

©

ª

h1 = x ∈ Rn: pt x ≥ k

donde p es un vector columna diferente de cero y k es un escalar. El otro semiespacio es el conjunto de puntos h2 = x ∈ Rn: pt x ≤ k donde h 1 ∪ h2 =

n

R

Punto extremo: Consideremos una región factible en programación lineal. Un punto extremo xˆ ∈ Rn es la intersección de n hiperplanos linealmente independientes. Es decir, xˆ es un punto extremo si hace activa n restricciones linealmente independientes. Si xˆ fuera también factible, entonces se dice que es un punto extremo factible. Punto Extremo No Degenerado: Un punto extremo no degenerado x˙ se dice no degenerado, si exactamente hace activa n restricciones. Cuando más de n restricciones son activas con respecto a x, ˙ entonces el punto se denomina extremo degenerado. Solución Básica: Consideremos el problema de programación lineal en la forma estándar definida en (1.1) donde existe n variables y m restricciones, con m < n, además asumiremos inicialmente que el rango de A es m. Minimizar z = ct x Ax = b

(1.6) (1.7)

5

1.1. EL PROBLEMA DE PROGRAMACIÓN LINEAL

x≥0

(1.8)

Una solución básica n−dimensional se obtiene eligiendo desde A una submatriz cuadrada B de rango m, a esta submatriz cuadrada e inversible se le llama matriz básica. Reordenando las columnas de A, si fuera necesario, y también las componentes de x, particionamos la matriz A y el vector de x de modo que se mantenga la compatibilidad: A = [B x =

N ]

" # xB xN

(1.9)

(1.10)

Al vector xN se le denomina vector no básico y a sus n − m componentes se les llama variables no básicas. y a sus m componentes se les llama variables básicas. Reemplazando (1.9) y (1.10) en (1.7) y despejando xB en función de x N , desde [B

N ]

" # xB xN

= b

obtenemos BxB B −1 BxB xB

+ NxN = b + B −1NxN = B −1 b + B −1NxN = B −1 b

xB = B −1 b − B −1 NxN

(1.11)

Obsérvese que si damos valores arbitrarios a las variables del vector xN , y evaluamos xB según (1.11), el vector x =

" # xB xN

satisface automáticamente (1.7), pero aún no

satisface (1.8). Si hacemos xN = 0 y calculamos nuevamente xB según (1.11), el vector x =

" #" #" # xB xN

=

xB 0

=

B −1b 0

además será un punto extremo, pues n restricciones linealmente independientes son activas en x. Esto es, n − m restricciones activas a causa de xN = 0 y m restricciones activas a causa de Ax = b. Este último vector x , calculando de ese modo, se denomina solución básica. Es decir, en programación lineal, una solución básica es un punto extremo.

Solución Básica Factible (SBF ): Haciendo x N = 0 y calculando x B según (1.11), si x B = B −1 b ≥ 0, entonces la solución básica x =

" #" # xB xN

=

B −1 b 0

6


será además factible, pues satisface (1.7) pero ahora también (1.8). En este caso, la solución básica se denomina solución básica factible. En programación lineal, una solución básica factible es un extremo factible.

Solución Básica Factible No degenerada: Es una solución básica factible donde x B > 0. Observe que una solución básica factible no degenerada hace exactamente n restricciones activas, por lo que es un punto extremo no degenerado. Cuando una o más variables básicas toman el valor de cero, la solución básica factible se denomina degenerada. Note que una solución básica factible degenerada hace más de n restricciones activas, lo que significa que es un punto extremo degenerado. Dirección Extrema: Una dirección extrema de un conjunto convexo es una dirección factible, la cual no puede ser representada como una combinación lineal positiva de dos direcciones factibles distintas (linealmente independientes). Por lo tanto, en programación lineal podemos decir que el conjunto de direcciones extremas generan todas las direcciones factibles. 6 0, d es una dirección de Dirección de Decrecimiento: Dado x factible y d ∈ Rn , d = decrecimiento partiendo de x, si existe  > 0 tal que, para todo α ∈ h0, ] ct (x + αd) < ct x

Para el caso de minimización, una dirección es de decrecimiento si, y sólo si, ctd < 0. Obviamente, si para un x factible no existe direcciones de decrecimiento d, tal que el punto factible, entonces una solución óptima. x + αd sea x es 6 0, obtenido de la Claramente, podemos veri ficar que el vector d = −∇(ct x) = −c = función objetivo, es una dirección de decrecimiento. Esto se debe a que. ct d = ct (−c) = − kck2 < 0

El vector −c no sólo es un vector de decrecimiento, sino que es el vector que proporciona ¯ el máximo decrecimiento. Es decir, si consideramos otra dirección de decrecimiento d, tal que k−ck = d¯ , entonces partiendo de un punto factible x y desplazándonos en las ¯ en otras palabras direcciones −c y d¯, el punto x − αc es mejor que el punto x + αd,

°°

¯ ct (x − αd) ≤ ct (x + αd)

Función Convexa: Si K conjunto convexo no vacío y f : K ⊂ Rn dice que f es convexa si:

→ R es

f (λx1 + (1 − λ)x2) ≤ λf (x1 ) + (1 − λ)f (x2)

una función, se

7

1.2. EL MÉTODO SIMPLEX f

λ f(x1)+(1- λ )f(x2)

f(x1 )

f(x2 )

x1

x = λ x1+(1- λ )x2 λ

Figura 1.3: función convexa

1.2. El Método Simplex Existen muchas variantes del método Simplex, aunque todas ellas se basan en la misma idea central. En esta sección se describe una de tales versiones. El método Simplex inicia con una solución básica factible, y a cada paso, se desplaza hacia una solución básica factible que es un mejor candidato que la anterior, al menos no peor, esto se repite hasta una solución básica factible óptima. Dada un problema de Programación lineal en su forma estándar Minimizar z = ct x

Ax = b

(1.12)

x ≥ 0

de rango m , (m < n), b ∈ Rm y c, x ∈ Rn Supongamos que tenemos una base B y una matriz N, ambas asociadas a la matriz A. Además las variables básicas y no básicas x B y x N , respectivamente. A ∈ Rm

n

×

x =

" # xB xN

, A = [B

N ] y c =

" # cB cN

, esto en (1.12) obtenemos

xB = B −1 b − B −1 NxN

= B −1b −

X

B−1 a j x j

(1.13)

j ∈R

= ¯b −

X

y j x j

j ∈R

donde R es el conjunto de los índices de las variables no básicas, donde a j representa las

8


j −columna de A. Además y j = B −1 a j y ¯b = B −1 b, si hacemos xN = 0, el punto x se convierte

en la solución básica con valor objetivo t B

t N

z0 = [c

c ]

" # B −1 b 0

= c tB B −1 b

Usando (1.13), la función objetivo en (1.12) puede ser expresado por z = ct x = ctB xB + ctN xN = ctB

Ã

B −1 b −

X

B −1 a j x j

j ∈R

= c tB B −1 b −

X ÃX X¡

c j x j

ctB B −1 a j x j +

= ctB B −1 b −

ctB B −1 a j −

= ctB B −1 b −

c j x j

(1.14)

j ∈R

j ∈R

c j

x j

j ∈R

ctB B −1 a j − c j x j

j ∈R

= z0 −

+

j ∈R

j ∈R

X

!X X X! ¢

(z j − c j ) x j

j ∈R

donde z j = ctB B −1 a j , j ∈ R, se observa que z = z0 cuando x es una solución básica ( xN = 0). Usando (1.13) y (1.14), el problema (1.12) puede ahora ser visto como Minimizar z = z0 −

X

(z j − c j ) x j

j ∈R

xB +

X

y j x j = ¯b

(1.15)

j ∈R

x j ≥ 0, j ∈ R, xB ≥ 0

de (1.15). los coeficientes z j − c j de las variables no básicas x j , j reducidos.

∈ R. se

denominan costos

1.2.1. Criterio de Optimalidad El método simplex reconoce que la actual solución básica factible es óptima, cuando todos los costos reducidos a las variables no básicos son menores o iguales a cero, es decir. z j − c j ≤ 0,

9

1.3. FINALIZACIÓN ∀ j ∈ R

1.2.2. Criterio de la Razón Si existe k ∈ R tal que zk − ck > 0, incrementando el valor de la variable no básica xk obtendremos una reducción en la función objetivo, que es justamente lo que deseamos en un problema de minimización. Pero, esto puede ocasionar que alguna de las variables básicas se torne negativa (infactible). Cuando la actual solución básica factible no es óptima, el método simplex sólo modifica una variable no básica en cada paso. Para eso, calcula el costo reducido zk − ck = m´ ax {z j − c : j ∈ R} > 0, elige la variable xk para incrementarla desde su nivel cero y hace x j = 0, ∀ j ∈ R − {k} . Ahora, seleccionada x k y haciendo x j = 0, la función objetivo y las restricciones de (1.15) pueden ser vistas del siguiente modo:

z = z0 − (zk − ck )xk

⎡ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣

xB1 xB2

.. .

xBr

.. .

xBm

⎤ ⎡ ⎥⎥ ⎢⎢ ⎥⎥ ⎢⎢ ⎥⎥ = ⎢⎢ ⎥⎥ ⎢⎢ ⎥⎥ ⎢⎢ ⎦ ⎣

¯b1 ¯b2

.. .

¯br

.. .

¯bm

⎤ ⎡ ⎥⎥ ⎢⎢ ⎥⎥ ⎢⎢ ⎥⎥ − ⎢⎢ ⎥⎥ ⎢⎢ ⎥⎥ ⎢⎢ ⎦ ⎣

y1,k y2,k

.. .

⎤ ⎥⎥ ⎥⎥ ⎥⎥ x ⎥⎥ ⎥⎥ ⎦

k

yr,k

.. .

ym,k

(1.16)

(1.17)

para evitar que alguna variable básica se torne negativa, el valor de xk debe ser elegido mediante el criterio de la razón ¯br xk = = m´ın yr,k

½

¯bi : y i,k > 0, yi,k

¾

i = 1,...,m

1.3. Finalización El simplex finaliza reconociendo que: Una única solución óptima: si z j − c j < 0 Existen soluciones óptimas alternativas: Si z j − c j ≤ 0, ∀ j ∈ R, pero alguna variable no básica, digamos xq , q ∈ R, tiene su correspondiente costo reducido z q − cq = 0, El problema tiene solución óptima ilimitada: Dada una solución básica factible, si zk − ck > 0 y yk ≤ 0, para alguna variable no básica xk , entonces la solución óptima, y por

consecuencia, el valor objetivo óptimo, serán ilimitados.

10


1.3.1. El Algoritmo Simplex Elegir una matriz básica B tal que x =

" # B −1 b 0

sea una solución básica factible.

Paso 1: Calcular x B = B −1 b = ¯b. hacer x N = 0 y calcular z = ctB xB . Paso 2: Calcular w = ctB B −1 y hallar zk − ck = max {zi − c j } , donde z j − c j = wa j − c j j ∈R

1. Si zk − ck ≤ 0, detenerse, la SBF óptima es x =

" # B −1 b 0

2. Caso contrario, ir al paso 3.

Paso 3: Calcular y k = B −1 ak 1. Si yk ≤ 0, detenerse. Existe soluciones óptimas ilimitadas 2. Caso contrario, ir al paso 4.

Paso 4: Calcular

¯br xk = = m´ın yr,k

½

¯bi : y i,k > 0, yi,k

¾

i = 1,...,m

1. Actualizar la base B. cambiando aBr por a k . 2. Actualizar R. cambiando k por el índice B r . 3. Volver al paso 1.

1.4. Teoremas de convergencia Teorema 1.1 La región factible P = {x ∈ Rn : Ax = b, x ≥ 0} , de soluciones de un problema de programación lineal, es un conjunto convexo

Prueba. Sea P la región factible del problema (1.1), sea x1, x2 ∈ P arbitrarios, definamos x = λx1 + (1 − λ)x2 con λ ∈ [0, 1] , queremos probar que x ∈ P y x ≥ 0. i) Como Ax1 = b y x 1 ≥ 0, Ax2 = b y x 2 ≥ 0, entonces. Ax = A(λx1 + (1 − λ)x2 ) = λAx1 + (1 − λ)Ax2 = λb + (1 − λ)b = b

ii) Ahora, x 1 ≥ 0, x2 ≥ 0. El hecho de que λ ∈ [0, 1] implica que λ ≥ 0 y que (1 − λ) ≥ 0, de aquí se tiene que λx1 ≥ 0, (1 − λ)x2 ≥ 0, por lo que x = λx1 + (1 − λ)x2 ≥ 0. Por lo tanto P es un conjunto convexo.

11

1.4. TEOREMAS DE CONVERGENCIA

Definición 1.1 Sea P ⊂

n

un conjunto convexo y f : siguiente formulación matemática dada por. R

n

R

→ R una función convexa. La

Minimizar f (x) x ∈ P

(1.18)

se denomina problema de programación convexa.

El problema de programación convexa tiene importantes propiedades que lo distingue de problemas más generales, uno de los resultados más interesantes está dado en el siguiente teorema.

Teorema 1.2 En un problema de programación convexa. a) Todo minimizador local es minimizador global. b) El conjunto de minimizadores es convexo. c) Si f es estrictamente convexa, no puede haber más de un mininizador.

Prueba. Para mostrar (a) supongamos que x∗ sea un minimizador local pero no global del = x∗ , tal que problema de programación convexa dada en (1.18), entonces existe un x ∈ P , x 6 f (x) < f (x∗ ). Ahora, para λ ∈ [0, 1] , consideremos la combinación convexa xλ = (1−λ)x∗ +λx, por la convexidad de P , tenemos que xλ ∈ P , y por la convexidad de f f (xλ )

= f ((1 − λ)x∗ + λx) ≤ (1 − λ)f (x∗) + λf (x) = f (x∗ ) + λ(f (x) − f (x∗ )) < f (x∗ )

Ahora, para λ > 0 suficientemente próximo de 0, xλ se torna arbitrariamente próximo pero diferente de x∗ , pero además f (xλ) < f (x∗ ). Esto significa que para cualquier vecindad V (x∗ , ε) es posible encontrar un xλ ∈ V (x∗ , ε) tal que f (xλ ) < f (x∗ ), lo cual contradice el hecho que x ∗ es un minimizador local. Para mostrar (b) nosotros llamamos de S al conjunto de todos los minimizadores del problema de programación convexa. sea x, y ∈ S, entonces f (x) = f (y) ≤ f (λx + (1 − λ)y), ∀ λ ∈ [0, 1]

12


Por la convexidad de f tenemos f (λx + (1 − λ)y) ≤ λf (x) + (1 − λ)f (y) = f (y) + λ(f (x) − f (y))

(f (x) = f (y))

= f (y) + λ(0) = f (y), ∀ λ ∈ [0, 1]

Entonces f (x) = f (y) ≤ f (λx + (1 − λ)y) ≤ f (x) = f (y), ∀ λ ∈ [0, 1] , este último significa que λx + (1 − λ)y también es un minimizador y por tanto λx + (1 − λ)y ∈ S, de donde S es convexo. Para la parte (c), supongamos que existe x, y ∈ S con x = 6 y, entonces para λ ∈ [0, 1] se tiene que f (x) = f (y) ≤ f (λx + (1 − λ)y), debido a que x, y son minimizadores globales. Ahora, como f es estrictamente convexa se verifica f (λx + (1 − λ)y) < f (x) = f (y)

con lo que tenemos la contradicción deseada y la prueba está concluida.

Teorema 1.3 El problema de programación lineal en la forma estándar es un problema de programación convexa.

Prueba. Dado el siguiente problema de programación lineal Minimizar f (x) = ct x sujeto a

Ax = b x ≥ 0

debido a la linealidad de f , satisface la siguiente condición ct (λx + (1 − λ)y) = λct x + (1 − λ)ct y,

n ∀x, y ∈ R ,

λ ∈ [0, 1]

Luego, f es una función convexa. Falta sólo veri ficar que la región de factibilidad es un conjunto convexo. Para esto, de finamos Ω = {x ∈ Rn : Ax = b, x ≥ 0} . Dado x1 , x2 ∈ Ω, luego Ax1 = b, x1 ≥ 0, Ax2 = b, x2 ≥ 0. Sea x = λx1 + (1 − λ)x2 y λ ∈ [0, 1] Ax = A(λx1 + (1 − λ)x2 ) = λAx1 + Ax2 − λAx2 = b

Además λx1 + (1 − λ)x2 ≥ 0

Por lo tanto x ∈ Ω, y Ω es convexo.

13


Consideremos ahora el problema de programación lineal en la forma estándar Minimizar c t x

x ∈ P

donde P = {x ∈ Rn : Ax = b, x ≥ 0} , b ∈ Rm , A ∈ Rm

n

×

(1.19)

es de rango m, con m < n.

Teorema 1.4 Un extremo factible x de P es una solución básica factible para (1.19). Inversamente, si x es una solución básica factible de (1.19), entonces x es un punto extremo factible de P.

Prueba. Como x es un punto extremo de P , existe n hiperplanos linealmente independientes pasando por x. Ahora, como x es factible, satisface Ax = b, de donde ya existen m restricciones linealmente independientes y activas en x, pues hemos asumido que el rango de la matriz A es m. Naturalmente, los restantes n − m hiperplanos linealmente independientes pasando por x deben ser del tipo xi = 0, donde xi es una componente del vector x. Denotemos tales componentes xi como parte del vector xN , y las columnas i de A asociadas a cada xi por la matriz N . Mientras que las demás componentes x j de x la hacemos parte del vector xB , y las respectivas columnas j de A por la matriz B . Reordenando, si fuera necesario, de modo que no se altere el problema (1.19), particionamos x de A de modo tal que: x =

" # xB xN

A = [B N ]

Por hipótesis, el vector x pertenece a P y satisface Ax = b y x ≥ 0. Es decir, satisface el sistema de n × n : BxB + NxN = b (1.20) IxN = 0

donde I ∈ R(n−m) (n−m) es la matriz identidad. El sistema (1.20) implica que xN = 0 y Bm m inversible, debido a la independencia lineal de las ecuaciones. Luego, x B = B −1b. como x es factible, se tiene que B −1b > 0, esto significa que ×

×

x =

es una solución básica factible.

" # B −1 b 0

14


Inversamente, si x es una solución básica factible, considerando (1.20), x satisface Ax = b, es decir, hace m restricciones linealmente independientes activas, pues B −1 existe. Además, n − m restricciones linealmente independientes son activas debido a xN = 0, lo cual nos dice que x, definido de ese modo, hace n restricciones linealmente independientes activas. Podemos concluir con respecto al problema dado en (1.19) que: La colección de puntos extremos factibles corresponde a la colección de soluciones básicas factibles. Ambos son no vacíos desde que la región factible sea no vacía. Para resolver el problema de programación lineal en la forma estándar, basta analizar los puntos estremos del poliedro de factibilidad, pues si una solución óptima existe, entonces existirá un extremo óptimo. Además, sólo tenemos que veri ficar que tal punto extremo es un minimizador local.

Ejemplo 1.1 Considere el siguiente problema de programación lineal Minimizar x1 + x2 4x1 − x2 ≥ 2 2x1 + 3x2 ≥ 8 −3x1 + 5x2 ≤ 15 x1 − x2 ≤ 3 10x1 + 7x2 ≤ 70 x1 , x2 ≥ 0

(1.21)

El grá fi co de la región de factibilidad y una curva de nivel de la función objetivo están representadas en la fi gura adjunta

5

x2

4

3

2

1

x1 0

1

2

3

4

5

Figura 1.4: Región factible y una curva de nivel


15

Vemos que la curva de nivel asociada a ctx = 0 es la recta que pasa por el origen, note 1 que ∇f (x1 , x2 ) = c = . Además, se sabe que ∇f (x1 , x2 ) es un vector que apunta en 1 la dirección donde la función objetivo más crece. Por lo tanto, −c estará apuntando hacia la dirección en la cual la función objetivo alcanza su mejor decrecimiento.

" #

Para minimizar el problema grá ficamente, tenemos que partir de un punto factible y desplazarnos en la dirección −c sobre R, pues en esa dirección la función objetivo decrecerá, pero además el punto al cual nos desplacemos debe ser también factible. También podemos observar del gráfico que la curva de nivel asociado a ct x = 3, es la recta que pasa por el punto [1 2]t . El único modo de mantener la factibilidad y al mismo tiempo obtener el menor valor de la función objetivo c t x, está relacionado en este punto [1 2]t , pues en ese punto la función objetivo vale 3, si continuamos desplazándonos en la dirección −c caeremos en infactibilidad.

Figura 1.5: Curvas de nivel de f sobre R Podemos concluir entonces que la solución óptima es x∗ = [1 2]t y el valor objetivo óptimo es ct x∗ = x∗1 + x∗2 = 1 + 2 = 3

16


del gráfico se puede observar que la solución óptima es un extremo, un vértice del poliedro factible. Llevando a su forma estándar, para ello introduciendo las variables auxiliares, tenemos Minimizar x 1 + x2 + 0x3 + 0x4 + 0x5 + 0x6 + 0x7 4x1 − x2 − x3 = 2 2x1 + 3x2 − x4 = 8 −3x1 + 5x2 + x5 = 15 x1 − x2 + x6 = 3 10x1 + 7x2 + x7 = 70 x1, x2 , x3, x4 , x5 , x6 , x7 ≥ 0

Los datos a introducir en el computador son

⎡ 4 ⎢⎢ 2 ⎢ −3 A = ⎢ ⎢⎢ ⎣ 1

−1 −1

3 5 −1 10 7

0 0 0 0

0 −1 0 0 0

0 0 1 0 0

0 0 0 1 0

0 0 0 0 1

⎤ ⎥⎥ ⎥⎥ , ⎥⎥ ⎦

⎡2⎤ ⎢⎢ 8 ⎥⎥ ⎢ 15 ⎥⎥ , b = ⎢ ⎢⎢ ⎥⎥ ⎣3⎦ 70

⎡1⎤ ⎢⎢ 1 ⎥⎥ ⎢⎢ 0 ⎥⎥ ⎢⎢ ⎥⎥ c = ⎢ 0 ⎥ ⎢⎢ 0 ⎥⎥ ⎢⎢ ⎥⎥ ⎣0⎦ 0

colB = [1 2 3 4 5];

Donde colB : vector de índices básicos iniciales. (SIMPLEX (A,b,c,colB) se encuentra en el apéndice), llamando SIMP LEX (A,b,c,colB) nos otorga:

x= 1.00 2.00 0 0 8.00 4.00 46.00 iter= 3.00 z= 3.00

1.5. DUALIDAD Y CONDICIONES DE OPTIMALIDAD

17

después de 3 iteraciones principales ejecutadas por el programa, [z,x,iter] = SIMP LEX (A,b,c,co Ahora proyectando x 3 sobre el x1 x2 -espacio, podemos ver lo que sucedió en el problema original. Así, la solución aproximada de (1.21) es

" #" # x1 x2

=

1.00 2.00

el valor óptimo es: z =3.00, y la solucion óptima es x = [1 2]t

1.5. Dualidad y Condiciones de Optimalidad Por lo general, la dualidad y las condiciones de optimalidad son usadas para la construcción de métodos que son variaciones del Simplex, tal como veremos en el capítulo 4, el método Primal-Dual de Puntos Interiores para programación lineal. Uno de los aspectos teóricos más importantes asociado a la programación lineal es la dualidad . Mediante esto, es posible caracterizar el problema de programación lineal desde otro punto de vista.

1.5.1. Dualidad Dado el problema de programación lineal en la forma estándar Minimizar c t x Ax = b x≥0

(1.22)

Definiendo el dual de un problema de programación lineal en la forma estándar Maximizar bt y At y + z = c z≥0

(1.23)

donde z ∈ Rn . Dado un problema de programación lineal, denominado problema primal , existe otro problema de programación lineal, denominado problema dual , intimamente relacionado con él. Se dice que ambos problemas son mutuamente duales.

18


Problema de Minimización

Problema de Maximización [coe fi cientes ]t número de restricciones número de variables

⇔

[coe fi cientes ] número de variables número de restricciones

⇔ ⇔

Variable

Restricción

≥0

⇔

≤

≤0

⇔

≥

Irrestricta Restricción

⇔

=

≥

⇔

≥0

≤

⇔

≤0

=

⇔

Irrestricta

Variable

Relación variable - restricción entre el primal y su dual

Ejemplo 1.2 Primal

Dual

Maximizar 3x1 + 5x2 x1 ≤ 4 2x2 ≤ 12 3x1 + 2x2 ≤ 18 x1 , x2 ≥ 0

Minimizar 4y1 + 12y2 + 18y3 y1 + 2y3 ≥ 3 2y2 + 2y3 ≥ 5 y1 , y2, y3 ≥ 0 .

⇔

Primal

Dual

Minimizar 3x1 + 5x2 − 4x3 x1 + 2x2 − 4x3 ≥ 6 3x1 − 2x2 + x3 = −10 x1 ≤ 0, x2 ≥ 0, x3 irrestricta

Maximizar 6y1 − 10y y1 + 3y2 ≥ 3 2y1 − 2y2 ≤ 5 −4y1 + y2 = −4 y1 ≥ 0, y2 irrestricto

.

⇔

1.5.2. Propiedades Estas propiedades en [YZH] página 34, se describen como proposiciones y en [LDG] como Lema y corolario

Lema 1.1 Si x, y son soluciones factibles del problema de minimización y maximización, respectivamente, entonces bt y ≤ ct x

19

1.5. DUALIDAD Y CONDICIONES DE OPTIMALIDAD

Prueba. Consideremos la forma canónica de dualidad, sean x¯ e y¯ soluciones factibles del problema de minimización y maximización, respectivamente. Veamos que A¯ x ≥ b,

x¯ ≥ 0

(1.24)

At y¯ ≤ c,

y¯ ≥ 0

(1.25)

Además Multiplicando por y¯ ≥ 0 en (1.24) y por x¯ ≥ 0 en (1.25), obtenemos y¯t A¯ x ≥ y¯t b = b ty¯

y x ¯t At y¯ ≤ x¯t c = c tx¯

de donde b t y¯ ≤ y¯t A¯x ≤ x¯t At y¯ ≤ x¯tc ≤ ct x¯

Corolario 1.1 Si x∗ e y∗ son soluciones factibles del problema de minimización y maximización, respectivamente, tal que ct x∗ = bt y∗ , entonces x∗ e y∗ son soluciones óptimas de sus respectivos problemas asociados.

Prueba. Supongamos que x∗ no sea una solución óptima del problema de minimización, entonces ct xˆ < ct x∗ para algún ˆx factible de este problema. Usando el lema anterior, tendríamos que b t y∗ ≤ ct xˆ < ct x∗ , y en consecuencia b t y∗ 6 = c t x∗

1.5.3. Condiciones de Optimalidad de Karush-Kuhn-Tucker (KKT ) Ya que en este trabajo hemos dado más importancia al problema de programación lineal en su forma estándar. Consideremos el problema de programación lineal en esta forma, al cual denominaremos primal: Minimizar c t x (1.26) Ax = b x≥0

Por lo que estas condiciones KKT para el problema (1.26) están de finidas por. Si x ∈ Rn es una solución óptima del problema de programación lineal estándar si, y sólo si, existen vectores w ∈ Rm y v ∈ Rn tales que: Ax = b, At w + v = c,

(Factibilidad Primal)

x≥0

w irrestricto, vt x = 0

v ≥0

(Factibilidad Dual) (Complementariedad)

20


Estas condiciones de Karush Kunh Tucker representan condiciones necesarias y su ficientes para un minimizador global del problema (1.26). Estas tres condiciones de óptimalidad tienen suma importancia, tanto teórico como en la práctica, pues nos permiten enfocar el problema de programación lineal desde otro punto de vista. Más aún, juegan un rol importante en la construcción del método de puntos interiores primal-dual. Dada una solución óptima del problema (1.26) x∗ ∈ Rn , entonces, existen vectores v∗ ∈ Rn y w∗ ∈ Rm tales que (x∗ , v∗ , w∗ ) satisface estas condiciones. Inversamente, si existen x ¯, ¯v y w, ¯ tal que (¯ x, v¯, w) ¯ satisfacen también las condiciones, entonces ¯ x es una solución óptima del problema (1.26).

1.6. El Método de Newton, Lagrange y Penalización Interna En esta sección describiremos el método de Newton , en el estudio de sistemas de ecuaciones no lineales. Este método tiene un lugar importante dentro de la optimización, los cuales serán empleados en la construcción del método Primal-Dual para programación lineal. El método de Newton. Este método numérico es aplicado para resolver sistemas de ecua-

ciones no lineales y es aplicado en la resolución de problemas de optimización sin la presencia de restricciones. El método de Lagrange . Este método numérico es aplicado frecuentemente para resolver

problemas de optimización con restricciones de igualdad. El método de Penalización interna . Este método numérico es aplicado frecuentemente para

resolver problemas de optimización con restricciones de desigualdad. los cuales nos permitirán construir el método de Puntos interiores Primal dual para programación lineal.

1.6.1. Método de Newton Este método numérico es aplicado con mayor frecuencia en la resolución de problemas de optimización sin restricciones. y será modi ficado antes de poderlo utilizar en puntos lejanos de la solución, añadiendo en un inicio una longitud de paso α. Recordemos inicialmente el método de Newton para hallar un cero de una función de una variable, es decir resolver f (x) = 0 con primera derivada contínua f ∈ C 1 (R). Dado un punto inicial xk , con f (xk ) 6 = 0 para todo k = 0, 1, 2, 3, ..., . 0

1.6. EL MÉTODO DE NEWTON, LAGRANGE Y PENALIZACIÓN INTERNA

21

Partiendo de un punto xk , como en la fi gura adjunta, trazamos la recta tangente a la curva y = f (x) recordemos la ecuación punto pendiente y − y0 = m(x − x0 ) en el punto inicial (xk , f (xk ).La cual es f (x) − f (xk ) = f (xk )(x − xk ) despejando f (x) tenemos 0

(Recta tangente a la curva)

f (x) = f (xk ) + f (xk )(x − xk ) 0

El próximo punto,xk+1 , está definido como la solución de f (x) = 0 0 = f (xk ) + f (xk )(xk+1 − xk ) 0

que corta al eje x 1 en el punto x k+1 = xk − f (1x ) f (xk ), se observa que debe ser f (xk ) 6 =0 0

0

k

Por lo que el método consiste en construir una sucesión de finida por k+1

x

f (xk ) = x − f (xk ) k

0

k = 0, 1, 2, 3,...

¯ ¯

cuando f (xk ) < ε, donde ε es un parámetro de precisión, entonces x k será una buena aproximación a un cero de f : I ⊂ R −→ R

Figura 1.6: Newton para hallar raices de ecuaciones f (x) = 0 Pero el método de Newton no sólo está orientado hacia ecuaciones involucrando una función real de variable real, sino que su verdadero potencial surge cuando es utilizado en la resolución de sistemas de ecuaciones no lineales, uno de los problemas de cálculo fundamental en ciencia e ingenieria.

22


Sea x ∈ Rn , F : Rn −→ Rn , es decir hay n variables y n restricciones la cuál acepta primera derivada contínua. F ∈ C 1(Rn ):

⎡x ⎢⎢ x x = ⎢ ⎢⎣ ...

1 2

xn

⎤ ⎡ f (x) ⎤ ⎡ f (x , x , ...x ) ⎤ ⎥⎥ ⎢⎢ f (x) ⎥⎥ ⎢⎢ f (x , x , ...x ) ⎥⎥ ⎥⎥ , F (x) = ⎢⎢ .. ⎥⎥ ≡ ⎢⎢ ⎥⎥ .. . ⎦ ⎣ . ⎦ ⎣ ⎦ 1

1

1

2

n

2

2

1

2

n

f n (x)

f n(x1 , x2 ,...xn )

Las funciones f i : Rn −→ R f i ∈ C 1 (Rn ), para i = 1, 2,...,n son las funciones componentes. La ecuación F (x) = 0 representa un sistema de n ecuaciones con n variables: f 1 (x) = 0 f 2 (x) = 0

.. .

f n(x) = 0

El Jacobiano de F en el punto x, denotado por F (x) o J (x) es definido como la matriz cuyas componentes (i, j) son: 0

∂f i (x), ∂x j

1 ≤ i, j ≤ n

explícitamente

F (x) = ∇F (x)t 0

⎛ ⎜⎜ = J (x) = ⎜ ⎜⎝

∂f 1 ∂x 1 ∂f 2 ∂x 1

∂f 1 ∂x 2 ∂f 2 ∂x 2

· ·· · ··

∂f 1 ∂x n ∂f 2 ∂x n

∂f n ∂x 1

∂f n ∂x 2

· ··

∂f n ∂x n

.. .

.. .

...

.. .

⎞ ⎟⎟ ⎟⎟ ∈ ⎠

n×n

R

El método de Newton es de caracter iterativo y está basado en la aproximación lineal de la función F en torno al punto actual xk F (x) ≈ Lk (x) = F (xk ) + J (xk )(x − xk )

siguiendo el principio descrito anteriormente, el siguiente punto x k+1 es una solución de Lk (x) = 0

Si J (xk ) es no singular de (1.27), entonces tiene solución única, esto es 0 = F (xk ) + J (xk )(xk+1 − xk )

(1.27)


23

despejando el punto actual x k+1 xk+1 = xk −

1 F (xk ) k J (x )

xk+1 = xk − J −1(xk )F (xk )

(1.28)

Empezando de un punto inicial x0. Si x0 está suficientemente cerca a la solución x∗ satisfaciendo F (x∗ ) = 0, entonces bajo las condiciones apropiadas, podemos esperar que la sucesión converga a la solución x ∗ ; es decir l´ım xk = x∗

k→∞

Obtener xk+1 directamente de (1.28) debería ser evitado para fines computacionales, ya que calcular J −1 (xk ) es considerado costoso. Entonces una i-teración de Newton consiste en resolver un sistema lineal: (1.29) J (xk )dk = −F (xk ) (1.30)

xk+1 = xk + dk

lo que nos indica (1.30) es que una vez conocido el paso de newton dx , es posible calcular xk+1 . Aunque el método de Newton es muy atractivo en función de sus propiedades de convergencia cerca de la solución, ha de modi ficarse antes de poderlo utilizar en puntos más lejados de la solución, una primera modi ficación es añadir una longitud de paso ( α).a la dirección de newton. Vale recalcar que el método de newton elige como dirección de búsqueda dx = −J (xk )−1 F (xk )

Es importante notar que la dirección dx no se puede calcular si J (xk ) es una matriz singular. Incluso en el caso de que no lo fuese, dx no sería necesariamente una dirección de descenso cuando la matriz J (xk ) no es definida positiva. La programación convexa es un caso importante para que la dirección de newton sea de descenso. Resumiendo este método más popular para resolver sistemas de ecuaciones no lineales. Dada una función F n : Rn 7→ Rn y un sistema F n (x) = 0. El Método de Newton consiste en repetir consecutivamente, para k = 0,1,2,.., hasta que F n (xk ) < ε, donde ε > 0, la siguiente regla: k+1

x

° °

F n (xk ) = x − F n (xk ) k

0

Note que xk denota un vector de aproximación de la solución del sistema de ecuaciones F n(x) = 0, además, x 0 es el vector inicial y F n(xk ) es el Jacobiano de F n en x k . 0

24


El siguiente programa en Matlab resuelve el sistema F (x) = 0, si prec representa el parámetro de presición, J = F (x) como el Jacobiano del sistema de ecuaciones. y x0 un punto inicial. 0

function x=newton(x) iter=0; prec=0.000001; while norm(Fn(x))>prec iter=iter+1; x=x-inv(Jn(x))*Fn(x); fprintf(’iter= %2i; x = [ %7.4f %7.4f]\n’, iter,x); if iter>1000 error(’parece que newton no converge’); end end

Veamos algunos ejemplos

Ejemplo 1.3 Se quiere resolver el sistema de ecuaciones no lineales −x + 2y − 4 = 0

(x − 6)2 − y + 4 = 0 Ejecutando Newton, con los siguientes puntos iniciales de fi namos el sistema de ecuaciones y el Jacobiano

0

1 y 2

" # " #" "#

F (x) =

F n (xk ) = J (x) =

" #" #

−x1 + 2x2 − 4

(x1 − 6)2 − x2 + 2

df dx1 dg dx1

df dx2 dg dx2

=

6 tenemos, pero primero 10

;

−1

2 2(x1 − 6) −1

#

Después de ingresar la función y la matríz Jacobiana en el coputador, ejecutamos el programa 1 newton, primero con el punto inicial con un parámetro de precisión 0,000001 tenemos: 2


25

>> x = newton([1 2] ) 0

iter =

1

;

x

=

[

3.3333

3.6667]

iter =

2

;

x

=

[ 4.2667

4.1333]

iter =

3

;

x

=

[ 4.4863

4.2431]

iter =

4

;

x

=

[

4.4999

4.2500]

iter =

5

;

x

=

[

4.5000

4.2500]

x= 4.5000 4.2500

Converge a la raíz

"

4,5000 4,2500

#

proyectando en el plano x1 x2 podemos observar todo los puntos

de iteración, desde el punto inicial, hasta la raíz

Figura 1.7: Newton con punto inicial x0 = [1 2]t

26


Ahora ejecutando newton con el otro punto inicial

"# 9 3

tenemos:

>> x = newton([9 3] ) 0

iter = iter = iter = iter =

1 2 3 4

; ; ; ;

x x x x

= [8.1818 = [8.0086 = [8.0000 = [8.0000

6.0909] 6.0043] 6.0000] 6.0000]

x = 8,0000 6,0000

que también converge a otra raíz

"

8,0000 6,0000

iteración.hacia la raíz en el grá fico adjunto

#

, y también podemos observar los puntos de

Figura 1.8: Newton con punto inicial x0 = [9 3]t

Problema 1.1 Se quiere resolver el sistema de ecuaciones no lineales 5x21 + 6x1 x2 + 5x22 − 4x1 + 4x2 − 4 = 0 x21 + x22 − 1 = 0 Ejecutando el programa newton con punto inicial

" # 1 1

.


Sean

"

F (x) =

J (x) =

"

df dx1 dg dx1

5x21 + 6x1 x2 + 5x22 − 4x1 + 4x2 x21 + x22 − 1 df dx2 dg dx2

#" =

−4

27

#

10x1 + 6x2 − 4 6x1 + 10x2 + 4 2x1 2x2

#

ejecutando el programa de newton, se puede hallar una raíz en cinco iteraciones. >> x = newton([1 1] ) 0

iter= 1 ; x = [1.2500

0.2500]

iter= 2 ; x = [0.9917

0.2917]

iter= 3 ; x = [0.9575

0.2904]

iter= 4 ; x = [0.9569

0.2903]

iter= 5 ; x = [0.9569

0.2903]

x= 0.9569 0.2903

Ejemplo 1.4 Considere el sistema de ecuaciones no lineales x22 − (x1 − 1)3 = 0 ln(1 + x22 ) + e−x1 − 3x1 + 2x2 = 0 Observe que: F (x) =

"

x22

3

− (x1 − 1) = 0

ln(1 + x22 ) + e−x1 − 3x1 + 2x2 = 0

y J (x) =

"

2

−3(x1 − 1) −e−x1 − 3

2x1 2x1 +2 1+x2 2

#

#

Al ejecutar Newton con punto inicial x 0 = [1 1]t y con un parámetro de precisión ε = 0,000001 nos otorga:

x(1) = x(2) = x(3) =

0.57273607817503 0.50000000000000 0.15438003883215 -0.05711721521426 0.41714604428871 0.33025949727142

28


(4)

=

x(5) = x(6) = x(7) = x(8) =

3.40954932017067 4.48252254203493 3.30977931047506 3.60789101090839 3.41941663803305 3.75489556515174 3.41233877698485 3.74673397345497 3.41230437719520 3.74669341650677

Ejemplo 1.5 Se requiere resolver el sistema x31 + x32 − 2x23 + x24 + 3 = 0 x1 x2 − x3 x4 + 2 = 0 x1 + x22 + x3 + x4 = 0 x1 x3 + x2 x4 = 0 Al ejecutar Newton con punto inicial x0 = [1 2 3 4]t y con un parámetro de precisión ε = 0,000001 nos otorga:

-1.68820905394305

x(15) =

1.40864292542700 -0.78048646281991 0.48442062541582

Al ejecutar Newton con punto inicial x 0 = [−1 −1  = 0,000001 nos otorga:

− 1 − 1]t

y con un parámetro de precisión

0.30104099749255

x(6) =

-1.49799619587565 -1.00745590938951 -1.53757769093631

1.6.2. Método de Lagrange Este método numérico transforma un problema de optimización con restricciones de igualdad en un problema irrestricto.


29

Si intentamos resolver el problema Minimizar f (x) Sujeto a: gi (x) = 0, i = 1,...,m

He aquí tenemos un problema de optimización con restricciones de igualdad, formando la función Lagrangiana tenemos. m

L(x, y) = f (x) −

X

yi gi (x)

i=1

y entonces teniendo esta función irrestricta nos toca minimizar la función L(x, y), para ello podemos utilizar el método de Newton expuesta anteriormente. ∂L ∂f = (x) − ∂x j ∂x j

m

X

yi

i=1

∂L = −gi (x) = 0 ∂x i

∂g i (x) = 0 ∂x j

j = 1,...,n

i = 1,...,m

Este sistema de ecuaciones constituye las condiciones de optimalidad de primer orden para el problema de minimización con restricciones de igualdad ( Karush Kuhn Tucker ), las cuales son necesarias para un minimizador local. Pero cuando este problema es un problema de programación convexa, es decir, que tanto la función objetivo como el conjunto de finido por las restricciones son convexos, entonces las condiciones son además su ficientes y un minimizador local es global. O sea, resolver el último sistema nos otorga un minimizador global del problema de optimización con restricciones de igualdad.

Figura 1.9: Curvas de niveles,

∇(f ) y ∇g(x)

30


Ilustración de las condiciones de Karush-Kuhn-Tucker para el caso de una restricción de igualdad y dos variables, y también podemos observar que el mínimo se alcanza en un punto en el que el gradiente de la función objetivo y el de la restricción son linealmente dependientes. Esto es lo que representa las condiciones de optimalidad de primer orden. Recordando además que las variables de Lagrange son irrestrictas en el signo.

Ejemplo 1.6 f (x1 , x2 ) = x1 + 2x2

Minimizar s.a

g(x1 , x2) = x21 − 6x1 + x22 − 8x2 + 20 = 0

formando la función Lagrangiana

¡

L(x1 , x2 , λ) = x1 + 2x2 + λ x21 − 6x1 + x22 − 8x2 + 20

entonces el sistema a resolver es:

¢

∂L = 1 + (2x1 − 6)λ = 0 ∂x 1 ∂L = 2 + (2x2 − 8)λ = 0 ∂x 2 ∂L = x21 − 6x1 + x22 − 8x2 + 20 = 0 ∂λ

donde λ es la variable lagrangiana asociada a la restricción del problema. De finimos

⎡ 1 + (2x ⎢ F (x, λ) = ⎣ 2 + (2x

1

− 6)λ

2

− 8)λ

x21 − 6x1 + x22 − 8x2 + 20

⎤ ⎥⎦

y calculando su respectivo Jacobiano

⎡ ⎢ J (x, λ) = ⎣

⎤ −6 ⎥ −8 ⎦

2λ 0 2x1 0 2λ 2x2 2x1 − 6 2x2 − 8 0

Para resolver este problema haremos una modi ficación al programa de newton, e implementaremos en el lenguaje de programación (Matlab) a la que denominamos newton_2. Teniendo en cuenta que x es la variable, y λ que es la variable lagrangiana irrestricta en signo, que lo


31

denotaremos por y. function [x,iter]=newton_2(x,y) n=length(x);m=length(y);iter=0; while norm(Fn(x,y))>0.0001 d=-inv(Jn(x,y))*Fn(x,y); dx=d(1:n); dy=d(n+1:n+m); a=-1/min(min(dx./x),-1); x=x+0.95*a*dx; y=y+0.95*a*dy; iter=iter+1; end a es el criterio de la razón

Ejecutando el programa de newton_2, para un punto inicial de [1 2 ]t y variable lagrangiana λ = 3

>> x = [1 2 ]t , y = 3 ,

Obtenemos (1)

x

=

x(2) = x(3) = x(4) = x(5) = x(6) =

" " " " " "

1.4354 2.2771 1.6967 2.1546 1.9563 1.9981 1.9979 1.9994 1.9999 2.0000 2.0000 2.0000

# # # # # #

(1)

h h h h h h

,

y

= 1.0406

,

y(2) = 0.5294

,

y(3) = 0.4964

,

y(4) = 0.4997

,

y(5) = 0.5000

,

y(6) = 0.5000

i i i i i i

Después de 6 iteraciones principales ejecutadas por el programa, proyectando x(6) sobres x1 x2 −espacio, podemos ver lo que sucedió en el problema original. Asi como la solución es

" #" x1 x2

=

2.0000 2.0000

#

32


Figura 1.10: Lagrangiano Del gráfico podemos observar que el mínimo se alcanza en un punto en el que, gradiente de la función objetivo y el de la restricción son linealmente dependientes.

1.6.3. Método de Barrera Logarímica Este método numérico transforma un problema de optimización con restricciones de desigualdad en un problema irrestricto. Dada un problema de la forma. Minimizar f (x) Sujeto a: x j ≥ 0 añadiendo la función barrera

n

−μ

P

(1.31)

j = 1, 2, ...n

log(x j ) a la función objetivo. El problema equivalente e

j =1

irrestricto asociado a (1.31) está dado por n

X

Minimizar B(x | μ) = f (x) − μ

log(x j )

(1.32)

j =1

mientras que si μ −→ 0, entonces x(μ) −→ x∗ . Añadiendo log(x j ) a la función objetivo ocasionará que la función objetivo se incrementará in finitamente cuando x j −→ 0 (Cuando x j se aproxima a la frontera de Ω). Pero la solución óptima está justamente sobre el vértice de la región factible, por eso se hace necesario la introducción del parámetro de penalización μ ∈ R+ que equilibre la verdadera condición de la función objetiva contra el de la función barrera.


33

El conjunto de minimizadores x(μ) es llamado trayectoria central . La siguiente figura ad junta muestra el comportamiento de la función logarítmica cuando decrece el parámetro μ. Se puede observar que a medida que este parámetro tiende a cero, la función se acerca cada vez más a los ejes.

Figura 1.11: El parámetro μ 1 < μ2 La siguiente figura muestra el comportamiento de la barrera logarítmica en relación a varias restricciones cuando μ −→ 0. A medida que esto ocurre la pared que forma la barrera logarítmica se asemeja cada vez más a las paredes del polítopo, de modo que en el límite, coincidirán.

Figura 1.12: Barrera logarítmica con várias restricciones

Ejemplo 1.7 Considere el siguiente problema. Minimizar (x1 + 2)2 + (x2 − 3)2 x1, x2 ≥ 0

(1.33)

El minimizador irrestricto debería ser (−2,3), pero el minimizador para (1.33) es (x∗1 , x∗2 ) =

34


(0,3) como se ve el siguiente grá fico

Figura 1.13: Punto óptimo en [0 3]t Para resolver el problema (1.33) aplicamos barrera logarítmica para μ > 0 tenemos B(x | μ) = (x1 + 2)2 + (x2 − 3)2 − μ (log(x1 ) + log(x2 ))

Para minimizar B(x | μ), derivamos parcialmente con respecto a x 1 y x 2 ∂B μ = 2(x1 + 2) − =0 ∂x 1 x1 ∂B μ = 2(x2 − 3) − =0 ∂x 2 x2

Resolviendo el sistema obtenemos x1 (μ) = x2 (μ) =

r r

μ +1−1 2 μ 9 3 + + 2 4 2

Observamos además: l´ım x1 (μ) = 0

μ−→0

l´ım x2 (μ) = 3

μ−→0

n

X

l´ım μ

μ−→0

j =1

log(x j ) = 0

(1.34)


En (1.34), si hacemos μ = 18 345 obtenemos los puntos x1 (μ) = 94,7784 x2 (μ) = 97,2849

el cual esta muy alejado de x ∗ . y si disminuimos μ, como por ejemplo, μ = 100, tenemos x1 (μ) = 6,1414 x2 (μ) = 8,7284

el cual se acerca a x ∗ = (0,3)

500 400 300 200 100 0 -100 -200 -300 -400 4

3

2

1

0 0

0.5

1

1.5

2

2.5

3

3.5

Figura 1.14: Con μ = 100 Ahora para un μ = 0,01, entonces tendremos que los nuevos puntos x1 (μ) = 0,0025 x2 (μ) = 3,0017

4

35

36


se encuentran muy próximos a la solución óptima x ∗ .= (0,3)

40

35

30

25

20

15

10

5

0 4 3 2 1 0

0

1

0.5

1.5

2

2.5

3

3.5

4

Figura 1.15: Con μ = 000, 1

Ejemplo 1.8 dada el problema

¡ ¢

2

Minimizar x1 − 25 + 2 (x2 − 2)2 3x1 + 2x2 ≤ 6

(1.35)

podemos observar que el minimizador irrestricto debería ser ( 52 , 2) aplicando barrera logarítmica al problema (1.35) para μ > 0 se tiene

µ ¶

5 B(x | μ) = x1 − 2

2

+ 2 (x2 − 2)2 − μ log(6 − 3x1 − 2x2)

los puntos estacionarios de B(x | μ) en función de μ son: ∂B 5 3μ = 2(x1 − ) + =0 ∂x 1 2 6 − 3x1 − 2x2 ∂B 2μ = 2(x2 − 2) + =0 ∂x 2 6 − 3x1 − 2x2 y resolviendo este sistema de ecuaciones obtenemos: x1 (μ) = x2 (μ) =

77 − 77 −

p p

99(11 + 8μ) 44 11(11 + 8μ) 44

(1.36)


37

donde ( x1 (μ), x2 (μ)) se aproximan a ( 1, 32 ) a medida que μ −→ 0.

Figura 1.16: Punto Óptimo en [1 3/2]t En general, no podemos dar una solución directa para (1.33) y (1.35) como en los dos ejemplos anteriores, pero podemos usar el siguiente algoritmo: Dados x0 > 0, μ0 > 0, ε > 0 y k = 0

Paso 1. Hallar x k (μk ), el minimzador de B(x | μ) Paso 2. Parar si μk < ε. Caso contrario, elegir un μk+1, tal que 0 < μk+1 < μk

Paso 3. Hacer k = k + 1 y volver al paso 1. Entonces x k (μk ) −→ x∗ cuando μ k −→ 0 El parámetro de penalización μk debe aproximarse iterativamente hacia cero, en general consiste en hacer μ 0 = 1 y μk+1 = μ10 para k = 1, 2,... En el paso 1, para hallar x(μ) usamos el método de Newton para resolver el sistema de n ecuaciones con n variables: k

∂B ∂f (x) μ (x | μ) = = 0, − ∂x j ∂x j x j

j = 1,...,n

38


Capítulo 2 Puntos Interior para programación lineal Desde su publicación por George B Dantzig en 1947 el método simplex ha demostrado sobradamente ser altamente e ficaz para resolver todo tipo de problemas de programación lineal, el hecho de que en determinadas circuns-tancias su complejidad sea exponencial ha motivado en los últimos años un elevado número de intentos, tanto teóricos como prácticos, de obtener otros procedimientos con mejor convergencia computacional. El primero de gran transcendencia teórica es debido a L. G. Khachiyan en 1979. Este autor recopilando una serie de trabajos e ideas de los autores rusos Shor, Yudin y Nemirovskii sobre un método basado en la generación de una sucesión de elipsoides para resolver problemas de programación no lineal, los aplicó e hizo extensivo a programación lineal dando lugar al conocido como método de las elipsoides . Su principal ventaja teórica sobre el método simplex radica en que resuelve problemas de programación lineal en tiempo polinómico; es decir, posee una convergencia computacional teórica polinómica. Las implementaciones prácticas, sin embargo, han demostrado su desventaja res-pecto al método simplex ya que, por un lado, el número de iteraciones que requiere para llegar al óptimo es muy grande y, por otro, el trabajo que implica cada una de ellas es mucho más costoso que el requerido por el método simplex, pues no se pueden aplicar las técnicas de matrices dispersas tan ca-racterísticas de las implementaciones en ordenador de este último. En el año de 1 984, N. K. Karmarkar, de los laboratorios Bell de la compañía AT&T, propuso un nuevo algoritmo de complejidad polinómica para resolver pro-blemas de programación lineal. Lo contrario que el método de las elipsoides, las implementaciones prácticas a que ha dado lugar en los últimos años los presentan como un gran rival del método simplex en cualquier tipo de pro-blemas y especialmente en los de gran tamaño que surgen de aplicaciones del mundo real. El método propuesto por Karmarkar es bastante diferente del simplex. Considera el problema lineal dentro de una estructura geométrica de tipo simplicial y una estrategia de movimien39

40

CAPÍTULO 2. PUNTOS INTERIOR PARA PROGRAMACIÓN LINEAL

tos para obtener los distintos puntos del proceso iterativo en el interior del polítopo o poliedro que definen las condiciones. Este polítopo se transforma de tal manera que el punto del proceso iterativo en el que esté en cada procedimiento sea el centro del polítopo o poliedro en el espacio transformado. Este nuevo concepto o enfoque de tratar el problema es decir, llegar al óptimo a través de puntos interiores de la región factible, ha motivado numerosas nuevas ideas dentro del apartado genérico de algoritmos o métodos de puntos interiores, tanto para programación lineal como no lineal.

Método de Punto interior para Programación lineal: Se denominan métodos de punto interior precisamente porque los puntos generados por estos algoritmos se hallan en el interior de la región factible. Esta es una clara diferencia respecto del método simplex, el cuál avanza por la frontera de dicha región moviéndose de un punto extremo a otro. En la figura adjunta se muestra la trayectoria seguida para alcanzar el punto óptimo x∗ desde el punto inicial x0 .

a) El método del simplex b) El método de punto interior

Figura 2.1: Simplex y Punto interio Los procedimientos numéricos de puntos interiores para resolver problemas de programación lineal basan su estrategia en. Encontrar un punto de partida factible en el interior de la región factible del problema Definir una dirección de movimiento tal que, conservando la factibilidad del problema, moviéndose a lo largo de ella, se avance hasta un punto en el que se deduzca el valor de la función objetivo.

41 Cuándo pararse. Es decir, cuántas veces se debe realizar la operación descrita en el punto 2 y cómo identificar que se ha alcanzado el óptimo del problema. La dirección de movimiento d calculada en cada iteración del algoritmo no puede ser una dirección cualquiera. Concretamente, d debe verificar dos condiciones para ser considerada una buena dirección : Una que preserve la factibilidad del nuevo punto y la otra debe mejorar el valor de la función objetivo, la dirección de máxima pendiente: en este caso −c , por ser lineal la función objetivo.

Desventajas del Método Simplex Puede realizar un número exponencial de iteraciones, pues se desplaza por los vértices de un poliedro. Para ciertos problemas en espacios de dimensión n , puede tomar alrededor de 2 n iteraciones Implementación computacional sofisticada Peligro de ciclaje ante degeneración Ineficiencia para problemas de gran tamaño.

Ventajas y Desventajas del MPIPD Desventajas No converge a un vértice como lo hace el Simplex

Ventajas Rápido, aproximadamente nL iteraciones Fácil de implementar.

Teorema 2.1 Sea K ⊂

abierto y convexo, f : K → 7 R, f ∈ C 1 (K ). Entonces f es convexa si, y sólo si f (y) ≥ f (x) + ∇f (x)t (y − x), para todo x, y ∈ K. R

n

Prueba. Sea f convexa como en la hipótesis del teorema, x, y ∈ K, λ ∈ [0, 1] . Luego, f (λy + (1 − λ)x) ≤ λf (y) + (1 − λ)f (x). Por lo tanto, f (x + λ(y − x)) − f (x) ≤ λ(f (y) − f (x))

Entonces

f (x + λ(y − x)) − f (x) ≤ f (y) − f (x) λ→0 λ l´ım

42


Luego, ∇f (x)t (y − x) ≤ f (y) − f (x)

Por lo que tenemos para todo x, y ∈ K.

f (x) + ∇f (x)t (y − x) ≤ f (y)

Reciprocamente, si f (y)

≥ f (x) + ∇f (x)t (y − x), para

todo x, y λy + (1 − λ)x, ∀ λ ∈ [0, 1] y también y = x o y = y tenemos:

∈ K, al

hacer x = zλ =

f (x) ≥ f (zλ ) + ∇f (zλ )t (x − zλ ) f (y) ≥ f (zλ ) + ∇f (zλ )t (y − zλ )

Multiplicando por (1 − λ) a la primera y por λ al segundo, y sumando se tiene (1 − λ)f (x) + λf (y) ≥ (1 − λ)(f (zλ ) + ∇f (zλ )t (x − zλ )) +λ(f (zλ ) + ∇f (zλ )t (y − zλ )) = f (zλ ) + ∇f (zλ )t (x − zλ − λx + λzλ + λy − λzλ ) = f (zλ ) + ∇f (zλ )t (λy + (1 − λ)x − zλ ) = f ((1 − λ)x + λy).

Teorema 2.2 Si el problema de minimización con restricciones de igualdad y desigualdad es un problema de programación convexa, y en x ∗ son válidas las condiciones generales de KKT, entonces x∗ es minimizador global.

6 x∗ . Si Prueba. Definamos Ω = {x ∈ Rn : h(x) = 0, g(x) ≤ 0} y tomamos x ∈ Ω con x = λ ∈ Rn y μ ∈ R p son los multiplicadores lagrangianos respectivamente, entonces m

∗

∇f (x ) +

X

λi ∇hi (x∗ ) + μi ∇gi (x∗ ) = 0

(2.1)

i=1

h(x∗ ) = 0 μi gi (x∗) = 0, μi ≥ 0, gi (x∗) ≤ 0,

(2.2)

i = 1,...,p

i = 1,...,p i = 1,...,p

(2.3) (2.4)

(2.5)

43

2.1. PUNTOS INTERIOR BASADO EN BARRERA LOGARÍTMICA

Ahora, (2.2), (2.4) y (2.5) implican m

f (x) ≥ f (x) +

P

X

λi hi (x) +

i=1

X

μi gi (x)

(2.6)

i=1

ya que hi(x) = 0 para i = 1, ..., m, g i(x) ≤ 0, i = 1,...,p, y vale (2.4). Aplicando la desigualdad del teorema anterior sobre f, hi y gi , vemos que f (x) ≥ f (x∗ ) + ∇f (x∗)t (x − x∗ ) hi (x) = hi (x∗ ) + ∇hi (x∗ )t (x − x∗ ) i = 1,...,m gi (x) ≥ gi (x∗ ) + ∇gi (x∗ )t (x − x∗ ) i = 1,...,p

Reemplazando en (2.6) m

∗ t

∗

∗

f (x) ≥ f (x ) + ∇f (x ) (x − x ) +

X

λi (hi (x∗) + ∇hi (x∗ )t (x − x∗ ))

i=1

P

X

μi (gi (x∗ ) + ∇gi (x∗ )t (x − x∗)).

+

i=1

Usando las condiciones (2.1) - (2.5) vemos que f (x) ≥ f (x∗ ). Por lo tanto, x∗ es minimizador global de un problema convexo.

2.1. Puntos Interior Basado en Barrera Logarítmica Se denomina métodos de puntos interiores precisamente porque los puntos generados por estos algoritmos se hallan en el interior de la región factible. Esta es una clara diferencia respecto del método Simplex, el cual avanza por la frontera de dicha región moviéndose de un punto extremo a otro Un punto interior para problemas de programación lineal en la forma estándar en un vector que satisface las restricciones de igualdad, pero sus componentes deben ser estrictamente positivas. Si Ω = {x ∈ Rn : Ax = b, x ≥ 0} denota la región de factibilidad para un problema de programación lineal en la forna estándar, entonces ,nosotro decimos que el punto x˙ es un punto interior de Ω, o que se encuentra en el interior de Ω, si x˙ ∈ Ω0 = {x ∈ Rn : Ax = b, x ≥ 0} . La construcción de este método está soportada por tres técnicas numéricas muy utilizadas en optimización: Penalización interna, Lagrange y el método de Newton. Un Newton para minimización irrestricta, y los métodos de Penalización interna y de Lagrange para convertir problemas con restricciones en problemas irrestrictos..

44


2.1.1. Construcción del método Consideremos los problemas Primal y Dual para programación lineal en forma estándar Minimizar ctx Sujeto a:

(P)

Ax = b x≥0

Donde c, x ∈ Rn, b ∈ Rm y A ∈ Rm

n

×

(m < n)

Maximizar bty Sujeto a:

(D)

At y + z = c z≥0

Donde z

∈ Rn , y ∈ Rm .irrestricta.

Pero debemos tener en cuenta lo siguiente

H1. El conjunto P = {x ∈ Rn : Ax = b,

x > 0} no vacío

H2. El conjunto T = {y ∈ Rm , z ∈ Rn : At y + z = c,

z > 0} no vacío

H3. La matriz A es de rango m En estas condiciones, de acuerdo con el teorema de la adualidad, los problemas (P) y (D) tienen solución óptima, coincidiendo los valores de sus funciones objetivo Además el conjunto de soluciones óptimas de ambos problemas están acotados.. La estratégia de este método resuelve ambos problemas en forma simultánea. Para x factible y (y, z), el duality gap (intervalo de dualidad) duality gap = ct x − bt y = xt c − xt At y

(Ax = b)

= xt (c − At y) = xt z

(c − At y = z)

≥ 0

(x, z ≥ 0)

duality gap es cerrado en el óptimo x∗ ctx∗ − bt y ∗ = (x∗ )t z ∗ = 0

45


por lo que las condiciones de complementaridad son: x j∗ z j∗ = 0,

j = 1, 2,...,n

Las condiciones de no negatividad serán manejadas por el método de Ba-rrera logaritmica, y las restricciones de igualdad por Lagrange. La función irrestricta resultante por el método de Newton.

Introduciendo primeramente la función Barrera a los (P) y (D) Minimizar ct x − μ

n

P

log(x j )

(2.7)

j =1

Ax − b = 0 n

P P

Maximizar b y + μ t

log(z j )

(2.8)

j =1

At y + z − c = 0

Vale recalcar que las barreras,

n

−

P

log(x j ) y +

j =1

n

log(z j ), no permiten que los puntos in-

j =1

teriores de (P) y (D) se aproximen a la frontera de la región factible. Las barreras al estar presentes en la función objetivo del problema de minimización aumentan la función objetivo infinitamente a medida que el punto x se acerca a la frontera. Pero en programación lineal el óptimo x ∗ esta en un extremo, y aparentemente la barrera nos impediría de llegar a x ∗ . Para contrarrestar este efecto usamos un parámetro de penalización equilibrante μ > 0.

A continuación formando las funciones lagrangianas correspondiente de (2.7) y (2.8) tenemos: n t LP (x, y) = c x − μ log(x j ) − y t(Ax − b) (2.9)

X X j =1

n

t

LD (x, y) = b y +

log(z j ) − xt (Aty + z − c)

(2.10)

j =1

Calculando las derivadas parciales respectivas de (2.9) y (2.10), tenemos: ∂L P μ = c− − At y = 0 ∂x x j ∂L P = b − Ax = 0 ∂y

(2.11) (2.12)

46

CAPÍTULO 2. PUNTOS INTERIOR PARA PROGRAMACIÓN LINEAL ∂L D = At y + z − c = 0 ∂x ∂L D = b − Ax = 0 ∂y ∂L D μ = − x = 0; −→ x j z j = μ ∂z z j

(2.13)

(2.14)

(2.15)

De (2.13) despejamos c = At y + z y sustituimos en (2.11) se tiene At y + z −

μ t −A y = 0 x j x j z j = μ

j = 1,...,n

estas expresiones después de simpli ficar las expresiones algebraicas tenemo Ax = b At y + z = c x j z j = μ,

(2.16) j = 1,...n

Se observa que cuando μ = 0, las condiciones de OPtimalidad son satisfechas, pero por cuestiones numéricas debemos fi jar un μ > 0 y aproximarlo iterativamente hacia cero

Dados los puntos x, z ∈ Rn. definamos e = [1 1 1 ... 1]t ∈ Rn (n-vector de unos) y

⎡x ⎢⎢ 0 X = ⎢ ⎢⎣ ...

0 x2

· ·· · ··

0 0

0

0

· ··

xn

⎡ ⎢⎢ =⎢ ⎢⎣

1

0

0 0

xn

1

.. .

...

.. .

⎤ ⎥⎥ ⎥⎥ ⎦

⎡z ⎢⎢ 0 Z = ⎢ ⎢⎣ ...

1

y

0 z2

· ·· · ··

0 0

0

· ··

zn

.. .

0

...

.. .

⎤ ⎥⎥ ⎥⎥ ⎦

notando que

−1

X

0

.. .

x2

· ·· · ··

0

0

· ··

x1

1

.. .

...

.. .

1

⎤ ⎥⎥ ⎥⎥ ⎦

y

−1

Z

⎡ ⎢⎢ =⎢ ⎢⎣

1 z1

0

0

.. .

z2

· ·· · ··

0

0

· ··

1

.. .

...

0 0

zn

.. .

1

⎤ ⎥⎥ ⎥⎥ ⎦

Es decir, las condiciones de optimalidad simultáneas, para un mismo μ de los problemas (P) y (D) están dadas por: Ax = b At y + z = c XZ e = μe

(2.17)

47


Si la solución del problema del sistema de ecuaciones lineales (2.17) se designa mediante [x(μ), y(μ), z(μ)]t , para cada μ > 0. El duality gap es g(μ) = ct x(μ) − bt y(μ) = xt (μ)c − xt (μ)Ay(μ)

(Ax(μ) = b)

= xt (μ) (c − Ay(μ))

¡

= xt (μ)z(μ)

¢

c − Aty(μ) = z(μ)

= μ

Al tender μ → 0, g (μ) converge a cero, lo que implica que x(μ) y y(μ) converge a la solución de los problemas (P) y (D), respectivamente.

μ es una aproximación del intervalo de dualidad, y debemos aclarar, que si su valor de-

crece muy rápidamente se puede perder convergencia. El intervalo de dualidad puede incluso aumentar. Si su valor disminuye muy lentamente se requiere muchas iteraciones.

Pero, por cuestiones numéricas, nosotros debemos comenzar fi jando un μ > 0 inicial y aproximarlo iterativamente hacia cero.

Definamos la trayectoria central como: C = {(xμ , yμ , zμ ); μ > 0}

La trayectoria central C es un arco de puntos estrctamente factibles, que juegan un rol importante en la teoria de algoritmos Primal-Dual. Para cada μ > 0 se tiene los puntos (xμ, yμ , zμ ) ∈ C.

Otra manera de definir C es

⎛ ⎜ F (x , y , z ) = ⎝ μ

μ

μ

0 0 XZe -μe

⎞ ⎟⎠ ,

(xμ , zμ ) > 0

48


Se muestra de manera intuitiva la trayectoria central

*

x 0 x µ Figura 2.2: Trayectoria central El método de puntos interiores Primal-Dual, encuentra soluciones Óptimas (x∗, y∗ , z∗ ), aplicando variantes del método de Newton a las tres condiciones de igualdad (2.17) y modi ficando las direcciones de búsqueda y longitud de paso, de manera que las desigualdades (x, z) ≥ 0 sean satisfechas estrictamente en cada iteración.

Cálculo de la dirección de movimiento (Paso de Newton) El método Primal-Dual se basa en tomar como direcciones de movimiento las de Newton, resultantes de resolver el sistema (2.17). Dada los puntos iniciales x0 , z 0 ∈ Rn, y0 ∈ Rm, donde x0 > 0 y z 0 > 0. Fijemos adecuadamente un valor inicial para el parámetro de penalización μ = μ0 . El sistema (2.17) lo podemos respresentar por F (x0 , y0 , z 0 ) = 0.

⎡ Ax − b ⎤ ⎢ ⎥ F (x , y , z ) = ⎣ A y + z − c ⎦ = 0 0

0

0

0

t 0

0

XZe − μ0e

La matriz Jacobiana de F (x0 , y0 , z 0) es

⎡A ⎢ J (x , y , z ) = ⎣ 0 0

0

0

Z

0

0

At

I X

0

⎤ ⎥⎦

(2.18)

49


para realizar una iteración de Newton, resolvamos el sistema J (x0 , y0 , z 0)dw = −F (x0, y0 , z 0 ) donde dw = [dx,dy,dz]t es el paso de Newton.

⎡A ⎢⎣ 0

Z

0

0

At

I X

0

⎤ ⎡ dx ⎤ ⎥⎦ ⎢⎣ dy ⎥⎦

⎡ Ax − b ⎤ ⎢ ⎥ −⎣ A y + z − c ⎦ XZe − μ e ⎡ b − Ax ⎤ ⎢⎣ c − A y − z ⎥⎦ 0

=

dz

t 0

0

0

0

=

t 0

0

μ0 e − XZe

definamos los vectores residuales Primal y Dual r p = b − Ax0 rd = c − At y0 − z 0 rc = μ0 e − XZe

Es decir, se resuelve el sistema 0

⎡A ⎢⎣ 0

At

Z

Obsérvamos que si x

k

0

0

£ ¤

∈ P y yk , z k

t

⎤ ⎡ dx ⎤ ⎡ r ⎤ ⎥⎢ ⎥ ⎢ ⎥ I ⎦ ⎣ dy ⎦ = ⎣ r ⎦ p

(2.19)

d

X

dz

rc

∈ T. k = 0, 1, 2, ...., entonces r p = 0 y r d = 0.

El sistema (2.19) puede ser resuelto indirectamente, y puede ser visto como Adx = r p

(2.20)

At dy + dz = rd

(2.21)

Zdx + Xdz = rc

(2.22)

multiplicando (2.21) por X y restando con (2.22), obtenemos XA t dy − Zdx = X rd − rc

(2.23)

multiplicando la ecuación (2.23) por AZ −1 e involucrando (2.20) llegamos a un sistema cuadrado para dy esto es AZ −1 XA t dy = AZ −1 (Xr d − rc ) + r p

calculada dy se puede calcular fácilmente dz y dx,de (2.21) y (2.22), en conjunto tenemos.

50

CAPÍTULO 2. PUNTOS INTERIOR PARA PROGRAMACIÓN LINEAL dy = (AZ −1XAt )−1 (AZ −1 (Xrd − rc ) + r p) dz = rd − At dy dx = Z −1 (rc − Xdz)

(2.24)

Amplitud de movimiento (Longitud de paso) El algoritmo básico seguido por el método de Puntos interiores consiste en, comenzando por xo , generar una sucesión de puntos { xk } que nos lleve hasta el punto óptimo x ∗ Cada nuevo punto se obtendrá del anterior, siguiendo el procedimiento iterativo. x(k+1) = x(k) + βαdx

Mediante (2.24) tenemos calculada el paso de Newton (dx,dy,dz)t . Ahora usamos dx como una dirección en el x−espacio y (dy,dz) como una dirección en el (y, z)−espacio. Seguidamente debemos elegir la longitud de los pasos en los dos espacios de modo que podamos desplazarnos prudentemente, manteniendo x > 0, z > 0, esto es posible utilizando: αP =

m´ın 1≤ j ≤n

αD =

m´ın 1≤ j ≤n

½ ½

−x j

: dx j −z j : dz j

¾ ¾

dx j < 0 dz j < 0

Ahora podemos dezplazarnos prudentemente, el siguiente punto (x1, y1 , z 1 ) será calculado mediante x1 = x0 + βα P dx y 1 = y0 + βα D dy z 1 = z 0 + βα D dz

donde β ∈ h0, 1i es un parámetro que nos permite controlar que las posibles soluciones no lleguen a la frontera de la región factible, ni mucho menos salgan de la región. En la práctica se obtuvieron buenos resultados usándo β = 0.95 y/o β = 0.98, esto completa una iteración de Newton. La técnica sólo resuelve el sistema (2.17) con respecto a un μ0 , para obtener una aproximación a la solución óptima, deberíamos reducir μ0 , en cada iteración. Los puntos (x1 , y1 , z 1 )t son utilizados como iniciales para la segunda iteración de newton, t al repetir esto, el método genera una sucesión xk , yk , z k que converge sólo hacia una solución del sistema (2.17).

n¡

Ajuste del Parámetro Penalizador (μ)

¢o


51

Ahora debemos establecer un criterio para determinar cuando el punto actual x(k) está lo suficientemente próximo a x∗. Penalización sugiere la disminución iterativa de μ con la finalidad de obtener una solución para el problema de programación lineal estándar, y nada impide la disminución del parámetro de penalización en cada iteración de newton, dado μ k , el siguiente parámetro penalizador es elegido mediante:

¯

t k

t k

¯

c x −b y , n2

μk+1 =

k = 0, 1, 2, ..

donde n es el número de variables, esto produce una reducción substancial en μ a cada paso de newton, las iteraciones continúan hasta que la diferencia

¯

t k

t k

c x −b y

¯

sea suficientemente pequeña, esto se consigue controlando el error relativo

¯

¯

ct xk − bt yk <ε 1 + |bt yk |

¯ ¯

(2.25)

de (2.25), el denominador 1 + bt yk nos permite mantener una proporción acerca de la aproximación entre los valores objetivos primal y dual, se le sumó 1 para evitar problemas de precisión cuando el costo óptimo se acerca a cero. Por ejemplo, si para cada cierta iteración k, ctxk = 100000001 y bt y k = 100000000, el error absoluto es ct xk − bty k = 1. Sin embargo, el x −b y error relativo c 1+ tendrá una aproximación con una precisión de 8 dígitos, con lo cual ya b y es razonable parar el algorímo. Donde ε > 0 es la precisión deseada. En la práctica se utiliza valores ε = [10−6 , 10−8 ] . t k

t k

t k

2.1.2. El Algoritmo Inicialización: Dados y ∈ Rm y x, z ∈ Rn, donde x, z > 0 . Definir la precisión requerida ε > 0, el n−vector columna e = [1 1..,1]t y el valor inicial de μ > 0. c x( ) −b y ( ) Paso 1: Mientras | 1+|b y( )| | < ε, no satisface, HACER t

k

t

t

k

k

Paso 2: Definir las matrices diagonales X = diag(x) y Z = diag(z)

⎧⎪ r = b − Ax ⎨ = c − A y − z ⎪⎩ rr = μ e − XZ ⎧⎪ dy = (AZ XA ) (AZ ⎨ = r − A dy ⎪⎩ dz dx = Z (r − Xdz) p

Paso 3: Calcular

t

d c

0

−1

Paso 4: Calcular

d

−1

t −1

t

c

−1

(Xrd − rc ) + r p )

52


Paso 5: Calcular

αP = min

1≤ j ≤n

αD = min

1≤ j ≤n

n n

−xj

o o

: dx j < 0

dxj −zj

: dz j < 0

dzj

⎧⎪ x ← x + (0.95)(α ⎨ y + (0.95)(α ⎪⎩ yz ← ← z + (0 .95)(α

)dx D )dy D )dz P

Paso 6: Asignar

Paso 7: Hacer

μk+1 =

|ct x−bt y| n2

y volver al paso 1.

2.1.3. Implementación y Experimentos La implemetación está hecha de un modo particular y con fines experimentales. hemos utilizado Matlab 6.0, esto se justi fica debido a la interface grá fica y matricial que Matlab otorga

Método de puntos interiores Primal-Dual function [x,v,iter]=MPInterior(A,b,c) [m,n]=size(A); x=ones(n,1); y=ones(m,1); z=x; u=100; e=x iter=0; while abs(c’*x-b’*y)/(1+abs(b’*y)) > 0.000001; X=diag(x); Z=diag(z); rp=b-A*x; rd=c-A’*y+z; rc=u*e-X*Z*e; dy=(A*inv(Z)*X*A’)\(A*inv(Z)*(X*rd-rc)+rp); dz=rd-A’*dy; dx=Z\(rc-X*dz); aP=-1/min(min(dx./x),-1); aD=-1/min(min(dz./z),-1); x=x+0.95*aP*dx;

y=y+0.95*aD*dy; z=z+0.95*aD*dz;

u=abs(c‘*x-b’*y)/n^2;

iter=iter+1;

end v=c’*x;

El programa resuelve problemas de programación lineal en la forma estándar . Minimizar c t x Ax = b x≥0

debemos advertir que αP =aP y αD =aD, ambos ejecutan el criterio de la razón. El parámetro con que inicimos es μ =100, para reducir μ a cada iteración se usó u = abs(c0 ∗ x − b0 ∗ y)/nˆ2.

53


Con un parámetro de precisión de ε = 0,000001, y los puntos iniciales x = [1 1 . . .1]t, y = [0 0 . . .0]t y z = [1 1 . . .1]t , en la fila siete se tiene dx=Z\(rc-X*dz); que es equivalente a Zdx=(rc-X*dz). Una vez ingresada la matriz A y los vectores b y c, llamamos MPInterior

Ejemplo 2.1 Considere el siguiente ejemplo de programación lineal Minimizar x1 + x2 + x3 x1 + 2x2 + 3x3 ≥ 10 2x1 + x2 − x3 ≥ 20 x3 ≥ 4 x1 , x2 , x3 ≥ 0 llevando a su forma estándar tenemos Minimizar x1 + x2 + x3 x1 2x1

+2x2 +x2

+3x3 −x4 −x3 −x5 x3 −x6 x1 , x2 , x3, x4, x5 , x6 ≥ 0

⎛1 ⎜ A = ⎝ 2

2 3 -1 0 0 1 -1 0 -1 0 0 0 1 0 0 -1

⎞ ⎟⎠

,

⎛ 10 ⎞ ⎜ ⎟ b = ⎝ 20 ⎠

= 10 = 20 =4

,

4

⎛ ⎞ ⎜⎜ 11 ⎟⎟ ⎜⎜ ⎟⎟ ⎜1⎟ c = ⎜ ⎟ ⎜⎜ 0 ⎟⎟ ⎜⎝ 0 ⎟⎠ 0

La solución factible para este problema dado por (puntos iniciales):

x(0)

⎛ ⎞ ⎜⎜ 305 ⎟⎟ ⎜⎜ ⎟⎟ ⎜ 15 ⎟⎟ , =⎜ ⎜⎜ 75 ⎟⎟ ⎜⎝ 30 ⎟⎠ 11

y (0)

⎡ 1/4 ⎤ ⎢ ⎥ = ⎣ 1/4 ⎦ , 1/4

z (0)

⎛ ⎞ ⎜⎜ 1/4 ⎟⎟ 1/4 ⎜⎜ ⎟⎟ ⎜ 1/4 ⎟⎟ =⎜ ⎜⎜ 1/4 ⎟⎟ ⎜⎝ 1/4 ⎟⎠ 1/4

con estos vectores iniciales se procede a aplicar el algoritmo para una i-teración. Iteración 1

dado el vector de solución, x(0) y z (0), y las matrices, X y Z, para la primera iteración

54


están disponibles a través de:

⎡ ⎢⎢ 300 ⎢⎢ ⎢0 X = ⎢ ⎢⎢ 0 ⎢⎣ 0

0 5

0

0 0 0

0 0 0 0

0 0 0 0 0

0 15 0 0 75 0 0 0 30 0 0 0 0 11

0

⎡ ⎢⎢ 1/4 ⎢⎢ 0 ⎢ 0 Z = ⎢ ⎢⎢ 0 ⎢⎣ 0

0 0

0

0 0

1/4

0 0 0 0

0 0 0

1/4

0 0 0

0 0 0 0

1/4

0 0

⎤ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎦

1/4

0 0 0 0 0

0

1/4

⎤ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎦

⎡ ⎤ ⎢⎢ 11 ⎥⎥ ⎢⎢ ⎥⎥ ⎢1⎥ e = ⎢ ⎥ ⎢⎢ 1 ⎥⎥ ⎢⎣ 1 ⎥⎦ 1

Las funciones objetivo del primal y dual, así como la brecha de dualidad ( duality gap ) y la barrera de parámetros μ para estos vectores de inicio están dados por ct x = 50.00 ; bt y = 8.50;

MIENTRAS

¯

¯

ct x − bt y = 41.50

|ct x − bt y| = 4,3684 < ε, HACER 1 + |bt y|

calculando dy =

¡

AZ −1 XA t

dz = rd − At dy

¢

−1

dx = Z −1 (rc − Xdz)

(AZ −1 (Xr d − rc ) + r p )


55

Pues r p = b − Ax(0)

⎡0⎤ ⎢ ⎥ =⎣ 0 ⎦ 0

rd = c − At y(0) − z (0)

⎡ ⎤ ⎢⎢ 00 ⎥⎥ ⎢⎢ ⎥⎥ ⎢0⎥ =⎢ ⎥ ⎢⎢ 0 ⎥⎥ ⎢⎣ 0 ⎥⎦ 0

tenemos

⎡ -2.83 ⎤ ⎢ ⎥ dy = ⎣ 0.31 ⎦ 5.11

;

⎡ ⎤ ⎢⎢ 2.21 ⎥⎥ 5.35 ⎢⎢ ⎥⎥ ⎢ 3.68 ⎥⎥ dz = ⎢ ⎢⎢ -2.83 ⎥⎥ ⎢⎣ 0.31 ⎥⎦ 5.11

;

⎡ ⎤ ⎢⎢ 104.59 ⎥⎥ 288.06 ⎢⎢ ⎥⎥ ⎢ 164.17 ⎥⎥ dx = ⎢ ⎢⎢ 1173.23 ⎥⎥ ⎢⎣ 333.06 ⎥⎦ 164.17

Usando el criterio de la razón, se encuentra el tamaño del paso del vector primal

½

¾

x j (k) α p = min : ∀ dx j (k) < 0, 1 ≤ j ≤ n −dx j (k) α p = m´ın {1.00} = 1.00

de igual manera para el vector dual.

½

¾

z j (k) αd = min : ∀ dz j (k) < 0, 1 ≤ j ≤ n −dz j (k)

αd = m´ın {0.09} = 0.09

Usando un factor de tamaño de paso del 95%, para las nuevas iteraciones al final de la primera iteración, para el vector primal x, el vector dual y, y el vector reducción de costo z

56


se obtiene

x(1)

y(1)

⎛ ⎞ ⎡ ⎤ ⎛ ⎞ 129.36 ⎜ ⎟⎟ ⎜⎜305⎟⎟ ⎢⎢ 104.59 ⎥ ⎥ ⎜ 288.06 278.66 ⎜⎜ ⎟⎟ ⎢⎢ ⎥⎥ ⎜⎜ ⎟⎟ ⎜15⎟ ⎢ 164.17 ⎥⎥ = ⎜⎜ 170.97 ⎟⎟ = ⎜ ⎟ + (0.95) (1) ⎢ ⎜⎜75⎟⎟ ⎢⎢ 1173.23 ⎥⎥ ⎜⎜1189.57⎟⎟ ⎜⎝30⎟⎠ ⎢⎣ 333.06 ⎥⎦ ⎜⎝ 346.41 ⎟⎠ 11 164.17 166.97 ⎛0.25⎞ ⎡ -2.83 ⎤ ⎛0.01⎞ ⎜ ⎟ ⎢ ⎥ ⎜ ⎟ = ⎝0.25⎠ + (0.95) (0.09) ⎣ 0.31 ⎦ = ⎝0.28⎠ 0.25

z (1)

5.11

0.68

⎛ ⎞ ⎡ ⎤ ⎛ ⎞ 104.59 0.44 ⎜⎜0.25 ⎟ ⎢ ⎥ ⎜ ⎟⎟ ⎟ ⎢ ⎥ ⎜ 0.25 288.06 0.70 ⎜⎜ ⎟⎟ ⎢⎢ ⎥⎥ ⎜⎜ ⎟⎟ ⎟⎟ + (0.98) (0.09) ⎢⎢ 164.17 ⎥⎥ = ⎜⎜0.56⎟⎟ =⎜ ⎜⎜0.25 ⎢⎢ 1173.23 ⎥⎥ ⎜⎜0.01⎟⎟ ⎜⎜0.25⎟⎟⎟ ⎢⎣ 333.06 ⎥⎦ ⎜⎝0.28⎟⎠ ⎝0.25⎠ 0.25

164.17

0.68

reduciendo el valor del parámetro u ct x =578.98 , bt y = 8.36 μ(1) = 0,1 ∗

|ct x−bt y| n2

=0.1∗ 570.62 36 =1.59

Con las nuevas iteraciones disponibles para el vector de costo reducido z, y el vector dual y, un nuevo ciclo de iteración se puede empezar. Las primeras diez iteraciones se muestran en las siguientes tablas. El resumen para el vector dual y , y el valor objetivo del dual se muestra en la tabla 1; para el vector solución del primal y el valor objetivo del primal se muestra en la tabla 2; y para el vector de costo reducido z, y el parámetro de barrera, μ, se muestra en la tabla 3.

57


iteración

y1

y3 Objetivo

y2

1

0.01

0.28

0.68

8.36

2

0.01

0.28

0.69

8.53

3

0.01

0.37

0.84

10.89

4

0.01

0.49

1.23

14.79

5

0.00

0.50

1.44

15.77

6

0.00

0.50

1.50

15.98

7

0.00

0.50

1.50

16.00

8

0.00

0.50

1.50

16.00

9

0.00

0.50

1.50

16.00

10

0.00

0.50

1.50

16.00

Tabla 1

iter.

objt.

x1

x2

x3

x4

x5

x6

1

129.36

278.66

170.97

1189.57

346.41

166.97

578.98

2

14.88

21.41

13.84

89.24

17.32

9.84

50.13

3

8.99

8.02

5.13

30.42

0.87

1.13

22.14

4

10.00

4.22

4.06

20.62

0.17

0.06

18.28

5

11.78

0.47

4.00

14.04

0.00

0.00

16.02

6

11.99

0.00

4.00

14.00

0.00

0.00

16.00

7

12.00

0.00

4.00

14.00

0.00

0.00

16.00

8

12.00

0.00

4.00

14.00

0.00

0.00

16.00

9

12.00

0.00

4.00

14.00

0.00

0.00

16.00

10

12.00

0.00

4.00

14.00

0.00

0.00

16.00

Tabla 2

58


iter

z1

z2

z3

z4

z5

z6

μ

1

0.44

0.70

0.56

0.01

0.28

0.68

1.59

2

0.42

0.69

0.55

0.01

0.28

0.69

0.12

3

0.25

0.60

0.49

0.01

0.37

0.84

0.03

4

0.01

0.49

0.23

0.01

0.49

1.23

0.01

5

0.00

0.49

0.04

0.00

0.50

1.44

0.00

6

0.00

0.50

0.00

0.00

0.50

1.50

0.00

7

0.00

0.50

0.00

0.00

0.50

1.50

0.00

8

0.00

0.50

0.00

0.00

0.50

1.50

0.00

9

0.00

0.50

0.00

0.00

0.50

1.50

0.00

10

0.00

0.50

0.00

0.00

0.500

1.50

0.00

Tabla 3

de la tabla 2, después de 10 iteraciones se tiene el punto óptimo y el valor óptimo x∗ = [12 0 4 14 0 0]

0

v = 16

Ejemplo 2.2 Considere el siguiente problema Minimizar −2x1 − 7x2 s.a 4x1 + 5x2 ≤ 40 2x1 + x2 ≥ 8 2x1 + 5x2 ≥ 20 x1 , x2 ≥ 0

59

2.1. PUNTOS INTERIOR BASADO EN BARRERA LOGARÍTMICA su forma estándar, dada por Minimizar −2x1 − 7x2 + 0x3 − 0x4 − 0x5 s.a 4x1 + 5x2 + x3 = 40 2x1 + x2 =8 −x4 2x1 + 5x2 −x5 = 20 x1, x2 , x3 , x4 , x5 ≥ 0 y su correspondiente Dual del problema es Maximizar s.a

40y1 + 8y2 + 20y3 4y1 + 2y2 + 2y3 + z1 5y1 + y2 + 5y3 +z2 y1 +z3 +z4 −y2 +z5 −y3 z1 , z2 , z3 , z4 , z5 ≥ 0

Comencemos tomando un punto inicial factitible

x0

= −2 = −7 =0 =0 =0

⎡5⎤ ⎢⎢ 3 ⎥⎥ ⎢ 5 ⎥⎥ , =⎢ ⎢⎢ ⎥⎥ ⎣5⎦ 5

y0

⎡ -4 ⎤ ⎢ ⎥ =⎣ 1 ⎦yz 1

0

⎡ 10 ⎤ ⎢⎢ 7 ⎥⎥ ⎢ 4 ⎥⎥ en esta condición también se puede observar que r =⎢ ⎢⎢ ⎥⎥ ⎣1⎦

p

=0y

1

rd = 0.

A partir de esta solución factible, la figura adjunta muestra la sucesión generada por el

60


algoritmo hasta llegar a la solución óptima factible, el punto [0 8]t .

8

7

6

5

4

3

2

1

0 0

1

2

3

4

5

6

7

8

Punto inicial factible [5 3]t

iter 1 2 3 4 5 6

x0 x1 x2 x3 x4 x5 x6

= = = = = = =

5

3

2.35 6.07 2.27 6.12 1.18 7.02 0.12 7.90 0.01 7.99 0.00 8.00

Observamos ahora que, dada otro punto inicial infactible x 0

⎡1⎤ ⎢⎢ 1 ⎥⎥ ⎢ 1 ⎥⎥ , =⎢ ⎢⎢ ⎥⎥ ⎣1⎦ 1

9

10

61


y0

⎡1⎤ ⎢ ⎥ = ⎣ 1 ⎦ y z 1

0

⎡1⎤ ⎢⎢ 1 ⎥⎥ ⎢ 1 ⎥⎥ , Aquí =⎢ ⎢⎢ ⎥⎥ ⎣1⎦

= 0 y r d = 6 0, pero en el punto óptima se observó que r p 6

1 estos residuos son r p = 0 y r d = 0.

Se observó aquí, a pesar que el punto inicial es infactible [1 1]t , también el algoritmo hace converge a la misma solución óptima del problema la cual es [0 8]t , veamos la figura adjunta

8

7

6

5

4

3

2

1

0 0

1

2

3

4

5

6

7

punto inicial infactible [1 1] t

8

9

10

62

CAPÍTULO 2. PUNTOS INTERIOR PARA PROGRAMACIÓN LINEAL iter 1 2 3 4 5 6 7 8

x0 x1 x2 x3 x4 x5 x6 x7 x8

= = = = = = = = =

1

1

0.93 1.77 3.14 5.16 2.99 5.53 1.13 7.07 0.09 7.93 0.01 8.00 0.00 8.00 0.00 8.00

ser puede observar que x0 y x1 corresponden a puntos aún infactible, mientras que xk , k = 2, 3, 4, ..., están asociados a puntos dentro de la región de factibilidad. La solución óptima en 0 el punto x ∗ = 8

" #

Capítulo 3 Experimentos Computacionales Al ejecutar el programa con los ejemplos siguientes y ver su comportamiento de estos, hemos usado el sistema operativo Windows-XP (Pentium-III), con un PC compatible y usando procesador celeron(R) cpu 2.40 GHz, con 256 MB de memoria RAM. Usamos el Matlab 6.0 para windows, Matlab es el nombre abreaviado de "MATrix LABoratory". MATLAB es un programa para realizar cálculos numéricos con vectores y matrices. Como caso particular puede también trabajar con números escalares, tanto reales como complejos. Una de las capacidades más atractivas es la de realizar una amplia variedad de grá ficos de dos y tres dimensiones. Matlab tiene también un lenguaje de programación propio. Cada resultado obtenido fué comparado con la ejecución del algoritmo Simplex, utilizando la técnica M-grande.

Experimentos Experimento 1: Un carpintero hace sillas y mesas solamente. Las utilidades proyectadas para los dos productos son S/. 20 por silla y S/. 30 por mesa respectivamente. La demanda proyectada es 400 sillas y 100 mesas, cada silla requiere 2 pies cúbicos de madera, mientras que cada mesa requiere 4 pies cúbicos. El carpintero tiene un total de 1 000 pies cúbicos de madera en stock. ¿Cuántas sillas y mesas debe construir para aumentear al máximo sus utilidades? Sean x1 el número de sillas x2 el número de mesas

estas son las variables desconocidas, para éste problema. 63

64

CAPÍTULO 3. EXPERIMENTOS COMPUTACIONALES

El carpintero quiere aumentar al máximo sus utilidades, esto es 20x1 + 30x2

sujeto a las restricciones que:

El importe total de madera para los dos productos no puede exceder los 1 000 pies cúbicos disponibles Los números de sillas y mesas a ser hechas, no deben exceder las demandas,

El número de sillas y mesas necesitan ser no negativos. por lo que tenemos: Maximizar 20x1 + 30x2 2x1 + 4x2 ≤ 1 000 0 ≤ x1 ≤ 400 0 ≤ x2 ≤ 100

añadiendo las variables auxiliares x 3 , x4, x5 y llevando a su forma estandar, tenemos. Minimizar

− 20x1 − 30x2 + 0x3 + 0x4 + 0x5

2x1 + 4x2 + x3 + 0x4 + 0x5 = 1 000 x1 + 0x2 + 0x3 + x4 + 0x5 = 400 0x1 + x2 + 0x3 + 0x4 + x5 = 100 x1 , x2, x3 , x4 , x5 ≥ 0

Al introducir los datos en el computador

⎡2 ⎢ A = ⎣ 1

4 1 0 0 0 0 1 0 0 1 0 0 1

⎤ ⎥⎦

⎡ -20 ⎤ ⎡ 1000 ⎤ ⎢⎢ -30 ⎥⎥ ⎢ ⎥ ⎢⎢ 0 ⎥⎥⎥ b = ⎣ 400 ⎦ ; c = ⎢ ⎢⎣ 0 ⎥⎦ 100 0

65 y ejecutando el programa de punto interior se obtuvo

⎡ 400.00 ⎤ ⎢⎢ 50.00 ⎥⎥ ⎢ 0.00 ⎥⎥ x = ⎢ ⎢⎢ ⎥⎥ ⎣ 0.00 ⎦ 50.00

con valor Óptimo v = -9499.99, como el problema es de maximización v = 9499.99

Experimento 2 Supongamos que un carpintero fabrica tres tipos de puertas en diferentes medidas: pequeñas, medianas y grandes. La utilidad que obtiene por cada una de ellas es de 100, 130 y 150 nuevos soles, respectivamente. El tiempo de elaboración es de 5 h, 7 h y 10 h por cada tipo de puerta en el orden en que se menciona. Además, el consumo em m 2 de madera es de 2,4, 2,8 y 3,2 para cada una de las puertas, en el mismo orden antes mencionado. Finalmente, el carpintero estima que puede vender a lo mucho 3 puertas pequeñas, 5 medianas y 5 grandes semanalmente. Por último, él dispone de 40h a la semana y de 30m2 de madera por semana. Encuentre el modelo de programación lineal que maximice la utilidad del carpintero.

Solución: Sean x1 : numéro de puertas pequeñas x2 : numéro de puertas medianas x3 : numéro de puertas grandes

má x 100x1 + 130x2 + 150x3 5x1 + 7x2 + 10x3 ≤ 40 2,4x1 + 2,8x2 + 3,2x3 ≤ 30 x1 ≤ 3 x2 ≤ 5 x3 ≤ 5 x1 , x2 , x3 ≥ 0

añadiendo las variables auxiliares x3 , x4 , x5 ,x6 y llevando a su forma estandar, tenemos. − m´ın − 100x1 − 130x2 − 150x3 + 0x4 + 0x5 + 0x6 + 0x7 + 0x8

66

CAPÍTULO 3. EXPERIMENTOS COMPUTACIONALES 5x1 + 7x2 + 10x3 + x4 = 40 2,4x1 + 2,8x2 + 3,2x3 + x5 = 30 x1 + x6 = 3 x2 + x7 = 5 x3 + x8 = 5 x1 , x2 , x3 , x4 , x5 , x6, x7 , x8 ≥ 0

⎡ 5 ⎢⎢ 2,4 ⎢ A = ⎢ ⎢⎢ 1 ⎣ 0 0

7 10 1 0 0 2,8 3,2 0 1 0 0 0 0 0 1 1 0 0 0 0 0 1 0 0 0

0 0 0 1 0

0 0 0 0 1

⎤ ⎥⎥ ⎥⎥ ; ⎥⎥ ⎦

⎡ 40 ⎤ ⎢⎢ 30 ⎥⎥ ⎢ ⎥ ⎢ b = ⎢ 3 ⎥ ; ⎢⎣ 5 ⎥⎥⎦ 5

⎡ ⎤ ⎢⎢ −−100 ⎥⎥ 130 ⎢⎢ ⎥⎥ ⎢⎢ −150 ⎥⎥ ⎢ 0 ⎥ ⎢ ⎥⎥ c = ⎢ ⎢⎢ 0 ⎥⎥ ⎢⎢ 0 ⎥⎥ ⎢⎣ 0 ⎥⎦ 0

y ejecutando el programa de punto interior se obtuvo

x(11)

⎡ ⎤ ⎢⎢ 34 ⎥⎥ ⎢⎢ ⎥⎥ ⎢⎢ 0 ⎥⎥ ⎢ 0 ⎥⎥ =⎢ ⎢⎢ 13 ⎥⎥ ⎢⎢ 0 ⎥⎥ ⎢⎢ ⎥⎥ ⎣ 1⎦ 5

Experimento 3 Una compañia produce 2 tipos de fertilizantes: HI-Fosfato y LO-Fosfato. Tres materias primas son utilizadas para manufacturar estes fertilizantes de acuerdo con la tabla. Cantidad disponible Materia Prima HI-fosfato LO-fosfato (En toneladas) 1

2

1

1500

2

1

1

1200

3

1

0

500

Lucro líquido

S/. 15

S/. 10

Observamos en las columnas HI-Fosfato y LO-Fosfato indican las toneladas de materia prima que son necesarios para producir una tonelada del respectivo fertilizante. Por ejemplo,

67 necesitamos de 2 toneladas de la materia 1 para producir 1 tonelada de HI-Fosfato. Se desea saber cuanto de fertilizante debemos producir para obtener una máxima ganancia.

Solución: Sean x1 :

toneladas de HI-fosfato producido

x2 :

toneladas de LO-fosfato producido

má x 15x1 + 10x2 2x1 + x2 ≤ 1500 x1 + x2 ≤ 1200 x1 ≤ 500 x1 , x2 ≥ 0

añadiendo las variables auxiliares x3 , x4 , x5 y llevando a su forma estandar, tenemos. − m´ın −15x1 − 10x2 + 0x3 + 0x4 + 0x5

2x1 + x2 + x3 + 0x4 + 0x5 = 1500 x1 + x2 + 0x3 + x4 + 0x5 = 1200 x1 + 0x2 + 0x3 + 0x4 + x5 = 500 x1 , x2 , x3 , x4 , x5 ≥ 0

⎡2 ⎢ A = ⎣ 1

1 1 0 1 0 1 1 0 0 0

⎡ −15 ⎤ ⎤ ⎡ 1500 ⎤ ⎢⎢ −10 ⎥⎥ 0 ⎢ 0 ⎥⎥ ⎥ ⎢ ⎥ 0 ⎦ ; b = ⎣ 1200 ⎦ ; c = ⎢ ⎢⎢ ⎥⎥ 1 500 ⎣ 0 ⎦ 0

y ejecutando el programa de punto interior se obtuvo

x(7)

⎡ 300 ⎤ ⎢⎢ 900 ⎥⎥ ⎢ 0 ⎥⎥ =⎢ ⎢⎢ ⎥⎥ ⎣ 0 ⎦ 200

Experimento 4 (Un problema de plani ficación) Supongamos que tenemos que escribir un trabajo en un tiempo máximo de 40 horas, el trabajo consta de 200 páginas de texto, 62 tablas y 38 figuras. Para realizar este trabajo

68


fueron consultados 4 trabajadores diferentes (tipeadores), los detalles de la rapidez, tarifa y disponibilidad de cada trabajador son dados en la tabla adjunta. ¿Cuál será la estrategia para distribuir el trabajo entre los cuatro trabajadores de modo que el costo sea lo mínimo posible. Tipo

Traba jador 1

Traba jador 2

Traba jador 3

Traba jador 4

texto

8

0.50

10

0.55

9

0.60

14

0.65

pág/h

sol/pág

pág/h

sol/pág

pág/h

sol/pág

pág/h

sol/pág

12

0.60

11

0.60

10

0.50

11

0.55

tab/h

so l/tab

tab/ h

sol/tab

tab/h

sol/tab

ta b/h

sol/ tab

5

1.20

2

0.90

3

1.00

2

0.80

fi g/h

sol/fig

fi g/h

sol/fig

fi g/h

sol/fig

fig/h

sol/fig

tabla

fi gura

Disponibilidad

14 horas

16 horas

12 horas

10 horas

Para formular matemáticamente este problema, es decir para construir un modelo que represente este problema, debemos primeramente de finir las variables, esto es fundamental, pues la complejidad de la formulación depende de esta de finición. Debemos aclarar que no existe una fórmula explícita para el proceso de modelación, sino que depende de la práctica, sentido común y habilidad.

La variable. Sea xi,j la cantidad de unidades de trabajo del tipo i a ser realiza por el trabajador j . Claramente, i = 1, 2, 3 y j = 1, 2, 3, 4. La función objetivo. Note que 0,5x1,1 +0,6x2,1 +1,2x3,1 representa el costo que demanda el trabajador 1 se le es encargado realizar x1,1 páginas de texto, x2,1 tablas y x3,1 figuras. El mismo razonamiento vale para el trabajador 2, el trabajador 3 y el cuarto trabajador. Asi, el costo total debido a los cuatro trabajadores estará dado por 0,5x1,1 + 0,6x2,1 + 1,2x3,1 + 0,55x1,2 + 0,6x2,2 + 0,9x3,2 +

|

{z | {z

Trabajador 1

} | } |

{z {z

Trabajador 2

} }

+0,6x1,3 + 0,5x2,3 + x3,3 + 0,65x1,4 + 0,55x2,4 + 0,8x3,4 Trabajador 3

Trabajador 4

Esta función (objetivo) debemos minimizar tanto como sea posible.

Las restricciones. Como son 200 páginas de texto, entonces el total de páginas de texto digitadas por los cuatro trabajadores debe satisfacer x1,1 + x1,2 + x1,3 + x1,4 = 200

Como son 62 tablas, la suma total de tablas realizada por los cuatro trabajadores debe satisfacer x2,1 + x2,2 + x2,3 + x2,4 = 62

69 Por otro lado, como son 38 figuras, el total de figuras hechas por los cuatro trabajadores debe satisfacer x3,1 + x3,2 + x3,3 + x3,4 = 38

Además, cada trabajador no debe exceder su disponibilidad. Así, si el primer trabajador escribe 8 páginas de texto en una hora, él realizará una página en 1/8 de hora, entonces digitar x1,1 páginas de texto le tomará 81 x1,1 horas (el mismo razonamiento vale para las tablas y figuras), como su disponibilidad es 14 horas, debemos tener 1 1 1 x1,1 + x2,1 + x3,1 ≤ 14 8 12 5

Análogamente para el segundo trabajador 1 1 1 x1,2 + x2,2 + x3,2 ≤ 16 10 11 2

Para el tercer trabajador

1 1 1 x1,3 + x2,3 + x3,3 ≤ 12 9 10 3

Para el cuarto trabajador

1 1 1 x1,4 + x2,4 + x3,4 ≤ 10 14 11 2

Finalmente, como la ejecución del trabajo no debe exceder las 40 horas, entonces la suma de todas las horas utilizadas por los cuatro trabajadores no debe exceder esta cantidad, es decir 1 1 1 1 1 1 x1,1 + x2,1 + x3,1 + x1,2 + x2,2 + x3,2 + 8 12 5 10 11 2

| |

{z {z

} | } |

xi,j ≥ 0,

i = 1, 2, 3

Trabajador 1

{z {z

Trabajador 2

} }

1 1 1 1 1 1 + x1,3 + x2,3 + x3,3 + x1,4 + x2,4 + x3,4 ≤ 40 9 10 3 14 11 2 Trabajador 3

Trabajador 4

Las condiciones de no negatividad. Observe que las variables no pueden ser negativas, es decir j = 1, 2, 3, 4

La formulación matemática. La función objetivo, las restricciones y las condiciones de no negatividad, producen un modelo particular de problema de programación lineal, una de las representaciones que usaremos en en este ejemplo corresponde a:

Minimizar 0,5x1,1 + 0,6x2,1 + 1,2x3,1 + 0,55x1,2 + 0,6x2,2 + 0,9x3,2 + 0,6x1,3 + 0,5x2,3 + x3,3 + 0,65x1,4 + 0,55x2,4 + 0,8x3,4

70


x1,1 + x1,2 + x1,3 + x1,4 = 200 x2,1 + x2,2 + x2,3 + x2,4 = 62 x3,1 + x3,2 + x3,3 + x3,4 = 38 1 1 x + 12 x2,1 + 51 x3,1 ≤ 14 8 1,1 1 1 x + 11 x2,2 + 21 x3,2 ≤ 16 10 1,2 1 1 x + 10 x2,3 + 31 x3,3 ≤ 12 9 1,3 1 1 x + 11 x2,4 + 21 x3,4 ≤ 10 14 1,4 1 1 1 1 x + 12 x2,1 + 51 x3,1 + 10 x1,2 + 11 x2,2 + 21 x3,2 + 8 1,1 1 1 1 + 19 x1,3 + 10 x2,3 + 31 x3,3 + 14 x1,4 + 11 x2,4 + 21 x3,4 ≤ 40

xi,j ≥ 0,

i = 1, 2, 3

j = 1, 2, 3, 4

Llevando a su forma estándar se resolvió este problema, tenemos:

x(13)

⎡ 88 ⎤ ⎢⎢ 0 ⎥⎥ ⎢⎢ 15 ⎥⎥ ⎢⎢ ⎥⎥ ⎢⎢ 112 ⎥⎥ ⎢⎢ 0 ⎥⎥ ⎢⎢ 0 ⎥⎥ ⎢⎢ 0 ⎥⎥ ⎢⎢ ⎥⎥ ⎢⎢ 62 ⎥⎥ =⎢ ⎢⎢ 17,4 ⎥⎥⎥ ⎢⎢ 0 ⎥⎥ ⎢⎢ 0 ⎥⎥ ⎢⎢ 5,6 ⎥⎥ ⎢⎢ 0 ⎥⎥ ⎢⎢ ⎥⎥ ⎢⎢ 4,8 ⎥⎥ ⎢⎢ 0 ⎥⎥ ⎢⎣ 7,2 ⎥⎦ 0

(3.1)

71 interpretando estos resultados x1,1 = 88 :

Al trabajador 1 se le encarga 88 páginas de texto

x3,1 = 15 :

Al trabajador 1 se le encarga 15 figuras

x1,2 = 112 : Al trabajador 2 se le encarga 112 páginas de texto x2,3 = 62 :

Al trabajador 3 se le encarga 62 tablas

x3,3 = 17,4 : Al trabajador 3 se le encarga 17 figuras x3,4 = 5,6 : Al trabajador 4 se le encarga 6 figuras

Experimento 5 ( un problema de alimentación ) Una fábrica agrícola produce comida para ganado, oveja, y pollos. Esto se hace mezclando los siguientes ingredientes básicos: maíz, cal, soya, y harina de pescado. Esos ingredientes contienen los siguientes nutrientes: vitaminas, proteínas, calcio, y grasa. Los contenidos de los nutrientes en cada kilogramo de ingredientes son resumidos en la tabla. Nutrientes Ingredientes Vitaminas Proteinas Calcio Grasa Maíz 8 10 6 8 Cal 6 5 10 6 Soya 10 12 6 6 Harina de pescado 4 8 6 9 La fábrica debe producir 10, 6 y 8 toneladas de comida para ganado, comida para oveja, y comida para pollo. Debido al desabastecimiento, una suma limitada de los ingredientes está disponible, a saber: 6 toneladas de maíz, 10 toneladas de cal, 4 toneladas de soya, y 5 toneladas de harina de pescado. Los precios por kilogramo de ingredientes son respectivamente 0.20, 0.12, 0.24 y 0.12 dólares. La cantidad de nutrientes mínimo y máximo que son permitidos por kilogramo de comida de ganado, comida de oveja, y comida de pollo, están resumidos a seguir: Nutrientes producto

vitaminas

proteinas

calcio

Grasa

minimo

máximo

minimo

máximo

minimo

máximo

minimo

máximo

co mida-ganado

6

∞

6

∞

7

∞

4

∞

comida-oveja

6

∞

6

∞

6

∞

4

∞

comida-pollo

4

∞

6

∞

6

∞

4

∞

72


El objetivo de esta fábrica es adoptar una política la cual sea la más óptima posible. Es decir, que el costo de producción sea minimizado, determinando la cantidad de kilogramos de los ingredientes básicos para alimentar el ganado, ovino y pollos, y cumplir además con todas las restricciones.

Solución: El modelo matemático correspondiente es m´ın z = 0,20x11 + 0,20x12 + 0,20x13 + 0,12x21 + 0,12x22 + 0,12x23 + 0,24x31 0,24x32 + 0,24x33 + 0,12x41 + 0,12x42 + 0,12x43

Sujeto a: x11 + x21 + x31 + x41 = 1000 x12 + x22 + x32 + x42 = 6000 x13 + x23 + x33 + x43 = 8000 x11 + x12 + x13 ≤ 6000 x21 + x22 + x23 ≤ 10000 x31 + x32 + x33 ≤ 4000 x41 + x42 + x43 ≤ 5000 2x11 + 0x21 + 4x31 − 2x41 ≥ 0 4x11 − x21 + 6x31 + 2x41 ≥ 0 −x11 + 3x21 − x31 − x41 ≥ 0 4x11 + 2x21 + 2x31 + 5x41 ≥ 0 2x12 + 0x22 + 4x32 − 2x42 ≥ 0 4x12 − x22 + 6x32 + 2x42 ≥ 0 2x12 + 4x22 + 0x32 + 0x42 ≥ 0 4x12 + 2x22 + 2x32 + 5x42 ≥ 0 4x13 + 2x23 + 6x33 + 0x43 ≥ 0 x13 − 2x23 + 6x33 + 2x43 ≥ 0 0x13 + 4x23 + 0x33 + 0x43 ≥ 0 4x13 + 2x23 + 2x33 + 5x43 ≥ 0

Cuya solución óptima aproximada fué x∗ = [3894,2 4273,65 1313,187 518,6 662,6 4077,8240 873,3 386,256 1443,17 1648,5 813,5 4094,8]t

Interpretando x1,1 = 3894,2. x1,2 = 4273,65 y x1,3 = 1313,187 son los kilogramos de maíz para ser incluidos en los alimentos para ganado, ovejas y pollos, respectivamente. Los demás valores se pueden interpretar análogamente.

Observación 3.1 A partir del experimento 6, para generar problemas inmensos se construyó

73 un programa en Matlab, a lo que llamamos generaPL_1 al ingresar las dimensiones m y n, el programa produce un ejemplar de programación lineal.

function [A,b,c]=generaPL_1(m,n) A=round(100*rand(m,n)-100*rand(m,n)); x=round(100*rand(n,1)); b=A*x; c=round(100*rand(n,1)-100*rand(n,1));

Experimento 6: Aquí se generó un problema de programación lineal con 90 restricciones y 110 variables, resumiendo de estos resultados se tiene el siguiente cuadro. m

n

iteraciones

tiempo (seg.)

Valor Objetivo

Simplex

90

110

164

8.03

1597.36

MPIPD

90

110

15

0.437

1597.40

Aquí se observa también que el algoritmo Simplex necesita aproximadamente 11 veces más de iteraciones, que las iteraciones de MPIPD, el cual solo realizó 15 iteraciones. En cuanto al tiempo de ejecución, la diferencia es de 7.59 segundos

Experimento 7: Aquí se generó un problema de programación lineal con 220 restricciones y 320 variables y se obtuvo los siguientes resultados, m

n

iteraciones

tiempo (seg.)

Valor Objetivo

Simplex

220

320

762

692.28

-185790

MPIPD

220

320

21

8.656

-185791.78

En este problema se muestra la e ficacia del MPIPD respecto al Simplex, en cuanto al tiempo de ejecución y las respectivas iteraciones

74


Conclusiones En este trabajo se construyó, detalladamente, el método más e ficiente en la actualidad que nos permite resolver el problema de programación lineal a gran dimensión: El Método de puntos interiores Primal-Dual para programación lineal. Podemos concluir entonces que:

1. El método de puntos interiores Primal-Dual es relativamente nuevo, y constituye una alternativa al método Simplex. 2. Los métodos de puntos interiores son mecanismos que abordan el problema de un modo diferente: crean una sucesión de puntos por el interior de la región de factibilidad, la cual converge hacia un óptimo. 3. El método de puntos interiores Primal-Dual está basado en tres estrategias numéricas, El método de Newton, Lagrange y Penalización interna. 4. El método de puntos interiores Primal-Dual es de fácil implementación computacional. En este caso usamos el Matlab. 5. En la parte experimental, cuando el Método es sometido a problemas de grandes dimensiones, resultó tener un buen desempeño.

75

76

CAPÍTULO CAPÍTULO 3. EXPERIMENTOS EXPERIMENTOS COMPUT COMPUTACIONALES ACIONALES

.

Apéndice A

Simplex: Una Implementación en Matlab del algoritmo simplex function function [z,x,iter [z,x,iter] ] = SIMPLEX(A SIMPLEX(A,b,c, ,b,c,colB colB) ) format format compac compact; t; format format bank; bank; [m,n]=size(A);rmin=inf;iter=0; parar=’no’; while while parar==’n parar==’no’ o’ iter=iter+1; Bv=[ Bv=[ ];cBv= ];cBv=[ [ ];Nv=[ ];Nv=[ ];cNv= ];cNv=[ [ ];Rv=[ ];Rv=[ ]; for j=1:m j=1:m Bv=[Bv,A(:,colB(j))]; cBv=[cBv;c(colB(j))]; end for j=1:n j=1:n repite=max(colB==j); if repite==0 repite==0 Nv=[Nv,A(:,j)]; cNv=[cNv;c(j)]; Rv=[Rv;j]; end end if abs(det(B abs(det(Bv)) v))<1.e-9 error(’Mat error(’Matriz riz B probablem probablemente ente sea singular singular’); ’); elseif elseif inv(Bv)*b inv(Bv)*b>-1.e-9 B=Bv;N=Nv;cB=cBv;cN=cNv;R=Rv; else error( error(’B ’B no define define SBF’); SBF’); end

77 Zj_Cj=cB’*inv(B)*N-cN’; b_=inv(B)*b; z0=cB’*inv(B)*b; [valor,k]=max(Zj_Cj); yk=inv(B)*N(:,k); x=zeros(n,1); x(colB)=b_; if valor valor<=1.e-9 parar=’si parar=’si’; ’; z=z0; z=z0; elseif elseif yk<=1.e-9 parar=’si parar=’si’; ’; z=-inf; z=-inf; else rmin=inf; for i=1:m i=1:m if yk(i) yk(i)>1.e-9 1.e-9 & b_(i)/ b_(i)/yk( yk(i) i)
end end

78

CAPÍTULO CAPÍTULO 3. EXPERIMENTOS EXPERIMENTOS COMPUT COMPUTACIONALES ACIONALES

problema 1.6 Minimizar

x1 + 2x 2x2 x21 − 6x1 + x22 − 8x2 + 20 = 0

function F = func (x,y)

function J = Jacobian(x,y)

F = [1+(2*x(1)-6)*y

J=[2*y 0 2*x(1)-6

2+(2*x(2)-8)*y

0 2*y 2*x(2)-8

x(1)^2 x(1)^2+x( +x(2)^ 2)^2-6 2-6*x( *x(1)1)-8*x 8*x(2) (2)+20 +20]; ];

2*x(1) 2*x(1)-6 -6 2*x(2) 2*x(2)-8 -8 0];

Los valores iniciales >> x = [1 2 ]t, y = 3,

Al resolver en el computador se observa después de 6 iteraciones : x(6) =

"

2.0000 2.0000

#

79

Apéndice B Notación En este apéndice se lista la notación principal usada en este trabajo. A x(k) x∗ xi X, Z k n L SBF SB F MPIPD PL KK T

Matriz m × n Diferentes puntos obtenidos por el algoritmo Denotarán puntos óptimos o punto solución Componente Componente de vectores vectores Matrices Se reserva reserva para indicar el número de iteraciones Número de variables variables Tamaño de un ejemplar Solución básica factible Método de puntos interiores Primal-Dual Programación Programación Lineal Condiciones de optimalidad de Karush-Kuhn-T Karush-Kuhn-Tuck ucker er

80


LIBRO Optimizaciòn Lineal

Recommend Documents