Práctica 6: PROGRAMACIÓN LINEAL CON R

1 Introducción

Un problema de programación lineal (PPL) es una situación compleja en la que el responsable debe tomar decisiones del tipo:

• elegir cantidades de producción para maximizar beneficios,
• decidir la ruta de un viaje para minimizar costes,
• determinar un reparto de tareas para minimizar el tiempo de finalización,
• etc.

de modo que se cumpla con ciertas limitaciones, como:

• no superar la materia prima, el personal y el almacenaje disponibles,
• pasar necesariamente por una serie de puntos comprometidos,
• no exceder la capacidad de los encargados,
• etc.

y de manera que, tanto la función objetivo (a maximizar o minimizar), como las limitaciones se pueden escribir como ecuaciones lineales en las variables de decisión.

Una forma estándar en la que todo PPL puede escribirse es:

\[ \begin{array}{rrrrrrrrc} \text{max } z(\vec{x}) = &c_1 x_1 & + & c_2 x_2 & + & \cdots & + & c_n x_n & &\\ \text{s.a} & & & & & & & & & \\ &a_{11} x_1 & + &a_{12} x_2 & + &\cdots & + &a_{1n} x_n & = &b_1 \\ &a_{21} x_1 & + & a_{22} x_2 & + & \cdots & + & a_{2n} x_n & = & b_2 \\ & \vdots & \vdots &\vdots & \vdots &\ddots & \vdots &\vdots & \vdots &\vdots \\ &a_{m1} x_1 & +& a_{m2} x_2 & + &\cdots & +& a_{mn} x_n & =& b_m \\ & x_1, & & x_2, & & \ldots, & & x_n & \geq & 0 \end{array} \]

donde todas las restricciones quedan como ecuaciones lineales.

El conjunto de valores que cumplen las restricciones se llama región factible y se denota por \(F\).

La resolución “manual” de un PPL pasa por:

• Modelizarlo:
  • identificando las variables de decisión,
  • escribiendo la función objetivo,
  • escribiendo las restricciones.
• Pasarlo a forma estándar:
  • Transformando variables no restringidas como resta de dos no negativas,
  • Transformando inecuaciones en ecuaciones:
    • Sumando variables de holgura,
    • restando variables de exceso,
    • cambiando signo si el término independiente \(b_j\) sale negativo
• Aplicar el algoritmo SIMPLEX:
  • Añadir variables artificiales (si hace falta) para la tabla simplex inicial.
  • Aplicar el algoritmo hasta su detención.

2 El algoritmo SIMPLEX

Algunos detalles de la región factible \(F\) son:

• Sus elementos son candidatos a solución del problema, ya que cumplen las restricciones, y se llaman soluciones factibles.
• \(F\) es cerrado, y por eso:
  • Si \(F\) es acotado, entonces es seguro que la función objetivo alcanza su óptimo en un punto de \(F\).
  • Si \(F\) no es acotado, es posible que el óptimo sea \(+\infty\) (en ese caso sería un problema no acotado).
• \(F\) es convexo, y por eso:
  • Tiene unos puntos especiales, llamados puntos extremos.
  • Si \(F\) es acotado, los puntos no extremos de \(F\) siempre se pueden ubicar “entre unos cuantos de los puntos extremos”.
  • Si se alcanza el óptimo en un punto de \(F\), es seguro que hay un punto extremo de \(F\) que también alcanza el óptimo.

El algoritmo Simplex (George B. Dantzig, 1947) se basa en esto último. Mediante sencillas operaciones:

1. Sale de un punto extremo de \(F\) (tabla simplex inicial).
2. Es capaz de saltar de un punto extremo de \(F\) a otro adyacente, siempre que se mejore la función objetivo (actualizar tabla simplex).
3. Se detiene si saltar a cualquier otro punto extemo adyacente no mejora la función objetivo.

Existen muchos programas para la aplicación del algoritmo simplex, y algunos permiten tratar el PPL sin estar en forma estándar, lo cual quita trabajo al usuario.

Este algoritmo no está exento de inestabilidad numérica, ya que en la actualización de la tabla simplex, al saltar de un extremo a otro, se calculan divisiones con un elemento pivote como denominador. Si se trata de un valor cercano a 0 la exactitud del algoritmo peligra.

3 Función simplex() del package boot

Sirve para problemas de tamaño no muy grande (leer Nota en la ayuda help(simplex)). Instalar y cargar package previamente.

• Función: simplex(a, A1, b1, A2, b2, A3, b3, maxi, n.iter, eps)
• Argumentos:
  • a: coeficientes de la función objetivo
  • A1: matriz de restricciones de tipo \(\leq\) (sólo la parte de coeficientes)
  • b1: vector de términos independientes de las restricciones de tipo \(\leq\) (correspondientes a la matriz A1)
  • A2: idem con restricciones de tipo \(\geq\)
  • b2: idem correspondientes a la matriz A2
  • A3: idem con restricciones de tipo \(=\)
  • b3: idem correspondientes a la matriz A3
  • maxi: Poner a TRUE si se desea maximizar la función objetivo
  • n.iter: número máximo de iteraciones en cada fase (usa método de dos fases)
  • eps: tolerancia de coma flotante
• Devuelve: un objeto de la clase simplex, que es una lista con componentes (entre otras):
  • $soln: solución óptima, si se alcanza
  • $solved: código de terminación:
    • 1 indica solución óptima encontrada
    • 0 indica que se ha alcanzado el número máximo de iteraciones sin acabar
    • -1 indica que no existe solución factible
  • $value: valor de la función objetivo en la solución $soln
  • $val.aux: valor de las variables artificiales en el óptimo
  • etc.

Ejemplo 1: Resuelve el PPL expresado como: \[ \begin{array}{rl} \text{max} &x_1 + 5x_2 - x_3 \\ \text{s.a} & \\ & x_1 + x_2 + x_3 \geq 100 \\ & 5x_1 - x_2 + 8x_3 \leq 500 \\ & - x_1 + x_2 + 2x_3 = 0 \\ & 2x_1 + x_2 + 12x_3 \leq 800 \\ & x_1, x_2, x_3 \geq 0 \end{array} \]

Ejemplo 2: Resuelve el PPL expresado como: \[ \begin{array}{rl} \text{max} &15 x_1 + 10x_2\\ \text{s.a} & \\ & x_2 \leq 50 \\ & 1.5 x_1 - x_2 \geq 20 \\ & x_1, x_2 \geq 0 \end{array} \]

Ejemplo 3: Resuelve el PPL expresado como: \[ \begin{array}{rl} \text{min} &50 x_1 -20 x_2\\ \text{s.a} &\\ & - 1.5 x_2 \leq 15 \\ & x_1 + x_2 \geq 10 \\ & x_1, x_2 \geq 0 \end{array} \]

4 Función lp() del package lpSolve

Se trara de otra función para resolver PPLs, pero más completa pues admite variables enteras (programación entera) y variables binarias (programación 0-1). Para usarla hay que instalar y cargar el package previamente.

• Función: lp (direction, objective.in, const.mat, const.dir, const.rhs, int.vec, binary.vec)
• Argumentos:
  • direction: por defecto vale 'min'. Cambiar a 'max' si se quiere maximizar la función objetivo
  • objective.in: vector con coeficientes de la función objetivo
  • const.mat: matriz de coeficientes de las restricciones
  • const.dir: vector de cadenas de texto indicando la dirección de cada restricción ("<=", "=," o ">=")
  • const.rhs: vector de términos independientes de las restricciones
  • int.vec: Vector indicando qué variables son enteras
  • binary.vec: Vector indicando qué variables son 0-1
• Devuelve: un objeto de la clase lp, que es una lista con componentes (entre otras):
  • $objval: valor de la función objetivo en el óptimo
  • $solution: solución óptima
  • $status: indicador numérico:
    • 0: éxito
    • 2: no hay solución factible

Ejemplo 4: Resuelve los ejemplos 1, 2 y 3 por medio de la función lp.

4 Ejercicios evaluables

set.seed(20200507) # PON TU FECHA DE NACIMIENTO AÑOMESDIA
misproblemas = sort(sample(x=1:6, size=2))

Ejercicio 1: De la lista que figura más abajo, borra todos los ejercicios excepto los números 2 y 5, y resuélvelos por duplicado (usando ambas funciones de R para cada PPL), comentando el resultado alcanzado según:

• Solución óptima: justificar, según salida de R, escribir la solución óptima (variables) y valor de la función objetivo.
• Problema no acotado: justificar, según salida de R.
• Problema infactible: justificar, según salida de R.

LISTA DE PROBLEMAS (hacer sólo los 2 indicados en el enunciado):

1. \(\left\{\begin{array}{rl} \text{max} &x_1 -2x_2 -3x_3 -x_4\\ \text{s.a} & \\ & x_1 - x_2 - 2x_3 - x_4 \leq 4 \\ & 2x_1 + x_3 - 4x_4 \leq 2 \\ & -2x_1 + x_2 + x_4 \leq 1 \\ & x_1, x_2, x_3, x_4 \geq 0 \end{array} \right.\)

# AQUÍ EL CÓDIGO

Y AQUÍ LOS COMENTARIOS SOBRE LA SOLUCIÓN QUE TE HA DADO, O NO PUNTUARÁ

---

2. \(\left\{\begin{array}{rl} \text{min} &3y_1 -y_2 + 2y_3\\ \text{s.a} & \\ & 2y_1 - y_2 + y_3 \geq -1 \\ & y_1 + 2y_3 \geq 2 \\ & -7y_1 + 4y_2 - 6y_3 \geq 1 \\ & y_1, y_2, y_3 \geq 0 \end{array} \right.\)

# AQUÍ EL CÓDIGO

Y AQUÍ LOS COMENTARIOS SOBRE LA SOLUCIÓN QUE TE HA DADO, O NO PUNTUARÁ

---

3. \(\left\{\begin{array}{rl} \text{max} &-x_1 -x_2 + 2x_3\\ \text{s.a} & \\ & -3x_1 + 3x_2 + x_3 \leq 3 \\ & 2x_1 - x_2 - 2x_3 \leq 1 \\ & -x_1 + x_3 \leq 1 \\ & x_1, x_2, x_3 \geq 0 \end{array} \right.\)

# AQUÍ EL CÓDIGO

Y AQUÍ LOS COMENTARIOS SOBRE LA SOLUCIÓN QUE TE HA DADO, O NO PUNTUARÁ

---

4. \(\left\{\begin{array}{rl} \text{min} &5y_1 -2y_2 -y_3\\ \text{s.a} & \\ & -2y_1 + 3y_3 \geq -1 \\ & 2y_1 - y_2 + y_3 \geq 1 \\ & 3y_1 + 2y_2 - y_3 \geq 0 \\ & y_1, y_2, y_3 \geq 0 \end{array} \right.\)

# AQUÍ EL CÓDIGO

Y AQUÍ LOS COMENTARIOS SOBRE LA SOLUCIÓN QUE TE HA DADO, O NO PUNTUARÁ

---

5. \(\left\{\begin{array}{rl} \text{min} &-2y_2 +y_3 \\ \text{s.a} & \\ & -y_1 - 2y_2 \geq -3 \\ & 4y_1 + y_2 + 7y_3 \geq -1 \\ & 2y_1 - 3y_2 + y_3 \geq -5 \\ & y_1, y_2, y_3 \geq 0 \end{array} \right.\)

# AQUÍ EL CÓDIGO

Y AQUÍ LOS COMENTARIOS SOBRE LA SOLUCIÓN QUE TE HA DADO, O NO PUNTUARÁ

---

6. \(\left\{\begin{array}{rl} \text{max} &3x_1 + 4x_2 + 5x_3 \\ \text{s.a} & \\ & x_1 + 2x_2 + 2x_3 \leq 1 \\ & -3x_1 + x_3 \leq -1 \\ & -2x_1 - x_2 \leq -1 \\ & x_1, x_2, x_3 \geq 0 \end{array} \right.\)

# AQUÍ EL CÓDIGO

Y AQUÍ LOS COMENTARIOS SOBRE LA SOLUCIÓN QUE TE HA DADO, O NO PUNTUARÁ

Ejercicio 2: Inicia la programación del algoritmo simplex en lenguaje R, en una versión lo más simplificada posible. Instrucciones:

• Solo valdrá para:
  • Problemas de programación lineal (no entera ni 0-1)
  • En forma estándar
    • \(\max z = \vec{c} \vec{\mathbf{x}}\)
    • s.a. \(A \vec{\mathbf{x}} = \vec{b}\) (donde \(\vec{b}\) será de coeficientes no negativos)
  • Donde la matriz tecnológica \(A\) tendrá incluida la submatriz identidad (en las últimas columnas)
• La función debe devolver la solución óptima (si la hay), o una frase que diga “problema no acotado” o “problema infactible”.
• Se facilita el principio de la programación:

simplex0 = function(c, A, b) {
  # comprobaciones necesarias
  stopifnot( is.numeric(c), is.numeric(b), is.matrix(A))
  n = length(c); m = length(b)
  stopifnot(b>=0, n>=m, dim(A)[1]==m, dim(A)[2]==n)
  ident = diag(x=m)
  stopifnot( identical(A[,(m+1):n], ident) )
  base = (m+1):n # columnas base inicial
  x.base = b # primera sol básica
  c.base = c[base] # coeficientes de las var. básicas
  Y = A # la matriz de coeficientes y_kj
  zj.minus.cj = apply(X=c.base * Y, MARGIN=2, FUN='sum') - c # los zj menos cj
  z.base = sum(c.base*x.base)
  # a partir de aquí te toca a ti, suerte
}

Comprueba con un ejemplo completo de clase.