Hide navigation sidebar
Hide table of contents sidebar
Back to top
Toggle table of contents sidebar

# facility.R

    # Copyright 2025, Gurobi Optimization, LLC
    #
    # Facility location: a company currently ships its product from 5 plants
    # to 4 warehouses. It is considering closing some plants to reduce
    # costs. What plant(s) should the company close, in order to minimize
    # transportation and fixed costs?
    #
    # Note that this example uses lists instead of dictionaries.  Since
    # it does not work with sparse data, lists are a reasonable option.
    #
    # Based on an example from Frontline Systems:
    #   http://www.solver.com/disfacility.htm
    # Used with permission.

    library(Matrix)
    library(gurobi)

    # define primitive data
    nPlants     <- 5
    nWarehouses <- 4
    # Warehouse demand in thousands of units
    Demand      <- c(15, 18, 14, 20)
    # Plant capacity in thousands of units
    Capacity    <- c(20, 22, 17, 19, 18)
    # Fixed costs for each plant
    FixedCosts  <- c( 12000, 15000, 17000, 13000, 16000)
    # Transportation costs per thousand units
    TransCosts  <- c(4000, 2000, 3000, 2500, 4500,
                     2500, 2600, 3400, 3000, 4000,
                     1200, 1800, 2600, 4100, 3000,
                     2200, 2600, 3100, 3700, 3200 )

    flowidx <- function(w, p) {nPlants * (w-1) + p}

    # Build model
    model <- list()
    model$modelname <- 'facility'
    model$modelsense <- 'min'

    # initialize data for variables
    model$lb       <- 0
    model$ub       <- c(rep(1, nPlants),   rep(Inf, nPlants * nWarehouses))
    model$vtype    <- c(rep('B', nPlants), rep('C', nPlants * nWarehouses))
    model$obj      <- c(FixedCosts, TransCosts)
    model$varnames <- c(paste0(rep('Open',nPlants),1:nPlants),
                        sprintf('Trans%d,%d',
                                c(mapply(rep,1:nWarehouses,nPlants)),
                                1:nPlants))

    # build production constraint matrix
    A1 <- spMatrix(nPlants, nPlants, i = c(1:nPlants), j = (1:nPlants), x = -Capacity)
    A2 <- spMatrix(nPlants, nPlants * nWarehouses,
                   i = c(mapply(rep, 1:nPlants, nWarehouses)),
                   j = mapply(flowidx,1:nWarehouses,c(mapply(rep,1:nPlants,nWarehouses))),
                   x = rep(1, nWarehouses * nPlants))
    A3 <- spMatrix(nWarehouses, nPlants)
    A4 <- spMatrix(nWarehouses, nPlants * nWarehouses,
                   i = c(mapply(rep, 1:nWarehouses, nPlants)),
                   j = mapply(flowidx,c(mapply(rep,1:nWarehouses,nPlants)),1:nPlants),
                   x = rep(1, nPlants * nWarehouses))
    model$A           <- rbind(cbind(A1, A2), cbind(A3, A4))
    model$rhs         <- c(rep(0, nPlants),   Demand)
    model$sense       <- c(rep('<', nPlants), rep('=', nWarehouses))
    model$constrnames <- c(sprintf('Capacity%d',1:nPlants),
                           sprintf('Demand%d',1:nWarehouses))

    # Save model
    gurobi_write(model,'facilityR.lp')

    # Guess at the starting point: close the plant with the highest fixed
    # costs; open all others first open all plants
    model$start <- c(rep(1,nPlants),rep(NA, nPlants * nWarehouses))

    # find most expensive plant, and close it in mipstart
    cat('Initial guess:\n')
    worstidx <- which.max(FixedCosts)
    model$start[worstidx] <- 0
    cat('Closing plant',worstidx,'\n')

    # set parameters
    params <- list()
    params$method <- 2

    # Optimize
    res <- gurobi(model, params)

    # Print solution
    if (res$status == 'OPTIMAL') {
      cat('\nTotal Costs:',res$objval,'\nsolution:\n')
      cat('Facilities:', model$varnames[which(res$x[1:nPlants]>1e-5)], '\n')
      active <- nPlants + which(res$x[(nPlants+1):(nPlants*(nWarehouses+1))] > 1e-5)
      cat('Flows: ')
      cat(sprintf('%s=%g ',model$varnames[active], res$x[active]), '\n')
      rm(active)
    } else {
      cat('No solution\n')
    }

    # Clear space
    rm(res, model, params, A1, A2, A3, A4)

Help and Feedback

Press `/` or click to search this documentation | Gurobi website | Gurobi community forum
Go to other Gurobi documentation | Gurobi download | Contact support