预处理1

2026-01-31 13:43:21 +08:00
parent 84fb77d060
commit 8f767b45b6
3 changed files with 202 additions and 0 deletions
--- a/A题/分析/框架2/TTE预测与不确定性量化.py
+++ b/A题/分析/框架2/TTE预测与不确定性量化.py
@@ -0,0 +1,256 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import seaborn as sns
+from scipy.stats import spearmanr
+import os
+import time
+
+# ==========================================
+# 1. 配置与初始化
+# ==========================================
+def configure_environment():
+    """
+    配置绘图环境，解决中文乱码和负号显示问题，
+    并设置符合学术规范的绘图风格。
+    """
+    # 设置字体为 Times New Roman (英文) 或 SimHei (中文兼容)
+    # 为了美赛(MCM)标准，主要使用英文标签，但配置好中文支持以防万一
+    plt.rcParams['font.family'] = 'sans-serif'
+    plt.rcParams['font.sans-serif'] = ['SimHei', 'Arial'] # 优先使用黑体显示中文
+    plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示为方块的问题
+    
+    # 学术风格配置
+    plt.style.use('seaborn-v0_8-whitegrid')
+    plt.rcParams.update({
+        'font.size': 12,
+        'axes.labelsize': 14,
+        'axes.titlesize': 16,
+        'xtick.labelsize': 12,
+        'ytick.labelsize': 12,
+        'legend.fontsize': 12,
+        'figure.figsize': (10, 6),
+        'figure.dpi': 150
+    })
+
+    # 确保输出目录存在
+    if not os.path.exists('output_q2'):
+        os.makedirs('output_q2')
+
+# ==========================================
+# 2. 核心物理模型类
+# ==========================================
+class StochasticBatteryModel:
+    def __init__(self):
+        # 电池物理参数
+        self.Q_design = 4000 / 1000.0  # 4000 mAh -> 4.0 Ah
+        self.R_int_base = 0.15         # 内阻 (Ohm)
+        
+        # Shepherd 模型参数 (OCV 曲线)
+        self.E0 = 3.4
+        self.K = 0.02
+        self.A = 0.5
+        self.B = 15.0
+
+    def get_capacity_correction(self, temp_c):
+        """
+        根据温度修正有效容量 (Arrhenius 效应)
+        Temp_c: 环境温度 (摄氏度)
+        """
+        temp_k = temp_c + 273.15
+        ref_k = 298.15 # 25°C
+        # 温度越低，容量越小；温度越高，容量略微增加但有限制
+        factor = np.exp(0.8 * (1 - ref_k / temp_k))
+        return np.clip(factor, 0.4, 1.1) # 限制修正因子范围
+
+    def get_ocv(self, soc):
+        """计算开路电压 (Open Circuit Voltage)"""
+        soc = np.clip(soc, 0.001, 1.0)
+        # Shepherd Model + 线性项
+        term1 = self.E0
+        term2 = -self.K / soc
+        term3 = -self.R_int_base * 0.1 # 简化极化项
+        term4 = self.A * np.exp(-self.B * (1 - soc))
+        term5 = 0.7 * soc # 线性部分
+        v_ocv = 3.0 + term5 + term4 + term2 
+        return np.clip(v_ocv, 2.8, 4.35)
+
+    def calculate_power_drain(self, base_power, signal_dbm, noise_scale=0.1):
+        """
+        计算瞬时功率，包含信号强度的非线性影响和随机噪声
+        P_total = P_base + P_signal + Noise
+        """
+        # 1. 信号功耗模型：信号越弱(越负)，功耗呈指数上升
+        # 假设 -60dBm 为基准，每下降 10dBm 功耗显著增加
+        # 归一化信号强度：将 -120dBm 到 -50dBm 映射到 0-1 之间
+        sig_norm = np.clip((-signal_dbm - 50) / 70, 0, 1)
+        p_net = 2.0 * (sig_norm ** 3) # 三次幂关系，模拟弱信号时的急剧恶化
+        
+        # 2. 随机噪声 (模拟 CPU 动态调频)
+        noise = np.random.normal(0, noise_scale)
+        
+        return max(0.1, base_power + p_net + noise)
+
+# ==========================================
+# 3. 蒙特卡洛模拟引擎
+# ==========================================
+def run_monte_carlo_simulation(n_simulations=2000):
+    """
+    执行 N 次蒙特卡洛模拟
+    """
+    print(f"Starting Monte Carlo Simulation with {n_simulations} samples...")
+    
+    model = StochasticBatteryModel()
+    results = []
+    
+    # 定义时间步长 (秒)
+    dt = 60.0  # 1分钟一步，平衡精度与速度
+    max_time = 48 * 3600 # 最大模拟 48 小时
+    
+    start_time = time.time()
+
+    for i in range(n_simulations):
+        # --- A. 随机生成初始条件 (输入分布) ---
+        
+        # 1. 初始电量 SOC0: 均匀分布 U(0.2, 1.0)
+        # 模拟用户在不同电量下开始使用手机
+        soc_0 = np.random.uniform(0.2, 1.0)
+        
+        # 2. 环境温度 Temp: 正态分布 N(25, 8)，截断在 [-10, 45]
+        temp_c = np.random.normal(25, 8)
+        temp_c = np.clip(temp_c, -10, 45)
+        
+        # 3. 平均信号强度 Signal: 偏态分布
+        # 大部分时候信号较好 (-80dBm), 偶尔很差 (-110dBm)
+        signal_dbm = np.random.triangular(-120, -85, -50)
+        
+        # 4. 用户行为模式 (基准功率)
+        # 混合高斯分布：待机为主(0.5W)，偶尔重度使用(4W)
+        user_type = np.random.choice(['Light', 'Medium', 'Heavy'], p=[0.4, 0.4, 0.2])
+        if user_type == 'Light':
+            base_power_avg = 0.5 # 待机、阅读
+        elif user_type == 'Medium':
+            base_power_avg = 1.5 # 视频、社交
+        else:
+            base_power_avg = 4.0 # 游戏
+            
+        # --- B. 执行单次时间步进模拟 (Euler Integration) ---
+        soc = soc_0
+        t = 0
+        q_eff = model.Q_design * model.get_capacity_correction(temp_c)
+        
+        while soc > 0.02 and t < max_time: # 截止电压设为 2%
+            # 计算当前功率 (加入随机性)
+            p_inst = model.calculate_power_drain(base_power_avg, signal_dbm, noise_scale=0.2)
+            
+            # 计算端电压
+            v_term = model.get_ocv(soc) - (p_inst / 3.7) * model.R_int_base
+            v_term = max(v_term, 3.0) # 防止电压过低
+            
+            # 计算电流 I = P / V
+            current = p_inst / v_term
+            
+            # 更新 SOC: dS = -I * dt / Q
+            d_soc = -current * dt / (q_eff * 3600)
+            soc += d_soc
+            t += dt
+            
+        # --- C. 记录结果 ---
+        tte_hours = t / 3600.0
+        results.append({
+            'Initial_SOC': soc_0,
+            'Temperature_C': temp_c,
+            'Signal_dBm': signal_dbm,
+            'User_Profile': user_type,
+            'Base_Power': base_power_avg,
+            'TTE_Hours': tte_hours
+        })
+        
+        if (i+1) % 500 == 0:
+            print(f"  Processed {i+1}/{n_simulations} simulations...")
+
+    print(f"Simulation completed in {time.time() - start_time:.2f} seconds.")
+    return pd.DataFrame(results)
+
+# ==========================================
+# 4. 结果分析与可视化
+# ==========================================
+def analyze_and_plot(df):
+    """
+    对模拟结果进行统计分析和绘图
+    """
+    print("\n=== Analysis Report ===")
+    print(df.describe())
+    
+    # --- 图1: TTE 分布直方图 (Probability Distribution) ---
+    plt.figure(figsize=(10, 6))
+    sns.histplot(data=df, x='TTE_Hours', hue='User_Profile', element="step", stat="density", common_norm=False, palette='viridis')
+    plt.title('Distribution of Time-to-Empty (TTE) by User Profile', fontweight='bold')
+    plt.xlabel('Time to Empty (Hours)')
+    plt.ylabel('Probability Density')
+    plt.grid(True, alpha=0.3)
+    plt.savefig('output_q2/fig1_tte_distribution.png')
+    plt.show()
+    
+    # --- 图2: 初始电量 vs TTE (非线性关系展示) ---
+    plt.figure(figsize=(10, 6))
+    # 使用散点图，颜色映射温度
+    sc = plt.scatter(df['Initial_SOC']*100, df['TTE_Hours'], c=df['Temperature_C'], cmap='coolwarm', alpha=0.6, s=20)
+    plt.colorbar(sc, label='Temperature (°C)')
+    plt.title('Impact of Initial SOC and Temperature on TTE', fontweight='bold')
+    plt.xlabel('Initial State of Charge (%)')
+    plt.ylabel('Time to Empty (Hours)')
+    plt.grid(True, linestyle='--', alpha=0.5)
+    
+    # 添加拟合曲线 (展示非线性)
+    # 简单的多项式拟合用于视觉引导
+    z = np.polyfit(df['Initial_SOC']*100, df['TTE_Hours'], 2)
+    p = np.poly1d(z)
+    x_range = np.linspace(20, 100, 100)
+    plt.plot(x_range, p(x_range), 'k--', linewidth=2, label='Trend Line')
+    plt.legend()
+    plt.savefig('output_q2/fig2_soc_vs_tte.png')
+    plt.show()
+    
+    # --- 图3: 全局灵敏度分析 (Tornado Plot 替代品) ---
+    # 计算 Spearman 相关系数
+    corr_cols = ['Initial_SOC', 'Temperature_C', 'Signal_dBm', 'Base_Power']
+    corr_matrix = df[corr_cols + ['TTE_Hours']].corr(method='spearman')
+    sensitivity = corr_matrix['TTE_Hours'].drop('TTE_Hours').sort_values()
+    
+    plt.figure(figsize=(10, 5))
+    colors = ['red' if x < 0 else 'green' for x in sensitivity.values]
+    sensitivity.plot(kind='barh', color=colors, alpha=0.8)
+    plt.title('Sensitivity Analysis: Correlation with TTE', fontweight='bold')
+    plt.xlabel('Spearman Correlation Coefficient')
+    plt.axvline(0, color='black', linewidth=0.8)
+    plt.grid(axis='x', linestyle='--', alpha=0.5)
+    
+    # 添加数值标签
+    for index, value in enumerate(sensitivity):
+        plt.text(value, index, f' {value:.2f}', va='center', fontsize=10, fontweight='bold')
+        
+    plt.tight_layout()
+    plt.savefig('output_q2/fig3_sensitivity.png')
+    plt.show()
+
+    print("\nKey Insights:")
+    print(f"1. Base Power Correlation: {sensitivity['Base_Power']:.2f} (Dominant negative factor)")
+    print(f"2. Initial SOC Correlation: {sensitivity['Initial_SOC']:.2f} (Dominant positive factor)")
+    print(f"3. Signal Strength Correlation: {sensitivity['Signal_dBm']:.2f} (Significant environmental factor)")
+
+# ==========================================
+# 5. 主程序入口
+# ==========================================
+if __name__ == "__main__":
+    # 1. 配置环境
+    configure_environment()
+    
+    # 2. 运行模拟
+    # 模拟 3000 次以获得平滑的分布
+    simulation_data = run_monte_carlo_simulation(n_simulations=3000)
+    
+    # 3. 分析与绘图
+    analyze_and_plot(simulation_data)
+    
+    print("\nAll tasks completed. Results saved in 'output_q2' directory.")