CLRS-Lean
Chapters 1–35 · CLRS Fourth Edition
Chapter 29. Linear Programming

Chapter 29 — Linear Programming

CLRS, fourth edition · Lean 4 formalization

The proofs below use the models and assumptions described in the scope and implementation notes.

Imports

import CLRSLean.FourthEdition.Chapter_29.Section_29_1_Standard_And_Slack_Forms.Definitions
import CLRSLean.FourthEdition.Chapter_29.Section_29_1_Standard_And_Slack_Forms.SlackVariables
import CLRSLean.FourthEdition.Chapter_29.Section_29_1_Standard_And_Slack_Forms.Equivalence
import CLRSLean.FourthEdition.Chapter_29.Section_29_1_Standard_And_Slack_Forms.Normalization

29.1. Standard and Slack Forms

The represented foundation defines standard-form maximization programs over finite real matrices, and converts the general linear-programming form to standard form while preserving feasibility and the objective. Subsequent child modules add slack variables and prove their exact feasibility equivalence.

The canonical Chapter 29 guide additionally imports SolverWrapper, which connects this normalization to the initialized SIMPLEX development in the supplementary online material. Its basic/nonbasic dictionary model and GeneralLP.solve_complete certify infeasible, optimal or unbounded outcomes. This section is the formulation layer of that larger development; it does not claim a polynomial SIMPLEX bound or executable exact-real choices.

Implementation details

The split proof layers remain available outside the main sidebar:

              namespace CLRSnamespace Chapter29end Chapter29end CLRS

Definitions and proofs

CLRSLean.FourthEdition.Chapter_29.Section_29_1_Standard_And_Slack_Forms.Definitions

29.1 Standard-form linear programs

This module defines the finite real matrix model used by Chapter 29. A standard-form program maximizes cᵀx subject to Ax ≤ b and 0 ≤ x.

Main declarations:

IsNonnegative: pointwise nonnegativity of a finite vector.
StandardLP: coefficients, bounds, and objective coefficients.
StandardLP.IsFeasible: primal standard-form feasibility.
StandardLP.objective: the value cᵀx.

Downstream layers:

Slack-variable equivalence is proved in later Section 29.1 modules.
Basic/nonbasic dictionaries and SIMPLEX are developed in Sections 29.3--29.5.

              namespace CLRSnamespace Chapter29open Matrix

A finite real vector is nonnegative when every coordinate is nonnegative.

def IsNonnegative {n : ℕ} (x : Fin n → ℝ) : Prop :=
  ∀ j, 0 ≤ x j

A maximization linear program in CLRS standard form: maximize cᵀx subject to Ax ≤ b and 0 ≤ x.

The constraint coefficient matrix.

The constraint right-hand side.

The objective coefficient vector.

              structure StandardLP (m n : ℕ) where  A : Matrix (Fin m) (Fin n) ℝ  b : Fin m → ℝ  c : Fin n → ℝ

namespace StandardLP

A vector is primal feasible when it is nonnegative and satisfies every row inequality of the standard-form program.

def IsFeasible {m n : ℕ} (P : StandardLP m n) (x : Fin n → ℝ) : Prop :=
  IsNonnegative x ∧ ∀ i, (P.A *ᵥ x) i ≤ P.b i

The objective value cᵀx of a standard-form assignment.

def objective {m n : ℕ} (P : StandardLP m n) (x : Fin n → ℝ) : ℝ :=
  P.c ⬝ᵥ x

              end StandardLPend Chapter29end CLRS

CLRSLean.FourthEdition.Chapter_29.Section_29_1_Standard_And_Slack_Forms.SlackVariables

29.1 Slack-variable construction

The canonical slack vector for an assignment x is b - Ax. Primal feasibility makes this vector nonnegative and converts every row inequality into an equality.

Main results:

slack_nonnegative_of_feasible.
slack_equation.
slackExtension_of_feasible.

              namespace CLRSnamespace Chapter29open Matrixnamespace StandardLP

The canonical slack vector b - Ax.

def slack {m n : ℕ} (P : StandardLP m n) (x : Fin n → ℝ) : Fin m → ℝ :=
  fun i => P.b i - (P.A *ᵥ x) i

A nonnegative slack extension satisfies Ax + s = b coordinatewise.

def IsSlackExtension {m n : ℕ} (P : StandardLP m n)
    (x : Fin n → ℝ) (s : Fin m → ℝ) : Prop :=
  IsNonnegative x ∧ IsNonnegative s ∧
    ∀ i, (P.A *ᵥ x) i + s i = P.b i

A feasible assignment has a nonnegative canonical slack vector.

theorem slack_nonnegative_of_feasible {m n : ℕ} {P : StandardLP m n}
    {x : Fin n → ℝ} (hx : P.IsFeasible x) :
    IsNonnegative (P.slack x) := by
  intro i
  exact sub_nonneg.mpr (hx.2 i)

The canonical slack vector satisfies Ax + slack(x) = b.

theorem slack_equation {m n : ℕ} (P : StandardLP m n) (x : Fin n → ℝ) :
    ∀ i, (P.A *ᵥ x) i + P.slack x i = P.b i := by
  intro i
  simp [slack]

Every primal-feasible assignment extends canonically to a nonnegative equality-form assignment.

theorem slackExtension_of_feasible {m n : ℕ} {P : StandardLP m n}
    {x : Fin n → ℝ} (hx : P.IsFeasible x) :
    P.IsSlackExtension x (P.slack x) := by
  exact ⟨hx.1, slack_nonnegative_of_feasible hx, P.slack_equation x⟩

              end StandardLPend Chapter29end CLRS

CLRSLean.FourthEdition.Chapter_29.Section_29_1_Standard_And_Slack_Forms.Equivalence

29.1 Standard/slack feasibility equivalence

This module proves the exact semantic bridge Ax ≤ b ↔ ∃ s ≥ 0, Ax + s = b for nonnegative decision variables. The slack vector is uniquely determined by the decision assignment.

Main results:

isFeasible_iff_exists_slackExtension.
slackExtension_eq_slack.
existsUnique_slackExtension_iff.

              namespace CLRSnamespace Chapter29namespace StandardLP

Eliminating nonnegative slack variables recovers primal feasibility.

theorem feasible_of_slackExtension {m n : ℕ} {P : StandardLP m n}
    {x : Fin n → ℝ} {s : Fin m → ℝ} (hxs : P.IsSlackExtension x s) :
    P.IsFeasible x := by
  refine ⟨hxs.1, ?_⟩
  intro i
  have hs : 0 ≤ s i := hxs.2.1 i
  have heq := hxs.2.2 i
  linarith

Standard-form feasibility is equivalent to the existence of a nonnegative slack vector satisfying the equality system.

theorem isFeasible_iff_exists_slackExtension {m n : ℕ}
    (P : StandardLP m n) {x : Fin n → ℝ} :
    P.IsFeasible x ↔ ∃ s, P.IsSlackExtension x s := by
  constructor
  · intro hx
    exact ⟨P.slack x, slackExtension_of_feasible hx⟩
  · rintro ⟨s, hxs⟩
    exact feasible_of_slackExtension hxs

Every slack extension equals the canonical vector b - Ax.

theorem slackExtension_eq_slack {m n : ℕ} {P : StandardLP m n}
    {x : Fin n → ℝ} {s : Fin m → ℝ} (hxs : P.IsSlackExtension x s) :
    s = P.slack x := by
  funext i
  have heq := hxs.2.2 i
  simp only [slack]
  linarith

A standard-form assignment is feasible exactly when it has a unique nonnegative slack extension.

theorem existsUnique_slackExtension_iff {m n : ℕ}
    (P : StandardLP m n) {x : Fin n → ℝ} :
    P.IsFeasible x ↔ ∃! s, P.IsSlackExtension x s := by
  constructor
  · intro hx
    refine ⟨P.slack x, slackExtension_of_feasible hx, ?_⟩
    intro s hxs
    exact slackExtension_eq_slack hxs
  · rintro ⟨s, hxs, _⟩
    exact feasible_of_slackExtension hxs

              end StandardLPend Chapter29end CLRS

CLRSLean.FourthEdition.Chapter_29.Section_29_1_Standard_And_Slack_Forms.Normalization

29.1 General-form normalization

CLRS §29.1 lists the general form of a linear program (equations (29.16)-- (29.19)): the objective may be maximized or minimized, each constraint may be =, ≥, or ≤, and each variable is either nonnegative or free. This module converts any such program into a finite standard-form StandardLP (maximize a linear objective subject to Ax ≤ b, x ≥ 0) while preserving feasibility and the objective value up to the maximization/minimization sign.

Main results:

GeneralLP: the general-form program representation.
GeneralLP.toStandardLP: the normalization.
GeneralLP.lift / GeneralLP.proj: the variable expansion and its inverse.
GeneralLP.feasible_iff_lift: feasibility is exactly preserved.
GeneralLP.objective_lift: the objective is preserved up to sign.
GeneralLP.solve / GeneralLP.solve_complete: a canonical main-text solver wrapper around the initialized SIMPLEX.

Each original variable x j is expanded to two nonnegative variables x j⁺ and x j⁻ with x j = x j⁺ - x j⁻; for a nonnegative variable an extra constraint forces x j⁻ = 0.

              namespace CLRSnamespace Chapter29open Matrixopen scoped BigOperators

The relation of a general-form constraint row.

inductive ConstraintRel where
  | le | eq | ge
  deriving DecidableEq, Repr

A general-form linear program: optimize cᵀx over constraint rows of mixed relation and nonnegative or free variables (CLRS (29.16)--(29.19)).

Number of variables.

Number of constraints.

true when the objective is maximized, false when minimized.

Objective coefficients.

The relation of each constraint row.

Constraint coefficient matrix.

Constraint right-hand side.

true when variable j is free (unrestricted).

              structure GeneralLP where  n : ℕ  m : ℕ  maximize : Bool  c : Fin n → ℝ  rel : Fin m → ConstraintRel  A : Matrix (Fin m) (Fin n) ℝ  b : Fin m → ℝ  free : Fin n → Bool

namespace GeneralLPvariable (G : GeneralLP)

An assignment is feasible when every nonnegative variable is nonnegative and every constraint holds under its relation.

def IsFeasible (x : Fin G.n → ℝ) : Prop :=
  (∀ j, ¬ G.free j → 0 ≤ x j) ∧
    ∀ i, match G.rel i with
      | ConstraintRel.le => (G.A *ᵥ x) i ≤ G.b i
      | ConstraintRel.eq => (G.A *ᵥ x) i = G.b i
      | ConstraintRel.ge => G.b i ≤ (G.A *ᵥ x) i

The objective cᵀx (before the maximization/minimization sign).

def objective (x : Fin G.n → ℝ) : ℝ :=
  G.c ⬝ᵥ x

The sign +1 for maximization and -1 for minimization, so that objectiveSign · objective is always the quantity being maximized.

def objectiveSign : ℝ :=
  if G.maximize then 1 else -1

An assignment is optimal when it is feasible and its signed objective dominates every other feasible assignment.

def IsOptimal (x : Fin G.n → ℝ) : Prop :=
  G.IsFeasible x ∧
    ∀ z, G.IsFeasible z → G.objectiveSign * G.objective z ≤ G.objectiveSign * G.objective x

The signed objective is unbounded above on feasible assignments.

def IsUnbounded : Prop :=
  ∀ M : ℝ, ∃ x, G.IsFeasible x ∧ M < G.objectiveSign * G.objective x

Positive part of the expanded variable j.

abbrev posIndex (j : Fin G.n) : Fin (G.n + G.n) :=
  Fin.castAdd G.n j

Negative part of the expanded variable j.

abbrev negIndex (j : Fin G.n) : Fin (G.n + G.n) :=
  Fin.natAdd G.n j

Upper form of constraint row i.

abbrev upperIndex (i : Fin G.m) : Fin ((G.m + G.m) + G.n) :=
  Fin.castAdd G.n (Fin.castAdd G.m i)

Lower form of constraint row i.

abbrev lowerIndex (i : Fin G.m) : Fin ((G.m + G.m) + G.n) :=
  Fin.castAdd G.n (Fin.natAdd G.m i)

The nonnegativity-lock row for variable j.

abbrev lockIndex (j : Fin G.n) : Fin ((G.m + G.m) + G.n) :=
  Fin.natAdd (G.m + G.m) j

The normalized constraint matrix: each equality contributes an upper and lower inequality, each ≥ contributes a lower inequality, and each nonnegative variable contributes a lock row forcing its negative part to zero.

def normalizedA : Matrix (Fin ((G.m + G.m) + G.n)) (Fin (G.n + G.n)) ℝ :=
  fun i =>
    Fin.addCases (motive := fun _ => Fin (G.n + G.n) → ℝ)
      (fun ij =>
        Fin.addCases (motive := fun _ => Fin (G.n + G.n) → ℝ)
          (fun iU =>
            fun j =>
              Fin.addCases (motive := fun _ => ℝ)
                (fun jP => if G.rel iU = ConstraintRel.ge then 0 else G.A iU jP)
                (fun jN => if G.rel iU = ConstraintRel.ge then 0 else -(G.A iU jN))
                j)
          (fun iL =>
            fun j =>
              Fin.addCases (motive := fun _ => ℝ)
                (fun jP => if G.rel iL = ConstraintRel.le then 0 else -(G.A iL jP))
                (fun jN => if G.rel iL = ConstraintRel.le then 0 else G.A iL jN)
                j)
          ij)
      (fun jL =>
        fun j =>
          Fin.addCases (motive := fun _ => ℝ)
            (fun _ => 0)
            (fun jN => if G.free jL then 0 else (if jN = jL then 1 else 0))
            j)
      i

The normalized right-hand side, mirroring normalizedA.

def normalizedB : Fin ((G.m + G.m) + G.n) → ℝ :=
  fun i =>
    Fin.addCases (motive := fun _ => ℝ)
      (fun ij =>
        Fin.addCases (motive := fun _ => ℝ)
          (fun iU => if G.rel iU = ConstraintRel.ge then 0 else G.b iU)
          (fun iL => if G.rel iL = ConstraintRel.le then 0 else -(G.b iL))
          ij)
      (fun _ => 0)
      i

The normalized objective: sign · c on positive parts and -sign · c on negative parts.

def normalizedC : Fin (G.n + G.n) → ℝ :=
  fun j =>
    Fin.addCases (motive := fun _ => ℝ)
      (fun jP => G.objectiveSign * G.c jP)
      (fun jN => -(G.objectiveSign * G.c jN))
      j

The normalized standard-form program.

abbrev toStandardLP : StandardLP ((G.m + G.m) + G.n) (G.n + G.n) :=
  { A := G.normalizedA, b := G.normalizedB, c := G.normalizedC }

Expand an assignment to the nonnegative positive/negative parts.

def lift (x : Fin G.n → ℝ) : Fin (G.n + G.n) → ℝ :=
  fun j =>
    Fin.addCases (motive := fun _ => ℝ)
      (fun jP => max (x jP) 0)
      (fun jN => max (-(x jN)) 0)
      j

Collapse an expanded assignment back to the signed difference.

def proj (x' : Fin (G.n + G.n) → ℝ) : Fin G.n → ℝ :=
  fun j => x' (G.posIndex j) - x' (G.negIndex j)

max a 0 - max (-a) 0 = a.


lemma max_sub_max_neg (a : ℝ) : max a 0 - max (-a) 0 = a := by
  rcases le_total 0 a with h | h
  · have hneg : -a ≤ 0 := by linarith
    simp [max_eq_left h, max_eq_right hneg]
  · have hneg : 0 ≤ -a := by linarith
    have ha : a ≤ 0 := by linarith
    simp [max_eq_right ha, max_eq_left hneg]

Index reductions

              lemma lift_pos (x : Fin G.n → ℝ) (j : Fin G.n) :
    G.lift x (G.posIndex j) = max (x j) 0 := by
  simp [lift, posIndex, Fin.addCases, Fin.castAdd, Fin.castLT]lemma lift_neg (x : Fin G.n → ℝ) (j : Fin G.n) :
    G.lift x (G.negIndex j) = max (-(x j)) 0 := by
  simp [lift, negIndex, Fin.addCases, Fin.natAdd]

lemma normalizedA_upper_pos (iU : Fin G.m) (jP : Fin G.n) :
    G.normalizedA (G.upperIndex iU) (G.posIndex jP) =
      if G.rel iU = ConstraintRel.ge then 0 else G.A iU jP := by
  have hlt : ↑iU < G.m + G.m := by omega
  simp [normalizedA, upperIndex, posIndex, Fin.addCases, Fin.castAdd, Fin.castLT, hlt]

lemma normalizedA_upper_neg (iU : Fin G.m) (jN : Fin G.n) :
    G.normalizedA (G.upperIndex iU) (G.negIndex jN) =
      if G.rel iU = ConstraintRel.ge then 0 else -(G.A iU jN) := by
  have hlt : ↑iU < G.m + G.m := by omega
  simp [normalizedA, upperIndex, negIndex, Fin.addCases, Fin.castAdd, Fin.castLT, Fin.natAdd, hlt]lemma normalizedA_lower_pos (iL : Fin G.m) (jP : Fin G.n) :
    G.normalizedA (G.lowerIndex iL) (G.posIndex jP) =
      if G.rel iL = ConstraintRel.le then 0 else -(G.A iL jP) := by
  simp [normalizedA, lowerIndex, posIndex, Fin.addCases, Fin.castAdd, Fin.castLT, Fin.natAdd]lemma normalizedA_lower_neg (iL : Fin G.m) (jN : Fin G.n) :
    G.normalizedA (G.lowerIndex iL) (G.negIndex jN) =
      if G.rel iL = ConstraintRel.le then 0 else G.A iL jN := by
  simp [normalizedA, lowerIndex, negIndex, Fin.addCases, Fin.castAdd, Fin.castLT, Fin.natAdd]

lemma normalizedA_lock_pos (jL jP : Fin G.n) :
    G.normalizedA (G.lockIndex jL) (G.posIndex jP) = 0 := by
  simp [normalizedA, lockIndex, posIndex, Fin.addCases, Fin.castAdd, This simp argument is unused:
  Fin.castLT

Hint: Omit it from the simp argument list.
  simp [normalizedA, lockIndex, posIndex, Fin.addCases, Fin.castAdd, Fin.c̵a̵s̵t̵L̵T̵,̵ ̵F̵i̵n.̵n̵atAdd]

Note: This linter can be disabled with `set_option linter.unusedSimpArgs false`Fin.castLT, Fin.natAdd]

lemma normalizedA_lock_neg (jL jN : Fin G.n) :
    G.normalizedA (G.lockIndex jL) (G.negIndex jN) =
      if G.free jL then 0 else (if jN = jL then 1 else 0) := by
  simp [normalizedA, lockIndex, negIndex, Fin.addCases, This simp argument is unused:
  Fin.castAdd

Hint: Omit it from the simp argument list.
  simp [normalizedA, lockIndex, negIndex, Fin.addCases, Fin.castA̵d̵d̵,̵ ̵F̵i̵n̵.̵c̵a̵s̵t̵LT, Fin.natAdd]

Note: This linter can be disabled with `set_option linter.unusedSimpArgs false`Fin.castAdd, This simp argument is unused:
  Fin.castLT

Hint: Omit it from the simp argument list.
  simp [normalizedA, lockIndex, negIndex, Fin.addCases, Fin.castAdd, Fin.c̵a̵s̵t̵L̵T̵,̵ ̵F̵i̵n.̵n̵atAdd]

Note: This linter can be disabled with `set_option linter.unusedSimpArgs false`Fin.castLT, Fin.natAdd]

lemma normalizedB_upper (iU : Fin G.m) :
    G.normalizedB (G.upperIndex iU) = if G.rel iU = ConstraintRel.ge then 0 else G.b iU := by
  have hlt : ↑iU < G.m + G.m := by omega
  simp [normalizedB, upperIndex, Fin.addCases, Fin.castAdd, Fin.castLT, hlt]lemma normalizedB_lower (iL : Fin G.m) :
    G.normalizedB (G.lowerIndex iL) = if G.rel iL = ConstraintRel.le then 0 else -(G.b iL) := by
  simp [normalizedB, lowerIndex, Fin.addCases, Fin.castAdd, Fin.castLT, Fin.natAdd]lemma normalizedB_lock (jL : Fin G.n) : G.normalizedB (G.lockIndex jL) = 0 := by
  simp [normalizedB, lockIndex, Fin.addCases, Fin.natAdd]lemma normalizedC_pos (jP : Fin G.n) :
    G.normalizedC (G.posIndex jP) = G.objectiveSign * G.c jP := by
  simp [normalizedC, posIndex, Fin.addCases, Fin.castAdd, Fin.castLT]lemma normalizedC_neg (jN : Fin G.n) :
    G.normalizedC (G.negIndex jN) = -(G.objectiveSign * G.c jN) := by
  simp [normalizedC, negIndex, Fin.addCases, Fin.natAdd]

Row-vector identities

              

lemma normalizedA_mulVec_upper (x' : Fin (G.n + G.n) → ℝ) (iU : Fin G.m)
    (hge : G.rel iU ≠ ConstraintRel.ge) :
    (G.normalizedA *ᵥ x') (G.upperIndex iU) =
      ∑ j : Fin G.n, G.A iU j * (x' (G.posIndex j) - x' (G.negIndex j)) := by
  simp only [Matrix.mulVec, dotProduct]
  rw [Fin.sum_univ_add]
  have h1 : (∑ j : Fin G.n, G.normalizedA (G.upperIndex iU) (G.posIndex j) * x' (G.posIndex j)) =
      ∑ j : Fin G.n, G.A iU j * x' (G.posIndex j) := by
    apply Finset.sum_congr rfl
    intro j hj
    rw [normalizedA_upper_pos, if_neg hge]
  have h2 : (∑ j : Fin G.n, G.normalizedA (G.upperIndex iU) (G.negIndex j) * x' (G.negIndex j)) =
      ∑ j : Fin G.n, (-(G.A iU j)) * x' (G.negIndex j) := by
    apply Finset.sum_congr rfl
    intro j hj
    rw [normalizedA_upper_neg, if_neg hge]
  rw [h1, h2]
  rw [← Finset.sum_add_distrib]
  apply Finset.sum_congr rfl
  intro j hj
  ring

lemma normalizedA_mulVec_upper_ge (x' : Fin (G.n + G.n) → ℝ) (iU : Fin G.m)
    (hge : G.rel iU = ConstraintRel.ge) :
    (G.normalizedA *ᵥ x') (G.upperIndex iU) = 0 := by
  simp only [Matrix.mulVec, dotProduct]
  rw [Fin.sum_univ_add]
  have h1 : (∑ j : Fin G.n, G.normalizedA (G.upperIndex iU) (G.posIndex j) * x' (G.posIndex j)) = 0 := by
    apply Finset.sum_eq_zero
    intro j hj
    rw [normalizedA_upper_pos, if_pos hge, zero_mul]
  have h2 : (∑ j : Fin G.n, G.normalizedA (G.upperIndex iU) (G.negIndex j) * x' (G.negIndex j)) = 0 := by
    apply Finset.sum_eq_zero
    intro j hj
    rw [normalizedA_upper_neg, if_pos hge, zero_mul]
  rw [h1, h2]
  simp

lemma normalizedA_mulVec_lower (x' : Fin (G.n + G.n) → ℝ) (iL : Fin G.m)
    (hle : G.rel iL ≠ ConstraintRel.le) :
    (G.normalizedA *ᵥ x') (G.lowerIndex iL) =
      -(∑ j : Fin G.n, G.A iL j * (x' (G.posIndex j) - x' (G.negIndex j))) := by
  simp only [Matrix.mulVec, dotProduct]
  rw [Fin.sum_univ_add]
  have h1 : (∑ j : Fin G.n, G.normalizedA (G.lowerIndex iL) (G.posIndex j) * x' (G.posIndex j)) =
      ∑ j : Fin G.n, -(G.A iL j * x' (G.posIndex j)) := by
    apply Finset.sum_congr rfl
    intro j hj
    rw [normalizedA_lower_pos, if_neg hle, neg_mul]
  have h2 : (∑ j : Fin G.n, G.normalizedA (G.lowerIndex iL) (G.negIndex j) * x' (G.negIndex j)) =
      ∑ j : Fin G.n, G.A iL j * x' (G.negIndex j) := by
    apply Finset.sum_congr rfl
    intro j hj
    rw [normalizedA_lower_neg, if_neg hle]
  rw [h1, h2]
  rw [Finset.sum_neg_distrib]
  simp only [mul_sub]
  rw [Finset.sum_sub_distrib]
  ring

lemma normalizedA_mulVec_lower_le (x' : Fin (G.n + G.n) → ℝ) (iL : Fin G.m)
    (hle : G.rel iL = ConstraintRel.le) :
    (G.normalizedA *ᵥ x') (G.lowerIndex iL) = 0 := by
  simp only [Matrix.mulVec, dotProduct]
  rw [Fin.sum_univ_add]
  have h1 : (∑ j : Fin G.n, G.normalizedA (G.lowerIndex iL) (G.posIndex j) * x' (G.posIndex j)) = 0 := by
    apply Finset.sum_eq_zero
    intro j hj
    rw [normalizedA_lower_pos, if_pos hle, zero_mul]
  have h2 : (∑ j : Fin G.n, G.normalizedA (G.lowerIndex iL) (G.negIndex j) * x' (G.negIndex j)) = 0 := by
    apply Finset.sum_eq_zero
    intro j hj
    rw [normalizedA_lower_neg, if_pos hle, zero_mul]
  rw [h1, h2]
  simp

lemma normalizedA_mulVec_lock (x' : Fin (G.n + G.n) → ℝ) (jL : Fin G.n) :
    (G.normalizedA *ᵥ x') (G.lockIndex jL) = if G.free jL then 0 else x' (G.negIndex jL) := by
  simp only [Matrix.mulVec, dotProduct]
  rw [Fin.sum_univ_add]
  have h1 : (∑ j : Fin G.n, G.normalizedA (G.lockIndex jL) (G.posIndex j) * x' (G.posIndex j)) = 0 := by
    apply Finset.sum_eq_zero
    intro j hj
    rw [normalizedA_lock_pos, zero_mul]
  rw [h1]
  simp only [zero_add]
  by_cases hfree : G.free jL
  · rw [if_pos hfree]
    apply Finset.sum_eq_zero
    intro j hj
    rw [normalizedA_lock_neg, if_pos hfree, zero_mul]
  · rw [if_neg hfree]
    have hsum : (∑ j : Fin G.n, (if j = jL then 1 else 0) * x' (G.negIndex j)) = x' (G.negIndex jL) := by
      rw [Finset.sum_eq_single jL]
      · simp
      · intro j hj hne
        simp [hne]
      · intro h
        exact False.elim (h (Finset.mem_univ jL))
    calc
      (∑ j : Fin G.n, G.normalizedA (G.lockIndex jL) (G.negIndex j) * x' (G.negIndex j))
          = ∑ j : Fin G.n, (if G.free jL then 0 else (if j = jL then 1 else 0)) * x' (G.negIndex j) := by
              apply Finset.sum_congr rfl
              intro j hj
              rw [normalizedA_lock_neg]
      _ = ∑ j : Fin G.n, (if j = jL then 1 else 0) * x' (G.negIndex j) := by
              simp [hfree]
      _ = x' (G.negIndex jL) := hsum

Feasibility

An expanded assignment is coordinatewise nonnegative.


theorem lift_nonnegative (x : Fin G.n → ℝ) : IsNonnegative (G.lift x) := by
  intro j
  exact Fin.addCases (motive := fun j => 0 ≤ G.lift x j)
    (fun jP => by rw [G.lift_pos]; exact le_max_right _ _)
    (fun jN => by rw [G.lift_neg]; exact le_max_right _ _)
    j

The upper form of a constraint holds for the expanded assignment.


lemma feasible_constraint_upper {x : Fin G.n → ℝ} (hx : G.IsFeasible x) (iU : Fin G.m) :
    (G.normalizedA *ᵥ G.lift x) (G.upperIndex iU) ≤ G.normalizedB (G.upperIndex iU) := by
  by_cases hge : G.rel iU = ConstraintRel.ge
  · simp [G.normalizedA_mulVec_upper_ge (G.lift x) iU hge, G.normalizedB_upper, hge]
  · rw [G.normalizedA_mulVec_upper (G.lift x) iU hge]
    rw [G.normalizedB_upper]
    simp only [hge, if_false]
    have hsum : (∑ j : Fin G.n, G.A iU j * (G.lift x (G.posIndex j) - G.lift x (G.negIndex j))) =
        (G.A *ᵥ x) iU := by
      simp only [Matrix.mulVec, dotProduct]
      apply Finset.sum_congr rfl
      intro j hj
      rw [G.lift_pos, G.lift_neg, max_sub_max_neg]
    rw [hsum]
    cases hrel : G.rel iU with
    | le => exact by simpa [hrel] using hx.2 iU
    | eq => exact le_of_eq (by simpa [hrel] using hx.2 iU)
    | ge => exact False.elim (hge hrel)

The lower form of a constraint holds for the expanded assignment.


lemma feasible_constraint_lower {x : Fin G.n → ℝ} (hx : G.IsFeasible x) (iL : Fin G.m) :
    (G.normalizedA *ᵥ G.lift x) (G.lowerIndex iL) ≤ G.normalizedB (G.lowerIndex iL) := by
  by_cases hle : G.rel iL = ConstraintRel.le
  · simp [G.normalizedA_mulVec_lower_le (G.lift x) iL hle, G.normalizedB_lower, hle]
  · rw [G.normalizedA_mulVec_lower (G.lift x) iL hle]
    rw [G.normalizedB_lower]
    simp only [hle, if_false]
    have hsum : (∑ j : Fin G.n, G.A iL j * (G.lift x (G.posIndex j) - G.lift x (G.negIndex j))) =
        (G.A *ᵥ x) iL := by
      simp only [Matrix.mulVec, dotProduct]
      apply Finset.sum_congr rfl
      intro j hj
      rw [G.lift_pos, G.lift_neg, max_sub_max_neg]
    rw [hsum]
    cases hrel : G.rel iL with
    | le => exact False.elim (hle hrel)
    | eq => exact le_of_eq (by simpa [hrel] using hx.2 iL)
    | ge => exact neg_le_neg (by simpa [hrel] using hx.2 iL)

The nonnegativity-lock row holds for the expanded assignment.


lemma feasible_constraint_lock {x : Fin G.n → ℝ} (hx : G.IsFeasible x) (jL : Fin G.n) :
    (G.normalizedA *ᵥ G.lift x) (G.lockIndex jL) ≤ G.normalizedB (G.lockIndex jL) := by
  rw [G.normalizedA_mulVec_lock, G.normalizedB_lock]
  by_cases hfree : G.free jL
  · simp [hfree]
  · rw [if_neg hfree]
    have hnn : 0 ≤ x jL := hx.1 jL hfree
    rw [G.lift_neg]
    have hneg : -(x jL) ≤ 0 := by linarith
    rw [max_eq_right hneg]

A feasible assignment expands to a feasible standard-form assignment.

theorem normalized_feasible_of_feasible {x : Fin G.n → ℝ} (hx : G.IsFeasible x) :
    (G.toStandardLP).IsFeasible (G.lift x) := by
  refine ⟨G.lift_nonnegative x, ?_⟩
  intro i
  exact Fin.addCases (motive := fun i => (G.normalizedA *ᵥ G.lift x) i ≤ G.normalizedB i)
    (fun ij =>
      Fin.addCases (motive := fun ij => (G.normalizedA *ᵥ G.lift x) (Fin.castAdd G.n ij) ≤ G.normalizedB (Fin.castAdd G.n ij))
        (fun iU => G.feasible_constraint_upper hx iU)
        (fun iL => G.feasible_constraint_lower hx iL)
        ij)
    (fun jL => G.feasible_constraint_lock hx jL)
    i

Every standard-form-feasible expanded assignment collapses to a feasible general-form assignment.


theorem feasible_of_normalized_feasible {x' : Fin (G.n + G.n) → ℝ}
    (hx' : (G.toStandardLP).IsFeasible x') : G.IsFeasible (G.proj x') := by
  refine ⟨?_, ?_⟩
  · intro j hjfree
    have hlock : (G.normalizedA *ᵥ x') (G.lockIndex j) ≤ G.normalizedB (G.lockIndex j) :=
      hx'.2 (G.lockIndex j)
    rw [G.normalizedA_mulVec_lock, G.normalizedB_lock] at hlock
    by_cases hfree : G.free j
    · exfalso
      exact hjfree hfree
    · rw [if_neg hfree] at hlock
      have hneg_nonneg := hx'.1 (G.negIndex j)
      have hneg_zero : x' (G.negIndex j) = 0 := le_antisymm hlock hneg_nonneg
      unfold proj
      rw [hneg_zero]
      simpa using hx'.1 (G.posIndex j)
  · intro i
    cases hrel : G.rel i with
    | le =>
        have hne : G.rel i ≠ ConstraintRel.ge := by
          intro hc
          rw [hrel] at hc
          cases hc
        have hupper : (G.normalizedA *ᵥ x') (G.upperIndex i) ≤ G.normalizedB (G.upperIndex i) :=
          hx'.2 (G.upperIndex i)
        rw [G.normalizedA_mulVec_upper x' i hne, G.normalizedB_upper] at hupper
        simp only [hrel, This simp argument is unused:
  if_false

Hint: Omit it from the simp argument list.
  simp only [hrel,̵ ̵i̵f̵_̵f̵a̵l̵s̵e̵] at hupper

Note: This linter can be disabled with `set_option linter.unusedSimpArgs false`if_false] at hupper
        simpa [Matrix.mulVec, dotProduct, proj] using hupper
    | eq =>
        have hne_ge : G.rel i ≠ ConstraintRel.ge := by
          intro hc
          rw [hrel] at hc
          cases hc
        have hne_le : G.rel i ≠ ConstraintRel.le := by
          intro hc
          rw [hrel] at hc
          cases hc
        have hupper : (G.normalizedA *ᵥ x') (G.upperIndex i) ≤ G.normalizedB (G.upperIndex i) :=
          hx'.2 (G.upperIndex i)
        have hlower : (G.normalizedA *ᵥ x') (G.lowerIndex i) ≤ G.normalizedB (G.lowerIndex i) :=
          hx'.2 (G.lowerIndex i)
        rw [G.normalizedA_mulVec_upper x' i hne_ge, G.normalizedB_upper] at hupper
        rw [G.normalizedA_mulVec_lower x' i hne_le, G.normalizedB_lower] at hlower
        simp only [hrel, This simp argument is unused:
  if_false

Hint: Omit it from the simp argument list.
  simp only [hrel,̵ ̵i̵f̵_̵f̵a̵l̵s̵e̵] at hupper hlower

Note: This linter can be disabled with `set_option linter.unusedSimpArgs false`if_false] at hupper hlower
        have hle : (G.A *ᵥ G.proj x') i ≤ G.b i := by
          simpa [Matrix.mulVec, dotProduct, proj] using hupper
        have hge : G.b i ≤ (G.A *ᵥ G.proj x') i := by
          have h := neg_le_neg_iff.mp hlower
          simpa [Matrix.mulVec, dotProduct, proj] using h
        exact le_antisymm hle hge
    | ge =>
        have hne : G.rel i ≠ ConstraintRel.le := by
          intro hc
          rw [hrel] at hc
          cases hc
        have hlower : (G.normalizedA *ᵥ x') (G.lowerIndex i) ≤ G.normalizedB (G.lowerIndex i) :=
          hx'.2 (G.lowerIndex i)
        rw [G.normalizedA_mulVec_lower x' i hne, G.normalizedB_lower] at hlower
        simp only [hrel, This simp argument is unused:
  if_false

Hint: Omit it from the simp argument list.
  simp only [hrel,̵ ̵i̵f̵_̵f̵a̵l̵s̵e̵] at hlower

Note: This linter can be disabled with `set_option linter.unusedSimpArgs false`if_false] at hlower
        have h := neg_le_neg_iff.mp hlower
        simpa [Matrix.mulVec, dotProduct, proj] using h

The expansion collapses to the original assignment.


theorem proj_lift (x : Fin G.n → ℝ) : G.proj (G.lift x) = x := by
  funext j
  unfold proj
  rw [G.lift_pos, G.lift_neg]
  exact max_sub_max_neg (x j)

Feasibility is exactly preserved under the expansion.

theorem feasible_iff_lift (x : Fin G.n → ℝ) :
    G.IsFeasible x ↔ (G.toStandardLP).IsFeasible (G.lift x) := by
  constructor
  · exact G.normalized_feasible_of_feasible
  · intro h
    simpa [G.proj_lift x] using G.feasible_of_normalized_feasible h

Existence of a feasible assignment is preserved under the normalization.

theorem feasible_iff_exists :
    (∃ x : Fin G.n → ℝ, G.IsFeasible x) ↔
      (∃ x' : Fin (G.n + G.n) → ℝ, (G.toStandardLP).IsFeasible x') := by
  constructor
  · rintro ⟨x, hx⟩
    exact ⟨G.lift x, G.normalized_feasible_of_feasible hx⟩
  · rintro ⟨x', hx'⟩
    exact ⟨G.proj x', G.feasible_of_normalized_feasible hx'⟩

Objective preservation

The normalized objective of an arbitrary expanded assignment is the signed objective of its collapse.


theorem objective_eq_sign_objective_proj (x' : Fin (G.n + G.n) → ℝ) :
    (G.toStandardLP).objective x' = G.objectiveSign * G.objective (G.proj x') := by
  simp only [StandardLP.objective, objective, proj, dotProduct]
  rw [Fin.sum_univ_add]
  simp only [normalizedC_pos, normalizedC_neg]
  calc
    (∑ jP : Fin G.n, (G.objectiveSign * G.c jP) * x' (G.posIndex jP))
        + ∑ jN : Fin G.n, (-(G.objectiveSign * G.c jN)) * x' (G.negIndex jN)
        = (∑ jP : Fin G.n, G.objectiveSign * (G.c jP * x' (G.posIndex jP)))
          + ∑ jN : Fin G.n, (-(G.objectiveSign)) * (G.c jN * x' (G.negIndex jN)) := by
            congr 1
            · apply Finset.sum_congr rfl
              intro j hj
              ring
            · apply Finset.sum_congr rfl
              intro j hj
              ring
    _ = G.objectiveSign * (∑ jP : Fin G.n, G.c jP * x' (G.posIndex jP))
          + (-(G.objectiveSign)) * (∑ jN : Fin G.n, G.c jN * x' (G.negIndex jN)) := by
            rw [← Finset.mul_sum, ← Finset.mul_sum]
    _ = G.objectiveSign * (∑ j : Fin G.n, G.c j * (x' (G.posIndex j) - x' (G.negIndex j))) := by
            simp only [mul_sub]
            rw [Finset.sum_sub_distrib]
            ring

The normalized objective of an expanded assignment is the signed general objective.


theorem objective_lift (x : Fin G.n → ℝ) :
    (G.toStandardLP).objective (G.lift x) = G.objectiveSign * G.objective x := by
  rw [G.objective_eq_sign_objective_proj]
  congr 1
  rw [G.proj_lift x]

              end GeneralLPend Chapter29end CLRS

CLRSLean.FourthEdition.Chapter_29.Section_29_1_Standard_And_Slack_Forms.SolverWrapper

29.1 Canonical main-text solver wrapper

This module combines the general-form normalization of the previous module with the initialized SIMPLEX solver (online material) to obtain a canonical main-text solver for general-form programs. The result type exposes exactly the three textbook outcomes.

Main results:

GeneralLP.Result: infeasible, optimal, or unbounded.
GeneralLP.solve: normalize, solve, and translate back.
GeneralLP.solve_complete: the wrapper always certifies one outcome.

              namespace CLRSnamespace Chapter29namespace GeneralLPvariable (G : GeneralLP)

Certified outcomes of the normalized solver, mapped back to the general program.

inductive Result where
  | infeasible : (¬ ∃ x : Fin G.n → ℝ, G.IsFeasible x) → Result
  | optimal : (x : Fin G.n → ℝ) → G.IsOptimal x → Result
  | unbounded : G.IsUnbounded → Result

A canonical main-text solver wrapper: normalize, run the initialized SIMPLEX, and translate the outcome back to the general program.


noncomputable def solve : Result G := by
  classical
  match G.toStandardLP.initializedSimplex with
  | StandardLP.InitializedSimplexResult.infeasible h =>
      exact .infeasible (by
        intro hx
        exact h (G.feasible_iff_exists.mp hx))
  | StandardLP.InitializedSimplexResult.optimal x' hx' =>
      exact .optimal (G.proj x') (by
        refine ⟨G.feasible_of_normalized_feasible hx'.1, ?_⟩
        intro z hz
        have hz' : (G.toStandardLP).IsFeasible (G.lift z) := G.normalized_feasible_of_feasible hz
        have hbound := hx'.2 (G.lift z) hz'
        rw [G.objective_lift z, G.objective_eq_sign_objective_proj] at hbound
        exact hbound)
  | StandardLP.InitializedSimplexResult.unbounded h =>
      exact .unbounded (by
        intro M
        obtain ⟨x', hx', hM⟩ := h M
        refine ⟨G.proj x', G.feasible_of_normalized_feasible hx', ?_⟩
        rw [G.objective_eq_sign_objective_proj] at hM
        exact hM)

The wrapper always certifies infeasibility, an optimum, or unboundedness.

theorem solve_complete :
    (¬ ∃ x : Fin G.n → ℝ, G.IsFeasible x) ∨
      (∃ x : Fin G.n → ℝ, G.IsOptimal x) ∨ G.IsUnbounded := by
  cases G.solve with
  | infeasible h => exact Or.inl h
  | optimal x hx => exact Or.inr (Or.inl ⟨x, hx⟩)
  | unbounded h => exact Or.inr (Or.inr h)

              end GeneralLPend Chapter29end CLRS

Imports

import CLRSLean.FourthEdition.Chapter_29.Section_29_2_Formulating_Problems_As_Linear_Programs.NetworkFlow
import CLRSLean.FourthEdition.Chapter_29.Section_29_2_Formulating_Problems_As_Linear_Programs.ShortestPath
import CLRSLean.FourthEdition.Chapter_29.Section_29_2_Formulating_Problems_As_Linear_Programs.MaximumFlow
import CLRSLean.FourthEdition.Chapter_29.Section_29_2_Formulating_Problems_As_Linear_Programs.MinimumCostFlow
import CLRSLean.FourthEdition.Chapter_29.Section_29_2_Formulating_Problems_As_Linear_Programs.MulticommodityFlow

29.2. Formulating Problems as Linear Programs

This section contains the four textbook formulations: shortest path, maximum flow, minimum-cost flow, and multicommodity flow.

Implementation details

              namespace CLRSnamespace Chapter29end Chapter29end CLRS

Definitions and proofs

CLRSLean.FourthEdition.Chapter_29.Section_29_2_Formulating_Problems_As_Linear_Programs.NetworkFlow

29.2: Common finite-network definitions

CLRS writes its flow linear programs with one nonnegative variable f u v for every ordered pair of vertices. A missing edge is represented by capacity zero. This file records that shared finite-network vocabulary, together with the finite-standard-form encoding helpers used by every §29.2 formulation.

              namespace CLRSnamespace Chapter29open Finsetopen scoped BigOperators

A finite directed capacitated network. Capacities are total; assigning capacity zero to nonedges gives exactly the convention used in CLRS §29.2.

structure FlowNetwork (V : Type*) [Fintype V] where
  source : V
  sink : V
  source_ne_sink : source ≠ sink
  capacity : V → V → ℝ
  capacity_nonnegative : ∀ u v, 0 ≤ capacity u v

namespace FlowNetworkvariable {V : Type*} [Fintype V] (N : FlowNetwork V)

Total flow leaving a vertex.

def outflow (f : V → V → ℝ) (u : V) : ℝ := ∑ v, f u v

Total flow entering a vertex.

def inflow (f : V → V → ℝ) (u : V) : ℝ := ∑ v, f v u

Net flow leaving a vertex.

def netOutflow (f : V → V → ℝ) (u : V) : ℝ := outflow f u - inflow f u

Flow conservation at one vertex.

def ConservesAt (f : V → V → ℝ) (u : V) : Prop := inflow f u = outflow f u

The nonnegativity and capacity inequalities shared by the flow LPs.

def IsCapacityFeasible (f : V → V → ℝ) : Prop :=
  (∀ u v, 0 ≤ f u v) ∧ ∀ u v, f u v ≤ N.capacity u v

A feasible single-commodity flow: capacity constraints everywhere and conservation away from the source and sink.

def IsFlow (f : V → V → ℝ) : Prop :=
  N.IsCapacityFeasible f ∧
    ∀ u, u ≠ N.source → u ≠ N.sink → ConservesAt f u

end FlowNetwork

Finite standard-form encoding helpers

These helpers reindex semantic vectors over a finite type into the Fin-indexed vectors used by StandardLP, and prove the sum identities that the assignment bridges rely on.

namespace FinEncodingvariable {ι : Type*} [Fintype ι]

Lift a semantic vector over a finite type to the Fin-indexed vector of the same coordinates.

noncomputable def lift (x : ι → ℝ) : Fin (Fintype.card ι) → ℝ :=
  x ∘ (Fintype.equivFin ι).symm

Project a Fin-indexed vector back to a semantic vector.

noncomputable def proj (x : Fin (Fintype.card ι) → ℝ) : ι → ℝ :=
  x ∘ (Fintype.equivFin ι)

Lifting then projecting recovers the original vector.

theorem proj_lift (x : ι → ℝ) (a : ι) : proj (lift x) a = x a := by
  simp [proj, lift]

A finite sum over Fin (card ι) reindexes to a sum over ι.

theorem sum_reindex {M : Type*} [AddCommMonoid M] (f : Fin (Fintype.card ι) → M) :
    (∑ i : Fin (Fintype.card ι), f i) = ∑ a : ι, f (Fintype.equivFin ι a) := by
  exact (Equiv.sum_comp (Fintype.equivFin ι) f).symm

The indicator sum that collapses an ι-indicator to its selected coordinate.


theorem sum_indicator [DecidableEq ι] (a : ι) (x : ι → ℝ) :
    (∑ j : Fin (Fintype.card ι),
        (if (Fintype.equivFin ι).symm j = a then 1 else 0) * lift x j) = x a := by
  calc
    (∑ j : Fin (Fintype.card ι),
        (if (Fintype.equivFin ι).symm j = a then 1 else 0) * lift x j)
        = ∑ b : ι, (if b = a then 1 else 0) * x b := by
            rw [sum_reindex]
            simp [lift]
    _ = x a := by
            rw [Finset.sum_eq_single a]
            · simp
            · intro b _ hb; simp [hb]
            · intro h; exact False.elim (h (Finset.mem_univ a))

A direct Fin-indicator sum (no reindexing): the coefficient is nonzero exactly at the selected coordinate.


theorem sum_indicator_fin {n : ℕ} (e : Fin n) (x : Fin n → ℝ) :
    (∑ j : Fin n, (if j = e then 1 else 0) * x j) = x e := by
  rw [Finset.sum_eq_single e]
  · simp
  · intro b _ hb; simp [hb]
  · intro h; exact False.elim (h (Finset.mem_univ e))

Fin.addCases reduces on a left (castAdd) index.

theorem addCases_castAdd {m n : ℕ} {C : Type u} (left : Fin m → C) (right : Fin n → C) (e : Fin m) :
    Fin.addCases left right (Fin.castAdd n e) = left e := by
  simp [Fin.addCases]

Fin.addCases reduces on a right (natAdd) index.

theorem addCases_natAdd {m n : ℕ} {C : Type u} (left : Fin m → C) (right : Fin n → C) (v : Fin n) :
    Fin.addCases left right (Fin.natAdd m v) = right v := by
  simp [Fin.addCases]

Lift a binary vector to the Fin-indexed vector over ordered pairs.

noncomputable def lift₂ {α β : Type*} [Fintype α] [Fintype β] (x : α → β → ℝ) :
    Fin (Fintype.card (α × β)) → ℝ :=
  lift (Function.uncurry x)

Project a Fin-indexed vector over ordered pairs back to a binary vector.

noncomputable def proj₂ {α β : Type*} [Fintype α] [Fintype β]
    (x : Fin (Fintype.card (α × β)) → ℝ) : α → β → ℝ :=
  Function.curry (proj x)

Sum of an indicator on the first coordinate over all ordered pairs.


theorem sum_indicator_fst {α β : Type*} [Fintype α] [Fintype β] [DecidableEq α]
    (a : α) (x : α → β → ℝ) :
    (∑ j : Fin (Fintype.card (α × β)),
        (if ((Fintype.equivFin (α × β)).symm j).1 = a then 1 else 0) * lift₂ x j) =
      ∑ b : β, x a b := by
  calc
    (∑ j : Fin (Fintype.card (α × β)),
        (if ((Fintype.equivFin (α × β)).symm j).1 = a then 1 else 0) * lift₂ x j)
        = ∑ p : α × β, (if p.1 = a then 1 else 0) * x p.1 p.2 := by
            rw [sum_reindex]
            simp [lift₂, lift, Function.uncurry]
    _ = ∑ b : β, x a b := by
            rw [Fintype.sum_prod_type]
            rw [Finset.sum_eq_single a]
            · simp
            · intro c _ hc; simp [hc]
            · intro h; exact False.elim (h (Finset.mem_univ a))

Sum of an indicator on the second coordinate over all ordered pairs.


theorem sum_indicator_snd {α β : Type*} [Fintype α] [Fintype β] [DecidableEq β]
    (a : β) (x : α → β → ℝ) :
    (∑ j : Fin (Fintype.card (α × β)),
        (if ((Fintype.equivFin (α × β)).symm j).2 = a then 1 else 0) * lift₂ x j) =
      ∑ b : α, x b a := by
  calc
    (∑ j : Fin (Fintype.card (α × β)),
        (if ((Fintype.equivFin (α × β)).symm j).2 = a then 1 else 0) * lift₂ x j)
        = ∑ p : α × β, (if p.2 = a then 1 else 0) * x p.1 p.2 := by
            rw [sum_reindex]
            simp [lift₂, lift, Function.uncurry]
    _ = ∑ b : α, x b a := by
            rw [Fintype.sum_prod_type]
            rw [Finset.sum_comm]
            rw [Finset.sum_eq_single a]
            · simp
            · intro c _ hc; simp [hc]
            · intro h; exact False.elim (h (Finset.mem_univ a))

              end FinEncodingend Chapter29end CLRS

CLRSLean.FourthEdition.Chapter_29.Section_29_2_Formulating_Problems_As_Linear_Programs.ShortestPath

29.2: Shortest paths as a linear program

For a source s and target t, CLRS maximizes d t, subject to d s = 0 and d v ≤ d u + w u v on every edge. Thus every feasible d t is a lower bound on every s-to-t walk; an attained bound is optimal.

This module also records the finite standard-form encoding: the general-form program is reindexed into a concrete StandardLP (one variable per vertex), and feasibility and the objective value are proved to be preserved.

              namespace CLRSnamespace Chapter29namespace ShortestPathLPopen Chapter24open Matrixopen scoped BigOperatorsvariable {V : Type*} [Fintype V] [DecidableEq V]

Feasibility constraints of the CLRS shortest-path LP.

def IsFeasible (G : WeightedGraph V) (s : V) (d : V → ℝ) : Prop :=
  d s = 0 ∧ ∀ u v, (u, v) ∈ G.edges → d v ≤ d u + G.w u v

Optimality for the maximization objective d t.

def IsOptimal (G : WeightedGraph V) (s t : V) (d : V → ℝ) : Prop :=
  IsFeasible G s d ∧ ∀ e, IsFeasible G s e → e t ≤ d t

Every feasible potential is a lower bound on the weight of every walk from the source. This is the core correctness property of the formulation.

theorem feasible_le_walkWeight {G : WeightedGraph V} {s t : V} {d : V → ℝ}
    (hd : IsFeasible G s d) (p : List V) (hp : G.IsWalkFrom s t p) :
    d t ≤ Chapter24.WeightedGraph.walkWeight G.w p := by
  have h := Chapter24.WeightedGraph.le_add_walkWeight_of_potential G d hd.2 s t p hp
  simpa [hd.1] using h

If a feasible potential is attained by an actual source-to-target walk, then it solves the shortest-path LP.

theorem optimal_of_attained_walk {G : WeightedGraph V} {s t : V} {d : V → ℝ}
    (hd : IsFeasible G s d) (p : List V) (hp : G.IsWalkFrom s t p)
    (hweight : Chapter24.WeightedGraph.walkWeight G.w p = d t) : IsOptimal G s t d := by
  refine ⟨hd, ?_⟩
  intro e he
  calc
    e t ≤ Chapter24.WeightedGraph.walkWeight G.w p := feasible_le_walkWeight he p hp
    _ = d t := hweight

Finite standard-form encoding

The type of directed edges that are actually present in the graph.

abbrev Edge (G : WeightedGraph V) := {e : V × V // e ∈ G.edges}

The general-form program for the shortest-path LP: maximize d t subject to d s = 0 and d v - d u ≤ w u v on every edge. All vertices are free variables.

noncomputable def toGeneralLP (G : WeightedGraph V) (s t : V) : GeneralLP where
  n := Fintype.card V
  m := Fintype.card (Edge G) + 1
  maximize := true
  c := fun j => if (Fintype.equivFin V).symm j = t then 1 else 0
  rel := fun i =>
    Fin.cases
      (ConstraintRel.eq)
      (fun _ => ConstraintRel.le)
      i
  A := fun i =>
    Fin.cases
      (fun j => if (Fintype.equivFin V).symm j = s then 1 else 0)
      (fun e j =>
        let uv := (Fintype.equivFin (Edge G)).symm e
        (if (Fintype.equivFin V).symm j = uv.1.2 then 1 else 0) -
          (if (Fintype.equivFin V).symm j = uv.1.1 then 1 else 0))
      i
  b := fun i =>
    Fin.cases
      (0 : ℝ)
      (fun e =>
        let uv := (Fintype.equivFin (Edge G)).symm e
        G.w uv.1.1 uv.1.2)
      i
  free := fun _ => true

The index of the source-equality constraint.

abbrev sourceIndex (G : WeightedGraph V) : Fin (Fintype.card (Edge G) + 1) := 0

The index of the constraint attached to a present edge.

abbrev edgeIndex (G : WeightedGraph V) (e : Fin (Fintype.card (Edge G))) :
    Fin (Fintype.card (Edge G) + 1) :=
  e.succ

              lemma A_source (G : WeightedGraph V) (s t : V) (j : Fin (Fintype.card V)) :
    (toGeneralLP G s t).A (sourceIndex G) j =
      if (Fintype.equivFin V).symm j = s then 1 else 0 := by
  simp [toGeneralLP, sourceIndex]lemma b_source (G : WeightedGraph V) (s t : V) :
    (toGeneralLP G s t).b (sourceIndex G) = 0 := by
  simp [toGeneralLP, sourceIndex]lemma rel_source (G : WeightedGraph V) (s t : V) :
    (toGeneralLP G s t).rel (sourceIndex G) = ConstraintRel.eq := by
  simp [toGeneralLP, sourceIndex]lemma A_edge (G : WeightedGraph V) (s t : V) (e : Fin (Fintype.card (Edge G)))
    (j : Fin (Fintype.card V)) :
    (toGeneralLP G s t).A (edgeIndex G e) j =
      (if (Fintype.equivFin V).symm j = ((Fintype.equivFin (Edge G)).symm e).1.2 then 1 else 0) -
        (if (Fintype.equivFin V).symm j = ((Fintype.equivFin (Edge G)).symm e).1.1 then 1 else 0) := by
  simp [toGeneralLP, edgeIndex]lemma b_edge (G : WeightedGraph V) (s t : V) (e : Fin (Fintype.card (Edge G))) :
    (toGeneralLP G s t).b (edgeIndex G e) =
      G.w ((Fintype.equivFin (Edge G)).symm e).1.1 ((Fintype.equivFin (Edge G)).symm e).1.2 := by
  simp [toGeneralLP, edgeIndex]lemma rel_edge (G : WeightedGraph V) (s t : V) (e : Fin (Fintype.card (Edge G))) :
    (toGeneralLP G s t).rel (edgeIndex G e) = ConstraintRel.le := by
  simp [toGeneralLP, edgeIndex]

The source row of the encoded constraint matrix reads out d s.

lemma mulVec_source (G : WeightedGraph V) (s t : V) (d : V → ℝ) :
    ((toGeneralLP G s t).A *ᵥ FinEncoding.lift d) (sourceIndex G) = d s := by
  calc
    ((toGeneralLP G s t).A *ᵥ FinEncoding.lift d) (sourceIndex G)
        = ∑ j : Fin (Fintype.card V), (toGeneralLP G s t).A (sourceIndex G) j * FinEncoding.lift d j := by rfl
    _ = ∑ j : Fin (Fintype.card V),
          (if (Fintype.equivFin V).symm j = s then 1 else 0) * FinEncoding.lift d j := by
          simp only [A_source]
    _ = d s := FinEncoding.sum_indicator s d

An edge row of the encoded constraint matrix reads out d v - d u.


lemma mulVec_edge (G : WeightedGraph V) (s t : V) (e : Fin (Fintype.card (Edge G)))
    (d : V → ℝ) :
    ((toGeneralLP G s t).A *ᵥ FinEncoding.lift d) (edgeIndex G e) =
      d ((Fintype.equivFin (Edge G)).symm e).1.2 - d ((Fintype.equivFin (Edge G)).symm e).1.1 := by
  calc
    ((toGeneralLP G s t).A *ᵥ FinEncoding.lift d) (edgeIndex G e)
        = ∑ j : Fin (Fintype.card V), (toGeneralLP G s t).A (edgeIndex G e) j * FinEncoding.lift d j := by rfl
    _ = ∑ j : Fin (Fintype.card V),
          ((if (Fintype.equivFin V).symm j = ((Fintype.equivFin (Edge G)).symm e).1.2 then 1 else 0) -
             (if (Fintype.equivFin V).symm j = ((Fintype.equivFin (Edge G)).symm e).1.1 then 1 else 0)) *
            FinEncoding.lift d j := by
          simp only [A_edge]
    _ = (∑ j : Fin (Fintype.card V),
          (if (Fintype.equivFin V).symm j = ((Fintype.equivFin (Edge G)).symm e).1.2 then 1 else 0) * FinEncoding.lift d j) -
          (∑ j : Fin (Fintype.card V),
          (if (Fintype.equivFin V).symm j = ((Fintype.equivFin (Edge G)).symm e).1.1 then 1 else 0) * FinEncoding.lift d j) := by
          simp only [sub_mul, Finset.sum_sub_distrib]
    _ = d ((Fintype.equivFin (Edge G)).symm e).1.2 - d ((Fintype.equivFin (Edge G)).symm e).1.1 := by
          rw [FinEncoding.sum_indicator, FinEncoding.sum_indicator]

Semantic feasibility is exactly encoded general-form feasibility.

Try `simp at hfree` instead of `simpa using hfree`

Note: This linter can be disabled with `set_option linter.unnecessarySimpa false`
theorem feasible_iff_generalLP (G : WeightedGraph V) (s t : V) (d : V → ℝ) :
    IsFeasible G s d ↔ (toGeneralLP G s t).IsFeasible (FinEncoding.lift d) := by
  unfold IsFeasible GeneralLP.IsFeasible
  constructor
  · intro hd
    refine ⟨?_, ?_⟩
    · intro j hfree
      exfalso
      Try `simp at hfree` instead of `simpa using hfree`

Note: This linter can be disabled with `set_option linter.unnecessarySimpa false`simpa [toGeneralLP] using hfree
    · intro i
      exact Fin.cases (motive := fun i => match (toGeneralLP G s t).rel i with
          | ConstraintRel.le => ((toGeneralLP G s t).A *ᵥ FinEncoding.lift d) i ≤ (toGeneralLP G s t).b i
          | ConstraintRel.eq => ((toGeneralLP G s t).A *ᵥ FinEncoding.lift d) i = (toGeneralLP G s t).b i
          | ConstraintRel.ge => (toGeneralLP G s t).b i ≤ ((toGeneralLP G s t).A *ᵥ FinEncoding.lift d) i)
        (by
          rw [rel_source, mulVec_source, b_source]
          exact hd.1)
        (fun e => by
          rw [rel_edge, mulVec_edge, b_edge]
          have h := hd.2 ((Fintype.equivFin (Edge G)).symm e).1.1 ((Fintype.equivFin (Edge G)).symm e).1.2 ((Fintype.equivFin (Edge G)).symm e).2
          linarith)
        i
  · intro h
    refine ⟨?_, ?_⟩
    · have hsource := h.2 (sourceIndex G)
      rw [rel_source] at hsource
      rw [mulVec_source, b_source] at hsource
      exact hsource
    · intro u v huv
      let e : Edge G := ⟨(u, v), huv⟩
      have hedge := h.2 (edgeIndex G (Fintype.equivFin (Edge G) e))
      rw [rel_edge] at hedge
      rw [mulVec_edge, b_edge] at hedge
      have hedge' : d v - d u ≤ G.w u v := by simpa [e] using hedge
      linarith

The general-form objective reads out d t.


lemma objective_generalLP (G : WeightedGraph V) (s t : V) (d : V → ℝ) :
    (toGeneralLP G s t).objective (FinEncoding.lift d) = d t := by
  dsimp [GeneralLP.objective, toGeneralLP]
  simp only [dotProduct]
  rw [FinEncoding.sum_reindex]
  simp only [FinEncoding.lift]
  rw [Finset.sum_eq_single t]
  · simp
  · intro v _ hv; simp [hv]
  · intro h; exact False.elim (h (Finset.mem_univ t))

The finite standard-form program obtained by normalizing the shortest-path general program.

noncomputable def toStandardLP (G : WeightedGraph V) (s t : V) :=
  (toGeneralLP G s t).toStandardLP

The expansion of a potential to the normalized standard-form variable vector.

noncomputable def fullLift (G : WeightedGraph V) (s t : V) (d : V → ℝ) :
    Fin (Fintype.card V + Fintype.card V) → ℝ :=
  (toGeneralLP G s t).lift (FinEncoding.lift d)

The normalized standard-form objective is exactly the shortest-path objective d t.


theorem objective_toStandardLP (G : WeightedGraph V) (s t : V) (d : V → ℝ) :
    (toStandardLP G s t).objective (fullLift G s t d) = d t := by
  unfold toStandardLP fullLift
  rw [GeneralLP.objective_lift]
  rw [objective_generalLP]
  simp [toGeneralLP, GeneralLP.objectiveSign]

Feasibility is exactly preserved under the finite standard-form encoding.


theorem feasible_iff_toStandardLP (G : WeightedGraph V) (s t : V) (d : V → ℝ) :
    IsFeasible G s d ↔ (toStandardLP G s t).IsFeasible (fullLift G s t d) := by
  unfold toStandardLP fullLift
  rw [← GeneralLP.feasible_iff_lift]
  exact feasible_iff_generalLP G s t d

              end ShortestPathLPend Chapter29end CLRS

CLRSLean.FourthEdition.Chapter_29.Section_29_2_Formulating_Problems_As_Linear_Programs.MaximumFlow

29.2: Maximum flow as a linear program

The variables are the gross directed flows f u v. The constraints are nonnegativity, capacity, and conservation at every vertex other than s,t; the objective is the net flow leaving s, exactly as in CLRS §29.2.

This module also records the finite standard-form encoding: one nonnegative variable per ordered pair, a capacity inequality per ordered pair, and a conservation equality per internal vertex, reindexed into a concrete StandardLP.

              namespace CLRSnamespace Chapter29namespace MaximumFlowLPopen Finsetopen Matrixopen scoped BigOperatorsvariable {V : Type*} [Fintype V] (N : FlowNetwork V)

The constraints of the maximum-flow LP.

def IsFeasible (f : V → V → ℝ) : Prop :=
  (∀ u v, 0 ≤ f u v) ∧
  (∀ u v, f u v ≤ N.capacity u v) ∧
  ∀ u, u ≠ N.source → u ≠ N.sink → FlowNetwork.ConservesAt f u

The textbook maximum-flow objective: source outflow minus source inflow.

def objective (f : V → V → ℝ) : ℝ := FlowNetwork.netOutflow f N.source

A feasible flow maximizing the source outflow.

def IsOptimal (f : V → V → ℝ) : Prop :=
  IsFeasible N f ∧ ∀ g, IsFeasible N g → objective N g ≤ objective N f

The displayed LP constraints are precisely the usual finite-network flow constraints.

theorem isFeasible_iff (f : V → V → ℝ) : IsFeasible N f ↔ N.IsFlow f := by
  simp only [IsFeasible, FlowNetwork.IsFlow, FlowNetwork.IsCapacityFeasible]
  tauto

The LP optimum is precisely a maximum flow under the same objective.

theorem isOptimal_iff (f : V → V → ℝ) :
    IsOptimal N f ↔
      N.IsFlow f ∧ ∀ g, N.IsFlow g →
        FlowNetwork.netOutflow g N.source ≤ FlowNetwork.netOutflow f N.source := by
  simp only [IsOptimal, objective, isFeasible_iff]

Finite standard-form encoding

section Encodingvariable [DecidableEq V]

The internal vertices at which flow is conserved.

abbrev Internal (N : FlowNetwork V) := {u : V // u ≠ N.source ∧ u ≠ N.sink}

The general-form program for the maximum-flow LP: one nonnegative variable per ordered pair, a capacity inequality per ordered pair, and a conservation equality per internal vertex.

noncomputable def toGeneralLP (N : FlowNetwork V) : GeneralLP where
  n := Fintype.card (V × V)
  m := Fintype.card (V × V) + Fintype.card (Internal N)
  maximize := true
  c := fun j =>
    (if ((Fintype.equivFin (V × V)).symm j).1 = N.source then 1 else 0) -
      (if ((Fintype.equivFin (V × V)).symm j).2 = N.source then 1 else 0)
  rel := fun i =>
    Fin.addCases (motive := fun _ => ConstraintRel)
      (fun _ => ConstraintRel.le)
      (fun _ => ConstraintRel.eq)
      i
  A := fun i =>
    Fin.addCases (motive := fun _ => Fin (Fintype.card (V × V)) → ℝ)
      (fun e j => if j = e then 1 else 0)
      (fun u j =>
        let v := ((Fintype.equivFin (Internal N)).symm u).1
        (if ((Fintype.equivFin (V × V)).symm j).1 = v then 1 else 0) -
          (if ((Fintype.equivFin (V × V)).symm j).2 = v then 1 else 0))
      i
  b := fun i =>
    Fin.addCases (motive := fun _ => ℝ)
      (fun e => N.capacity ((Fintype.equivFin (V × V)).symm e).1 ((Fintype.equivFin (V × V)).symm e).2)
      (fun _ => 0)
      i
  free := fun _ => false

The index of the capacity constraint for a directed pair.

abbrev capIndex (N : FlowNetwork V) (e : Fin (Fintype.card (V × V))) :
    Fin (Fintype.card (V × V) + Fintype.card (Internal N)) :=
  Fin.castAdd (Fintype.card (Internal N)) e

The index of the conservation constraint for an internal vertex.

abbrev consIndex (N : FlowNetwork V) (u : Fin (Fintype.card (Internal N))) :
    Fin (Fintype.card (V × V) + Fintype.card (Internal N)) :=
  Fin.natAdd (Fintype.card (V × V)) u

              

lemma A_cap (N : FlowNetwork V) (e : Fin (Fintype.card (V × V))) (j : Fin (Fintype.card (V × V))) :
    (toGeneralLP N).A (capIndex N e) j = if j = e then 1 else 0 := by
  dsimp [toGeneralLP, capIndex]
  rw [FinEncoding.addCases_castAdd]

lemma b_cap (N : FlowNetwork V) (e : Fin (Fintype.card (V × V))) :
    (toGeneralLP N).b (capIndex N e) =
      N.capacity ((Fintype.equivFin (V × V)).symm e).1 ((Fintype.equivFin (V × V)).symm e).2 := by
  dsimp [toGeneralLP, capIndex]
  rw [FinEncoding.addCases_castAdd]

lemma rel_cap (N : FlowNetwork V) (e : Fin (Fintype.card (V × V))) :
    (toGeneralLP N).rel (capIndex N e) = ConstraintRel.le := by
  dsimp [toGeneralLP, capIndex]
  rw [FinEncoding.addCases_castAdd]

lemma A_cons (N : FlowNetwork V) (u : Fin (Fintype.card (Internal N))) (j : Fin (Fintype.card (V × V))) :
    (toGeneralLP N).A (consIndex N u) j =
      (if ((Fintype.equivFin (V × V)).symm j).1 = ((Fintype.equivFin (Internal N)).symm u).1 then 1 else 0) -
        (if ((Fintype.equivFin (V × V)).symm j).2 = ((Fintype.equivFin (Internal N)).symm u).1 then 1 else 0) := by
  dsimp [toGeneralLP, consIndex]
  rw [FinEncoding.addCases_natAdd]

lemma b_cons (N : FlowNetwork V) (u : Fin (Fintype.card (Internal N))) :
    (toGeneralLP N).b (consIndex N u) = 0 := by
  dsimp [toGeneralLP, consIndex]
  rw [FinEncoding.addCases_natAdd]

lemma rel_cons (N : FlowNetwork V) (u : Fin (Fintype.card (Internal N))) :
    (toGeneralLP N).rel (consIndex N u) = ConstraintRel.eq := by
  dsimp [toGeneralLP, consIndex]
  rw [FinEncoding.addCases_natAdd]

A capacity row reads out the gross flow on its directed pair.

lemma mulVec_cap (N : FlowNetwork V) (e : Fin (Fintype.card (V × V))) (f : V → V → ℝ) :
    ((toGeneralLP N).A *ᵥ FinEncoding.lift₂ f) (capIndex N e) = FinEncoding.lift₂ f e := by
  calc
    ((toGeneralLP N).A *ᵥ FinEncoding.lift₂ f) (capIndex N e)
        = ∑ j : Fin (Fintype.card (V × V)), (toGeneralLP N).A (capIndex N e) j * FinEncoding.lift₂ f j := by rfl
    _ = ∑ j : Fin (Fintype.card (V × V)), (if j = e then 1 else 0) * FinEncoding.lift₂ f j := by
          simp only [A_cap]
    _ = FinEncoding.lift₂ f e := FinEncoding.sum_indicator_fin e (FinEncoding.lift₂ f)

A conservation row reads out the net outflow at its vertex.


lemma mulVec_cons (N : FlowNetwork V) (u : Fin (Fintype.card (Internal N))) (f : V → V → ℝ) :
    ((toGeneralLP N).A *ᵥ FinEncoding.lift₂ f) (consIndex N u) =
      FlowNetwork.outflow f ((Fintype.equivFin (Internal N)).symm u).1 -
        FlowNetwork.inflow f ((Fintype.equivFin (Internal N)).symm u).1 := by
  calc
    ((toGeneralLP N).A *ᵥ FinEncoding.lift₂ f) (consIndex N u)
        = ∑ j : Fin (Fintype.card (V × V)), (toGeneralLP N).A (consIndex N u) j * FinEncoding.lift₂ f j := by rfl
    _ = ∑ j : Fin (Fintype.card (V × V)),
          ((if ((Fintype.equivFin (V × V)).symm j).1 = ((Fintype.equivFin (Internal N)).symm u).1 then 1 else 0) -
             (if ((Fintype.equivFin (V × V)).symm j).2 = ((Fintype.equivFin (Internal N)).symm u).1 then 1 else 0)) *
            FinEncoding.lift₂ f j := by
          simp only [A_cons]
    _ = (∑ j : Fin (Fintype.card (V × V)),
          (if ((Fintype.equivFin (V × V)).symm j).1 = ((Fintype.equivFin (Internal N)).symm u).1 then 1 else 0) * FinEncoding.lift₂ f j) -
          (∑ j : Fin (Fintype.card (V × V)),
          (if ((Fintype.equivFin (V × V)).symm j).2 = ((Fintype.equivFin (Internal N)).symm u).1 then 1 else 0) * FinEncoding.lift₂ f j) := by
          simp only [sub_mul, Finset.sum_sub_distrib]
    _ = (∑ b : V, f ((Fintype.equivFin (Internal N)).symm u).1 b) -
          (∑ b : V, f b ((Fintype.equivFin (Internal N)).symm u).1) := by
          rw [FinEncoding.sum_indicator_fst, FinEncoding.sum_indicator_snd]
    _ = FlowNetwork.outflow f ((Fintype.equivFin (Internal N)).symm u).1 -
          FlowNetwork.inflow f ((Fintype.equivFin (Internal N)).symm u).1 := rfl

Semantic feasibility is exactly encoded general-form feasibility.


theorem feasible_iff_generalLP (N : FlowNetwork V) (f : V → V → ℝ) :
    IsFeasible N f ↔ (toGeneralLP N).IsFeasible (FinEncoding.lift₂ f) := by
  unfold IsFeasible GeneralLP.IsFeasible
  constructor
  · intro hd
    refine ⟨?_, ?_⟩
    · intro j hfree
      simpa [FinEncoding.lift₂, FinEncoding.lift, Function.uncurry] using
        (hd.1 ((Fintype.equivFin (V × V)).symm j).1 ((Fintype.equivFin (V × V)).symm j).2)
    · intro i
      exact Fin.addCases (motive := fun i => match (toGeneralLP N).rel i with
          | ConstraintRel.le => ((toGeneralLP N).A *ᵥ FinEncoding.lift₂ f) i ≤ (toGeneralLP N).b i
          | ConstraintRel.eq => ((toGeneralLP N).A *ᵥ FinEncoding.lift₂ f) i = (toGeneralLP N).b i
          | ConstraintRel.ge => (toGeneralLP N).b i ≤ ((toGeneralLP N).A *ᵥ FinEncoding.lift₂ f) i)
        (fun e => by
          rw [rel_cap, mulVec_cap, b_cap]
          simpa [FinEncoding.lift₂, FinEncoding.lift, Function.uncurry] using
            (hd.2.1 ((Fintype.equivFin (V × V)).symm e).1 ((Fintype.equivFin (V × V)).symm e).2))
        (fun u => by
          rw [rel_cons, mulVec_cons, b_cons]
          have hc := hd.2.2 ((Fintype.equivFin (Internal N)).symm u).1
            ((Fintype.equivFin (Internal N)).symm u).2.1
            ((Fintype.equivFin (Internal N)).symm u).2.2
          unfold FlowNetwork.ConservesAt at hc
          linarith)
        i
  · intro h
    refine ⟨?_, ?_, ?_⟩
    · intro u v
      have hn := h.1 (Fintype.equivFin (V × V) (u, v))
      simpa [FinEncoding.lift₂, FinEncoding.lift, Function.uncurry] using (hn (by simp [toGeneralLP]))
    · intro u v
      have hc := h.2 (capIndex N (Fintype.equivFin (V × V) (u, v)))
      rw [rel_cap] at hc
      rw [mulVec_cap, b_cap] at hc
      simpa [FinEncoding.lift₂, FinEncoding.lift, Function.uncurry] using hc
    · intro u hu hw
      let e : Internal N := ⟨u, ⟨hu, hw⟩⟩
      have hcons := h.2 (consIndex N (Fintype.equivFin (Internal N) e))
      rw [rel_cons] at hcons
      rw [mulVec_cons, b_cons] at hcons
      have hcons' : FlowNetwork.outflow f u - FlowNetwork.inflow f u = 0 := by
        simpa [e] using hcons
      unfold FlowNetwork.ConservesAt
      linarith

The general-form objective is the textbook maximum-flow objective.


lemma objective_generalLP (N : FlowNetwork V) (f : V → V → ℝ) :
    (toGeneralLP N).objective (FinEncoding.lift₂ f) = objective N f := by
  dsimp [GeneralLP.objective, toGeneralLP, objective, FlowNetwork.netOutflow, FlowNetwork.outflow, FlowNetwork.inflow]
  simp only [dotProduct, sub_mul, Finset.sum_sub_distrib]
  rw [FinEncoding.sum_indicator_fst, FinEncoding.sum_indicator_snd]

The finite standard-form program obtained by normalizing the maximum-flow general program.

noncomputable def toStandardLP (N : FlowNetwork V) :=
  (toGeneralLP N).toStandardLP

The expansion of a gross-flow vector to the normalized standard-form variable vector.

noncomputable def fullLift (N : FlowNetwork V) (f : V → V → ℝ) :
    Fin (Fintype.card (V × V) + Fintype.card (V × V)) → ℝ :=
  (toGeneralLP N).lift (FinEncoding.lift₂ f)

The normalized standard-form objective is the maximum-flow objective.


theorem objective_toStandardLP (N : FlowNetwork V) (f : V → V → ℝ) :
    (toStandardLP N).objective (fullLift N f) = objective N f := by
  unfold toStandardLP fullLift
  rw [GeneralLP.objective_lift]
  rw [objective_generalLP]
  simp [toGeneralLP, GeneralLP.objectiveSign]

Feasibility is exactly preserved under the finite standard-form encoding.


theorem feasible_iff_toStandardLP (N : FlowNetwork V) (f : V → V → ℝ) :
    IsFeasible N f ↔ (toStandardLP N).IsFeasible (fullLift N f) := by
  unfold toStandardLP fullLift
  rw [← GeneralLP.feasible_iff_lift]
  exact feasible_iff_generalLP N f

              end Encodingend MaximumFlowLPend Chapter29end CLRS

CLRSLean.FourthEdition.Chapter_29.Section_29_2_Formulating_Problems_As_Linear_Programs.MinimumCostFlow

29.2: Minimum-cost flow as a linear program

In addition to capacities, a network assigns a unit cost to every directed pair. A demand-d flow has net source outflow d; minimizing the sum of unit cost times flow gives the CLRS minimum-cost-flow formulation.

This module also records the finite standard-form encoding: one nonnegative variable per ordered pair, a capacity inequality per ordered pair, a conservation equality per internal vertex, and one source-demand equality, reindexed into a concrete StandardLP.

              namespace CLRSnamespace Chapter29open Finsetopen Matrixopen scoped BigOperators

A capacitated network with a real unit cost on each directed pair.

structure CostedFlowNetwork (V : Type*) [Fintype V] extends FlowNetwork V where
  cost : V → V → ℝ

namespace MinimumCostFlowLPvariable {V : Type*} [Fintype V] (N : CostedFlowNetwork V) (demand : ℝ)

The capacity, conservation, and exact-demand constraints.

def IsFeasible (f : V → V → ℝ) : Prop :=
  N.toFlowNetwork.IsFlow f ∧
    FlowNetwork.netOutflow f N.source = demand

Total cost Σ_u Σ_v a_uv f_uv.

def objective (f : V → V → ℝ) : ℝ :=
  ∑ u, ∑ v, N.cost u v * f u v

A demand flow of minimum total cost.

def IsOptimal (f : V → V → ℝ) : Prop :=
  IsFeasible N demand f ∧
    ∀ g, IsFeasible N demand g → objective N f ≤ objective N g

The LP constraints are exactly capacity-feasible flow plus the prescribed source demand.

theorem isFeasible_iff (f : V → V → ℝ) :
    IsFeasible N demand f ↔
      N.toFlowNetwork.IsFlow f ∧
        FlowNetwork.netOutflow f N.source = demand := by
  rfl

The minimization specification agrees exactly with minimum-cost demand flow.

theorem isOptimal_iff (f : V → V → ℝ) :
    IsOptimal N demand f ↔
      (N.toFlowNetwork.IsFlow f ∧
        FlowNetwork.netOutflow f N.source = demand) ∧
      ∀ g, (N.toFlowNetwork.IsFlow g ∧
        FlowNetwork.netOutflow g N.source = demand) →
        (∑ u, ∑ v, N.cost u v * f u v) ≤
          ∑ u, ∑ v, N.cost u v * g u v := by
  rfl

Finite standard-form encoding

section Encodingvariable [DecidableEq V]

The internal vertices at which flow is conserved.

abbrev Internal (N : CostedFlowNetwork V) := {u : V // u ≠ N.source ∧ u ≠ N.sink}

The general-form program for the minimum-cost-flow LP: one nonnegative variable per ordered pair, a capacity inequality per ordered pair, a conservation equality per internal vertex, and a source-demand equality. The objective is minimized (cost is the sum of unit cost times flow).

noncomputable def toGeneralLP (N : CostedFlowNetwork V) (demand : ℝ) : GeneralLP where
  n := Fintype.card (V × V)
  m := Fintype.card (V × V) + Fintype.card (Internal N) + 1
  maximize := false
  c := fun j => N.cost ((Fintype.equivFin (V × V)).symm j).1 ((Fintype.equivFin (V × V)).symm j).2
  rel := fun i =>
    Fin.cases
      (ConstraintRel.eq)
      (fun i => Fin.addCases (motive := fun _ => ConstraintRel)
        (fun _ => ConstraintRel.le)
        (fun _ => ConstraintRel.eq)
        i)
      i
  A := fun i =>
    Fin.cases
      (fun j =>
        (if ((Fintype.equivFin (V × V)).symm j).1 = N.source then 1 else 0) -
          (if ((Fintype.equivFin (V × V)).symm j).2 = N.source then 1 else 0))
      (fun i j => Fin.addCases (motive := fun _ => ℝ)
        (fun e => if j = e then 1 else 0)
        (fun u =>
          let v := ((Fintype.equivFin (Internal N)).symm u).1
          (if ((Fintype.equivFin (V × V)).symm j).1 = v then 1 else 0) -
            (if ((Fintype.equivFin (V × V)).symm j).2 = v then 1 else 0))
        i)
      i
  b := fun i =>
    Fin.cases
      (demand)
      (fun i => Fin.addCases (motive := fun _ => ℝ)
        (fun e => N.capacity ((Fintype.equivFin (V × V)).symm e).1 ((Fintype.equivFin (V × V)).symm e).2)
        (fun _ => 0)
        i)
      i
  free := fun _ => false

The index of the source-demand equality constraint.

abbrev demandIndex (N : CostedFlowNetwork V) :
    Fin (Fintype.card (V × V) + Fintype.card (Internal N) + 1) := 0

The index of the capacity constraint for a directed pair.

abbrev capIndex (N : CostedFlowNetwork V) (e : Fin (Fintype.card (V × V))) :
    Fin (Fintype.card (V × V) + Fintype.card (Internal N) + 1) :=
  (Fin.castAdd (Fintype.card (Internal N)) e).succ

The index of the conservation constraint for an internal vertex.

abbrev consIndex (N : CostedFlowNetwork V) (u : Fin (Fintype.card (Internal N))) :
    Fin (Fintype.card (V × V) + Fintype.card (Internal N) + 1) :=
  (Fin.natAdd (Fintype.card (V × V)) u).succ

              lemma A_demand (N : CostedFlowNetwork V) (demand : ℝ) (j : Fin (Fintype.card (V × V))) :
    (toGeneralLP N demand).A (demandIndex N) j =
      (if ((Fintype.equivFin (V × V)).symm j).1 = N.source then 1 else 0) -
        (if ((Fintype.equivFin (V × V)).symm j).2 = N.source then 1 else 0) := by
  simp [toGeneralLP, demandIndex]lemma b_demand (N : CostedFlowNetwork V) (demand : ℝ) :
    (toGeneralLP N demand).b (demandIndex N) = demand := by
  simp [toGeneralLP, demandIndex]lemma rel_demand (N : CostedFlowNetwork V) (demand : ℝ) :
    (toGeneralLP N demand).rel (demandIndex N) = ConstraintRel.eq := by
  simp [toGeneralLP, demandIndex]

lemma A_cap (N : CostedFlowNetwork V) (demand : ℝ) (e : Fin (Fintype.card (V × V))) (j : Fin (Fintype.card (V × V))) :
    (toGeneralLP N demand).A (capIndex N e) j = if j = e then 1 else 0 := by
  dsimp [toGeneralLP, capIndex]
  rw [FinEncoding.addCases_castAdd]

lemma b_cap (N : CostedFlowNetwork V) (demand : ℝ) (e : Fin (Fintype.card (V × V))) :
    (toGeneralLP N demand).b (capIndex N e) =
      N.capacity ((Fintype.equivFin (V × V)).symm e).1 ((Fintype.equivFin (V × V)).symm e).2 := by
  dsimp [toGeneralLP, capIndex]
  rw [FinEncoding.addCases_castAdd]

lemma rel_cap (N : CostedFlowNetwork V) (demand : ℝ) (e : Fin (Fintype.card (V × V))) :
    (toGeneralLP N demand).rel (capIndex N e) = ConstraintRel.le := by
  dsimp [toGeneralLP, capIndex]
  rw [FinEncoding.addCases_castAdd]

lemma A_cons (N : CostedFlowNetwork V) (demand : ℝ) (u : Fin (Fintype.card (Internal N))) (j : Fin (Fintype.card (V × V))) :
    (toGeneralLP N demand).A (consIndex N u) j =
      (if ((Fintype.equivFin (V × V)).symm j).1 = ((Fintype.equivFin (Internal N)).symm u).1 then 1 else 0) -
        (if ((Fintype.equivFin (V × V)).symm j).2 = ((Fintype.equivFin (Internal N)).symm u).1 then 1 else 0) := by
  dsimp [toGeneralLP, consIndex]
  rw [FinEncoding.addCases_natAdd]

lemma b_cons (N : CostedFlowNetwork V) (demand : ℝ) (u : Fin (Fintype.card (Internal N))) :
    (toGeneralLP N demand).b (consIndex N u) = 0 := by
  dsimp [toGeneralLP, consIndex]
  rw [FinEncoding.addCases_natAdd]

lemma rel_cons (N : CostedFlowNetwork V) (demand : ℝ) (u : Fin (Fintype.card (Internal N))) :
    (toGeneralLP N demand).rel (consIndex N u) = ConstraintRel.eq := by
  dsimp [toGeneralLP, consIndex]
  rw [FinEncoding.addCases_natAdd]

The demand row reads out the source net outflow.


lemma mulVec_demand (N : CostedFlowNetwork V) (demand : ℝ) (f : V → V → ℝ) :
    ((toGeneralLP N demand).A *ᵥ FinEncoding.lift₂ f) (demandIndex N) =
      FlowNetwork.netOutflow f N.source := by
  calc
    ((toGeneralLP N demand).A *ᵥ FinEncoding.lift₂ f) (demandIndex N)
        = ∑ j : Fin (Fintype.card (V × V)), (toGeneralLP N demand).A (demandIndex N) j * FinEncoding.lift₂ f j := by rfl
    _ = ∑ j : Fin (Fintype.card (V × V)),
          ((if ((Fintype.equivFin (V × V)).symm j).1 = N.source then 1 else 0) -
             (if ((Fintype.equivFin (V × V)).symm j).2 = N.source then 1 else 0)) *
            FinEncoding.lift₂ f j := by
          simp only [A_demand]
    _ = (∑ j : Fin (Fintype.card (V × V)),
          (if ((Fintype.equivFin (V × V)).symm j).1 = N.source then 1 else 0) * FinEncoding.lift₂ f j) -
          (∑ j : Fin (Fintype.card (V × V)),
          (if ((Fintype.equivFin (V × V)).symm j).2 = N.source then 1 else 0) * FinEncoding.lift₂ f j) := by
          simp only [sub_mul, Finset.sum_sub_distrib]
    _ = (∑ b : V, f N.source b) - (∑ b : V, f b N.source) := by
          rw [FinEncoding.sum_indicator_fst, FinEncoding.sum_indicator_snd]
    _ = FlowNetwork.netOutflow f N.source := rfl

A capacity row reads out the gross flow on its directed pair.

lemma mulVec_cap (N : CostedFlowNetwork V) (demand : ℝ) (e : Fin (Fintype.card (V × V))) (f : V → V → ℝ) :
    ((toGeneralLP N demand).A *ᵥ FinEncoding.lift₂ f) (capIndex N e) = FinEncoding.lift₂ f e := by
  calc
    ((toGeneralLP N demand).A *ᵥ FinEncoding.lift₂ f) (capIndex N e)
        = ∑ j : Fin (Fintype.card (V × V)), (toGeneralLP N demand).A (capIndex N e) j * FinEncoding.lift₂ f j := by rfl
    _ = ∑ j : Fin (Fintype.card (V × V)), (if j = e then 1 else 0) * FinEncoding.lift₂ f j := by
          simp only [A_cap]
    _ = FinEncoding.lift₂ f e := FinEncoding.sum_indicator_fin e (FinEncoding.lift₂ f)

A conservation row reads out the net outflow at its vertex.


lemma mulVec_cons (N : CostedFlowNetwork V) (demand : ℝ) (u : Fin (Fintype.card (Internal N))) (f : V → V → ℝ) :
    ((toGeneralLP N demand).A *ᵥ FinEncoding.lift₂ f) (consIndex N u) =
      FlowNetwork.outflow f ((Fintype.equivFin (Internal N)).symm u).1 -
        FlowNetwork.inflow f ((Fintype.equivFin (Internal N)).symm u).1 := by
  calc
    ((toGeneralLP N demand).A *ᵥ FinEncoding.lift₂ f) (consIndex N u)
        = ∑ j : Fin (Fintype.card (V × V)), (toGeneralLP N demand).A (consIndex N u) j * FinEncoding.lift₂ f j := by rfl
    _ = ∑ j : Fin (Fintype.card (V × V)),
          ((if ((Fintype.equivFin (V × V)).symm j).1 = ((Fintype.equivFin (Internal N)).symm u).1 then 1 else 0) -
             (if ((Fintype.equivFin (V × V)).symm j).2 = ((Fintype.equivFin (Internal N)).symm u).1 then 1 else 0)) *
            FinEncoding.lift₂ f j := by
          simp only [A_cons]
    _ = (∑ j : Fin (Fintype.card (V × V)),
          (if ((Fintype.equivFin (V × V)).symm j).1 = ((Fintype.equivFin (Internal N)).symm u).1 then 1 else 0) * FinEncoding.lift₂ f j) -
          (∑ j : Fin (Fintype.card (V × V)),
          (if ((Fintype.equivFin (V × V)).symm j).2 = ((Fintype.equivFin (Internal N)).symm u).1 then 1 else 0) * FinEncoding.lift₂ f j) := by
          simp only [sub_mul, Finset.sum_sub_distrib]
    _ = (∑ b : V, f ((Fintype.equivFin (Internal N)).symm u).1 b) -
          (∑ b : V, f b ((Fintype.equivFin (Internal N)).symm u).1) := by
          rw [FinEncoding.sum_indicator_fst, FinEncoding.sum_indicator_snd]
    _ = FlowNetwork.outflow f ((Fintype.equivFin (Internal N)).symm u).1 -
          FlowNetwork.inflow f ((Fintype.equivFin (Internal N)).symm u).1 := rfl

Semantic feasibility is exactly encoded general-form feasibility.


theorem feasible_iff_generalLP (N : CostedFlowNetwork V) (demand : ℝ) (f : V → V → ℝ) :
    IsFeasible N demand f ↔ (toGeneralLP N demand).IsFeasible (FinEncoding.lift₂ f) := by
  unfold IsFeasible GeneralLP.IsFeasible FlowNetwork.IsFlow FlowNetwork.IsCapacityFeasible
  constructor
  · intro hd
    refine ⟨?_, ?_⟩
    · intro j hfree
      simpa [FinEncoding.lift₂, FinEncoding.lift, Function.uncurry] using
        (hd.1.1.1 ((Fintype.equivFin (V × V)).symm j).1 ((Fintype.equivFin (V × V)).symm j).2)
    · intro i
      exact Fin.cases (motive := fun i => match (toGeneralLP N demand).rel i with
          | ConstraintRel.le => ((toGeneralLP N demand).A *ᵥ FinEncoding.lift₂ f) i ≤ (toGeneralLP N demand).b i
          | ConstraintRel.eq => ((toGeneralLP N demand).A *ᵥ FinEncoding.lift₂ f) i = (toGeneralLP N demand).b i
          | ConstraintRel.ge => (toGeneralLP N demand).b i ≤ ((toGeneralLP N demand).A *ᵥ FinEncoding.lift₂ f) i)
        (by
          rw [rel_demand, mulVec_demand, b_demand]
          exact hd.2)
        (fun i => by
          exact Fin.addCases (motive := fun i => match (toGeneralLP N demand).rel i.succ with
              | ConstraintRel.le => ((toGeneralLP N demand).A *ᵥ FinEncoding.lift₂ f) i.succ ≤ (toGeneralLP N demand).b i.succ
              | ConstraintRel.eq => ((toGeneralLP N demand).A *ᵥ FinEncoding.lift₂ f) i.succ = (toGeneralLP N demand).b i.succ
              | ConstraintRel.ge => (toGeneralLP N demand).b i.succ ≤ ((toGeneralLP N demand).A *ᵥ FinEncoding.lift₂ f) i.succ)
            (fun e => by
              rw [rel_cap, mulVec_cap, b_cap]
              simpa [FinEncoding.lift₂, FinEncoding.lift, Function.uncurry] using
                (hd.1.1.2 ((Fintype.equivFin (V × V)).symm e).1 ((Fintype.equivFin (V × V)).symm e).2))
            (fun u => by
              rw [rel_cons, mulVec_cons, b_cons]
              have hc := hd.1.2 ((Fintype.equivFin (Internal N)).symm u).1
                ((Fintype.equivFin (Internal N)).symm u).2.1
                ((Fintype.equivFin (Internal N)).symm u).2.2
              unfold FlowNetwork.ConservesAt at hc
              linarith)
            i)
        i
  · intro h
    refine ⟨⟨⟨?_, ?_⟩, ?_⟩, ?_⟩
    · intro u v
      have hn := h.1 (Fintype.equivFin (V × V) (u, v))
      simpa [FinEncoding.lift₂, FinEncoding.lift, Function.uncurry] using (hn (by simp [toGeneralLP]))
    · intro u v
      have hc := h.2 (capIndex N (Fintype.equivFin (V × V) (u, v)))
      rw [rel_cap] at hc
      rw [mulVec_cap, b_cap] at hc
      simpa [FinEncoding.lift₂, FinEncoding.lift, Function.uncurry] using hc
    · intro u hu hw
      let e : Internal N := ⟨u, ⟨hu, hw⟩⟩
      have hcons := h.2 (consIndex N (Fintype.equivFin (Internal N) e))
      rw [rel_cons] at hcons
      rw [mulVec_cons, b_cons] at hcons
      have hcons' : FlowNetwork.outflow f u - FlowNetwork.inflow f u = 0 := by
        simpa [e] using hcons
      unfold FlowNetwork.ConservesAt
      linarith
    · have hdemand := h.2 (demandIndex N)
      rw [rel_demand] at hdemand
      rw [mulVec_demand, b_demand] at hdemand
      exact hdemand

The general-form objective is the total cost (minimized).


lemma objective_generalLP (N : CostedFlowNetwork V) (demand : ℝ) (f : V → V → ℝ) :
    (toGeneralLP N demand).objective (FinEncoding.lift₂ f) = objective N f := by
  dsimp [GeneralLP.objective, toGeneralLP, objective]
  simp only [dotProduct]
  rw [FinEncoding.sum_reindex]
  simp [FinEncoding.lift₂, FinEncoding.lift, Function.uncurry]
  rw [Fintype.sum_prod_type]

The finite standard-form program obtained by normalizing the minimum-cost general program.

noncomputable def toStandardLP (N : CostedFlowNetwork V) (demand : ℝ) :=
  (toGeneralLP N demand).toStandardLP

The expansion of a gross-flow vector to the normalized standard-form variable vector.

noncomputable def fullLift (N : CostedFlowNetwork V) (demand : ℝ) (f : V → V → ℝ) :
    Fin (Fintype.card (V × V) + Fintype.card (V × V)) → ℝ :=
  (toGeneralLP N demand).lift (FinEncoding.lift₂ f)

The normalized standard-form objective equals the signed total cost.


theorem objective_toStandardLP (N : CostedFlowNetwork V) (demand : ℝ) (f : V → V → ℝ) :
    (toStandardLP N demand).objective (fullLift N demand f) = -objective N f := by
  unfold toStandardLP fullLift
  rw [GeneralLP.objective_lift]
  rw [objective_generalLP]
  simp [toGeneralLP, GeneralLP.objectiveSign]

Feasibility is exactly preserved under the finite standard-form encoding.


theorem feasible_iff_toStandardLP (N : CostedFlowNetwork V) (demand : ℝ) (f : V → V → ℝ) :
    IsFeasible N demand f ↔ (toStandardLP N demand).IsFeasible (fullLift N demand f) := by
  unfold toStandardLP fullLift
  rw [← GeneralLP.feasible_iff_lift]
  exact feasible_iff_generalLP N demand f

              end Encodingend MinimumCostFlowLPend Chapter29end CLRS

CLRSLean.FourthEdition.Chapter_29.Section_29_2_Formulating_Problems_As_Linear_Programs.MulticommodityFlow

29.2: Multicommodity flow as a linear program

Each commodity has its own source, sink, demand, and nonnegative gross flow. All commodities share the edge capacities. The optional cost objective also records Exercise 29.2-7's minimum-cost multicommodity formulation.

This module also records the finite standard-form encoding: one nonnegative variable per commodity/ordered-pair, a shared capacity inequality per ordered pair, a conservation equality per internal commodity/vertex, and a demand equality per commodity, reindexed into a concrete StandardLP.

              namespace CLRSnamespace Chapter29open Finsetopen Matrixopen scoped BigOperators

Source, sink, and nonnegative demand of one commodity.

structure Commodity (V : Type*) where
  source : V
  sink : V
  source_ne_sink : source ≠ sink
  demand : ℝ
  demand_nonnegative : 0 ≤ demand

namespace MulticommodityFlowLPvariable {V K : Type*} [Fintype V] [Fintype K]

Total flow of all commodities on one directed pair.

def aggregate (f : K → V → V → ℝ) (u v : V) : ℝ := ∑ i, f i u v

The CLRS multicommodity feasibility constraints.

def IsFeasible (N : FlowNetwork V) (commodity : K → Commodity V)
    (f : K → V → V → ℝ) : Prop :=
  (∀ i u v, 0 ≤ f i u v) ∧
  (∀ u v, aggregate f u v ≤ N.capacity u v) ∧
  (∀ i u, u ≠ (commodity i).source → u ≠ (commodity i).sink →
    FlowNetwork.ConservesAt (f i) u) ∧
  ∀ i, FlowNetwork.netOutflow (f i) (commodity i).source = (commodity i).demand

An expanded statement of the displayed multicommodity LP.

theorem isFeasible_iff (N : FlowNetwork V) (commodity : K → Commodity V)
    (f : K → V → V → ℝ) :
    IsFeasible N commodity f ↔
      (∀ i u v, 0 ≤ f i u v) ∧
      (∀ u v, (∑ i, f i u v) ≤ N.capacity u v) ∧
      (∀ i u, u ≠ (commodity i).source → u ≠ (commodity i).sink →
        FlowNetwork.inflow (f i) u = FlowNetwork.outflow (f i) u) ∧
      ∀ i, FlowNetwork.outflow (f i) (commodity i).source -
          FlowNetwork.inflow (f i) (commodity i).source = (commodity i).demand := by
  rfl

Total cost of the aggregate multicommodity flow.

def cost (unitCost : V → V → ℝ) (f : K → V → V → ℝ) : ℝ :=
  ∑ u, ∑ v, unitCost u v * aggregate f u v

Minimum-cost feasible multicommodity routing.

def IsMinimumCost (N : FlowNetwork V) (commodity : K → Commodity V)
    (unitCost : V → V → ℝ) (f : K → V → V → ℝ) : Prop :=
  IsFeasible N commodity f ∧
    ∀ g, IsFeasible N commodity g → cost unitCost f ≤ cost unitCost g

Expanded optimality statement for minimum-cost multicommodity flow.

theorem isMinimumCost_iff (N : FlowNetwork V) (commodity : K → Commodity V)
    (unitCost : V → V → ℝ) (f : K → V → V → ℝ) :
    IsMinimumCost N commodity unitCost f ↔
      IsFeasible N commodity f ∧
      ∀ g, IsFeasible N commodity g →
        (∑ u, ∑ v, unitCost u v * (∑ i, f i u v)) ≤
          ∑ u, ∑ v, unitCost u v * (∑ i, g i u v) := by
  rfl

Finite standard-form encoding

section Encodingvariable [DecidableEq V] [DecidableEq K]

Lift a commodity-indexed flow to the Fin-indexed vector over K × (V × V).

noncomputable def lift₃ (f : K → V → V → ℝ) : Fin (Fintype.card (K × (V × V))) → ℝ :=
  FinEncoding.lift (fun p : K × (V × V) => f p.1 p.2.1 p.2.2)

The internal commodity/vertex pairs at which flow is conserved.

abbrev Internal (Variable name `N` is not explicitly referenced.

The binding can be removed (if unused) or named `_` (if used implicitly).

Note: This linter can be disabled with `set_option linter.unusedVariables false`N : FlowNetwork V) (commodity : K → Commodity V) :=
  {p : K × V // p.2 ≠ (commodity p.1).source ∧ p.2 ≠ (commodity p.1).sink}

Sum of an indicator on the commodity and the source coordinate (outflow).


lemma sum_indicator_fst_fst (i : K) (u : V) (f : K → V → V → ℝ) :
    (∑ j : Fin (Fintype.card (K × (V × V))),
      (if ((Fintype.equivFin (K × (V × V))).symm j).1 = i ∧
           ((Fintype.equivFin (K × (V × V))).symm j).2.1 = u then 1 else 0) * lift₃ f j) =
      ∑ b : V, f i u b := by
  rw [FinEncoding.sum_reindex]
  simp [lift₃, FinEncoding.lift]
  rw [Fintype.sum_prod_type]
  rw [Finset.sum_eq_single i]
  · rw [Fintype.sum_prod_type]
    rw [Finset.sum_eq_single u]
    · simp
    · intro a _ ha; simp [ha]
    · intro h; exact False.elim (h (Finset.mem_univ u))
  · intro k _ hk; simp [hk]
  · intro h; exact False.elim (h (Finset.mem_univ i))

Sum of an indicator on the commodity and the target coordinate (inflow).


lemma sum_indicator_fst_snd (i : K) (u : V) (f : K → V → V → ℝ) :
    (∑ j : Fin (Fintype.card (K × (V × V))),
      (if ((Fintype.equivFin (K × (V × V))).symm j).1 = i ∧
           ((Fintype.equivFin (K × (V × V))).symm j).2.2 = u then 1 else 0) * lift₃ f j) =
      ∑ a : V, f i a u := by
  rw [FinEncoding.sum_reindex]
  simp [lift₃, FinEncoding.lift]
  rw [Fintype.sum_prod_type]
  rw [Finset.sum_eq_single i]
  · rw [Fintype.sum_prod_type]
    rw [Finset.sum_comm]
    rw [Finset.sum_eq_single u]
    · simp
    · intro a _ ha; simp [ha]
    · intro h; exact False.elim (h (Finset.mem_univ u))
  · intro k _ hk; simp [hk]
  · intro h; exact False.elim (h (Finset.mem_univ i))

Sum of an indicator on the directed pair (shared capacity).

automatically included section variable(s) unused in theorem `CLRS.Chapter29.MulticommodityFlowLP.sum_indicator_pair`:
  [DecidableEq K]
consider restructuring your `variable` declarations so that the variables are not in scope or explicitly omit them:
  omit [DecidableEq K] in theorem ...

Note: This linter can be disabled with `set_option linter.unusedSectionVars false`automatically included section variable(s) unused in theorem `CLRS.Chapter29.MulticommodityFlowLP.sum_indicator_pair`:
  [DecidableEq K]
consider restructuring your `variable` declarations so that the variables are not in scope or explicitly omit them:
  omit [DecidableEq K] in theorem ...

Note: This linter can be disabled with `set_option linter.unusedSectionVars false`automatically included section variable(s) unused in theorem `CLRS.Chapter29.MulticommodityFlowLP.sum_indicator_pair`:
  [DecidableEq K]
consider restructuring your `variable` declarations so that the variables are not in scope or explicitly omit them:
  omit [DecidableEq K] in theorem ...

Note: This linter can be disabled with `set_option linter.unusedSectionVars false`automatically included section variable(s) unused in theorem `CLRS.Chapter29.MulticommodityFlowLP.sum_indicator_pair`:
  [DecidableEq K]
consider restructuring your `variable` declarations so that the variables are not in scope or explicitly omit them:
  omit [DecidableEq K] in theorem ...

Note: This linter can be disabled with `set_option linter.unusedSectionVars false`automatically included section variable(s) unused in theorem `CLRS.Chapter29.MulticommodityFlowLP.sum_indicator_pair`:
  [DecidableEq K]
consider restructuring your `variable` declarations so that the variables are not in scope or explicitly omit them:
  omit [DecidableEq K] in theorem ...

Note: This linter can be disabled with `set_option linter.unusedSectionVars false`automatically included section variable(s) unused in theorem `CLRS.Chapter29.MulticommodityFlowLP.sum_indicator_pair`:
  [DecidableEq K]
consider restructuring your `variable` declarations so that the variables are not in scope or explicitly omit them:
  omit [DecidableEq K] in theorem ...

Note: This linter can be disabled with `set_option linter.unusedSectionVars false`automatically included section variable(s) unused in theorem `CLRS.Chapter29.MulticommodityFlowLP.sum_indicator_pair`:
  [DecidableEq K]
consider restructuring your `variable` declarations so that the variables are not in scope or explicitly omit them:
  omit [DecidableEq K] in theorem ...

Note: This linter can be disabled with `set_option linter.unusedSectionVars false`
automatically included section variable(s) unused in theorem `CLRS.Chapter29.MulticommodityFlowLP.sum_indicator_pair`:
  [DecidableEq K]
consider restructuring your `variable` declarations so that the variables are not in scope or explicitly omit them:
  omit [DecidableEq K] in theorem ...

Note: This linter can be disabled with `set_option linter.unusedSectionVars false`lemma sum_indicator_pair (u v : V) (f : K → V → V → ℝ) :
    (∑ j : Fin (Fintype.card (K × (V × V))),
      (if ((Fintype.equivFin (K × (V × V))).symm j).2 = (u, v) then 1 else 0) * lift₃ f j) =
      ∑ i : K, f i u v := by
  rw [FinEncoding.sum_reindex]
  simp [lift₃, FinEncoding.lift]
  rw [Fintype.sum_prod_type]
  simp [Finset.sum_ite_eq']

The general-form program for the multicommodity-flow LP.

noncomputable def toGeneralLP (N : FlowNetwork V) (commodity : K → Commodity V)
    (unitCost : V → V → ℝ) : GeneralLP where
  n := Fintype.card (K × (V × V))
  m := Fintype.card (V × V) + Fintype.card (Internal N commodity) + Fintype.card K
  maximize := false
  c := fun j => unitCost ((Fintype.equivFin (K × (V × V))).symm j).2.1 ((Fintype.equivFin (K × (V × V))).symm j).2.2
  rel := fun i =>
    Fin.addCases (motive := fun _ => ConstraintRel)
      (fun ij => Fin.addCases (motive := fun _ => ConstraintRel)
        (fun _ => ConstraintRel.le)
        (fun _ => ConstraintRel.eq)
        ij)
      (fun _ => ConstraintRel.eq)
      i
  A := fun i =>
    Fin.addCases (motive := fun _ => Fin (Fintype.card (K × (V × V))) → ℝ)
      (fun ij => Fin.addCases (motive := fun _ => Fin (Fintype.card (K × (V × V))) → ℝ)
        (fun e j => if ((Fintype.equivFin (K × (V × V))).symm j).2 = (Fintype.equivFin (V × V)).symm e then 1 else 0)
        (fun p j =>
          let iu := (Fintype.equivFin (Internal N commodity)).symm p
          (if ((Fintype.equivFin (K × (V × V))).symm j).1 = iu.1.1 ∧
               ((Fintype.equivFin (K × (V × V))).symm j).2.1 = iu.1.2 then 1 else 0) -
            (if ((Fintype.equivFin (K × (V × V))).symm j).1 = iu.1.1 ∧
                 ((Fintype.equivFin (K × (V × V))).symm j).2.2 = iu.1.2 then 1 else 0))
        ij)
      (fun i j =>
        let ci := (Fintype.equivFin K).symm i
        (if ((Fintype.equivFin (K × (V × V))).symm j).1 = ci ∧
             ((Fintype.equivFin (K × (V × V))).symm j).2.1 = (commodity ci).source then 1 else 0) -
          (if ((Fintype.equivFin (K × (V × V))).symm j).1 = ci ∧
               ((Fintype.equivFin (K × (V × V))).symm j).2.2 = (commodity ci).source then 1 else 0))
      i
  b := fun i =>
    Fin.addCases (motive := fun _ => ℝ)
      (fun ij => Fin.addCases (motive := fun _ => ℝ)
        (fun e => N.capacity ((Fintype.equivFin (V × V)).symm e).1 ((Fintype.equivFin (V × V)).symm e).2)
        (fun _ => 0)
        ij)
      (fun i => (commodity ((Fintype.equivFin K).symm i)).demand)
      i
  free := fun _ => false

The index of the shared-capacity constraint for a directed pair.

abbrev capIndex (N : FlowNetwork V) (commodity : K → Commodity V) (e : Fin (Fintype.card (V × V))) :
    Fin (Fintype.card (V × V) + Fintype.card (Internal N commodity) + Fintype.card K) :=
  Fin.castAdd (Fintype.card K) (Fin.castAdd (Fintype.card (Internal N commodity)) e)

The index of the conservation constraint for an internal commodity/vertex.

abbrev consIndex (N : FlowNetwork V) (commodity : K → Commodity V) (p : Fin (Fintype.card (Internal N commodity))) :
    Fin (Fintype.card (V × V) + Fintype.card (Internal N commodity) + Fintype.card K) :=
  Fin.castAdd (Fintype.card K) (Fin.natAdd (Fintype.card (V × V)) p)

The index of the demand constraint for a commodity.

abbrev demandIndex (N : FlowNetwork V) (commodity : K → Commodity V) (i : Fin (Fintype.card K)) :
    Fin (Fintype.card (V × V) + Fintype.card (Internal N commodity) + Fintype.card K) :=
  Fin.natAdd (Fintype.card (V × V) + Fintype.card (Internal N commodity)) i

              

lemma rel_cap (N : FlowNetwork V) (commodity : K → Commodity V) (e : Fin (Fintype.card (V × V))) :
    (toGeneralLP N commodity (fun _ _ => 0)).rel (capIndex N commodity e) = ConstraintRel.le := by
  dsimp [toGeneralLP, capIndex]
  rw [FinEncoding.addCases_castAdd]
  rw [FinEncoding.addCases_castAdd]

lemma b_cap (N : FlowNetwork V) (commodity : K → Commodity V) (e : Fin (Fintype.card (V × V))) :
    (toGeneralLP N commodity (fun _ _ => 0)).b (capIndex N commodity e) =
      N.capacity ((Fintype.equivFin (V × V)).symm e).1 ((Fintype.equivFin (V × V)).symm e).2 := by
  dsimp [toGeneralLP, capIndex]
  rw [FinEncoding.addCases_castAdd]
  rw [FinEncoding.addCases_castAdd]

lemma A_cap (N : FlowNetwork V) (commodity : K → Commodity V) (e : Fin (Fintype.card (V × V)))
    (j : Fin (Fintype.card (K × (V × V)))) :
    (toGeneralLP N commodity (fun _ _ => 0)).A (capIndex N commodity e) j =
      if ((Fintype.equivFin (K × (V × V))).symm j).2 = (Fintype.equivFin (V × V)).symm e then 1 else 0 := by
  dsimp [toGeneralLP, capIndex]
  rw [FinEncoding.addCases_castAdd]
  rw [FinEncoding.addCases_castAdd]

lemma rel_cons (N : FlowNetwork V) (commodity : K → Commodity V) (p : Fin (Fintype.card (Internal N commodity))) :
    (toGeneralLP N commodity (fun _ _ => 0)).rel (consIndex N commodity p) = ConstraintRel.eq := by
  dsimp [toGeneralLP, consIndex]
  rw [FinEncoding.addCases_castAdd]
  rw [FinEncoding.addCases_natAdd]

lemma b_cons (N : FlowNetwork V) (commodity : K → Commodity V) (p : Fin (Fintype.card (Internal N commodity))) :
    (toGeneralLP N commodity (fun _ _ => 0)).b (consIndex N commodity p) = 0 := by
  dsimp [toGeneralLP, consIndex]
  rw [FinEncoding.addCases_castAdd]
  rw [FinEncoding.addCases_natAdd]

lemma A_cons (N : FlowNetwork V) (commodity : K → Commodity V) (p : Fin (Fintype.card (Internal N commodity)))
    (j : Fin (Fintype.card (K × (V × V)))) :
    (toGeneralLP N commodity (fun _ _ => 0)).A (consIndex N commodity p) j =
      (if ((Fintype.equivFin (K × (V × V))).symm j).1 = ((Fintype.equivFin (Internal N commodity)).symm p).1.1 ∧
           ((Fintype.equivFin (K × (V × V))).symm j).2.1 = ((Fintype.equivFin (Internal N commodity)).symm p).1.2 then 1 else 0) -
        (if ((Fintype.equivFin (K × (V × V))).symm j).1 = ((Fintype.equivFin (Internal N commodity)).symm p).1.1 ∧
             ((Fintype.equivFin (K × (V × V))).symm j).2.2 = ((Fintype.equivFin (Internal N commodity)).symm p).1.2 then 1 else 0) := by
  dsimp [toGeneralLP, consIndex]
  rw [FinEncoding.addCases_castAdd]
  rw [FinEncoding.addCases_natAdd]

lemma rel_demand (N : FlowNetwork V) (commodity : K → Commodity V) (i : Fin (Fintype.card K)) :
    (toGeneralLP N commodity (fun _ _ => 0)).rel (demandIndex N commodity i) = ConstraintRel.eq := by
  dsimp [toGeneralLP, demandIndex]
  rw [FinEncoding.addCases_natAdd]

lemma b_demand (N : FlowNetwork V) (commodity : K → Commodity V) (i : Fin (Fintype.card K)) :
    (toGeneralLP N commodity (fun _ _ => 0)).b (demandIndex N commodity i) =
      (commodity ((Fintype.equivFin K).symm i)).demand := by
  dsimp [toGeneralLP, demandIndex]
  rw [FinEncoding.addCases_natAdd]

lemma A_demand (N : FlowNetwork V) (commodity : K → Commodity V) (i : Fin (Fintype.card K))
    (j : Fin (Fintype.card (K × (V × V)))) :
    (toGeneralLP N commodity (fun _ _ => 0)).A (demandIndex N commodity i) j =
      (if ((Fintype.equivFin (K × (V × V))).symm j).1 = (Fintype.equivFin K).symm i ∧
           ((Fintype.equivFin (K × (V × V))).symm j).2.1 = (commodity ((Fintype.equivFin K).symm i)).source then 1 else 0) -
        (if ((Fintype.equivFin (K × (V × V))).symm j).1 = (Fintype.equivFin K).symm i ∧
             ((Fintype.equivFin (K × (V × V))).symm j).2.2 = (commodity ((Fintype.equivFin K).symm i)).source then 1 else 0) := by
  dsimp [toGeneralLP, demandIndex]
  rw [FinEncoding.addCases_natAdd]

A capacity row reads out the aggregate flow on its directed pair.

lemma mulVec_cap (N : FlowNetwork V) (commodity : K → Commodity V) (e : Fin (Fintype.card (V × V)))
    (f : K → V → V → ℝ) :
    ((toGeneralLP N commodity (fun _ _ => 0)).A *ᵥ lift₃ f) (capIndex N commodity e) =
      aggregate f ((Fintype.equivFin (V × V)).symm e).1 ((Fintype.equivFin (V × V)).symm e).2 := by
  calc
    ((toGeneralLP N commodity (fun _ _ => 0)).A *ᵥ lift₃ f) (capIndex N commodity e)
        = ∑ j : Fin (Fintype.card (K × (V × V))), (toGeneralLP N commodity (fun _ _ => 0)).A (capIndex N commodity e) j * lift₃ f j := by rfl
    _ = ∑ j : Fin (Fintype.card (K × (V × V))),
          (if ((Fintype.equivFin (K × (V × V))).symm j).2 = (Fintype.equivFin (V × V)).symm e then 1 else 0) * lift₃ f j := by
          simp only [A_cap]
    _ = ∑ i : K, f i ((Fintype.equivFin (V × V)).symm e).1 ((Fintype.equivFin (V × V)).symm e).2 := by
          simpa using (sum_indicator_pair ((Fintype.equivFin (V × V)).symm e).1 ((Fintype.equivFin (V × V)).symm e).2 f)
    _ = aggregate f ((Fintype.equivFin (V × V)).symm e).1 ((Fintype.equivFin (V × V)).symm e).2 := rfl

A conservation row reads out the net outflow of one commodity at its internal vertex.


lemma mulVec_cons (N : FlowNetwork V) (commodity : K → Commodity V) (p : Fin (Fintype.card (Internal N commodity)))
    (f : K → V → V → ℝ) :
    ((toGeneralLP N commodity (fun _ _ => 0)).A *ᵥ lift₃ f) (consIndex N commodity p) =
      FlowNetwork.outflow (f ((Fintype.equivFin (Internal N commodity)).symm p).1.1) ((Fintype.equivFin (Internal N commodity)).symm p).1.2 -
        FlowNetwork.inflow (f ((Fintype.equivFin (Internal N commodity)).symm p).1.1) ((Fintype.equivFin (Internal N commodity)).symm p).1.2 := by
  calc
    ((toGeneralLP N commodity (fun _ _ => 0)).A *ᵥ lift₃ f) (consIndex N commodity p)
        = ∑ j : Fin (Fintype.card (K × (V × V))), (toGeneralLP N commodity (fun _ _ => 0)).A (consIndex N commodity p) j * lift₃ f j := by rfl
    _ = ∑ j : Fin (Fintype.card (K × (V × V))),
          ((if ((Fintype.equivFin (K × (V × V))).symm j).1 = ((Fintype.equivFin (Internal N commodity)).symm p).1.1 ∧
                ((Fintype.equivFin (K × (V × V))).symm j).2.1 = ((Fintype.equivFin (Internal N commodity)).symm p).1.2 then 1 else 0) -
             (if ((Fintype.equivFin (K × (V × V))).symm j).1 = ((Fintype.equivFin (Internal N commodity)).symm p).1.1 ∧
                  ((Fintype.equivFin (K × (V × V))).symm j).2.2 = ((Fintype.equivFin (Internal N commodity)).symm p).1.2 then 1 else 0)) *
            lift₃ f j := by
          simp only [A_cons]
    _ = (∑ j : Fin (Fintype.card (K × (V × V))),
          (if ((Fintype.equivFin (K × (V × V))).symm j).1 = ((Fintype.equivFin (Internal N commodity)).symm p).1.1 ∧
               ((Fintype.equivFin (K × (V × V))).symm j).2.1 = ((Fintype.equivFin (Internal N commodity)).symm p).1.2 then 1 else 0) * lift₃ f j) -
          (∑ j : Fin (Fintype.card (K × (V × V))),
          (if ((Fintype.equivFin (K × (V × V))).symm j).1 = ((Fintype.equivFin (Internal N commodity)).symm p).1.1 ∧
               ((Fintype.equivFin (K × (V × V))).symm j).2.2 = ((Fintype.equivFin (Internal N commodity)).symm p).1.2 then 1 else 0) * lift₃ f j) := by
          simp only [sub_mul, Finset.sum_sub_distrib]
    _ = (∑ b : V, f ((Fintype.equivFin (Internal N commodity)).symm p).1.1 ((Fintype.equivFin (Internal N commodity)).symm p).1.2 b) -
          (∑ a : V, f ((Fintype.equivFin (Internal N commodity)).symm p).1.1 a ((Fintype.equivFin (Internal N commodity)).symm p).1.2) := by
          rw [sum_indicator_fst_fst, sum_indicator_fst_snd]
    _ = FlowNetwork.outflow (f ((Fintype.equivFin (Internal N commodity)).symm p).1.1) ((Fintype.equivFin (Internal N commodity)).symm p).1.2 -
          FlowNetwork.inflow (f ((Fintype.equivFin (Internal N commodity)).symm p).1.1) ((Fintype.equivFin (Internal N commodity)).symm p).1.2 := rfl

A demand row reads out the net outflow of one commodity at its source.


lemma mulVec_demand (N : FlowNetwork V) (commodity : K → Commodity V) (i : Fin (Fintype.card K))
    (f : K → V → V → ℝ) :
    ((toGeneralLP N commodity (fun _ _ => 0)).A *ᵥ lift₃ f) (demandIndex N commodity i) =
      FlowNetwork.netOutflow (f ((Fintype.equivFin K).symm i)) (commodity ((Fintype.equivFin K).symm i)).source := by
  calc
    ((toGeneralLP N commodity (fun _ _ => 0)).A *ᵥ lift₃ f) (demandIndex N commodity i)
        = ∑ j : Fin (Fintype.card (K × (V × V))), (toGeneralLP N commodity (fun _ _ => 0)).A (demandIndex N commodity i) j * lift₃ f j := by rfl
    _ = ∑ j : Fin (Fintype.card (K × (V × V))),
          ((if ((Fintype.equivFin (K × (V × V))).symm j).1 = (Fintype.equivFin K).symm i ∧
                ((Fintype.equivFin (K × (V × V))).symm j).2.1 = (commodity ((Fintype.equivFin K).symm i)).source then 1 else 0) -
             (if ((Fintype.equivFin (K × (V × V))).symm j).1 = (Fintype.equivFin K).symm i ∧
                  ((Fintype.equivFin (K × (V × V))).symm j).2.2 = (commodity ((Fintype.equivFin K).symm i)).source then 1 else 0)) *
            lift₃ f j := by
          simp only [A_demand]
    _ = (∑ j : Fin (Fintype.card (K × (V × V))),
          (if ((Fintype.equivFin (K × (V × V))).symm j).1 = (Fintype.equivFin K).symm i ∧
               ((Fintype.equivFin (K × (V × V))).symm j).2.1 = (commodity ((Fintype.equivFin K).symm i)).source then 1 else 0) * lift₃ f j) -
          (∑ j : Fin (Fintype.card (K × (V × V))),
          (if ((Fintype.equivFin (K × (V × V))).symm j).1 = (Fintype.equivFin K).symm i ∧
               ((Fintype.equivFin (K × (V × V))).symm j).2.2 = (commodity ((Fintype.equivFin K).symm i)).source then 1 else 0) * lift₃ f j) := by
          simp only [sub_mul, Finset.sum_sub_distrib]
    _ = (∑ b : V, f ((Fintype.equivFin K).symm i) (commodity ((Fintype.equivFin K).symm i)).source b) -
          (∑ a : V, f ((Fintype.equivFin K).symm i) a (commodity ((Fintype.equivFin K).symm i)).source) := by
          rw [sum_indicator_fst_fst, sum_indicator_fst_snd]
    _ = FlowNetwork.netOutflow (f ((Fintype.equivFin K).symm i)) (commodity ((Fintype.equivFin K).symm i)).source := rfl

Semantic feasibility is exactly encoded general-form feasibility.


theorem feasible_iff_generalLP (N : FlowNetwork V) (commodity : K → Commodity V)
    (f : K → V → V → ℝ) :
    IsFeasible N commodity f ↔ (toGeneralLP N commodity (fun _ _ => 0)).IsFeasible (lift₃ f) := by
  unfold IsFeasible GeneralLP.IsFeasible
  constructor
  · intro hd
    refine ⟨?_, ?_⟩
    · intro j hfree
      simpa [lift₃, FinEncoding.lift] using
        (hd.1 ((Fintype.equivFin (K × (V × V))).symm j).1
          ((Fintype.equivFin (K × (V × V))).symm j).2.1
          ((Fintype.equivFin (K × (V × V))).symm j).2.2)
    · intro i
      exact Fin.addCases (motive := fun i => match (toGeneralLP N commodity (fun _ _ => 0)).rel i with
          | ConstraintRel.le => ((toGeneralLP N commodity (fun _ _ => 0)).A *ᵥ lift₃ f) i ≤ (toGeneralLP N commodity (fun _ _ => 0)).b i
          | ConstraintRel.eq => ((toGeneralLP N commodity (fun _ _ => 0)).A *ᵥ lift₃ f) i = (toGeneralLP N commodity (fun _ _ => 0)).b i
          | ConstraintRel.ge => (toGeneralLP N commodity (fun _ _ => 0)).b i ≤ ((toGeneralLP N commodity (fun _ _ => 0)).A *ᵥ lift₃ f) i)
        (fun ij =>
          Fin.addCases (motive := fun ij => match (toGeneralLP N commodity (fun _ _ => 0)).rel (Fin.castAdd (Fintype.card K) ij) with
              | ConstraintRel.le => ((toGeneralLP N commodity (fun _ _ => 0)).A *ᵥ lift₃ f) (Fin.castAdd (Fintype.card K) ij) ≤ (toGeneralLP N commodity (fun _ _ => 0)).b (Fin.castAdd (Fintype.card K) ij)
              | ConstraintRel.eq => ((toGeneralLP N commodity (fun _ _ => 0)).A *ᵥ lift₃ f) (Fin.castAdd (Fintype.card K) ij) = (toGeneralLP N commodity (fun _ _ => 0)).b (Fin.castAdd (Fintype.card K) ij)
              | ConstraintRel.ge => (toGeneralLP N commodity (fun _ _ => 0)).b (Fin.castAdd (Fintype.card K) ij) ≤ ((toGeneralLP N commodity (fun _ _ => 0)).A *ᵥ lift₃ f) (Fin.castAdd (Fintype.card K) ij))
            (fun e => by
              rw [rel_cap, mulVec_cap, b_cap]
              simpa [aggregate] using
                (hd.2.1 ((Fintype.equivFin (V × V)).symm e).1 ((Fintype.equivFin (V × V)).symm e).2))
            (fun p => by
              rw [rel_cons, mulVec_cons, b_cons]
              have hc := hd.2.2.1 ((Fintype.equivFin (Internal N commodity)).symm p).1.1
                ((Fintype.equivFin (Internal N commodity)).symm p).1.2
                ((Fintype.equivFin (Internal N commodity)).symm p).2.1
                ((Fintype.equivFin (Internal N commodity)).symm p).2.2
              unfold FlowNetwork.ConservesAt at hc
              linarith)
            ij)
        (fun i => by
          rw [rel_demand, mulVec_demand, b_demand]
          exact hd.2.2.2 ((Fintype.equivFin K).symm i))
        i
  · intro h
    refine ⟨?_, ?_, ?_, ?_⟩
    · intro i u v
      have hn := h.1 (Fintype.equivFin (K × (V × V)) (i, (u, v)))
      simpa [lift₃, FinEncoding.lift] using (hn (by simp [toGeneralLP]))
    · intro u v
      have hc := h.2 (capIndex N commodity (Fintype.equivFin (V × V) (u, v)))
      rw [rel_cap] at hc
      rw [mulVec_cap, b_cap] at hc
      simpa [aggregate] using hc
    · intro i u hu hw
      let e : Internal N commodity := ⟨(i, u), ⟨hu, hw⟩⟩
      have hcons := h.2 (consIndex N commodity (Fintype.equivFin (Internal N commodity) e))
      rw [rel_cons] at hcons
      rw [mulVec_cons, b_cons] at hcons
      have hcons' : FlowNetwork.outflow (f i) u - FlowNetwork.inflow (f i) u = 0 := by
        simpa [e] using hcons
      unfold FlowNetwork.ConservesAt
      linarith
    · intro i
      have hdemand := h.2 (demandIndex N commodity (Fintype.equivFin K i))
      rw [rel_demand] at hdemand
      rw [mulVec_demand, b_demand] at hdemand
      simpa using hdemand

The general-form objective is the aggregate minimum cost.


lemma objective_generalLP (N : FlowNetwork V) (commodity : K → Commodity V)
    (unitCost : V → V → ℝ) (f : K → V → V → ℝ) :
    (toGeneralLP N commodity unitCost).objective (lift₃ f) = cost unitCost f := by
  dsimp [GeneralLP.objective, toGeneralLP, cost, aggregate]
  simp only [dotProduct]
  rw [FinEncoding.sum_reindex]
  simp [lift₃, FinEncoding.lift]
  rw [Fintype.sum_prod_type]
  rw [Finset.sum_comm]
  rw [Fintype.sum_prod_type]
  simp [Finset.mul_sum]

The finite standard-form program obtained by normalizing the multicommodity general program.

noncomputable def toStandardLP (N : FlowNetwork V) (commodity : K → Commodity V)
    (unitCost : V → V → ℝ) :=
  (toGeneralLP N commodity unitCost).toStandardLP

The expansion of a commodity flow to the normalized standard-form variable vector.

noncomputable def fullLift (N : FlowNetwork V) (commodity : K → Commodity V)
    (unitCost : V → V → ℝ) (f : K → V → V → ℝ) :
    Fin (Fintype.card (K × (V × V)) + Fintype.card (K × (V × V))) → ℝ :=
  (toGeneralLP N commodity unitCost).lift (lift₃ f)

The normalized standard-form objective equals the signed aggregate cost.


theorem objective_toStandardLP (N : FlowNetwork V) (commodity : K → Commodity V)
    (unitCost : V → V → ℝ) (f : K → V → V → ℝ) :
    (toStandardLP N commodity unitCost).objective (fullLift N commodity unitCost f) = -cost unitCost f := by
  unfold toStandardLP fullLift
  rw [GeneralLP.objective_lift]
  rw [objective_generalLP]
  simp [toGeneralLP, GeneralLP.objectiveSign]

Feasibility is exactly preserved under the finite standard-form encoding.


theorem feasible_iff_toStandardLP (N : FlowNetwork V) (commodity : K → Commodity V)
    (f : K → V → V → ℝ) :
    IsFeasible N commodity f ↔ (toStandardLP N commodity (fun _ _ => 0)).IsFeasible (fullLift N commodity (fun _ _ => 0) f) := by
  unfold toStandardLP fullLift
  rw [← GeneralLP.feasible_iff_lift]
  exact feasible_iff_generalLP N commodity f

              end Encodingend MulticommodityFlowLPend Chapter29end CLRS

Imports

import CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.Definitions
import CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.WeakDuality
import CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.Optimality
import CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.ComplementarySlackness
import CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.TerminalCertificate
import CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.DictionaryBridge
import CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.StrongDuality
import CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.ComplementarySlacknessTheorem

29.3. Duality

The represented layer proves weak duality, extracts a dual certificate from a terminal SIMPLEX dictionary, derives strong duality, and proves both directions of complementary slackness.

Implementation details

The split proof layers remain available outside the main sidebar:

              namespace CLRSnamespace Chapter29end Chapter29end CLRS

Definitions and proofs

CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.Definitions

29.3 Dual linear programs

For a primal maximization program max cᵀx with Ax ≤ b and x ≥ 0, a dual assignment satisfies y ≥ 0 and Aᵀy ≥ c; its objective is bᵀy.

Main declarations:

StandardLP.IsDualFeasible.
StandardLP.dualObjective.

Downstream layers:

Weak duality is proved in the next module.
Strong duality and complementary slackness are proved by the later Section 29.3 modules together with the initialized solver (online material).

              namespace CLRSnamespace Chapter29open Matrixnamespace StandardLP

A nonnegative vector satisfying c ≤ Aᵀy is dual feasible.

def IsDualFeasible {m n : ℕ} (P : StandardLP m n) (y : Fin m → ℝ) : Prop :=
  IsNonnegative y ∧ ∀ j, P.c j ≤ (P.A.transpose *ᵥ y) j

The dual objective value bᵀy.

def dualObjective {m n : ℕ} (P : StandardLP m n) (y : Fin m → ℝ) : ℝ :=
  P.b ⬝ᵥ y

namespace IsDualFeasible

A dual-feasible assignment is coordinatewise nonnegative.

theorem nonnegative {m n : ℕ} {P : StandardLP m n} {y : Fin m → ℝ}
    (hy : P.IsDualFeasible y) : IsNonnegative y :=
  hy.1

A dual-feasible assignment bounds each primal objective coefficient by the corresponding coordinate of Aᵀy.

theorem coefficient_le {m n : ℕ} {P : StandardLP m n} {y : Fin m → ℝ}
    (hy : P.IsDualFeasible y) :
    ∀ j, P.c j ≤ (P.A.transpose *ᵥ y) j :=
  hy.2

              end IsDualFeasibleend StandardLPend Chapter29end CLRS

CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.WeakDuality

29.3 Weak duality

This module proves CLRS Theorem 29.8. For every primal-feasible x and dual-feasible y, the primal objective is bounded by the dual objective: cᵀx ≤ bᵀy.

The proof follows the CLRS calculation cᵀx ≤ (Aᵀy)ᵀx = yᵀAx ≤ yᵀb.

Main result:

StandardLP.weak_duality: CLRS Theorem 29.8.

Downstream layers:

Later Section 29.3 modules prove strong duality (Theorem 29.9) and complementary slackness (Theorem 29.10).

              namespace CLRSnamespace Chapter29open Matrixnamespace StandardLP

Increasing the left factor of a dot product preserves order when the right factor is coordinatewise nonnegative.

theorem dotProduct_mono_right_of_nonnegative {n : ℕ}
    {a b x : Fin n → ℝ} (hab : ∀ i, a i ≤ b i) (hx : IsNonnegative x) :
    a ⬝ᵥ x ≤ b ⬝ᵥ x := by
  simp only [dotProduct]
  exact Finset.sum_le_sum fun i _ =>
    mul_le_mul_of_nonneg_right (hab i) (hx i)

Increasing the right factor of a dot product preserves order when the left factor is coordinatewise nonnegative.

theorem dotProduct_mono_left_of_nonnegative {n : ℕ}
    {a b y : Fin n → ℝ} (hab : ∀ i, a i ≤ b i) (hy : IsNonnegative y) :
    y ⬝ᵥ a ≤ y ⬝ᵥ b := by
  simp only [dotProduct]
  exact Finset.sum_le_sum fun i _ =>
    mul_le_mul_of_nonneg_left (hab i) (hy i)

Moving a matrix transpose across a finite dot product exchanges the order of summation: (Aᵀy)ᵀx = yᵀ(Ax).

theorem transpose_mulVec_dotProduct {m n : ℕ}
    (A : Matrix (Fin m) (Fin n) ℝ) (y : Fin m → ℝ) (x : Fin n → ℝ) :
    (A.transpose *ᵥ y) ⬝ᵥ x = y ⬝ᵥ (A *ᵥ x) := by
  simp only [dotProduct, Matrix.mulVec, Matrix.transpose_apply]
  simp_rw [Finset.sum_mul, Finset.mul_sum]
  conv_lhs => rw [Finset.sum_comm]
  apply Finset.sum_congr rfl
  intro i _
  apply Finset.sum_congr rfl
  intro j _
  ring

CLRS Theorem 29.8 (weak duality). Every primal-feasible objective value is at most every dual-feasible objective value.

theorem weak_duality {m n : ℕ} (P : StandardLP m n)
    {x : Fin n → ℝ} {y : Fin m → ℝ}
    (hx : P.IsFeasible x) (hy : P.IsDualFeasible y) :
    P.objective x ≤ P.dualObjective y := by
  calc
    P.objective x = P.c ⬝ᵥ x := rfl
    _ ≤ (P.A.transpose *ᵥ y) ⬝ᵥ x :=
      dotProduct_mono_right_of_nonnegative hy.coefficient_le hx.1
    _ = y ⬝ᵥ (P.A *ᵥ x) := transpose_mulVec_dotProduct P.A y x
    _ ≤ y ⬝ᵥ P.b :=
      dotProduct_mono_left_of_nonnegative hx.2 hy.nonnegative
    _ = P.dualObjective y := by
      simp only [dualObjective, dotProduct]
      apply Finset.sum_congr rfl
      intro i _
      ring

              end StandardLPend Chapter29end CLRS

CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.Optimality

29.3 Primal and dual optimality specifications

These predicates state the terminal contracts used by strong duality, complementary slackness, and the full initialized SIMPLEX solver.

              namespace CLRSnamespace Chapter29namespace StandardLP

A feasible primal assignment dominating every other feasible assignment.

def IsOptimal (P : StandardLP m n) (x : Fin n → ℝ) : Prop :=
  P.IsFeasible x ∧ ∀ z, P.IsFeasible z → P.objective z ≤ P.objective x

A feasible dual assignment no worse than every other dual assignment.

def IsDualOptimal (P : StandardLP m n) (y : Fin m → ℝ) : Prop :=
  P.IsDualFeasible y ∧
    ∀ z, P.IsDualFeasible z → P.dualObjective y ≤ P.dualObjective z

The primal objective exceeds every real bound on feasible assignments.

def IsUnbounded (P : StandardLP m n) : Prop :=
  ∀ M : ℝ, ∃ x, P.IsFeasible x ∧ M < P.objective x

              end StandardLPend Chapter29end CLRS

CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.ComplementarySlackness

29.3 Complementary slackness

The duality gap splits exactly into primal-slack and dual-slack products. For feasible assignments all products are nonnegative, so equality of the two objectives is equivalent to the textbook complementary-slackness equations.

              namespace CLRSnamespace Chapter29open Matrixopen scoped BigOperatorsnamespace StandardLP

Slack in primal constraint row i.

def primalSlack (P : StandardLP m n) (x : Fin n → ℝ) (i : Fin m) : ℝ :=
  P.b i - (P.A *ᵥ x) i

Slack in dual constraint column j.

def dualSlack (P : StandardLP m n) (y : Fin m → ℝ) (j : Fin n) : ℝ :=
  (P.A.transpose *ᵥ y) j - P.c j

The textbook complementary-slackness equations.

def ComplementarySlackness (P : StandardLP m n)
    (x : Fin n → ℝ) (y : Fin m → ℝ) : Prop :=
  (∀ i, y i * P.primalSlack x i = 0) ∧
    ∀ j, x j * P.dualSlack y j = 0

Exact decomposition of the duality gap into complementary-slackness products.


theorem dualityGap_eq_slackSums (P : StandardLP m n)
    (x : Fin n → ℝ) (y : Fin m → ℝ) :
    P.dualObjective y - P.objective x =
      (∑ i, y i * P.primalSlack x i) +
        ∑ j, x j * P.dualSlack y j := by
  simp only [dualObjective, objective, primalSlack, dualSlack]
  simp_rw [mul_sub]
  rw [Finset.sum_sub_distrib, Finset.sum_sub_distrib]
  change P.b ⬝ᵥ y - P.c ⬝ᵥ x =
    (y ⬝ᵥ P.b - y ⬝ᵥ (P.A *ᵥ x)) +
      (x ⬝ᵥ (P.A.transpose *ᵥ y) - x ⬝ᵥ P.c)
  rw [dotProduct_comm P.b y, dotProduct_comm x (P.A.transpose *ᵥ y),
    dotProduct_comm x P.c, transpose_mulVec_dotProduct]
  ring

For feasible primal and dual assignments, complementary slackness holds exactly when their objective values agree.


theorem complementarySlackness_iff_objective_eq (P : StandardLP m n)
    {x : Fin n → ℝ} {y : Fin m → ℝ}
    (hx : P.IsFeasible x) (hy : P.IsDualFeasible y) :
    P.ComplementarySlackness x y ↔
      P.objective x = P.dualObjective y := by
  have hpnonneg : ∀ i, 0 ≤ y i * P.primalSlack x i := by
    intro i
    exact mul_nonneg (hy.1 i) (sub_nonneg.mpr (hx.2 i))
  have hdnonneg : ∀ j, 0 ≤ x j * P.dualSlack y j := by
    intro j
    exact mul_nonneg (hx.1 j) (sub_nonneg.mpr (hy.2 j))
  constructor
  · rintro ⟨hp, hd⟩
    have hpsum : (∑ i, y i * P.primalSlack x i) = 0 := by
      apply Finset.sum_eq_zero
      intro i _
      exact hp i
    have hdsum : (∑ j, x j * P.dualSlack y j) = 0 := by
      apply Finset.sum_eq_zero
      intro j _
      exact hd j
    have hgap := P.dualityGap_eq_slackSums x y
    rw [hpsum, hdsum, add_zero] at hgap
    linarith
  · intro hobj
    have hgap := P.dualityGap_eq_slackSums x y
    have htotal :
        (∑ i, y i * P.primalSlack x i) +
          ∑ j, x j * P.dualSlack y j = 0 := by
      linarith
    have hpsum_nonneg : 0 ≤ ∑ i, y i * P.primalSlack x i :=
      Finset.sum_nonneg fun i _ => hpnonneg i
    have hdsum_nonneg : 0 ≤ ∑ j, x j * P.dualSlack y j :=
      Finset.sum_nonneg fun j _ => hdnonneg j
    have hpsum : (∑ i, y i * P.primalSlack x i) = 0 := by
      linarith
    have hdsum : (∑ j, x j * P.dualSlack y j) = 0 := by
      linarith
    have hpzero : (fun i => y i * P.primalSlack x i) = 0 :=
      (Fintype.sum_eq_zero_iff_of_nonneg hpnonneg).mp hpsum
    have hdzero : (fun j => x j * P.dualSlack y j) = 0 :=
      (Fintype.sum_eq_zero_iff_of_nonneg hdnonneg).mp hdsum
    exact ⟨fun i => congrFun hpzero i, fun j => congrFun hdzero j⟩

Complementary slackness certifies primal optimality.

theorem optimal_of_complementarySlackness (P : StandardLP m n)
    {x : Fin n → ℝ} {y : Fin m → ℝ}
    (hx : P.IsFeasible x) (hy : P.IsDualFeasible y)
    (hcs : P.ComplementarySlackness x y) : P.IsOptimal x := by
  have heq := (P.complementarySlackness_iff_objective_eq hx hy).1 hcs
  refine ⟨hx, ?_⟩
  intro z hz
  calc
    P.objective z ≤ P.dualObjective y := P.weak_duality hz hy
    _ = P.objective x := heq.symm

Complementary slackness certifies dual optimality.

theorem dualOptimal_of_complementarySlackness (P : StandardLP m n)
    {x : Fin n → ℝ} {y : Fin m → ℝ}
    (hx : P.IsFeasible x) (hy : P.IsDualFeasible y)
    (hcs : P.ComplementarySlackness x y) : P.IsDualOptimal y := by
  have heq := (P.complementarySlackness_iff_objective_eq hx hy).1 hcs
  refine ⟨hy, ?_⟩
  intro z hz
  calc
    P.dualObjective y = P.objective x := heq.symm
    _ ≤ P.dualObjective z := P.weak_duality hx hz

              end StandardLPend Chapter29end CLRS

CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.TerminalCertificate

29.3 Dual certificates from terminal dictionaries

When every reduced cost is nonpositive, the negated objective coefficients of the original slack variables are the textbook dual variables. Dictionary equivalence proves both dual feasibility and equality of objective values.

              namespace CLRSnamespace Chapter29open Matrixopen scoped BigOperatorsnamespace Dictionary

Nonpositive reduced costs make the stable objective coefficient of every variable nonpositive; basic variables have coefficient zero.

theorem objectiveCoeff_nonpos_of_reducedCosts (D : Dictionary m n)
    (hc : ∀ j, D.c j ≤ 0) (q : LPVar m n) :
    D.objectiveCoeff q ≤ 0 := by
  rcases D.exists_basic_or_nonbasic q with ⟨i, rfl⟩ | ⟨j, rfl⟩
  · simp
  · simpa using hc j

The shadow-price vector read from a terminal dictionary.

def dualCertificate (D : Dictionary m n) : Fin m → ℝ :=
  fun i => -D.objectiveCoeff (.inr i)

The shadow prices of a terminal dictionary form a feasible solution of the dual of the original standard-form program.


theorem dualCertificate_isDualFeasible (P : StandardLP m n)
    (D : Dictionary m n) (hEq : P.initialDictionary.Equivalent D)
    (hc : ∀ j, D.c j ≤ 0) :
    P.IsDualFeasible D.dualCertificate := by
  constructor
  · intro i
    exact neg_nonneg.mpr
      (D.objectiveCoeff_nonpos_of_reducedCosts hc (.inr i))
  · intro j
    have hid := hEq.entering_coefficient_identity j
    have horiginal : D.objectiveCoeff (.inl j) ≤ 0 :=
      D.objectiveCoeff_nonpos_of_reducedCosts hc (.inl j)
    have hsum :
        (P.A.transpose *ᵥ D.dualCertificate) j =
          -(∑ i, D.objectiveCoeff (.inr i) * P.A i j) := by
      change (∑ i, P.A i j * (-D.objectiveCoeff (.inr i))) = _
      rw [← Finset.sum_neg_distrib]
      apply Finset.sum_congr rfl
      intro i _
      ring
    change P.c j = D.objectiveCoeff (.inl j) -
      ∑ i, D.objectiveCoeff (.inr i) * P.A i j at hid
    rw [hsum]
    linarith

The objective of the terminal dual certificate equals the terminal dictionary's objective constant.


theorem dualCertificate_objective_eq_v (P : StandardLP m n)
    (D : Dictionary m n) (hEq : P.initialDictionary.Equivalent D) :
    P.dualObjective D.dualCertificate = D.v := by
  let D₀ := P.initialDictionary
  have hobj := hEq.2 D₀.basicAssignment D₀.basicAssignment_satisfies
  have hsum :
      (∑ q, D.objectiveCoeff q * D₀.basicAssignment q) =
        ∑ i, D.objectiveCoeff (.inr i) * P.b i := by
    rw [D₀.sum_eq_sum_labels, Fintype.sum_sum_type]
    change
      (∑ i, D.objectiveCoeff (D₀.basicVar i) *
          D₀.basicAssignment (D₀.basicVar i)) +
        (∑ j, D.objectiveCoeff (D₀.nonbasicVar j) *
          D₀.basicAssignment (D₀.nonbasicVar j)) = _
    simp only [Dictionary.basicAssignment_basicVar,
      Dictionary.basicAssignment_nonbasicVar, mul_zero,
      Finset.sum_const_zero, add_zero]
    change (∑ i, D.objectiveCoeff (.inr i) * P.b i) = _
    rfl
  have hzero :
      0 = D.v + ∑ i, D.objectiveCoeff (.inr i) * P.b i := by
    calc
      0 = D₀.objectiveRhs D₀.basicAssignment := by
        rw [D₀.objectiveRhs_basicAssignment]
        rfl
      _ = D.objectiveRhs D₀.basicAssignment := hobj
      _ = D.v + ∑ q, D.objectiveCoeff q * D₀.basicAssignment q :=
        D.objectiveRhs_eq_fullSum D₀.basicAssignment
      _ = D.v + ∑ i, D.objectiveCoeff (.inr i) * P.b i := by
        rw [hsum]
  change (∑ i, P.b i * (-D.objectiveCoeff (.inr i))) = D.v
  have hneg :
      (∑ i, P.b i * (-D.objectiveCoeff (.inr i))) =
        -(∑ i, D.objectiveCoeff (.inr i) * P.b i) := by
    rw [← Finset.sum_neg_distrib]
    apply Finset.sum_congr rfl
    intro i _
    ring
  rw [hneg]
  linarith

              end Dictionaryend Chapter29end CLRS

CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.DictionaryBridge

29.3 Bridge from dictionary assignments to the primal program

The initial dictionary uses one stable assignment containing both original and slack variables. This module projects that assignment back to the standard-form primal variables and transports optimality and unboundedness.

              namespace CLRSnamespace Chapter29namespace Dictionary

Original-variable coordinates of a complete dictionary assignment.

def assignmentOriginal (z : LPVar m n → ℝ) : Fin n → ℝ :=
  fun j => z (.inl j)

Slack-variable coordinates of a complete dictionary assignment.

def assignmentSlack (z : LPVar m n → ℝ) : Fin m → ℝ :=
  fun i => z (.inr i)

@[simp] theorem assignmentOriginal_combinedAssignment
    (x : Fin n → ℝ) (s : Fin m → ℝ) :
    assignmentOriginal (StandardLP.combinedAssignment x s) = x := by
  funext j
  rfl

@[simp] theorem assignmentSlack_combinedAssignment
    (x : Fin n → ℝ) (s : Fin m → ℝ) :
    assignmentSlack (StandardLP.combinedAssignment x s) = s := by
  funext i
  rfl

Splitting and recombining a stable assignment is the identity.

theorem combinedAssignment_parts (z : LPVar m n → ℝ) :
    StandardLP.combinedAssignment (assignmentOriginal z)
      (assignmentSlack z) = z := by
  funext q
  cases q <;> rfl

A nonnegative assignment satisfying the initial dictionary projects to a feasible assignment of the original standard-form program.


theorem original_feasible_of_initialDictionary (P : StandardLP m n)
    {z : LPVar m n → ℝ} (hznonneg : IsNonnegativeAssignment z)
    (hzsat : P.initialDictionary.Satisfies z) :
    P.IsFeasible (assignmentOriginal z) := by
  apply StandardLP.feasible_of_slackExtension
  refine ⟨fun j => hznonneg (.inl j), fun i => hznonneg (.inr i), ?_⟩
  apply (P.initialDictionary_satisfies_iff
    (assignmentOriginal z) (assignmentSlack z)).1
  rw [combinedAssignment_parts]
  exact hzsat

On every complete assignment, the initial dictionary objective is the standard-form objective of its original-variable projection.


theorem initialDictionary_objectiveRhs_eq_objective_original
    (P : StandardLP m n) (z : LPVar m n → ℝ) :
    P.initialDictionary.objectiveRhs z =
      P.objective (assignmentOriginal z) := by
  rw [← combinedAssignment_parts z]
  exact P.initialDictionary_objectiveRhs _ _

Dictionary optimality for the initial representation yields primal optimality for the original standard-form program.

theorem initialDictionary_optimal_to_standardLP (P : StandardLP m n)
    {z : LPVar m n → ℝ}
    (hz : P.initialDictionary.IsOptimalAssignment z) :
    P.IsOptimal (assignmentOriginal z) := by
  refine ⟨original_feasible_of_initialDictionary P hz.1 hz.2.1, ?_⟩
  intro x hx
  let w := StandardLP.combinedAssignment x (P.slack x)
  have hwext := P.slackExtension_of_feasible hx
  have hwnonneg : IsNonnegativeAssignment w :=
    (StandardLP.combinedAssignment_nonnegative_iff x (P.slack x)).2
      ⟨hwext.1, hwext.2.1⟩
  have hwsat : P.initialDictionary.Satisfies w :=
    P.initialDictionary_satisfies_of_slackExtension hwext
  have hle := hz.2.2 w hwnonneg hwsat
  calc
    P.objective x = P.initialDictionary.objectiveRhs w := by
      exact (P.initialDictionary_objectiveRhs x (P.slack x)).symm
    _ ≤ P.initialDictionary.objectiveRhs z := hle
    _ = P.objective (assignmentOriginal z) :=
      initialDictionary_objectiveRhs_eq_objective_original P z

Dictionary unboundedness for the initial representation yields unboundedness of the original standard-form program.


theorem initialDictionary_unbounded_to_standardLP (P : StandardLP m n)
    (h : P.initialDictionary.IsUnbounded) : P.IsUnbounded := by
  intro M
  obtain ⟨z, hznonneg, hzsat, hzobj⟩ := h M
  refine ⟨assignmentOriginal z,
    original_feasible_of_initialDictionary P hznonneg hzsat, ?_⟩
  rw [← initialDictionary_objectiveRhs_eq_objective_original P z]
  exact hzobj

              end Dictionaryend Chapter29end CLRS

CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.StrongDuality

29.3 Strong duality from finite SIMPLEX

For an initially basic-feasible standard-form program, finite Bland-SIMPLEX either produces an unbounded primal ray or a terminal dictionary. In the terminal case its basic assignment and shadow prices are primal/dual optimal and have the same objective value.

              namespace CLRSnamespace Chapter29namespace StandardLP

Strong duality from any basic-feasible dictionary equivalent to the program's initial dictionary.


theorem strongDuality_or_unbounded_of_equivalent_isBasicFeasible
    (P : StandardLP m n) (D : Dictionary m n)
    (hEq₀ : P.initialDictionary.Equivalent D)
    (hD : D.IsBasicFeasible) :
    P.IsUnbounded ∨
      ∃ x y, P.IsOptimal x ∧ P.IsDualOptimal y ∧
        P.objective x = P.dualObjective y := by
  let fuel := Dictionary.basisCount m n
  cases hrun : D.simplexRun fuel with
  | optimal terminal hc =>
      have hEq : P.initialDictionary.Equivalent terminal :=
        hEq₀.trans (by
          simpa [hrun, Dictionary.SimplexRunResult.terminalDictionary] using
            D.simplexRun_equivalent fuel)
      have hterminal : terminal.IsBasicFeasible := by
        simpa [hrun, Dictionary.SimplexRunResult.terminalDictionary] using
          D.simplexRun_isBasicFeasible fuel hD
      have hdictOptimal :
          P.initialDictionary.IsOptimalAssignment terminal.basicAssignment :=
        hEq.isOptimalAssignment
          (terminal.basicAssignment_optimal_of_reducedCosts_nonpos
            hterminal hc)
      let x := Dictionary.assignmentOriginal terminal.basicAssignment
      let y := terminal.dualCertificate
      have hx : P.IsOptimal x := by
        exact Dictionary.initialDictionary_optimal_to_standardLP P
          hdictOptimal
      have hyfeasible : P.IsDualFeasible y := by
        exact terminal.dualCertificate_isDualFeasible P hEq hc
      have hterminalSatD₀ :
          P.initialDictionary.Satisfies terminal.basicAssignment :=
        (hEq.1 terminal.basicAssignment).2
          terminal.basicAssignment_satisfies
      have hprimalValue : P.objective x = terminal.v := by
        calc
          P.objective x =
              P.initialDictionary.objectiveRhs terminal.basicAssignment :=
            (Dictionary.initialDictionary_objectiveRhs_eq_objective_original
              P terminal.basicAssignment).symm
          _ = terminal.objectiveRhs terminal.basicAssignment :=
            hEq.2 terminal.basicAssignment hterminalSatD₀
          _ = terminal.v := terminal.objectiveRhs_basicAssignment
      have hdualValue : P.dualObjective y = terminal.v :=
        terminal.dualCertificate_objective_eq_v P hEq
      have hvalue : P.objective x = P.dualObjective y :=
        hprimalValue.trans hdualValue.symm
      have hcs : P.ComplementarySlackness x y :=
        (P.complementarySlackness_iff_objective_eq hx.1 hyfeasible).2 hvalue
      have hy : P.IsDualOptimal y :=
        P.dualOptimal_of_complementarySlackness hx.1 hyfeasible hcs
      exact Or.inr ⟨x, y, hx, hy, hvalue⟩
  | unbounded terminal entering he ha =>
      have hEq : P.initialDictionary.Equivalent terminal :=
        hEq₀.trans (by
          simpa [hrun, Dictionary.SimplexRunResult.terminalDictionary] using
            D.simplexRun_equivalent fuel)
      have hterminal : terminal.IsBasicFeasible := by
        simpa [hrun, Dictionary.SimplexRunResult.terminalDictionary] using
          D.simplexRun_isBasicFeasible fuel hD
      have hunbounded : P.initialDictionary.IsUnbounded :=
        hEq.isUnbounded
          (terminal.unbounded_of_entering_column hterminal entering he.1 ha)
      exact Or.inl
        (Dictionary.initialDictionary_unbounded_to_standardLP P hunbounded)
  | exhausted terminal =>
      exact False.elim (D.simplexRun_basisCount_not_exhausted hD (by
        simp [fuel, hrun, Dictionary.SimplexRunResult.IsExhausted]))

Strong duality, with the alternative unbounded outcome made explicit, for programs whose initial slack dictionary is basic feasible.

theorem strongDuality_or_unbounded_of_initialDictionary_isBasicFeasible
    (P : StandardLP m n) (hP : P.initialDictionary.IsBasicFeasible) :
    P.IsUnbounded ∨
      ∃ x y, P.IsOptimal x ∧ P.IsDualOptimal y ∧
        P.objective x = P.dualObjective y :=
  P.strongDuality_or_unbounded_of_equivalent_isBasicFeasible
    P.initialDictionary (Dictionary.Equivalent.refl _) hP

If the primal is not unbounded, finite SIMPLEX yields primal and dual optima with equal objective values.

theorem strongDuality_of_initialDictionary_isBasicFeasible
    (P : StandardLP m n) (hP : P.initialDictionary.IsBasicFeasible)
    (hbounded : ¬P.IsUnbounded) :
    ∃ x y, P.IsOptimal x ∧ P.IsDualOptimal y ∧
      P.objective x = P.dualObjective y := by
  rcases P.strongDuality_or_unbounded_of_initialDictionary_isBasicFeasible hP
      with hunbounded | hopt
  · exact False.elim (hbounded hunbounded)
  · exact hopt

              end StandardLPend Chapter29end CLRS

CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.ComplementarySlacknessTheorem

29.3 The complementary-slackness theorem

This closes the converse direction of the textbook theorem: for feasible primal and dual assignments, simultaneous optimality is equivalent to the complementary-slackness equations.

              namespace CLRSnamespace Chapter29namespace StandardLP

Any dual-feasible assignment gives a finite upper bound on all primal feasible objective values.

theorem not_isUnbounded_of_isDualFeasible (P : StandardLP m n)
    {y : Fin m → ℝ} (hy : P.IsDualFeasible y) : ¬P.IsUnbounded := by
  intro hunbounded
  obtain ⟨x, hx, hlarge⟩ := hunbounded (P.dualObjective y)
  have hweak := P.weak_duality hx hy
  linarith

CLRS complementary-slackness theorem for an initially basic-feasible program: feasible primal and dual assignments are both optimal exactly when their complementary products vanish.


theorem complementarySlackness_iff_optimal_of_initialDictionary_isBasicFeasible
    (P : StandardLP m n) (hP : P.initialDictionary.IsBasicFeasible)
    {x : Fin n → ℝ} {y : Fin m → ℝ}
    (hx : P.IsFeasible x) (hy : P.IsDualFeasible y) :
    P.ComplementarySlackness x y ↔
      P.IsOptimal x ∧ P.IsDualOptimal y := by
  constructor
  · intro hcs
    exact ⟨P.optimal_of_complementarySlackness hx hy hcs,
      P.dualOptimal_of_complementarySlackness hx hy hcs⟩
  · rintro ⟨hxoptimal, hyoptimal⟩
    obtain ⟨x₀, y₀, hx₀, hy₀, hvalue₀⟩ :=
      P.strongDuality_of_initialDictionary_isBasicFeasible hP
        (P.not_isUnbounded_of_isDualFeasible hy)
    have hprimal : P.objective x = P.objective x₀ :=
      le_antisymm (hx₀.2 x hx) (hxoptimal.2 x₀ hx₀.1)
    have hdual : P.dualObjective y = P.dualObjective y₀ :=
      le_antisymm (hyoptimal.2 y₀ hy₀.1) (hy₀.2 y hy)
    have hvalue : P.objective x = P.dualObjective y := by
      calc
        P.objective x = P.objective x₀ := hprimal
        _ = P.dualObjective y₀ := hvalue₀
        _ = P.dualObjective y := hdual.symm
    exact (P.complementarySlackness_iff_objective_eq hx hy).2 hvalue

              end StandardLPend Chapter29end CLRS

CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.GeneralStrongDuality

29.3 General strong duality

Phase II supplies the equivalent basic-feasible dictionary required by the finite-SIMPLEX strong-duality proof. Projecting both primal and dual locked coordinates yields the unrestricted textbook theorems for the original LP.

This module is the canonical fourth-edition owner of the general (no initial-basis hypothesis) strong-duality and complementary-slackness theorems. The phase-I/II initialization machinery it builds on remains online material (reachable through CLRSLean.OnlineMaterial); only the resulting duality declarations live here.

              namespace CLRSnamespace Chapter29namespace StandardLP

Every feasible standard-form program is either unbounded or has primal and dual optima with equal values.


theorem strongDuality_or_unbounded_of_feasible (P : StandardLP m n)
    (hfeasible : ∃ x, P.IsFeasible x) :
    P.IsUnbounded ∨
      ∃ x y, P.IsOptimal x ∧ P.IsDualOptimal y ∧
        P.objective x = P.dualObjective y := by
  have hlocked :=
    P.lockedAuxiliary.strongDuality_or_unbounded_of_equivalent_isBasicFeasible
      P.phaseTwoStart P.phaseTwoStart_equivalent_lockedAuxiliary
      (P.phaseTwoStart_isBasicFeasible hfeasible)
  rcases hlocked with hunbounded | ⟨z, y, hz, hy, hvalue⟩
  · exact Or.inl (P.lockedAuxiliary_unbounded_to_original hunbounded)
  · let x₀ := auxiliaryTail z
    let y₀ := lockedDualTail y
    have hx₀ : P.IsOptimal x₀ :=
      P.lockedAuxiliary_optimal_to_original hz
    have hy₀feasible : P.IsDualFeasible y₀ :=
      P.lockedAuxiliary_dualFeasible_to_original hy.1
    have hvalue₀ : P.objective x₀ = P.dualObjective y₀ := by
      calc
        P.objective x₀ = P.lockedAuxiliary.objective z :=
          (P.lockedAuxiliary_objective z).symm
        _ = P.lockedAuxiliary.dualObjective y := hvalue
        _ = P.dualObjective y₀ := P.lockedAuxiliary_dualObjective y
    have hy₀ : P.IsDualOptimal y₀ := by
      refine ⟨hy₀feasible, ?_⟩
      intro y' hy'
      calc
        P.dualObjective y₀ = P.objective x₀ := hvalue₀.symm
        _ ≤ P.dualObjective y' := P.weak_duality hx₀.1 hy'
    exact Or.inr ⟨x₀, y₀, hx₀, hy₀, hvalue₀⟩

General strong duality: a feasible, bounded primal program has primal and dual optima with equal objective values.

theorem strongDuality (P : StandardLP m n)
    (hfeasible : ∃ x, P.IsFeasible x) (hbounded : ¬P.IsUnbounded) :
    ∃ x y, P.IsOptimal x ∧ P.IsDualOptimal y ∧
      P.objective x = P.dualObjective y := by
  rcases P.strongDuality_or_unbounded_of_feasible hfeasible with
    hunbounded | hopt
  · exact False.elim (hbounded hunbounded)
  · exact hopt

An attained primal optimum rules out primal unboundedness.

theorem not_isUnbounded_of_isOptimal (P : StandardLP m n)
    {x : Fin n → ℝ} (hx : P.IsOptimal x) : ¬P.IsUnbounded := by
  intro hunbounded
  obtain ⟨z, hz, hlarge⟩ := hunbounded (P.objective x)
  exact (not_lt_of_ge (hx.2 z hz)) hlarge

Textbook strong-duality form for a given attained primal optimum.

theorem strongDuality_of_isOptimal (P : StandardLP m n)
    {x : Fin n → ℝ} (hx : P.IsOptimal x) :
    ∃ y, P.IsDualOptimal y ∧ P.objective x = P.dualObjective y := by
  obtain ⟨x₀, y, hx₀, hy, hvalue₀⟩ :=
    P.strongDuality ⟨x, hx.1⟩ (P.not_isUnbounded_of_isOptimal hx)
  have hprimal : P.objective x = P.objective x₀ :=
    le_antisymm (hx₀.2 x hx.1) (hx.2 x₀ hx₀.1)
  exact ⟨y, hy, hprimal.trans hvalue₀⟩

CLRS complementary-slackness theorem without an initial-basis hypothesis: feasible primal and dual assignments are simultaneously optimal exactly when all complementary products vanish.

theorem complementarySlackness_iff_optimal (P : StandardLP m n)
    {x : Fin n → ℝ} {y : Fin m → ℝ}
    (hx : P.IsFeasible x) (hy : P.IsDualFeasible y) :
    P.ComplementarySlackness x y ↔
      P.IsOptimal x ∧ P.IsDualOptimal y := by
  constructor
  · intro hcs
    exact ⟨P.optimal_of_complementarySlackness hx hy hcs,
      P.dualOptimal_of_complementarySlackness hx hy hcs⟩
  · rintro ⟨hxoptimal, hyoptimal⟩
    obtain ⟨y₀, hy₀, hvalue₀⟩ := P.strongDuality_of_isOptimal hxoptimal
    have hdual : P.dualObjective y = P.dualObjective y₀ :=
      le_antisymm (hyoptimal.2 y₀ hy₀.1) (hy₀.2 y hy)
    have hvalue : P.objective x = P.dualObjective y :=
      hvalue₀.trans hdual.symm
    exact (P.complementarySlackness_iff_objective_eq hx hy).2 hvalue

              end StandardLPend Chapter29end CLRS

Scope and implementation notes

Imports

import CLRSLean.Chapter_29
import CLRSLean.FourthEdition.Chapter_29.Section_29_1_Standard_And_Slack_Forms
import CLRSLean.FourthEdition.Chapter_29.Section_29_2_Formulating_Problems_As_Linear_Programs
import CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality
import CLRSLean.FourthEdition.Chapter_29.Section_29_3_Duality.GeneralStrongDuality
import CLRSLean.FourthEdition.Chapter_29.Section_29_1_Standard_And_Slack_Forms.SolverWrapper

Current source

Sections 29.1--29.3 are native fourth-edition sections (standard and slack forms, formulating problems as linear programs, and duality), imported directly from Section 29.1, Section 29.2, and Section 29.3. The duality development builds on the online simplex machinery (legacy Section 29.3), and the difference-constraints bridge imports the fourth-edition single-source-shortest-paths sources (Chapter 22). Declarations keep their current namespaces; the third-edition-numbered imports CLRSLean.Chapter_29 and CLRSLean.Chapter_29.Section_29_* forward to these sources.

Implementation details

The supporting implementation pages remain available outside the main sidebar:

Coverage boundary

The native sections supply the represented fourth-edition linear-programming sections (§29.1 formulations and algorithms, §29.2 formulating problems, §29.3 duality with Theorems 29.8--29.10). The detailed simplex algorithm (legacy Section 29.3) and the initial basic feasible solution (legacy Section 29.5) are retained as supplementary online material (reachable through CLRSLean.OnlineMaterial). The canonical guide already imports SolverWrapper, whose normalization/initialized-solver composition certifies all three outcomes for general-form programs. The construction uses classical real-valued choices; no polynomial SIMPLEX or machine-runtime bound is claimed.

See docs/clrs-fourth-edition-map.csv for the section-level mapping and docs/migrations/clrs4.md for compatibility and deprecation policy.

CLRS, fourth edition · Chapter 29 of 35