diff --git a/p4/__pycache__/environmental_impact.cpython-313.pyc b/p4/__pycache__/environmental_impact.cpython-313.pyc index d826e35..df9b05d 100644 Binary files a/p4/__pycache__/environmental_impact.cpython-313.pyc and b/p4/__pycache__/environmental_impact.cpython-313.pyc differ diff --git a/p4/environmental_comparison.png b/p4/environmental_comparison.png index 7bb4668..a3d16ed 100644 Binary files a/p4/environmental_comparison.png and b/p4/environmental_comparison.png differ diff --git a/p4/environmental_impact.py b/p4/environmental_impact.py index d56e5da..b399603 100644 --- a/p4/environmental_impact.py +++ b/p4/environmental_impact.py @@ -490,99 +490,326 @@ def optimize_for_environment( # ============== 可视化 ============== def plot_environmental_comparison(save_dir: str = '/Volumes/Files/code/mm/20260130_b/p4'): - """绘制环境影响对比图""" + """绘制环境影响对比图（2张子图版本，马卡龙色系论文风格）""" df = analyze_all_scenarios() - fig, axes = plt.subplots(2, 2, figsize=(14, 12)) + # 马卡龙色系 + 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.6 + width = 0.35 - # ========== 图1: CO2排放分解 ========== - ax1 = axes[0, 0] + # ========== 左图: 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, co2_combustion, width, label='Combustion', color='#FF6B6B') - bars2 = ax1.bar(x, co2_production, width, bottom=co2_combustion, label='Fuel Production', color='#4ECDC4') - bars3 = ax1.bar(x, co2_electricity, width, bottom=co2_combustion+co2_production, - label='Electricity', color='#45B7D1') + 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)', fontsize=12) - ax1.set_title('CO₂ Emissions Breakdown by Scenario', fontsize=14) - ax1.set_xticks(x) - ax1.set_xticklabels(scenarios, rotation=15, ha='right') - ax1.legend() - ax1.grid(True, alpha=0.3, axis='y') + 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, total + 50, f'{total:.0f}', ha='center', fontsize=10, fontweight='bold') + ax1.text(i - width/2, total + 300, f'{total:.0f}', ha='center', fontsize=9, fontweight='bold') - # ========== 图2: 平流层水蒸气 ========== - ax2 = axes[0, 1] + # H2O数据（右侧柱子） + ax1_twin = ax1.twinx() h2o_surface = df['H2O Total (Mt)'].values - df['H2O Stratosphere (Mt)'].values h2o_strat = df['H2O Stratosphere (Mt)'].values - bars1 = ax2.bar(x, h2o_surface, width, label='Troposphere', color='#96CEB4') - bars2 = ax2.bar(x, h2o_strat, width, bottom=h2o_surface, label='Stratosphere', color='#D4A5A5') + 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) - ax2.set_ylabel('H₂O Emissions (Mt)', fontsize=12) - ax2.set_title('Water Vapor Emissions\n(Stratospheric emissions have climate impact)', fontsize=14) - ax2.set_xticks(x) - ax2.set_xticklabels(scenarios, rotation=15, ha='right') - ax2.legend() - ax2.grid(True, alpha=0.3, axis='y') + ax1_twin.set_ylabel('H₂O Emissions (Mt)', color='#7B68A0') + ax1_twin.tick_params(axis='y', labelcolor='#7B68A0') - # ========== 图3: 年均CO2排放 ========== - ax3 = axes[1, 0] + 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 - colors = ['#FF6B6B', '#4ECDC4', '#45B7D1', '#96CEB4'] - - bars = ax3.bar(x, annual_co2, width, color=colors) - - # 添加参考线：全球年排放量的比例 - global_annual_co2 = 37000 # Mt/yr (2023年全球CO2排放) - for i, val in enumerate(annual_co2): - pct = val / global_annual_co2 * 100 - ax3.text(i, val + 1, f'{val:.1f}\n({pct:.2f}%)', ha='center', fontsize=9) - - ax3.set_ylabel('Annual CO₂ (Mt/yr)', fontsize=12) - ax3.set_title('Annual CO₂ Emissions\n(% of global annual emissions shown)', fontsize=14) - ax3.set_xticks(x) - ax3.set_xticklabels(scenarios, rotation=15, ha='right') - ax3.grid(True, alpha=0.3, axis='y') - - # ========== 图4: 能效对比 (CO2/ton) ========== - ax4 = axes[1, 1] - co2_per_ton = df['CO2/ton (kg)'].values - bars = ax4.bar(x, co2_per_ton, width, color=colors) + # 年均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) - ax4.set_ylabel('CO₂ per ton payload (kg)', fontsize=12) - ax4.set_title('Carbon Intensity: CO₂ Emissions per Ton Payload', fontsize=14) - ax4.set_xticks(x) - ax4.set_xticklabels(scenarios, rotation=15, ha='right') - ax4.grid(True, alpha=0.3, axis='y') + ax2.set_ylabel('Annual CO₂ (Mt/yr)', color='#333333') + ax2.tick_params(axis='y', labelcolor='#333333') - for bar, val in zip(bars, co2_per_ton): - ax4.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 5, - f'{val:.0f}', ha='center', fontsize=10, fontweight='bold') + # 添加年均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=150, bbox_inches='tight') + 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图""" diff --git a/p4/marginal_benefit.png b/p4/marginal_benefit.png new file mode 100644 index 0000000..cd83ea6 Binary files /dev/null and b/p4/marginal_benefit.png differ diff --git a/p4/radar_chart.png b/p4/radar_chart.png new file mode 100644 index 0000000..6075cd5 Binary files /dev/null and b/p4/radar_chart.png differ