Optimization

Consigna 1

Nunca quemes un núcleo de tu computador

  • Usa monitor de temperatura
  • R-Cran tiene herramientas muy avanzadas para esto
# library(multicore)
# load the package
# x1=1:5;x2=5:10
# create 2 objects
# p1=parallel(factorial(x1)) # run in parallel
# p2=parallel(factorial(x2)) # run in parallel
# collect(list(p1,p2))
# collect results

# $‘8995‘ [1]  1  2  6 24 120
# $‘8996‘ [1]  120  720  5040  40320  362880 3628800

Primeros pasos en R

R se puede usar como una calculadora avanzada en el escritorio de windows linux o android.

3+5
[1] 8
a<-3 
b<-5 
a+b
[1] 8
sin(pi/2)
[1] 1
sqrt(9)
[1] 3

Filosofía Matricial

  • Estructura de Datos
  • Variables Indexadas
  • Definición de Funciones

Estructura de Datos

  • Variables como a <-3 pueden ser valores numéricos, pero R los interpreta como matrices. a en este caso es una matriz de 1x1 con el valor guardado 3 .

  • Se emplea la flecha y no el signo igual para asignar a efectos de no confundir con la operación de comparación que sí usa =.

  • En caso de if condicional se usa == para preguntar por igualdad

Tipos de Objetos Matriz

  • vectores (matriz de una fila o columna)
  • listas (matriz con caracteres alfabéticos)
  • data-frame (matriz con elementos numéricos y alfabéticos) cargados en memoria

Vectores

car_name <-c("Honda","BMW","Ferrari")
car_color =c("Black","Blue","Red")
car_cc =c(2000,3400,4000)

car_name
[1] "Honda"   "BMW"     "Ferrari"
class(car_name)
[1] "character"
class(car_cc)
[1] "numeric"

Listas

Son más versátiles que los vectores, pero pueden tener elementos numéricos y alfabéticos que provienen de vectores disimiles. Piense en ellos como vector de vectores.

cars <-list(name =c("Honda","BMW","Ferrari"),
color =c("Black","Blue","Red"),
cc =c(2000,3400,4000))

cars
$name
[1] "Honda"   "BMW"     "Ferrari"

$color
[1] "Black" "Blue"  "Red"  

$cc
[1] 2000 3400 4000

Usar el comando class

Matrices

Son estructuras como las listas, pero limitadas sólo a datos numéricos. Pueden venir de aparatos que emiten constantemente valores (caja registradora) o fuentes de web (Servicio Meteorológico), incluso pueden venir del teclado usando la función scan()

mdat <-matrix(c(1,2,3, 11,12,13), nrow =2, ncol =3, byrow =TRUE,
dimnames =list(c("Fila1", "Fila2"),
c("Col1", "Col2", "Client3")))
mdat
      Col1 Col2 Client3
Fila1    1    2       3
Fila2   11   12      13

DataFrame

Los dataframe extienden las matrices con la capacidad adicional de contener tipos de datos heterogéneos. En un dataframe, puede almacenar variables de caracteres, numéricas y de factores en diferentes columnas del mismo df. En casi todas las tareas de análisis de datos, con filas y columnas de datos, df se presenta como una opción natural para almacenar los datos y maneja mucho mejor la memoria y recursos del hardware

L3 <-LETTERS[1:3]
fac <-sample(L3, 10, replace =TRUE)
df <-data.frame(x =1, y =1:10, fac = fac)

DataFrame 2

class(df$x)
[1] "numeric"
class(df$y)
[1] "integer"
class(df$fac)
[1] "factor"

Nueva forma de invocar vectores = variable $ columna

Variables Indexadas o Subseting

R tiene uno de los operadores de subconjunto más avanzados, potentes y rápidos en comparación con cualquier otro lenguaje de programación. Es poderoso hasta el punto de que, salvo en algunos casos que discutiremos en la siguiente sección, no existe una construcción de bucle for o while sino que se trabaja todo maticialmente, aunque R explícitamente proporciona estas palabras reservadas si es necesario.

Subseting 2

Aunque es muy poderoso, sintácticamente podría llegar a ser una pesadilla o podría aparecer un error grave si no se presta atención cuidadosa al colocar el número requerido de paréntesis, corchetes y comas. Los operadores [, [[y $ se usan para subconjuntos, dependiendo de qué estructura de datos contiene los datos. También es posible combinar el subconjunto con la asignación para realizar una función realmente complicada con muy pocas líneas de código.

Subseting en Vectores

Para los vectores, la subconjunto podría hacerse haciendo referencia al índice respectivo de los elementos almacenados en un vector. Por ejemplo, car name [c (1,2)] devolverá los elementos almacenados en el índice 1 y 2, y car name [-2] devolverá todos los elementos excepto el segundo. También es posible utilizar operadores binarios para indicar al vector que recupere o no un elemento.

carname <-c ("Honda", "BMW", "Ferrari")
#Selecciona el 1 ° y el 2 ° elemento del vector
carname [c (1,2)]
[1] "Honda" "BMW"  

Subseting Vectores 2

#Seleciona todo, salve el 2do elemento 
carname[-2]
[1] "Honda"   "Ferrari"

Otra forma de seleccionar el 2do elemento

carname[c(FALSE,TRUE,FALSE)]
[1] "BMW"

Subseting en Listas

En listas es similar al subseting en un vector; sin embargo, dado que una lista es una colección de muchos vectores, se debe usar corchetes dobles para recuperar un elemento de la lista. Por ejemplo, cars \( [2] \) recupera el segundo vector completo de la lista y \( cars [[c (2,1)]] \) recupera el primer elemento del segundo vector.

cars <-list(name =c("Honda","BMW","Ferrari"),
color =c("Black","Blue","Red"),
cc =c(2000,3400,4000))

cars[2]$color
[1] "Black" "Blue"  "Red"  
cars[[c(2,1)]]
[1] "Black"

Subseting en Matrices

Las matrices tienen un subconjunto similar como vectores. Sin embargo, en lugar de especificar un índice para recuperar los datos, necesitamos dos índices aquí: uno que indique la fila y el otro para la columna.

Por ejemplo, mdat [1: 2,] recupera todas las columnas de las primeras dos filas, mientras que \( mdat [1: 2, "C.1"] \) recupera las primeras dos filas y C.1 column.mn.

Subseting en Matrices 2

mdat <-matrix(c(1,2,3, 11,12,13), nrow =2, ncol =3,
byrow =TRUE, dimnames =list(c("F.1", "F.2"),
c("C.1", "C.2", "C.3")))

Subseting en Matrices 3

mdat[1:2,] #Select first two rows and all columns
    C.1 C.2 C.3
F.1   1   2   3
F.2  11  12  13
mean(mdat[1,])
[1] 2
mdat[2,3]
[1] 13

Sumario o información sumaria

summary(cars)
      Length Class  Mode     
name  3      -none- character
color 3      -none- character
cc    3      -none- numeric  

Cuidado ! , existe un dataset llamado cars. Si definimos una variable con el nombre de un dataset podemos tener problemas de confusión.

Ejercicio

  • Generar 500 número aleatorios con distribución normal, media = 50 , desvio estandard = 12

  • Ver media, calcular varianza de la muestra, moda, mediana y cuantiles

Cargar Paquetes

Como se ha señalado existen más de veinte mil paquetes oficiales a la fecha y cada paquete suele venir con un set de datos de prueba que sirven para ver como trabaja el modelo.

  • ejecute el comando data() para ver que tiene disponible

  • ejecute data(package = “datasets”) para ver el contenido del datos del paquete datasets

  • tienes que instalar antes el paquete con:

  • install.packages(“dataset”)

  • y cargarlo con

  • library(dataset)

Ver paquetes instalados

  • Siquieres saber que paquetes tienes instalados ejecuta:
  • Elimina head para ver todos
head(installed.packages())
         Package    LibPath                                         
abind    "abind"    "/home/rpalma/R/x86_64-pc-linux-gnu-library/3.4"
acepack  "acepack"  "/home/rpalma/R/x86_64-pc-linux-gnu-library/3.4"
adabag   "adabag"   "/home/rpalma/R/x86_64-pc-linux-gnu-library/3.4"
alabama  "alabama"  "/home/rpalma/R/x86_64-pc-linux-gnu-library/3.4"
alluvial "alluvial" "/home/rpalma/R/x86_64-pc-linux-gnu-library/3.4"
askpass  "askpass"  "/home/rpalma/R/x86_64-pc-linux-gnu-library/3.4"
         Version    Priority Depends                            
abind    "1.4-5"    NA       "R (>= 1.5.0)"                     
acepack  "1.4.1"    NA       NA                                 
adabag   "4.2"      NA       "rpart, caret, foreach, doParallel"
alabama  "2015.3-1" NA       "R (>= 2.10.1), numDeriv"          
alluvial "0.1-2"    NA       NA                                 
askpass  "1.1"      NA       NA                                 
         Imports          LinkingTo
abind    "methods, utils" NA       
acepack  NA               NA       
adabag   NA               NA       
alabama  NA               NA       
alluvial NA               NA       
askpass  "sys (>= 2.1)"   NA       
         Suggests                                                Enhances
abind    NA                                                      NA      
acepack  "testthat"                                              NA      
adabag   "mlbench"                                               NA      
alabama  NA                                                      NA      
alluvial "devtools, testthat, reshape2, knitr, rmarkdown, dplyr" NA      
askpass  "testthat"                                              NA      
         License              License_is_FOSS License_restricts_use
abind    "LGPL (>= 2)"        NA              NA                   
acepack  "MIT + file LICENSE" NA              NA                   
adabag   "GPL (>= 2)"         NA              NA                   
alabama  "GPL (>= 2)"         NA              NA                   
alluvial "MIT + file LICENSE" NA              NA                   
askpass  "MIT + file LICENSE" NA              NA                   
         OS_type MD5sum NeedsCompilation Built  
abind    NA      NA     "no"             "3.4.4"
acepack  NA      NA     "yes"            "3.4.4"
adabag   NA      NA     "no"             "3.4.4"
alabama  NA      NA     "no"             "3.4.4"
alluvial NA      NA     "no"             "3.4.4"
askpass  NA      NA     "yes"            "3.4.4"

Ver descripción de paquetes

Si esto te da error debes instalar el paquete alluvial.

packageDescription("alluvial")
Package: alluvial
Type: Package
Title: Alluvial Diagrams
Version: 0.1-2
Date: 2016-09-09
Authors@R: c( person("Michal", "Bojanowski", role=c("aut", "cre"),
       email="michal2992@gmail.com"), person("Robin", "Edwards",
       role="aut", email="robin.edwards@ucl.ac.uk") )
Description: Creating alluvial diagrams (also known as parallel
       sets plots) for multivariate and time series-like data.
URL: https://github.com/mbojan/alluvial
BugReports: https://github.com/mbojan/alluvial/issues
Suggests: devtools, testthat, reshape2, knitr, rmarkdown, dplyr
License: MIT + file LICENSE
LazyLoad: yes
LazyData: yes
VignetteBuilder: knitr
RoxygenNote: 5.0.1
NeedsCompilation: no
Packaged: 2016-09-09 09:58:05 UTC; mbojan
Author: Michal Bojanowski [aut, cre], Robin Edwards [aut]
Maintainer: Michal Bojanowski <michal2992@gmail.com>
Repository: CRAN
Date/Publication: 2016-09-09 13:08:51
Built: R 3.4.4; ; 2019-04-05 02:31:51 UTC; unix

-- File: /home/rpalma/R/x86_64-pc-linux-gnu-library/3.4/alluvial/Meta/package.rds 

Ayuda de uso de paquete

Para obtener ayuda de un paquete ejecutar

help(package = "alluvial")

Ver ventana de ayuda abajo a la derecha …

Muestreo

Generar un muestreo entre 10 individuos

Primero generamos los individuos

x <- sample(1:10)
x
 [1]  9  1  7  4 10  3  6  8  5  2

¿Qué valor tienen las muestras 4 a 8?

match(c(4,8),x)
[1] 4 8

Elegir individuos entre el cuatro y el octavo

match(c(4,8),x)
[1] 4 8

Muestreo con repetición

Generamos un muestreo de 30 muestras con 7 individuos y repetición

x <- sample(1:7,30,replace=TRUE)
x
 [1] 2 4 6 5 2 7 3 3 2 2 4 2 5 6 5 2 3 7 1 1 6 5 4 1 5 3 1 6 1 2

¿En que casos aparecieron los individuos 2 o 4?

which(x %in% c(2,4))
 [1]  1  2  5  9 10 11 12 16 23 30

Importar DATOS

SECCIÓN II

Captura de datos desde el teclado

Generaremos un vector llamado productividad obtenido de un estudio realizado sobre 25 empresas.

productividad <- scan()

Tipear estos números. Al final de cada número presione entre e ingrese el que sigue. En el último caso tipear dos veces enter.

Se puede ejecutar este comando desde el teclado para cargar un archivo local

library(car)
library(readr)
partners <- read.table("Socios.csv",header=TRUE,sep=";")
partners2 <- read-xlsx ("D:/Datos/Socios2.xls")

Usar getwd() y setwd()

Leer datos de la web

library(car)
# Buscamos la tabla http://ceal.fing.uncu.edu.ar/r-cran/Socios.csv
partners <-read.table("http://themys.sid.uncu.edu.ar/rpalma/R-cran/Socios.csv",header=TRUE,sep=";")

names(partners)
[1] "TIR"    "VAN"    "Invent" "Asset"  "Type"  
head(partners)
  TIR VAN Invent Asset Type
1 5.1 3.5    1.4   0.2   Co
2 4.9 3.0    1.4   0.2   Co
3 4.7 3.2    1.3   0.2   Co
4 4.6 3.1    1.5   0.2   Co
5 5.0 3.6    1.4   0.2   Co
6 5.4 3.9    1.7   0.4   Co

¿Cuántas columnas tiene la tabla partners?

length(partners)
[1] 5

¿Cuantas filas tiene?

Podemos saberlo examinando la longitud de cualquier columna

length(partners$VAN)
[1] 150

Muestreo Estadístico

Para ello debemos pasar la tabla a memoria (attach)

attach(partners)
str(partners)
'data.frame':   150 obs. of  5 variables:
 $ TIR   : num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
 $ VAN   : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
 $ Invent: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
 $ Asset : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
 $ Type  : Factor w/ 3 levels "Co","Da","Pi": 1 1 1 1 1 1 1 1 1 1 ...
  • ¿Qué hace attach(partners)?
  • ¿ Qué hace str(partners)?

Separaración población muestra

Separaremos a la población en una muestra que tiene el 30% y reservaremos el resto (70%) de la población para otros cálculos o verificación de modelos

muestra_2 <- sample(2, nrow(partners), replace=TRUE, prob=c(0.7, 0.3))
  Ejecutar View muestra 2
muestra_2
  [1] 2 2 1 1 1 2 2 1 1 1 1 1 1 2 2 1 2 1 1 1 1 1 1 2 2 1 2 1 2 1 2 1 1 1 1
 [36] 1 1 1 1 1 1 2 1 2 2 1 1 1 2 1 2 1 2 2 1 1 1 1 2 2 1 2 1 1 1 2 1 1 1 1
 [71] 2 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 2 1 1 1 1 1
[106] 1 2 2 1 1 1 2 1 2 2 1 1 2 1 1 1 1 2 2 1 2 2 1 1 2 1 1 1 1 2 1 1 1 1 1
[141] 2 2 1 1 2 1 1 2 2 1

Mostrar los datos de los casos usados en la muestra

head( partners [muestra_2==1,])
   TIR VAN Invent Asset Type
3  4.7 3.2    1.3   0.2   Co
4  4.6 3.1    1.5   0.2   Co
5  5.0 3.6    1.4   0.2   Co
8  5.0 3.4    1.5   0.2   Co
9  4.4 2.9    1.4   0.2   Co
10 4.9 3.1    1.5   0.1   Co

Hoja de Ruta

SUPER EJERCICIO

*Tenemos un centro de manufactura automatizada que tiene siete clientes que le envían trabajos a sus tres CNC.

*Según el BOM cada operación ha sido identificada y se conoce que secuencia se debe seguir en cada máquina.

Hoja de Ruta 2

  • Se ha identificado la secuencia de máquinas a seguir por los pedidos de cada cliente.

  • Todos los pedidos pasan por las tres máquina sin importar en que orden lo hagan.

  • Se conocen los tiempos de mecanizado y las demoras que se generarán en el transporte desde uno a otro más el tiempo de espera estimado en la cola de entrada del centro de mecanizado.

Generar las hojas de ruta

  • El tiempo de mecanizado sigue una ley beta (pert) de promedio 0.6 hs con apertura 2 y excentricidad 0.1 + 0.3

  • En tanto el transit time más la espera se han aproximado con una weibull de promedio 0.5 min con spread de 1

  • Graficar histogramas de ambos tiempos para el período de trabajo (130 operaciones) que en promedio se realizan en un día.

Preguntas

  • ¿Cual es el tiempo máximo de manufactura? wtime
  • Cual es el tiempo máximo de traslado / espera? qtiem histograma
  • Realizar los mismos histogramas con 20 cajas
  • Realizar el gráfico de densidad
  • Determinar las horas a facturar a cada cliente
  • Hacer un análisis sumario del tiempo que requiere cada cliente
  • ¿Que probabilidad tengo de que el tiempo de manufactura sea menor o igual a 2 horas?

Generación simulada hoja de ruta

x <- sample(1:7,130,replace=TRUE)
HojaRuta <- list (cliente = LETTERS[x] , 
                  op_1 = round(runif(130,1,4)),
                  op_2 = round(runif(130,1,4)),
                  op_3= round(runif(130,1,4)) ,
                  wtime= 0.3+rbeta(130,0.6,2,0.1),
                  queuet= (1+rweibull(130,0.5,1))/60 )
summary(HojaRuta$wtime)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.3004  0.3667  0.5036  0.5509  0.7381  1.2679 

Tiempo máximo de manufactura wtime

Convertiremos la lista en matriz

 matrizHR <- matrix(unlist(HojaRuta), ncol = 6, byrow = FALSE)
  View matrizHR
library(nnet)
which.is.max(matrizHR) 
[1] 109
matrizHR[74, ]
[1] "F"                 "2"                 "3"                
[4] "3"                 "0.550351039095638" "0.118797098914528"
  ver el comando class matrizHR

Densidad de la función

Trnasit and queue time

density( as.numeric(matrizHR[ ,5]))

Call:
    density.default(x = as.numeric(matrizHR[, 5]))

Data: as.numeric(matrizHR[, 5]) (130 obs.); Bandwidth 'bw' = 0.0749

       x                 y           
 Min.   :0.07566   Min.   :0.000475  
 1st Qu.:0.42990   1st Qu.:0.120869  
 Median :0.78414   Median :0.454183  
 Mean   :0.78414   Mean   :0.704941  
 3rd Qu.:1.13838   3rd Qu.:1.107017  
 Max.   :1.49262   Max.   :2.127150  

Gráficas de Tiempos

par(mfrow=c(2,1))
hist( as.numeric(matrizHR[ ,5]),
      breaks = 25, 
      main = "Histograma aplanado",
      sub = "Tiempo de manufactura wtime",
      ylab= "Repeticiones",
      xlab="Tiempo H:mm" )

plot(density( as.numeric(matrizHR[ ,5])))

plot of chunk unnamed-chunk-36

Trucos

Para resolver estos puntos se recomienda inspirarse en este ejemplo

 x=seq(0,35,length=350)
 x1=seq(0,13,length=130)
 x2=seq(13,35,length=(350-130))
 y=dnorm(x,mean=14.7,sd=5.7)
 y1=dnorm(x1,mean=14.7,sd=5.7)
 y2=dnorm(x2,mean=14.7,sd=5.7)

 plot(x,y,type="l", lwd=2, col="blue") 
 polygon(c(0,x1,13),c(0,y1,0),col=rgb(1, 0.9, 0,0.5) )

plot of chunk unnamed-chunk-37

Distribución Normal

  • R tiene una colección importante de distribuciones de probabilidades.
  • Sin lugar a dudas la Distribución Normal es la que más se usa (acepta)
  • Es posible obtener ayuda sobre el comando que trabaja con esta distribución con
?pnorm

Mire en la ventana de ayuda lo que aparece y explique que hacen los comando dnorm, pnorm, qnorm y rnorm.

Ejemplos de Normales

n <- 100
x <- rnorm(n)
par(mfrow=c(1,2), las=1)
for(i in 1:8) {
  y <- i*x + rnorm(n)
  plot(x, y, main=i)
}

plot of chunk unnamed-chunk-39plot of chunk unnamed-chunk-39plot of chunk unnamed-chunk-39plot of chunk unnamed-chunk-39

Uso de tablas

n <- 100
w <- rnorm(n)
y <- 2*w + rnorm(n)
out <- lm(y ~ w)
library(knitr)
kable(summary(out)$coef, digits=2)
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.14 0.1 -1.46 0.15
w 1.97 0.1 19.84 0.00
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.1 0.10 1.02 0.31
w 2.1 0.09 23.40 0.00

Otras Distribuciones

Distrib. Comandos
Beta pbeta qbeta dbeta rbeta
Binomial pbinom qbinom dbinom rbinom
Cauchy pcauchy qcauchy dcauchy rcauchy
Chi-Square pchisq qchisq dchisq rchisq
Exponential pexp qexp dexp rexp
F pf qf df rf
Gamma pgamma qgamma dgamma rgamma
Geometric pgeom qgeom dgeom rgeom

Otras distribuciones (2)

Distrib. Comandos
Hypergeometr phyper qhyper dhyper rhyper
Logistic plogis qlogis dlogis rlogis
Log Normal plnorm qlnorm dlnorm rlnorm
Neg.Binomial pnbinom qnbinom dnbinom rnbinom
Normal pnorm qnorm dnorm rnorm
Poisson ppois qpois dpois rpois
Student t pt qt dt rt
Studentized Range ptukey qtukey dtukey rtukey

Otras distribuciones (3)

Distrib. Comandos
Uniform punif qunif dunif runif
Weibull pweibull qweibull dweibull rweibull
Wilcox Rank pwilcox qwilcox dwilcox rwilcox
WilcoxSigned psignrank qsignrank dsignrank rsignrank

Timeline

Las variables vinculadas al tiempo (fecha y hora) siempre son un problema a la hora de operar con ellas. Ver ejemplo de excel.

La biblioteca timevis provee muchas instacias para solucionar estos problemas

library(timevis)

data1 <- data.frame(
  id      = 1:4,
  content = c("Inicio", "Reunión_1",
              "Rev.Papers", "Modelado"),
  start   = c("2018-01-10", "2018-01-11",
              "2018-01-20", "2018-02-14 15:00:00"),
  end     = c(NA, NA, "2018-02-04", NA)
)

Timeview 2

Ejecutar desde la línea de comando

timevis(data1)

Ver resultado en ventana de ayuda

Prolblemas de Optimización

SECCIÓN III

Clasificación

  • Optimización Continua
  • Optimización Convexa
  • Optimización Convexa Lineal'' o “Programación Lineal''
  • Programación no lineal
  • Optimización Combinatoria
  • Optimización Mixta
  • Optimización binaria mixta

  • Clasificación base a la naturaleza probabilística

Clasificación de métodos de resolución

  • Resolución mediante cálculo
  • Resolución mediante técnicas de búsquedas
  • Resolución mediante técnicas de convergencia de soluciones

Resolución mediante cálculo

Los métodos de resolución por cálculo apelan al cálculo de derivadas para determinar para qué valores del dominio la función presenta un máximo o un mínimo. Son métodos de optimización muy robustos, pero que requieren mucho esfuerzo de cómputo y que tanto la función objetivo como las restricciones presenten determinadas condiciones (por ejemplo, que estén definidas, que sean continuas, etc.)

Resolución mediante técnicas de búsquedas:

Dentro de este apartado, podemos encontrar un gran abanico de técnicas, desde el viejo método de prueba y error hasta las modernas técnicas de programación matemática. En forma genérica, consisten en el siguiente algoritmo: Seleccionar una solución inicial y hacerla la solución actual:

Algoritmo de búsqueda

1) Hacer Mejor Solución = Solución Actual

2) Buscar n soluciones cercanas a la Solución Actual

3) Para cada una de las n soluciones cercanas hacer

  a. Si el valor de la función objetivo de la solución a verificar es mayor (o menor) al valor generado por la solución actual, hacer Mejor Solución = Solución Evaluada

  b. Si Mejor Solución = Solución Actual, Finalizar del procedimiento, en caso contrario Hacer Solución Actual = Mejor Solución y volver a 2)

Introducción a la Programación Lineal

Los problemas de Programación Lineal son problemas de optimización convexa en el cual tanto la función objetivo como las restricciones son funciones lineales. Si bien su definición parece un tanto restrictiva, en realidad muchos de los problemas estudiados en el mundo real entran dentro de esta categoría, por ejemplo:

  • Asignación de rutas de transporte
  • Balance de líneas de producción
  • Asignación de procesos a equipos de procesamiento
  • Definición de mix de productos
  • Definición de mezcla de componentes para productos
  • Nivelación de recursos

Caracterización

Un problema de programación lineal se caracteriza porque su conjunto de soluciones factibles forma un politopo (o sea, la generalización de un polígono para n dimensiones). Por su propia naturaleza, al estar los gradientes definidos en una solo dirección, los puntos extremos se encuentran en la frontera del politopo, más específicamente en algunos de sus vértices. Esto permitió el desarrollo del muy eficiente Método Simplex para su resolución, que se aprovecha de esta particularidad para buscar soluciones solo en dichos puntos.

Cómo resolver problemas de PL con R-Cran}

Dentro de los paquetes “estándares'' que se distribuyen con R, tenemos el paquete “boot”. Entre otras, este paquete tiene una función denominada simplex que se encarga de aplicar el algoritmo de dicho nombre. Para usarla, primero carguemos el paquete:

library(boot)

La función “simplex” se encarga de resolver problemas de programación lineal. Permite únicamente el uso de solo variables continuas y fuerza las soluciones no negativas. Veamos un ejemplo de uso:

Ejemplo:

problema1

problema1

Resolución

Este problema se puede resolver en R-cran invocando la biblioteca anterior como sigue:

coefFuncionZ <- c(1,2,3)
matrizRestriccionesMenorOIgual <-
matrix(c(0,2,3,5,0,3),2,3,byrow=TRUE)
vectorTermIndep <- c(100,100)
mySimplex <- simplex(coefFuncionZ, A1 =
matrizRestriccionesMenorOIgual, b1= vectorTermIndep,
maxi=TRUE )
mySimplex

Linear Programming Results

Call : simplex(a = coefFuncionZ, A1 = matrizRestriccionesMenorOIgual, 
    b1 = vectorTermIndep, maxi = TRUE)

Maximization Problem with Objective Function Coefficients
x1 x2 x3 
 1  2  3 


Optimal solution has the following values
x1 x2 x3 
20 50  0 
The optimal value of the objective  function is 120.

Deventajas

simplex() es una función compleja para utilizar debido al gran número de parámetros que requiere. Veamos su ayuda:

?simplex

Simplex como objeto

Como objeto hecho y derecho, podemos acceder a los atributos del objeto generado por simplex:

class(mySimplex)
[1] "simplex"
  • mySimplex$soln # Soluciones
  • mySimplex$soln[2] # Solución de la segunda variable
  • mySimplex$value # Valor de la función objetivo
  • mySimplex$slack # Holgura de cada restricción menor igual
  • mySimplex$surplus # Excedente de cada restricción mayor igual

Atributos

Este objeto tiene, además, otros atributos relativos al cálculo:

mySimplex$a # Coeficientes de Z con variables básicas 
[1] 0.0 0.0 0.6 1.0 0.2
#             con coeficientes igual a cero
mySimplex$basic # Los índices de las variables básicas 
[1] 2 1
#                 en la solución encontrada
mySimplex$A # Matriz de restricciones expresada en 
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    1 -1.5 -0.5  0.0
[2,]    1    0 -0.6  0.0 -0.2
#             términos de variables no básicas

Otros paquetes de PL: lpSolve

Vamos a trabajar con otro paquete para resolver problemas de optimización lineal. Como no es estándar, tenemos que descargarlo desde CRAN:

  install.packages("lpSolve")

Procedimiento

En RStudio podría haber ido al menú “Tools”, seleccionar el ítem “Install Packages” y escribir “lpSolve” en el cuadro de texto:

Una vez instalado, tenemos que cargarlo al workspace 38 actual con la función “library”:

library(lpSolve)

Característica

lpSolve encapsula el acceso al programa también llamado lpSolve .

El mismo tiene licencias LGPL y está disponible en forma de bibliotecas.

El paquete lpSolve implementa algunas de las funcionalidades del programa lpSolve.

¿Se acuerdan cómo funcionaba la función “simplex” del paquete “boot”?. Teníamos tres matrices de coeficientes de restricciones y tres vectores de términos independientes, uno por cada tipo de restricción ( \( \leq \), \( = \), \( \geq \)). En este paquete las cosas funcionan un poco más simples.

Bueno, resolvamos un problema sencillo con lpSolve.

Z función objetivo a optimizar

Primero, definamos los coeficientes de la función lineal objetivo Z:

problem.coef.obj <- c(1, 9, 3)

Como vemos, es una función de tres variables. Supongamos que tenemos dos restricciones, una del tipo \( \geq \) y otra del tipo \( \leq \):

Restricciones

problem.coef.restr <- matrix (c(1, 2, 3, 3, 2, 2), nrow=2,
  byrow=TRUE)

        problem.coef.restr
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    3    2    2

Restricciones

En una única matriz defino mis dos restricciones.

Como el parámetro “byrow” vale “TRUE”, tenemos la primera restricción de coeficientes

    (1, 2, 3) 

y la segunda de coeficientes

    (3,-2,2).

Signo de las restricciones

Y ahora defino los signos de mis restricciones en un vector del tipo “character”:

problem.dir.restr <- c(">=", "<=")
problem.dir.restr
[1] ">=" "<="

Lado derecho de la ecuación

Seguimos con los términos independientes:

problem.term.restr <- c(1, 15)
problem.term.restr
[1]  1 15

Maximicemos la función objetivo:

lp ("max", problem.coef.obj, problem.coef.restr,
problem.dir.restr, problem.term.restr)
Success: the objective function is 67.5 

Ventajas

En cuanto a sintaxis, mucho mejor que usar la función simplex. Además, el algoritmo simplex que implementa la función lp() es mucho más veloz. También a diferencia de solver en Excel tiene capacidad de administrar problemas mucho más grandes, incluso en máquinas de pequeño porte sin provocar desbordamientos de memoria “la temida pantalla azul de windows”

Como ya sabemos, todo en R es un objeto:

mylp <- lp ("max", problem.coef.obj, problem.coef.restr,
problem.dir.restr, problem.term.restr)
mylp
Success: the objective function is 67.5 

Y como objeto, podemos acceder a sus miembros. Primero, las soluciones del problema:

mylp$solution
[1] 0.0 7.5 0.0

También el valor de Z alcanzado:

mylp$objval
[1] 67.5

Sensibilidad

¿Queremos hacer análisis de sensibilidad?, bueno, hagamos valer TRUE el parámetro “compute.sens” y volvamos a resolver:

mylp <- lp ("max", problem.coef.obj, problem.coef.restr,
problem.dir.restr, problem.term.restr, compute.sens=TRUE)

Sensibilidad

¿Queremos hacer análisis de sensibilidad?, bueno, hagamos valer TRUE el parámetro “compute.sens”:

mylp$sens.coef.from
[1] -1e+30  3e+00 -1e+30

Y accedamos a los atributos “sens.coef.from” y “sens.coef.from” para ver la sensibilidad de los coeficientes de la función objetivo:

mylp$sens.coef.from
[1] -1e+30  3e+00 -1e+30

“compute.sens=TRUE” nos permite hacer análisis de dualidad también:

Sensibilidad

mylp$duals
[1]   0.0   4.5 -12.5   0.0  -6.0

Donde sí tengo “n” variables y “m” restricciones, los primeros m valores corresponden a los valores duales de las restricciones y los siguientes n valores a los duales de las actuales variables de decisión. Puedo ver la sensibilidad de los duales también:

Dualidad

mylp$duals.from
[1] -1.0e+30  1.0e+00 -1.0e+30 -1.0e+30 -1.4e+01
mylp$duals.to
[1] 1.0e+30 1.0e+30 5.0e+00 1.0e+30 7.5e+00

No me devuelve las holguras, pero las puedo calcular a mano. Primero devuelvo la matriz de las restricciones:

Slack Constraints

mylp$constraints
              [,1] [,2]
                 1    3
                 2    2
                 3    2
const.dir.num    2    1
const.rhs        1   15

La matriz se devuelve transpuesta (es la forma en que lpSolve trabaja internamente). La primera columna es la primera restricción, la segunda es la segunda restricción. La última fila representa a los términos independientes, la anteúltima son referencias al sentido de las restricciones (2 es \( \leq \), 1 es \( \geq \) y 3 \( = \)) y los primeros valores los coeficientes de las restricciones. Entonces, podemos hacer:

Constraints

sum(mylp$constraints[1:3,1] * mylp$solution) -
mylp$constraints[5,1]
const.rhs 
       14 

Constraints

sum(mylp$constraints[1:3,2] * mylp$solution) -
mylp$constraints[5,2]
    const.rhs 
-1.776357e-15 

Problema de Almacenes y Clientes

El paquete lpSolve tiene algunas funcionalidades adicionales. Viene con dos funciones que implementan los algoritmos de asignación y transporte. Veamos el caso de transporte. Definamos la matriz de costo:

matrizCosto <- matrix (10000, 8, 5)
matrizCosto
      [,1]  [,2]  [,3]  [,4]  [,5]
[1,] 10000 10000 10000 10000 10000
[2,] 10000 10000 10000 10000 10000
[3,] 10000 10000 10000 10000 10000
[4,] 10000 10000 10000 10000 10000
[5,] 10000 10000 10000 10000 10000
[6,] 10000 10000 10000 10000 10000
[7,] 10000 10000 10000 10000 10000
[8,] 10000 10000 10000 10000 10000

Este problema podría enunciarse del siguiente modo: Tengo 8 almacenes (filas) que tienen que abastecer a 5 Clientes. Genero una matriz con los costos de transportar desde cada almacén a cada cliente. En principio he tomado un valor muy alto para ese costo, con el objeto de poner luego costos reales y que al minimizar el costo total jamás se tome ese camino de costo alto.

En esta matriz cada elemento [i,j] nos dice cual es costo unitario de transporte desde i a j. La modificaremos un poco para hacerla más real:

Alternativa (cargar costos a mano)

Podría usar scan() para cargar los datos a mano, pero quiero que vean como puedo hacerlo por código. Incluso puedo armar un bucle for / while para hacer que el costo crezca según una ley dada

matrizCosto[4,1]<- matrizCosto[-4,5]<- 0
matrizCosto
      [,1]  [,2]  [,3]  [,4]  [,5]
[1,] 10000 10000 10000 10000     0
[2,] 10000 10000 10000 10000     0
[3,] 10000 10000 10000 10000     0
[4,]     0 10000 10000 10000 10000
[5,] 10000 10000 10000 10000     0
[6,] 10000 10000 10000 10000     0
[7,] 10000 10000 10000 10000     0
[8,] 10000 10000 10000 10000     0
matrizCosto[1,2]<- matrizCosto[2,3]<- matrizCosto[3,4] <- 7
matrizCosto[1,3]<- matrizCosto[2,4]<- 7.7
matrizCosto
      [,1]  [,2]    [,3]    [,4]  [,5]
[1,] 10000     7     7.7 10000.0     0
[2,] 10000 10000     7.0     7.7     0
[3,] 10000 10000 10000.0     7.0     0
[4,]     0 10000 10000.0 10000.0 10000
[5,] 10000 10000 10000.0 10000.0     0
[6,] 10000 10000 10000.0 10000.0     0
[7,] 10000 10000 10000.0 10000.0     0
[8,] 10000 10000 10000.0 10000.0     0
matrizCosto[5,1]<- matrizCosto[7,3]<- 8
matrizCosto
      [,1]  [,2]    [,3]    [,4]  [,5]
[1,] 10000     7     7.7 10000.0     0
[2,] 10000 10000     7.0     7.7     0
[3,] 10000 10000 10000.0     7.0     0
[4,]     0 10000 10000.0 10000.0 10000
[5,]     8 10000 10000.0 10000.0     0
[6,] 10000 10000 10000.0 10000.0     0
[7,] 10000 10000     8.0 10000.0     0
[8,] 10000 10000 10000.0 10000.0     0

Modifiquemos aún más

matrizCosto[1,4] <- 8.4
matrizCosto[6,2] <- 9
matrizCosto[8,4] <- 10
matrizCosto[4,2:4] <- c(.7, 1.4, 2.1)
matrizCosto
      [,1]  [,2]    [,3]    [,4]  [,5]
[1,] 10000 7e+00     7.7     8.4     0
[2,] 10000 1e+04     7.0     7.7     0
[3,] 10000 1e+04 10000.0     7.0     0
[4,]     0 7e-01     1.4     2.1 10000
[5,]     8 1e+04 10000.0 10000.0     0
[6,] 10000 9e+00 10000.0 10000.0     0
[7,] 10000 1e+04     8.0 10000.0     0
[8,] 10000 1e+04 10000.0    10.0     0

Restricciones

Armemos el vector de restricciones. Estos vectores indican que en ninguno caso la suma de las filas debe ser superior a la cantidad que tiene cada almacén. Por otro lado la suma de las columnas tienen que ser exactamente lo que cada cliente pide.

row.dir <- rep ("<", 8)
row.term <- c(200, 300, 350, 200, 100, 50, 100, 150)
col.dir <- rep (">", 5)
col.term <- c(250, 100, 400, 500, 200)
  • La primer línea dice que todos los 8 almacenes tiene una restricción del tipo “menor que”. Fijarse como con el comando rep repito la dirección de la restricción

  • La segunda línea señala la cantidad que cada almacén tiene

  • La tercer es la restricción que dice que no podemos enviar más que lo que el cliente pidió

  • Finalmente tenemos un vector que dice cuanto han pedido los clientes.

Para resolver el problema Corramos hora el comando lp.transport():

Resolver problema de transporte

myTransport <- lp.transport(matrizCosto, "min", row.dir,
row.term, col.dir, col.term)
myTransport
Success: the objective function is 7790 

Soluciones

Veamos cuanta carga se ha asignado para obtener este costo mínimo

myTransport$solution
     [,1] [,2] [,3] [,4] [,5]
[1,]    0  100  100    0    0
[2,]    0    0  200  100    0
[3,]    0    0    0  350    0
[4,]  200    0    0    0    0
[5,]   50    0    0    0   50
[6,]    0    0    0    0   50
[7,]    0    0  100    0    0
[8,]    0    0    0   50  100

Muy “inteligentemente”, R nos devuelve una matriz de asignación de orígenes a destinos, en vez de una lista de valores de variables de decisión.

El problema de asignación

Definir un problema de asignación es más fácil, ya que al ser todas las restricciones \( \leq \) con variables binarias, lpSolve me automatiza la carga, por lo cual solo me ocupo de la matriz de costo:

Este problema de asignación es típico en casos de problemas de pallets sobre piso en posiciones flotantes, también se usa para diseño de layouts, casos simples de ruteo, etc. En la biografiara existe un método llamado “Método Húngaro” pero la heurística de lpsolve es mucho más poderosa.

Matriz de Costos

Generaremos la matriz de costos de asignación, que debe interpretarse así. Cada elemento A(i,j) representa el costo de poner el elemento que se señala en la fila i en la posición de piso j.

matrizCostoAsignacion <- matrix (c(2, 7, 7, 2, 7, 7, 3, 2,
7, 2, 8, 10, 1, 9, 8, 2), 4, 4)
matrizCostoAsignacion
     [,1] [,2] [,3] [,4]
[1,]    2    7    7    1
[2,]    7    7    2    9
[3,]    7    3    8    8
[4,]    2    2   10    2

Caracteristicas

  • En el método Húngaro, de debía proporcionar una solución inicial.

  • Esta matriz solución llamada H era una matriz diagonal unitaria.

  • Se establecían restricciones del tipo:

  • La suma de toda las filas debe ser igual a 1

  • La suma de todas las columnas debe ser igual a 1

    Esto se puede interpretar como “un lugar para cada cosa” y “cada cosa en su lugar”

Ventajas de lpsolve

Sin embargo no es necesario usar esto en lpsolve.

Se invoca el problema como de asignación y se resuelve.

myAssign <- lp.assign(matrizCostoAsignacion, "min")
myAssign
Success: the objective function is 8 
  • Este es el costo mínimo.

  • Luego el plan de uso de las posiciones o asignaciones de productos al almacén serán las que se exponen a continuación.

  • La cantidad de problemas que se pueden resolver por este método es realmente impresionante.

  • Este es uno de los terrenos en los que R-Cran corre con ventajas frente a otros productos.

  • En la Universidad Nacional de Cuyo existía un problema serio de aisgnación de aulas por franja horaria, con vaientntes que complicaban más aún el problema, tales como la necesidad para algunos profesores śolo en algunas clases al año de alta conectividad en el aula o necesidad de un proyectos de mayor resolución de pantalla que los que se utilizan a diario.

  • Este método ha demostrado ser muy eficaz para resolverlo.

Solución encontrada

myAssign$solution
     [,1] [,2] [,3] [,4]
[1,]    0    0    0    1
[2,]    0    0    1    0
[3,]    0    1    0    0
[4,]    1    0    0    0

Otras alternativas de Optimización con R

En los ejemplos anteriores hemos visto la importancia de modelar (antes que resolver) problemas

Hemos visto también los distintos tipoos de probleas de optimización y las limitaciones que podemos encontrar al trabajarlos con hojas de calculo. En este apartado nos sentraremos un una serie de bilbiotecas (library) que pueden ayudarnos cuando la potencia de las hojas de cálculo no nos dejan tener certeza de los resultados obtenidos.

Otros paquetes de Programación Lineal: lpSolveAPI

lpSolve no es la única implementación de “lpSolve” para R, tenemos también lpSolveAPI:

#install.packages("lpSolveAPI")
library(lpSolveAPI)

Los modelos con lpSolveAPI se construyen de otra forma. Para empezar, no son objetos R hechos y derechos, así que es muy parecido su uso a como se usaría desde C.

Con lpSolveAPI comienzo así:

Uso de lpsolveapi

lprec <- make.lp(0, 4)
lprec
Model name: 
            C1    C2    C3    C4    
Minimize     0     0     0     0    
Kind       Std   Std   Std   Std    
Type      Real  Real  Real  Real    
Upper      Inf   Inf   Inf   Inf    
Lower        0     0     0     0    

Con esto creo un objeto “modelo lineal” con 4 variables de decisión. Este es un objeto vacio, pero que tiene la estructura para que podamos cargar nuestro modelo (o problema a resolver) Definamos ahora los coeficientes:

Coeficientes

set.objfn(lprec, c(1, 3, 6.24, 0.1))
lprec
Model name: 
            C1    C2    C3    C4    
Minimize     1     3  6.24   0.1    
Kind       Std   Std   Std   Std    
Type      Real  Real  Real  Real    
Upper      Inf   Inf   Inf   Inf    
Lower        0     0     0     0    

Restricciones

Agrego las restricciones:

add.constraint(lprec, c(0, 78.26, 0, 2.9), ">=", 92.3)
add.constraint(lprec, c(0.24, 0, 11.31, 0), "<=", 14.8)
add.constraint(lprec, c(12.68, 0, 0.08, 0.9), ">=", 4)
lprec
Model name: 
             C1     C2     C3     C4          
Minimize      1      3   6.24    0.1          
R1            0  78.26      0    2.9  >=  92.3
R2         0.24      0  11.31      0  <=  14.8
R3        12.68      0   0.08    0.9  >=     4
Kind        Std    Std    Std    Std          
Type       Real   Real   Real   Real          
Upper       Inf    Inf    Inf    Inf          
Lower         0      0      0      0          

En caso de restricciones cuya única función sea acotar una única variable, las agrego separadamente como “Bounds”: Esto suele ser casi siempre una limitación del tipo presupuestaria (no puedo superar cierto gasto, cierno número de camiones, etc)

Limites

Esto no se puede hacer con los otros métodos ya vistos

set.bounds(lprec, lower = c(28.6, 18), columns = c(1, 4))
set.bounds(lprec, upper = 48.98, columns = 4)
lprec
Model name: 
             C1     C2     C3     C4          
Minimize      1      3   6.24    0.1          
R1            0  78.26      0    2.9  >=  92.3
R2         0.24      0  11.31      0  <=  14.8
R3        12.68      0   0.08    0.9  >=     4
Kind        Std    Std    Std    Std          
Type       Real   Real   Real   Real          
Upper       Inf    Inf    Inf  48.98          
Lower      28.6      0      0     18          

Normalmente el modelo como variables puede ser claro para quien modela el problema, pero puede ser confuso para quien toma decisiones. Tengo ponerle nombre:

Nombre de las variables de decisión

rowNames <- c("Restricción_01", "Restricción_02", "Restricción_03")
colNames <- c("Var_X1", "Var_X2", "Var_X3", "Var-X4")
dimnames(lprec) <- list(rowNames, colNames)
lprec
Model name: 
                Var_X1  Var_X2  Var_X3  Var-X4          
Minimize             1       3    6.24     0.1          
Restricción_01       0   78.26       0     2.9  >=  92.3
Restricción_02    0.24       0   11.31       0  <=  14.8
Restricción_03   12.68       0    0.08     0.9  >=     4
Kind               Std     Std     Std     Std          
Type              Real    Real    Real    Real          
Upper              Inf     Inf     Inf   48.98          
Lower             28.6       0       0      18          

Resolución

Listo, ahora podemos pasar a resolverlo: Sorprende la facilidad con la que lo hacemos, solamente llamamos a un comando invocando el objeto creado.

solve(lprec)
[1] 0

Si deseo ver el valor optimal de Z (valor alcanzado por la función a optimizar) invoco el comando siguiente.

get.objective(lprec)
[1] 31.78276

Valores de las variables

Obtener las soluciones o valores de la variables en el óptimo:

get.variables(lprec)
[1] 28.60000  0.00000  0.00000 31.82759

Valor de las restricciones asumidos o proporcionados

get.constraints(lprec)
[1]  92.3000   6.8640 391.2928

Programación Lieal Entera

Resulta comprensible el por qué es más conveniente esta bilioteca que las otras (nativas) que trae el paquete R.

Trnasformar un problema de LP en ILP (Integer Linear Programin) resulta un juego de niños.

Supongamos que tenemos el mismo problema anterior y queremos que la variable 4 sea etntera. Esta variable puede represetar el número de containers que tendré en inventario, la cantidad de choferes a asignar a cada ruta u otra variable que me interese manter entera.

Problema entero

set.type(lprec, 4, "integer")
lprec
Model name: 
                Var_X1  Var_X2  Var_X3  Var-X4          
Minimize             1       3    6.24     0.1          
Restricción_01       0   78.26       0     2.9  >=  92.3
Restricción_02    0.24       0   11.31       0  <=  14.8
Restricción_03   12.68       0    0.08     0.9  >=     4
Kind               Std     Std     Std     Std          
Type              Real    Real    Real     Int          
Upper              Inf     Inf     Inf   48.98          
Lower             28.6       0       0      18          

Soluciones enteras

Veamos ahora las souciones obtenidas y coparemos con el caso anterior.

solve(lprec)
[1] 0
get.objective(lprec)
[1] 31.792
get.variables(lprec)
[1] 28.60000000  0.03066701  0.00000000 31.00000000
get.constraints(lprec)
[1]  92.300   6.864 390.548

Otro Problema de Asignación limitado

Este tipo de problema se llama bounded, quiere decir que además de la matrz de unos y ceros, tenemos variables que solamente pueden tomar valores entre ciertos límites que pueden ser superiores o inferiores o ambos.

Podemos intuir que siguiendo la síntaxis anterior sería posible declarar una variable binaria (en lugar de entera). Esto, tal como se vio en los primeros ejercicios, no lleva a tener la base para armar un problema de asignación.

Bouded Problem

Usando el ejemplo anterior veamos como declaro una variable binaria.

Si la segunda variable debiera ser binaria:

set.type(lprec, 2, "binary")
lprec
Model name: 
                Var_X1  Var_X2  Var_X3  Var-X4          
Minimize             1       3    6.24     0.1          
Restricción_01       0   78.26       0     2.9  >=  92.3
Restricción_02    0.24       0   11.31       0  <=  14.8
Restricción_03   12.68       0    0.08     0.9  >=     4
Kind               Std     Std     Std     Std          
Type              Real     Int    Real     Int          
Upper              Inf       1     Inf   48.98          
Lower             28.6       0       0      18          

Soluciones

Luego repasemos las soluciones Z, valores de X y restricciones:

solve(lprec)
[1] 0
get.objective(lprec)
[1] 31.8
get.variables(lprec)
[1] 28.6  0.0  0.0 32.0
get.constraints(lprec)
[1]  92.800   6.864 391.448

Algoritmo Genético

Problema tipo PnakSack

El algoritmo genético (GA) es una técnica de optimización basada en la búsqueda basada en los principios de genética y selección natural. Con frecuencia se utiliza para encontrar soluciones óptimas o casi óptimas a problemas difíciles que de otra manera podrían llevar una vida para resolverse. Se utiliza con frecuencia para resolver problemas de optimización, cálculo variacional, simulación y análisis de escenarios, así como en la investigación, y en el aprendizaje automático (Deep Learning).

GA

GA se basa en la generción de una población, los individuos en esta población (a menudo representdos como cromosomas) tienen un estado dado. Una vez que se genera la población, el estado de estos individuos se evalúa y se clasifica en su valor o desempeño. Se debe crear una función de desempeña para que a partir de estado cromosómico se obtenga una cifra de preformance. Las mejores candidatos son tomadas y cruzadas -en la esperanza de generar una mejor descendencia, para formar la nueva población. En algunos casos se conservan los mejores individuos de la población con el fin de garantizar “buenos candidatos” en la nueva generación (esto se llama elitismo, entendido como procentaje de padres que prevalecen frente a hijos).

GA

Para explicar el ejemplo se utilizarar una versión simplificada de la problema mochila (pnacksack).

Te proponen pasar un mes en el desierto. Partirás llevando una mochila y tu supervivencia depende de lo que cargues, sin embargo, el peso máximo que puede transportar es de 20 kilogramos. Tienes una serie de artículos de supervivencia disponibles, cada uno con su propio número de “puntos de supervivencia”. El objetivo es maximizar el número de puntos de supervivencia.

GA

En la siguiente tabla se muestran los elementos que puede elegir.

Ga

Resolución con R

library(genalg)
library(ggplot2)

Estructura de datos

dataset <- data.frame(item = c("pocketknife", "beans", "potatoes", "unions", 
    "sleeping bag", "rope", "compass"), survivalpoints = c(10, 20, 15, 2, 30, 
    10, 30), weight = c(1, 5, 10, 1, 7, 5, 1))
weightlimit <- 20

Formato

Veamos como ha quedado conformado el set de datos

dataset
          item survivalpoints weight
1  pocketknife             10      1
2        beans             20      5
3     potatoes             15     10
4       unions              2      1
5 sleeping bag             30      7
6         rope             10      5
7      compass             30      1

Configuración de la Función de Evaluación

Antes de crear el modelo tenemos que configurar una función de evaluación. La función de evaluación evaluará a los diferentes individuos (cromosomas) de la población sobre el valor de su configuración genética.

Por ejemplo, un individuo puede tener la siguiente configuración de genes: 1001100. Esto representa algo así como que llevo los itmes señalados con 1 en el cromosoma, tomando como referencia la columna 2 del dataset. Podemos elegir a esta como una solución inicial, independientemente de cumplir o no con la restricción de 20 kg.

Evaluación ...

Cada número de esta cadena binaria representa si desea o no llevar un elemento con usted. Un valor de 1 se refiere a poner el artículo específico en la mochila mientras que un 0 se refiere a dejar el elemento en casa. Dado el ejemplo de configuración de genes, tomaríamos los siguientes elementos;

Podemos comprobar a qué cantidad de puntos de supervivencia de esta configuración. Esto le dará un valor a la configuración génica de un cromosoma dado. Esto es exactamente lo que hace la función de evaluación.

Codoficación Cromosomica

chromosome = c(1, 0, 0, 1, 1, 0, 0)
dataset[chromosome == 1, ]
          item survivalpoints weight
1  pocketknife             10      1
4       unions              2      1
5 sleeping bag             30      7

Operatoria

Mutiplico cada uno y cadad cero del cromosoma por el survivalpoint

cat(chromosome %*% dataset$survivalpoints)
42

Crearemos una función de desempeño para no tener que repetir en pada paso los cálculo.

Funcicón para nuestro problema de mochila.

En GAs, tenemos una uestra o una población inicial de posibles soluciones para el problema dado, no nos interesa si cumplen o no con la restricción. Estas soluciones se someten a la recombinación y mutación (como en la genética natural), produciendo nuevos niños, y el proceso se repite a lo largo de varias generaciones. A cada individuo (o solución candidata) se le asigna un valor de aptitud (basado en su valor de función objetivo) y a los individuos más aptos se les da una mayor probabilidad de aparearse y producir más individuos hijos. Esto está en consonancia con la teoría darwiniana de la “supervivencia del más apto'' (que en realidad es una falacia).

El más apto

El algoritmo genalg intenta optimizar siempre hacia el valor mínimo. Por lo tanto, el valor se calcula por la función que construimos, y se multiplica por-1 si se quiere maximizar. Una configuración que lleva a exceder la restricción de peso devuelve un valor de 0 (también se puede dar un valor más alto o distinto). Esto último se conoce como tunnear el algoritmo.

evalFunc <- function(x) {
    current_solution_survivalpoints <- x %*% dataset$survivalpoints
    current_solution_weight <- x %*% dataset$weight

    if (current_solution_weight > weightlimit) 
        return(0) else return(-current_solution_survivalpoints)
}

nota: observe que en R el símbolo asterisco se utiliza para señalar multiplicación. En este caso si lo utilizo tendría el producto vectorial de dos matrices y obtendría error si no coinciden la cantidad de filas de la primerca con el número de columnas de la segunda. Para poder obtener la multiplicación de cada valor de una matriz (o como en este caso un vector) por su homólogo en la otra entidad se ha adoptado el signo asterísco pero rodeado por el de porcentaje. Esto hace el producti de las componentes una a una.

Corrida del modelo

A continuación, elegimos el número de iteraciones, diseñamos y ejecutamos el modelo.

iter = 100
GAmodel <- rbga.bin(size = 7, popSize = 200, iters = iter, 
      mutationChance = 0.01,  elitism = T, evalFunc = evalFunc)
      plot(GAmodel)

plot of chunk unnamed-chunk-93

Interpretación gráfica

Generamos un histograma de la población elitista que mejor resuleve el problema

plot(GAmodel,type = "hist")

plot of chunk unnamed-chunk-94

Vemos que es una población de comportamento discreto.

Encontrar la mejor solución

Qué significa la mejor solución?

solution = c(1, 1,0, 1, 1, 1, 1)
dataset[solution == 1, ] 
          item survivalpoints weight
1  pocketknife             10      1
2        beans             20      5
4       unions              2      1
5 sleeping bag             30      7
6         rope             10      5
7      compass             30      1

Solución

La mejor solución se encuentra en 1111101. Esto nos lleva a tomar los siguientes artículos con nosotros en nuestro viaje a la naturaleza.

cat(paste(solution %*% dataset$survivalpoints, "/", sum(dataset$survivalpoints)))
102 / 117

Ejercicios para la evaluación

SECCIÓN 4

  • Este ejercicio debe ser resuletos con R-Cran primero con con lpsolve_api y luego con Algoritmo Genético

  • Se deberá medir la performance de tiempo que casa solución demanda

    En caso de no querer resolverlo se puede presentar un caso de estudio redactado en formato paper, que no supere las 10 páginas, en el que se modela y resuelve usando R-Studio con alguno de los métodos vistos.

Plan de producción

    Plantear este ejercicio en R

Consideremos la siguiente situación: Una pequeña empresa vende dos productos, denominados Producto 1 y Producto 2. Cada tonelada del Producto 1 consume 30 horas de trabajo, y cada tonelada del producto 2 consume 20 horas de trabajo. El negocio tiene un máximo de 2.700 horas de trabajo durante el período considerado. En cuanto a las horas de máquina, cada tonelada de Productos 1 y 2 consume 5 y 10 horas de máquina respectivamente. Hay 850 horas de máquina disponibles.

Plan de producción

Cada tonelada del Producto 1 produce 20 Me de ganancia, mientras que el Producto 2 produce 60 Me por cada tonelada vendida. Por razones técnicas, la empresa debe producir un mínimo de 95 toneladas en total entre ambos productos. Necesitamos saber cuántas toneladas de Producto 1 y 2 deben ser producidas para maximizar el beneficio total.

Esta situación es apta para ser modelada como un modelo PL. Primero, necesitamos definir las variables de decisión. En este caso tenemos:

  • P 1 número de toneladas producidas y vendidas del Producto 1
  • P 2 número de toneladas producidas y vendidas del producto 2

Coeficientes

Los coeficientes de costo de estas variables son 20 y 60, respectivamente. Por lo tanto, la función objetivo se define multiplicando cada variable por su Coeficiente de coste correspondiente.

Las limitaciones de este LP son:

  • Una limitación de Horas de Trabajo que haga que la cantidad total de horas de trabajo utilizadas en el Producto 1 y Producto 2, que sea igual a 30P 1 + 20P 2, sea inferior o igual a 2.700 horas.

  • Una restricción similar Horas de máquina que hacen que el total de horas de máquina 5P 1 + 10P 2 sea menor o igual que 850.

  • Una restricción PM que hace que el total de unidades producidas y vendidas P1 + P2 son mayores o iguales que 95.

Condición de no negatividad

Juntando todo esto, y considerando que las variables de decisión son no negativas, el LP que maximiza el beneficio es:

Problema

Resolver problemas difíciles

En Ciencias de la computación, hay un gran conjunto de problemas, que son NP-Hard. Lo que esto significa esencialmente es que, incluso los sistemas de computación más poderosos toman un tiempo muy largo (¡ incluso años!) para resolver ese problema. En tal escenario, GAs demuestra ser una herramienta eficiente para proporcionar soluciones utilizables casi óptimas en un corto período de tiempo.

aaa borrar

For more details on authoring R presentations please visit https://support.rstudio.com/hc/en-us/articles/200486468.

  • Bullet 1
  • Bullet 2
  • Bullet 3

Slide With Code

summary(cars)
      Length Class  Mode     
name  3      -none- character
color 3      -none- character
cc    3      -none- numeric  

Slide With Plot

Error in xy.coords(x, y, xlabel, ylabel, log) : 
  'x' is a list, but does not have components 'x' and 'y'