domingo, 26 de mayo de 2013

Using R to visualize geo optimization algorithms
 
May 26, 2013
By JohannesSch

http://www.r-bloggers.com/

(This article was first published on progRamming for the fun of it , and kindly contributed to R-bloggers)

site optimization is the process of finding an optimal location for a plant or a warehouse to minimize transportation costs and duration. A simple model only consists of one good and no restrictions regarding transportation capacities or delivery time. The optimizing algorithms are often hard to understand. Fortunately, R is a great tool to make them more comprehensible.

The basic math
A basic contineous optimization problem is to deliver goods to $n$ customers. Every customer has a individual demand $a_j$ and is approached directly from our warehouse. The single routes have a length of $d_j$. The objective function to be minimized is:
$$min(\sum^{n}_{j=1} a_j \cdot d_j)$$
How can on calculate distances? Typical distance measures for flat surfaces most of you are probably familiar with are Taxicab distance and Euclidean distance. We will calculate spherical distances to get a better solution for large routes: $$\\d_{a,b} = arccos[sin(x_a)\cdot sin(x_b)+cos(x_a)\cdot cos(x_b)\cdot cos(y_a-y_b)] \cdot R\\$$ R denotes the Earth radius: ~ 6,370 km. Keep in mind:
    ...To use geographical coordinates in the form of decimal degrees
    ... To ransform them to radian units (by multiplying pi/360)

How to optimize
We will use an iterative way for optimizing:
    set an initial warehouse location: $x_{start}$ and $y_{start}$
    for every iteration the current warehouse location $(x_w, y_w)$ is calculated by: $$x_w=\frac{\sum_{j=1}^{n}{a_j \cdot \overline{x}/d_j}}{\sum_{j=1}^{n}{a_j/d_j}}  \quad \quad y_w=\frac{\sum_{j=1}^{n}{a_j \cdot \overline{y}/d_j}}{\sum_{j=1}^{n}{a_j/d_j}}$$
Results

Red circles are customer locations. Their surface area represents the individual demands. The single iterations start with dark blue circled and end with the final warehouse location, colored in green.
Click here to download the data file. Below the code:

require(maps)    #if required install with install.packages("maps")
require(mapdata) #If required install with install.packages("mapdata")
distxy <- br="" calculate="" coorda="" coordb="" custom="" distances="" earth="" function="" his="" is="" our="" to="">  {
    dist <- acos="" br="" coorda="" coordb="" cos="" pi="" sin="">    return(dist)
  }
quartz() # Use windows() when running R on a windows machine
data <- br="" data="" koord.csv="" read.csv="" read="" required="">lim <- border="" br="" c="" plotting="" set="">lim2 <- 55="" br="" c=""># Set iteration start point
xstart <- 45="" br="">ystart <- -110="" br="">number_it <- 20="" br="">iteration_colors <- alpha="1.0," end="0.7)<br" number_it="" rainbow="" start="0.33,">data2 <- data.frame="" x="data[,3]," y="data[,4])<br">xstart_o <- 45="" br="">ystart_o <- -110="" br=""># plot map and customers
it <- data.frame="" x="0," y="0)<br">costs <- br="" c=""># optimization iterations
for (i in 1:number_it)
{
  distance <- 1="" apply="" br="" calculate="" coordb="c(xstart,ystart))" current="" data2="" distances="" distxy="">  nenner <- br="" data="" denominator="" distance="">  zahler_x <- br="" data2="" data="" distance="" nominator_y="">  zahler_y <- br="" data2="" data="" distance="" nominator_y="">  xstart <- br="" nenner="" sum="" zahler_x="">  ystart <- br="" nenner="" sum="" zahler_y="">  it[i,1] <- br="" current="" save="" x_coordinate="" xstart="">  it[i,2] <- br="" current="" save="" y_coordinate="" ystart="">  costs<-append br="" calculate="" costs="" current="" data="" distance="" of="" solution="" sum="">}
# Now plot the results
map("state", lwd=0.3, mar=c(0,0,0,0), col="#ffe895", border=1, interior=TRUE,fill=TRUE, bg="#93c3f2", resolution=0, , xlim = lim, ylim= lim2,)
points(data[,4], data[,3], col="darkred", pch=21,bg="red", cex = sqrt(data[,5]/3)*0.1)
points(ystart_o, xstart_o, col="darkblue", pch=16, cex=1) # starting point
points(it[,2], it[,1],  col = iteration_colors[number_it:1], pch=16)

No hay comentarios: