import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import minimize

# Define the second-order system
def system(state, t, Kp, Ki, Kd, setpoint):
    x, x_dot, integral = state
    error = setpoint - x
    u = Kp * error + Ki * integral + Kd * (-x_dot)
    m, c, k = 1.0, 0.2, 1.0  # Mass, damping, spring constant
    x_ddot = (u - c * x_dot - k * x) / m
    return [x_dot, x_ddot, error]

# Simulate the system for given PID parameters
def simulate_pid(params, t, setpoint=1.0):
    Kp, Ki, Kd = params
    initial_state = [0.0, 0.0, 0.0]  # Initial position, velocity, integral
    states = odeint(system, initial_state, t, args=(Kp, Ki, Kd, setpoint))
    return states[:, 0], states[:, 2]  # Return position and integral of error

# Objective function for optimization (minimize ITAE)
def objective(params, t, setpoint=1.0):
    y, integral = simulate_pid(params, t, setpoint)
    error = setpoint - y
    itae = np.trapz(t * np.abs(error), t)
    return itae

# Optimization and simulation
t = np.linspace(0, 10, 1000)
initial_params = [1.0, 0.1, 0.1]  # Initial guess for Kp, Ki, Kd
bounds = [(0, 10), (0, 10), (0, 10)]  # Bounds for Kp, Ki, Kd
result = minimize(objective, initial_params, args=(t,), bounds=bounds)

# Get optimized parameters
Kp_opt, Ki_opt, Kd_opt = result.x
y_opt, _ = simulate_pid(result.x, t)

# Simulate with initial parameters for comparison
y_init, _ = simulate_pid(initial_params, t)

# Plot results
plt.figure(figsize=(10, 6))
plt.plot(t, y_opt, label=f'Optimized: Kp={Kp_opt:.2f}, Ki={Ki_opt:.2f}, Kd={Kd_opt:.2f}')
plt.plot(t, y_init, label=f'Initial: Kp={initial_params[0]}, Ki={initial_params[1]}, Kd={initial_params[2]}')
plt.plot(t, [1.0]*len(t), 'k--', label='Setpoint')
plt.xlabel('Time (s)')
plt.ylabel('Output')
plt.title('PID Controller Response: Initial vs Optimized')
plt.legend()
plt.grid(True)
plt.show()

print(f"Optimized Parameters: Kp={Kp_opt:.2f}, Ki={Ki_opt:.2f}, Kd={Kd_opt:.2f}")
print(f"ITAE (Optimized): {objective(result.x, t):.2f}")
print(f"ITAE (Initial): {objective(initial_params, t):.2f}")

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar

# Parameters
F = 1.0        # Bioavailability
D = 500        # Dose in mg
ka = 1.5       # Absorption rate constant (1/hr)
ke = 0.3       # Elimination rate constant (1/hr)
Vd = 40        # Volume of distribution (L)
t_end = 72     # Total simulation time (hr)
dose_count = 10

# Target therapeutic range
C_min = 10     # Minimum effective concentration (mg/L)
C_max = 20     # Maximum safe concentration (mg/L)

def concentration_single_dose(t):
    """Concentration after a single oral dose at time t"""
    return (F * D * ka) / (Vd * (ka - ke)) * (np.exp(-ke * t) - np.exp(-ka * t))

def concentration_multiple_doses(tau):
    """Simulate multiple doses given at interval tau"""
    t = np.linspace(0, t_end, 1000)
    C = np.zeros_like(t)
    for n in range(dose_count):
        tn = n * tau
        C += concentration_single_dose(np.clip(t - tn, 0, None))
    return t, C

def objective(tau):
    """Objective: Penalize time outside therapeutic range"""
    t, C = concentration_multiple_doses(tau)
    penalty = np.sum((C < C_min) * (C_min - C)**2 + (C > C_max) * (C - C_max)**2)
    return penalty

# Optimization
result = minimize_scalar(objective, bounds=(4, 24), method='bounded')
optimal_tau = result.x
print(f"Optimal dosing interval: {optimal_tau:.2f} hours")

# Final simulation with optimal tau
t, C = concentration_multiple_doses(optimal_tau)

plt.figure(figsize=(12, 5))
plt.plot(t, C, label='Drug Concentration')
plt.axhline(C_min, color='green', linestyle='--', label='Min Effective (10 mg/L)')
plt.axhline(C_max, color='red', linestyle='--', label='Max Safe (20 mg/L)')
plt.title(f"Optimized Drug Dosing (Interval: {optimal_tau:.2f} hr)")
plt.xlabel('Time (hours)')
plt.ylabel('Concentration (mg/L)')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Set random seed for reproducibility
np.random.seed(42)

# Depot location
depot = np.array([50, 50])

# Generate random delivery locations
n_deliveries = 10
locations = np.random.rand(n_deliveries, 2) * 100  # 100x100 area

# Calculate energy usage for each round trip (Euclidean distance * 2)
def calculate_total_energy(depot, locations):
    distances = np.linalg.norm(locations - depot, axis=1)
    return np.sum(distances * 2)

total_energy = calculate_total_energy(depot, locations)

print(f"Total round-trip energy usage (arbitrary units): {total_energy:.2f}")

# Plot depot and delivery locations
plt.figure(figsize=(8, 8))
plt.scatter(*depot, color='red', label='Depot', s=100)
plt.scatter(locations[:, 0], locations[:, 1], label='Deliveries', s=60)

# Draw lines for round trips
for loc in locations:
    plt.plot([depot[0], loc[0]], [depot[1], loc[1]], color='gray', linestyle='--', alpha=0.6)

plt.title('Electric Delivery Routes')
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.grid(True)
plt.legend()
plt.show()

from sklearn.cluster import KMeans

# Assume we use 3 trucks and cluster the deliveries
kmeans = KMeans(n_clusters=3, random_state=42)
kmeans.fit(locations)
labels = kmeans.labels_
centroids = kmeans.cluster_centers_

# Recalculate energy: depot to centroid and back (simulating drop-off and collection hub)
cluster_energy = np.sum([np.linalg.norm(c - depot) * 2 for c in centroids])

print(f"Energy usage with clustered deliveries: {cluster_energy:.2f}")

plt.bar(['Individual Trips', 'Clustered'], [total_energy, cluster_energy], color=['blue', 'green'])
plt.title('Energy Usage Comparison')
plt.ylabel('Energy Units')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import linprog
import pandas as pd
from matplotlib.patches import Rectangle
import warnings
warnings.filterwarnings('ignore')

# Set random seed for reproducibility
np.random.seed(42)

# Problem Setup: 10 potential locations for renewable energy installation
n_locations = 10
location_names = [f'Site_{i+1}' for i in range(n_locations)]

# Generate realistic data for each location
# Wind speed data (m/s) - affects wind turbine efficiency
wind_speeds = np.random.uniform(4.5, 12.0, n_locations)

# Solar irradiance data (kWh/m²/day) - affects solar panel efficiency  
solar_irradiance = np.random.uniform(3.5, 6.5, n_locations)

# Wind turbine specifications
turbine_capacity = 2.5  # MW per turbine
turbine_cost = 3.0      # Million $ per turbine
turbine_maintenance = 0.1  # Million $ per year per turbine

# Solar panel specifications  
solar_capacity_per_unit = 1.0  # MW per solar unit
solar_cost = 1.5              # Million $ per solar unit
solar_maintenance = 0.05      # Million $ per year per solar unit

# Grid connection costs (varies by location)
grid_connection_cost = np.random.uniform(0.2, 0.8, n_locations)

# Calculate energy output coefficients
# Wind energy output: P = 0.5 * ρ * A * v³ * Cp (simplified)
# We'll use a simplified efficiency curve
def wind_efficiency(wind_speed):
    """Calculate wind turbine efficiency based on wind speed"""
    if wind_speed < 3:
        return 0
    elif wind_speed < 12:
        return 0.35 * (wind_speed / 12) ** 3
    else:
        return 0.35

# Solar efficiency is more linear with irradiance
def solar_efficiency(irradiance):
    """Calculate solar panel efficiency based on irradiance"""
    return 0.18 * (irradiance / 5.0)

# Calculate annual energy output for each location (GWh/year)
wind_output = np.array([wind_efficiency(ws) * turbine_capacity * 8760 / 1000 
                       for ws in wind_speeds])
solar_output = np.array([solar_efficiency(si) * solar_capacity_per_unit * 8760 / 1000 
                        for si in solar_irradiance])

# Total costs per installation (including 5-year maintenance)
wind_total_cost = turbine_cost + 5 * turbine_maintenance + grid_connection_cost
solar_total_cost = solar_cost + 5 * solar_maintenance + grid_connection_cost

# Budget constraint
total_budget = 50.0  # Million $

print("=== Renewable Energy Optimization Problem ===")
print(f"Budget: ${total_budget} million")
print(f"Number of potential sites: {n_locations}")
print(f"Wind turbine capacity: {turbine_capacity} MW")
print(f"Solar unit capacity: {solar_capacity_per_unit} MW")
print("\n=== Location Data ===")

# Create comprehensive data table
data_table = pd.DataFrame({
    'Location': location_names,
    'Wind_Speed_ms': np.round(wind_speeds, 2),
    'Solar_Irradiance_kWh_m2_day': np.round(solar_irradiance, 2),
    'Wind_Output_GWh_year': np.round(wind_output, 2),
    'Solar_Output_GWh_year': np.round(solar_output, 2),
    'Wind_Cost_Million': np.round(wind_total_cost, 2),
    'Solar_Cost_Million': np.round(solar_total_cost, 2),
    'Grid_Connection_Million': np.round(grid_connection_cost, 2)
})

print(data_table.to_string(index=False))

# Mathematical Formulation
print("\n=== Mathematical Formulation ===")
print("Objective Function:")
print("Maximize: Σ(wind_output[i] * x_wind[i] + solar_output[i] * x_solar[i])")
print("\nSubject to:")
print("Σ(wind_cost[i] * x_wind[i] + solar_cost[i] * x_solar[i]) ≤ Budget")
print("x_wind[i], x_solar[i] ∈ {0, 1} for all i")

# Since scipy.optimize.linprog doesn't handle integer programming directly,
# we'll use a relaxed LP and then apply rounding heuristics
# For exact integer solutions, specialized solvers like PuLP would be better

# Prepare linear programming formulation
# Variables: [x_wind_0, x_wind_1, ..., x_wind_9, x_solar_0, x_solar_1, ..., x_solar_9]
c = np.concatenate([-wind_output, -solar_output])  # Negative for maximization
A_ub = np.concatenate([wind_total_cost, solar_total_cost]).reshape(1, -1)
b_ub = np.array([total_budget])

# Bounds: 0 ≤ x ≤ 1 for relaxed problem
bounds = [(0, 1) for _ in range(2 * n_locations)]

# Solve relaxed linear program
print("\n=== Solving Optimization Problem ===")
result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')

if result.success:
    print("Optimization successful!")
    
    # Extract solutions
    x_wind_relaxed = result.x[:n_locations]
    x_solar_relaxed = result.x[n_locations:]
    
    # Apply rounding heuristic for integer solution
    # Sort by efficiency ratio (output/cost) and select greedily
    wind_efficiency_ratio = wind_output / wind_total_cost
    solar_efficiency_ratio = solar_output / solar_total_cost
    
    # Create candidate list with efficiency ratios
    candidates = []
    for i in range(n_locations):
        candidates.append(('wind', i, wind_efficiency_ratio[i], wind_total_cost[i], wind_output[i]))
        candidates.append(('solar', i, solar_efficiency_ratio[i], solar_total_cost[i], solar_output[i]))
    
    # Sort by efficiency ratio (descending)
    candidates.sort(key=lambda x: x[2], reverse=True)
    
    # Greedy selection
    selected_wind = np.zeros(n_locations, dtype=int)
    selected_solar = np.zeros(n_locations, dtype=int)
    remaining_budget = total_budget
    total_output = 0
    
    print("\n=== Greedy Selection Process ===")
    for energy_type, location, ratio, cost, output in candidates:
        if cost <= remaining_budget:
            if energy_type == 'wind':
                selected_wind[location] = 1
            else:
                selected_solar[location] = 1
            remaining_budget -= cost
            total_output += output
            print(f"Selected: {energy_type.upper()} at {location_names[location]} "
                  f"(Efficiency: {ratio:.3f}, Cost: ${cost:.2f}M, Output: {output:.2f} GWh/year)")
    
    print(f"\nRemaining budget: ${remaining_budget:.2f} million")
    print(f"Total annual energy output: {total_output:.2f} GWh/year")
    
    # Calculate total investment
    total_investment = total_budget - remaining_budget
    print(f"Total investment: ${total_investment:.2f} million")
    
else:
    print("Optimization failed!")
    selected_wind = np.zeros(n_locations, dtype=int)
    selected_solar = np.zeros(n_locations, dtype=int)

# Create comprehensive results analysis
print("\n=== Final Solution Analysis ===")
solution_df = pd.DataFrame({
    'Location': location_names,
    'Wind_Selected': selected_wind,
    'Solar_Selected': selected_solar,
    'Wind_Speed': np.round(wind_speeds, 2),
    'Solar_Irradiance': np.round(solar_irradiance, 2),
    'Total_Output_GWh': np.round(
        selected_wind * wind_output + selected_solar * solar_output, 2),
    'Total_Cost_Million': np.round(
        selected_wind * wind_total_cost + selected_solar * solar_total_cost, 2)
})

print(solution_df.to_string(index=False))

# Visualization Section
plt.style.use('seaborn-v0_8')
fig = plt.figure(figsize=(20, 15))

# 1. Location Overview Map-style Plot
ax1 = plt.subplot(2, 3, 1)
scatter_wind = ax1.scatter([i for i in range(n_locations) if selected_wind[i]], 
                         [wind_speeds[i] for i in range(n_locations) if selected_wind[i]],
                         s=[wind_output[i]*50 for i in range(n_locations) if selected_wind[i]], 
                         c='blue', alpha=0.7, label='Wind Turbines', marker='^')

scatter_solar = ax1.scatter([i for i in range(n_locations) if selected_solar[i]], 
                          [solar_irradiance[i] for i in range(n_locations) if selected_solar[i]],
                          s=[solar_output[i]*50 for i in range(n_locations) if selected_solar[i]], 
                          c='orange', alpha=0.7, label='Solar Panels', marker='s')

# Add unselected locations in gray
unselected_idx = [i for i in range(n_locations) 
                 if not selected_wind[i] and not selected_solar[i]]
if unselected_idx:
    ax1.scatter(unselected_idx, [wind_speeds[i] for i in unselected_idx],
               c='gray', alpha=0.3, s=30, label='Unselected')

ax1.set_xlabel('Location Index')
ax1.set_ylabel('Resource Quality')
ax1.set_title('Selected Renewable Energy Sites\n(Bubble size = Annual Output)')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 2. Cost vs Output Analysis
ax2 = plt.subplot(2, 3, 2)
wind_costs = [wind_total_cost[i] if selected_wind[i] else 0 for i in range(n_locations)]
solar_costs = [solar_total_cost[i] if selected_solar[i] else 0 for i in range(n_locations)]
wind_outputs = [wind_output[i] if selected_wind[i] else 0 for i in range(n_locations)]
solar_outputs = [solar_output[i] if selected_solar[i] else 0 for i in range(n_locations)]

x_pos = np.arange(n_locations)
width = 0.35

bars1 = ax2.bar(x_pos - width/2, wind_costs, width, label='Wind Cost', color='lightblue', alpha=0.7)
bars2 = ax2.bar(x_pos + width/2, solar_costs, width, label='Solar Cost', color='lightyellow', alpha=0.7)

ax2_twin = ax2.twinx()
line1 = ax2_twin.plot(x_pos, wind_outputs, 'bo-', label='Wind Output', linewidth=2, markersize=6)
line2 = ax2_twin.plot(x_pos, solar_outputs, 'ro-', label='Solar Output', linewidth=2, markersize=6)

ax2.set_xlabel('Location')
ax2.set_ylabel('Cost (Million $)', color='black')
ax2_twin.set_ylabel('Output (GWh/year)', color='black')
ax2.set_title('Cost vs Output by Location')
ax2.set_xticks(x_pos)
ax2.set_xticklabels([f'S{i+1}' for i in range(n_locations)], rotation=45)

# Combine legends
lines1, labels1 = ax2.get_legend_handles_labels()
lines2, labels2 = ax2_twin.get_legend_handles_labels()
ax2.legend(lines1 + lines2, labels1 + labels2, loc='upper left')
ax2.grid(True, alpha=0.3)

# 3. Resource Quality Distribution
ax3 = plt.subplot(2, 3, 3)
selected_locations = [i for i in range(n_locations) 
                     if selected_wind[i] or selected_solar[i]]
unselected_locations = [i for i in range(n_locations) 
                       if not (selected_wind[i] or selected_solar[i])]

if selected_locations:
    ax3.scatter([wind_speeds[i] for i in selected_locations], 
               [solar_irradiance[i] for i in selected_locations],
               c='green', s=100, alpha=0.7, label='Selected Sites')

if unselected_locations:
    ax3.scatter([wind_speeds[i] for i in unselected_locations], 
               [solar_irradiance[i] for i in unselected_locations],
               c='red', s=100, alpha=0.7, label='Unselected Sites')

ax3.set_xlabel('Wind Speed (m/s)')
ax3.set_ylabel('Solar Irradiance (kWh/m²/day)')
ax3.set_title('Resource Quality: Selected vs Unselected')
ax3.legend()
ax3.grid(True, alpha=0.3)

# 4. Budget Utilization
ax4 = plt.subplot(2, 3, 4)
budget_categories = ['Wind\nInvestment', 'Solar\nInvestment', 'Remaining\nBudget']
budget_values = [
    np.sum(selected_wind * wind_total_cost),
    np.sum(selected_solar * solar_total_cost),
    remaining_budget
]
colors = ['skyblue', 'gold', 'lightgray']

wedges, texts, autotexts = ax4.pie(budget_values, labels=budget_categories, 
                                  autopct='%1.1f%%', startangle=90, colors=colors)
ax4.set_title(f'Budget Utilization\n(Total: ${total_budget}M)')

# 5. Efficiency Comparison
ax5 = plt.subplot(2, 3, 5)
efficiency_data = {
    'Technology': ['Wind Turbines', 'Solar Panels'],
    'Avg_Efficiency_Ratio': [
        np.mean([wind_efficiency_ratio[i] for i in range(n_locations) if selected_wind[i]]) 
        if np.sum(selected_wind) > 0 else 0,
        np.mean([solar_efficiency_ratio[i] for i in range(n_locations) if selected_solar[i]]) 
        if np.sum(selected_solar) > 0 else 0
    ],
    'Count': [np.sum(selected_wind), np.sum(selected_solar)]
}

bars = ax5.bar(efficiency_data['Technology'], efficiency_data['Avg_Efficiency_Ratio'], 
               color=['lightblue', 'lightyellow'], alpha=0.8)

# Add count labels on bars
for bar, count in zip(bars, efficiency_data['Count']):
    height = bar.get_height()
    ax5.text(bar.get_x() + bar.get_width()/2., height + 0.01,
             f'Count: {count}', ha='center', va='bottom')

ax5.set_ylabel('Average Efficiency Ratio (GWh/$M)')
ax5.set_title('Technology Efficiency Comparison')
ax5.grid(True, alpha=0.3)

# 6. Timeline Projection
ax6 = plt.subplot(2, 3, 6)
years = np.arange(1, 11)  # 10-year projection
annual_output = total_output
cumulative_output = annual_output * years
cumulative_investment = total_investment + np.arange(1, 11) * 0.1 * total_investment  # 10% annual maintenance

ax6.plot(years, cumulative_output, 'g-', linewidth=3, label='Cumulative Energy Output (GWh)')
ax6_twin = ax6.twinx()
ax6_twin.plot(years, cumulative_investment, 'r--', linewidth=2, label='Cumulative Investment ($M)')

ax6.set_xlabel('Years')
ax6.set_ylabel('Cumulative Energy (GWh)', color='green')
ax6_twin.set_ylabel('Cumulative Cost ($M)', color='red')
ax6.set_title('10-Year Energy & Investment Projection')
ax6.grid(True, alpha=0.3)

# Combine legends
lines1, labels1 = ax6.get_legend_handles_labels()
lines2, labels2 = ax6_twin.get_legend_handles_labels()
ax6.legend(lines1 + lines2, labels1 + labels2, loc='upper left')

plt.tight_layout()
plt.show()

# Mathematical Summary
print("\n" + "="*60)
print("MATHEMATICAL FORMULATION SUMMARY")
print("="*60)

print("\nObjective Function (Maximization):")
print("$$\\max \\sum_{i=1}^{n} \\left( E_{wind}^i \\cdot x_{wind}^i + E_{solar}^i \\cdot x_{solar}^i \\right)$$")

print("\nWhere:")
print("- $E_{wind}^i$: Annual wind energy output at location $i$")
print("- $E_{solar}^i$: Annual solar energy output at location $i$")  
print("- $x_{wind}^i$: Binary decision variable for wind installation at location $i$")
print("- $x_{solar}^i$: Binary decision variable for solar installation at location $i$")

print("\nConstraints:")
print("Budget Constraint:")
print("$$\\sum_{i=1}^{n} \\left( C_{wind}^i \\cdot x_{wind}^i + C_{solar}^i \\cdot x_{solar}^i \\right) \\leq B$$")

print("\nBinary Constraints:")
print("$$x_{wind}^i, x_{solar}^i \\in \\{0,1\\} \\quad \\forall i \\in \\{1,2,...,n\\}$$")

print(f"\nProblem Instance:")
print(f"- $n = {n_locations}$ (number of locations)")
print(f"- $B = {total_budget}$ million USD (budget)")
print(f"- Total variables: {2*n_locations} (binary)")

print(f"\nSolution Quality:")
print(f"- Total annual output: {total_output:.2f} GWh/year")
print(f"- Budget utilization: {(total_investment/total_budget)*100:.1f}%")
print(f"- Wind installations: {np.sum(selected_wind)}")
print(f"- Solar installations: {np.sum(selected_solar)}")
print(f"- Average ROI: {total_output/total_investment:.2f} GWh per million USD")

print("\n" + "="*60)
print("OPTIMIZATION COMPLETE")
print("="*60)

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# Given constants
P = 1000          # Load in N
L = 2.0           # Length in meters
rho = 7850        # Density of steel (kg/m^3)
g = 9.81          # Gravity (m/s^2)
sigma_yield = 250e6  # Yield stress in Pa
b = 0.05          # Fixed width in meters

# Objective function: weight of the beam
def weight(x):
    h = x[0]
    volume = b * h * L
    return rho * g * volume  # Weight = mass * gravity

# Constraint: maximum bending stress must be less than yield
def stress_constraint(x):
    h = x[0]
    sigma_max = 6 * P * L / (b * h**2)
    return sigma_yield - sigma_max

# Bounds and constraints
bounds = [(0.01, 0.5)]  # height h should be between 1cm and 50cm
constraints = [{'type': 'ineq', 'fun': stress_constraint}]

# Initial guess
x0 = [0.1]

# Run optimization
result = minimize(weight, x0, method='SLSQP', bounds=bounds, constraints=constraints)

# Extract result
h_opt = result.x[0]
w_opt = weight([h_opt])

print(f"Optimal height h = {h_opt:.4f} m")
print(f"Minimum weight = {w_opt:.2f} N")

def weight(x):
    h = x[0]
    volume = b * h * L
    return rho * g * volume

def stress_constraint(x):
    h = x[0]
    sigma_max = 6 * P * L / (b * h**2)
    return sigma_yield - sigma_max

h_vals = np.linspace(0.01, 0.5, 100)
weights = rho * g * b * h_vals * L
stresses = 6 * P * L / (b * h_vals**2)

fig, ax1 = plt.subplots()

ax1.set_xlabel('Height h (m)')
ax1.set_ylabel('Weight (N)', color='tab:blue')
ax1.plot(h_vals, weights, label='Weight', color='tab:blue')
ax1.tick_params(axis='y', labelcolor='tab:blue')

ax2 = ax1.twinx()
ax2.set_ylabel('Max Bending Stress (Pa)', color='tab:red')
ax2.plot(h_vals, stresses, label='Stress', color='tab:red')
ax2.axhline(sigma_yield, color='gray', linestyle='--', label='Yield Stress')
ax2.tick_params(axis='y', labelcolor='tab:red')

plt.title("Beam Weight and Stress vs Height")
fig.tight_layout()
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import random
from collections import defaultdict
import seaborn as sns
from scipy.optimize import differential_evolution
import warnings
warnings.filterwarnings('ignore')

# Set random seed for reproducibility
np.random.seed(42)
random.seed(42)

class HPProteinFolder:
    """
    HP Model Protein Folding Optimizer
    
    This class implements a simplified protein folding model where:
    - H represents hydrophobic amino acids
    - P represents polar amino acids
    - Goal: Maximize hydrophobic contacts while avoiding collisions
    """
    
    def __init__(self, sequence):
        """
        Initialize the protein folder with an HP sequence
        
        Args:
            sequence (str): String of 'H' and 'P' characters representing the protein
        """
        self.sequence = sequence
        self.length = len(sequence)
        self.moves = [(0, 1), (1, 0), (0, -1), (-1, 0)]  # Up, Right, Down, Left
        self.move_names = ['U', 'R', 'D', 'L']
        
    def decode_folding(self, angles):
        """
        Convert angle representation to 2D coordinates
        
        Args:
            angles (list): List of angles (0-3) representing folding directions
            
        Returns:
            tuple: (coordinates, valid_folding_flag)
        """
        if len(angles) != self.length - 1:
            return None, False
            
        coordinates = [(0, 0)]  # Start at origin
        current_pos = (0, 0)
        occupied = {(0, 0)}
        
        for i, angle in enumerate(angles):
            # Convert continuous angle to discrete move
            move_idx = int(angle) % 4
            dx, dy = self.moves[move_idx]
            new_pos = (current_pos[0] + dx, current_pos[1] + dy)
            
            # Check for collision (self-intersection)
            if new_pos in occupied:
                return None, False
                
            coordinates.append(new_pos)
            occupied.add(new_pos)
            current_pos = new_pos
            
        return coordinates, True
    
    def calculate_energy(self, coordinates):
        """
        Calculate the energy of a folding configuration
        
        Energy = -1 * (number of H-H contacts that are adjacent in space but not in sequence)
        
        Args:
            coordinates (list): List of (x, y) coordinates for each amino acid
            
        Returns:
            float: Energy value (lower is better)
        """
        if coordinates is None:
            return float('inf')  # Invalid folding has infinite energy
            
        energy = 0
        
        # Check all pairs of amino acids
        for i in range(self.length):
            for j in range(i + 2, self.length):  # Skip adjacent in sequence
                # Only count H-H interactions
                if self.sequence[i] == 'H' and self.sequence[j] == 'H':
                    # Check if they are adjacent in 2D space
                    dist = abs(coordinates[i][0] - coordinates[j][0]) + abs(coordinates[i][1] - coordinates[j][1])
                    if dist == 1:  # Manhattan distance = 1 means adjacent
                        energy -= 1  # Favorable interaction
                        
        return energy
    
    def objective_function(self, angles):
        """
        Objective function for optimization
        
        Args:
            angles (array): Array of folding angles
            
        Returns:
            float: Energy to minimize
        """
        coordinates, valid = self.decode_folding(angles)
        if not valid:
            return 1000  # Heavy penalty for invalid folding
            
        return self.calculate_energy(coordinates)
    
    def optimize_folding(self, method='differential_evolution', max_iterations=1000):
        """
        Optimize protein folding using specified optimization method
        
        Args:
            method (str): Optimization method to use
            max_iterations (int): Maximum number of iterations
            
        Returns:
            tuple: (best_angles, best_energy, optimization_history)
        """
        bounds = [(0, 3.99) for _ in range(self.length - 1)]
        history = []
        
        def callback(xk, convergence=None):
            energy = self.objective_function(xk)
            history.append(energy)
            return False
        
        if method == 'differential_evolution':
            result = differential_evolution(
                self.objective_function,
                bounds,
                maxiter=max_iterations // 10,
                popsize=15,
                callback=callback,
                seed=42
            )
            best_angles = result.x
            best_energy = result.fun
            
        return best_angles, best_energy, history
    
    def monte_carlo_optimization(self, max_iterations=10000, temperature=1.0, cooling_rate=0.995):
        """
        Monte Carlo optimization with simulated annealing
        
        Args:
            max_iterations (int): Maximum number of iterations
            temperature (float): Initial temperature
            cooling_rate (float): Temperature cooling rate
            
        Returns:
            tuple: (best_angles, best_energy, history)
        """
        current_angles = [random.uniform(0, 3.99) for _ in range(self.length - 1)]
        current_energy = self.objective_function(current_angles)
        
        best_angles = current_angles.copy()
        best_energy = current_energy
        history = [current_energy]
        
        temp = temperature
        
        for iteration in range(max_iterations):
            # Generate neighbor solution
            new_angles = current_angles.copy()
            idx = random.randint(0, len(new_angles) - 1)
            new_angles[idx] = random.uniform(0, 3.99)
            
            new_energy = self.objective_function(new_angles)
            
            # Accept or reject based on Metropolis criterion
            if new_energy < current_energy or random.random() < np.exp(-(new_energy - current_energy) / temp):
                current_angles = new_angles
                current_energy = new_energy
                
                if current_energy < best_energy:
                    best_angles = current_angles.copy()
                    best_energy = current_energy
            
            history.append(best_energy)
            temp *= cooling_rate
            
            if iteration % 1000 == 0:
                print(f"Iteration {iteration}: Best Energy = {best_energy:.3f}, Temp = {temp:.4f}")
        
        return best_angles, best_energy, history

# Example protein sequences for testing
sequences = {
    'simple': 'HPHPPHHPHH',  # 10 amino acids
    'medium': 'HPHPPHHPHHPPHPPHHPHH',  # 20 amino acids  
    'complex': 'HHPPHPHPHPPHPPHPPHPHPHHPPHPPH'  # 27 amino acids
}

print("=== Protein Folding Optimization Results ===\n")

# Let's work with the medium complexity sequence
sequence = sequences['medium']
print(f"Protein sequence: {sequence}")
print(f"Length: {len(sequence)} amino acids")
print(f"Hydrophobic residues: {sequence.count('H')}")
print(f"Polar residues: {sequence.count('P')}\n")

# Initialize the protein folder
folder = HPProteinFolder(sequence)

# Optimize using differential evolution
print("Running Differential Evolution optimization...")
de_angles, de_energy, de_history = folder.optimize_folding(method='differential_evolution')
de_coords, de_valid = folder.decode_folding(de_angles)

print(f"DE - Best energy: {de_energy:.3f}")
print(f"DE - Valid folding: {de_valid}\n")

# Optimize using Monte Carlo with simulated annealing
print("Running Monte Carlo with Simulated Annealing...")
mc_angles, mc_energy, mc_history = folder.monte_carlo_optimization(max_iterations=5000)
mc_coords, mc_valid = folder.decode_folding(mc_angles)

print(f"MC - Best energy: {mc_energy:.3f}")
print(f"MC - Valid folding: {mc_valid}\n")

# Analysis and visualization
def analyze_folding(coordinates, sequence, title):
    """
    Analyze and display folding statistics
    """
    if coordinates is None:
        print(f"{title}: Invalid folding")
        return
        
    # Calculate statistics
    h_positions = [i for i, aa in enumerate(sequence) if aa == 'H']
    p_positions = [i for i, aa in enumerate(sequence) if aa == 'P']
    
    # Count H-H contacts
    hh_contacts = 0
    for i in range(len(sequence)):
        for j in range(i + 2, len(sequence)):
            if sequence[i] == 'H' and sequence[j] == 'H':
                dist = abs(coordinates[i][0] - coordinates[j][0]) + abs(coordinates[i][1] - coordinates[j][1])
                if dist == 1:
                    hh_contacts += 1
    
    # Calculate compactness (radius of gyration)
    center_x = np.mean([coord[0] for coord in coordinates])
    center_y = np.mean([coord[1] for coord in coordinates])
    rg = np.sqrt(np.mean([(coord[0] - center_x)**2 + (coord[1] - center_y)**2 for coord in coordinates]))
    
    print(f"{title} Analysis:")
    print(f"  H-H contacts: {hh_contacts}")
    print(f"  Energy: {-hh_contacts}")
    print(f"  Radius of gyration: {rg:.3f}")
    print(f"  Span X: {max(coord[0] for coord in coordinates) - min(coord[0] for coord in coordinates)}")
    print(f"  Span Y: {max(coord[1] for coord in coordinates) - min(coord[1] for coord in coordinates)}")
    print()

# Analyze both solutions
analyze_folding(de_coords, sequence, "Differential Evolution")
analyze_folding(mc_coords, sequence, "Monte Carlo")

# Create comprehensive visualization
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('Protein Folding Optimization Results', fontsize=16, fontweight='bold')

# Plot 1: DE Folding Structure
if de_coords and de_valid:
    ax = axes[0, 0]
    x_coords = [coord[0] for coord in de_coords]
    y_coords = [coord[1] for coord in de_coords]
    
    # Plot backbone
    ax.plot(x_coords, y_coords, 'k-', linewidth=2, alpha=0.7, label='Backbone')
    
    # Plot amino acids with different colors
    for i, (x, y) in enumerate(de_coords):
        color = 'red' if sequence[i] == 'H' else 'blue'
        marker = 'o' if sequence[i] == 'H' else 's'
        ax.scatter(x, y, c=color, s=150, marker=marker, 
                  edgecolors='black', linewidth=1, zorder=5)
        ax.annotate(f'{i}', (x, y), xytext=(5, 5), textcoords='offset points', 
                   fontsize=8, ha='left')
    
    ax.set_title(f'Differential Evolution\nEnergy: {de_energy:.1f}', fontweight='bold')
    ax.set_xlabel('X coordinate')
    ax.set_ylabel('Y coordinate')
    ax.grid(True, alpha=0.3)
    ax.legend(['Backbone', 'Hydrophobic (H)', 'Polar (P)'])
    ax.set_aspect('equal')

# Plot 2: MC Folding Structure  
if mc_coords and mc_valid:
    ax = axes[0, 1]
    x_coords = [coord[0] for coord in mc_coords]
    y_coords = [coord[1] for coord in mc_coords]
    
    # Plot backbone
    ax.plot(x_coords, y_coords, 'k-', linewidth=2, alpha=0.7)
    
    # Plot amino acids
    for i, (x, y) in enumerate(mc_coords):
        color = 'red' if sequence[i] == 'H' else 'blue'
        marker = 'o' if sequence[i] == 'H' else 's'
        ax.scatter(x, y, c=color, s=150, marker=marker, 
                  edgecolors='black', linewidth=1, zorder=5)
        ax.annotate(f'{i}', (x, y), xytext=(5, 5), textcoords='offset points', 
                   fontsize=8, ha='left')
    
    ax.set_title(f'Monte Carlo\nEnergy: {mc_energy:.1f}', fontweight='bold')
    ax.set_xlabel('X coordinate')
    ax.set_ylabel('Y coordinate')
    ax.grid(True, alpha=0.3)
    ax.set_aspect('equal')

# Plot 3: Convergence Comparison
ax = axes[0, 2]
if len(de_history) > 0:
    ax.plot(range(len(de_history)), de_history, 'g-', linewidth=2, label='Differential Evolution')
if len(mc_history) > 0:
    # Sample MC history for better visualization
    sample_indices = np.linspace(0, len(mc_history)-1, min(200, len(mc_history)), dtype=int)
    sampled_history = [mc_history[i] for i in sample_indices]
    ax.plot(sample_indices, sampled_history, 'r-', linewidth=2, label='Monte Carlo')

ax.set_title('Optimization Convergence', fontweight='bold')
ax.set_xlabel('Iteration')
ax.set_ylabel('Best Energy Found')
ax.legend()
ax.grid(True, alpha=0.3)

# Plot 4: Energy Landscape Analysis
ax = axes[1, 0]
# Generate random foldings to show energy distribution
random_energies = []
for _ in range(1000):
    random_angles = [random.uniform(0, 3.99) for _ in range(len(sequence) - 1)]
    energy = folder.objective_function(random_angles)
    if energy < 100:  # Only valid foldings
        random_energies.append(energy)

if random_energies:
    ax.hist(random_energies, bins=30, alpha=0.7, color='skyblue', edgecolor='black')
    ax.axvline(de_energy, color='green', linestyle='--', linewidth=2, label=f'DE: {de_energy:.1f}')
    ax.axvline(mc_energy, color='red', linestyle='--', linewidth=2, label=f'MC: {mc_energy:.1f}')
    ax.set_title('Energy Distribution\n(Random vs Optimized)', fontweight='bold')
    ax.set_xlabel('Energy')
    ax.set_ylabel('Frequency')
    ax.legend()
    ax.grid(True, alpha=0.3)

# Plot 5: Sequence Analysis
ax = axes[1, 1]
positions = list(range(len(sequence)))
colors = ['red' if aa == 'H' else 'blue' for aa in sequence]
bars = ax.bar(positions, [1]*len(sequence), color=colors, alpha=0.7, edgecolor='black')

ax.set_title('Protein Sequence', fontweight='bold')
ax.set_xlabel('Position')
ax.set_ylabel('Amino Acid Type')
ax.set_xticks(positions[::2])
ax.set_xticklabels(positions[::2])
ax.set_yticks([0.5, 1])
ax.set_yticklabels(['', 'H/P'])

# Add legend
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor='red', alpha=0.7, label='Hydrophobic (H)'),
                  Patch(facecolor='blue', alpha=0.7, label='Polar (P)')]
ax.legend(handles=legend_elements)

# Plot 6: Contact Map
ax = axes[1, 2]
contact_matrix = np.zeros((len(sequence), len(sequence)))

# Fill contact matrix for the best solution (MC in this case)
best_coords = mc_coords if mc_valid else de_coords
if best_coords:
    for i in range(len(sequence)):
        for j in range(len(sequence)):
            if i != j:
                dist = abs(best_coords[i][0] - best_coords[j][0]) + abs(best_coords[i][1] - best_coords[j][1])
                if dist == 1:  # Adjacent in space
                    contact_matrix[i, j] = 1

im = ax.imshow(contact_matrix, cmap='RdYlBu_r', aspect='equal')
ax.set_title('Contact Map\n(Best Solution)', fontweight='bold')
ax.set_xlabel('Residue Index')
ax.set_ylabel('Residue Index')

# Add colorbar
cbar = plt.colorbar(im, ax=ax, shrink=0.8)
cbar.set_label('Contact')

plt.tight_layout()
plt.show()

# Performance comparison
print("=== Performance Comparison ===")
print(f"Differential Evolution:")
print(f"  Final Energy: {de_energy:.3f}")
print(f"  Convergence Steps: {len(de_history)}")

print(f"Monte Carlo + Simulated Annealing:")
print(f"  Final Energy: {mc_energy:.3f}")
print(f"  Total Steps: {len(mc_history)}")

# Determine the winner
if de_valid and mc_valid:
    if de_energy < mc_energy:
        print(f"\n🏆 Winner: Differential Evolution (Energy: {de_energy:.3f})")
    elif mc_energy < de_energy:
        print(f"\n🏆 Winner: Monte Carlo (Energy: {mc_energy:.3f})")
    else:
        print(f"\n🤝 Tie! Both methods achieved energy: {de_energy:.3f}")
else:
    print(f"\n⚠️  Some optimizations produced invalid foldings")

print("\n=== Mathematical Analysis ===")
print("Energy Function:")
print("E = -∑(i,j) δ_HiHj * contact(i,j)")
print("where δ_HiHj = 1 if both residues i,j are hydrophobic and adjacent")
print("Lower energy indicates better folding (more H-H contacts)")

def calculate_energy(self, coordinates):
    energy = 0
    for i in range(self.length):
        for j in range(i + 2, self.length):  # Skip adjacent in sequence
            if self.sequence[i] == 'H' and self.sequence[j] == 'H':
                dist = abs(coordinates[i][0] - coordinates[j][0]) + abs(coordinates[i][1] - coordinates[j][1])
                if dist == 1:  # Manhattan distance = 1 means adjacent
                    energy -= 1  # Favorable interaction
    return energy

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar

# Constants
CYCLE_TIME = 60  # Total signal cycle time (seconds)
lambda_NS = 0.6  # Arrival rate per second from North-South
lambda_EW = 0.4  # Arrival rate per second from East-West

# Cost function to minimize: total waiting time
def total_waiting_time(g):
    if g <= 0 or g >= CYCLE_TIME:
        return np.inf  # avoid invalid green times
    red_NS = CYCLE_TIME - g
    red_EW = g
    wait_NS = (lambda_NS * red_NS ** 2) / (2 * g)
    wait_EW = (lambda_EW * red_EW ** 2) / (2 * (CYCLE_TIME - g))
    return wait_NS + wait_EW

# Find the optimal green time for NS direction
result = minimize_scalar(total_waiting_time, bounds=(1, CYCLE_TIME - 1), method='bounded')

optimal_g = result.x
min_wait = result.fun

print(f"✅ Optimal green time for NS: {optimal_g:.2f} seconds")
print(f"⏱️ Minimum total waiting time: {min_wait:.2f} vehicle-seconds")

green_times = np.linspace(1, CYCLE_TIME - 1, 100)
waiting_times = [total_waiting_time(g) for g in green_times]

plt.figure(figsize=(10, 6))
plt.plot(green_times, waiting_times, label='Total Waiting Time', color='blue')
plt.axvline(optimal_g, color='red', linestyle='--', label=f'Optimal g = {optimal_g:.2f}s')
plt.title('🚦 Traffic Signal Optimization: Total Waiting Time vs Green Time')
plt.xlabel('Green Time for NS (seconds)')
plt.ylabel('Total Waiting Time (vehicle-seconds)')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import linprog
import warnings
warnings.filterwarnings('ignore')

# Set up the problem parameters
class BasketballOptimizer:
    def __init__(self):
        # Player data: [Name, Offensive Rating, Defensive Rating, Fatigue Factor, Max Minutes]
        self.players_data = [
            ['LeBron James', 118, 110, 0.8, 42],
            ['Stephen Curry', 125, 105, 1.2, 40],
            ['Kawhi Leonard', 115, 115, 1.0, 35],
            ['Giannis Antetokounmpo', 120, 112, 0.9, 38],
            ['Kevin Durant', 122, 108, 1.1, 36],
            ['Luka Doncic', 116, 104, 1.3, 37],
            ['Jayson Tatum', 114, 109, 1.0, 39],
            ['Joel Embiid', 119, 113, 1.4, 32],
            ['Nikola Jokic', 121, 106, 0.7, 36],
            ['Jimmy Butler', 113, 111, 0.9, 35]
        ]
        
        # Create DataFrame
        self.df = pd.DataFrame(self.players_data, 
                              columns=['Player', 'Off_Rating', 'Def_Rating', 'Fatigue_Factor', 'Max_Minutes'])
        
        # Game parameters
        self.game_duration = 48  # minutes
        self.time_periods = 12   # 4-minute periods
        self.players_on_court = 5
        
    def calculate_performance_matrix(self):
        """Calculate performance for each player in each time period"""
        n_players = len(self.df)
        performance_matrix = np.zeros((n_players, self.time_periods))
        
        for i in range(n_players):
            base_performance = (self.df.iloc[i]['Off_Rating'] + self.df.iloc[i]['Def_Rating']) / 2
            fatigue_factor = self.df.iloc[i]['Fatigue_Factor']
            
            for t in range(self.time_periods):
                # Performance decreases with fatigue over time
                fatigue_penalty = fatigue_factor * (t + 1) * 0.5
                performance_matrix[i, t] = max(base_performance - fatigue_penalty, 50)
                
        return performance_matrix
    
    def solve_optimization(self):
        """Solve the player rotation optimization problem"""
        n_players = len(self.df)
        n_periods = self.time_periods
        
        # Get performance matrix
        performance_matrix = self.calculate_performance_matrix()
        
        # Create decision variables: x[i,t] for player i in period t
        # We'll use a simple greedy approach for demonstration
        rotation_schedule = np.zeros((n_players, n_periods))
        player_minutes = np.zeros(n_players)
        
        # For each time period, select top 5 available players
        for t in range(n_periods):
            # Calculate current performance considering fatigue and availability
            current_scores = []
            for i in range(n_players):
                if player_minutes[i] < self.df.iloc[i]['Max_Minutes']:
                    current_scores.append((performance_matrix[i, t], i))
                else:
                    current_scores.append((0, i))
            
            # Sort by performance and select top 5
            current_scores.sort(reverse=True)
            selected_players = [idx for _, idx in current_scores[:self.players_on_court]]
            
            # Update rotation schedule and minutes
            for player_idx in selected_players:
                rotation_schedule[player_idx, t] = 1
                player_minutes[player_idx] += 4  # 4 minutes per period
        
        return rotation_schedule, performance_matrix
    
    def analyze_results(self, rotation_schedule, performance_matrix):
        """Analyze optimization results"""
        results = {}
        
        # Calculate total minutes per player
        total_minutes = np.sum(rotation_schedule, axis=1) * 4
        
        # Calculate total performance
        total_performance = np.sum(rotation_schedule * performance_matrix)
        
        # Create results dataframe
        results_df = self.df.copy()
        results_df['Minutes_Played'] = total_minutes
        results_df['Utilization'] = (total_minutes / results_df['Max_Minutes'] * 100).round(1)
        
        # Calculate performance contribution
        performance_contribution = np.sum(rotation_schedule * performance_matrix, axis=1)
        results_df['Performance_Contribution'] = performance_contribution.round(1)
        
        results['summary'] = results_df
        results['total_performance'] = total_performance
        results['rotation_matrix'] = rotation_schedule
        
        return results
    
    def create_visualizations(self, results):
        """Create comprehensive visualizations"""
        fig, axes = plt.subplots(2, 2, figsize=(16, 12))
        
        # 1. Player Minutes Distribution
        ax1 = axes[0, 0]
        players = results['summary']['Player']
        minutes = results['summary']['Minutes_Played']
        max_minutes = results['summary']['Max_Minutes']
        
        x_pos = np.arange(len(players))
        bars1 = ax1.bar(x_pos, minutes, alpha=0.7, label='Minutes Played', color='skyblue')
        bars2 = ax1.bar(x_pos, max_minutes, alpha=0.3, label='Max Minutes', color='red')
        
        ax1.set_xlabel('Players')
        ax1.set_ylabel('Minutes')
        ax1.set_title('Player Minutes: Played vs Maximum Available')
        ax1.set_xticks(x_pos)
        ax1.set_xticklabels([name.split()[-1] for name in players], rotation=45)
        ax1.legend()
        ax1.grid(True, alpha=0.3)
        
        # 2. Performance Contribution
        ax2 = axes[0, 1]
        performance = results['summary']['Performance_Contribution']
        bars = ax2.bar(range(len(players)), performance, color='lightgreen', alpha=0.7)
        
        ax2.set_xlabel('Players')
        ax2.set_ylabel('Performance Contribution')
        ax2.set_title('Individual Player Performance Contribution')
        ax2.set_xticks(range(len(players)))
        ax2.set_xticklabels([name.split()[-1] for name in players], rotation=45)
        ax2.grid(True, alpha=0.3)
        
        # Add value labels on bars
        for bar, value in zip(bars, performance):
            ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1,
                    f'{value:.0f}', ha='center', va='bottom')
        
        # 3. Rotation Schedule Heatmap
        ax3 = axes[1, 0]
        rotation_matrix = results['rotation_matrix']
        
        # Create heatmap
        im = ax3.imshow(rotation_matrix, cmap='YlOrRd', aspect='auto')
        ax3.set_xlabel('Time Periods (4-min each)')
        ax3.set_ylabel('Players')
        ax3.set_title('Player Rotation Schedule')
        ax3.set_xticks(range(self.time_periods))
        ax3.set_xticklabels([f'P{i+1}' for i in range(self.time_periods)])
        ax3.set_yticks(range(len(players)))
        ax3.set_yticklabels([name.split()[-1] for name in players])
        
        # Add colorbar
        plt.colorbar(im, ax=ax3, label='Playing (1) / Resting (0)')
        
        # 4. Utilization Rate
        ax4 = axes[1, 1]
        utilization = results['summary']['Utilization']
        colors = ['red' if x > 90 else 'orange' if x > 75 else 'green' for x in utilization]
        
        bars = ax4.bar(range(len(players)), utilization, color=colors, alpha=0.7)
        ax4.axhline(y=100, color='red', linestyle='--', alpha=0.5, label='Max Capacity')
        ax4.axhline(y=75, color='orange', linestyle='--', alpha=0.5, label='High Utilization')
        
        ax4.set_xlabel('Players')
        ax4.set_ylabel('Utilization Rate (%)')
        ax4.set_title('Player Utilization Rate')
        ax4.set_xticks(range(len(players)))
        ax4.set_xticklabels([name.split()[-1] for name in players], rotation=45)
        ax4.legend()
        ax4.grid(True, alpha=0.3)
        
        # Add percentage labels
        for bar, value in zip(bars, utilization):
            ax4.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1,
                    f'{value:.1f}%', ha='center', va='bottom')
        
        plt.tight_layout()
        plt.show()
        
        return fig
    
    def generate_insights(self, results):
        """Generate strategic insights from the optimization"""
        summary = results['summary']
        
        print("=== BASKETBALL TEAM ROTATION OPTIMIZATION RESULTS ===\n")
        
        print("📊 PERFORMANCE SUMMARY:")
        print(f"Total Team Performance Score: {results['total_performance']:.1f}")
        print(f"Average Player Utilization: {summary['Utilization'].mean():.1f}%")
        print(f"Players at >90% Utilization: {sum(summary['Utilization'] > 90)}")
        
        print("\n🏀 TOP PERFORMERS:")
        top_performers = summary.nlargest(3, 'Performance_Contribution')
        for idx, (_, player) in enumerate(top_performers.iterrows()):
            print(f"{idx+1}. {player['Player']}: {player['Performance_Contribution']:.1f} points "
                  f"({player['Minutes_Played']} min, {player['Utilization']:.1f}% utilization)")
        
        print("\n⚠️  STRATEGIC RECOMMENDATIONS:")
        
        # High utilization players
        high_util = summary[summary['Utilization'] > 85]
        if len(high_util) > 0:
            print("• High Utilization Alert:")
            for _, player in high_util.iterrows():
                print(f"  - {player['Player']}: {player['Utilization']:.1f}% utilization - consider rest periods")
        
        # Underutilized players
        low_util = summary[summary['Utilization'] < 60]
        if len(low_util) > 0:
            print("• Underutilized Players:")
            for _, player in low_util.iterrows():
                print(f"  - {player['Player']}: {player['Utilization']:.1f}% utilization - opportunity for more minutes")
        
        # Performance efficiency
        summary['Efficiency'] = summary['Performance_Contribution'] / summary['Minutes_Played']
        top_efficiency = summary.nlargest(3, 'Efficiency')
        print("\n💡 MOST EFFICIENT PLAYERS (Performance per Minute):")
        for idx, (_, player) in enumerate(top_efficiency.iterrows()):
            print(f"{idx+1}. {player['Player']}: {player['Efficiency']:.2f} points/minute")
        
        return summary

# Execute the optimization
def main():
    # Initialize optimizer
    optimizer = BasketballOptimizer()
    
    print("🏀 Basketball Team Rotation Optimization")
    print("=" * 50)
    
    # Display initial team data
    print("\n📋 TEAM ROSTER:")
    print(optimizer.df.to_string(index=False))
    
    # Solve optimization
    print("\n🔄 Solving optimization problem...")
    rotation_schedule, performance_matrix = optimizer.solve_optimization()
    
    # Analyze results
    results = optimizer.analyze_results(rotation_schedule, performance_matrix)
    
    # Generate insights
    insights = optimizer.generate_insights(results)
    
    # Create visualizations
    print("\n📈 Generating visualizations...")
    fig = optimizer.create_visualizations(results)
    
    # Additional analysis: Performance over time
    plt.figure(figsize=(14, 8))
    
    # Plot team performance over time periods
    plt.subplot(2, 1, 1)
    period_performance = np.sum(rotation_schedule * performance_matrix, axis=0)
    periods = [f'Q{(i//3)+1}-{(i%3)*4+4}min' for i in range(optimizer.time_periods)]
    
    plt.plot(periods, period_performance, marker='o', linewidth=2, markersize=8, color='blue')
    plt.title('Team Performance Throughout the Game', fontsize=14, fontweight='bold')
    plt.xlabel('Game Period')
    plt.ylabel('Team Performance Score')
    plt.grid(True, alpha=0.3)
    plt.xticks(rotation=45)
    
    # Plot active players per period
    plt.subplot(2, 1, 2)
    active_players = np.sum(rotation_schedule, axis=0)
    plt.bar(periods, active_players, alpha=0.7, color='lightcoral')
    plt.title('Number of Active Players per Period', fontsize=14, fontweight='bold')
    plt.xlabel('Game Period')
    plt.ylabel('Active Players')
    plt.ylim(0, 6)
    plt.grid(True, alpha=0.3)
    plt.xticks(rotation=45)
    
    # Add horizontal line at 5 players
    plt.axhline(y=5, color='red', linestyle='--', alpha=0.7, label='Required Players')
    plt.legend()
    
    plt.tight_layout()
    plt.show()
    
    return optimizer, results

# Run the main function
if __name__ == "__main__":
    optimizer, results = main()