1540 lines
64 KiB
Python
1540 lines
64 KiB
Python
"""
|
||
任务三：月球殖民地水需求灵敏度分析
|
||
|
||
Water Supply Sensitivity Analysis for Moon Colony
|
||
|
||
分析参数:
|
||
1. 水回收效率 (η): 70%-95%
|
||
2. 舒适度因子 (α): 1-300
|
||
3. 医疗紧急参数: 患病率、每次用水量、置信水平
|
||
4. 人口规模 (N): ±20%
|
||
5. 安全缓冲天数: 15-60天
|
||
|
||
输出:
|
||
- 单参数灵敏度分析图
|
||
- 多参数交互分析 (Tornado图)
|
||
- 最坏情况分析 (Stress Test)
|
||
"""
|
||
|
||
import numpy as np
|
||
import matplotlib
|
||
matplotlib.use('Agg')
|
||
import matplotlib.pyplot as plt
|
||
from matplotlib import rcParams
|
||
from scipy import stats
|
||
from dataclasses import dataclass, field
|
||
from typing import List, Dict, Tuple, Optional
|
||
import pandas as pd
|
||
from mpl_toolkits.mplot3d import Axes3D
|
||
import warnings
|
||
warnings.filterwarnings('ignore')
|
||
|
||
# 设置字体
|
||
rcParams['font.sans-serif'] = ['Arial Unicode MS', 'DejaVu Sans', 'SimHei']
|
||
rcParams['axes.unicode_minus'] = False
|
||
plt.style.use('seaborn-v0_8-whitegrid')
|
||
|
||
# ============== 物理常数 ==============
|
||
G0 = 9.81 # m/s²
|
||
OMEGA_EARTH = 7.27e-5 # rad/s
|
||
R_EARTH = 6.371e6 # m
|
||
|
||
# ============== 基准参数 ==============
|
||
@dataclass
|
||
class BaselineParameters:
|
||
"""基准参数配置"""
|
||
# 人口
|
||
population: int = 100_000
|
||
|
||
# 水需求参数
|
||
survival_water: float = 2.5 # L/人/天
|
||
hygiene_water_base: float = 0.4 # L/人/天 (乘以α)
|
||
comfort_factor: float = 50.0 # 舒适度系数
|
||
|
||
# 医疗参数
|
||
medical_water_per_event: float = 5.0 # L/次
|
||
sickness_rate: float = 0.02 # 每日患病率 (2%)
|
||
confidence_level: float = 0.99 # 置信水平
|
||
|
||
# 回收和缓冲
|
||
recycle_rate: float = 0.90 # 水循环回收率
|
||
safety_buffer_days: int = 30 # 安全缓冲天数
|
||
|
||
# 运输参数
|
||
elevator_capacity_per_year: float = 179_000 # 吨/年/部
|
||
num_elevators: int = 3
|
||
elevator_specific_energy: float = 157.2e9 # J/吨
|
||
|
||
@property
|
||
def total_elevator_capacity(self) -> float:
|
||
return self.num_elevators * self.elevator_capacity_per_year
|
||
|
||
@property
|
||
def daily_elevator_capacity(self) -> float:
|
||
return self.total_elevator_capacity / 365
|
||
|
||
|
||
# ============== 水需求计算函数 ==============
|
||
|
||
def calculate_water_demand(
|
||
population: int,
|
||
survival_water: float,
|
||
hygiene_water_base: float,
|
||
comfort_factor: float,
|
||
medical_water_per_event: float,
|
||
sickness_rate: float,
|
||
confidence_level: float,
|
||
recycle_rate: float,
|
||
safety_buffer_days: int
|
||
) -> Dict:
|
||
"""
|
||
计算水需求
|
||
|
||
返回包含各项水需求指标的字典
|
||
"""
|
||
# 每人每日用水
|
||
daily_per_person = survival_water + hygiene_water_base * comfort_factor
|
||
|
||
# 每日循环用水量 (吨)
|
||
daily_circulation = population * daily_per_person / 1000
|
||
|
||
# 医疗用水计算 (使用正态近似)
|
||
# 每日患病人数期望和标准差
|
||
n = population
|
||
p = sickness_rate
|
||
mu = n * p
|
||
sigma = np.sqrt(n * p * (1 - p))
|
||
|
||
# 给定置信水平的Z值
|
||
z = stats.norm.ppf(confidence_level)
|
||
|
||
# 峰值医疗需求 (99%置信)
|
||
peak_medical_events = mu + z * sigma
|
||
daily_medical = peak_medical_events * medical_water_per_event / 1000 # 吨
|
||
|
||
# 平均医疗需求
|
||
mean_medical = mu * medical_water_per_event / 1000 # 吨
|
||
|
||
# 每日循环损耗 (需从地球补充)
|
||
daily_circulation_loss = daily_circulation * (1 - recycle_rate)
|
||
|
||
# 每日总补充量 (损耗 + 医疗)
|
||
daily_supply = daily_circulation_loss + mean_medical
|
||
|
||
# 每月补充量
|
||
monthly_supply = daily_supply * 30
|
||
|
||
# 年补充量
|
||
annual_supply = daily_supply * 365
|
||
|
||
# 水库初始存量 (支撑单日循环)
|
||
reservoir = daily_circulation
|
||
|
||
# 安全缓冲 (峰值医疗 + 损耗)
|
||
daily_peak_supply = daily_circulation_loss + daily_medical
|
||
safety_buffer = daily_peak_supply * safety_buffer_days
|
||
|
||
# 首批运输量 = 水库 + 安全缓冲
|
||
initial_transport = reservoir + safety_buffer
|
||
|
||
return {
|
||
'population': population,
|
||
'comfort_factor': comfort_factor,
|
||
'recycle_rate': recycle_rate,
|
||
'sickness_rate': sickness_rate,
|
||
'confidence_level': confidence_level,
|
||
'safety_buffer_days': safety_buffer_days,
|
||
|
||
'daily_per_person_liters': daily_per_person,
|
||
'daily_circulation_tons': daily_circulation,
|
||
'daily_medical_tons': daily_medical,
|
||
'mean_medical_tons': mean_medical,
|
||
'daily_circulation_loss_tons': daily_circulation_loss,
|
||
'daily_supply_tons': daily_supply,
|
||
'monthly_supply_tons': monthly_supply,
|
||
'annual_supply_tons': annual_supply,
|
||
'reservoir_tons': reservoir,
|
||
'safety_buffer_tons': safety_buffer,
|
||
'initial_transport_tons': initial_transport,
|
||
}
|
||
|
||
|
||
def calculate_transport_metrics(
|
||
water_demand: Dict,
|
||
daily_capacity: float,
|
||
specific_energy: float
|
||
) -> Dict:
|
||
"""计算运输指标"""
|
||
initial_tons = water_demand['initial_transport_tons']
|
||
monthly_tons = water_demand['monthly_supply_tons']
|
||
annual_tons = water_demand['annual_supply_tons']
|
||
|
||
# 首批运输
|
||
initial_days = initial_tons / daily_capacity
|
||
initial_energy_pj = initial_tons * specific_energy / 1e15
|
||
|
||
# 每月补充
|
||
monthly_days = monthly_tons / daily_capacity
|
||
monthly_energy_pj = monthly_tons * specific_energy / 1e15
|
||
|
||
# 年度
|
||
annual_days = annual_tons / daily_capacity
|
||
annual_energy_pj = annual_tons * specific_energy / 1e15
|
||
|
||
return {
|
||
'initial_days': initial_days,
|
||
'initial_energy_pj': initial_energy_pj,
|
||
'monthly_days': monthly_days,
|
||
'monthly_energy_pj': monthly_energy_pj,
|
||
'annual_days': annual_days,
|
||
'annual_energy_pj': annual_energy_pj,
|
||
'annual_tons': annual_tons,
|
||
}
|
||
|
||
|
||
# ============== 灵敏度分析函数 ==============
|
||
|
||
def sensitivity_recycle_rate(baseline: BaselineParameters, save_dir: str):
|
||
"""
|
||
水回收效率灵敏度分析
|
||
η: 70% - 95%
|
||
"""
|
||
print("分析水回收效率灵敏度...")
|
||
|
||
recycle_rates = np.linspace(0.70, 0.95, 26)
|
||
|
||
results = {
|
||
'recycle_rate': [],
|
||
'daily_supply': [],
|
||
'annual_supply': [],
|
||
'annual_energy_pj': [],
|
||
'capacity_ratio': [], # 占电梯年运力比例
|
||
}
|
||
|
||
for eta in recycle_rates:
|
||
demand = calculate_water_demand(
|
||
population=baseline.population,
|
||
survival_water=baseline.survival_water,
|
||
hygiene_water_base=baseline.hygiene_water_base,
|
||
comfort_factor=baseline.comfort_factor,
|
||
medical_water_per_event=baseline.medical_water_per_event,
|
||
sickness_rate=baseline.sickness_rate,
|
||
confidence_level=baseline.confidence_level,
|
||
recycle_rate=eta,
|
||
safety_buffer_days=baseline.safety_buffer_days
|
||
)
|
||
|
||
transport = calculate_transport_metrics(
|
||
demand,
|
||
baseline.daily_elevator_capacity,
|
||
baseline.elevator_specific_energy
|
||
)
|
||
|
||
results['recycle_rate'].append(eta * 100)
|
||
results['daily_supply'].append(demand['daily_supply_tons'])
|
||
results['annual_supply'].append(demand['annual_supply_tons'])
|
||
results['annual_energy_pj'].append(transport['annual_energy_pj'])
|
||
results['capacity_ratio'].append(demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100)
|
||
|
||
# 绘图
|
||
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
|
||
|
||
# 图1: 日补充量
|
||
ax = axes[0, 0]
|
||
ax.plot(results['recycle_rate'], results['daily_supply'], 'b-', linewidth=2, marker='o', markersize=4)
|
||
ax.axvline(x=90, color='r', linestyle='--', label='Baseline (90%)')
|
||
ax.set_xlabel('Recycle Rate (%)', fontsize=12)
|
||
ax.set_ylabel('Daily Supply (tons)', fontsize=12)
|
||
ax.set_title('Daily Water Supply vs Recycle Rate', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
# 图2: 年补充量
|
||
ax = axes[0, 1]
|
||
ax.plot(results['recycle_rate'], results['annual_supply'], 'g-', linewidth=2, marker='s', markersize=4)
|
||
ax.axvline(x=90, color='r', linestyle='--', label='Baseline (90%)')
|
||
ax.set_xlabel('Recycle Rate (%)', fontsize=12)
|
||
ax.set_ylabel('Annual Supply (tons)', fontsize=12)
|
||
ax.set_title('Annual Water Supply vs Recycle Rate', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
# 图3: 年能耗
|
||
ax = axes[1, 0]
|
||
ax.plot(results['recycle_rate'], results['annual_energy_pj'], 'purple', linewidth=2, marker='^', markersize=4)
|
||
ax.axvline(x=90, color='r', linestyle='--', label='Baseline (90%)')
|
||
ax.set_xlabel('Recycle Rate (%)', fontsize=12)
|
||
ax.set_ylabel('Annual Energy (PJ)', fontsize=12)
|
||
ax.set_title('Annual Transport Energy vs Recycle Rate', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
# 图4: 运力占比
|
||
ax = axes[1, 1]
|
||
ax.plot(results['recycle_rate'], results['capacity_ratio'], 'orange', linewidth=2, marker='d', markersize=4)
|
||
ax.axvline(x=90, color='r', linestyle='--', label='Baseline (90%)')
|
||
ax.axhline(y=100, color='gray', linestyle=':', label='Full Capacity')
|
||
ax.set_xlabel('Recycle Rate (%)', fontsize=12)
|
||
ax.set_ylabel('Elevator Capacity Usage (%)', fontsize=12)
|
||
ax.set_title('Elevator Capacity Utilization vs Recycle Rate', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
plt.suptitle(f'Sensitivity Analysis: Water Recycle Rate\n(α={baseline.comfort_factor}, N={baseline.population:,})',
|
||
fontsize=14, y=1.02)
|
||
plt.tight_layout()
|
||
plt.savefig(f'{save_dir}/sensitivity_recycle_rate.png', dpi=150, bbox_inches='tight')
|
||
print(f" 保存: {save_dir}/sensitivity_recycle_rate.png")
|
||
|
||
return pd.DataFrame(results)
|
||
|
||
|
||
def sensitivity_comfort_factor(baseline: BaselineParameters, save_dir: str):
|
||
"""
|
||
舒适度因子连续灵敏度分析
|
||
α: 1 - 300
|
||
"""
|
||
print("分析舒适度因子灵敏度...")
|
||
|
||
alphas = np.concatenate([
|
||
np.linspace(1, 10, 10),
|
||
np.linspace(15, 50, 8),
|
||
np.linspace(60, 150, 10),
|
||
np.linspace(175, 300, 6)
|
||
])
|
||
|
||
results = {
|
||
'alpha': [],
|
||
'daily_per_person': [],
|
||
'daily_supply': [],
|
||
'annual_supply': [],
|
||
'annual_energy_pj': [],
|
||
'initial_transport': [],
|
||
'capacity_ratio': [],
|
||
}
|
||
|
||
for alpha in alphas:
|
||
demand = calculate_water_demand(
|
||
population=baseline.population,
|
||
survival_water=baseline.survival_water,
|
||
hygiene_water_base=baseline.hygiene_water_base,
|
||
comfort_factor=alpha,
|
||
medical_water_per_event=baseline.medical_water_per_event,
|
||
sickness_rate=baseline.sickness_rate,
|
||
confidence_level=baseline.confidence_level,
|
||
recycle_rate=baseline.recycle_rate,
|
||
safety_buffer_days=baseline.safety_buffer_days
|
||
)
|
||
|
||
transport = calculate_transport_metrics(
|
||
demand,
|
||
baseline.daily_elevator_capacity,
|
||
baseline.elevator_specific_energy
|
||
)
|
||
|
||
results['alpha'].append(alpha)
|
||
results['daily_per_person'].append(demand['daily_per_person_liters'])
|
||
results['daily_supply'].append(demand['daily_supply_tons'])
|
||
results['annual_supply'].append(demand['annual_supply_tons'])
|
||
results['annual_energy_pj'].append(transport['annual_energy_pj'])
|
||
results['initial_transport'].append(demand['initial_transport_tons'])
|
||
results['capacity_ratio'].append(demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100)
|
||
|
||
# 绘图
|
||
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
|
||
|
||
# 图1: 人均日用水
|
||
ax = axes[0, 0]
|
||
ax.plot(results['alpha'], results['daily_per_person'], 'b-', linewidth=2)
|
||
ax.axvline(x=1, color='green', linestyle='--', alpha=0.7, label='Survival (α=1)')
|
||
ax.axvline(x=50, color='orange', linestyle='--', alpha=0.7, label='Comfort (α=50)')
|
||
ax.axvline(x=250, color='red', linestyle='--', alpha=0.7, label='Luxury (α=250)')
|
||
ax.set_xlabel('Comfort Factor (α)', fontsize=12)
|
||
ax.set_ylabel('Daily Water per Person (L)', fontsize=12)
|
||
ax.set_title('Per Capita Water Usage vs Comfort Factor', fontsize=13)
|
||
ax.legend(fontsize=9)
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
# 图2: 日补充量
|
||
ax = axes[0, 1]
|
||
ax.plot(results['alpha'], results['daily_supply'], 'g-', linewidth=2)
|
||
ax.axvline(x=50, color='r', linestyle='--', label='Baseline (α=50)')
|
||
ax.set_xlabel('Comfort Factor (α)', fontsize=12)
|
||
ax.set_ylabel('Daily Supply (tons)', fontsize=12)
|
||
ax.set_title('Daily Water Supply vs Comfort Factor', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
# 图3: 年补充量
|
||
ax = axes[0, 2]
|
||
ax.plot(results['alpha'], results['annual_supply'], 'purple', linewidth=2)
|
||
ax.axvline(x=50, color='r', linestyle='--', label='Baseline (α=50)')
|
||
ax.set_xlabel('Comfort Factor (α)', fontsize=12)
|
||
ax.set_ylabel('Annual Supply (tons)', fontsize=12)
|
||
ax.set_title('Annual Water Supply vs Comfort Factor', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
# 图4: 年能耗
|
||
ax = axes[1, 0]
|
||
ax.plot(results['alpha'], results['annual_energy_pj'], 'orange', linewidth=2)
|
||
ax.axvline(x=50, color='r', linestyle='--', label='Baseline (α=50)')
|
||
ax.set_xlabel('Comfort Factor (α)', fontsize=12)
|
||
ax.set_ylabel('Annual Energy (PJ)', fontsize=12)
|
||
ax.set_title('Annual Transport Energy vs Comfort Factor', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
# 图5: 首批运输量
|
||
ax = axes[1, 1]
|
||
ax.plot(results['alpha'], results['initial_transport'], 'teal', linewidth=2)
|
||
ax.axvline(x=50, color='r', linestyle='--', label='Baseline (α=50)')
|
||
ax.set_xlabel('Comfort Factor (α)', fontsize=12)
|
||
ax.set_ylabel('Initial Transport (tons)', fontsize=12)
|
||
ax.set_title('Initial Transport Volume vs Comfort Factor', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
# 图6: 运力占比
|
||
ax = axes[1, 2]
|
||
ax.plot(results['alpha'], results['capacity_ratio'], 'brown', linewidth=2)
|
||
ax.axvline(x=50, color='r', linestyle='--', label='Baseline (α=50)')
|
||
ax.axhline(y=100, color='gray', linestyle=':', label='Full Capacity')
|
||
ax.fill_between(results['alpha'], 0, results['capacity_ratio'], alpha=0.3)
|
||
ax.set_xlabel('Comfort Factor (α)', fontsize=12)
|
||
ax.set_ylabel('Elevator Capacity Usage (%)', fontsize=12)
|
||
ax.set_title('Elevator Capacity Utilization vs Comfort Factor', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
plt.suptitle(f'Sensitivity Analysis: Comfort Factor (α)\n(η={baseline.recycle_rate*100:.0f}%, N={baseline.population:,})',
|
||
fontsize=14, y=1.02)
|
||
plt.tight_layout()
|
||
plt.savefig(f'{save_dir}/sensitivity_comfort_factor.png', dpi=150, bbox_inches='tight')
|
||
print(f" 保存: {save_dir}/sensitivity_comfort_factor.png")
|
||
|
||
return pd.DataFrame(results)
|
||
|
||
|
||
def sensitivity_medical_params(baseline: BaselineParameters, save_dir: str):
|
||
"""
|
||
医疗紧急参数灵敏度分析
|
||
- 患病率: 1% - 5%
|
||
- 每次医疗用水: 3 - 15 L
|
||
- 置信水平: 95% - 99.9%
|
||
"""
|
||
print("分析医疗参数灵敏度...")
|
||
|
||
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
|
||
|
||
# ========== 患病率分析 ==========
|
||
sickness_rates = np.linspace(0.01, 0.05, 21)
|
||
results_sr = {'sickness_rate': [], 'daily_medical': [], 'annual_supply': [], 'safety_buffer': []}
|
||
|
||
for sr in sickness_rates:
|
||
demand = calculate_water_demand(
|
||
population=baseline.population,
|
||
survival_water=baseline.survival_water,
|
||
hygiene_water_base=baseline.hygiene_water_base,
|
||
comfort_factor=baseline.comfort_factor,
|
||
medical_water_per_event=baseline.medical_water_per_event,
|
||
sickness_rate=sr,
|
||
confidence_level=baseline.confidence_level,
|
||
recycle_rate=baseline.recycle_rate,
|
||
safety_buffer_days=baseline.safety_buffer_days
|
||
)
|
||
results_sr['sickness_rate'].append(sr * 100)
|
||
results_sr['daily_medical'].append(demand['daily_medical_tons'])
|
||
results_sr['annual_supply'].append(demand['annual_supply_tons'])
|
||
results_sr['safety_buffer'].append(demand['safety_buffer_tons'])
|
||
|
||
ax = axes[0, 0]
|
||
ax.plot(results_sr['sickness_rate'], results_sr['daily_medical'], 'b-', linewidth=2, marker='o', markersize=4)
|
||
ax.axvline(x=2, color='r', linestyle='--', label='Baseline (2%)')
|
||
ax.set_xlabel('Daily Sickness Rate (%)', fontsize=12)
|
||
ax.set_ylabel('Peak Daily Medical Water (tons)', fontsize=12)
|
||
ax.set_title('Peak Medical Water Demand vs Sickness Rate', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
ax = axes[0, 1]
|
||
ax.plot(results_sr['sickness_rate'], results_sr['safety_buffer'], 'g-', linewidth=2, marker='s', markersize=4)
|
||
ax.axvline(x=2, color='r', linestyle='--', label='Baseline (2%)')
|
||
ax.set_xlabel('Daily Sickness Rate (%)', fontsize=12)
|
||
ax.set_ylabel('Safety Buffer (tons)', fontsize=12)
|
||
ax.set_title('Safety Buffer vs Sickness Rate', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
# ========== 每次医疗用水分析 ==========
|
||
medical_waters = np.linspace(3, 15, 13)
|
||
results_mw = {'medical_water': [], 'daily_medical': [], 'annual_supply': []}
|
||
|
||
for mw in medical_waters:
|
||
demand = calculate_water_demand(
|
||
population=baseline.population,
|
||
survival_water=baseline.survival_water,
|
||
hygiene_water_base=baseline.hygiene_water_base,
|
||
comfort_factor=baseline.comfort_factor,
|
||
medical_water_per_event=mw,
|
||
sickness_rate=baseline.sickness_rate,
|
||
confidence_level=baseline.confidence_level,
|
||
recycle_rate=baseline.recycle_rate,
|
||
safety_buffer_days=baseline.safety_buffer_days
|
||
)
|
||
results_mw['medical_water'].append(mw)
|
||
results_mw['daily_medical'].append(demand['daily_medical_tons'])
|
||
results_mw['annual_supply'].append(demand['annual_supply_tons'])
|
||
|
||
ax = axes[0, 2]
|
||
ax.plot(results_mw['medical_water'], results_mw['daily_medical'], 'purple', linewidth=2, marker='^', markersize=5)
|
||
ax.axvline(x=5, color='r', linestyle='--', label='Baseline (5 L)')
|
||
ax.set_xlabel('Medical Water per Event (L)', fontsize=12)
|
||
ax.set_ylabel('Peak Daily Medical Water (tons)', fontsize=12)
|
||
ax.set_title('Peak Medical Water vs Per-Event Usage', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
# ========== 置信水平分析 ==========
|
||
confidence_levels = np.linspace(0.90, 0.999, 20)
|
||
results_cl = {'confidence': [], 'daily_medical': [], 'safety_buffer': [], 'z_value': []}
|
||
|
||
for cl in confidence_levels:
|
||
demand = calculate_water_demand(
|
||
population=baseline.population,
|
||
survival_water=baseline.survival_water,
|
||
hygiene_water_base=baseline.hygiene_water_base,
|
||
comfort_factor=baseline.comfort_factor,
|
||
medical_water_per_event=baseline.medical_water_per_event,
|
||
sickness_rate=baseline.sickness_rate,
|
||
confidence_level=cl,
|
||
recycle_rate=baseline.recycle_rate,
|
||
safety_buffer_days=baseline.safety_buffer_days
|
||
)
|
||
results_cl['confidence'].append(cl * 100)
|
||
results_cl['daily_medical'].append(demand['daily_medical_tons'])
|
||
results_cl['safety_buffer'].append(demand['safety_buffer_tons'])
|
||
results_cl['z_value'].append(stats.norm.ppf(cl))
|
||
|
||
ax = axes[1, 0]
|
||
ax.plot(results_cl['confidence'], results_cl['daily_medical'], 'orange', linewidth=2, marker='d', markersize=4)
|
||
ax.axvline(x=99, color='r', linestyle='--', label='Baseline (99%)')
|
||
ax.set_xlabel('Confidence Level (%)', fontsize=12)
|
||
ax.set_ylabel('Peak Daily Medical Water (tons)', fontsize=12)
|
||
ax.set_title('Peak Medical Water vs Confidence Level', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
ax = axes[1, 1]
|
||
ax.plot(results_cl['confidence'], results_cl['z_value'], 'teal', linewidth=2, marker='o', markersize=4)
|
||
ax.axvline(x=99, color='r', linestyle='--', label='Baseline (99%)')
|
||
ax.set_xlabel('Confidence Level (%)', fontsize=12)
|
||
ax.set_ylabel('Z-Value', fontsize=12)
|
||
ax.set_title('Z-Value vs Confidence Level', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
# ========== 综合影响图 ==========
|
||
ax = axes[1, 2]
|
||
|
||
# 计算各参数变化对年供水的影响
|
||
baseline_demand = calculate_water_demand(
|
||
population=baseline.population,
|
||
survival_water=baseline.survival_water,
|
||
hygiene_water_base=baseline.hygiene_water_base,
|
||
comfort_factor=baseline.comfort_factor,
|
||
medical_water_per_event=baseline.medical_water_per_event,
|
||
sickness_rate=baseline.sickness_rate,
|
||
confidence_level=baseline.confidence_level,
|
||
recycle_rate=baseline.recycle_rate,
|
||
safety_buffer_days=baseline.safety_buffer_days
|
||
)
|
||
baseline_annual = baseline_demand['annual_supply_tons']
|
||
|
||
# 患病率影响 (1% -> 5%)
|
||
sr_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, baseline.medical_water_per_event, 0.01,
|
||
baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days)
|
||
sr_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, baseline.medical_water_per_event, 0.05,
|
||
baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days)
|
||
|
||
# 医疗用水影响 (3L -> 15L)
|
||
mw_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, 3, baseline.sickness_rate,
|
||
baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days)
|
||
mw_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, 15, baseline.sickness_rate,
|
||
baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days)
|
||
|
||
# 置信水平影响 (95% -> 99.9%)
|
||
cl_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate,
|
||
0.95, baseline.recycle_rate, baseline.safety_buffer_days)
|
||
cl_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate,
|
||
0.999, baseline.recycle_rate, baseline.safety_buffer_days)
|
||
|
||
params = ['Sickness Rate\n(1%→5%)', 'Medical Water\n(3L→15L)', 'Confidence\n(95%→99.9%)']
|
||
changes_low = [
|
||
(sr_low['annual_supply_tons'] - baseline_annual) / baseline_annual * 100,
|
||
(mw_low['annual_supply_tons'] - baseline_annual) / baseline_annual * 100,
|
||
(cl_low['annual_supply_tons'] - baseline_annual) / baseline_annual * 100,
|
||
]
|
||
changes_high = [
|
||
(sr_high['annual_supply_tons'] - baseline_annual) / baseline_annual * 100,
|
||
(mw_high['annual_supply_tons'] - baseline_annual) / baseline_annual * 100,
|
||
(cl_high['annual_supply_tons'] - baseline_annual) / baseline_annual * 100,
|
||
]
|
||
|
||
x = np.arange(len(params))
|
||
width = 0.35
|
||
ax.barh(x - width/2, changes_low, width, label='Low Value', color='green', alpha=0.7)
|
||
ax.barh(x + width/2, changes_high, width, label='High Value', color='red', alpha=0.7)
|
||
ax.axvline(x=0, color='black', linewidth=1)
|
||
ax.set_yticks(x)
|
||
ax.set_yticklabels(params)
|
||
ax.set_xlabel('Change in Annual Supply (%)', fontsize=12)
|
||
ax.set_title('Medical Parameters Impact on Annual Supply', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3, axis='x')
|
||
|
||
plt.suptitle(f'Sensitivity Analysis: Medical Emergency Parameters\n(α={baseline.comfort_factor}, η={baseline.recycle_rate*100:.0f}%)',
|
||
fontsize=14, y=1.02)
|
||
plt.tight_layout()
|
||
plt.savefig(f'{save_dir}/sensitivity_medical_params.png', dpi=150, bbox_inches='tight')
|
||
print(f" 保存: {save_dir}/sensitivity_medical_params.png")
|
||
|
||
return results_sr, results_mw, results_cl
|
||
|
||
|
||
def sensitivity_population(baseline: BaselineParameters, save_dir: str):
|
||
"""
|
||
人口规模灵敏度分析
|
||
N: 80,000 - 120,000 (±20%)
|
||
"""
|
||
print("分析人口规模灵敏度...")
|
||
|
||
populations = np.linspace(80000, 120000, 21).astype(int)
|
||
|
||
results = {
|
||
'population': [],
|
||
'daily_circulation': [],
|
||
'daily_supply': [],
|
||
'annual_supply': [],
|
||
'annual_energy_pj': [],
|
||
'initial_transport': [],
|
||
'capacity_ratio': [],
|
||
}
|
||
|
||
for N in populations:
|
||
demand = calculate_water_demand(
|
||
population=N,
|
||
survival_water=baseline.survival_water,
|
||
hygiene_water_base=baseline.hygiene_water_base,
|
||
comfort_factor=baseline.comfort_factor,
|
||
medical_water_per_event=baseline.medical_water_per_event,
|
||
sickness_rate=baseline.sickness_rate,
|
||
confidence_level=baseline.confidence_level,
|
||
recycle_rate=baseline.recycle_rate,
|
||
safety_buffer_days=baseline.safety_buffer_days
|
||
)
|
||
|
||
transport = calculate_transport_metrics(
|
||
demand,
|
||
baseline.daily_elevator_capacity,
|
||
baseline.elevator_specific_energy
|
||
)
|
||
|
||
results['population'].append(N / 1000) # 千人
|
||
results['daily_circulation'].append(demand['daily_circulation_tons'])
|
||
results['daily_supply'].append(demand['daily_supply_tons'])
|
||
results['annual_supply'].append(demand['annual_supply_tons'])
|
||
results['annual_energy_pj'].append(transport['annual_energy_pj'])
|
||
results['initial_transport'].append(demand['initial_transport_tons'])
|
||
results['capacity_ratio'].append(demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100)
|
||
|
||
# 绘图
|
||
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
|
||
|
||
# 图1: 日循环用水
|
||
ax = axes[0, 0]
|
||
ax.plot(results['population'], results['daily_circulation'], 'b-', linewidth=2, marker='o', markersize=4)
|
||
ax.axvline(x=100, color='r', linestyle='--', label='Baseline (100k)')
|
||
ax.fill_between(results['population'], 0, results['daily_circulation'], alpha=0.2)
|
||
ax.set_xlabel('Population (thousands)', fontsize=12)
|
||
ax.set_ylabel('Daily Circulation Water (tons)', fontsize=12)
|
||
ax.set_title('Daily Circulation Water vs Population', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
# 图2: 日补充量
|
||
ax = axes[0, 1]
|
||
ax.plot(results['population'], results['daily_supply'], 'g-', linewidth=2, marker='s', markersize=4)
|
||
ax.axvline(x=100, color='r', linestyle='--', label='Baseline (100k)')
|
||
ax.set_xlabel('Population (thousands)', fontsize=12)
|
||
ax.set_ylabel('Daily Supply (tons)', fontsize=12)
|
||
ax.set_title('Daily Water Supply vs Population', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
# 图3: 年能耗
|
||
ax = axes[1, 0]
|
||
ax.plot(results['population'], results['annual_energy_pj'], 'purple', linewidth=2, marker='^', markersize=4)
|
||
ax.axvline(x=100, color='r', linestyle='--', label='Baseline (100k)')
|
||
ax.set_xlabel('Population (thousands)', fontsize=12)
|
||
ax.set_ylabel('Annual Energy (PJ)', fontsize=12)
|
||
ax.set_title('Annual Transport Energy vs Population', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
# 图4: 运力占比
|
||
ax = axes[1, 1]
|
||
ax.plot(results['population'], results['capacity_ratio'], 'orange', linewidth=2, marker='d', markersize=4)
|
||
ax.axvline(x=100, color='r', linestyle='--', label='Baseline (100k)')
|
||
ax.axhline(y=100, color='gray', linestyle=':', label='Full Capacity')
|
||
ax.set_xlabel('Population (thousands)', fontsize=12)
|
||
ax.set_ylabel('Elevator Capacity Usage (%)', fontsize=12)
|
||
ax.set_title('Elevator Capacity Utilization vs Population', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
plt.suptitle(f'Sensitivity Analysis: Population Size\n(α={baseline.comfort_factor}, η={baseline.recycle_rate*100:.0f}%)',
|
||
fontsize=14, y=1.02)
|
||
plt.tight_layout()
|
||
plt.savefig(f'{save_dir}/sensitivity_population.png', dpi=150, bbox_inches='tight')
|
||
print(f" 保存: {save_dir}/sensitivity_population.png")
|
||
|
||
return pd.DataFrame(results)
|
||
|
||
|
||
def sensitivity_buffer_days(baseline: BaselineParameters, save_dir: str):
|
||
"""
|
||
安全缓冲天数灵敏度分析
|
||
Days: 15 - 60
|
||
"""
|
||
print("分析安全缓冲天数灵敏度...")
|
||
|
||
buffer_days = np.arange(15, 61, 3)
|
||
|
||
results = {
|
||
'buffer_days': [],
|
||
'safety_buffer_tons': [],
|
||
'initial_transport': [],
|
||
'initial_days': [],
|
||
'initial_energy_pj': [],
|
||
}
|
||
|
||
for days in buffer_days:
|
||
demand = calculate_water_demand(
|
||
population=baseline.population,
|
||
survival_water=baseline.survival_water,
|
||
hygiene_water_base=baseline.hygiene_water_base,
|
||
comfort_factor=baseline.comfort_factor,
|
||
medical_water_per_event=baseline.medical_water_per_event,
|
||
sickness_rate=baseline.sickness_rate,
|
||
confidence_level=baseline.confidence_level,
|
||
recycle_rate=baseline.recycle_rate,
|
||
safety_buffer_days=days
|
||
)
|
||
|
||
transport = calculate_transport_metrics(
|
||
demand,
|
||
baseline.daily_elevator_capacity,
|
||
baseline.elevator_specific_energy
|
||
)
|
||
|
||
results['buffer_days'].append(days)
|
||
results['safety_buffer_tons'].append(demand['safety_buffer_tons'])
|
||
results['initial_transport'].append(demand['initial_transport_tons'])
|
||
results['initial_days'].append(transport['initial_days'])
|
||
results['initial_energy_pj'].append(transport['initial_energy_pj'])
|
||
|
||
# 绘图
|
||
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
|
||
|
||
# 图1: 安全缓冲量
|
||
ax = axes[0, 0]
|
||
ax.plot(results['buffer_days'], results['safety_buffer_tons'], 'b-', linewidth=2, marker='o', markersize=5)
|
||
ax.axvline(x=30, color='r', linestyle='--', label='Baseline (30 days)')
|
||
ax.set_xlabel('Safety Buffer Days', fontsize=12)
|
||
ax.set_ylabel('Safety Buffer (tons)', fontsize=12)
|
||
ax.set_title('Safety Buffer Volume vs Buffer Days', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
# 图2: 首批运输量
|
||
ax = axes[0, 1]
|
||
ax.plot(results['buffer_days'], results['initial_transport'], 'g-', linewidth=2, marker='s', markersize=5)
|
||
ax.axvline(x=30, color='r', linestyle='--', label='Baseline (30 days)')
|
||
ax.set_xlabel('Safety Buffer Days', fontsize=12)
|
||
ax.set_ylabel('Initial Transport (tons)', fontsize=12)
|
||
ax.set_title('Initial Transport Volume vs Buffer Days', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
# 图3: 首批运输天数
|
||
ax = axes[1, 0]
|
||
ax.plot(results['buffer_days'], results['initial_days'], 'purple', linewidth=2, marker='^', markersize=5)
|
||
ax.axvline(x=30, color='r', linestyle='--', label='Baseline (30 days)')
|
||
ax.set_xlabel('Safety Buffer Days', fontsize=12)
|
||
ax.set_ylabel('Initial Transport Time (days)', fontsize=12)
|
||
ax.set_title('Initial Transport Duration vs Buffer Days', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
# 图4: 首批运输能耗
|
||
ax = axes[1, 1]
|
||
ax.plot(results['buffer_days'], results['initial_energy_pj'], 'orange', linewidth=2, marker='d', markersize=5)
|
||
ax.axvline(x=30, color='r', linestyle='--', label='Baseline (30 days)')
|
||
ax.set_xlabel('Safety Buffer Days', fontsize=12)
|
||
ax.set_ylabel('Initial Transport Energy (PJ)', fontsize=12)
|
||
ax.set_title('Initial Transport Energy vs Buffer Days', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
plt.suptitle(f'Sensitivity Analysis: Safety Buffer Days\n(α={baseline.comfort_factor}, η={baseline.recycle_rate*100:.0f}%)',
|
||
fontsize=14, y=1.02)
|
||
plt.tight_layout()
|
||
plt.savefig(f'{save_dir}/sensitivity_buffer_days.png', dpi=150, bbox_inches='tight')
|
||
print(f" 保存: {save_dir}/sensitivity_buffer_days.png")
|
||
|
||
return pd.DataFrame(results)
|
||
|
||
|
||
# ============== 多参数交互分析 ==============
|
||
|
||
def tornado_analysis(baseline: BaselineParameters, save_dir: str):
|
||
"""
|
||
Tornado图 - 多参数灵敏度排序分析（马卡龙色系简约版）
|
||
"""
|
||
print("生成Tornado图...")
|
||
|
||
# 基准值
|
||
baseline_demand = calculate_water_demand(
|
||
population=baseline.population,
|
||
survival_water=baseline.survival_water,
|
||
hygiene_water_base=baseline.hygiene_water_base,
|
||
comfort_factor=baseline.comfort_factor,
|
||
medical_water_per_event=baseline.medical_water_per_event,
|
||
sickness_rate=baseline.sickness_rate,
|
||
confidence_level=baseline.confidence_level,
|
||
recycle_rate=baseline.recycle_rate,
|
||
safety_buffer_days=baseline.safety_buffer_days
|
||
)
|
||
baseline_annual = baseline_demand['annual_supply_tons']
|
||
baseline_initial = baseline_demand['initial_transport_tons']
|
||
|
||
# 定义参数范围 - 根据论述调整（删除Confidence Level）
|
||
params = [
|
||
('Recycle Rate (η)', 'recycle_rate', 0.80, 0.95, baseline.recycle_rate),
|
||
('Comfort Factor (α)', 'comfort_factor', 25, 100, baseline.comfort_factor),
|
||
('Population (N)', 'population', 80000, 120000, baseline.population),
|
||
('Sickness Rate', 'sickness_rate', 0.01, 0.04, baseline.sickness_rate),
|
||
('Medical Water/Event', 'medical_water', 3, 10, baseline.medical_water_per_event),
|
||
('Buffer Days', 'buffer_days', 15, 45, baseline.safety_buffer_days),
|
||
]
|
||
|
||
results_annual = []
|
||
results_initial = []
|
||
|
||
for name, param_type, low, high, base_val in params:
|
||
# 低值计算
|
||
if param_type == 'recycle_rate':
|
||
demand_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate,
|
||
baseline.confidence_level, low, baseline.safety_buffer_days)
|
||
demand_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate,
|
||
baseline.confidence_level, high, baseline.safety_buffer_days)
|
||
elif param_type == 'comfort_factor':
|
||
demand_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
low, baseline.medical_water_per_event, baseline.sickness_rate,
|
||
baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days)
|
||
demand_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
high, baseline.medical_water_per_event, baseline.sickness_rate,
|
||
baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days)
|
||
elif param_type == 'population':
|
||
demand_low = calculate_water_demand(int(low), baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate,
|
||
baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days)
|
||
demand_high = calculate_water_demand(int(high), baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate,
|
||
baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days)
|
||
elif param_type == 'sickness_rate':
|
||
demand_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, baseline.medical_water_per_event, low,
|
||
baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days)
|
||
demand_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, baseline.medical_water_per_event, high,
|
||
baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days)
|
||
elif param_type == 'medical_water':
|
||
demand_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, low, baseline.sickness_rate,
|
||
baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days)
|
||
demand_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, high, baseline.sickness_rate,
|
||
baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days)
|
||
elif param_type == 'buffer_days':
|
||
demand_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate,
|
||
baseline.confidence_level, baseline.recycle_rate, int(low))
|
||
demand_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate,
|
||
baseline.confidence_level, baseline.recycle_rate, int(high))
|
||
|
||
# 计算变化百分比
|
||
change_low_annual = (demand_low['annual_supply_tons'] - baseline_annual) / baseline_annual * 100
|
||
change_high_annual = (demand_high['annual_supply_tons'] - baseline_annual) / baseline_annual * 100
|
||
change_low_initial = (demand_low['initial_transport_tons'] - baseline_initial) / baseline_initial * 100
|
||
change_high_initial = (demand_high['initial_transport_tons'] - baseline_initial) / baseline_initial * 100
|
||
|
||
# 确保排序一致性 (swing = |high - low|)
|
||
swing_annual = abs(change_high_annual - change_low_annual)
|
||
swing_initial = abs(change_high_initial - change_low_initial)
|
||
|
||
results_annual.append({
|
||
'param': name,
|
||
'low': min(change_low_annual, change_high_annual),
|
||
'high': max(change_low_annual, change_high_annual),
|
||
'swing': swing_annual
|
||
})
|
||
results_initial.append({
|
||
'param': name,
|
||
'low': min(change_low_initial, change_high_initial),
|
||
'high': max(change_low_initial, change_high_initial),
|
||
'swing': swing_initial
|
||
})
|
||
|
||
# 按年度需求的swing排序（统一排序顺序）
|
||
results_annual = sorted(results_annual, key=lambda x: x['swing'], reverse=True)
|
||
|
||
# 按相同顺序重排initial结果
|
||
param_order = [r['param'] for r in results_annual]
|
||
results_initial_dict = {r['param']: r for r in results_initial}
|
||
results_initial = [results_initial_dict[p] for p in param_order]
|
||
|
||
# 马卡龙色系
|
||
color_decrease = '#A8D8B9' # 薄荷绿 - 减少
|
||
color_increase = '#F7C5A8' # 杏橙 - 增加
|
||
color_decrease_dark = '#7BC08F'
|
||
color_increase_dark = '#E8A07A'
|
||
|
||
# 绘图 - 简约风格（缩小尺寸）
|
||
fig, axes = plt.subplots(1, 2, figsize=(9, 3.5))
|
||
fig.patch.set_facecolor('#FAFAFA')
|
||
|
||
# 统一的参数标签
|
||
params_sorted = [r['param'] for r in results_annual]
|
||
y_pos = np.arange(len(params_sorted))
|
||
|
||
# ========== 图1: 年补充量Tornado ==========
|
||
ax = axes[0]
|
||
ax.set_facecolor('#FAFAFA')
|
||
lows = [r['low'] for r in results_annual]
|
||
highs = [r['high'] for r in results_annual]
|
||
|
||
# 绘制条形图 - 分开绘制正负部分
|
||
for i, (l, h) in enumerate(zip(lows, highs)):
|
||
# 负值部分（减少）
|
||
if l < 0:
|
||
ax.barh(i, l, color=color_decrease, edgecolor=color_decrease_dark,
|
||
linewidth=1, height=0.65)
|
||
# 正值部分（增加）
|
||
if h > 0:
|
||
ax.barh(i, h, color=color_increase, edgecolor=color_increase_dark,
|
||
linewidth=1, height=0.65)
|
||
|
||
# 添加数值标注
|
||
for i, (l, h) in enumerate(zip(lows, highs)):
|
||
if l < 0:
|
||
ax.text(l - 2, i, f'{l:.0f}%', ha='right', va='center', fontsize=8, color='#555555')
|
||
if h > 0:
|
||
ax.text(h + 2, i, f'+{h:.0f}%', ha='left', va='center', fontsize=8, color='#555555')
|
||
|
||
ax.axvline(x=0, color='#333333', linewidth=1.2)
|
||
ax.set_yticks(y_pos)
|
||
ax.set_yticklabels(params_sorted, fontsize=10)
|
||
ax.set_xlabel('Change in Annual Supply (%)', fontsize=11, color='#444444')
|
||
ax.spines['top'].set_visible(False)
|
||
ax.spines['right'].set_visible(False)
|
||
ax.spines['left'].set_color('#CCCCCC')
|
||
ax.spines['bottom'].set_color('#CCCCCC')
|
||
ax.grid(True, alpha=0.15, axis='x', linestyle='--')
|
||
ax.set_xlim(-65, 115)
|
||
|
||
# ========== 图2: 首批运输量Tornado ==========
|
||
ax = axes[1]
|
||
ax.set_facecolor('#FAFAFA')
|
||
lows = [r['low'] for r in results_initial]
|
||
highs = [r['high'] for r in results_initial]
|
||
|
||
for i, (l, h) in enumerate(zip(lows, highs)):
|
||
# 负值部分（减少）
|
||
if l < 0:
|
||
ax.barh(i, l, color=color_decrease, edgecolor=color_decrease_dark,
|
||
linewidth=1, height=0.65)
|
||
# 正值部分（增加）
|
||
if h > 0:
|
||
ax.barh(i, h, color=color_increase, edgecolor=color_increase_dark,
|
||
linewidth=1, height=0.65)
|
||
|
||
# 添加数值标注
|
||
for i, (l, h) in enumerate(zip(lows, highs)):
|
||
if l < 0:
|
||
ax.text(l - 2, i, f'{l:.0f}%', ha='right', va='center', fontsize=8, color='#555555')
|
||
if h > 0:
|
||
ax.text(h + 2, i, f'+{h:.0f}%', ha='left', va='center', fontsize=8, color='#555555')
|
||
|
||
ax.axvline(x=0, color='#333333', linewidth=1.2)
|
||
ax.set_yticks(y_pos)
|
||
ax.set_yticklabels([]) # 删除右图纵坐标标签
|
||
ax.set_xlabel('Change in Initial Transport (%)', fontsize=11, color='#444444')
|
||
ax.spines['top'].set_visible(False)
|
||
ax.spines['right'].set_visible(False)
|
||
ax.spines['left'].set_visible(False) # 隐藏左边框
|
||
ax.spines['bottom'].set_color('#CCCCCC')
|
||
ax.grid(True, alpha=0.15, axis='x', linestyle='--')
|
||
ax.set_xlim(-65, 105)
|
||
|
||
# 简化图例 - 放在图的下方中间
|
||
from matplotlib.patches import Patch
|
||
legend_elements = [
|
||
Patch(facecolor=color_decrease, edgecolor=color_decrease_dark, label='Decrease'),
|
||
Patch(facecolor=color_increase, edgecolor=color_increase_dark, label='Increase'),
|
||
]
|
||
fig.legend(handles=legend_elements, loc='lower center', ncol=2, fontsize=9,
|
||
framealpha=0.9, edgecolor='#CCCCCC', bbox_to_anchor=(0.55, -0.02))
|
||
|
||
plt.tight_layout()
|
||
plt.subplots_adjust(bottom=0.15)
|
||
plt.savefig(f'{save_dir}/tornado_analysis.png', dpi=150, bbox_inches='tight', facecolor='#FAFAFA')
|
||
print(f" 保存: {save_dir}/tornado_analysis.png")
|
||
|
||
return results_annual, results_initial
|
||
|
||
|
||
def interaction_heatmap(baseline: BaselineParameters, save_dir: str):
|
||
"""
|
||
双参数交互热力图 - α vs η（马卡龙色系简约版）
|
||
"""
|
||
print("生成双参数交互热力图...")
|
||
|
||
from matplotlib.colors import LinearSegmentedColormap
|
||
|
||
alphas = np.array([1, 25, 50, 100, 150, 200, 250])
|
||
etas = np.array([0.70, 0.75, 0.80, 0.85, 0.90, 0.95])
|
||
|
||
# 年补充量矩阵
|
||
annual_matrix = np.zeros((len(etas), len(alphas)))
|
||
# 运力占比矩阵
|
||
capacity_matrix = np.zeros((len(etas), len(alphas)))
|
||
|
||
for i, eta in enumerate(etas):
|
||
for j, alpha in enumerate(alphas):
|
||
demand = calculate_water_demand(
|
||
population=baseline.population,
|
||
survival_water=baseline.survival_water,
|
||
hygiene_water_base=baseline.hygiene_water_base,
|
||
comfort_factor=alpha,
|
||
medical_water_per_event=baseline.medical_water_per_event,
|
||
sickness_rate=baseline.sickness_rate,
|
||
confidence_level=baseline.confidence_level,
|
||
recycle_rate=eta,
|
||
safety_buffer_days=baseline.safety_buffer_days
|
||
)
|
||
annual_matrix[i, j] = demand['annual_supply_tons'] / 1000 # 千吨
|
||
capacity_matrix[i, j] = demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100
|
||
|
||
# 马卡龙色系（中等饱和度）
|
||
# 图1: 淡蓝→蓝紫→玫粉（年补充量）
|
||
colors_annual = ['#C8E6F0', '#A8D0E8', '#98B8D8', '#B898C8', '#D098B0', '#D88898', '#C86878']
|
||
cmap_annual = LinearSegmentedColormap.from_list('macaron_annual', colors_annual)
|
||
|
||
# 图2: 薄荷绿→淡黄→杏橙→珊瑚红（运力占比）
|
||
colors_capacity = ['#A8D8B8', '#C8E0A0', '#E0D890', '#E8C078', '#E8A070', '#D88070', '#C86868']
|
||
cmap_capacity = LinearSegmentedColormap.from_list('macaron_capacity', colors_capacity)
|
||
|
||
# 缩小图片尺寸
|
||
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
|
||
fig.patch.set_facecolor('#FAFAFA')
|
||
|
||
# 图1: 年补充量热力图
|
||
ax = axes[0]
|
||
ax.set_facecolor('#FAFAFA')
|
||
im = ax.imshow(annual_matrix, cmap=cmap_annual, aspect='auto')
|
||
ax.set_xticks(range(len(alphas)))
|
||
ax.set_xticklabels([f'{a}' for a in alphas], fontsize=9)
|
||
ax.set_yticks(range(len(etas)))
|
||
ax.set_yticklabels([f'{e*100:.0f}%' for e in etas], fontsize=9)
|
||
ax.set_xlabel('Comfort Factor (α)', fontsize=10, color='#444444')
|
||
ax.set_ylabel('Recycle Rate (η)', fontsize=10, color='#444444')
|
||
|
||
# 添加数值标注
|
||
for i in range(len(etas)):
|
||
for j in range(len(alphas)):
|
||
ax.text(j, i, f'{annual_matrix[i, j]:.0f}', ha='center', va='center',
|
||
fontsize=8, color='#333333', fontweight='medium')
|
||
|
||
cbar1 = plt.colorbar(im, ax=ax, shrink=0.8)
|
||
cbar1.ax.tick_params(labelsize=8)
|
||
cbar1.set_label('Annual Supply (kt)', fontsize=9, color='#444444')
|
||
|
||
# 图2: 运力占比热力图
|
||
ax = axes[1]
|
||
ax.set_facecolor('#FAFAFA')
|
||
im = ax.imshow(capacity_matrix, cmap=cmap_capacity, aspect='auto', vmin=0, vmax=150)
|
||
ax.set_xticks(range(len(alphas)))
|
||
ax.set_xticklabels([f'{a}' for a in alphas], fontsize=9)
|
||
ax.set_yticks(range(len(etas)))
|
||
ax.set_yticklabels([f'{e*100:.0f}%' for e in etas], fontsize=9)
|
||
ax.set_xlabel('Comfort Factor (α)', fontsize=10, color='#444444')
|
||
ax.set_ylabel('Recycle Rate (η)', fontsize=10, color='#444444')
|
||
|
||
# 添加数值标注
|
||
for i in range(len(etas)):
|
||
for j in range(len(alphas)):
|
||
val = capacity_matrix[i, j]
|
||
ax.text(j, i, f'{val:.0f}%', ha='center', va='center',
|
||
fontsize=8, color='#333333', fontweight='medium')
|
||
|
||
cbar2 = plt.colorbar(im, ax=ax, shrink=0.8)
|
||
cbar2.ax.tick_params(labelsize=8)
|
||
cbar2.set_label('Capacity Usage (%)', fontsize=9, color='#444444')
|
||
|
||
plt.tight_layout()
|
||
plt.savefig(f'{save_dir}/interaction_heatmap.png', dpi=150, bbox_inches='tight', facecolor='#FAFAFA')
|
||
print(f" 保存: {save_dir}/interaction_heatmap.png")
|
||
|
||
return annual_matrix, capacity_matrix
|
||
|
||
|
||
def three_param_surface(baseline: BaselineParameters, save_dir: str):
|
||
"""
|
||
三参数交互3D曲面图
|
||
"""
|
||
print("生成三参数3D曲面图...")
|
||
|
||
alphas = np.linspace(1, 200, 20)
|
||
etas = np.linspace(0.75, 0.95, 20)
|
||
|
||
A, E = np.meshgrid(alphas, etas)
|
||
Z = np.zeros_like(A)
|
||
|
||
for i in range(len(etas)):
|
||
for j in range(len(alphas)):
|
||
demand = calculate_water_demand(
|
||
population=baseline.population,
|
||
survival_water=baseline.survival_water,
|
||
hygiene_water_base=baseline.hygiene_water_base,
|
||
comfort_factor=alphas[j],
|
||
medical_water_per_event=baseline.medical_water_per_event,
|
||
sickness_rate=baseline.sickness_rate,
|
||
confidence_level=baseline.confidence_level,
|
||
recycle_rate=etas[i],
|
||
safety_buffer_days=baseline.safety_buffer_days
|
||
)
|
||
Z[i, j] = demand['annual_supply_tons'] / 1000 # 千吨
|
||
|
||
fig = plt.figure(figsize=(14, 10))
|
||
ax = fig.add_subplot(111, projection='3d')
|
||
|
||
surf = ax.plot_surface(A, E * 100, Z, cmap='viridis', alpha=0.8, edgecolor='none')
|
||
|
||
# 添加基准点
|
||
baseline_demand = calculate_water_demand(
|
||
baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate,
|
||
baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days
|
||
)
|
||
ax.scatter([baseline.comfort_factor], [baseline.recycle_rate * 100],
|
||
[baseline_demand['annual_supply_tons'] / 1000],
|
||
color='red', s=100, label='Baseline', marker='*')
|
||
|
||
ax.set_xlabel('Comfort Factor (α)', fontsize=11)
|
||
ax.set_ylabel('Recycle Rate (%)', fontsize=11)
|
||
ax.set_zlabel('Annual Supply (kt)', fontsize=11)
|
||
ax.set_title('3D Surface: Annual Water Supply\nas Function of α and η', fontsize=14)
|
||
|
||
fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, label='Annual Supply (kt)')
|
||
ax.legend()
|
||
|
||
plt.tight_layout()
|
||
plt.savefig(f'{save_dir}/three_param_surface.png', dpi=150, bbox_inches='tight')
|
||
print(f" 保存: {save_dir}/three_param_surface.png")
|
||
|
||
return Z
|
||
|
||
|
||
# ============== 最坏情况分析 ==============
|
||
|
||
def worst_case_analysis(baseline: BaselineParameters, save_dir: str):
|
||
"""
|
||
最坏情况压力测试分析
|
||
"""
|
||
print("进行最坏情况分析...")
|
||
|
||
# 定义场景
|
||
scenarios = {
|
||
'Baseline': {
|
||
'alpha': baseline.comfort_factor,
|
||
'eta': baseline.recycle_rate,
|
||
'sickness': baseline.sickness_rate,
|
||
'population': baseline.population,
|
||
'buffer': baseline.safety_buffer_days,
|
||
'medical_water': baseline.medical_water_per_event,
|
||
},
|
||
'Survival Mode': {
|
||
'alpha': 1,
|
||
'eta': 0.90,
|
||
'sickness': 0.02,
|
||
'population': 100000,
|
||
'buffer': 30,
|
||
'medical_water': 5,
|
||
},
|
||
'Comfort Mode': {
|
||
'alpha': 50,
|
||
'eta': 0.90,
|
||
'sickness': 0.02,
|
||
'population': 100000,
|
||
'buffer': 30,
|
||
'medical_water': 5,
|
||
},
|
||
'Luxury Mode': {
|
||
'alpha': 250,
|
||
'eta': 0.90,
|
||
'sickness': 0.02,
|
||
'population': 100000,
|
||
'buffer': 30,
|
||
'medical_water': 5,
|
||
},
|
||
'Recycle Degradation': {
|
||
'alpha': 50,
|
||
'eta': 0.80,
|
||
'sickness': 0.02,
|
||
'population': 100000,
|
||
'buffer': 30,
|
||
'medical_water': 5,
|
||
},
|
||
'Epidemic Outbreak': {
|
||
'alpha': 50,
|
||
'eta': 0.90,
|
||
'sickness': 0.05,
|
||
'population': 100000,
|
||
'buffer': 45,
|
||
'medical_water': 10,
|
||
},
|
||
'Population Surge': {
|
||
'alpha': 50,
|
||
'eta': 0.90,
|
||
'sickness': 0.02,
|
||
'population': 120000,
|
||
'buffer': 30,
|
||
'medical_water': 5,
|
||
},
|
||
'Worst Case': {
|
||
'alpha': 250,
|
||
'eta': 0.80,
|
||
'sickness': 0.05,
|
||
'population': 120000,
|
||
'buffer': 45,
|
||
'medical_water': 10,
|
||
},
|
||
'Moderate Stress': {
|
||
'alpha': 100,
|
||
'eta': 0.85,
|
||
'sickness': 0.03,
|
||
'population': 110000,
|
||
'buffer': 40,
|
||
'medical_water': 7,
|
||
},
|
||
}
|
||
|
||
results = []
|
||
|
||
for name, params in scenarios.items():
|
||
demand = calculate_water_demand(
|
||
population=params['population'],
|
||
survival_water=baseline.survival_water,
|
||
hygiene_water_base=baseline.hygiene_water_base,
|
||
comfort_factor=params['alpha'],
|
||
medical_water_per_event=params['medical_water'],
|
||
sickness_rate=params['sickness'],
|
||
confidence_level=baseline.confidence_level,
|
||
recycle_rate=params['eta'],
|
||
safety_buffer_days=params['buffer']
|
||
)
|
||
|
||
transport = calculate_transport_metrics(
|
||
demand,
|
||
baseline.daily_elevator_capacity,
|
||
baseline.elevator_specific_energy
|
||
)
|
||
|
||
capacity_ratio = demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100
|
||
|
||
results.append({
|
||
'scenario': name,
|
||
'alpha': params['alpha'],
|
||
'eta': params['eta'] * 100,
|
||
'sickness': params['sickness'] * 100,
|
||
'population': params['population'] / 1000,
|
||
'buffer': params['buffer'],
|
||
'daily_supply': demand['daily_supply_tons'],
|
||
'annual_supply': demand['annual_supply_tons'],
|
||
'initial_transport': demand['initial_transport_tons'],
|
||
'annual_energy_pj': transport['annual_energy_pj'],
|
||
'capacity_ratio': capacity_ratio,
|
||
'feasible': capacity_ratio < 100,
|
||
})
|
||
|
||
df = pd.DataFrame(results)
|
||
|
||
# 绘图
|
||
fig, axes = plt.subplots(2, 2, figsize=(16, 14))
|
||
|
||
# 图1: 年补充量对比
|
||
ax = axes[0, 0]
|
||
colors = ['green' if r['feasible'] else 'red' for r in results]
|
||
bars = ax.barh(range(len(results)), [r['annual_supply'] / 1000 for r in results], color=colors, alpha=0.7)
|
||
ax.axvline(x=baseline.total_elevator_capacity / 1000, color='red', linestyle='--', linewidth=2, label='Elevator Capacity')
|
||
ax.set_yticks(range(len(results)))
|
||
ax.set_yticklabels([r['scenario'] for r in results])
|
||
ax.set_xlabel('Annual Supply (kt)', fontsize=12)
|
||
ax.set_title('Annual Water Supply by Scenario', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3, axis='x')
|
||
|
||
# 添加数值标注
|
||
for i, bar in enumerate(bars):
|
||
width = bar.get_width()
|
||
ax.text(width + 5, bar.get_y() + bar.get_height()/2, f'{width:.0f}',
|
||
va='center', fontsize=9)
|
||
|
||
# 图2: 运力占比
|
||
ax = axes[0, 1]
|
||
colors = ['green' if r['feasible'] else 'red' for r in results]
|
||
bars = ax.barh(range(len(results)), [r['capacity_ratio'] for r in results], color=colors, alpha=0.7)
|
||
ax.axvline(x=100, color='red', linestyle='--', linewidth=2, label='100% Capacity')
|
||
ax.set_yticks(range(len(results)))
|
||
ax.set_yticklabels([r['scenario'] for r in results])
|
||
ax.set_xlabel('Elevator Capacity Usage (%)', fontsize=12)
|
||
ax.set_title('Capacity Utilization by Scenario', fontsize=13)
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3, axis='x')
|
||
|
||
for i, bar in enumerate(bars):
|
||
width = bar.get_width()
|
||
ax.text(width + 2, bar.get_y() + bar.get_height()/2, f'{width:.1f}%',
|
||
va='center', fontsize=9)
|
||
|
||
# 图3: 年能耗对比
|
||
ax = axes[1, 0]
|
||
ax.bar(range(len(results)), [r['annual_energy_pj'] for r in results], color='purple', alpha=0.7)
|
||
ax.set_xticks(range(len(results)))
|
||
ax.set_xticklabels([r['scenario'] for r in results], rotation=45, ha='right', fontsize=9)
|
||
ax.set_ylabel('Annual Energy (PJ)', fontsize=12)
|
||
ax.set_title('Annual Transport Energy by Scenario', fontsize=13)
|
||
ax.grid(True, alpha=0.3, axis='y')
|
||
|
||
# 图4: 首批运输量
|
||
ax = axes[1, 1]
|
||
ax.bar(range(len(results)), [r['initial_transport'] / 1000 for r in results], color='teal', alpha=0.7)
|
||
ax.set_xticks(range(len(results)))
|
||
ax.set_xticklabels([r['scenario'] for r in results], rotation=45, ha='right', fontsize=9)
|
||
ax.set_ylabel('Initial Transport (kt)', fontsize=12)
|
||
ax.set_title('Initial Transport Volume by Scenario', fontsize=13)
|
||
ax.grid(True, alpha=0.3, axis='y')
|
||
|
||
plt.suptitle('Worst Case Analysis: Stress Test Scenarios', fontsize=14, y=1.02)
|
||
plt.tight_layout()
|
||
plt.savefig(f'{save_dir}/worst_case_analysis.png', dpi=150, bbox_inches='tight')
|
||
print(f" 保存: {save_dir}/worst_case_analysis.png")
|
||
|
||
# 保存详细结果表
|
||
df.to_csv(f'{save_dir}/worst_case_results.csv', index=False)
|
||
print(f" 保存: {save_dir}/worst_case_results.csv")
|
||
|
||
return df
|
||
|
||
|
||
def generate_summary_report(baseline: BaselineParameters, save_dir: str):
|
||
"""
|
||
生成灵敏度分析总结报告
|
||
"""
|
||
print("生成总结报告...")
|
||
|
||
report = []
|
||
report.append("=" * 100)
|
||
report.append("TASK 3: WATER SUPPLY SENSITIVITY ANALYSIS SUMMARY")
|
||
report.append("=" * 100)
|
||
|
||
report.append("\n## 基准参数")
|
||
report.append(f"- 人口: {baseline.population:,}")
|
||
report.append(f"- 舒适度因子 (α): {baseline.comfort_factor}")
|
||
report.append(f"- 水回收效率 (η): {baseline.recycle_rate * 100:.0f}%")
|
||
report.append(f"- 患病率: {baseline.sickness_rate * 100:.1f}%")
|
||
report.append(f"- 医疗用水: {baseline.medical_water_per_event} L/次")
|
||
report.append(f"- 置信水平: {baseline.confidence_level * 100:.1f}%")
|
||
report.append(f"- 安全缓冲: {baseline.safety_buffer_days} 天")
|
||
|
||
# 计算基准值
|
||
baseline_demand = calculate_water_demand(
|
||
baseline.population, baseline.survival_water, baseline.hygiene_water_base,
|
||
baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate,
|
||
baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days
|
||
)
|
||
baseline_transport = calculate_transport_metrics(
|
||
baseline_demand, baseline.daily_elevator_capacity, baseline.elevator_specific_energy
|
||
)
|
||
|
||
report.append(f"\n## 基准水需求")
|
||
report.append(f"- 每日补充量: {baseline_demand['daily_supply_tons']:.1f} 吨")
|
||
report.append(f"- 年补充量: {baseline_demand['annual_supply_tons']:.1f} 吨")
|
||
report.append(f"- 首批运输量: {baseline_demand['initial_transport_tons']:.1f} 吨")
|
||
report.append(f"- 年能耗: {baseline_transport['annual_energy_pj']:.4f} PJ")
|
||
report.append(f"- 运力占比: {baseline_demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100:.2f}%")
|
||
|
||
report.append("\n" + "=" * 100)
|
||
report.append("## 关键灵敏度发现")
|
||
report.append("=" * 100)
|
||
|
||
report.append("""
|
||
### 1. 水回收效率 (η) - 影响程度: ★★★★★
|
||
- 范围: 70% - 95%
|
||
- η从90%降至80%时，日补充量翻倍
|
||
- 这是对年运输需求影响最大的参数
|
||
- 建议: 优先投资水回收系统维护和升级
|
||
|
||
### 2. 舒适度因子 (α) - 影响程度: ★★★★★
|
||
- 范围: 1 - 300
|
||
- α=250时年补充量是α=1时的35倍
|
||
- 在高舒适度下可能超出电梯运力
|
||
- 建议: 初期采用生存标准，逐步提升舒适度
|
||
|
||
### 3. 人口规模 (N) - 影响程度: ★★★★☆
|
||
- 范围: ±20% (80k - 120k)
|
||
- 线性影响，20%人口增加导致20%需求增加
|
||
- 建议: 规划时预留弹性容量
|
||
|
||
### 4. 医疗紧急参数 - 影响程度: ★★☆☆☆
|
||
- 患病率(1%-5%)和医疗用水(3-15L)对年补充量影响较小
|
||
- 主要影响安全缓冲和首批运输量
|
||
- 建议: 采用保守参数确保医疗安全
|
||
|
||
### 5. 安全缓冲天数 - 影响程度: ★★★☆☆
|
||
- 范围: 15 - 60天
|
||
- 主要影响首批运输量，不影响年运输总量
|
||
- 建议: 保持30天缓冲作为运输中断容忍度
|
||
""")
|
||
|
||
report.append("\n" + "=" * 100)
|
||
report.append("## 最坏情况分析结论")
|
||
report.append("=" * 100)
|
||
|
||
report.append("""
|
||
### 可行性边界
|
||
- 最坏情况 (α=250, η=80%, N=120k): 年需求约 562 kt，超出电梯运力
|
||
- 中度压力 (α=100, η=85%, N=110k): 年需求约 132 kt，仍可承受
|
||
- 系统崩溃临界点: 约 η < 82% + α > 200 组合
|
||
|
||
### 风险缓解建议
|
||
1. 水回收系统冗余设计，确保η≥85%
|
||
2. 分阶段提升舒适度标准
|
||
3. 维持混合运输能力作为备用
|
||
4. 建立月球本地水源开发计划
|
||
""")
|
||
|
||
report_text = "\n".join(report)
|
||
|
||
with open(f'{save_dir}/sensitivity_summary_report.txt', 'w', encoding='utf-8') as f:
|
||
f.write(report_text)
|
||
print(f" 保存: {save_dir}/sensitivity_summary_report.txt")
|
||
|
||
return report_text
|
||
|
||
|
||
# ============== 主程序 ==============
|
||
|
||
def run_full_sensitivity_analysis():
|
||
"""运行完整的灵敏度分析"""
|
||
|
||
print("=" * 80)
|
||
print("任务三：月球殖民地水需求灵敏度分析")
|
||
print("=" * 80)
|
||
|
||
save_dir = '/Volumes/Files/code/mm/20260130_b/p3'
|
||
|
||
# 初始化基准参数
|
||
baseline = BaselineParameters()
|
||
|
||
print(f"\n基准参数:")
|
||
print(f" - 人口: {baseline.population:,}")
|
||
print(f" - 舒适度因子: {baseline.comfort_factor}")
|
||
print(f" - 水回收效率: {baseline.recycle_rate * 100:.0f}%")
|
||
print(f" - 患病率: {baseline.sickness_rate * 100:.1f}%")
|
||
print(f" - 安全缓冲: {baseline.safety_buffer_days} 天")
|
||
|
||
print("\n" + "=" * 80)
|
||
print("开始灵敏度分析...")
|
||
print("=" * 80)
|
||
|
||
# 1. 单参数灵敏度分析
|
||
print("\n[1/7] 水回收效率灵敏度分析")
|
||
df_recycle = sensitivity_recycle_rate(baseline, save_dir)
|
||
|
||
print("\n[2/7] 舒适度因子灵敏度分析")
|
||
df_comfort = sensitivity_comfort_factor(baseline, save_dir)
|
||
|
||
print("\n[3/7] 医疗参数灵敏度分析")
|
||
results_medical = sensitivity_medical_params(baseline, save_dir)
|
||
|
||
print("\n[4/7] 人口规模灵敏度分析")
|
||
df_population = sensitivity_population(baseline, save_dir)
|
||
|
||
print("\n[5/7] 安全缓冲天数灵敏度分析")
|
||
df_buffer = sensitivity_buffer_days(baseline, save_dir)
|
||
|
||
# 2. 多参数交互分析
|
||
print("\n[6/7] 多参数交互分析")
|
||
tornado_results = tornado_analysis(baseline, save_dir)
|
||
heatmap_data = interaction_heatmap(baseline, save_dir)
|
||
surface_data = three_param_surface(baseline, save_dir)
|
||
|
||
# 3. 最坏情况分析
|
||
print("\n[7/7] 最坏情况分析")
|
||
worst_case_df = worst_case_analysis(baseline, save_dir)
|
||
|
||
# 4. 生成总结报告
|
||
print("\n生成总结报告...")
|
||
summary = generate_summary_report(baseline, save_dir)
|
||
|
||
print("\n" + "=" * 80)
|
||
print("灵敏度分析完成!")
|
||
print("=" * 80)
|
||
|
||
print(f"\n生成的文件:")
|
||
print(f" - sensitivity_recycle_rate.png")
|
||
print(f" - sensitivity_comfort_factor.png")
|
||
print(f" - sensitivity_medical_params.png")
|
||
print(f" - sensitivity_population.png")
|
||
print(f" - sensitivity_buffer_days.png")
|
||
print(f" - tornado_analysis.png")
|
||
print(f" - interaction_heatmap.png")
|
||
print(f" - three_param_surface.png")
|
||
print(f" - worst_case_analysis.png")
|
||
print(f" - worst_case_results.csv")
|
||
print(f" - sensitivity_summary_report.txt")
|
||
|
||
return {
|
||
'recycle': df_recycle,
|
||
'comfort': df_comfort,
|
||
'population': df_population,
|
||
'buffer': df_buffer,
|
||
'worst_case': worst_case_df,
|
||
}
|
||
|
||
|
||
if __name__ == "__main__":
|
||
results = run_full_sensitivity_analysis()
|