698 lines
27 KiB
Python
698 lines
27 KiB
Python
# -*- coding: utf-8 -*-
|
||
"""
|
||
Problem 3: Sensitivity Analysis and Assumption Testing Figures
|
||
All data sourced from 整合输出.md (frozen model specification)
|
||
"""
|
||
|
||
import matplotlib.pyplot as plt
|
||
import numpy as np
|
||
import seaborn as sns
|
||
from math import pi
|
||
import os
|
||
|
||
# ==========================================
|
||
# Style Configuration
|
||
# ==========================================
|
||
plt.rcParams['font.family'] = 'Arial'
|
||
plt.rcParams['font.size'] = 11
|
||
plt.rcParams['axes.labelsize'] = 12
|
||
plt.rcParams['axes.titlesize'] = 13
|
||
plt.rcParams['figure.dpi'] = 150
|
||
|
||
colors = ['#2ecc71', '#3498db', '#9b59b6', '#e74c3c', '#f39c12', '#1abc9c', '#34495e']
|
||
|
||
os.makedirs('figures', exist_ok=True)
|
||
|
||
def save_fig(name):
|
||
plt.tight_layout()
|
||
plt.savefig(f'figures/{name}.png', dpi=300, bbox_inches='tight')
|
||
plt.close()
|
||
print(f"Saved {name}.png")
|
||
|
||
|
||
# ==========================================
|
||
# Data from 整合输出.md (FROZEN)
|
||
# ==========================================
|
||
|
||
# SOBOL_TABLE_v1
|
||
SOBOL_TABLE = {
|
||
'k_L': {'S_i': 0.412, 'ST_i': 0.445},
|
||
'k_C': {'S_i': 0.285, 'ST_i': 0.312},
|
||
'kappa': {'S_i': 0.164, 'ST_i': 0.198},
|
||
'k_N': {'S_i': 0.042, 'ST_i': 0.065},
|
||
'R_ref': {'S_i': 0.021, 'ST_i': 0.048},
|
||
'alpha_Q': {'S_i': 0.011, 'ST_i': 0.032},
|
||
}
|
||
|
||
# COMPUTE_LOG_v1
|
||
SOBOL_COMPUTE = {
|
||
'N_base': 512,
|
||
'D': 6,
|
||
'N_evals_total': 4096,
|
||
'seed': 20260201,
|
||
'sampling_scheme': 'Saltelli',
|
||
}
|
||
|
||
# SCENARIO_TTE_TABLE_v1
|
||
SCENARIO_TABLE = {
|
||
'S0': {'desc': 'Baseline', 'TTE': 4.60, 'delta': 0.00, 'reason': 'SOC_ZERO'},
|
||
'S1': {'desc': 'Brightness Reduced (0.5x)', 'TTE': 5.82, 'delta': 1.22, 'reason': 'SOC_ZERO'},
|
||
'S2': {'desc': 'CPU Reduced (0.5x)', 'TTE': 5.45, 'delta': 0.85, 'reason': 'SOC_ZERO'},
|
||
'S3': {'desc': 'Network Reduced (0.5x)', 'TTE': 4.92, 'delta': 0.32, 'reason': 'SOC_ZERO'},
|
||
'S4': {'desc': 'Poor Signal (Psi=0.2)', 'TTE': 2.78, 'delta': -1.82, 'reason': 'SOC_ZERO'},
|
||
'S5': {'desc': 'Cold Ambient (0C)', 'TTE': 3.15, 'delta': -1.45, 'reason': 'V_CUTOFF'},
|
||
'S6': {'desc': 'Hot Ambient (40C)', 'TTE': 4.98, 'delta': 0.38, 'reason': 'SOC_ZERO'},
|
||
'S7': {'desc': 'Background Cut (0.5x)', 'TTE': 4.74, 'delta': 0.14, 'reason': 'SOC_ZERO'},
|
||
}
|
||
|
||
# MECH_SIGNATURES_v1
|
||
MECH_SIGNATURES = {
|
||
'S0': {'avg_P': 3.22, 'max_I': 1.54, 'min_Delta': 8.15, 'avg_R0': 0.108, 'avg_Qeff': 4.00},
|
||
'S4': {'avg_P': 5.32, 'max_I': 2.45, 'min_Delta': 3.82, 'avg_R0': 0.112, 'avg_Qeff': 4.00},
|
||
'S5': {'avg_P': 3.28, 'max_I': 1.92, 'min_Delta': 0.85, 'avg_R0': 0.235, 'avg_Qeff': 3.52},
|
||
}
|
||
|
||
# UQ_SUMMARY_v1
|
||
UQ_SUMMARY = {
|
||
'mean': 4.6021,
|
||
'std': 0.0542,
|
||
'p10': 4.5314,
|
||
'p50': 4.6018,
|
||
'p90': 4.6725,
|
||
'CI95_low': 4.5959,
|
||
'CI95_high': 4.6083,
|
||
}
|
||
|
||
# REPRODUCIBILITY_v1
|
||
REPRODUCIBILITY = {
|
||
'seed': 20260201,
|
||
'M': 300,
|
||
'theta': 1/600,
|
||
'sigma': 0.02,
|
||
'dt': 1.0,
|
||
}
|
||
|
||
# Assumption Robustness Data (from document analysis)
|
||
ASSUMPTION_ROBUSTNESS = {
|
||
'Baseline': {'TTE': 4.60, 'delta': 0.0, 'delta_pct': 0.0},
|
||
'OCV Linear→Poly': {'TTE': 4.68, 'delta': 0.08, 'delta_pct': 1.7},
|
||
'CPL→CC': {'TTE': 5.12, 'delta': 0.52, 'delta_pct': 11.3},
|
||
'Lumped→FEM': {'TTE': 4.48, 'delta': -0.12, 'delta_pct': -2.6},
|
||
'Signal Exp→Lin': {'TTE': 5.49, 'delta': 0.89, 'delta_pct': 19.3},
|
||
'OU θ Change': {'TTE': 4.63, 'delta': 0.03, 'delta_pct': 0.7},
|
||
}
|
||
|
||
# Physical Decoupling Data
|
||
DECOUPLING_DATA = {
|
||
'Full Model (All Coupling)': {'TTE': 4.60, 'delta_pct': 0.0},
|
||
'No Thermal (dT/dt=0)': {'TTE': 4.85, 'delta_pct': 5.4},
|
||
'No CPL (I=const)': {'TTE': 5.12, 'delta_pct': 11.3},
|
||
'No Signal (Psi=0.9)': {'TTE': 6.42, 'delta_pct': 39.6},
|
||
'No R(T) (R0=const)': {'TTE': 4.73, 'delta_pct': 2.8},
|
||
'No Coupling (Linear)': {'TTE': 7.21, 'delta_pct': 56.7},
|
||
}
|
||
|
||
# Extreme Scenario Data
|
||
EXTREME_SCENARIOS = {
|
||
'E0: Baseline': {'TTE': 4.60, 'delta_pct': 0, 'confidence': 5},
|
||
'E1: Arctic (-10°C)': {'TTE': 1.68, 'delta_pct': -63.5, 'confidence': 3},
|
||
'E2: Basement (Weak Signal)': {'TTE': 1.85, 'delta_pct': -59.8, 'confidence': 4},
|
||
'E3: Desert (50°C)': {'TTE': 3.78, 'delta_pct': -17.8, 'confidence': 3},
|
||
'E4: Perfect Storm': {'TTE': 0.92, 'delta_pct': -80.0, 'confidence': 2},
|
||
'E5: Aged Battery (SOH=70%)': {'TTE': 3.22, 'delta_pct': -30.0, 'confidence': 4},
|
||
'E6: Gaming Max': {'TTE': 1.54, 'delta_pct': -66.5, 'confidence': 3},
|
||
}
|
||
|
||
# Usage Fluctuation Data (from OU process analysis)
|
||
FLUCTUATION_DATA = {
|
||
'Low (σ=0.01)': {'mean': 4.60, 'std': 0.027, 'p10': 4.57, 'p90': 4.63, 'cv': 0.59},
|
||
'Baseline (σ=0.02)': {'mean': 4.60, 'std': 0.054, 'p10': 4.53, 'p90': 4.67, 'cv': 1.18},
|
||
'High (σ=0.04)': {'mean': 4.60, 'std': 0.108, 'p10': 4.46, 'p90': 4.74, 'cv': 2.35},
|
||
'Extreme (σ=0.08)': {'mean': 4.60, 'std': 0.215, 'p10': 4.32, 'p90': 4.88, 'cv': 4.67},
|
||
}
|
||
|
||
# Load Model Comparison (CPL vs CC vs CR)
|
||
LOAD_MODEL_COMPARISON = {
|
||
'CPL (Ours)': {'TTE': 4.60, 'end_current': 1.01, 'current_increase': 46},
|
||
'CC (Constant Current)': {'TTE': 5.12, 'end_current': 0.69, 'current_increase': 0},
|
||
'CR (Constant Resistance)': {'TTE': 5.38, 'end_current': 0.59, 'current_increase': -14},
|
||
}
|
||
|
||
# Signal Mapping Validation Data
|
||
SIGNAL_MAPPING = {
|
||
0.9: {'measured': 0.78, 'exp_model': 0.80, 'lin_model': 0.82},
|
||
0.5: {'measured': 2.15, 'exp_model': 2.18, 'lin_model': 1.52},
|
||
0.2: {'measured': 5.32, 'exp_model': 5.28, 'lin_model': 2.22},
|
||
0.1: {'measured': 8.15, 'exp_model': 8.21, 'lin_model': 2.92},
|
||
}
|
||
|
||
# Amplification Factor by SOC
|
||
AMPLIFICATION_BY_SOC = {
|
||
'1.0-0.7 (Early)': 1.8,
|
||
'0.7-0.4 (Mid)': 2.7,
|
||
'0.4-0.2 (Late)': 3.5,
|
||
'0.2-0 (Critical)': 4.2,
|
||
}
|
||
|
||
|
||
# ==========================================
|
||
# Fig 14: Sobol Sensitivity Indices
|
||
# ==========================================
|
||
def plot_fig14():
|
||
"""
|
||
Sobol sensitivity indices from SOBOL_TABLE_v1
|
||
Shows first-order (S_i) and total-order (ST_i) indices
|
||
"""
|
||
params = list(SOBOL_TABLE.keys())
|
||
param_labels = [r'$k_L$ (Screen)', r'$k_C$ (CPU)', r'$\kappa$ (Signal)',
|
||
r'$k_N$ (Network)', r'$R_{ref}$ (Resistance)', r'$\alpha_Q$ (Capacity)']
|
||
S1 = [SOBOL_TABLE[p]['S_i'] for p in params]
|
||
ST = [SOBOL_TABLE[p]['ST_i'] for p in params]
|
||
interaction = [st - s for st, s in zip(ST, S1)]
|
||
|
||
fig, ax = plt.subplots(figsize=(11, 6))
|
||
|
||
x = np.arange(len(params))
|
||
width = 0.35
|
||
|
||
bars1 = ax.bar(x - width/2, S1, width, label=r'First-order $S_i$', color=colors[1], alpha=0.8, edgecolor='black')
|
||
bars2 = ax.bar(x + width/2, ST, width, label=r'Total-order $ST_i$', color=colors[4], alpha=0.8, edgecolor='black')
|
||
|
||
# Add interaction values as text
|
||
for i, (xi, sti, inter) in enumerate(zip(x, ST, interaction)):
|
||
ax.annotate(f'Int={inter:.3f}', xy=(xi + width/2, sti + 0.01),
|
||
ha='center', fontsize=9, color='gray')
|
||
|
||
ax.set_ylabel('Sobol Index')
|
||
ax.set_title(f'Global Sensitivity Analysis: Sobol Indices (N={SOBOL_COMPUTE["N_evals_total"]})')
|
||
ax.set_xticks(x)
|
||
ax.set_xticklabels(param_labels, rotation=15, ha='right')
|
||
ax.legend(loc='upper right')
|
||
ax.set_ylim(0, 0.55)
|
||
|
||
# Add cumulative line on secondary axis
|
||
ax2 = ax.twinx()
|
||
cumsum = np.cumsum(ST) / np.sum(ST) * 100
|
||
ax2.plot(x, cumsum, 'o--', color='red', alpha=0.7, linewidth=2, markersize=8, label='Cumulative %')
|
||
ax2.set_ylabel('Cumulative Contribution (%)', color='red')
|
||
ax2.tick_params(axis='y', labelcolor='red')
|
||
ax2.set_ylim(0, 110)
|
||
|
||
# Mark 75% threshold
|
||
ax2.axhline(75, color='red', linestyle=':', alpha=0.5)
|
||
ax2.text(2.5, 77, '75% (Top 3 params)', fontsize=9, color='red')
|
||
|
||
# Stats box
|
||
stats_text = (f"SOBOL_TABLE_v1:\n"
|
||
f"Top contributor: k_L (ST={SOBOL_TABLE['k_L']['ST_i']:.3f})\n"
|
||
f"Top 2 total: {SOBOL_TABLE['k_L']['ST_i'] + SOBOL_TABLE['k_C']['ST_i']:.1%}\n"
|
||
f"Max interaction: kappa ({max(interaction):.3f})")
|
||
ax.text(0.02, 0.98, stats_text, transform=ax.transAxes, fontsize=9,
|
||
verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.9))
|
||
|
||
save_fig('fig14_sobol_indices')
|
||
|
||
|
||
# ==========================================
|
||
# Fig 15: Assumption Robustness Waterfall
|
||
# ==========================================
|
||
def plot_fig15():
|
||
"""
|
||
Bar chart showing impact of changing modeling assumptions on TTE.
|
||
MCM O-Award style: clean, professional, data-driven visualization.
|
||
Data from ASSUMPTION_ROBUSTNESS
|
||
"""
|
||
# Extract data (skip baseline for comparison)
|
||
assumptions_full = list(ASSUMPTION_ROBUSTNESS.keys())
|
||
assumptions = assumptions_full[1:] # Exclude 'Baseline'
|
||
tte_vals = [ASSUMPTION_ROBUSTNESS[a]['TTE'] for a in assumptions]
|
||
delta_pcts = [ASSUMPTION_ROBUSTNESS[a]['delta_pct'] for a in assumptions]
|
||
|
||
baseline = 4.60
|
||
|
||
fig, ax = plt.subplots(figsize=(10, 6))
|
||
|
||
x = np.arange(len(assumptions))
|
||
|
||
# Vibrant color coding: bright green for robust, orange for moderate, bright red for critical
|
||
bar_colors = []
|
||
for dp in delta_pcts:
|
||
if abs(dp) > 10:
|
||
bar_colors.append('#e74c3c') # Red - Critical
|
||
elif abs(dp) > 5:
|
||
bar_colors.append('#f39c12') # Orange - Moderate
|
||
else:
|
||
bar_colors.append('#27ae60') # Green - Robust
|
||
|
||
# Create bars
|
||
bars = ax.bar(x, tte_vals, width=0.6, color=bar_colors, alpha=0.9,
|
||
edgecolor='#2C3E50', linewidth=1.5)
|
||
|
||
# Baseline reference line
|
||
ax.axhline(baseline, color='#2c3e50', linestyle='--', linewidth=2.5,
|
||
label=f'Baseline CPL Model ({baseline:.2f}h)', zorder=1)
|
||
|
||
# Add value labels ABOVE bars (not overlapping)
|
||
for i, (bar, tte, dp) in enumerate(zip(bars, tte_vals, delta_pcts)):
|
||
height = bar.get_height()
|
||
# Position label above the bar
|
||
y_pos = height + 0.12
|
||
sign = '+' if dp > 0 else ''
|
||
label = f'{tte:.2f}h\n({sign}{dp:.1f}%)'
|
||
ax.text(bar.get_x() + bar.get_width()/2, y_pos, label,
|
||
ha='center', va='bottom', fontsize=10, fontweight='bold')
|
||
|
||
# Styling
|
||
ax.set_xticks(x)
|
||
ax.set_xticklabels(assumptions, rotation=20, ha='right', fontsize=10)
|
||
ax.set_ylabel('Time to Exhaustion (hours)', fontsize=11, fontweight='bold')
|
||
ax.set_xlabel('Modeling Assumption Change', fontsize=11, fontweight='bold')
|
||
ax.set_title('Assumption Robustness Analysis: Impact on TTE Prediction',
|
||
fontsize=13, fontweight='bold', pad=15)
|
||
ax.set_ylim(0, 6.8)
|
||
|
||
# Add legend with color explanation - positioned at upper right, outside plot area
|
||
from matplotlib.patches import Patch
|
||
legend_elements = [
|
||
plt.Line2D([0], [0], color='#2c3e50', linestyle='--', linewidth=2.5, label=f'Baseline ({baseline:.2f}h)'),
|
||
Patch(facecolor='#27ae60', edgecolor='#2C3E50', label='Robust (|Δ| < 5%)'),
|
||
Patch(facecolor='#f39c12', edgecolor='#2C3E50', label='Moderate (5-10%)'),
|
||
Patch(facecolor='#e74c3c', edgecolor='#2C3E50', label='Critical (|Δ| > 10%)')
|
||
]
|
||
ax.legend(handles=legend_elements, loc='upper right', fontsize=9,
|
||
framealpha=0.95, edgecolor='gray', bbox_to_anchor=(0.99, 0.99))
|
||
|
||
# Grid for readability
|
||
ax.yaxis.grid(True, linestyle=':', alpha=0.6)
|
||
ax.set_axisbelow(True)
|
||
|
||
plt.tight_layout()
|
||
save_fig('fig15_assumption_robustness')
|
||
|
||
|
||
# ==========================================
|
||
# Fig 16: Physical Coupling Decoupling Analysis
|
||
# ==========================================
|
||
def plot_fig16():
|
||
"""
|
||
Shows impact of disabling each physical coupling mechanism.
|
||
Data from DECOUPLING_DATA
|
||
"""
|
||
experiments = list(DECOUPLING_DATA.keys())
|
||
tte_vals = [DECOUPLING_DATA[e]['TTE'] for e in experiments]
|
||
delta_pcts = [DECOUPLING_DATA[e]['delta_pct'] for e in experiments]
|
||
|
||
fig, ax = plt.subplots(figsize=(10, 6))
|
||
|
||
# Color by magnitude
|
||
bar_colors = [colors[0] if d == 0 else (colors[3] if d > 10 else colors[1]) for d in delta_pcts]
|
||
|
||
bars = ax.barh(experiments, tte_vals, color=bar_colors, alpha=0.8, edgecolor='black', height=0.6)
|
||
|
||
ax.axvline(4.60, color='red', linestyle='--', linewidth=2, label='Full Model (4.60h)')
|
||
|
||
# Add percentage labels
|
||
for bar, pct in zip(bars, delta_pcts):
|
||
width = bar.get_width()
|
||
label = f'+{pct:.1f}%' if pct > 0 else ('Baseline' if pct == 0 else f'{pct:.1f}%')
|
||
ax.text(width + 0.15, bar.get_y() + bar.get_height()/2, label,
|
||
va='center', fontsize=10, fontweight='bold')
|
||
|
||
ax.set_xlabel('TTE (hours)')
|
||
ax.set_title('Physical Coupling Decoupling Analysis: TTE Overestimation When Ignoring Feedback')
|
||
ax.set_xlim(0, 8)
|
||
ax.legend(loc='lower right')
|
||
ax.grid(axis='x', alpha=0.3)
|
||
|
||
# Key insight box
|
||
insight_text = ("Key Insight:\n"
|
||
"Signal-Power coupling: +39.6%\n"
|
||
"CPL feedback: +11.3%\n"
|
||
"All coupling off: +56.7%")
|
||
ax.text(0.98, 0.02, insight_text, transform=ax.transAxes, fontsize=9,
|
||
ha='right', va='bottom', bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.9))
|
||
|
||
save_fig('fig16_decoupling')
|
||
|
||
|
||
# ==========================================
|
||
# Fig 17: Extreme Scenario Stress Testing
|
||
# ==========================================
|
||
def plot_fig17():
|
||
"""
|
||
Bar chart showing TTE under extreme conditions.
|
||
Data from EXTREME_SCENARIOS
|
||
"""
|
||
scenarios = list(EXTREME_SCENARIOS.keys())
|
||
tte_vals = [EXTREME_SCENARIOS[s]['TTE'] for s in scenarios]
|
||
delta_pcts = [EXTREME_SCENARIOS[s]['delta_pct'] for s in scenarios]
|
||
confidences = [EXTREME_SCENARIOS[s]['confidence'] for s in scenarios]
|
||
|
||
fig, ax = plt.subplots(figsize=(11, 6))
|
||
|
||
# Color by severity (green to red)
|
||
norm_delta = [(d - min(delta_pcts)) / (max(delta_pcts) - min(delta_pcts) + 1e-6) for d in delta_pcts]
|
||
bar_colors = plt.cm.RdYlGn_r(norm_delta)
|
||
|
||
bars = ax.bar(scenarios, tte_vals, color=bar_colors, edgecolor='black', alpha=0.85)
|
||
|
||
ax.axhline(4.60, color='gray', linestyle='--', alpha=0.7, linewidth=2, label='Baseline (4.60h)')
|
||
|
||
# Add labels with confidence stars
|
||
for bar, d, c in zip(bars, delta_pcts, confidences):
|
||
height = bar.get_height()
|
||
stars = '*' * c
|
||
label = f'{d:+.0f}%\n[{stars}]' if d != 0 else f'Baseline\n[{stars}]'
|
||
ax.text(bar.get_x() + bar.get_width()/2, height + 0.15, label,
|
||
ha='center', va='bottom', fontsize=9)
|
||
|
||
ax.set_ylabel('TTE (hours)')
|
||
ax.set_title('Extreme Scenario Stress Testing (Confidence: * to *****)')
|
||
ax.set_ylim(0, 6)
|
||
plt.xticks(rotation=30, ha='right')
|
||
ax.legend(loc='upper right')
|
||
|
||
# Highlight perfect storm
|
||
ax.annotate('Perfect Storm\n(DANGER)', xy=(4, 0.92), xytext=(5, 2),
|
||
arrowprops=dict(arrowstyle='->', color='red', lw=2),
|
||
fontsize=10, color='red', fontweight='bold')
|
||
|
||
save_fig('fig17_extreme_scenarios')
|
||
|
||
|
||
# ==========================================
|
||
# Fig 18: Usage Pattern Fluctuation Impact
|
||
# ==========================================
|
||
def plot_fig18():
|
||
"""
|
||
Shows how usage pattern volatility affects TTE uncertainty.
|
||
Data from FLUCTUATION_DATA
|
||
"""
|
||
labels = list(FLUCTUATION_DATA.keys())
|
||
means = [FLUCTUATION_DATA[l]['mean'] for l in labels]
|
||
stds = [FLUCTUATION_DATA[l]['std'] for l in labels]
|
||
p10s = [FLUCTUATION_DATA[l]['p10'] for l in labels]
|
||
p90s = [FLUCTUATION_DATA[l]['p90'] for l in labels]
|
||
cvs = [FLUCTUATION_DATA[l]['cv'] for l in labels]
|
||
|
||
fig, ax = plt.subplots(figsize=(9, 6))
|
||
|
||
x = np.arange(len(labels))
|
||
bar_colors = [colors[0], colors[1], colors[4], colors[3]]
|
||
|
||
# Plot 90% confidence intervals
|
||
for i, (m, low, high, c) in enumerate(zip(means, p10s, p90s, bar_colors)):
|
||
ax.plot([i, i], [low, high], color='black', linewidth=3, alpha=0.4)
|
||
ax.plot(i, m, 'o', markersize=14, color=c, markeredgecolor='black', markeredgewidth=1.5)
|
||
ax.plot([i-0.15, i+0.15], [low, low], 'k-', linewidth=2)
|
||
ax.plot([i-0.15, i+0.15], [high, high], 'k-', linewidth=2)
|
||
|
||
# Add width and CV annotations
|
||
for i, (low, high, cv) in enumerate(zip(p10s, p90s, cvs)):
|
||
width = high - low
|
||
ax.annotate(f'Width: {width:.2f}h\nCV: {cv:.1f}%',
|
||
xy=(i, high + 0.03), ha='center', fontsize=9,
|
||
bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
|
||
|
||
ax.set_xticks(x)
|
||
ax.set_xticklabels(labels)
|
||
ax.set_ylabel('TTE (hours)')
|
||
ax.set_title('Usage Pattern Fluctuation: Impact on TTE Uncertainty (90% CI)')
|
||
ax.set_ylim(4.1, 5.1)
|
||
ax.axhline(4.60, color='gray', linestyle='--', alpha=0.5, label='Mean TTE')
|
||
ax.legend(loc='upper left')
|
||
ax.grid(axis='y', alpha=0.3)
|
||
|
||
# Robustness conclusion
|
||
ax.text(0.98, 0.02, 'Robust: CV < 2.5% for reasonable volatility',
|
||
transform=ax.transAxes, ha='right', va='bottom', fontsize=10, fontweight='bold',
|
||
bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.8))
|
||
|
||
save_fig('fig18_fluctuation')
|
||
|
||
|
||
# ==========================================
|
||
# Fig 19: CPL vs CC vs CR Load Model Comparison
|
||
# ==========================================
|
||
def plot_fig19():
|
||
"""
|
||
Compares three load models: CPL, CC, CR.
|
||
Shows TTE prediction differences and end-of-life current behavior.
|
||
MCM O-Award style: clear comparison with proper annotations.
|
||
"""
|
||
models = list(LOAD_MODEL_COMPARISON.keys())
|
||
tte_vals = [LOAD_MODEL_COMPARISON[m]['TTE'] for m in models]
|
||
currents = [LOAD_MODEL_COMPARISON[m]['end_current'] for m in models]
|
||
|
||
# Baseline values
|
||
baseline_tte = 4.60 # CPL model TTE
|
||
nominal_current = 0.69 # Nominal operating current
|
||
|
||
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5.5))
|
||
|
||
# Distinct colors for each model
|
||
bar_colors = ['#27ae60', '#3498db', '#9b59b6'] # Green, Blue, Purple
|
||
|
||
# ===== Left Panel: TTE Comparison =====
|
||
x1 = np.arange(len(models))
|
||
bars1 = ax1.bar(x1, tte_vals, width=0.6, color=bar_colors, alpha=0.85,
|
||
edgecolor='#2c3e50', linewidth=1.5)
|
||
|
||
# Baseline reference line (CPL model)
|
||
ax1.axhline(baseline_tte, color='#e74c3c', linestyle='--', linewidth=2,
|
||
label=f'CPL Baseline ({baseline_tte:.2f}h)')
|
||
|
||
# Add value labels ABOVE bars
|
||
for i, (bar, t) in enumerate(zip(bars1, tte_vals)):
|
||
height = bar.get_height()
|
||
delta_pct = (t - baseline_tte) / baseline_tte * 100
|
||
if abs(delta_pct) < 0.1:
|
||
label = f'{t:.2f}h\n(Baseline)'
|
||
else:
|
||
label = f'{t:.2f}h\n({delta_pct:+.1f}%)'
|
||
ax1.text(bar.get_x() + bar.get_width()/2, height + 0.12, label,
|
||
ha='center', va='bottom', fontsize=10, fontweight='bold')
|
||
|
||
ax1.set_xticks(x1)
|
||
ax1.set_xticklabels(models, fontsize=10)
|
||
ax1.set_ylabel('Time to Exhaustion (hours)', fontsize=11, fontweight='bold')
|
||
ax1.set_title('(a) TTE Prediction by Load Model', fontsize=12, fontweight='bold')
|
||
ax1.set_ylim(0, 6.5)
|
||
ax1.legend(loc='upper left', fontsize=9)
|
||
ax1.yaxis.grid(True, linestyle=':', alpha=0.5)
|
||
ax1.set_axisbelow(True)
|
||
|
||
# ===== Right Panel: End-of-Life Current =====
|
||
x2 = np.arange(len(models))
|
||
bars2 = ax2.bar(x2, currents, width=0.6, color=bar_colors, alpha=0.85,
|
||
edgecolor='#2c3e50', linewidth=1.5)
|
||
|
||
# Nominal current reference line
|
||
ax2.axhline(nominal_current, color='#7f8c8d', linestyle='--', linewidth=2,
|
||
label=f'Nominal Current ({nominal_current:.2f}A)')
|
||
|
||
# Measured range shading (28-45% increase from literature)
|
||
measured_low = nominal_current * 1.28
|
||
measured_high = nominal_current * 1.45
|
||
ax2.axhspan(measured_low, measured_high, alpha=0.25, color='#27ae60',
|
||
label=f'Measured Range ({measured_low:.2f}-{measured_high:.2f}A)')
|
||
|
||
# Add value labels ABOVE bars
|
||
for i, (bar, c) in enumerate(zip(bars2, currents)):
|
||
height = bar.get_height()
|
||
delta_pct = (c - nominal_current) / nominal_current * 100
|
||
label = f'{c:.2f}A\n({delta_pct:+.0f}%)'
|
||
ax2.text(bar.get_x() + bar.get_width()/2, height + 0.03, label,
|
||
ha='center', va='bottom', fontsize=10, fontweight='bold')
|
||
|
||
ax2.set_xticks(x2)
|
||
ax2.set_xticklabels(models, fontsize=10)
|
||
ax2.set_ylabel('Current at SOC=0.1 (A)', fontsize=11, fontweight='bold')
|
||
ax2.set_title('(b) End-of-Life Current Surge', fontsize=12, fontweight='bold')
|
||
ax2.set_ylim(0, 1.3)
|
||
ax2.legend(loc='upper right', fontsize=9)
|
||
ax2.yaxis.grid(True, linestyle=':', alpha=0.5)
|
||
ax2.set_axisbelow(True)
|
||
|
||
# Add key insight annotation on right panel
|
||
ax2.annotate('CPL captures\nreal-world surge',
|
||
xy=(0, 1.01), xytext=(0.8, 1.15),
|
||
fontsize=9, ha='center',
|
||
arrowprops=dict(arrowstyle='->', color='#27ae60', lw=1.5),
|
||
bbox=dict(boxstyle='round,pad=0.3', facecolor='#d5f5e3', edgecolor='#27ae60'))
|
||
|
||
plt.tight_layout()
|
||
save_fig('fig19_cpl_comparison')
|
||
|
||
|
||
# ==========================================
|
||
# Fig 20: Signal-Power Mapping Validation
|
||
# ==========================================
|
||
def plot_fig20():
|
||
"""
|
||
Validates exponential signal-power mapping against literature.
|
||
Shows linear model failure at low signal quality.
|
||
"""
|
||
psi_vals = list(SIGNAL_MAPPING.keys())
|
||
measured = [SIGNAL_MAPPING[p]['measured'] for p in psi_vals]
|
||
exp_model = [SIGNAL_MAPPING[p]['exp_model'] for p in psi_vals]
|
||
lin_model = [SIGNAL_MAPPING[p]['lin_model'] for p in psi_vals]
|
||
|
||
fig, ax = plt.subplots(figsize=(9, 6))
|
||
|
||
x = np.arange(len(psi_vals))
|
||
width = 0.25
|
||
|
||
ax.bar(x - width, measured, width, label='Measured (Literature)', color='gray', alpha=0.7, edgecolor='black')
|
||
ax.bar(x, exp_model, width, label='Exponential Model (Ours)', color=colors[0], alpha=0.8, edgecolor='black')
|
||
ax.bar(x + width, lin_model, width, label='Linear Model', color=colors[3], alpha=0.8, edgecolor='black')
|
||
|
||
ax.set_xlabel('Signal Quality Ψ')
|
||
ax.set_ylabel('Network Power (W)')
|
||
ax.set_title('Signal-Power Mapping: Model Validation')
|
||
ax.set_xticks(x)
|
||
ax.set_xticklabels([f'Ψ={p}' for p in psi_vals])
|
||
ax.legend(loc='upper left')
|
||
|
||
# Calculate and annotate errors
|
||
for i, (m, l) in enumerate(zip(measured, lin_model)):
|
||
err = (l - m) / m * 100
|
||
if abs(err) > 20:
|
||
ax.annotate(f'Error: {err:.0f}%',
|
||
xy=(i + width, l), xytext=(i + width + 0.3, l + 0.8),
|
||
arrowprops=dict(arrowstyle='->', color='red', lw=1.5),
|
||
color='red', fontsize=9, fontweight='bold')
|
||
|
||
# Add conclusion box
|
||
ax.text(0.98, 0.98, 'Linear model fails at Ψ < 0.3\n(Error > 50%)',
|
||
transform=ax.transAxes, ha='right', va='top', fontsize=10,
|
||
bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.9))
|
||
|
||
save_fig('fig20_signal_validation')
|
||
|
||
|
||
# ==========================================
|
||
# Fig 21: 3D Sensitivity Framework Radar
|
||
# ==========================================
|
||
def plot_fig21():
|
||
"""
|
||
Radar chart summarizing sensitivity across six dimensions.
|
||
Higher = more sensitive / less robust (risk indicator)
|
||
"""
|
||
categories = ['Temperature\nSensitivity', 'Signal\nSensitivity', 'Load\nSensitivity',
|
||
'Assumption\nRobustness', 'Fluctuation\nRobustness', 'Extreme\nResilience']
|
||
|
||
# Values: higher = worse (more sensitive or less robust)
|
||
# Based on analysis: Temp=4.8 (63.5% impact), Signal=4.5 (59.8%), Load=3.2,
|
||
# Assumption=4.0 (robust to most), Fluctuation=4.5 (CV<2.5%), Extreme=2.5 (E4 problematic)
|
||
values = [4.8, 4.5, 3.2, 4.0, 4.5, 2.5]
|
||
|
||
N = len(categories)
|
||
angles = [n / float(N) * 2 * pi for n in range(N)]
|
||
values_plot = values + values[:1]
|
||
angles += angles[:1]
|
||
|
||
fig, ax = plt.subplots(figsize=(8, 8), subplot_kw=dict(polar=True))
|
||
|
||
ax.plot(angles, values_plot, 'o-', linewidth=2, color=colors[2], markersize=8)
|
||
ax.fill(angles, values_plot, alpha=0.25, color=colors[2])
|
||
|
||
# Risk zones
|
||
theta = np.linspace(0, 2*pi, 100)
|
||
ax.fill_between(theta, 4, 5, alpha=0.15, color='red')
|
||
ax.fill_between(theta, 0, 2, alpha=0.15, color='green')
|
||
|
||
ax.set_xticks(angles[:-1])
|
||
ax.set_xticklabels(categories, size=10)
|
||
ax.set_ylim(0, 5)
|
||
ax.set_yticks([1, 2, 3, 4, 5])
|
||
ax.set_title('3D Sensitivity Framework Summary\n(Higher = More Sensitive / Less Robust)', pad=20)
|
||
|
||
# Add legend
|
||
from matplotlib.patches import Patch
|
||
legend_elements = [Patch(facecolor='red', alpha=0.3, label='High Risk Zone (>4)'),
|
||
Patch(facecolor='green', alpha=0.3, label='Safe Zone (<2)')]
|
||
ax.legend(handles=legend_elements, loc='upper right', bbox_to_anchor=(1.3, 1.0))
|
||
|
||
save_fig('fig21_framework_radar')
|
||
|
||
|
||
# ==========================================
|
||
# Fig 22: Fluctuation Amplification by SOC
|
||
# ==========================================
|
||
def plot_fig22():
|
||
"""
|
||
Shows how fluctuation amplification increases at low SOC.
|
||
CPL feedback causes "last 20% drops fast" phenomenon.
|
||
"""
|
||
soc_ranges = list(AMPLIFICATION_BY_SOC.keys())
|
||
beta_vals = list(AMPLIFICATION_BY_SOC.values())
|
||
|
||
fig, ax = plt.subplots(figsize=(9, 5))
|
||
|
||
# Color gradient (green to red)
|
||
colors_grad = plt.cm.Reds(np.linspace(0.3, 0.9, len(beta_vals)))
|
||
|
||
bars = ax.bar(soc_ranges, beta_vals, color=colors_grad, alpha=0.85, edgecolor='black')
|
||
|
||
ax.axhline(1.0, color='gray', linestyle='--', alpha=0.6, label='No Amplification (β=1)')
|
||
ax.set_ylabel('Amplification Factor β')
|
||
ax.set_xlabel('SOC Range')
|
||
ax.set_title('Fluctuation Amplification by Battery State: Why "Last 20% Drops Fast"')
|
||
ax.set_ylim(0, 5)
|
||
ax.legend(loc='upper left')
|
||
|
||
for bar, b in zip(bars, beta_vals):
|
||
ax.text(bar.get_x() + bar.get_width()/2, b + 0.15, f'{b:.1f}×',
|
||
ha='center', fontweight='bold', fontsize=11)
|
||
|
||
# Physical explanation
|
||
ax.text(0.98, 0.98, 'Physics: I=P/V causes\npositive feedback at low V',
|
||
transform=ax.transAxes, ha='right', va='top', fontsize=10,
|
||
bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.9))
|
||
|
||
ax.grid(axis='y', alpha=0.3)
|
||
|
||
save_fig('fig22_amplification')
|
||
|
||
|
||
# ==========================================
|
||
# Main Execution
|
||
# ==========================================
|
||
if __name__ == "__main__":
|
||
print("=" * 60)
|
||
print("Problem 3 Figure Generation (Data from 整合输出.md)")
|
||
print("=" * 60)
|
||
|
||
print("\n[1/9] Fig 14: Sobol Sensitivity Indices...")
|
||
plot_fig14()
|
||
|
||
print("[2/9] Fig 15: Assumption Robustness Waterfall...")
|
||
plot_fig15()
|
||
|
||
print("[3/9] Fig 16: Physical Decoupling Analysis...")
|
||
plot_fig16()
|
||
|
||
print("[4/9] Fig 17: Extreme Scenario Testing...")
|
||
plot_fig17()
|
||
|
||
print("[5/9] Fig 18: Usage Pattern Fluctuation...")
|
||
plot_fig18()
|
||
|
||
print("[6/9] Fig 19: CPL vs CC vs CR Comparison...")
|
||
plot_fig19()
|
||
|
||
print("[7/9] Fig 20: Signal Mapping Validation...")
|
||
plot_fig20()
|
||
|
||
print("[8/9] Fig 21: 3D Framework Radar...")
|
||
plot_fig21()
|
||
|
||
print("[9/9] Fig 22: SOC Amplification...")
|
||
plot_fig22()
|
||
|
||
print("\n" + "=" * 60)
|
||
print("All Problem 3 figures generated successfully!")
|
||
print("Output directory: figures/")
|
||
print("=" * 60)
|