diff --git a/p1/data_window_analysis.py b/p1/data_window_analysis.py index 03cc32f..6625a5d 100644 --- a/p1/data_window_analysis.py +++ b/p1/data_window_analysis.py @@ -1,351 +1,351 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Data Window Sensitivity Analysis for Richards Model - -Analyzes how different historical data windows affect the K (carrying capacity) -estimate in the Richards growth model. - -Windows analyzed: -- 1957-2025: All available data (including Cold War era) -- 1990-2025: Post-Cold War era -- 2000-2025: Commercial space age -- 2010-2025: SpaceX era -- 2015-2025: Recent rapid growth -""" - -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 - - -# ============================================================ -# Richards Growth Model -# ============================================================ - -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) - - -# ============================================================ -# Data Loading -# ============================================================ - -def load_data(filepath="rocket_launch_counts.csv"): - """Load and preprocess 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 - - -# ============================================================ -# Fit Richards Model (Unconstrained) -# ============================================================ - -def fit_richards_unconstrained(data, base_year=1957): - """Fit Richards model without constraining K""" - years = data["year"].values - launches = data["launches"].values - t = (years - base_year).astype(float) - - # Initial parameters - p0 = [5000.0, 0.08, 80.0, 2.0] - - # Wide bounds - let data determine K - bounds = ([500, 0.005, 10, 0.2], [100000, 1.0, 200, 10.0]) - - try: - popt, pcov = curve_fit(richards, t, launches, p0=p0, bounds=bounds, maxfev=100000) - perr = np.sqrt(np.diag(pcov)) - y_pred = richards(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) - - K, r, t0, v = popt - K_err = perr[0] - - return { - "success": True, - "K": K, - "K_err": K_err, - "r": r, - "t0": t0, - "v": v, - "r_squared": r_squared, - "n_points": len(data), - "years": years, - "launches": launches, - "params": popt, - "y_pred": y_pred, - } - except Exception as e: - return {"success": False, "error": str(e)} - - -# ============================================================ -# Main Analysis -# ============================================================ - -def analyze_data_windows(df): - """Analyze K estimates for different data windows""" - - # Define windows - windows = { - "1957-2025 (全部68年)": (1957, 2025), - "1990-2025 (35年)": (1990, 2025), - "2000-2025 (25年)": (2000, 2025), - "2010-2025 (15年)": (2010, 2025), - "2015-2025 (10年)": (2015, 2025), - } - - results = {} - - print("=" * 80) - print("分析不同历史数据窗口对 K (carrying capacity) 估计的影响") - print("=" * 80) - - for name, (start, end) in windows.items(): - data = df[(df["year"] >= start) & (df["year"] <= end)].copy() - result = fit_richards_unconstrained(data) - - if result["success"]: - result["start"] = start - result["end"] = end - result["name"] = name - results[name] = result - - print(f"\n--- {name} ---") - print(f" 数据点数: {result['n_points']}") - print(f" K = {result['K']:.0f} (carrying capacity)") - print(f" r = {result['r']:.4f} (growth rate)") - print(f" R² = {result['r_squared']:.4f}") - else: - print(f"\n--- {name} ---") - print(f" 拟合失败: {result['error']}") - - return results - - -# ============================================================ -# Visualization -# ============================================================ - -def plot_window_analysis(df, results, save_path="launch_capacity_window_analysis.png"): - """Generate comprehensive window analysis plot""" - - fig, axes = plt.subplots(2, 2, figsize=(14, 11)) - - names = list(results.keys()) - physical_max = 3650 - - # ========== Plot 1: K Values Comparison ========== - ax1 = axes[0, 0] - - K_vals = [results[n]["K"] for n in names] - colors = ["#E74C3C" if k < physical_max else "#27AE60" for k in K_vals] - - y_pos = range(len(names)) - bars = ax1.barh(y_pos, K_vals, color=colors, alpha=0.7, edgecolor="black") - ax1.axvline(physical_max, color="blue", ls="--", lw=2.5, - label=f"Physical Max ({physical_max})") - - ax1.set_yticks(y_pos) - ax1.set_yticklabels(names, fontsize=10) - ax1.set_xlabel("K (Carrying Capacity, launches/year)", fontsize=11) - ax1.set_title("K Estimates by Data Window\n(Unconstrained Richards Model)", fontsize=12) - ax1.legend(loc="lower right", fontsize=10) - ax1.grid(True, alpha=0.3, axis="x") - - for bar, k in zip(bars, K_vals): - ax1.text(k + 200, bar.get_y() + bar.get_height()/2, - f"{k:.0f}", va="center", fontsize=10, fontweight="bold") - - # ========== Plot 2: R² Comparison ========== - ax2 = axes[0, 1] - - r2_vals = [results[n]["r_squared"] for n in names] - colors2 = ["#27AE60" if r2 > 0.9 else "#F39C12" if r2 > 0.7 else "#E74C3C" - for r2 in r2_vals] - - bars2 = ax2.barh(y_pos, r2_vals, color=colors2, alpha=0.7, edgecolor="black") - ax2.axvline(0.9, color="green", ls="--", lw=1.5, alpha=0.7, label="R²=0.9 (Good fit)") - ax2.axvline(0.7, color="orange", ls=":", lw=1.5, alpha=0.7, label="R²=0.7 (Acceptable)") - - ax2.set_yticks(y_pos) - ax2.set_yticklabels(names, fontsize=10) - ax2.set_xlabel("R² (Goodness of Fit)", fontsize=11) - ax2.set_title("Model Fit Quality by Data Window", fontsize=12) - ax2.set_xlim(0, 1.1) - ax2.legend(loc="lower right", fontsize=9) - ax2.grid(True, alpha=0.3, axis="x") - - for bar, r2 in zip(bars2, r2_vals): - ax2.text(r2 + 0.02, bar.get_y() + bar.get_height()/2, - f"{r2:.3f}", va="center", fontsize=10, fontweight="bold") - - # ========== Plot 3: Historical Data with All Model Fits ========== - ax3 = axes[1, 0] - - # Plot all historical data - ax3.scatter(df["year"], df["launches"], color="gray", s=30, alpha=0.5, - label="Historical Data", zorder=1) - - # Plot each model prediction - colors_line = ["#E74C3C", "#F39C12", "#3498DB", "#27AE60", "#9B59B6"] - for i, (name, res) in enumerate(results.items()): - start = res["start"] - years_proj = np.arange(start, 2080) - t_proj = years_proj - 1957 - pred = richards(t_proj, *res["params"]) - - label = f'{name.split(" ")[0]}: K={res["K"]:.0f}' - ax3.plot(years_proj, pred, color=colors_line[i], lw=2, alpha=0.8, - label=label, zorder=2) - - ax3.axhline(physical_max, color="black", ls=":", lw=2, - label=f"Physical Max ({physical_max})") - - ax3.set_xlabel("Year", fontsize=11) - ax3.set_ylabel("Annual Launches", fontsize=11) - ax3.set_title("Richards Model Predictions\n(Different Data Windows)", fontsize=12) - ax3.legend(loc="upper left", fontsize=8) - ax3.grid(True, alpha=0.3) - ax3.set_xlim(1955, 2080) - ax3.set_ylim(0, 15000) - - # ========== Plot 4: Summary Table ========== - ax4 = axes[1, 1] - ax4.axis("off") - - # Build summary text - summary = """ -┌────────────────────────────────────────────────────────────────────────┐ -│ 数据窗口对 K 估计的影响分析 │ -├────────────────────────────────────────────────────────────────────────┤ -│ │ -│ 核心发现: │ -│ ───────── │ -│ 1. 数据窗口选择对 K 估计影响巨大 │ -│ • 全部数据: K ≈ 500 │ -│ • 近15年: K ≈ 12,000 │ -│ │ -│ 2. 1957-2025 全部数据: │ -│ • R² = 0.08 (极差拟合) │ -│ • 冷战时期数据干扰，模型误判为"已饱和" │ -│ • 不适合用于预测未来增长 │ -│ │ -│ 3. 2010-2025 近15年数据: │ -│ • R² = 0.90 (良好拟合) │ -│ • 准确捕捉商业航天时代的增长趋势 │ -│ • 预测 K ≈ 12,000 >> 物理上限 (3,650) │ -│ │ -│ 结论: │ -│ ───── │ -│ • 选择 2010-2025 数据窗口是合理的 (R² = 0.90, 反映当前趋势) │ -│ • 数据驱动 K ≈ 12,000 反映增长潜力 (需要 ~33 个发射站) │ -│ • 但物理约束 (10站 × 365天 = 3,650) 才是真正的上限 │ -│ • 因此采用 1 次/站/天 作为保守估计是合理的 │ -│ │ -│ 物理上限: 3,650 次/年 (10个发射站 × 365天) │ -│ │ -└────────────────────────────────────────────────────────────────────────┘ -""" - - ax4.text(0.02, 0.98, summary, transform=ax4.transAxes, fontsize=10, - verticalalignment="top", family="monospace", - bbox=dict(boxstyle="round", facecolor="lightyellow", alpha=0.9)) - - plt.suptitle("Data Window Sensitivity Analysis for Richards Model", - fontsize=14, fontweight="bold", y=1.02) - plt.tight_layout() - plt.savefig(save_path, dpi=150, bbox_inches="tight") - plt.close() - print(f"\nPlot saved: {save_path}") - - -# ============================================================ -# Print Detailed Table -# ============================================================ - -def print_comparison_table(results): - """Print detailed comparison table""" - - physical_max = 3650 - - print("\n" + "=" * 80) - print("详细对比表") - print("=" * 80) - - header = f"{'数据窗口':<22} | {'K值':>8} | {'R²':>7} | {'数据点':>5} | {'K/物理上限':>10} | {'说明':<18}" - print(header) - print("-" * 22 + "-+-" + "-" * 8 + "-+-" + "-" * 7 + "-+-" + "-" * 5 + "-+-" + "-" * 10 + "-+-" + "-" * 18) - - for name, res in results.items(): - k = res["K"] - r2 = res["r_squared"] - n = res["n_points"] - ratio = k / physical_max - - if r2 < 0.5: - note = "拟合差，不可用" - elif k < physical_max: - note = "低于物理上限" - elif ratio < 2: - note = "接近物理上限" - else: - note = "远超物理上限" - - print(f"{name:<22} | {k:>8.0f} | {r2:>7.4f} | {n:>5} | {ratio:>10.2f}x | {note:<18}") - - print("\n物理上限: 3,650 次/年 (10个发射站 × 每站每天1次 × 365天)") - - -# ============================================================ -# Main -# ============================================================ - -def main(): - print("Loading data...") - df = load_data() - print(f"Data loaded: {len(df)} records from {df['year'].min()} to {df['year'].max()}") - - # Analyze different windows - results = analyze_data_windows(df) - - # Print comparison table - print_comparison_table(results) - - # Generate visualization - print("\n" + "=" * 80) - print("Generating visualization...") - print("=" * 80) - plot_window_analysis(df, results) - - print("\n" + "=" * 80) - print("分析完成!") - print("=" * 80) - - return results - - -if __name__ == "__main__": - results = main() +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Data Window Sensitivity Analysis for Richards Model + +Analyzes how different historical data windows affect the K (carrying capacity) +estimate in the Richards growth model. + +Windows analyzed: +- 1957-2025: All available data (including Cold War era) +- 1990-2025: Post-Cold War era +- 2000-2025: Commercial space age +- 2010-2025: SpaceX era +- 2015-2025: Recent rapid growth +""" + +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 + + +# ============================================================ +# Richards Growth Model +# ============================================================ + +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) + + +# ============================================================ +# Data Loading +# ============================================================ + +def load_data(filepath="rocket_launch_counts.csv"): + """Load and preprocess 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 + + +# ============================================================ +# Fit Richards Model (Unconstrained) +# ============================================================ + +def fit_richards_unconstrained(data, base_year=1957): + """Fit Richards model without constraining K""" + years = data["year"].values + launches = data["launches"].values + t = (years - base_year).astype(float) + + # Initial parameters + p0 = [5000.0, 0.08, 80.0, 2.0] + + # Wide bounds - let data determine K + bounds = ([500, 0.005, 10, 0.2], [100000, 1.0, 200, 10.0]) + + try: + popt, pcov = curve_fit(richards, t, launches, p0=p0, bounds=bounds, maxfev=100000) + perr = np.sqrt(np.diag(pcov)) + y_pred = richards(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) + + K, r, t0, v = popt + K_err = perr[0] + + return { + "success": True, + "K": K, + "K_err": K_err, + "r": r, + "t0": t0, + "v": v, + "r_squared": r_squared, + "n_points": len(data), + "years": years, + "launches": launches, + "params": popt, + "y_pred": y_pred, + } + except Exception as e: + return {"success": False, "error": str(e)} + + +# ============================================================ +# Main Analysis +# ============================================================ + +def analyze_data_windows(df): + """Analyze K estimates for different data windows""" + + # Define windows + windows = { + "1957-2025 (全部68年)": (1957, 2025), + "1990-2025 (35年)": (1990, 2025), + "2000-2025 (25年)": (2000, 2025), + "2010-2025 (15年)": (2010, 2025), + "2015-2025 (10年)": (2015, 2025), + } + + results = {} + + print("=" * 80) + print("分析不同历史数据窗口对 K (carrying capacity) 估计的影响") + print("=" * 80) + + for name, (start, end) in windows.items(): + data = df[(df["year"] >= start) & (df["year"] <= end)].copy() + result = fit_richards_unconstrained(data) + + if result["success"]: + result["start"] = start + result["end"] = end + result["name"] = name + results[name] = result + + print(f"\n--- {name} ---") + print(f" 数据点数: {result['n_points']}") + print(f" K = {result['K']:.0f} (carrying capacity)") + print(f" r = {result['r']:.4f} (growth rate)") + print(f" R² = {result['r_squared']:.4f}") + else: + print(f"\n--- {name} ---") + print(f" 拟合失败: {result['error']}") + + return results + + +# ============================================================ +# Visualization +# ============================================================ + +def plot_window_analysis(df, results, save_path="launch_capacity_window_analysis.png"): + """Generate comprehensive window analysis plot""" + + fig, axes = plt.subplots(2, 2, figsize=(14, 11)) + + names = list(results.keys()) + physical_max = 3650 + + # ========== Plot 1: K Values Comparison ========== + ax1 = axes[0, 0] + + K_vals = [results[n]["K"] for n in names] + colors = ["#E74C3C" if k < physical_max else "#27AE60" for k in K_vals] + + y_pos = range(len(names)) + bars = ax1.barh(y_pos, K_vals, color=colors, alpha=0.7, edgecolor="black") + ax1.axvline(physical_max, color="blue", ls="--", lw=2.5, + label=f"Physical Max ({physical_max})") + + ax1.set_yticks(y_pos) + ax1.set_yticklabels(names, fontsize=10) + ax1.set_xlabel("K (Carrying Capacity, launches/year)", fontsize=11) + ax1.set_title("K Estimates by Data Window\n(Unconstrained Richards Model)", fontsize=12) + ax1.legend(loc="lower right", fontsize=10) + ax1.grid(True, alpha=0.3, axis="x") + + for bar, k in zip(bars, K_vals): + ax1.text(k + 200, bar.get_y() + bar.get_height()/2, + f"{k:.0f}", va="center", fontsize=10, fontweight="bold") + + # ========== Plot 2: R² Comparison ========== + ax2 = axes[0, 1] + + r2_vals = [results[n]["r_squared"] for n in names] + colors2 = ["#27AE60" if r2 > 0.9 else "#F39C12" if r2 > 0.7 else "#E74C3C" + for r2 in r2_vals] + + bars2 = ax2.barh(y_pos, r2_vals, color=colors2, alpha=0.7, edgecolor="black") + ax2.axvline(0.9, color="green", ls="--", lw=1.5, alpha=0.7, label="R²=0.9 (Good fit)") + ax2.axvline(0.7, color="orange", ls=":", lw=1.5, alpha=0.7, label="R²=0.7 (Acceptable)") + + ax2.set_yticks(y_pos) + ax2.set_yticklabels(names, fontsize=10) + ax2.set_xlabel("R² (Goodness of Fit)", fontsize=11) + ax2.set_title("Model Fit Quality by Data Window", fontsize=12) + ax2.set_xlim(0, 1.1) + ax2.legend(loc="lower right", fontsize=9) + ax2.grid(True, alpha=0.3, axis="x") + + for bar, r2 in zip(bars2, r2_vals): + ax2.text(r2 + 0.02, bar.get_y() + bar.get_height()/2, + f"{r2:.3f}", va="center", fontsize=10, fontweight="bold") + + # ========== Plot 3: Historical Data with All Model Fits ========== + ax3 = axes[1, 0] + + # Plot all historical data + ax3.scatter(df["year"], df["launches"], color="gray", s=30, alpha=0.5, + label="Historical Data", zorder=1) + + # Plot each model prediction + colors_line = ["#E74C3C", "#F39C12", "#3498DB", "#27AE60", "#9B59B6"] + for i, (name, res) in enumerate(results.items()): + start = res["start"] + years_proj = np.arange(start, 2080) + t_proj = years_proj - 1957 + pred = richards(t_proj, *res["params"]) + + label = f'{name.split(" ")[0]}: K={res["K"]:.0f}' + ax3.plot(years_proj, pred, color=colors_line[i], lw=2, alpha=0.8, + label=label, zorder=2) + + ax3.axhline(physical_max, color="black", ls=":", lw=2, + label=f"Physical Max ({physical_max})") + + ax3.set_xlabel("Year", fontsize=11) + ax3.set_ylabel("Annual Launches", fontsize=11) + ax3.set_title("Richards Model Predictions\n(Different Data Windows)", fontsize=12) + ax3.legend(loc="upper left", fontsize=8) + ax3.grid(True, alpha=0.3) + ax3.set_xlim(1955, 2080) + ax3.set_ylim(0, 15000) + + # ========== Plot 4: Summary Table ========== + ax4 = axes[1, 1] + ax4.axis("off") + + # Build summary text + summary = """ +┌────────────────────────────────────────────────────────────────────────┐ +│ 数据窗口对 K 估计的影响分析 │ +├────────────────────────────────────────────────────────────────────────┤ +│ │ +│ 核心发现: │ +│ ───────── │ +│ 1. 数据窗口选择对 K 估计影响巨大 │ +│ • 全部数据: K ≈ 500 │ +│ • 近15年: K ≈ 12,000 │ +│ │ +│ 2. 1957-2025 全部数据: │ +│ • R² = 0.08 (极差拟合) │ +│ • 冷战时期数据干扰，模型误判为"已饱和" │ +│ • 不适合用于预测未来增长 │ +│ │ +│ 3. 2010-2025 近15年数据: │ +│ • R² = 0.90 (良好拟合) │ +│ • 准确捕捉商业航天时代的增长趋势 │ +│ • 预测 K ≈ 12,000 >> 物理上限 (3,650) │ +│ │ +│ 结论: │ +│ ───── │ +│ • 选择 2010-2025 数据窗口是合理的 (R² = 0.90, 反映当前趋势) │ +│ • 数据驱动 K ≈ 12,000 反映增长潜力 (需要 ~33 个发射站) │ +│ • 但物理约束 (10站 × 365天 = 3,650) 才是真正的上限 │ +│ • 因此采用 1 次/站/天 作为保守估计是合理的 │ +│ │ +│ 物理上限: 3,650 次/年 (10个发射站 × 365天) │ +│ │ +└────────────────────────────────────────────────────────────────────────┘ +""" + + ax4.text(0.02, 0.98, summary, transform=ax4.transAxes, fontsize=10, + verticalalignment="top", family="monospace", + bbox=dict(boxstyle="round", facecolor="lightyellow", alpha=0.9)) + + plt.suptitle("Data Window Sensitivity Analysis for Richards Model", + fontsize=14, fontweight="bold", y=1.02) + plt.tight_layout() + plt.savefig(save_path, dpi=150, bbox_inches="tight") + plt.close() + print(f"\nPlot saved: {save_path}") + + +# ============================================================ +# Print Detailed Table +# ============================================================ + +def print_comparison_table(results): + """Print detailed comparison table""" + + physical_max = 3650 + + print("\n" + "=" * 80) + print("详细对比表") + print("=" * 80) + + header = f"{'数据窗口':<22} | {'K值':>8} | {'R²':>7} | {'数据点':>5} | {'K/物理上限':>10} | {'说明':<18}" + print(header) + print("-" * 22 + "-+-" + "-" * 8 + "-+-" + "-" * 7 + "-+-" + "-" * 5 + "-+-" + "-" * 10 + "-+-" + "-" * 18) + + for name, res in results.items(): + k = res["K"] + r2 = res["r_squared"] + n = res["n_points"] + ratio = k / physical_max + + if r2 < 0.5: + note = "拟合差，不可用" + elif k < physical_max: + note = "低于物理上限" + elif ratio < 2: + note = "接近物理上限" + else: + note = "远超物理上限" + + print(f"{name:<22} | {k:>8.0f} | {r2:>7.4f} | {n:>5} | {ratio:>10.2f}x | {note:<18}") + + print("\n物理上限: 3,650 次/年 (10个发射站 × 每站每天1次 × 365天)") + + +# ============================================================ +# Main +# ============================================================ + +def main(): + print("Loading data...") + df = load_data() + print(f"Data loaded: {len(df)} records from {df['year'].min()} to {df['year'].max()}") + + # Analyze different windows + results = analyze_data_windows(df) + + # Print comparison table + print_comparison_table(results) + + # Generate visualization + print("\n" + "=" * 80) + print("Generating visualization...") + print("=" * 80) + plot_window_analysis(df, results) + + print("\n" + "=" * 80) + print("分析完成!") + print("=" * 80) + + return results + + +if __name__ == "__main__": + results = main() diff --git a/p1/find_optimal_window.py b/p1/find_optimal_window.py index 5acb8d3..fe9f329 100644 --- a/p1/find_optimal_window.py +++ b/p1/find_optimal_window.py @@ -1,230 +1,230 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Find Optimal Data Window for K ≈ 3650 - -Enumerate all possible starting years to find which data window -produces K closest to the physical limit (3650 launches/year). -""" - -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 - - -# Richards model -def richards(t, K, r, t0, v): - 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"): - 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 fit_richards(data, base_year=1957): - """Fit unconstrained Richards model""" - years = data["year"].values - launches = data["launches"].values - t = (years - base_year).astype(float) - - p0 = [5000.0, 0.08, 80.0, 2.0] - bounds = ([500, 0.005, 10, 0.2], [100000, 1.0, 200, 10.0]) - - try: - popt, pcov = curve_fit(richards, t, launches, p0=p0, bounds=bounds, maxfev=100000) - y_pred = richards(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) if ss_tot > 0 else 0 - - return { - "success": True, - "K": popt[0], - "r": popt[1], - "t0": popt[2], - "v": popt[3], - "r_squared": r_squared, - "n_points": len(data), - } - except: - return {"success": False} - - -def main(): - print("=" * 80) - print("枚举起始年份，寻找 K ≈ 3650 的数据窗口") - print("=" * 80) - - df = load_data() - print(f"数据范围: {df['year'].min()} - {df['year'].max()}\n") - - target_K = 3650 - end_year = 2025 - - # Enumerate all starting years - results = [] - - for start_year in range(1957, 2020): - data = df[(df["year"] >= start_year) & (df["year"] <= end_year)].copy() - - if len(data) < 5: # Need at least 5 data points - continue - - result = fit_richards(data) - - if result["success"]: - diff = abs(result["K"] - target_K) - results.append({ - "start_year": start_year, - "end_year": end_year, - "years_span": end_year - start_year + 1, - "n_points": result["n_points"], - "K": result["K"], - "r_squared": result["r_squared"], - "diff_from_target": diff, - "ratio": result["K"] / target_K, - }) - - # Sort by difference from target - results_df = pd.DataFrame(results) - results_df = results_df.sort_values("diff_from_target") - - # Print full table - print("完整枚举结果表（按与目标值 3650 的差距排序）:") - print("-" * 80) - print(f"{'起始年份':>8} | {'时间跨度':>8} | {'数据点':>6} | {'K值':>10} | {'R²':>8} | {'K/3650':>8} | {'差距':>10}") - print("-" * 80) - - for _, row in results_df.iterrows(): - marker = "★" if row["diff_from_target"] < 500 else "" - print(f"{int(row['start_year']):>8} | {int(row['years_span']):>6}年 | {int(row['n_points']):>6} | " - f"{row['K']:>10.0f} | {row['r_squared']:>8.4f} | {row['ratio']:>8.2f}x | {row['diff_from_target']:>10.0f} {marker}") - - # Find closest matches - print("\n" + "=" * 80) - print("最接近 K = 3650 的数据窗口 (Top 10)") - print("=" * 80) - - top10 = results_df.head(10) - print(f"\n{'排名':>4} | {'起始年份':>8} | {'K值':>10} | {'R²':>8} | {'差距':>10} | {'评价':<20}") - print("-" * 75) - - for i, (_, row) in enumerate(top10.iterrows(), 1): - if row["r_squared"] < 0.5: - comment = "❌ 拟合差，不可信" - elif row["r_squared"] < 0.8: - comment = "⚠️ 拟合一般" - else: - comment = "✅ 拟合良好" - - print(f"{i:>4} | {int(row['start_year']):>8} | {row['K']:>10.0f} | {row['r_squared']:>8.4f} | " - f"{row['diff_from_target']:>10.0f} | {comment:<20}") - - # Analysis - print("\n" + "=" * 80) - print("分析结论") - print("=" * 80) - - # Find best with good R² - good_fit = results_df[results_df["r_squared"] >= 0.8] - if len(good_fit) > 0: - best_good = good_fit.iloc[0] - print(f"\n在 R² ≥ 0.8 的条件下，最接近 K=3650 的窗口:") - print(f" 起始年份: {int(best_good['start_year'])}") - print(f" K = {best_good['K']:.0f}") - print(f" R² = {best_good['r_squared']:.4f}") - print(f" 差距: {best_good['diff_from_target']:.0f}") - - # Summary - print("\n关键发现:") - print("-" * 40) - - # Check if any good fit gives K near 3650 - near_target = results_df[(results_df["diff_from_target"] < 1000) & (results_df["r_squared"] >= 0.7)] - - if len(near_target) == 0: - print(" ⚠️ 没有任何数据窗口能在良好拟合(R²≥0.7)的条件下得到 K≈3650") - print(" ⚠️ 所有良好拟合的窗口都给出 K >> 3650 或 K << 3650") - print("\n 这说明:") - print(" • K=3650 不是数据自然支持的结论") - print(" • K=3650 来自物理约束，而非统计预测") - print(" • 论文中应该明确说明这是'物理上限'而非'数据预测'") - else: - print(f" 找到 {len(near_target)} 个窗口使 K 接近 3650:") - for _, row in near_target.iterrows(): - print(f" {int(row['start_year'])}-2025: K={row['K']:.0f}, R²={row['r_squared']:.4f}") - - # Generate visualization - print("\n" + "=" * 80) - print("生成可视化...") - print("=" * 80) - - fig, axes = plt.subplots(1, 2, figsize=(14, 6)) - - # Plot 1: K vs Start Year - ax1 = axes[0] - colors = ['#27AE60' if r2 >= 0.8 else '#F39C12' if r2 >= 0.5 else '#E74C3C' - for r2 in results_df["r_squared"]] - - ax1.scatter(results_df["start_year"], results_df["K"], c=colors, s=60, alpha=0.7, edgecolor='black') - ax1.axhline(3650, color='blue', ls='--', lw=2, label='Target K=3650') - ax1.axhline(3650*0.9, color='blue', ls=':', lw=1, alpha=0.5) - ax1.axhline(3650*1.1, color='blue', ls=':', lw=1, alpha=0.5) - - ax1.set_xlabel("Starting Year", fontsize=11) - ax1.set_ylabel("K (Carrying Capacity)", fontsize=11) - ax1.set_title("K vs Starting Year\n(Color: Green=R²≥0.8, Yellow=R²≥0.5, Red=R²<0.5)", fontsize=12) - ax1.legend() - ax1.grid(True, alpha=0.3) - ax1.set_xlim(1955, 2020) - - # Plot 2: R² vs Start Year - ax2 = axes[1] - ax2.scatter(results_df["start_year"], results_df["r_squared"], c=colors, s=60, alpha=0.7, edgecolor='black') - ax2.axhline(0.8, color='green', ls='--', lw=1.5, label='R²=0.8 (Good)') - ax2.axhline(0.5, color='orange', ls=':', lw=1.5, label='R²=0.5 (Acceptable)') - - ax2.set_xlabel("Starting Year", fontsize=11) - ax2.set_ylabel("R² (Goodness of Fit)", fontsize=11) - ax2.set_title("Model Fit Quality vs Starting Year", fontsize=12) - ax2.legend() - ax2.grid(True, alpha=0.3) - ax2.set_xlim(1955, 2020) - ax2.set_ylim(-0.1, 1.1) - - plt.tight_layout() - plt.savefig("find_optimal_window.png", dpi=150, bbox_inches='tight') - plt.close() - print("图表已保存: find_optimal_window.png") - - # Save to CSV - results_df.to_csv("window_enumeration.csv", index=False) - print("数据已保存: window_enumeration.csv") - - print("\n" + "=" * 80) - print("分析完成!") - print("=" * 80) - - return results_df - - -if __name__ == "__main__": - results = main() +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Find Optimal Data Window for K ≈ 3650 + +Enumerate all possible starting years to find which data window +produces K closest to the physical limit (3650 launches/year). +""" + +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 + + +# Richards model +def richards(t, K, r, t0, v): + 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"): + 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 fit_richards(data, base_year=1957): + """Fit unconstrained Richards model""" + years = data["year"].values + launches = data["launches"].values + t = (years - base_year).astype(float) + + p0 = [5000.0, 0.08, 80.0, 2.0] + bounds = ([500, 0.005, 10, 0.2], [100000, 1.0, 200, 10.0]) + + try: + popt, pcov = curve_fit(richards, t, launches, p0=p0, bounds=bounds, maxfev=100000) + y_pred = richards(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) if ss_tot > 0 else 0 + + return { + "success": True, + "K": popt[0], + "r": popt[1], + "t0": popt[2], + "v": popt[3], + "r_squared": r_squared, + "n_points": len(data), + } + except: + return {"success": False} + + +def main(): + print("=" * 80) + print("枚举起始年份，寻找 K ≈ 3650 的数据窗口") + print("=" * 80) + + df = load_data() + print(f"数据范围: {df['year'].min()} - {df['year'].max()}\n") + + target_K = 3650 + end_year = 2025 + + # Enumerate all starting years + results = [] + + for start_year in range(1957, 2020): + data = df[(df["year"] >= start_year) & (df["year"] <= end_year)].copy() + + if len(data) < 5: # Need at least 5 data points + continue + + result = fit_richards(data) + + if result["success"]: + diff = abs(result["K"] - target_K) + results.append({ + "start_year": start_year, + "end_year": end_year, + "years_span": end_year - start_year + 1, + "n_points": result["n_points"], + "K": result["K"], + "r_squared": result["r_squared"], + "diff_from_target": diff, + "ratio": result["K"] / target_K, + }) + + # Sort by difference from target + results_df = pd.DataFrame(results) + results_df = results_df.sort_values("diff_from_target") + + # Print full table + print("完整枚举结果表（按与目标值 3650 的差距排序）:") + print("-" * 80) + print(f"{'起始年份':>8} | {'时间跨度':>8} | {'数据点':>6} | {'K值':>10} | {'R²':>8} | {'K/3650':>8} | {'差距':>10}") + print("-" * 80) + + for _, row in results_df.iterrows(): + marker = "★" if row["diff_from_target"] < 500 else "" + print(f"{int(row['start_year']):>8} | {int(row['years_span']):>6}年 | {int(row['n_points']):>6} | " + f"{row['K']:>10.0f} | {row['r_squared']:>8.4f} | {row['ratio']:>8.2f}x | {row['diff_from_target']:>10.0f} {marker}") + + # Find closest matches + print("\n" + "=" * 80) + print("最接近 K = 3650 的数据窗口 (Top 10)") + print("=" * 80) + + top10 = results_df.head(10) + print(f"\n{'排名':>4} | {'起始年份':>8} | {'K值':>10} | {'R²':>8} | {'差距':>10} | {'评价':<20}") + print("-" * 75) + + for i, (_, row) in enumerate(top10.iterrows(), 1): + if row["r_squared"] < 0.5: + comment = "❌ 拟合差，不可信" + elif row["r_squared"] < 0.8: + comment = "⚠️ 拟合一般" + else: + comment = "✅ 拟合良好" + + print(f"{i:>4} | {int(row['start_year']):>8} | {row['K']:>10.0f} | {row['r_squared']:>8.4f} | " + f"{row['diff_from_target']:>10.0f} | {comment:<20}") + + # Analysis + print("\n" + "=" * 80) + print("分析结论") + print("=" * 80) + + # Find best with good R² + good_fit = results_df[results_df["r_squared"] >= 0.8] + if len(good_fit) > 0: + best_good = good_fit.iloc[0] + print(f"\n在 R² ≥ 0.8 的条件下，最接近 K=3650 的窗口:") + print(f" 起始年份: {int(best_good['start_year'])}") + print(f" K = {best_good['K']:.0f}") + print(f" R² = {best_good['r_squared']:.4f}") + print(f" 差距: {best_good['diff_from_target']:.0f}") + + # Summary + print("\n关键发现:") + print("-" * 40) + + # Check if any good fit gives K near 3650 + near_target = results_df[(results_df["diff_from_target"] < 1000) & (results_df["r_squared"] >= 0.7)] + + if len(near_target) == 0: + print(" ⚠️ 没有任何数据窗口能在良好拟合(R²≥0.7)的条件下得到 K≈3650") + print(" ⚠️ 所有良好拟合的窗口都给出 K >> 3650 或 K << 3650") + print("\n 这说明:") + print(" • K=3650 不是数据自然支持的结论") + print(" • K=3650 来自物理约束，而非统计预测") + print(" • 论文中应该明确说明这是'物理上限'而非'数据预测'") + else: + print(f" 找到 {len(near_target)} 个窗口使 K 接近 3650:") + for _, row in near_target.iterrows(): + print(f" {int(row['start_year'])}-2025: K={row['K']:.0f}, R²={row['r_squared']:.4f}") + + # Generate visualization + print("\n" + "=" * 80) + print("生成可视化...") + print("=" * 80) + + fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + + # Plot 1: K vs Start Year + ax1 = axes[0] + colors = ['#27AE60' if r2 >= 0.8 else '#F39C12' if r2 >= 0.5 else '#E74C3C' + for r2 in results_df["r_squared"]] + + ax1.scatter(results_df["start_year"], results_df["K"], c=colors, s=60, alpha=0.7, edgecolor='black') + ax1.axhline(3650, color='blue', ls='--', lw=2, label='Target K=3650') + ax1.axhline(3650*0.9, color='blue', ls=':', lw=1, alpha=0.5) + ax1.axhline(3650*1.1, color='blue', ls=':', lw=1, alpha=0.5) + + ax1.set_xlabel("Starting Year", fontsize=11) + ax1.set_ylabel("K (Carrying Capacity)", fontsize=11) + ax1.set_title("K vs Starting Year\n(Color: Green=R²≥0.8, Yellow=R²≥0.5, Red=R²<0.5)", fontsize=12) + ax1.legend() + ax1.grid(True, alpha=0.3) + ax1.set_xlim(1955, 2020) + + # Plot 2: R² vs Start Year + ax2 = axes[1] + ax2.scatter(results_df["start_year"], results_df["r_squared"], c=colors, s=60, alpha=0.7, edgecolor='black') + ax2.axhline(0.8, color='green', ls='--', lw=1.5, label='R²=0.8 (Good)') + ax2.axhline(0.5, color='orange', ls=':', lw=1.5, label='R²=0.5 (Acceptable)') + + ax2.set_xlabel("Starting Year", fontsize=11) + ax2.set_ylabel("R² (Goodness of Fit)", fontsize=11) + ax2.set_title("Model Fit Quality vs Starting Year", fontsize=12) + ax2.legend() + ax2.grid(True, alpha=0.3) + ax2.set_xlim(1955, 2020) + ax2.set_ylim(-0.1, 1.1) + + plt.tight_layout() + plt.savefig("find_optimal_window.png", dpi=150, bbox_inches='tight') + plt.close() + print("图表已保存: find_optimal_window.png") + + # Save to CSV + results_df.to_csv("window_enumeration.csv", index=False) + print("数据已保存: window_enumeration.csv") + + print("\n" + "=" * 80) + print("分析完成!") + print("=" * 80) + + return results_df + + +if __name__ == "__main__": + results = main() diff --git a/p1/lambda_cost_analysis.py b/p1/lambda_cost_analysis.py index 9233ca0..4ead3b9 100644 --- a/p1/lambda_cost_analysis.py +++ b/p1/lambda_cost_analysis.py @@ -1,844 +1,844 @@ -""" -Lambda Time Cost Analysis for Moon Colony Logistics - -This module introduces the time opportunity cost (λ) to find the optimal -operating point on the Energy-Time trade-off curve. - -Objective Function: J = E_total + λ × T -where: - - E_total: Total energy consumption (PJ) - - T: Construction timeline (years) - - λ: Time opportunity cost (PJ/year) - -The optimal point satisfies: dE/dT = -λ (marginal energy saving = time cost) -""" - -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 - -# 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 - -# Mission parameters -TOTAL_PAYLOAD = 100e6 # 100 million metric tons - -# Space Elevator parameters -NUM_ELEVATORS = 3 -ELEVATOR_CAPACITY_PER_YEAR = 179000 # metric tons per elevator per year -TOTAL_ELEVATOR_CAPACITY = NUM_ELEVATORS * ELEVATOR_CAPACITY_PER_YEAR # 537,000 tons/year -ELEVATOR_SPECIFIC_ENERGY = 157.2e9 # J per metric ton (157.2 MJ/kg × 1000) - -# Rocket parameters - LOX/CH4 (Raptor-class) -PAYLOAD_PER_LAUNCH = 125 # metric tons per launch -ISP = 363 # Specific impulse (s) -SPECIFIC_FUEL_ENERGY = 11.9e6 # J/kg -ALPHA = 0.10 # Structural coefficient -NUM_STAGES = 3 -DELTA_V_BASE = 13300 # m/s (LEO + TLI + LOI) - - -# ============== Launch Site Definition ============== -@dataclass -class LaunchSite: - name: str - short_name: str - latitude: float - max_launches_per_day: int = 1 - - @property - def abs_latitude(self) -> float: - return abs(self.latitude) - - @property - def delta_v_loss(self) -> float: - v_equator = OMEGA_EARTH * R_EARTH - v_site = OMEGA_EARTH * R_EARTH * np.cos(np.radians(self.abs_latitude)) - return v_equator - v_site - - @property - def total_delta_v(self) -> float: - return DELTA_V_BASE + self.delta_v_loss - - -LAUNCH_SITES = sorted([ - LaunchSite("Kourou (French Guiana)", "Kourou", 5.2), - LaunchSite("Satish Dhawan (India)", "SDSC", 13.7), - LaunchSite("Boca Chica (Texas)", "Texas", 26.0), - LaunchSite("Cape Canaveral (Florida)", "Florida", 28.5), - LaunchSite("Vandenberg (California)", "California", 34.7), - LaunchSite("Wallops (Virginia)", "Virginia", 37.8), - LaunchSite("Taiyuan (China)", "Taiyuan", 38.8), - LaunchSite("Mahia (New Zealand)", "Mahia", 39.3), - LaunchSite("Baikonur (Kazakhstan)", "Baikonur", 45.6), - LaunchSite("Kodiak (Alaska)", "Alaska", 57.4), -], key=lambda x: x.abs_latitude) - - -# ============== Core Calculation Functions ============== - -def fuel_ratio_multistage(delta_v: float) -> float: - """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 rocket_energy_per_ton(site: LaunchSite) -> float: - """Energy consumption per ton of payload for rocket launch (J/ton)""" - k = fuel_ratio_multistage(site.total_delta_v) - fuel_per_ton = k * 1000 # kg fuel per metric ton payload - return fuel_per_ton * SPECIFIC_FUEL_ENERGY - - -def calculate_scenario(completion_years: float) -> Optional[Dict]: - """ - Calculate optimal scenario for given completion timeline - (elevator priority + low-latitude rockets) - """ - # Space elevator transport - elevator_payload = min(TOTAL_ELEVATOR_CAPACITY * completion_years, TOTAL_PAYLOAD) - elevator_energy = elevator_payload * 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, - 'elevator_energy_PJ': elevator_energy / 1e15, - 'rocket_energy_PJ': 0, - 'total_energy_PJ': elevator_energy / 1e15, - 'rocket_launches': 0, - 'sites_used': 0, - 'elevator_fraction': 1.0, - } - - # Rocket launches needed - rocket_launches_needed = int(np.ceil(remaining_payload / PAYLOAD_PER_LAUNCH)) - - # Allocate by latitude priority - days_available = completion_years * 365 - max_launches_per_site = int(days_available) - - # Check feasibility - total_rocket_capacity = len(LAUNCH_SITES) * max_launches_per_site * PAYLOAD_PER_LAUNCH - if remaining_payload > total_rocket_capacity: - return None - - rocket_energy = 0 - sites_used = 0 - remaining_launches = rocket_launches_needed - - for site in LAUNCH_SITES: - if remaining_launches <= 0: - break - allocated = min(remaining_launches, max_launches_per_site) - rocket_energy += rocket_energy_per_ton(site) * PAYLOAD_PER_LAUNCH * allocated - remaining_launches -= allocated - if allocated > 0: - sites_used += 1 - - rocket_payload = rocket_launches_needed * PAYLOAD_PER_LAUNCH - total_energy = elevator_energy + rocket_energy - - return { - 'years': completion_years, - 'elevator_payload': elevator_payload, - 'rocket_payload': rocket_payload, - 'elevator_energy_PJ': elevator_energy / 1e15, - 'rocket_energy_PJ': rocket_energy / 1e15, - 'total_energy_PJ': total_energy / 1e15, - 'rocket_launches': rocket_launches_needed, - 'sites_used': sites_used, - 'elevator_fraction': elevator_payload / TOTAL_PAYLOAD, - } - - -# ============== Generate Trade-off Curve ============== - -def generate_tradeoff_curve( - year_min: float = 100, - year_max: float = 250, - num_points: int = 500 -) -> pd.DataFrame: - """Generate Energy-Time trade-off curve data""" - years_range = np.linspace(year_min, year_max, num_points) - - results = [] - for years in years_range: - scenario = calculate_scenario(years) - if scenario is not None: - results.append({ - 'years': years, - 'energy_PJ': scenario['total_energy_PJ'], - 'elevator_fraction': scenario['elevator_fraction'], - 'sites_used': scenario['sites_used'], - 'rocket_launches': scenario['rocket_launches'], - }) - - return pd.DataFrame(results) - - -# ============== Lambda Cost Analysis ============== - -def calculate_total_cost(df: pd.DataFrame, lambda_cost: float) -> np.ndarray: - """ - Calculate total cost J = E + λ × T - - Args: - df: Trade-off curve data - lambda_cost: Time opportunity cost (PJ/year) - - Returns: - Total cost array - """ - return df['energy_PJ'].values + lambda_cost * df['years'].values - - -def find_optimal_point(df: pd.DataFrame, lambda_cost: float) -> Dict: - """ - Find optimal point that minimizes J = E + λ × T - - Args: - df: Trade-off curve data - lambda_cost: Time opportunity cost (PJ/year) - - Returns: - Optimal point information - """ - total_cost = calculate_total_cost(df, lambda_cost) - opt_idx = np.argmin(total_cost) - - return { - 'index': opt_idx, - 'years': df['years'].iloc[opt_idx], - 'energy_PJ': df['energy_PJ'].iloc[opt_idx], - 'total_cost': total_cost[opt_idx], - 'elevator_fraction': df['elevator_fraction'].iloc[opt_idx], - 'lambda': lambda_cost, - } - - -def calculate_marginal_energy_saving(df: pd.DataFrame) -> np.ndarray: - """ - Calculate marginal energy saving rate: -dE/dT (PJ/year) - - This represents how much energy is saved per additional year of timeline. - """ - years = df['years'].values - energy = df['energy_PJ'].values - - # Use central difference for interior points - marginal = -np.gradient(energy, years) - - return marginal - - -def sensitivity_analysis( - df: pd.DataFrame, - lambda_range: np.ndarray -) -> pd.DataFrame: - """ - Perform sensitivity analysis on λ parameter - - Args: - df: Trade-off curve data - lambda_range: Array of λ values to test - - Returns: - DataFrame with optimal points for each λ - """ - results = [] - for lam in lambda_range: - opt = find_optimal_point(df, lam) - results.append({ - 'lambda_PJ_per_year': lam, - 'optimal_years': opt['years'], - 'optimal_energy_PJ': opt['energy_PJ'], - 'total_cost_PJ': opt['total_cost'], - 'elevator_fraction': opt['elevator_fraction'], - }) - - return pd.DataFrame(results) - - -# ============== Visualization Functions ============== - -def plot_lambda_analysis( - df: pd.DataFrame, - save_path: str = '/Volumes/Files/code/mm/20260130_b/p1/lambda_cost_analysis.png' -): - """ - Comprehensive visualization of λ time cost analysis - Focus on critical range λ = 400-600 PJ/year - """ - fig, axes = plt.subplots(2, 2, figsize=(14, 12)) - - years = df['years'].values - energy = df['energy_PJ'].values - - # Key boundaries - T_min = years.min() # ~100.7 years (fastest) - T_elev = TOTAL_PAYLOAD / TOTAL_ELEVATOR_CAPACITY # ~186 years (elevator only) - E_min = energy[years >= T_elev].min() if any(years >= T_elev) else energy.min() - - # ========== Plot 1: Trade-off Curve with λ Iso-cost Lines ========== - ax1 = axes[0, 0] - - # Main curve - ax1.plot(years, energy, 'b-', linewidth=2.5, label='Energy-Time Trade-off') - - # Mark key points - ax1.axvline(x=T_elev, color='green', linestyle='--', alpha=0.7, label=f'Elevator-only: {T_elev:.1f} years') - ax1.axvline(x=T_min, color='red', linestyle='--', alpha=0.7, label=f'Minimum time: {T_min:.1f} years') - - # Draw iso-cost lines for different λ (focus on 400-600 range) - lambda_values = [420, 480, 500, 550] - colors = ['#2ecc71', '#e74c3c', '#9b59b6', '#3498db'] - - for lam, color in zip(lambda_values, colors): - opt = find_optimal_point(df, lam) - # Iso-cost line: E + λT = const → E = const - λT - T_line = np.linspace(80, 220, 100) - E_line = opt['total_cost'] - lam * T_line - valid = (E_line > 0) & (E_line < 70000) - ax1.plot(T_line[valid], E_line[valid], '--', color=color, alpha=0.6, linewidth=1.5) - ax1.plot(opt['years'], opt['energy_PJ'], 'o', color=color, markersize=12, - markeredgecolor='black', markeredgewidth=1.5, - label=f'λ={lam}: T={opt["years"]:.1f}y, E={opt["energy_PJ"]:.0f}PJ') - - ax1.set_xlabel('Construction Timeline T (years)', fontsize=12) - ax1.set_ylabel('Total Energy E (PJ)', fontsize=12) - ax1.set_title('Energy-Time Trade-off with λ Iso-cost Lines (λ=400-600)\n$J = E + λT$', fontsize=13) - ax1.legend(loc='upper right', fontsize=9) - ax1.grid(True, alpha=0.3) - ax1.set_xlim(95, 200) - ax1.set_ylim(10000, 65000) - - # ========== Plot 2: Marginal Energy Saving Rate ========== - ax2 = axes[0, 1] - - marginal = calculate_marginal_energy_saving(df) - - ax2.plot(years, marginal, 'r-', linewidth=2, label='Marginal Energy Saving -dE/dT') - ax2.axhline(y=0, color='black', linestyle='-', alpha=0.3) - - # Mark critical λ values in 400-600 range - for lam, color in zip([450, 480, 500, 550], ['#2ecc71', '#e74c3c', '#9b59b6', '#3498db']): - ax2.axhline(y=lam, color=color, linestyle='--', alpha=0.7, label=f'λ = {lam} PJ/year') - # Find intersection - intersect_idx = np.argmin(np.abs(marginal - lam)) - ax2.plot(years[intersect_idx], marginal[intersect_idx], 'o', color=color, markersize=8) - - ax2.set_xlabel('Construction Timeline T (years)', fontsize=12) - ax2.set_ylabel('Marginal Energy Saving -dE/dT (PJ/year)', fontsize=12) - ax2.set_title('Marginal Analysis: Optimal Condition is -dE/dT = λ', fontsize=13) - ax2.legend(loc='upper right', fontsize=9) - ax2.grid(True, alpha=0.3) - ax2.set_xlim(95, 200) - ax2.set_ylim(350, 620) - - # ========== Plot 3: Sensitivity Analysis (400-600 range) ========== - ax3 = axes[1, 0] - - lambda_range = np.linspace(400, 600, 200) - sensitivity_df = sensitivity_analysis(df, lambda_range) - - # Plot optimal years vs λ - ax3.plot(sensitivity_df['lambda_PJ_per_year'], sensitivity_df['optimal_years'], - 'b-', linewidth=2.5, label='Optimal Timeline T*') - - # Mark key regions - ax3.axhline(y=T_elev, color='green', linestyle='--', alpha=0.7, label=f'Elevator-only: {T_elev:.1f}y') - ax3.axhline(y=T_min, color='red', linestyle='--', alpha=0.7, label=f'Min timeline: {T_min:.1f}y') - ax3.axhline(y=139, color='orange', linestyle=':', linewidth=2, alpha=0.8, label='Original knee: 139y') - - # Find λ corresponding to 139 years - idx_139 = np.argmin(np.abs(sensitivity_df['optimal_years'] - 139)) - if idx_139 < len(sensitivity_df): - lambda_139 = sensitivity_df['lambda_PJ_per_year'].iloc[idx_139] - ax3.axvline(x=lambda_139, color='orange', linestyle=':', linewidth=2, alpha=0.6) - ax3.scatter([lambda_139], [139], s=150, c='orange', marker='*', zorder=5, - edgecolors='black', linewidths=1.5, label=f'λ≈{lambda_139:.0f} for T*=139y') - - # Mark critical transition - ax3.axvline(x=480, color='purple', linestyle='--', alpha=0.5) - ax3.text(482, 175, 'Critical\nTransition', fontsize=9, color='purple') - - ax3.set_xlabel('Time Opportunity Cost λ (PJ/year)', fontsize=12) - ax3.set_ylabel('Optimal Timeline T* (years)', fontsize=12) - ax3.set_title('Sensitivity Analysis: Optimal Timeline vs λ (400-600 range)', fontsize=13) - ax3.legend(loc='upper right', fontsize=9) - ax3.grid(True, alpha=0.3) - ax3.set_xlim(400, 600) - ax3.set_ylim(95, 195) - - # ========== Plot 4: Total Cost Curves (400-600 range) ========== - ax4 = axes[1, 1] - - for lam, color in [(450, '#2ecc71'), (480, '#e74c3c'), (500, '#9b59b6'), (550, '#3498db')]: - total_cost = calculate_total_cost(df, lam) - ax4.plot(years, total_cost / 1000, color=color, linestyle='-', - linewidth=2, label=f'λ={lam} PJ/year') - - # Mark minimum - opt = find_optimal_point(df, lam) - ax4.plot(opt['years'], opt['total_cost'] / 1000, 'o', color=color, markersize=10, - markeredgecolor='black', markeredgewidth=1.5) - - ax4.set_xlabel('Construction Timeline T (years)', fontsize=12) - ax4.set_ylabel('Total Cost J = E + λT (×10³ PJ)', fontsize=12) - ax4.set_title('Total Cost Function for λ = 400-600 PJ/year', fontsize=13) - ax4.legend(loc='upper right', fontsize=10) - 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"Lambda analysis plot saved to: {save_path}") - - return fig - - -def plot_decision_recommendation( - df: pd.DataFrame, - save_path: str = '/Volumes/Files/code/mm/20260130_b/p1/lambda_decision_map.png' -): - """ - Decision map showing optimal choice based on λ preference - """ - fig, ax = plt.subplots(figsize=(12, 8)) - - years = df['years'].values - energy = df['energy_PJ'].values - - # Main trade-off curve - ax.plot(years, energy, 'b-', linewidth=3, label='Feasible Trade-off Curve') - - # Key boundaries - T_elev = TOTAL_PAYLOAD / TOTAL_ELEVATOR_CAPACITY - T_min = years.min() - - # Shade decision regions - # Region 1: Low λ (cost priority) → longer timeline - ax.axvspan(160, 190, alpha=0.2, color='green', label='Low λ (<150): Cost Priority') - # Region 2: Medium λ (balanced) → middle ground - ax.axvspan(130, 160, alpha=0.2, color='yellow', label='Medium λ (150-250): Balanced') - # Region 3: High λ (time priority) → shorter timeline - ax.axvspan(100, 130, alpha=0.2, color='red', label='High λ (>250): Time Priority') - - # Mark specific strategy points - strategies = [ - {'name': 'A: Cost-Prioritized', 'years': 186, 'lambda': 0, 'color': 'green'}, - {'name': 'C: Balanced (λ≈200)', 'years': 139, 'lambda': 200, 'color': 'orange'}, - {'name': 'B: Time-Prioritized', 'years': 101, 'lambda': 500, 'color': 'red'}, - ] - - for s in strategies: - idx = np.argmin(np.abs(years - s['years'])) - ax.plot(years[idx], energy[idx], 'o', color=s['color'], markersize=15, - markeredgecolor='black', markeredgewidth=2) - ax.annotate(s['name'], (years[idx], energy[idx]), - textcoords="offset points", xytext=(10, 10), fontsize=11, - fontweight='bold', color=s['color']) - - # Add decision guidance text - textstr = '\n'.join([ - 'Decision Guidance:', - '─────────────────', - 'λ < 150 PJ/year → Strategy A', - ' (Long-term cost efficiency)', - '', - '150 ≤ λ ≤ 250 → Strategy C', - ' (Balanced trade-off)', - '', - 'λ > 250 PJ/year → Strategy B', - ' (Time-critical mission)', - ]) - props = dict(boxstyle='round', facecolor='wheat', alpha=0.8) - ax.text(0.02, 0.98, textstr, transform=ax.transAxes, fontsize=10, - verticalalignment='top', bbox=props, family='monospace') - - ax.set_xlabel('Construction Timeline T (years)', fontsize=13) - ax.set_ylabel('Total Energy E (PJ)', fontsize=13) - ax.set_title('Decision Map: Optimal Strategy Selection Based on Time Opportunity Cost λ', fontsize=14) - ax.legend(loc='upper right', fontsize=10) - ax.grid(True, alpha=0.3) - ax.set_xlim(95, 200) - ax.set_ylim(10000, 65000) - - plt.tight_layout() - plt.savefig(save_path, dpi=150, bbox_inches='tight') - print(f"Decision map saved to: {save_path}") - - return fig - - -# ============== Critical Point Analysis ============== - -def analyze_curve_structure(df: pd.DataFrame) -> Dict: - """ - Analyze the structure of the trade-off curve to understand - why certain λ values lead to specific optimal points. - """ - years = df['years'].values - energy = df['energy_PJ'].values - - # Calculate marginal rates - marginal = -np.gradient(energy, years) - - # Find the elevator-only boundary - T_elev = TOTAL_PAYLOAD / TOTAL_ELEVATOR_CAPACITY - idx_elev = np.argmin(np.abs(years - T_elev)) - - # Analyze marginal rate distribution - # Region 1: T > T_elev (pure elevator, marginal ≈ 0) - # Region 2: T < T_elev (need rockets, marginal > 0) - - # Find critical points where marginal rate changes significantly - marginal_near_elev = marginal[idx_elev - 5:idx_elev + 5] - marginal_at_139 = marginal[np.argmin(np.abs(years - 139))] - marginal_at_101 = marginal[np.argmin(np.abs(years - 101))] - - return { - 'T_elev': T_elev, - 'idx_elev': idx_elev, - 'marginal_at_elev': marginal[idx_elev], - 'marginal_at_139': marginal_at_139, - 'marginal_at_101': marginal_at_101, - 'energy_at_elev': energy[idx_elev], - 'energy_at_139': energy[np.argmin(np.abs(years - 139))], - 'energy_at_101': energy[np.argmin(np.abs(years - 101))], - } - - -def find_critical_lambda(df: pd.DataFrame) -> Dict: - """ - Find the critical λ values where optimal point transitions occur. - """ - years = df['years'].values - energy = df['energy_PJ'].values - - # Dense lambda scan - lambda_range = np.linspace(1, 1000, 2000) - - transitions = [] - prev_years = None - - for lam in lambda_range: - opt = find_optimal_point(df, lam) - if prev_years is not None and abs(opt['years'] - prev_years) > 5: - transitions.append({ - 'lambda': lam, - 'from_years': prev_years, - 'to_years': opt['years'], - }) - prev_years = opt['years'] - - return transitions - - -def plot_comprehensive_analysis( - df: pd.DataFrame, - save_path: str = '/Volumes/Files/code/mm/20260130_b/p1/lambda_comprehensive.png' -): - """ - Comprehensive plot showing: - 1. Trade-off curve with marginal rates - 2. Critical λ transitions - 3. Decision regions - """ - fig = plt.figure(figsize=(16, 14)) - - years = df['years'].values - energy = df['energy_PJ'].values - marginal = calculate_marginal_energy_saving(df) - - T_elev = TOTAL_PAYLOAD / TOTAL_ELEVATOR_CAPACITY - T_min = years.min() - - # ========== Plot 1: Trade-off curve with annotations ========== - ax1 = fig.add_subplot(2, 2, 1) - - ax1.plot(years, energy / 1000, 'b-', linewidth=2.5, label='Trade-off Curve') - - # Mark key points - key_points = [ - (T_min, 'Minimum Time\n(~101 years)', 'red'), - (139, 'Original Knee\n(139 years)', 'orange'), - (T_elev, 'Elevator-Only\n(~186 years)', 'green'), - ] - - for t, label, color in key_points: - idx = np.argmin(np.abs(years - t)) - ax1.plot(years[idx], energy[idx] / 1000, 'o', color=color, markersize=12, - markeredgecolor='black', markeredgewidth=2, zorder=5) - ax1.annotate(label, (years[idx], energy[idx] / 1000), - textcoords="offset points", xytext=(10, 10), fontsize=10, - color=color, fontweight='bold') - - ax1.set_xlabel('Construction Timeline T (years)', fontsize=12) - ax1.set_ylabel('Total Energy E (×10³ PJ)', fontsize=12) - ax1.set_title('Energy-Time Trade-off Curve', fontsize=14) - ax1.grid(True, alpha=0.3) - ax1.set_xlim(95, 200) - ax1.legend(loc='upper right') - - # ========== Plot 2: Marginal Rate Analysis ========== - ax2 = fig.add_subplot(2, 2, 2) - - ax2.plot(years, marginal, 'r-', linewidth=2, label='Marginal Saving Rate') - ax2.fill_between(years, 0, marginal, alpha=0.3, color='red') - - # Mark critical thresholds - ax2.axhline(y=480, color='purple', linestyle='--', linewidth=2, - label='Critical λ ≈ 480 PJ/year') - ax2.axvline(x=T_elev, color='green', linestyle=':', alpha=0.7) - - # Annotate the jump at elevator boundary - ax2.annotate('Marginal rate jumps\nat elevator capacity limit', - xy=(T_elev, marginal[np.argmin(np.abs(years - T_elev))]), - xytext=(150, 400), fontsize=10, - arrowprops=dict(arrowstyle='->', color='black')) - - ax2.set_xlabel('Construction Timeline T (years)', fontsize=12) - ax2.set_ylabel('Marginal Energy Saving -dE/dT (PJ/year)', fontsize=12) - ax2.set_title('Marginal Analysis: Why 186 years is optimal for most λ', fontsize=14) - ax2.grid(True, alpha=0.3) - ax2.set_xlim(95, 200) - ax2.set_ylim(-50, 600) - ax2.legend(loc='upper right') - - # ========== Plot 3: λ Sensitivity with Phase Transitions (400-600 focus) ========== - ax3 = fig.add_subplot(2, 2, 3) - - lambda_range = np.linspace(400, 600, 200) - opt_years = [] - for lam in lambda_range: - opt = find_optimal_point(df, lam) - opt_years.append(opt['years']) - - ax3.plot(lambda_range, opt_years, 'b-', linewidth=2.5) - - # Shade phases - ax3.axhspan(180, 190, alpha=0.3, color='green', label='Phase 1: Elevator-Only (λ < 480)') - ax3.axhspan(100, 145, alpha=0.3, color='red', label='Phase 2: Hybrid (λ > 480)') - - # Mark critical λ - ax3.axvline(x=480, color='purple', linestyle='--', linewidth=2) - ax3.annotate('Critical Transition\nλ ≈ 480 PJ/year', - xy=(480, 160), fontsize=11, color='purple', fontweight='bold', - ha='center') - - # Mark 139 year point - ax3.axhline(y=139, color='orange', linestyle=':', linewidth=2, alpha=0.8) - ax3.scatter([500], [139], s=200, c='orange', marker='*', zorder=5, - edgecolors='black', linewidths=2, label='T*=139y at λ≈500') - - ax3.set_xlabel('Time Opportunity Cost λ (PJ/year)', fontsize=12) - ax3.set_ylabel('Optimal Timeline T* (years)', fontsize=12) - ax3.set_title('Phase Transition in Optimal Strategy (λ=400-600)', fontsize=14) - ax3.grid(True, alpha=0.3) - ax3.set_xlim(400, 600) - ax3.set_ylim(95, 195) - ax3.legend(loc='upper right') - - # ========== Plot 4: Decision Summary Table ========== - ax4 = fig.add_subplot(2, 2, 4) - ax4.axis('off') - - # Create summary table - summary_text = """ - ┌─────────────────────────────────────────────────────────────────┐ - │ KEY FINDINGS FROM λ COST ANALYSIS │ - ├─────────────────────────────────────────────────────────────────┤ - │ │ - │ 1. CURVE STRUCTURE │ - │ • Sharp discontinuity at T = 186 years (elevator capacity) │ - │ • Marginal rate jumps from ~0 to ~480 PJ/year at boundary │ - │ │ - │ 2. OPTIMAL POINT SELECTION │ - │ • λ < 480 PJ/year → T* = 186 years (elevator-only) │ - │ • λ ≈ 500 PJ/year → T* = 139 years (original knee) │ - │ • λ > 600 PJ/year → T* → 101 years (time-priority) │ - │ │ - │ 3. IMPLICATIONS FOR THE 139-YEAR KNEE POINT │ - │ • Requires implicit assumption: λ ≈ 500 PJ/year │ - │ • This means: 1 year delay costs ~500 PJ opportunity cost │ - │ • Equivalent to: ~0.5% of total elevator energy per year │ - │ │ - │ 4. RECOMMENDATION │ - │ • Either justify λ ≈ 500 PJ/year with physical reasoning │ - │ • Or acknowledge 186 years as cost-optimal baseline │ - │ • Present 139 years as a time-constrained scenario │ - │ │ - └─────────────────────────────────────────────────────────────────┘ - """ - - ax4.text(0.05, 0.95, summary_text, transform=ax4.transAxes, - fontsize=11, verticalalignment='top', family='monospace', - bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8)) - - plt.tight_layout() - plt.savefig(save_path, dpi=150, bbox_inches='tight') - print(f"Comprehensive analysis saved to: {save_path}") - - return fig - - -# ============== Main Analysis ============== - -def main(): - print("=" * 70) - print("Lambda Time Cost Analysis for Moon Colony Logistics") - print("=" * 70) - - # Generate trade-off curve - print("\n[1] Generating Energy-Time trade-off curve...") - df = generate_tradeoff_curve(year_min=100, year_max=220, num_points=500) - print(f" Generated {len(df)} data points") - - # Key statistics - T_min = df['years'].min() - T_max = df['years'].max() - E_min = df['energy_PJ'].min() - E_max = df['energy_PJ'].max() - T_elev = TOTAL_PAYLOAD / TOTAL_ELEVATOR_CAPACITY - - print(f"\n[2] Trade-off Curve Statistics:") - print(f" Timeline range: {T_min:.1f} - {T_max:.1f} years") - print(f" Energy range: {E_min:.0f} - {E_max:.0f} PJ") - print(f" Elevator-only timeline: {T_elev:.1f} years") - - # Curve structure analysis - print("\n[3] Curve Structure Analysis:") - structure = analyze_curve_structure(df) - print(f" At elevator boundary (T={structure['T_elev']:.1f}y):") - print(f" - Energy: {structure['energy_at_elev']:.0f} PJ") - print(f" - Marginal rate: {structure['marginal_at_elev']:.1f} PJ/year") - print(f" At T=139 years:") - print(f" - Energy: {structure['energy_at_139']:.0f} PJ") - print(f" - Marginal rate: {structure['marginal_at_139']:.1f} PJ/year") - print(f" At T=101 years:") - print(f" - Energy: {structure['energy_at_101']:.0f} PJ") - print(f" - Marginal rate: {structure['marginal_at_101']:.1f} PJ/year") - - # Find critical transitions - print("\n[4] Critical λ Transitions:") - transitions = find_critical_lambda(df) - for t in transitions: - print(f" λ ≈ {t['lambda']:.0f} PJ/year: T* jumps from {t['from_years']:.0f}y to {t['to_years']:.0f}y") - - # Sensitivity analysis - print("\n[5] Sensitivity Analysis Results:") - print(" " + "-" * 55) - print(f" {'λ (PJ/year)':<15} {'Optimal T (years)':<20} {'Energy (PJ)':<15}") - print(" " + "-" * 55) - - test_lambdas = [100, 200, 300, 400, 450, 480, 500, 550, 600] - sensitivity_results = [] - - for lam in test_lambdas: - opt = find_optimal_point(df, lam) - print(f" {lam:<15.0f} {opt['years']:<20.1f} {opt['energy_PJ']:<15.0f}") - sensitivity_results.append(opt) - - # Find λ that gives 139 years - print("\n[6] Original Knee Point Analysis (139 years):") - lambda_for_139 = None - for lam in np.linspace(400, 600, 1000): - opt = find_optimal_point(df, lam) - if abs(opt['years'] - 139) < 1: - lambda_for_139 = lam - print(f" λ ≈ {lam:.0f} PJ/year corresponds to T* ≈ 139 years") - break - - if lambda_for_139: - print(f"\n INTERPRETATION:") - print(f" To justify 139-year knee point, one must argue that:") - print(f" • Each year of delay costs ~{lambda_for_139:.0f} PJ in opportunity") - print(f" • This is {lambda_for_139/E_min*100:.1f}% of minimum total energy per year") - print(f" • Over 47 years (139→186), this amounts to {lambda_for_139*47:.0f} PJ") - - # Generate plots - print("\n[7] Generating visualization plots...") - plot_lambda_analysis(df) - plot_decision_recommendation(df) - plot_comprehensive_analysis(df) - - # Save sensitivity data - sensitivity_df = sensitivity_analysis(df, np.linspace(50, 600, 120)) - sensitivity_df.to_csv('/Volumes/Files/code/mm/20260130_b/p1/lambda_sensitivity.csv', index=False) - print(" Data saved to: lambda_sensitivity.csv") - - # Paper modification recommendations - print("\n" + "=" * 70) - print("RECOMMENDATIONS FOR PAPER MODIFICATION") - print("=" * 70) - print(""" - 1. REFRAME THE OPTIMIZATION PROBLEM - Current: "Multi-objective Pareto optimization" - Suggested: "Constrained optimization with time-cost trade-off" - - 2. INTRODUCE λ AS DECISION PARAMETER - Add equation: J = E_total + λ × T - where λ represents "time opportunity cost" (PJ/year) - - 3. JUSTIFY THE KNEE POINT SELECTION - Option A: Argue λ ≈ 480-500 PJ/year is reasonable because: - - Delayed lunar resource extraction - - Extended Earth-side operational costs - - Strategic/geopolitical value of early completion - - Option B: Present multiple scenarios: - - "Cost-optimal" (λ < 480): T* = 186 years - - "Balanced" (λ ≈ 500): T* = 139 years - - "Time-critical" (λ > 600): T* → 101 years - - 4. ADD SENSITIVITY ANALYSIS FIGURE - Show how optimal T* changes with λ (Figure generated) - - 5. ACKNOWLEDGE THE DISCONTINUITY - Note that the trade-off curve has a sharp transition at - T = 186 years due to elevator capacity constraints - """) - - print("=" * 70) - print("Analysis Complete!") - print("=" * 70) - - return df, sensitivity_results - - -if __name__ == "__main__": - df, results = main() +""" +Lambda Time Cost Analysis for Moon Colony Logistics + +This module introduces the time opportunity cost (λ) to find the optimal +operating point on the Energy-Time trade-off curve. + +Objective Function: J = E_total + λ × T +where: + - E_total: Total energy consumption (PJ) + - T: Construction timeline (years) + - λ: Time opportunity cost (PJ/year) + +The optimal point satisfies: dE/dT = -λ (marginal energy saving = time cost) +""" + +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 + +# 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 + +# Mission parameters +TOTAL_PAYLOAD = 100e6 # 100 million metric tons + +# Space Elevator parameters +NUM_ELEVATORS = 3 +ELEVATOR_CAPACITY_PER_YEAR = 179000 # metric tons per elevator per year +TOTAL_ELEVATOR_CAPACITY = NUM_ELEVATORS * ELEVATOR_CAPACITY_PER_YEAR # 537,000 tons/year +ELEVATOR_SPECIFIC_ENERGY = 157.2e9 # J per metric ton (157.2 MJ/kg × 1000) + +# Rocket parameters - LOX/CH4 (Raptor-class) +PAYLOAD_PER_LAUNCH = 125 # metric tons per launch +ISP = 363 # Specific impulse (s) +SPECIFIC_FUEL_ENERGY = 11.9e6 # J/kg +ALPHA = 0.10 # Structural coefficient +NUM_STAGES = 3 +DELTA_V_BASE = 13300 # m/s (LEO + TLI + LOI) + + +# ============== Launch Site Definition ============== +@dataclass +class LaunchSite: + name: str + short_name: str + latitude: float + max_launches_per_day: int = 1 + + @property + def abs_latitude(self) -> float: + return abs(self.latitude) + + @property + def delta_v_loss(self) -> float: + v_equator = OMEGA_EARTH * R_EARTH + v_site = OMEGA_EARTH * R_EARTH * np.cos(np.radians(self.abs_latitude)) + return v_equator - v_site + + @property + def total_delta_v(self) -> float: + return DELTA_V_BASE + self.delta_v_loss + + +LAUNCH_SITES = sorted([ + LaunchSite("Kourou (French Guiana)", "Kourou", 5.2), + LaunchSite("Satish Dhawan (India)", "SDSC", 13.7), + LaunchSite("Boca Chica (Texas)", "Texas", 26.0), + LaunchSite("Cape Canaveral (Florida)", "Florida", 28.5), + LaunchSite("Vandenberg (California)", "California", 34.7), + LaunchSite("Wallops (Virginia)", "Virginia", 37.8), + LaunchSite("Taiyuan (China)", "Taiyuan", 38.8), + LaunchSite("Mahia (New Zealand)", "Mahia", 39.3), + LaunchSite("Baikonur (Kazakhstan)", "Baikonur", 45.6), + LaunchSite("Kodiak (Alaska)", "Alaska", 57.4), +], key=lambda x: x.abs_latitude) + + +# ============== Core Calculation Functions ============== + +def fuel_ratio_multistage(delta_v: float) -> float: + """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 rocket_energy_per_ton(site: LaunchSite) -> float: + """Energy consumption per ton of payload for rocket launch (J/ton)""" + k = fuel_ratio_multistage(site.total_delta_v) + fuel_per_ton = k * 1000 # kg fuel per metric ton payload + return fuel_per_ton * SPECIFIC_FUEL_ENERGY + + +def calculate_scenario(completion_years: float) -> Optional[Dict]: + """ + Calculate optimal scenario for given completion timeline + (elevator priority + low-latitude rockets) + """ + # Space elevator transport + elevator_payload = min(TOTAL_ELEVATOR_CAPACITY * completion_years, TOTAL_PAYLOAD) + elevator_energy = elevator_payload * 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, + 'elevator_energy_PJ': elevator_energy / 1e15, + 'rocket_energy_PJ': 0, + 'total_energy_PJ': elevator_energy / 1e15, + 'rocket_launches': 0, + 'sites_used': 0, + 'elevator_fraction': 1.0, + } + + # Rocket launches needed + rocket_launches_needed = int(np.ceil(remaining_payload / PAYLOAD_PER_LAUNCH)) + + # Allocate by latitude priority + days_available = completion_years * 365 + max_launches_per_site = int(days_available) + + # Check feasibility + total_rocket_capacity = len(LAUNCH_SITES) * max_launches_per_site * PAYLOAD_PER_LAUNCH + if remaining_payload > total_rocket_capacity: + return None + + rocket_energy = 0 + sites_used = 0 + remaining_launches = rocket_launches_needed + + for site in LAUNCH_SITES: + if remaining_launches <= 0: + break + allocated = min(remaining_launches, max_launches_per_site) + rocket_energy += rocket_energy_per_ton(site) * PAYLOAD_PER_LAUNCH * allocated + remaining_launches -= allocated + if allocated > 0: + sites_used += 1 + + rocket_payload = rocket_launches_needed * PAYLOAD_PER_LAUNCH + total_energy = elevator_energy + rocket_energy + + return { + 'years': completion_years, + 'elevator_payload': elevator_payload, + 'rocket_payload': rocket_payload, + 'elevator_energy_PJ': elevator_energy / 1e15, + 'rocket_energy_PJ': rocket_energy / 1e15, + 'total_energy_PJ': total_energy / 1e15, + 'rocket_launches': rocket_launches_needed, + 'sites_used': sites_used, + 'elevator_fraction': elevator_payload / TOTAL_PAYLOAD, + } + + +# ============== Generate Trade-off Curve ============== + +def generate_tradeoff_curve( + year_min: float = 100, + year_max: float = 250, + num_points: int = 500 +) -> pd.DataFrame: + """Generate Energy-Time trade-off curve data""" + years_range = np.linspace(year_min, year_max, num_points) + + results = [] + for years in years_range: + scenario = calculate_scenario(years) + if scenario is not None: + results.append({ + 'years': years, + 'energy_PJ': scenario['total_energy_PJ'], + 'elevator_fraction': scenario['elevator_fraction'], + 'sites_used': scenario['sites_used'], + 'rocket_launches': scenario['rocket_launches'], + }) + + return pd.DataFrame(results) + + +# ============== Lambda Cost Analysis ============== + +def calculate_total_cost(df: pd.DataFrame, lambda_cost: float) -> np.ndarray: + """ + Calculate total cost J = E + λ × T + + Args: + df: Trade-off curve data + lambda_cost: Time opportunity cost (PJ/year) + + Returns: + Total cost array + """ + return df['energy_PJ'].values + lambda_cost * df['years'].values + + +def find_optimal_point(df: pd.DataFrame, lambda_cost: float) -> Dict: + """ + Find optimal point that minimizes J = E + λ × T + + Args: + df: Trade-off curve data + lambda_cost: Time opportunity cost (PJ/year) + + Returns: + Optimal point information + """ + total_cost = calculate_total_cost(df, lambda_cost) + opt_idx = np.argmin(total_cost) + + return { + 'index': opt_idx, + 'years': df['years'].iloc[opt_idx], + 'energy_PJ': df['energy_PJ'].iloc[opt_idx], + 'total_cost': total_cost[opt_idx], + 'elevator_fraction': df['elevator_fraction'].iloc[opt_idx], + 'lambda': lambda_cost, + } + + +def calculate_marginal_energy_saving(df: pd.DataFrame) -> np.ndarray: + """ + Calculate marginal energy saving rate: -dE/dT (PJ/year) + + This represents how much energy is saved per additional year of timeline. + """ + years = df['years'].values + energy = df['energy_PJ'].values + + # Use central difference for interior points + marginal = -np.gradient(energy, years) + + return marginal + + +def sensitivity_analysis( + df: pd.DataFrame, + lambda_range: np.ndarray +) -> pd.DataFrame: + """ + Perform sensitivity analysis on λ parameter + + Args: + df: Trade-off curve data + lambda_range: Array of λ values to test + + Returns: + DataFrame with optimal points for each λ + """ + results = [] + for lam in lambda_range: + opt = find_optimal_point(df, lam) + results.append({ + 'lambda_PJ_per_year': lam, + 'optimal_years': opt['years'], + 'optimal_energy_PJ': opt['energy_PJ'], + 'total_cost_PJ': opt['total_cost'], + 'elevator_fraction': opt['elevator_fraction'], + }) + + return pd.DataFrame(results) + + +# ============== Visualization Functions ============== + +def plot_lambda_analysis( + df: pd.DataFrame, + save_path: str = '/Volumes/Files/code/mm/20260130_b/p1/lambda_cost_analysis.png' +): + """ + Comprehensive visualization of λ time cost analysis + Focus on critical range λ = 400-600 PJ/year + """ + fig, axes = plt.subplots(2, 2, figsize=(14, 12)) + + years = df['years'].values + energy = df['energy_PJ'].values + + # Key boundaries + T_min = years.min() # ~100.7 years (fastest) + T_elev = TOTAL_PAYLOAD / TOTAL_ELEVATOR_CAPACITY # ~186 years (elevator only) + E_min = energy[years >= T_elev].min() if any(years >= T_elev) else energy.min() + + # ========== Plot 1: Trade-off Curve with λ Iso-cost Lines ========== + ax1 = axes[0, 0] + + # Main curve + ax1.plot(years, energy, 'b-', linewidth=2.5, label='Energy-Time Trade-off') + + # Mark key points + ax1.axvline(x=T_elev, color='green', linestyle='--', alpha=0.7, label=f'Elevator-only: {T_elev:.1f} years') + ax1.axvline(x=T_min, color='red', linestyle='--', alpha=0.7, label=f'Minimum time: {T_min:.1f} years') + + # Draw iso-cost lines for different λ (focus on 400-600 range) + lambda_values = [420, 480, 500, 550] + colors = ['#2ecc71', '#e74c3c', '#9b59b6', '#3498db'] + + for lam, color in zip(lambda_values, colors): + opt = find_optimal_point(df, lam) + # Iso-cost line: E + λT = const → E = const - λT + T_line = np.linspace(80, 220, 100) + E_line = opt['total_cost'] - lam * T_line + valid = (E_line > 0) & (E_line < 70000) + ax1.plot(T_line[valid], E_line[valid], '--', color=color, alpha=0.6, linewidth=1.5) + ax1.plot(opt['years'], opt['energy_PJ'], 'o', color=color, markersize=12, + markeredgecolor='black', markeredgewidth=1.5, + label=f'λ={lam}: T={opt["years"]:.1f}y, E={opt["energy_PJ"]:.0f}PJ') + + ax1.set_xlabel('Construction Timeline T (years)', fontsize=12) + ax1.set_ylabel('Total Energy E (PJ)', fontsize=12) + ax1.set_title('Energy-Time Trade-off with λ Iso-cost Lines (λ=400-600)\n$J = E + λT$', fontsize=13) + ax1.legend(loc='upper right', fontsize=9) + ax1.grid(True, alpha=0.3) + ax1.set_xlim(95, 200) + ax1.set_ylim(10000, 65000) + + # ========== Plot 2: Marginal Energy Saving Rate ========== + ax2 = axes[0, 1] + + marginal = calculate_marginal_energy_saving(df) + + ax2.plot(years, marginal, 'r-', linewidth=2, label='Marginal Energy Saving -dE/dT') + ax2.axhline(y=0, color='black', linestyle='-', alpha=0.3) + + # Mark critical λ values in 400-600 range + for lam, color in zip([450, 480, 500, 550], ['#2ecc71', '#e74c3c', '#9b59b6', '#3498db']): + ax2.axhline(y=lam, color=color, linestyle='--', alpha=0.7, label=f'λ = {lam} PJ/year') + # Find intersection + intersect_idx = np.argmin(np.abs(marginal - lam)) + ax2.plot(years[intersect_idx], marginal[intersect_idx], 'o', color=color, markersize=8) + + ax2.set_xlabel('Construction Timeline T (years)', fontsize=12) + ax2.set_ylabel('Marginal Energy Saving -dE/dT (PJ/year)', fontsize=12) + ax2.set_title('Marginal Analysis: Optimal Condition is -dE/dT = λ', fontsize=13) + ax2.legend(loc='upper right', fontsize=9) + ax2.grid(True, alpha=0.3) + ax2.set_xlim(95, 200) + ax2.set_ylim(350, 620) + + # ========== Plot 3: Sensitivity Analysis (400-600 range) ========== + ax3 = axes[1, 0] + + lambda_range = np.linspace(400, 600, 200) + sensitivity_df = sensitivity_analysis(df, lambda_range) + + # Plot optimal years vs λ + ax3.plot(sensitivity_df['lambda_PJ_per_year'], sensitivity_df['optimal_years'], + 'b-', linewidth=2.5, label='Optimal Timeline T*') + + # Mark key regions + ax3.axhline(y=T_elev, color='green', linestyle='--', alpha=0.7, label=f'Elevator-only: {T_elev:.1f}y') + ax3.axhline(y=T_min, color='red', linestyle='--', alpha=0.7, label=f'Min timeline: {T_min:.1f}y') + ax3.axhline(y=139, color='orange', linestyle=':', linewidth=2, alpha=0.8, label='Original knee: 139y') + + # Find λ corresponding to 139 years + idx_139 = np.argmin(np.abs(sensitivity_df['optimal_years'] - 139)) + if idx_139 < len(sensitivity_df): + lambda_139 = sensitivity_df['lambda_PJ_per_year'].iloc[idx_139] + ax3.axvline(x=lambda_139, color='orange', linestyle=':', linewidth=2, alpha=0.6) + ax3.scatter([lambda_139], [139], s=150, c='orange', marker='*', zorder=5, + edgecolors='black', linewidths=1.5, label=f'λ≈{lambda_139:.0f} for T*=139y') + + # Mark critical transition + ax3.axvline(x=480, color='purple', linestyle='--', alpha=0.5) + ax3.text(482, 175, 'Critical\nTransition', fontsize=9, color='purple') + + ax3.set_xlabel('Time Opportunity Cost λ (PJ/year)', fontsize=12) + ax3.set_ylabel('Optimal Timeline T* (years)', fontsize=12) + ax3.set_title('Sensitivity Analysis: Optimal Timeline vs λ (400-600 range)', fontsize=13) + ax3.legend(loc='upper right', fontsize=9) + ax3.grid(True, alpha=0.3) + ax3.set_xlim(400, 600) + ax3.set_ylim(95, 195) + + # ========== Plot 4: Total Cost Curves (400-600 range) ========== + ax4 = axes[1, 1] + + for lam, color in [(450, '#2ecc71'), (480, '#e74c3c'), (500, '#9b59b6'), (550, '#3498db')]: + total_cost = calculate_total_cost(df, lam) + ax4.plot(years, total_cost / 1000, color=color, linestyle='-', + linewidth=2, label=f'λ={lam} PJ/year') + + # Mark minimum + opt = find_optimal_point(df, lam) + ax4.plot(opt['years'], opt['total_cost'] / 1000, 'o', color=color, markersize=10, + markeredgecolor='black', markeredgewidth=1.5) + + ax4.set_xlabel('Construction Timeline T (years)', fontsize=12) + ax4.set_ylabel('Total Cost J = E + λT (×10³ PJ)', fontsize=12) + ax4.set_title('Total Cost Function for λ = 400-600 PJ/year', fontsize=13) + ax4.legend(loc='upper right', fontsize=10) + 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"Lambda analysis plot saved to: {save_path}") + + return fig + + +def plot_decision_recommendation( + df: pd.DataFrame, + save_path: str = '/Volumes/Files/code/mm/20260130_b/p1/lambda_decision_map.png' +): + """ + Decision map showing optimal choice based on λ preference + """ + fig, ax = plt.subplots(figsize=(12, 8)) + + years = df['years'].values + energy = df['energy_PJ'].values + + # Main trade-off curve + ax.plot(years, energy, 'b-', linewidth=3, label='Feasible Trade-off Curve') + + # Key boundaries + T_elev = TOTAL_PAYLOAD / TOTAL_ELEVATOR_CAPACITY + T_min = years.min() + + # Shade decision regions + # Region 1: Low λ (cost priority) → longer timeline + ax.axvspan(160, 190, alpha=0.2, color='green', label='Low λ (<150): Cost Priority') + # Region 2: Medium λ (balanced) → middle ground + ax.axvspan(130, 160, alpha=0.2, color='yellow', label='Medium λ (150-250): Balanced') + # Region 3: High λ (time priority) → shorter timeline + ax.axvspan(100, 130, alpha=0.2, color='red', label='High λ (>250): Time Priority') + + # Mark specific strategy points + strategies = [ + {'name': 'A: Cost-Prioritized', 'years': 186, 'lambda': 0, 'color': 'green'}, + {'name': 'C: Balanced (λ≈200)', 'years': 139, 'lambda': 200, 'color': 'orange'}, + {'name': 'B: Time-Prioritized', 'years': 101, 'lambda': 500, 'color': 'red'}, + ] + + for s in strategies: + idx = np.argmin(np.abs(years - s['years'])) + ax.plot(years[idx], energy[idx], 'o', color=s['color'], markersize=15, + markeredgecolor='black', markeredgewidth=2) + ax.annotate(s['name'], (years[idx], energy[idx]), + textcoords="offset points", xytext=(10, 10), fontsize=11, + fontweight='bold', color=s['color']) + + # Add decision guidance text + textstr = '\n'.join([ + 'Decision Guidance:', + '─────────────────', + 'λ < 150 PJ/year → Strategy A', + ' (Long-term cost efficiency)', + '', + '150 ≤ λ ≤ 250 → Strategy C', + ' (Balanced trade-off)', + '', + 'λ > 250 PJ/year → Strategy B', + ' (Time-critical mission)', + ]) + props = dict(boxstyle='round', facecolor='wheat', alpha=0.8) + ax.text(0.02, 0.98, textstr, transform=ax.transAxes, fontsize=10, + verticalalignment='top', bbox=props, family='monospace') + + ax.set_xlabel('Construction Timeline T (years)', fontsize=13) + ax.set_ylabel('Total Energy E (PJ)', fontsize=13) + ax.set_title('Decision Map: Optimal Strategy Selection Based on Time Opportunity Cost λ', fontsize=14) + ax.legend(loc='upper right', fontsize=10) + ax.grid(True, alpha=0.3) + ax.set_xlim(95, 200) + ax.set_ylim(10000, 65000) + + plt.tight_layout() + plt.savefig(save_path, dpi=150, bbox_inches='tight') + print(f"Decision map saved to: {save_path}") + + return fig + + +# ============== Critical Point Analysis ============== + +def analyze_curve_structure(df: pd.DataFrame) -> Dict: + """ + Analyze the structure of the trade-off curve to understand + why certain λ values lead to specific optimal points. + """ + years = df['years'].values + energy = df['energy_PJ'].values + + # Calculate marginal rates + marginal = -np.gradient(energy, years) + + # Find the elevator-only boundary + T_elev = TOTAL_PAYLOAD / TOTAL_ELEVATOR_CAPACITY + idx_elev = np.argmin(np.abs(years - T_elev)) + + # Analyze marginal rate distribution + # Region 1: T > T_elev (pure elevator, marginal ≈ 0) + # Region 2: T < T_elev (need rockets, marginal > 0) + + # Find critical points where marginal rate changes significantly + marginal_near_elev = marginal[idx_elev - 5:idx_elev + 5] + marginal_at_139 = marginal[np.argmin(np.abs(years - 139))] + marginal_at_101 = marginal[np.argmin(np.abs(years - 101))] + + return { + 'T_elev': T_elev, + 'idx_elev': idx_elev, + 'marginal_at_elev': marginal[idx_elev], + 'marginal_at_139': marginal_at_139, + 'marginal_at_101': marginal_at_101, + 'energy_at_elev': energy[idx_elev], + 'energy_at_139': energy[np.argmin(np.abs(years - 139))], + 'energy_at_101': energy[np.argmin(np.abs(years - 101))], + } + + +def find_critical_lambda(df: pd.DataFrame) -> Dict: + """ + Find the critical λ values where optimal point transitions occur. + """ + years = df['years'].values + energy = df['energy_PJ'].values + + # Dense lambda scan + lambda_range = np.linspace(1, 1000, 2000) + + transitions = [] + prev_years = None + + for lam in lambda_range: + opt = find_optimal_point(df, lam) + if prev_years is not None and abs(opt['years'] - prev_years) > 5: + transitions.append({ + 'lambda': lam, + 'from_years': prev_years, + 'to_years': opt['years'], + }) + prev_years = opt['years'] + + return transitions + + +def plot_comprehensive_analysis( + df: pd.DataFrame, + save_path: str = '/Volumes/Files/code/mm/20260130_b/p1/lambda_comprehensive.png' +): + """ + Comprehensive plot showing: + 1. Trade-off curve with marginal rates + 2. Critical λ transitions + 3. Decision regions + """ + fig = plt.figure(figsize=(16, 14)) + + years = df['years'].values + energy = df['energy_PJ'].values + marginal = calculate_marginal_energy_saving(df) + + T_elev = TOTAL_PAYLOAD / TOTAL_ELEVATOR_CAPACITY + T_min = years.min() + + # ========== Plot 1: Trade-off curve with annotations ========== + ax1 = fig.add_subplot(2, 2, 1) + + ax1.plot(years, energy / 1000, 'b-', linewidth=2.5, label='Trade-off Curve') + + # Mark key points + key_points = [ + (T_min, 'Minimum Time\n(~101 years)', 'red'), + (139, 'Original Knee\n(139 years)', 'orange'), + (T_elev, 'Elevator-Only\n(~186 years)', 'green'), + ] + + for t, label, color in key_points: + idx = np.argmin(np.abs(years - t)) + ax1.plot(years[idx], energy[idx] / 1000, 'o', color=color, markersize=12, + markeredgecolor='black', markeredgewidth=2, zorder=5) + ax1.annotate(label, (years[idx], energy[idx] / 1000), + textcoords="offset points", xytext=(10, 10), fontsize=10, + color=color, fontweight='bold') + + ax1.set_xlabel('Construction Timeline T (years)', fontsize=12) + ax1.set_ylabel('Total Energy E (×10³ PJ)', fontsize=12) + ax1.set_title('Energy-Time Trade-off Curve', fontsize=14) + ax1.grid(True, alpha=0.3) + ax1.set_xlim(95, 200) + ax1.legend(loc='upper right') + + # ========== Plot 2: Marginal Rate Analysis ========== + ax2 = fig.add_subplot(2, 2, 2) + + ax2.plot(years, marginal, 'r-', linewidth=2, label='Marginal Saving Rate') + ax2.fill_between(years, 0, marginal, alpha=0.3, color='red') + + # Mark critical thresholds + ax2.axhline(y=480, color='purple', linestyle='--', linewidth=2, + label='Critical λ ≈ 480 PJ/year') + ax2.axvline(x=T_elev, color='green', linestyle=':', alpha=0.7) + + # Annotate the jump at elevator boundary + ax2.annotate('Marginal rate jumps\nat elevator capacity limit', + xy=(T_elev, marginal[np.argmin(np.abs(years - T_elev))]), + xytext=(150, 400), fontsize=10, + arrowprops=dict(arrowstyle='->', color='black')) + + ax2.set_xlabel('Construction Timeline T (years)', fontsize=12) + ax2.set_ylabel('Marginal Energy Saving -dE/dT (PJ/year)', fontsize=12) + ax2.set_title('Marginal Analysis: Why 186 years is optimal for most λ', fontsize=14) + ax2.grid(True, alpha=0.3) + ax2.set_xlim(95, 200) + ax2.set_ylim(-50, 600) + ax2.legend(loc='upper right') + + # ========== Plot 3: λ Sensitivity with Phase Transitions (400-600 focus) ========== + ax3 = fig.add_subplot(2, 2, 3) + + lambda_range = np.linspace(400, 600, 200) + opt_years = [] + for lam in lambda_range: + opt = find_optimal_point(df, lam) + opt_years.append(opt['years']) + + ax3.plot(lambda_range, opt_years, 'b-', linewidth=2.5) + + # Shade phases + ax3.axhspan(180, 190, alpha=0.3, color='green', label='Phase 1: Elevator-Only (λ < 480)') + ax3.axhspan(100, 145, alpha=0.3, color='red', label='Phase 2: Hybrid (λ > 480)') + + # Mark critical λ + ax3.axvline(x=480, color='purple', linestyle='--', linewidth=2) + ax3.annotate('Critical Transition\nλ ≈ 480 PJ/year', + xy=(480, 160), fontsize=11, color='purple', fontweight='bold', + ha='center') + + # Mark 139 year point + ax3.axhline(y=139, color='orange', linestyle=':', linewidth=2, alpha=0.8) + ax3.scatter([500], [139], s=200, c='orange', marker='*', zorder=5, + edgecolors='black', linewidths=2, label='T*=139y at λ≈500') + + ax3.set_xlabel('Time Opportunity Cost λ (PJ/year)', fontsize=12) + ax3.set_ylabel('Optimal Timeline T* (years)', fontsize=12) + ax3.set_title('Phase Transition in Optimal Strategy (λ=400-600)', fontsize=14) + ax3.grid(True, alpha=0.3) + ax3.set_xlim(400, 600) + ax3.set_ylim(95, 195) + ax3.legend(loc='upper right') + + # ========== Plot 4: Decision Summary Table ========== + ax4 = fig.add_subplot(2, 2, 4) + ax4.axis('off') + + # Create summary table + summary_text = """ + ┌─────────────────────────────────────────────────────────────────┐ + │ KEY FINDINGS FROM λ COST ANALYSIS │ + ├─────────────────────────────────────────────────────────────────┤ + │ │ + │ 1. CURVE STRUCTURE │ + │ • Sharp discontinuity at T = 186 years (elevator capacity) │ + │ • Marginal rate jumps from ~0 to ~480 PJ/year at boundary │ + │ │ + │ 2. OPTIMAL POINT SELECTION │ + │ • λ < 480 PJ/year → T* = 186 years (elevator-only) │ + │ • λ ≈ 500 PJ/year → T* = 139 years (original knee) │ + │ • λ > 600 PJ/year → T* → 101 years (time-priority) │ + │ │ + │ 3. IMPLICATIONS FOR THE 139-YEAR KNEE POINT │ + │ • Requires implicit assumption: λ ≈ 500 PJ/year │ + │ • This means: 1 year delay costs ~500 PJ opportunity cost │ + │ • Equivalent to: ~0.5% of total elevator energy per year │ + │ │ + │ 4. RECOMMENDATION │ + │ • Either justify λ ≈ 500 PJ/year with physical reasoning │ + │ • Or acknowledge 186 years as cost-optimal baseline │ + │ • Present 139 years as a time-constrained scenario │ + │ │ + └─────────────────────────────────────────────────────────────────┘ + """ + + ax4.text(0.05, 0.95, summary_text, transform=ax4.transAxes, + fontsize=11, verticalalignment='top', family='monospace', + bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8)) + + plt.tight_layout() + plt.savefig(save_path, dpi=150, bbox_inches='tight') + print(f"Comprehensive analysis saved to: {save_path}") + + return fig + + +# ============== Main Analysis ============== + +def main(): + print("=" * 70) + print("Lambda Time Cost Analysis for Moon Colony Logistics") + print("=" * 70) + + # Generate trade-off curve + print("\n[1] Generating Energy-Time trade-off curve...") + df = generate_tradeoff_curve(year_min=100, year_max=220, num_points=500) + print(f" Generated {len(df)} data points") + + # Key statistics + T_min = df['years'].min() + T_max = df['years'].max() + E_min = df['energy_PJ'].min() + E_max = df['energy_PJ'].max() + T_elev = TOTAL_PAYLOAD / TOTAL_ELEVATOR_CAPACITY + + print(f"\n[2] Trade-off Curve Statistics:") + print(f" Timeline range: {T_min:.1f} - {T_max:.1f} years") + print(f" Energy range: {E_min:.0f} - {E_max:.0f} PJ") + print(f" Elevator-only timeline: {T_elev:.1f} years") + + # Curve structure analysis + print("\n[3] Curve Structure Analysis:") + structure = analyze_curve_structure(df) + print(f" At elevator boundary (T={structure['T_elev']:.1f}y):") + print(f" - Energy: {structure['energy_at_elev']:.0f} PJ") + print(f" - Marginal rate: {structure['marginal_at_elev']:.1f} PJ/year") + print(f" At T=139 years:") + print(f" - Energy: {structure['energy_at_139']:.0f} PJ") + print(f" - Marginal rate: {structure['marginal_at_139']:.1f} PJ/year") + print(f" At T=101 years:") + print(f" - Energy: {structure['energy_at_101']:.0f} PJ") + print(f" - Marginal rate: {structure['marginal_at_101']:.1f} PJ/year") + + # Find critical transitions + print("\n[4] Critical λ Transitions:") + transitions = find_critical_lambda(df) + for t in transitions: + print(f" λ ≈ {t['lambda']:.0f} PJ/year: T* jumps from {t['from_years']:.0f}y to {t['to_years']:.0f}y") + + # Sensitivity analysis + print("\n[5] Sensitivity Analysis Results:") + print(" " + "-" * 55) + print(f" {'λ (PJ/year)':<15} {'Optimal T (years)':<20} {'Energy (PJ)':<15}") + print(" " + "-" * 55) + + test_lambdas = [100, 200, 300, 400, 450, 480, 500, 550, 600] + sensitivity_results = [] + + for lam in test_lambdas: + opt = find_optimal_point(df, lam) + print(f" {lam:<15.0f} {opt['years']:<20.1f} {opt['energy_PJ']:<15.0f}") + sensitivity_results.append(opt) + + # Find λ that gives 139 years + print("\n[6] Original Knee Point Analysis (139 years):") + lambda_for_139 = None + for lam in np.linspace(400, 600, 1000): + opt = find_optimal_point(df, lam) + if abs(opt['years'] - 139) < 1: + lambda_for_139 = lam + print(f" λ ≈ {lam:.0f} PJ/year corresponds to T* ≈ 139 years") + break + + if lambda_for_139: + print(f"\n INTERPRETATION:") + print(f" To justify 139-year knee point, one must argue that:") + print(f" • Each year of delay costs ~{lambda_for_139:.0f} PJ in opportunity") + print(f" • This is {lambda_for_139/E_min*100:.1f}% of minimum total energy per year") + print(f" • Over 47 years (139→186), this amounts to {lambda_for_139*47:.0f} PJ") + + # Generate plots + print("\n[7] Generating visualization plots...") + plot_lambda_analysis(df) + plot_decision_recommendation(df) + plot_comprehensive_analysis(df) + + # Save sensitivity data + sensitivity_df = sensitivity_analysis(df, np.linspace(50, 600, 120)) + sensitivity_df.to_csv('/Volumes/Files/code/mm/20260130_b/p1/lambda_sensitivity.csv', index=False) + print(" Data saved to: lambda_sensitivity.csv") + + # Paper modification recommendations + print("\n" + "=" * 70) + print("RECOMMENDATIONS FOR PAPER MODIFICATION") + print("=" * 70) + print(""" + 1. REFRAME THE OPTIMIZATION PROBLEM + Current: "Multi-objective Pareto optimization" + Suggested: "Constrained optimization with time-cost trade-off" + + 2. INTRODUCE λ AS DECISION PARAMETER + Add equation: J = E_total + λ × T + where λ represents "time opportunity cost" (PJ/year) + + 3. JUSTIFY THE KNEE POINT SELECTION + Option A: Argue λ ≈ 480-500 PJ/year is reasonable because: + - Delayed lunar resource extraction + - Extended Earth-side operational costs + - Strategic/geopolitical value of early completion + + Option B: Present multiple scenarios: + - "Cost-optimal" (λ < 480): T* = 186 years + - "Balanced" (λ ≈ 500): T* = 139 years + - "Time-critical" (λ > 600): T* → 101 years + + 4. ADD SENSITIVITY ANALYSIS FIGURE + Show how optimal T* changes with λ (Figure generated) + + 5. ACKNOWLEDGE THE DISCONTINUITY + Note that the trade-off curve has a sharp transition at + T = 186 years due to elevator capacity constraints + """) + + print("=" * 70) + print("Analysis Complete!") + print("=" * 70) + + return df, sensitivity_results + + +if __name__ == "__main__": + df, results = main() diff --git a/p1/richards_curve_1984.py b/p1/richards_curve_1984.py index 5dec4d9..d4c246a 100644 --- a/p1/richards_curve_1984.py +++ b/p1/richards_curve_1984.py @@ -1,183 +1,183 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Richards S-Curve Fit for 1984-2025 Data - -Fits and visualizes the Richards growth model to rocket launch data -starting from 1984. -""" - -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) - - N(t) = K / (1 + exp(-r*(t - t0)))^(1/v) - """ - 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 fit_richards(years, launches, base_year=1984): - """Fit Richards model to data""" - t = (years - base_year).astype(float) - - # Initial parameters and bounds (unconstrained K) - p0 = [5000.0, 0.1, 40.0, 2.0] - bounds = ([500, 0.005, 10, 0.2], [100000, 1.0, 150, 10.0]) - - popt, pcov = curve_fit(richards, t, launches, p0=p0, bounds=bounds, maxfev=100000) - perr = np.sqrt(np.diag(pcov)) - - # Calculate R² - y_pred = richards(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) - - return popt, perr, r_squared - - -def main(): - print("=" * 60) - print("Richards S-Curve Fit (1984-2025)") - print("=" * 60) - - # Load data - df = load_data() - - # Filter 1984-2025 - start_year = 1984 - data = df[(df["year"] >= start_year) & (df["year"] <= 2025)].copy() - years = data["year"].values - launches = data["launches"].values - - print(f"Data range: {start_year} - 2025") - print(f"Data points: {len(data)}") - - # Fit model - popt, perr, r_squared = fit_richards(years, launches, base_year=start_year) - K, r, t0, v = popt - - print(f"\nFitted Parameters:") - print(f" K (carrying capacity) = {K:.0f} launches/year") - 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}") - - # Physical limit - physical_max = 3650 - print(f"\nPhysical limit: {physical_max} (10 sites × 365 days)") - print(f"K / Physical limit = {K/physical_max:.2f}x") - - # ========== Create Visualization ========== - fig, ax = plt.subplots(figsize=(12, 7)) - - # Historical data points - ax.scatter(years, launches, color='#2C3E50', s=80, alpha=0.8, - label='Historical Data (1984-2025)', zorder=3, edgecolor='white', linewidth=0.5) - - # Generate smooth S-curve - years_smooth = np.linspace(start_year, 2100, 500) - t_smooth = years_smooth - start_year - pred_smooth = richards(t_smooth, *popt) - - # S-curve prediction - ax.plot(years_smooth, pred_smooth, color='#27AE60', lw=3, - label=f'Richards Model (K={K:.0f}, R²={r_squared:.3f})', zorder=2) - - # Confidence band (approximate using parameter errors) - K_low = max(500, K - 2*perr[0]) - K_high = K + 2*perr[0] - pred_low = richards(t_smooth, K_low, r, t0, v) - pred_high = richards(t_smooth, K_high, r, t0, v) - ax.fill_between(years_smooth, pred_low, pred_high, color='#27AE60', alpha=0.15, - label='95% Confidence Band') - - # Physical limit line - ax.axhline(physical_max, color='#E74C3C', ls='--', lw=2.5, - label=f'Physical Limit: {physical_max} (1/site/day)') - - # Mark inflection point - inflection_year = start_year + t0 - inflection_value = K / (v + 1) ** (1/v) - ax.scatter([inflection_year], [inflection_value], color='#F39C12', s=150, - marker='*', zorder=5, label=f'Inflection Point ({inflection_year:.0f})') - - # Mark key years - ax.axvline(2025, color='gray', ls=':', lw=1.5, alpha=0.7) - ax.axvline(2050, color='#3498DB', ls=':', lw=2, alpha=0.8) - ax.text(2026, K*0.95, '2025\n(Now)', fontsize=10, color='gray') - ax.text(2051, K*0.85, '2050\n(Target)', fontsize=10, color='#3498DB') - - # Prediction for 2050 - t_2050 = 2050 - start_year - pred_2050 = richards(t_2050, *popt) - ax.scatter([2050], [pred_2050], color='#3498DB', s=100, marker='D', zorder=4) - ax.annotate(f'{pred_2050:.0f}', xy=(2050, pred_2050), xytext=(2055, pred_2050-300), - fontsize=10, color='#3498DB', fontweight='bold') - - # Formatting - ax.set_xlabel('Year', fontsize=12) - ax.set_ylabel('Annual Launches', fontsize=12) - ax.set_title('Richards S-Curve Model Fit (1984-2025 Data)\nRocket Launch Capacity Projection', - fontsize=14, fontweight='bold') - ax.legend(loc='upper left', fontsize=10) - ax.grid(True, alpha=0.3) - ax.set_xlim(1982, 2100) - ax.set_ylim(0, K * 1.15) - - # Add text box with model info - textstr = f'''Model: N(t) = K / (1 + e^(-r(t-t₀)))^(1/v) - -Data Window: {start_year}-2025 ({len(data)} points) -K = {K:.0f} launches/year -r = {r:.4f} -t₀ = {start_year + t0:.1f} -v = {v:.3f} -R² = {r_squared:.4f} - -Note: Physical limit (3,650) shown as -dashed red line''' - - props = dict(boxstyle='round', facecolor='wheat', alpha=0.8) - ax.text(0.98, 0.35, textstr, transform=ax.transAxes, fontsize=9, - verticalalignment='top', horizontalalignment='right', bbox=props, family='monospace') - - plt.tight_layout() - plt.savefig('richards_curve_1984.png', dpi=150, bbox_inches='tight') - plt.close() - - print("\nPlot saved: richards_curve_1984.png") - print("=" * 60) - - -if __name__ == "__main__": - main() +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Richards S-Curve Fit for 1984-2025 Data + +Fits and visualizes the Richards growth model to rocket launch data +starting from 1984. +""" + +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) + + N(t) = K / (1 + exp(-r*(t - t0)))^(1/v) + """ + 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 fit_richards(years, launches, base_year=1984): + """Fit Richards model to data""" + t = (years - base_year).astype(float) + + # Initial parameters and bounds (unconstrained K) + p0 = [5000.0, 0.1, 40.0, 2.0] + bounds = ([500, 0.005, 10, 0.2], [100000, 1.0, 150, 10.0]) + + popt, pcov = curve_fit(richards, t, launches, p0=p0, bounds=bounds, maxfev=100000) + perr = np.sqrt(np.diag(pcov)) + + # Calculate R² + y_pred = richards(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) + + return popt, perr, r_squared + + +def main(): + print("=" * 60) + print("Richards S-Curve Fit (1984-2025)") + print("=" * 60) + + # Load data + df = load_data() + + # Filter 1984-2025 + start_year = 1984 + data = df[(df["year"] >= start_year) & (df["year"] <= 2025)].copy() + years = data["year"].values + launches = data["launches"].values + + print(f"Data range: {start_year} - 2025") + print(f"Data points: {len(data)}") + + # Fit model + popt, perr, r_squared = fit_richards(years, launches, base_year=start_year) + K, r, t0, v = popt + + print(f"\nFitted Parameters:") + print(f" K (carrying capacity) = {K:.0f} launches/year") + 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}") + + # Physical limit + physical_max = 3650 + print(f"\nPhysical limit: {physical_max} (10 sites × 365 days)") + print(f"K / Physical limit = {K/physical_max:.2f}x") + + # ========== Create Visualization ========== + fig, ax = plt.subplots(figsize=(12, 7)) + + # Historical data points + ax.scatter(years, launches, color='#2C3E50', s=80, alpha=0.8, + label='Historical Data (1984-2025)', zorder=3, edgecolor='white', linewidth=0.5) + + # Generate smooth S-curve + years_smooth = np.linspace(start_year, 2100, 500) + t_smooth = years_smooth - start_year + pred_smooth = richards(t_smooth, *popt) + + # S-curve prediction + ax.plot(years_smooth, pred_smooth, color='#27AE60', lw=3, + label=f'Richards Model (K={K:.0f}, R²={r_squared:.3f})', zorder=2) + + # Confidence band (approximate using parameter errors) + K_low = max(500, K - 2*perr[0]) + K_high = K + 2*perr[0] + pred_low = richards(t_smooth, K_low, r, t0, v) + pred_high = richards(t_smooth, K_high, r, t0, v) + ax.fill_between(years_smooth, pred_low, pred_high, color='#27AE60', alpha=0.15, + label='95% Confidence Band') + + # Physical limit line + ax.axhline(physical_max, color='#E74C3C', ls='--', lw=2.5, + label=f'Physical Limit: {physical_max} (1/site/day)') + + # Mark inflection point + inflection_year = start_year + t0 + inflection_value = K / (v + 1) ** (1/v) + ax.scatter([inflection_year], [inflection_value], color='#F39C12', s=150, + marker='*', zorder=5, label=f'Inflection Point ({inflection_year:.0f})') + + # Mark key years + ax.axvline(2025, color='gray', ls=':', lw=1.5, alpha=0.7) + ax.axvline(2050, color='#3498DB', ls=':', lw=2, alpha=0.8) + ax.text(2026, K*0.95, '2025\n(Now)', fontsize=10, color='gray') + ax.text(2051, K*0.85, '2050\n(Target)', fontsize=10, color='#3498DB') + + # Prediction for 2050 + t_2050 = 2050 - start_year + pred_2050 = richards(t_2050, *popt) + ax.scatter([2050], [pred_2050], color='#3498DB', s=100, marker='D', zorder=4) + ax.annotate(f'{pred_2050:.0f}', xy=(2050, pred_2050), xytext=(2055, pred_2050-300), + fontsize=10, color='#3498DB', fontweight='bold') + + # Formatting + ax.set_xlabel('Year', fontsize=12) + ax.set_ylabel('Annual Launches', fontsize=12) + ax.set_title('Richards S-Curve Model Fit (1984-2025 Data)\nRocket Launch Capacity Projection', + fontsize=14, fontweight='bold') + ax.legend(loc='upper left', fontsize=10) + ax.grid(True, alpha=0.3) + ax.set_xlim(1982, 2100) + ax.set_ylim(0, K * 1.15) + + # Add text box with model info + textstr = f'''Model: N(t) = K / (1 + e^(-r(t-t₀)))^(1/v) + +Data Window: {start_year}-2025 ({len(data)} points) +K = {K:.0f} launches/year +r = {r:.4f} +t₀ = {start_year + t0:.1f} +v = {v:.3f} +R² = {r_squared:.4f} + +Note: Physical limit (3,650) shown as +dashed red line''' + + props = dict(boxstyle='round', facecolor='wheat', alpha=0.8) + ax.text(0.98, 0.35, textstr, transform=ax.transAxes, fontsize=9, + verticalalignment='top', horizontalalignment='right', bbox=props, family='monospace') + + plt.tight_layout() + plt.savefig('richards_curve_1984.png', dpi=150, bbox_inches='tight') + plt.close() + + print("\nPlot saved: richards_curve_1984.png") + print("=" * 60) + + +if __name__ == "__main__": + main() diff --git a/p3/water_sensitivity_analysis.py b/p3/water_sensitivity_analysis.py index 2032be5..0af7c1b 100644 --- a/p3/water_sensitivity_analysis.py +++ b/p3/water_sensitivity_analysis.py @@ -1,1500 +1,1500 @@ -""" -任务三：月球殖民地水需求灵敏度分析 - -Water Supply Sensitivity Analysis for Moon Colony - -分析参数: -1. 水回收效率 (η): 70%-95% -2. 舒适度因子 (α): 1-300 -3. 医疗紧急参数: 患病率、每次用水量、置信水平 -4. 人口规模 (N): ±20% -5. 安全缓冲天数: 15-60天 - -输出: -- 单参数灵敏度分析图 -- 多参数交互分析 (Tornado图) -- 最坏情况分析 (Stress Test) -""" - -import numpy as np -import matplotlib -matplotlib.use('Agg') -import matplotlib.pyplot as plt -from matplotlib import rcParams -from scipy import stats -from dataclasses import dataclass, field -from typing import List, Dict, Tuple, Optional -import pandas as pd -from mpl_toolkits.mplot3d import Axes3D -import warnings -warnings.filterwarnings('ignore') - -# 设置字体 -rcParams['font.sans-serif'] = ['Arial Unicode MS', 'DejaVu Sans', 'SimHei'] -rcParams['axes.unicode_minus'] = False -plt.style.use('seaborn-v0_8-whitegrid') - -# ============== 物理常数 ============== -G0 = 9.81 # m/s² -OMEGA_EARTH = 7.27e-5 # rad/s -R_EARTH = 6.371e6 # m - -# ============== 基准参数 ============== -@dataclass -class BaselineParameters: - """基准参数配置""" - # 人口 - population: int = 100_000 - - # 水需求参数 - survival_water: float = 2.5 # L/人/天 - hygiene_water_base: float = 0.4 # L/人/天 (乘以α) - comfort_factor: float = 50.0 # 舒适度系数 - - # 医疗参数 - medical_water_per_event: float = 5.0 # L/次 - sickness_rate: float = 0.02 # 每日患病率 (2%) - confidence_level: float = 0.99 # 置信水平 - - # 回收和缓冲 - recycle_rate: float = 0.90 # 水循环回收率 - safety_buffer_days: int = 30 # 安全缓冲天数 - - # 运输参数 - elevator_capacity_per_year: float = 179_000 # 吨/年/部 - num_elevators: int = 3 - elevator_specific_energy: float = 157.2e9 # J/吨 - - @property - def total_elevator_capacity(self) -> float: - return self.num_elevators * self.elevator_capacity_per_year - - @property - def daily_elevator_capacity(self) -> float: - return self.total_elevator_capacity / 365 - - -# ============== 水需求计算函数 ============== - -def calculate_water_demand( - population: int, - survival_water: float, - hygiene_water_base: float, - comfort_factor: float, - medical_water_per_event: float, - sickness_rate: float, - confidence_level: float, - recycle_rate: float, - safety_buffer_days: int -) -> Dict: - """ - 计算水需求 - - 返回包含各项水需求指标的字典 - """ - # 每人每日用水 - daily_per_person = survival_water + hygiene_water_base * comfort_factor - - # 每日循环用水量 (吨) - daily_circulation = population * daily_per_person / 1000 - - # 医疗用水计算 (使用正态近似) - # 每日患病人数期望和标准差 - n = population - p = sickness_rate - mu = n * p - sigma = np.sqrt(n * p * (1 - p)) - - # 给定置信水平的Z值 - z = stats.norm.ppf(confidence_level) - - # 峰值医疗需求 (99%置信) - peak_medical_events = mu + z * sigma - daily_medical = peak_medical_events * medical_water_per_event / 1000 # 吨 - - # 平均医疗需求 - mean_medical = mu * medical_water_per_event / 1000 # 吨 - - # 每日循环损耗 (需从地球补充) - daily_circulation_loss = daily_circulation * (1 - recycle_rate) - - # 每日总补充量 (损耗 + 医疗) - daily_supply = daily_circulation_loss + mean_medical - - # 每月补充量 - monthly_supply = daily_supply * 30 - - # 年补充量 - annual_supply = daily_supply * 365 - - # 水库初始存量 (支撑单日循环) - reservoir = daily_circulation - - # 安全缓冲 (峰值医疗 + 损耗) - daily_peak_supply = daily_circulation_loss + daily_medical - safety_buffer = daily_peak_supply * safety_buffer_days - - # 首批运输量 = 水库 + 安全缓冲 - initial_transport = reservoir + safety_buffer - - return { - 'population': population, - 'comfort_factor': comfort_factor, - 'recycle_rate': recycle_rate, - 'sickness_rate': sickness_rate, - 'confidence_level': confidence_level, - 'safety_buffer_days': safety_buffer_days, - - 'daily_per_person_liters': daily_per_person, - 'daily_circulation_tons': daily_circulation, - 'daily_medical_tons': daily_medical, - 'mean_medical_tons': mean_medical, - 'daily_circulation_loss_tons': daily_circulation_loss, - 'daily_supply_tons': daily_supply, - 'monthly_supply_tons': monthly_supply, - 'annual_supply_tons': annual_supply, - 'reservoir_tons': reservoir, - 'safety_buffer_tons': safety_buffer, - 'initial_transport_tons': initial_transport, - } - - -def calculate_transport_metrics( - water_demand: Dict, - daily_capacity: float, - specific_energy: float -) -> Dict: - """计算运输指标""" - initial_tons = water_demand['initial_transport_tons'] - monthly_tons = water_demand['monthly_supply_tons'] - annual_tons = water_demand['annual_supply_tons'] - - # 首批运输 - initial_days = initial_tons / daily_capacity - initial_energy_pj = initial_tons * specific_energy / 1e15 - - # 每月补充 - monthly_days = monthly_tons / daily_capacity - monthly_energy_pj = monthly_tons * specific_energy / 1e15 - - # 年度 - annual_days = annual_tons / daily_capacity - annual_energy_pj = annual_tons * specific_energy / 1e15 - - return { - 'initial_days': initial_days, - 'initial_energy_pj': initial_energy_pj, - 'monthly_days': monthly_days, - 'monthly_energy_pj': monthly_energy_pj, - 'annual_days': annual_days, - 'annual_energy_pj': annual_energy_pj, - 'annual_tons': annual_tons, - } - - -# ============== 灵敏度分析函数 ============== - -def sensitivity_recycle_rate(baseline: BaselineParameters, save_dir: str): - """ - 水回收效率灵敏度分析 - η: 70% - 95% - """ - print("分析水回收效率灵敏度...") - - recycle_rates = np.linspace(0.70, 0.95, 26) - - results = { - 'recycle_rate': [], - 'daily_supply': [], - 'annual_supply': [], - 'annual_energy_pj': [], - 'capacity_ratio': [], # 占电梯年运力比例 - } - - for eta in recycle_rates: - demand = calculate_water_demand( - population=baseline.population, - survival_water=baseline.survival_water, - hygiene_water_base=baseline.hygiene_water_base, - comfort_factor=baseline.comfort_factor, - medical_water_per_event=baseline.medical_water_per_event, - sickness_rate=baseline.sickness_rate, - confidence_level=baseline.confidence_level, - recycle_rate=eta, - safety_buffer_days=baseline.safety_buffer_days - ) - - transport = calculate_transport_metrics( - demand, - baseline.daily_elevator_capacity, - baseline.elevator_specific_energy - ) - - results['recycle_rate'].append(eta * 100) - results['daily_supply'].append(demand['daily_supply_tons']) - results['annual_supply'].append(demand['annual_supply_tons']) - results['annual_energy_pj'].append(transport['annual_energy_pj']) - results['capacity_ratio'].append(demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100) - - # 绘图 - fig, axes = plt.subplots(2, 2, figsize=(14, 12)) - - # 图1: 日补充量 - ax = axes[0, 0] - ax.plot(results['recycle_rate'], results['daily_supply'], 'b-', linewidth=2, marker='o', markersize=4) - ax.axvline(x=90, color='r', linestyle='--', label='Baseline (90%)') - ax.set_xlabel('Recycle Rate (%)', fontsize=12) - ax.set_ylabel('Daily Supply (tons)', fontsize=12) - ax.set_title('Daily Water Supply vs Recycle Rate', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - # 图2: 年补充量 - ax = axes[0, 1] - ax.plot(results['recycle_rate'], results['annual_supply'], 'g-', linewidth=2, marker='s', markersize=4) - ax.axvline(x=90, color='r', linestyle='--', label='Baseline (90%)') - ax.set_xlabel('Recycle Rate (%)', fontsize=12) - ax.set_ylabel('Annual Supply (tons)', fontsize=12) - ax.set_title('Annual Water Supply vs Recycle Rate', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - # 图3: 年能耗 - ax = axes[1, 0] - ax.plot(results['recycle_rate'], results['annual_energy_pj'], 'purple', linewidth=2, marker='^', markersize=4) - ax.axvline(x=90, color='r', linestyle='--', label='Baseline (90%)') - ax.set_xlabel('Recycle Rate (%)', fontsize=12) - ax.set_ylabel('Annual Energy (PJ)', fontsize=12) - ax.set_title('Annual Transport Energy vs Recycle Rate', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - # 图4: 运力占比 - ax = axes[1, 1] - ax.plot(results['recycle_rate'], results['capacity_ratio'], 'orange', linewidth=2, marker='d', markersize=4) - ax.axvline(x=90, color='r', linestyle='--', label='Baseline (90%)') - ax.axhline(y=100, color='gray', linestyle=':', label='Full Capacity') - ax.set_xlabel('Recycle Rate (%)', fontsize=12) - ax.set_ylabel('Elevator Capacity Usage (%)', fontsize=12) - ax.set_title('Elevator Capacity Utilization vs Recycle Rate', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - plt.suptitle(f'Sensitivity Analysis: Water Recycle Rate\n(α={baseline.comfort_factor}, N={baseline.population:,})', - fontsize=14, y=1.02) - plt.tight_layout() - plt.savefig(f'{save_dir}/sensitivity_recycle_rate.png', dpi=150, bbox_inches='tight') - print(f" 保存: {save_dir}/sensitivity_recycle_rate.png") - - return pd.DataFrame(results) - - -def sensitivity_comfort_factor(baseline: BaselineParameters, save_dir: str): - """ - 舒适度因子连续灵敏度分析 - α: 1 - 300 - """ - print("分析舒适度因子灵敏度...") - - alphas = np.concatenate([ - np.linspace(1, 10, 10), - np.linspace(15, 50, 8), - np.linspace(60, 150, 10), - np.linspace(175, 300, 6) - ]) - - results = { - 'alpha': [], - 'daily_per_person': [], - 'daily_supply': [], - 'annual_supply': [], - 'annual_energy_pj': [], - 'initial_transport': [], - 'capacity_ratio': [], - } - - for alpha in alphas: - demand = calculate_water_demand( - population=baseline.population, - survival_water=baseline.survival_water, - hygiene_water_base=baseline.hygiene_water_base, - comfort_factor=alpha, - medical_water_per_event=baseline.medical_water_per_event, - sickness_rate=baseline.sickness_rate, - confidence_level=baseline.confidence_level, - recycle_rate=baseline.recycle_rate, - safety_buffer_days=baseline.safety_buffer_days - ) - - transport = calculate_transport_metrics( - demand, - baseline.daily_elevator_capacity, - baseline.elevator_specific_energy - ) - - results['alpha'].append(alpha) - results['daily_per_person'].append(demand['daily_per_person_liters']) - results['daily_supply'].append(demand['daily_supply_tons']) - results['annual_supply'].append(demand['annual_supply_tons']) - results['annual_energy_pj'].append(transport['annual_energy_pj']) - results['initial_transport'].append(demand['initial_transport_tons']) - results['capacity_ratio'].append(demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100) - - # 绘图 - fig, axes = plt.subplots(2, 3, figsize=(18, 12)) - - # 图1: 人均日用水 - ax = axes[0, 0] - ax.plot(results['alpha'], results['daily_per_person'], 'b-', linewidth=2) - ax.axvline(x=1, color='green', linestyle='--', alpha=0.7, label='Survival (α=1)') - ax.axvline(x=50, color='orange', linestyle='--', alpha=0.7, label='Comfort (α=50)') - ax.axvline(x=250, color='red', linestyle='--', alpha=0.7, label='Luxury (α=250)') - ax.set_xlabel('Comfort Factor (α)', fontsize=12) - ax.set_ylabel('Daily Water per Person (L)', fontsize=12) - ax.set_title('Per Capita Water Usage vs Comfort Factor', fontsize=13) - ax.legend(fontsize=9) - ax.grid(True, alpha=0.3) - - # 图2: 日补充量 - ax = axes[0, 1] - ax.plot(results['alpha'], results['daily_supply'], 'g-', linewidth=2) - ax.axvline(x=50, color='r', linestyle='--', label='Baseline (α=50)') - ax.set_xlabel('Comfort Factor (α)', fontsize=12) - ax.set_ylabel('Daily Supply (tons)', fontsize=12) - ax.set_title('Daily Water Supply vs Comfort Factor', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - # 图3: 年补充量 - ax = axes[0, 2] - ax.plot(results['alpha'], results['annual_supply'], 'purple', linewidth=2) - ax.axvline(x=50, color='r', linestyle='--', label='Baseline (α=50)') - ax.set_xlabel('Comfort Factor (α)', fontsize=12) - ax.set_ylabel('Annual Supply (tons)', fontsize=12) - ax.set_title('Annual Water Supply vs Comfort Factor', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - # 图4: 年能耗 - ax = axes[1, 0] - ax.plot(results['alpha'], results['annual_energy_pj'], 'orange', linewidth=2) - ax.axvline(x=50, color='r', linestyle='--', label='Baseline (α=50)') - ax.set_xlabel('Comfort Factor (α)', fontsize=12) - ax.set_ylabel('Annual Energy (PJ)', fontsize=12) - ax.set_title('Annual Transport Energy vs Comfort Factor', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - # 图5: 首批运输量 - ax = axes[1, 1] - ax.plot(results['alpha'], results['initial_transport'], 'teal', linewidth=2) - ax.axvline(x=50, color='r', linestyle='--', label='Baseline (α=50)') - ax.set_xlabel('Comfort Factor (α)', fontsize=12) - ax.set_ylabel('Initial Transport (tons)', fontsize=12) - ax.set_title('Initial Transport Volume vs Comfort Factor', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - # 图6: 运力占比 - ax = axes[1, 2] - ax.plot(results['alpha'], results['capacity_ratio'], 'brown', linewidth=2) - ax.axvline(x=50, color='r', linestyle='--', label='Baseline (α=50)') - ax.axhline(y=100, color='gray', linestyle=':', label='Full Capacity') - ax.fill_between(results['alpha'], 0, results['capacity_ratio'], alpha=0.3) - ax.set_xlabel('Comfort Factor (α)', fontsize=12) - ax.set_ylabel('Elevator Capacity Usage (%)', fontsize=12) - ax.set_title('Elevator Capacity Utilization vs Comfort Factor', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - plt.suptitle(f'Sensitivity Analysis: Comfort Factor (α)\n(η={baseline.recycle_rate*100:.0f}%, N={baseline.population:,})', - fontsize=14, y=1.02) - plt.tight_layout() - plt.savefig(f'{save_dir}/sensitivity_comfort_factor.png', dpi=150, bbox_inches='tight') - print(f" 保存: {save_dir}/sensitivity_comfort_factor.png") - - return pd.DataFrame(results) - - -def sensitivity_medical_params(baseline: BaselineParameters, save_dir: str): - """ - 医疗紧急参数灵敏度分析 - - 患病率: 1% - 5% - - 每次医疗用水: 3 - 15 L - - 置信水平: 95% - 99.9% - """ - print("分析医疗参数灵敏度...") - - fig, axes = plt.subplots(2, 3, figsize=(18, 12)) - - # ========== 患病率分析 ========== - sickness_rates = np.linspace(0.01, 0.05, 21) - results_sr = {'sickness_rate': [], 'daily_medical': [], 'annual_supply': [], 'safety_buffer': []} - - for sr in sickness_rates: - demand = calculate_water_demand( - population=baseline.population, - survival_water=baseline.survival_water, - hygiene_water_base=baseline.hygiene_water_base, - comfort_factor=baseline.comfort_factor, - medical_water_per_event=baseline.medical_water_per_event, - sickness_rate=sr, - confidence_level=baseline.confidence_level, - recycle_rate=baseline.recycle_rate, - safety_buffer_days=baseline.safety_buffer_days - ) - results_sr['sickness_rate'].append(sr * 100) - results_sr['daily_medical'].append(demand['daily_medical_tons']) - results_sr['annual_supply'].append(demand['annual_supply_tons']) - results_sr['safety_buffer'].append(demand['safety_buffer_tons']) - - ax = axes[0, 0] - ax.plot(results_sr['sickness_rate'], results_sr['daily_medical'], 'b-', linewidth=2, marker='o', markersize=4) - ax.axvline(x=2, color='r', linestyle='--', label='Baseline (2%)') - ax.set_xlabel('Daily Sickness Rate (%)', fontsize=12) - ax.set_ylabel('Peak Daily Medical Water (tons)', fontsize=12) - ax.set_title('Peak Medical Water Demand vs Sickness Rate', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - ax = axes[0, 1] - ax.plot(results_sr['sickness_rate'], results_sr['safety_buffer'], 'g-', linewidth=2, marker='s', markersize=4) - ax.axvline(x=2, color='r', linestyle='--', label='Baseline (2%)') - ax.set_xlabel('Daily Sickness Rate (%)', fontsize=12) - ax.set_ylabel('Safety Buffer (tons)', fontsize=12) - ax.set_title('Safety Buffer vs Sickness Rate', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - # ========== 每次医疗用水分析 ========== - medical_waters = np.linspace(3, 15, 13) - results_mw = {'medical_water': [], 'daily_medical': [], 'annual_supply': []} - - for mw in medical_waters: - demand = calculate_water_demand( - population=baseline.population, - survival_water=baseline.survival_water, - hygiene_water_base=baseline.hygiene_water_base, - comfort_factor=baseline.comfort_factor, - medical_water_per_event=mw, - sickness_rate=baseline.sickness_rate, - confidence_level=baseline.confidence_level, - recycle_rate=baseline.recycle_rate, - safety_buffer_days=baseline.safety_buffer_days - ) - results_mw['medical_water'].append(mw) - results_mw['daily_medical'].append(demand['daily_medical_tons']) - results_mw['annual_supply'].append(demand['annual_supply_tons']) - - ax = axes[0, 2] - ax.plot(results_mw['medical_water'], results_mw['daily_medical'], 'purple', linewidth=2, marker='^', markersize=5) - ax.axvline(x=5, color='r', linestyle='--', label='Baseline (5 L)') - ax.set_xlabel('Medical Water per Event (L)', fontsize=12) - ax.set_ylabel('Peak Daily Medical Water (tons)', fontsize=12) - ax.set_title('Peak Medical Water vs Per-Event Usage', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - # ========== 置信水平分析 ========== - confidence_levels = np.linspace(0.90, 0.999, 20) - results_cl = {'confidence': [], 'daily_medical': [], 'safety_buffer': [], 'z_value': []} - - for cl in confidence_levels: - demand = calculate_water_demand( - population=baseline.population, - survival_water=baseline.survival_water, - hygiene_water_base=baseline.hygiene_water_base, - comfort_factor=baseline.comfort_factor, - medical_water_per_event=baseline.medical_water_per_event, - sickness_rate=baseline.sickness_rate, - confidence_level=cl, - recycle_rate=baseline.recycle_rate, - safety_buffer_days=baseline.safety_buffer_days - ) - results_cl['confidence'].append(cl * 100) - results_cl['daily_medical'].append(demand['daily_medical_tons']) - results_cl['safety_buffer'].append(demand['safety_buffer_tons']) - results_cl['z_value'].append(stats.norm.ppf(cl)) - - ax = axes[1, 0] - ax.plot(results_cl['confidence'], results_cl['daily_medical'], 'orange', linewidth=2, marker='d', markersize=4) - ax.axvline(x=99, color='r', linestyle='--', label='Baseline (99%)') - ax.set_xlabel('Confidence Level (%)', fontsize=12) - ax.set_ylabel('Peak Daily Medical Water (tons)', fontsize=12) - ax.set_title('Peak Medical Water vs Confidence Level', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - ax = axes[1, 1] - ax.plot(results_cl['confidence'], results_cl['z_value'], 'teal', linewidth=2, marker='o', markersize=4) - ax.axvline(x=99, color='r', linestyle='--', label='Baseline (99%)') - ax.set_xlabel('Confidence Level (%)', fontsize=12) - ax.set_ylabel('Z-Value', fontsize=12) - ax.set_title('Z-Value vs Confidence Level', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - # ========== 综合影响图 ========== - ax = axes[1, 2] - - # 计算各参数变化对年供水的影响 - baseline_demand = calculate_water_demand( - population=baseline.population, - survival_water=baseline.survival_water, - hygiene_water_base=baseline.hygiene_water_base, - comfort_factor=baseline.comfort_factor, - medical_water_per_event=baseline.medical_water_per_event, - sickness_rate=baseline.sickness_rate, - confidence_level=baseline.confidence_level, - recycle_rate=baseline.recycle_rate, - safety_buffer_days=baseline.safety_buffer_days - ) - baseline_annual = baseline_demand['annual_supply_tons'] - - # 患病率影响 (1% -> 5%) - sr_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, baseline.medical_water_per_event, 0.01, - baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) - sr_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, baseline.medical_water_per_event, 0.05, - baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) - - # 医疗用水影响 (3L -> 15L) - mw_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, 3, baseline.sickness_rate, - baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) - mw_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, 15, baseline.sickness_rate, - baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) - - # 置信水平影响 (95% -> 99.9%) - cl_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, - 0.95, baseline.recycle_rate, baseline.safety_buffer_days) - cl_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, - 0.999, baseline.recycle_rate, baseline.safety_buffer_days) - - params = ['Sickness Rate\n(1%→5%)', 'Medical Water\n(3L→15L)', 'Confidence\n(95%→99.9%)'] - changes_low = [ - (sr_low['annual_supply_tons'] - baseline_annual) / baseline_annual * 100, - (mw_low['annual_supply_tons'] - baseline_annual) / baseline_annual * 100, - (cl_low['annual_supply_tons'] - baseline_annual) / baseline_annual * 100, - ] - changes_high = [ - (sr_high['annual_supply_tons'] - baseline_annual) / baseline_annual * 100, - (mw_high['annual_supply_tons'] - baseline_annual) / baseline_annual * 100, - (cl_high['annual_supply_tons'] - baseline_annual) / baseline_annual * 100, - ] - - x = np.arange(len(params)) - width = 0.35 - ax.barh(x - width/2, changes_low, width, label='Low Value', color='green', alpha=0.7) - ax.barh(x + width/2, changes_high, width, label='High Value', color='red', alpha=0.7) - ax.axvline(x=0, color='black', linewidth=1) - ax.set_yticks(x) - ax.set_yticklabels(params) - ax.set_xlabel('Change in Annual Supply (%)', fontsize=12) - ax.set_title('Medical Parameters Impact on Annual Supply', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3, axis='x') - - plt.suptitle(f'Sensitivity Analysis: Medical Emergency Parameters\n(α={baseline.comfort_factor}, η={baseline.recycle_rate*100:.0f}%)', - fontsize=14, y=1.02) - plt.tight_layout() - plt.savefig(f'{save_dir}/sensitivity_medical_params.png', dpi=150, bbox_inches='tight') - print(f" 保存: {save_dir}/sensitivity_medical_params.png") - - return results_sr, results_mw, results_cl - - -def sensitivity_population(baseline: BaselineParameters, save_dir: str): - """ - 人口规模灵敏度分析 - N: 80,000 - 120,000 (±20%) - """ - print("分析人口规模灵敏度...") - - populations = np.linspace(80000, 120000, 21).astype(int) - - results = { - 'population': [], - 'daily_circulation': [], - 'daily_supply': [], - 'annual_supply': [], - 'annual_energy_pj': [], - 'initial_transport': [], - 'capacity_ratio': [], - } - - for N in populations: - demand = calculate_water_demand( - population=N, - survival_water=baseline.survival_water, - hygiene_water_base=baseline.hygiene_water_base, - comfort_factor=baseline.comfort_factor, - medical_water_per_event=baseline.medical_water_per_event, - sickness_rate=baseline.sickness_rate, - confidence_level=baseline.confidence_level, - recycle_rate=baseline.recycle_rate, - safety_buffer_days=baseline.safety_buffer_days - ) - - transport = calculate_transport_metrics( - demand, - baseline.daily_elevator_capacity, - baseline.elevator_specific_energy - ) - - results['population'].append(N / 1000) # 千人 - results['daily_circulation'].append(demand['daily_circulation_tons']) - results['daily_supply'].append(demand['daily_supply_tons']) - results['annual_supply'].append(demand['annual_supply_tons']) - results['annual_energy_pj'].append(transport['annual_energy_pj']) - results['initial_transport'].append(demand['initial_transport_tons']) - results['capacity_ratio'].append(demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100) - - # 绘图 - fig, axes = plt.subplots(2, 2, figsize=(14, 12)) - - # 图1: 日循环用水 - ax = axes[0, 0] - ax.plot(results['population'], results['daily_circulation'], 'b-', linewidth=2, marker='o', markersize=4) - ax.axvline(x=100, color='r', linestyle='--', label='Baseline (100k)') - ax.fill_between(results['population'], 0, results['daily_circulation'], alpha=0.2) - ax.set_xlabel('Population (thousands)', fontsize=12) - ax.set_ylabel('Daily Circulation Water (tons)', fontsize=12) - ax.set_title('Daily Circulation Water vs Population', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - # 图2: 日补充量 - ax = axes[0, 1] - ax.plot(results['population'], results['daily_supply'], 'g-', linewidth=2, marker='s', markersize=4) - ax.axvline(x=100, color='r', linestyle='--', label='Baseline (100k)') - ax.set_xlabel('Population (thousands)', fontsize=12) - ax.set_ylabel('Daily Supply (tons)', fontsize=12) - ax.set_title('Daily Water Supply vs Population', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - # 图3: 年能耗 - ax = axes[1, 0] - ax.plot(results['population'], results['annual_energy_pj'], 'purple', linewidth=2, marker='^', markersize=4) - ax.axvline(x=100, color='r', linestyle='--', label='Baseline (100k)') - ax.set_xlabel('Population (thousands)', fontsize=12) - ax.set_ylabel('Annual Energy (PJ)', fontsize=12) - ax.set_title('Annual Transport Energy vs Population', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - # 图4: 运力占比 - ax = axes[1, 1] - ax.plot(results['population'], results['capacity_ratio'], 'orange', linewidth=2, marker='d', markersize=4) - ax.axvline(x=100, color='r', linestyle='--', label='Baseline (100k)') - ax.axhline(y=100, color='gray', linestyle=':', label='Full Capacity') - ax.set_xlabel('Population (thousands)', fontsize=12) - ax.set_ylabel('Elevator Capacity Usage (%)', fontsize=12) - ax.set_title('Elevator Capacity Utilization vs Population', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - plt.suptitle(f'Sensitivity Analysis: Population Size\n(α={baseline.comfort_factor}, η={baseline.recycle_rate*100:.0f}%)', - fontsize=14, y=1.02) - plt.tight_layout() - plt.savefig(f'{save_dir}/sensitivity_population.png', dpi=150, bbox_inches='tight') - print(f" 保存: {save_dir}/sensitivity_population.png") - - return pd.DataFrame(results) - - -def sensitivity_buffer_days(baseline: BaselineParameters, save_dir: str): - """ - 安全缓冲天数灵敏度分析 - Days: 15 - 60 - """ - print("分析安全缓冲天数灵敏度...") - - buffer_days = np.arange(15, 61, 3) - - results = { - 'buffer_days': [], - 'safety_buffer_tons': [], - 'initial_transport': [], - 'initial_days': [], - 'initial_energy_pj': [], - } - - for days in buffer_days: - demand = calculate_water_demand( - population=baseline.population, - survival_water=baseline.survival_water, - hygiene_water_base=baseline.hygiene_water_base, - comfort_factor=baseline.comfort_factor, - medical_water_per_event=baseline.medical_water_per_event, - sickness_rate=baseline.sickness_rate, - confidence_level=baseline.confidence_level, - recycle_rate=baseline.recycle_rate, - safety_buffer_days=days - ) - - transport = calculate_transport_metrics( - demand, - baseline.daily_elevator_capacity, - baseline.elevator_specific_energy - ) - - results['buffer_days'].append(days) - results['safety_buffer_tons'].append(demand['safety_buffer_tons']) - results['initial_transport'].append(demand['initial_transport_tons']) - results['initial_days'].append(transport['initial_days']) - results['initial_energy_pj'].append(transport['initial_energy_pj']) - - # 绘图 - fig, axes = plt.subplots(2, 2, figsize=(14, 12)) - - # 图1: 安全缓冲量 - ax = axes[0, 0] - ax.plot(results['buffer_days'], results['safety_buffer_tons'], 'b-', linewidth=2, marker='o', markersize=5) - ax.axvline(x=30, color='r', linestyle='--', label='Baseline (30 days)') - ax.set_xlabel('Safety Buffer Days', fontsize=12) - ax.set_ylabel('Safety Buffer (tons)', fontsize=12) - ax.set_title('Safety Buffer Volume vs Buffer Days', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - # 图2: 首批运输量 - ax = axes[0, 1] - ax.plot(results['buffer_days'], results['initial_transport'], 'g-', linewidth=2, marker='s', markersize=5) - ax.axvline(x=30, color='r', linestyle='--', label='Baseline (30 days)') - ax.set_xlabel('Safety Buffer Days', fontsize=12) - ax.set_ylabel('Initial Transport (tons)', fontsize=12) - ax.set_title('Initial Transport Volume vs Buffer Days', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - # 图3: 首批运输天数 - ax = axes[1, 0] - ax.plot(results['buffer_days'], results['initial_days'], 'purple', linewidth=2, marker='^', markersize=5) - ax.axvline(x=30, color='r', linestyle='--', label='Baseline (30 days)') - ax.set_xlabel('Safety Buffer Days', fontsize=12) - ax.set_ylabel('Initial Transport Time (days)', fontsize=12) - ax.set_title('Initial Transport Duration vs Buffer Days', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - # 图4: 首批运输能耗 - ax = axes[1, 1] - ax.plot(results['buffer_days'], results['initial_energy_pj'], 'orange', linewidth=2, marker='d', markersize=5) - ax.axvline(x=30, color='r', linestyle='--', label='Baseline (30 days)') - ax.set_xlabel('Safety Buffer Days', fontsize=12) - ax.set_ylabel('Initial Transport Energy (PJ)', fontsize=12) - ax.set_title('Initial Transport Energy vs Buffer Days', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3) - - plt.suptitle(f'Sensitivity Analysis: Safety Buffer Days\n(α={baseline.comfort_factor}, η={baseline.recycle_rate*100:.0f}%)', - fontsize=14, y=1.02) - plt.tight_layout() - plt.savefig(f'{save_dir}/sensitivity_buffer_days.png', dpi=150, bbox_inches='tight') - print(f" 保存: {save_dir}/sensitivity_buffer_days.png") - - return pd.DataFrame(results) - - -# ============== 多参数交互分析 ============== - -def tornado_analysis(baseline: BaselineParameters, save_dir: str): - """ - Tornado图 - 多参数灵敏度排序分析 - """ - print("生成Tornado图...") - - # 基准值 - baseline_demand = calculate_water_demand( - population=baseline.population, - survival_water=baseline.survival_water, - hygiene_water_base=baseline.hygiene_water_base, - comfort_factor=baseline.comfort_factor, - medical_water_per_event=baseline.medical_water_per_event, - sickness_rate=baseline.sickness_rate, - confidence_level=baseline.confidence_level, - recycle_rate=baseline.recycle_rate, - safety_buffer_days=baseline.safety_buffer_days - ) - baseline_annual = baseline_demand['annual_supply_tons'] - baseline_initial = baseline_demand['initial_transport_tons'] - - # 定义参数范围 - params = [ - ('Recycle Rate (η)', 'recycle_rate', 0.80, 0.95, baseline.recycle_rate), - ('Comfort Factor (α)', 'comfort_factor', 25, 100, baseline.comfort_factor), - ('Population (N)', 'population', 80000, 120000, baseline.population), - ('Sickness Rate', 'sickness_rate', 0.01, 0.04, baseline.sickness_rate), - ('Medical Water/Event', 'medical_water', 3, 10, baseline.medical_water_per_event), - ('Buffer Days', 'buffer_days', 15, 45, baseline.safety_buffer_days), - ('Confidence Level', 'confidence', 0.95, 0.999, baseline.confidence_level), - ] - - results_annual = [] - results_initial = [] - - for name, param_type, low, high, base_val in params: - # 低值计算 - if param_type == 'recycle_rate': - demand_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, - baseline.confidence_level, low, baseline.safety_buffer_days) - demand_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, - baseline.confidence_level, high, baseline.safety_buffer_days) - elif param_type == 'comfort_factor': - demand_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - low, baseline.medical_water_per_event, baseline.sickness_rate, - baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) - demand_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - high, baseline.medical_water_per_event, baseline.sickness_rate, - baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) - elif param_type == 'population': - demand_low = calculate_water_demand(int(low), baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, - baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) - demand_high = calculate_water_demand(int(high), baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, - baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) - elif param_type == 'sickness_rate': - demand_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, baseline.medical_water_per_event, low, - baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) - demand_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, baseline.medical_water_per_event, high, - baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) - elif param_type == 'medical_water': - demand_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, low, baseline.sickness_rate, - baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) - demand_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, high, baseline.sickness_rate, - baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) - elif param_type == 'buffer_days': - demand_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, - baseline.confidence_level, baseline.recycle_rate, int(low)) - demand_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, - baseline.confidence_level, baseline.recycle_rate, int(high)) - elif param_type == 'confidence': - demand_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, - low, baseline.recycle_rate, baseline.safety_buffer_days) - demand_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, - high, baseline.recycle_rate, baseline.safety_buffer_days) - - # 计算变化百分比 - change_low_annual = (demand_low['annual_supply_tons'] - baseline_annual) / baseline_annual * 100 - change_high_annual = (demand_high['annual_supply_tons'] - baseline_annual) / baseline_annual * 100 - change_low_initial = (demand_low['initial_transport_tons'] - baseline_initial) / baseline_initial * 100 - change_high_initial = (demand_high['initial_transport_tons'] - baseline_initial) / baseline_initial * 100 - - # 确保排序一致性 (swing = |high - low|) - swing_annual = abs(change_high_annual - change_low_annual) - swing_initial = abs(change_high_initial - change_low_initial) - - results_annual.append({ - 'param': name, - 'low': min(change_low_annual, change_high_annual), - 'high': max(change_low_annual, change_high_annual), - 'swing': swing_annual - }) - results_initial.append({ - 'param': name, - 'low': min(change_low_initial, change_high_initial), - 'high': max(change_low_initial, change_high_initial), - 'swing': swing_initial - }) - - # 按swing排序 - results_annual = sorted(results_annual, key=lambda x: x['swing'], reverse=True) - results_initial = sorted(results_initial, key=lambda x: x['swing'], reverse=True) - - # 绘图 - fig, axes = plt.subplots(1, 2, figsize=(16, 8)) - - # 图1: 年补充量Tornado - ax = axes[0] - y_pos = np.arange(len(results_annual)) - params_sorted = [r['param'] for r in results_annual] - lows = [r['low'] for r in results_annual] - highs = [r['high'] for r in results_annual] - - # 绘制条形图 - for i, (l, h) in enumerate(zip(lows, highs)): - if l < 0: - ax.barh(i, l, color='green', alpha=0.7, height=0.6) - if h > 0: - ax.barh(i, h, color='red', alpha=0.7, height=0.6) - if l >= 0: - ax.barh(i, l, left=0, color='green', alpha=0.7, height=0.6) - if h <= 0: - ax.barh(i, h, left=0, color='red', alpha=0.7, height=0.6) - # 如果跨越0，绘制完整范围 - if l < 0 < h: - pass # 已经绘制 - elif l >= 0: - ax.barh(i, h - l, left=l, color='orange', alpha=0.7, height=0.6) - elif h <= 0: - ax.barh(i, h - l, left=l, color='orange', alpha=0.7, height=0.6) - - ax.axvline(x=0, color='black', linewidth=1.5) - ax.set_yticks(y_pos) - ax.set_yticklabels(params_sorted) - ax.set_xlabel('Change in Annual Supply (%)', fontsize=12) - ax.set_title('Tornado Diagram: Annual Water Supply', fontsize=14) - ax.grid(True, alpha=0.3, axis='x') - - # 图2: 首批运输量Tornado - ax = axes[1] - params_sorted = [r['param'] for r in results_initial] - lows = [r['low'] for r in results_initial] - highs = [r['high'] for r in results_initial] - - for i, (l, h) in enumerate(zip(lows, highs)): - width = h - l - left = l - color = 'steelblue' if width > 0 else 'gray' - ax.barh(i, width, left=left, color=color, alpha=0.7, height=0.6) - - ax.axvline(x=0, color='black', linewidth=1.5) - ax.set_yticks(y_pos) - ax.set_yticklabels(params_sorted) - ax.set_xlabel('Change in Initial Transport (%)', fontsize=12) - ax.set_title('Tornado Diagram: Initial Transport Volume', fontsize=14) - ax.grid(True, alpha=0.3, axis='x') - - # 添加图例说明 - from matplotlib.patches import Patch - legend_elements = [ - Patch(facecolor='green', alpha=0.7, label='Decrease'), - Patch(facecolor='red', alpha=0.7, label='Increase'), - ] - axes[0].legend(handles=legend_elements, loc='lower right') - - plt.suptitle('Multi-Parameter Sensitivity Analysis (Tornado Diagram)\nBaseline: α=50, η=90%, N=100k', - fontsize=14, y=1.02) - plt.tight_layout() - plt.savefig(f'{save_dir}/tornado_analysis.png', dpi=150, bbox_inches='tight') - print(f" 保存: {save_dir}/tornado_analysis.png") - - return results_annual, results_initial - - -def interaction_heatmap(baseline: BaselineParameters, save_dir: str): - """ - 双参数交互热力图 - α vs η - """ - print("生成双参数交互热力图...") - - alphas = np.array([1, 25, 50, 100, 150, 200, 250]) - etas = np.array([0.70, 0.75, 0.80, 0.85, 0.90, 0.95]) - - # 年补充量矩阵 - annual_matrix = np.zeros((len(etas), len(alphas))) - # 运力占比矩阵 - capacity_matrix = np.zeros((len(etas), len(alphas))) - - for i, eta in enumerate(etas): - for j, alpha in enumerate(alphas): - demand = calculate_water_demand( - population=baseline.population, - survival_water=baseline.survival_water, - hygiene_water_base=baseline.hygiene_water_base, - comfort_factor=alpha, - medical_water_per_event=baseline.medical_water_per_event, - sickness_rate=baseline.sickness_rate, - confidence_level=baseline.confidence_level, - recycle_rate=eta, - safety_buffer_days=baseline.safety_buffer_days - ) - annual_matrix[i, j] = demand['annual_supply_tons'] / 1000 # 千吨 - capacity_matrix[i, j] = demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100 - - fig, axes = plt.subplots(1, 2, figsize=(16, 7)) - - # 图1: 年补充量热力图 - ax = axes[0] - im = ax.imshow(annual_matrix, cmap='YlOrRd', aspect='auto') - ax.set_xticks(range(len(alphas))) - ax.set_xticklabels([f'{a}' for a in alphas]) - ax.set_yticks(range(len(etas))) - ax.set_yticklabels([f'{e*100:.0f}%' for e in etas]) - ax.set_xlabel('Comfort Factor (α)', fontsize=12) - ax.set_ylabel('Recycle Rate (η)', fontsize=12) - ax.set_title('Annual Water Supply (kt)', fontsize=13) - - # 添加数值标注 - for i in range(len(etas)): - for j in range(len(alphas)): - text_color = 'white' if annual_matrix[i, j] > annual_matrix.max() * 0.6 else 'black' - ax.text(j, i, f'{annual_matrix[i, j]:.0f}', ha='center', va='center', fontsize=9, color=text_color) - - plt.colorbar(im, ax=ax, label='Annual Supply (kt)') - - # 图2: 运力占比热力图 - ax = axes[1] - im = ax.imshow(capacity_matrix, cmap='RdYlGn_r', aspect='auto', vmin=0, vmax=150) - ax.set_xticks(range(len(alphas))) - ax.set_xticklabels([f'{a}' for a in alphas]) - ax.set_yticks(range(len(etas))) - ax.set_yticklabels([f'{e*100:.0f}%' for e in etas]) - ax.set_xlabel('Comfort Factor (α)', fontsize=12) - ax.set_ylabel('Recycle Rate (η)', fontsize=12) - ax.set_title('Elevator Capacity Usage (%)', fontsize=13) - - # 添加数值标注 - for i in range(len(etas)): - for j in range(len(alphas)): - val = capacity_matrix[i, j] - text_color = 'white' if val > 60 else 'black' - ax.text(j, i, f'{val:.0f}%', ha='center', va='center', fontsize=9, color=text_color) - - # 添加100%等高线 - cbar = plt.colorbar(im, ax=ax, label='Capacity Usage (%)') - - plt.suptitle('Two-Parameter Interaction Analysis: α vs η\n(N=100,000, Buffer=30 days)', - fontsize=14, y=1.02) - plt.tight_layout() - plt.savefig(f'{save_dir}/interaction_heatmap.png', dpi=150, bbox_inches='tight') - print(f" 保存: {save_dir}/interaction_heatmap.png") - - return annual_matrix, capacity_matrix - - -def three_param_surface(baseline: BaselineParameters, save_dir: str): - """ - 三参数交互3D曲面图 - """ - print("生成三参数3D曲面图...") - - alphas = np.linspace(1, 200, 20) - etas = np.linspace(0.75, 0.95, 20) - - A, E = np.meshgrid(alphas, etas) - Z = np.zeros_like(A) - - for i in range(len(etas)): - for j in range(len(alphas)): - demand = calculate_water_demand( - population=baseline.population, - survival_water=baseline.survival_water, - hygiene_water_base=baseline.hygiene_water_base, - comfort_factor=alphas[j], - medical_water_per_event=baseline.medical_water_per_event, - sickness_rate=baseline.sickness_rate, - confidence_level=baseline.confidence_level, - recycle_rate=etas[i], - safety_buffer_days=baseline.safety_buffer_days - ) - Z[i, j] = demand['annual_supply_tons'] / 1000 # 千吨 - - fig = plt.figure(figsize=(14, 10)) - ax = fig.add_subplot(111, projection='3d') - - surf = ax.plot_surface(A, E * 100, Z, cmap='viridis', alpha=0.8, edgecolor='none') - - # 添加基准点 - baseline_demand = calculate_water_demand( - baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, - baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days - ) - ax.scatter([baseline.comfort_factor], [baseline.recycle_rate * 100], - [baseline_demand['annual_supply_tons'] / 1000], - color='red', s=100, label='Baseline', marker='*') - - ax.set_xlabel('Comfort Factor (α)', fontsize=11) - ax.set_ylabel('Recycle Rate (%)', fontsize=11) - ax.set_zlabel('Annual Supply (kt)', fontsize=11) - ax.set_title('3D Surface: Annual Water Supply\nas Function of α and η', fontsize=14) - - fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, label='Annual Supply (kt)') - ax.legend() - - plt.tight_layout() - plt.savefig(f'{save_dir}/three_param_surface.png', dpi=150, bbox_inches='tight') - print(f" 保存: {save_dir}/three_param_surface.png") - - return Z - - -# ============== 最坏情况分析 ============== - -def worst_case_analysis(baseline: BaselineParameters, save_dir: str): - """ - 最坏情况压力测试分析 - """ - print("进行最坏情况分析...") - - # 定义场景 - scenarios = { - 'Baseline': { - 'alpha': baseline.comfort_factor, - 'eta': baseline.recycle_rate, - 'sickness': baseline.sickness_rate, - 'population': baseline.population, - 'buffer': baseline.safety_buffer_days, - 'medical_water': baseline.medical_water_per_event, - }, - 'Survival Mode': { - 'alpha': 1, - 'eta': 0.90, - 'sickness': 0.02, - 'population': 100000, - 'buffer': 30, - 'medical_water': 5, - }, - 'Comfort Mode': { - 'alpha': 50, - 'eta': 0.90, - 'sickness': 0.02, - 'population': 100000, - 'buffer': 30, - 'medical_water': 5, - }, - 'Luxury Mode': { - 'alpha': 250, - 'eta': 0.90, - 'sickness': 0.02, - 'population': 100000, - 'buffer': 30, - 'medical_water': 5, - }, - 'Recycle Degradation': { - 'alpha': 50, - 'eta': 0.80, - 'sickness': 0.02, - 'population': 100000, - 'buffer': 30, - 'medical_water': 5, - }, - 'Epidemic Outbreak': { - 'alpha': 50, - 'eta': 0.90, - 'sickness': 0.05, - 'population': 100000, - 'buffer': 45, - 'medical_water': 10, - }, - 'Population Surge': { - 'alpha': 50, - 'eta': 0.90, - 'sickness': 0.02, - 'population': 120000, - 'buffer': 30, - 'medical_water': 5, - }, - 'Worst Case': { - 'alpha': 250, - 'eta': 0.80, - 'sickness': 0.05, - 'population': 120000, - 'buffer': 45, - 'medical_water': 10, - }, - 'Moderate Stress': { - 'alpha': 100, - 'eta': 0.85, - 'sickness': 0.03, - 'population': 110000, - 'buffer': 40, - 'medical_water': 7, - }, - } - - results = [] - - for name, params in scenarios.items(): - demand = calculate_water_demand( - population=params['population'], - survival_water=baseline.survival_water, - hygiene_water_base=baseline.hygiene_water_base, - comfort_factor=params['alpha'], - medical_water_per_event=params['medical_water'], - sickness_rate=params['sickness'], - confidence_level=baseline.confidence_level, - recycle_rate=params['eta'], - safety_buffer_days=params['buffer'] - ) - - transport = calculate_transport_metrics( - demand, - baseline.daily_elevator_capacity, - baseline.elevator_specific_energy - ) - - capacity_ratio = demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100 - - results.append({ - 'scenario': name, - 'alpha': params['alpha'], - 'eta': params['eta'] * 100, - 'sickness': params['sickness'] * 100, - 'population': params['population'] / 1000, - 'buffer': params['buffer'], - 'daily_supply': demand['daily_supply_tons'], - 'annual_supply': demand['annual_supply_tons'], - 'initial_transport': demand['initial_transport_tons'], - 'annual_energy_pj': transport['annual_energy_pj'], - 'capacity_ratio': capacity_ratio, - 'feasible': capacity_ratio < 100, - }) - - df = pd.DataFrame(results) - - # 绘图 - fig, axes = plt.subplots(2, 2, figsize=(16, 14)) - - # 图1: 年补充量对比 - ax = axes[0, 0] - colors = ['green' if r['feasible'] else 'red' for r in results] - bars = ax.barh(range(len(results)), [r['annual_supply'] / 1000 for r in results], color=colors, alpha=0.7) - ax.axvline(x=baseline.total_elevator_capacity / 1000, color='red', linestyle='--', linewidth=2, label='Elevator Capacity') - ax.set_yticks(range(len(results))) - ax.set_yticklabels([r['scenario'] for r in results]) - ax.set_xlabel('Annual Supply (kt)', fontsize=12) - ax.set_title('Annual Water Supply by Scenario', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3, axis='x') - - # 添加数值标注 - for i, bar in enumerate(bars): - width = bar.get_width() - ax.text(width + 5, bar.get_y() + bar.get_height()/2, f'{width:.0f}', - va='center', fontsize=9) - - # 图2: 运力占比 - ax = axes[0, 1] - colors = ['green' if r['feasible'] else 'red' for r in results] - bars = ax.barh(range(len(results)), [r['capacity_ratio'] for r in results], color=colors, alpha=0.7) - ax.axvline(x=100, color='red', linestyle='--', linewidth=2, label='100% Capacity') - ax.set_yticks(range(len(results))) - ax.set_yticklabels([r['scenario'] for r in results]) - ax.set_xlabel('Elevator Capacity Usage (%)', fontsize=12) - ax.set_title('Capacity Utilization by Scenario', fontsize=13) - ax.legend() - ax.grid(True, alpha=0.3, axis='x') - - for i, bar in enumerate(bars): - width = bar.get_width() - ax.text(width + 2, bar.get_y() + bar.get_height()/2, f'{width:.1f}%', - va='center', fontsize=9) - - # 图3: 年能耗对比 - ax = axes[1, 0] - ax.bar(range(len(results)), [r['annual_energy_pj'] for r in results], color='purple', alpha=0.7) - ax.set_xticks(range(len(results))) - ax.set_xticklabels([r['scenario'] for r in results], rotation=45, ha='right', fontsize=9) - ax.set_ylabel('Annual Energy (PJ)', fontsize=12) - ax.set_title('Annual Transport Energy by Scenario', fontsize=13) - ax.grid(True, alpha=0.3, axis='y') - - # 图4: 首批运输量 - ax = axes[1, 1] - ax.bar(range(len(results)), [r['initial_transport'] / 1000 for r in results], color='teal', alpha=0.7) - ax.set_xticks(range(len(results))) - ax.set_xticklabels([r['scenario'] for r in results], rotation=45, ha='right', fontsize=9) - ax.set_ylabel('Initial Transport (kt)', fontsize=12) - ax.set_title('Initial Transport Volume by Scenario', fontsize=13) - ax.grid(True, alpha=0.3, axis='y') - - plt.suptitle('Worst Case Analysis: Stress Test Scenarios', fontsize=14, y=1.02) - plt.tight_layout() - plt.savefig(f'{save_dir}/worst_case_analysis.png', dpi=150, bbox_inches='tight') - print(f" 保存: {save_dir}/worst_case_analysis.png") - - # 保存详细结果表 - df.to_csv(f'{save_dir}/worst_case_results.csv', index=False) - print(f" 保存: {save_dir}/worst_case_results.csv") - - return df - - -def generate_summary_report(baseline: BaselineParameters, save_dir: str): - """ - 生成灵敏度分析总结报告 - """ - print("生成总结报告...") - - report = [] - report.append("=" * 100) - report.append("TASK 3: WATER SUPPLY SENSITIVITY ANALYSIS SUMMARY") - report.append("=" * 100) - - report.append("\n## 基准参数") - report.append(f"- 人口: {baseline.population:,}") - report.append(f"- 舒适度因子 (α): {baseline.comfort_factor}") - report.append(f"- 水回收效率 (η): {baseline.recycle_rate * 100:.0f}%") - report.append(f"- 患病率: {baseline.sickness_rate * 100:.1f}%") - report.append(f"- 医疗用水: {baseline.medical_water_per_event} L/次") - report.append(f"- 置信水平: {baseline.confidence_level * 100:.1f}%") - report.append(f"- 安全缓冲: {baseline.safety_buffer_days} 天") - - # 计算基准值 - baseline_demand = calculate_water_demand( - baseline.population, baseline.survival_water, baseline.hygiene_water_base, - baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, - baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days - ) - baseline_transport = calculate_transport_metrics( - baseline_demand, baseline.daily_elevator_capacity, baseline.elevator_specific_energy - ) - - report.append(f"\n## 基准水需求") - report.append(f"- 每日补充量: {baseline_demand['daily_supply_tons']:.1f} 吨") - report.append(f"- 年补充量: {baseline_demand['annual_supply_tons']:.1f} 吨") - report.append(f"- 首批运输量: {baseline_demand['initial_transport_tons']:.1f} 吨") - report.append(f"- 年能耗: {baseline_transport['annual_energy_pj']:.4f} PJ") - report.append(f"- 运力占比: {baseline_demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100:.2f}%") - - report.append("\n" + "=" * 100) - report.append("## 关键灵敏度发现") - report.append("=" * 100) - - report.append(""" -### 1. 水回收效率 (η) - 影响程度: ★★★★★ -- 范围: 70% - 95% -- η从90%降至80%时，日补充量翻倍 -- 这是对年运输需求影响最大的参数 -- 建议: 优先投资水回收系统维护和升级 - -### 2. 舒适度因子 (α) - 影响程度: ★★★★★ -- 范围: 1 - 300 -- α=250时年补充量是α=1时的35倍 -- 在高舒适度下可能超出电梯运力 -- 建议: 初期采用生存标准，逐步提升舒适度 - -### 3. 人口规模 (N) - 影响程度: ★★★★☆ -- 范围: ±20% (80k - 120k) -- 线性影响，20%人口增加导致20%需求增加 -- 建议: 规划时预留弹性容量 - -### 4. 医疗紧急参数 - 影响程度: ★★☆☆☆ -- 患病率(1%-5%)和医疗用水(3-15L)对年补充量影响较小 -- 主要影响安全缓冲和首批运输量 -- 建议: 采用保守参数确保医疗安全 - -### 5. 安全缓冲天数 - 影响程度: ★★★☆☆ -- 范围: 15 - 60天 -- 主要影响首批运输量，不影响年运输总量 -- 建议: 保持30天缓冲作为运输中断容忍度 -""") - - report.append("\n" + "=" * 100) - report.append("## 最坏情况分析结论") - report.append("=" * 100) - - report.append(""" -### 可行性边界 -- 最坏情况 (α=250, η=80%, N=120k): 年需求约 562 kt，超出电梯运力 -- 中度压力 (α=100, η=85%, N=110k): 年需求约 132 kt，仍可承受 -- 系统崩溃临界点: 约 η < 82% + α > 200 组合 - -### 风险缓解建议 -1. 水回收系统冗余设计，确保η≥85% -2. 分阶段提升舒适度标准 -3. 维持混合运输能力作为备用 -4. 建立月球本地水源开发计划 -""") - - report_text = "\n".join(report) - - with open(f'{save_dir}/sensitivity_summary_report.txt', 'w', encoding='utf-8') as f: - f.write(report_text) - print(f" 保存: {save_dir}/sensitivity_summary_report.txt") - - return report_text - - -# ============== 主程序 ============== - -def run_full_sensitivity_analysis(): - """运行完整的灵敏度分析""" - - print("=" * 80) - print("任务三：月球殖民地水需求灵敏度分析") - print("=" * 80) - - save_dir = '/Volumes/Files/code/mm/20260130_b/p3' - - # 初始化基准参数 - baseline = BaselineParameters() - - print(f"\n基准参数:") - print(f" - 人口: {baseline.population:,}") - print(f" - 舒适度因子: {baseline.comfort_factor}") - print(f" - 水回收效率: {baseline.recycle_rate * 100:.0f}%") - print(f" - 患病率: {baseline.sickness_rate * 100:.1f}%") - print(f" - 安全缓冲: {baseline.safety_buffer_days} 天") - - print("\n" + "=" * 80) - print("开始灵敏度分析...") - print("=" * 80) - - # 1. 单参数灵敏度分析 - print("\n[1/7] 水回收效率灵敏度分析") - df_recycle = sensitivity_recycle_rate(baseline, save_dir) - - print("\n[2/7] 舒适度因子灵敏度分析") - df_comfort = sensitivity_comfort_factor(baseline, save_dir) - - print("\n[3/7] 医疗参数灵敏度分析") - results_medical = sensitivity_medical_params(baseline, save_dir) - - print("\n[4/7] 人口规模灵敏度分析") - df_population = sensitivity_population(baseline, save_dir) - - print("\n[5/7] 安全缓冲天数灵敏度分析") - df_buffer = sensitivity_buffer_days(baseline, save_dir) - - # 2. 多参数交互分析 - print("\n[6/7] 多参数交互分析") - tornado_results = tornado_analysis(baseline, save_dir) - heatmap_data = interaction_heatmap(baseline, save_dir) - surface_data = three_param_surface(baseline, save_dir) - - # 3. 最坏情况分析 - print("\n[7/7] 最坏情况分析") - worst_case_df = worst_case_analysis(baseline, save_dir) - - # 4. 生成总结报告 - print("\n生成总结报告...") - summary = generate_summary_report(baseline, save_dir) - - print("\n" + "=" * 80) - print("灵敏度分析完成!") - print("=" * 80) - - print(f"\n生成的文件:") - print(f" - sensitivity_recycle_rate.png") - print(f" - sensitivity_comfort_factor.png") - print(f" - sensitivity_medical_params.png") - print(f" - sensitivity_population.png") - print(f" - sensitivity_buffer_days.png") - print(f" - tornado_analysis.png") - print(f" - interaction_heatmap.png") - print(f" - three_param_surface.png") - print(f" - worst_case_analysis.png") - print(f" - worst_case_results.csv") - print(f" - sensitivity_summary_report.txt") - - return { - 'recycle': df_recycle, - 'comfort': df_comfort, - 'population': df_population, - 'buffer': df_buffer, - 'worst_case': worst_case_df, - } - - -if __name__ == "__main__": - results = run_full_sensitivity_analysis() +""" +任务三：月球殖民地水需求灵敏度分析 + +Water Supply Sensitivity Analysis for Moon Colony + +分析参数: +1. 水回收效率 (η): 70%-95% +2. 舒适度因子 (α): 1-300 +3. 医疗紧急参数: 患病率、每次用水量、置信水平 +4. 人口规模 (N): ±20% +5. 安全缓冲天数: 15-60天 + +输出: +- 单参数灵敏度分析图 +- 多参数交互分析 (Tornado图) +- 最坏情况分析 (Stress Test) +""" + +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from matplotlib import rcParams +from scipy import stats +from dataclasses import dataclass, field +from typing import List, Dict, Tuple, Optional +import pandas as pd +from mpl_toolkits.mplot3d import Axes3D +import warnings +warnings.filterwarnings('ignore') + +# 设置字体 +rcParams['font.sans-serif'] = ['Arial Unicode MS', 'DejaVu Sans', 'SimHei'] +rcParams['axes.unicode_minus'] = False +plt.style.use('seaborn-v0_8-whitegrid') + +# ============== 物理常数 ============== +G0 = 9.81 # m/s² +OMEGA_EARTH = 7.27e-5 # rad/s +R_EARTH = 6.371e6 # m + +# ============== 基准参数 ============== +@dataclass +class BaselineParameters: + """基准参数配置""" + # 人口 + population: int = 100_000 + + # 水需求参数 + survival_water: float = 2.5 # L/人/天 + hygiene_water_base: float = 0.4 # L/人/天 (乘以α) + comfort_factor: float = 50.0 # 舒适度系数 + + # 医疗参数 + medical_water_per_event: float = 5.0 # L/次 + sickness_rate: float = 0.02 # 每日患病率 (2%) + confidence_level: float = 0.99 # 置信水平 + + # 回收和缓冲 + recycle_rate: float = 0.90 # 水循环回收率 + safety_buffer_days: int = 30 # 安全缓冲天数 + + # 运输参数 + elevator_capacity_per_year: float = 179_000 # 吨/年/部 + num_elevators: int = 3 + elevator_specific_energy: float = 157.2e9 # J/吨 + + @property + def total_elevator_capacity(self) -> float: + return self.num_elevators * self.elevator_capacity_per_year + + @property + def daily_elevator_capacity(self) -> float: + return self.total_elevator_capacity / 365 + + +# ============== 水需求计算函数 ============== + +def calculate_water_demand( + population: int, + survival_water: float, + hygiene_water_base: float, + comfort_factor: float, + medical_water_per_event: float, + sickness_rate: float, + confidence_level: float, + recycle_rate: float, + safety_buffer_days: int +) -> Dict: + """ + 计算水需求 + + 返回包含各项水需求指标的字典 + """ + # 每人每日用水 + daily_per_person = survival_water + hygiene_water_base * comfort_factor + + # 每日循环用水量 (吨) + daily_circulation = population * daily_per_person / 1000 + + # 医疗用水计算 (使用正态近似) + # 每日患病人数期望和标准差 + n = population + p = sickness_rate + mu = n * p + sigma = np.sqrt(n * p * (1 - p)) + + # 给定置信水平的Z值 + z = stats.norm.ppf(confidence_level) + + # 峰值医疗需求 (99%置信) + peak_medical_events = mu + z * sigma + daily_medical = peak_medical_events * medical_water_per_event / 1000 # 吨 + + # 平均医疗需求 + mean_medical = mu * medical_water_per_event / 1000 # 吨 + + # 每日循环损耗 (需从地球补充) + daily_circulation_loss = daily_circulation * (1 - recycle_rate) + + # 每日总补充量 (损耗 + 医疗) + daily_supply = daily_circulation_loss + mean_medical + + # 每月补充量 + monthly_supply = daily_supply * 30 + + # 年补充量 + annual_supply = daily_supply * 365 + + # 水库初始存量 (支撑单日循环) + reservoir = daily_circulation + + # 安全缓冲 (峰值医疗 + 损耗) + daily_peak_supply = daily_circulation_loss + daily_medical + safety_buffer = daily_peak_supply * safety_buffer_days + + # 首批运输量 = 水库 + 安全缓冲 + initial_transport = reservoir + safety_buffer + + return { + 'population': population, + 'comfort_factor': comfort_factor, + 'recycle_rate': recycle_rate, + 'sickness_rate': sickness_rate, + 'confidence_level': confidence_level, + 'safety_buffer_days': safety_buffer_days, + + 'daily_per_person_liters': daily_per_person, + 'daily_circulation_tons': daily_circulation, + 'daily_medical_tons': daily_medical, + 'mean_medical_tons': mean_medical, + 'daily_circulation_loss_tons': daily_circulation_loss, + 'daily_supply_tons': daily_supply, + 'monthly_supply_tons': monthly_supply, + 'annual_supply_tons': annual_supply, + 'reservoir_tons': reservoir, + 'safety_buffer_tons': safety_buffer, + 'initial_transport_tons': initial_transport, + } + + +def calculate_transport_metrics( + water_demand: Dict, + daily_capacity: float, + specific_energy: float +) -> Dict: + """计算运输指标""" + initial_tons = water_demand['initial_transport_tons'] + monthly_tons = water_demand['monthly_supply_tons'] + annual_tons = water_demand['annual_supply_tons'] + + # 首批运输 + initial_days = initial_tons / daily_capacity + initial_energy_pj = initial_tons * specific_energy / 1e15 + + # 每月补充 + monthly_days = monthly_tons / daily_capacity + monthly_energy_pj = monthly_tons * specific_energy / 1e15 + + # 年度 + annual_days = annual_tons / daily_capacity + annual_energy_pj = annual_tons * specific_energy / 1e15 + + return { + 'initial_days': initial_days, + 'initial_energy_pj': initial_energy_pj, + 'monthly_days': monthly_days, + 'monthly_energy_pj': monthly_energy_pj, + 'annual_days': annual_days, + 'annual_energy_pj': annual_energy_pj, + 'annual_tons': annual_tons, + } + + +# ============== 灵敏度分析函数 ============== + +def sensitivity_recycle_rate(baseline: BaselineParameters, save_dir: str): + """ + 水回收效率灵敏度分析 + η: 70% - 95% + """ + print("分析水回收效率灵敏度...") + + recycle_rates = np.linspace(0.70, 0.95, 26) + + results = { + 'recycle_rate': [], + 'daily_supply': [], + 'annual_supply': [], + 'annual_energy_pj': [], + 'capacity_ratio': [], # 占电梯年运力比例 + } + + for eta in recycle_rates: + demand = calculate_water_demand( + population=baseline.population, + survival_water=baseline.survival_water, + hygiene_water_base=baseline.hygiene_water_base, + comfort_factor=baseline.comfort_factor, + medical_water_per_event=baseline.medical_water_per_event, + sickness_rate=baseline.sickness_rate, + confidence_level=baseline.confidence_level, + recycle_rate=eta, + safety_buffer_days=baseline.safety_buffer_days + ) + + transport = calculate_transport_metrics( + demand, + baseline.daily_elevator_capacity, + baseline.elevator_specific_energy + ) + + results['recycle_rate'].append(eta * 100) + results['daily_supply'].append(demand['daily_supply_tons']) + results['annual_supply'].append(demand['annual_supply_tons']) + results['annual_energy_pj'].append(transport['annual_energy_pj']) + results['capacity_ratio'].append(demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100) + + # 绘图 + fig, axes = plt.subplots(2, 2, figsize=(14, 12)) + + # 图1: 日补充量 + ax = axes[0, 0] + ax.plot(results['recycle_rate'], results['daily_supply'], 'b-', linewidth=2, marker='o', markersize=4) + ax.axvline(x=90, color='r', linestyle='--', label='Baseline (90%)') + ax.set_xlabel('Recycle Rate (%)', fontsize=12) + ax.set_ylabel('Daily Supply (tons)', fontsize=12) + ax.set_title('Daily Water Supply vs Recycle Rate', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + # 图2: 年补充量 + ax = axes[0, 1] + ax.plot(results['recycle_rate'], results['annual_supply'], 'g-', linewidth=2, marker='s', markersize=4) + ax.axvline(x=90, color='r', linestyle='--', label='Baseline (90%)') + ax.set_xlabel('Recycle Rate (%)', fontsize=12) + ax.set_ylabel('Annual Supply (tons)', fontsize=12) + ax.set_title('Annual Water Supply vs Recycle Rate', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + # 图3: 年能耗 + ax = axes[1, 0] + ax.plot(results['recycle_rate'], results['annual_energy_pj'], 'purple', linewidth=2, marker='^', markersize=4) + ax.axvline(x=90, color='r', linestyle='--', label='Baseline (90%)') + ax.set_xlabel('Recycle Rate (%)', fontsize=12) + ax.set_ylabel('Annual Energy (PJ)', fontsize=12) + ax.set_title('Annual Transport Energy vs Recycle Rate', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + # 图4: 运力占比 + ax = axes[1, 1] + ax.plot(results['recycle_rate'], results['capacity_ratio'], 'orange', linewidth=2, marker='d', markersize=4) + ax.axvline(x=90, color='r', linestyle='--', label='Baseline (90%)') + ax.axhline(y=100, color='gray', linestyle=':', label='Full Capacity') + ax.set_xlabel('Recycle Rate (%)', fontsize=12) + ax.set_ylabel('Elevator Capacity Usage (%)', fontsize=12) + ax.set_title('Elevator Capacity Utilization vs Recycle Rate', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + plt.suptitle(f'Sensitivity Analysis: Water Recycle Rate\n(α={baseline.comfort_factor}, N={baseline.population:,})', + fontsize=14, y=1.02) + plt.tight_layout() + plt.savefig(f'{save_dir}/sensitivity_recycle_rate.png', dpi=150, bbox_inches='tight') + print(f" 保存: {save_dir}/sensitivity_recycle_rate.png") + + return pd.DataFrame(results) + + +def sensitivity_comfort_factor(baseline: BaselineParameters, save_dir: str): + """ + 舒适度因子连续灵敏度分析 + α: 1 - 300 + """ + print("分析舒适度因子灵敏度...") + + alphas = np.concatenate([ + np.linspace(1, 10, 10), + np.linspace(15, 50, 8), + np.linspace(60, 150, 10), + np.linspace(175, 300, 6) + ]) + + results = { + 'alpha': [], + 'daily_per_person': [], + 'daily_supply': [], + 'annual_supply': [], + 'annual_energy_pj': [], + 'initial_transport': [], + 'capacity_ratio': [], + } + + for alpha in alphas: + demand = calculate_water_demand( + population=baseline.population, + survival_water=baseline.survival_water, + hygiene_water_base=baseline.hygiene_water_base, + comfort_factor=alpha, + medical_water_per_event=baseline.medical_water_per_event, + sickness_rate=baseline.sickness_rate, + confidence_level=baseline.confidence_level, + recycle_rate=baseline.recycle_rate, + safety_buffer_days=baseline.safety_buffer_days + ) + + transport = calculate_transport_metrics( + demand, + baseline.daily_elevator_capacity, + baseline.elevator_specific_energy + ) + + results['alpha'].append(alpha) + results['daily_per_person'].append(demand['daily_per_person_liters']) + results['daily_supply'].append(demand['daily_supply_tons']) + results['annual_supply'].append(demand['annual_supply_tons']) + results['annual_energy_pj'].append(transport['annual_energy_pj']) + results['initial_transport'].append(demand['initial_transport_tons']) + results['capacity_ratio'].append(demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100) + + # 绘图 + fig, axes = plt.subplots(2, 3, figsize=(18, 12)) + + # 图1: 人均日用水 + ax = axes[0, 0] + ax.plot(results['alpha'], results['daily_per_person'], 'b-', linewidth=2) + ax.axvline(x=1, color='green', linestyle='--', alpha=0.7, label='Survival (α=1)') + ax.axvline(x=50, color='orange', linestyle='--', alpha=0.7, label='Comfort (α=50)') + ax.axvline(x=250, color='red', linestyle='--', alpha=0.7, label='Luxury (α=250)') + ax.set_xlabel('Comfort Factor (α)', fontsize=12) + ax.set_ylabel('Daily Water per Person (L)', fontsize=12) + ax.set_title('Per Capita Water Usage vs Comfort Factor', fontsize=13) + ax.legend(fontsize=9) + ax.grid(True, alpha=0.3) + + # 图2: 日补充量 + ax = axes[0, 1] + ax.plot(results['alpha'], results['daily_supply'], 'g-', linewidth=2) + ax.axvline(x=50, color='r', linestyle='--', label='Baseline (α=50)') + ax.set_xlabel('Comfort Factor (α)', fontsize=12) + ax.set_ylabel('Daily Supply (tons)', fontsize=12) + ax.set_title('Daily Water Supply vs Comfort Factor', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + # 图3: 年补充量 + ax = axes[0, 2] + ax.plot(results['alpha'], results['annual_supply'], 'purple', linewidth=2) + ax.axvline(x=50, color='r', linestyle='--', label='Baseline (α=50)') + ax.set_xlabel('Comfort Factor (α)', fontsize=12) + ax.set_ylabel('Annual Supply (tons)', fontsize=12) + ax.set_title('Annual Water Supply vs Comfort Factor', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + # 图4: 年能耗 + ax = axes[1, 0] + ax.plot(results['alpha'], results['annual_energy_pj'], 'orange', linewidth=2) + ax.axvline(x=50, color='r', linestyle='--', label='Baseline (α=50)') + ax.set_xlabel('Comfort Factor (α)', fontsize=12) + ax.set_ylabel('Annual Energy (PJ)', fontsize=12) + ax.set_title('Annual Transport Energy vs Comfort Factor', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + # 图5: 首批运输量 + ax = axes[1, 1] + ax.plot(results['alpha'], results['initial_transport'], 'teal', linewidth=2) + ax.axvline(x=50, color='r', linestyle='--', label='Baseline (α=50)') + ax.set_xlabel('Comfort Factor (α)', fontsize=12) + ax.set_ylabel('Initial Transport (tons)', fontsize=12) + ax.set_title('Initial Transport Volume vs Comfort Factor', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + # 图6: 运力占比 + ax = axes[1, 2] + ax.plot(results['alpha'], results['capacity_ratio'], 'brown', linewidth=2) + ax.axvline(x=50, color='r', linestyle='--', label='Baseline (α=50)') + ax.axhline(y=100, color='gray', linestyle=':', label='Full Capacity') + ax.fill_between(results['alpha'], 0, results['capacity_ratio'], alpha=0.3) + ax.set_xlabel('Comfort Factor (α)', fontsize=12) + ax.set_ylabel('Elevator Capacity Usage (%)', fontsize=12) + ax.set_title('Elevator Capacity Utilization vs Comfort Factor', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + plt.suptitle(f'Sensitivity Analysis: Comfort Factor (α)\n(η={baseline.recycle_rate*100:.0f}%, N={baseline.population:,})', + fontsize=14, y=1.02) + plt.tight_layout() + plt.savefig(f'{save_dir}/sensitivity_comfort_factor.png', dpi=150, bbox_inches='tight') + print(f" 保存: {save_dir}/sensitivity_comfort_factor.png") + + return pd.DataFrame(results) + + +def sensitivity_medical_params(baseline: BaselineParameters, save_dir: str): + """ + 医疗紧急参数灵敏度分析 + - 患病率: 1% - 5% + - 每次医疗用水: 3 - 15 L + - 置信水平: 95% - 99.9% + """ + print("分析医疗参数灵敏度...") + + fig, axes = plt.subplots(2, 3, figsize=(18, 12)) + + # ========== 患病率分析 ========== + sickness_rates = np.linspace(0.01, 0.05, 21) + results_sr = {'sickness_rate': [], 'daily_medical': [], 'annual_supply': [], 'safety_buffer': []} + + for sr in sickness_rates: + demand = calculate_water_demand( + population=baseline.population, + survival_water=baseline.survival_water, + hygiene_water_base=baseline.hygiene_water_base, + comfort_factor=baseline.comfort_factor, + medical_water_per_event=baseline.medical_water_per_event, + sickness_rate=sr, + confidence_level=baseline.confidence_level, + recycle_rate=baseline.recycle_rate, + safety_buffer_days=baseline.safety_buffer_days + ) + results_sr['sickness_rate'].append(sr * 100) + results_sr['daily_medical'].append(demand['daily_medical_tons']) + results_sr['annual_supply'].append(demand['annual_supply_tons']) + results_sr['safety_buffer'].append(demand['safety_buffer_tons']) + + ax = axes[0, 0] + ax.plot(results_sr['sickness_rate'], results_sr['daily_medical'], 'b-', linewidth=2, marker='o', markersize=4) + ax.axvline(x=2, color='r', linestyle='--', label='Baseline (2%)') + ax.set_xlabel('Daily Sickness Rate (%)', fontsize=12) + ax.set_ylabel('Peak Daily Medical Water (tons)', fontsize=12) + ax.set_title('Peak Medical Water Demand vs Sickness Rate', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + ax = axes[0, 1] + ax.plot(results_sr['sickness_rate'], results_sr['safety_buffer'], 'g-', linewidth=2, marker='s', markersize=4) + ax.axvline(x=2, color='r', linestyle='--', label='Baseline (2%)') + ax.set_xlabel('Daily Sickness Rate (%)', fontsize=12) + ax.set_ylabel('Safety Buffer (tons)', fontsize=12) + ax.set_title('Safety Buffer vs Sickness Rate', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + # ========== 每次医疗用水分析 ========== + medical_waters = np.linspace(3, 15, 13) + results_mw = {'medical_water': [], 'daily_medical': [], 'annual_supply': []} + + for mw in medical_waters: + demand = calculate_water_demand( + population=baseline.population, + survival_water=baseline.survival_water, + hygiene_water_base=baseline.hygiene_water_base, + comfort_factor=baseline.comfort_factor, + medical_water_per_event=mw, + sickness_rate=baseline.sickness_rate, + confidence_level=baseline.confidence_level, + recycle_rate=baseline.recycle_rate, + safety_buffer_days=baseline.safety_buffer_days + ) + results_mw['medical_water'].append(mw) + results_mw['daily_medical'].append(demand['daily_medical_tons']) + results_mw['annual_supply'].append(demand['annual_supply_tons']) + + ax = axes[0, 2] + ax.plot(results_mw['medical_water'], results_mw['daily_medical'], 'purple', linewidth=2, marker='^', markersize=5) + ax.axvline(x=5, color='r', linestyle='--', label='Baseline (5 L)') + ax.set_xlabel('Medical Water per Event (L)', fontsize=12) + ax.set_ylabel('Peak Daily Medical Water (tons)', fontsize=12) + ax.set_title('Peak Medical Water vs Per-Event Usage', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + # ========== 置信水平分析 ========== + confidence_levels = np.linspace(0.90, 0.999, 20) + results_cl = {'confidence': [], 'daily_medical': [], 'safety_buffer': [], 'z_value': []} + + for cl in confidence_levels: + demand = calculate_water_demand( + population=baseline.population, + survival_water=baseline.survival_water, + hygiene_water_base=baseline.hygiene_water_base, + comfort_factor=baseline.comfort_factor, + medical_water_per_event=baseline.medical_water_per_event, + sickness_rate=baseline.sickness_rate, + confidence_level=cl, + recycle_rate=baseline.recycle_rate, + safety_buffer_days=baseline.safety_buffer_days + ) + results_cl['confidence'].append(cl * 100) + results_cl['daily_medical'].append(demand['daily_medical_tons']) + results_cl['safety_buffer'].append(demand['safety_buffer_tons']) + results_cl['z_value'].append(stats.norm.ppf(cl)) + + ax = axes[1, 0] + ax.plot(results_cl['confidence'], results_cl['daily_medical'], 'orange', linewidth=2, marker='d', markersize=4) + ax.axvline(x=99, color='r', linestyle='--', label='Baseline (99%)') + ax.set_xlabel('Confidence Level (%)', fontsize=12) + ax.set_ylabel('Peak Daily Medical Water (tons)', fontsize=12) + ax.set_title('Peak Medical Water vs Confidence Level', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + ax = axes[1, 1] + ax.plot(results_cl['confidence'], results_cl['z_value'], 'teal', linewidth=2, marker='o', markersize=4) + ax.axvline(x=99, color='r', linestyle='--', label='Baseline (99%)') + ax.set_xlabel('Confidence Level (%)', fontsize=12) + ax.set_ylabel('Z-Value', fontsize=12) + ax.set_title('Z-Value vs Confidence Level', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + # ========== 综合影响图 ========== + ax = axes[1, 2] + + # 计算各参数变化对年供水的影响 + baseline_demand = calculate_water_demand( + population=baseline.population, + survival_water=baseline.survival_water, + hygiene_water_base=baseline.hygiene_water_base, + comfort_factor=baseline.comfort_factor, + medical_water_per_event=baseline.medical_water_per_event, + sickness_rate=baseline.sickness_rate, + confidence_level=baseline.confidence_level, + recycle_rate=baseline.recycle_rate, + safety_buffer_days=baseline.safety_buffer_days + ) + baseline_annual = baseline_demand['annual_supply_tons'] + + # 患病率影响 (1% -> 5%) + sr_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, baseline.medical_water_per_event, 0.01, + baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) + sr_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, baseline.medical_water_per_event, 0.05, + baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) + + # 医疗用水影响 (3L -> 15L) + mw_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, 3, baseline.sickness_rate, + baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) + mw_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, 15, baseline.sickness_rate, + baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) + + # 置信水平影响 (95% -> 99.9%) + cl_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, + 0.95, baseline.recycle_rate, baseline.safety_buffer_days) + cl_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, + 0.999, baseline.recycle_rate, baseline.safety_buffer_days) + + params = ['Sickness Rate\n(1%→5%)', 'Medical Water\n(3L→15L)', 'Confidence\n(95%→99.9%)'] + changes_low = [ + (sr_low['annual_supply_tons'] - baseline_annual) / baseline_annual * 100, + (mw_low['annual_supply_tons'] - baseline_annual) / baseline_annual * 100, + (cl_low['annual_supply_tons'] - baseline_annual) / baseline_annual * 100, + ] + changes_high = [ + (sr_high['annual_supply_tons'] - baseline_annual) / baseline_annual * 100, + (mw_high['annual_supply_tons'] - baseline_annual) / baseline_annual * 100, + (cl_high['annual_supply_tons'] - baseline_annual) / baseline_annual * 100, + ] + + x = np.arange(len(params)) + width = 0.35 + ax.barh(x - width/2, changes_low, width, label='Low Value', color='green', alpha=0.7) + ax.barh(x + width/2, changes_high, width, label='High Value', color='red', alpha=0.7) + ax.axvline(x=0, color='black', linewidth=1) + ax.set_yticks(x) + ax.set_yticklabels(params) + ax.set_xlabel('Change in Annual Supply (%)', fontsize=12) + ax.set_title('Medical Parameters Impact on Annual Supply', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3, axis='x') + + plt.suptitle(f'Sensitivity Analysis: Medical Emergency Parameters\n(α={baseline.comfort_factor}, η={baseline.recycle_rate*100:.0f}%)', + fontsize=14, y=1.02) + plt.tight_layout() + plt.savefig(f'{save_dir}/sensitivity_medical_params.png', dpi=150, bbox_inches='tight') + print(f" 保存: {save_dir}/sensitivity_medical_params.png") + + return results_sr, results_mw, results_cl + + +def sensitivity_population(baseline: BaselineParameters, save_dir: str): + """ + 人口规模灵敏度分析 + N: 80,000 - 120,000 (±20%) + """ + print("分析人口规模灵敏度...") + + populations = np.linspace(80000, 120000, 21).astype(int) + + results = { + 'population': [], + 'daily_circulation': [], + 'daily_supply': [], + 'annual_supply': [], + 'annual_energy_pj': [], + 'initial_transport': [], + 'capacity_ratio': [], + } + + for N in populations: + demand = calculate_water_demand( + population=N, + survival_water=baseline.survival_water, + hygiene_water_base=baseline.hygiene_water_base, + comfort_factor=baseline.comfort_factor, + medical_water_per_event=baseline.medical_water_per_event, + sickness_rate=baseline.sickness_rate, + confidence_level=baseline.confidence_level, + recycle_rate=baseline.recycle_rate, + safety_buffer_days=baseline.safety_buffer_days + ) + + transport = calculate_transport_metrics( + demand, + baseline.daily_elevator_capacity, + baseline.elevator_specific_energy + ) + + results['population'].append(N / 1000) # 千人 + results['daily_circulation'].append(demand['daily_circulation_tons']) + results['daily_supply'].append(demand['daily_supply_tons']) + results['annual_supply'].append(demand['annual_supply_tons']) + results['annual_energy_pj'].append(transport['annual_energy_pj']) + results['initial_transport'].append(demand['initial_transport_tons']) + results['capacity_ratio'].append(demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100) + + # 绘图 + fig, axes = plt.subplots(2, 2, figsize=(14, 12)) + + # 图1: 日循环用水 + ax = axes[0, 0] + ax.plot(results['population'], results['daily_circulation'], 'b-', linewidth=2, marker='o', markersize=4) + ax.axvline(x=100, color='r', linestyle='--', label='Baseline (100k)') + ax.fill_between(results['population'], 0, results['daily_circulation'], alpha=0.2) + ax.set_xlabel('Population (thousands)', fontsize=12) + ax.set_ylabel('Daily Circulation Water (tons)', fontsize=12) + ax.set_title('Daily Circulation Water vs Population', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + # 图2: 日补充量 + ax = axes[0, 1] + ax.plot(results['population'], results['daily_supply'], 'g-', linewidth=2, marker='s', markersize=4) + ax.axvline(x=100, color='r', linestyle='--', label='Baseline (100k)') + ax.set_xlabel('Population (thousands)', fontsize=12) + ax.set_ylabel('Daily Supply (tons)', fontsize=12) + ax.set_title('Daily Water Supply vs Population', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + # 图3: 年能耗 + ax = axes[1, 0] + ax.plot(results['population'], results['annual_energy_pj'], 'purple', linewidth=2, marker='^', markersize=4) + ax.axvline(x=100, color='r', linestyle='--', label='Baseline (100k)') + ax.set_xlabel('Population (thousands)', fontsize=12) + ax.set_ylabel('Annual Energy (PJ)', fontsize=12) + ax.set_title('Annual Transport Energy vs Population', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + # 图4: 运力占比 + ax = axes[1, 1] + ax.plot(results['population'], results['capacity_ratio'], 'orange', linewidth=2, marker='d', markersize=4) + ax.axvline(x=100, color='r', linestyle='--', label='Baseline (100k)') + ax.axhline(y=100, color='gray', linestyle=':', label='Full Capacity') + ax.set_xlabel('Population (thousands)', fontsize=12) + ax.set_ylabel('Elevator Capacity Usage (%)', fontsize=12) + ax.set_title('Elevator Capacity Utilization vs Population', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + plt.suptitle(f'Sensitivity Analysis: Population Size\n(α={baseline.comfort_factor}, η={baseline.recycle_rate*100:.0f}%)', + fontsize=14, y=1.02) + plt.tight_layout() + plt.savefig(f'{save_dir}/sensitivity_population.png', dpi=150, bbox_inches='tight') + print(f" 保存: {save_dir}/sensitivity_population.png") + + return pd.DataFrame(results) + + +def sensitivity_buffer_days(baseline: BaselineParameters, save_dir: str): + """ + 安全缓冲天数灵敏度分析 + Days: 15 - 60 + """ + print("分析安全缓冲天数灵敏度...") + + buffer_days = np.arange(15, 61, 3) + + results = { + 'buffer_days': [], + 'safety_buffer_tons': [], + 'initial_transport': [], + 'initial_days': [], + 'initial_energy_pj': [], + } + + for days in buffer_days: + demand = calculate_water_demand( + population=baseline.population, + survival_water=baseline.survival_water, + hygiene_water_base=baseline.hygiene_water_base, + comfort_factor=baseline.comfort_factor, + medical_water_per_event=baseline.medical_water_per_event, + sickness_rate=baseline.sickness_rate, + confidence_level=baseline.confidence_level, + recycle_rate=baseline.recycle_rate, + safety_buffer_days=days + ) + + transport = calculate_transport_metrics( + demand, + baseline.daily_elevator_capacity, + baseline.elevator_specific_energy + ) + + results['buffer_days'].append(days) + results['safety_buffer_tons'].append(demand['safety_buffer_tons']) + results['initial_transport'].append(demand['initial_transport_tons']) + results['initial_days'].append(transport['initial_days']) + results['initial_energy_pj'].append(transport['initial_energy_pj']) + + # 绘图 + fig, axes = plt.subplots(2, 2, figsize=(14, 12)) + + # 图1: 安全缓冲量 + ax = axes[0, 0] + ax.plot(results['buffer_days'], results['safety_buffer_tons'], 'b-', linewidth=2, marker='o', markersize=5) + ax.axvline(x=30, color='r', linestyle='--', label='Baseline (30 days)') + ax.set_xlabel('Safety Buffer Days', fontsize=12) + ax.set_ylabel('Safety Buffer (tons)', fontsize=12) + ax.set_title('Safety Buffer Volume vs Buffer Days', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + # 图2: 首批运输量 + ax = axes[0, 1] + ax.plot(results['buffer_days'], results['initial_transport'], 'g-', linewidth=2, marker='s', markersize=5) + ax.axvline(x=30, color='r', linestyle='--', label='Baseline (30 days)') + ax.set_xlabel('Safety Buffer Days', fontsize=12) + ax.set_ylabel('Initial Transport (tons)', fontsize=12) + ax.set_title('Initial Transport Volume vs Buffer Days', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + # 图3: 首批运输天数 + ax = axes[1, 0] + ax.plot(results['buffer_days'], results['initial_days'], 'purple', linewidth=2, marker='^', markersize=5) + ax.axvline(x=30, color='r', linestyle='--', label='Baseline (30 days)') + ax.set_xlabel('Safety Buffer Days', fontsize=12) + ax.set_ylabel('Initial Transport Time (days)', fontsize=12) + ax.set_title('Initial Transport Duration vs Buffer Days', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + # 图4: 首批运输能耗 + ax = axes[1, 1] + ax.plot(results['buffer_days'], results['initial_energy_pj'], 'orange', linewidth=2, marker='d', markersize=5) + ax.axvline(x=30, color='r', linestyle='--', label='Baseline (30 days)') + ax.set_xlabel('Safety Buffer Days', fontsize=12) + ax.set_ylabel('Initial Transport Energy (PJ)', fontsize=12) + ax.set_title('Initial Transport Energy vs Buffer Days', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3) + + plt.suptitle(f'Sensitivity Analysis: Safety Buffer Days\n(α={baseline.comfort_factor}, η={baseline.recycle_rate*100:.0f}%)', + fontsize=14, y=1.02) + plt.tight_layout() + plt.savefig(f'{save_dir}/sensitivity_buffer_days.png', dpi=150, bbox_inches='tight') + print(f" 保存: {save_dir}/sensitivity_buffer_days.png") + + return pd.DataFrame(results) + + +# ============== 多参数交互分析 ============== + +def tornado_analysis(baseline: BaselineParameters, save_dir: str): + """ + Tornado图 - 多参数灵敏度排序分析 + """ + print("生成Tornado图...") + + # 基准值 + baseline_demand = calculate_water_demand( + population=baseline.population, + survival_water=baseline.survival_water, + hygiene_water_base=baseline.hygiene_water_base, + comfort_factor=baseline.comfort_factor, + medical_water_per_event=baseline.medical_water_per_event, + sickness_rate=baseline.sickness_rate, + confidence_level=baseline.confidence_level, + recycle_rate=baseline.recycle_rate, + safety_buffer_days=baseline.safety_buffer_days + ) + baseline_annual = baseline_demand['annual_supply_tons'] + baseline_initial = baseline_demand['initial_transport_tons'] + + # 定义参数范围 + params = [ + ('Recycle Rate (η)', 'recycle_rate', 0.80, 0.95, baseline.recycle_rate), + ('Comfort Factor (α)', 'comfort_factor', 25, 100, baseline.comfort_factor), + ('Population (N)', 'population', 80000, 120000, baseline.population), + ('Sickness Rate', 'sickness_rate', 0.01, 0.04, baseline.sickness_rate), + ('Medical Water/Event', 'medical_water', 3, 10, baseline.medical_water_per_event), + ('Buffer Days', 'buffer_days', 15, 45, baseline.safety_buffer_days), + ('Confidence Level', 'confidence', 0.95, 0.999, baseline.confidence_level), + ] + + results_annual = [] + results_initial = [] + + for name, param_type, low, high, base_val in params: + # 低值计算 + if param_type == 'recycle_rate': + demand_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, + baseline.confidence_level, low, baseline.safety_buffer_days) + demand_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, + baseline.confidence_level, high, baseline.safety_buffer_days) + elif param_type == 'comfort_factor': + demand_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + low, baseline.medical_water_per_event, baseline.sickness_rate, + baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) + demand_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + high, baseline.medical_water_per_event, baseline.sickness_rate, + baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) + elif param_type == 'population': + demand_low = calculate_water_demand(int(low), baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, + baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) + demand_high = calculate_water_demand(int(high), baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, + baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) + elif param_type == 'sickness_rate': + demand_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, baseline.medical_water_per_event, low, + baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) + demand_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, baseline.medical_water_per_event, high, + baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) + elif param_type == 'medical_water': + demand_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, low, baseline.sickness_rate, + baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) + demand_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, high, baseline.sickness_rate, + baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days) + elif param_type == 'buffer_days': + demand_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, + baseline.confidence_level, baseline.recycle_rate, int(low)) + demand_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, + baseline.confidence_level, baseline.recycle_rate, int(high)) + elif param_type == 'confidence': + demand_low = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, + low, baseline.recycle_rate, baseline.safety_buffer_days) + demand_high = calculate_water_demand(baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, + high, baseline.recycle_rate, baseline.safety_buffer_days) + + # 计算变化百分比 + change_low_annual = (demand_low['annual_supply_tons'] - baseline_annual) / baseline_annual * 100 + change_high_annual = (demand_high['annual_supply_tons'] - baseline_annual) / baseline_annual * 100 + change_low_initial = (demand_low['initial_transport_tons'] - baseline_initial) / baseline_initial * 100 + change_high_initial = (demand_high['initial_transport_tons'] - baseline_initial) / baseline_initial * 100 + + # 确保排序一致性 (swing = |high - low|) + swing_annual = abs(change_high_annual - change_low_annual) + swing_initial = abs(change_high_initial - change_low_initial) + + results_annual.append({ + 'param': name, + 'low': min(change_low_annual, change_high_annual), + 'high': max(change_low_annual, change_high_annual), + 'swing': swing_annual + }) + results_initial.append({ + 'param': name, + 'low': min(change_low_initial, change_high_initial), + 'high': max(change_low_initial, change_high_initial), + 'swing': swing_initial + }) + + # 按swing排序 + results_annual = sorted(results_annual, key=lambda x: x['swing'], reverse=True) + results_initial = sorted(results_initial, key=lambda x: x['swing'], reverse=True) + + # 绘图 + fig, axes = plt.subplots(1, 2, figsize=(16, 8)) + + # 图1: 年补充量Tornado + ax = axes[0] + y_pos = np.arange(len(results_annual)) + params_sorted = [r['param'] for r in results_annual] + lows = [r['low'] for r in results_annual] + highs = [r['high'] for r in results_annual] + + # 绘制条形图 + for i, (l, h) in enumerate(zip(lows, highs)): + if l < 0: + ax.barh(i, l, color='green', alpha=0.7, height=0.6) + if h > 0: + ax.barh(i, h, color='red', alpha=0.7, height=0.6) + if l >= 0: + ax.barh(i, l, left=0, color='green', alpha=0.7, height=0.6) + if h <= 0: + ax.barh(i, h, left=0, color='red', alpha=0.7, height=0.6) + # 如果跨越0，绘制完整范围 + if l < 0 < h: + pass # 已经绘制 + elif l >= 0: + ax.barh(i, h - l, left=l, color='orange', alpha=0.7, height=0.6) + elif h <= 0: + ax.barh(i, h - l, left=l, color='orange', alpha=0.7, height=0.6) + + ax.axvline(x=0, color='black', linewidth=1.5) + ax.set_yticks(y_pos) + ax.set_yticklabels(params_sorted) + ax.set_xlabel('Change in Annual Supply (%)', fontsize=12) + ax.set_title('Tornado Diagram: Annual Water Supply', fontsize=14) + ax.grid(True, alpha=0.3, axis='x') + + # 图2: 首批运输量Tornado + ax = axes[1] + params_sorted = [r['param'] for r in results_initial] + lows = [r['low'] for r in results_initial] + highs = [r['high'] for r in results_initial] + + for i, (l, h) in enumerate(zip(lows, highs)): + width = h - l + left = l + color = 'steelblue' if width > 0 else 'gray' + ax.barh(i, width, left=left, color=color, alpha=0.7, height=0.6) + + ax.axvline(x=0, color='black', linewidth=1.5) + ax.set_yticks(y_pos) + ax.set_yticklabels(params_sorted) + ax.set_xlabel('Change in Initial Transport (%)', fontsize=12) + ax.set_title('Tornado Diagram: Initial Transport Volume', fontsize=14) + ax.grid(True, alpha=0.3, axis='x') + + # 添加图例说明 + from matplotlib.patches import Patch + legend_elements = [ + Patch(facecolor='green', alpha=0.7, label='Decrease'), + Patch(facecolor='red', alpha=0.7, label='Increase'), + ] + axes[0].legend(handles=legend_elements, loc='lower right') + + plt.suptitle('Multi-Parameter Sensitivity Analysis (Tornado Diagram)\nBaseline: α=50, η=90%, N=100k', + fontsize=14, y=1.02) + plt.tight_layout() + plt.savefig(f'{save_dir}/tornado_analysis.png', dpi=150, bbox_inches='tight') + print(f" 保存: {save_dir}/tornado_analysis.png") + + return results_annual, results_initial + + +def interaction_heatmap(baseline: BaselineParameters, save_dir: str): + """ + 双参数交互热力图 - α vs η + """ + print("生成双参数交互热力图...") + + alphas = np.array([1, 25, 50, 100, 150, 200, 250]) + etas = np.array([0.70, 0.75, 0.80, 0.85, 0.90, 0.95]) + + # 年补充量矩阵 + annual_matrix = np.zeros((len(etas), len(alphas))) + # 运力占比矩阵 + capacity_matrix = np.zeros((len(etas), len(alphas))) + + for i, eta in enumerate(etas): + for j, alpha in enumerate(alphas): + demand = calculate_water_demand( + population=baseline.population, + survival_water=baseline.survival_water, + hygiene_water_base=baseline.hygiene_water_base, + comfort_factor=alpha, + medical_water_per_event=baseline.medical_water_per_event, + sickness_rate=baseline.sickness_rate, + confidence_level=baseline.confidence_level, + recycle_rate=eta, + safety_buffer_days=baseline.safety_buffer_days + ) + annual_matrix[i, j] = demand['annual_supply_tons'] / 1000 # 千吨 + capacity_matrix[i, j] = demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100 + + fig, axes = plt.subplots(1, 2, figsize=(16, 7)) + + # 图1: 年补充量热力图 + ax = axes[0] + im = ax.imshow(annual_matrix, cmap='YlOrRd', aspect='auto') + ax.set_xticks(range(len(alphas))) + ax.set_xticklabels([f'{a}' for a in alphas]) + ax.set_yticks(range(len(etas))) + ax.set_yticklabels([f'{e*100:.0f}%' for e in etas]) + ax.set_xlabel('Comfort Factor (α)', fontsize=12) + ax.set_ylabel('Recycle Rate (η)', fontsize=12) + ax.set_title('Annual Water Supply (kt)', fontsize=13) + + # 添加数值标注 + for i in range(len(etas)): + for j in range(len(alphas)): + text_color = 'white' if annual_matrix[i, j] > annual_matrix.max() * 0.6 else 'black' + ax.text(j, i, f'{annual_matrix[i, j]:.0f}', ha='center', va='center', fontsize=9, color=text_color) + + plt.colorbar(im, ax=ax, label='Annual Supply (kt)') + + # 图2: 运力占比热力图 + ax = axes[1] + im = ax.imshow(capacity_matrix, cmap='RdYlGn_r', aspect='auto', vmin=0, vmax=150) + ax.set_xticks(range(len(alphas))) + ax.set_xticklabels([f'{a}' for a in alphas]) + ax.set_yticks(range(len(etas))) + ax.set_yticklabels([f'{e*100:.0f}%' for e in etas]) + ax.set_xlabel('Comfort Factor (α)', fontsize=12) + ax.set_ylabel('Recycle Rate (η)', fontsize=12) + ax.set_title('Elevator Capacity Usage (%)', fontsize=13) + + # 添加数值标注 + for i in range(len(etas)): + for j in range(len(alphas)): + val = capacity_matrix[i, j] + text_color = 'white' if val > 60 else 'black' + ax.text(j, i, f'{val:.0f}%', ha='center', va='center', fontsize=9, color=text_color) + + # 添加100%等高线 + cbar = plt.colorbar(im, ax=ax, label='Capacity Usage (%)') + + plt.suptitle('Two-Parameter Interaction Analysis: α vs η\n(N=100,000, Buffer=30 days)', + fontsize=14, y=1.02) + plt.tight_layout() + plt.savefig(f'{save_dir}/interaction_heatmap.png', dpi=150, bbox_inches='tight') + print(f" 保存: {save_dir}/interaction_heatmap.png") + + return annual_matrix, capacity_matrix + + +def three_param_surface(baseline: BaselineParameters, save_dir: str): + """ + 三参数交互3D曲面图 + """ + print("生成三参数3D曲面图...") + + alphas = np.linspace(1, 200, 20) + etas = np.linspace(0.75, 0.95, 20) + + A, E = np.meshgrid(alphas, etas) + Z = np.zeros_like(A) + + for i in range(len(etas)): + for j in range(len(alphas)): + demand = calculate_water_demand( + population=baseline.population, + survival_water=baseline.survival_water, + hygiene_water_base=baseline.hygiene_water_base, + comfort_factor=alphas[j], + medical_water_per_event=baseline.medical_water_per_event, + sickness_rate=baseline.sickness_rate, + confidence_level=baseline.confidence_level, + recycle_rate=etas[i], + safety_buffer_days=baseline.safety_buffer_days + ) + Z[i, j] = demand['annual_supply_tons'] / 1000 # 千吨 + + fig = plt.figure(figsize=(14, 10)) + ax = fig.add_subplot(111, projection='3d') + + surf = ax.plot_surface(A, E * 100, Z, cmap='viridis', alpha=0.8, edgecolor='none') + + # 添加基准点 + baseline_demand = calculate_water_demand( + baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, + baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days + ) + ax.scatter([baseline.comfort_factor], [baseline.recycle_rate * 100], + [baseline_demand['annual_supply_tons'] / 1000], + color='red', s=100, label='Baseline', marker='*') + + ax.set_xlabel('Comfort Factor (α)', fontsize=11) + ax.set_ylabel('Recycle Rate (%)', fontsize=11) + ax.set_zlabel('Annual Supply (kt)', fontsize=11) + ax.set_title('3D Surface: Annual Water Supply\nas Function of α and η', fontsize=14) + + fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, label='Annual Supply (kt)') + ax.legend() + + plt.tight_layout() + plt.savefig(f'{save_dir}/three_param_surface.png', dpi=150, bbox_inches='tight') + print(f" 保存: {save_dir}/three_param_surface.png") + + return Z + + +# ============== 最坏情况分析 ============== + +def worst_case_analysis(baseline: BaselineParameters, save_dir: str): + """ + 最坏情况压力测试分析 + """ + print("进行最坏情况分析...") + + # 定义场景 + scenarios = { + 'Baseline': { + 'alpha': baseline.comfort_factor, + 'eta': baseline.recycle_rate, + 'sickness': baseline.sickness_rate, + 'population': baseline.population, + 'buffer': baseline.safety_buffer_days, + 'medical_water': baseline.medical_water_per_event, + }, + 'Survival Mode': { + 'alpha': 1, + 'eta': 0.90, + 'sickness': 0.02, + 'population': 100000, + 'buffer': 30, + 'medical_water': 5, + }, + 'Comfort Mode': { + 'alpha': 50, + 'eta': 0.90, + 'sickness': 0.02, + 'population': 100000, + 'buffer': 30, + 'medical_water': 5, + }, + 'Luxury Mode': { + 'alpha': 250, + 'eta': 0.90, + 'sickness': 0.02, + 'population': 100000, + 'buffer': 30, + 'medical_water': 5, + }, + 'Recycle Degradation': { + 'alpha': 50, + 'eta': 0.80, + 'sickness': 0.02, + 'population': 100000, + 'buffer': 30, + 'medical_water': 5, + }, + 'Epidemic Outbreak': { + 'alpha': 50, + 'eta': 0.90, + 'sickness': 0.05, + 'population': 100000, + 'buffer': 45, + 'medical_water': 10, + }, + 'Population Surge': { + 'alpha': 50, + 'eta': 0.90, + 'sickness': 0.02, + 'population': 120000, + 'buffer': 30, + 'medical_water': 5, + }, + 'Worst Case': { + 'alpha': 250, + 'eta': 0.80, + 'sickness': 0.05, + 'population': 120000, + 'buffer': 45, + 'medical_water': 10, + }, + 'Moderate Stress': { + 'alpha': 100, + 'eta': 0.85, + 'sickness': 0.03, + 'population': 110000, + 'buffer': 40, + 'medical_water': 7, + }, + } + + results = [] + + for name, params in scenarios.items(): + demand = calculate_water_demand( + population=params['population'], + survival_water=baseline.survival_water, + hygiene_water_base=baseline.hygiene_water_base, + comfort_factor=params['alpha'], + medical_water_per_event=params['medical_water'], + sickness_rate=params['sickness'], + confidence_level=baseline.confidence_level, + recycle_rate=params['eta'], + safety_buffer_days=params['buffer'] + ) + + transport = calculate_transport_metrics( + demand, + baseline.daily_elevator_capacity, + baseline.elevator_specific_energy + ) + + capacity_ratio = demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100 + + results.append({ + 'scenario': name, + 'alpha': params['alpha'], + 'eta': params['eta'] * 100, + 'sickness': params['sickness'] * 100, + 'population': params['population'] / 1000, + 'buffer': params['buffer'], + 'daily_supply': demand['daily_supply_tons'], + 'annual_supply': demand['annual_supply_tons'], + 'initial_transport': demand['initial_transport_tons'], + 'annual_energy_pj': transport['annual_energy_pj'], + 'capacity_ratio': capacity_ratio, + 'feasible': capacity_ratio < 100, + }) + + df = pd.DataFrame(results) + + # 绘图 + fig, axes = plt.subplots(2, 2, figsize=(16, 14)) + + # 图1: 年补充量对比 + ax = axes[0, 0] + colors = ['green' if r['feasible'] else 'red' for r in results] + bars = ax.barh(range(len(results)), [r['annual_supply'] / 1000 for r in results], color=colors, alpha=0.7) + ax.axvline(x=baseline.total_elevator_capacity / 1000, color='red', linestyle='--', linewidth=2, label='Elevator Capacity') + ax.set_yticks(range(len(results))) + ax.set_yticklabels([r['scenario'] for r in results]) + ax.set_xlabel('Annual Supply (kt)', fontsize=12) + ax.set_title('Annual Water Supply by Scenario', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3, axis='x') + + # 添加数值标注 + for i, bar in enumerate(bars): + width = bar.get_width() + ax.text(width + 5, bar.get_y() + bar.get_height()/2, f'{width:.0f}', + va='center', fontsize=9) + + # 图2: 运力占比 + ax = axes[0, 1] + colors = ['green' if r['feasible'] else 'red' for r in results] + bars = ax.barh(range(len(results)), [r['capacity_ratio'] for r in results], color=colors, alpha=0.7) + ax.axvline(x=100, color='red', linestyle='--', linewidth=2, label='100% Capacity') + ax.set_yticks(range(len(results))) + ax.set_yticklabels([r['scenario'] for r in results]) + ax.set_xlabel('Elevator Capacity Usage (%)', fontsize=12) + ax.set_title('Capacity Utilization by Scenario', fontsize=13) + ax.legend() + ax.grid(True, alpha=0.3, axis='x') + + for i, bar in enumerate(bars): + width = bar.get_width() + ax.text(width + 2, bar.get_y() + bar.get_height()/2, f'{width:.1f}%', + va='center', fontsize=9) + + # 图3: 年能耗对比 + ax = axes[1, 0] + ax.bar(range(len(results)), [r['annual_energy_pj'] for r in results], color='purple', alpha=0.7) + ax.set_xticks(range(len(results))) + ax.set_xticklabels([r['scenario'] for r in results], rotation=45, ha='right', fontsize=9) + ax.set_ylabel('Annual Energy (PJ)', fontsize=12) + ax.set_title('Annual Transport Energy by Scenario', fontsize=13) + ax.grid(True, alpha=0.3, axis='y') + + # 图4: 首批运输量 + ax = axes[1, 1] + ax.bar(range(len(results)), [r['initial_transport'] / 1000 for r in results], color='teal', alpha=0.7) + ax.set_xticks(range(len(results))) + ax.set_xticklabels([r['scenario'] for r in results], rotation=45, ha='right', fontsize=9) + ax.set_ylabel('Initial Transport (kt)', fontsize=12) + ax.set_title('Initial Transport Volume by Scenario', fontsize=13) + ax.grid(True, alpha=0.3, axis='y') + + plt.suptitle('Worst Case Analysis: Stress Test Scenarios', fontsize=14, y=1.02) + plt.tight_layout() + plt.savefig(f'{save_dir}/worst_case_analysis.png', dpi=150, bbox_inches='tight') + print(f" 保存: {save_dir}/worst_case_analysis.png") + + # 保存详细结果表 + df.to_csv(f'{save_dir}/worst_case_results.csv', index=False) + print(f" 保存: {save_dir}/worst_case_results.csv") + + return df + + +def generate_summary_report(baseline: BaselineParameters, save_dir: str): + """ + 生成灵敏度分析总结报告 + """ + print("生成总结报告...") + + report = [] + report.append("=" * 100) + report.append("TASK 3: WATER SUPPLY SENSITIVITY ANALYSIS SUMMARY") + report.append("=" * 100) + + report.append("\n## 基准参数") + report.append(f"- 人口: {baseline.population:,}") + report.append(f"- 舒适度因子 (α): {baseline.comfort_factor}") + report.append(f"- 水回收效率 (η): {baseline.recycle_rate * 100:.0f}%") + report.append(f"- 患病率: {baseline.sickness_rate * 100:.1f}%") + report.append(f"- 医疗用水: {baseline.medical_water_per_event} L/次") + report.append(f"- 置信水平: {baseline.confidence_level * 100:.1f}%") + report.append(f"- 安全缓冲: {baseline.safety_buffer_days} 天") + + # 计算基准值 + baseline_demand = calculate_water_demand( + baseline.population, baseline.survival_water, baseline.hygiene_water_base, + baseline.comfort_factor, baseline.medical_water_per_event, baseline.sickness_rate, + baseline.confidence_level, baseline.recycle_rate, baseline.safety_buffer_days + ) + baseline_transport = calculate_transport_metrics( + baseline_demand, baseline.daily_elevator_capacity, baseline.elevator_specific_energy + ) + + report.append(f"\n## 基准水需求") + report.append(f"- 每日补充量: {baseline_demand['daily_supply_tons']:.1f} 吨") + report.append(f"- 年补充量: {baseline_demand['annual_supply_tons']:.1f} 吨") + report.append(f"- 首批运输量: {baseline_demand['initial_transport_tons']:.1f} 吨") + report.append(f"- 年能耗: {baseline_transport['annual_energy_pj']:.4f} PJ") + report.append(f"- 运力占比: {baseline_demand['annual_supply_tons'] / baseline.total_elevator_capacity * 100:.2f}%") + + report.append("\n" + "=" * 100) + report.append("## 关键灵敏度发现") + report.append("=" * 100) + + report.append(""" +### 1. 水回收效率 (η) - 影响程度: ★★★★★ +- 范围: 70% - 95% +- η从90%降至80%时，日补充量翻倍 +- 这是对年运输需求影响最大的参数 +- 建议: 优先投资水回收系统维护和升级 + +### 2. 舒适度因子 (α) - 影响程度: ★★★★★ +- 范围: 1 - 300 +- α=250时年补充量是α=1时的35倍 +- 在高舒适度下可能超出电梯运力 +- 建议: 初期采用生存标准，逐步提升舒适度 + +### 3. 人口规模 (N) - 影响程度: ★★★★☆ +- 范围: ±20% (80k - 120k) +- 线性影响，20%人口增加导致20%需求增加 +- 建议: 规划时预留弹性容量 + +### 4. 医疗紧急参数 - 影响程度: ★★☆☆☆ +- 患病率(1%-5%)和医疗用水(3-15L)对年补充量影响较小 +- 主要影响安全缓冲和首批运输量 +- 建议: 采用保守参数确保医疗安全 + +### 5. 安全缓冲天数 - 影响程度: ★★★☆☆ +- 范围: 15 - 60天 +- 主要影响首批运输量，不影响年运输总量 +- 建议: 保持30天缓冲作为运输中断容忍度 +""") + + report.append("\n" + "=" * 100) + report.append("## 最坏情况分析结论") + report.append("=" * 100) + + report.append(""" +### 可行性边界 +- 最坏情况 (α=250, η=80%, N=120k): 年需求约 562 kt，超出电梯运力 +- 中度压力 (α=100, η=85%, N=110k): 年需求约 132 kt，仍可承受 +- 系统崩溃临界点: 约 η < 82% + α > 200 组合 + +### 风险缓解建议 +1. 水回收系统冗余设计，确保η≥85% +2. 分阶段提升舒适度标准 +3. 维持混合运输能力作为备用 +4. 建立月球本地水源开发计划 +""") + + report_text = "\n".join(report) + + with open(f'{save_dir}/sensitivity_summary_report.txt', 'w', encoding='utf-8') as f: + f.write(report_text) + print(f" 保存: {save_dir}/sensitivity_summary_report.txt") + + return report_text + + +# ============== 主程序 ============== + +def run_full_sensitivity_analysis(): + """运行完整的灵敏度分析""" + + print("=" * 80) + print("任务三：月球殖民地水需求灵敏度分析") + print("=" * 80) + + save_dir = '/Volumes/Files/code/mm/20260130_b/p3' + + # 初始化基准参数 + baseline = BaselineParameters() + + print(f"\n基准参数:") + print(f" - 人口: {baseline.population:,}") + print(f" - 舒适度因子: {baseline.comfort_factor}") + print(f" - 水回收效率: {baseline.recycle_rate * 100:.0f}%") + print(f" - 患病率: {baseline.sickness_rate * 100:.1f}%") + print(f" - 安全缓冲: {baseline.safety_buffer_days} 天") + + print("\n" + "=" * 80) + print("开始灵敏度分析...") + print("=" * 80) + + # 1. 单参数灵敏度分析 + print("\n[1/7] 水回收效率灵敏度分析") + df_recycle = sensitivity_recycle_rate(baseline, save_dir) + + print("\n[2/7] 舒适度因子灵敏度分析") + df_comfort = sensitivity_comfort_factor(baseline, save_dir) + + print("\n[3/7] 医疗参数灵敏度分析") + results_medical = sensitivity_medical_params(baseline, save_dir) + + print("\n[4/7] 人口规模灵敏度分析") + df_population = sensitivity_population(baseline, save_dir) + + print("\n[5/7] 安全缓冲天数灵敏度分析") + df_buffer = sensitivity_buffer_days(baseline, save_dir) + + # 2. 多参数交互分析 + print("\n[6/7] 多参数交互分析") + tornado_results = tornado_analysis(baseline, save_dir) + heatmap_data = interaction_heatmap(baseline, save_dir) + surface_data = three_param_surface(baseline, save_dir) + + # 3. 最坏情况分析 + print("\n[7/7] 最坏情况分析") + worst_case_df = worst_case_analysis(baseline, save_dir) + + # 4. 生成总结报告 + print("\n生成总结报告...") + summary = generate_summary_report(baseline, save_dir) + + print("\n" + "=" * 80) + print("灵敏度分析完成!") + print("=" * 80) + + print(f"\n生成的文件:") + print(f" - sensitivity_recycle_rate.png") + print(f" - sensitivity_comfort_factor.png") + print(f" - sensitivity_medical_params.png") + print(f" - sensitivity_population.png") + print(f" - sensitivity_buffer_days.png") + print(f" - tornado_analysis.png") + print(f" - interaction_heatmap.png") + print(f" - three_param_surface.png") + print(f" - worst_case_analysis.png") + print(f" - worst_case_results.csv") + print(f" - sensitivity_summary_report.txt") + + return { + 'recycle': df_recycle, + 'comfort': df_comfort, + 'population': df_population, + 'buffer': df_buffer, + 'worst_case': worst_case_df, + } + + +if __name__ == "__main__": + results = run_full_sensitivity_analysis() diff --git a/plane/__pycache__/cost_analysis.cpython-314.pyc b/plane/__pycache__/cost_analysis.cpython-314.pyc new file mode 100644 index 0000000..50d86d5 Binary files /dev/null and b/plane/__pycache__/cost_analysis.cpython-314.pyc differ diff --git a/plane/__pycache__/cost_analysis.cpython-38.pyc b/plane/__pycache__/cost_analysis.cpython-38.pyc new file mode 100644 index 0000000..11c4d04 Binary files /dev/null and b/plane/__pycache__/cost_analysis.cpython-38.pyc differ diff --git a/plane/cost_analysis.py b/plane/cost_analysis.py index 44d53f2..7517d0c 100644 --- a/plane/cost_analysis.py +++ b/plane/cost_analysis.py @@ -556,6 +556,143 @@ def plot_comprehensive_summary(data, save_path='aviation_cost_summary.png'): return fig +# ============== 图表6: 燃油占比趋势（简洁版，中低饱和度）============== + +def plot_fuel_share_trend_simple(data, save_path='fuel_share_trend.png'): + """绘制燃油成本占总成本的比例趋势 - 简洁版，无标题，中低饱和度色调""" + + years = data['years'] + fuel_cents = data['fuel_cents'] + total_cents = data['total_cents'] + + # 计算燃油占比 + fuel_share = (fuel_cents / total_cents) * 100 + + # 创建图表 - 更小但保持2:1比例 + fig, ax1 = plt.subplots(figsize=(10, 5)) + + # 中低饱和度色调 + color1 = '#B87373' # 柔和的红色/玫瑰色 + color2 = '#6B8FAF' # 柔和的蓝色 + color3 = '#7AAE7A' # 柔和的绿色 + + # 绘制燃油占比（左轴） + ax1.fill_between(years, 0, fuel_share, alpha=0.25, color=color1) + line1, = ax1.plot(years, fuel_share, 'o-', color=color1, linewidth=2, + markersize=6, label='Fuel Share (%)') + ax1.set_xlabel('Year', fontsize=11) + ax1.set_ylabel('Fuel as % of Total Operating Cost', fontsize=11, color=color1) + ax1.tick_params(axis='y', labelcolor=color1) + ax1.set_ylim(0, 45) + + # 添加平均线 + avg_share = np.mean(fuel_share) + ax1.axhline(y=avg_share, color=color1, linestyle='--', linewidth=1.2, alpha=0.6) + ax1.text(years[-1]+0.5, avg_share+1, f'Avg: {avg_share:.1f}%', + fontsize=10, color=color1, fontweight='bold') + + # 右轴：绝对成本 + ax2 = ax1.twinx() + line2, = ax2.plot(years, total_cents, 's--', color=color2, linewidth=1.8, + markersize=5, alpha=0.75, label='Total CASM (cents)') + line3, = ax2.plot(years, fuel_cents, '^--', color=color3, linewidth=1.8, + markersize=5, alpha=0.75, label='Fuel CASM (cents)') + ax2.set_ylabel('Cost per ASM (cents)', fontsize=11, color=color2) + ax2.tick_params(axis='y', labelcolor=color2) + + # 标注关键年份 + # 2008年油价高峰 + idx_2008 = np.where(years == 2008)[0][0] + ax1.annotate(f'2008 Oil Crisis\n{fuel_share[idx_2008]:.1f}%', + xy=(2008, fuel_share[idx_2008]), + xytext=(2005, fuel_share[idx_2008]+8), + fontsize=9, ha='center', + arrowprops=dict(arrowstyle='->', color='#888888'), + bbox=dict(boxstyle='round,pad=0.3', facecolor='#FFF5E6', alpha=0.8)) + + # 2015年油价低谷 + idx_2015 = np.where(years == 2015)[0][0] + ax1.annotate(f'2015 Oil Crash\n{fuel_share[idx_2015]:.1f}%', + xy=(2015, fuel_share[idx_2015]), + xytext=(2017, fuel_share[idx_2015]-8), + fontsize=9, ha='center', + arrowprops=dict(arrowstyle='->', color='#888888'), + bbox=dict(boxstyle='round,pad=0.3', facecolor='#E6F3FA', alpha=0.8)) + + # 图例 + lines = [line1, line2, line3] + labels = [l.get_label() for l in lines] + ax1.legend(lines, labels, loc='upper left', fontsize=10) + + # 添加数据来源 + fig.text(0.99, 0.01, 'Data Source: MIT Airline Data Project / US DOT Form 41', + fontsize=8, ha='right', style='italic', alpha=0.6) + + plt.tight_layout() + plt.savefig(save_path, dpi=150, bbox_inches='tight') + print(f"图表已保存: {save_path}") + + return fig + + +# ============== 图表7: 燃油-总成本相关性散点图（单独版）============== + +def plot_fuel_total_correlation_scatter(data, save_path='fuel_total_correlation.png'): + """绘制燃油成本与总成本的相关性散点图 - 单独版，无标题，中低饱和度色调""" + + years = data['years'] + fuel = data['fuel_cents'] + total = data['total_cents'] + + # 计算相关系数 + correlation = np.corrcoef(fuel, total)[0, 1] + + # 线性回归 + coeffs = np.polyfit(fuel, total, 1) + poly = np.poly1d(coeffs) + fuel_range = np.linspace(fuel.min(), fuel.max(), 100) + + # 创建图表 - 与fuel_share_trend类似的尺寸比例 + fig, ax = plt.subplots(figsize=(10, 5)) + + # 中低饱和度色调 - 与fuel_share_trend类似但不完全相同 + scatter_cmap = 'YlGnBu' # 柔和的黄-绿-蓝渐变色 + line_color = '#9B7B6B' # 柔和的棕红色 + text_bg_color = '#D4E8D4' # 柔和的淡绿色 + + # 散点图 + scatter = ax.scatter(fuel, total, c=years, cmap=scatter_cmap, s=90, alpha=0.85, + edgecolors='white', linewidth=1.2) + + # 回归线 + ax.plot(fuel_range, poly(fuel_range), '--', color=line_color, linewidth=2, + label=f'Linear Fit (R²={correlation**2:.3f})') + + # 颜色条 + cbar = plt.colorbar(scatter, ax=ax) + cbar.set_label('Year', fontsize=10) + + ax.set_xlabel('Fuel Cost per ASM (cents)', fontsize=11) + ax.set_ylabel('Total Cost per ASM (cents)', fontsize=11) + ax.legend(fontsize=10, loc='lower right') + ax.grid(True, alpha=0.3) + + # 添加相关性解释 + ax.text(0.05, 0.95, f'Pearson r = {correlation:.3f}\nEnergy drives ~{correlation**2*100:.0f}% of cost variance', + transform=ax.transAxes, fontsize=10, verticalalignment='top', + bbox=dict(boxstyle='round', facecolor=text_bg_color, alpha=0.8)) + + # 添加数据来源 + fig.text(0.99, 0.01, 'Data Source: MIT Airline Data Project / US DOT Form 41', + fontsize=8, ha='right', style='italic', alpha=0.6) + + plt.tight_layout() + plt.savefig(save_path, dpi=150, bbox_inches='tight') + print(f"图表已保存: {save_path}") + + return fig + + # ============== 主程序 ============== if __name__ == "__main__": diff --git a/plane/fuel_share_trend.png b/plane/fuel_share_trend.png index 0479843..f9395d0 100644 Binary files a/plane/fuel_share_trend.png and b/plane/fuel_share_trend.png differ diff --git a/plane/fuel_total_correlation.png b/plane/fuel_total_correlation.png new file mode 100644 index 0000000..e554568 Binary files /dev/null and b/plane/fuel_total_correlation.png differ diff --git a/prob/p_B.pdf b/prob/p_B.pdf index e66cc47..4a71604 100644 Binary files a/prob/p_B.pdf and b/prob/p_B.pdf differ