diff --git a/p1/__pycache__/combined_scenario.cpython-313.pyc b/p1/__pycache__/combined_scenario.cpython-313.pyc new file mode 100644 index 0000000..83b5917 Binary files /dev/null and b/p1/__pycache__/combined_scenario.cpython-313.pyc differ diff --git a/p1/__pycache__/pareto_optimization.cpython-313.pyc b/p1/__pycache__/pareto_optimization.cpython-313.pyc new file mode 100644 index 0000000..fdd6a6d Binary files /dev/null and b/p1/__pycache__/pareto_optimization.cpython-313.pyc differ diff --git a/p1/__pycache__/sensitivity_analysis.cpython-313.pyc b/p1/__pycache__/sensitivity_analysis.cpython-313.pyc new file mode 100644 index 0000000..4ff62ca Binary files /dev/null and b/p1/__pycache__/sensitivity_analysis.cpython-313.pyc differ diff --git a/p1/combined_scenario.py b/p1/combined_scenario.py index 6f8ac3f..93d774d 100644 --- a/p1/combined_scenario.py +++ b/p1/combined_scenario.py @@ -515,7 +515,7 @@ def print_summary(df: pd.DataFrame): def plot_three_scenarios_comparison( year_min: float = 100, - year_max: float = 300, + year_max: float = 250, save_path: str = '/Volumes/Files/code/mm/20260130_b/three_scenarios_comparison.png' ): """ @@ -579,74 +579,45 @@ def plot_three_scenarios_comparison( else: combined_energy.append(np.nan) - # 创建图表 - fig, ax = plt.subplots(figsize=(14, 8)) + # 创建图表 - 使用0.618黄金比例，缩小尺寸以放大字体 + fig_width = 8 + fig_height = fig_width * 0.618 + fig, ax = plt.subplots(figsize=(fig_width, fig_height)) + + # 中偏低饱和度配色 + color_rocket = '#B85450' # 暗砖红色 + color_elevator = '#5B8266' # 灰绿色 + color_combined = '#4A6FA5' # 灰蓝色 # 绘制三条折线 - ax.plot(years, rocket_energy, 'r-', linewidth=2.5, label='Rocket Only', marker='', markersize=4) - ax.plot(years, elevator_energy, 'g-', linewidth=2.5, label='Elevator Only', marker='', markersize=4) - ax.plot(years, combined_energy, 'b-', linewidth=2.5, label='Combined (Elevator + Rocket)', marker='', markersize=4) + ax.plot(years, rocket_energy, color=color_rocket, linewidth=2.5, label='Rocket Only', marker='', markersize=4) + ax.plot(years, elevator_energy, color=color_elevator, linewidth=2.5, label='Elevator Only', marker='', markersize=4) + ax.plot(years, combined_energy, color=color_combined, linewidth=2.5, label='Combined (Elevator + Rocket)', marker='', markersize=4) # 标记关键时间点 - ax.axvline(x=combined_min_years, color='blue', linestyle=':', alpha=0.7, linewidth=1.5) - ax.axvline(x=elevator_only_min_years, color='green', linestyle=':', alpha=0.7, linewidth=1.5) - ax.axvline(x=rocket_only_min_years, color='red', linestyle=':', alpha=0.7, linewidth=1.5) + ax.axvline(x=combined_min_years, color=color_combined, linestyle=':', alpha=0.7, linewidth=1.5) + ax.axvline(x=elevator_only_min_years, color=color_elevator, linestyle=':', alpha=0.7, linewidth=1.5) + ax.axvline(x=rocket_only_min_years, color=color_rocket, linestyle=':', alpha=0.7, linewidth=1.5) - # 添加标注 - y_pos = ax.get_ylim()[1] * 0.95 + # 添加标注 - 使用低饱和度背景色，避免遮挡曲线 + y_max = ax.get_ylim()[1] ax.annotate(f'Combined\nMin: {combined_min_years:.0f}y', - xy=(combined_min_years, y_pos), fontsize=9, ha='center', - bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.8)) + xy=(combined_min_years, y_max * 0.55), fontsize=9, ha='center', + bbox=dict(boxstyle='round', facecolor='#C5D5E8', edgecolor=color_combined, alpha=0.9)) ax.annotate(f'Elevator\nMin: {elevator_only_min_years:.0f}y', - xy=(elevator_only_min_years, y_pos * 0.85), fontsize=9, ha='center', - bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.8)) + xy=(elevator_only_min_years, y_max * 0.72), fontsize=9, ha='center', + bbox=dict(boxstyle='round', facecolor='#C8DBC8', edgecolor=color_elevator, alpha=0.9)) ax.annotate(f'Rocket\nMin: {rocket_only_min_years:.0f}y', - xy=(rocket_only_min_years, y_pos * 0.75), fontsize=9, ha='center', - bbox=dict(boxstyle='round', facecolor='lightsalmon', alpha=0.8)) + xy=(rocket_only_min_years, y_max * 0.55), fontsize=9, ha='center', + bbox=dict(boxstyle='round', facecolor='#E8C8C5', edgecolor=color_rocket, alpha=0.9)) - # 填充区域表示节能效果 - # Combined比Rocket节省的能量 - rocket_arr = np.array(rocket_energy) - combined_arr = np.array(combined_energy) - valid_mask = ~np.isnan(rocket_arr) & ~np.isnan(combined_arr) - ax.fill_between(years[valid_mask], combined_arr[valid_mask], rocket_arr[valid_mask], - alpha=0.2, color='purple', label='Energy saved by Combined vs Rocket') - - ax.set_xlabel('Completion Time (years)', fontsize=13) - ax.set_ylabel('Total Energy Consumption (PJ)', fontsize=13) - ax.set_title('Moon Colony Construction: Energy vs Completion Time\n' - '(100 Million Metric Tons Payload Comparison)', fontsize=15) - ax.legend(loc='upper right', fontsize=11) + ax.set_xlabel('Completion Time (years)', fontsize=11) + ax.set_ylabel('Total Energy Consumption (PJ)', fontsize=11) + ax.legend(loc='upper left', fontsize=9) ax.grid(True, alpha=0.3) ax.set_xlim(year_min, year_max) - # 添加数据表格 - # 选择几个关键年份的数据 - key_years = [110, 150, 186, 220, 250] - table_data = [] - for ky in key_years: - idx = np.argmin(np.abs(years - ky)) - table_data.append([ - f'{ky}', - f'{combined_energy[idx]:.0f}' if not np.isnan(combined_energy[idx]) else 'N/A', - f'{elevator_energy[idx]:.0f}' if not np.isnan(elevator_energy[idx]) else 'N/A', - f'{rocket_energy[idx]:.0f}' if not np.isnan(rocket_energy[idx]) else 'N/A', - ]) - - # 在图表下方添加表格 - table = ax.table( - cellText=table_data, - colLabels=['Year', 'Combined (PJ)', 'Elevator (PJ)', 'Rocket (PJ)'], - loc='bottom', - bbox=[0.15, -0.25, 0.7, 0.18], - cellLoc='center' - ) - table.auto_set_font_size(False) - table.set_fontsize(10) - table.scale(1, 1.5) - - # 调整布局给表格留空间 - plt.subplots_adjust(bottom=0.25) + plt.tight_layout() plt.savefig(save_path, dpi=150, bbox_inches='tight') print(f"\n三方案对比图已保存至: {save_path}") @@ -700,6 +671,6 @@ if __name__ == "__main__": print(f"{name:<25} {data['years']:>12.0f} {data['energy_PJ']:>15.0f} {savings:>17.0f}%") print("=" * 90) - # 生成三方案对比折线图 (100-300年) - print("\n\n正在生成三方案对比折线图 (100-300年)...") - plot_three_scenarios_comparison(year_min=100, year_max=300) + # 生成三方案对比折线图 (100-250年) + print("\n\n正在生成三方案对比折线图 (100-250年)...") + plot_three_scenarios_comparison(year_min=100, year_max=250) diff --git a/p1/latitude_effects.png b/p1/latitude_effects.png index 9655f7d..2dea388 100644 Binary files a/p1/latitude_effects.png and b/p1/latitude_effects.png differ diff --git a/p1/marginal_benefit.png b/p1/marginal_benefit.png new file mode 100644 index 0000000..5b9a960 Binary files /dev/null and b/p1/marginal_benefit.png differ diff --git a/p1/pareto_optimization.py b/p1/pareto_optimization.py index 579c28a..594e247 100644 --- a/p1/pareto_optimization.py +++ b/p1/pareto_optimization.py @@ -796,6 +796,60 @@ def plot_combined_range_analysis( return fig +def plot_marginal_benefit_only( + df: pd.DataFrame, + knee_analysis: Dict, + save_path: str = '/Volumes/Files/code/mm/20260130_b/p1/marginal_benefit.png' +): + """ + 仅绘制边际收益图（右下角图的简化版） + - 长宽比 0.618 + - 中低饱和度色系 + - 无标题 + - 无 Elevator share 线 + - Y轴范围 450-600 + """ + # 使用0.618黄金比例，缩小尺寸以放大字体 + fig_width = 8 + fig_height = fig_width * 0.618 + fig, ax = plt.subplots(figsize=(fig_width, fig_height)) + + T = df['years'].values + E = df['energy_PJ'].values + + # 边际节省 + dE_dT = np.gradient(E, T) + marginal_savings = -dE_dT # 每多1年节省的能量 + + # 中低饱和度色系 + color_main = '#4A6FA5' # 灰蓝色 + color_fill = '#A8C0D8' # 浅灰蓝 + color_vline = '#B85450' # 暗砖红 + + # 边际收益曲线 + ax.plot(T, marginal_savings, color=color_main, linewidth=2.5, label='Marginal Saving (PJ/year)') + ax.fill_between(T, marginal_savings, alpha=0.25, color=color_fill) + + # 标记推荐点 + rec = knee_analysis['recommended'] + ax.axvline(x=rec['years'], color=color_vline, linestyle=':', linewidth=2, + label=f'Recommended: {rec["years"]:.0f}y') + + ax.set_xlabel('Completion Time (years)', fontsize=11) + ax.set_ylabel('Marginal Energy Saving (PJ/year)', fontsize=11) + ax.legend(loc='upper right', fontsize=9) + ax.grid(True, alpha=0.3) + + # 设置Y轴范围 450-600 + ax.set_ylim(450, 600) + + plt.tight_layout() + plt.savefig(save_path, dpi=150, bbox_inches='tight') + print(f"边际收益图已保存至: {save_path}") + + return fig + + def plot_combined_decision( df: pd.DataFrame, knee_analysis: Dict, diff --git a/p1/richards_curve_2010.py b/p1/richards_curve_2010.py index c1b4102..c2b649d 100644 --- a/p1/richards_curve_2010.py +++ b/p1/richards_curve_2010.py @@ -103,43 +103,49 @@ def main(): print(f"K / Physical limit = {K_fixed/physical_max:.2f}x") # ========== Create Visualization ========== - fig, ax = plt.subplots(figsize=(12, 7)) + # 缩小图片尺寸（比例不变），使字体相对更大 + fig, ax = plt.subplots(figsize=(8, 4.67)) + + # 中低饱和度配色 + color_data = '#5D6D7E' # 灰蓝色 - 数据点 + color_curve = '#52796F' # 暗绿色 - S曲线 + color_target = '#7B9EA8' # 灰蓝绿色 - 2050标记 # Historical data points - ax.scatter(years, launches, color='#2C3E50', s=100, alpha=0.9, + ax.scatter(years, launches, color=color_data, s=80, alpha=0.85, label='Historical Data (2010-2025)', zorder=4, edgecolor='white', linewidth=1) # Generate smooth S-curve - years_smooth = np.linspace(start_year, 2080, 500) + years_smooth = np.linspace(start_year, 2055, 500) t_smooth = years_smooth - start_year pred_smooth = richards(t_smooth, K_fixed, r, t0, v) # S-curve prediction - ax.plot(years_smooth, pred_smooth, color='#27AE60', lw=3, + ax.plot(years_smooth, pred_smooth, color=color_curve, lw=2.5, label=f'Richards Model (K={K_fixed}, R²={r_squared:.3f})', zorder=2) # K=4298 saturation line - ax.axhline(K_fixed, color='#27AE60', ls=':', lw=2, alpha=0.7, + ax.axhline(K_fixed, color=color_curve, ls=':', lw=1.5, alpha=0.6, label=f'K = {K_fixed}') # Mark 2050 line only - ax.axvline(2050, color='#3498DB', ls=':', lw=2, alpha=0.8) - ax.text(2051, K_fixed*0.85, '2050\n(Target)', fontsize=10, color='#3498DB') + ax.axvline(2050, color=color_target, ls=':', lw=1.5, alpha=0.7) + ax.text(2050.5, K_fixed*0.83, '2050\n(Target)', fontsize=9, color=color_target) # Only show 2050 prediction point t_2050 = 2050 - start_year pred_2050 = richards(t_2050, K_fixed, r, t0, v) - ax.scatter([2050], [pred_2050], color='#3498DB', s=80, marker='D', zorder=4) + ax.scatter([2050], [pred_2050], color=color_target, s=60, marker='D', zorder=4) ax.annotate(f'{pred_2050:.0f}', xy=(2050, pred_2050), - xytext=(2051, pred_2050+150), - fontsize=9, color='#3498DB', fontweight='bold') + xytext=(2050.5, pred_2050+180), + fontsize=9, color=color_target, fontweight='bold') # Formatting - ax.set_xlabel('Year', fontsize=12) - ax.set_ylabel('Annual Launches', fontsize=12) - ax.legend(loc='upper left', fontsize=10) - ax.grid(True, alpha=0.3) - ax.set_xlim(2008, 2080) + ax.set_xlabel('Year', fontsize=11) + ax.set_ylabel('Annual Launches', fontsize=11) + ax.legend(loc='upper left', fontsize=9) + ax.grid(True, alpha=0.25) + ax.set_xlim(2010, 2055) ax.set_ylim(0, K_fixed * 1.15) plt.tight_layout() diff --git a/p1/richards_curve_2010_K4298.png b/p1/richards_curve_2010_K4298.png index 6840316..4c2ad33 100644 Binary files a/p1/richards_curve_2010_K4298.png and b/p1/richards_curve_2010_K4298.png differ diff --git a/p1/sensitivity_analysis.py b/p1/sensitivity_analysis.py index f4d4938..9e72063 100644 --- a/p1/sensitivity_analysis.py +++ b/p1/sensitivity_analysis.py @@ -585,6 +585,96 @@ def plot_sensitivity_analysis( return fig +def plot_tornado_chart_standalone( + payload_df: pd.DataFrame, + elevator_df: pd.DataFrame, + frequency_df: pd.DataFrame, + alpha_df: pd.DataFrame, + save_path: str = '/Volumes/Files/code/mm/20260130_b/p1/tornado_chart.png' +): + """ + 单独绘制 Tornado Chart（敏感性汇总图） + - XY轴反转（垂直柱状图） + - 长宽比 0.618 + - 无标题 + - 有图例 + - 低饱和度色系 + """ + # 使用0.618黄金比例 + fig_width = 8 + fig_height = fig_width * 0.618 + fig, ax = plt.subplots(figsize=(fig_width, fig_height)) + + # 计算敏感性指标 + baseline_T = 139 + + sensitivities = [ + ('Payload\n(100→150t)', + (payload_df['optimal_T_lambda504'].iloc[0] - baseline_T) / baseline_T * 100, + (payload_df['optimal_T_lambda504'].iloc[-1] - baseline_T) / baseline_T * 100), + ('Elevator Cap.\n(80%→120%)', + (elevator_df['optimal_T_lambda504'].iloc[0] - baseline_T) / baseline_T * 100, + (elevator_df['optimal_T_lambda504'].iloc[-1] - baseline_T) / baseline_T * 100), + ('Launch Freq.\n(0.5→2/day)', + (frequency_df['optimal_T_lambda504'].iloc[0] - baseline_T) / baseline_T * 100, + (frequency_df['optimal_T_lambda504'].iloc[-1] - baseline_T) / baseline_T * 100), + ('Struct. Coef.\n(0.06→0.14)', + (alpha_df['optimal_T_lambda504'].iloc[0] - baseline_T) / baseline_T * 100, + (alpha_df['optimal_T_lambda504'].iloc[-1] - baseline_T) / baseline_T * 100), + ] + + # 低饱和度配色 + color_positive = '#7A9CBF' # 灰蓝色（参数增加导致T*增加） + color_negative = '#C4816B' # 灰橙红色（参数增加导致T*减少） + + x_pos = np.arange(len(sensitivities)) + bar_width = 0.6 + + # 绘制柱状图（XY反转，使用垂直柱状图） + for i, (name, low, high) in enumerate(sensitivities): + bar_color = color_positive if high > low else color_negative + bar_height = high - low + bar_bottom = min(low, high) + + ax.bar(i, bar_height, bottom=bar_bottom, width=bar_width, + color=bar_color, alpha=0.8, edgecolor='#555555', linewidth=0.8) + + # 标记端点 + ax.plot(i, low, 'o', color='#444444', markersize=7, zorder=5) + ax.plot(i, high, '^', color='#444444', markersize=7, zorder=5) + ax.plot([i, i], [low, high], '-', color='#444444', linewidth=1.5, zorder=4) + + # 基准线 + ax.axhline(y=0, color='#333333', linestyle='-', linewidth=1.2) + + # 设置轴 + ax.set_xticks(x_pos) + ax.set_xticklabels([s[0] for s in sensitivities], fontsize=10) + ax.set_ylabel('Change in Optimal Timeline T* (%)', fontsize=11) + ax.grid(True, alpha=0.3, axis='y') + + # 添加图例 + from matplotlib.patches import Patch + from matplotlib.lines import Line2D + legend_elements = [ + Patch(facecolor=color_negative, edgecolor='#555555', alpha=0.8, + label='Parameter ↑ → T* ↓'), + Patch(facecolor=color_positive, edgecolor='#555555', alpha=0.8, + label='Parameter ↑ → T* ↑'), + Line2D([0], [0], marker='o', color='#444444', linestyle='None', + markersize=6, label='Low value'), + Line2D([0], [0], marker='^', color='#444444', linestyle='None', + markersize=6, label='High value'), + ] + ax.legend(handles=legend_elements, loc='upper right', fontsize=9, framealpha=0.9) + + plt.tight_layout() + plt.savefig(save_path, dpi=150, bbox_inches='tight') + print(f"Tornado chart saved to: {save_path}") + + return fig + + def plot_tradeoff_comparison( save_path: str = '/Volumes/Files/code/mm/20260130_b/p1/sensitivity_tradeoff.png' ): diff --git a/p1/specific_energy_comparison.py b/p1/specific_energy_comparison.py index e4b38da..632764e 100644 --- a/p1/specific_energy_comparison.py +++ b/p1/specific_energy_comparison.py @@ -79,29 +79,29 @@ class LaunchSite: # 预定义的发射场 LAUNCH_SITES = { # 赤道参考 - '赤道': LaunchSite("Equator (参考)", 0.0), + '赤道': LaunchSite("Equator", 0.0), # 法属圭亚那 (接近赤道) - 'Kourou': LaunchSite("Kourou (French Guiana)", 5.2), + 'Kourou': LaunchSite("French Guiana", 5.2), # 印度 - 'SDSC': LaunchSite("Satish Dhawan (India)", 13.7), + 'SDSC': LaunchSite("Satish Dhawan", 13.7), # 美国 - 'Texas': LaunchSite("Boca Chica (Texas)", 26.0), - 'Florida': LaunchSite("Cape Canaveral (Florida)", 28.5), - 'California': LaunchSite("Vandenberg (California)", 34.7), - 'Virginia': LaunchSite("Wallops (Virginia)", 37.8), - 'Alaska': LaunchSite("Kodiak (Alaska)", 57.4), + 'Texas': LaunchSite("Texas", 26.0), + 'Florida': LaunchSite("Florida", 28.5), + 'California': LaunchSite("California", 34.7), + 'Virginia': LaunchSite("Virginia", 37.8), + 'Alaska': LaunchSite("Alaska", 57.4), # 中国 - 'Taiyuan': LaunchSite("Taiyuan (China)", 38.8), + 'Taiyuan': LaunchSite("Taiyuan", 38.8), # 新西兰 - 'Mahia': LaunchSite("Mahia (New Zealand)", -39.3), # 南半球 + 'Mahia': LaunchSite("Mahia", -39.3), # 南半球 # 哈萨克斯坦 - 'Baikonur': LaunchSite("Baikonur (Kazakhstan)", 45.6), + 'Baikonur': LaunchSite("Kazakhstan", 45.6), } @@ -535,8 +535,12 @@ def plot_latitude_effects( mission_plot = mission title_suffix = " (Target: Equatorial Orbit)" - # Golden ratio aspect: width=10, height=10*0.618 - fig, ax = plt.subplots(figsize=(10, 10 * 0.618)) + # Golden ratio aspect: 缩小尺寸使字体相对更大 + fig, ax = plt.subplots(figsize=(8, 8 * 0.618)) + + # 中低饱和度配色 + color_curve = '#52796F' # 暗绿色 - 曲线 + color_point = '#8E7B9A' # 灰紫色 - 数据点 # Continuous latitude range latitudes = np.linspace(0, 65, 100) @@ -549,29 +553,32 @@ def plot_latitude_effects( fuel_ratios.append(result['fuel_ratio']) # Plot continuous curve - ax.plot(latitudes, fuel_ratios, 'g-', linewidth=2, label='Fuel Ratio vs Latitude') + ax.plot(latitudes, fuel_ratios, color=color_curve, linewidth=2, label='Fuel Ratio vs Latitude') # Mark launch sites with custom offsets to avoid overlap with curve # Offsets: (x, y) in points - positive y moves text up, positive x moves right + # 纬度参考: Equator(0), French Guiana(5.2), Satish Dhawan(13.7), Texas(26), + # Florida(28.5), California(34.7), Virginia(37.8), Taiyuan(38.8), + # Mahia(39.3), Kazakhstan(45.6), Alaska(57.4) label_offsets = { '赤道': (5, -15), # Below curve - 'Kourou': (5, 8), # Above curve - 'SDSC': (5, -15), # Below curve - 'Texas': (-50, 8), # Left and above - 'Florida': (5, 8), # Above curve - 'California': (-60, -8), # Left - 'Virginia': (5, 8), # Above curve - 'Taiyuan': (-45, -12), # Left and below - 'Mahia': (5, 8), # Above curve - 'Alaska': (5, -15), # Below curve - 'Baikonur': (5, 8), # Above curve + 'Kourou': (-75, 15), # Left above (French Guiana) + 'SDSC': (8, -15), # Below right (Satish Dhawan) + 'Texas': (-38, -15), # Left below + 'Florida': (8, -15), # Right below + 'California': (8, -15), # Right below + 'Virginia': (8, 12), # Right above + 'Taiyuan': (-52, 12), # Left above (避开California) + 'Mahia': (8, -18), # Right below + 'Alaska': (8, 8), # Right above + 'Baikonur': (8, 8), # Right above (Kazakhstan) } for name, site in LAUNCH_SITES.items(): abs_lat = abs(site.latitude) site_temp = LaunchSite(site.name, abs_lat) result = ground_launch_specific_energy_with_latitude(engine, mission_plot, site_temp, constants) - ax.plot(abs_lat, result['fuel_ratio'], 'mo', markersize=8) + ax.plot(abs_lat, result['fuel_ratio'], 'o', color=color_point, markersize=8) # Simplified label label = site.name.split('(')[0].strip() # Get custom offset for this site, default to (5, 8) @@ -583,7 +590,7 @@ def plot_latitude_effects( ax.set_ylabel('Fuel / Payload Mass Ratio', fontsize=12) # ax.set_title(f'Fuel Requirement vs Launch Latitude{title_suffix}', fontsize=13) ax.grid(True, alpha=0.3) - ax.set_xlim(0, 65) + ax.set_xlim(-2, 65) plt.tight_layout() plt.savefig(save_path, dpi=150, bbox_inches='tight') diff --git a/p1/three_scenarios_comparison.png b/p1/three_scenarios_comparison.png index ff8b9f3..f69bfc7 100644 Binary files a/p1/three_scenarios_comparison.png and b/p1/three_scenarios_comparison.png differ diff --git a/p1/tornado_chart.png b/p1/tornado_chart.png new file mode 100644 index 0000000..7c57808 Binary files /dev/null and b/p1/tornado_chart.png differ diff --git a/three_scenarios_comparison.png b/three_scenarios_comparison.png index 33f0e56..f69bfc7 100644 Binary files a/three_scenarios_comparison.png and b/three_scenarios_comparison.png differ