""" Task 4: Environmental Impact Analysis for Moon Colony Construction 使用液氧甲烷(LOX/CH4)燃料的环境影响评估 分析不同建设方案对地球环境的影响，并给出最小化环境影响的模型调整策略。 """ import numpy as np import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt from matplotlib import rcParams from dataclasses import dataclass from typing import List, Dict, Tuple import pandas as pd # 设置字体 rcParams['font.sans-serif'] = ['Arial Unicode MS', 'DejaVu Sans', 'SimHei'] rcParams['axes.unicode_minus'] = False # ============== 物理常数 ============== G0 = 9.81 # m/s² OMEGA_EARTH = 7.27e-5 # rad/s R_EARTH = 6.371e6 # m # ============== 任务参数 (与任务一一致) ============== TOTAL_PAYLOAD = 100e6 # 100 million metric tons # ============== 太空电梯参数 ============== NUM_ELEVATORS = 3 ELEVATOR_CAPACITY_PER_YEAR = 179000 # metric tons per elevator per year TOTAL_ELEVATOR_CAPACITY = NUM_ELEVATORS * ELEVATOR_CAPACITY_PER_YEAR # 537,000 tons/year ELEVATOR_SPECIFIC_ENERGY = 157.2e9 # J per metric ton (157.2 MJ/kg) # 电梯能量分解 ELEVATOR_LIFT_ENERGY = 57.0e9 # J/ton (电力提升能) ELEVATOR_TOP_ROCKET_ENERGY = 100.2e9 # J/ton (顶端火箭推进) # ============== 液氧甲烷燃料参数 (LOX/CH4) ============== @dataclass class MethaneFuelParams: """液氧甲烷燃料参数""" name: str = "LOX/CH4" isp: float = 360 # 比冲 (秒), Raptor发动机级别 specific_energy: float = 12.9e6 # 燃料比能量 (J/kg) ox_fuel_ratio: float = 3.5 # 氧燃比 O2:CH4 @property def exhaust_velocity(self) -> float: """排气速度 (m/s)""" return self.isp * G0 @property def ch4_mass_fraction(self) -> float: """甲烷在燃料中的质量分数""" return 1 / (1 + self.ox_fuel_ratio) # ≈ 0.222 @property def lox_mass_fraction(self) -> float: """液氧在燃料中的质量分数""" return self.ox_fuel_ratio / (1 + self.ox_fuel_ratio) # ≈ 0.778 # ============== 排放因子 (基于化学计量) ============== @dataclass class EmissionFactors: """排放因子 (kg排放物/kg燃料)""" # 燃烧排放: CH4 + 2O2 → CO2 + 2H2O # 每kg CH4产生: 44/16 = 2.75 kg CO2, 36/16 = 2.25 kg H2O fuel_params: MethaneFuelParams = None def __post_init__(self): if self.fuel_params is None: self.fuel_params = MethaneFuelParams() @property def co2_per_kg_fuel(self) -> float: """每kg燃料的CO2排放 (kg)""" return self.fuel_params.ch4_mass_fraction * 2.75 # ≈ 0.611 kg CO2/kg燃料 @property def h2o_per_kg_fuel(self) -> float: """每kg燃料的H2O排放 (kg)""" return self.fuel_params.ch4_mass_fraction * 2.25 # ≈ 0.500 kg H2O/kg燃料 # 燃料生产排放 ch4_production_co2: float = 2.5 # kg CO2/kg CH4 (天然气制甲烷，含开采和加工) lox_production_kwh: float = 0.5 # kWh/kg LOX (空气分离) # 电网碳强度 grid_carbon_intensity: float = 0.475 # kg CO2/kWh (全球平均, IEA 2023) grid_carbon_intensity_renewable: float = 0.05 # kg CO2/kWh (可再生能源) # 平流层排放比例 stratosphere_emission_fraction: float = 0.40 # 40%的排放进入平流层 stratosphere_h2o_residence_time: float = 3.0 # 年 (平流层水蒸气滞留时间) # ============== 发射场定义 ============== @dataclass class LaunchSite: name: str short_name: str latitude: float max_launches_per_day: int = 1 @property def abs_latitude(self) -> float: return abs(self.latitude) @property def rotation_velocity(self) -> float: return OMEGA_EARTH * R_EARTH * np.cos(np.radians(self.abs_latitude)) @property def delta_v_loss(self) -> float: v_equator = OMEGA_EARTH * R_EARTH return v_equator - self.rotation_velocity # 10个发射场 (按纬度排序) LAUNCH_SITES = sorted([ LaunchSite("Kourou (French Guiana)", "Kourou", 5.2), LaunchSite("Satish Dhawan (India)", "SDSC", 13.7), LaunchSite("Boca Chica (Texas)", "Texas", 26.0), LaunchSite("Cape Canaveral (Florida)", "Florida", 28.5), LaunchSite("Vandenberg (California)", "California", 34.7), LaunchSite("Wallops (Virginia)", "Virginia", 37.8), LaunchSite("Taiyuan (China)", "Taiyuan", 38.8), LaunchSite("Mahia (New Zealand)", "Mahia", 39.3), LaunchSite("Baikonur (Kazakhstan)", "Baikonur", 45.6), LaunchSite("Kodiak (Alaska)", "Alaska", 57.4), ], key=lambda x: x.abs_latitude) # ============== 火箭参数 ============== PAYLOAD_PER_LAUNCH = 125 # metric tons per launch ALPHA = 0.10 # 结构系数 NUM_STAGES = 3 DELTA_V_BASE = 13300 # m/s (赤道发射到月球) FUEL_PARAMS = MethaneFuelParams() EMISSION_FACTORS = EmissionFactors(fuel_params=FUEL_PARAMS) # ============== 核心计算函数 ============== def fuel_ratio_multistage(delta_v: float, fuel_params: MethaneFuelParams = FUEL_PARAMS) -> float: """多级火箭燃料/载荷比""" ve = fuel_params.exhaust_velocity delta_v_per_stage = delta_v / NUM_STAGES R_stage = np.exp(delta_v_per_stage / ve) denominator = 1 - ALPHA * (R_stage - 1) if denominator <= 0: return np.inf k_stage = (R_stage - 1) / denominator total_fuel_ratio = 0 remaining_ratio = 1.0 for _ in range(NUM_STAGES): fuel_this_stage = remaining_ratio * k_stage total_fuel_ratio += fuel_this_stage remaining_ratio *= (1 + k_stage * (1 + ALPHA)) return total_fuel_ratio def rocket_fuel_per_ton(site: LaunchSite) -> float: """每吨载荷的燃料消耗 (kg燃料/ton载荷)""" delta_v = DELTA_V_BASE + site.delta_v_loss k = fuel_ratio_multistage(delta_v) return k * 1000 # kg fuel per metric ton payload def rocket_energy_per_ton(site: LaunchSite) -> float: """火箭发射每吨载荷的能量消耗 (J/ton)""" return rocket_fuel_per_ton(site) * FUEL_PARAMS.specific_energy # ============== 环境影响计算 ============== @dataclass class EnvironmentalImpact: """环境影响结果""" scenario_name: str completion_years: float total_payload_tons: float # 能量消耗 total_energy_PJ: float rocket_energy_PJ: float elevator_energy_PJ: float # 燃料消耗 total_fuel_Mt: float # 百万吨 rocket_fuel_Mt: float elevator_top_rocket_fuel_Mt: float # 排放量 co2_combustion_Mt: float # 燃烧直接排放 co2_fuel_production_Mt: float # 燃料生产排放 co2_electricity_Mt: float # 电力消耗排放 h2o_total_Mt: float h2o_stratosphere_Mt: float # 平流层水蒸气 # 发射统计 rocket_launches: int @property def co2_total_Mt(self) -> float: """总CO2排放 (百万吨)""" return self.co2_combustion_Mt + self.co2_fuel_production_Mt + self.co2_electricity_Mt @property def co2_per_ton_payload(self) -> float: """每吨载荷的CO2排放 (kg)""" return self.co2_total_Mt * 1e9 / self.total_payload_tons @property def annual_co2_Mt(self) -> float: """年均CO2排放 (百万吨/年)""" return self.co2_total_Mt / self.completion_years if self.completion_years > 0 else 0 def calculate_scenario_impact( scenario_name: str, completion_years: float, elevator_fraction: float, rocket_sites_used: List[LaunchSite], emission_factors: EmissionFactors = EMISSION_FACTORS, renewable_energy_fraction: float = 0.0 ) -> EnvironmentalImpact: """ 计算单个情景的环境影响 参数: scenario_name: 情景名称 completion_years: 完成年限 elevator_fraction: 电梯运输比例 (0-1) rocket_sites_used: 使用的火箭发射场列表 emission_factors: 排放因子 renewable_energy_fraction: 可再生能源占比 (0-1) """ # 载荷分配 elevator_payload = TOTAL_PAYLOAD * elevator_fraction rocket_payload = TOTAL_PAYLOAD * (1 - elevator_fraction) # ===== 电梯部分 ===== elevator_energy = elevator_payload * ELEVATOR_SPECIFIC_ENERGY # J elevator_lift_energy = elevator_payload * ELEVATOR_LIFT_ENERGY # J (电力) elevator_top_rocket_energy = elevator_payload * ELEVATOR_TOP_ROCKET_ENERGY # J (化学能) # 电梯顶端火箭燃料 (使用较小的燃料比，因为已在高轨道) # 顶端ΔV约1800 m/s (LOI + 轨道调整) elevator_top_delta_v = 1800 elevator_top_fuel_ratio = fuel_ratio_multistage(elevator_top_delta_v) elevator_top_rocket_fuel = elevator_payload * 1000 * elevator_top_fuel_ratio # kg # ===== 火箭部分 ===== if rocket_payload > 0 and len(rocket_sites_used) > 0: launches_per_site = int(np.ceil(rocket_payload / PAYLOAD_PER_LAUNCH / len(rocket_sites_used))) total_launches = launches_per_site * len(rocket_sites_used) # 按发射场计算燃料消耗 (简化：均匀分配) rocket_fuel_total = 0 rocket_energy_total = 0 for site in rocket_sites_used: site_payload = rocket_payload / len(rocket_sites_used) site_fuel = site_payload * 1000 * fuel_ratio_multistage(DELTA_V_BASE + site.delta_v_loss) rocket_fuel_total += site_fuel rocket_energy_total += site_fuel * FUEL_PARAMS.specific_energy else: total_launches = 0 rocket_fuel_total = 0 rocket_energy_total = 0 # ===== 总燃料消耗 ===== total_fuel = rocket_fuel_total + elevator_top_rocket_fuel # kg # ===== 排放计算 ===== # 1. 燃烧直接排放 co2_combustion = total_fuel * emission_factors.co2_per_kg_fuel # kg h2o_combustion = total_fuel * emission_factors.h2o_per_kg_fuel # kg # 2. 燃料生产排放 ch4_mass = total_fuel * FUEL_PARAMS.ch4_mass_fraction lox_mass = total_fuel * FUEL_PARAMS.lox_mass_fraction co2_ch4_production = ch4_mass * emission_factors.ch4_production_co2 co2_lox_production = lox_mass * emission_factors.lox_production_kwh * emission_factors.grid_carbon_intensity co2_fuel_production = co2_ch4_production + co2_lox_production # 3. 电梯电力消耗排放 elevator_electricity_kwh = elevator_lift_energy / 3.6e6 # J -> kWh grid_intensity = (emission_factors.grid_carbon_intensity * (1 - renewable_energy_fraction) + emission_factors.grid_carbon_intensity_renewable * renewable_energy_fraction) co2_electricity = elevator_electricity_kwh * grid_intensity # 4. 平流层水蒸气 (仅火箭发射部分) # 电梯顶端释放点已在GEO以上，不经过平流层 h2o_rocket = rocket_fuel_total * emission_factors.h2o_per_kg_fuel h2o_stratosphere = h2o_rocket * emission_factors.stratosphere_emission_fraction return EnvironmentalImpact( scenario_name=scenario_name, completion_years=completion_years, total_payload_tons=TOTAL_PAYLOAD, total_energy_PJ=(rocket_energy_total + elevator_energy) / 1e15, rocket_energy_PJ=rocket_energy_total / 1e15, elevator_energy_PJ=elevator_energy / 1e15, total_fuel_Mt=total_fuel / 1e9, rocket_fuel_Mt=rocket_fuel_total / 1e9, elevator_top_rocket_fuel_Mt=elevator_top_rocket_fuel / 1e9, co2_combustion_Mt=co2_combustion / 1e9, co2_fuel_production_Mt=co2_fuel_production / 1e9, co2_electricity_Mt=co2_electricity / 1e9, h2o_total_Mt=h2o_combustion / 1e9, h2o_stratosphere_Mt=h2o_stratosphere / 1e9, rocket_launches=total_launches ) # ============== 情景定义 ============== def define_scenarios() -> List[Dict]: """ 定义四个基准情景 (基于任务一结果) """ scenarios = [] # 情景1: 纯火箭 (全部10个发射场) rocket_only_years = TOTAL_PAYLOAD / (len(LAUNCH_SITES) * 365 * PAYLOAD_PER_LAUNCH) scenarios.append({ 'name': 'Rocket Only', 'completion_years': rocket_only_years, # ~219年 'elevator_fraction': 0.0, 'rocket_sites': LAUNCH_SITES.copy(), 'description': '纯火箭方案：10个发射场满负荷' }) # 情景2: 纯电梯 elevator_only_years = TOTAL_PAYLOAD / TOTAL_ELEVATOR_CAPACITY scenarios.append({ 'name': 'Elevator Only', 'completion_years': elevator_only_years, # ~186年 'elevator_fraction': 1.0, 'rocket_sites': [], 'description': '纯电梯方案：3个电梯满负荷' }) # 情景3: 混合方案 (最短时间) combined_capacity = TOTAL_ELEVATOR_CAPACITY + len(LAUNCH_SITES) * 365 * PAYLOAD_PER_LAUNCH combined_min_years = TOTAL_PAYLOAD / combined_capacity elevator_payload_min = TOTAL_ELEVATOR_CAPACITY * combined_min_years elevator_fraction_min = elevator_payload_min / TOTAL_PAYLOAD scenarios.append({ 'name': 'Combined (Min Time)', 'completion_years': combined_min_years, # ~101年 'elevator_fraction': elevator_fraction_min, # ~54% 'rocket_sites': LAUNCH_SITES.copy(), 'description': '混合方案(最短时间)：电梯+全部火箭站点' }) # 情景4: 混合方案 (膝点优化) knee_years = 139 # 从任务一Pareto分析 elevator_payload_knee = min(TOTAL_ELEVATOR_CAPACITY * knee_years, TOTAL_PAYLOAD) elevator_fraction_knee = elevator_payload_knee / TOTAL_PAYLOAD # 膝点时只需要部分低纬度发射场 scenarios.append({ 'name': 'Combined (Balanced)', 'completion_years': knee_years, 'elevator_fraction': elevator_fraction_knee, # ~74.6% 'rocket_sites': LAUNCH_SITES[:5], # 只使用5个低纬度站点 'description': '混合方案(膝点)：电梯优先+低纬度火箭' }) return scenarios def analyze_all_scenarios(renewable_fraction: float = 0.0) -> pd.DataFrame: """分析所有情景""" scenarios = define_scenarios() results = [] for scenario in scenarios: impact = calculate_scenario_impact( scenario_name=scenario['name'], completion_years=scenario['completion_years'], elevator_fraction=scenario['elevator_fraction'], rocket_sites_used=scenario['rocket_sites'], renewable_energy_fraction=renewable_fraction ) results.append({ 'Scenario': impact.scenario_name, 'Years': impact.completion_years, 'Elevator %': impact.elevator_energy_PJ / impact.total_energy_PJ * 100 if impact.total_energy_PJ > 0 else 0, 'Energy (PJ)': impact.total_energy_PJ, 'Fuel (Mt)': impact.total_fuel_Mt, 'CO2 Total (Mt)': impact.co2_total_Mt, 'CO2 Combustion (Mt)': impact.co2_combustion_Mt, 'CO2 Production (Mt)': impact.co2_fuel_production_Mt, 'CO2 Electricity (Mt)': impact.co2_electricity_Mt, 'H2O Total (Mt)': impact.h2o_total_Mt, 'H2O Stratosphere (Mt)': impact.h2o_stratosphere_Mt, 'Annual CO2 (Mt/yr)': impact.annual_co2_Mt, 'CO2/ton (kg)': impact.co2_per_ton_payload, 'Launches': impact.rocket_launches, }) return pd.DataFrame(results) # ============== 模型优化：环境约束优化 ============== def optimize_for_environment( max_annual_co2_Mt: float = 50.0, min_completion_years: float = 100.0, renewable_fraction: float = 0.5 ) -> Dict: """ 带环境约束的优化模型 目标: 最小化 总CO2排放 约束: - 完成时间 >= min_completion_years - 年均CO2 <= max_annual_co2_Mt 决策变量: - 电梯运输比例 (0-1) - 使用的火箭站点数量 (0-10) """ best_result = None best_co2 = float('inf') # 网格搜索 for elevator_frac in np.linspace(0, 1, 21): for n_sites in range(0, 11): # 计算完成时间 elevator_capacity = TOTAL_ELEVATOR_CAPACITY rocket_capacity = n_sites * 365 * PAYLOAD_PER_LAUNCH if n_sites > 0 else 0 elevator_payload = TOTAL_PAYLOAD * elevator_frac rocket_payload = TOTAL_PAYLOAD * (1 - elevator_frac) if elevator_frac == 1.0: years = elevator_payload / elevator_capacity if elevator_capacity > 0 else float('inf') elif elevator_frac == 0.0: years = rocket_payload / rocket_capacity if rocket_capacity > 0 else float('inf') else: # 电梯和火箭并行 elevator_years = elevator_payload / elevator_capacity rocket_years = rocket_payload / rocket_capacity if rocket_capacity > 0 else float('inf') years = max(elevator_years, rocket_years) if years < min_completion_years: continue # 计算环境影响 sites_used = LAUNCH_SITES[:n_sites] if n_sites > 0 else [] impact = calculate_scenario_impact( scenario_name='Optimization', completion_years=years, elevator_fraction=elevator_frac, rocket_sites_used=sites_used, renewable_energy_fraction=renewable_fraction ) # 检查约束 if impact.annual_co2_Mt > max_annual_co2_Mt: continue # 更新最优解 if impact.co2_total_Mt < best_co2: best_co2 = impact.co2_total_Mt best_result = { 'elevator_fraction': elevator_frac, 'n_sites': n_sites, 'years': years, 'impact': impact } return best_result # ============== 可视化 ============== def plot_environmental_comparison(save_dir: str = '/Volumes/Files/code/mm/20260130_b/p4'): """绘制环境影响对比图（2张子图版本，马卡龙色系论文风格）""" df = analyze_all_scenarios() # 马卡龙色系 MACARON = { 'pink': '#F8B4C0', # 樱花粉 'mint': '#A8E6CF', # 薄荷绿 'lavender': '#C5A3FF', # 薰衣草紫 'peach': '#FFCBA4', # 蜜桃橙 'sky': '#A8D8EA', # 天空蓝 'cream': '#FFEAA7', # 奶油黄 'rose': '#DDA0DD', # 玫瑰粉 'sage': '#B5EAD7', # 鼠尾草绿 } # 论文风格设置 plt.rcParams.update({ 'font.size': 11, 'axes.labelsize': 12, 'axes.titlesize': 13, 'xtick.labelsize': 10, 'ytick.labelsize': 10, 'legend.fontsize': 9, 'axes.linewidth': 0.8, 'axes.edgecolor': '#666666', }) fig, axes = plt.subplots(1, 2, figsize=(10, 4.5)) fig.patch.set_facecolor('white') scenarios = df['Scenario'].tolist() x = np.arange(len(scenarios)) width = 0.35 # ========== 左图: CO2排放分解 + H2O排放 ========== ax1 = axes[0] ax1.set_facecolor('white') # CO2数据（左侧柱子） co2_combustion = df['CO2 Combustion (Mt)'].values co2_production = df['CO2 Production (Mt)'].values co2_electricity = df['CO2 Electricity (Mt)'].values bars1 = ax1.bar(x - width/2, co2_combustion, width, label='Combustion', color=MACARON['pink'], edgecolor='#666666', linewidth=0.5) bars2 = ax1.bar(x - width/2, co2_production, width, bottom=co2_combustion, label='Fuel Prod.', color=MACARON['mint'], edgecolor='#666666', linewidth=0.5) bars3 = ax1.bar(x - width/2, co2_electricity, width, bottom=co2_combustion+co2_production, label='Electricity', color=MACARON['sky'], edgecolor='#666666', linewidth=0.5) ax1.set_ylabel('CO₂ Emissions (Mt)', color='#333333') ax1.tick_params(axis='y', labelcolor='#333333') # 添加CO2总量标签 totals = df['CO2 Total (Mt)'].values for i, total in enumerate(totals): ax1.text(i - width/2, total + 300, f'{total:.0f}', ha='center', fontsize=9, fontweight='bold') # H2O数据（右侧柱子） ax1_twin = ax1.twinx() h2o_surface = df['H2O Total (Mt)'].values - df['H2O Stratosphere (Mt)'].values h2o_strat = df['H2O Stratosphere (Mt)'].values bars4 = ax1_twin.bar(x + width/2, h2o_surface, width, label='Troposphere', color=MACARON['sage'], edgecolor='#666666', linewidth=0.5) bars5 = ax1_twin.bar(x + width/2, h2o_strat, width, bottom=h2o_surface, label='Stratosphere', color=MACARON['rose'], edgecolor='#666666', linewidth=0.5) ax1_twin.set_ylabel('H₂O Emissions (Mt)', color='#7B68A0') ax1_twin.tick_params(axis='y', labelcolor='#7B68A0') ax1.set_title('(a) CO₂ & H₂O Emissions', fontweight='bold', pad=10) ax1.set_xticks(x) ax1.set_xticklabels(scenarios, rotation=20, ha='right') ax1.grid(True, alpha=0.2, axis='y', linestyle='--') ax1.set_axisbelow(True) # 合并图例 lines1, labels1 = ax1.get_legend_handles_labels() lines2, labels2 = ax1_twin.get_legend_handles_labels() ax1.legend(lines1 + lines2, labels1 + labels2, loc='upper right', framealpha=0.95, edgecolor='#cccccc', ncol=2) # ========== 右图: 年均CO2排放 + 碳强度 ========== ax2 = axes[1] ax2.set_facecolor('white') annual_co2 = df['Annual CO2 (Mt/yr)'].values co2_per_ton = df['CO2/ton (kg)'].values # 年均CO2（左侧柱子）- 使用渐变马卡龙色 colors_annual = [MACARON['pink'], MACARON['mint'], MACARON['sky'], MACARON['sage']] bars_annual = ax2.bar(x - width/2, annual_co2, width, color=colors_annual, edgecolor='#666666', linewidth=0.5) ax2.set_ylabel('Annual CO₂ (Mt/yr)', color='#333333') ax2.tick_params(axis='y', labelcolor='#333333') # 添加年均CO2标签 global_annual_co2 = 37000 for i, val in enumerate(annual_co2): pct = val / global_annual_co2 * 100 ax2.text(i - width/2, val + 2, f'{val:.1f}\n({pct:.2f}%)', ha='center', fontsize=8) # 碳强度（右侧柱子） ax2_twin = ax2.twinx() bars_intensity = ax2_twin.bar(x + width/2, co2_per_ton, width, color=MACARON['cream'], edgecolor='#666666', linewidth=0.5) ax2_twin.set_ylabel('CO₂ per ton payload (kg)', color='#B8860B') ax2_twin.tick_params(axis='y', labelcolor='#B8860B') # 添加碳强度标签 for i, val in enumerate(co2_per_ton): ax2_twin.text(i + width/2, val + 3000, f'{val:.0f}', ha='center', fontsize=9, fontweight='bold') ax2.set_title('(b) Annual Emissions & Carbon Intensity', fontweight='bold', pad=10) ax2.set_xticks(x) ax2.set_xticklabels(scenarios, rotation=20, ha='right') ax2.grid(True, alpha=0.2, axis='y', linestyle='--') ax2.set_axisbelow(True) # 图例说明 from matplotlib.patches import Patch legend_elements = [Patch(facecolor=MACARON['pink'], edgecolor='#666666', label='Annual CO₂'), Patch(facecolor=MACARON['cream'], edgecolor='#666666', label='Carbon Intensity')] ax2.legend(handles=legend_elements, loc='upper right', framealpha=0.95, edgecolor='#cccccc') plt.tight_layout() plt.savefig(f'{save_dir}/environmental_comparison.png', dpi=300, bbox_inches='tight', facecolor='white', edgecolor='none') print(f"环境影响对比图已保存至: {save_dir}/environmental_comparison.png") return df def plot_radar_chart(save_dir: str = '/Volumes/Files/code/mm/20260130_b/p4'): """单独绘制雷达图（马卡龙色系，无标题，论文风格）""" df = analyze_all_scenarios() # 马卡龙色系 MACARON = { 'pink': '#F8B4C0', # 樱花粉 'mint': '#A8E6CF', # 薄荷绿 'lavender': '#C5A3FF', # 薰衣草紫 'sky': '#A8D8EA', # 天空蓝 } colors = [MACARON['pink'], MACARON['mint'], MACARON['sky'], MACARON['lavender']] # 论文风格设置 plt.rcParams.update({ 'font.size': 11, 'axes.labelsize': 12, 'xtick.labelsize': 11, 'ytick.labelsize': 10, 'legend.fontsize': 10, }) fig, ax = plt.subplots(figsize=(5, 5), subplot_kw=dict(projection='polar')) fig.patch.set_facecolor('white') categories = ['Time\n(inv)', 'Energy\n(inv)', 'CO₂\n(inv)', 'Launches\n(inv)', 'Strat. H₂O\n(inv)'] N = len(categories) # 归一化数据 (取逆以使"更好"在外围) max_vals = df[['Years', 'Energy (PJ)', 'CO2 Total (Mt)', 'Launches', 'H2O Stratosphere (Mt)']].max() angles = np.linspace(0, 2*np.pi, N, endpoint=False).tolist() angles += angles[:1] for i, (idx, row) in enumerate(df.iterrows()): values = [ 1 - row['Years'] / max_vals['Years'], 1 - row['Energy (PJ)'] / max_vals['Energy (PJ)'], 1 - row['CO2 Total (Mt)'] / max_vals['CO2 Total (Mt)'], 1 - row['Launches'] / max_vals['Launches'] if max_vals['Launches'] > 0 else 1, 1 - row['H2O Stratosphere (Mt)'] / max_vals['H2O Stratosphere (Mt)'] if max_vals['H2O Stratosphere (Mt)'] > 0 else 1 ] values += values[:1] ax.plot(angles, values, 'o-', linewidth=2, label=row['Scenario'], color=colors[i], markersize=5) ax.fill(angles, values, alpha=0.2, color=colors[i]) ax.set_xticks(angles[:-1]) ax.set_xticklabels(categories, fontsize=10) ax.set_ylim(0, 1) # 设置径向刻度 ax.set_yticks([0.2, 0.4, 0.6, 0.8, 1.0]) ax.set_yticklabels(['0.2', '0.4', '0.6', '0.8', '1.0'], fontsize=8, color='#666666') # 图例放在右侧 ax.legend(loc='upper left', bbox_to_anchor=(1.05, 1.0), framealpha=0.95, edgecolor='#cccccc') # 网格样式 ax.grid(True, alpha=0.3, linestyle='--') ax.spines['polar'].set_color('#666666') ax.spines['polar'].set_linewidth(0.8) plt.tight_layout() plt.savefig(f'{save_dir}/radar_chart.png', dpi=300, bbox_inches='tight', facecolor='white', edgecolor='none') print(f"雷达图已保存至: {save_dir}/radar_chart.png") return df def plot_marginal_benefit(save_dir: str = '/Volumes/Files/code/mm/20260130_b/p4'): """单独绘制边际收益分析图（马卡龙色系，无标题，论文风格）""" # 马卡龙色系 MACARON = { 'sky': '#A8D8EA', # 天空蓝 'pink': '#F8B4C0', # 樱花粉 'mint': '#A8E6CF', # 薄荷绿 } # 论文风格设置 plt.rcParams.update({ 'font.size': 11, 'axes.labelsize': 12, 'xtick.labelsize': 10, 'ytick.labelsize': 10, 'legend.fontsize': 10, 'axes.linewidth': 0.8, }) # 生成Pareto数据 years_range = np.linspace(101, 200, 100) results = [] for years in years_range: # 计算电梯运输比例 elevator_payload = min(TOTAL_ELEVATOR_CAPACITY * years, TOTAL_PAYLOAD) elevator_frac = elevator_payload / TOTAL_PAYLOAD # 确定使用的发射场 rocket_payload = TOTAL_PAYLOAD - elevator_payload if rocket_payload > 0: rocket_capacity_per_site = 365 * years * PAYLOAD_PER_LAUNCH n_sites_needed = int(np.ceil(rocket_payload / rocket_capacity_per_site)) n_sites_needed = min(n_sites_needed, 10) sites = LAUNCH_SITES[:n_sites_needed] else: sites = [] impact = calculate_scenario_impact( scenario_name=f'{years:.0f}yr', completion_years=years, elevator_fraction=elevator_frac, rocket_sites_used=sites ) results.append({ 'Years': years, 'Energy (PJ)': impact.total_energy_PJ, 'CO2 (Mt)': impact.co2_total_Mt, }) df = pd.DataFrame(results) # 计算边际减少率 T = df['Years'].values E = df['Energy (PJ)'].values C = df['CO2 (Mt)'].values dE_dT = np.gradient(E, T) dC_dT = np.gradient(C, T) energy_reduction_rate = -dE_dT # 每多1年节省的能量 co2_reduction_rate = -dC_dT # 每多1年减少的CO2 # 绘图（黄金比例） fig, ax = plt.subplots(figsize=(8, 8*0.618)) fig.patch.set_facecolor('white') ax.set_facecolor('white') # 能量减少率 ax.plot(T, energy_reduction_rate, color=MACARON['sky'], linewidth=2, label='Energy reduction rate') ax.fill_between(T, energy_reduction_rate, alpha=0.2, color=MACARON['sky']) # CO2减少率（使用第二Y轴） ax2 = ax.twinx() ax2.plot(T, co2_reduction_rate, color=MACARON['pink'], linewidth=2, label='CO₂ reduction rate') ax2.fill_between(T, co2_reduction_rate, alpha=0.2, color=MACARON['pink']) # 推荐点 recommended_year = 186 ax.axvline(x=recommended_year, color=MACARON['mint'], linestyle='--', linewidth=2, label=f'Recommended: {recommended_year}yr') ax.set_xlabel('Completion Time (years)') ax.set_ylabel('Energy Reduction (PJ/yr)', color='#5B9BD5') ax.tick_params(axis='y', labelcolor='#5B9BD5') ax2.set_ylabel('CO₂ Reduction (Mt/yr)', color='#E57373') ax2.tick_params(axis='y', labelcolor='#E57373') ax.grid(True, alpha=0.2, linestyle='--') ax.set_axisbelow(True) # 合并图例 lines1, labels1 = ax.get_legend_handles_labels() lines2, labels2 = ax2.get_legend_handles_labels() ax.legend(lines1 + lines2, labels1 + labels2, loc='upper right', framealpha=0.95, edgecolor='#cccccc') plt.tight_layout() plt.savefig(f'{save_dir}/marginal_benefit.png', dpi=300, bbox_inches='tight', facecolor='white', edgecolor='none') print(f"边际收益图已保存至: {save_dir}/marginal_benefit.png") return df def plot_pareto_with_environment(save_dir: str = '/Volumes/Files/code/mm/20260130_b/p4'): """绘制时间-能量-环境三目标Pareto图""" fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # 生成Pareto数据 years_range = np.linspace(101, 220, 50) results = [] for years in years_range: # 计算电梯运输比例 elevator_payload = min(TOTAL_ELEVATOR_CAPACITY * years, TOTAL_PAYLOAD) elevator_frac = elevator_payload / TOTAL_PAYLOAD # 确定使用的发射场 rocket_payload = TOTAL_PAYLOAD - elevator_payload if rocket_payload > 0: rocket_capacity_per_site = 365 * years * PAYLOAD_PER_LAUNCH n_sites_needed = int(np.ceil(rocket_payload / rocket_capacity_per_site)) n_sites_needed = min(n_sites_needed, 10) sites = LAUNCH_SITES[:n_sites_needed] else: sites = [] impact = calculate_scenario_impact( scenario_name=f'{years:.0f}yr', completion_years=years, elevator_fraction=elevator_frac, rocket_sites_used=sites ) results.append({ 'Years': years, 'Energy (PJ)': impact.total_energy_PJ, 'CO2 (Mt)': impact.co2_total_Mt, 'Elevator %': elevator_frac * 100 }) df = pd.DataFrame(results) # ========== 图1: 时间 vs CO2 ========== ax1 = axes[0] scatter = ax1.scatter(df['Years'], df['CO2 (Mt)'], c=df['Elevator %'], cmap='RdYlGn', s=50, alpha=0.8) plt.colorbar(scatter, ax=ax1, label='Elevator %') ax1.set_xlabel('Completion Time (years)', fontsize=12) ax1.set_ylabel('Total CO₂ Emissions (Mt)', fontsize=12) ax1.set_title('Pareto Front: Time vs CO₂ Emissions\n(Color: Elevator Usage)', fontsize=14) ax1.grid(True, alpha=0.3) # 标记关键点 ax1.axvline(x=139, color='red', linestyle='--', alpha=0.7, label='Balanced (139yr)') ax1.axvline(x=186, color='green', linestyle='--', alpha=0.7, label='Elevator Only (186yr)') ax1.legend() # ========== 图2: 能量 vs CO2 ========== ax2 = axes[1] scatter = ax2.scatter(df['Energy (PJ)'], df['CO2 (Mt)'], c=df['Years'], cmap='viridis', s=50, alpha=0.8) plt.colorbar(scatter, ax=ax2, label='Years') ax2.set_xlabel('Total Energy (PJ)', fontsize=12) ax2.set_ylabel('Total CO₂ Emissions (Mt)', fontsize=12) ax2.set_title('Pareto Front: Energy vs CO₂ Emissions\n(Color: Completion Time)', fontsize=14) ax2.grid(True, alpha=0.3) plt.tight_layout() plt.savefig(f'{save_dir}/pareto_environmental.png', dpi=150, bbox_inches='tight') print(f"Pareto环境分析图已保存至: {save_dir}/pareto_environmental.png") return df def plot_mitigation_strategies(save_dir: str = '/Volumes/Files/code/mm/20260130_b/p4'): """绘制减排策略对比图""" # 基准情景 (混合膝点) baseline = calculate_scenario_impact( scenario_name='Baseline', completion_years=139, elevator_fraction=0.746, rocket_sites_used=LAUNCH_SITES[:5], renewable_energy_fraction=0.0 ) # 策略1: 50%可再生能源 strategy1 = calculate_scenario_impact( scenario_name='50% Renewable', completion_years=139, elevator_fraction=0.746, rocket_sites_used=LAUNCH_SITES[:5], renewable_energy_fraction=0.5 ) # 策略2: 100%可再生能源 strategy2 = calculate_scenario_impact( scenario_name='100% Renewable', completion_years=139, elevator_fraction=0.746, rocket_sites_used=LAUNCH_SITES[:5], renewable_energy_fraction=1.0 ) # 策略3: 延长工期到纯电梯 strategy3 = calculate_scenario_impact( scenario_name='Elevator Only', completion_years=186, elevator_fraction=1.0, rocket_sites_used=[], renewable_energy_fraction=0.0 ) # 策略4: 纯电梯 + 100%可再生 strategy4 = calculate_scenario_impact( scenario_name='Elevator + 100% Renewable', completion_years=186, elevator_fraction=1.0, rocket_sites_used=[], renewable_energy_fraction=1.0 ) strategies = [baseline, strategy1, strategy2, strategy3, strategy4] fig, ax = plt.subplots(figsize=(12, 7)) names = [s.scenario_name for s in strategies] co2_vals = [s.co2_total_Mt for s in strategies] reductions = [(baseline.co2_total_Mt - s.co2_total_Mt) / baseline.co2_total_Mt * 100 for s in strategies] x = np.arange(len(names)) colors = ['#FF6B6B', '#FFE66D', '#4ECDC4', '#45B7D1', '#96CEB4'] bars = ax.bar(x, co2_vals, color=colors) ax.set_ylabel('Total CO₂ Emissions (Mt)', fontsize=12) ax.set_title('CO₂ Reduction Strategies Comparison\n(Balanced baseline: 139 years)', fontsize=14) ax.set_xticks(x) ax.set_xticklabels(names, rotation=15, ha='right') ax.grid(True, alpha=0.3, axis='y') # 添加数值和减排比例 for i, (bar, val, red) in enumerate(zip(bars, co2_vals, reductions)): ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 20, f'{val:.0f} Mt\n({red:+.1f}%)', ha='center', fontsize=10) # 添加基准线 ax.axhline(y=baseline.co2_total_Mt, color='red', linestyle='--', alpha=0.5, label=f'Baseline: {baseline.co2_total_Mt:.0f} Mt') ax.legend() plt.tight_layout() plt.savefig(f'{save_dir}/mitigation_strategies.png', dpi=150, bbox_inches='tight') print(f"减排策略对比图已保存至: {save_dir}/mitigation_strategies.png") def plot_comprehensive_summary(save_dir: str = '/Volumes/Files/code/mm/20260130_b/p4'): """绘制综合汇总图""" df = analyze_all_scenarios() fig = plt.figure(figsize=(16, 12)) # 创建4x4网格 gs = fig.add_gridspec(3, 3, hspace=0.3, wspace=0.3) scenarios = df['Scenario'].tolist() # ========== 左上: 雷达图 ========== ax1 = fig.add_subplot(gs[0, 0], projection='polar') categories = ['Time\n(inv)', 'Energy\n(inv)', 'CO2\n(inv)', 'Launches\n(inv)', 'Strat H2O\n(inv)'] N = len(categories) # 归一化数据 (取逆以使"更好"在外围) max_vals = df[['Years', 'Energy (PJ)', 'CO2 Total (Mt)', 'Launches', 'H2O Stratosphere (Mt)']].max() angles = np.linspace(0, 2*np.pi, N, endpoint=False).tolist() angles += angles[:1] colors = ['#FF6B6B', '#4ECDC4', '#45B7D1', '#96CEB4'] for i, (idx, row) in enumerate(df.iterrows()): values = [ 1 - row['Years'] / max_vals['Years'], 1 - row['Energy (PJ)'] / max_vals['Energy (PJ)'], 1 - row['CO2 Total (Mt)'] / max_vals['CO2 Total (Mt)'], 1 - row['Launches'] / max_vals['Launches'] if max_vals['Launches'] > 0 else 1, 1 - row['H2O Stratosphere (Mt)'] / max_vals['H2O Stratosphere (Mt)'] if max_vals['H2O Stratosphere (Mt)'] > 0 else 1 ] values += values[:1] ax1.plot(angles, values, 'o-', linewidth=2, label=row['Scenario'], color=colors[i]) ax1.fill(angles, values, alpha=0.15, color=colors[i]) ax1.set_xticks(angles[:-1]) ax1.set_xticklabels(categories, fontsize=9) ax1.set_title('Multi-Criteria Comparison\n(outer = better)', fontsize=12, pad=20) ax1.legend(loc='upper right', bbox_to_anchor=(1.3, 1.0), fontsize=8) # ========== 右上: 汇总表格 ========== ax2 = fig.add_subplot(gs[0, 1:]) ax2.axis('off') table_data = [] for idx, row in df.iterrows(): table_data.append([ row['Scenario'], f"{row['Years']:.0f}", f"{row['Energy (PJ)']:.0f}", f"{row['CO2 Total (Mt)']:.0f}", f"{row['H2O Stratosphere (Mt)']:.1f}", f"{row['Annual CO2 (Mt/yr)']:.1f}", ]) table = ax2.table( cellText=table_data, colLabels=['Scenario', 'Years', 'Energy(PJ)', 'CO2(Mt)', 'Strat.H2O(Mt)', 'Annual CO2'], loc='center', cellLoc='center' ) table.auto_set_font_size(False) table.set_fontsize(10) table.scale(1.2, 1.8) # 高亮最优值 for i in range(len(df)): for j in range(1, 6): cell = table[(i+1, j)] if j == 1 and df.iloc[i]['Years'] == df['Years'].min(): cell.set_facecolor('#90EE90') elif j == 2 and df.iloc[i]['Energy (PJ)'] == df['Energy (PJ)'].min(): cell.set_facecolor('#90EE90') elif j == 3 and df.iloc[i]['CO2 Total (Mt)'] == df['CO2 Total (Mt)'].min(): cell.set_facecolor('#90EE90') ax2.set_title('Scenario Summary Table\n(Green = Best in Category)', fontsize=12, pad=10) # ========== 中间行: CO2分解 + 时间线 ========== ax3 = fig.add_subplot(gs[1, :2]) x = np.arange(len(scenarios)) width = 0.25 ax3.bar(x - width, df['CO2 Combustion (Mt)'], width, label='Combustion', color='#FF6B6B') ax3.bar(x, df['CO2 Production (Mt)'], width, label='Fuel Production', color='#4ECDC4') ax3.bar(x + width, df['CO2 Electricity (Mt)'], width, label='Electricity', color='#45B7D1') ax3.set_ylabel('CO₂ (Mt)', fontsize=11) ax3.set_title('CO₂ Emissions by Source', fontsize=12) ax3.set_xticks(x) ax3.set_xticklabels(scenarios, rotation=15, ha='right') ax3.legend() ax3.grid(True, alpha=0.3, axis='y') # ========== 右侧: 时间-CO2权衡 ========== ax4 = fig.add_subplot(gs[1, 2]) ax4.scatter(df['Years'], df['CO2 Total (Mt)'], s=200, c=colors, edgecolors='black', linewidth=2) for i, row in df.iterrows(): ax4.annotate(row['Scenario'].replace(' ', '\n'), (row['Years'], row['CO2 Total (Mt)']), textcoords="offset points", xytext=(0, 15), ha='center', fontsize=8) ax4.set_xlabel('Completion Time (years)', fontsize=11) ax4.set_ylabel('Total CO₂ (Mt)', fontsize=11) ax4.set_title('Time-CO₂ Trade-off', fontsize=12) ax4.grid(True, alpha=0.3) # ========== 底行: 结论文字 ========== ax5 = fig.add_subplot(gs[2, :]) ax5.axis('off') conclusions = """ KEY FINDINGS (LOX/CH4 Fuel): 1. Rocket Only: Highest CO₂ emissions ({:.0f} Mt) due to large fuel consumption ({:.1f} Mt fuel). - Annual emissions: {:.1f} Mt/yr (≈{:.3f}% of global annual emissions) - Stratospheric H₂O injection: {:.1f} Mt (potential climate impact) 2. Elevator Only: Lowest CO₂ emissions ({:.0f} Mt), but longest timeline ({:.0f} years). - CO₂ reduction vs Rocket: {:.0f}% (primary savings from avoiding rocket fuel production) - No stratospheric H₂O injection from main transport 3. Combined (Balanced): Best trade-off at 139 years with {:.0f} Mt CO₂. - 74.6% payload via elevator, 25.4% via low-latitude rockets - Reasonable timeline with significant environmental benefit 4. RECOMMENDATIONS for Minimizing Environmental Impact: a) Prioritize space elevator transport (saves ~{:.0f}% CO₂ vs pure rocket) b) Use renewable energy for elevator operation (reduces electricity emissions) c) If rockets needed, prefer low-latitude sites (lower fuel per kg) d) Consider carbon capture for CH4 production phase """.format( df.loc[df['Scenario']=='Rocket Only', 'CO2 Total (Mt)'].values[0], df.loc[df['Scenario']=='Rocket Only', 'Fuel (Mt)'].values[0], df.loc[df['Scenario']=='Rocket Only', 'Annual CO2 (Mt/yr)'].values[0], df.loc[df['Scenario']=='Rocket Only', 'Annual CO2 (Mt/yr)'].values[0] / 37000 * 100, df.loc[df['Scenario']=='Rocket Only', 'H2O Stratosphere (Mt)'].values[0], df.loc[df['Scenario']=='Elevator Only', 'CO2 Total (Mt)'].values[0], df.loc[df['Scenario']=='Elevator Only', 'Years'].values[0], (1 - df.loc[df['Scenario']=='Elevator Only', 'CO2 Total (Mt)'].values[0] / df.loc[df['Scenario']=='Rocket Only', 'CO2 Total (Mt)'].values[0]) * 100, df.loc[df['Scenario']=='Combined (Balanced)', 'CO2 Total (Mt)'].values[0], (1 - df.loc[df['Scenario']=='Elevator Only', 'CO2 Total (Mt)'].values[0] / df.loc[df['Scenario']=='Rocket Only', 'CO2 Total (Mt)'].values[0]) * 100 ) ax5.text(0.02, 0.95, conclusions, transform=ax5.transAxes, fontsize=10, verticalalignment='top', fontfamily='monospace', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5)) plt.suptitle('Moon Colony Environmental Impact Analysis\n' '(100 Million Tons Payload, LOX/CH4 Fuel)', fontsize=14, y=0.98) plt.savefig(f'{save_dir}/comprehensive_summary.png', dpi=150, bbox_inches='tight') print(f"综合汇总图已保存至: {save_dir}/comprehensive_summary.png") # ============== 报告生成 ============== def generate_report(save_dir: str = '/Volumes/Files/code/mm/20260130_b/p4'): """生成分析报告""" df = analyze_all_scenarios() # 优化结果 opt_result = optimize_for_environment( max_annual_co2_Mt=50.0, min_completion_years=120.0, renewable_fraction=0.5 ) report = [] report.append("=" * 80) report.append("TASK 4: ENVIRONMENTAL IMPACT ANALYSIS") report.append("Moon Colony Construction - 100 Million Metric Tons") report.append("Fuel: LOX/CH4 (Liquid Oxygen / Methane)") report.append("=" * 80) report.append("\n" + "=" * 80) report.append("1. FUEL PARAMETERS (LOX/CH4)") report.append("=" * 80) report.append(f" Specific Impulse (Isp): {FUEL_PARAMS.isp} s") report.append(f" Exhaust Velocity: {FUEL_PARAMS.exhaust_velocity:.0f} m/s") report.append(f" Specific Energy: {FUEL_PARAMS.specific_energy/1e6:.1f} MJ/kg") report.append(f" O2/CH4 Mass Ratio: {FUEL_PARAMS.ox_fuel_ratio}:1") report.append(f" CH4 Mass Fraction: {FUEL_PARAMS.ch4_mass_fraction*100:.1f}%") report.append(f"\n Combustion Reaction: CH4 + 2O2 → CO2 + 2H2O") report.append(f" CO2 per kg fuel: {EMISSION_FACTORS.co2_per_kg_fuel:.3f} kg") report.append(f" H2O per kg fuel: {EMISSION_FACTORS.h2o_per_kg_fuel:.3f} kg") report.append("\n" + "=" * 80) report.append("2. SCENARIO COMPARISON") report.append("=" * 80) for idx, row in df.iterrows(): report.append(f"\n {row['Scenario']}:") report.append(f" Completion Time: {row['Years']:.0f} years") report.append(f" Total Energy: {row['Energy (PJ)']:.0f} PJ") report.append(f" Total Fuel: {row['Fuel (Mt)']:.2f} Mt") report.append(f" CO2 Emissions: {row['CO2 Total (Mt)']:.0f} Mt") report.append(f" - Combustion: {row['CO2 Combustion (Mt)']:.0f} Mt") report.append(f" - Fuel Production: {row['CO2 Production (Mt)']:.0f} Mt") report.append(f" - Electricity: {row['CO2 Electricity (Mt)']:.1f} Mt") report.append(f" H2O Emissions: {row['H2O Total (Mt)']:.2f} Mt") report.append(f" - Stratospheric: {row['H2O Stratosphere (Mt)']:.2f} Mt") report.append(f" Annual CO2: {row['Annual CO2 (Mt/yr)']:.2f} Mt/yr ({row['Annual CO2 (Mt/yr)']/370:.3f}% of global)") report.append(f" CO2 per ton payload: {row['CO2/ton (kg)']:.0f} kg/ton") report.append("\n" + "=" * 80) report.append("3. ENVIRONMENTAL IMPACT RANKING") report.append("=" * 80) df_sorted = df.sort_values('CO2 Total (Mt)') report.append("\n By Total CO2 (lowest first):") for i, (idx, row) in enumerate(df_sorted.iterrows(), 1): report.append(f" {i}. {row['Scenario']}: {row['CO2 Total (Mt)']:.0f} Mt") rocket_co2 = df.loc[df['Scenario']=='Rocket Only', 'CO2 Total (Mt)'].values[0] report.append("\n CO2 Reduction vs Rocket Only:") for idx, row in df.iterrows(): if row['Scenario'] != 'Rocket Only': reduction = (rocket_co2 - row['CO2 Total (Mt)']) / rocket_co2 * 100 report.append(f" {row['Scenario']}: -{reduction:.1f}%") report.append("\n" + "=" * 80) report.append("4. OPTIMIZATION FOR MINIMUM ENVIRONMENTAL IMPACT") report.append("=" * 80) if opt_result: report.append(f"\n Optimization Constraints:") report.append(f" - Minimum completion time: 120 years") report.append(f" - Maximum annual CO2: 50 Mt/yr") report.append(f" - Renewable energy fraction: 50%") report.append(f"\n Optimal Solution:") report.append(f" - Elevator fraction: {opt_result['elevator_fraction']*100:.1f}%") report.append(f" - Rocket sites used: {opt_result['n_sites']}") report.append(f" - Completion time: {opt_result['years']:.0f} years") report.append(f" - Total CO2: {opt_result['impact'].co2_total_Mt:.0f} Mt") report.append(f" - Annual CO2: {opt_result['impact'].annual_co2_Mt:.2f} Mt/yr") report.append("\n" + "=" * 80) report.append("5. MODEL ADJUSTMENT RECOMMENDATIONS") report.append("=" * 80) report.append(""" To minimize environmental impact, adjust the model as follows: A. OBJECTIVE FUNCTION MODIFICATION: Original: min(α×Time + β×Energy) Modified: min(α×Time + β×Energy + γ×CO2_total + δ×H2O_stratosphere) Where γ and δ are weights for environmental costs. B. RECOMMENDED STRATEGIES (by effectiveness): 1. Prioritize Elevator Transport (CO2 reduction: ~70-80%) - Elevator has no direct combustion emissions - Top rocket stage emissions much smaller than full rocket launch 2. Use Renewable Energy for Elevator (CO2 reduction: 5-10%) - Replace grid electricity with solar/wind - Reduces electricity-related emissions 3. Low-Latitude Launch Sites Only (CO2 reduction: 3-5%) - Lower ΔV requirements = less fuel - Kourou (5.2°) vs Alaska (57.4°): ~10% fuel difference 4. Green Methane Production (CO2 reduction: 15-20%) - Use bio-methane or synthetic methane from CO2 capture - Reduces fuel production emissions 5. Extend Timeline (CO2 reduction: proportional) - More time → more elevator use → less rocket use - Trade-off: slower colony establishment C. QUANTIFIED IMPACT: - Baseline (Rocket Only): {:.0f} Mt CO2 - Best Case (Elevator + 100% Renewable): {:.0f} Mt CO2 - Maximum Reduction: {:.1f}% """.format( df.loc[df['Scenario']=='Rocket Only', 'CO2 Total (Mt)'].values[0], calculate_scenario_impact('Best', 186, 1.0, [], renewable_energy_fraction=1.0).co2_total_Mt, (1 - calculate_scenario_impact('Best', 186, 1.0, [], renewable_energy_fraction=1.0).co2_total_Mt / df.loc[df['Scenario']=='Rocket Only', 'CO2 Total (Mt)'].values[0]) * 100 )) report.append("\n" + "=" * 80) report.append("6. STRATOSPHERIC IMPACT ANALYSIS") report.append("=" * 80) report.append(f""" Stratospheric Water Vapor Injection (Rocket Only scenario): - Total H2O emitted: {df.loc[df['Scenario']=='Rocket Only', 'H2O Total (Mt)'].values[0]:.2f} Mt - Stratospheric fraction: {EMISSION_FACTORS.stratosphere_emission_fraction*100:.0f}% - Stratospheric H2O: {df.loc[df['Scenario']=='Rocket Only', 'H2O Stratosphere (Mt)'].values[0]:.2f} Mt - Residence time: {EMISSION_FACTORS.stratosphere_h2o_residence_time:.0f} years Potential Climate Effects: - Stratospheric H2O has greenhouse effect (radiative forcing) - May affect ozone chemistry - Elevator-only scenario eliminates this impact entirely """) report.append("=" * 80) # 保存报告 report_text = '\n'.join(report) with open(f'{save_dir}/environmental_report.txt', 'w', encoding='utf-8') as f: f.write(report_text) print(f"分析报告已保存至: {save_dir}/environmental_report.txt") # 保存CSV df.to_csv(f'{save_dir}/scenario_comparison.csv', index=False) print(f"数据表已保存至: {save_dir}/scenario_comparison.csv") return report_text # ============== 主程序 ============== if __name__ == "__main__": print("=" * 80) print("TASK 4: ENVIRONMENTAL IMPACT ANALYSIS") print("Moon Colony Construction with LOX/CH4 Fuel") print("=" * 80) save_dir = '/Volumes/Files/code/mm/20260130_b/p4' # 1. 基础情景对比 print("\n1. Analyzing baseline scenarios...") df = plot_environmental_comparison(save_dir) print("\nScenario Comparison:") print(df.to_string(index=False)) # 2. Pareto分析 print("\n2. Generating Pareto analysis...") plot_pareto_with_environment(save_dir) # 3. 减排策略对比 print("\n3. Analyzing mitigation strategies...") plot_mitigation_strategies(save_dir) # 4. 综合汇总 print("\n4. Generating comprehensive summary...") plot_comprehensive_summary(save_dir) # 5. 生成报告 print("\n5. Generating report...") generate_report(save_dir) print("\n" + "=" * 80) print("Analysis complete! All outputs saved to:", save_dir) print("=" * 80)