---
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: <https://www.scipopt.org/>
- 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).
