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)")

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)