Capacity optimization of Supply Chain

import numpy as np
import scipy.optimize as opt
from scipy.stats import uniform
import matplotlib.pyplot as plt

# Define parameters
a = 100  # Upper bound of uniform demand distribution
p = 70  # Retail price
m = 5  # Unit cost
B = 10  # Scaling factor for cost function
pi_res = 20  # Supplier's reservation profit

# Define cost functions
c_linear = lambda K: B * K
c_quadratic = lambda K: B * K**2

# Demand distribution functions
F = lambda x: uniform.cdf(x, 0, a)
f = lambda x: uniform.pdf(x, 0, a)

# Expected demand function
E_D = a / 2
E_min = lambda K: (a - K) ** 2 / (2 * a) + K

# Centralized model optimization
def centralized_model(cost_function):
    obj_func = lambda K: (p - m) * E_min(K) - cost_function(K)
    res = opt.minimize_scalar(obj_func, bounds=(0, a), method='bounded')
    return res.x

K_c_linear = centralized_model(c_linear)
K_c_quadratic = centralized_model(c_quadratic)

# Wholesale price contract
def wholesale_price_contract(cost_function):
    obj_func = lambda K: (p - m) * E_min(K) - cost_function(K)
    res = opt.minimize_scalar(obj_func, bounds=(0, a), method='bounded')
    return res.x

K_w_linear = wholesale_price_contract(c_linear)
K_w_quadratic = wholesale_price_contract(c_quadratic)

# Linear contract
def linear_contract(K_c, cost_function):
    w = p
    t = - (p - m) * (E_D - np.trapz([(x - K_c) * f(x) for x in np.linspace(K_c, a, 100)], dx=(a-K_c)/100)) + cost_function(K_c) + pi_res
    return w, t

w_linear, t_linear = linear_contract(K_c_linear, c_linear)

# Cost-sharing contract
def cost_sharing_contract(K_c, cost_function):
    w = m + pi_res / ((p - m) * E_min(K_c) - cost_function(K_c))
    sigma = pi_res / ((p - m) * E_min(K_c) - cost_function(K_c))
    return w, sigma

w_cost_sharing, sigma_cost_sharing = cost_sharing_contract(K_c_linear, c_linear)

# Print results
print(f"Centralized Model: K_c_linear = {K_c_linear}, K_c_quadratic = {K_c_quadratic}")
print(f"Wholesale Price Contract: K_w_linear = {K_w_linear}, K_w_quadratic = {K_w_quadratic}")
print(f"Linear Contract: w = {w_linear}, t = {t_linear}")
print(f"Cost Sharing Contract: w = {w_cost_sharing}, sigma = {sigma_cost_sharing}")

# Plot efficiency comparison
labels = ['Centralized', 'Wholesale', 'Linear', 'Cost Sharing']
k_values = [K_c_linear, K_w_linear, K_c_linear, K_c_linear]

plt.figure(figsize=(8, 5))
plt.bar(labels, k_values, color=['blue', 'red', 'green', 'purple'])
plt.xlabel('Contract Type')
plt.ylabel('Optimal Capacity (K)')
plt.title('Optimal Capacity Under Different Contracts')
plt.show()