Manufacturing Defect Detection

import pandas as pd
import numpy as np
from scipy import stats

df = pd.DataFrame({
    "batch_id":    [f"B{i:03d}" for i in range(1, 21)],
    "dimension_mm":[24.98,25.02,25.11,24.94,25.28,25.04,25.18,25.31,
                    25.07,24.96,25.42,25.08,24.91,25.38,25.19,25.51,
                    25.12,25.44,25.03,25.29],
    "temp_c":      [182.4,183.1,184.8,182.0,191.2,183.4,185.2,190.8,
                    183.8,182.6,193.4,184.1,182.2,192.8,186.4,195.2,
                    184.8,193.1,183.2,191.8],
    "pressure_bar":[6.2,6.1,6.4,6.0,7.1,6.2,6.5,7.0,
                    6.3,6.1,7.3,6.3,6.0,7.2,6.6,7.6,
                    6.4,7.4,6.2,7.0],
    "vibration_g": [0.42,0.38,0.51,0.39,0.84,0.44,0.58,0.81,
                    0.46,0.40,0.92,0.47,0.38,0.88,0.62,1.04,
                    0.52,0.94,0.41,0.82],
    "speed_rpm":   [2840,2855,2830,2860,2780,2848,2822,2792,
                    2844,2858,2768,2842,2862,2774,2818,2752,
                    2836,2762,2852,2786],
    "defects":     [4,3,6,3,14,4,7,13,5,3,18,5,3,16,8,22,6,19,4,12],
    "units":       [250]*20
})

# ── Defect rate ───────────────────────────────────────────────────────────────
df["defect_rate"] = (df["defects"] / df["units"] * 100).round(2)
df["batch_num"]   = range(1, 21)

# ── X-bar control chart: establish limits from baseline (batches 1-10) ────────
BASELINE = 10   # first N batches used to estimate control limits
baseline = df.iloc[:BASELINE]

xbar   = baseline["dimension_mm"].mean()
sigma  = baseline["dimension_mm"].std(ddof=1)
UCL    = xbar + 3 * sigma
LCL    = xbar - 3 * sigma
UCL_2s = xbar + 2 * sigma   # 2σ zone for Western Electric rules
LCL_2s = xbar - 2 * sigma
UCL_1s = xbar + 1 * sigma
LCL_1s = xbar - 1 * sigma

print(f"X-bar Control Chart Parameters (baseline batches 1–{BASELINE}):")
print(f"  Process mean (X̄):  {xbar:.4f} mm")
print(f"  Process std  (σ):   {sigma:.4f} mm")
print(f"  UCL (+3σ):          {UCL:.4f} mm")
print(f"  LCL (-3σ):          {LCL:.4f} mm")
print(f"  Engineering USL:    25.20 mm  (target ± 0.20)")
print(f"  Engineering LSL:    24.80 mm")

# ── Flag out-of-control points ────────────────────────────────────────────────
df["z_score_dim"]   = ((df["dimension_mm"] - xbar) / sigma).round(3)
df["ooc_3sigma"]    = ((df["dimension_mm"] > UCL) | (df["dimension_mm"] < LCL)).astype(int)

print(f"\nControl chart — all batches:")
print(df[["batch_id","dimension_mm","z_score_dim","defect_rate","ooc_3sigma"]].to_string(index=False))
print(f"\nOut-of-control batches (3σ violations): {df['ooc_3sigma'].sum()}")

USL = 25.20   # engineering upper specification limit
LSL = 24.80   # engineering lower specification limit
x   = df["dimension_mm"].values

# ── Western Electric Rules ────────────────────────────────────────────────────
violations = []

# Rule 1: any point beyond 3σ (already computed — add WE label)
for i, row in df[df["ooc_3sigma"]==1].iterrows():
    violations.append((row["batch_id"], "WE Rule 1", "Point beyond 3σ control limit"))

# Rule 2: 6 consecutive points steadily increasing or decreasing (trend)
for i in range(5, len(x)):
    window = x[i-5:i+1]
    if all(window[j] < window[j+1] for j in range(5)):
        violations.append((df.iloc[i]["batch_id"], "WE Rule 2", "6 consecutive increasing points"))
    elif all(window[j] > window[j+1] for j in range(5)):
        violations.append((df.iloc[i]["batch_id"], "WE Rule 2", "6 consecutive decreasing points"))

# Rule 3: 2 of 3 consecutive points beyond 2σ on same side
for i in range(2, len(x)):
    window_z = df["z_score_dim"].values[i-2:i+1]
    if sum(z > 2 for z in window_z) >= 2:
        violations.append((df.iloc[i]["batch_id"], "WE Rule 3", "2 of 3 beyond +2σ"))
    if sum(z < -2 for z in window_z) >= 2:
        violations.append((df.iloc[i]["batch_id"], "WE Rule 3", "2 of 3 beyond -2σ"))

# Rule 4: 8 consecutive points on same side of centreline
for i in range(7, len(x)):
    window_z = df["z_score_dim"].values[i-7:i+1]
    if all(z > 0 for z in window_z):
        violations.append((df.iloc[i]["batch_id"], "WE Rule 4", "8 consecutive above centreline"))
    elif all(z < 0 for z in window_z):
        violations.append((df.iloc[i]["batch_id"], "WE Rule 4", "8 consecutive below centreline"))

viol_df = pd.DataFrame(violations, columns=["batch_id","rule","description"])
print(f"Western Electric Rule Violations ({len(viol_df)} total):")
print(viol_df.to_string(index=False))

# ── Process Capability Indices ────────────────────────────────────────────────
# Use baseline sigma (stable process estimate)
Cp  = (USL - LSL) / (6 * sigma)
Cpu = (USL - xbar) / (3 * sigma)
Cpl = (xbar - LSL) / (3 * sigma)
Cpk = min(Cpu, Cpl)

print(f"\nProcess Capability (baseline parameters):")
print(f"  Cp  = {Cp:.3f}  ({'Capable' if Cp >= 1.33 else 'Marginal' if Cp >= 1.0 else 'Not capable'})")
print(f"  Cpk = {Cpk:.3f}  ({'Meets standard' if Cpk >= 1.33 else 'Below automotive standard (1.33)'})")
print(f"  Cpu = {Cpu:.3f}  (upper half capability)")
print(f"  Cpl = {Cpl:.3f}  (lower half capability)")
print(f"\nProcess centring: mean = {xbar:.4f}, target = 25.0000")
print(f"  Off-target by: {(xbar - 25.0)*1000:.1f} µm (positive = shifted high)")

sensor_cols = ["temp_c","pressure_bar","vibration_g","speed_rpm"]

# ── Baseline stats from first 10 batches ─────────────────────────────────────
base_stats = baseline[sensor_cols].agg(["mean","std"])

# ── Z-scores for all batches ──────────────────────────────────────────────────
for col in sensor_cols:
    mu  = base_stats.loc["mean", col]
    std = base_stats.loc["std",  col]
    df[f"z_{col}"] = ((df[col] - mu) / std).round(3)
    df[f"anom_{col}"] = (df[f"z_{col}"].abs() > 3).astype(int)

z_cols   = [f"z_{c}" for c in sensor_cols]
anom_cols= [f"anom_{c}" for c in sensor_cols]
df["total_anomalies"] = df[anom_cols].sum(axis=1)

print("Sensor Z-scores and anomaly flags:")
display_cols = ["batch_id","defect_rate"] + z_cols + ["total_anomalies"]
print(df[display_cols].to_string(index=False))

# ── Co-occurrence: sensor anomalies vs high-defect batches ────────────────────
HIGH_DEFECT_THRESHOLD = 5.0   # % defect rate
df["high_defect"] = (df["defect_rate"] >= HIGH_DEFECT_THRESHOLD).astype(int)

print(f"\nSensor anomaly co-occurrence with high-defect batches (>{HIGH_DEFECT_THRESHOLD}%):")
for col in sensor_cols:
    anom_col = f"anom_{col}"
    both   = ((df[anom_col]==1) & (df["high_defect"]==1)).sum()
    anom_n = df[anom_col].sum()
    hd_n   = df["high_defect"].sum()
    prec   = both / anom_n if anom_n > 0 else 0
    recall = both / hd_n   if hd_n  > 0 else 0
    print(f"  {col:<16}  anomalies={anom_n}  co-occur={both}  "
          f"precision={prec:.0%}  recall={recall:.0%}")

print("Parameter–Defect Rate Correlation Analysis:")
print(f"{'Parameter':<16} {'Pearson r':>10} {'Spearman ρ':>11} {'p-value':>9} {'Sig':>5} {'Rank'}")
print("─" * 62)

corr_results = []
for col in sensor_cols:
    r_p,  p_p  = stats.pearsonr(df[col],  df["defect_rate"])
    r_sp, p_sp = stats.spearmanr(df[col], df["defect_rate"])
    corr_results.append({
        "param":    col,
        "pearson":  round(r_p,  4),
        "spearman": round(r_sp, 4),
        "p_val":    round(p_p,  4),
        "abs_r":    abs(r_p)
    })

corr_df = pd.DataFrame(corr_results).sort_values("abs_r", ascending=False).reset_index(drop=True)
corr_df["rank"] = range(1, len(corr_df)+1)

for _, row in corr_df.iterrows():
    sig = "***" if row["p_val"]<0.001 else "**" if row["p_val"]<0.01 else "*" if row["p_val"]<0.05 else "ns"
    direction = "↑" if row["pearson"] > 0 else "↓"
    print(f"  {row['param']:<14} {row['pearson']:>+10.4f} {row['spearman']:>+11.4f} "
          f"{row['p_val']:>9.4f} {sig:>5}   #{row['rank']} {direction}")

top_param = corr_df.iloc[0]["param"]
top_r     = corr_df.iloc[0]["pearson"]
print(f"\nPrimary root cause candidate: {top_param}")
print(f"  Pearson r = {top_r:.4f} — {'strong' if abs(top_r)>0.7 else 'moderate'} positive correlation with defect rate")

# Linear relationship: defect rate predicted by top parameter
slope_coef, intercept_coef, r_val, p_val, se = stats.linregress(df[top_param], df["defect_rate"])
print(f"\nSimple regression: defect_rate = {slope_coef:.3f} × {top_param} + ({intercept_coef:.3f})")
print(f"  R² = {r_val**2:.4f}  |  p = {p_val:.4f}")

UNIT_COST        = 42.50    # manufacturing cost per unit £
SCRAP_MULTIPLIER = 1.0      # scrapped parts = full unit cost lost
REWORK_SHARE     = 0.40     # 40% of defects can be reworked at 30% of unit cost
REWORK_COST      = 0.30
BATCHES_PER_DAY  = 8
DAYS_PER_YEAR    = 250

# Current vs baseline defect rates
baseline_rate = df.iloc[:BASELINE]["defect_rate"].mean() / 100
current_rate  = df["defect_rate"].mean() / 100

print(f"Defect rate comparison:")
print(f"  Baseline (B001–B010):  {baseline_rate*100:.2f}%")
print(f"  Current  (all 20):     {current_rate*100:.2f}%")
print(f"  Increase:              +{(current_rate-baseline_rate)*100:.2f}pp")

# Cost per batch
units_per_batch = 250
def batch_cost(rate):
    defects     = units_per_batch * rate
    scrap_n     = defects * (1 - REWORK_SHARE)
    rework_n    = defects * REWORK_SHARE
    return scrap_n * UNIT_COST * SCRAP_MULTIPLIER + rework_n * UNIT_COST * REWORK_COST

cost_baseline = batch_cost(baseline_rate)
cost_current  = batch_cost(current_rate)
cost_increase = cost_current - cost_baseline

annual_cost_current  = cost_current  * BATCHES_PER_DAY * DAYS_PER_YEAR
annual_cost_baseline = cost_baseline * BATCHES_PER_DAY * DAYS_PER_YEAR
annual_saving        = annual_cost_current - annual_cost_baseline

print(f"\nCost per batch (250 units × £{UNIT_COST}):")
print(f"  Baseline defect cost:  £{cost_baseline:,.2f}")
print(f"  Current  defect cost:  £{cost_current:,.2f}")
print(f"  Excess cost per batch: £{cost_increase:,.2f}")

print(f"\nAnnualised ({BATCHES_PER_DAY} batches/day × {DAYS_PER_YEAR} days):")
print(f"  Baseline annual cost:  £{annual_cost_baseline:,.0f}")
print(f"  Current  annual cost:  £{annual_cost_current:,.0f}")
print(f"  Annual excess cost:    £{annual_saving:,.0f}")
print(f"\nIf process restored to baseline: saving of £{annual_saving:,.0f}/year")

# ── Concise SPC dashboard ─────────────────────────────────────────────────────
print("=" * 60)
print("  LINE 3 SPC BRIEF — AXON PRECISION COMPONENTS")
print("=" * 60)

print(f"\n  PROCESS STATE:        OUT OF CONTROL")
print(f"  Control chart:        4 batches beyond 3σ UCL")
print(f"  WE rule violations:   7 (Rules 1, 2, 3 triggered)")
print(f"  Trend detected:       6 consecutive increasing batches (B010–B015)")

print(f"\n  CAPABILITY:")
print(f"  Cp  = {Cp:.3f}   (automotive standard: 1.33)")
print(f"  Cpk = {Cpk:.3f}   (FAIL — process NOT capable)")
print(f"  Mean shifted +39µm above target (25.039 vs 25.000)")

print(f"\n  ROOT CAUSE (ranked by correlation with defect rate):")
for _, row in corr_df.iterrows():
    direction = "HIGH" if row["pearson"] > 0 else "LOW"
    print(f"    #{row['rank']}  {row['param']:<16} r={row['pearson']:+.4f}  {direction}")

print(f"\n  MOST LIKELY MECHANISM:")
print(f"    Bearing wear → vibration ↑ → thermal load ↑ → temp/pressure ↑")
print(f"    → spindle instability → RPM ↓ → dimension drift → defects ↑")

print(f"\n  FINANCIAL IMPACT:")
print(f"    Excess defect cost:   £{annual_saving:,.0f}/year")
print(f"    Bearing replacement:  ~£10,000 (payback < 3 weeks)")

print(f"\n  RECOMMENDED ACTIONS:")
print(f"    1. STOP Line 3 before Monday's run")
print(f"    2. Inspect and replace main spindle bearing")
print(f"    3. Re-establish baseline with 10 qualification batches")
print(f"    4. Set vibration alarm threshold at baseline +2σ ({0.40 + 2*0.066:.3f} g)")
print(f"    5. Re-calculate Cpk — must reach 1.33 before full production resumes")
print("=" * 60)