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.")