Linear_Programs

Linear Programs

A scenario linear program (LP) minimises a linear cost subject to linear constraints that depend on an uncertain parameter $\delta \in \mathbb{R}^q$. The distribution of $\delta$ is unknown; only $N$ independent samples $\delta_1, \dots, \delta_N$ (the scenarios) are available. Scen-Opt solves the program built from these scenarios and then reports the complexity $k$ and the risk bounds $[\underline{\varepsilon}, \overline{\varepsilon}]$ of the solution. See Scenario Theory for what these quantities mean.

This page covers the LP formulation, its inputs in the web interface, in the program JSON file and in Python, and a worked example. For a first walk-through of the interface itself, see the Web Interface example; for Python usage, see the Python example.

Formulation

The most general LP is the relaxed–regularized LP:

$$\begin{aligned} \min_{x \in \mathbb{R}^d,\ \zeta_1,\dots,\zeta_N \ge 0} \quad & c^\top x + \tau \lVert x - \bar{x} \rVert_p + \rho \sum_{i=1}^{N} \zeta_i \\\ \text{subject to} \quad & A(\delta_i)\, x + b(\delta_i) \le \zeta_i, \qquad i = 1, \dots, N, \\\ & G\, x + h \le 0, \end{aligned}$$

with $\tau \ge 0$ and $\rho \ge 0$. Inequalities between vectors hold entry-wise.

Symbol	Meaning
$x \in \mathbb{R}^d$	decision vector
$c \in \mathbb{R}^d$	cost vector
$A(\delta) \in \mathbb{R}^{\mathfrak{m} \times d}$, $b(\delta) \in \mathbb{R}^{\mathfrak{m}}$	scenario constraints: each scenario $\delta_i$ adds the $\mathfrak{m}$ rows $A(\delta_i)x + b(\delta_i) \le \zeta_i$
$G \in \mathbb{R}^{\mathfrak{n} \times d}$, $h \in \mathbb{R}^{\mathfrak{n}}$	hard constraints, which do not depend on $\delta$ ($\mathfrak{n} = 0$ means none)
$\zeta_i \ge 0$	slack of scenario $i$: one scalar per scenario, shared by all $\mathfrak{m}$ rows of that scenario
$\rho \ge 0$	penalty on the slacks (relaxation)
$\tau \ge 0$, $\bar{x} \in \mathbb{R}^d$, $p$	weight, reference point and norm of the regularization term

The hard constraints are always enforced. They are never part of the support list, so they do not count towards the complexity $k$. Each scenario counts as one constraint, whatever the value of $\mathfrak{m}$. As a result, $k$ counts scenarios, not rows.

The four variants

Each variant is a special case of the formulation above. The hard constraints $Gx + h \le 0$ are present in all four.

Variant	`option` value	Objective	Scenario constraints
Robust	`robust`	$c^\top x$	$A(\delta_i)x + b(\delta_i) \le 0$
Regularized	`regularization`	$c^\top x + \tau \lVert x - \bar{x} \rVert_p$	$A(\delta_i)x + b(\delta_i) \le 0$
Relaxed	`relaxation`	$c^\top x + \rho \sum_i \zeta_i$	$A(\delta_i)x + b(\delta_i) \le \zeta_i$, $\zeta_i \ge 0$
Relaxed–regularized	`regularization-relaxation`	$c^\top x + \tau \lVert x - \bar{x} \rVert_p + \rho \sum_i \zeta_i$	$A(\delta_i)x + b(\delta_i) \le \zeta_i$, $\zeta_i \ge 0$

Robust: $\tau = 0$, $\rho = 0$ and $\zeta_i = 0$ for all $i$. Every scenario constraint must hold.
Relaxed: $\tau = 0$. Scenario constraints may be violated, and each violation is penalised by $\rho \zeta_i$.
Regularized: $\rho = 0$ and $\zeta_i = 0$. The solution is pulled towards $\bar{x}$.

How the variant is selected:

Web interface: the Select Option dropdown sets the variant, and it shows only the parameters that variant uses ($\tau$, $\bar{x}$ and $p$ for regularization; $\rho$ for relaxation). Changing the option clears $\tau$, $\rho$, $p$ and $\bar{x}$.
Relaxed options with $\rho = 0$: the slacks are kept even when $\rho = 0$. The scenario constraints can then be violated at no cost, and an LP is usually unbounded. The error column then shows LP did not solve to optimality. Status: unbounded, .... Use $\rho > 0$.
Python: solve_lp has no option argument. Slack variables are added when include_slack=True. With the default include_slack=None, they are added whenever rho != 0. The regularization term is active whenever tau > 0.

Regularization norm

You can set the norm $p$ in the web field p, in the JSON key p, or with the Python argument norm_type. If you leave it empty, $p = 2$. The web field accepts a single number $\ge 1$, inf or fro.

For LPs (and QPs), the decision variable is stored as a $d \times 1$ column, and CVXPY evaluates the norm as a matrix norm. For a column vector, $p = 1$, $2$ and inf give the usual vector norms, and fro equals the 2-norm. Other orders, for example $p = 3$, fail with the error Unsupported matrix norm. In Python, use norm_type=1, 2, np.inf or 'fro'.

Inputs

Problem data

Symbol	Shape	Web interface	Program JSON key	`solve_lp` argument
$d$	integer $\ge 1$	dimension field d	`n_x`	inferred from `A_d`
$\mathfrak{m}$	integer $\ge 1$	dimension field 𝔪	`rows_A`	inferred from `A_d`
$\mathfrak{n}$	integer $\ge 0$	dimension field 𝔫	`rows_G`	inferred from `G`
$c$	$d \times 1$	Edit c	`c`	`c`: array `(d, 1)` or `(d,)`
$A(\delta)$	$\mathfrak{m} \times d$, may use `delta[k]`	Edit A(δ) (symbolic mode only)	`A_d`	`A_d`: callable, `delta -> (m, d)` array
$b(\delta)$	$\mathfrak{m} \times 1$, may use `delta[k]`	Edit b(δ) (symbolic mode only)	`b_d`	`b_d`: callable, `delta -> (m, 1)` array
$G$	$\mathfrak{n} \times d$, numeric	Edit G (enabled when 𝔫 > 0)	`G`	`G`: array `(n, d)`, or `np.array([])`
$h$	$\mathfrak{n} \times 1$, numeric	Edit h (enabled when 𝔫 > 0)	`h`	`h`: array `(n, 1)`, or `np.array([])`
scenarios	see input modes	Scenarios upload box	separate file	`deltas`: 2-D array `(N, q)`

Parameters

Symbol	Web interface	Program JSON key	Python	Default
$\tau$	τ ≥ 0 (comma-separated list allowed)	`tau`	`tau` (float)	0
$\bar{x}$	Edit x̄ ($d \times 1$)	`x_ref`	`x_ref`: scalar (applied to every entry) or `d` entries	0
$p$	p	`p`	`norm_type`	2
$\rho$	ρ ≥ 0 (comma-separated list allowed)	`rho`	`rho` (float)	0
$\beta$	Confidence ... where β	`confidence`	`beta` argument of `quantify_risk(k, N, beta)`	none; choose $0 < \beta < 1$, e.g. `1e-06`
solver	Solver dropdown	not stored	`solver`, e.g. `"CLARABEL"`	CVXPY's choice (`None`)

Notes for Python callers

deltas must be 2-D, with shape (N, q), even when $q = 1$, because the expressions index delta[k]. For example, use np.loadtxt(path, delimiter=",", ndmin=2).
b_d(delta) must return a column of shape (m, 1), and h must have shape (n, 1). A 1-D array is broadcast against the (m, 1) product $A(\delta_i)x$ and silently changes the constraints. For example, a 1-D b_d in the example below returns a wrong interval.
Hard constraints are used only when both G and h are non-empty.
x_ref must contain either 1 or d entries; any other size raises a ValueError.
solve_lp(deltas, A_d, b_d, G, h, c, tau=0.0, x_ref=np.array([0.0]), rho=0.0, norm_type=2, solver=None, include_slack=None) returns (x, zeta, cost, N, complexity, constraints, degeneracy):
- x has shape (d, 1).
- zeta has shape (N,) and is all zeros without slacks.
- cost is the optimal value of the full objective, including the $\tau$ and $\rho$ terms. The web interface reports only $c^\top x^\ast$.
- constraints holds one CVXPY constraint per scenario.
Pass complexity and N to quantify_risk to get the risk bounds. Implementation details are on the LP.py and Risk.py pages.

Writing scenario-dependent entries

In symbolic mode, each entry of $A(\delta)$ and $b(\delta)$ is a string expression in the components delta[0], delta[1], …, delta[q-1] of a scenario. Indexing starts at 0: delta[0] is the first column of the scenario file. Expressions are evaluated by a restricted parser, which allows:

numbers, + - * / // % **, unary minus and parentheses;
delta[k];
functions and constants of Python's math module, written with the prefix, e.g. math.sin(delta[1]), math.exp(-delta[0]) or math.pi * delta[2].

Anything else is rejected. For example, a bare sin(delta[1]) gives the error disallowed name: sin, and so does abs(...) (use math.fabs(...) instead). The matrices $c$, $G$ and $h$ must not contain delta.

In the program JSON, write every matrix entry as a string, constants included ("0", "1.5", "delta[0]"), as the benchmark files do. Plain JSON numbers are not accepted for c, A_d, b_d, G or h in v1.0.

Symbolic and numeric input modes

Use the Symbolic / Numeric toggle in the web interface, or the JSON key mode, to pick how scenario data are supplied.

Symbolic (default). $A(\delta)$ and $b(\delta)$ are entered as expressions. The scenario file has one row per scenario and $q$ columns, delta[0], ..., delta[q-1], so it is an $N \times q$ table.
Numeric. The Edit A(δ) and Edit b(δ) buttons are hidden. Instead, row $i$ of the scenario file holds the already evaluated $A_i = A(\delta_i)$ and $b_i = b(\delta_i)$: first $A_i$ flattened row by row ($A_i[0,0], A_i[0,1], \dots, A_i[0,d-1], A_i[1,0], \dots$), then the $\mathfrak{m}$ entries of $b_i$. Each row therefore has $\mathfrak{m} d + \mathfrak{m}$ numbers. The server checks this and otherwise reports Each scenario row must have ... elements. Because there is no $A(\delta)$ to infer them from, the dimensions $d$ and $\mathfrak{m}$ must be set, through the dimension fields or the JSON keys n_x/c and rows_A.

Both modes solve the same program. The example below gives identical results in either mode. Scenario files are uploaded as CSV without a header row, or in another format accepted by the upload box (see the Web Interface example).

Parameter sweeps

The $\rho$ and $\tau$ fields accept comma-separated lists, for example 0.02, 0.1, 1; $p$ cannot be a list. The tool solves one program for each $(\tau, \rho)$ pair, $n_\tau n_\rho$ programs in total. Each has its own complexity and risk bounds, computed with the β you entered.

By the union bound, all of these bounds hold together with confidence at least $1 - \beta n_\tau n_\rho$. The interface shows this formula next to the confidence field and its value in the results table. The results also include $\rho$ and $\tau$ graphs (see In-built Graphics). In Python, loop over the values yourself.

Program JSON

A program file holds everything except the scenarios and the solver. Load it with Upload Program, next to the LP/QP/SDP tabs: either choose the file or paste its content, then press Load. After a solve, the Download Program button saves the current problem in the same format.

Key	Required	Content
`type`	yes	`"LP"`
`mode`	no	`"symbolic"` (default) or `"numeric"`
`option`	no	`"robust"`, `"regularization"`, `"relaxation"` or `"regularization-relaxation"`. If absent, it is inferred from the first value of `rho` and `tau`: both > 0 gives `regularization-relaxation`, only `rho` > 0 gives `relaxation`, only `tau` > 0 gives `regularization`, otherwise `robust`.
`c`	yes	$d \times 1$ array of strings
`A_d`, `b_d`	symbolic mode	$\mathfrak{m} \times d$ and $\mathfrak{m} \times 1$ arrays of strings
`G`, `h`	no (both or neither)	$\mathfrak{n} \times d$ and $\mathfrak{n} \times 1$ arrays of strings
`n_x`, `rows_A`, `rows_G`	`rows_A` in numeric mode	integers $d$, $\mathfrak{m}$, $\mathfrak{n}$. When `c`, `A_d` or `G` are present, the dimensions are taken from their shapes.
`rho`, `tau`	no	strings, e.g. `"0"` or `"0.1, 1, 10"`
`confidence`	no	β as a string, e.g. `"1e-06"`
`p`	no	`"1"`, `"2"`, `"inf"` or `"fro"`
`x_ref`	no	$d \times 1$ array, e.g. `[["0"], ["0"]]`

Optional fields that a problem does not use can simply be omitted. For the robust option, leave tau and rho at "0" or omit them.

Worked example: smallest enclosing interval

Benchmark: benchmarks/LP_half_width_2d.

The goal is the smallest interval $[o - w,\ o + w]$ that contains $N = 100$ points $\delta_i$ sampled from $[0, 1]$. The decision vector is $x = [o,\ w]^\top$ (centre, half-width), and the cost $c = [0,\ 1]^\top$ minimises the half-width. Each point adds two rows:

$$A(\delta_i) = \begin{bmatrix} -1 & -1 \\ 1 & -1 \end{bmatrix}, \qquad b(\delta_i) = \begin{bmatrix} \delta_i \\ -\delta_i \end{bmatrix},$$

so that $A(\delta_i)x + b(\delta_i) \le 0$ reads $\delta_i \le o + w$ and $o - w \le \delta_i$. There are no hard constraints ($\mathfrak{n} = 0$). The program is robust ($\tau = \rho = 0$), and $\beta = 10^{-6}$.

Program file

data/program_symbolic.json (whitespace reformatted):

{
  "type": "LP",
  "A_d": [["-1", "-1"], ["1", "-1"]],
  "b_d": [["delta[0]"], ["-delta[0]"]],
  "c": [["0"], ["1"]],
  "rho": "0",
  "tau": "0",
  "confidence": "1e-06",
  "mode": "symbolic",
  "n_x": 2,
  "rows_A": 2,
  "rows_G": 0
}

The scenario file data/scenarios.csv has 100 rows with one number each (delta[0]).

In the web interface

Click Upload Program, choose data/program_symbolic.json and press Load. The LP tab opens with option Robust, $d = 2$, 𝔪 = 2 and 𝔫 = 0.
Under Scenarios, upload data/scenarios.csv.
Choose a solver (CLARABEL was used for the numbers below) and press Solve.

The results table shows:

optimal cost 0.4983;
optimal x = (0.5006, 0.4983);
complexity 2;
risk bounds [0, 0.2085];
confidence 0.999999;
no degeneracy.

In Python

Run the snippet from the root of the Scen-Opt repository, so that src and benchmarks/ can be found. It reads the same JSON file and converts it with the parser that the web interface uses.

import json
import numpy as np
from src.LP import solve_lp
from src.Risk import quantify_risk
from src.parsing import generate_matrix_function, generate_matrix

data = "benchmarks/LP_half_width_2d/data/"
prog = json.load(open(data + "program_symbolic.json"))

# Turn the JSON entries into Python objects, exactly as the web interface does
A_d = generate_matrix_function(json.dumps(prog["A_d"]))  # callable: delta -> (2, 2) array
b_d = generate_matrix_function(json.dumps(prog["b_d"]))  # callable: delta -> (2, 1) array
c = generate_matrix(json.dumps(prog["c"]))               # (2, 1) array
G, h = np.array([]), np.array([])                        # no hard constraints

deltas = np.loadtxt(data + "scenarios.csv", delimiter=",", ndmin=2)  # shape (100, 1)

# Robust LP: tau = 0 and rho = 0 (the defaults)
x, zeta, cost, N, k, constraints, degenerate = solve_lp(
    deltas, A_d, b_d, G, h, c, solver="CLARABEL")
eps_L, eps_U = quantify_risk(k, N, beta=1e-6)

print("x* = [center, half-width]:", x.ravel())
print("cost:", cost)
print("N =", N, "| k =", k, "| degenerate:", degenerate)
print("risk bounds:", float(eps_L), float(eps_U))

Output (CVXPY 1.7.2, CLARABEL):

Discarded 1 spuriously screened constraint(s) (strictly slack, near-threshold duals); remaining support scenarios are irreducible: non-degenerate.
x* = [center, half-width]: [0.5006347  0.49830997]
cost: 0.49830996866324617
N = 100 | k = 2 | degenerate: False
risk bounds: 0.0 0.20852406037505716

The first line is an informational message from the support-list computation (see Miscellaneous.py). The JSON is not required: you can replace the A_d, b_d and c lines with hand-written Python objects and get the same result:

A_d = lambda delta: np.array([[-1.0, -1.0], [1.0, -1.0]])
b_d = lambda delta: np.array([[delta[0]], [-delta[0]]])   # (2, 1), not 1-D
c = np.array([[0.0], [1.0]])

Reading the result

Solution: the interval is $[0.0023,\ 0.9989]$.
Complexity: the support list consists of the smallest sample (0.0023) and the largest (0.9989), so $k = 2$.
Risk: with confidence at least $1 - 10^{-6}$, the probability that a new point drawn from the same distribution falls outside the interval is at most $\overline{\varepsilon} = 0.2085$. The lower bound is $\underline{\varepsilon} = 0$.
MOSEK: the stored MOSEK result in results/metrics.json agrees.

Numeric mode

For numeric mode, use data/program_numeric.json, which has no A_d/b_d and declares "rows_A": 2, together with data/scenarios_numeric.csv. Each row of that file is $[A_i[0,:] \mid A_i[1,:] \mid b_i]$. The first row, for $\delta_1 \approx 0.013984$, is

-1, -1, 1, -1, 0.013984..., -0.013984...

and the result is the same as in symbolic mode.

Relaxing the constraints

Select option Relaxation and enter ρ = 0.02, 0.1, 1. The tool then solves three relaxed LPs (web interface, CLARABEL, β = 10⁻⁶):

$\rho$	half-width $w$	violated points	$k$	$[\underline{\varepsilon}, \overline{\varepsilon}]$
0.02	0.2737	50	50	[0.2367, 0.7638]
0.1	0.4644	10	10	[0, 0.3432]
1	0.4983	0	2	[0, 0.2085]

A small $\rho$ makes it cheap to leave points outside, so the interval shrinks and the risk grows. The reported confidence is $1 - 3 \times 10^{-6} = 0.999997$.

Other LP benchmarks

Each benchmark's README describes the problem, its matrices and the expected results.

LP_growth_bound_12d: data-driven growth bound for a finite abstraction of 3-D vehicle dynamics, from 3127 sampled transitions. Robust LP with $d = 12$; degeneracy is detected.
LP_inventory_11d: daily ordering of 5 perishable products under uncertain yield, demand and packing, with budget and warehouse limits ($d = 11$, $N = 500$).
LP_portfolio_cvar_13d: CVaR (95%) portfolio over 8 ETFs, using historical returns from Yahoo Finance ($d = 13$, $N = 1255$).
LP_power_dispatch_120d: 24-hour dispatch of 3 thermal generators with uncertain wind and solar output ($d = 120$, 330 hard constraints).

Each of these benchmarks ships data/program_symbolic.json with data/scenarios.csv, and data/program_numeric.json with data/scenarios_numeric.csv.

Linear_Programs

Linear Programs

Formulation

The four variants

Regularization norm

Inputs

Problem data

Parameters

Notes for Python callers

Writing scenario-dependent entries

Symbolic and numeric input modes

Parameter sweeps

Program JSON

Worked example: smallest enclosing interval

Program file

In the web interface

In Python

Reading the result

Numeric mode

Relaxing the constraints

Other LP benchmarks

See also

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Using The Tool

Starting Up

Convex Optimization

Graphical Results

Getting into the Weeds

Scenario Theory

Implementation Details

Clone this wiki locally