Linear Programming

import numpy as np

def simplex(c, A, b):
    """Simplex method for: min c^T x, s.t. Ax <= b, x >= 0."""
    m, n = A.shape
    # Add slack variables
    M = np.hstack([A, np.eye(m)])
    c_aug = np.concatenate([c, np.zeros(m)])
    
    # Initial BFS: slack variables are basic
    basis = list(range(n, n + m))
    non_basis = list(range(n))
    
    while True:
        # Compute reduced costs
        B = M[:, basis]
        N = M[:, non_basis]
        c_B = c_aug[basis]
        c_N = c_aug[non_basis]
        
        y = np.linalg.solve(B.T, c_B)
        reduced_costs = c_N - N.T @ y
        
        if np.all(reduced_costs >= -1e-10):
            break  # Optimal
        
        # Enter: most negative reduced cost
        enter_idx = np.argmin(reduced_costs)
        enter_var = non_basis[enter_idx]
        
        # Direction vector
        d = np.linalg.solve(B, M[:, enter_var])
        
        if np.all(d <= 1e-10):
            raise ValueError("LP is unbounded")
        
        # Leave: minimum ratio test
        x_B = np.linalg.solve(B, b)
        ratios = np.where(d > 1e-10, x_B / d, np.inf)
        leave_idx = np.argmin(ratios)
        leave_var = basis[leave_idx]
        
        # Update basis
        basis[leave_idx] = enter_var
        non_basis[enter_idx] = leave_var
    
    x = np.zeros(n + m)
    x_B = np.linalg.solve(M[:, basis], b)
    for i, v in enumerate(basis):
        x[v] = x_B[i]
    
    return x[:n], c_aug @ x

# Example
c = np.array([-1, -2])
A = np.array([[1, 1], [1, -1], [-1, 1]])
b = np.array([4, 2, 2])
x_opt, obj = simplex(c, A, b)
print(f"Optimal: x={x_opt}, obj={-obj}")

from scipy.optimize import linprog
import numpy as np

def transportation_problem(costs, supply, demand):
    """Solve transportation problem using LP."""
    m, n = costs.shape
    # Decision variables: x_{ij} flattened as x[i*n + j]
    c = costs.flatten()
    
    # Supply constraints: sum over j for each i
    A_eq_supply = np.zeros((m, m * n))
    for i in range(m):
        A_eq_supply[i, i*n:(i+1)*n] = 1
    
    # Demand constraints: sum over i for each j
    A_eq_demand = np.zeros((n, m * n))
    for j in range(n):
        A_eq_demand[j, j::n] = 1
    
    A_eq = np.vstack([A_eq_supply, A_eq_demand])
    b_eq = np.concatenate([supply, demand])
    
    result = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=[(0, None)] * (m*n), method='highs')
    
    if result.success:
        return result.x.reshape(m, n), result.fun
    return None, None

costs = np.array([[8, 6, 10, 9], [9, 12, 13, 7], [14, 9, 16, 5]])
supply = np.array([35, 50, 40])
demand = np.array([45, 20, 30, 30])

shipment, total_cost = transportation_problem(costs, supply, demand)
print(f"Optimal shipments:\n{shipment}")
print(f"Total cost: {total_cost}")

from scipy.optimize import linprog
import numpy as np

# Example: Maximize z = 3x + 2y
# Subject to: x + y <= 4, x <= 2, y <= 3, x,y >= 0
c = [-3, -2]  # Negate for maximization
A_ub = [[1, 1], [1, 0], [0, 1]]
b_ub = [4, 2, 3]

result = linprog(c, A_ub=A_ub, b_ub=b_ub, method='highs')
print(f"Optimal x: {result.x}")
print(f"Maximum z: {-result.fun}")

# Equality constraints example
# min 5x + 4y, s.t. 6x + 4y = 24, x + 2y = 6, x,y >= 0
c2 = [5, 4]
A_eq = [[6, 4], [1, 2]]
b_eq = [24, 6]

result2 = linprog(c2, A_eq=A_eq, b_eq=b_eq, bounds=[(0, None), (0, None)])
print(f"Optimal: x={result2.x}, cost={result2.fun}")

from pulp import *

# Define the problem
prob = LpProblem("Production_Planning", LpMaximize)

# Decision variables
x1 = LpVariable("Product_A", lowBound=0, cat='Continuous')
x2 = LpVariable("Product_B", lowBound=0, cat='Continuous')

# Objective function
prob += 3 * x1 + 2 * x2, "Total_Profit"

# Constraints
prob += x1 + x2 <= 4, "Resource_1"
prob += x1 <= 2, "Resource_2"
prob += x2 <= 3, "Resource_3"

# Solve
prob.solve(PULP_CBC_CMD(msg=0))

print(f"Status: {LpStatus[prob.status]}")
print(f"Product A = {x1.varValue}")
print(f"Product B = {x2.varValue}")
print(f"Total Profit = {value(prob.objective)}")

# Mixed-Integer example
prob2 = LpProblem("ILP_Example", LpMaximize)
x = LpVariable("x", lowBound=0, cat='Integer')
y = LpVariable("y", lowBound=0, cat='Continuous')

prob2 += 5 * x + 4 * y
prob2 += 2 * x + y <= 8
prob2 += x + 3 * y <= 9
prob2 += x <= 3

prob2.solve(PULP_CBC_CMD(msg=0))
print(f"ILP: x={x.varValue}, y={y.varValue}, obj={value(prob2.objective)}")

import numpy as np
from scipy.optimize import linprog

# Lasso regression as LP: min ||y - Xw||_1
# Reformulated: min sum(t_i), s.t. y - Xw <= t, -(y - Xw) <= t
np.random.seed(42)
X = np.random.randn(50, 5)
true_w = np.array([1, 0, 0, 0, 0], dtype=float)
y = X @ true_w + 0.1 * np.random.randn(50)

n, p = X.shape
# Variables: [w (p), t (n)]
c = np.concatenate([np.zeros(p), np.ones(n)])

# Constraints: y - Xw <= t  =>  -Xw - t <= -y
#             -(y - Xw) <= t  =>  Xw - t <= y
A_ub = np.vstack([
    np.hstack([-X, -np.eye(n)]),
    np.hstack([X, -np.eye(n)])
])
b_ub = np.concatenate([-y, y])

result = linprog(c, A_ub=A_ub, b_ub=b_ub, method='highs')
w_lasso = result.x[:p]
print(f"Lasso coefficients: {w_lasso.round(3)}")

from scipy.optimize import linprog

# Linear risk model: minimize tracking error
# min sum |x_i - w_i| subject to budget and sector constraints
n_assets = 5
target = np.array([0.2, 0.2, 0.2, 0.2, 0.2])

# Variables: [x (n), s_plus (n), s_minus (n)]
c = np.concatenate([np.zeros(n_assets), np.ones(n_assets), np.ones(n_assets)])

# x - target = s+ - s-
# x - s+ + s- = target  (equality)
A_eq = np.hstack([
    np.eye(n_assets),
    -np.eye(n_assets),
    np.eye(n_assets)
])
b_eq = target

# Budget: sum(x) = 1
A_budget = np.hstack([np.ones((1, n_assets)), np.zeros((1, 2 * n_assets))])
b_budget = [1.0]

A_eq_full = np.vstack([A_eq, A_budget])
b_eq_full = np.concatenate([b_eq, b_budget])

result = linprog(c, A_eq=A_eq_full, b_eq=b_eq_full,
                 bounds=[(0, 1)] * n_assets + [(0, None)] * (2 * n_assets),
                 method='highs')
print(f"Optimal weights: {result.x[:n_assets].round(3)}")