--- title: "Solving Optimization Problems with SCIP" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Solving Optimization Problems with SCIP} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") library(scip) ``` ## Introduction The `scip` package provides an R interface to [SCIP](https://www.scipopt.org/) (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. This vignette demonstrates the package through examples adapted from the [SCIP documentation and example suite](https://www.scipopt.org/doc/html/EXAMPLES.php) (Copyright 2002--2026 Zuse Institute Berlin, Apache-2.0 license). 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$$ ```{r production-lp} 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 -res$objval # negate back to get the maximized profit res$x ``` ## 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\}$$ ```{r knapsack} 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 items_packed <- which(res$x == 1) items_packed total_weight <- sum(c(7, 2, 7, 5, 1, 3)[items_packed]) total_weight ``` ## 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$ ```{r nqueens-function} 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) } ``` ```{r nqueens-solve} result <- solve_nqueens(8) result$status result$n_queens result$positions ``` Visualize the board (`.` = empty, `Q` = queen): ```{r nqueens-board} board_str <- apply(result$board, 1, function(row) { paste(ifelse(row == 1, "Q", "."), collapse = " ") }) cat(board_str, sep = "\n") ``` ## 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. ```{r circlepacking} 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") sol <- scip_get_solution(m) W_opt <- sol$x[w_idx] cat(sprintf("Minimum width: %.3f (height fixed at %.1f)\n", W_opt, H)) 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]])) } ``` ## 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$$ ```{r facility-location} ## 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") cat("Total cost:", sol$objval, "\n") open_warehouses <- which(round(sol$x[y_idx]) == 1) cat("Open warehouses:", open_warehouses, "\n") shipments <- matrix(sol$x[x_idx], nrow = p) cat("\nShipment plan (rows=warehouses, cols=customers):\n") rownames(shipments) <- paste0("W", 1:p) colnames(shipments) <- paste0("C", 1:q) print(round(shipments, 1)) ``` ## 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$. ```{r indicator} 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") cat("Profit:", sol$objval, "\n") 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)) } ``` ## Solver Controls Use `scip_control()` with the one-shot interface, or `scip_set_param()` with the model-building interface, to tune solver behavior. ```{r controls} ## 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](https://www.scipopt.org/doc/html/PARAMETERS.php) for the complete list. ## References - SCIP Optimization Suite: - Bestuzheva, K., Besançon, M., Chen, W.-K., Chmiela, A., Donkiewicz, T., van Doornmalen, J., et al. (2021). The SCIP Optimization Suite 8.0. *ZIB-Report* 21-41, Zuse Institute Berlin. - Examples in this vignette are adapted from the SCIP example suite (Copyright 2002--2026 Zuse Institute Berlin, Apache-2.0 license), available at `examples/` in the [SCIP source](https://github.com/scipopt/scip).