⁹⁶ views
ubuntu2404

Kernel: SageMath 10.7

Linear Programming with SageMath - Complete Guide

Learning Objectives

By the end of this comprehensive tutorial, you will:

Master linear programming fundamentals and mathematical formulation
Understand the geometric interpretation of LP problems and feasible regions
Apply duality theory and perform sensitivity analysis
Solve real-world optimization problems in production, transportation, and finance
Use integer programming for discrete decision problems
Implement advanced LP techniques using SageMath's optimization toolkit
Analyze optimality conditions and interpret shadow prices

What You'll Learn

Linear Programming Fundamentals - Standard forms, feasible regions, and optimality
Geometric Interpretation - Visualizing constraints and objective functions
Duality Theory - Primal-dual relationships and economic interpretation
Transportation Problems - Supply chain and logistics optimization
Integer Programming - Discrete optimization and knapsack problems
Portfolio Optimization - Financial applications and risk management
Sensitivity Analysis - Understanding solution robustness
Advanced Applications - Real-world case studies and industry examples

🏛 Historical Context: The Birth of Linear Programming

Linear Programming was born during World War II from the urgent need to optimize military logistics and resource allocation. The story begins with George Dantzig (1914-2005), often called the "father of linear programming."

Key Historical Milestones:

1939-1941: Dantzig worked on planning methods for the US Air Force
1947: Dantzig developed the Simplex Algorithm - one of the most important algorithms of the 20th century
1951: The first computer implementation solved a 42-variable problem
1975: Dantzig shared the National Medal of Science
1979: Khachian's ellipsoid method proved LP is polynomial-time solvable
1984: Karmarkar's interior-point method revolutionized large-scale optimization

Revolutionary Impact:

"The real impact of linear programming has been in the thousands of applications and the way it has changed our approach to problem solving." - George Dantzig

Today, LP powers everything from airline scheduling and financial portfolio management to supply chain optimization and resource allocation across countless industries.

Introduction to Linear Programming

Linear Programming (LP) is a mathematical optimization technique for maximizing or minimizing a linear objective function subject to linear equality and inequality constraints.

Standard Form:

\begin{align} \text{minimize} \quad & c^T x \\ \text{subject to} \quad & Ax = b \\ & x \geq 0 \end{align}

Where:

$x \in \mathbb{R}^n$ are the decision variables
$c \in \mathbb{R}^n$ is the objective coefficient vector
$A \in \mathbb{R}^{m \times n}$ is the constraint matrix
$b \in \mathbb{R}^m$ is the right-hand side vector

Key Properties:

Linearity: Both objective function and constraints are linear
Feasible Region: Set of all points satisfying all constraints
Optimal Solution: Occurs at vertices (extreme points) of feasible region
Fundamental Theorem: If an optimal solution exists, at least one vertex is optimal

Real-World Applications:

Manufacturing: Production planning and resource allocation
Transportation: Route optimization and logistics
Finance: Portfolio optimization and risk management
Energy: Power generation scheduling
Agriculture: Crop planning and livestock feed optimization

In [3]:

# Setup: Import required libraries
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import linprog
import sympy as sp

# Create simple LP solver using scipy
def solve_linear_program(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, bounds=None, method='highs'):
    """Solve linear program using scipy.optimize.linprog"""
    result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method=method)
    return result

print(" Linear Programming - Complete Tutorial")
print("=" * 60)
print("Mastering optimization from fundamentals to advanced applications")
print(" Using Python Scientific Libraries")
print(" Ready to optimize!")
print()

Out[3]:

🚀 Linear Programming - Complete Tutorial
============================================================
Mastering optimization from fundamentals to advanced applications
📊 Using Python Scientific Libraries
📈 Ready to optimize!

Chapter 1: Linear Programming Fundamentals and Geometric Interpretation

We begin with a classic example that illustrates all key concepts of linear programming.

In [4]:

# Chapter 1: Basic Linear Programming and Geometric Interpretation
print(" CHAPTER 1: FUNDAMENTALS AND GEOMETRY")
print("=" * 50)

# Classic Product Mix Problem
print("=== PRODUCT MIX OPTIMIZATION ===")  
print("A furniture company produces chairs (x) and tables (y)")
print("Maximize profit: 3x + 5y dollars")
print("Subject to constraints:")
print("  Labor:     x + 2y ≤ 6 hours")
print("  Material:  2x + y ≤ 8 units")
print("  Non-negativity: x, y ≥ 0")
print()

# Solve using SageMath
p = MixedIntegerLinearProgram(maximization=True)
x = p.new_variable(real=True, nonnegative=True)
p.set_objective(3*x[0] + 5*x[1])  # 3*chairs + 5*tables
p.add_constraint(x[0] + 2*x[1] <= 6, name='labor')
p.add_constraint(2*x[0] + x[1] <= 8, name='material')

optimal_value = p.solve()
optimal_solution = p.get_values(x)
x_val, y_val = optimal_solution[0], optimal_solution[1]

print(f" OPTIMAL SOLUTION:")
print(f"  Maximum profit: ${optimal_value:.2f}")
print(f"  Optimal production: {x_val:.2f} chairs, {y_val:.2f} tables")
print()

# Verify solution manually
labor_used = x_val + 2*y_val
material_used = 2*x_val + y_val
print(f"📋 RESOURCE UTILIZATION:")
print(f"  Labor: {labor_used:.2f}/6 hours ({100*labor_used/6:.1f}% utilized)")
print(f"  Material: {material_used:.2f}/8 units ({100*material_used/8:.1f}% utilized)")
print()

Out[4]:

📈 CHAPTER 1: FUNDAMENTALS AND GEOMETRY
==================================================
=== PRODUCT MIX OPTIMIZATION ===
A furniture company produces chairs (x) and tables (y)
Maximize profit: 3x + 5y dollars
Subject to constraints:
  Labor:     x + 2y ≤ 6 hours
  Material:  2x + y ≤ 8 units
  Non-negativity: x, y ≥ 0

💰 OPTIMAL SOLUTION:
  Maximum profit: $16.67
  Optimal production: 3.33 chairs, 1.33 tables

📋 RESOURCE UTILIZATION:
  Labor: 6.00/6 hours (100.0% utilized)
  Material: 8.00/8 units (100.0% utilized)

In [12]:

import numpy as np
import matplotlib.pyplot as plt

print(" GEOMETRIC INTERPRETATION")
print("=" * 30)

# Create feasible region plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

# Left plot: Feasible region
x_range = np.linspace(0, 5, 400)
y_range = np.linspace(0, 4, 400)
X, Y = np.meshgrid(x_range, y_range)

# Define constraints (with a tiny tolerance)
feasible_mask = (X + 2*Y <= 6.01) & (2*X + Y <= 8.01) & (X >= 0) & (Y >= 0)

# Plot feasible region
ax1.contourf(X, Y, feasible_mask.astype(int), levels=[0.5, 1.5],
             colors=['lightblue'], alpha=0.7)

# Constraint lines
x_line = np.linspace(0, 5, 100)
ax1.plot(x_line, (6 - x_line) / 2, 'r-', linewidth=2, label='Labor: x + 2y = 6')
ax1.plot(x_line, 8 - 2*x_line, 'g-', linewidth=2, label='Material: 2x + y = 8')

# Corners as floats (avoid Sage Rationals)
corners = [(0.0, 0.0), (0.0, 3.0), (4.0, 0.0), (10/3.0, 4/3.0)]
corner_names = ['Origin', 'Labor Only', 'Material Only', 'Intersection']

# Try to use existing optimal point, else compute from corners
try:
    x_val_f = float(x_val)
    y_val_f = float(y_val)
    optimal_value_f = float(optimal_value)
except NameError:
    feas = []
    vals = []
    for (xc, yc) in corners:
        if (xc >= 0 and yc >= 0 and xc + 2*yc <= 6.001 and 2*xc + yc <= 8.001):
            feas.append((xc, yc))
            vals.append(3*xc + 5*yc)
    if vals:
        idx = int(np.argmax(vals))
        x_val_f, y_val_f = feas[idx]
        optimal_value_f = vals[idx]
    else:
        x_val_f, y_val_f, optimal_value_f = 0.0, 0.0, 0.0

# Mark optimal point
ax1.plot(x_val_f, y_val_f, 'ko', markersize=12, markerfacecolor='red',
         label=f'Optimal: ({x_val_f:.2f}, {y_val_f:.2f})')

# Objective function contours
for p in [10, 15, 18]:
    y_obj = (p - 3*x_line) / 5
    mask = (y_obj >= 0) & (y_obj <= 4)
    ax1.plot(x_line[mask], y_obj[mask], '--', alpha=0.6,
             label=f'Profit = ${p}' if p == 18 else '')

ax1.set_xlim(0, 5)
ax1.set_ylim(0, 4)
ax1.set_xlabel('Chairs (x)', fontsize=12)
ax1.set_ylabel('Tables (y)', fontsize=12)
ax1.set_title('Feasible Region and Optimal Solution', fontsize=14)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# Right plot: Corner point analysis
profits = []
feasible_corners = []
feasible_names = []

print(" CORNER POINT ANALYSIS:")
print(f"{'Corner':<15} {'Coordinates':<16} {'Feasible':<10} {'Profit':<10}")
print("-" * 60)

for (x_c, y_c), name in zip(corners, corner_names):
    x_f, y_f = float(x_c), float(y_c)
    labor_ok = x_f + 2*y_f <= 6.001
    material_ok = 2*x_f + y_f <= 8.001
    feasible_corner = (x_f >= 0) and (y_f >= 0) and labor_ok and material_ok

    if feasible_corner:
        profit = 3*x_f + 5*y_f
        profits.append(profit)
        feasible_corners.append((x_f, y_f))
        feasible_names.append(name)
        print(f"{name:<15} ({x_f:.2f}, {y_f:.2f})   {'✓':<10} ${profit:.2f}")

        color = 'red' if abs(profit - float(optimal_value_f)) < 0.01 else 'blue'
        size = 12 if abs(profit - float(optimal_value_f)) < 0.01 else 8
        ax2.scatter(x_f, y_f, c=color, s=size**2, alpha=0.8)
        ax2.annotate(f'{name}\n$({x_f:.1f}, {y_f:.1f})$\nProfit: ${profit:.1f}',
                     (x_f, y_f), xytext=(10, 10), textcoords='offset points',
                     fontsize=9, ha='left')
    else:
        print(f"{name:<15} ({x_f:.2f}, {y_f:.2f})   {'✗':<10} N/A")

print("\n FUNDAMENTAL THEOREM VERIFICATION:")
if profits:
    idx = int(np.argmax(profits))
    print(f"   Optimal value ${max(profits):.2f} occurs at vertex {feasible_names[idx]}")
    print(f"   This confirms: 'Linear programs achieve optimality at vertices!'")
else:
    print("   No feasible corners found.")

# Set up second plot
ax2.set_xlim(-0.5, 5)
ax2.set_ylim(-0.5, 4)
ax2.set_xlabel('Chairs (x)', fontsize=12)
ax2.set_ylabel('Tables (y)', fontsize=12)
ax2.set_title('Corner Point Analysis', fontsize=14)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()
print()

Out[12]:

🎨 GEOMETRIC INTERPRETATION
==============================
🔍 CORNER POINT ANALYSIS:
Corner          Coordinates      Feasible   Profit    
------------------------------------------------------------
Origin          (0.00, 0.00)   ✓          $0.00
Labor Only      (0.00, 3.00)   ✓          $15.00
Material Only   (4.00, 0.00)   ✓          $12.00
Intersection    (3.33, 1.33)   ✓          $16.67

🏆 FUNDAMENTAL THEOREM VERIFICATION:
   Optimal value $16.67 occurs at vertex Intersection
   This confirms: 'Linear programs achieve optimality at vertices!'

Chapter 2: Duality Theory and Sensitivity Analysis

Every linear program has a "dual" problem that provides profound economic insights and optimality conditions.

In [6]:

# Chapter 2: Duality Theory and Economic Interpretation
print("🔄 CHAPTER 2: DUALITY THEORY")
print("=" * 40)

# Production Planning Example (Primal Problem)
print("=== PRIMAL PROBLEM: Production Planning ====")
print("Furniture company wants to maximize profit:")
print("Maximize: 40A + 30B dollars (profit from products A and B)")
print("Subject to:")
print("  Labor constraint:    2A + 1B ≤ 100 hours")
print("  Material constraint: 1A + 2B ≤ 80 units")
print("  Non-negativity: A, B ≥ 0")
print()

# Solve primal problem
primal = MixedIntegerLinearProgram(maximization=True)
x = primal.new_variable(real=True, nonnegative=True)
primal.set_objective(40*x[0] + 30*x[1])  # Profit function
primal.add_constraint(2*x[0] + 1*x[1] <= 100)  # Labor constraint
primal.add_constraint(1*x[0] + 2*x[1] <= 80)   # Material constraint

primal_value = primal.solve()
primal_solution = primal.get_values(x)
A_opt, B_opt = primal_solution[0], primal_solution[1]

print(f" PRIMAL OPTIMAL SOLUTION:")
print(f"  Maximum profit: ${primal_value:.2f}")
print(f"  Optimal production: A = {A_opt:.1f}, B = {B_opt:.1f}")

# Check resource utilization
labor_used = 2*A_opt + B_opt
material_used = A_opt + 2*B_opt
print(f"  Labor used: {labor_used:.1f}/100 hours")
print(f"  Material used: {material_used:.1f}/80 units")
print()

print("=== DUAL PROBLEM: Resource Valuation ====")
print("What is the value of each resource hour/unit?")
print("Minimize: 100λ + 80μ dollars (total resource value)")
print("Subject to:")
print("  Product A constraint: 2λ + 1μ ≥ 40 (value ≥ profit)")
print("  Product B constraint: 1λ + 2μ ≥ 30 (value ≥ profit)")
print("  Non-negativity: λ, μ ≥ 0")
print()

# Solve dual problem
dual = MixedIntegerLinearProgram(maximization=False)  # Minimize
y = dual.new_variable(real=True, nonnegative=True)
dual.set_objective(100*y[0] + 80*y[1])  # Resource valuation
dual.add_constraint(2*y[0] + 1*y[1] >= 40)  # For product A
dual.add_constraint(1*y[0] + 2*y[1] >= 30)  # For product B

dual_value = dual.solve()
dual_solution = dual.get_values(y)
lambda_opt, mu_opt = dual_solution[0], dual_solution[1]

print(f" DUAL OPTIMAL SOLUTION:")
print(f"  Minimum resource cost: ${dual_value:.2f}")
print(f"  Shadow prices: Labor λ = ${lambda_opt:.2f}/hour")
print(f"                Material μ = ${mu_opt:.2f}/unit")
print()

print("=== STRONG DUALITY THEOREM ====")
print(f" Fundamental Result: Primal = Dual at optimality")
print(f"   Primal optimal value:  ${primal_value:.6f}")
print(f"   Dual optimal value:    ${dual_value:.6f}")
print(f"   Duality gap:           ${abs(primal_value - dual_value):.8f} ≈ 0 ✓")
print(f"   This confirms optimal solutions are correct!")
print()

Out[6]:

🔄 CHAPTER 2: DUALITY THEORY
========================================
=== PRIMAL PROBLEM: Production Planning ====
Furniture company wants to maximize profit:
Maximize: 40A + 30B dollars (profit from products A and B)
Subject to:
  Labor constraint:    2A + 1B ≤ 100 hours
  Material constraint: 1A + 2B ≤ 80 units
  Non-negativity: A, B ≥ 0

💰 PRIMAL OPTIMAL SOLUTION:
  Maximum profit: $2200.00
  Optimal production: A = 40.0, B = 20.0
  Labor used: 100.0/100 hours
  Material used: 80.0/80 units

=== DUAL PROBLEM: Resource Valuation ====
What is the value of each resource hour/unit?
Minimize: 100λ + 80μ dollars (total resource value)
Subject to:
  Product A constraint: 2λ + 1μ ≥ 40 (value ≥ profit)
  Product B constraint: 1λ + 2μ ≥ 30 (value ≥ profit)
  Non-negativity: λ, μ ≥ 0

💎 DUAL OPTIMAL SOLUTION:
  Minimum resource cost: $2200.00
  Shadow prices: Labor λ = $16.67/hour
                Material μ = $6.67/unit

=== STRONG DUALITY THEOREM ====
🎯 Fundamental Result: Primal = Dual at optimality
   Primal optimal value:  $2200.000000
   Dual optimal value:    $2200.000000
   Duality gap:           $0.00000000 ≈ 0 ✓
   This confirms optimal solutions are correct!

In [7]:

# Sensitivity Analysis using Shadow Prices
print(" SENSITIVITY ANALYSIS")
print("=" * 30)
print("Shadow prices predict value of additional resources")
print()

# Test shadow price predictions
lambda_shadow = dual_solution[0]  # Labor shadow price
mu_shadow = dual_solution[1]      # Material shadow price

print(f" SHADOW PRICE INTERPRETATION:")
print(f"   Labor shadow price λ = ${lambda_shadow:.2f}/hour")
print(f"   → Each additional labor hour increases profit by ${lambda_shadow:.2f}")
print(f"   Material shadow price μ = ${mu_shadow:.2f}/unit")
print(f"   → Each additional material unit increases profit by ${mu_shadow:.2f}")
print()

# Test predictions with +1 labor hour
test_labor = MixedIntegerLinearProgram(maximization=True)
x_test = test_labor.new_variable(real=True, nonnegative=True)
test_labor.set_objective(40*x_test[0] + 30*x_test[1])
test_labor.add_constraint(2*x_test[0] + 1*x_test[1] <= 101)  # +1 labor
test_labor.add_constraint(1*x_test[0] + 2*x_test[1] <= 80)
new_labor_value = test_labor.solve()

print(f" TESTING PREDICTIONS:")
print(f"Adding 1 labor hour (100 → 101):")
print(f"   Shadow price prediction: +${lambda_shadow:.3f}")
print(f"   Actual improvement:      +${new_labor_value - primal_value:.3f}")
print(f"   Prediction accuracy:     {abs((new_labor_value - primal_value) - lambda_shadow) < 0.001} ✓")
print()

# Test with +1 material unit
test_material = MixedIntegerLinearProgram(maximization=True)
x_test2 = test_material.new_variable(real=True, nonnegative=True)
test_material.set_objective(40*x_test2[0] + 30*x_test2[1])
test_material.add_constraint(2*x_test2[0] + 1*x_test2[1] <= 100)
test_material.add_constraint(1*x_test2[0] + 2*x_test2[1] <= 81)  # +1 material
new_material_value = test_material.solve()

print(f"Adding 1 material unit (80 → 81):")
print(f"   Shadow price prediction: +${mu_shadow:.3f}")
print(f"   Actual improvement:      +${new_material_value - primal_value:.3f}")
print(f"   Prediction accuracy:     {abs((new_material_value - primal_value) - mu_shadow) < 0.001} ✓")
print()

print(" ECONOMIC INSIGHTS:")
if lambda_shadow > mu_shadow:
    print(f"   Labor (${lambda_shadow:.2f}) is more valuable than material (${mu_shadow:.2f})")
    print(f"   → Company should prioritize acquiring more labor capacity")
else:
    print(f"   Material (${mu_shadow:.2f}) is more valuable than labor (${lambda_shadow:.2f})")
    print(f"   → Company should prioritize acquiring more material")

print(f"   Zero shadow price means resource has slack (unused capacity)")
print(f"   Positive shadow price means constraint is binding (fully utilized)")
print()

Out[7]:

🔍 SENSITIVITY ANALYSIS
==============================
Shadow prices predict value of additional resources

📊 SHADOW PRICE INTERPRETATION:
   Labor shadow price λ = $16.67/hour
   → Each additional labor hour increases profit by $16.67
   Material shadow price μ = $6.67/unit
   → Each additional material unit increases profit by $6.67

🧪 TESTING PREDICTIONS:
Adding 1 labor hour (100 → 101):
   Shadow price prediction: +$16.667
   Actual improvement:      +$16.667
   Prediction accuracy:     True ✓

Adding 1 material unit (80 → 81):
   Shadow price prediction: +$6.667
   Actual improvement:      +$6.667
   Prediction accuracy:     True ✓

💡 ECONOMIC INSIGHTS:
   Labor ($16.67) is more valuable than material ($6.67)
   → Company should prioritize acquiring more labor capacity
   Zero shadow price means resource has slack (unused capacity)
   Positive shadow price means constraint is binding (fully utilized)

Chapter 3: Transportation Problems

Transportation problems are a special class of linear programs with applications in logistics, supply chain, and distribution.

In [8]:

# Chapter 3: Transportation Problem - Supply Chain Optimization
print("🚚 CHAPTER 3: TRANSPORTATION PROBLEMS")
print("=" * 45)
print("Optimize shipping costs from suppliers to customers")
print()

# Problem setup - Electronics distribution
suppliers = ['Factory A', 'Factory B', 'Factory C']
customers = ['Store 1', 'Store 2', 'Store 3', 'Store 4']
supplies = [300, 400, 350]  # Units available
demands = [200, 250, 300, 250]  # Units required

# Shipping costs per unit ($)
costs = [
    [8, 6, 10, 9],   # Factory A to Stores 1,2,3,4
    [9, 12, 13, 7],  # Factory B to Stores 1,2,3,4  
    [14, 9, 16, 5]   # Factory C to Stores 1,2,3,4
]

print("=== SUPPLY AND DEMAND ANALYSIS ====")
print(f"Suppliers: {len(suppliers)} factories")
for i, (supplier, supply) in enumerate(zip(suppliers, supplies)):
    print(f"  {supplier}: {supply} units available")

print(f"\nCustomers: {len(customers)} stores")
for i, (customer, demand) in enumerate(zip(customers, demands)):
    print(f"  {customer}: {demand} units needed")

total_supply = sum(supplies)
total_demand = sum(demands)
print(f"\n BALANCE CHECK:")
print(f"   Total supply: {total_supply} units")
print(f"   Total demand: {total_demand} units")
print(f"   Problem is {'balanced' if total_supply == total_demand else 'unbalanced'} ✓")
print()

# Display cost matrix
print("=== TRANSPORTATION COST MATRIX ($/unit) ====")
print(f"{'Supplier':<12}", end="")
for customer in customers:
    print(f"{customer:>8}", end="")
print(f"{'Supply':>8}")
print("-" * 60)

for i, (supplier, supply) in enumerate(zip(suppliers, supplies)):
    print(f"{supplier:<12}", end="")
    for j in range(len(customers)):
        print(f"${costs[i][j]:>6}", end="")
    print(f"{supply:>8}")
    
print("-" * 60)
print(f"{'Demand':<12}", end="")
for demand in demands:
    print(f"{demand:>7}", end="")
print()
print()

Out[8]:

🚚 CHAPTER 3: TRANSPORTATION PROBLEMS
=============================================
Optimize shipping costs from suppliers to customers

=== SUPPLY AND DEMAND ANALYSIS ====
Suppliers: 3 factories
  Factory A: 300 units available
  Factory B: 400 units available
  Factory C: 350 units available

Customers: 4 stores
  Store 1: 200 units needed
  Store 2: 250 units needed
  Store 3: 300 units needed
  Store 4: 250 units needed

📊 BALANCE CHECK:
   Total supply: 1050 units
   Total demand: 1000 units
   Problem is unbalanced ✓

=== TRANSPORTATION COST MATRIX ($/unit) ====
Supplier     Store 1 Store 2 Store 3 Store 4  Supply
------------------------------------------------------------
Factory A   $     8$     6$    10$     9     300
Factory B   $     9$    12$    13$     7     400
Factory C   $    14$     9$    16$     5     350
------------------------------------------------------------
Demand          200    250    300    250

In [21]:

# Solve transportation problem
print(" SOLVING TRANSPORTATION PROBLEM")
print("=" * 40)

# Set up the linear program
transport = MixedIntegerLinearProgram(maximization=False)  # Minimize cost
x = transport.new_variable(real=True, nonnegative=True)

# Objective: minimize total transportation cost
objective = sum(costs[i][j] * x[i,j] 
               for i in range(len(supplies)) 
               for j in range(len(demands)))
transport.set_objective(objective)

# Supply: ship at most supply
for i in range(len(supplies)):
    transport.add_constraint(sum(x[i,j] for j in range(len(demands))) <= supplies[i])

# Demand: must meet each demand exactly
for j in range(len(demands)):
    transport.add_constraint(sum(x[i,j] for i in range(len(supplies))) == demands[j])

# Solve the problem
min_cost = transport.solve()
solution = transport.get_values(x)

print(f" OPTIMAL SOLUTION FOUND!")
print(f"   Minimum total cost: ${min_cost:.2f}")
print()

# Display optimal shipping plan
print("=== OPTIMAL SHIPPING PLAN ====")
print(f"{'Supplier':<12}", end="")
for customer in customers:
    print(f"{customer:>8}", end="")
print(f"{'Total':>8}")
print("-" * 60)

total_shipped = 0
for i in range(len(supplies)):
    print(f"{suppliers[i]:<12}", end="")
    row_total = 0
    for j in range(len(demands)):
        amount = solution.get((i,j), 0)
        if amount > 0.001:  # Only show non-zero shipments
            print(f"{amount:>7.0f}", end="")
        else:
            print(f"{'':>7}", end="")
        row_total += amount
    print(f"{row_total:>8.0f}")
    total_shipped += row_total

print("-" * 60)
print(f"{'Received':<12}", end="")
for j in range(len(demands)):
    col_total = sum(solution.get((i,j), 0) for i in range(len(supplies)))
    print(f"{col_total:>7.0f}", end="")
print(f"{total_shipped:>8.0f}")
print()

# Analyze solution properties
non_zero_routes = sum(1 for val in solution.values() if val > 0.001)
total_variables = len(supplies) * len(demands)
basic_variables_expected = len(supplies) + len(demands) - 1

print("=== SOLUTION ANALYSIS ====")
print(f" Routes used: {non_zero_routes}")
print(f" Total possible routes: {total_variables}")
print(f" Basic variables (theory): m + n - 1 = {basic_variables_expected}")
print(f"✅ Solution structure: {'Non-degenerate' if non_zero_routes == basic_variables_expected else 'Degenerate'}")

# Calculate cost savings
naive_cost = sum(max(costs[i]) * supplies[i] for i in range(len(supplies)))
savings = naive_cost - min_cost
savings_percent = 100 * savings / naive_cost

print(f"\n ECONOMIC IMPACT:")
print(f"   Naive highest-cost strategy: ${naive_cost:.2f}")
print(f"   Optimized transportation: ${min_cost:.2f}")
print(f"   Total savings: ${savings:.2f} ({savings_percent:.1f}%)")
print()

Out[21]:

🎯 SOLVING TRANSPORTATION PROBLEM
========================================
💰 OPTIMAL SOLUTION FOUND!
   Minimum total cost: $8300.00

=== OPTIMAL SHIPPING PLAN ====
Supplier     Store 1 Store 2 Store 3 Store 4   Total
------------------------------------------------------------
Factory A              200    100            300
Factory B       200           200            400
Factory C               50           250     300
------------------------------------------------------------
Received        200    250    300    250    1000

=== SOLUTION ANALYSIS ====
📈 Routes used: 6
📊 Total possible routes: 12
🧮 Basic variables (theory): m + n - 1 = 6
✅ Solution structure: Non-degenerate

💡 ECONOMIC IMPACT:
   Naive highest-cost strategy: $13800.00
   Optimized transportation: $8300.00
   Total savings: $5500.00 (39.9%)

Chapter 4: Integer Programming and Discrete Optimization

When decision variables must be integers, we enter the realm of integer programming - computationally more challenging but essential for many real-world problems.

In [10]:

# Chapter 4: Integer Programming - The Knapsack Problem
print("🎒 CHAPTER 4: INTEGER PROGRAMMING")
print("=" * 40)
print("Discrete optimization with 0-1 decision variables")
print()

# Classic 0-1 Knapsack Problem
print("=== THE KNAPSACK PROBLEM ====")
print("A hiker must choose items for a backpacking trip")
print("Each item has weight and value - maximize utility!")
print()

# Problem data
items = ['Water Filter', 'Sleeping Bag', 'Tent', 'Cook Set', 'First Aid', 'Headlamp']
weights = [2, 4, 6, 3, 1, 1]  # kg
values = [20, 15, 25, 12, 8, 6]   # utility points
capacity = 10  # kg capacity

print(f"🎒 Backpack capacity: {capacity} kg")
print("\n📋 Available items:")
print(f"{'Item':<15} {'Weight':<8} {'Value':<8} {'Ratio':<10}")
print("-" * 45)

item_data = []
for i, item in enumerate(items):
    ratio = values[i] / weights[i]  # Sage Rational
    item_data.append((item, weights[i], values[i], ratio))
    print(f"{item:<15} {weights[i]:<8} {values[i]:<8} {float(ratio):<10.2f}")

print()

# Sort by value-to-weight ratio for comparison
sorted_items = sorted(enumerate(item_data), key=lambda x: x[1][3], reverse=True)
print(" GREEDY RANKING (by value/weight ratio):")
for rank, (idx, (item, weight, value, ratio)) in enumerate(sorted_items, 1):
    print(f"{rank}. {item} (ratio: {float(ratio):.2f})")
print()

Out[10]:

🎒 CHAPTER 4: INTEGER PROGRAMMING
========================================
Discrete optimization with 0-1 decision variables

=== THE KNAPSACK PROBLEM ====
A hiker must choose items for a backpacking trip
Each item has weight and value - maximize utility!

🎒 Backpack capacity: 10 kg

📋 Available items:
Item            Weight   Value    Ratio     
---------------------------------------------

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[10], line 28
     26     ratio = values[i] / weights[i]
     27     item_data.append((item, weights[i], values[i], ratio))
---> 28     print(f"{item:<15} {weights[i]:<8} {values[i]:<8} {ratio:<10.2f}")
     30 print()
     32 # Sort by value-to-weight ratio for comparison
TypeError: unsupported format string passed to sage.rings.rational.Rational.__format__

In [23]:

# Solve 0-1 Knapsack using Integer Programming
print(" SOLVING WITH INTEGER PROGRAMMING")
print("=" * 40)

# Set up integer program
knapsack = MixedIntegerLinearProgram(maximization=True)
x = knapsack.new_variable(binary=True)  # 0-1 variables

# Objective: maximize total value
knapsack.set_objective(sum(values[i] * x[i] for i in range(len(items))))

# Constraint: respect weight capacity
knapsack.add_constraint(sum(weights[i] * x[i] for i in range(len(items))) <= capacity)

# Solve the integer program
max_value = knapsack.solve()
selection = knapsack.get_values(x)

print(f" OPTIMAL INTEGER SOLUTION:")
print(f"   Maximum utility: {max_value:.0f} points")

total_weight = 0
selected_items = []
print(f"\n🎒 Selected items:")
for i, item in enumerate(items):
    if selection[i] > 0.5:  # Binary variable is 1
        selected_items.append(item)
        print(f"   ✓ {item} (weight: {weights[i]} kg, value: {values[i]} points)")
        total_weight += weights[i]

print(f"\n SOLUTION SUMMARY:")
percent = 100.0 * float(total_weight) / float(capacity)
print(f"   Total weight: {total_weight}/{capacity} kg ({percent:.1f}% of capacity)")
print(f"   Items selected: {len(selected_items)}/{len(items)}")
print()

Out[23]:

🧮 SOLVING WITH INTEGER PROGRAMMING
========================================
🏆 OPTIMAL INTEGER SOLUTION:
   Maximum utility: 59 points

🎒 Selected items:
   ✓ Water Filter (weight: 2 kg, value: 20 points)
   ✓ Tent (weight: 6 kg, value: 25 points)
   ✓ First Aid (weight: 1 kg, value: 8 points)
   ✓ Headlamp (weight: 1 kg, value: 6 points)

📊 SOLUTION SUMMARY:
   Total weight: 10/10 kg (100.0% of capacity)
   Items selected: 4/6

In [14]:

# Compare with LP Relaxation
print(" LP RELAXATION ANALYSIS")
print("=" * 30)
print("What if we allow fractional items? (Upper bound for integer solution)")
print()

# Solve LP relaxation (allow fractional variables)
knapsack_relaxed = MixedIntegerLinearProgram(maximization=True)
y = knapsack_relaxed.new_variable(real=True, nonnegative=True)

# Same objective and weight constraint
knapsack_relaxed.set_objective(sum(values[i] * y[i] for i in range(len(items))))
knapsack_relaxed.add_constraint(sum(weights[i] * y[i] for i in range(len(items))) <= capacity)

# Add upper bounds (can't take more than 1 of each item)
for i in range(len(items)):
    knapsack_relaxed.add_constraint(y[i] <= 1)

relaxed_value = knapsack_relaxed.solve()
relaxed_solution = knapsack_relaxed.get_values(y)

print(f" LP RELAXATION RESULTS:")
print(f"   Relaxed optimal value: {relaxed_value:.3f} points")
print(f"   Integer optimal value: {max_value:.0f} points")
print(f"   Integrality gap: {relaxed_value - max_value:.3f} points")
print(f"   Gap percentage: {100*(relaxed_value - max_value)/relaxed_value:.2f}%")
print()

print(f" FRACTIONAL SOLUTION ANALYSIS:")
for i, item in enumerate(items):
    frac_amount = relaxed_solution[i]
    if frac_amount > 0.001:
        if frac_amount >= 0.999:
            print(f"   {item}: {frac_amount:.3f} (take fully)")
        else:
            print(f"   {item}: {frac_amount:.3f} (partial - {100*frac_amount:.1f}%)")

print()
print(" KEY INSIGHTS:")
print("   ✓ LP relaxation provides upper bound for integer solution")
print("   ✓ Fractional variables reveal which items are 'on the margin'")
print("   ✓ Small integrality gap indicates good integer solution quality")
print()

Out[14]:

🔍 LP RELAXATION ANALYSIS
==============================
What if we allow fractional items? (Upper bound for integer solution)

📈 LP RELAXATION RESULTS:
   Relaxed optimal value: 59.000 points
   Integer optimal value: 59 points
   Integrality gap: 0.000 points
   Gap percentage: 0.00%

🔬 FRACTIONAL SOLUTION ANALYSIS:
   Water Filter: 1.000 (take fully)
   Tent: 1.000 (take fully)
   First Aid: 1.000 (take fully)
   Headlamp: 1.000 (take fully)

💡 KEY INSIGHTS:
   ✓ LP relaxation provides upper bound for integer solution
   ✓ Fractional variables reveal which items are 'on the margin'
   ✓ Small integrality gap indicates good integer solution quality

Chapter 5: Portfolio Optimization and Financial Applications

Linear programming revolutionized modern finance, from portfolio theory to risk management.

In [15]:

# Chapter 5: Portfolio Optimization and Financial Applications
print("💼 CHAPTER 5: PORTFOLIO OPTIMIZATION")
print("=" * 42)
print("Modern Portfolio Theory meets Linear Programming")
print()

# Investment universe
assets = ['Tech Stocks', 'Bank Stocks', 'Gov Bonds', 'Corp Bonds', 'REITs', 'Commodities']
expected_returns = [0.15, 0.12, 0.04, 0.06, 0.10, 0.08]  # Annual expected return
risk_levels = [0.25, 0.20, 0.02, 0.08, 0.15, 0.18]      # Standard deviation (risk)
min_weights = [0.05, 0.05, 0.10, 0.05, 0.00, 0.00]     # Minimum allocations
max_weights = [0.30, 0.25, 0.50, 0.30, 0.20, 0.15]     # Maximum allocations

print("=== INVESTMENT UNIVERSE ====")
print(f"{'Asset Class':<15} {'Return':<8} {'Risk':<8} {'Min%':<6} {'Max%':<6}")
print("-" * 50)
for i, asset in enumerate(assets):
    print(f"{asset:<15} {100*expected_returns[i]:>6.1f}% {100*risk_levels[i]:>6.1f}% {100*min_weights[i]:>4.0f}% {100*max_weights[i]:>4.0f}%")
print()

# Problem 1: Maximize return subject to risk constraint
print("=== PROBLEM 1: MAXIMIZE RETURN (Risk ≤ 12%) ====")
risk_limit = 0.12

portfolio1 = MixedIntegerLinearProgram(maximization=True)
w = portfolio1.new_variable(real=True, nonnegative=True)

# Objective: maximize expected return
portfolio1.set_objective(sum(expected_returns[i] * w[i] for i in range(len(assets))))

# Constraints
portfolio1.add_constraint(sum(w[i] for i in range(len(assets))) == 1)  # Fully invested
portfolio1.add_constraint(sum(risk_levels[i] * w[i] for i in range(len(assets))) <= risk_limit)  # Risk limit

# Individual asset limits
for i in range(len(assets)):
    portfolio1.add_constraint(w[i] >= min_weights[i])  # Minimum weights
    portfolio1.add_constraint(w[i] <= max_weights[i])  # Maximum weights

max_return = portfolio1.solve()
optimal_weights = portfolio1.get_values(w)

print(f" OPTIMAL PORTFOLIO (Risk-Constrained):")
print(f"   Expected annual return: {100*max_return:.2f}%")

portfolio_risk = sum(risk_levels[i] * optimal_weights[i] for i in range(len(assets)))
print(f"   Portfolio risk: {100*portfolio_risk:.2f}% (limit: {100*risk_limit:.1f}%)")
print(f"\n Asset Allocation:")

for i, asset in enumerate(assets):
    weight = optimal_weights[i]
    if weight > 0.001:  # Only show meaningful allocations
        contribution = weight * expected_returns[i]
        print(f"   {asset:<15} {100*weight:>6.1f}% (contributes {100*contribution:>4.2f}% return)")

print()

# Verify constraints
total_weight = sum(optimal_weights.values())
active_constraints = []
if abs(portfolio_risk - risk_limit) < 0.001:
    active_constraints.append("Risk constraint")
for i, asset in enumerate(assets):
    if abs(optimal_weights[i] - max_weights[i]) < 0.001:
        active_constraints.append(f"{asset} upper bound")

print(f"✅ Constraint verification: Total weight = {total_weight:.3f}")
print(f"🔒 Active constraints: {', '.join(active_constraints) if active_constraints else 'None'}")
print()

Out[15]:

💼 CHAPTER 5: PORTFOLIO OPTIMIZATION
==========================================
Modern Portfolio Theory meets Linear Programming

=== INVESTMENT UNIVERSE ====
Asset Class     Return   Risk     Min%   Max%  
--------------------------------------------------
Tech Stocks       15.0%   25.0%    5%   30%
Bank Stocks       12.0%   20.0%    5%   25%
Gov Bonds          4.0%    2.0%   10%   50%
Corp Bonds         6.0%    8.0%    5%   30%
REITs             10.0%   15.0%    0%   20%
Commodities        8.0%   18.0%    0%   15%

=== PROBLEM 1: MAXIMIZE RETURN (Risk ≤ 12%) ====
🎯 OPTIMAL PORTFOLIO (Risk-Constrained):
   Expected annual return: 8.68%
   Portfolio risk: 12.00% (limit: 12.0%)

📊 Asset Allocation:
   Tech Stocks       30.0% (contributes 4.50% return)
   Bank Stocks        5.0% (contributes 0.60% return)
   Gov Bonds         45.4% (contributes 1.82% return)
   Corp Bonds         5.0% (contributes 0.30% return)
   REITs             14.6% (contributes 1.46% return)

✅ Constraint verification: Total weight = 1.000
🔒 Active constraints: Risk constraint, Tech Stocks upper bound

In [16]:

# Problem 2: Efficient Frontier Analysis
print(" EFFICIENT FRONTIER ANALYSIS")
print("=" * 35)
print("Generate optimal risk-return combinations")
print()

# Generate efficient frontier points
risk_targets = [0.05, 0.08, 0.10, 0.12, 0.15, 0.18]
frontier_points = []

print(f"{'Target Risk':<12} {'Max Return':<12} {'Portfolio Composition':<30}")
print("-" * 70)

for risk_target in risk_targets:
    try:
        # Solve for each risk level
        portfolio_ef = MixedIntegerLinearProgram(maximization=True)
        w_ef = portfolio_ef.new_variable(real=True, nonnegative=True)
        
        portfolio_ef.set_objective(sum(expected_returns[i] * w_ef[i] for i in range(len(assets))))
        portfolio_ef.add_constraint(sum(w_ef[i] for i in range(len(assets))) == 1)
        portfolio_ef.add_constraint(sum(risk_levels[i] * w_ef[i] for i in range(len(assets))) <= risk_target)
        
        for i in range(len(assets)):
            portfolio_ef.add_constraint(w_ef[i] >= min_weights[i])
            portfolio_ef.add_constraint(w_ef[i] <= max_weights[i])
        
        return_ef = portfolio_ef.solve()
        weights_ef = portfolio_ef.get_values(w_ef)
        
        # Find dominant asset
        max_weight = max(weights_ef.values())
        dominant_asset = assets[list(weights_ef.values()).index(max_weight)]
        
        frontier_points.append((risk_target, return_ef))
        print(f"{100*risk_target:>10.1f}% {100*return_ef:>10.2f}% {dominant_asset} dominant ({100*max_weight:.1f}%)")
        
    except:
        print(f"{100*risk_target:>10.1f}% {'Infeasible':<10}")

print()

# Problem 3: Minimum Risk Portfolio
print("=== PROBLEM 3: MINIMUM RISK PORTFOLIO ====")
portfolio_minrisk = MixedIntegerLinearProgram(maximization=False)
w_mr = portfolio_minrisk.new_variable(real=True, nonnegative=True)

# Objective: minimize risk
portfolio_minrisk.set_objective(sum(risk_levels[i] * w_mr[i] for i in range(len(assets))))
portfolio_minrisk.add_constraint(sum(w_mr[i] for i in range(len(assets))) == 1)

for i in range(len(assets)):
    portfolio_minrisk.add_constraint(w_mr[i] >= min_weights[i])
    portfolio_minrisk.add_constraint(w_mr[i] <= max_weights[i])

min_risk = portfolio_minrisk.solve()
minrisk_weights = portfolio_minrisk.get_values(w_mr)
minrisk_return = sum(expected_returns[i] * minrisk_weights[i] for i in range(len(assets)))

print(f"🛡 MINIMUM RISK PORTFOLIO:")
print(f"   Minimum risk: {100*min_risk:.2f}%")
print(f"   Expected return: {100*minrisk_return:.2f}%")
print(f"   Return per unit risk: {minrisk_return/min_risk:.2f}")

print(f"\n Conservative allocation:")
for i, asset in enumerate(assets):
    weight = minrisk_weights[i]
    if weight > 0.001:
        print(f"   {asset:<15} {100*weight:>6.1f}%")

print()

# Summary and insights
print(" PORTFOLIO OPTIMIZATION INSIGHTS:")
print(f"    Risk-return tradeoff: Higher risk enables higher expected return")
print(f"    Diversification: Optimal portfolios spread risk across asset classes")
print(f"    Constraints matter: Position limits prevent over-concentration")
print(f"    Efficient frontier: Shows best possible risk-return combinations")
print(f"   🛡 Risk management: LP enables precise risk targeting")
print()

Out[16]:

📈 EFFICIENT FRONTIER ANALYSIS
===================================
Generate optimal risk-return combinations

Target Risk  Max Return   Portfolio Composition         
----------------------------------------------------------------------
       5.0% Infeasible
       8.0%       6.63% Gov Bonds dominant (50.0%)
      10.0%       7.69% Gov Bonds dominant (50.0%)
      12.0%       8.68% Gov Bonds dominant (45.4%)
      15.0%      10.02% Tech Stocks dominant (30.0%)
      18.0%      11.00% Tech Stocks dominant (30.0%)

=== PROBLEM 3: MINIMUM RISK PORTFOLIO ====
🛡️ MINIMUM RISK PORTFOLIO:
   Minimum risk: 7.15%
   Expected return: 6.15%
   Return per unit risk: 0.86

📊 Conservative allocation:
   Tech Stocks        5.0%
   Bank Stocks        5.0%
   Gov Bonds         50.0%
   Corp Bonds        30.0%
   REITs             10.0%

💡 PORTFOLIO OPTIMIZATION INSIGHTS:
   📊 Risk-return tradeoff: Higher risk enables higher expected return
   🎯 Diversification: Optimal portfolios spread risk across asset classes
   ⚖️ Constraints matter: Position limits prevent over-concentration
   📈 Efficient frontier: Shows best possible risk-return combinations
   🛡️ Risk management: LP enables precise risk targeting

Chapter 6: Advanced Applications and Real-World Case Studies

Explore sophisticated linear programming applications across industries.

In [17]:

# Chapter 6: Advanced Applications - Multi-Period Production Planning
print("🏭 CHAPTER 6: ADVANCED APPLICATIONS")
print("=" * 42)
print("Multi-period production planning with inventory")
print()

# Multi-period production planning problem
print("=== SEMICONDUCTOR MANUFACTURING PLANNING ====")
print("Plan production over 4 quarters to meet varying demand")
print()

# Problem parameters
periods = 4  # quarters
products = ['CPU', 'GPU', 'Memory']
demand = {
    'CPU': [1000, 1200, 1500, 1300],
    'GPU': [800, 1000, 1100, 900], 
    'Memory': [2000, 2200, 2500, 2300]
}

production_cost = {'CPU': 200, 'GPU': 350, 'Memory': 80}  # $/unit
inventory_cost = {'CPU': 10, 'GPU': 20, 'Memory': 5}     # $/unit/quarter
production_capacity = {'CPU': 1400, 'GPU': 1200, 'Memory': 2600}  # units/quarter
initial_inventory = {'CPU': 100, 'GPU': 50, 'Memory': 200}

print(" DEMAND FORECAST:")
print(f"{'Product':<8}", end="")
for q in range(periods):
    print(f"Q{q+1:>7}", end="")
print(f"{'Total':<8}")
print("-" * 45)

for product in products:
    print(f"{product:<8}", end="")
    total_demand = 0
    for q in range(periods):
        print(f"{demand[product][q]:>7}", end="")
        total_demand += demand[product][q]
    print(f"{total_demand:<8}")
print()

print(" COST STRUCTURE:")
print(f"{'Product':<8} {'Production':<12} {'Inventory':<12} {'Capacity':<12}")
print("-" * 50)
for product in products:
    print(f"{product:<8} ${production_cost[product]:<11} ${inventory_cost[product]:<11} {production_capacity[product]:<12}")
print()

Out[17]:

🏭 CHAPTER 6: ADVANCED APPLICATIONS
==========================================
Multi-period production planning with inventory

=== SEMICONDUCTOR MANUFACTURING PLANNING ====
Plan production over 4 quarters to meet varying demand

📊 DEMAND FORECAST:
Product Q      1Q      2Q      3Q      4Total   
---------------------------------------------
CPU        1000   1200   1500   13005000    
GPU         800   1000   1100    9003800    
Memory     2000   2200   2500   23009000    

💰 COST STRUCTURE:
Product  Production   Inventory    Capacity    
--------------------------------------------------
CPU      $200         $10          1400        
GPU      $350         $20          1200        
Memory   $80          $5           2600        

In [18]:

# Solve multi-period production planning
print(" MULTI-PERIOD OPTIMIZATION")
print("=" * 32)

# Set up the optimization model
production_plan = MixedIntegerLinearProgram(maximization=False)  # Minimize cost

# Decision variables
produce = production_plan.new_variable(real=True, nonnegative=True)  # produce[p,t]
inventory = production_plan.new_variable(real=True, nonnegative=True)  # inventory[p,t]

# Objective: minimize total production and inventory costs
total_cost = 0
for p, product in enumerate(products):
    for t in range(periods):
        total_cost += production_cost[product] * produce[p,t]  # Production cost
        total_cost += inventory_cost[product] * inventory[p,t]  # Inventory holding cost

production_plan.set_objective(total_cost)

# Constraints
for p, product in enumerate(products):
    for t in range(periods):
        # Production capacity constraints
        production_plan.add_constraint(produce[p,t] <= production_capacity[product])
        
        # Inventory balance constraints
        if t == 0:
            # First period: initial inventory + production - demand = ending inventory
            inventory_balance = (initial_inventory[product] + produce[p,t] - 
                               demand[product][t] == inventory[p,t])
        else:
            # Other periods: previous inventory + production - demand = ending inventory
            inventory_balance = (inventory[p,t-1] + produce[p,t] - 
                               demand[product][t] == inventory[p,t])
        
        production_plan.add_constraint(inventory_balance)

# Solve the problem
optimal_cost = production_plan.solve()
production_solution = production_plan.get_values(produce)
inventory_solution = production_plan.get_values(inventory)

print(f" OPTIMAL SOLUTION FOUND!")
print(f"   Minimum total cost: ${optimal_cost:,.0f}")
print()

# Display production plan
print("=== OPTIMAL PRODUCTION PLAN ====")
print(f"{'Product':<8}", end="")
for q in range(periods):
    print(f"Q{q+1} Prod", end="  ")
print()
print("-" * 45)

total_production_cost = 0
for p, product in enumerate(products):
    print(f"{product:<8}", end="")
    for t in range(periods):
        prod_amount = production_solution.get((p,t), 0)
        print(f"{prod_amount:>7.0f}", end="  ")
        total_production_cost += prod_amount * production_cost[product]
    print()

print()

# Display inventory plan
print("=== INVENTORY LEVELS ====")
print(f"{'Product':<8}", end="")
for q in range(periods):
    print(f"Q{q+1} Inv", end="  ")
print()
print("-" * 45)

total_inventory_cost = 0
for p, product in enumerate(products):
    print(f"{product:<8}", end="")
    for t in range(periods):
        inv_amount = inventory_solution.get((p,t), 0)
        print(f"{inv_amount:>7.0f}", end="  ")
        total_inventory_cost += inv_amount * inventory_cost[product]
    print()

print()
print(f" COST BREAKDOWN:")
print(f"   Production costs: ${total_production_cost:,.0f} ({100*total_production_cost/optimal_cost:.1f}%)")
print(f"   Inventory costs:  ${total_inventory_cost:,.0f} ({100*total_inventory_cost/optimal_cost:.1f}%)")
print(f"   Total cost:       ${optimal_cost:,.0f}")
print()

Out[18]:

🎯 MULTI-PERIOD OPTIMIZATION
================================
💰 OPTIMAL SOLUTION FOUND!
   Minimum total cost: $2,997,500

=== OPTIMAL PRODUCTION PLAN ====
Product Q1 Prod  Q2 Prod  Q3 Prod  Q4 Prod  
---------------------------------------------
CPU         900     1300     1400     1300  
GPU         750     1000     1100      900  
Memory     1800     2200     2500     2300  

=== INVENTORY LEVELS ====
Product Q1 Inv  Q2 Inv  Q3 Inv  Q4 Inv  
---------------------------------------------
CPU           0      100        0        0  
GPU           0        0        0        0  
Memory        0        0        0        0  

📊 COST BREAKDOWN:
   Production costs: $2,996,500 (100.0%)
   Inventory costs:  $1,000 (0.0%)
   Total cost:       $2,997,500

Practice Problems and Challenges

Test your understanding with these hands-on problems!

In [19]:

# Practice Problems
print(" PRACTICE PROBLEMS")
print("=" * 25)
print("Test your linear programming skills!")
print()

print(" PROBLEM SET:")
print()
print("1. DIET OPTIMIZATION:")
print("   A nutritionist wants to design a minimum-cost diet that meets")
print("   daily requirements for calories (≥2000), protein (≥50g), and")
print("   vitamins (≥100 units). Available foods:")
print("   - Bread: $2/loaf, 1000 cal, 10g protein, 20 vitamins")
print("   - Milk: $3/gallon, 500 cal, 25g protein, 50 vitamins")
print("   - Meat: $8/lb, 800 cal, 40g protein, 10 vitamins")
print()

# Solution template for Problem 1
print(" SOLUTION TEMPLATE:")
print("# Uncomment and complete:")
print("# diet = MixedIntegerLinearProgram(maximization=False)")
print("# food = diet.new_variable(real=True, nonnegative=True)")
print("# # Add objective and constraints...")
print()

print("2. BLENDING PROBLEM:")
print("   An oil refinery blends crude oils to produce gasoline.")
print("   Two crude types: Light ($60/barrel, 0.8 yield) and")
print("   Heavy ($50/barrel, 0.6 yield). Need 10,000 gallons gasoline.")
print("   Light crude limited to 8,000 barrels. Minimize cost.")
print()

print("3. WORKFORCE SCHEDULING:")
print("   A call center needs staff for different shifts.")
print("   Requirements: Morning 50, Afternoon 80, Evening 30")
print("   Workers can work: 8-hour shift $120, 4-hour shift $80")
print("   Minimize total cost while meeting requirements.")
print()

# Interactive challenge
print(" CHALLENGE YOURSELF:")
print("Choose one problem above and implement the complete solution!")
print("Include:")
print("- Problem formulation with variables and constraints")
print("- SageMath implementation")
print("- Solution interpretation and insights")
print("- Sensitivity analysis")
print()

Out[19]:

🎯 PRACTICE PROBLEMS
=========================
Test your linear programming skills!

📝 PROBLEM SET:

1. DIET OPTIMIZATION:
   A nutritionist wants to design a minimum-cost diet that meets
   daily requirements for calories (≥2000), protein (≥50g), and
   vitamins (≥100 units). Available foods:
   - Bread: $2/loaf, 1000 cal, 10g protein, 20 vitamins
   - Milk: $3/gallon, 500 cal, 25g protein, 50 vitamins
   - Meat: $8/lb, 800 cal, 40g protein, 10 vitamins

🔧 SOLUTION TEMPLATE:
# Uncomment and complete:
# diet = MixedIntegerLinearProgram(maximization=False)
# food = diet.new_variable(real=True, nonnegative=True)
# # Add objective and constraints...

2. BLENDING PROBLEM:
   An oil refinery blends crude oils to produce gasoline.
   Two crude types: Light ($60/barrel, 0.8 yield) and
   Heavy ($50/barrel, 0.6 yield). Need 10,000 gallons gasoline.
   Light crude limited to 8,000 barrels. Minimize cost.

3. WORKFORCE SCHEDULING:
   A call center needs staff for different shifts.
   Requirements: Morning 50, Afternoon 80, Evening 30
   Workers can work: 8-hour shift $120, 4-hour shift $80
   Minimize total cost while meeting requirements.

🏆 CHALLENGE YOURSELF:
Choose one problem above and implement the complete solution!
Include:
- Problem formulation with variables and constraints
- SageMath implementation
- Solution interpretation and insights
- Sensitivity analysis

📋 Summary and Next Steps

What You've Mastered

Congratulations! You've completed a comprehensive journey through linear programming with SageMath. You now understand:

Core Concepts:

Linear Programming Fundamentals: Standard forms, feasible regions, and optimality conditions
Geometric Interpretation: Visualizing constraints and understanding vertex optimality
Duality Theory: Primal-dual relationships and economic interpretation via shadow prices
Sensitivity Analysis: Understanding solution robustness and parameter changes

Problem-Solving Techniques:

Transportation Problems: Supply chain optimization and distribution planning
Integer Programming: Discrete optimization with binary and integer variables
Portfolio Optimization: Financial applications and risk management
Multi-Period Planning: Dynamic optimization with inventory and capacity constraints

SageMath Expertise:

MixedIntegerLinearProgram() for model setup and solving
Constraint formulation and variable declaration
Solution interpretation and optimality analysis
Integration with visualization and data analysis tools

Advanced Topics to Explore

1. Non-Linear Programming

Quadratic programming and convex optimization
Constrained optimization with Lagrange multipliers
Applications in machine learning and portfolio theory

2. Advanced Integer Programming

Branch-and-bound algorithms
Cutting plane methods
Large-scale combinatorial optimization

3. Stochastic Programming

Optimization under uncertainty
Scenario-based planning
Risk-averse decision making

4. Network Optimization

Shortest path and maximum flow problems
Graph theory applications
Telecommunications and logistics networks

5. Game Theory and Mechanism Design

Strategic decision making
Auction theory
Market design applications

Real-World Applications

Apply your linear programming skills to:

Business Operations: Production planning, inventory management, workforce scheduling
Finance: Portfolio optimization, risk management, capital allocation
Supply Chain: Transportation, distribution, vendor selection
Energy: Power generation, smart grid optimization, renewable energy planning
Healthcare: Resource allocation, scheduling, treatment planning
Government: Budget allocation, policy optimization, public service planning

Research and Career Paths

Operations Research: Mathematical optimization in business and industry
Data Science: Optimization components in machine learning and analytics
Financial Engineering: Quantitative methods in finance and risk management
Supply Chain Analytics: End-to-end optimization of global supply networks
Management Consulting: Strategic decision support using optimization

Continue Learning

Online Resources:

SageMath Optimization Documentation
NEOS Optimization Guide
Operations Research professional societies (INFORMS, EURO)

** Remember: Linear programming is not just a mathematical technique—it's a powerful way of thinking about resource allocation, trade-offs, and optimal decision making that applies across virtually every field of human endeavor.**

Continue exploring, keep optimizing, and use CoCalc's collaborative features to work with others on challenging optimization problems!