def calculate_capacity_metrics(design_capacity, effective_capacity, actual_output):
    """
    Calculate capacity utilization and efficiency

    Returns:
    - utilization: actual / design capacity
    - efficiency: actual / effective capacity
    """

    utilization = (actual_output / design_capacity) * 100
    efficiency = (actual_output / effective_capacity) * 100

    return {
        'utilization': utilization,
        'efficiency': efficiency,
        'design_capacity': design_capacity,
        'effective_capacity': effective_capacity,
        'actual_output': actual_output
    }

# Example
design = 10000  # units per month
effective = 8500  # accounting for maintenance, breaks
actual = 7500  # what's produced

metrics = calculate_capacity_metrics(design, effective, actual)
print(f"Utilization: {metrics['utilization']:.1f}%")  # 75%
print(f"Efficiency: {metrics['efficiency']:.1f}%")    # 88.2%

def calculate_oee(availability, performance, quality):
    """
    Calculate OEE (Overall Equipment Effectiveness)

    Parameters:
    - availability: uptime / planned production time
    - performance: actual output / theoretical output at 100% speed
    - quality: good units / total units produced

    World-class OEE: > 85%
    """

    oee = availability * performance * quality * 100

    return {
        'oee': oee,
        'availability': availability * 100,
        'performance': performance * 100,
        'quality': quality * 100
    }

# Example
availability = 0.90  # 90% uptime
performance = 0.85   # 85% of theoretical speed
quality = 0.95       # 95% good units

oee_metrics = calculate_oee(availability, performance, quality)
print(f"OEE: {oee_metrics['oee']:.1f}%")  # 72.7%

import pandas as pd
import numpy as np

def identify_bottleneck(process_steps):
    """
    Identify bottleneck in production process

    Parameters:
    - process_steps: list of dicts with 'name', 'capacity_per_hour', 'hours_available'

    Returns bottleneck step and throughput
    """

    df = pd.DataFrame(process_steps)

    # Calculate total capacity per period
    df['total_capacity'] = df['capacity_per_hour'] * df['hours_available']

    # Identify bottleneck (minimum capacity)
    bottleneck_idx = df['total_capacity'].idxmin()
    bottleneck = df.loc[bottleneck_idx]

    # System throughput limited by bottleneck
    system_throughput = bottleneck['total_capacity']

    # Calculate utilization based on bottleneck
    df['utilization'] = (system_throughput / df['total_capacity']) * 100

    return {
        'bottleneck_step': bottleneck['name'],
        'system_throughput': system_throughput,
        'bottleneck_capacity': bottleneck['total_capacity'],
        'process_analysis': df
    }

# Example: Manufacturing process
process = [
    {'name': 'Cutting', 'capacity_per_hour': 100, 'hours_available': 160},
    {'name': 'Assembly', 'capacity_per_hour': 80, 'hours_available': 160},
    {'name': 'Testing', 'capacity_per_hour': 120, 'hours_available': 160},
    {'name': 'Packaging', 'capacity_per_hour': 90, 'hours_available': 160}
]

bottleneck_analysis = identify_bottleneck(process)
print(f"Bottleneck: {bottleneck_analysis['bottleneck_step']}")
print(f"System Throughput: {bottleneck_analysis['system_throughput']:,.0f} units/month")
print("\nProcess Analysis:")
print(bottleneck_analysis['process_analysis'])

class DrumBufferRope:
    """
    Theory of Constraints scheduling method

    - Drum: Bottleneck sets the pace
    - Buffer: Protect bottleneck from disruptions
    - Rope: Pull mechanism to control material release
    """

    def __init__(self, bottleneck_capacity, buffer_time_days=3):
        self.bottleneck_capacity = bottleneck_capacity
        self.buffer_time = buffer_time_days

    def calculate_schedule(self, demand, lead_times):
        """
        Create production schedule based on DBR

        Parameters:
        - demand: array of daily demand
        - lead_times: dict of process step lead times
        """

        schedule = []

        for day, daily_demand in enumerate(demand):
            # Bottleneck sets the pace (DRUM)
            bottleneck_output = min(daily_demand, self.bottleneck_capacity)

            # Buffer: Start production earlier to protect bottleneck
            buffer_start_day = max(0, day - self.buffer_time)

            # Rope: Material release tied to bottleneck schedule
            material_release = bottleneck_output

            schedule.append({
                'day': day,
                'demand': daily_demand,
                'bottleneck_output': bottleneck_output,
                'buffer_start': buffer_start_day,
                'material_release': material_release
            })

        return pd.DataFrame(schedule)

# Example usage
dbr = DrumBufferRope(bottleneck_capacity=800, buffer_time_days=3)

# Daily demand for next 10 days
demand = np.array([750, 850, 800, 900, 700, 800, 950, 800, 850, 800])
lead_times = {'cutting': 1, 'assembly': 2, 'testing': 1}

schedule = dbr.calculate_schedule(demand, lead_times)
print(schedule)

from pulp import *
import pandas as pd
import numpy as np

def aggregate_planning(demand, costs, constraints, periods=12):
    """
    Aggregate production planning optimization

    Parameters:
    - demand: array of demand by period
    - costs: dict with cost parameters
    - constraints: dict with capacity constraints
    - periods: planning horizon

    Returns optimal plan
    """

    # Create problem
    prob = LpProblem("Aggregate_Planning", LpMinimize)

    # Decision variables
    P = LpVariable.dicts("Production", range(periods), lowBound=0)
    W = LpVariable.dicts("Workforce", range(periods), lowBound=0, cat='Integer')
    O = LpVariable.dicts("Overtime", range(periods), lowBound=0)
    I = LpVariable.dicts("Inventory", range(periods), lowBound=0)
    H = LpVariable.dicts("Hire", range(periods), lowBound=0, cat='Integer')
    F = LpVariable.dicts("Fire", range(periods), lowBound=0, cat='Integer')
    B = LpVariable.dicts("Backorder", range(periods), lowBound=0)
    S = LpVariable.dicts("Subcontract", range(periods), lowBound=0)

    # Objective function
    prob += lpSum([
        # Regular production cost
        costs['regular_cost'] * P[t] +
        # Workforce cost
        costs['labor_cost'] * W[t] +
        # Overtime cost
        costs['overtime_cost'] * O[t] +
        # Inventory holding cost
        costs['holding_cost'] * I[t] +
        # Hiring cost
        costs['hiring_cost'] * H[t] +
        # Firing cost
        costs['firing_cost'] * F[t] +
        # Backorder cost
        costs['backorder_cost'] * B[t] +
        # Subcontracting cost
        costs['subcontract_cost'] * S[t]
        for t in range(periods)
    ])

    # Constraints

    # Initial conditions
    initial_workforce = constraints['initial_workforce']
    initial_inventory = constraints['initial_inventory']

    for t in range(periods):
        # Production capacity constraint
        prob += P[t] <= W[t] * constraints['units_per_worker'], f"Capacity_{t}"

        # Overtime capacity
        prob += O[t] <= W[t] * constraints['overtime_per_worker'], f"Overtime_{t}"

        # Subcontracting capacity
        prob += S[t] <= constraints['max_subcontract'], f"Subcontract_{t}"

        # Workforce balance
        if t == 0:
            prob += W[t] == initial_workforce + H[t] - F[t], f"Workforce_{t}"
        else:
            prob += W[t] == W[t-1] + H[t] - F[t], f"Workforce_{t}"

        # Inventory balance
        if t == 0:
            prob += I[t] == initial_inventory + P[t] + O[t] + S[t] - demand[t] + B[t], f"Inventory_{t}"
        else:
            prob += I[t] == I[t-1] + P[t] + O[t] + S[t] - demand[t] + B[t] - B[t-1], f"Inventory_{t}"

        # Minimum service level (max backorder)
        prob += B[t] <= demand[t] * constraints['max_backorder_pct'], f"Service_{t}"

    # Solve
    prob.solve(PULP_CBC_CMD(msg=0))

    # Extract results
    results = {
        'status': LpStatus[prob.status],
        'total_cost': value(prob.objective),
        'production': [P[t].varValue for t in range(periods)],
        'workforce': [W[t].varValue for t in range(periods)],
        'overtime': [O[t].varValue for t in range(periods)],
        'inventory': [I[t].varValue for t in range(periods)],
        'hired': [H[t].varValue for t in range(periods)],
        'fired': [F[t].varValue for t in range(periods)],
        'backorders': [B[t].varValue for t in range(periods)],
        'subcontract': [S[t].varValue for t in range(periods)]
    }

    # Create summary DataFrame
    df = pd.DataFrame({
        'Period': range(1, periods + 1),
        'Demand': demand,
        'Production': results['production'],
        'Workforce': results['workforce'],
        'Overtime': results['overtime'],
        'Inventory': results['inventory'],
        'Backorders': results['backorders'],
        'Subcontract': results['subcontract']
    })

    results['summary'] = df

    return results

# Example usage
demand = np.array([1000, 1200, 1500, 1800, 2000, 1800,
                   1600, 1400, 1300, 1200, 1100, 1000])

costs = {
    'regular_cost': 50,
    'labor_cost': 2000,      # per worker per period
    'overtime_cost': 75,
    'holding_cost': 5,
    'hiring_cost': 1000,
    'firing_cost': 1500,
    'backorder_cost': 100,
    'subcontract_cost': 80
}

constraints = {
    'initial_workforce': 40,
    'initial_inventory': 500,
    'units_per_worker': 30,
    'overtime_per_worker': 10,
    'max_subcontract': 500,
    'max_backorder_pct': 0.10
}

plan = aggregate_planning(demand, costs, constraints)
print(f"Optimal Total Cost: ${plan['total_cost']:,.0f}")
print("\nProduction Plan:")
print(plan['summary'])

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

class CapacityRequirementsPlanning:
    """
    CRP - Calculate and analyze capacity requirements
    based on production schedule and routings
    """

    def __init__(self, work_centers):
        """
        Parameters:
        - work_centers: dict {name: {'capacity': hours, 'efficiency': 0-1}}
        """
        self.work_centers = work_centers

    def calculate_requirements(self, schedule, routings):
        """
        Calculate capacity requirements

        Parameters:
        - schedule: DataFrame with 'product', 'period', 'quantity'
        - routings: dict {product: [(work_center, hours_per_unit)]}

        Returns DataFrame with requirements by work center and period
        """

        requirements = []

        for idx, row in schedule.iterrows():
            product = row['product']
            period = row['period']
            quantity = row['quantity']

            # Get routing for product
            routing = routings.get(product, [])

            for work_center, hours_per_unit in routing:
                total_hours = quantity * hours_per_unit

                # Adjust for efficiency
                efficiency = self.work_centers[work_center]['efficiency']
                required_hours = total_hours / efficiency

                requirements.append({
                    'period': period,
                    'work_center': work_center,
                    'product': product,
                    'quantity': quantity,
                    'hours_required': required_hours
                })

        df = pd.DataFrame(requirements)

        # Aggregate by work center and period
        summary = df.groupby(['period', 'work_center'])['hours_required'].sum().reset_index()

        # Add capacity and utilization
        summary['capacity'] = summary['work_center'].map(
            lambda wc: self.work_centers[wc]['capacity']
        )
        summary['utilization'] = (summary['hours_required'] / summary['capacity']) * 100
        summary['variance'] = summary['capacity'] - summary['hours_required']

        return summary

    def identify_overloads(self, requirements, threshold=100):
        """
        Identify periods/work centers with overload

        Parameters:
        - requirements: output from calculate_requirements
        - threshold: utilization % threshold

        Returns overloaded resources
        """

        overloads = requirements[requirements['utilization'] > threshold].copy()
        overloads = overloads.sort_values(['period', 'work_center'])

        return overloads

    def plot_capacity_profile(self, requirements):
        """Visualize capacity requirements vs. available"""

        work_centers = requirements['work_center'].unique()

        fig, axes = plt.subplots(len(work_centers), 1,
                                figsize=(12, 4 * len(work_centers)),
                                squeeze=False)

        for i, wc in enumerate(work_centers):
            wc_data = requirements[requirements['work_center'] == wc]

            ax = axes[i, 0]

            # Plot capacity line
            ax.axhline(y=wc_data['capacity'].iloc[0],
                      color='green', linestyle='--',
                      linewidth=2, label='Capacity')

            # Plot requirements
            ax.bar(wc_data['period'], wc_data['hours_required'],
                  alpha=0.7, label='Required')

            # Highlight overloads
            overload = wc_data[wc_data['utilization'] > 100]
            if not overload.empty:
                ax.bar(overload['period'], overload['hours_required'],
                      color='red', alpha=0.7, label='Overload')

            ax.set_title(f'{wc} - Capacity Profile')
            ax.set_xlabel('Period')
            ax.set_ylabel('Hours')
            ax.legend()
            ax.grid(True, alpha=0.3)

        plt.tight_layout()
        return fig

# Example usage
work_centers = {
    'Cutting': {'capacity': 160, 'efficiency': 0.90},
    'Welding': {'capacity': 160, 'efficiency': 0.85},
    'Assembly': {'capacity': 160, 'efficiency': 0.92},
    'Inspection': {'capacity': 160, 'efficiency': 0.95}
}

crp = CapacityRequirementsPlanning(work_centers)

# Production schedule
schedule = pd.DataFrame({
    'product': ['A', 'A', 'B', 'B', 'C', 'C'] * 3,
    'period': [1, 2, 1, 2, 1, 2] * 3,
    'quantity': [100, 120, 80, 90, 60, 70] * 3
})

# Routings: hours per unit at each work center
routings = {
    'A': [('Cutting', 0.5), ('Welding', 0.8), ('Assembly', 1.0), ('Inspection', 0.3)],
    'B': [('Cutting', 0.6), ('Assembly', 1.2), ('Inspection', 0.4)],
    'C': [('Cutting', 0.4), ('Welding', 1.0), ('Assembly', 0.8), ('Inspection', 0.2)]
}

# Calculate requirements
requirements = crp.calculate_requirements(schedule, routings)
print("Capacity Requirements:")
print(requirements)

# Identify overloads
overloads = crp.identify_overloads(requirements)
if not overloads.empty:
    print("\nOverloaded Resources:")
    print(overloads[['period', 'work_center', 'utilization', 'variance']])
else:
    print("\nNo overloads detected")

# Plot
crp.plot_capacity_profile(requirements)

import numpy as np

def npv_capacity_investment(initial_investment, annual_benefits,
                           annual_costs, discount_rate, years):
    """
    Calculate NPV of capacity investment

    Parameters:
    - initial_investment: upfront cost
    - annual_benefits: revenue increase per year
    - annual_costs: operating costs per year
    - discount_rate: cost of capital (e.g., 0.10 for 10%)
    - years: investment horizon
    """

    cash_flows = [-initial_investment]

    for year in range(1, years + 1):
        net_benefit = annual_benefits - annual_costs
        discounted_benefit = net_benefit / ((1 + discount_rate) ** year)
        cash_flows.append(discounted_benefit)

    npv = sum(cash_flows)

    # Calculate IRR (Internal Rate of Return)
    irr = np.irr([-initial_investment] + [annual_benefits - annual_costs] * years)

    # Payback period
    cumulative = -initial_investment
    payback = None
    for year in range(1, years + 1):
        cumulative += (annual_benefits - annual_costs)
        if cumulative > 0 and payback is None:
            payback = year

    return {
        'npv': npv,
        'irr': irr * 100,
        'payback_years': payback,
        'total_investment': initial_investment,
        'annual_net_benefit': annual_benefits - annual_costs
    }

# Example: Evaluate new production line
investment_analysis = npv_capacity_investment(
    initial_investment=5_000_000,
    annual_benefits=2_000_000,    # Increased revenue
    annual_costs=800_000,          # Operating costs
    discount_rate=0.12,            # 12% cost of capital
    years=10
)

print("Investment Analysis:")
print(f"  NPV: ${investment_analysis['npv']:,.0f}")
print(f"  IRR: {investment_analysis['irr']:.1f}%")
print(f"  Payback: {investment_analysis['payback_years']} years")

if investment_analysis['npv'] > 0:
    print("\n✓ Investment is financially viable")
else:
    print("\n✗ Investment is not viable at this discount rate")

import matplotlib.pyplot as plt
import numpy as np

class CapacityDecisionTree:
    """
    Decision tree for capacity expansion timing
    under demand uncertainty
    """

    def __init__(self):
        self.scenarios = []

    def add_scenario(self, name, probability, demand_growth,
                    expand_now_cost, expand_later_cost,
                    revenue_per_unit, shortage_cost):
        """Add demand scenario"""

        self.scenarios.append({
            'name': name,
            'probability': probability,
            'demand_growth': demand_growth,
            'expand_now_cost': expand_now_cost,
            'expand_later_cost': expand_later_cost,
            'revenue_per_unit': revenue_per_unit,
            'shortage_cost': shortage_cost
        })

    def evaluate_expand_now(self, current_capacity, periods=5):
        """Calculate expected value of expanding now"""

        ev = 0

        for scenario in self.scenarios:
            # Cost of expanding now
            cost = -scenario['expand_now_cost']

            # Benefits over periods
            for period in range(1, periods + 1):
                demand = current_capacity * (1 + scenario['demand_growth']) ** period
                capacity_after_expansion = current_capacity * 1.5  # Assume 50% expansion

                # Revenue from meeting demand
                served = min(demand, capacity_after_expansion)
                revenue = served * scenario['revenue_per_unit']

                # Discount
                discounted_revenue = revenue / (1.10 ** period)
                cost += discounted_revenue

            # Weight by probability
            ev += cost * scenario['probability']

        return ev

    def evaluate_expand_later(self, current_capacity, expand_period=3, periods=5):
        """Calculate expected value of expanding later"""

        ev = 0

        for scenario in self.scenarios:
            cost = 0

            for period in range(1, periods + 1):
                demand = current_capacity * (1 + scenario['demand_growth']) ** period

                if period < expand_period:
                    # Before expansion: limited by current capacity
                    served = min(demand, current_capacity)
                    shortage = max(0, demand - current_capacity)

                    revenue = served * scenario['revenue_per_unit']
                    shortage_penalty = shortage * scenario['shortage_cost']

                    discounted_value = (revenue - shortage_penalty) / (1.10 ** period)

                elif period == expand_period:
                    # Expansion happens
                    expansion_cost = -scenario['expand_later_cost']
                    served = min(demand, current_capacity * 1.5)
                    revenue = served * scenario['revenue_per_unit']

                    discounted_value = (expansion_cost + revenue) / (1.10 ** period)

                else:
                    # After expansion
                    served = min(demand, current_capacity * 1.5)
                    revenue = served * scenario['revenue_per_unit']

                    discounted_value = revenue / (1.10 ** period)

                cost += discounted_value

            # Weight by probability
            ev += cost * scenario['probability']

        return ev

    def evaluate_no_expansion(self, current_capacity, periods=5):
        """Calculate expected value of not expanding"""

        ev = 0

        for scenario in self.scenarios:
            cost = 0

            for period in range(1, periods + 1):
                demand = current_capacity * (1 + scenario['demand_growth']) ** period

                # Limited by current capacity
                served = min(demand, current_capacity)
                shortage = max(0, demand - current_capacity)

                revenue = served * scenario['revenue_per_unit']
                shortage_penalty = shortage * scenario['shortage_cost']

                discounted_value = (revenue - shortage_penalty) / (1.10 ** period)
                cost += discounted_value

            # Weight by probability
            ev += cost * scenario['probability']

        return ev

    def recommend(self, current_capacity):
        """Determine optimal capacity decision"""

        ev_now = self.evaluate_expand_now(current_capacity)
        ev_later = self.evaluate_expand_later(current_capacity)
        ev_no = self.evaluate_no_expansion(current_capacity)

        results = {
            'expand_now': ev_now,
            'expand_later': ev_later,
            'no_expansion': ev_no
        }

        best_option = max(results, key=results.get)

        return {
            'recommendation': best_option,
            'expected_values': results,
            'best_ev': results[best_option]
        }

# Example usage
dt = CapacityDecisionTree()

# Add demand scenarios
dt.add_scenario(
    name='High Growth',
    probability=0.30,
    demand_growth=0.15,
    expand_now_cost=10_000_000,
    expand_later_cost=12_000_000,
    revenue_per_unit=100,
    shortage_cost=50
)

dt.add_scenario(
    name='Moderate Growth',
    probability=0.50,
    demand_growth=0.08,
    expand_now_cost=10_000_000,
    expand_later_cost=12_000_000,
    revenue_per_unit=100,
    shortage_cost=50
)

dt.add_scenario(
    name='Low Growth',
    probability=0.20,
    demand_growth=0.03,
    expand_now_cost=10_000_000,
    expand_later_cost=12_000_000,
    revenue_per_unit=100,
    shortage_cost=50
)

# Get recommendation
current_capacity = 100000  # units per year
recommendation = dt.recommend(current_capacity)

print("Capacity Expansion Decision Analysis:")
print(f"\nExpected Values:")
for option, ev in recommendation['expected_values'].items():
    print(f"  {option}: ${ev:,.0f}")

print(f"\n✓ Recommendation: {recommendation['recommendation'].upper()}")
print(f"  Expected Value: ${recommendation['best_ev']:,.0f}")

def value_of_flexibility(demand_scenarios, fixed_capacity,
                        flexible_capacity, flexibility_cost):
    """
    Calculate value of flexible capacity using real options approach

    Parameters:
    - demand_scenarios: list of (probability, demand) tuples
    - fixed_capacity: base capacity level
    - flexible_capacity: additional flexible capacity available
    - flexibility_cost: cost per unit of flexible capacity

    Returns value of flexibility option
    """

    # Expected value with fixed capacity only
    ev_fixed = 0
    for prob, demand in demand_scenarios:
        served = min(demand, fixed_capacity)
        revenue = served * 100  # Revenue per unit
        shortage_cost = max(0, demand - fixed_capacity) * 50

        ev_fixed += prob * (revenue - shortage_cost)

    # Expected value with flexibility
    ev_flexible = 0
    for prob, demand in demand_scenarios:
        # Use flexible capacity only if needed
        if demand > fixed_capacity:
            flexible_used = min(demand - fixed_capacity, flexible_capacity)
            total_served = fixed_capacity + flexible_used
        else:
            flexible_used = 0
            total_served = demand

        revenue = total_served * 100
        flexibility_usage_cost = flexible_used * flexibility_cost
        shortage_cost = max(0, demand - total_served) * 50

        ev_flexible += prob * (revenue - flexibility_usage_cost - shortage_cost)

    value_of_flexibility = ev_flexible - ev_fixed

    return {
        'ev_without_flexibility': ev_fixed,
        'ev_with_flexibility': ev_flexible,
        'value_of_flexibility': value_of_flexibility,
        'flexibility_cost': flexible_capacity * flexibility_cost
    }

# Example
demand_scenarios = [
    (0.20, 8000),   # Low demand
    (0.50, 10000),  # Expected demand
    (0.30, 13000)   # High demand
]

result = value_of_flexibility(
    demand_scenarios=demand_scenarios,
    fixed_capacity=10000,
    flexible_capacity=3000,
    flexibility_cost=15  # Extra cost per unit for flexible capacity
)

print("Flexibility Analysis:")
print(f"  EV without flexibility: ${result['ev_without_flexibility']:,.0f}")
print(f"  EV with flexibility: ${result['ev_with_flexibility']:,.0f}")
print(f"  Value of flexibility: ${result['value_of_flexibility']:,.0f}")

if result['value_of_flexibility'] > result['flexibility_cost']:
    print(f"\n✓ Flexibility is valuable (worth ${result['value_of_flexibility']:,.0f})")
else:
    print(f"\n✗ Flexibility not cost-effective")