import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_bvp

# Parameters
alpha = 1.0   # prey birth rate
beta = 0.1    # predation rate
delta = 0.075 # predator reproduction rate
gamma = 1.5   # predator death rate
rho = 0.01    # control penalty
x_d, y_d = 30, 5  # desired population levels

T = 10
n_points = 100
t = np.linspace(0, T, n_points)

def odes(t, Y):
    x, y, lam1, lam2 = Y
    u = -lam1 / (2 * rho)
    dx = alpha * x - beta * x * y + u
    dy = delta * x * y - gamma * y
    dlam1 = - (2 * (x - x_d) + lam1 * (alpha - beta * y) + lam2 * delta * y)
    dlam2 = - (2 * (y - y_d) + lam1 * (-beta * x) + lam2 * (delta * x - gamma))
    return np.vstack([dx, dy, dlam1, dlam2])

def bc(Ya, Yb):
    return [Ya[0] - 40, Ya[1] - 9,  # initial conditions
            Yb[2], Yb[3]]          # final costate free

# Initial guess for solution
Y_guess = np.zeros((4, t.size))
Y_guess[0] = 40
Y_guess[1] = 9

sol = solve_bvp(odes, bc, t, Y_guess)

# Extract solutions
t_plot = sol.x
x, y, lam1, lam2 = sol.y
u = -lam1 / (2 * rho)

plt.figure(figsize=(14, 6))

# Population plots
plt.subplot(1, 2, 1)
plt.plot(t_plot, x, label='Prey (x)', lw=2)
plt.plot(t_plot, y, label='Predator (y)', lw=2)
plt.axhline(x_d, color='gray', ls='--', label='x desired')
plt.axhline(y_d, color='gray', ls='-.', label='y desired')
plt.xlabel('Time')
plt.ylabel('Population')
plt.title('Population Dynamics')
plt.legend()
plt.grid(True)

# Control input plot
plt.subplot(1, 2, 2)
plt.plot(t_plot, u, label='Control input u(t)', color='purple')
plt.xlabel('Time')
plt.ylabel('Control effort')
plt.title('Optimal Control Over Time')
plt.grid(True)
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
import warnings
warnings.filterwarnings('ignore')

# Material and geometric properties
E = 200e9  # Young's modulus (Pa)
rho = 7850  # Density (kg/m³)
sigma_allow = 200e6  # Allowable stress (Pa)
L = 1.0  # Beam length (m)
b = 0.05  # Beam width (m)
F = 1000  # Applied force (N)

# Discretization
n_elements = 20
x = np.linspace(0, L, n_elements + 1)
dx = L / n_elements

print("=== Cantilever Beam Shape Optimization ===")
print(f"Beam length: {L} m")
print(f"Applied force: {F} N")
print(f"Allowable stress: {sigma_allow/1e6:.0f} MPa")
print(f"Material: Steel (E = {E/1e9:.0f} GPa, ρ = {rho} kg/m³)")
print(f"Beam width: {b*1000:.0f} mm")
print(f"Number of elements: {n_elements}")

def bending_moment(x_pos):
    """Calculate bending moment at position x"""
    return F * (L - x_pos)

def second_moment_area(height):
    """Calculate second moment of area for rectangular cross-section"""
    return b * height**3 / 12

def stress_at_position(x_pos, height):
    """Calculate maximum stress at position x"""
    M = bending_moment(x_pos)
    I = second_moment_area(height)
    return M * height / 2 / I

def objective_function(h):
    """Objective function: minimize weight"""
    # Weight = density × volume
    # Volume = width × ∫(height dx) ≈ width × Σ(height × dx)
    weight = rho * b * np.sum(h * dx)
    return weight

def constraint_function(h):
    """Constraint function: stress constraint"""
    constraints = []
    
    # Stress constraint at each node
    for i, x_pos in enumerate(x):
        if h[i] > 0:  # Avoid division by zero
            stress = stress_at_position(x_pos, h[i])
            # Constraint: stress - allowable_stress ≤ 0
            constraints.append(stress - sigma_allow)
    
    return np.array(constraints)

def solve_optimization():
    """Solve the shape optimization problem"""
    
    # Initial guess: uniform height
    h0 = np.ones(n_elements + 1) * 0.1  # 100 mm initial height
    
    # Bounds: minimum height 10 mm, maximum height 500 mm
    bounds = [(0.01, 0.5) for _ in range(n_elements + 1)]
    
    # Constraint dictionary for scipy.optimize
    constraints = {'type': 'ineq', 
                   'fun': lambda h: -constraint_function(h)}  # Convert ≤ 0 to ≥ 0
    
    print("\n=== Optimization Process ===")
    print("Starting optimization...")
    print(f"Initial weight: {objective_function(h0):.2f} kg")
    
    # Solve optimization problem
    result = minimize(objective_function, h0, 
                     method='SLSQP', 
                     bounds=bounds, 
                     constraints=constraints,
                     options={'maxiter': 1000, 'ftol': 1e-9})
    
    return result

# Solve the optimization problem
result = solve_optimization()

# Extract optimal solution
h_optimal = result.x
weight_optimal = result.fun

print(f"\nOptimization {'converged' if result.success else 'failed'}")
print(f"Final weight: {weight_optimal:.2f} kg")
print(f"Weight reduction: {((objective_function(np.ones(n_elements + 1) * 0.1) - weight_optimal) / objective_function(np.ones(n_elements + 1) * 0.1) * 100):.1f}%")

# Calculate stress distribution for optimal shape
stress_optimal = np.zeros(n_elements + 1)
moment_optimal = np.zeros(n_elements + 1)

for i, x_pos in enumerate(x):
    moment_optimal[i] = bending_moment(x_pos)
    if h_optimal[i] > 0:
        stress_optimal[i] = stress_at_position(x_pos, h_optimal[i])

print(f"Maximum stress: {np.max(stress_optimal)/1e6:.1f} MPa")
print(f"Stress ratio: {np.max(stress_optimal)/sigma_allow:.3f}")

# Theoretical optimal solution (for comparison)
# For a cantilever beam, the optimal shape follows: h(x) ∝ √(L-x)
x_theory = np.linspace(0, L, 100)
# Scale factor to match stress constraint at the fixed end
scale_factor = np.sqrt(12 * F * L / (b * sigma_allow))
h_theory = scale_factor * np.sqrt(L - x_theory)
h_theory[h_theory < 0.01] = 0.01  # Minimum thickness constraint

print(f"\n=== Comparison with Theory ===")
print(f"Theoretical optimal shape: h(x) ∝ √(L-x)")
print(f"Scale factor: {scale_factor:.3f}")

# Create comprehensive plots
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))

# Plot 1: Beam shape comparison
ax1.plot(x, h_optimal*1000, 'b-o', linewidth=2, markersize=4, label='Optimized shape')
ax1.plot(x_theory, h_theory*1000, 'r--', linewidth=2, label='Theoretical optimal')
ax1.fill_between(x, 0, h_optimal*1000, alpha=0.3, color='blue')
ax1.set_xlabel('Position along beam (m)')
ax1.set_ylabel('Height (mm)')
ax1.set_title('Optimal Beam Shape')
ax1.grid(True, alpha=0.3)
ax1.legend()
ax1.set_ylim(0, max(h_optimal)*1000*1.1)

# Plot 2: Bending moment diagram
ax2.plot(x, moment_optimal, 'g-', linewidth=2, label='Bending moment')
ax2.fill_between(x, 0, moment_optimal, alpha=0.3, color='green')
ax2.set_xlabel('Position along beam (m)')
ax2.set_ylabel('Bending Moment (N⋅m)')
ax2.set_title('Bending Moment Distribution')
ax2.grid(True, alpha=0.3)
ax2.legend()

# Plot 3: Stress distribution
ax3.plot(x, stress_optimal/1e6, 'r-o', linewidth=2, markersize=4, label='Stress')
ax3.axhline(y=sigma_allow/1e6, color='k', linestyle='--', linewidth=2, label='Allowable stress')
ax3.fill_between(x, 0, stress_optimal/1e6, alpha=0.3, color='red')
ax3.set_xlabel('Position along beam (m)')
ax3.set_ylabel('Stress (MPa)')
ax3.set_title('Stress Distribution')
ax3.grid(True, alpha=0.3)
ax3.legend()
ax3.set_ylim(0, max(stress_optimal/1e6)*1.1)

# Plot 4: Convergence and comparison metrics
uniform_height = 0.1  # 100 mm uniform beam
uniform_weight = rho * b * L * uniform_height
weights = [uniform_weight, weight_optimal]
labels = ['Uniform\nBeam', 'Optimized\nBeam']
colors = ['lightcoral', 'lightblue']

bars = ax4.bar(labels, weights, color=colors, alpha=0.7, edgecolor='black')
ax4.set_ylabel('Weight (kg)')
ax4.set_title('Weight Comparison')
ax4.grid(True, alpha=0.3, axis='y')

# Add value labels on bars
for bar, weight in zip(bars, weights):
    height = bar.get_height()
    ax4.text(bar.get_x() + bar.get_width()/2., height + height*0.01,
             f'{weight:.2f} kg', ha='center', va='bottom', fontweight='bold')

plt.tight_layout()
plt.show()

# Additional analysis
print(f"\n=== Detailed Analysis ===")
print(f"Uniform beam weight: {uniform_weight:.2f} kg")
print(f"Optimized beam weight: {weight_optimal:.2f} kg")
print(f"Material savings: {uniform_weight - weight_optimal:.2f} kg")
print(f"Efficiency gain: {(uniform_weight - weight_optimal)/uniform_weight*100:.1f}%")

# Print height distribution
print(f"\n=== Height Distribution (mm) ===")
for i in range(0, len(x), 5):  # Print every 5th point
    print(f"x = {x[i]:.2f} m: h = {h_optimal[i]*1000:.1f} mm")

# Verify constraint satisfaction
max_stress_ratio = np.max(stress_optimal) / sigma_allow
print(f"\n=== Constraint Verification ===")
print(f"Maximum stress ratio: {max_stress_ratio:.4f}")
print(f"Constraint satisfied: {'Yes' if max_stress_ratio <= 1.0 else 'No'}")

# Calculate theoretical weight for comparison
weight_theory = rho * b * np.trapz(h_theory, x_theory)
print(f"Theoretical optimal weight: {weight_theory:.2f} kg")
print(f"Numerical vs theoretical difference: {abs(weight_optimal - weight_theory)/weight_theory*100:.1f}%")

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar, minimize
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D

# Set up the problem parameters
class ReactorOptimization:
    def __init__(self):
        # Reaction parameters
        self.k = 0.5  # reaction rate constant (1/min)
        self.C_A0 = 10.0  # inlet concentration (mol/L)
        self.q = 100.0  # volumetric flow rate (L/min)
        
        # Cost parameters
        self.alpha = 0.1  # cost coefficient for reactor volume ($/L·min)
        self.beta = 2.0   # cost coefficient for raw material ($/mol)
        
    def concentration_steady_state(self, V):
        """Calculate steady-state concentration of A"""
        tau = V / self.q  # residence time
        C_A = self.C_A0 / (1 + self.k * tau)
        return C_A
    
    def conversion(self, V):
        """Calculate conversion of reactant A"""
        C_A = self.concentration_steady_state(V)
        X_A = (self.C_A0 - C_A) / self.C_A0
        return X_A
    
    def total_cost(self, V):
        """Calculate total cost function"""
        C_A = self.concentration_steady_state(V)
        operating_cost = self.alpha * V
        material_cost = self.beta * self.q * C_A
        return operating_cost + material_cost
    
    def optimize_reactor_volume(self):
        """Find optimal reactor volume"""
        # Define bounds for reactor volume (reasonable range)
        bounds = (1, 1000)  # L
        
        # Minimize the total cost
        result = minimize_scalar(self.total_cost, bounds=bounds, method='bounded')
        
        return result

# Create reactor optimization instance
reactor = ReactorOptimization()

# Solve the optimization problem
optimization_result = reactor.optimize_reactor_volume()
optimal_volume = optimization_result.x
optimal_cost = optimization_result.fun

print("=== Reactor Design Optimization Results ===")
print(f"Optimal reactor volume: {optimal_volume:.2f} L")
print(f"Minimum total cost: ${optimal_cost:.2f}/min")
print(f"Optimal residence time: {optimal_volume/reactor.q:.3f} min")
print(f"Steady-state concentration: {reactor.concentration_steady_state(optimal_volume):.3f} mol/L")
print(f"Conversion at optimal conditions: {reactor.conversion(optimal_volume):.3f}")

# Generate data for visualization
V_range = np.linspace(1, 500, 1000)
C_A_values = [reactor.concentration_steady_state(V) for V in V_range]
X_A_values = [reactor.conversion(V) for V in V_range]
cost_values = [reactor.total_cost(V) for V in V_range]
operating_cost_values = [reactor.alpha * V for V in V_range]
material_cost_values = [reactor.beta * reactor.q * reactor.concentration_steady_state(V) for V in V_range]

# Create comprehensive visualization
fig = plt.figure(figsize=(15, 12))

# Plot 1: Concentration vs Volume
plt.subplot(2, 3, 1)
plt.plot(V_range, C_A_values, 'b-', linewidth=2, label='$C_A$')
plt.axvline(optimal_volume, color='r', linestyle='--', alpha=0.7, label=f'Optimal V = {optimal_volume:.1f} L')
plt.xlabel('Reactor Volume (L)')
plt.ylabel('Concentration (mol/L)')
plt.title('Steady-State Concentration vs Volume')
plt.grid(True, alpha=0.3)
plt.legend()

# Plot 2: Conversion vs Volume
plt.subplot(2, 3, 2)
plt.plot(V_range, X_A_values, 'g-', linewidth=2, label='Conversion')
plt.axvline(optimal_volume, color='r', linestyle='--', alpha=0.7, label=f'Optimal V = {optimal_volume:.1f} L')
plt.xlabel('Reactor Volume (L)')
plt.ylabel('Conversion')
plt.title('Conversion vs Volume')
plt.grid(True, alpha=0.3)
plt.legend()

# Plot 3: Cost Analysis
plt.subplot(2, 3, 3)
plt.plot(V_range, cost_values, 'k-', linewidth=2, label='Total Cost')
plt.plot(V_range, operating_cost_values, 'b--', linewidth=1.5, label='Operating Cost')
plt.plot(V_range, material_cost_values, 'r--', linewidth=1.5, label='Material Cost')
plt.axvline(optimal_volume, color='r', linestyle='--', alpha=0.7)
plt.axhline(optimal_cost, color='r', linestyle='--', alpha=0.7)
plt.plot(optimal_volume, optimal_cost, 'ro', markersize=8, label=f'Optimum: ${optimal_cost:.1f}/min')
plt.xlabel('Reactor Volume (L)')
plt.ylabel('Cost ($/min)')
plt.title('Cost Analysis')
plt.grid(True, alpha=0.3)
plt.legend()

# Plot 4: Residence Time vs Conversion
tau_range = V_range / reactor.q
plt.subplot(2, 3, 4)
plt.plot(tau_range, X_A_values, 'purple', linewidth=2)
plt.axvline(optimal_volume/reactor.q, color='r', linestyle='--', alpha=0.7, 
           label=f'Optimal τ = {optimal_volume/reactor.q:.3f} min')
plt.xlabel('Residence Time (min)')
plt.ylabel('Conversion')
plt.title('Conversion vs Residence Time')
plt.grid(True, alpha=0.3)
plt.legend()

# Plot 5: Sensitivity Analysis - Effect of k
plt.subplot(2, 3, 5)
k_values = np.linspace(0.1, 1.0, 5)
colors = plt.cm.viridis(np.linspace(0, 1, len(k_values)))

for i, k_val in enumerate(k_values):
    reactor_temp = ReactorOptimization()
    reactor_temp.k = k_val
    cost_temp = [reactor_temp.total_cost(V) for V in V_range]
    plt.plot(V_range, cost_temp, color=colors[i], linewidth=2, 
             label=f'k = {k_val:.1f} min⁻¹')

plt.xlabel('Reactor Volume (L)')
plt.ylabel('Total Cost ($/min)')
plt.title('Sensitivity to Reaction Rate Constant')
plt.grid(True, alpha=0.3)
plt.legend()

# Plot 6: 3D Surface Plot
plt.subplot(2, 3, 6)
ax = fig.add_subplot(2, 3, 6, projection='3d')

# Create meshgrid for 3D plot
V_mesh = np.linspace(50, 300, 30)
k_mesh = np.linspace(0.2, 1.0, 30)
V_grid, k_grid = np.meshgrid(V_mesh, k_mesh)
Z_grid = np.zeros_like(V_grid)

for i in range(len(k_mesh)):
    for j in range(len(V_mesh)):
        reactor_temp = ReactorOptimization()
        reactor_temp.k = k_grid[i, j]
        Z_grid[i, j] = reactor_temp.total_cost(V_grid[i, j])

surf = ax.plot_surface(V_grid, k_grid, Z_grid, cmap='viridis', alpha=0.8)
ax.set_xlabel('Volume (L)')
ax.set_ylabel('Rate Constant (1/min)')
ax.set_zlabel('Total Cost ($/min)')
ax.set_title('Cost Surface')

plt.tight_layout()
plt.show()

# Create detailed results table
results_df = pd.DataFrame({
    'Parameter': ['Optimal Volume (L)', 'Optimal Residence Time (min)', 
                  'Steady-State Concentration (mol/L)', 'Conversion', 
                  'Minimum Total Cost ($/min)', 'Operating Cost ($/min)', 
                  'Material Cost ($/min)'],
    'Value': [optimal_volume, optimal_volume/reactor.q, 
              reactor.concentration_steady_state(optimal_volume),
              reactor.conversion(optimal_volume), optimal_cost,
              reactor.alpha * optimal_volume,
              reactor.beta * reactor.q * reactor.concentration_steady_state(optimal_volume)]
})

print("\n=== Detailed Results Table ===")
print(results_df.to_string(index=False, float_format='%.3f'))

# Perform parametric study
print("\n=== Parametric Study: Effect of Flow Rate ===")
q_values = [50, 75, 100, 125, 150]
parametric_results = []

for q_val in q_values:
    reactor_param = ReactorOptimization()
    reactor_param.q = q_val
    result = reactor_param.optimize_reactor_volume()
    parametric_results.append({
        'Flow Rate (L/min)': q_val,
        'Optimal Volume (L)': result.x,
        'Optimal Residence Time (min)': result.x/q_val,
        'Conversion': reactor_param.conversion(result.x),
        'Min Cost ($/min)': result.fun
    })

parametric_df = pd.DataFrame(parametric_results)
print(parametric_df.to_string(index=False, float_format='%.3f'))

# Calculate economic metrics
print("\n=== Economic Analysis ===")
annual_cost = optimal_cost * 60 * 24 * 365  # Convert to $/year
print(f"Annual operating cost: ${annual_cost:,.0f}/year")

# Calculate payback period assuming capital cost
capital_cost = 5000 * optimal_volume  # Assume $5000/L capital cost
print(f"Estimated capital cost: ${capital_cost:,.0f}")

# Calculate production rate
production_rate = reactor.q * reactor.conversion(optimal_volume) * reactor.C_A0  # mol/min
annual_production = production_rate * 60 * 24 * 365  # mol/year
print(f"Annual production rate: {annual_production:,.0f} mol/year")
print(f"Specific production rate: {production_rate/optimal_volume:.3f} mol/(L·min)")

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linprog
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns

# Set up the optimization problem
print("=== Petroleum Refinery Optimization Problem ===")
print("Maximizing profit from gasoline, diesel, and fuel oil production")
print()

# Objective function coefficients (profits per barrel)
# We use negative values because linprog minimizes by default
c = [-80, -60, -40]  # Gasoline: $80/barrel, Diesel: $60/barrel, Fuel Oil: $40/barrel

# Inequality constraint matrix A and bounds b (Ax <= b)
A = [
    [0.3, 0.2, 0.1],  # CDU capacity constraint
    [0.4, 0.3, 0.1],  # CCU capacity constraint  
    [1, 0, 0],        # Gasoline demand constraint
    [0, 1, 0],        # Diesel demand constraint
    [0, 0, 1]         # Fuel oil demand constraint
]

b = [2000, 1800, 4000, 3000, 2000]

# Variable bounds (non-negativity)
x_bounds = [(0, None), (0, None), (0, None)]

# Solve the linear programming problem
print("Solving the optimization problem...")
result = linprog(c, A_ub=A, b_ub=b, bounds=x_bounds, method='highs')

print("\n=== OPTIMIZATION RESULTS ===")
if result.success:
    gasoline_prod = result.x[0]
    diesel_prod = result.x[1]
    fuel_oil_prod = result.x[2]
    max_profit = -result.fun  # Convert back to positive (maximization)
    
    print(f"Optimal Production Plan:")
    print(f"  Gasoline: {gasoline_prod:.2f} barrels/day")
    print(f"  Diesel: {diesel_prod:.2f} barrels/day") 
    print(f"  Fuel Oil: {fuel_oil_prod:.2f} barrels/day")
    print(f"  Maximum Daily Profit: ${max_profit:.2f}")
    
    # Calculate resource utilization
    cdu_usage = 0.3*gasoline_prod + 0.2*diesel_prod + 0.1*fuel_oil_prod
    ccu_usage = 0.4*gasoline_prod + 0.3*diesel_prod + 0.1*fuel_oil_prod
    
    print(f"\nResource Utilization:")
    print(f"  CDU Usage: {cdu_usage:.2f}/2000 ({cdu_usage/2000*100:.1f}%)")
    print(f"  CCU Usage: {ccu_usage:.2f}/1800 ({ccu_usage/1800*100:.1f}%)")
    
    # Store results for visualization
    production_data = {
        'Product': ['Gasoline', 'Diesel', 'Fuel Oil'],
        'Production (barrels/day)': [gasoline_prod, diesel_prod, fuel_oil_prod],
        'Unit Profit ($/barrel)': [80, 60, 40],
        'Total Revenue ($)': [gasoline_prod*80, diesel_prod*60, fuel_oil_prod*40]
    }
    
    utilization_data = {
        'Resource': ['CDU', 'CCU'],
        'Usage': [cdu_usage, ccu_usage],
        'Capacity': [2000, 1800],
        'Utilization (%)': [cdu_usage/2000*100, ccu_usage/1800*100]
    }
    
else:
    print("Optimization failed!")
    print(result.message)

# Sensitivity Analysis
print("\n=== SENSITIVITY ANALYSIS ===")
print("Analyzing how profit changes with different gasoline prices...")

gasoline_prices = np.linspace(60, 100, 20)
profits = []

for price in gasoline_prices:
    c_temp = [-price, -60, -40]
    result_temp = linprog(c_temp, A_ub=A, b_ub=b, bounds=x_bounds, method='highs')
    if result_temp.success:
        profits.append(-result_temp.fun)
    else:
        profits.append(0)

sensitivity_data = pd.DataFrame({
    'Gasoline_Price': gasoline_prices,
    'Max_Profit': profits
})

print("Sensitivity analysis complete. Results will be visualized in graphs.")

# Create comprehensive visualizations
fig = plt.figure(figsize=(20, 15))

# 1. Production Mix Pie Chart
ax1 = plt.subplot(2, 3, 1)
products = ['Gasoline', 'Diesel', 'Fuel Oil']
production_values = [gasoline_prod, diesel_prod, fuel_oil_prod]
colors = ['#FF6B6B', '#4ECDC4', '#45B7D1']
plt.pie(production_values, labels=products, autopct='%1.1f%%', colors=colors, startangle=90)
plt.title('Optimal Production Mix\n(Barrels per Day)', fontsize=14, fontweight='bold')

# 2. Revenue Contribution Bar Chart
ax2 = plt.subplot(2, 3, 2)
revenues = [gasoline_prod*80, diesel_prod*60, fuel_oil_prod*40]
bars = plt.bar(products, revenues, color=colors, alpha=0.8)
plt.title('Revenue Contribution by Product', fontsize=14, fontweight='bold')
plt.ylabel('Revenue ($)')
plt.xticks(rotation=45)
# Add value labels on bars
for bar, revenue in zip(bars, revenues):
    plt.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1000,
             f'${revenue:,.0f}', ha='center', va='bottom', fontweight='bold')

# 3. Resource Utilization
ax3 = plt.subplot(2, 3, 3)
resources = ['CDU', 'CCU']
usage = [cdu_usage, ccu_usage]
capacity = [2000, 1800]
utilization_pct = [u/c*100 for u, c in zip(usage, capacity)]

x_pos = np.arange(len(resources))
bars1 = plt.bar(x_pos - 0.2, usage, 0.4, label='Current Usage', color='#FF6B6B', alpha=0.8)
bars2 = plt.bar(x_pos + 0.2, capacity, 0.4, label='Total Capacity', color='#4ECDC4', alpha=0.8)

plt.title('Resource Utilization', fontsize=14, fontweight='bold')
plt.ylabel('Capacity (units)')
plt.xlabel('Processing Units')
plt.xticks(x_pos, resources)
plt.legend()

# Add percentage labels
for i, (usage_val, util_pct) in enumerate(zip(usage, utilization_pct)):
    plt.text(i - 0.2, usage_val + 50, f'{util_pct:.1f}%', 
             ha='center', va='bottom', fontweight='bold')

# 4. Sensitivity Analysis
ax4 = plt.subplot(2, 3, 4)
plt.plot(gasoline_prices, profits, 'o-', color='#45B7D1', linewidth=3, markersize=6)
plt.title('Profit Sensitivity to Gasoline Price', fontsize=14, fontweight='bold')
plt.xlabel('Gasoline Price ($/barrel)')
plt.ylabel('Maximum Daily Profit ($)')
plt.grid(True, alpha=0.3)
plt.axvline(x=80, color='red', linestyle='--', alpha=0.7, label='Current Price')
plt.legend()

# 5. Constraint Analysis (Feasible Region Visualization for 2D case)
ax5 = plt.subplot(2, 3, 5)
# For visualization, we'll fix fuel oil production and show gasoline vs diesel
x1_range = np.linspace(0, 5000, 100)
x2_range = np.linspace(0, 4000, 100)

# CDU constraint: 0.3*x1 + 0.2*x2 <= 2000 (assuming x3=0 for visualization)
x2_cdu = (2000 - 0.3*x1_range) / 0.2
# CCU constraint: 0.4*x1 + 0.3*x2 <= 1800 (assuming x3=0 for visualization)
x2_ccu = (1800 - 0.4*x1_range) / 0.3

plt.plot(x1_range, x2_cdu, label='CDU Constraint', color='red', linewidth=2)
plt.plot(x1_range, x2_ccu, label='CCU Constraint', color='blue', linewidth=2)
plt.axhline(y=3000, color='green', linestyle='--', label='Diesel Demand Limit')
plt.axvline(x=4000, color='orange', linestyle='--', label='Gasoline Demand Limit')

# Plot optimal point (projected to 2D)
plt.plot(gasoline_prod, diesel_prod, 'ro', markersize=10, label='Optimal Solution')

plt.xlim(0, 5000)
plt.ylim(0, 4000)
plt.xlabel('Gasoline Production (barrels/day)')
plt.ylabel('Diesel Production (barrels/day)')
plt.title('Feasible Region (2D Projection)', fontsize=14, fontweight='bold')
plt.legend()
plt.grid(True, alpha=0.3)

# 6. Profit Breakdown
ax6 = plt.subplot(2, 3, 6)
profit_breakdown = [gasoline_prod*80, diesel_prod*60, fuel_oil_prod*40]
cumulative_profit = np.cumsum([0] + profit_breakdown)

for i, (product, profit) in enumerate(zip(products, profit_breakdown)):
    plt.barh(0, profit, left=cumulative_profit[i], height=0.5, 
             color=colors[i], alpha=0.8, label=f'{product}: ${profit:,.0f}')

plt.title('Profit Breakdown by Product', fontsize=14, fontweight='bold')
plt.xlabel('Cumulative Profit ($)')
plt.yticks([])
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')

plt.tight_layout()
plt.show()

# Print detailed analysis
print("\n=== DETAILED ANALYSIS ===")
print(f"1. PRODUCTION EFFICIENCY:")
print(f"   - Total production: {sum(production_values):,.0f} barrels/day")
print(f"   - Gasoline dominates at {gasoline_prod/sum(production_values)*100:.1f}% of total production")
print(f"   - This reflects gasoline's higher profit margin (${80}/barrel)")

print(f"\n2. RESOURCE BOTTLENECKS:")
if ccu_usage/1800 > cdu_usage/2000:
    print(f"   - CCU is the limiting factor at {ccu_usage/1800*100:.1f}% utilization")
    print(f"   - CDU has spare capacity at {cdu_usage/2000*100:.1f}% utilization")
    print(f"   - Consider expanding CCU capacity for further optimization")
else:
    print(f"   - CDU is the limiting factor at {cdu_usage/2000*100:.1f}% utilization")
    print(f"   - CCU has spare capacity at {ccu_usage/1800*100:.1f}% utilization")
    print(f"   - Consider expanding CDU capacity for further optimization")

print(f"\n3. MARKET POSITION:")
print(f"   - Gasoline: Using {gasoline_prod/4000*100:.1f}% of market demand")
print(f"   - Diesel: Using {diesel_prod/3000*100:.1f}% of market demand") 
print(f"   - Fuel Oil: Using {fuel_oil_prod/2000*100:.1f}% of market demand")

print(f"\n4. PROFITABILITY INSIGHTS:")
print(f"   - Revenue per barrel (weighted avg): ${max_profit/sum(production_values):.2f}")
print(f"   - Gasoline contributes {gasoline_prod*80/max_profit*100:.1f}% of total profit")
print(f"   - High gasoline focus aligns with profit maximization strategy")

print(f"\n5. SENSITIVITY INSIGHTS:")
max_profit_idx = np.argmax(profits)
optimal_gas_price = gasoline_prices[max_profit_idx]
print(f"   - Current gasoline price ($80/barrel) vs optimal range")
print(f"   - Profit increases linearly with gasoline price in current range")
print(f"   - At $100/barrel gasoline price, profit would be ${profits[-1]:,.2f}")

print("\n=== OPTIMIZATION COMPLETE ===")
print("The refinery should focus on maximizing gasoline production while")
print("efficiently utilizing both processing units to achieve optimal profitability.")

result = linprog(c, A_ub=A, b_ub=b, bounds=x_bounds, method='highs')

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

# Constants
n = 1.5  # Refractive index of lens material (e.g., glass)
L = 100  # Total system length in mm
target_focal_length = 50  # Target focal length for the whole system

def lens_focal_length(R1, R2, n=1.5):
    return 1 / ((n - 1) * (1/R1 - 1/R2))

def system_focal_length(f1, f2, d):
    return 1 / (1/f1 + 1/f2 - d/(f1 * f2))

# Aberration model (simplified): Assume spot size increases with deviation from target focal length
def spot_radius(f_sys, target_f=50):
    return (f_sys - target_f)**2  # Quadratic penalty

# Optimization function: curvatures (R1_1, R2_1, R1_2, R2_2)
def objective(x):
    R1_1, R2_1, R1_2, R2_2 = x
    try:
        f1 = lens_focal_length(R1_1, R2_1)
        f2 = lens_focal_length(R1_2, R2_2)
        f_sys = system_focal_length(f1, f2, L)
        return spot_radius(f_sys)
    except:
        return 1e6  # Return a large number if division by zero or invalid geometry

# Initial guess for radii (in mm)
x0 = [30, -30, 30, -30]

# Bounds to avoid very small radii
bounds = [(-100, -10), (10, 100), (-100, -10), (10, 100)]

result = minimize(objective, x0, bounds=bounds)

# Extract result
R1_1_opt, R2_1_opt, R1_2_opt, R2_2_opt = result.x
f1_opt = lens_focal_length(R1_1_opt, R2_1_opt)
f2_opt = lens_focal_length(R1_2_opt, R2_2_opt)
f_sys_opt = system_focal_length(f1_opt, f2_opt, L)
spot_opt = spot_radius(f_sys_opt)

print(f"Optimized Curvatures: {result.x}")
print(f"Optimized Focal Lengths: f1 = {f1_opt:.2f} mm, f2 = {f2_opt:.2f} mm")
print(f"System Focal Length: {f_sys_opt:.2f} mm")
print(f"RMS Spot Radius Penalty: {spot_opt:.4f}")

from mpl_toolkits.mplot3d import Axes3D

# Grid scan for visualization
R_range = np.linspace(10, 100, 50)
f_sys_map = np.zeros((len(R_range), len(R_range)))

for i, R1 in enumerate(R_range):
    for j, R2 in enumerate(R_range):
        f1 = lens_focal_length(R1, -R1)
        f2 = lens_focal_length(R2, -R2)
        try:
            f_sys = system_focal_length(f1, f2, L)
        except:
            f_sys = np.nan
        f_sys_map[i, j] = f_sys

R1_grid, R2_grid = np.meshgrid(R_range, R_range)

fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(R1_grid, R2_grid, f_sys_map, cmap='viridis', edgecolor='none')
ax.set_xlabel("R1 (mm)")
ax.set_ylabel("R2 (mm)")
ax.set_zlabel("System Focal Length (mm)")
ax.set_title("System Focal Length vs. Lens Curvatures")
plt.colorbar(surf, shrink=0.5, aspect=10)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Rastrigin function
def rastrigin(x):
    return 10 * len(x) + sum([(xi**2 - 10 * np.cos(2 * np.pi * xi)) for xi in x])

# PSO implementation
def pso(n_particles, dimensions, bounds, max_iter, w=0.5, c1=1.5, c2=1.5):
    # Initialize particles and velocities
    particles = np.random.uniform(bounds[0], bounds[1], (n_particles, dimensions))
    velocities = np.random.uniform(-1, 1, (n_particles, dimensions))
    p_best = particles.copy()
    p_best_scores = np.array([rastrigin(p) for p in p_best])
    g_best = p_best[np.argmin(p_best_scores)]
    g_best_score = min(p_best_scores)
    
    # Store history for visualization
    history = [particles.copy()]
    
    # Main PSO loop
    for _ in range(max_iter):
        # Update velocities
        r1, r2 = np.random.random((2, n_particles, 1))
        velocities = (w * velocities +
                     c1 * r1 * (p_best - particles) +
                     c2 * r2 * (g_best - particles))
        
        # Update positions
        particles += velocities
        particles = np.clip(particles, bounds[0], bounds[1])
        
        # Evaluate fitness
        scores = np.array([rastrigin(p) for p in particles])
        
        # Update personal and global best
        improved = scores < p_best_scores
        p_best[improved] = particles[improved]
        p_best_scores[improved] = scores[improved]
        if min(scores) < g_best_score:
            g_best = particles[np.argmin(scores)]
            g_best_score = min(scores)
        
        history.append(particles.copy())
    
    return g_best, g_best_score, history

# Parameters
n_particles = 30
dimensions = 2
bounds = (-5.12, 5.12)
max_iter = 100

# Run PSO
best_position, best_score, history = pso(n_particles, dimensions, bounds, max_iter)

# Visualization
# Contour plot
x = np.linspace(bounds[0], bounds[1], 100)
y = np.linspace(bounds[0], bounds[1], 100)
X, Y = np.meshgrid(x, y)
Z = np.array([[rastrigin([x_, y_]) for x_, y_ in zip(x_row, y_row)] for x_row, y_row in zip(X, Y)])

fig = plt.figure(figsize=(12, 5))

# Contour plot with particle movement
ax1 = fig.add_subplot(121)
ax1.contour(X, Y, Z, levels=20)
ax1.scatter(history[0][:, 0], history[0][:, 1], c='blue', label='Initial Particles', alpha=0.5)
ax1.scatter(history[-1][:, 0], history[-1][:, 1], c='red', label='Final Particles', alpha=0.5)
ax1.scatter(best_position[0], best_position[1], c='green', marker='*', s=200, label='Global Best')
ax1.set_xlabel('X')
ax1.set_ylabel('Y')
ax1.set_title('Particle Movement')
ax1.legend()

# 3D Surface plot
ax2 = fig.add_subplot(122, projection='3d')
ax2.plot_surface(X, Y, Z, cmap='viridis', alpha=0.8)
ax2.scatter(best_position[0], best_position[1], best_score, c='red', s=100, label='Global Best')
ax2.set_xlabel('X')
ax2.set_ylabel('Y')
ax2.set_zlabel('f(x, y)')
ax2.set_title('Rastrigin Function')
ax2.legend()

plt.tight_layout()
plt.show()

print(f"Best position: {best_position}")
print(f"Best score: {best_score}")

import numpy as np
import matplotlib.pyplot as plt

# Constants
f_target = 1000  # 1 kHz
C_E12_series = np.array([1.0, 1.2, 1.5, 1.8, 2.2, 2.7, 3.3, 3.9, 4.7, 5.6, 6.8, 8.2]) * 1e-9  # nF to F
multipliers = [1, 10, 100]  # extend to 1nF to 100nF
C_values = np.array([c * m for m in multipliers for c in C_E12_series])

# Calculate R for each C
R_values = 1 / (2 * np.pi * f_target * C_values)

# Filter valid R values within range
R_min, R_max = 1e3, 1e6  # ohms
valid_indices = np.where((R_values >= R_min) & (R_values <= R_max))[0]

# Select the optimal design
optimal_index = valid_indices[np.argmax(R_values[valid_indices])]
optimal_C = C_values[optimal_index]
optimal_R = R_values[optimal_index]

# Display result
print(f"Optimal C: {optimal_C*1e9:.2f} nF")
print(f"Optimal R: {optimal_R/1e3:.2f} kΩ")
print(f"Cutoff frequency: {1 / (2 * np.pi * optimal_R * optimal_C):.2f} Hz")

# Plotting
plt.figure(figsize=(10, 5))
plt.plot(C_values * 1e9, R_values / 1e3, marker='o', label='R vs C')
plt.axhline(R_max / 1e3, color='red', linestyle='--', label='R max limit')
plt.axhline(R_min / 1e3, color='green', linestyle='--', label='R min limit')
plt.scatter(optimal_C * 1e9, optimal_R / 1e3, color='black', label='Optimal Design', zorder=5)
plt.xlabel('Capacitance (nF)')
plt.ylabel('Resistance (kΩ)')
plt.title('RC Filter Design Optimization')
plt.grid(True)
plt.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import torch
import torch.nn as nn
import torch.optim as optim

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

# Generate synthetic data
x = np.linspace(0, 1, 100)
y_true = np.sin(2 * np.pi * x)
y_noisy = y_true + 0.1 * np.random.randn(*x.shape)

# Convert to PyTorch tensors
x_tensor = torch.tensor(x, dtype=torch.float32).unsqueeze(1)
y_tensor = torch.tensor(y_noisy, dtype=torch.float32).unsqueeze(1)

# Define a simple feedforward neural network
class Net(nn.Module):
    def __init__(self):
        super(Net, self).__init__()
        self.net = nn.Sequential(
            nn.Linear(1, 32),
            nn.ReLU(),
            nn.Linear(32, 32),
            nn.ReLU(),
            nn.Linear(32, 1)
        )

    def forward(self, x):
        return self.net(x)

model = Net()

# Define loss and optimizer with momentum
criterion = nn.MSELoss()
optimizer = optim.SGD(model.parameters(), lr=0.1, momentum=0.9)

# Training loop
epochs = 500
losses = []

for epoch in range(epochs):
    model.train()
    optimizer.zero_grad()
    outputs = model(x_tensor)
    loss = criterion(outputs, y_tensor)
    loss.backward()
    optimizer.step()
    losses.append(loss.item())

# Predictions
model.eval()
with torch.no_grad():
    predictions = model(x_tensor).numpy()

# Plot loss curve
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.plot(losses)
plt.title("Training Loss")
plt.xlabel("Epoch")
plt.ylabel("MSE Loss")
plt.grid(True)

# Plot predictions vs true function
plt.subplot(1, 2, 2)
plt.scatter(x, y_noisy, label="Noisy Data", alpha=0.6)
plt.plot(x, y_true, label="True Function", linestyle="dashed")
plt.plot(x, predictions, label="NN Prediction", color="red")
plt.title("Function Approximation")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

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

# Define power plant characteristics
plants = {
    'Coal': {
        'capacity_min': 200,  # MW
        'capacity_max': 800,  # MW
        'marginal_cost': 25,  # $/MWh
        'fixed_cost': 5000,   # $/h when operating
        'ramp_rate': 50,      # MW/h maximum change
        'efficiency': 0.35    # 35% efficiency
    },
    'Gas': {
        'capacity_min': 100,  # MW
        'capacity_max': 500,  # MW
        'marginal_cost': 45,  # $/MWh
        'fixed_cost': 2000,   # $/h when operating
        'ramp_rate': 200,     # MW/h maximum change
        'efficiency': 0.45    # 45% efficiency
    },
    'Peaker': {
        'capacity_min': 0,    # MW (can be completely shut down)
        'capacity_max': 200,  # MW
        'marginal_cost': 80,  # $/MWh
        'fixed_cost': 500,    # $/h when operating
        'ramp_rate': 150,     # MW/h maximum change
        'efficiency': 0.38    # 38% efficiency
    }
}

# Generate realistic 24-hour demand profile (typical summer day)
hours = np.arange(24)
base_demand = 800  # MW base load
peak_demand = 1200  # MW peak demand

# Create demand curve with morning and evening peaks
demand_profile = base_demand + (peak_demand - base_demand) * (
    0.3 * np.sin(2 * np.pi * (hours - 6) / 24) ** 2 +  # Morning peak
    0.7 * np.sin(2 * np.pi * (hours - 18) / 24) ** 2    # Evening peak
)

# Add some realistic noise
np.random.seed(42)
demand_profile += np.random.normal(0, 20, 24)
demand_profile = np.maximum(demand_profile, 600)  # Minimum demand floor

print("24-Hour Electricity Demand Profile:")
for i, demand in enumerate(demand_profile):
    print(f"Hour {i:2d}: {demand:6.1f} MW")

# Cost function for each plant
def plant_cost(power, plant_name):
    """Calculate total cost for a plant given power output"""
    if power <= 0:
        return 0
    
    plant = plants[plant_name]
    # Quadratic cost function: Fixed cost + marginal cost * power + quadratic term
    variable_cost = plant['marginal_cost'] * power + 0.01 * power**2
    fixed_cost = plant['fixed_cost'] if power > 0 else 0
    
    return fixed_cost + variable_cost

# Objective function for optimization
def total_cost(x):
    """Calculate total system cost for all plants across all hours"""
    # Reshape optimization variable: x = [coal_powers, gas_powers, peaker_powers]
    n_hours = 24
    coal_powers = x[0:n_hours]
    gas_powers = x[n_hours:2*n_hours]
    peaker_powers = x[2*n_hours:3*n_hours]
    
    total = 0
    for t in range(n_hours):
        total += plant_cost(coal_powers[t], 'Coal')
        total += plant_cost(gas_powers[t], 'Gas')
        total += plant_cost(peaker_powers[t], 'Peaker')
    
    return total

# Constraint functions
def power_balance_constraint(x):
    """Ensure power generation meets demand at each hour"""
    n_hours = 24
    coal_powers = x[0:n_hours]
    gas_powers = x[n_hours:2*n_hours]
    peaker_powers = x[2*n_hours:3*n_hours]
    
    # Return difference between generation and demand (should be zero)
    return coal_powers + gas_powers + peaker_powers - demand_profile

def ramp_rate_constraints(x):
    """Ensure power plants don't change output too quickly"""
    n_hours = 24
    constraints = []
    
    # Check ramp rates for each plant type
    plant_names = ['Coal', 'Gas', 'Peaker']
    for i, plant_name in enumerate(plant_names):
        powers = x[i*n_hours:(i+1)*n_hours]
        ramp_limit = plants[plant_name]['ramp_rate']
        
        for t in range(1, n_hours):
            # Constraint: |P(t) - P(t-1)| <= ramp_rate
            # Implemented as two constraints: P(t) - P(t-1) <= ramp_rate and P(t-1) - P(t) <= ramp_rate
            constraints.append(ramp_limit - (powers[t] - powers[t-1]))
            constraints.append(ramp_limit - (powers[t-1] - powers[t]))
    
    return np.array(constraints)

# Set up bounds for each variable (capacity constraints)
bounds = []
for plant_name in ['Coal', 'Gas', 'Peaker']:
    for _ in range(24):  # 24 hours
        bounds.append((plants[plant_name]['capacity_min'], 
                      plants[plant_name]['capacity_max']))

# Initial guess: distribute demand proportionally
initial_guess = []
total_capacity = sum(plants[p]['capacity_max'] for p in plants.keys())
for plant_name in ['Coal', 'Gas', 'Peaker']:
    plant_share = plants[plant_name]['capacity_max'] / total_capacity
    for demand in demand_profile:
        initial_power = demand * plant_share
        # Ensure initial guess is within bounds
        initial_power = max(plants[plant_name]['capacity_min'], 
                          min(initial_power, plants[plant_name]['capacity_max']))
        initial_guess.append(initial_power)

# Define constraints for the optimizer
constraints = [
    {'type': 'eq', 'fun': power_balance_constraint},
    {'type': 'ineq', 'fun': ramp_rate_constraints}
]

print("\nSolving economic dispatch optimization...")
print("This may take a moment...")

# Solve the optimization problem
result = minimize(
    total_cost,
    initial_guess,
    method='SLSQP',
    bounds=bounds,
    constraints=constraints,
    options={'maxiter': 1000, 'ftol': 1e-6}
)

if not result.success:
    print(f"Optimization failed: {result.message}")
else:
    print(f"Optimization successful!")
    print(f"Total daily cost: ${result.fun:,.2f}")

# Extract results
optimal_solution = result.x
n_hours = 24
coal_schedule = optimal_solution[0:n_hours]
gas_schedule = optimal_solution[n_hours:2*n_hours]
peaker_schedule = optimal_solution[2*n_hours:3*n_hours]

# Create results DataFrame
results_df = pd.DataFrame({
    'Hour': hours,
    'Demand_MW': demand_profile,
    'Coal_MW': coal_schedule,
    'Gas_MW': gas_schedule,
    'Peaker_MW': peaker_schedule,
    'Total_Generation_MW': coal_schedule + gas_schedule + peaker_schedule
})

# Calculate hourly costs
hourly_costs = []
for t in range(24):
    hour_cost = (plant_cost(coal_schedule[t], 'Coal') + 
                 plant_cost(gas_schedule[t], 'Gas') + 
                 plant_cost(peaker_schedule[t], 'Peaker'))
    hourly_costs.append(hour_cost)

results_df['Hourly_Cost_$'] = hourly_costs

print("\nOptimal Generation Schedule:")
print(results_df.round(1))

# Calculate summary statistics
total_generation = {
    'Coal': np.sum(coal_schedule),
    'Gas': np.sum(gas_schedule),
    'Peaker': np.sum(peaker_schedule)
}

print(f"\nDaily Generation Summary:")
for plant, total in total_generation.items():
    percentage = 100 * total / np.sum(demand_profile)
    avg_cost = plants[plant]['marginal_cost']
    print(f"{plant:7s}: {total:6.1f} MWh ({percentage:4.1f}%) - Avg Cost: ${avg_cost}/MWh")

print(f"\nTotal Daily Demand: {np.sum(demand_profile):6.1f} MWh")
print(f"Total Daily Cost: ${result.fun:,.2f}")
print(f"Average Cost per MWh: ${result.fun/np.sum(demand_profile):.2f}")

# Plotting results
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))

# Plot 1: Generation Schedule
ax1.plot(hours, demand_profile, 'k-', linewidth=3, label='Demand', marker='o')
ax1.plot(hours, coal_schedule, 'brown', linewidth=2, label='Coal', marker='s')
ax1.plot(hours, gas_schedule, 'blue', linewidth=2, label='Natural Gas', marker='^')
ax1.plot(hours, peaker_schedule, 'red', linewidth=2, label='Peaker', marker='d')

ax1.set_xlabel('Hour of Day')
ax1.set_ylabel('Power (MW)')
ax1.set_title('Optimal Generation Schedule vs Demand')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_xticks(range(0, 24, 2))

# Plot 2: Stacked Generation
ax2.fill_between(hours, 0, coal_schedule, alpha=0.7, color='brown', label='Coal')
ax2.fill_between(hours, coal_schedule, coal_schedule + gas_schedule, 
                alpha=0.7, color='blue', label='Natural Gas')
ax2.fill_between(hours, coal_schedule + gas_schedule, 
                coal_schedule + gas_schedule + peaker_schedule,
                alpha=0.7, color='red', label='Peaker')
ax2.plot(hours, demand_profile, 'k-', linewidth=2, label='Demand')

ax2.set_xlabel('Hour of Day')
ax2.set_ylabel('Power (MW)')
ax2.set_title('Generation Mix (Stacked)')
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.set_xticks(range(0, 24, 2))

# Plot 3: Hourly Costs
ax3.bar(hours, hourly_costs, alpha=0.7, color='green')
ax3.set_xlabel('Hour of Day')
ax3.set_ylabel('Hourly Cost ($)')
ax3.set_title('Hourly Generation Costs')
ax3.grid(True, alpha=0.3)
ax3.set_xticks(range(0, 24, 2))

# Plot 4: Generation Mix Pie Chart
total_gen = [total_generation[plant] for plant in ['Coal', 'Gas', 'Peaker']]
colors = ['brown', 'blue', 'red']
ax4.pie(total_gen, labels=['Coal', 'Natural Gas', 'Peaker'], 
        colors=colors, autopct='%1.1f%%', startangle=90)
ax4.set_title('Daily Generation Mix')

plt.tight_layout()
plt.show()

# Additional analysis: Marginal costs
print(f"\nMarginal Cost Analysis:")
print(f"Coal marginal cost: ${plants['Coal']['marginal_cost']}/MWh")
print(f"Gas marginal cost: ${plants['Gas']['marginal_cost']}/MWh") 
print(f"Peaker marginal cost: ${plants['Peaker']['marginal_cost']}/MWh")

# Calculate capacity factors
print(f"\nCapacity Factors:")
for plant_name in ['Coal', 'Gas', 'Peaker']:
    if plant_name == 'Coal':
        generation = coal_schedule
    elif plant_name == 'Gas':
        generation = gas_schedule
    else:
        generation = peaker_schedule
    
    max_capacity = plants[plant_name]['capacity_max']
    avg_generation = np.mean(generation)
    capacity_factor = avg_generation / max_capacity * 100
    print(f"{plant_name:7s}: {capacity_factor:5.1f}% (Avg: {avg_generation:5.1f} MW)")