The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

Linear programming with Clp

Clp is the linear programming code of the COIN-OR project: a simplex implementation with a barrier alternative, presolve, and warm starts. This package binds its callable library.

clp_version()$version
#> [1] "1.17.0"

Solving a problem in one call

A small product mix problem. Two products, three resource constraints, maximise the contribution.

A <- rbind(material = c(120, 210),
           labour   = c(110,  30),
           capacity = c(  1,   1))
b <- c(15000, 4000, 75)

fit <- clp_solve(c(143, 60), A, "<=", b, max = TRUE,
                 col_names = c("x", "y"), row_names = rownames(A))
fit
#> <clp_solution>
#>   status:    0 (optimal)
#>   objective: 6315.625 [max]
#>   solution:  21.875 53.125
#>   iterations: 2

The solution, the shadow prices and the reduced costs come back together.

fit$solution
#>      x      y 
#> 21.875 53.125
fit$duals
#> material   labour capacity 
#>   0.0000   1.0375  28.8750
fit$reduced_costs
#> [1] 0 0

A zero dual marks a resource that is not binding: here the material constraint has slack, while labour and capacity are tight.

data.frame(row = rownames(A), activity = fit$row_activity, limit = b,
           dual = fit$duals)
#>               row activity limit    dual
#> material material 13781.25 15000  0.0000
#> labour     labour  4000.00  4000  1.0375
#> capacity capacity    75.00    75 28.8750

Ranged constraints and free variables

dir and rhs cover one-sided rows. For a row bounded on both sides, give row_lower and row_upper instead. Variable bounds default to [0, Inf) and are set with lower and upper; -Inf makes a variable free.

clp_solve(c(1, 1), rbind(c(1, 1)),
          row_lower = 2, row_upper = 5,
          lower = -Inf, upper = Inf)$objval
#> [1] 2

Sparse input

Real models are sparse. Any of a Matrix sparse matrix, a slam::simple_triplet_matrix or plain triplets can be passed straight in. A sparse matrix goes to Clp as sparse: only the non-zero entries are passed, and the matrix is never expanded into a full grid of mostly zeros first, which for a large model is the difference between fitting in memory and not.

set.seed(1)
big <- Matrix::rsparsematrix(200, 500, density = 0.01)
res <- clp_solve(rep(1, 500), big, ">=", rep(-1, 200), lower = 0, upper = 10)
res$status_message
#> [1] "optimal"

Keeping a model

clp_solve() builds a model, solves it and throws it away. When a problem is solved repeatedly with small changes – a parametric study, a column generation loop, a rolling horizon – build the model once and re-solve it, so Clp can start from the basis it already has.

model <- clp_model()
clp_set_log_level(model, 0)
clp_load_problem(model, ncols = 2, nrows = 3,
                 start = c(0L, 3L, 6L),
                 index = c(0L, 1L, 2L, 0L, 1L, 2L),
                 value = c(120, 110, 1, 210, 30, 1),
                 obj   = c(-143, -60),
                 rowub = b)
clp_initial_solve(model)
#> [1] 0
c(objective = clp_objective_value(model), iterations = clp_iterations(model))
#>  objective iterations 
#>  -6315.625      2.000

The problem is loaded column by column in compressed sparse column form: start says where each column begins, index holds 0-based row positions and value the coefficients. Because the low level functions follow the C API, positions are 0-based here, as in the Clp documentation. Maximisation is expressed by negating the objective, or by clp_set_optimization_direction(model, -1).

Warm starts

Keep the basis, change the model, hand the basis back:

basis <- clp_status_array(model)

clp_set_row_upper(model, c(15000, 4000, 70))
clp_copyin_status(model, basis)
clp_dual_simplex(model)
#> [1] 0

c(objective = clp_objective_value(model), iterations = clp_iterations(model))
#>  objective iterations 
#>   -6171.25       0.00

The re-solve takes no iterations at all: the old basis is still optimal for the tightened problem.

Choosing the algorithm

for (alg in c("auto", "primal", "dual", "barrier")) {
    fit <- clp_solve(c(143, 60), A, "<=", b, max = TRUE,
                     control = clp_control(algorithm = alg))
    cat(sprintf("%-8s %.4f\n", alg, fit$objval))
}
#> auto     6315.6250
#> primal   6315.6250
#> dual     6315.6250
#> barrier  6315.6250

clp_control() also carries the tolerances, the iteration and time limits, the scaling mode and whether to presolve. For finer control over presolve there is a clp_options() object with the individual transformations (clp_options_set_do_dupcol() and friends), passed to clp_initial_solve_with_options().

MPS files

path <- system.file("extdata", "productmix.mps", package = "coinclp")
from_file <- clp_model()
clp_set_log_level(from_file, 0)
clp_read_mps(from_file, path)
clp_col_names(from_file)
#> [1] "x" "y"
clp_initial_solve(from_file)
#> [1] 0
clp_objective_value(from_file)
#> [1] -6315.625

Writing works the other way with clp_write_mps(). Clp only gained an MPS writer in its C API after the 1.17 series, so on a Clp without one – the version Rtools ships, for instance – the package writes the file itself. clp_features() says which entry points the current build has.

clp_features()
#>          write_mps modify_coefficient          set_names maximum_iterations 
#>              FALSE              FALSE              FALSE               TRUE

Coming from clpAPI

The archived clpAPI package is reproduced function for function, so code written against it runs here unchanged.

lp <- initProbCLP()
setLogLevelCLP(lp, 0)
loadProblemCLP(lp, 2, 3, c(0, 3, 6), c(0, 1, 2, 0, 1, 2),
               c(120, 110, 1, 210, 30, 1),
               lb = c(0, 0), ub = c(1e30, 1e30), obj_coef = c(143, 60),
               rlb = rep(-1e30, 3), rub = b)
setObjDirCLP(lp, -1)
solveInitialCLP(lp)
#> [1] 0
status_codeCLP(getSolStatusCLP(lp))
#> [1] "solution is optimal"
getObjValCLP(lp)
#> [1] 6315.625
delProbCLP(lp)

Note that clpAPI used 1e30 for an infinite bound, which is Clp’s own convention and what clp_inf() returns. The clp_solve() interface accepts Inf and converts it.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.