1. Foundations of Operations Research

Operations Research (OR) is the discipline of using mathematical models and algorithms to make better decisions about complex systems. This chapter covers where the field comes from, how a real-world problem becomes a mathematical model, and the linear algebra and convexity tools that the rest of the guide relies on. Linear programming, network flows, integer programming and the later chapters all build on these three ideas: decision variables, feasible regions and optimality.

<a id="11-introduction-to-operations-research"></a>

1.1 Introduction to Operations Research

This subchapter introduces OR as a way of thinking: what problems it solves, how models are built, and where it is used today.

History and Origins of OR

Why it matters. Knowing where OR came from explains its character. It started as a practical, interdisciplinary effort to answer urgent questions with data and mathematics, and that is still its core. It was never just pure mathematics.

Theory. The term operational research was coined in Britain in the late 1930s. Scientists were asked to study how newly invented radar systems could be integrated into air defence. The question was operational ("how do we use this system?"), not technical ("how do we build it?"), and that difference gave the field its name.

During World War II, mixed teams of physicists, mathematicians, biologists and statisticians studied military operations. The most famous was the group led by physicist Patrick Blackett, nicknamed "Blackett's Circus". Typical questions:

• At what depth should anti-submarine depth charges detonate?
• How large should merchant convoys be?
• Where should aircraft be inspected to reduce losses?

The convoy question is a classic example of the OR mindset. Intuition said small convoys were safer because they were harder to spot. The data showed that losses depended mainly on the convoy's perimeter, which is what submarines attack and escorts defend, while the number of ships grows with its area. Larger convoys were therefore safer per ship.

After the war, OR spread to industry and government, helped by two developments: new algorithms and electronic computers.

Period	Milestone	Key people
1939	Mathematical methods for production planning (early linear programming)	Leonid Kantorovich
1940s	Military operational research in the UK and US	Patrick Blackett, Philip Morse
1944	Theory of games (foundation of game theory and duality ideas)	John von Neumann, Oskar Morgenstern
1947	Simplex method for linear programming	George Dantzig
1950s	Dynamic programming; network flows; queueing applications	Richard Bellman; L. R. Ford Jr., D. R. Fulkerson
1958	Cutting planes for integer programming	Ralph Gomory
1970s	NP-completeness frames which problems are "hard"	Stephen Cook, Richard Karp
1979-1984	Polynomial-time LP algorithms (ellipsoid, interior point)	Leonid Khachiyan, Narendra Karmarkar
1990s-today	Industrial MIP solvers, metaheuristics, links with machine learning	Many

 1930s        1940s           1950s              1970s            1980s            today
   |            |               |                  |                |                |
 radar      WWII teams     simplex (1947)      complexity      interior point    large-scale
 studies    convoys,       dynamic prog.,      theory          methods           MIP solvers,
            depth charges  network flows                                         ML + OR
   |____________|_______________|__________________|________________|________________|
     military origin      ->   algorithms + computers   ->   theory   ->   industrial practice

The convoy argument fits in a few lines of code:

import math

# Model a convoy as a circle of radius r. Ships fill the area,
# escorts and submarine attacks act on the perimeter.
for r in [1, 2, 4, 8]:
    ships = math.pi * r**2          # number of ships grows with the area
    perimeter = 2 * math.pi * r     # exposed boundary grows only linearly
    exposure_per_ship = perimeter / ships   # = 2 / r
    print(f"radius={r}: ships={ships:6.1f}  exposure per ship={exposure_per_ship:.3f}")

# radius=1: ships=   3.1  exposure per ship=2.000
# radius=8: ships= 201.1  exposure per ship=0.250   <- 8x safer per ship

Key takeaways

• OR began in the late 1930s and WWII as the scientific study of operations: how to use systems, not how to build them.
• Its early success came from interdisciplinary teams, data analysis and simple quantitative models, often ones that went against intuition.
• The simplex method (1947), dynamic programming, network flows and computers turned OR from wartime analysis into a general decision-making discipline.
• Modern OR combines decades of algorithmic theory with industrial-strength solvers.

Practice

1. (Beginner) List three WWII problems that early OR teams studied and, for each, name the decision being made.
2. (Intermediate) Extend the convoy model: suppose each convoy needs one escort per 10 units of perimeter. Write a function that returns escorts needed per ship as a function of radius, and explain what it implies about convoy size.
3. (Interview) "The simplex method was published in 1947, yet LP only became widely used in industry years later. Why?" Hint: think about what solving a model with hundreds of variables required at the time.

The OR Modeling Process

Why it matters. Solving the equations is usually the easiest part of an OR project. Most failures come from solving the wrong problem, using bad data, or producing an answer nobody implements. A disciplined process protects against all three.

Theory. An OR model is a deliberately simplified representation of reality, much like a map. A map is useful because it leaves out detail: a subway map that showed every street would be useless for choosing a train. Likewise, a model keeps only the elements that drive the decision. As statistician George Box put it: "All models are wrong, but some are useful."

The OR process is a loop, not a straight line:

        +-----------------------------+
        |  1. Problem definition      |  What decision? Whose objective? Which limits?
        +--------------+--------------+
                       v
        +-----------------------------+
        |  2. Data collection         |  Costs, capacities, demands, times
        +--------------+--------------+
                       v
        +-----------------------------+
        |  3. Model formulation       |  Variables, objective, constraints
        +--------------+--------------+
                       v
        +-----------------------------+
        |  4. Solution                |  Algorithm / solver / heuristic
        +--------------+--------------+
                       v
        +-----------------------------+       does not match reality?
        |  5. Validation              | ------------------------------+
        +--------------+--------------+                               |
                       v                                              |
        +-----------------------------+                               |
        |  6. Implementation          |                               |
        +--------------+--------------+                               |
                       v                                              |
        +-----------------------------+                               |
        |  7. Monitoring & update     | ---> back to step 1 or 3 <----+
        +-----------------------------+

Step	Key question	Typical pitfall
Problem definition	What exactly are we deciding, and why?	Optimizing a symptom instead of the real problem
Data collection	Are the numbers reliable and current?	Using averages that hide variability
Formulation	Does the model capture the essential trade-offs?	Too much detail (unsolvable) or too little (useless)
Solution	Which method fits the model's structure?	Using a general solver where structure could be exploited
Validation	Does the model reproduce known past outcomes?	Skipping it because "the math is correct"
Implementation	Will people and systems adopt the result?	Ignoring organizational constraints
Monitoring	Is the solution still good as conditions change?	Treating the model as a one-time deliverable

Example: a small bakery. Steps 1 to 4 on a toy problem:

1. Problem: a bakery wants to choose how many batches of bread and cake to produce daily to maximize profit.
2. Data: bread earns 3 and cake earns 5 (hundreds of EUR per batch). Bread uses 1 h of the bread oven (4 h available). Cake uses 2 h of the cake oven (12 h available). Both share a mixer: 3 h per bread batch and 2 h per cake batch (18 h available).
3. Formulation: with $x_1$ = bread batches and $x_2$ = cake batches:

$$ \begin{aligned} \max \quad & 3x_1 + 5x_2 \ \text{s.t.} \quad & x_1 \le 4 \ & 2x_2 \le 12 \ & 3x_1 + 2x_2 \le 18 \ & x_1, x_2 \ge 0 \end{aligned} $$

4. Solution:

from scipy.optimize import linprog

# Decision variables: x1 = bread batches, x2 = cake batches
profit = [3, 5]                 # profit per batch (hundreds of EUR)

# linprog MINIMIZES, so we negate the profit to maximize it
c = [-p for p in profit]

A_ub = [[1, 0],                 # bread oven:  1 h per bread batch      <= 4 h
        [0, 2],                 # cake oven:   2 h per cake batch       <= 12 h
        [3, 2]]                 # shared mixer: 3 h bread, 2 h cake     <= 18 h
b_ub = [4, 12, 18]

res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=[(0, None), (0, None)])
print("plan   :", res.x)          # [2. 6.]  -> 2 bread, 6 cake batches
print("profit :", -res.fun)       # 36.0     -> undo the sign flip

5. Validation: would the bakery manager accept 2 bread and 6 cake batches? If cakes go unsold, the hidden assumption "everything produced is sold" is wrong. Add a demand constraint and go back to step 3.

This bakery model comes back several times in this chapter.

Key takeaways

• A model is a purposeful simplification. Its value is measured by the decisions it improves, not by its realism.
• The process is iterative: validation and monitoring regularly send you back to earlier steps.
• Problem definition and data collection usually take more effort than solving.
• An unimplemented optimal solution has zero value, so plan for adoption from the start.

Practice

1. (Beginner) For the bakery, name two pieces of data that would be hard to measure accurately in practice and explain how errors in them could change the plan.
2. (Intermediate) Add a market constraint "at most 5 cake batches can be sold per day" to the code above, re-solve, and describe how the plan changes and why.
3. (Interview) "You built a model whose optimal solution looks great on paper, but operations staff keep overriding it. How do you proceed?" Hint: the problem is often in steps 1, 2 or 5, not in the solver.

Decision Problems and Decision Variables

Why it matters. Every optimization model starts by answering the question "what can we actually control?". The answer, the decision variables, determines how easy the rest of the model is to write, how large it is, and how fast it can be solved.

Theory. A decision problem asks you to choose one alternative from a set of possibilities to best achieve a goal while respecting restrictions. Every OR model has four ingredients:

Ingredient	Meaning	Bakery example
Decision variables	Quantities the decision maker controls	$x_1, x_2$: batches of bread and cake
Parameters	Given data, outside the decision maker's control	Profits 3 and 5; capacities 4, 12, 18
Objective function	Measure of quality of a decision	$3x_1 + 5x_2$
Constraints	Rules any acceptable decision must satisfy	Oven and mixer hours

Think of decision variables as the knobs on a control panel and parameters as the weather. You turn the knobs, but you must accept the weather.

Decision variables come in three basic types:

Type	Domain	Example
Continuous	$x \in \mathbb{R}$, often $x \ge 0$	Liters of fuel to blend
Integer	$x \in \mathbb{Z}$	Number of trucks to buy
Binary	$x \in {0, 1}$	Open warehouse $j$ (1) or not (0)

Choosing variables well. The same problem can often be modeled with different variables, and the choice matters a lot. Classic example: staff scheduling where each nurse works 5 consecutive days.

• Option A: $y_t$ = number of nurses on duty on day $t$. This is easy to state, but the rule "5 consecutive days" cannot be expressed with $y_t$ alone.
• Option B: $x_t$ = number of nurses who start their 5-day shift on day $t$. Then the number on duty on day $t$ is the sum of $x_s$ over the 5 start days $s$ that cover $t$, and the consecutive-days rule is built in.

Option B works because its variables match the real decision: when each shift begins.

Size of the decision space. With $n$ binary variables there are $2^n$ candidate decisions. Brute-force enumeration works for tiny problems only:

from itertools import product

# Five candidate projects: (cost, value). Budget = 10.
projects = [(4, 7), (3, 5), (5, 8), (2, 3), (6, 8)]
budget = 10

best_value, best_choice, evaluated = -1, None, 0
# Each decision variable x_j is binary: 1 = fund project j, 0 = skip it
for x in product([0, 1], repeat=len(projects)):
    evaluated += 1
    cost = sum(xj * c for xj, (c, _) in zip(x, projects))
    if cost > budget:            # constraint violated -> infeasible choice
        continue
    value = sum(xj * v for xj, (_, v) in zip(x, projects))
    if value > best_value:
        best_value, best_choice = value, x

print("evaluated", evaluated, "combinations")      # 32 = 2^5
print("best choice:", best_choice, "value:", best_value)
# best choice: (0, 1, 1, 1, 0) value: 16  -> fund projects 2, 3 and 4

With 60 projects there would be $2^{60} \approx 10^{18}$ combinations. This explosion is why OR develops smarter methods instead of enumerating.

Key takeaways

• A decision problem has four parts: variables (controllable), parameters (given), an objective, and constraints.
• Variables can be continuous, integer or binary, and the type largely decides which algorithms apply.
• Good variable choices make constraints simple. If a constraint is awkward to write, try different variables.
• The number of candidate solutions grows exponentially with the number of discrete variables, so enumeration does not scale.

Practice

1. (Beginner) A delivery company decides how many vans to lease (vans cost 500/month) and which of 8 depots to keep open. Identify the decision variables, their types, and two parameters.
2. (Intermediate) Write out the nurse-scheduling constraint for Wednesday using Option B variables ($x_{\text{Mon}}, \dots, x_{\text{Sun}}$) if at least 12 nurses must be on duty and shifts wrap around the week.
3. (Interview) "How would you choose decision variables to model which route a single delivery truck takes through 20 customers?" Hint: consider a binary variable for each ordered pair of locations, and ask what it means for that variable to equal 1.

Deterministic vs Stochastic Models

Why it matters. Real data is uncertain: demand, travel times, prices and failures all vary. Treating uncertain quantities as fixed numbers is convenient, but it can lead to plans that look optimal on paper and perform worse in practice.

Theory. A deterministic model assumes every parameter is known exactly. A stochastic model treats some parameters as random variables with known (or estimated) probability distributions.

Aspect	Deterministic	Stochastic
Parameters	Fixed numbers	Random variables / scenarios
Output	A single optimal decision and value	A decision with a distribution of outcomes
Typical models	LP, IP, network flows, deterministic DP	Queueing, Markov chains, stochastic programming, simulation
Computation	Usually easier	Usually harder (expectations, many scenarios)
When appropriate	Low variability, short horizon, planning	High variability, risk matters

Decision theory also distinguishes between:

• Risk: outcomes are uncertain but their probabilities are known (dice, well-measured demand).
• Uncertainty (ambiguity): even the probabilities are unknown (a new product launch). This motivates robust and worst-case approaches (Chapters 7 and 9).

The flaw of averages. Replacing a random parameter with its mean and solving a deterministic model has a hidden cost. When the objective is nonlinear in the random parameter, as $\min(Q, D)$ is in demand $D$, the profit "at average demand" is not the average profit.

Example: newsvendor. A kiosk buys newspapers at 2 and sells them at 10. Unsold papers are worthless, and daily demand is uniform between 50 and 150 (mean 100).

import numpy as np

rng = np.random.default_rng(seed=42)
price, cost = 10, 2                  # sell at 10, buy at 2, leftovers are worthless
demand = rng.uniform(50, 150, size=100_000)   # uncertain demand scenarios

def expected_profit(q):
    sold = np.minimum(q, demand)     # cannot sell more than demand or stock
    return np.mean(price * sold - cost * q)

# Deterministic model: replace demand by its mean (100) and order exactly that
q_det = 100
print("deterministic promise :", price * q_det - cost * q_det)   # 800
print("deterministic reality :", round(expected_profit(q_det), 1)) # ~675

# Stochastic model: search the order quantity that maximizes expected profit
qs = np.arange(50, 151)
q_sto = qs[np.argmax([expected_profit(q) for q in qs])]
print("stochastic order      :", q_sto)                           # 130
print("stochastic reality    :", round(expected_profit(q_sto), 1)) # ~720

The deterministic model promises 800, but on average its plan earns about 675. The stochastic model orders more (130) because a lost sale costs 8 while an unsold paper costs only 2, and it earns about 720 on average. Chapter 10 derives this "order more than the mean" result in closed form.

Key takeaways

• Deterministic models assume perfect information. Stochastic models represent uncertain parameters with probability distributions.
• Plugging averages into a model can overestimate performance and recommend the wrong decision (the flaw of averages).
• Stochastic models are more realistic but more expensive to solve. Use them when variability is large and its consequences are asymmetric.
• "Risk" means known probabilities. "Uncertainty" means the probabilities themselves are unknown.

Practice

1. (Beginner) Classify each as better suited to a deterministic or stochastic model, and justify: (a) cutting steel bars to fixed lengths from known orders, (b) staffing a call center, (c) scheduling next year's university exams.
2. (Intermediate) Modify the newsvendor code so unsold papers can be returned for 1 each. Predict whether the optimal order goes up or down, then verify.
3. (Interview) "Our planning team uses forecasted average demand in an optimization model, and actual costs are consistently higher than planned. Why might that be?" Hint: compare $f(\mathbb{E}[D])$ with $\mathbb{E}[f(D)]$ for a nonlinear $f$.

Applications in Industry and Computing

Why it matters. OR is everywhere, usually out of sight: in the price of an airline ticket, the route of a parcel, the placement of a virtual machine, the order in which a database joins tables. Recognizing the OR model behind a practical problem is the first step toward solving it well.

Theory. Most applications reduce to a small set of standard model families. Learn the families and you can recognize them in very different settings.

Domain	Problem	Underlying OR model	Chapter
Logistics	Delivery routes for a fleet	Vehicle routing (integer programming)	4
Manufacturing	What to produce each week	Production planning (LP / MIP)	2, 4
Airlines	Assign crews to flights	Set covering / partitioning	4
Airlines, hotels	Seat or room pricing over time	Revenue management (DP, stochastic models)	6, 9
Energy	Which power plants to switch on each hour	Unit commitment (MIP)	4
Finance	Portfolio allocation	Quadratic programming	5
Healthcare	Operating room scheduling	Scheduling / IP	10
Telecom	Network capacity design	Multicommodity flows	3
Computer networks	Packet routing (e.g., OSPF)	Shortest paths (Dijkstra)	3
Compilers	Register allocation	Graph coloring	4
Databases	Join order in query planning	Dynamic programming	6
Cloud computing	Place VMs/containers on servers	Bin packing	4, 11
Operating systems	CPU and I/O scheduling	Queueing theory, scheduling	9, 10
Machine learning	Train a model, fit an SVM	Nonlinear / convex optimization	5

Example: VM placement as bin packing. A cloud provider must place VMs (each requesting some vCPUs) on hosts with 16 vCPUs each, using as few hosts as possible. This is the bin packing problem, which is NP-hard (Chapter 4). A fast greedy heuristic, First Fit Decreasing, usually performs very well:

def first_fit_decreasing(vm_sizes, capacity):
    """Place VMs (vCPU demand) on as few hosts as possible (heuristic)."""
    hosts = []                                   # each host = list of VM sizes
    for size in sorted(vm_sizes, reverse=True):  # largest VMs first
        for host in hosts:
            if sum(host) + size <= capacity:     # first host with room wins
                host.append(size)
                break
        else:                                    # no host fits -> open a new one
            hosts.append([size])
    return hosts

vms = [8, 11, 8, 3, 6, 3, 8, 2, 10]              # vCPUs requested by each VM
hosts = first_fit_decreasing(vms, capacity=16)
for i, h in enumerate(hosts, 1):
    print(f"host {i}: {h} -> {sum(h)}/16 vCPU")
# host 1: [11, 3, 2] -> 16/16 vCPU
# host 2: [10, 6] -> 16/16 vCPU
# host 3: [8, 8] -> 16/16 vCPU
# host 4: [8, 3] -> 11/16 vCPU

# Lower bound: total demand / capacity, rounded up
print("lower bound:", -(-sum(vms) // 16), "hosts")   # 4 -> the heuristic is optimal here

The lower bound of $\lceil 59 / 16 \rceil = 4$ hosts proves that this solution is optimal. Comparing a heuristic against a bound is a pattern used throughout OR.

Key takeaways

• A handful of model families (LP, IP, network flows, DP, queueing, scheduling) cover most applications.
• Computing is full of OR problems: routing, compilation, query planning, scheduling and resource allocation in the cloud.
• Recognizing the underlying model lets you reuse known algorithms, complexity results and solvers.
• Bounds tell you how good a heuristic solution is, even when the true optimum is unknown.

Practice

1. (Beginner) For each of these, name the OR model family: (a) finding the fastest route in a navigation app, (b) choosing which 3 of 10 cities get a new warehouse, (c) estimating average wait time at a supermarket checkout.
2. (Intermediate) Implement Best Fit Decreasing (place each VM on the host with the least remaining room that still fits) and compare it with First Fit Decreasing on random instances of 50 VMs.
3. (Interview) "How would you formulate placing containers on servers so that the number of active servers is minimized, as an integer program?" Hint: use a binary $y_j$ for "server $j$ is used" and a binary $x_{ij}$ for "container $i$ runs on server $j$", then link them through the capacity constraint.

<a id="12-mathematical-prerequisites"></a>

1.2 Mathematical Prerequisites

This subchapter reviews the linear algebra and convex analysis that OR algorithms are built on, always pointing out why each concept matters for optimization.

Vectors and Matrices

Why it matters. Real models have thousands of variables and constraints. Vectors and matrices let us write them compactly, as $\max, c^\top x$ subject to $Ax \le b$, and manipulate them efficiently in code.

Theory. A vector $x \in \mathbb{R}^n$ is an ordered list of $n$ numbers. In OR, a vector usually describes a plan: $x = (x_1, \dots, x_n)$, where $x_j$ is the amount of activity $j$. A matrix $A \in \mathbb{R}^{m \times n}$ is a table of numbers with $m$ rows and $n$ columns. In an LP, $a_{ij}$ is the amount of resource $i$ consumed by one unit of activity $j$.

Core operations:

Operation	Notation	Meaning in OR
Dot product	$c^\top x = \sum_j c_j x_j$	Total profit/cost of plan $x$
Matrix-vector product	$Ax$	Resource usage: entry $i$ is $\sum_j a_{ij} x_j$
Transpose	$A^\top$	Swaps rows and columns (used in duality, Chapter 2)
Identity	$I$	$Ix = x$; appears when slack variables are added
Inverse	$A^{-1}$ (square, nonsingular)	Solves $Ax = b$ as $x = A^{-1}b$
Componentwise inequality	$x \le y$ iff $x_i \le y_i$ for all $i$	$Ax \le b$ packs $m$ constraints into one line

Two ways to read $A$:

                     product 1  product 2          <- columns = activities
                   +----------+----------+
   bread oven      |    1     |    0     |   <= 4
   cake oven       |    0     |    2     |   <= 12        rows = resources
   mixer           |    3     |    2     |   <= 18
                   +----------+----------+
                        ^
                        |
           column j = "recipe" of activity j: what ONE unit of it consumes
           row i    = constraint i: how all activities together use resource i

import numpy as np

A = np.array([[1, 0],      # row i = how much of resource i each product uses
              [0, 2],
              [3, 2]])
b = np.array([4, 12, 18])  # resource availability
c = np.array([3, 5])       # profit per unit of each product
x = np.array([2, 6])       # a candidate production plan

usage = A @ x              # matrix-vector product: resource consumption
print("usage      :", usage)                      # [ 2 12 18]
print("feasible?  :", np.all(usage <= b) and np.all(x >= 0))   # True
print("leftover   :", b - usage)                  # [2 0 0] -> slack per resource
print("profit c.x :", c @ x)                      # 36 -> dot product = objective
print("column 2   :", A[:, 1])                    # [0 2 2] -> recipe of one cake batch

Key takeaways

• Vectors represent plans. Matrices represent how activities consume resources.
• $Ax \le b$ is $m$ inequalities written as one; $c^\top x$ is the objective value.
• Rows of $A$ are constraints, columns are activities. Both views are used constantly (the column view in the simplex method, the row view in duality).
• NumPy's @ operator performs both matrix-vector and dot products.

Practice

1. (Beginner) For the bakery data, compute by hand $Ax$ for $x = (4, 3)$ and decide whether the plan is feasible.
2. (Intermediate) Write a NumPy function check_plan(A, b, x) that returns the list of violated constraints and by how much each is violated.
3. (Interview) "In the LP $\max, c^\top x,\ Ax \le b$, what do a row and a column of $A$ represent, and what would adding a new column mean for the business?" Hint: a new column is a new activity, such as a new product.

Systems of Linear Equations

Why it matters. The corners of an LP's feasible region, where optimal solutions live, are exactly the solutions of systems of linear equations formed by "active" constraints. Solving linear systems is the computational heart of the simplex method.

Theory. A linear system $Ax = b$ with $m$ equations and $n$ unknowns asks for all $x$ satisfying every equation at once. Geometrically, each equation in $\mathbb{R}^2$ is a line (in $\mathbb{R}^n$, a hyperplane), and a solution is a point lying on all of them.

There are exactly three possible outcomes:

  Unique solution         No solution              Infinitely many
  (lines cross)           (parallel lines)         (same line)

      \   /                 ------                   ======
       \ /                                           (two equations
        X  <- solution      ------                   describe one line)
       / \

Case	Algebraic signature	Example
Unique	Square, $\det A \neq 0$	$x + y = 10,\ x - y = 2$ gives $(6, 4)$
None	Contradictory equations	$x + y = 10,\ 2x + 2y = 7$
Infinite	Fewer independent equations than unknowns	$x + y = 10$ alone

The OR case: underdetermined systems. When inequality constraints are turned into equations by adding slack variables (Chapter 2), we typically get $m$ equations with $n > m$ unknowns, which has infinitely many solutions. The simplex method picks out special ones: set $n - m$ variables to zero and solve the remaining $m \times m$ system. These are the basic solutions, and they correspond to corners of the feasible region.

import numpy as np

# Square system with a unique solution
A = np.array([[1.0, 1.0],
              [1.0, -1.0]])
b = np.array([10.0, 2.0])
print("unique solution:", np.linalg.solve(A, b))      # [6. 4.]

# Parallel lines: x + y = 10 and 2x + 2y = 7 -> singular matrix
A_bad = np.array([[1.0, 1.0],
                  [2.0, 2.0]])
try:
    np.linalg.solve(A_bad, np.array([10.0, 7.0]))
except np.linalg.LinAlgError as e:
    print("no unique solution:", e)                   # Singular matrix

# Underdetermined: 3x1 + 2x2 + s = 18 (one equation, three unknowns).
# Fix n - m = 2 variables to zero and solve for the remaining one:
row = np.array([3.0, 2.0, 1.0])
for keep in range(3):
    x = np.zeros(3)
    x[keep] = 18.0 / row[keep]          # solve the 1x1 system
    print("basic variable", ["x1", "x2", "s"][keep], "->", x)
# basic variable x1 -> [6. 0. 0.]
# basic variable x2 -> [0. 9. 0.]
# basic variable s  -> [ 0.  0. 18.]

Key takeaways

• A linear system has a unique solution, no solution, or infinitely many solutions. There is no other possibility.
• Each equation is a hyperplane, and solving the system means intersecting them.
• LPs in equality form are underdetermined ($n > m$). Basic solutions come from fixing $n - m$ variables to zero.
• In code, prefer np.linalg.solve to computing an explicit inverse. It is faster and numerically safer.

Practice

1. (Beginner) Solve $2x + y = 7,\ x - y = 2$ by hand, then verify with NumPy.
2. (Intermediate) For the system $x_1 + x_2 + s_1 = 4,\ x_1 + 3x_2 + s_2 = 6$ (2 equations, 4 unknowns), list all basic solutions obtained by setting two variables to zero. Which of them have all components $\ge 0$?
3. (Interview) "Why is computing $A^{-1}b$ explicitly a bad idea for solving $Ax = b$ in production code?" Hint: think about operation count and how rounding errors are amplified.

Gaussian Elimination

Why it matters. Gaussian elimination is the standard algorithm for solving linear systems, and a pivot in the simplex method is exactly a Gauss-Jordan elimination step. Understanding it by hand makes the simplex tableau (Chapter 2) straightforward.

Theory. The idea is to transform the system into an equivalent one that is trivial to solve, much like tidying a messy room one shelf at a time so that each shelf, once done, is never touched again. Three elementary row operations never change the solution set:

1. Swap two rows.
2. Multiply a row by a nonzero constant.
3. Add a multiple of one row to another row.

The algorithm has two phases:

• Forward elimination: use row operations to create zeros below the diagonal, reaching an upper-triangular (row echelon) form.
• Back substitution: solve the last equation first, then substitute upward.

Gauss-Jordan elimination continues until the matrix becomes the identity (reduced row echelon form), so the solution can be read off directly.

Worked example. Solve

$$ \begin{aligned} 2x + y - z &= 8 \ -3x - y + 2z &= -11 \ -2x + y + 2z &= -3 \end{aligned} $$

Augmented matrix [A | b]:

  [  2   1  -1 |   8 ]
  [ -3  -1   2 | -11 ]
  [ -2   1   2 |  -3 ]

R2 <- R2 + (3/2) R1          R3 <- R3 + R1
  [  2   1    -1   |  8 ]
  [  0  1/2  1/2   |  1 ]
  [  0   2    1    |  5 ]

R3 <- R3 - 4 R2
  [  2   1    -1   |  8 ]
  [  0  1/2  1/2   |  1 ]     <- upper triangular: forward phase done
  [  0   0    -1   |  1 ]

Back substitution:
  -z = 1                 ->  z = -1
  (1/2)y + (1/2)z = 1    ->  y = 3
  2x + y - z = 8         ->  x = 2

Partial pivoting. If a pivot is zero, elimination fails. If it is tiny, rounding errors explode. Partial pivoting swaps in the row with the largest absolute value in the current column before eliminating. Every serious implementation does this.

import numpy as np

def gauss_solve(A, b):
    """Solve Ax = b with Gaussian elimination and partial pivoting."""
    A = A.astype(float).copy()
    b = b.astype(float).copy()
    n = len(b)

    # Forward elimination: create zeros below the diagonal
    for k in range(n - 1):
        p = k + np.argmax(np.abs(A[k:, k]))   # row with largest pivot candidate
        if np.isclose(A[p, k], 0.0):
            raise ValueError("matrix is singular")
        A[[k, p]], b[[k, p]] = A[[p, k]], b[[p, k]]   # swap rows k and p
        for i in range(k + 1, n):
            m = A[i, k] / A[k, k]             # multiplier
            A[i, k:] -= m * A[k, k:]          # R_i <- R_i - m * R_k
            b[i] -= m * b[k]

    # Back substitution: solve from the last equation upward
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (b[i] - A[i, i + 1:] @ x[i + 1:]) / A[i, i]
    return x

A = np.array([[ 2,  1, -1],
              [-3, -1,  2],
              [-2,  1,  2]])
b = np.array([8, -11, -3])
print(gauss_solve(A, b))          # [ 2.  3. -1.]

Property	Value
Cost for an $n \times n$ system	About $\tfrac{2}{3}n^3$ floating-point operations, i.e. $O(n^3)$
Back substitution	$O(n^2)$
Stability	Good in practice with partial pivoting
Relation to LU factorization	Forward elimination computes $PA = LU$, reusable for many right-hand sides

Key takeaways

• Elementary row operations preserve the solution set. Gaussian elimination applies them systematically.
• Forward elimination produces a triangular system, and back substitution solves it.
• Partial pivoting is essential for numerical stability.
• A simplex pivot is a Gauss-Jordan step: it makes one column a unit vector.

Practice

1. (Beginner) Solve $x + 2y = 5,\ 3x + 4y = 11$ using row operations on the augmented matrix. Show each step.
2. (Intermediate) Extend gauss_solve into a gauss_jordan function that reduces $[A \mid b]$ to $[I \mid x]$, and print the matrix after each pivot.
3. (Interview) "What goes wrong if you run Gaussian elimination without pivoting on $\begin{bmatrix} 10^{-20} & 1 \ 1 & 1 \end{bmatrix} x = \begin{bmatrix} 1 \ 2 \end{bmatrix}$ in double precision?" Hint: compute the multiplier and look at what happens to $1 - 10^{20}$ in floating point.

Rank and Linear Independence

Why it matters. Rank tells you whether constraints are redundant and whether a system can be solved. In LP theory, a basis is a set of linearly independent columns of $A$, so independence is what makes basic solutions well-defined.

Theory. A linear combination of vectors $v_1, \dots, v_k$ is $\alpha_1 v_1 + \dots + \alpha_k v_k$. The vectors are linearly independent if the only combination giving the zero vector is the trivial one ($\alpha_1 = \dots = \alpha_k = 0$). Equivalently, no vector is a combination of the others: each brings genuinely new information.

Analogy: three recipes are dependent if one is just "2 portions of recipe A plus 3 portions of recipe B". It adds no new possibility.

• A basis of a space is a set of independent vectors that spans it.
• The rank of $A$ is the maximum number of independent columns, which always equals the maximum number of independent rows.
• Rank-nullity: for $A \in \mathbb{R}^{m \times n}$, $\operatorname{rank}(A) + \dim(\text{null space}) = n$.

Rouche-Capelli theorem (solvability):

Condition	Result
$\operatorname{rank}(A) \neq \operatorname{rank}([A \mid b])$	No solution
$\operatorname{rank}(A) = \operatorname{rank}([A \mid b]) = n$	Unique solution
$\operatorname{rank}(A) = \operatorname{rank}([A \mid b]) < n$	Infinitely many solutions ($n - \text{rank}$ free parameters)

In LP. We usually assume $A \in \mathbb{R}^{m \times n}$ has full row rank ($\operatorname{rank} A = m$), meaning no redundant equality constraints. A basis is then a choice of $m$ linearly independent columns, which gives an invertible $m \times m$ matrix $B$. There are at most $\binom{n}{m}$ such choices, which is finite but grows fast.

import numpy as np
from math import comb

v1 = np.array([1, 0, 2])
v2 = np.array([0, 1, 1])
v3 = np.array([2, 3, 7])           # v3 = 2*v1 + 3*v2  -> dependent

M = np.column_stack([v1, v2, v3])
print("rank of [v1 v2 v3]:", np.linalg.matrix_rank(M))   # 2, not 3

# Rouche-Capelli: compare rank(A) with rank([A | b])
A = M
for b in (np.array([3, 4, 10]), np.array([1, 1, 1])):
    rA  = np.linalg.matrix_rank(A)
    rAb = np.linalg.matrix_rank(np.column_stack([A, b]))
    status = "solvable" if rA == rAb else "inconsistent"
    print(f"b={b}: rank(A)={rA}, rank([A|b])={rAb} -> {status}")
# b=[ 3  4 10]: rank(A)=2, rank([A|b])=2 -> solvable      (b = 3*v1 + 4*v2)
# b=[1 1 1]:    rank(A)=2, rank([A|b])=3 -> inconsistent

# How many candidate bases does an LP with m=3 rows, n=6 columns have?
print("C(6, 3) =", comb(6, 3))     # 20

Key takeaways

• Vectors are independent when none can be built from the others.
• Rank is the number of independent rows (equivalently, columns). It measures how much non-redundant information $A$ contains.
• $Ax = b$ is solvable iff $\operatorname{rank}(A) = \operatorname{rank}([A \mid b])$, and the solution is unique iff that rank also equals $n$.
• An LP basis is $m$ independent columns of $A$. The number of bases, and hence of basic solutions, is finite.

Practice

1. (Beginner) Are $(1, 2)$ and $(2, 4)$ linearly independent? What about $(1, 2)$ and $(2, 3)$? Justify with a determinant.
2. (Intermediate) Build a $3 \times 5$ matrix with rank 2 and write code that finds which constraint (row) is redundant by removing rows one at a time and checking the rank.
3. (Interview) "Why do LP textbooks assume the constraint matrix has full row rank, and what can a solver do if it does not?" Hint: consider what a dependent equality row means: it is either redundant or it makes the system inconsistent.

Convex Sets

Why it matters. Convexity separates "easy" optimization problems from "hard" ones. If the feasible region is convex (and the objective is suitable), any local optimum is a global optimum. Every LP feasible region is convex.

Theory. A set $C$ is convex if, for any two points $x, y \in C$, the entire line segment between them also lies in $C$:

$$ x, y \in C,\ \lambda \in [0, 1] \implies \lambda x + (1 - \lambda) y \in C $$

Analogy: a room is convex if, standing anywhere inside, you can see every other point without looking around a corner.

     Convex                        Not convex

    .--------.                    .--------.
   /          \                  /    __    \
  |  x------y  |                |  x-/  \-y  |    <- segment leaves the set
   \          /                  \  |    |  /
    '--------'                    '-'    '-'

Important convex sets:

Set	Definition	Convex?
Hyperplane	${x : a^\top x = \beta}$	Yes
Halfspace	${x : a^\top x \le \beta}$	Yes
Ball	${x : \lVert x - x_0 \rVert \le r}$	Yes
Intersection of convex sets	$C_1 \cap C_2$	Always yes
Union of convex sets	$C_1 \cup C_2$	Not in general
Integer points	${x \in \mathbb{Z}^n : Ax \le b}$	No

The intersection rule is the key fact for OR. An LP feasible region is an intersection of halfspaces, so it is always convex. Integer constraints destroy convexity, which is the root of why integer programming is harder (Chapter 4).

A convex combination of points $x_1, \dots, x_k$ is $\sum_i \lambda_i x_i$ with $\lambda_i \ge 0$ and $\sum_i \lambda_i = 1$, a weighted average. The convex hull of a set is the set of all convex combinations of its points: the shape a rubber band would take if stretched around them.

import numpy as np

rng = np.random.default_rng(0)

def in_disk(p):        # convex: unit disk
    return p[0]**2 + p[1]**2 <= 1

def in_ring(p):        # NOT convex: annulus 0.5 <= r <= 1
    r2 = p[0]**2 + p[1]**2
    return 0.25 <= r2 <= 1

def looks_convex(member, trials=20_000):
    """Heuristic test: sample point pairs in the set, check their segment."""
    pts = [p for p in rng.uniform(-1, 1, size=(4000, 2)) if member(p)]
    for _ in range(trials):
        x, y = pts[rng.integers(len(pts))], pts[rng.integers(len(pts))]
        lam = rng.uniform()
        if not member(lam * x + (1 - lam) * y):   # convex combination escapes
            return False, (x, y, lam)
    return True, None

print("disk convex?", looks_convex(in_disk)[0])   # True
print("ring convex?", looks_convex(in_ring)[0])   # False

Sampling can only disprove convexity. Proving it requires the definition, or rules such as "an intersection of convex sets is convex".

Key takeaways

• Convex means the segment between any two points of the set stays in the set.
• Halfspaces, hyperplanes and balls are convex, and any intersection of convex sets is convex.
• Every LP feasible region is convex. Integer restrictions make regions non-convex.
• Convex combinations are weighted averages, and the convex hull is the smallest convex set containing given points.

Practice

1. (Beginner) Which are convex: (a) ${(x, y) : x + y \le 3,\ x \ge 0}$, (b) ${(x, y) : xy \ge 1}$, (c) ${(x, y) : xy \ge 1,\ x > 0}$? Sketch each.
2. (Intermediate) Prove from the definition that the intersection of two convex sets is convex.
3. (Interview) "Why does adding the constraint $x \in {0, 1}$ make an otherwise easy optimization problem potentially hard?" Hint: is ${0, 1}$ convex? What does that do to the "local optimum = global optimum" guarantee?

Convex and Concave Functions

Why it matters. Convex sets describe feasible regions, and convex functions describe "well-behaved" objectives and constraints. Minimizing a convex function over a convex set is the broadest class of problems that can be solved reliably at scale (Chapter 5).

Theory. A function $f$ is convex if the chord between any two points on its graph lies on or above the graph:

$$ f(\lambda x + (1 - \lambda) y) \le \lambda f(x) + (1 - \lambda) f(y), \quad \lambda \in [0, 1] $$

A function is concave if $-f$ is convex: the chord lies below the graph. Linear (affine) functions are both convex and concave, which is why LP sits comfortably inside convex optimization.

   Convex f: chord above graph            Concave f: chord below graph

   f |  *                *                f |        .-''''-.
     |   \      chord   /                   |      *----------*  <- chord
     |    *------------*                    |    /              \
     |     \          /                     |   /                \
     |      '-.____.-'                      |  /                  \
     +---------------------- x              +----------------------- x

Ways to recognize convexity:

Test	Applies when	Convex if
Definition (chord)	Always	Inequality above holds
Second derivative	$f: \mathbb{R} \to \mathbb{R}$, twice differentiable	$f''(x) \ge 0$ everywhere
Hessian	$f: \mathbb{R}^n \to \mathbb{R}$, twice differentiable	$\nabla^2 f(x)$ positive semidefinite (all eigenvalues $\ge 0$)
Composition rules	Built from known pieces	Nonnegative sums, pointwise max, composition with affine maps preserve convexity

Convex examples	Concave examples	Neither
$x^2$, $e^x$, $\lvert x \rvert$, $-\log x$, any norm $\lVert x \rVert$	$\log x$, $\sqrt{x}$, $-x^2$, $\min(x_1, x_2)$	$x^3$, $\sin x$, $x_1 x_2$

Why convexity matters for constraints. If $g$ is convex, the sublevel set ${x : g(x) \le 0}$ is a convex set. So a problem "minimize convex $f$ subject to convex $g_i(x) \le 0$ and affine equalities" has a convex feasible region and a convex objective, and every local minimum is global.

import numpy as np

def classify(H):
    """Classify a quadratic form by the eigenvalues of its (constant) Hessian."""
    eig = np.linalg.eigvalsh(H)              # symmetric -> real eigenvalues
    if np.all(eig >= 0):
        return eig, "convex (PSD)"
    if np.all(eig <= 0):
        return eig, "concave (NSD)"
    return eig, "neither (indefinite)"

# f(x, y) = x^2 + x*y + y^2  -> Hessian [[2, 1], [1, 2]]
print(classify(np.array([[2.0, 1.0], [1.0, 2.0]])))    # eigenvalues [1, 3] -> convex
# g(x, y) = x^2 - y^2        -> Hessian [[2, 0], [0, -2]]
print(classify(np.array([[2.0, 0.0], [0.0, -2.0]])))   # [-2, 2] -> indefinite (saddle)

# Check the definition for f(t) = exp(t)
x, y, lam = -1.0, 2.0, 0.3
lhs = np.exp(lam * x + (1 - lam) * y)          # f at the mixed point
rhs = lam * np.exp(x) + (1 - lam) * np.exp(y)  # mix of the f values (chord)
print(f"f(mix) = {lhs:.3f} <= chord = {rhs:.3f}: {lhs <= rhs}")   # 3.004 <= 5.283

Key takeaways

• Convex: chord above the graph. Concave: chord below. Linear: both.
• For twice-differentiable functions, check $f'' \ge 0$ (1D) or a positive semidefinite Hessian ($n$-D).
• Sums with nonnegative weights, maxima, and affine compositions of convex functions are convex.
• Minimizing a convex function (or maximizing a concave one) over a convex set has no spurious local optima.

Practice

1. (Beginner) Use the second derivative to classify $f(x) = x^4$, $g(x) = \log x$ (for $x > 0$) and $h(x) = x^3$.
2. (Intermediate) Show that $f(x) = \max(2x + 1,\ -x + 4,\ 3)$ is convex, and explain why piecewise-linear convex costs can be modeled in an LP.
3. (Interview) "A colleague says: 'My loss function is non-convex, so gradient descent is useless.' How do you respond?" Hint: distinguish between finding a good solution and certifying a global optimum.

Polyhedra and Polytopes

Why it matters. The feasible region of every linear program is a polyhedron, and if an LP has an optimal solution (and its region has at least one corner), one of the corners is optimal. That single fact turns a search over infinitely many points into a search over finitely many corners, which is the idea behind the simplex method.

Theory.

• A polyhedron is the intersection of finitely many halfspaces: $P = {x \in \mathbb{R}^n : Ax \le b}$. It is convex and may be unbounded.
• A polytope is a bounded polyhedron. By the Minkowski-Weyl theorem, a polytope is also the convex hull of finitely many points (its vertices), so it can be described either by inequalities or by corners.
• An extreme point (vertex) of $P$ is a point that is not a convex combination of two other distinct points of $P$. Algebraically, $x \in P$ is a vertex iff it makes $n$ linearly independent constraints active (hold with equality).
• Faces: vertices (dimension 0), edges (dimension 1), ..., facets (dimension $n-1$).

 Bakery feasible region (polytope):

  x2
   |
 6 *---------*  (2,6)          constraints:
   |##########\                  x1 <= 4
   |###########\                 2 x2 <= 12   (i.e. x2 <= 6)
   |############\                3 x1 + 2 x2 <= 18
 3 +#############* (4,3)         x1, x2 >= 0
   |#############|
   |#############|
   *-------------*-----------> x1
 (0,0)         (4,0)          top-left corner: (0,6)

 5 vertices: (0,0), (4,0), (4,3), (2,6), (0,6)

Brute-force vertex enumeration. In $\mathbb{R}^2$, every vertex is the intersection of 2 active constraints. Try every pair, solve the $2 \times 2$ system, keep the feasible points, and evaluate the objective:

import numpy as np
from itertools import combinations

# Bakery polyhedron written uniformly as A x <= b (nonnegativity as -x <= 0)
A = np.array([[ 1,  0],    # x1 <= 4
              [ 0,  2],    # 2 x2 <= 12
              [ 3,  2],    # 3 x1 + 2 x2 <= 18
              [-1,  0],    # x1 >= 0
              [ 0, -1]],   # x2 >= 0
             dtype=float)
b = np.array([4, 12, 18, 0, 0], dtype=float)
c = np.array([3, 5])

vertices = []
for rows in combinations(range(len(b)), 2):       # n = 2 active constraints
    sub_A, sub_b = A[list(rows)], b[list(rows)]
    if abs(np.linalg.det(sub_A)) < 1e-12:          # parallel lines: no vertex
        continue
    x = np.linalg.solve(sub_A, sub_b) + 0.0        # intersection (+0.0 cleans -0.0)
    if np.all(A @ x <= b + 1e-9):                  # keep it only if feasible
        vertices.append(x)

for v in vertices:
    print(f"vertex {v} -> objective {c @ v:.0f}")
# vertex [4. 3.] -> objective 27
# vertex [4. 0.] -> objective 12
# vertex [2. 6.] -> objective 36   <- best corner = the LP optimum found earlier
# vertex [0. 6.] -> objective 30
# vertex [0. 0.] -> objective 0

With $m$ constraints in $\mathbb{R}^n$ this approach checks up to $\binom{m}{n}$ systems, which becomes impossible for real models. The simplex method instead walks from vertex to adjacent vertex, always improving the objective.

Concept	Inequality view	Point view
Polyhedron	${x : Ax \le b}$	Convex hull of points + cone of directions
Polytope	Bounded ${x : Ax \le b}$	Convex hull of its vertices
Vertex	$n$ independent active constraints	Not a mix of two other points

Key takeaways

• A polyhedron is a finite intersection of halfspaces; a polytope is a bounded polyhedron.
• Polytopes have two equivalent descriptions: inequalities, or the convex hull of their vertices.
• A vertex is where $n$ linearly independent constraints are active.
• A linear objective over a polytope is always optimized at a vertex. This is the foundation of the simplex method.

Practice

1. (Beginner) Sketch ${x \ge 0 : x_1 + x_2 \le 5,\ x_1 - x_2 \le 1}$ and list its vertices.
2. (Intermediate) Modify the enumeration code to also print which constraints are active at each vertex, then add the constraint $x_1 + x_2 \le 7$ and observe which vertices change.
3. (Interview) "Why is it enough for an LP solver to look only at vertices, and when might an LP have optimal solutions that are not vertices?" Hint: consider an objective parallel to an edge of the polytope.

<a id="13-optimization-problem-formulation"></a>

1.3 Optimization Problem Formulation

This subchapter formalizes what an optimization problem is, introduces the vocabulary used in the rest of the guide, and explains the assumptions that every model makes.

Every problem in this guide can be written in the general form:

$$ \begin{aligned} \min_{x} \quad & f(x) \ \text{s.t.} \quad & g_i(x) \le 0, \quad i = 1, \dots, m \ & h_k(x) = 0, \quad k = 1, \dots, p \ & x \in X \quad (\text{e.g. } \mathbb{R}^n,\ \mathbb{Z}^n,\ {0,1}^n) \end{aligned} $$

Objective Functions

Why it matters. The objective function encodes what "better" means. Two models with identical constraints but different objectives can recommend very different decisions, so choosing the objective is a modeling decision with real consequences.

Theory. The objective function $f: \mathbb{R}^n \to \mathbb{R}$ assigns a single number to every decision $x$. We either minimize it (cost, time, distance, risk) or maximize it (profit, throughput, reliability). The two are interchangeable:

$$ \max_x f(x) = -\min_x \big(-f(x)\big) $$

The optimal point $x^*$ is the same; only the sign of the optimal value flips. Most solver APIs minimize by default, which is why the bakery example negates its profits.

Common objective shapes:

Shape	Form	Example	Notes
Linear	$c^\top x$	Total profit	LP / IP
Quadratic	$x^\top Q x + c^\top x$	Portfolio variance	QP (convex if $Q \succeq 0$)
Piecewise linear convex	$\max_i (a_i^\top x + b_i)$	Overtime cost tiers	Linearizable
Minimax	$\max_i f_i(x)$	Worst-case delivery time, makespan	Linearizable with an extra variable
Sum of absolute values	$\sum_i \lvert a_i^\top x - b_i \rvert$	Robust regression	Linearizable
Multiple objectives	$(f_1(x), \dots, f_k(x))$	Cost vs emissions	Chapter 7

Linearization trick for minimax. $\min_x \max_i f_i(x)$ is not linear, but it is equivalent to

$$ \min_{x, t}\ t \quad \text{s.t.}\quad t \ge f_i(x)\ \ \forall i $$

At the optimum, $t$ is pushed down onto the largest $f_i$. Similarly, $\lvert u \rvert$ can be replaced by $t$ with $t \ge u$ and $t \ge -u$.

Example: the objective changes the answer. Choose a depot location $x$ on a road serving customers at positions 1, 4 and 9. Compare "minimize the worst distance" with "minimize the total distance":

from scipy.optimize import linprog

a = [1, 4, 9]          # positions of three customers on a road
n = len(a)

# (1) Minimize the MAXIMUM distance: variables [x, t], minimize t
#     t >= x - a_i   ->   x - t <= a_i
#     t >= a_i - x   ->  -x - t <= -a_i
A_ub = [[1, -1] for ai in a] + [[-1, -1] for ai in a]
b_ub = a + [-ai for ai in a]
res = linprog(c=[0, 1], A_ub=A_ub, b_ub=b_ub, bounds=[(None, None), (0, None)])
print("minimax : x = %.1f, worst distance = %.1f" % (res.x[0], res.fun))
# minimax : x = 5.0, worst distance = 4.0      <- midpoint of the extremes

# (2) Minimize the SUM of distances: variables [x, t1, t2, t3], minimize sum t_i
A_ub, b_ub = [], []
for i, ai in enumerate(a):
    e = [0] * n; e[i] = -1                    # coefficient -1 on t_i
    A_ub.append([ 1] + e); b_ub.append( ai)   #  x - t_i <=  a_i
    A_ub.append([-1] + e); b_ub.append(-ai)   # -x - t_i <= -a_i
res = linprog(c=[0] + [1] * n, A_ub=A_ub, b_ub=b_ub,
              bounds=[(None, None)] + [(0, None)] * n)
print("min-sum : x = %.1f, total distance = %.1f" % (res.x[0], res.fun))
# min-sum : x = 4.0, total distance = 8.0      <- the median

The minimax objective is "fair" to the farthest customer. The sum objective is "efficient" for the group. Neither is wrong: they encode different values.

Key takeaways

• The objective maps each decision to one number to minimize or maximize; $\max f = -\min(-f)$.
• The choice of objective is a value judgment (efficiency, fairness, risk) and can change the optimal decision.
• Minimax and absolute-value objectives can be linearized with auxiliary variables.
• Keep units consistent: adding euros to hours in one objective is meaningless without explicit weights.

Practice

1. (Beginner) Rewrite $\max\ 4x_1 - 2x_2$ as a minimization problem. Is the optimal $x$ the same? Is the optimal value the same?
2. (Intermediate) Add a fourth customer at position 10 to the code. What happens to the min-sum solution? Explain why it is no longer unique.
3. (Interview) "A delivery company asks you to minimize total driving distance, and drivers complain that some routes are extremely long. How would you change the objective?" Hint: consider a minimax term, or a weighted combination of total and maximum route length.

Constraints and Feasible Regions

Why it matters. Constraints encode physical limits, business rules and logic. Their combined effect, the feasible region, determines whether a problem has a solution at all, and its shape determines how hard the problem is to solve.

Theory. A constraint is a condition that every acceptable decision must satisfy. A decision $x$ is feasible if it satisfies all constraints, and the feasible region (feasible set) $F$ is the set of all feasible decisions.

Constraint type	Form	Example
Inequality	$g(x) \le 0$	Mixer hours $3x_1 + 2x_2 \le 18$
Equality	$h(x) = 0$	Flow in = flow out at a warehouse
Bound	$l \le x \le u$	$0 \le x_j \le 100$
Integrality	$x_j \in \mathbb{Z}$	Whole trucks
Logical	"if $y = 1$ then $x \le 50$"	Use machine only if it is set up (Chapter 4)

Vocabulary used throughout OR:

• Binding (active) constraint at $x$: holds with equality, so there is no slack left.
• Non-binding constraint: holds strictly, with slack $= b_i - a_i^\top x > 0$.
• Redundant constraint: removing it does not change $F$.
• Hard constraint: must be satisfied (physics, law). Soft constraint: desirable, may be violated at a penalty, modeled by adding a deviation variable $s \ge 0$ and a cost on $s$ in the objective.

A problem falls into exactly one of three cases:

Status	Meaning
Infeasible	$F = \emptyset$: constraints contradict each other
Unbounded	$F \neq \emptyset$ but the objective can improve forever
Optimal	A best feasible solution exists (it may be one of many)

 Bakery: which constraints bind at the optimum (2,6)?

  x2
 6 *=========*  (2,6)  <- cake oven  2 x2 <= 12  BINDING (12 = 12)
   |########/ \
   |#######/   \_ mixer 3 x1 + 2 x2 <= 18        BINDING (18 = 18)
   |######/
 3 +#####* (4,3)
   |#####|  <- bread oven x1 <= 4                SLACK 2  (2 < 4)
   *-----*------> x1

import numpy as np

A = np.array([[1, 0], [0, 2], [3, 2]])
b = np.array([4, 12, 18])
names = ["bread oven", "cake oven", "mixer"]

def report(x):
    x = np.asarray(x)
    slack = b - A @ x                       # how much of each resource is left
    feasible = np.all(slack >= 0) and np.all(x >= 0)
    print(f"x = {x}: {'feasible' if feasible else 'INFEASIBLE'}")
    for name, s in zip(names, slack):
        tag = "binding" if np.isclose(s, 0) else ("violated" if s < 0 else "slack")
        print(f"   {name:<10} slack = {s:>3}  ({tag})")

report([2, 6])     # the optimal plan
report([4, 4])     # tempting, but uses 20 mixer hours
# x = [2 6]: feasible
#    bread oven slack =   2  (slack)
#    cake oven  slack =   0  (binding)
#    mixer      slack =   0  (binding)
# x = [4 4]: INFEASIBLE
#    bread oven slack =   0  (binding)
#    cake oven  slack =   4  (slack)
#    mixer      slack =  -2  (violated)

Binding constraints are the bottlenecks. Only extra cake-oven or mixer hours could increase profit, while extra bread-oven time is worthless at the current plan. Chapter 2 makes this precise with shadow prices.

Key takeaways

• The feasible region is the set of decisions satisfying every constraint, and its shape drives difficulty.
• Binding constraints are the bottlenecks at a solution, and non-binding ones have slack.
• Every problem is infeasible, unbounded, or has an optimum.
• Soft constraints are modeled with penalized deviation variables, hard constraints directly.

Practice

1. (Beginner) For the bakery, is $x = (3, 4)$ feasible? Which constraints bind?
2. (Intermediate) Write a two-variable LP that is infeasible and one that is unbounded. Run both through linprog and record the status and message fields.
3. (Interview) "A scheduling model returns 'infeasible' after a business user adds 40 new rules. How would you find the cause?" Hint: look up the idea of an Irreducible Infeasible Subsystem (IIS), or relax constraints as soft ones with penalties and see which ones get violated.

Local vs Global Optima

Why it matters. Most algorithms only look at the neighborhood of the current point, so they find local optima. Knowing when a local optimum is guaranteed to be global tells you whether you can trust a solver's answer or need further search.

Theory. For a minimization problem with feasible set $F$:

• $x^$ is a global minimum if $f(x^) \le f(x)$ for all $x \in F$.
• $x^$ is a local minimum if there is a radius $\varepsilon > 0$ such that $f(x^) \le f(x)$ for all $x \in F$ with $\lVert x - x^* \rVert < \varepsilon$.
• Strict versions use $<$ for $x \neq x^*$.

Analogy: a hiker in thick fog who always walks downhill will stop at the bottom of a valley, with no way of knowing whether a deeper valley lies beyond the next ridge.

  f(x) = x^4 - 3x^2 + x

  f |
  6 +                                          *
    |  *                                      /
  2 +   \                                    /
    |    \          .-'-.  <- local max     /
  0 +-----\--------/-----\-----------------/-------
    |      \      /       \               /
 -1 +       \    /         '-.__*__.-'         <- local min  x ~ +1.13, f ~ -1.07
    |        \  /
 -3.5        '*'                               <- GLOBAL min x ~ -1.30, f ~ -3.51
    +----+------+------+------+------+------+----> x
        -2    -1.3     0     0.17   1.13    2

def f(x):
    return x**4 - 3 * x**2 + x

def grad(x):
    return 4 * x**3 - 6 * x + 1

def gradient_descent(x, step=0.01, iters=2000):
    for _ in range(iters):
        x -= step * grad(x)      # always walk downhill
    return x

for start in [-2.0, -0.5, 0.5, 2.0]:
    x = gradient_descent(start)
    print(f"start {start:+.1f} -> x* = {x:+.4f}, f(x*) = {f(x):+.4f}")
# start -2.0 -> x* = -1.3008, f(x*) = -3.5139   <- global
# start -0.5 -> x* = -1.3008, f(x*) = -3.5139   <- global
# start +0.5 -> x* = +1.1309, f(x*) = -1.0702   <- stuck in a local minimum
# start +2.0 -> x* = +1.1309, f(x*) = -1.0702   <- stuck in a local minimum

When is local also global?

Situation	Local = global?
Convex $f$ minimized over a convex set	Yes: the fundamental theorem of convex optimization
Linear program	Yes (a special case of the above)
Non-convex objective (e.g., the function above)	No
Non-convex feasible set (integers, disjunctions, $x_1 x_2 = 1$)	No

Existence (Weierstrass theorem). A continuous function on a nonempty, closed and bounded (compact) set always attains a global minimum and maximum. Without compactness an optimum may not exist: $\min e^{-x}$ over $x \ge 0$ approaches 0 but never reaches it.

Strategies for non-convex problems: multi-start local search (as above), branch and bound for global optimality (Chapter 4), and metaheuristics that escape local optima (Chapter 11).

Key takeaways

• A local optimum is best within a neighborhood. A global optimum is best overall.
• Descent methods find local optima, and the starting point determines which one.
• For convex problems, including all LPs, every local optimum is global.
• Weierstrass: continuous objective plus nonempty compact feasible set guarantees an optimum exists.

Practice

1. (Beginner) Find all stationary points of $f(x) = x^3 - 3x$ and classify each as a local minimum, local maximum, or neither. Does $f$ have a global minimum on $\mathbb{R}$? On $[-3, 3]$?
2. (Intermediate) Run the gradient descent code with starts on a fine grid in $[-2, 2]$ and find the "watershed" point where the outcome switches between the two minima. How does it relate to the local maximum?
3. (Interview) "Your optimizer returns different answers every time you run it with random starts. What does that tell you about the problem, and what would you do?" Hint: think about non-convexity, and about lower bounds that could certify how far from optimal a solution is.

Classification of Optimization Problems

Why it matters. The class a problem belongs to determines which algorithms and solvers apply, whether a global optimum can be certified, and roughly how large an instance you can solve. Classifying a model correctly is the first step after formulating it.

Theory. Problems are classified along several independent dimensions:

Dimension	Options
Variables	Continuous / integer / binary / mixed
Functions	Linear / quadratic / general nonlinear
Convexity	Convex / non-convex
Constraints	Unconstrained / constrained
Data	Deterministic / stochastic
Objectives	Single / multiple
Time	Static / dynamic (multi-stage)

The major problem classes, all covered later in this guide:

Class	Form	Typical method	Difficulty	Chapter
LP	Linear $f, g$; continuous $x$	Simplex, interior point	Polynomial	2
Network flow	LP on a graph	Specialized combinatorial algorithms	Polynomial, very fast	3
ILP / MILP	Linear + integer variables	Branch and bound/cut	NP-hard in general	4
Convex NLP / QP	Convex $f$, convex $g_i$	Interior point, gradient methods	Efficiently solvable	5
Non-convex NLP	General nonlinear	Local methods, global B&B	Hard in general	5
Dynamic programming	Sequential decisions	Bellman recursion	Depends on state space size	6
Stochastic programming	Random parameters	Scenario decomposition	Hard (scale)	9

                      Are all variables continuous?
                       /                        \
                     yes                         no -> integer/binary present
                     /                                 |
        Are f and all g_i linear?               linear functions?
           /           \                        /             \
         yes            no                    yes              no
         /               \                    /                 \
       LP          convex f and g_i?       MILP / ILP          MINLP
   (Chapter 2)       /         \          (Chapter 4)       (hardest)
                   yes          no
                   /             \
            Convex NLP      Non-convex NLP
           (Chapter 5)     (local optima!)

Example: the same data in two classes. Integrality alone can change both the answer and the difficulty:

import numpy as np
from scipy.optimize import linprog, milp, LinearConstraint, Bounds

# max 5 x1 + 8 x2   s.t.  x1 + x2 <= 6,  5 x1 + 9 x2 <= 45,  x >= 0
c = np.array([-5, -8])                     # negate: both solvers minimize
A = np.array([[1, 1], [5, 9]])
b = np.array([6, 45])

# Continuous variables -> Linear Program
lp = linprog(c, A_ub=A, b_ub=b, bounds=[(0, None)] * 2)
print("LP   :", lp.x, "value", -lp.fun)          # [2.25 3.75] value 41.25

# Same data, integer variables -> Integer Linear Program
ip = milp(c, constraints=LinearConstraint(A, ub=b),
          integrality=np.ones(2),          # 1 = variable must be integer
          bounds=Bounds(0, np.inf))
print("ILP  :", ip.x, "value", -ip.fun)          # [0. 5.] value 40.0

# Naive rounding of the LP solution
for x in ([2, 4], [2, 3]):
    ok = np.all(A @ x <= b)
    print("round", x, "feasible" if ok else "infeasible", "value", -c @ x)
# round [2, 4] infeasible value 42     <- rounding up breaks a constraint
# round [2, 3] feasible value 34       <- rounding down is far from optimal (40)

The integer optimum $(0, 5)$ is nowhere near the LP optimum $(2.25, 3.75)$. Rounding is not a substitute for solving the integer problem, which is why Chapter 4 exists.

Key takeaways

• Classify by variable type, linearity, convexity, uncertainty, number of objectives and time structure.
• LP and convex problems are efficiently solvable with global optimality guarantees. Integer and non-convex problems are hard in general.
• The class determines the solver: linprog for LP, milp for MILP, NLP solvers for nonlinear models.
• Rounding an LP solution can be infeasible or badly suboptimal for the integer problem.

Practice

1. (Beginner) Classify: (a) $\min x_1^2 + x_2^2$ s.t. $x_1 + x_2 \ge 1$; (b) $\max 3x + 2y$ s.t. $x + y \le 4$, $x, y \in {0, 1, 2, \dots}$; (c) $\min \sin(x) + x^2/10$.
2. (Intermediate) Change the right-hand side 45 in the code to values from 40 to 54 and tabulate the LP value, the ILP value, and their difference (the integrality gap).
3. (Interview) "Given a new optimization problem, how do you decide which solver to use?" Hint: walk through the classification dimensions in order: variable types, linearity, convexity, size.

Modeling Assumptions and Limitations

Why it matters. Every model simplifies reality, and every simplification is an assumption that can fail. Knowing the assumptions behind a model tells you when to trust its output, when to choose a richer model class, and how to explain its limits to decision makers.

Theory. Linear programming rests on four classical assumptions. Each one, when it fails, leads to a different model class:

Assumption	Meaning	Violated when...	Remedy
Proportionality	Contribution of $x_j$ is $c_j x_j$, so doubling activity doubles cost/profit/resource use	Economies of scale, quantity discounts, setup costs	Piecewise linear models, fixed-charge MIP, NLP
Additivity	Total effect = sum of individual effects	Products interact (complements, cannibalization, shared setups)	Nonlinear terms, interaction variables
Divisibility	Variables can take fractional values	Trucks, people, yes/no decisions	Integer programming (Chapter 4)
Certainty	All parameters known exactly	Demand, prices, durations vary	Sensitivity analysis, stochastic and robust models (Chapters 2, 9)

Broader limitations apply to all OR models:

• Data quality: garbage in, garbage out. An optimizer will exploit any error in the data.
• Scope (model boundary): optimizing one department can hurt the whole company (suboptimization).
• Objective misspecification: "When a measure becomes a target, it ceases to be a good measure" (Goodhart's law). Minimizing cost alone may destroy service quality.
• Static view: a single-period model ignores how today's decision affects tomorrow (Chapter 6).
• Unmodeled human factors: acceptance, fairness and preferences that are hard to quantify.
• Computational limits: realistic models may be too large to solve exactly, forcing heuristics (Chapter 11).

Example: certainty and sensitivity. Suppose the bread profit of 3 is only an estimate. How sensitive is the bakery's plan?

from scipy.optimize import linprog

A_ub = [[1, 0], [0, 2], [3, 2]]
b_ub = [4, 12, 18]

# Certainty assumption: what if the bread profit estimate is wrong?
for bread_profit in [2, 4, 6, 7, 8, 10]:
    res = linprog([-bread_profit, -5], A_ub=A_ub, b_ub=b_ub)
    print(f"bread profit {bread_profit:>2}: plan {res.x}, profit {-res.fun:.0f}")
# bread profit  2: plan [2. 6.], profit 34
# bread profit  4: plan [2. 6.], profit 38
# bread profit  6: plan [2. 6.], profit 42
# bread profit  7: plan [2. 6.], profit 44
# bread profit  8: plan [4. 3.], profit 47   <- the optimal vertex jumps
# bread profit 10: plan [4. 3.], profit 55

The plan stays at $(2, 6)$ for a wide range of bread profits and then jumps to the adjacent vertex $(4, 3)$ once the profit exceeds 7.5. The threshold comes from geometry: $(2, 6)$ stays optimal while the objective ratio $c_1 / c_2 = c_1 / 5$ does not exceed the mixer constraint's coefficient ratio $3/2$. At exactly $c_1 = 7.5$ the objective is parallel to the mixer edge and every point on that edge is optimal. The plan is therefore robust to estimation errors within that range, and Chapter 2 shows how to compute such ranges exactly without re-solving.

Key takeaways

• LP assumes proportionality, additivity, divisibility and certainty. Each failure points to a richer model class.
• All models are limited by data quality, scope, objective choice and computational tractability.
• Sensitivity analysis tests how much the solution depends on uncertain parameters.
• Optimal solutions often stay unchanged over a range of parameter values and then jump to another vertex.

Practice

1. (Beginner) For each scenario, name the LP assumption that is violated: (a) the first 100 units cost 5 each, further units cost 4; (b) the number of aircraft to buy; (c) selling both coffee and milk increases the sales of each; (d) next month's demand is a forecast.
2. (Intermediate) Repeat the sensitivity experiment for the cake profit (keeping bread at 3) and find the range in which $(2, 6)$ remains optimal. Verify your answer geometrically using the slopes of the binding constraints.
3. (Interview) "Your model's optimal plan was rejected by the operations managers as 'unrealistic'. What could have gone wrong, and how would you investigate?" Hint: check each assumption (proportionality, additivity, divisibility, certainty), then the data and the model boundary, then ask which unwritten rules the managers apply.

<a id="2-linear-programming"></a>