diff --git a/p1/richards_curve_2010.py b/p1/richards_curve_2010.py index c2b649d..4e2e0b4 100644 --- a/p1/richards_curve_2010.py +++ b/p1/richards_curve_2010.py @@ -1,160 +1,160 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Richards S-Curve Fit for 2010-2025 Data with K=4298 - -Fits Richards model to 2010-2025 launch data with carrying capacity -constrained to K=4298 (close to physical limit of 3650). -""" - -import pandas as pd -import numpy as np -from scipy.optimize import curve_fit -import matplotlib -matplotlib.use('Agg') -import matplotlib.pyplot as plt -import warnings -warnings.filterwarnings('ignore') - -plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei', 'DejaVu Sans'] -plt.rcParams['axes.unicode_minus'] = False - - -def richards(t, K, r, t0, v): - """Richards curve (generalized logistic model)""" - exp_term = np.exp(-r * (t - t0)) - exp_term = np.clip(exp_term, 1e-10, 1e10) - return K / np.power(1 + exp_term, 1/v) - - -def richards_fixed_K(t, r, t0, v): - """Richards curve with fixed K=4298""" - K = 4298 - exp_term = np.exp(-r * (t - t0)) - exp_term = np.clip(exp_term, 1e-10, 1e10) - return K / np.power(1 + exp_term, 1/v) - - -def load_data(filepath="rocket_launch_counts.csv"): - """Load rocket launch data""" - df = pd.read_csv(filepath) - df = df.rename(columns={"YDate": "year", "Total": "launches"}) - df["year"] = pd.to_numeric(df["year"], errors="coerce") - df["launches"] = pd.to_numeric(df["launches"], errors="coerce") - df = df.dropna(subset=["year", "launches"]) - df = df[(df["year"] >= 1957) & (df["year"] <= 2025)] - df = df.astype({"year": int, "launches": int}) - df = df.sort_values("year").reset_index(drop=True) - return df - - -def main(): - print("=" * 60) - print("Richards S-Curve Fit (2010-2025 Data, K=4298)") - print("=" * 60) - - # Load data - df = load_data() - - # Filter 2010-2025 - start_year = 2010 - end_year = 2025 - data = df[(df["year"] >= start_year) & (df["year"] <= end_year)].copy() - years = data["year"].values - launches = data["launches"].values - - print(f"Data range: {start_year} - {end_year}") - print(f"Data points: {len(data)}") - print(f"\nHistorical data:") - for y, l in zip(years, launches): - print(f" {y}: {l} launches") - - # Fit with fixed K=4298 - K_fixed = 4298 - t = (years - start_year).astype(float) - - p0 = [0.2, 20.0, 2.0] # r, t0, v - bounds = ([0.01, 5, 0.5], [1.0, 100, 10.0]) - - try: - popt, pcov = curve_fit(richards_fixed_K, t, launches, p0=p0, bounds=bounds, maxfev=100000) - r, t0, v = popt - - # Calculate R² - y_pred = richards_fixed_K(t, *popt) - ss_res = np.sum((launches - y_pred) ** 2) - ss_tot = np.sum((launches - np.mean(launches)) ** 2) - r_squared = 1 - (ss_res / ss_tot) - - print(f"\nFitted Parameters (K fixed at {K_fixed}):") - print(f" K (carrying capacity) = {K_fixed} launches/year (FIXED)") - print(f" r (growth rate) = {r:.4f}") - print(f" t0 (inflection point) = {start_year + t0:.1f}") - print(f" v (shape parameter) = {v:.3f}") - print(f" R² = {r_squared:.4f}") - - except Exception as e: - print(f"Fitting error: {e}") - return - - # Physical limit - physical_max = 3650 - print(f"\nPhysical limit: {physical_max} (10 sites × 365 days)") - print(f"K / Physical limit = {K_fixed/physical_max:.2f}x") - - # ========== Create Visualization ========== - # 缩小图片尺寸（比例不变），使字体相对更大 - 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=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, 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=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=color_curve, ls=':', lw=1.5, alpha=0.6, - label=f'K = {K_fixed}') - - # Mark 2050 line only - 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=color_target, s=60, marker='D', zorder=4) - ax.annotate(f'{pred_2050:.0f}', xy=(2050, pred_2050), - xytext=(2050.5, pred_2050+180), - fontsize=9, color=color_target, fontweight='bold') - - # Formatting - 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() - plt.savefig('richards_curve_2010_K4298.png', dpi=150, bbox_inches='tight') - plt.close() - - print("\nPlot saved: richards_curve_2010_K4298.png") - print("=" * 60) - - -if __name__ == "__main__": - main() +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Richards S-Curve Fit for 2010-2025 Data with K=4298 + +Fits Richards model to 2010-2025 launch data with carrying capacity +constrained to K=4298 (close to physical limit of 3650). +""" + +import pandas as pd +import numpy as np +from scipy.optimize import curve_fit +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +import warnings +warnings.filterwarnings('ignore') + +plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei', 'DejaVu Sans'] +plt.rcParams['axes.unicode_minus'] = False + + +def richards(t, K, r, t0, v): + """Richards curve (generalized logistic model)""" + exp_term = np.exp(-r * (t - t0)) + exp_term = np.clip(exp_term, 1e-10, 1e10) + return K / np.power(1 + exp_term, 1/v) + + +def richards_fixed_K(t, r, t0, v): + """Richards curve with fixed K=4298""" + K = 4298 + exp_term = np.exp(-r * (t - t0)) + exp_term = np.clip(exp_term, 1e-10, 1e10) + return K / np.power(1 + exp_term, 1/v) + + +def load_data(filepath="rocket_launch_counts.csv"): + """Load rocket launch data""" + df = pd.read_csv(filepath) + df = df.rename(columns={"YDate": "year", "Total": "launches"}) + df["year"] = pd.to_numeric(df["year"], errors="coerce") + df["launches"] = pd.to_numeric(df["launches"], errors="coerce") + df = df.dropna(subset=["year", "launches"]) + df = df[(df["year"] >= 1957) & (df["year"] <= 2025)] + df = df.astype({"year": int, "launches": int}) + df = df.sort_values("year").reset_index(drop=True) + return df + + +def main(): + print("=" * 60) + print("Richards S-Curve Fit (2010-2025 Data, K=4298)") + print("=" * 60) + + # Load data + df = load_data() + + # Filter 2010-2025 + start_year = 2010 + end_year = 2025 + data = df[(df["year"] >= start_year) & (df["year"] <= end_year)].copy() + years = data["year"].values + launches = data["launches"].values + + print(f"Data range: {start_year} - {end_year}") + print(f"Data points: {len(data)}") + print(f"\nHistorical data:") + for y, l in zip(years, launches): + print(f" {y}: {l} launches") + + # Fit with fixed K=4298 + K_fixed = 4298 + t = (years - start_year).astype(float) + + p0 = [0.2, 20.0, 2.0] # r, t0, v + bounds = ([0.01, 5, 0.5], [1.0, 100, 10.0]) + + try: + popt, pcov = curve_fit(richards_fixed_K, t, launches, p0=p0, bounds=bounds, maxfev=100000) + r, t0, v = popt + + # Calculate R² + y_pred = richards_fixed_K(t, *popt) + ss_res = np.sum((launches - y_pred) ** 2) + ss_tot = np.sum((launches - np.mean(launches)) ** 2) + r_squared = 1 - (ss_res / ss_tot) + + print(f"\nFitted Parameters (K fixed at {K_fixed}):") + print(f" K (carrying capacity) = {K_fixed} launches/year (FIXED)") + print(f" r (growth rate) = {r:.4f}") + print(f" t0 (inflection point) = {start_year + t0:.1f}") + print(f" v (shape parameter) = {v:.3f}") + print(f" R² = {r_squared:.4f}") + + except Exception as e: + print(f"Fitting error: {e}") + return + + # Physical limit + physical_max = 3650 + print(f"\nPhysical limit: {physical_max} (10 sites × 365 days)") + print(f"K / Physical limit = {K_fixed/physical_max:.2f}x") + + # ========== Create Visualization ========== + # 缩小图片尺寸（比例不变），使字体相对更大 + 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=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, 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=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=color_curve, ls=':', lw=1.5, alpha=0.6, + label=f'K = {K_fixed}') + + # Mark 2050 line only + 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=color_target, s=60, marker='D', zorder=4) + ax.annotate(f'{pred_2050:.0f}', xy=(2050, pred_2050), + xytext=(2050.5, pred_2050+180), + fontsize=9, color=color_target, fontweight='bold') + + # Formatting + 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() + plt.savefig('richards_curve_2010_K4298.png', dpi=150, bbox_inches='tight') + plt.close() + + print("\nPlot saved: richards_curve_2010_K4298.png") + print("=" * 60) + + +if __name__ == "__main__": + main() diff --git a/p1/sensitivity_analysis.py b/p1/sensitivity_analysis.py index 9e72063..d9b5052 100644 --- a/p1/sensitivity_analysis.py +++ b/p1/sensitivity_analysis.py @@ -1,853 +1,853 @@ -""" -Task 1 Sensitivity Analysis for Moon Colony Logistics - -This module performs comprehensive sensitivity analysis on key model parameters: -1. Rocket payload capacity (100-150 tons, as specified in problem) -2. Elevator capacity (±20%) -3. Launch frequency (0.5-2 launches/day) -4. Engine technology (Isp variations) -5. Structural coefficient α - -Author: MCM 2026 Team -""" - -import numpy as np -import matplotlib -matplotlib.use('Agg') -import matplotlib.pyplot as plt -from matplotlib import rcParams -import pandas as pd -from typing import Dict, List, Tuple, Optional -from dataclasses import dataclass -import warnings -warnings.filterwarnings('ignore') - -# Font settings -rcParams['font.sans-serif'] = ['Arial Unicode MS', 'DejaVu Sans', 'SimHei'] -rcParams['axes.unicode_minus'] = False - -# ============== Physical Constants ============== -G0 = 9.81 # m/s² -OMEGA_EARTH = 7.27e-5 # rad/s -R_EARTH = 6.371e6 # m - -# Baseline Mission parameters -TOTAL_PAYLOAD = 100e6 # 100 million metric tons - -# Baseline Space Elevator parameters -BASELINE_NUM_ELEVATORS = 3 -BASELINE_ELEVATOR_CAPACITY_PER_YEAR = 179000 # metric tons per elevator per year -BASELINE_ELEVATOR_SPECIFIC_ENERGY = 157.2e9 # J per metric ton - -# Baseline Rocket parameters -BASELINE_PAYLOAD_PER_LAUNCH = 125 # metric tons (range: 100-150) -BASELINE_ISP = 363 # seconds (LOX/CH4) -BASELINE_SPECIFIC_FUEL_ENERGY = 11.9e6 # J/kg -BASELINE_ALPHA = 0.10 # structural coefficient -BASELINE_NUM_STAGES = 3 -BASELINE_DELTA_V = 13300 # m/s - -# Baseline Launch frequency -BASELINE_LAUNCHES_PER_DAY = 1 # per site - - -# ============== Launch Sites ============== -@dataclass -class LaunchSite: - name: str - short_name: str - latitude: float - - @property - def abs_latitude(self) -> float: - return abs(self.latitude) - - def delta_v_loss(self, omega=OMEGA_EARTH, r=R_EARTH) -> float: - v_equator = omega * r - v_site = omega * r * np.cos(np.radians(self.abs_latitude)) - return v_equator - v_site - - -LAUNCH_SITES = [ - LaunchSite("Kourou", "Kourou", 5.2), - LaunchSite("SDSC", "SDSC", 13.7), - LaunchSite("Texas", "Texas", 26.0), - LaunchSite("Florida", "Florida", 28.5), - LaunchSite("California", "California", 34.7), - LaunchSite("Virginia", "Virginia", 37.8), - LaunchSite("Taiyuan", "Taiyuan", 38.8), - LaunchSite("Mahia", "Mahia", 39.3), - LaunchSite("Baikonur", "Baikonur", 45.6), - LaunchSite("Alaska", "Alaska", 57.4), -] -LAUNCH_SITES = sorted(LAUNCH_SITES, key=lambda x: x.abs_latitude) -NUM_SITES = len(LAUNCH_SITES) - - -# ============== Core Calculation Functions ============== - -def fuel_ratio_multistage(delta_v: float, isp: float, alpha: float, num_stages: int) -> float: - """Calculate multi-stage rocket fuel/payload ratio""" - ve = isp * G0 - 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 calculate_scenario( - completion_years: float, - payload_per_launch: float = BASELINE_PAYLOAD_PER_LAUNCH, - elevator_capacity: float = BASELINE_NUM_ELEVATORS * BASELINE_ELEVATOR_CAPACITY_PER_YEAR, - launches_per_day: float = BASELINE_LAUNCHES_PER_DAY, - isp: float = BASELINE_ISP, - alpha: float = BASELINE_ALPHA, - num_stages: int = BASELINE_NUM_STAGES, - delta_v_base: float = BASELINE_DELTA_V, - specific_fuel_energy: float = BASELINE_SPECIFIC_FUEL_ENERGY, -) -> Optional[Dict]: - """Calculate scenario for given parameters""" - - # Space elevator transport - elevator_payload = min(elevator_capacity * completion_years, TOTAL_PAYLOAD) - elevator_energy = elevator_payload * BASELINE_ELEVATOR_SPECIFIC_ENERGY - - # Remaining payload for rockets - remaining_payload = TOTAL_PAYLOAD - elevator_payload - - if remaining_payload <= 0: - return { - 'years': completion_years, - 'elevator_payload': elevator_payload, - 'rocket_payload': 0, - 'total_energy_PJ': elevator_energy / 1e15, - 'rocket_launches': 0, - 'elevator_fraction': 1.0, - 'feasible': True, - } - - # Rocket launches needed - rocket_launches_needed = int(np.ceil(remaining_payload / payload_per_launch)) - - # Check feasibility - days_available = completion_years * 365 - max_launches_per_site = int(days_available * launches_per_day) - total_rocket_capacity = NUM_SITES * max_launches_per_site * payload_per_launch - - if remaining_payload > total_rocket_capacity: - return { - 'years': completion_years, - 'feasible': False, - 'shortage': remaining_payload - total_rocket_capacity, - } - - # Calculate rocket energy (low-latitude priority) - rocket_energy = 0 - remaining_launches = rocket_launches_needed - - for site in LAUNCH_SITES: - if remaining_launches <= 0: - break - allocated = min(remaining_launches, max_launches_per_site) - - # Calculate fuel ratio for this site - total_dv = delta_v_base + site.delta_v_loss() - k = fuel_ratio_multistage(total_dv, isp, alpha, num_stages) - fuel_per_ton = k * 1000 # kg fuel per metric ton payload - energy_per_ton = fuel_per_ton * specific_fuel_energy - - rocket_energy += energy_per_ton * payload_per_launch * allocated - remaining_launches -= allocated - - total_energy = elevator_energy + rocket_energy - - return { - 'years': completion_years, - 'elevator_payload': elevator_payload, - 'rocket_payload': remaining_payload, - 'total_energy_PJ': total_energy / 1e15, - 'rocket_launches': rocket_launches_needed, - 'elevator_fraction': elevator_payload / TOTAL_PAYLOAD, - 'feasible': True, - } - - -def find_minimum_timeline( - payload_per_launch: float = BASELINE_PAYLOAD_PER_LAUNCH, - elevator_capacity: float = BASELINE_NUM_ELEVATORS * BASELINE_ELEVATOR_CAPACITY_PER_YEAR, - launches_per_day: float = BASELINE_LAUNCHES_PER_DAY, - **kwargs -) -> float: - """Find minimum feasible timeline""" - for years in np.linspace(50, 300, 500): - result = calculate_scenario( - years, payload_per_launch, elevator_capacity, launches_per_day, **kwargs - ) - if result and result.get('feasible', False): - return years - return 300 # fallback - - -def generate_tradeoff_curve( - year_min: float = None, - year_max: float = 250, - num_points: int = 200, - **kwargs -) -> pd.DataFrame: - """Generate trade-off curve for given parameters""" - - if year_min is None: - year_min = find_minimum_timeline(**kwargs) - - years_range = np.linspace(year_min, year_max, num_points) - - results = [] - for years in years_range: - scenario = calculate_scenario(years, **kwargs) - if scenario and scenario.get('feasible', False): - results.append({ - 'years': years, - 'energy_PJ': scenario['total_energy_PJ'], - 'elevator_fraction': scenario['elevator_fraction'], - 'rocket_launches': scenario.get('rocket_launches', 0), - }) - - return pd.DataFrame(results) - - -def find_optimal_point(df: pd.DataFrame, lambda_cost: float) -> Dict: - """Find optimal point for given λ""" - if df.empty: - return {'years': np.nan, 'energy_PJ': np.nan} - - total_cost = df['energy_PJ'].values + lambda_cost * df['years'].values - opt_idx = np.argmin(total_cost) - - return { - 'years': df['years'].iloc[opt_idx], - 'energy_PJ': df['energy_PJ'].iloc[opt_idx], - 'elevator_fraction': df['elevator_fraction'].iloc[opt_idx], - } - - -# ============== Sensitivity Analysis Functions ============== - -def sensitivity_payload_capacity(): - """Analyze sensitivity to rocket payload capacity (100-150 tons)""" - print("\n[1] Analyzing Rocket Payload Capacity Sensitivity (100-150 tons)...") - - payload_range = np.linspace(100, 150, 11) - results = [] - - for payload in payload_range: - elevator_cap = BASELINE_NUM_ELEVATORS * BASELINE_ELEVATOR_CAPACITY_PER_YEAR - t_min = find_minimum_timeline(payload_per_launch=payload, elevator_capacity=elevator_cap) - t_elev = TOTAL_PAYLOAD / elevator_cap - - df = generate_tradeoff_curve( - year_min=t_min, year_max=250, - payload_per_launch=payload, elevator_capacity=elevator_cap - ) - - if not df.empty: - opt_504 = find_optimal_point(df, 504) - energy_at_139 = df[df['years'] >= 139]['energy_PJ'].iloc[0] if any(df['years'] >= 139) else np.nan - - results.append({ - 'payload_tons': payload, - 'min_timeline': t_min, - 'elevator_only_timeline': t_elev, - 'optimal_T_lambda504': opt_504['years'], - 'optimal_E_lambda504': opt_504['energy_PJ'], - 'energy_at_139y': energy_at_139, - }) - - return pd.DataFrame(results) - - -def sensitivity_elevator_capacity(): - """Analyze sensitivity to elevator capacity (±20%)""" - print("\n[2] Analyzing Elevator Capacity Sensitivity (±20%)...") - - baseline_cap = BASELINE_NUM_ELEVATORS * BASELINE_ELEVATOR_CAPACITY_PER_YEAR - capacity_factors = np.linspace(0.8, 1.2, 9) - results = [] - - for factor in capacity_factors: - elevator_cap = baseline_cap * factor - t_min = find_minimum_timeline(elevator_capacity=elevator_cap) - t_elev = TOTAL_PAYLOAD / elevator_cap - - df = generate_tradeoff_curve( - year_min=t_min, year_max=300, - elevator_capacity=elevator_cap - ) - - if not df.empty: - # Find critical λ - lambda_crit = None - for lam in np.linspace(300, 700, 200): - opt = find_optimal_point(df, lam) - if opt['years'] < t_elev - 5: - lambda_crit = lam - break - - opt_504 = find_optimal_point(df, 504) - - results.append({ - 'capacity_factor': factor, - 'capacity_tons_per_year': elevator_cap, - 'min_timeline': t_min, - 'elevator_only_timeline': t_elev, - 'critical_lambda': lambda_crit, - 'optimal_T_lambda504': opt_504['years'], - 'optimal_E_lambda504': opt_504['energy_PJ'], - }) - - return pd.DataFrame(results) - - -def sensitivity_launch_frequency(): - """Analyze sensitivity to launch frequency""" - print("\n[3] Analyzing Launch Frequency Sensitivity (0.5-2 launches/day)...") - - frequency_range = np.linspace(0.5, 2.0, 7) - results = [] - - for freq in frequency_range: - elevator_cap = BASELINE_NUM_ELEVATORS * BASELINE_ELEVATOR_CAPACITY_PER_YEAR - t_min = find_minimum_timeline(launches_per_day=freq, elevator_capacity=elevator_cap) - t_elev = TOTAL_PAYLOAD / elevator_cap - - df = generate_tradeoff_curve( - year_min=t_min, year_max=250, - launches_per_day=freq, elevator_capacity=elevator_cap - ) - - if not df.empty: - opt_504 = find_optimal_point(df, 504) - - results.append({ - 'launches_per_day': freq, - 'min_timeline': t_min, - 'elevator_only_timeline': t_elev, - 'optimal_T_lambda504': opt_504['years'], - 'optimal_E_lambda504': opt_504['energy_PJ'], - 'rocket_capacity_tons_per_year': freq * 365 * NUM_SITES * BASELINE_PAYLOAD_PER_LAUNCH, - }) - - return pd.DataFrame(results) - - -def sensitivity_engine_technology(): - """Analyze sensitivity to engine technology (Isp)""" - print("\n[4] Analyzing Engine Technology Sensitivity (Isp)...") - - engines = [ - ('Solid', 280, 5.0e6), - ('LOX/Kerosene', 350, 10.3e6), - ('LOX/Methane', 363, 11.9e6), - ('LOX/Hydrogen', 450, 15.5e6), - ] - - results = [] - - for name, isp, specific_energy in engines: - elevator_cap = BASELINE_NUM_ELEVATORS * BASELINE_ELEVATOR_CAPACITY_PER_YEAR - t_min = find_minimum_timeline( - isp=isp, specific_fuel_energy=specific_energy, elevator_capacity=elevator_cap - ) - - df = generate_tradeoff_curve( - year_min=t_min, year_max=250, - isp=isp, specific_fuel_energy=specific_energy, elevator_capacity=elevator_cap - ) - - if not df.empty: - opt_504 = find_optimal_point(df, 504) - - # Calculate fuel ratio at equator - k = fuel_ratio_multistage(BASELINE_DELTA_V, isp, BASELINE_ALPHA, BASELINE_NUM_STAGES) - - results.append({ - 'engine_type': name, - 'isp_seconds': isp, - 'specific_energy_MJ_kg': specific_energy / 1e6, - 'fuel_ratio': k, - 'min_timeline': t_min, - 'optimal_T_lambda504': opt_504['years'], - 'optimal_E_lambda504': opt_504['energy_PJ'], - }) - - return pd.DataFrame(results) - - -def sensitivity_structural_coefficient(): - """Analyze sensitivity to structural coefficient α""" - print("\n[5] Analyzing Structural Coefficient Sensitivity (α)...") - - alpha_range = np.linspace(0.06, 0.14, 9) - results = [] - - for alpha in alpha_range: - elevator_cap = BASELINE_NUM_ELEVATORS * BASELINE_ELEVATOR_CAPACITY_PER_YEAR - t_min = find_minimum_timeline(alpha=alpha, elevator_capacity=elevator_cap) - - df = generate_tradeoff_curve( - year_min=t_min, year_max=250, - alpha=alpha, elevator_capacity=elevator_cap - ) - - if not df.empty: - opt_504 = find_optimal_point(df, 504) - k = fuel_ratio_multistage(BASELINE_DELTA_V, BASELINE_ISP, alpha, BASELINE_NUM_STAGES) - - results.append({ - 'alpha': alpha, - 'fuel_ratio': k, - 'min_timeline': t_min, - 'optimal_T_lambda504': opt_504['years'], - 'optimal_E_lambda504': opt_504['energy_PJ'], - }) - - return pd.DataFrame(results) - - -# ============== Visualization ============== - -def plot_sensitivity_analysis( - payload_df: pd.DataFrame, - elevator_df: pd.DataFrame, - frequency_df: pd.DataFrame, - engine_df: pd.DataFrame, - alpha_df: pd.DataFrame, - save_path: str = '/Volumes/Files/code/mm/20260130_b/p1/sensitivity_analysis.png' -): - """Create comprehensive sensitivity analysis visualization""" - - fig, axes = plt.subplots(2, 3, figsize=(18, 12)) - - # ========== Plot 1: Rocket Payload Capacity ========== - ax1 = axes[0, 0] - - ax1.plot(payload_df['payload_tons'], payload_df['min_timeline'], - 'b-o', linewidth=2, markersize=8, label='Minimum Timeline') - ax1.plot(payload_df['payload_tons'], payload_df['optimal_T_lambda504'], - 'r-s', linewidth=2, markersize=8, label='Optimal T* (λ=504)') - - ax1.axhline(y=186.2, color='green', linestyle='--', alpha=0.7, label='Elevator-only') - ax1.axvline(x=125, color='gray', linestyle=':', alpha=0.7, label='Baseline (125t)') - - ax1.fill_betweenx([80, 220], 100, 150, alpha=0.1, color='blue', label='Problem range') - - ax1.set_xlabel('Rocket Payload Capacity (metric tons)', fontsize=11) - ax1.set_ylabel('Timeline (years)', fontsize=11) - ax1.set_title('(a) Sensitivity to Rocket Payload Capacity\n(Problem specifies 100-150 tons)', fontsize=12) - ax1.legend(loc='upper right', fontsize=9) - ax1.grid(True, alpha=0.3) - ax1.set_xlim(95, 155) - ax1.set_ylim(80, 200) - - # ========== Plot 2: Elevator Capacity ========== - ax2 = axes[0, 1] - - ax2.plot(elevator_df['capacity_factor'] * 100, elevator_df['elevator_only_timeline'], - 'g-o', linewidth=2, markersize=8, label='Elevator-only Timeline') - ax2.plot(elevator_df['capacity_factor'] * 100, elevator_df['optimal_T_lambda504'], - 'r-s', linewidth=2, markersize=8, label='Optimal T* (λ=504)') - ax2.plot(elevator_df['capacity_factor'] * 100, elevator_df['min_timeline'], - 'b-^', linewidth=2, markersize=8, label='Minimum Timeline') - - ax2.axvline(x=100, color='gray', linestyle=':', alpha=0.7, label='Baseline') - - ax2.set_xlabel('Elevator Capacity (% of baseline)', fontsize=11) - ax2.set_ylabel('Timeline (years)', fontsize=11) - ax2.set_title('(b) Sensitivity to Elevator Capacity (±20%)', fontsize=12) - ax2.legend(loc='upper right', fontsize=9) - ax2.grid(True, alpha=0.3) - ax2.set_xlim(75, 125) - - # ========== Plot 3: Launch Frequency ========== - ax3 = axes[0, 2] - - ax3.plot(frequency_df['launches_per_day'], frequency_df['min_timeline'], - 'b-o', linewidth=2, markersize=8, label='Minimum Timeline') - ax3.plot(frequency_df['launches_per_day'], frequency_df['optimal_T_lambda504'], - 'r-s', linewidth=2, markersize=8, label='Optimal T* (λ=504)') - - ax3.axhline(y=186.2, color='green', linestyle='--', alpha=0.7, label='Elevator-only') - ax3.axvline(x=1.0, color='gray', linestyle=':', alpha=0.7, label='Baseline (1/day)') - - ax3.set_xlabel('Launch Frequency (launches/site/day)', fontsize=11) - ax3.set_ylabel('Timeline (years)', fontsize=11) - ax3.set_title('(c) Sensitivity to Launch Frequency', fontsize=12) - ax3.legend(loc='upper right', fontsize=9) - ax3.grid(True, alpha=0.3) - - # ========== Plot 4: Engine Technology ========== - ax4 = axes[1, 0] - - x_pos = np.arange(len(engine_df)) - width = 0.35 - - bars1 = ax4.bar(x_pos - width/2, engine_df['optimal_E_lambda504'] / 1000, width, - label='Energy at T* (×10³ PJ)', color='steelblue', alpha=0.8) - - ax4_twin = ax4.twinx() - bars2 = ax4_twin.bar(x_pos + width/2, engine_df['fuel_ratio'], width, - label='Fuel Ratio k', color='coral', alpha=0.8) - - ax4.set_xlabel('Engine Type', fontsize=11) - ax4.set_ylabel('Energy (×10³ PJ)', fontsize=11, color='steelblue') - ax4_twin.set_ylabel('Fuel Ratio k', fontsize=11, color='coral') - ax4.set_title('(d) Sensitivity to Engine Technology', fontsize=12) - ax4.set_xticks(x_pos) - ax4.set_xticklabels(engine_df['engine_type'], rotation=15) - ax4.legend(loc='upper left', fontsize=9) - ax4_twin.legend(loc='upper right', fontsize=9) - ax4.grid(True, alpha=0.3, axis='y') - - # ========== Plot 5: Structural Coefficient ========== - ax5 = axes[1, 1] - - ax5.plot(alpha_df['alpha'], alpha_df['fuel_ratio'], - 'g-o', linewidth=2, markersize=8, label='Fuel Ratio k') - - ax5_twin = ax5.twinx() - ax5_twin.plot(alpha_df['alpha'], alpha_df['optimal_E_lambda504'] / 1000, - 'r-s', linewidth=2, markersize=8, label='Energy at T* (×10³ PJ)') - - ax5.axvline(x=0.10, color='gray', linestyle=':', alpha=0.7, label='Baseline') - - ax5.set_xlabel('Structural Coefficient α', fontsize=11) - ax5.set_ylabel('Fuel Ratio k', fontsize=11, color='green') - ax5_twin.set_ylabel('Energy (×10³ PJ)', fontsize=11, color='red') - ax5.set_title('(e) Sensitivity to Structural Coefficient α', fontsize=12) - ax5.legend(loc='upper left', fontsize=9) - ax5_twin.legend(loc='upper right', fontsize=9) - ax5.grid(True, alpha=0.3) - - # ========== Plot 6: Summary Tornado Chart ========== - ax6 = axes[1, 2] - - # Calculate sensitivity indices (% change in T* for parameter variation) - baseline_T = 139 # baseline optimal timeline - - 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), - ] - - y_pos = np.arange(len(sensitivities)) - - for i, (name, low, high) in enumerate(sensitivities): - ax6.barh(i, high - low, left=min(low, high), height=0.6, - color='steelblue' if high > low else 'coral', alpha=0.7) - ax6.plot([low, high], [i, i], 'k-', linewidth=2) - ax6.plot(low, i, 'ko', markersize=8) - ax6.plot(high, i, 'k^', markersize=8) - - ax6.axvline(x=0, color='black', linestyle='-', linewidth=1) - ax6.set_yticks(y_pos) - ax6.set_yticklabels([s[0] for s in sensitivities]) - ax6.set_xlabel('Change in Optimal Timeline T* (%)', fontsize=11) - ax6.set_title('(f) Sensitivity Summary (Tornado Chart)', fontsize=12) - ax6.grid(True, alpha=0.3, axis='x') - - plt.tight_layout() - plt.savefig(save_path, dpi=150, bbox_inches='tight') - print(f"\nSensitivity analysis plot saved to: {save_path}") - - 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' -): - """Plot trade-off curves for different parameter values""" - - fig, axes = plt.subplots(2, 2, figsize=(14, 12)) - - elevator_cap = BASELINE_NUM_ELEVATORS * BASELINE_ELEVATOR_CAPACITY_PER_YEAR - - # ========== Plot 1: Different Payload Capacities ========== - ax1 = axes[0, 0] - - for payload, color in [(100, 'red'), (125, 'blue'), (150, 'green')]: - df = generate_tradeoff_curve(payload_per_launch=payload, elevator_capacity=elevator_cap) - if not df.empty: - ax1.plot(df['years'], df['energy_PJ'] / 1000, color=color, linewidth=2, - label=f'{payload} tons') - - ax1.set_xlabel('Timeline (years)', fontsize=11) - ax1.set_ylabel('Energy (×10³ PJ)', fontsize=11) - ax1.set_title('Trade-off Curves: Rocket Payload Capacity', fontsize=12) - ax1.legend(title='Payload') - ax1.grid(True, alpha=0.3) - ax1.set_xlim(95, 200) - - # ========== Plot 2: Different Elevator Capacities ========== - ax2 = axes[0, 1] - - for factor, color in [(0.8, 'red'), (1.0, 'blue'), (1.2, 'green')]: - df = generate_tradeoff_curve(elevator_capacity=elevator_cap * factor) - if not df.empty: - ax2.plot(df['years'], df['energy_PJ'] / 1000, color=color, linewidth=2, - label=f'{factor*100:.0f}%') - - ax2.set_xlabel('Timeline (years)', fontsize=11) - ax2.set_ylabel('Energy (×10³ PJ)', fontsize=11) - ax2.set_title('Trade-off Curves: Elevator Capacity', fontsize=12) - ax2.legend(title='Capacity') - ax2.grid(True, alpha=0.3) - ax2.set_xlim(95, 250) - - # ========== Plot 3: Different Launch Frequencies ========== - ax3 = axes[1, 0] - - for freq, color in [(0.5, 'red'), (1.0, 'blue'), (2.0, 'green')]: - df = generate_tradeoff_curve(launches_per_day=freq, elevator_capacity=elevator_cap) - if not df.empty: - ax3.plot(df['years'], df['energy_PJ'] / 1000, color=color, linewidth=2, - label=f'{freq}/day') - - ax3.set_xlabel('Timeline (years)', fontsize=11) - ax3.set_ylabel('Energy (×10³ PJ)', fontsize=11) - ax3.set_title('Trade-off Curves: Launch Frequency', fontsize=12) - ax3.legend(title='Frequency') - ax3.grid(True, alpha=0.3) - ax3.set_xlim(45, 200) - - # ========== Plot 4: Different Engine Types ========== - ax4 = axes[1, 1] - - engines = [ - ('LOX/Kerosene', 350, 10.3e6, 'red'), - ('LOX/Methane', 363, 11.9e6, 'blue'), - ('LOX/Hydrogen', 450, 15.5e6, 'green'), - ] - - for name, isp, se, color in engines: - df = generate_tradeoff_curve(isp=isp, specific_fuel_energy=se, elevator_capacity=elevator_cap) - if not df.empty: - ax4.plot(df['years'], df['energy_PJ'] / 1000, color=color, linewidth=2, - label=f'{name}') - - ax4.set_xlabel('Timeline (years)', fontsize=11) - ax4.set_ylabel('Energy (×10³ PJ)', fontsize=11) - ax4.set_title('Trade-off Curves: Engine Technology', fontsize=12) - ax4.legend(title='Engine') - ax4.grid(True, alpha=0.3) - ax4.set_xlim(95, 200) - - plt.tight_layout() - plt.savefig(save_path, dpi=150, bbox_inches='tight') - print(f"Trade-off comparison plot saved to: {save_path}") - - return fig - - -# ============== Main ============== - -def main(): - print("=" * 70) - print("Task 1 Sensitivity Analysis for Moon Colony Logistics") - print("=" * 70) - - # Run all sensitivity analyses - payload_df = sensitivity_payload_capacity() - elevator_df = sensitivity_elevator_capacity() - frequency_df = sensitivity_launch_frequency() - engine_df = sensitivity_engine_technology() - alpha_df = sensitivity_structural_coefficient() - - # Print summary tables - print("\n" + "=" * 70) - print("SENSITIVITY ANALYSIS RESULTS") - print("=" * 70) - - print("\n[1] Rocket Payload Capacity (Problem range: 100-150 tons)") - print(payload_df.to_string(index=False)) - - print("\n[2] Elevator Capacity (±20%)") - print(elevator_df.to_string(index=False)) - - print("\n[3] Launch Frequency") - print(frequency_df.to_string(index=False)) - - print("\n[4] Engine Technology") - print(engine_df.to_string(index=False)) - - print("\n[5] Structural Coefficient α") - print(alpha_df.to_string(index=False)) - - # Generate plots - print("\n" + "=" * 70) - print("Generating Visualization...") - print("=" * 70) - - plot_sensitivity_analysis(payload_df, elevator_df, frequency_df, engine_df, alpha_df) - plot_tradeoff_comparison() - - # Save data - payload_df.to_csv('/Volumes/Files/code/mm/20260130_b/p1/sensitivity_payload.csv', index=False) - elevator_df.to_csv('/Volumes/Files/code/mm/20260130_b/p1/sensitivity_elevator.csv', index=False) - frequency_df.to_csv('/Volumes/Files/code/mm/20260130_b/p1/sensitivity_frequency.csv', index=False) - engine_df.to_csv('/Volumes/Files/code/mm/20260130_b/p1/sensitivity_engine.csv', index=False) - alpha_df.to_csv('/Volumes/Files/code/mm/20260130_b/p1/sensitivity_alpha.csv', index=False) - print("\nData saved to CSV files.") - - # Summary - print("\n" + "=" * 70) - print("KEY FINDINGS") - print("=" * 70) - print(""" - 1. ROCKET PAYLOAD CAPACITY (100-150 tons): - - Higher capacity reduces minimum timeline significantly - - Limited impact on energy when elevator dominates - - Baseline 125 tons is a reasonable middle estimate - - 2. ELEVATOR CAPACITY (±20%): - - Most sensitive parameter for timeline - - ±20% capacity changes elevator-only timeline by ~37 years - - Critical for long-term planning - - 3. LAUNCH FREQUENCY: - - Doubling frequency halves minimum rocket-only timeline - - Essential for time-constrained scenarios - - Less impact when elevator provides majority of transport - - 4. ENGINE TECHNOLOGY: - - Primarily affects energy, not timeline - - LOX/Hydrogen saves ~30% energy vs LOX/Methane - - But LOX/Methane has practical advantages (storage, cost) - - 5. STRUCTURAL COEFFICIENT α: - - Lower α (better technology) improves fuel efficiency - - Moderate impact on total energy - - Progress from 0.10 to 0.08 saves ~8% energy - """) - - print("=" * 70) - print("Analysis Complete!") - print("=" * 70) - - return payload_df, elevator_df, frequency_df, engine_df, alpha_df - - -if __name__ == "__main__": - results = main() +""" +Task 1 Sensitivity Analysis for Moon Colony Logistics + +This module performs comprehensive sensitivity analysis on key model parameters: +1. Rocket payload capacity (100-150 tons, as specified in problem) +2. Elevator capacity (±20%) +3. Launch frequency (0.5-2 launches/day) +4. Engine technology (Isp variations) +5. Structural coefficient α + +Author: MCM 2026 Team +""" + +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from matplotlib import rcParams +import pandas as pd +from typing import Dict, List, Tuple, Optional +from dataclasses import dataclass +import warnings +warnings.filterwarnings('ignore') + +# Font settings +rcParams['font.sans-serif'] = ['Arial Unicode MS', 'DejaVu Sans', 'SimHei'] +rcParams['axes.unicode_minus'] = False + +# ============== Physical Constants ============== +G0 = 9.81 # m/s² +OMEGA_EARTH = 7.27e-5 # rad/s +R_EARTH = 6.371e6 # m + +# Baseline Mission parameters +TOTAL_PAYLOAD = 100e6 # 100 million metric tons + +# Baseline Space Elevator parameters +BASELINE_NUM_ELEVATORS = 3 +BASELINE_ELEVATOR_CAPACITY_PER_YEAR = 179000 # metric tons per elevator per year +BASELINE_ELEVATOR_SPECIFIC_ENERGY = 157.2e9 # J per metric ton + +# Baseline Rocket parameters +BASELINE_PAYLOAD_PER_LAUNCH = 125 # metric tons (range: 100-150) +BASELINE_ISP = 363 # seconds (LOX/CH4) +BASELINE_SPECIFIC_FUEL_ENERGY = 11.9e6 # J/kg +BASELINE_ALPHA = 0.10 # structural coefficient +BASELINE_NUM_STAGES = 3 +BASELINE_DELTA_V = 13300 # m/s + +# Baseline Launch frequency +BASELINE_LAUNCHES_PER_DAY = 1 # per site + + +# ============== Launch Sites ============== +@dataclass +class LaunchSite: + name: str + short_name: str + latitude: float + + @property + def abs_latitude(self) -> float: + return abs(self.latitude) + + def delta_v_loss(self, omega=OMEGA_EARTH, r=R_EARTH) -> float: + v_equator = omega * r + v_site = omega * r * np.cos(np.radians(self.abs_latitude)) + return v_equator - v_site + + +LAUNCH_SITES = [ + LaunchSite("Kourou", "Kourou", 5.2), + LaunchSite("SDSC", "SDSC", 13.7), + LaunchSite("Texas", "Texas", 26.0), + LaunchSite("Florida", "Florida", 28.5), + LaunchSite("California", "California", 34.7), + LaunchSite("Virginia", "Virginia", 37.8), + LaunchSite("Taiyuan", "Taiyuan", 38.8), + LaunchSite("Mahia", "Mahia", 39.3), + LaunchSite("Baikonur", "Baikonur", 45.6), + LaunchSite("Alaska", "Alaska", 57.4), +] +LAUNCH_SITES = sorted(LAUNCH_SITES, key=lambda x: x.abs_latitude) +NUM_SITES = len(LAUNCH_SITES) + + +# ============== Core Calculation Functions ============== + +def fuel_ratio_multistage(delta_v: float, isp: float, alpha: float, num_stages: int) -> float: + """Calculate multi-stage rocket fuel/payload ratio""" + ve = isp * G0 + 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 calculate_scenario( + completion_years: float, + payload_per_launch: float = BASELINE_PAYLOAD_PER_LAUNCH, + elevator_capacity: float = BASELINE_NUM_ELEVATORS * BASELINE_ELEVATOR_CAPACITY_PER_YEAR, + launches_per_day: float = BASELINE_LAUNCHES_PER_DAY, + isp: float = BASELINE_ISP, + alpha: float = BASELINE_ALPHA, + num_stages: int = BASELINE_NUM_STAGES, + delta_v_base: float = BASELINE_DELTA_V, + specific_fuel_energy: float = BASELINE_SPECIFIC_FUEL_ENERGY, +) -> Optional[Dict]: + """Calculate scenario for given parameters""" + + # Space elevator transport + elevator_payload = min(elevator_capacity * completion_years, TOTAL_PAYLOAD) + elevator_energy = elevator_payload * BASELINE_ELEVATOR_SPECIFIC_ENERGY + + # Remaining payload for rockets + remaining_payload = TOTAL_PAYLOAD - elevator_payload + + if remaining_payload <= 0: + return { + 'years': completion_years, + 'elevator_payload': elevator_payload, + 'rocket_payload': 0, + 'total_energy_PJ': elevator_energy / 1e15, + 'rocket_launches': 0, + 'elevator_fraction': 1.0, + 'feasible': True, + } + + # Rocket launches needed + rocket_launches_needed = int(np.ceil(remaining_payload / payload_per_launch)) + + # Check feasibility + days_available = completion_years * 365 + max_launches_per_site = int(days_available * launches_per_day) + total_rocket_capacity = NUM_SITES * max_launches_per_site * payload_per_launch + + if remaining_payload > total_rocket_capacity: + return { + 'years': completion_years, + 'feasible': False, + 'shortage': remaining_payload - total_rocket_capacity, + } + + # Calculate rocket energy (low-latitude priority) + rocket_energy = 0 + remaining_launches = rocket_launches_needed + + for site in LAUNCH_SITES: + if remaining_launches <= 0: + break + allocated = min(remaining_launches, max_launches_per_site) + + # Calculate fuel ratio for this site + total_dv = delta_v_base + site.delta_v_loss() + k = fuel_ratio_multistage(total_dv, isp, alpha, num_stages) + fuel_per_ton = k * 1000 # kg fuel per metric ton payload + energy_per_ton = fuel_per_ton * specific_fuel_energy + + rocket_energy += energy_per_ton * payload_per_launch * allocated + remaining_launches -= allocated + + total_energy = elevator_energy + rocket_energy + + return { + 'years': completion_years, + 'elevator_payload': elevator_payload, + 'rocket_payload': remaining_payload, + 'total_energy_PJ': total_energy / 1e15, + 'rocket_launches': rocket_launches_needed, + 'elevator_fraction': elevator_payload / TOTAL_PAYLOAD, + 'feasible': True, + } + + +def find_minimum_timeline( + payload_per_launch: float = BASELINE_PAYLOAD_PER_LAUNCH, + elevator_capacity: float = BASELINE_NUM_ELEVATORS * BASELINE_ELEVATOR_CAPACITY_PER_YEAR, + launches_per_day: float = BASELINE_LAUNCHES_PER_DAY, + **kwargs +) -> float: + """Find minimum feasible timeline""" + for years in np.linspace(50, 300, 500): + result = calculate_scenario( + years, payload_per_launch, elevator_capacity, launches_per_day, **kwargs + ) + if result and result.get('feasible', False): + return years + return 300 # fallback + + +def generate_tradeoff_curve( + year_min: float = None, + year_max: float = 250, + num_points: int = 200, + **kwargs +) -> pd.DataFrame: + """Generate trade-off curve for given parameters""" + + if year_min is None: + year_min = find_minimum_timeline(**kwargs) + + years_range = np.linspace(year_min, year_max, num_points) + + results = [] + for years in years_range: + scenario = calculate_scenario(years, **kwargs) + if scenario and scenario.get('feasible', False): + results.append({ + 'years': years, + 'energy_PJ': scenario['total_energy_PJ'], + 'elevator_fraction': scenario['elevator_fraction'], + 'rocket_launches': scenario.get('rocket_launches', 0), + }) + + return pd.DataFrame(results) + + +def find_optimal_point(df: pd.DataFrame, lambda_cost: float) -> Dict: + """Find optimal point for given λ""" + if df.empty: + return {'years': np.nan, 'energy_PJ': np.nan} + + total_cost = df['energy_PJ'].values + lambda_cost * df['years'].values + opt_idx = np.argmin(total_cost) + + return { + 'years': df['years'].iloc[opt_idx], + 'energy_PJ': df['energy_PJ'].iloc[opt_idx], + 'elevator_fraction': df['elevator_fraction'].iloc[opt_idx], + } + + +# ============== Sensitivity Analysis Functions ============== + +def sensitivity_payload_capacity(): + """Analyze sensitivity to rocket payload capacity (100-150 tons)""" + print("\n[1] Analyzing Rocket Payload Capacity Sensitivity (100-150 tons)...") + + payload_range = np.linspace(100, 150, 11) + results = [] + + for payload in payload_range: + elevator_cap = BASELINE_NUM_ELEVATORS * BASELINE_ELEVATOR_CAPACITY_PER_YEAR + t_min = find_minimum_timeline(payload_per_launch=payload, elevator_capacity=elevator_cap) + t_elev = TOTAL_PAYLOAD / elevator_cap + + df = generate_tradeoff_curve( + year_min=t_min, year_max=250, + payload_per_launch=payload, elevator_capacity=elevator_cap + ) + + if not df.empty: + opt_504 = find_optimal_point(df, 504) + energy_at_139 = df[df['years'] >= 139]['energy_PJ'].iloc[0] if any(df['years'] >= 139) else np.nan + + results.append({ + 'payload_tons': payload, + 'min_timeline': t_min, + 'elevator_only_timeline': t_elev, + 'optimal_T_lambda504': opt_504['years'], + 'optimal_E_lambda504': opt_504['energy_PJ'], + 'energy_at_139y': energy_at_139, + }) + + return pd.DataFrame(results) + + +def sensitivity_elevator_capacity(): + """Analyze sensitivity to elevator capacity (±20%)""" + print("\n[2] Analyzing Elevator Capacity Sensitivity (±20%)...") + + baseline_cap = BASELINE_NUM_ELEVATORS * BASELINE_ELEVATOR_CAPACITY_PER_YEAR + capacity_factors = np.linspace(0.8, 1.2, 9) + results = [] + + for factor in capacity_factors: + elevator_cap = baseline_cap * factor + t_min = find_minimum_timeline(elevator_capacity=elevator_cap) + t_elev = TOTAL_PAYLOAD / elevator_cap + + df = generate_tradeoff_curve( + year_min=t_min, year_max=300, + elevator_capacity=elevator_cap + ) + + if not df.empty: + # Find critical λ + lambda_crit = None + for lam in np.linspace(300, 700, 200): + opt = find_optimal_point(df, lam) + if opt['years'] < t_elev - 5: + lambda_crit = lam + break + + opt_504 = find_optimal_point(df, 504) + + results.append({ + 'capacity_factor': factor, + 'capacity_tons_per_year': elevator_cap, + 'min_timeline': t_min, + 'elevator_only_timeline': t_elev, + 'critical_lambda': lambda_crit, + 'optimal_T_lambda504': opt_504['years'], + 'optimal_E_lambda504': opt_504['energy_PJ'], + }) + + return pd.DataFrame(results) + + +def sensitivity_launch_frequency(): + """Analyze sensitivity to launch frequency""" + print("\n[3] Analyzing Launch Frequency Sensitivity (0.5-2 launches/day)...") + + frequency_range = np.linspace(0.5, 2.0, 7) + results = [] + + for freq in frequency_range: + elevator_cap = BASELINE_NUM_ELEVATORS * BASELINE_ELEVATOR_CAPACITY_PER_YEAR + t_min = find_minimum_timeline(launches_per_day=freq, elevator_capacity=elevator_cap) + t_elev = TOTAL_PAYLOAD / elevator_cap + + df = generate_tradeoff_curve( + year_min=t_min, year_max=250, + launches_per_day=freq, elevator_capacity=elevator_cap + ) + + if not df.empty: + opt_504 = find_optimal_point(df, 504) + + results.append({ + 'launches_per_day': freq, + 'min_timeline': t_min, + 'elevator_only_timeline': t_elev, + 'optimal_T_lambda504': opt_504['years'], + 'optimal_E_lambda504': opt_504['energy_PJ'], + 'rocket_capacity_tons_per_year': freq * 365 * NUM_SITES * BASELINE_PAYLOAD_PER_LAUNCH, + }) + + return pd.DataFrame(results) + + +def sensitivity_engine_technology(): + """Analyze sensitivity to engine technology (Isp)""" + print("\n[4] Analyzing Engine Technology Sensitivity (Isp)...") + + engines = [ + ('Solid', 280, 5.0e6), + ('LOX/Kerosene', 350, 10.3e6), + ('LOX/Methane', 363, 11.9e6), + ('LOX/Hydrogen', 450, 15.5e6), + ] + + results = [] + + for name, isp, specific_energy in engines: + elevator_cap = BASELINE_NUM_ELEVATORS * BASELINE_ELEVATOR_CAPACITY_PER_YEAR + t_min = find_minimum_timeline( + isp=isp, specific_fuel_energy=specific_energy, elevator_capacity=elevator_cap + ) + + df = generate_tradeoff_curve( + year_min=t_min, year_max=250, + isp=isp, specific_fuel_energy=specific_energy, elevator_capacity=elevator_cap + ) + + if not df.empty: + opt_504 = find_optimal_point(df, 504) + + # Calculate fuel ratio at equator + k = fuel_ratio_multistage(BASELINE_DELTA_V, isp, BASELINE_ALPHA, BASELINE_NUM_STAGES) + + results.append({ + 'engine_type': name, + 'isp_seconds': isp, + 'specific_energy_MJ_kg': specific_energy / 1e6, + 'fuel_ratio': k, + 'min_timeline': t_min, + 'optimal_T_lambda504': opt_504['years'], + 'optimal_E_lambda504': opt_504['energy_PJ'], + }) + + return pd.DataFrame(results) + + +def sensitivity_structural_coefficient(): + """Analyze sensitivity to structural coefficient α""" + print("\n[5] Analyzing Structural Coefficient Sensitivity (α)...") + + alpha_range = np.linspace(0.06, 0.14, 9) + results = [] + + for alpha in alpha_range: + elevator_cap = BASELINE_NUM_ELEVATORS * BASELINE_ELEVATOR_CAPACITY_PER_YEAR + t_min = find_minimum_timeline(alpha=alpha, elevator_capacity=elevator_cap) + + df = generate_tradeoff_curve( + year_min=t_min, year_max=250, + alpha=alpha, elevator_capacity=elevator_cap + ) + + if not df.empty: + opt_504 = find_optimal_point(df, 504) + k = fuel_ratio_multistage(BASELINE_DELTA_V, BASELINE_ISP, alpha, BASELINE_NUM_STAGES) + + results.append({ + 'alpha': alpha, + 'fuel_ratio': k, + 'min_timeline': t_min, + 'optimal_T_lambda504': opt_504['years'], + 'optimal_E_lambda504': opt_504['energy_PJ'], + }) + + return pd.DataFrame(results) + + +# ============== Visualization ============== + +def plot_sensitivity_analysis( + payload_df: pd.DataFrame, + elevator_df: pd.DataFrame, + frequency_df: pd.DataFrame, + engine_df: pd.DataFrame, + alpha_df: pd.DataFrame, + save_path: str = '/Volumes/Files/code/mm/20260130_b/p1/sensitivity_analysis.png' +): + """Create comprehensive sensitivity analysis visualization""" + + fig, axes = plt.subplots(2, 3, figsize=(18, 12)) + + # ========== Plot 1: Rocket Payload Capacity ========== + ax1 = axes[0, 0] + + ax1.plot(payload_df['payload_tons'], payload_df['min_timeline'], + 'b-o', linewidth=2, markersize=8, label='Minimum Timeline') + ax1.plot(payload_df['payload_tons'], payload_df['optimal_T_lambda504'], + 'r-s', linewidth=2, markersize=8, label='Optimal T* (λ=504)') + + ax1.axhline(y=186.2, color='green', linestyle='--', alpha=0.7, label='Elevator-only') + ax1.axvline(x=125, color='gray', linestyle=':', alpha=0.7, label='Baseline (125t)') + + ax1.fill_betweenx([80, 220], 100, 150, alpha=0.1, color='blue', label='Problem range') + + ax1.set_xlabel('Rocket Payload Capacity (metric tons)', fontsize=11) + ax1.set_ylabel('Timeline (years)', fontsize=11) + ax1.set_title('(a) Sensitivity to Rocket Payload Capacity\n(Problem specifies 100-150 tons)', fontsize=12) + ax1.legend(loc='upper right', fontsize=9) + ax1.grid(True, alpha=0.3) + ax1.set_xlim(95, 155) + ax1.set_ylim(80, 200) + + # ========== Plot 2: Elevator Capacity ========== + ax2 = axes[0, 1] + + ax2.plot(elevator_df['capacity_factor'] * 100, elevator_df['elevator_only_timeline'], + 'g-o', linewidth=2, markersize=8, label='Elevator-only Timeline') + ax2.plot(elevator_df['capacity_factor'] * 100, elevator_df['optimal_T_lambda504'], + 'r-s', linewidth=2, markersize=8, label='Optimal T* (λ=504)') + ax2.plot(elevator_df['capacity_factor'] * 100, elevator_df['min_timeline'], + 'b-^', linewidth=2, markersize=8, label='Minimum Timeline') + + ax2.axvline(x=100, color='gray', linestyle=':', alpha=0.7, label='Baseline') + + ax2.set_xlabel('Elevator Capacity (% of baseline)', fontsize=11) + ax2.set_ylabel('Timeline (years)', fontsize=11) + ax2.set_title('(b) Sensitivity to Elevator Capacity (±20%)', fontsize=12) + ax2.legend(loc='upper right', fontsize=9) + ax2.grid(True, alpha=0.3) + ax2.set_xlim(75, 125) + + # ========== Plot 3: Launch Frequency ========== + ax3 = axes[0, 2] + + ax3.plot(frequency_df['launches_per_day'], frequency_df['min_timeline'], + 'b-o', linewidth=2, markersize=8, label='Minimum Timeline') + ax3.plot(frequency_df['launches_per_day'], frequency_df['optimal_T_lambda504'], + 'r-s', linewidth=2, markersize=8, label='Optimal T* (λ=504)') + + ax3.axhline(y=186.2, color='green', linestyle='--', alpha=0.7, label='Elevator-only') + ax3.axvline(x=1.0, color='gray', linestyle=':', alpha=0.7, label='Baseline (1/day)') + + ax3.set_xlabel('Launch Frequency (launches/site/day)', fontsize=11) + ax3.set_ylabel('Timeline (years)', fontsize=11) + ax3.set_title('(c) Sensitivity to Launch Frequency', fontsize=12) + ax3.legend(loc='upper right', fontsize=9) + ax3.grid(True, alpha=0.3) + + # ========== Plot 4: Engine Technology ========== + ax4 = axes[1, 0] + + x_pos = np.arange(len(engine_df)) + width = 0.35 + + bars1 = ax4.bar(x_pos - width/2, engine_df['optimal_E_lambda504'] / 1000, width, + label='Energy at T* (×10³ PJ)', color='steelblue', alpha=0.8) + + ax4_twin = ax4.twinx() + bars2 = ax4_twin.bar(x_pos + width/2, engine_df['fuel_ratio'], width, + label='Fuel Ratio k', color='coral', alpha=0.8) + + ax4.set_xlabel('Engine Type', fontsize=11) + ax4.set_ylabel('Energy (×10³ PJ)', fontsize=11, color='steelblue') + ax4_twin.set_ylabel('Fuel Ratio k', fontsize=11, color='coral') + ax4.set_title('(d) Sensitivity to Engine Technology', fontsize=12) + ax4.set_xticks(x_pos) + ax4.set_xticklabels(engine_df['engine_type'], rotation=15) + ax4.legend(loc='upper left', fontsize=9) + ax4_twin.legend(loc='upper right', fontsize=9) + ax4.grid(True, alpha=0.3, axis='y') + + # ========== Plot 5: Structural Coefficient ========== + ax5 = axes[1, 1] + + ax5.plot(alpha_df['alpha'], alpha_df['fuel_ratio'], + 'g-o', linewidth=2, markersize=8, label='Fuel Ratio k') + + ax5_twin = ax5.twinx() + ax5_twin.plot(alpha_df['alpha'], alpha_df['optimal_E_lambda504'] / 1000, + 'r-s', linewidth=2, markersize=8, label='Energy at T* (×10³ PJ)') + + ax5.axvline(x=0.10, color='gray', linestyle=':', alpha=0.7, label='Baseline') + + ax5.set_xlabel('Structural Coefficient α', fontsize=11) + ax5.set_ylabel('Fuel Ratio k', fontsize=11, color='green') + ax5_twin.set_ylabel('Energy (×10³ PJ)', fontsize=11, color='red') + ax5.set_title('(e) Sensitivity to Structural Coefficient α', fontsize=12) + ax5.legend(loc='upper left', fontsize=9) + ax5_twin.legend(loc='upper right', fontsize=9) + ax5.grid(True, alpha=0.3) + + # ========== Plot 6: Summary Tornado Chart ========== + ax6 = axes[1, 2] + + # Calculate sensitivity indices (% change in T* for parameter variation) + baseline_T = 139 # baseline optimal timeline + + 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), + ] + + y_pos = np.arange(len(sensitivities)) + + for i, (name, low, high) in enumerate(sensitivities): + ax6.barh(i, high - low, left=min(low, high), height=0.6, + color='steelblue' if high > low else 'coral', alpha=0.7) + ax6.plot([low, high], [i, i], 'k-', linewidth=2) + ax6.plot(low, i, 'ko', markersize=8) + ax6.plot(high, i, 'k^', markersize=8) + + ax6.axvline(x=0, color='black', linestyle='-', linewidth=1) + ax6.set_yticks(y_pos) + ax6.set_yticklabels([s[0] for s in sensitivities]) + ax6.set_xlabel('Change in Optimal Timeline T* (%)', fontsize=11) + ax6.set_title('(f) Sensitivity Summary (Tornado Chart)', fontsize=12) + ax6.grid(True, alpha=0.3, axis='x') + + plt.tight_layout() + plt.savefig(save_path, dpi=150, bbox_inches='tight') + print(f"\nSensitivity analysis plot saved to: {save_path}") + + 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' +): + """Plot trade-off curves for different parameter values""" + + fig, axes = plt.subplots(2, 2, figsize=(14, 12)) + + elevator_cap = BASELINE_NUM_ELEVATORS * BASELINE_ELEVATOR_CAPACITY_PER_YEAR + + # ========== Plot 1: Different Payload Capacities ========== + ax1 = axes[0, 0] + + for payload, color in [(100, 'red'), (125, 'blue'), (150, 'green')]: + df = generate_tradeoff_curve(payload_per_launch=payload, elevator_capacity=elevator_cap) + if not df.empty: + ax1.plot(df['years'], df['energy_PJ'] / 1000, color=color, linewidth=2, + label=f'{payload} tons') + + ax1.set_xlabel('Timeline (years)', fontsize=11) + ax1.set_ylabel('Energy (×10³ PJ)', fontsize=11) + ax1.set_title('Trade-off Curves: Rocket Payload Capacity', fontsize=12) + ax1.legend(title='Payload') + ax1.grid(True, alpha=0.3) + ax1.set_xlim(95, 200) + + # ========== Plot 2: Different Elevator Capacities ========== + ax2 = axes[0, 1] + + for factor, color in [(0.8, 'red'), (1.0, 'blue'), (1.2, 'green')]: + df = generate_tradeoff_curve(elevator_capacity=elevator_cap * factor) + if not df.empty: + ax2.plot(df['years'], df['energy_PJ'] / 1000, color=color, linewidth=2, + label=f'{factor*100:.0f}%') + + ax2.set_xlabel('Timeline (years)', fontsize=11) + ax2.set_ylabel('Energy (×10³ PJ)', fontsize=11) + ax2.set_title('Trade-off Curves: Elevator Capacity', fontsize=12) + ax2.legend(title='Capacity') + ax2.grid(True, alpha=0.3) + ax2.set_xlim(95, 250) + + # ========== Plot 3: Different Launch Frequencies ========== + ax3 = axes[1, 0] + + for freq, color in [(0.5, 'red'), (1.0, 'blue'), (2.0, 'green')]: + df = generate_tradeoff_curve(launches_per_day=freq, elevator_capacity=elevator_cap) + if not df.empty: + ax3.plot(df['years'], df['energy_PJ'] / 1000, color=color, linewidth=2, + label=f'{freq}/day') + + ax3.set_xlabel('Timeline (years)', fontsize=11) + ax3.set_ylabel('Energy (×10³ PJ)', fontsize=11) + ax3.set_title('Trade-off Curves: Launch Frequency', fontsize=12) + ax3.legend(title='Frequency') + ax3.grid(True, alpha=0.3) + ax3.set_xlim(45, 200) + + # ========== Plot 4: Different Engine Types ========== + ax4 = axes[1, 1] + + engines = [ + ('LOX/Kerosene', 350, 10.3e6, 'red'), + ('LOX/Methane', 363, 11.9e6, 'blue'), + ('LOX/Hydrogen', 450, 15.5e6, 'green'), + ] + + for name, isp, se, color in engines: + df = generate_tradeoff_curve(isp=isp, specific_fuel_energy=se, elevator_capacity=elevator_cap) + if not df.empty: + ax4.plot(df['years'], df['energy_PJ'] / 1000, color=color, linewidth=2, + label=f'{name}') + + ax4.set_xlabel('Timeline (years)', fontsize=11) + ax4.set_ylabel('Energy (×10³ PJ)', fontsize=11) + ax4.set_title('Trade-off Curves: Engine Technology', fontsize=12) + ax4.legend(title='Engine') + ax4.grid(True, alpha=0.3) + ax4.set_xlim(95, 200) + + plt.tight_layout() + plt.savefig(save_path, dpi=150, bbox_inches='tight') + print(f"Trade-off comparison plot saved to: {save_path}") + + return fig + + +# ============== Main ============== + +def main(): + print("=" * 70) + print("Task 1 Sensitivity Analysis for Moon Colony Logistics") + print("=" * 70) + + # Run all sensitivity analyses + payload_df = sensitivity_payload_capacity() + elevator_df = sensitivity_elevator_capacity() + frequency_df = sensitivity_launch_frequency() + engine_df = sensitivity_engine_technology() + alpha_df = sensitivity_structural_coefficient() + + # Print summary tables + print("\n" + "=" * 70) + print("SENSITIVITY ANALYSIS RESULTS") + print("=" * 70) + + print("\n[1] Rocket Payload Capacity (Problem range: 100-150 tons)") + print(payload_df.to_string(index=False)) + + print("\n[2] Elevator Capacity (±20%)") + print(elevator_df.to_string(index=False)) + + print("\n[3] Launch Frequency") + print(frequency_df.to_string(index=False)) + + print("\n[4] Engine Technology") + print(engine_df.to_string(index=False)) + + print("\n[5] Structural Coefficient α") + print(alpha_df.to_string(index=False)) + + # Generate plots + print("\n" + "=" * 70) + print("Generating Visualization...") + print("=" * 70) + + plot_sensitivity_analysis(payload_df, elevator_df, frequency_df, engine_df, alpha_df) + plot_tradeoff_comparison() + + # Save data + payload_df.to_csv('/Volumes/Files/code/mm/20260130_b/p1/sensitivity_payload.csv', index=False) + elevator_df.to_csv('/Volumes/Files/code/mm/20260130_b/p1/sensitivity_elevator.csv', index=False) + frequency_df.to_csv('/Volumes/Files/code/mm/20260130_b/p1/sensitivity_frequency.csv', index=False) + engine_df.to_csv('/Volumes/Files/code/mm/20260130_b/p1/sensitivity_engine.csv', index=False) + alpha_df.to_csv('/Volumes/Files/code/mm/20260130_b/p1/sensitivity_alpha.csv', index=False) + print("\nData saved to CSV files.") + + # Summary + print("\n" + "=" * 70) + print("KEY FINDINGS") + print("=" * 70) + print(""" + 1. ROCKET PAYLOAD CAPACITY (100-150 tons): + - Higher capacity reduces minimum timeline significantly + - Limited impact on energy when elevator dominates + - Baseline 125 tons is a reasonable middle estimate + + 2. ELEVATOR CAPACITY (±20%): + - Most sensitive parameter for timeline + - ±20% capacity changes elevator-only timeline by ~37 years + - Critical for long-term planning + + 3. LAUNCH FREQUENCY: + - Doubling frequency halves minimum rocket-only timeline + - Essential for time-constrained scenarios + - Less impact when elevator provides majority of transport + + 4. ENGINE TECHNOLOGY: + - Primarily affects energy, not timeline + - LOX/Hydrogen saves ~30% energy vs LOX/Methane + - But LOX/Methane has practical advantages (storage, cost) + + 5. STRUCTURAL COEFFICIENT α: + - Lower α (better technology) improves fuel efficiency + - Moderate impact on total energy + - Progress from 0.10 to 0.08 saves ~8% energy + """) + + print("=" * 70) + print("Analysis Complete!") + print("=" * 70) + + return payload_df, elevator_df, frequency_df, engine_df, alpha_df + + +if __name__ == "__main__": + results = main() diff --git a/p2/completion_time_distribution.png b/p2/completion_time_distribution.png index 2b9751f..a95653e 100644 Binary files a/p2/completion_time_distribution.png and b/p2/completion_time_distribution.png differ diff --git a/p2/completion_time_distribution_paper.pdf b/p2/completion_time_distribution_paper.pdf new file mode 100644 index 0000000..23c79c0 Binary files /dev/null and b/p2/completion_time_distribution_paper.pdf differ diff --git a/p2/completion_time_distribution_paper.png b/p2/completion_time_distribution_paper.png new file mode 100644 index 0000000..5cf3ae7 Binary files /dev/null and b/p2/completion_time_distribution_paper.png differ diff --git a/p2/energy_distribution_paper.pdf b/p2/energy_distribution_paper.pdf new file mode 100644 index 0000000..08ba962 Binary files /dev/null and b/p2/energy_distribution_paper.pdf differ diff --git a/p2/energy_distribution_paper.png b/p2/energy_distribution_paper.png new file mode 100644 index 0000000..8c14e0c Binary files /dev/null and b/p2/energy_distribution_paper.png differ diff --git a/p2/plot_completion_time_paper.py b/p2/plot_completion_time_paper.py new file mode 100644 index 0000000..a9e51e2 --- /dev/null +++ b/p2/plot_completion_time_paper.py @@ -0,0 +1,89 @@ +""" +生成论文用的完成时间分布图 - 改进版 +""" + +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from matplotlib import rcParams +import pandas as pd + +# 设置字体 +rcParams['font.sans-serif'] = ['Arial', 'DejaVu Sans', 'Helvetica'] +rcParams['axes.unicode_minus'] = False +rcParams['font.size'] = 11 + +# 读取模拟结果数据 +df = pd.read_csv('/Volumes/Files/code/mm/20260130_b/p2/simulation_results.csv') + +# 中低饱和度配色方案 - 柔和学术风格 +colors = ['#7BAE7F', '#6A9ECF', '#9B8EC2'] # 柔和的绿、蓝、紫 + +# 方案名称 +scenario_labels = { + 'Scenario_A': 'Cost Priority', + 'Scenario_B': 'Time Priority', + 'Scenario_C': 'Balanced' +} + +# 创建图表 - 缩小尺寸以放大字体效果 +fig, axes = plt.subplots(1, 3, figsize=(10, 3.2)) + +for idx, (scenario_key, label) in enumerate(scenario_labels.items()): + ax = axes[idx] + + # 获取该方案的数据 + data = df[df['scenario'] == scenario_key]['completion_years'] + + # 计算统计量 + mean_val = np.mean(data) + p5 = np.percentile(data, 5) + p95 = np.percentile(data, 95) + + # 绘制柱状图 - 实心柱状图，增粗柱子 + ax.hist(data, bins=15, color=colors[idx], alpha=0.9, + edgecolor='white', linewidth=0.8, rwidth=0.9) + + # 设置子图标题 - 简洁的 (a) (b) (c) 标签 + ax.set_title(f'({chr(97+idx)}) {label}', fontsize=12, fontweight='normal', pad=8) + + # 设置坐标轴标签 + ax.set_xlabel('Completion Years', fontsize=11) + if idx == 0: + ax.set_ylabel('Frequency', fontsize=11) + + # 在图内右上角添加纯文本统计信息（无边框） + text_str = f'Mean: {mean_val:.1f}\n5%: {p5:.1f}\n95%: {p95:.1f}' + ax.text(0.97, 0.97, text_str, transform=ax.transAxes, fontsize=9, + verticalalignment='top', horizontalalignment='right', + bbox=dict(boxstyle='round,pad=0.3', facecolor='white', + edgecolor='none', alpha=0.85)) + + # 网格线 - 非常淡 + ax.grid(True, alpha=0.15, linestyle='-', linewidth=0.5) + + # 设置刻度字体大小 + ax.tick_params(axis='both', labelsize=10) + + # 中间图 (b) 的 x 轴使用整数刻度 + if idx == 1: + from matplotlib.ticker import MaxNLocator + ax.xaxis.set_major_locator(MaxNLocator(integer=True)) + + # 简化边框 + ax.spines['top'].set_visible(False) + ax.spines['right'].set_visible(False) + +# 调整布局 +plt.tight_layout() + +# 保存图片 +plt.savefig('/Volumes/Files/code/mm/20260130_b/p2/completion_time_distribution_paper.png', + dpi=200, bbox_inches='tight', facecolor='white') +plt.savefig('/Volumes/Files/code/mm/20260130_b/p2/completion_time_distribution_paper.pdf', + dpi=200, bbox_inches='tight', facecolor='white') + +print("图表已保存:") +print(" - completion_time_distribution_paper.png") +print(" - completion_time_distribution_paper.pdf") diff --git a/p2/plot_energy_distribution_paper.py b/p2/plot_energy_distribution_paper.py new file mode 100644 index 0000000..e339940 --- /dev/null +++ b/p2/plot_energy_distribution_paper.py @@ -0,0 +1,85 @@ +""" +生成论文用的能量分布图 - 改进版 +""" + +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from matplotlib import rcParams +from matplotlib.ticker import MaxNLocator +import pandas as pd + +# 设置字体 +rcParams['font.sans-serif'] = ['Arial', 'DejaVu Sans', 'Helvetica'] +rcParams['axes.unicode_minus'] = False +rcParams['font.size'] = 11 + +# 读取模拟结果数据 +df = pd.read_csv('/Volumes/Files/code/mm/20260130_b/p2/simulation_results.csv') + +# 中低饱和度配色方案 - 不同于之前的绿蓝紫，使用暖色调 +colors = ['#D4936A', '#6AACAC', '#C48BB8'] # 柔和的橙、青、玫瑰 + +# 方案名称 +scenario_labels = { + 'Scenario_A': 'Cost Priority', + 'Scenario_B': 'Time Priority', + 'Scenario_C': 'Balanced' +} + +# 创建图表 - 缩小尺寸以放大字体效果 +fig, axes = plt.subplots(1, 3, figsize=(10, 3.2)) + +for idx, (scenario_key, label) in enumerate(scenario_labels.items()): + ax = axes[idx] + + # 获取该方案的数据 + data = df[df['scenario'] == scenario_key]['total_energy_pj'] + + # 计算统计量 + mean_val = np.mean(data) + p5 = np.percentile(data, 5) + p95 = np.percentile(data, 95) + + # 绘制柱状图 - 实心柱状图，增粗柱子 + ax.hist(data, bins=15, color=colors[idx], alpha=0.9, + edgecolor='white', linewidth=0.8, rwidth=0.9) + + # 设置子图标题 - 简洁的 (a) (b) (c) 标签 + ax.set_title(f'({chr(97+idx)}) {label}', fontsize=12, fontweight='normal', pad=8) + + # 设置坐标轴标签 + ax.set_xlabel('Total Energy (PJ)', fontsize=11) + if idx == 0: + ax.set_ylabel('Frequency', fontsize=11) + + # 在图内右上角添加纯文本统计信息（无边框） + text_str = f'Mean: {mean_val:.0f}\n5%: {p5:.0f}\n95%: {p95:.0f}' + ax.text(0.97, 0.97, text_str, transform=ax.transAxes, fontsize=9, + verticalalignment='top', horizontalalignment='right', + bbox=dict(boxstyle='round,pad=0.3', facecolor='white', + edgecolor='none', alpha=0.85)) + + # 网格线 - 非常淡 + ax.grid(True, alpha=0.15, linestyle='-', linewidth=0.5) + + # 设置刻度字体大小 + ax.tick_params(axis='both', labelsize=10) + + # 简化边框 + ax.spines['top'].set_visible(False) + ax.spines['right'].set_visible(False) + +# 调整布局 +plt.tight_layout() + +# 保存图片 +plt.savefig('/Volumes/Files/code/mm/20260130_b/p2/energy_distribution_paper.png', + dpi=200, bbox_inches='tight', facecolor='white') +plt.savefig('/Volumes/Files/code/mm/20260130_b/p2/energy_distribution_paper.pdf', + dpi=200, bbox_inches='tight', facecolor='white') + +print("图表已保存:") +print(" - energy_distribution_paper.png") +print(" - energy_distribution_paper.pdf")