Equipment Criticality Assessment:
├─ Critical (Production stoppage, safety risk)
│   └─ Strategy: Predictive or rigorous preventive
├─ Important (Reduced capacity, quality impact)
│   └─ Strategy: Preventive with some condition monitoring
└─ Non-critical (Minimal impact)
    └─ Strategy: Reactive or basic preventive

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from datetime import datetime, timedelta

class ReliabilityAnalysis:
    """
    Equipment reliability analysis
    Calculate MTBF, MTTR, availability, and failure distributions
    """

    def __init__(self, failure_data):
        """
        failure_data: DataFrame with columns
        - equipment_id
        - failure_date
        - repair_date
        - time_between_failures (optional)
        - downtime_hours
        """
        self.data = failure_data.copy()

        # Calculate time between failures if not provided
        if 'time_between_failures' not in self.data.columns:
            self.data = self.data.sort_values(['equipment_id', 'failure_date'])
            self.data['time_between_failures'] = (
                self.data.groupby('equipment_id')['failure_date']
                .diff().dt.total_seconds() / 3600  # Convert to hours
            )

    def calculate_mtbf(self):
        """
        Calculate Mean Time Between Failures

        MTBF = Total operating time / Number of failures
        """

        # Remove first failure for each equipment (no prior TBF)
        valid_tbf = self.data[self.data['time_between_failures'].notna()]

        mtbf = valid_tbf['time_between_failures'].mean()
        mtbf_by_equipment = valid_tbf.groupby('equipment_id')['time_between_failures'].mean()

        return {
            'overall_mtbf_hours': mtbf,
            'overall_mtbf_days': mtbf / 24,
            'mtbf_by_equipment': mtbf_by_equipment.to_dict(),
            'std_dev': valid_tbf['time_between_failures'].std()
        }

    def calculate_mttr(self):
        """
        Calculate Mean Time To Repair

        MTTR = Total repair time / Number of repairs
        """

        mttr = self.data['downtime_hours'].mean()
        mttr_by_equipment = self.data.groupby('equipment_id')['downtime_hours'].mean()

        return {
            'overall_mttr_hours': mttr,
            'mttr_by_equipment': mttr_by_equipment.to_dict(),
            'std_dev': self.data['downtime_hours'].std()
        }

    def calculate_availability(self, operating_hours_per_year=8760):
        """
        Calculate equipment availability

        Availability = MTBF / (MTBF + MTTR)
        or
        Availability = Uptime / Total Time
        """

        mtbf_result = self.calculate_mtbf()
        mttr_result = self.calculate_mttr()

        mtbf = mtbf_result['overall_mtbf_hours']
        mttr = mttr_result['overall_mttr_hours']

        # Inherent availability (MTBF / (MTBF + MTTR))
        availability = mtbf / (mtbf + mttr) if (mtbf + mttr) > 0 else 0

        # Actual availability based on data
        total_downtime = self.data['downtime_hours'].sum()
        actual_availability = 1 - (total_downtime / operating_hours_per_year)

        return {
            'inherent_availability_pct': availability * 100,
            'actual_availability_pct': actual_availability * 100,
            'mtbf_hours': mtbf,
            'mttr_hours': mttr,
            'total_downtime_hours': total_downtime,
            'total_failures': len(self.data)
        }

    def weibull_analysis(self, time_to_failure_data):
        """
        Weibull distribution analysis for failure prediction

        Weibull parameters:
        - Beta (shape): failure pattern
          - Beta < 1: Infant mortality (decreasing failure rate)
          - Beta = 1: Random failures (constant failure rate)
          - Beta > 1: Wear-out failures (increasing failure rate)
        - Eta (scale): characteristic life

        Parameters:
        - time_to_failure_data: array of times to failure
        """

        # Fit Weibull distribution
        shape, loc, scale = stats.weibull_min.fit(time_to_failure_data, floc=0)

        # Beta = shape, Eta = scale
        beta = shape
        eta = scale

        # Determine failure pattern
        if beta < 1:
            pattern = "Infant Mortality (decreasing failure rate)"
        elif 0.9 <= beta <= 1.1:
            pattern = "Random Failures (constant failure rate)"
        else:
            pattern = "Wear-out (increasing failure rate)"

        # Calculate reliability at different times
        times = np.linspace(0, time_to_failure_data.max(), 100)
        reliability = stats.weibull_min.sf(times, shape, loc, scale)

        # Calculate MTTF (Mean Time To Failure)
        mttf = stats.weibull_min.mean(shape, loc, scale)

        return {
            'beta_shape': beta,
            'eta_scale': eta,
            'failure_pattern': pattern,
            'mttf': mttf,
            'times': times,
            'reliability': reliability
        }

    def plot_weibull(self, weibull_results):
        """Plot Weibull reliability curve"""

        fig, ax = plt.subplots(figsize=(10, 6))

        ax.plot(weibull_results['times'], weibull_results['reliability'] * 100,
               'b-', linewidth=2)
        ax.axhline(50, color='red', linestyle='--', linewidth=1, label='50% Reliability')

        ax.set_xlabel('Time (hours)', fontsize=12, fontweight='bold')
        ax.set_ylabel('Reliability (%)', fontsize=12, fontweight='bold')
        ax.set_title(f'Weibull Reliability Curve\nβ={weibull_results["beta_shape"]:.2f}, η={weibull_results["eta_scale"]:.0f}\nPattern: {weibull_results["failure_pattern"]}',
                    fontsize=14, fontweight='bold')
        ax.grid(True, alpha=0.3)
        ax.legend()

        plt.tight_layout()
        return fig

    def failure_rate_trend(self):
        """Analyze failure rate trends over time"""

        # Group failures by month
        self.data['month'] = pd.to_datetime(self.data['failure_date']).dt.to_period('M')
        monthly_failures = self.data.groupby('month').size()

        # Calculate trend
        x = np.arange(len(monthly_failures))
        y = monthly_failures.values

        if len(x) > 1:
            slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)

            if p_value < 0.05:
                if slope > 0:
                    trend = "Increasing (equipment degradation)"
                else:
                    trend = "Decreasing (reliability improvement)"
            else:
                trend = "Stable (no significant trend)"
        else:
            trend = "Insufficient data"
            slope = 0

        return {
            'monthly_failures': monthly_failures.to_dict(),
            'trend': trend,
            'slope': slope,
            'avg_failures_per_month': monthly_failures.mean()
        }

failure_dates.append(fail_date)
repair_dates.append(fail_date + timedelta(hours=downtime))
equipment_ids.append(equipment_id)
downtime_hours.append(downtime)

fig = reliability.plot_weibull(weibull)
plt.show()

---

class PreventiveMaintenanceOptimizer:
    """
    Optimize preventive maintenance intervals
    Balance maintenance cost vs. failure cost
    """

    def __init__(self, mtbf, mttr, pm_cost, failure_cost, pm_duration_hours):
        """
        Parameters:
        - mtbf: Mean Time Between Failures (hours)
        - mttr: Mean Time To Repair (hours)
        - pm_cost: Cost of preventive maintenance ($)
        - failure_cost: Cost of failure (downtime + repair) ($)
        - pm_duration_hours: Duration of PM activity (hours)
        """
        self.mtbf = mtbf
        self.mttr = mttr
        self.pm_cost = pm_cost
        self.failure_cost = failure_cost
        self.pm_duration = pm_duration_hours

    def calculate_optimal_interval(self, interval_range=(100, 2000, 50)):
        """
        Find optimal PM interval that minimizes total cost

        Total Cost = PM Cost + Failure Cost
        """

        intervals = np.arange(*interval_range)
        results = []

        for T in intervals:
            # Expected number of PMs per year (8760 hours)
            n_pms = 8760 / T

            # Probability of failure before PM
            # Assuming exponential distribution
            failure_prob = 1 - np.exp(-T / self.mtbf)

            # Expected failures per year
            n_failures = n_pms * failure_prob

            # Total cost
            annual_pm_cost = n_pms * self.pm_cost
            annual_failure_cost = n_failures * self.failure_cost
            total_cost = annual_pm_cost + annual_failure_cost

            # Downtime
            pm_downtime = n_pms * self.pm_duration
            failure_downtime = n_failures * self.mttr
            total_downtime = pm_downtime + failure_downtime

            # Availability
            availability = 1 - (total_downtime / 8760)

            results.append({
                'interval_hours': T,
                'n_pms_per_year': n_pms,
                'n_failures_per_year': n_failures,
                'annual_pm_cost': annual_pm_cost,
                'annual_failure_cost': annual_failure_cost,
                'total_annual_cost': total_cost,
                'total_downtime_hours': total_downtime,
                'availability_pct': availability * 100
            })

        df = pd.DataFrame(results)

        # Find optimal
        optimal_idx = df['total_annual_cost'].idxmin()
        optimal = df.loc[optimal_idx]

        return {
            'optimal_interval_hours': optimal['interval_hours'],
            'optimal_interval_days': optimal['interval_hours'] / 24,
            'total_annual_cost': optimal['total_annual_cost'],
            'availability': optimal['availability_pct'],
            'analysis': df
        }

    def plot_optimization(self, optimization_results):
        """Plot cost optimization curve"""

        df = optimization_results['analysis']

        fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))

        # Cost curves
        ax1.plot(df['interval_hours'], df['annual_pm_cost'], 'b-', linewidth=2, label='PM Cost')
        ax1.plot(df['interval_hours'], df['annual_failure_cost'], 'r-', linewidth=2, label='Failure Cost')
        ax1.plot(df['interval_hours'], df['total_annual_cost'], 'g-', linewidth=3, label='Total Cost')

        # Mark optimal
        optimal_interval = optimization_results['optimal_interval_hours']
        optimal_cost = optimization_results['total_annual_cost']
        ax1.plot(optimal_interval, optimal_cost, 'go', markersize=12, markeredgewidth=2,
                markerfacecolor='yellow', label='Optimal')

        ax1.set_xlabel('PM Interval (hours)', fontsize=12, fontweight='bold')
        ax1.set_ylabel('Annual Cost ($)', fontsize=12, fontweight='bold')
        ax1.set_title('Preventive Maintenance Interval Optimization', fontsize=14, fontweight='bold')
        ax1.legend()
        ax1.grid(True, alpha=0.3)

        # Availability curve
        ax2.plot(df['interval_hours'], df['availability_pct'], 'b-', linewidth=2)
        ax2.axhline(optimization_results['availability'], color='red', linestyle='--',
                   label=f"Optimal Availability: {optimization_results['availability']:.1f}%")

        ax2.set_xlabel('PM Interval (hours)', fontsize=12, fontweight='bold')
        ax2.set_ylabel('Availability (%)', fontsize=12, fontweight='bold')
        ax2.set_title('Equipment Availability vs. PM Interval', fontsize=12, fontweight='bold')
        ax2.legend()
        ax2.grid(True, alpha=0.3)

        plt.tight_layout()
        return fig

---

class PredictiveMaintenanceAnalyzer:
    """
    Predictive maintenance using condition monitoring data
    Anomaly detection and RUL (Remaining Useful Life) estimation
    """

    def __init__(self, sensor_data, sensor_name):
        """
        sensor_data: time series of sensor readings (e.g., vibration, temperature)
        sensor_name: description of sensor
        """
        self.data = np.array(sensor_data)
        self.sensor_name = sensor_name

    def detect_anomalies(self, threshold_sigma=3):
        """
        Detect anomalies using statistical methods
        Values beyond threshold_sigma standard deviations are flagged

        Returns anomaly indices and scores
        """

        mean = self.data.mean()
        std = self.data.std()

        # Z-score
        z_scores = np.abs((self.data - mean) / std)

        # Anomalies
        anomalies = z_scores > threshold_sigma
        anomaly_indices = np.where(anomalies)[0]

        return {
            'anomaly_count': anomalies.sum(),
            'anomaly_indices': anomaly_indices,
            'anomaly_values': self.data[anomaly_indices],
            'z_scores': z_scores,
            'mean': mean,
            'std': std,
            'threshold': threshold_sigma
        }

    def trend_analysis(self):
        """
        Analyze trend in sensor data
        Increasing trend may indicate degradation
        """

        x = np.arange(len(self.data))
        y = self.data

        # Linear regression
        slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)

        # Determine trend
        if p_value < 0.05:
            if slope > 0:
                trend = "Increasing (potential degradation)"
                concern_level = "High" if slope > std_err * 2 else "Medium"
            else:
                trend = "Decreasing (improving or cooling)"
                concern_level = "Low"
        else:
            trend = "Stable (no significant trend)"
            concern_level = "Low"

        return {
            'trend': trend,
            'slope': slope,
            'r_squared': r_value ** 2,
            'p_value': p_value,
            'concern_level': concern_level
        }

    def estimate_rul(self, failure_threshold):
        """
        Estimate Remaining Useful Life

        Assumes linear degradation to failure threshold

        Parameters:
        - failure_threshold: sensor value at which failure occurs

        Returns estimated RUL
        """

        trend = self.trend_analysis()

        if trend['slope'] <= 0:
            return {
                'rul_cycles': float('inf'),
                'message': 'No degradation detected - RUL cannot be estimated'
            }

        # Current value (most recent)
        current_value = self.data[-1]

        # Cycles until threshold
        rul_cycles = (failure_threshold - current_value) / trend['slope']

        # Check if already exceeded
        if rul_cycles < 0:
            rul_cycles = 0
            message = "Warning: Failure threshold already exceeded!"
        else:
            message = f"Estimated {rul_cycles:.0f} cycles remaining until failure threshold"

        return {
            'current_value': current_value,
            'failure_threshold': failure_threshold,
            'degradation_rate': trend['slope'],
            'rul_cycles': max(0, rul_cycles),
            'message': message
        }

    def maintenance_recommendation(self, rul_result, lead_time_cycles=50):
        """
        Generate maintenance recommendation based on RUL

        Parameters:
        - rul_result: output from estimate_rul()
        - lead_time_cycles: lead time needed to plan maintenance
        """

        rul = rul_result['rul_cycles']

        if rul == float('inf'):
            recommendation = "Continue monitoring - no immediate action needed"
            priority = "Low"
        elif rul <= 0:
            recommendation = "URGENT: Schedule immediate maintenance"
            priority = "Critical"
        elif rul < lead_time_cycles:
            recommendation = "Schedule maintenance soon - approaching failure threshold"
            priority = "High"
        elif rul < lead_time_cycles * 2:
            recommendation = "Plan maintenance within next planning cycle"
            priority = "Medium"
        else:
            recommendation = "Continue monitoring - sufficient remaining life"
            priority = "Low"

        return {
            'recommendation': recommendation,
            'priority': priority,
            'estimated_rul': rul,
            'lead_time': lead_time_cycles
        }

    def plot_condition_monitoring(self, anomaly_results, trend_results, rul_results=None):
        """Plot condition monitoring dashboard"""

        fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))

        # Time series with anomalies
        x = np.arange(len(self.data))
        ax1.plot(x, self.data, 'b-', linewidth=1, label='Sensor Reading')

        # Mark anomalies
        if len(anomaly_results['anomaly_indices']) > 0:
            ax1.plot(anomaly_results['anomaly_indices'],
                    anomaly_results['anomaly_values'],
                    'ro', markersize=8, label='Anomalies')

        # Mean and control limits
        mean = anomaly_results['mean']
        std = anomaly_results['std']
        threshold = anomaly_results['threshold']

        ax1.axhline(mean, color='green', linestyle='-', linewidth=2, label='Mean')
        ax1.axhline(mean + threshold * std, color='red', linestyle='--', linewidth=1.5, label='Upper Limit')
        ax1.axhline(mean - threshold * std, color='red', linestyle='--', linewidth=1.5, label='Lower Limit')

        # Trend line
        if trend_results['p_value'] < 0.05:
            trend_line = trend_results['slope'] * x + (self.data[0] - trend_results['slope'] * 0)
            ax1.plot(x, trend_line, 'orange', linestyle='--', linewidth=2, label='Trend')

        ax1.set_xlabel('Cycle / Time', fontsize=12, fontweight='bold')
        ax1.set_ylabel(f'{self.sensor_name}', fontsize=12, fontweight='bold')
        ax1.set_title(f'Condition Monitoring - {self.sensor_name}\nTrend: {trend_results["trend"]}',
                     fontsize=14, fontweight='bold')
        ax1.legend(loc='upper left')
        ax1.grid(True, alpha=0.3)

        # Z-scores (anomaly scores)
        ax2.plot(x, anomaly_results['z_scores'], 'b-', linewidth=1)
        ax2.axhline(threshold, color='red', linestyle='--', linewidth=2, label=f'Threshold ({threshold}σ)')
        ax2.fill_between(x, 0, threshold, alpha=0.2, color='green', label='Normal')
        ax2.fill_between(x, threshold, anomaly_results['z_scores'].max() + 1,
                        alpha=0.2, color='red', label='Anomaly Zone')

        ax2.set_xlabel('Cycle / Time', fontsize=12, fontweight='bold')
        ax2.set_ylabel('Anomaly Score (Z-score)', fontsize=12, fontweight='bold')
        ax2.set_title('Anomaly Detection', fontsize=12, fontweight='bold')
        ax2.legend()
        ax2.grid(True, alpha=0.3)

        plt.tight_layout()
        return fig

---

class OEEAnalysis:
    """
    Overall Equipment Effectiveness (OEE) calculation
    OEE = Availability × Performance × Quality
    """

    def __init__(self, shift_data):
        """
        shift_data: DataFrame with production shift information
        Columns:
        - planned_production_time (minutes)
        - downtime (minutes)
        - units_produced
        - ideal_cycle_time (minutes per unit)
        - defects
        """
        self.data = shift_data

    def calculate_oee(self):
        """
        Calculate OEE and its components

        Availability = Operating Time / Planned Production Time
        Performance = (Actual Output / Max Possible Output)
        Quality = Good Units / Total Units
        OEE = Availability × Performance × Quality
        """

        results = []

        for idx, row in self.data.iterrows():
            # Availability
            operating_time = row['planned_production_time'] - row['downtime']
            availability = operating_time / row['planned_production_time']

            # Performance
            ideal_time_for_output = row['units_produced'] * row['ideal_cycle_time']
            performance = ideal_time_for_output / operating_time if operating_time > 0 else 0

            # Quality
            good_units = row['units_produced'] - row['defects']
            quality = good_units / row['units_produced'] if row['units_produced'] > 0 else 0

            # OEE
            oee = availability * performance * quality

            results.append({
                'shift': idx,
                'availability': availability * 100,
                'performance': performance * 100,
                'quality': quality * 100,
                'oee': oee * 100,
                'operating_time': operating_time,
                'good_units': good_units
            })

        df_results = pd.DataFrame(results)

        # Overall averages
        overall = {
            'avg_availability': df_results['availability'].mean(),
            'avg_performance': df_results['performance'].mean(),
            'avg_quality': df_results['quality'].mean(),
            'avg_oee': df_results['oee'].mean(),
            'world_class_gap': 85 - df_results['oee'].mean(),  # World class = 85%+
            'results_by_shift': df_results
        }

        return overall

    def six_big_losses_analysis(self, breakdown_hours, setup_hours,
                                 small_stops_hours, speed_loss_pct,
                                 process_defects, reduced_yield_pct):
        """
        Analyze TPM Six Big Losses

        1. Breakdowns (Availability)
        2. Setup/Changeovers (Availability)
        3. Small Stops (Performance)
        4. Speed Loss (Performance)
        5. Process Defects (Quality)
        6. Reduced Yield (Quality)
        """

        total_time = self.data['planned_production_time'].sum() / 60  # Convert to hours

        losses = {
            'Breakdowns': {
                'hours': breakdown_hours,
                'pct_of_total': (breakdown_hours / total_time) * 100,
                'category': 'Availability'
            },
            'Setup/Changeover': {
                'hours': setup_hours,
                'pct_of_total': (setup_hours / total_time) * 100,
                'category': 'Availability'
            },
            'Small Stops': {
                'hours': small_stops_hours,
                'pct_of_total': (small_stops_hours / total_time) * 100,
                'category': 'Performance'
            },
            'Speed Loss': {
                'hours': total_time * (speed_loss_pct / 100),
                'pct_of_total': speed_loss_pct,
                'category': 'Performance'
            },
            'Process Defects': {
                'units': process_defects,
                'pct_of_total': (process_defects / self.data['units_produced'].sum()) * 100,
                'category': 'Quality'
            },
            'Reduced Yield': {
                'pct_of_total': reduced_yield_pct,
                'category': 'Quality'
            }
        }

        return pd.DataFrame(losses).T

    def plot_oee_dashboard(self, oee_results):
        """Create OEE dashboard visualization"""

        df = oee_results['results_by_shift']

        fig = plt.figure(figsize=(14, 10))
        gs = fig.add_gridspec(3, 2, hspace=0.3, wspace=0.3)

        # OEE over time
        ax1 = fig.add_subplot(gs[0, :])
        ax1.plot(df['shift'], df['oee'], 'bo-', linewidth=2, markersize=8, label='OEE')
        ax1.axhline(85, color='green', linestyle='--', linewidth=2, label='World Class (85%)')
        ax1.axhline(oee_results['avg_oee'], color='orange', linestyle='-', linewidth=2, label='Average')
        ax1.set_xlabel('Shift', fontsize=12, fontweight='bold')
        ax1.set_ylabel('OEE (%)', fontsize=12, fontweight='bold')
        ax1.set_title('Overall Equipment Effectiveness (OEE) by Shift', fontsize=14, fontweight='bold')
        ax1.legend()
        ax1.grid(True, alpha=0.3)

        # OEE components (average)
        ax2 = fig.add_subplot(gs[1, 0])
        components = ['Availability', 'Performance', 'Quality', 'OEE']
        values = [
            oee_results['avg_availability'],
            oee_results['avg_performance'],
            oee_results['avg_quality'],
            oee_results['avg_oee']
        ]
        colors = ['skyblue', 'lightgreen', 'lightcoral', 'gold']

        bars = ax2.bar(components, values, color=colors, edgecolor='black', linewidth=1.5)
        ax2.axhline(85, color='green', linestyle='--', linewidth=2, label='World Class')
        ax2.set_ylabel('Percentage (%)', fontsize=12, fontweight='bold')
        ax2.set_title('Average OEE Components', fontsize=12, fontweight='bold')
        ax2.set_ylim([0, 100])
        ax2.legend()
        ax2.grid(True, axis='y', alpha=0.3)

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

        # Component trends
        ax3 = fig.add_subplot(gs[1, 1])
        ax3.plot(df['shift'], df['availability'], 'o-', label='Availability', linewidth=2)
        ax3.plot(df['shift'], df['performance'], 's-', label='Performance', linewidth=2)
        ax3.plot(df['shift'], df['quality'], '^-', label='Quality', linewidth=2)
        ax3.set_xlabel('Shift', fontsize=12, fontweight='bold')
        ax3.set_ylabel('Percentage (%)', fontsize=12, fontweight='bold')
        ax3.set_title('OEE Component Trends', fontsize=12, fontweight='bold')
        ax3.legend()
        ax3.grid(True, alpha=0.3)

        # Performance metrics table
        ax4 = fig.add_subplot(gs[2, :])
        ax4.axis('tight')
        ax4.axis('off')

        table_data = [
            ['Metric', 'Value', 'World Class', 'Gap'],
            ['Availability', f"{oee_results['avg_availability']:.1f}%", '90%', f"{90 - oee_results['avg_availability']:.1f}%"],
            ['Performance', f"{oee_results['avg_performance']:.1f}%", '95%', f"{95 - oee_results['avg_performance']:.1f}%"],
            ['Quality', f"{oee_results['avg_quality']:.1f}%", '99%', f"{99 - oee_results['avg_quality']:.1f}%"],
            ['OEE', f"{oee_results['avg_oee']:.1f}%", '85%', f"{oee_results['world_class_gap']:.1f}%"]
        ]

        table = ax4.table(cellText=table_data, cellLoc='center', loc='center',
                         colWidths=[0.3, 0.2, 0.25, 0.25])
        table.auto_set_font_size(False)
        table.set_fontsize(11)
        table.scale(1, 2)

        # Style header row
        for i in range(4):
            table[(0, i)].set_facecolor('#4CAF50')
            table[(0, i)].set_text_props(weight='bold', color='white')

        plt.suptitle('OEE Performance Dashboard', fontsize=16, fontweight='bold', y=0.98)

        return fig