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

# Set up the problem parameters
class ResourceDiplomacyOptimizer:
    def __init__(self):
        # Production capacity (barrels per day)
        self.total_capacity = 2_000_000
        
        # Trading partner characteristics
        self.partners = {
            'Country_A': {'price': 85.0, 'strategic_weight': 0.8, 'min_commitment': 300_000},
            'Country_B': {'price': 82.0, 'strategic_weight': 0.6, 'min_commitment': 200_000},
            'Country_C': {'price': 78.0, 'strategic_weight': 0.9, 'min_commitment': 400_000}
        }
        
        # Weights for multi-objective optimization
        self.economic_weight = 0.7  # Weight for economic returns
        self.diplomatic_weight = 0.3  # Weight for diplomatic benefits
        
    def objective_function(self, x):
        """
        Multi-objective function combining economic and diplomatic benefits
        x = [allocation_A, allocation_B, allocation_C]
        """
        allocation_A, allocation_B, allocation_C = x
        
        # Economic benefit calculation
        economic_benefit = (
            allocation_A * self.partners['Country_A']['price'] +
            allocation_B * self.partners['Country_B']['price'] +
            allocation_C * self.partners['Country_C']['price']
        )
        
        # Diplomatic benefit calculation (normalized strategic importance)
        diplomatic_benefit = (
            allocation_A * self.partners['Country_A']['strategic_weight'] +
            allocation_B * self.partners['Country_B']['strategic_weight'] +
            allocation_C * self.partners['Country_C']['strategic_weight']
        )
        
        # Combined objective (negative because we minimize)
        total_benefit = (
            self.economic_weight * economic_benefit +
            self.diplomatic_weight * diplomatic_benefit * 100_000  # Scale diplomatic benefits
        )
        
        return -total_benefit  # Negative for minimization
    
    def constraints(self):
        """Define constraints for the optimization problem"""
        constraints = []
        
        # Total capacity constraint
        constraints.append({
            'type': 'eq',
            'fun': lambda x: self.total_capacity - sum(x)
        })
        
        # Minimum commitment constraints
        constraints.append({
            'type': 'ineq',
            'fun': lambda x: x[0] - self.partners['Country_A']['min_commitment']
        })
        constraints.append({
            'type': 'ineq',
            'fun': lambda x: x[1] - self.partners['Country_B']['min_commitment']
        })
        constraints.append({
            'type': 'ineq',
            'fun': lambda x: x[2] - self.partners['Country_C']['min_commitment']
        })
        
        return constraints
    
    def solve_optimization(self):
        """Solve the resource allocation optimization problem"""
        # Initial guess (equal distribution after minimum commitments)
        remaining_capacity = self.total_capacity - sum(
            partner['min_commitment'] for partner in self.partners.values()
        )
        
        x0 = [
            self.partners['Country_A']['min_commitment'] + remaining_capacity/3,
            self.partners['Country_B']['min_commitment'] + remaining_capacity/3,
            self.partners['Country_C']['min_commitment'] + remaining_capacity/3
        ]
        
        # Bounds (non-negative allocations, maximum is total capacity)
        bounds = [(0, self.total_capacity) for _ in range(3)]
        
        # Solve optimization
        result = minimize(
            self.objective_function,
            x0,
            method='SLSQP',
            bounds=bounds,
            constraints=self.constraints()
        )
        
        return result
    
    def analyze_sensitivity(self, price_variations=None):
        """Analyze sensitivity to price changes"""
        if price_variations is None:
            price_variations = np.linspace(0.8, 1.2, 10)  # ±20% price variation
        
        results = []
        base_prices = [self.partners[country]['price'] for country in self.partners.keys()]
        
        for variation in price_variations:
            # Temporarily modify prices
            for country in self.partners.keys():
                self.partners[country]['price'] *= variation
            
            # Solve optimization
            result = self.solve_optimization()
            
            if result.success:
                allocation = result.x
                total_revenue = -result.fun
                results.append({
                    'price_multiplier': variation,
                    'country_A_allocation': allocation[0],
                    'country_B_allocation': allocation[1],
                    'country_C_allocation': allocation[2],
                    'total_benefit': total_revenue
                })
            
            # Restore original prices
            for i, country in enumerate(self.partners.keys()):
                self.partners[country]['price'] = base_prices[i]
        
        return pd.DataFrame(results)

# Initialize optimizer
optimizer = ResourceDiplomacyOptimizer()

# Solve the optimization problem
print("=== Resource Diplomacy Optimization Results ===")
result = optimizer.solve_optimization()

if result.success:
    optimal_allocation = result.x
    print(f"\nOptimization Status: {result.message}")
    print(f"\nOptimal Resource Allocation:")
    print(f"Country A: {optimal_allocation[0]:,.0f} barrels/day ({optimal_allocation[0]/optimizer.total_capacity*100:.1f}%)")
    print(f"Country B: {optimal_allocation[1]:,.0f} barrels/day ({optimal_allocation[1]/optimizer.total_capacity*100:.1f}%)")
    print(f"Country C: {optimal_allocation[2]:,.0f} barrels/day ({optimal_allocation[2]/optimizer.total_capacity*100:.1f}%)")
    print(f"\nTotal Allocation: {sum(optimal_allocation):,.0f} barrels/day")
    
    # Calculate economic and diplomatic benefits
    economic_benefit = sum(optimal_allocation[i] * list(optimizer.partners.values())[i]['price'] 
                          for i in range(3))
    diplomatic_benefit = sum(optimal_allocation[i] * list(optimizer.partners.values())[i]['strategic_weight'] 
                           for i in range(3))
    
    print(f"\nDaily Economic Benefit: ${economic_benefit:,.0f}")
    print(f"Diplomatic Benefit Index: {diplomatic_benefit:,.0f}")
    
else:
    print(f"Optimization failed: {result.message}")

# Perform sensitivity analysis
print("\n=== Conducting Sensitivity Analysis ===")
sensitivity_data = optimizer.analyze_sensitivity()

# Create comprehensive visualizations
plt.style.use('default')
fig = plt.figure(figsize=(20, 15))

# 1. Optimal Allocation Pie Chart
ax1 = plt.subplot(2, 3, 1)
countries = ['Country A', 'Country B', 'Country C']
colors = ['#FF6B6B', '#4ECDC4', '#45B7D1']
wedges, texts, autotexts = ax1.pie(optimal_allocation, labels=countries, autopct='%1.1f%%', 
                                   colors=colors, startangle=90, explode=(0.05, 0.05, 0.05))
ax1.set_title('Optimal Resource Allocation\n(Barrels per Day)', fontsize=14, fontweight='bold')

# 2. Economic vs Diplomatic Trade-off
ax2 = plt.subplot(2, 3, 2)
prices = [optimizer.partners[country]['price'] for country in optimizer.partners.keys()]
strategic_weights = [optimizer.partners[country]['strategic_weight'] for country in optimizer.partners.keys()]
scatter = ax2.scatter(prices, strategic_weights, s=[alloc/5000 for alloc in optimal_allocation], 
                     c=colors, alpha=0.7, edgecolors='black', linewidth=2)
ax2.set_xlabel('Oil Price ($/barrel)', fontsize=12)
ax2.set_ylabel('Strategic Importance', fontsize=12)
ax2.set_title('Price vs Strategic Importance\n(Bubble size = Allocation)', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

for i, country in enumerate(countries):
    ax2.annotate(country, (prices[i], strategic_weights[i]), 
                xytext=(5, 5), textcoords='offset points', fontsize=10)

# 3. Sensitivity Analysis - Allocation Changes
ax3 = plt.subplot(2, 3, 3)
ax3.plot(sensitivity_data['price_multiplier'], sensitivity_data['country_A_allocation']/1000, 
         'o-', label='Country A', color=colors[0], linewidth=2, markersize=6)
ax3.plot(sensitivity_data['price_multiplier'], sensitivity_data['country_B_allocation']/1000, 
         's-', label='Country B', color=colors[1], linewidth=2, markersize=6)
ax3.plot(sensitivity_data['price_multiplier'], sensitivity_data['country_C_allocation']/1000, 
         '^-', label='Country C', color=colors[2], linewidth=2, markersize=6)
ax3.set_xlabel('Price Multiplier', fontsize=12)
ax3.set_ylabel('Allocation (Thousands of barrels/day)', fontsize=12)
ax3.set_title('Sensitivity to Price Changes', fontsize=14, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)

# 4. Total Benefit vs Price Multiplier
ax4 = plt.subplot(2, 3, 4)
ax4.plot(sensitivity_data['price_multiplier'], sensitivity_data['total_benefit']/1_000_000, 
         'o-', color='green', linewidth=3, markersize=8)
ax4.set_xlabel('Price Multiplier', fontsize=12)
ax4.set_ylabel('Total Benefit (Millions $)', fontsize=12)
ax4.set_title('Total Economic Benefit\nvs Price Variation', fontsize=14, fontweight='bold')
ax4.grid(True, alpha=0.3)

# 5. Constraint Analysis
ax5 = plt.subplot(2, 3, 5)
min_commitments = [optimizer.partners[country]['min_commitment'] for country in optimizer.partners.keys()]
actual_allocations = optimal_allocation

x_pos = np.arange(len(countries))
width = 0.35

bars1 = ax5.bar(x_pos - width/2, [x/1000 for x in min_commitments], width, 
                label='Minimum Commitments', color='lightcoral', alpha=0.7)
bars2 = ax5.bar(x_pos + width/2, [x/1000 for x in actual_allocations], width, 
                label='Optimal Allocation', color='lightblue', alpha=0.7)

ax5.set_xlabel('Countries', fontsize=12)
ax5.set_ylabel('Allocation (Thousands of barrels/day)', fontsize=12)
ax5.set_title('Minimum Commitments vs\nOptimal Allocation', fontsize=14, fontweight='bold')
ax5.set_xticks(x_pos)
ax5.set_xticklabels(countries)
ax5.legend()
ax5.grid(True, alpha=0.3, axis='y')

# Add value labels on bars
for bar in bars1:
    height = bar.get_height()
    ax5.text(bar.get_x() + bar.get_width()/2., height + 10,
             f'{height:.0f}', ha='center', va='bottom', fontsize=9)
for bar in bars2:
    height = bar.get_height()
    ax5.text(bar.get_x() + bar.get_width()/2., height + 10,
             f'{height:.0f}', ha='center', va='bottom', fontsize=9)

# 6. Revenue Breakdown
ax6 = plt.subplot(2, 3, 6)
revenues = [optimal_allocation[i] * prices[i] for i in range(3)]
colors_revenue = ['#FF9999', '#66B2FF', '#99FF99']

bars = ax6.bar(countries, [r/1_000_000 for r in revenues], color=colors_revenue, alpha=0.8, edgecolor='black')
ax6.set_ylabel('Daily Revenue (Millions $)', fontsize=12)
ax6.set_title('Daily Revenue by Trading Partner', fontsize=14, fontweight='bold')
ax6.grid(True, alpha=0.3, axis='y')

# Add value labels on bars
for i, bar in enumerate(bars):
    height = bar.get_height()
    ax6.text(bar.get_x() + bar.get_width()/2., height + 0.5,
             f'${height:.1f}M', ha='center', va='bottom', fontsize=11, fontweight='bold')

plt.tight_layout()
plt.show()

# Summary statistics
print("\n=== Summary Statistics ===")
total_revenue = sum(revenues)
print(f"Total Daily Revenue: ${total_revenue:,.0f}")
print(f"Annual Revenue Projection: ${total_revenue * 365:,.0f}")

utilization_rates = [optimal_allocation[i] / optimizer.partners[list(optimizer.partners.keys())[i]]['min_commitment'] 
                    for i in range(3)]
print(f"\nUtilization above minimum commitments:")
for i, country in enumerate(countries):
    excess = optimal_allocation[i] - min_commitments[i]
    print(f"{country}: +{excess:,.0f} barrels/day ({(utilization_rates[i]-1)*100:.1f}% above minimum)")

# Pipeline and Maritime Transport Route Optimization
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 FancyBboxPatch
import warnings
warnings.filterwarnings('ignore')

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

class TransportOptimizer:
    def __init__(self, platforms, refineries):
        self.platforms = platforms
        self.refineries = refineries
        self.m = len(platforms)  # number of platforms
        self.n = len(refineries)  # number of refineries
        
        # Generate supply and demand data
        self.supply = np.array([p['capacity'] for p in platforms])
        self.demand = np.array([r['demand'] for r in refineries])
        
        # Generate cost matrices based on distances and transport modes
        self.pipeline_costs = self._generate_pipeline_costs()
        self.ship_costs = self._generate_ship_costs()
        
    def _calculate_distance(self, p1, p2):
        """Calculate Euclidean distance between two points"""
        return np.sqrt((p1[0] - p2[0])**2 + (p1[1] - p2[1])**2)
    
    def _generate_pipeline_costs(self):
        """Generate pipeline transportation costs based on distance"""
        costs = np.zeros((self.m, self.n))
        base_pipeline_cost = 15  # $/barrel per unit distance
        
        for i, platform in enumerate(self.platforms):
            for j, refinery in enumerate(self.refineries):
                distance = self._calculate_distance(platform['location'], refinery['location'])
                # Pipeline cost increases with distance, but has economies of scale
                costs[i, j] = base_pipeline_cost * distance * 0.8
                
                # Some routes may not have pipeline infrastructure
                if distance > 300:  # Long distances make pipelines impractical
                    costs[i, j] = float('inf')
                    
        return costs
    
    def _generate_ship_costs(self):
        """Generate maritime transportation costs"""
        costs = np.zeros((self.m, self.n))
        base_ship_cost = 8  # $/barrel per unit distance
        
        for i, platform in enumerate(self.platforms):
            for j, refinery in enumerate(self.refineries):
                distance = self._calculate_distance(platform['location'], refinery['location'])
                # Ship costs are more linear with distance but have fixed loading costs
                fixed_cost = 500  # Fixed loading/unloading cost per route
                variable_cost = base_ship_cost * distance
                costs[i, j] = fixed_cost / 1000 + variable_cost  # Normalize fixed cost
                
        return costs
    
    def optimize_routes(self):
        """Solve the transportation optimization problem using linear programming"""
        # Create cost vector for the linear program
        # Variables: [x_11, x_12, ..., x_mn, y_11, y_12, ..., y_mn]
        pipeline_costs_flat = self.pipeline_costs.flatten()
        ship_costs_flat = self.ship_costs.flatten()
        c = np.concatenate([pipeline_costs_flat, ship_costs_flat])
        
        # Replace infinite costs with a large number
        c[c == float('inf')] = 1e6
        
        # Constraint matrices
        A_eq = []
        b_eq = []
        
        # Supply constraints: sum over j of (x_ij + y_ij) <= S_i
        for i in range(self.m):
            constraint = np.zeros(2 * self.m * self.n)
            for j in range(self.n):
                # Pipeline variable x_ij
                constraint[i * self.n + j] = 1
                # Ship variable y_ij
                constraint[self.m * self.n + i * self.n + j] = 1
            A_eq.append(constraint)
            b_eq.append(self.supply[i])
        
        # Demand constraints: sum over i of (x_ij + y_ij) >= D_j
        # Convert to equality by adding slack variables would be complex,
        # so we'll use inequality constraints with bounds
        
        A_ub = []
        b_ub = []
        
        # Convert supply constraints to inequality (≤)
        for i in range(self.m):
            constraint = np.zeros(2 * self.m * self.n)
            for j in range(self.n):
                constraint[i * self.n + j] = 1
                constraint[self.m * self.n + i * self.n + j] = 1
            A_ub.append(constraint)
            b_ub.append(self.supply[i])
        
        # Demand constraints: -sum over i of (x_ij + y_ij) ≤ -D_j
        for j in range(self.n):
            constraint = np.zeros(2 * self.m * self.n)
            for i in range(self.m):
                constraint[i * self.n + j] = -1
                constraint[self.m * self.n + i * self.n + j] = -1
            A_ub.append(constraint)
            b_ub.append(-self.demand[j])
        
        # Bounds: all variables >= 0
        bounds = [(0, None) for _ in range(2 * self.m * self.n)]
        
        # Solve the linear program
        result = linprog(c, A_ub=np.array(A_ub), b_ub=np.array(b_ub), 
                        bounds=bounds, method='highs')
        
        if result.success:
            # Extract solution
            solution = result.x
            pipeline_solution = solution[:self.m * self.n].reshape((self.m, self.n))
            ship_solution = solution[self.m * self.n:].reshape((self.m, self.n))
            
            return {
                'pipeline_flows': pipeline_solution,
                'ship_flows': ship_solution,
                'total_cost': result.fun,
                'success': True
            }
        else:
            return {'success': False, 'message': result.message}

def create_visualization(optimizer, solution):
    """Create comprehensive visualizations of the optimization results"""
    
    fig, axes = plt.subplots(2, 2, figsize=(16, 12))
    fig.suptitle('Pipeline and Maritime Transport Route Optimization Results', 
                 fontsize=16, fontweight='bold')
    
    # 1. Geographic representation of routes
    ax1 = axes[0, 0]
    
    # Plot platforms and refineries
    platform_locs = np.array([p['location'] for p in optimizer.platforms])
    refinery_locs = np.array([r['location'] for r in optimizer.refineries])
    
    ax1.scatter(platform_locs[:, 0], platform_locs[:, 1], 
               c='red', s=200, marker='o', label='Oil Platforms', alpha=0.8)
    ax1.scatter(refinery_locs[:, 0], refinery_locs[:, 1], 
               c='blue', s=200, marker='s', label='Refineries', alpha=0.8)
    
    # Draw routes with thickness proportional to flow
    pipeline_flows = solution['pipeline_flows']
    ship_flows = solution['ship_flows']
    
    max_flow = max(np.max(pipeline_flows), np.max(ship_flows))
    
    for i in range(optimizer.m):
        for j in range(optimizer.n):
            if pipeline_flows[i, j] > 0.1:  # Only show significant flows
                thickness = (pipeline_flows[i, j] / max_flow) * 5
                ax1.plot([platform_locs[i, 0], refinery_locs[j, 0]], 
                        [platform_locs[i, 1], refinery_locs[j, 1]], 
                        'g-', linewidth=thickness, alpha=0.7, label='Pipeline' if i==0 and j==0 else "")
            
            if ship_flows[i, j] > 0.1:
                thickness = (ship_flows[i, j] / max_flow) * 5
                ax1.plot([platform_locs[i, 0], refinery_locs[j, 0]], 
                        [platform_locs[i, 1], refinery_locs[j, 1]], 
                        'b--', linewidth=thickness, alpha=0.7, label='Maritime' if i==0 and j==1 else "")
    
    ax1.set_xlabel('X Coordinate (km)')
    ax1.set_ylabel('Y Coordinate (km)')
    ax1.set_title('Geographic Route Visualization')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. Cost comparison heatmap
    ax2 = axes[0, 1]
    cost_comparison = optimizer.pipeline_costs + optimizer.ship_costs
    # Replace inf values for visualization
    cost_comparison[cost_comparison > 1000] = np.nan
    
    im2 = ax2.imshow(cost_comparison, cmap='YlOrRd', aspect='auto')
    ax2.set_xlabel('Refinery Index')
    ax2.set_ylabel('Platform Index')
    ax2.set_title('Total Transportation Cost Matrix\n($/barrel)')
    
    # Add colorbar
    plt.colorbar(im2, ax=ax2)
    
    # Add text annotations
    for i in range(optimizer.m):
        for j in range(optimizer.n):
            if not np.isnan(cost_comparison[i, j]):
                ax2.text(j, i, f'{cost_comparison[i, j]:.0f}', 
                        ha='center', va='center', fontsize=8)
    
    # 3. Flow distribution
    ax3 = axes[1, 0]
    
    # Create stacked bar chart for each platform-refinery pair
    x_labels = [f'P{i+1}→R{j+1}' for i in range(optimizer.m) for j in range(optimizer.n) 
                if pipeline_flows[i,j] > 0.1 or ship_flows[i,j] > 0.1]
    pipeline_values = [pipeline_flows[i,j] for i in range(optimizer.m) for j in range(optimizer.n) 
                      if pipeline_flows[i,j] > 0.1 or ship_flows[i,j] > 0.1]
    ship_values = [ship_flows[i,j] for i in range(optimizer.m) for j in range(optimizer.n) 
                  if pipeline_flows[i,j] > 0.1 or ship_flows[i,j] > 0.1]
    
    if x_labels:  # Only plot if there are active routes
        x_pos = np.arange(len(x_labels))
        ax3.bar(x_pos, pipeline_values, label='Pipeline', alpha=0.8, color='green')
        ax3.bar(x_pos, ship_values, bottom=pipeline_values, label='Maritime', alpha=0.8, color='blue')
        
        ax3.set_xlabel('Route (Platform → Refinery)')
        ax3.set_ylabel('Flow (1000 barrels/day)')
        ax3.set_title('Transportation Flow by Route and Mode')
        ax3.set_xticks(x_pos)
        ax3.set_xticklabels(x_labels, rotation=45)
        ax3.legend()
    
    # 4. Supply and demand analysis
    ax4 = axes[1, 1]
    
    platform_names = [f'Platform {i+1}' for i in range(optimizer.m)]
    refinery_names = [f'Refinery {j+1}' for j in range(optimizer.n)]
    
    total_supply = optimizer.supply
    total_outflow = np.sum(pipeline_flows + ship_flows, axis=1)
    total_demand = optimizer.demand
    total_inflow = np.sum(pipeline_flows + ship_flows, axis=0)
    
    x_supply = np.arange(len(platform_names))
    x_demand = np.arange(len(refinery_names)) + len(platform_names) + 1
    
    ax4.bar(x_supply, total_supply, alpha=0.7, label='Available Supply', color='orange')
    ax4.bar(x_supply, total_outflow, alpha=0.7, label='Actual Outflow', color='red')
    ax4.bar(x_demand, total_demand, alpha=0.7, label='Required Demand', color='lightblue')
    ax4.bar(x_demand, total_inflow, alpha=0.7, label='Actual Inflow', color='navy')
    
    # Combine labels
    all_labels = platform_names + [''] + refinery_names
    all_x = np.concatenate([x_supply, [len(platform_names)], x_demand])
    
    ax4.set_xticks(all_x)
    ax4.set_xticklabels(all_labels, rotation=45)
    ax4.set_ylabel('Volume (1000 barrels/day)')
    ax4.set_title('Supply and Demand Balance Analysis')
    ax4.legend()
    
    plt.tight_layout()
    plt.show()
    
    return fig

def print_detailed_results(optimizer, solution):
    """Print detailed analysis of the optimization results"""
    
    print("="*80)
    print("TRANSPORTATION OPTIMIZATION RESULTS")
    print("="*80)
    
    pipeline_flows = solution['pipeline_flows']
    ship_flows = solution['ship_flows']
    total_flows = pipeline_flows + ship_flows
    
    print(f"\nTotal Minimum Cost: ${solution['total_cost']:.2f}")
    print(f"Average Cost per Barrel: ${solution['total_cost']/np.sum(total_flows):.2f}")
    
    print("\nFLOW ALLOCATION SUMMARY:")
    print("-" * 50)
    
    total_pipeline = np.sum(pipeline_flows)
    total_ship = np.sum(ship_flows)
    total_volume = total_pipeline + total_ship
    
    print(f"Total Pipeline Transport: {total_pipeline:.1f} (1000 barrels/day)")
    print(f"Total Maritime Transport: {total_ship:.1f} (1000 barrels/day)")
    print(f"Pipeline Share: {(total_pipeline/total_volume)*100:.1f}%")
    print(f"Maritime Share: {(total_ship/total_volume)*100:.1f}%")
    
    print("\nROUTE-SPECIFIC ANALYSIS:")
    print("-" * 50)
    
    for i in range(optimizer.m):
        for j in range(optimizer.n):
            if total_flows[i, j] > 0.1:
                pipeline_cost = optimizer.pipeline_costs[i, j] * pipeline_flows[i, j]
                ship_cost = optimizer.ship_costs[i, j] * ship_flows[i, j]
                total_cost = pipeline_cost + ship_cost
                
                print(f"\nPlatform {i+1} → Refinery {j+1}:")
                print(f"  Pipeline: {pipeline_flows[i, j]:.1f} units @ ${optimizer.pipeline_costs[i, j]:.2f}/unit = ${pipeline_cost:.2f}")
                print(f"  Maritime: {ship_flows[i, j]:.1f} units @ ${optimizer.ship_costs[i, j]:.2f}/unit = ${ship_cost:.2f}")
                print(f"  Total Route Cost: ${total_cost:.2f}")
    
    print("\nCAPACITY UTILIZATION:")
    print("-" * 50)
    
    for i in range(optimizer.m):
        utilization = np.sum(total_flows[i, :]) / optimizer.supply[i] * 100
        print(f"Platform {i+1}: {utilization:.1f}% ({np.sum(total_flows[i, :]):.1f}/{optimizer.supply[i]:.1f})")
    
    print("\nDEMAND SATISFACTION:")
    print("-" * 50)
    
    for j in range(optimizer.n):
        satisfaction = np.sum(total_flows[:, j]) / optimizer.demand[j] * 100
        print(f"Refinery {j+1}: {satisfaction:.1f}% ({np.sum(total_flows[:, j]):.1f}/{optimizer.demand[j]:.1f})")

# Main execution
if __name__ == "__main__":
    # Define oil platforms with locations and capacities
    platforms = [
        {'location': (50, 100), 'capacity': 1500},   # Platform 1
        {'location': (150, 200), 'capacity': 2000},  # Platform 2
        {'location': (250, 150), 'capacity': 1200},  # Platform 3
        {'location': (100, 300), 'capacity': 1800},  # Platform 4
    ]
    
    # Define refineries with locations and demand
    refineries = [
        {'location': (200, 50), 'demand': 1800},     # Refinery 1
        {'location': (300, 250), 'demand': 1500},    # Refinery 2
        {'location': (400, 100), 'demand': 1200},    # Refinery 3
        {'location': (350, 300), 'demand': 1000},    # Refinery 4
    ]
    
    # Create optimizer instance
    optimizer = TransportOptimizer(platforms, refineries)
    
    # Display problem setup
    print("TRANSPORTATION OPTIMIZATION PROBLEM SETUP")
    print("="*60)
    print(f"Number of Oil Platforms: {len(platforms)}")
    print(f"Number of Refineries: {len(refineries)}")
    print(f"Total Supply Capacity: {sum(p['capacity'] for p in platforms):,} (1000 barrels/day)")
    print(f"Total Demand: {sum(r['demand'] for r in refineries):,} (1000 barrels/day)")
    
    # Solve optimization problem
    print("\nSolving optimization problem...")
    solution = optimizer.optimize_routes()
    
    if solution['success']:
        print("✓ Optimization completed successfully!")
        
        # Print detailed results
        print_detailed_results(optimizer, solution)
        
        # Create visualizations
        fig = create_visualization(optimizer, solution)
        
    else:
        print("✗ Optimization failed:", solution['message'])

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 style for better visualizations
plt.style.use('seaborn-v0_8')
sns.set_palette("husl")

print("=== Strategic Resource Allocation Optimization ===\n")

# Problem Definition
print("Problem: Maximize profit from manufacturing three products (A, B, C)")
print("Profit per unit: Product A = $50, Product B = $40, Product C = $60\n")

# Coefficients for the objective function (we negate for minimization)
c = [-50, -40, -60]  # Negative because linprog minimizes

# Inequality constraint matrix (A_ub * x <= b_ub)
A_ub = [
    [2, 3, 1],   # Labor constraint: 2x1 + 3x2 + x3 <= 100
    [1, 2, 3],   # Raw materials: x1 + 2x2 + 3x3 <= 80
    [3, 1, 2]    # Machine time: 3x1 + x2 + 2x3 <= 120
]

b_ub = [100, 80, 120]  # Right-hand side values

# Bounds for variables (all >= 0)
x_bounds = [(0, None), (0, None), (0, None)]

print("Constraints:")
print("Labor hours:     2x₁ + 3x₂ + x₃ ≤ 100")
print("Raw materials:   x₁ + 2x₂ + 3x₃ ≤ 80")
print("Machine time:    3x₁ + x₂ + 2x₃ ≤ 120")
print("Non-negativity:  x₁, x₂, x₃ ≥ 0\n")

# Solve the linear programming problem
result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=x_bounds, method='highs')

print("=== OPTIMIZATION RESULTS ===")
if result.success:
    x1, x2, x3 = result.x
    optimal_profit = -result.fun  # Convert back to positive (maximization)
    
    print(f"Optimal Solution Found!")
    print(f"Product A (x₁): {x1:.2f} units")
    print(f"Product B (x₂): {x2:.2f} units") 
    print(f"Product C (x₃): {x3:.2f} units")
    print(f"Maximum Profit: ${optimal_profit:.2f}")
    
    # Check resource utilization
    print(f"\n=== RESOURCE UTILIZATION ===")
    labor_used = 2*x1 + 3*x2 + x3
    materials_used = x1 + 2*x2 + 3*x3
    machine_used = 3*x1 + x2 + 2*x3
    
    print(f"Labor hours used: {labor_used:.2f} / 100 ({labor_used/100*100:.1f}%)")
    print(f"Raw materials used: {materials_used:.2f} / 80 ({materials_used/80*100:.1f}%)")
    print(f"Machine time used: {machine_used:.2f} / 120 ({machine_used/120*100:.1f}%)")
    
    # Identify binding constraints
    print(f"\n=== BINDING CONSTRAINTS ===")
    tolerances = [abs(labor_used - 100), abs(materials_used - 80), abs(machine_used - 120)]
    constraints = ['Labor', 'Raw Materials', 'Machine Time']
    
    for i, (constraint, tolerance) in enumerate(zip(constraints, tolerances)):
        if tolerance < 1e-6:
            print(f"✓ {constraint} constraint is BINDING (fully utilized)")
        else:
            available = [100, 80, 120][i]
            used = [labor_used, materials_used, machine_used][i]
            slack = available - used
            print(f"  {constraint} has slack of {slack:.2f} units")
    
else:
    print("Optimization failed!")
    print(result.message)

# Sensitivity Analysis
print(f"\n=== SENSITIVITY ANALYSIS ===")
print("Analyzing how profit changes with resource availability...")

# Vary each resource constraint
resource_variations = np.linspace(0.8, 1.2, 21)  # 80% to 120% of original
base_constraints = [100, 80, 120]
resource_names = ['Labor Hours', 'Raw Materials', 'Machine Time']

sensitivity_data = {}

for i, resource_name in enumerate(resource_names):
    profits = []
    for variation in resource_variations:
        # Create modified constraints
        modified_b_ub = base_constraints.copy()
        modified_b_ub[i] = base_constraints[i] * variation
        
        # Solve with modified constraint
        temp_result = linprog(c, A_ub=A_ub, b_ub=modified_b_ub, bounds=x_bounds, method='highs')
        
        if temp_result.success:
            profits.append(-temp_result.fun)
        else:
            profits.append(0)
    
    sensitivity_data[resource_name] = profits

# Create visualizations
fig = plt.figure(figsize=(16, 12))

# 1. Resource Utilization Chart
ax1 = fig.add_subplot(2, 3, 1)
resources = ['Labor', 'Materials', 'Machine']
utilization = [labor_used, materials_used, machine_used]
capacity = [100, 80, 120]
utilization_pct = [u/c*100 for u, c in zip(utilization, capacity)]

bars = ax1.bar(resources, utilization_pct, alpha=0.7, color=['#FF6B6B', '#4ECDC4', '#45B7D1'])
ax1.axhline(y=100, color='red', linestyle='--', alpha=0.7, label='Full Capacity')
ax1.set_ylabel('Utilization (%)')
ax1.set_title('Resource Utilization Analysis')
ax1.set_ylim(0, 120)

# Add value labels on bars
for bar, pct in zip(bars, utilization_pct):
    height = bar.get_height()
    ax1.text(bar.get_x() + bar.get_width()/2., height + 1,
             f'{pct:.1f}%', ha='center', va='bottom', fontweight='bold')

ax1.legend()
ax1.grid(True, alpha=0.3)

# 2. Product Mix Visualization
ax2 = fig.add_subplot(2, 3, 2)
products = ['Product A', 'Product B', 'Product C']
quantities = [x1, x2, x3]
profits_per_product = [q*p for q, p in zip(quantities, [50, 40, 60])]

colors = ['#FF6B6B', '#4ECDC4', '#45B7D1']
bars = ax2.bar(products, quantities, alpha=0.7, color=colors)
ax2.set_ylabel('Units Produced')
ax2.set_title('Optimal Product Mix')

# Add value labels
for bar, qty in zip(bars, quantities):
    height = bar.get_height()
    ax2.text(bar.get_x() + bar.get_width()/2., height + 0.5,
             f'{qty:.1f}', ha='center', va='bottom', fontweight='bold')

ax2.grid(True, alpha=0.3)

# 3. Profit Breakdown
ax3 = fig.add_subplot(2, 3, 3)
wedges, texts, autotexts = ax3.pie(profits_per_product, labels=products, autopct='%1.1f%%',
                                   colors=colors, startangle=90)
ax3.set_title('Profit Contribution by Product')

# Make percentage text bold
for autotext in autotexts:
    autotext.set_color('white')
    autotext.set_fontweight('bold')

# 4. Sensitivity Analysis - Labor
ax4 = fig.add_subplot(2, 3, 4)
labor_percentages = resource_variations * 100
ax4.plot(labor_percentages, sensitivity_data['Labor Hours'], 
         marker='o', linewidth=2, markersize=4, color='#FF6B6B')
ax4.axvline(x=100, color='red', linestyle='--', alpha=0.7, label='Current Level')
ax4.set_xlabel('Labor Hours (% of current)')
ax4.set_ylabel('Maximum Profit ($)')
ax4.set_title('Sensitivity: Labor Hours vs Profit')
ax4.grid(True, alpha=0.3)
ax4.legend()

# 5. Sensitivity Analysis - Raw Materials
ax5 = fig.add_subplot(2, 3, 5)
ax5.plot(resource_variations * 100, sensitivity_data['Raw Materials'],
         marker='s', linewidth=2, markersize=4, color='#4ECDC4')
ax5.axvline(x=100, color='red', linestyle='--', alpha=0.7, label='Current Level')
ax5.set_xlabel('Raw Materials (% of current)')
ax5.set_ylabel('Maximum Profit ($)')
ax5.set_title('Sensitivity: Raw Materials vs Profit')
ax5.grid(True, alpha=0.3)
ax5.legend()

# 6. Sensitivity Analysis - Machine Time
ax6 = fig.add_subplot(2, 3, 6)
ax6.plot(resource_variations * 100, sensitivity_data['Machine Time'],
         marker='^', linewidth=2, markersize=4, color='#45B7D1')
ax6.axvline(x=100, color='red', linestyle='--', alpha=0.7, label='Current Level')
ax6.set_xlabel('Machine Time (% of current)')
ax6.set_ylabel('Maximum Profit ($)')
ax6.set_title('Sensitivity: Machine Time vs Profit')
ax6.grid(True, alpha=0.3)
ax6.legend()

plt.tight_layout()
plt.show()

# Summary table
print(f"\n=== SUMMARY TABLE ===")
summary_df = pd.DataFrame({
    'Product': ['A', 'B', 'C'],
    'Optimal Quantity': [f"{x:.2f}" for x in [x1, x2, x3]],
    'Profit per Unit ($)': [50, 40, 60],
    'Total Profit ($)': [f"{p:.2f}" for p in profits_per_product],
    'Profit Share (%)': [f"{p/sum(profits_per_product)*100:.1f}" for p in profits_per_product]
})

print(summary_df.to_string(index=False))

# Resource efficiency analysis
print(f"\n=== RESOURCE EFFICIENCY ANALYSIS ===")
resource_efficiency_df = pd.DataFrame({
    'Resource': ['Labor Hours', 'Raw Materials', 'Machine Time'],
    'Available': [100, 80, 120],
    'Used': [f"{u:.2f}" for u in [labor_used, materials_used, machine_used]],
    'Utilization (%)': [f"{u:.1f}" for u in utilization_pct],
    'Slack': [f"{s:.2f}" for s in [100-labor_used, 80-materials_used, 120-machine_used]]
})

print(resource_efficiency_df.to_string(index=False))

print(f"\n=== BUSINESS INSIGHTS ===")
print("1. The optimal solution maximizes profit while efficiently using available resources")
print("2. Resource utilization analysis helps identify bottlenecks and expansion opportunities")
print("3. Sensitivity analysis reveals which resources have the highest impact on profitability")
print("4. This model can be extended to include additional constraints like storage, demand limits, etc.")

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

# Set up the optimization problem
print("=" * 60)
print("ENERGY SUPPLY DIVERSIFICATION OPTIMIZATION")
print("=" * 60)

# Define energy sources
sources = ['Coal', 'Natural Gas', 'Nuclear', 'Solar', 'Wind']
n_sources = len(sources)

# Cost per MWh ($/MWh)
costs = np.array([30, 45, 60, 80, 70])

# Emissions per MWh (tons CO2/MWh)
emissions = np.array([0.9, 0.5, 0.02, 0.0, 0.0])

# Capacity limits (MWh)
capacities = np.array([400, 350, 300, 200, 250])

# Total demand (MWh)
demand = 1000

# Maximum allowed emissions (tons CO2)
max_emissions = 400

print(f"Total Demand: {demand} MWh")
print(f"Maximum Emissions Allowed: {max_emissions} tons CO2")
print("\nEnergy Source Details:")
print("-" * 50)

# Create a summary table
source_data = pd.DataFrame({
    'Source': sources,
    'Cost ($/MWh)': costs,
    'Emissions (tons CO2/MWh)': emissions,
    'Capacity (MWh)': capacities
})
print(source_data.to_string(index=False))

# Solve the optimization problem using scipy.optimize.linprog
# Note: linprog minimizes, so we use costs directly

# Objective function coefficients (costs to minimize)
c = costs

# Inequality constraint matrix (emissions constraint)
# emissions * x <= max_emissions
A_ub = emissions.reshape(1, -1)
b_ub = np.array([max_emissions])

# Equality constraint matrix (demand constraint)
# sum(x) = demand
A_eq = np.ones((1, n_sources))
b_eq = np.array([demand])

# Bounds for each variable (0 <= x_i <= capacity_i)
bounds = [(0, cap) for cap in capacities]

print(f"\nSolving optimization problem...")
print(f"Minimize: {' + '.join([f'{c[i]:.0f}*x_{i+1}' for i in range(n_sources)])}")
print(f"Subject to:")
print(f"  Demand: {' + '.join([f'x_{i+1}' for i in range(n_sources)])} = {demand}")
print(f"  Emissions: {' + '.join([f'{emissions[i]:.2f}*x_{i+1}' for i in range(n_sources)])} ≤ {max_emissions}")
print(f"  Capacity: 0 ≤ x_i ≤ Cap_i for all i")

# Solve the linear programming problem
result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, 
                bounds=bounds, method='highs')

if result.success:
    optimal_solution = result.x
    optimal_cost = result.fun
    
    print(f"\nOptimization successful!")
    print(f"Minimum total cost: ${optimal_cost:,.2f}")
    
    # Calculate actual emissions
    actual_emissions = np.sum(emissions * optimal_solution)
    
    print(f"Total emissions: {actual_emissions:.2f} tons CO2")
    print(f"Emission constraint utilization: {actual_emissions/max_emissions*100:.1f}%")
    
    print(f"\nOptimal Energy Mix:")
    print("-" * 40)
    
    results_df = pd.DataFrame({
        'Source': sources,
        'Optimal Output (MWh)': optimal_solution,
        'Percentage (%)': optimal_solution / demand * 100,
        'Cost ($)': optimal_solution * costs,
        'Emissions (tons CO2)': optimal_solution * emissions
    })
    
    # Format the results nicely
    results_df['Optimal Output (MWh)'] = results_df['Optimal Output (MWh)'].round(2)
    results_df['Percentage (%)'] = results_df['Percentage (%)'].round(2)
    results_df['Cost ($)'] = results_df['Cost ($)'].round(2)
    results_df['Emissions (tons CO2)'] = results_df['Emissions (tons CO2)'].round(2)
    
    print(results_df.to_string(index=False))
    
else:
    print("Optimization failed!")
    print(f"Status: {result.message}")

# Sensitivity Analysis: What happens if we tighten emission constraints?
print(f"\n" + "=" * 60)
print("SENSITIVITY ANALYSIS")
print("=" * 60)

emission_limits = np.linspace(200, 600, 11)
costs_sensitivity = []
emission_usage = []

for limit in emission_limits:
    b_ub_temp = np.array([limit])
    result_temp = linprog(c, A_ub=A_ub, b_ub=b_ub_temp, A_eq=A_eq, b_eq=b_eq, 
                         bounds=bounds, method='highs')
    
    if result_temp.success:
        costs_sensitivity.append(result_temp.fun)
        actual_em = np.sum(emissions * result_temp.x)
        emission_usage.append(actual_em)
    else:
        costs_sensitivity.append(np.nan)
        emission_usage.append(np.nan)

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

# 1. Optimal Energy Mix (Pie Chart)
ax1 = plt.subplot(2, 3, 1)
non_zero_mask = optimal_solution > 0.1  # Only show sources with significant contribution
pie_data = optimal_solution[non_zero_mask]
pie_labels = [sources[i] for i in range(n_sources) if non_zero_mask[i]]
colors = ['#8B4513', '#FF6B35', '#4169E1', '#FFD700', '#32CD32']
pie_colors = [colors[i] for i in range(n_sources) if non_zero_mask[i]]

plt.pie(pie_data, labels=pie_labels, autopct='%1.1f%%', startangle=90, colors=pie_colors)
plt.title('Optimal Energy Mix Distribution', fontsize=14, fontweight='bold')

# 2. Cost vs Emissions Trade-off
ax2 = plt.subplot(2, 3, 2)
plt.plot(emission_usage, costs_sensitivity, 'b-o', linewidth=2, markersize=6)
plt.xlabel('Total Emissions (tons CO2)', fontsize=12)
plt.ylabel('Total Cost ($)', fontsize=12)
plt.title('Cost vs Emissions Trade-off Curve', fontsize=14, fontweight='bold')
plt.grid(True, alpha=0.3)
# Highlight our solution
opt_emissions = np.sum(emissions * optimal_solution)
plt.plot(opt_emissions, optimal_cost, 'ro', markersize=10, label=f'Optimal Solution\n({opt_emissions:.1f} tons, ${optimal_cost:,.0f})')
plt.legend()

# 3. Source Utilization vs Capacity
ax3 = plt.subplot(2, 3, 3)
utilization = (optimal_solution / capacities) * 100
bars = plt.bar(sources, utilization, color=colors, alpha=0.7, edgecolor='black')
plt.ylabel('Capacity Utilization (%)', fontsize=12)
plt.title('Capacity Utilization by Source', fontsize=14, fontweight='bold')
plt.xticks(rotation=45)
plt.ylim(0, 100)

# Add value labels on bars
for bar, util in zip(bars, utilization):
    height = bar.get_height()
    plt.text(bar.get_x() + bar.get_width()/2., height + 1,
             f'{util:.1f}%', ha='center', va='bottom', fontweight='bold')

# 4. Cost Breakdown
ax4 = plt.subplot(2, 3, 4)
cost_breakdown = optimal_solution * costs
bars = plt.bar(sources, cost_breakdown, color=colors, alpha=0.7, edgecolor='black')
plt.ylabel('Total Cost ($)', fontsize=12)
plt.title('Cost Breakdown by Source', fontsize=14, fontweight='bold')
plt.xticks(rotation=45)

# Add value labels on bars
for bar, cost in zip(bars, cost_breakdown):
    height = bar.get_height()
    if height > 0:
        plt.text(bar.get_x() + bar.get_width()/2., height + height*0.01,
                 f'${cost:,.0f}', ha='center', va='bottom', fontweight='bold')

# 5. Emissions Breakdown
ax5 = plt.subplot(2, 3, 5)
emission_breakdown = optimal_solution * emissions
bars = plt.bar(sources, emission_breakdown, color=colors, alpha=0.7, edgecolor='black')
plt.ylabel('Total Emissions (tons CO2)', fontsize=12)
plt.title('Emissions Breakdown by Source', fontsize=14, fontweight='bold')
plt.xticks(rotation=45)

# Add value labels on bars
for bar, em in zip(bars, emission_breakdown):
    height = bar.get_height()
    if height > 0:
        plt.text(bar.get_x() + bar.get_width()/2., height + height*0.01,
                 f'{em:.1f}', ha='center', va='bottom', fontweight='bold')

# 6. Comparison: Current vs Alternative Scenarios
ax6 = plt.subplot(2, 3, 6)

# Create comparison scenarios
scenarios = {
    'Optimal': optimal_solution,
    'Coal Heavy': np.array([400, 300, 100, 100, 100]),  # High coal usage
    'Renewable Focus': np.array([100, 200, 200, 200, 250])  # High renewable
}

scenario_costs = {}
scenario_emissions = {}

for name, mix in scenarios.items():
    if np.sum(mix) <= demand * 1.01:  # Allow small tolerance
        scenario_costs[name] = np.sum(mix * costs)
        scenario_emissions[name] = np.sum(mix * emissions)

x_pos = np.arange(len(scenario_costs))
cost_values = list(scenario_costs.values())
emission_values = list(scenario_emissions.values())

# Create dual y-axis plot
ax6_twin = ax6.twinx()

bars1 = ax6.bar([x - 0.2 for x in x_pos], cost_values, 0.4, 
                label='Cost', color='lightblue', alpha=0.7, edgecolor='black')
bars2 = ax6_twin.bar([x + 0.2 for x in x_pos], emission_values, 0.4, 
                     label='Emissions', color='lightcoral', alpha=0.7, edgecolor='black')

ax6.set_xlabel('Scenario', fontsize=12)
ax6.set_ylabel('Total Cost ($)', fontsize=12, color='blue')
ax6_twin.set_ylabel('Total Emissions (tons CO2)', fontsize=12, color='red')
ax6.set_title('Scenario Comparison', fontsize=14, fontweight='bold')
ax6.set_xticks(x_pos)
ax6.set_xticklabels(list(scenario_costs.keys()))

# Add value labels
for bar, value in zip(bars1, cost_values):
    height = bar.get_height()
    ax6.text(bar.get_x() + bar.get_width()/2., height + height*0.01,
             f'${value:,.0f}', ha='center', va='bottom', fontweight='bold', color='blue')

for bar, value in zip(bars2, emission_values):
    height = bar.get_height()
    ax6_twin.text(bar.get_x() + bar.get_width()/2., height + height*0.01,
                  f'{value:.0f}', ha='center', va='bottom', fontweight='bold', color='red')

# Add legends
ax6.legend(loc='upper left')
ax6_twin.legend(loc='upper right')

plt.tight_layout()
plt.show()

# Additional Analysis: Marginal Cost Analysis
print(f"\n" + "=" * 60)
print("MARGINAL COST ANALYSIS")
print("=" * 60)

# Calculate shadow prices (marginal costs)
print("Shadow Prices (Marginal Values):")
print(f"• Demand constraint: ${result.eqlin.marginals[0]:.2f}/MWh")
print(f"• Emission constraint: ${-result.ineqlin.marginals[0]:.2f}/ton CO2")

# Economic interpretation
print(f"\nEconomic Interpretation:")
print(f"• An additional 1 MWh of demand would increase costs by ${result.eqlin.marginals[0]:.2f}")
print(f"• Relaxing emission limit by 1 ton would save ${-result.ineqlin.marginals[0]:.2f}")

# Efficiency metrics
total_renewable = optimal_solution[3] + optimal_solution[4]  # Solar + Wind
renewable_percentage = (total_renewable / demand) * 100
clean_energy = optimal_solution[2] + optimal_solution[3] + optimal_solution[4]  # Nuclear + Solar + Wind
clean_percentage = (clean_energy / demand) * 100

print(f"\nSustainability Metrics:")
print(f"• Renewable energy share: {renewable_percentage:.1f}%")
print(f"• Clean energy share: {clean_percentage:.1f}%")
print(f"• Average cost per MWh: ${optimal_cost/demand:.2f}")
print(f"• Emissions intensity: {actual_emissions/demand:.3f} tons CO2/MWh")

# Risk analysis - what if one source fails?
print(f"\n" + "=" * 60)
print("RISK ANALYSIS: SOURCE FAILURE SCENARIOS")
print("=" * 60)

risk_analysis = []

for i, source in enumerate(sources):
    if optimal_solution[i] > 0:  # Only analyze sources we're using
        # Create modified capacity with one source unavailable
        modified_capacities = capacities.copy()
        modified_capacities[i] = 0
        
        # Modified bounds
        modified_bounds = [(0, cap) for cap in modified_capacities]
        
        # Solve with modified constraints
        result_risk = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, 
                             bounds=modified_bounds, method='highs')
        
        if result_risk.success:
            cost_increase = result_risk.fun - optimal_cost
            cost_increase_pct = (cost_increase / optimal_cost) * 100
            risk_analysis.append({
                'Failed Source': source,
                'Cost Increase ($)': cost_increase,
                'Cost Increase (%)': cost_increase_pct,
                'Feasible': 'Yes'
            })
        else:
            risk_analysis.append({
                'Failed Source': source,
                'Cost Increase ($)': 'N/A',
                'Cost Increase (%)': 'N/A',
                'Feasible': 'No'
            })

risk_df = pd.DataFrame(risk_analysis)
print(risk_df.to_string(index=False))

# Create risk visualization
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

# Risk analysis chart
feasible_risks = [r for r in risk_analysis if r['Feasible'] == 'Yes']
if feasible_risks:
    sources_at_risk = [r['Failed Source'] for r in feasible_risks]
    cost_increases = [r['Cost Increase (%)'] for r in feasible_risks]
    
    bars = ax1.bar(sources_at_risk, cost_increases, color='orange', alpha=0.7, edgecolor='black')
    ax1.set_ylabel('Cost Increase (%)', fontsize=12)
    ax1.set_title('Cost Impact of Source Failures', fontsize=14, fontweight='bold')
    ax1.set_xticklabels(sources_at_risk, rotation=45)
    
    # Add value labels
    for bar, increase in zip(bars, cost_increases):
        height = bar.get_height()
        ax1.text(bar.get_x() + bar.get_width()/2., height + height*0.01,
                 f'{increase:.1f}%', ha='center', va='bottom', fontweight='bold')

# Diversification index
ax2.pie(optimal_solution, labels=sources, autopct='%1.1f%%', startangle=90, colors=colors)
ax2.set_title('Energy Portfolio Diversification', fontsize=14, fontweight='bold')

plt.tight_layout()
plt.show()

print(f"\n" + "=" * 60)
print("SUMMARY AND RECOMMENDATIONS")
print("=" * 60)

print("Key Findings:")
print(f"1. Optimal mix achieves ${optimal_cost:,.0f} total cost")
print(f"2. {clean_percentage:.1f}% of energy comes from clean sources")
print(f"3. Emission constraint is {actual_emissions/max_emissions*100:.1f}% utilized")
print(f"4. Most vulnerable to {max(feasible_risks, key=lambda x: x['Cost Increase (%)'])['Failed Source'] if feasible_risks else 'N/A'} failure")

print(f"\nStrategic Recommendations:")
print(f"• Diversification reduces risk - no single source dominates")
print(f"• Emission constraints drive toward cleaner sources")
print(f"• Consider backup capacity for critical sources")
print(f"• Monitor fuel price volatility for cost-sensitive sources")

# Objective: minimize total cost
c = costs

# Constraint matrices
A_ub = emissions.reshape(1, -1)  # Emission constraint coefficients
b_ub = np.array([max_emissions])  # Emission limit

A_eq = np.ones((1, n_sources))   # Demand constraint (all sources sum to demand)
b_eq = np.array([demand])        # Total demand requirement