"""Synthetic fault-detection benchmark for two-channel beam monitoring. The simulation is an idealized single-degree-of-freedom beam model, not measured sensor or crack data. It evaluates a predeclared standardized fusion rule, temperature compensation residual, and a false-alarm budget. """ from __future__ import annotations import json import math import random from pathlib import Path ROOT = Path(__file__).resolve().parent SEED = 690304 WINDOW_RATE_HZ = 1 FALSE_ALARM_TARGET_PER_100_H = 0.5 HEALTHY_TEST_HOURS = 1000 HEALTHY_WINDOWS = HEALTHY_TEST_HOURS * 3600 * WINDOW_RATE_HZ DAMAGE_CASES = 20_000 BASE_STRAIN_MICROSTRAIN = 500.0 STRAIN_READ_NOISE_MICROSTRAIN = 11.0 THERMAL_EXPANSION_MICROSTRAIN_PER_C = 12.0 TEMP_COMPENSATION_RESIDUAL_SD_C = 0.5 MODAL_FREQUENCY_NOISE_REL = 0.01 def q(x: float) -> float: return 0.5 * math.erfc(x / math.sqrt(2)) def q_inverse(tail_probability: float) -> float: lo, hi = 0.0, 12.0 for _ in range(100): mid = (lo + hi) / 2 if q(mid) > tail_probability: lo = mid else: hi = mid return (lo + hi) / 2 def sensor_means(stiffness_retained: float) -> tuple[float, float]: # Constant-load beam assumption: strain varies inversely with stiffness; # first natural frequency varies as sqrt(stiffness). strain_delta = BASE_STRAIN_MICROSTRAIN * (1 / stiffness_retained - 1) strain_sigma = math.sqrt(STRAIN_READ_NOISE_MICROSTRAIN**2 + (THERMAL_EXPANSION_MICROSTRAIN_PER_C * TEMP_COMPENSATION_RESIDUAL_SD_C)**2) strain_z = strain_delta / strain_sigma frequency_shift = 1 - math.sqrt(stiffness_retained) frequency_z = frequency_shift / MODAL_FREQUENCY_NOISE_REL return strain_z, frequency_z def main() -> None: rng = random.Random(SEED) alpha = FALSE_ALARM_TARGET_PER_100_H / (100 * 3600 * WINDOW_RATE_HZ) cutoff = q_inverse(alpha) healthy_false = 0 # Equal, independent, standardized residuals after calibration/temperature correction. for _ in range(HEALTHY_WINDOWS): z_strain = rng.gauss(0.0, 1.0) z_modal = rng.gauss(0.0, 1.0) fused_score = (z_strain + z_modal) / math.sqrt(2) healthy_false += fused_score >= cutoff cases = [] for label, retained in [("10_percent_stiffness_loss", 0.90), ("4_percent_stiffness_loss", 0.96)]: mu_strain, mu_modal = sensor_means(retained) detected = 0 for _ in range(DAMAGE_CASES): z_strain = rng.gauss(mu_strain, 1.0) z_modal = rng.gauss(mu_modal, 1.0) detected += (z_strain + z_modal) / math.sqrt(2) >= cutoff cases.append({"damage_case": label, "stiffness_retained": retained, "standardized_strain_mean": mu_strain, "standardized_modal_shift_mean": mu_modal, "fused_recall": detected / DAMAGE_CASES, "detected": detected, "trials": DAMAGE_CASES, "target_recall_95_met": detected / DAMAGE_CASES >= 0.95}) false_rate_per_100h = healthy_false / HEALTHY_TEST_HOURS * 100 result = { "status": "synthetic software benchmark; no beam, sensor, or physical crack tested", "seed": SEED, "sampling": {"window_rate_hz": WINDOW_RATE_HZ, "healthy_test_hours": HEALTHY_TEST_HOURS, "healthy_windows": HEALTHY_WINDOWS, "damage_trials_each": DAMAGE_CASES}, "model": { "beam_assumption": "single-degree-of-freedom; fixed load; strain proportional to 1/stiffness; modal frequency proportional to sqrt(stiffness)", "nominal_strain_microstrain": BASE_STRAIN_MICROSTRAIN, "strain_read_noise_sd_microstrain": STRAIN_READ_NOISE_MICROSTRAIN, "thermal_expansion_microstrain_per_c": THERMAL_EXPANSION_MICROSTRAIN_PER_C, "temperature_compensation_residual_sd_c": TEMP_COMPENSATION_RESIDUAL_SD_C, "modal_frequency_relative_noise_sd": MODAL_FREQUENCY_NOISE_REL, "residual_assumption": "independent standard-normal strain/modal channels after calibration and thermal compensation", }, "classifier": { "rule": "(z_strain + z_modal_drop)/sqrt(2) >= one-sided Gaussian threshold", "false_alarm_budget_per_100h": FALSE_ALARM_TARGET_PER_100_H, "per_window_alpha": alpha, "standard_normal_threshold": cutoff, }, "healthy_holdout": {"false_alarms": healthy_false, "false_alarms_per_100h": false_rate_per_100h, "point_estimate_target_less_than_1_per100h_met": false_rate_per_100h < 1.0}, "damage_holdout": cases, "interpretation": [ "The 10% stiffness-loss scenario is an explicit synthetic test case, not a validated crack signature.", "The 4% case measures the chosen model's sensitivity boundary; it is not a universal minimum detectable crack.", "Real strain/modal channels can share temperature, load, mounting, timing, and power errors; independence may fail and invalidate the nominal false-alarm rate.", "A physical study needs calibrated sensors, independently characterized beam damage, blinded holdout specimens, and empirical environmental drift before using the stated thresholds.", ], } out = ROOT / "benchmark.json" out.write_text(json.dumps(result, ensure_ascii=False, indent=2) + "\n", encoding="utf-8") print(json.dumps(result, ensure_ascii=True, indent=2)) def export(): import csv import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt main() result=json.loads((ROOT/'benchmark.json').read_text()) def write(name,rows): with (ROOT/name).open('w',newline='',encoding='utf-8') as f: w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows) def wilson_low(k,n): z=1.6448536269514722;p=k/n;d=1+z*z/n return (p+z*z/(2*n)-z*math.sqrt(p*(1-p)/n+z*z/(4*n*n)))/d cutoff=result['classifier']['standard_normal_threshold'] curves=[] for i in range(151): loss=i/1000;a,b=sensor_means(1-loss) curves.append(dict(stiffness_loss_fraction=loss,mean_strain_z=a,mean_modal_z=b,analytical_detection_probability=q(cutoff-(a+b)/math.sqrt(2)))) correlation=[dict(correlation=r/100,healthy_score_variance=1+r/100,expected_false_alarms_per_100h=360000*q(cutoff/math.sqrt(1+r/100))) for r in range(91)] observed=[] for r in result['damage_holdout']: lo=wilson_low(r['detected'],r['trials']) observed.append(dict(stiffness_loss_fraction=round(1-r['stiffness_retained'],8),detected=r['detected'],trials=r['trials'],recall=r['fused_recall'],one_sided95_wilson_lower=lo,lower_bound_at_least95percent=lo>=.95)) # Exact Poisson one-sided upper confidence bound solves P(X<=observed|mu)=.05. count=result['healthy_holdout']['false_alarms'];lo=0.;hi=100. for _ in range(100): mid=(lo+hi)/2;cdf=math.exp(-mid)*sum(mid**j/math.factorial(j) for j in range(count+1)) if cdf>.05:lo=mid else:hi=mid upper=(lo+hi)/2/1000*100 result['healthy_holdout']['confidence_upper95_per100h']=upper result['healthy_holdout']['confidence_gate_passes']=upper<1 result['healthy_holdout']['gate_scope']='Point-estimate flag and one-sided 95% Poisson confidence gate are distinct; independent rare-event model only.' (ROOT/'benchmark.json').write_text(json.dumps(result,ensure_ascii=False,indent=2)+'\n',encoding='utf-8') write('detection-curve.csv',curves);write('correlation-sensitivity.csv',correlation);write('damage-observations.csv',observed) write('healthy-observations.csv',[dict(windows=HEALTHY_WINDOWS,hours=HEALTHY_TEST_HOURS,false_alarms=count,observed_per100h=count/10,poisson_upper95_per100h=upper,upper_bound_less_than1=upper<1)]) checks=dict(original_counts_match=count==5 and [r['detected'] for r in observed]==[19652,364],healthy_poisson_upper95_per100h=upper,healthy_confidence_gate_passes=upper<1,recall_confidence_bounds=observed) assert checks['original_counts_match'] (ROOT/'checks.json').write_text(json.dumps(checks,indent=2)+'\n') plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False}) fig,ax=plt.subplots(1,2,figsize=(10,4.4),layout='constrained') ax[0].plot([100*r['stiffness_loss_fraction'] for r in curves],[100*r['analytical_detection_probability'] for r in curves],color='#285ec7',label='Independent Gaussian model') for r in observed: x=100*r['stiffness_loss_fraction'];y=100*r['recall'];low=100*r['one_sided95_wilson_lower'] ax[0].plot(x,y,'o',color='#ad367f');ax[0].vlines(x,low,y,color='#ad367f',lw=2) ax[0].annotate(f"{r['detected']:,}/{r['trials']:,}",(x,y),xytext=(6,-18 if y>50 else 10),textcoords='offset points',fontsize=8) ax[0].axhline(95,ls=':',color='gray');ax[0].set(xlabel='Prescribed stiffness loss (%)',ylabel='Detection probability / recall (%)',title='A. Both damage scenarios matter',ylim=(-3,104));ax[0].legend(fontsize=8,loc='center left') ax[1].semilogy([r['correlation'] for r in correlation],[r['expected_false_alarms_per_100h'] for r in correlation],color='#285ec7') ax[1].axhline(.5,color='gray',ls=':',label='Nominal design budget: 0.5') ax[1].set(xlabel='Correlation between standardized channels',ylabel='Expected false alarms per 100 h',title='B. Keep the threshold; change correlation');ax[1].legend(fontsize=8) for a in ax:a.grid(alpha=.2) fig.suptitle('Synthetic structural monitoring • no beam or physical crack tested',fontsize=12) fig.supxlabel('Dots: original counts; downward stems: one-sided 95% Wilson lower bounds. Right panel is a new analytical sensitivity calculation.',fontsize=8) for ext in ('png','svg','pdf'):fig.savefig(ROOT/f'monitoring-study.{ext}',dpi=200,bbox_inches='tight') print(json.dumps(checks,indent=2)) if __name__=='__main__':export()