Prescriptive analytics: optimization for decisions

1. Optimization: the mathematical core

The last deck of the course opens by naming the engine behind everything seen so far. Mathematical optimization (aka mathematical programming) is the selection of a best element from some set of available alternatives. An optimization problem consists of minimizing (maximizing) a real function by choosing input values from within an allowed set — and, the deck adds, more generally, optimization includes finding best available values of some objective function, given a defined domain (or input), where objective functions and domains can be very diverse.

The deck's title slide carries an epigraph by Euler that says the whole programme in Latin: nihil omnino in mundo contingit, in quo non maximi minimiQve ratio quaepiam eluceat — nothing at all happens in the world in which some principle of maximum or minimum does not appear. It is the philosophical ancestor of everything this chapter models: even when the world looks like a forecast, what an operator actually needs is the best thing to do next.

Position this inside the course. The first chapter promised three payoffs — faster decisions, fault fixing, productivity — and drew the loop that closes this course: operations produce data, statistics describe it, models predict it, and prescriptive optimization turns the prediction into a decision that acts back on operations. Chapter 9 fitted the parameters of models; this chapter is the other half of Part IV: it takes the fitted models' outputs — a forecast of demand, a predicted failure probability, an estimated processing time — as input data for a decision problem, and solves that problem. The forecast is not the product; the decision is.

2. The map of optimization problems

Optimization is a landscape, not a single algorithm. The deck lists the major subfields in one slide — a taxonomy worth keeping whole:

… and many more, the deck says. Two remarks tie the taxonomy to this course. First, integer programming and combinatorial optimization are the subfields where most operational decisions live — you cannot send 3.7 trucks or assign half a nurse — and they are exactly the ones that will dominate the rest of this chapter. Second, stochastic and robust programming are the bridge back to chapters 2–7: a forecast is a random variable, and a decision model that consumes it either treats its uncertainty explicitly (stochastic) or insists the decision stay valid for a whole range of scenarios (robust).

3. Nonlinear, linear and integer linear optimization

The deck walks the taxonomy from the most general to the most rigid. Nonlinear optimization: the objective function or the constraints or both can have nonlinear components, and the problem can be constrained or unconstrained. The deck's one-liner example: min x² − 2x + 3, s.t. x² ≤ 9. The unconstrained minimum of the parabola is x = 1 (derivative 2x − 2 = 0), which lies inside the constraint x ∈ [−3, 3] — so the constraint is inactive and the answer is x = 1 with f = 2. That is the pattern to internalize: constraints only matter when they bite.

Linear optimization: both the objective function and the constraints are modelled by linear functions; constraints (there have to be some) are in the form of a system of equations or inequalities. The deck's example is the one the widget below draws:

min    −x₁ − x₂
s.t.   x₁ − x₂ ≤ 1.5
       −7x₁ + 5x₂ ≤ −17.5
       −x₂ ≤ −4.5
       x₁, x₂ ≥ 0

Rewritten as half-planes: x₂ ≥ x₁ − 1.5, x₂ ≤ 1.4·x₁ − 3.5, x₂ ≥ 4.5, plus nonnegativity. The objective moves in the (1, 1) direction — the north-east — and the feasible region is the wedge between the three half-planes.

Integer linear optimization: linear optimization problems where (some of) the variables are required to take on integer values. If all variables are integer: Integer Linear Programming (ILP); if some of the variables are continuous: Mixed Integer Linear Programming (MILP). The same example with "and integer" added: the optimal vertex of the continuous relaxation is generally not an integer point, and the optimum over the lattice can sit far from it — which is why the deck warns that integer programming is in general much more difficult than regular linear programming.

4. P, NP and the size of real instances

The deck's slide on complexity theory is two Venn diagrams: the world P = NP and the world P ≠ NP, with the NP-complete region inside NP but outside P. The deck does not take sides — it draws both and moves on, because either way the operational conclusion is the same. The next slide says it: the size of the instances met in real world applications rules out the possibility of solving them to optimality in an acceptable time. Nevertheless, these instances must be solved. Thus the need to look for suboptimal solutions, provided that they are of "acceptable quality" and that they can be found in "acceptable time".

The consequence is the working contract of the whole chapter: acceptable quality in acceptable time. Every method that follows — constructive heuristics, local search, metaheuristics — is a bid on that contract: it gives up the guarantee of optimality (which is unattainable anyway at real sizes) in exchange for solutions that are good enough, found fast enough. Chapter 9's black-box optimizers were exactly this in the parameter space; here the same machinery acts on decision spaces.

5. Combinatorial optimization: components, feasibility and cost

The deck then gives the formal vocabulary that all the rest of the chapter speaks. A combinatorial optimization problem is defined over a set C = {c₁, …, cₙ} of basic components. A solution of the problem is a subset S ⊂ C; F ⊆ 2^C is the subset of feasible solutions (a solution S is feasible iff S ∈ F). z: 2^C → ℝ is the cost function. The objective is to find a minimum cost feasible solution S°, i.e., to find S° ∈ F such that z(S°) ≤ z(S), ∀S ∈ F. Failing this, the algorithm anyway returns the best feasible solution found, S* ∈ F.

Three moves in this definition matter. First, a solution is a subset of components — the "decision variables" of section 3 become, in combinatorial language, the choice of which components to include. Second, feasibility is a separate layer: F is only a part of 2^C, and an algorithm that proposes S ∉ F has failed regardless of cost. Third — and this is the sentence the whole chapter is built on — the objective is S°, but the algorithm may fail to find it and then returns the best feasible solution found, S*. The "acceptable quality" of section 4 is here made formal: S* is the honest deliverable, and the quality gap z(S*) − z(S°) is what heuristics try to shrink.

6. The Travelling Salesman Problem

The canonical combinatorial problem gets the deck's full treatment. n locations given (home of salesperson and n−1 customers); [d_ij] distance matrix given, contains distances (Km, time, …) between each pair of locations. The TSP asks for the shortest tour starting home, visiting each customer once and returning home. Symmetric and asymmetric versions exist (depending on whether d_ij = d_ji).

Why it is the canonical case: the explosion is immediate. n nodes ⇒ (n−1)! possible tours (for the symmetric undirected case, half that — the direction of travel does not matter). The deck's table: n = 5 → 120; n = 10 → 3,628,800; n = 15 → 1.31·10¹²; n = 20 → 2.43·10¹⁴. "Difficult to say… which is the shortest tour" — indeed.

Then the deck does the formal translation that connects TSP to section 5: TSP is defined over a weighted (un)directed graph G=(V,E,D), where V is the set of vertices, E is the set of edges and D is the set of edge weights. The component set corresponds to E (C=E); F corresponds to the set of Hamiltonian cycles in G; z(S) is the sum of the weights associated to the edges belonging to the solution S. A tour is a subset of edges — components — that happens to form a cycle visiting every vertex once.

7. Vehicle routing and the generalized assignment problem

From one salesperson to a fleet. The Vehicle Routing Problem (VRP) is a classic optimization problem in logistics. Goal: determine optimal routes for a fleet of vehicles delivering goods to customers. The constraints multiply: vehicle capacities, and may also include delivery time windows, customer service times, driver working hours / route duration limits, split deliveries, pickup and delivery, …. The objective stays the same shape: minimize total distance, cost, or time. Input data: V = {0, …, n}: set of nodes (0 = depot, others = customers); E: set of edges; c_ij: cost (distance) from node i to node j; q_i: demand of customer i; Q: vehicle capacity; K: number of vehicles. "Mathematics gets more complicated…" — the deck's dry summary of moving from one vehicle to K.

The second model of the section is the one with a full math formulation on the slides: the Generalized Assignment Problem (GAP). It is modeled on a bipartite complete directed graph G(I ∪ J, A), where the vertices in I denote the potential facilities, the vertices in J the customers, and the arcs in A ⊂ I × J are associated with the material flows between the potential facilities and the demand points. Demand may be divisible or not. The notation: Q_i, i ∈ I, the capacity of facility i; q_ij, j ∈ J, the demand of customer j when serviced by facility i; x_ij a decision variable representing the (percentual) amount of product sent from site i to demand point j; c_ij the cost of transporting all products from site i to customer j. The deck's math model (MIP):

z(GAP) = min  ∑i∈I ∑j∈J cij xij                    (1)
subject to   ∑i∈I xij = 1        j ∈ J      (2)
             ∑j∈J qij xij ≤ Qi   i ∈ I      (3)
             xij ∈ {0,1}  (or xij ≥ 0)  i ∈ I, j ∈ J   (4)

Read constraint by constraint, as the deck does: variables x_ij determine customer allocations to facilities; the objective function (1) is the cost for transporting goods from facilities to customers; constraints (2) ensure that each customer request is satisfied; constraints (3) force the total request to a facility not to exceed the corresponding capacity; constraints (4) could be either x_ij ≥ 0 — clients' service can be split — or x_ij ∈ {0,1} — clients are singly assigned to a warehouse. This is a MILP exactly as section 3 defined it, and the widget below solves a small instance of it by hand with the tools of the next two sections.

8. Constructive heuristics

With the models in place, the algorithms begin. The first family builds a solution from nothing:

function Constructive():
1. Sort the components in C by increasing costs.
2. Set S* = ∅ and i = 1.
3. Repeat
        if (S* ∪ cᵢ is a feasible partial solution)
            then S* = S* ∪ cᵢ.
        i = i + 1.
   until S* ∈ F.

Two honest caveats from the deck. A constructive approach can yield optimal solutions for certain problems, e.g., the MST — the minimum spanning tree is exactly the case where the sorted-edge greedy (Kruskal's algorithm) is provably optimal. In other cases it could be unable to construct a feasible solution — the "until S* ∈ F" loop can simply run out of components. For the TSP the deck lists the classic constructive strategies: Farthest Insertion, the Nearest Neighbor, the Best Insertion or the Sweep algorithms.

The nearest neighbor is worked in full because its weakness is instructive:

function NearestNeighbor():
0. Start from any city; e.g. the first one
1. Find nearest neighbor (to the last city) not already visited
2. repeat 1 until all cities are visited.
   then connect last city with starting point

Ties can be broken arbitrarily.
Starting from different starting points gives different solutions!
Some will be bad, some will be good.

And the deck's punchline: GREEDY ⇒ the last connections tend to be long and improvement heuristics should be used afterwards. The greedy has perfect local information and zero foresight: the first choices look excellent, and the bill arrives at the end, when the remaining cities force long closing edges. That is why the constructive phase is only the first phase — the next section's local search exists to pay the bill.

9. Local search: neighborhoods, 2-opt and 3-opt

The improvement machinery starts with the neighborhood. The neighborhood of a solution S, N(S), is a subset of 2^C defined by a neighborhood function N: 2^C → 2^2C. Often only feasible solutions are considered, so neighborhood functions N: F → 2^F. The deck's warning: the specific function used has a deep impact on the algorithm performance and its choice is left to the algorithm designer — the neighborhood is the designer's first real decision.

function LocalSearch():
1. Generate an initial feasible solution S.
2. Find S' ∈ N(S), such that z(S') = min z(S^), ∀S^ ∈ N(S).
3. If z(S') < z(S) then S = S'
   go to step 2.
4. S* = S.

The update of a solution in step 3 is called a move from S to S' — and the deck notes it could be made to the first improving solution found rather than the best of the whole neighborhood (the difference between first-improvement and best-improvement descent). Two properties anchor the method: there are problems for which a local search approach guarantees to find an optimal solution (ex. the simplex algorithm) — LP is one giant local search over vertices — and for combinatorial problems it generally does not, because it stops at the first solution whose whole neighborhood is worse, the local optimum.

For the TSP, the deck names the two classic neighborhoods: 2-opt and 3-opt, which take a solution (a list of n vertices) and exhaustively swap the elements of all pairs or triplets of arcs in it — and the more sophisticated neighborhoods that produce the more effective heuristics, among which Lin and Kernighan [LK73]. The widget below runs first-improvement 2-opt exactly as the deck describes it: try all inversions of a subsequence; if one improves, apply it and start again from the beginning.

10. Metaheuristics: diversification versus intensification

The deck motivates its last big family with the failure mode just observed. Simple heuristics can perform very well, but can get trapped in poor quality solutions (local optima). The response, starting in the mid-1970s, is a class of algorithms that manage the trade-off between search diversification, when search is going on in bad regions of the search space, and intensification, aimed at finding the best solutions within the region being analyzed. These algorithms have been named metaheuristics.

The deck gives the two canonical definitions, worth keeping verbatim because they say the same thing with different emphasis:

The deck then lists the relevant family: iterated local search, simulated annealing, tabu search, GRASP, variable neighborhood search, ALNS, CMSA, ejection chains, genetic algorithms, evolution strategies, ant colony optimization, scatter search, particle swarm optimization, … — particle swarm already met in chapter 9. And the closing remark of the slide is the deck's own research program in one line: all relevant metaheuristics have been hybridized with MP [mathematical programming], generating matheuristics.

11. Iterated local search and simulated annealing

Iterated Local Search (ILS) is the simplest metaheuristic and the honest baseline. Repeat local searches starting from different initial solutions. The refinement that makes it work: initial solutions should not be random to increase effectiveness. Bespoke perturbations of previous solutions or of problem data work best. Instead of restarting blindly, ILS perturbs the current local optimum — a big jump in neighborhood space, not a random draw — and runs local search again:

function IteratedLocalSearch():
1. Generate an initial feasible solution x⁰
2. Apply a local search on x⁰ obtaining a local optimum x○
3. repeat
4.   Set x = Perturbation(x○, history)
5.   Set x = LocalSearch(x)
6.   Set x○ = AcceptanceCriterion(x○, x, history)
7.   if (z(x*) > z(x○)) then set x* = x○
8. until termination condition
9. return x*

Simulated Annealing (SA), a simple and widely used metaheuristic, replaces the "restart" strategy with a probabilistic escape: a worse move is accepted with a probability that decays as the search cools:

function SimulatedAnnealing(T):
1. Generate a feasible solution x; Set x* = x
2. Generate a feasible solution x′ ∈ N(x)
3. if (z(x′) < z(x)) then set x = x′
4. if (z(x*) > z(x)) then x* = x
5. else set x = x′ with probability p = exp(−(z(x′) − z(x))/kT)
6. if (annealing condition) then decrease T
7. if not (terminating condition) then go to 2
8. return x*

The deck's example is a TSP trace on the benchmark instance ulysses16: the plot shows z values over iterations, oscillating around a descending trend — the accepted uphill moves are what let the search climb out of one basin and into a deeper one. The widget below reproduces that trace live on a small TSP with the 2-opt neighborhood of section 9.

12. Tabu search

Where SA randomizes the uphill decision, Tabu Search (TS) [GL97] makes it deterministic and memory-driven. TS escapes from local minima by moving onto the best solution of the neighborhood at each iteration, even though it is worse than the current one. A memory structure called tabu list, TL, forbids to return to already explored solutions:

function TabuSearch():
1. Generate an initial feasible solution S,
   set S* = S and initialize TL = ∅
2. Find S′ ∈ N(S), such that
   z(S′) = min { z(S^), ∀S^ ∈ N(S), S^ ∉ TL }
3. S = S′, TL = TL ∪ {S}, if (z(S*) > z(S)) set S* = S
4. If not (end condition) go to step 2

Two details make it work. First, "best of the neighborhood even though it is worse": when no neighbor improves, TS still moves — to the least-bad one — which is how it walks out of a basin that 2-opt would call the answer. Second, the tabu list is what keeps that walk from immediately returning: the move that just made things worse is exactly the move a naive "go to the best neighbor" rule would immediately undo, so TS forbids the recently used moves for a while (a tenure). In practice the list stores features of moves rather than whole solutions, and an aspiration criterion — a tabu move that would beat the global best is allowed anyway — prevents the memory from blocking the best move available.

13. GRASP, variable neighborhood search and ALNS

Three more members of the family, each with a distinct idea. GRASP (Greedy Randomized Adaptive Search Procedure) restarts search from another promising region of the search space as soon as a local optimum is reached. GRASP consists in a multistart approach with a suggestion on how to construct the initial solutions — the construction itself is randomized:

function GRASP():
1. Build a solution S (=S* in the first iteration) by a constructive
   greedy randomized procedure based on candidate lists.
2. Apply a local search procedure starting from S and producing S′.
3. If z(S′) < z(S*) then S* = S′.
4. If not (end condition) go to step 1.

The trace the deck shows for GRASP on ulysses16 plots three lines: zcons (the cost of each constructed solution), zloc (the cost after the local search that follows each construction) and zbest (the best found so far). The pattern is the signature of multistart: zcons jumps around (randomized construction), zloc always lands below its zcons (local search does its job), and zbest descends in steps, one per restart.

Variable Neighborhood Search (VNS) explores a sequence of neighborhoods at each iteration (called a Variable Neighborhood Descent, VND), followed by a perturbation of the incumbent solution. The insight is the deck's best single sentence of the section: a local minimum with respect to one neighborhood is not necessarily a local minimum for other neighborhoods. Only a global optimum is a local optimum for any neighborhood. So VND applies the chosen neighborhoods in a fixed order, restarting from the beginning each time an improvement is found; VNS wraps it with a perturbation step (re-instantiating search on another seed) and an acceptance step (evaluating exploration versus exploitation).

function VariableNeighborhoodDescent(x):
1. Define Nₘ, k = 1, ..., k_max          // the neighborhoods
2. Set k = 1
3. repeat
4.   Set x′ = FindBestNeighbor(x, Nₘ)
5.   if (z(x′) < z(x)) then x = x′; x* = x; k = 1
6.   else k = k + 1
7. until k > k_max
8. return x*

function VariableNeighborhoodSearch():
1. Generate an initial feasible solution x
2. Set x* = x
3. repeat
4.   Set x′ = VariableNeighborhoodDescent(x)   // can be any local search
5.   if (z(x*) > z(x′)) then set x* = x′
6.   Set x = Perturb(x′)                      // perturbation
7.   Set x = Accept(x, history)                // acceptance
8. until termination condition
9. return x*

Finally ALNS — Adaptive Large Neighborhood Search — maintains sets of destroy and reconstruct operators, records the performance of each, and dynamically adjusts the weights used to choose (in probability or otherwise) which method to use at each iteration:

function ALNS():
// Given: an initial solution
//         a set of destroy methods Ω⁻ and associated weights w⁻
//         a set of repair methods Ω⁺ and associated weights w⁺
1. Repeat until stop-criterion:
      Choose a method from Ω⁻ and a method from Ω⁺
      Destroy part of solution
      Repair solution
      Accept or reject new solution
      Update destroy and repair weights w⁻ and w⁺
2. Return best solution encountered

ALNS is the metaheuristic as learner: the algorithm itself is fixed, and what adapts is which of its tools gets used, based on their recorded performance. The destroy/repair pair is a very large neighborhood — "destroy part of the solution, repair it" — searched adaptively.

14. The loop closes: prescriptive analytics in the exam project

This is the last deck, and it closes the loop drawn on the course index: operations produce data; descriptive statistics summarize it; predictive models extend it into the future; prescriptive optimization turns the forecast into decisions that act back on operations. The index plate's dashed vermilion line — "decision taken NOW acts back on operations" — is exactly what this chapter's models formalize.

Where chapter 9 ended with fitting the parameters of models, this chapter began where those models' outputs become inputs: the forecast of demand is the qᵢ of a GAP or the q_i of a VRP; the predicted failure probability of a machine decides whether it is in the maintenance schedule; the predicted service time decides how crews are rostered. The vocabulary of this chapter is what lets you name the model — LP, ILP, MILP, combinatorial — and then choose the method: exact when the instance is small (simplex, branch-and-bound, the closed forms of chapter 9), heuristic when it is not (constructive → local search → metaheuristic).

And one last thread to the beginning of the course. The first chapter's three payoffs — faster decisions, fault fixing, productivity — are the objective function of this chapter: every model here is a formal way of choosing the decision that delivers them. The Euler epigraph of slide 2 — nothing happens in the world without some principle of maximum or minimum — is, in the end, the one-line summary of the whole course: measurements describe the world (chapters 2–8), and optimization decides in it (chapters 9–10).

Check your understanding