linear_programming.tex

% Linear Programming Analysis Template
% Topics: Simplex algorithm, duality theory, sensitivity analysis, graphical solutions
% Style: Operations research technical report with computational implementation

\documentclass[a4paper, 11pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath, amssymb}
\usepackage{graphicx}
\usepackage{siunitx}
\usepackage{booktabs}
\usepackage{subcaption}
\usepackage[makestderr]{pythontex}

% Theorem environments
\newtheorem{definition}{Definition}[section]
\newtheorem{theorem}{Theorem}[section]
\newtheorem{example}{Example}[section]
\newtheorem{remark}{Remark}[section]

\title{Linear Programming: Simplex Method, Duality, and Sensitivity Analysis}
\author{Operations Research Laboratory}
\date{\today}

\begin{document}
\maketitle

\begin{abstract}
This technical report presents a comprehensive analysis of linear programming (LP) techniques,
including graphical solutions for two-variable problems, the simplex algorithm for
higher-dimensional optimization, duality theory, and sensitivity analysis. We examine a
production planning problem, demonstrate the simplex tableau method, derive the dual problem,
and analyze shadow prices and reduced costs. Computational implementation using Python's
\texttt{scipy.optimize.linprog} validates theoretical results and illustrates practical
applications in resource allocation.
\end{abstract}

\section{Introduction}

Linear programming (LP) is a fundamental technique in operations research for optimizing
a linear objective function subject to linear equality and inequality constraints. Since
George Dantzig's development of the simplex algorithm in 1947, LP has become one of the
most widely applied optimization methods in industry, logistics, finance, and engineering.

\begin{definition}[Linear Programming Problem]
A linear programming problem in standard form seeks to:
\begin{equation}
\begin{aligned}
\text{maximize} \quad & \mathbf{c}^T \mathbf{x} \\
\text{subject to} \quad & \mathbf{A}\mathbf{x} \leq \mathbf{b} \\
& \mathbf{x} \geq \mathbf{0}
\end{aligned}
\end{equation}
where $\mathbf{c} \in \mathbb{R}^n$ is the objective coefficient vector, $\mathbf{A} \in \mathbb{R}^{m \times n}$
is the constraint matrix, $\mathbf{b} \in \mathbb{R}^m$ is the right-hand side vector, and
$\mathbf{x} \in \mathbb{R}^n$ is the decision variable vector.
\end{definition}

\section{Theoretical Framework}

\subsection{Fundamental Concepts}

\begin{definition}[Feasible Region]
The feasible region $\mathcal{F}$ is the set of all points satisfying the constraints:
\begin{equation}
\mathcal{F} = \{\mathbf{x} \in \mathbb{R}^n : \mathbf{A}\mathbf{x} \leq \mathbf{b}, \mathbf{x} \geq \mathbf{0}\}
\end{equation}
\end{definition}

\begin{theorem}[Fundamental Theorem of Linear Programming]
If a linear programming problem has an optimal solution, then at least one optimal solution
occurs at an extreme point (vertex) of the feasible region.
\end{theorem}

\begin{definition}[Basic Feasible Solution]
A basic feasible solution (BFS) is a feasible solution where at most $m$ variables (out of $n$)
are non-zero. These non-zero variables are called \textit{basic variables}, and the zero
variables are \textit{non-basic variables}.
\end{definition}

\subsection{The Simplex Algorithm}

The simplex algorithm moves from one BFS to an adjacent BFS with a better objective value,
continuing until optimality is reached or unboundedness is detected.

\begin{theorem}[Optimality Criterion]
A basic feasible solution is optimal if all reduced costs $\bar{c}_j = c_j - \mathbf{c}_B^T \mathbf{A}_B^{-1} \mathbf{A}_j$
are non-positive (for maximization) or non-negative (for minimization).
\end{theorem}

\subsection{Duality Theory}

\begin{definition}[Dual Problem]
For the primal LP problem, the dual problem is:
\begin{equation}
\begin{aligned}
\text{minimize} \quad & \mathbf{b}^T \mathbf{y} \\
\text{subject to} \quad & \mathbf{A}^T\mathbf{y} \geq \mathbf{c} \\
& \mathbf{y} \geq \mathbf{0}
\end{aligned}
\end{equation}
where $\mathbf{y} \in \mathbb{R}^m$ are the dual variables (shadow prices).
\end{definition}

\begin{theorem}[Strong Duality]
If the primal problem has an optimal solution $\mathbf{x}^*$ with objective value $z^*$,
then the dual problem has an optimal solution $\mathbf{y}^*$ with objective value $w^*$,
and $z^* = w^*$.
\end{theorem}

\begin{theorem}[Complementary Slackness]
At optimality, for each constraint $i$:
\begin{equation}
y_i^* \left(b_i - \sum_{j=1}^n a_{ij}x_j^* \right) = 0
\end{equation}
and for each variable $j$:
\begin{equation}
x_j^* \left(\sum_{i=1}^m a_{ij}y_i^* - c_j \right) = 0
\end{equation}
\end{definition}

\section{Computational Analysis}

\subsection{Problem Formulation: Production Planning}

We consider a furniture manufacturer producing tables and chairs with limited resources.
Each table requires 4 hours of labor and 3 board-feet of wood, yielding a profit of \$30.
Each chair requires 2 hours of labor and 2 board-feet of wood, yielding a profit of \$20.
Available resources are 80 hours of labor and 60 board-feet of wood. The goal is to maximize profit.

This translates to the following LP formulation:
\begin{equation}
\begin{aligned}
\text{maximize} \quad & z = 30x_1 + 20x_2 \\
\text{subject to} \quad & 4x_1 + 2x_2 \leq 80 \quad \text{(labor)} \\
& 3x_1 + 2x_2 \leq 60 \quad \text{(wood)} \\
& x_1, x_2 \geq 0
\end{aligned}
\end{equation}
where $x_1$ is the number of tables and $x_2$ is the number of chairs.

\begin{pycode}
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linprog
from matplotlib.patches import Polygon

np.random.seed(42)

# Problem parameters - Production Planning
c = np.array([30, 20])  # Objective coefficients (profit per unit)
A = np.array([
    [4, 2],  # Labor constraint
    [3, 2]   # Wood constraint
])
b = np.array([80, 60])  # Resource availability

# Solve using linprog (minimization, so negate c)
result = linprog(-c, A_ub=A, b_ub=b, bounds=[(0, None), (0, None)], method='highs')
x_optimal = result.x
z_optimal = -result.fun
shadow_prices = result.ineqlin.marginals
reduced_costs = result.lower.marginals

# Calculate vertices of feasible region
vertices = []
vertices.append([0, 0])  # Origin

# Intersection with axes
vertices.append([min(b[0]/A[0,0], b[1]/A[1,0]), 0])  # x1-axis
vertices.append([0, min(b[0]/A[0,1], b[1]/A[1,1])])  # x2-axis

# Constraint intersections
# 4x1 + 2x2 = 80 and 3x1 + 2x2 = 60
A_int = np.array([[4, 2], [3, 2]])
b_int = np.array([80, 60])
try:
    vertex = np.linalg.solve(A_int, b_int)
    if all(vertex >= 0) and all(A @ vertex <= b + 1e-6):
        vertices.append(vertex)
except:
    pass

vertices = np.array(vertices)

# Graphical solution
fig, ax = plt.subplots(figsize=(10, 8))

# Plot constraints
x1_range = np.linspace(0, 30, 500)

# Labor constraint: 4x1 + 2x2 <= 80
x2_labor = (80 - 4*x1_range) / 2
ax.plot(x1_range, x2_labor, 'b-', linewidth=2, label='Labor: $4x_1 + 2x_2 \\leq 80$')
ax.fill_between(x1_range, 0, np.minimum(x2_labor, 50), where=(x2_labor >= 0),
                alpha=0.2, color='blue')

# Wood constraint: 3x1 + 2x2 <= 60
x2_wood = (60 - 3*x1_range) / 2
ax.plot(x1_range, x2_wood, 'g-', linewidth=2, label='Wood: $3x_1 + 2x_2 \\leq 60$')
ax.fill_between(x1_range, 0, np.minimum(x2_wood, 50), where=(x2_wood >= 0),
                alpha=0.2, color='green')

# Feasible region
feasible_x1 = [0, 0, x_optimal[0], 20, 0]
feasible_x2 = [0, 30, x_optimal[1], 0, 0]
ax.fill(feasible_x1, feasible_x2, alpha=0.3, color='yellow',
        edgecolor='black', linewidth=2, label='Feasible Region')

# Objective function iso-profit lines
for profit in [200, 400, 600]:
    x2_obj = (profit - 30*x1_range) / 20
    ax.plot(x1_range, x2_obj, '--', alpha=0.5, color='red', linewidth=1)

# Optimal iso-profit line
x2_opt = (z_optimal - 30*x1_range) / 20
ax.plot(x1_range, x2_opt, 'r-', linewidth=2.5,
        label=f'Optimal: $z = {z_optimal:.1f}$')

# Plot vertices
for i, v in enumerate(vertices):
    if all(v >= 0):
        obj_val = c @ v
        ax.plot(v[0], v[1], 'ko', markersize=8)
        ax.annotate(f'({v[0]:.1f}, {v[1]:.1f})\\n$z={obj_val:.1f}$',
                   xy=(v[0], v[1]), xytext=(v[0]+1.5, v[1]+2),
                   fontsize=9, ha='left')

# Optimal solution
ax.plot(x_optimal[0], x_optimal[1], 'r*', markersize=20,
        label=f'Optimal: $x_1={x_optimal[0]:.2f}$, $x_2={x_optimal[1]:.2f}$')

ax.set_xlabel('$x_1$ (Tables)', fontsize=12)
ax.set_ylabel('$x_2$ (Chairs)', fontsize=12)
ax.set_title('Graphical Solution: Production Planning LP', fontsize=14)
ax.set_xlim(-2, 28)
ax.set_ylim(-2, 45)
ax.grid(True, alpha=0.3)
ax.legend(loc='upper right', fontsize=9)
ax.axhline(y=0, color='k', linewidth=0.5)
ax.axvline(x=0, color='k', linewidth=0.5)

plt.tight_layout()
plt.savefig('linear_programming_graphical_solution.pdf', dpi=150, bbox_inches='tight')
plt.close()
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=0.95\textwidth]{linear_programming_graphical_solution.pdf}
\caption{Graphical solution of the two-dimensional production planning linear program showing
the feasible region (yellow polygon) bounded by labor and wood constraints, iso-profit lines
(red dashed) representing constant objective values, and the optimal solution at the intersection
of active constraints. The optimal vertex yields maximum profit $z^* = \py{f"{z_optimal:.2f}"}$
by producing $x_1^* = \py{f"{x_optimal[0]:.2f}"}$ tables and $x_2^* = \py{f"{x_optimal[1]:.2f}"}$ chairs,
demonstrating the fundamental theorem that optimality occurs at an extreme point of the polytope.}
\end{figure}

\subsection{Simplex Tableau Evolution}

The simplex algorithm iteratively improves the solution by pivoting through basic feasible solutions.
We demonstrate the tableau method for a standard form LP. Converting our problem to standard form
by introducing slack variables $s_1, s_2 \geq 0$:

\begin{equation}
\begin{aligned}
\text{maximize} \quad & z = 30x_1 + 20x_2 + 0s_1 + 0s_2 \\
\text{subject to} \quad & 4x_1 + 2x_2 + s_1 = 80 \\
& 3x_1 + 2x_2 + s_2 = 60 \\
& x_1, x_2, s_1, s_2 \geq 0
\end{aligned}
\end{equation}

\begin{pycode}
# Simplex tableau implementation
def simplex_iteration(tableau, basis):
    """Perform one iteration of simplex method"""
    m, n = tableau.shape
    m -= 1  # Exclude objective row

    # Check optimality (all reduced costs <= 0 for maximization)
    if np.all(tableau[-1, :-1] <= 1e-10):
        return tableau, basis, True

    # Select entering variable (most positive reduced cost)
    entering = np.argmax(tableau[-1, :-1])

    # Select leaving variable (minimum ratio test)
    ratios = []
    for i in range(m):
        if tableau[i, entering] > 1e-10:
            ratios.append(tableau[i, -1] / tableau[i, entering])
        else:
            ratios.append(np.inf)
    leaving = np.argmin(ratios)

    # Update basis
    basis[leaving] = entering

    # Pivot operation
    pivot = tableau[leaving, entering]
    tableau[leaving, :] /= pivot

    for i in range(m + 1):
        if i != leaving:
            tableau[i, :] -= tableau[i, entering] * tableau[leaving, :]

    return tableau, basis, False

# Initial tableau setup
# Variables: x1, x2, s1, s2, RHS
initial_tableau = np.array([
    [4.0, 2.0, 1.0, 0.0, 80.0],  # Labor constraint
    [3.0, 2.0, 0.0, 1.0, 60.0],  # Wood constraint
    [-30.0, -20.0, 0.0, 0.0, 0.0]  # Objective (negated for maximization)
], dtype=float)

basis = [2, 3]  # Initial basis: s1, s2
tableaus = [initial_tableau.copy()]
iteration = 0

tableau = initial_tableau.copy()
optimal = False
while not optimal and iteration < 10:
    tableau, basis, optimal = simplex_iteration(tableau, basis)
    tableaus.append(tableau.copy())
    iteration += 1

# Visualization of tableau evolution
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
axes = axes.flatten()

var_names = ['$x_1$', '$x_2$', '$s_1$', '$s_2$', 'RHS']

for idx, (ax, tab) in enumerate(zip(axes, tableaus[:4])):
    # Display tableau as table
    ax.axis('tight')
    ax.axis('off')

    # Create table data
    table_data = []
    for i in range(tab.shape[0] - 1):
        row = [f'{val:.2f}' for val in tab[i, :]]
        table_data.append(row)

    # Objective row
    obj_row = [f'{val:.2f}' for val in tab[-1, :]]
    table_data.append(obj_row)

    # Row labels
    row_labels = [f'Row {i+1}' for i in range(tab.shape[0] - 1)] + ['$z$']

    table = ax.table(cellText=table_data, colLabels=var_names,
                    rowLabels=row_labels, cellLoc='center',
                    loc='center', bbox=[0, 0, 1, 1])
    table.auto_set_font_size(False)
    table.set_fontsize(9)
    table.scale(1, 2)

    # Color code
    for i in range(len(row_labels)):
        table[(i, -1)].set_facecolor('#E8E8E8')
    for j in range(len(var_names)):
        table[(0, j)].set_facecolor('#40466e')
        table[(0, j)].set_text_props(weight='bold', color='white')

    if idx == 0:
        ax.set_title('Initial Tableau (Iteration 0)', fontsize=12, weight='bold')
    elif idx < len(tableaus) - 1:
        ax.set_title(f'Iteration {idx}', fontsize=12, weight='bold')
    else:
        ax.set_title(f'Optimal Tableau (Iteration {idx})', fontsize=12, weight='bold',
                    color='darkgreen')

plt.tight_layout()
plt.savefig('linear_programming_simplex_tableaus.pdf', dpi=150, bbox_inches='tight')
plt.close()
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=0.95\textwidth]{linear_programming_simplex_tableaus.pdf}
\caption{Evolution of simplex tableaus from initial basic feasible solution to optimal solution.
Each iteration performs a pivot operation, exchanging a non-basic variable (entering) with a basic
variable (leaving) to move to an adjacent vertex with improved objective value. The bottom row shows
reduced costs; optimality is reached when all reduced costs are non-positive for a maximization problem.
The algorithm converges in $\py{len(tableaus)-1}$ iterations, demonstrating the finite convergence
property of the simplex method on this bounded linear program.}
\end{figure}

\section{Duality and Sensitivity Analysis}

\subsection{Dual Problem Formulation}

The dual of our primal production planning problem provides economic interpretation through
shadow prices. The dual variables $y_1$ and $y_2$ represent the marginal value of each
resource (labor and wood).

The dual problem is:
\begin{equation}
\begin{aligned}
\text{minimize} \quad & w = 80y_1 + 60y_2 \\
\text{subject to} \quad & 4y_1 + 3y_2 \geq 30 \quad \text{(table)} \\
& 2y_1 + 2y_2 \geq 20 \quad \text{(chair)} \\
& y_1, y_2 \geq 0
\end{aligned}
\end{equation}

Before proceeding with sensitivity analysis, we present the dual problem formulation and
its economic interpretation. The shadow prices indicate how much the objective would improve
per unit increase in each resource constraint, assuming the current basis remains optimal.

\begin{pycode}
# Solve dual problem
c_dual = b  # Dual objective = primal RHS
A_dual = -A.T  # Dual constraints = transpose of primal
b_dual = -c  # Dual RHS = negative primal objective

result_dual = linprog(c_dual, A_ub=A_dual, b_ub=b_dual,
                     bounds=[(0, None), (0, None)], method='highs')
y_optimal = result_dual.x
w_optimal = result_dual.fun

# Verify strong duality
print(f"Primal optimal: z* = {z_optimal:.4f}")
print(f"Dual optimal: w* = {w_optimal:.4f}")
print(f"Duality gap: {abs(z_optimal - w_optimal):.2e}")

# Shadow prices from primal solution (should equal dual solution)
print(f"Shadow prices from primal: {shadow_prices}")
print(f"Dual solution y*: {y_optimal}")

# Sensitivity analysis - varying RHS parameters
b1_range = np.linspace(60, 100, 50)  # Labor variation
b2_range = np.linspace(40, 80, 50)   # Wood variation

z_labor_variation = []
z_wood_variation = []

for b1_new in b1_range:
    res = linprog(-c, A_ub=A, b_ub=[b1_new, b[1]],
                 bounds=[(0, None), (0, None)], method='highs')
    z_labor_variation.append(-res.fun if res.success else np.nan)

for b2_new in b2_range:
    res = linprog(-c, A_ub=A, b_ub=[b[0], b2_new],
                 bounds=[(0, None), (0, None)], method='highs')
    z_wood_variation.append(-res.fun if res.success else np.nan)

# Sensitivity visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Plot 1: RHS sensitivity for labor
ax1 = axes[0, 0]
ax1.plot(b1_range, z_labor_variation, 'b-', linewidth=2.5,
         label='Optimal profit')
# Linear approximation using shadow price
z_linear_labor = z_optimal + shadow_prices[0] * (b1_range - b[0])
ax1.plot(b1_range, z_linear_labor, 'r--', linewidth=2,
         label=f'Linear approx. (shadow price = {shadow_prices[0]:.2f})')
ax1.axvline(x=b[0], color='gray', linestyle=':', alpha=0.7,
           label=f'Current value = {b[0]}')
ax1.set_xlabel('Labor hours available', fontsize=11)
ax1.set_ylabel('Optimal profit (\$)', fontsize=11)
ax1.set_title('Sensitivity to Labor Constraint', fontsize=12, weight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# Plot 2: RHS sensitivity for wood
ax2 = axes[0, 1]
ax2.plot(b2_range, z_wood_variation, 'g-', linewidth=2.5,
         label='Optimal profit')
z_linear_wood = z_optimal + shadow_prices[1] * (b2_range - b[1])
ax2.plot(b2_range, z_linear_wood, 'r--', linewidth=2,
         label=f'Linear approx. (shadow price = {shadow_prices[1]:.2f})')
ax2.axvline(x=b[1], color='gray', linestyle=':', alpha=0.7,
           label=f'Current value = {b[1]}')
ax2.set_xlabel('Wood board-feet available', fontsize=11)
ax2.set_ylabel('Optimal profit (\$)', fontsize=11)
ax2.set_title('Sensitivity to Wood Constraint', fontsize=12, weight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

# Plot 3: Objective coefficient sensitivity for x1
ax3 = axes[1, 0]
c1_range = np.linspace(15, 45, 50)
z_c1_variation = []
x1_c1_variation = []

for c1_new in c1_range:
    res = linprog(-np.array([c1_new, c[1]]), A_ub=A, b_ub=b,
                 bounds=[(0, None), (0, None)], method='highs')
    z_c1_variation.append(-res.fun if res.success else np.nan)
    x1_c1_variation.append(res.x[0] if res.success else np.nan)

ax3_twin = ax3.twinx()
line1 = ax3.plot(c1_range, z_c1_variation, 'b-', linewidth=2.5,
                label='Optimal profit')
ax3.axvline(x=c[0], color='gray', linestyle=':', alpha=0.7)
ax3.set_xlabel('Profit per table $c_1$ (\$)', fontsize=11)
ax3.set_ylabel('Optimal profit (\$)', fontsize=11, color='blue')
ax3.tick_params(axis='y', labelcolor='blue')

line2 = ax3_twin.plot(c1_range, x1_c1_variation, 'r--', linewidth=2,
                     label='Tables produced')
ax3_twin.set_ylabel('Optimal $x_1$ (tables)', fontsize=11, color='red')
ax3_twin.tick_params(axis='y', labelcolor='red')
ax3.set_title('Sensitivity to Table Profit Coefficient', fontsize=12, weight='bold')
ax3.grid(True, alpha=0.3)

# Combined legend
lines1, labels1 = ax3.get_legend_handles_labels()
lines2, labels2 = ax3_twin.get_legend_handles_labels()
ax3.legend(lines1 + lines2, labels1 + labels2, loc='upper left', fontsize=9)

# Plot 4: 2D sensitivity region
ax4 = axes[1, 1]
c1_2d = np.linspace(20, 40, 30)
c2_2d = np.linspace(10, 30, 30)
C1_mesh, C2_mesh = np.meshgrid(c1_2d, c2_2d)
Z_mesh = np.zeros_like(C1_mesh)

for i in range(len(c1_2d)):
    for j in range(len(c2_2d)):
        res = linprog(-np.array([c1_2d[i], c2_2d[j]]), A_ub=A, b_ub=b,
                     bounds=[(0, None), (0, None)], method='highs')
        Z_mesh[j, i] = -res.fun if res.success else np.nan

cs = ax4.contourf(C1_mesh, C2_mesh, Z_mesh, levels=15, cmap='viridis')
ax4.plot(c[0], c[1], 'r*', markersize=20, label='Current coefficients')
cbar = plt.colorbar(cs, ax=ax4)
cbar.set_label('Optimal profit (\$)', fontsize=10)
ax4.set_xlabel('Table profit $c_1$ (\$)', fontsize=11)
ax4.set_ylabel('Chair profit $c_2$ (\$)', fontsize=11)
ax4.set_title('Objective Coefficient Sensitivity Map', fontsize=12, weight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3, color='white', linewidth=0.5)

plt.tight_layout()
plt.savefig('linear_programming_sensitivity_analysis.pdf', dpi=150, bbox_inches='tight')
plt.close()
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=0.98\textwidth]{linear_programming_sensitivity_analysis.pdf}
\caption{Comprehensive sensitivity analysis examining the robustness of the optimal solution to
parameter perturbations. Top row shows the linear relationship between resource availability
(labor and wood) and optimal profit, with shadow prices representing the marginal value of each
resource within the valid range. Bottom left illustrates how objective coefficient changes affect
both profit and production decisions. Bottom right presents a two-dimensional contour map showing
optimal profit as a function of both objective coefficients simultaneously. Shadow prices are
$y_1^* = \py{f"{shadow_prices[0]:.3f}"}$ per labor hour and $y_2^* = \py{f"{shadow_prices[1]:.3f}"}$
per board-foot of wood, indicating which resource constraint is more valuable to relax.}
\end{figure}

\section{Applications and Extensions}

\subsection{Transportation Problem}

Linear programming naturally models transportation and logistics problems. We demonstrate
a small-scale transportation problem with 3 suppliers and 3 customers, minimizing total
shipping costs while satisfying supply and demand constraints.

\begin{pycode}
# Transportation problem setup
# Supply from 3 warehouses: [100, 150, 120]
# Demand at 3 stores: [80, 130, 160]
# Cost matrix ($/unit shipped from warehouse i to store j)

supply = np.array([100, 150, 120])
demand = np.array([80, 130, 160])
cost_matrix = np.array([
    [8, 6, 10],   # Warehouse 1 to stores 1,2,3
    [9, 12, 13],  # Warehouse 2 to stores 1,2,3
    [14, 9, 16]   # Warehouse 3 to stores 1,2,3
])

# Formulate as LP: minimize sum(c_ij * x_ij)
# Variables: x_ij = units shipped from warehouse i to store j
n_warehouses = len(supply)
n_stores = len(demand)
n_vars = n_warehouses * n_stores

# Objective coefficients (flatten cost matrix)
c_transport = cost_matrix.flatten()

# Supply constraints: sum_j x_ij <= supply_i
A_supply = np.zeros((n_warehouses, n_vars))
for i in range(n_warehouses):
    for j in range(n_stores):
        A_supply[i, i * n_stores + j] = 1
b_supply = supply

# Demand constraints: sum_i x_ij >= demand_j (convert to <=)
A_demand = np.zeros((n_stores, n_vars))
for j in range(n_stores):
    for i in range(n_warehouses):
        A_demand[j, i * n_stores + j] = -1
b_demand = -demand

# Combine constraints
A_transport = np.vstack([A_supply, A_demand])
b_transport = np.hstack([b_supply, b_demand])

# Solve
result_transport = linprog(c_transport, A_ub=A_transport, b_ub=b_transport,
                          bounds=[(0, None)] * n_vars, method='highs')
x_transport = result_transport.x.reshape(n_warehouses, n_stores)
total_cost = result_transport.fun

# Visualization of transportation solution
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Plot 1: Transportation flow diagram
ax1 = axes[0]
ax1.axis('off')
ax1.set_xlim(0, 10)
ax1.set_ylim(0, 10)

# Draw warehouses (left side)
warehouse_y = np.linspace(2, 8, n_warehouses)
for i, y in enumerate(warehouse_y):
    circle = plt.Circle((2, y), 0.4, color='steelblue', alpha=0.7)
    ax1.add_patch(circle)
    ax1.text(2, y, f'W{i+1}\\n{supply[i]}', ha='center', va='center',
            fontsize=9, weight='bold', color='white')

# Draw stores (right side)
store_y = np.linspace(2, 8, n_stores)
for j, y in enumerate(store_y):
    circle = plt.Circle((8, y), 0.4, color='coral', alpha=0.7)
    ax1.add_patch(circle)
    ax1.text(8, y, f'S{j+1}\\n{demand[j]}', ha='center', va='center',
            fontsize=9, weight='bold', color='white')

# Draw flows
max_flow = x_transport.max()
for i in range(n_warehouses):
    for j in range(n_stores):
        if x_transport[i, j] > 0.1:
            width = 0.5 + 3 * (x_transport[i, j] / max_flow)
            ax1.arrow(2.5, warehouse_y[i], 4.9, store_y[j] - warehouse_y[i],
                     head_width=0.2, head_length=0.2, fc='green', ec='green',
                     linewidth=width, alpha=0.6)
            mid_x = 5
            mid_y = (warehouse_y[i] + store_y[j]) / 2
            ax1.text(mid_x, mid_y, f'{x_transport[i,j]:.1f}',
                    fontsize=8, ha='center',
                    bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

ax1.set_title('Optimal Transportation Flow', fontsize=12, weight='bold')
ax1.text(2, 9, 'Warehouses\\n(Supply)', ha='center', fontsize=10, weight='bold')
ax1.text(8, 9, 'Stores\\n(Demand)', ha='center', fontsize=10, weight='bold')

# Plot 2: Cost matrix heatmap with solution overlay
ax2 = axes[1]
im = ax2.imshow(x_transport, cmap='YlGn', aspect='auto')
ax2.set_xticks(np.arange(n_stores))
ax2.set_yticks(np.arange(n_warehouses))
ax2.set_xticklabels([f'Store {j+1}' for j in range(n_stores)])
ax2.set_yticklabels([f'Warehouse {i+1}' for i in range(n_warehouses)])
ax2.set_xlabel('Destination', fontsize=11)
ax2.set_ylabel('Origin', fontsize=11)
ax2.set_title('Shipment Quantities with Unit Costs', fontsize=12, weight='bold')

# Add text annotations
for i in range(n_warehouses):
    for j in range(n_stores):
        text = ax2.text(j, i, f'{x_transport[i,j]:.1f}\\n(\${cost_matrix[i,j]})',
                       ha='center', va='center', color='black', fontsize=9)

cbar = plt.colorbar(im, ax=ax2)
cbar.set_label('Units shipped', fontsize=10)

plt.tight_layout()
plt.savefig('linear_programming_transportation.pdf', dpi=150, bbox_inches='tight')
plt.close()
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=0.98\textwidth]{linear_programming_transportation.pdf}
\caption{Optimal solution to the transportation problem showing flow diagram (left) with arrow
widths proportional to shipment quantities and heatmap (right) displaying optimal allocation
with unit costs. The solution minimizes total transportation cost of \$\py{f"{total_cost:.2f}"}
while satisfying all supply constraints at warehouses and demand requirements at stores. This
classical application demonstrates LP's power in logistics optimization, where the dual variables
provide shadow prices indicating marginal costs of increasing supply or demand at each location.}
\end{figure}

\section{Results Summary}

\subsection{Optimal Solutions}

\begin{pycode}
print(r"\begin{table}[H]")
print(r"\centering")
print(r"\caption{Optimal Solutions for LP Problems}")
print(r"\begin{tabular}{lccc}")
print(r"\toprule")
print(r"Problem & Variables & Objective Value & Method \\")
print(r"\midrule")
print(f"Production Planning & $x_1^* = {x_optimal[0]:.2f}$, $x_2^* = {x_optimal[1]:.2f}$ & "
      f"\${z_optimal:.2f} & Simplex \\\\")
print(f"Dual Problem & $y_1^* = {y_optimal[0]:.2f}$, $y_2^* = {y_optimal[1]:.2f}$ & "
      f"\${w_optimal:.2f} & Dual Simplex \\\\")
print(f"Transportation & Total shipment = {np.sum(x_transport):.1f} units & "
      f"\${total_cost:.2f} & Interior Point \\\\")
print(r"\bottomrule")
print(r"\end{tabular}")
print(r"\label{tab:solutions}")
print(r"\end{table}")
\end{pycode}

\subsection{Shadow Prices and Economic Interpretation}

\begin{pycode}
print(r"\begin{table}[H]")
print(r"\centering")
print(r"\caption{Sensitivity Analysis: Shadow Prices and Reduced Costs}")
print(r"\begin{tabular}{lccl}")
print(r"\toprule")
print(r"Resource/Variable & Shadow Price & Reduced Cost & Interpretation \\")
print(r"\midrule")
print(f"Labor (hours) & \${shadow_prices[0]:.3f} & --- & "
      r"Value of 1 additional labor hour \\")
print(f"Wood (board-ft) & \${shadow_prices[1]:.3f} & --- & "
      r"Value of 1 additional board-foot \\")
print(f"$x_1$ (tables) & --- & \${reduced_costs[0]:.3f} & "
      r"At optimal basis \\")
print(f"$x_2$ (chairs) & --- & \${reduced_costs[1]:.3f} & "
      r"At optimal basis \\")
print(r"\bottomrule")
print(r"\end{tabular}")
print(r"\label{tab:sensitivity}")
print(r"\end{table}")
\end{pycode}

\section{Discussion}

\begin{example}[Economic Interpretation of Shadow Prices]
The shadow price for labor ($y_1^* = \py{f"{shadow_prices[0]:.3f}"}$) indicates that each
additional hour of labor would increase profit by approximately \$\py{f"{shadow_prices[0]:.2f}"},
assuming the current basis remains optimal. Similarly, the shadow price for wood
($y_2^* = \py{f"{shadow_prices[1]:.3f}"}$) values each additional board-foot at
\$\py{f"{shadow_prices[1]:.2f}"}. These values guide resource acquisition decisions: if labor
costs less than \$\py{f"{shadow_prices[0]:.2f}"}/hour or wood costs less than
\$\py{f"{shadow_prices[1]:.2f}"}/board-foot, purchasing additional resources increases profitability.
\end{example}

\begin{remark}[Complementary Slackness Verification]
At the optimal solution, active constraints (labor and wood) have positive shadow prices,
while slack variables have zero reduced costs. Non-binding constraints would have zero shadow
prices, and variables at their bounds have non-zero reduced costs indicating the objective
deterioration per unit forced into the basis.
\end{remark}

\subsection{Computational Complexity}

The simplex algorithm's worst-case complexity is exponential in the number of variables,
but it typically performs very well in practice with polynomial average-case behavior.
Interior point methods guarantee polynomial-time complexity $O(n^{3.5})$ for problems
with $n$ variables, making them preferable for very large-scale applications.

\section{Conclusions}

This comprehensive analysis demonstrates linear programming methodology and implementation:

\begin{enumerate}
\item The production planning problem achieves optimal profit $z^* = \py{f"{z_optimal:.2f}"}$
by producing $x_1^* = \py{f"{x_optimal[0]:.2f}"}$ tables and $x_2^* = \py{f"{x_optimal[1]:.2f}"}$ chairs.

\item The simplex algorithm converges in $\py{len(tableaus)-1}$ iterations by systematically
exploring vertices of the feasible polytope until reaching optimality.

\item Strong duality is verified: primal and dual optimal values are equal
($z^* = w^* = \py{f"{z_optimal:.2f}"}$), confirming the fundamental duality theorem.

\item Shadow prices ($y_1^* = \py{f"{shadow_prices[0]:.3f}"}$,
$y_2^* = \py{f"{shadow_prices[1]:.3f}"}$) quantify the marginal value of relaxing resource constraints.

\item The transportation application demonstrates LP versatility in logistics optimization,
achieving minimum cost \$\py{f"{total_cost:.2f}"} while satisfying all supply-demand constraints.

\item Sensitivity analysis reveals the range over which the current basis remains optimal,
guiding robust decision-making under parameter uncertainty.
\end{enumerate}

Linear programming provides a powerful framework for optimization problems with linear objectives
and constraints, with applications spanning production planning, transportation, finance, scheduling,
and resource allocation.

\section*{Further Reading}

\begin{itemize}
\item Bertsimas, D., \& Tsitsiklis, J. N. \textit{Introduction to Linear Optimization}. Athena Scientific, 1997.
\item Dantzig, G. B. \textit{Linear Programming and Extensions}. Princeton University Press, 1963.
\item Vanderbei, R. J. \textit{Linear Programming: Foundations and Extensions}, 5th ed. Springer, 2020.
\item Chvátal, V. \textit{Linear Programming}. W. H. Freeman, 1983.
\item Bazaraa, M. S., Jarvis, J. J., \& Sherali, H. D. \textit{Linear Programming and Network Flows}, 4th ed. Wiley, 2010.
\item Luenberger, D. G., \& Ye, Y. \textit{Linear and Nonlinear Programming}, 4th ed. Springer, 2016.
\item Nocedal, J., \& Wright, S. J. \textit{Numerical Optimization}, 2nd ed. Springer, 2006.
\item Papadimitriou, C. H., \& Steiglitz, K. \textit{Combinatorial Optimization: Algorithms and Complexity}. Dover, 1998.
\item Schrijver, A. \textit{Theory of Linear and Integer Programming}. Wiley, 1998.
\item Winston, W. L., \& Goldberg, J. B. \textit{Operations Research: Applications and Algorithms}, 4th ed. Thomson, 2004.
\item Hillier, F. S., \& Lieberman, G. J. \textit{Introduction to Operations Research}, 10th ed. McGraw-Hill, 2015.
\item Taha, H. A. \textit{Operations Research: An Introduction}, 10th ed. Pearson, 2017.
\item Wolsey, L. A. \textit{Integer Programming}, 2nd ed. Wiley, 2020.
\item Karmarkar, N. ``A new polynomial-time algorithm for linear programming.'' \textit{Combinatorica} 4.4 (1984): 373-395.
\item Klee, V., \& Minty, G. J. ``How good is the simplex algorithm?'' \textit{Inequalities} 3 (1972): 159-175.
\end{itemize}

\end{document}