Solving Optimization Problems with SCIP

Introduction

The scip package provides an R interface to SCIP (Solving Constraint Integer Programs), one of the fastest non-commercial solvers for mixed-integer programming (MIP) and mixed-integer nonlinear programming (MINLP). SCIP supports linear, quadratic, SOS, and indicator constraints with continuous, binary, and integer variables.

The package offers two interfaces:

One-shot interface (scip_solve): pass a constraint matrix and get a solution back in one call.
Model-building interface (scip_model, scip_add_var, etc.): build a model incrementally, useful for complex formulations.

Example 1: A Production Planning LP

A company produces two products. Each unit of Product 1 yields $5 profit and each unit of Product 2 yields $4 profit. Production is limited by labour (6 hours available) and materials (8 kg available). Product 1 requires 1 hour and 2 kg per unit; Product 2 requires 2 hours and 1 kg per unit.

\[\max\; 5x_1 + 4x_2\] subject to: \[x_1 + 2x_2 \le 6 \quad\text{(labour)}\] \[2x_1 + x_2 \le 8 \quad\text{(materials)}\] \[x_1, x_2 \ge 0\]

A <- matrix(c(1, 2,
              2, 1), nrow = 2, byrow = TRUE)
b <- c(6, 8)
## scip_solve minimizes, so negate the objective for maximization
res <- scip_solve(obj = c(-5, -4), A = A, b = b,
                  sense = c("<=", "<="),
                  control = list(verbose = FALSE))
res$status
#> [1] "optimal"
-res$objval  # negate back to get the maximized profit
#> [1] 22
res$x
#> [1] 3.333333 1.333333

Example 2: The Knapsack Problem

Adapted from SCIP’s solveknapsackexactly.c unit test (Marc Pfetsch, Gregor Hendel, Zuse Institute Berlin).

A hiker can carry at most 13 kg. Six items are available with weights 7, 2, 7, 5, 1, and 3 kg, each with equal value. Which items should the hiker pack to carry as many as possible?

\[\max\; \sum_{i=1}^{6} x_i\] subject to: \[7x_1 + 2x_2 + 7x_3 + 5x_4 + x_5 + 3x_6 \le 13\] \[x_i \in \{0,1\}\]

A <- matrix(c(7, 2, 7, 5, 1, 3), nrow = 1)
res <- scip_solve(obj = c(-1, -1, -1, -1, -1, -1),
                  A = A, b = 13, sense = "<=",
                  vtype = "B",
                  control = list(verbose = FALSE))
res$status
#> [1] "optimal"
items_packed <- which(res$x == 1)
items_packed
#> [1] 1 2 5 6
total_weight <- sum(c(7, 2, 7, 5, 1, 3)[items_packed])
total_weight
#> [1] 13

Example 3: N-Queens Problem

Adapted from SCIP’s Queens example (examples/Queens/src/queens.cpp, Cornelius Schwarz, University of Bayreuth).

Place $n$ queens on an $n \times n$ chessboard so that no two queens attack each other. This classic combinatorial problem is naturally modelled as a binary integer program.

Let $x_{i,j} \in \{0,1\}$ indicate whether a queen is placed at row $i$, column $j$. The constraints are:

Rows: exactly one queen per row, $\sum_j x_{i,j} = 1$
Columns: exactly one queen per column, $\sum_i x_{i,j} = 1$
Diagonals: at most one queen per diagonal, $\sum_{(i,j)\in D} x_{i,j} \le 1$

solve_nqueens <- function(n) {
    m <- scip_model(paste0(n, "-queens"))
    scip_set_objective_sense(m, "maximize")

    ## Create n*n binary variables x[i,j] with objective 1
    ## Variable (i,j) is stored at index (i-1)*n + j
    idx <- function(i, j) (i - 1L) * n + j
    scip_add_vars(m, obj = rep(1, n * n), lb = 0, ub = 1, vtype = "B",
                  names = sprintf("x%d_%d", rep(1:n, each = n), rep(1:n, n)))

    ## Row constraints: exactly one queen per row
    for (i in 1:n) {
        scip_add_linear_cons(m, vars = idx(i, 1:n), coefs = rep(1, n),
                             lhs = 1, rhs = 1, name = sprintf("row_%d", i))
    }

    ## Column constraints: exactly one queen per column
    for (j in 1:n) {
        scip_add_linear_cons(m, vars = idx(1:n, j), coefs = rep(1, n),
                             lhs = 1, rhs = 1, name = sprintf("col_%d", j))
    }

    ## Diagonal constraints: at most one queen per diagonal
    ## Down-right diagonals
    for (d in (-(n - 2)):(n - 2)) {
        rows <- max(1, 1 + d):min(n, n + d)
        cols <- rows - d
        if (length(rows) >= 2) {
            scip_add_linear_cons(m, vars = idx(rows, cols),
                                 coefs = rep(1, length(rows)),
                                 rhs = 1, name = sprintf("diag_down_%d", d))
        }
    }
    ## Down-left diagonals (anti-diagonals)
    for (s in 3:(2 * n - 1)) {
        rows <- max(1, s - n):min(n, s - 1)
        cols <- s - rows
        if (length(rows) >= 2) {
            scip_add_linear_cons(m, vars = idx(rows, cols),
                                 coefs = rep(1, length(rows)),
                                 rhs = 1, name = sprintf("diag_up_%d", s))
        }
    }

    scip_optimize(m)
    sol <- scip_get_solution(m)

    ## Extract queen positions
    x_mat <- matrix(round(sol$x), nrow = n, byrow = TRUE)
    positions <- which(x_mat == 1, arr.ind = TRUE)
    colnames(positions) <- c("row", "col")
    list(status = scip_get_status(m),
         n_queens = as.integer(sol$objval),
         positions = positions[order(positions[, "row"]), ],
         board = x_mat)
}

result <- solve_nqueens(8)
result$status
#> [1] "optimal"
result$n_queens
#> [1] 8
result$positions
#>      row col
#> [1,]   1   4
#> [2,]   2   6
#> [3,]   3   1
#> [4,]   4   5
#> [5,]   5   2
#> [6,]   6   8
#> [7,]   7   3
#> [8,]   8   7

Visualize the board (. = empty, Q = queen):

board_str <- apply(result$board, 1, function(row) {
    paste(ifelse(row == 1, "Q", "."), collapse = " ")
})
cat(board_str, sep = "\n")
#> . . . Q . . . .
#> . . . . . Q . .
#> Q . . . . . . .
#> . . . . Q . . .
#> . Q . . . . . .
#> . . . . . . . Q
#> . . Q . . . . .
#> . . . . . . Q .

Example 4: Circle Packing (Quadratic Constraints)

Adapted from SCIP’s CallableLibrary example (examples/CallableLibrary/src/circlepacking.c, Jose Salmeron and Stefan Vigerske, Zuse Institute Berlin).

Pack $n$ circles with given radii into a rectangle of minimum area. For each circle $i$ with radius $r_i$, find center coordinates $(x_i, y_i)$. The rectangle has width $W$ and height $H$.

Constraints:

Circles stay within the rectangle: $r_i \le x_i \le W - r_i$ and $r_i \le y_i \le H - r_i$
Circles do not overlap: $(x_i - x_j)^2 + (y_i - y_j)^2 \ge (r_i + r_j)^2$ for all $i < j$

Objective: Minimize $W \times H$ (area).

Since bilinear objectives are harder to handle, we use a simpler approach: minimize $W + H$ as a proxy and fix one dimension. Here we minimize the width given a fixed height, packing 3 circles.

radii <- c(1.0, 1.5, 0.8)
n <- length(radii)
H <- 5.0  # fixed height

m <- scip_model("circlepacking")

## Variables: x_i, y_i for each circle, plus W
x_idx <- integer(n)
y_idx <- integer(n)
for (i in seq_len(n)) {
    x_idx[i] <- scip_add_var(m, obj = 0, lb = radii[i], ub = 100,
                              name = sprintf("x%d", i))
    y_idx[i] <- scip_add_var(m, obj = 0, lb = radii[i], ub = H - radii[i],
                              name = sprintf("y%d", i))
}
w_idx <- scip_add_var(m, obj = 1, lb = 0, ub = 100, name = "W")

## x_i <= W - r_i  =>  x_i - W <= -r_i
for (i in seq_len(n)) {
    scip_add_linear_cons(m, vars = c(x_idx[i], w_idx), coefs = c(1, -1),
                         rhs = -radii[i],
                         name = sprintf("fit_x_%d", i))
}

## Non-overlap: (x_i - x_j)^2 + (y_i - y_j)^2 >= (r_i + r_j)^2
## Expand: x_i^2 - 2*x_i*x_j + x_j^2 + y_i^2 - 2*y_i*y_j + y_j^2 >= (r_i+r_j)^2
for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
        min_dist_sq <- (radii[i] + radii[j])^2
        scip_add_quadratic_cons(m,
            quadvars1 = c(x_idx[i], x_idx[i], x_idx[j], y_idx[i], y_idx[i], y_idx[j]),
            quadvars2 = c(x_idx[i], x_idx[j], x_idx[j], y_idx[i], y_idx[j], y_idx[j]),
            quadcoefs = c(1, -2, 1, 1, -2, 1),
            lhs = min_dist_sq,
            name = sprintf("nooverlap_%d_%d", i, j))
    }
}

scip_optimize(m)
cat("Status:", scip_get_status(m), "\n")
#> Status: optimal
sol <- scip_get_solution(m)
W_opt <- sol$x[w_idx]
cat(sprintf("Minimum width: %.3f (height fixed at %.1f)\n", W_opt, H))
#> Minimum width: 3.515 (height fixed at 5.0)
for (i in seq_len(n)) {
    cat(sprintf("  Circle %d (r=%.1f): center = (%.3f, %.3f)\n",
                i, radii[i], sol$x[x_idx[i]], sol$x[y_idx[i]]))
}
#>   Circle 1 (r=1.0): center = (1.000, 4.000)
#>   Circle 2 (r=1.5): center = (1.500, 1.500)
#>   Circle 3 (r=0.8): center = (2.715, 3.453)

Example 5: Facility Location (Mixed-Integer)

Adapted from SCIP’s SCFLP example (examples/SCFLP/doc/xternal_scflp.c, Zuse Institute Berlin).

A company must decide which of $p$ potential warehouse locations to open. Each warehouse $i$ has a fixed opening cost $f_i$ and a capacity $k_i$. Each of $q$ customers $j$ has a demand $d_j$ and a per-unit shipping cost $c_{ij}$ from warehouse $i$. Minimize total cost.

\[\min\; \sum_i f_i y_i + \sum_{i,j} c_{ij} x_{ij}\] subject to: \[\sum_i x_{ij} \ge d_j \quad\forall j \quad\text{(demand)}\] \[\sum_j x_{ij} \le k_i\, y_i \quad\forall i \quad\text{(capacity)}\] \[y_i \in \{0,1\},\; x_{ij} \ge 0\]

## Problem data
p <- 3  # warehouses
q <- 4  # customers
fixed_cost <- c(100, 150, 120)
capacity <- c(50, 60, 40)
demand <- c(20, 25, 15, 30)
## Shipping cost matrix (warehouses x customers)
ship_cost <- matrix(c(
    4, 8, 5, 6,
    6, 3, 7, 4,
    5, 5, 4, 8
), nrow = p, byrow = TRUE)

m <- scip_model("facility_location")

## Binary variables y_i: open warehouse i?
y_idx <- integer(p)
for (i in 1:p) {
    y_idx[i] <- scip_add_var(m, obj = fixed_cost[i], vtype = "B",
                              name = sprintf("y%d", i))
}

## Continuous variables x_ij: amount shipped from i to j
x_idx <- matrix(0L, nrow = p, ncol = q)
for (i in 1:p) {
    for (j in 1:q) {
        x_idx[i, j] <- scip_add_var(m, obj = ship_cost[i, j],
                                     lb = 0, ub = max(capacity),
                                     name = sprintf("x%d_%d", i, j))
    }
}

## Demand constraints: sum_i x_ij >= d_j
for (j in 1:q) {
    scip_add_linear_cons(m, vars = x_idx[, j], coefs = rep(1, p),
                         lhs = demand[j],
                         name = sprintf("demand_%d", j))
}

## Capacity constraints: sum_j x_ij - k_i * y_i <= 0
for (i in 1:p) {
    scip_add_linear_cons(m,
                         vars = c(x_idx[i, ], y_idx[i]),
                         coefs = c(rep(1, q), -capacity[i]),
                         rhs = 0,
                         name = sprintf("capacity_%d", i))
}

scip_optimize(m)
sol <- scip_get_solution(m)
cat("Status:", scip_get_status(m), "\n")
#> Status: optimal
cat("Total cost:", sol$objval, "\n")
#> Total cost: 600

open_warehouses <- which(round(sol$x[y_idx]) == 1)
cat("Open warehouses:", open_warehouses, "\n")
#> Open warehouses: 1 2

shipments <- matrix(sol$x[x_idx], nrow = p)
cat("\nShipment plan (rows=warehouses, cols=customers):\n")
#> 
#> Shipment plan (rows=warehouses, cols=customers):
rownames(shipments) <- paste0("W", 1:p)
colnames(shipments) <- paste0("C", 1:q)
print(round(shipments, 1))
#>    C1 C2 C3 C4
#> W1 20  0 15  0
#> W2  0 25  0 30
#> W3  0  0  0  0

Example 6: Indicator Constraints

Adapted from SCIP’s indicator test instances (check/instances/Indicator/, Zuse Institute Berlin).

Indicator constraints model logical implications: “if a binary variable is 1, then a linear constraint must hold.” They are widely used in scheduling, network design, and disjunctive programming.

Here we model a simple production problem with setup costs: a machine must be “turned on” ($z_i = 1$) before it can produce. Turning on a machine costs $50. Each unit produced on machine $i$ earns revenue $r_i$. Machine $i$ can produce at most $u_i$ units, but only if turned on.

\[\max\; \sum_i (r_i x_i - 50 z_i)\] subject to: \[z_i = 1 \implies x_i \le u_i \quad \text{(capacity when on)}\] \[z_i = 0 \implies x_i = 0 \quad \text{(nothing when off)}\] \[x_i \ge 0,\; z_i \in \{0,1\}\]

We model $z_i = 0 \implies x_i \le 0$ as an indicator constraint on the negated variable. Equivalently, we add the big-M constraint $x_i \le u_i z_i$.

revenue <- c(12, 8, 15)
max_prod <- c(10, 20, 8)
setup_cost <- 50

m <- scip_model("setup_production")
scip_set_objective_sense(m, "maximize")

z_idx <- integer(3)
x_idx <- integer(3)
for (i in 1:3) {
    z_idx[i] <- scip_add_var(m, obj = -setup_cost, vtype = "B",
                              name = sprintf("z%d", i))
    x_idx[i] <- scip_add_var(m, obj = revenue[i], lb = 0, ub = max_prod[i],
                              name = sprintf("x%d", i))
    ## x_i <= max_prod[i] * z_i  =>  x_i - max_prod[i] * z_i <= 0
    scip_add_linear_cons(m, vars = c(x_idx[i], z_idx[i]),
                         coefs = c(1, -max_prod[i]),
                         rhs = 0,
                         name = sprintf("link_%d", i))
}

scip_optimize(m)
sol <- scip_get_solution(m)
cat("Status:", scip_get_status(m), "\n")
#> Status: optimal
cat("Profit:", sol$objval, "\n")
#> Profit: 250
for (i in 1:3) {
    on <- round(sol$x[z_idx[i]])
    prod <- sol$x[x_idx[i]]
    cat(sprintf("  Machine %d: %s, produce %.0f units\n",
                i, if (on) "ON" else "OFF", prod))
}
#>   Machine 1: ON, produce 10 units
#>   Machine 2: ON, produce 20 units
#>   Machine 3: ON, produce 8 units

Solver Controls

Use scip_control() with the one-shot interface, or scip_set_param() with the model-building interface, to tune solver behavior.

## One-shot: time limit and gap tolerance
ctrl <- scip_control(verbose = FALSE, time_limit = 60, gap_limit = 0.01)

## Model-building: set SCIP parameters directly
m <- scip_model("tuning_example")
scip_set_param(m, "display/verblevel", 0L)   # suppress output
scip_set_param(m, "limits/time", 60.0)        # 60-second time limit
scip_set_param(m, "limits/gap", 0.01)         # 1% optimality gap

Common SCIP parameters:

Parameter	Type	Description
`display/verblevel`	int	Verbosity (0=silent, 5=max)
`limits/time`	real	Time limit in seconds
`limits/gap`	real	Relative MIP gap tolerance
`limits/nodes`	longint	Maximum B&B nodes
`limits/solutions`	int	Stop after this many solutions

See the SCIP parameter documentation for the complete list.