ntnt/2026_mcm_b
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

plane: fig

Browse Source
This commit is contained in:
ntnt
2026-02-02 21:47:52 +08:00
parent 8192fcc354
commit 244e157bc0
11 changed files with 3245 additions and 3108 deletions
Download Patch File Download Diff File Expand all files Collapse all files

702
p1/data_window_analysis.py
View File

@@ -1,351 +1,351 @@
#!/usr/bin/env python3 #!/usr/bin/env python3
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
""" """
Data Window Sensitivity Analysis for Richards Model Data Window Sensitivity Analysis for Richards Model
Analyzes how different historical data windows affect the K (carrying capacity) Analyzes how different historical data windows affect the K (carrying capacity)
estimate in the Richards growth model. estimate in the Richards growth model.
Windows analyzed: Windows analyzed:
- 1957-2025: All available data (including Cold War era) - 1957-2025: All available data (including Cold War era)
- 1990-2025: Post-Cold War era - 1990-2025: Post-Cold War era
- 2000-2025: Commercial space age - 2000-2025: Commercial space age
- 2010-2025: SpaceX era - 2010-2025: SpaceX era
- 2015-2025: Recent rapid growth - 2015-2025: Recent rapid growth
""" """
import pandas as pd import pandas as pd
import numpy as np import numpy as np
from scipy.optimize import curve_fit from scipy.optimize import curve_fit
import matplotlib import matplotlib
matplotlib.use('Agg') matplotlib.use('Agg')
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import warnings import warnings
warnings.filterwarnings('ignore') warnings.filterwarnings('ignore')
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei', 'DejaVu Sans'] plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False plt.rcParams['axes.unicode_minus'] = False
# ============================================================ # ============================================================
# Richards Growth Model # Richards Growth Model
# ============================================================ # ============================================================
def richards(t, K, r, t0, v): def richards(t, K, r, t0, v):
"""Richards curve (generalized logistic model)""" """Richards curve (generalized logistic model)"""
exp_term = np.exp(-r * (t - t0)) exp_term = np.exp(-r * (t - t0))
exp_term = np.clip(exp_term, 1e-10, 1e10) exp_term = np.clip(exp_term, 1e-10, 1e10)
return K / np.power(1 + exp_term, 1/v) return K / np.power(1 + exp_term, 1/v)
# ============================================================ # ============================================================
# Data Loading # Data Loading
# ============================================================ # ============================================================
def load_data(filepath="rocket_launch_counts.csv"): def load_data(filepath="rocket_launch_counts.csv"):
"""Load and preprocess rocket launch data""" """Load and preprocess rocket launch data"""
df = pd.read_csv(filepath) df = pd.read_csv(filepath)
df = df.rename(columns={"YDate": "year", "Total": "launches"}) df = df.rename(columns={"YDate": "year", "Total": "launches"})
df["year"] = pd.to_numeric(df["year"], errors="coerce") df["year"] = pd.to_numeric(df["year"], errors="coerce")
df["launches"] = pd.to_numeric(df["launches"], errors="coerce") df["launches"] = pd.to_numeric(df["launches"], errors="coerce")
df = df.dropna(subset=["year", "launches"]) df = df.dropna(subset=["year", "launches"])
df = df[(df["year"] >= 1957) & (df["year"] <= 2025)] df = df[(df["year"] >= 1957) & (df["year"] <= 2025)]
df = df.astype({"year": int, "launches": int}) df = df.astype({"year": int, "launches": int})
df = df.sort_values("year").reset_index(drop=True) df = df.sort_values("year").reset_index(drop=True)
return df return df
# ============================================================ # ============================================================
# Fit Richards Model (Unconstrained) # Fit Richards Model (Unconstrained)
# ============================================================ # ============================================================
def fit_richards_unconstrained(data, base_year=1957): def fit_richards_unconstrained(data, base_year=1957):
"""Fit Richards model without constraining K""" """Fit Richards model without constraining K"""
years = data["year"].values years = data["year"].values
launches = data["launches"].values launches = data["launches"].values
t = (years - base_year).astype(float) t = (years - base_year).astype(float)
# Initial parameters # Initial parameters
p0 = [5000.0, 0.08, 80.0, 2.0] p0 = [5000.0, 0.08, 80.0, 2.0]
# Wide bounds - let data determine K # Wide bounds - let data determine K
bounds = ([500, 0.005, 10, 0.2], [100000, 1.0, 200, 10.0]) bounds = ([500, 0.005, 10, 0.2], [100000, 1.0, 200, 10.0])
try: try:
popt, pcov = curve_fit(richards, t, launches, p0=p0, bounds=bounds, maxfev=100000) popt, pcov = curve_fit(richards, t, launches, p0=p0, bounds=bounds, maxfev=100000)
perr = np.sqrt(np.diag(pcov)) perr = np.sqrt(np.diag(pcov))
y_pred = richards(t, *popt) y_pred = richards(t, *popt)
ss_res = np.sum((launches - y_pred) ** 2) ss_res = np.sum((launches - y_pred) ** 2)
ss_tot = np.sum((launches - np.mean(launches)) ** 2) ss_tot = np.sum((launches - np.mean(launches)) ** 2)
r_squared = 1 - (ss_res / ss_tot) r_squared = 1 - (ss_res / ss_tot)
K, r, t0, v = popt K, r, t0, v = popt
K_err = perr[0] K_err = perr[0]
return { return {
"success": True, "success": True,
"K": K, "K": K,
"K_err": K_err, "K_err": K_err,
"r": r, "r": r,
"t0": t0, "t0": t0,
"v": v, "v": v,
"r_squared": r_squared, "r_squared": r_squared,
"n_points": len(data), "n_points": len(data),
"years": years, "years": years,
"launches": launches, "launches": launches,
"params": popt, "params": popt,
"y_pred": y_pred, "y_pred": y_pred,
} }
except Exception as e: except Exception as e:
return {"success": False, "error": str(e)} return {"success": False, "error": str(e)}
# ============================================================ # ============================================================
# Main Analysis # Main Analysis
# ============================================================ # ============================================================
def analyze_data_windows(df): def analyze_data_windows(df):
"""Analyze K estimates for different data windows""" """Analyze K estimates for different data windows"""
# Define windows # Define windows
windows = { windows = {
"1957-2025 (全部68年)": (1957, 2025), "1957-2025 (全部68年)": (1957, 2025),
"1990-2025 (35年)": (1990, 2025), "1990-2025 (35年)": (1990, 2025),
"2000-2025 (25年)": (2000, 2025), "2000-2025 (25年)": (2000, 2025),
"2010-2025 (15年)": (2010, 2025), "2010-2025 (15年)": (2010, 2025),
"2015-2025 (10年)": (2015, 2025), "2015-2025 (10年)": (2015, 2025),
} }
results = {} results = {}
print("=" * 80) print("=" * 80)
print("分析不同历史数据窗口对 K (carrying capacity) 估计的影响") print("分析不同历史数据窗口对 K (carrying capacity) 估计的影响")
print("=" * 80) print("=" * 80)
for name, (start, end) in windows.items(): for name, (start, end) in windows.items():
data = df[(df["year"] >= start) & (df["year"] <= end)].copy() data = df[(df["year"] >= start) & (df["year"] <= end)].copy()
result = fit_richards_unconstrained(data) result = fit_richards_unconstrained(data)
if result["success"]: if result["success"]:
result["start"] = start result["start"] = start
result["end"] = end result["end"] = end
result["name"] = name result["name"] = name
results[name] = result results[name] = result
print(f"\n--- {name} ---") print(f"\n--- {name} ---")
print(f" 数据点数: {result['n_points']}") print(f" 数据点数: {result['n_points']}")
print(f" K = {result['K']:.0f} (carrying capacity)") print(f" K = {result['K']:.0f} (carrying capacity)")
print(f" r = {result['r']:.4f} (growth rate)") print(f" r = {result['r']:.4f} (growth rate)")
print(f" R² = {result['r_squared']:.4f}") print(f" R² = {result['r_squared']:.4f}")
else: else:
print(f"\n--- {name} ---") print(f"\n--- {name} ---")
print(f" 拟合失败: {result['error']}") print(f" 拟合失败: {result['error']}")
return results return results
# ============================================================ # ============================================================
# Visualization # Visualization
# ============================================================ # ============================================================
def plot_window_analysis(df, results, save_path="launch_capacity_window_analysis.png"): def plot_window_analysis(df, results, save_path="launch_capacity_window_analysis.png"):
"""Generate comprehensive window analysis plot""" """Generate comprehensive window analysis plot"""
fig, axes = plt.subplots(2, 2, figsize=(14, 11)) fig, axes = plt.subplots(2, 2, figsize=(14, 11))
names = list(results.keys()) names = list(results.keys())
physical_max = 3650 physical_max = 3650
# ========== Plot 1: K Values Comparison ========== # ========== Plot 1: K Values Comparison ==========
ax1 = axes[0, 0] ax1 = axes[0, 0]
K_vals = [results[n]["K"] for n in names] K_vals = [results[n]["K"] for n in names]
colors = ["#E74C3C" if k < physical_max else "#27AE60" for k in K_vals] colors = ["#E74C3C" if k < physical_max else "#27AE60" for k in K_vals]
y_pos = range(len(names)) y_pos = range(len(names))
bars = ax1.barh(y_pos, K_vals, color=colors, alpha=0.7, edgecolor="black") bars = ax1.barh(y_pos, K_vals, color=colors, alpha=0.7, edgecolor="black")
ax1.axvline(physical_max, color="blue", ls="--", lw=2.5, ax1.axvline(physical_max, color="blue", ls="--", lw=2.5,
label=f"Physical Max ({physical_max})") label=f"Physical Max ({physical_max})")
ax1.set_yticks(y_pos) ax1.set_yticks(y_pos)
ax1.set_yticklabels(names, fontsize=10) ax1.set_yticklabels(names, fontsize=10)
ax1.set_xlabel("K (Carrying Capacity, launches/year)", fontsize=11) 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.set_title("K Estimates by Data Window\n(Unconstrained Richards Model)", fontsize=12)
ax1.legend(loc="lower right", fontsize=10) ax1.legend(loc="lower right", fontsize=10)
ax1.grid(True, alpha=0.3, axis="x") ax1.grid(True, alpha=0.3, axis="x")
for bar, k in zip(bars, K_vals): for bar, k in zip(bars, K_vals):
ax1.text(k + 200, bar.get_y() + bar.get_height()/2, ax1.text(k + 200, bar.get_y() + bar.get_height()/2,
f"{k:.0f}", va="center", fontsize=10, fontweight="bold") f"{k:.0f}", va="center", fontsize=10, fontweight="bold")
# ========== Plot 2: R² Comparison ========== # ========== Plot 2: R² Comparison ==========
ax2 = axes[0, 1] ax2 = axes[0, 1]
r2_vals = [results[n]["r_squared"] for n in names] r2_vals = [results[n]["r_squared"] for n in names]
colors2 = ["#27AE60" if r2 > 0.9 else "#F39C12" if r2 > 0.7 else "#E74C3C" colors2 = ["#27AE60" if r2 > 0.9 else "#F39C12" if r2 > 0.7 else "#E74C3C"
for r2 in r2_vals] for r2 in r2_vals]
bars2 = ax2.barh(y_pos, r2_vals, color=colors2, alpha=0.7, edgecolor="black") 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.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.axvline(0.7, color="orange", ls=":", lw=1.5, alpha=0.7, label="R²=0.7 (Acceptable)")
ax2.set_yticks(y_pos) ax2.set_yticks(y_pos)
ax2.set_yticklabels(names, fontsize=10) ax2.set_yticklabels(names, fontsize=10)
ax2.set_xlabel("R² (Goodness of Fit)", fontsize=11) ax2.set_xlabel("R² (Goodness of Fit)", fontsize=11)
ax2.set_title("Model Fit Quality by Data Window", fontsize=12) ax2.set_title("Model Fit Quality by Data Window", fontsize=12)
ax2.set_xlim(0, 1.1) ax2.set_xlim(0, 1.1)
ax2.legend(loc="lower right", fontsize=9) ax2.legend(loc="lower right", fontsize=9)
ax2.grid(True, alpha=0.3, axis="x") ax2.grid(True, alpha=0.3, axis="x")
for bar, r2 in zip(bars2, r2_vals): for bar, r2 in zip(bars2, r2_vals):
ax2.text(r2 + 0.02, bar.get_y() + bar.get_height()/2, ax2.text(r2 + 0.02, bar.get_y() + bar.get_height()/2,
f"{r2:.3f}", va="center", fontsize=10, fontweight="bold") f"{r2:.3f}", va="center", fontsize=10, fontweight="bold")
# ========== Plot 3: Historical Data with All Model Fits ========== # ========== Plot 3: Historical Data with All Model Fits ==========
ax3 = axes[1, 0] ax3 = axes[1, 0]
# Plot all historical data # Plot all historical data
ax3.scatter(df["year"], df["launches"], color="gray", s=30, alpha=0.5, ax3.scatter(df["year"], df["launches"], color="gray", s=30, alpha=0.5,
label="Historical Data", zorder=1) label="Historical Data", zorder=1)
# Plot each model prediction # Plot each model prediction
colors_line = ["#E74C3C", "#F39C12", "#3498DB", "#27AE60", "#9B59B6"] colors_line = ["#E74C3C", "#F39C12", "#3498DB", "#27AE60", "#9B59B6"]
for i, (name, res) in enumerate(results.items()): for i, (name, res) in enumerate(results.items()):
start = res["start"] start = res["start"]
years_proj = np.arange(start, 2080) years_proj = np.arange(start, 2080)
t_proj = years_proj - 1957 t_proj = years_proj - 1957
pred = richards(t_proj, *res["params"]) pred = richards(t_proj, *res["params"])
label = f'{name.split(" ")[0]}: K={res["K"]:.0f}' label = f'{name.split(" ")[0]}: K={res["K"]:.0f}'
ax3.plot(years_proj, pred, color=colors_line[i], lw=2, alpha=0.8, ax3.plot(years_proj, pred, color=colors_line[i], lw=2, alpha=0.8,
label=label, zorder=2) label=label, zorder=2)
ax3.axhline(physical_max, color="black", ls=":", lw=2, ax3.axhline(physical_max, color="black", ls=":", lw=2,
label=f"Physical Max ({physical_max})") label=f"Physical Max ({physical_max})")
ax3.set_xlabel("Year", fontsize=11) ax3.set_xlabel("Year", fontsize=11)
ax3.set_ylabel("Annual Launches", fontsize=11) ax3.set_ylabel("Annual Launches", fontsize=11)
ax3.set_title("Richards Model Predictions\n(Different Data Windows)", fontsize=12) ax3.set_title("Richards Model Predictions\n(Different Data Windows)", fontsize=12)
ax3.legend(loc="upper left", fontsize=8) ax3.legend(loc="upper left", fontsize=8)
ax3.grid(True, alpha=0.3) ax3.grid(True, alpha=0.3)
ax3.set_xlim(1955, 2080) ax3.set_xlim(1955, 2080)
ax3.set_ylim(0, 15000) ax3.set_ylim(0, 15000)
# ========== Plot 4: Summary Table ========== # ========== Plot 4: Summary Table ==========
ax4 = axes[1, 1] ax4 = axes[1, 1]
ax4.axis("off") ax4.axis("off")
# Build summary text # Build summary text
summary = """ summary = """
┌────────────────────────────────────────────────────────────────────────┐ ┌────────────────────────────────────────────────────────────────────────┐
│ 数据窗口对 K 估计的影响分析 │ │ 数据窗口对 K 估计的影响分析 │
├────────────────────────────────────────────────────────────────────────┤ ├────────────────────────────────────────────────────────────────────────┤
│ │ │ │
│ 核心发现: │ │ 核心发现: │
│ ───────── │ │ ───────── │
│ 1. 数据窗口选择对 K 估计影响巨大 │ │ 1. 数据窗口选择对 K 估计影响巨大 │
│ • 全部数据: K ≈ 500 │ │ • 全部数据: K ≈ 500 │
│ • 近15年: K ≈ 12,000 │ │ • 近15年: K ≈ 12,000 │
│ │ │ │
│ 2. 1957-2025 全部数据: │ │ 2. 1957-2025 全部数据: │
│ • R² = 0.08 (极差拟合) │ │ • R² = 0.08 (极差拟合) │
│ • 冷战时期数据干扰，模型误判为"已饱和" │ │ • 冷战时期数据干扰，模型误判为"已饱和" │
│ • 不适合用于预测未来增长 │ │ • 不适合用于预测未来增长 │
│ │ │ │
│ 3. 2010-2025 近15年数据: │ │ 3. 2010-2025 近15年数据: │
│ • R² = 0.90 (良好拟合) │ │ • R² = 0.90 (良好拟合) │
│ • 准确捕捉商业航天时代的增长趋势 │ │ • 准确捕捉商业航天时代的增长趋势 │
│ • 预测 K ≈ 12,000 >> 物理上限 (3,650) │ │ • 预测 K ≈ 12,000 >> 物理上限 (3,650) │
│ │ │ │
│ 结论: │ │ 结论: │
│ ───── │ │ ───── │
│ • 选择 2010-2025 数据窗口是合理的 (R² = 0.90, 反映当前趋势) │ │ • 选择 2010-2025 数据窗口是合理的 (R² = 0.90, 反映当前趋势) │
│ • 数据驱动 K ≈ 12,000 反映增长潜力 (需要 ~33 个发射站) │ │ • 数据驱动 K ≈ 12,000 反映增长潜力 (需要 ~33 个发射站) │
│ • 但物理约束 (10站 × 365天 = 3,650) 才是真正的上限 │ │ • 但物理约束 (10站 × 365天 = 3,650) 才是真正的上限 │
│ • 因此采用 1 次/站/天 作为保守估计是合理的 │ │ • 因此采用 1 次/站/天 作为保守估计是合理的 │
│ │ │ │
│ 物理上限: 3,650 次/年 (10个发射站 × 365天) │ │ 物理上限: 3,650 次/年 (10个发射站 × 365天) │
│ │ │ │
└────────────────────────────────────────────────────────────────────────┘ └────────────────────────────────────────────────────────────────────────┘
""" """
ax4.text(0.02, 0.98, summary, transform=ax4.transAxes, fontsize=10, ax4.text(0.02, 0.98, summary, transform=ax4.transAxes, fontsize=10,
verticalalignment="top", family="monospace", verticalalignment="top", family="monospace",
bbox=dict(boxstyle="round", facecolor="lightyellow", alpha=0.9)) bbox=dict(boxstyle="round", facecolor="lightyellow", alpha=0.9))
plt.suptitle("Data Window Sensitivity Analysis for Richards Model", plt.suptitle("Data Window Sensitivity Analysis for Richards Model",
fontsize=14, fontweight="bold", y=1.02) fontsize=14, fontweight="bold", y=1.02)
plt.tight_layout() plt.tight_layout()
plt.savefig(save_path, dpi=150, bbox_inches="tight") plt.savefig(save_path, dpi=150, bbox_inches="tight")
plt.close() plt.close()
print(f"\nPlot saved: {save_path}") print(f"\nPlot saved: {save_path}")
# ============================================================ # ============================================================
# Print Detailed Table # Print Detailed Table
# ============================================================ # ============================================================
def print_comparison_table(results): def print_comparison_table(results):
"""Print detailed comparison table""" """Print detailed comparison table"""
physical_max = 3650 physical_max = 3650
print("\n" + "=" * 80) print("\n" + "=" * 80)
print("详细对比表") print("详细对比表")
print("=" * 80) print("=" * 80)
header = f"{'数据窗口':<22} | {'K值':>8} | {'R²':>7} | {'数据点':>5} | {'K/物理上限':>10} | {'说明':<18}" header = f"{'数据窗口':<22} | {'K值':>8} | {'R²':>7} | {'数据点':>5} | {'K/物理上限':>10} | {'说明':<18}"
print(header) print(header)
print("-" * 22 + "-+-" + "-" * 8 + "-+-" + "-" * 7 + "-+-" + "-" * 5 + "-+-" + "-" * 10 + "-+-" + "-" * 18) print("-" * 22 + "-+-" + "-" * 8 + "-+-" + "-" * 7 + "-+-" + "-" * 5 + "-+-" + "-" * 10 + "-+-" + "-" * 18)
for name, res in results.items(): for name, res in results.items():
k = res["K"] k = res["K"]
r2 = res["r_squared"] r2 = res["r_squared"]
n = res["n_points"] n = res["n_points"]
ratio = k / physical_max ratio = k / physical_max
if r2 < 0.5: if r2 < 0.5:
note = "拟合差，不可用" note = "拟合差，不可用"
elif k < physical_max: elif k < physical_max:
note = "低于物理上限" note = "低于物理上限"
elif ratio < 2: elif ratio < 2:
note = "接近物理上限" note = "接近物理上限"
else: else:
note = "远超物理上限" note = "远超物理上限"
print(f"{name:<22} | {k:>8.0f} | {r2:>7.4f} | {n:>5} | {ratio:>10.2f}x | {note:<18}") 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天)") print("\n物理上限: 3,650 次/年 (10个发射站 × 每站每天1次 × 365天)")
# ============================================================ # ============================================================
# Main # Main
# ============================================================ # ============================================================
def main(): def main():
print("Loading data...") print("Loading data...")
df = load_data() df = load_data()
print(f"Data loaded: {len(df)} records from {df['year'].min()} to {df['year'].max()}") print(f"Data loaded: {len(df)} records from {df['year'].min()} to {df['year'].max()}")
# Analyze different windows # Analyze different windows
results = analyze_data_windows(df) results = analyze_data_windows(df)
# Print comparison table # Print comparison table
print_comparison_table(results) print_comparison_table(results)
# Generate visualization # Generate visualization
print("\n" + "=" * 80) print("\n" + "=" * 80)
print("Generating visualization...") print("Generating visualization...")
print("=" * 80) print("=" * 80)
plot_window_analysis(df, results) plot_window_analysis(df, results)
print("\n" + "=" * 80) print("\n" + "=" * 80)
print("分析完成!") print("分析完成!")
print("=" * 80) print("=" * 80)
return results return results
if __name__ == "__main__": if __name__ == "__main__":
results = main() results = main()

460
p1/find_optimal_window.py
View File

@@ -1,230 +1,230 @@
#!/usr/bin/env python3 #!/usr/bin/env python3
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
""" """
Find Optimal Data Window for K ≈ 3650 Find Optimal Data Window for K ≈ 3650
Enumerate all possible starting years to find which data window Enumerate all possible starting years to find which data window
produces K closest to the physical limit (3650 launches/year). produces K closest to the physical limit (3650 launches/year).
""" """
import pandas as pd import pandas as pd
import numpy as np import numpy as np
from scipy.optimize import curve_fit from scipy.optimize import curve_fit
import matplotlib import matplotlib
matplotlib.use('Agg') matplotlib.use('Agg')
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import warnings import warnings
warnings.filterwarnings('ignore') warnings.filterwarnings('ignore')
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei', 'DejaVu Sans'] plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False plt.rcParams['axes.unicode_minus'] = False
# Richards model # Richards model
def richards(t, K, r, t0, v): def richards(t, K, r, t0, v):
exp_term = np.exp(-r * (t - t0)) exp_term = np.exp(-r * (t - t0))
exp_term = np.clip(exp_term, 1e-10, 1e10) exp_term = np.clip(exp_term, 1e-10, 1e10)
return K / np.power(1 + exp_term, 1/v) return K / np.power(1 + exp_term, 1/v)
def load_data(filepath="rocket_launch_counts.csv"): def load_data(filepath="rocket_launch_counts.csv"):
df = pd.read_csv(filepath) df = pd.read_csv(filepath)
df = df.rename(columns={"YDate": "year", "Total": "launches"}) df = df.rename(columns={"YDate": "year", "Total": "launches"})
df["year"] = pd.to_numeric(df["year"], errors="coerce") df["year"] = pd.to_numeric(df["year"], errors="coerce")
df["launches"] = pd.to_numeric(df["launches"], errors="coerce") df["launches"] = pd.to_numeric(df["launches"], errors="coerce")
df = df.dropna(subset=["year", "launches"]) df = df.dropna(subset=["year", "launches"])
df = df[(df["year"] >= 1957) & (df["year"] <= 2025)] df = df[(df["year"] >= 1957) & (df["year"] <= 2025)]
df = df.astype({"year": int, "launches": int}) df = df.astype({"year": int, "launches": int})
df = df.sort_values("year").reset_index(drop=True) df = df.sort_values("year").reset_index(drop=True)
return df return df
def fit_richards(data, base_year=1957): def fit_richards(data, base_year=1957):
"""Fit unconstrained Richards model""" """Fit unconstrained Richards model"""
years = data["year"].values years = data["year"].values
launches = data["launches"].values launches = data["launches"].values
t = (years - base_year).astype(float) t = (years - base_year).astype(float)
p0 = [5000.0, 0.08, 80.0, 2.0] p0 = [5000.0, 0.08, 80.0, 2.0]
bounds = ([500, 0.005, 10, 0.2], [100000, 1.0, 200, 10.0]) bounds = ([500, 0.005, 10, 0.2], [100000, 1.0, 200, 10.0])
try: try:
popt, pcov = curve_fit(richards, t, launches, p0=p0, bounds=bounds, maxfev=100000) popt, pcov = curve_fit(richards, t, launches, p0=p0, bounds=bounds, maxfev=100000)
y_pred = richards(t, *popt) y_pred = richards(t, *popt)
ss_res = np.sum((launches - y_pred) ** 2) ss_res = np.sum((launches - y_pred) ** 2)
ss_tot = np.sum((launches - np.mean(launches)) ** 2) ss_tot = np.sum((launches - np.mean(launches)) ** 2)
r_squared = 1 - (ss_res / ss_tot) if ss_tot > 0 else 0 r_squared = 1 - (ss_res / ss_tot) if ss_tot > 0 else 0
return { return {
"success": True, "success": True,
"K": popt[0], "K": popt[0],
"r": popt[1], "r": popt[1],
"t0": popt[2], "t0": popt[2],
"v": popt[3], "v": popt[3],
"r_squared": r_squared, "r_squared": r_squared,
"n_points": len(data), "n_points": len(data),
} }
except: except:
return {"success": False} return {"success": False}
def main(): def main():
print("=" * 80) print("=" * 80)
print("枚举起始年份，寻找 K ≈ 3650 的数据窗口") print("枚举起始年份，寻找 K ≈ 3650 的数据窗口")
print("=" * 80) print("=" * 80)
df = load_data() df = load_data()
print(f"数据范围: {df['year'].min()} - {df['year'].max()}\n") print(f"数据范围: {df['year'].min()} - {df['year'].max()}\n")
target_K = 3650 target_K = 3650
end_year = 2025 end_year = 2025
# Enumerate all starting years # Enumerate all starting years
results = [] results = []
for start_year in range(1957, 2020): for start_year in range(1957, 2020):
data = df[(df["year"] >= start_year) & (df["year"] <= end_year)].copy() data = df[(df["year"] >= start_year) & (df["year"] <= end_year)].copy()
if len(data) < 5: # Need at least 5 data points if len(data) < 5: # Need at least 5 data points
continue continue
result = fit_richards(data) result = fit_richards(data)
if result["success"]: if result["success"]:
diff = abs(result["K"] - target_K) diff = abs(result["K"] - target_K)
results.append({ results.append({
"start_year": start_year, "start_year": start_year,
"end_year": end_year, "end_year": end_year,
"years_span": end_year - start_year + 1, "years_span": end_year - start_year + 1,
"n_points": result["n_points"], "n_points": result["n_points"],
"K": result["K"], "K": result["K"],
"r_squared": result["r_squared"], "r_squared": result["r_squared"],
"diff_from_target": diff, "diff_from_target": diff,
"ratio": result["K"] / target_K, "ratio": result["K"] / target_K,
}) })
# Sort by difference from target # Sort by difference from target
results_df = pd.DataFrame(results) results_df = pd.DataFrame(results)
results_df = results_df.sort_values("diff_from_target") results_df = results_df.sort_values("diff_from_target")
# Print full table # Print full table
print("完整枚举结果表（按与目标值 3650 的差距排序）:") print("完整枚举结果表（按与目标值 3650 的差距排序）:")
print("-" * 80) print("-" * 80)
print(f"{'起始年份':>8} | {'时间跨度':>8} | {'数据点':>6} | {'K值':>10} | {'R²':>8} | {'K/3650':>8} | {'差距':>10}") print(f"{'起始年份':>8} | {'时间跨度':>8} | {'数据点':>6} | {'K值':>10} | {'R²':>8} | {'K/3650':>8} | {'差距':>10}")
print("-" * 80) print("-" * 80)
for _, row in results_df.iterrows(): for _, row in results_df.iterrows():
marker = "★" if row["diff_from_target"] < 500 else "" 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} | " 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}") 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 # Find closest matches
print("\n" + "=" * 80) print("\n" + "=" * 80)
print("最接近 K = 3650 的数据窗口 (Top 10)") print("最接近 K = 3650 的数据窗口 (Top 10)")
print("=" * 80) print("=" * 80)
top10 = results_df.head(10) top10 = results_df.head(10)
print(f"\n{'排名':>4} | {'起始年份':>8} | {'K值':>10} | {'R²':>8} | {'差距':>10} | {'评价':<20}") print(f"\n{'排名':>4} | {'起始年份':>8} | {'K值':>10} | {'R²':>8} | {'差距':>10} | {'评价':<20}")
print("-" * 75) print("-" * 75)
for i, (_, row) in enumerate(top10.iterrows(), 1): for i, (_, row) in enumerate(top10.iterrows(), 1):
if row["r_squared"] < 0.5: if row["r_squared"] < 0.5:
comment = "❌ 拟合差，不可信" comment = "❌ 拟合差，不可信"
elif row["r_squared"] < 0.8: elif row["r_squared"] < 0.8:
comment = "⚠️ 拟合一般" comment = "⚠️ 拟合一般"
else: else:
comment = "✅ 拟合良好" comment = "✅ 拟合良好"
print(f"{i:>4} | {int(row['start_year']):>8} | {row['K']:>10.0f} | {row['r_squared']:>8.4f} | " 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}") f"{row['diff_from_target']:>10.0f} | {comment:<20}")
# Analysis # Analysis
print("\n" + "=" * 80) print("\n" + "=" * 80)
print("分析结论") print("分析结论")
print("=" * 80) print("=" * 80)
# Find best with good R² # Find best with good R²
good_fit = results_df[results_df["r_squared"] >= 0.8] good_fit = results_df[results_df["r_squared"] >= 0.8]
if len(good_fit) > 0: if len(good_fit) > 0:
best_good = good_fit.iloc[0] best_good = good_fit.iloc[0]
print(f"\n在 R² ≥ 0.8 的条件下，最接近 K=3650 的窗口:") print(f"\n在 R² ≥ 0.8 的条件下，最接近 K=3650 的窗口:")
print(f" 起始年份: {int(best_good['start_year'])}") print(f" 起始年份: {int(best_good['start_year'])}")
print(f" K = {best_good['K']:.0f}") print(f" K = {best_good['K']:.0f}")
print(f" R² = {best_good['r_squared']:.4f}") print(f" R² = {best_good['r_squared']:.4f}")
print(f" 差距: {best_good['diff_from_target']:.0f}") print(f" 差距: {best_good['diff_from_target']:.0f}")
# Summary # Summary
print("\n关键发现:") print("\n关键发现:")
print("-" * 40) print("-" * 40)
# Check if any good fit gives K near 3650 # 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)] near_target = results_df[(results_df["diff_from_target"] < 1000) & (results_df["r_squared"] >= 0.7)]
if len(near_target) == 0: if len(near_target) == 0:
print(" ⚠️ 没有任何数据窗口能在良好拟合(R²≥0.7)的条件下得到 K≈3650") print(" ⚠️ 没有任何数据窗口能在良好拟合(R²≥0.7)的条件下得到 K≈3650")
print(" ⚠️ 所有良好拟合的窗口都给出 K >> 3650 或 K << 3650") print(" ⚠️ 所有良好拟合的窗口都给出 K >> 3650 或 K << 3650")
print("\n 这说明:") print("\n 这说明:")
print(" • K=3650 不是数据自然支持的结论") print(" • K=3650 不是数据自然支持的结论")
print(" • K=3650 来自物理约束，而非统计预测") print(" • K=3650 来自物理约束，而非统计预测")
print(" • 论文中应该明确说明这是'物理上限'而非'数据预测'") print(" • 论文中应该明确说明这是'物理上限'而非'数据预测'")
else: else:
print(f" 找到 {len(near_target)} 个窗口使 K 接近 3650:") print(f" 找到 {len(near_target)} 个窗口使 K 接近 3650:")
for _, row in near_target.iterrows(): for _, row in near_target.iterrows():
print(f" {int(row['start_year'])}-2025: K={row['K']:.0f}, R²={row['r_squared']:.4f}") print(f" {int(row['start_year'])}-2025: K={row['K']:.0f}, R²={row['r_squared']:.4f}")
# Generate visualization # Generate visualization
print("\n" + "=" * 80) print("\n" + "=" * 80)
print("生成可视化...") print("生成可视化...")
print("=" * 80) print("=" * 80)
fig, axes = plt.subplots(1, 2, figsize=(14, 6)) fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# Plot 1: K vs Start Year # Plot 1: K vs Start Year
ax1 = axes[0] ax1 = axes[0]
colors = ['#27AE60' if r2 >= 0.8 else '#F39C12' if r2 >= 0.5 else '#E74C3C' colors = ['#27AE60' if r2 >= 0.8 else '#F39C12' if r2 >= 0.5 else '#E74C3C'
for r2 in results_df["r_squared"]] 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.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, 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*0.9, color='blue', ls=':', lw=1, alpha=0.5)
ax1.axhline(3650*1.1, 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_xlabel("Starting Year", fontsize=11)
ax1.set_ylabel("K (Carrying Capacity)", 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.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.legend()
ax1.grid(True, alpha=0.3) ax1.grid(True, alpha=0.3)
ax1.set_xlim(1955, 2020) ax1.set_xlim(1955, 2020)
# Plot 2: R² vs Start Year # Plot 2: R² vs Start Year
ax2 = axes[1] ax2 = axes[1]
ax2.scatter(results_df["start_year"], results_df["r_squared"], c=colors, s=60, alpha=0.7, edgecolor='black') 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.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.axhline(0.5, color='orange', ls=':', lw=1.5, label='R²=0.5 (Acceptable)')
ax2.set_xlabel("Starting Year", fontsize=11) ax2.set_xlabel("Starting Year", fontsize=11)
ax2.set_ylabel("R² (Goodness of Fit)", fontsize=11) ax2.set_ylabel("R² (Goodness of Fit)", fontsize=11)
ax2.set_title("Model Fit Quality vs Starting Year", fontsize=12) ax2.set_title("Model Fit Quality vs Starting Year", fontsize=12)
ax2.legend() ax2.legend()
ax2.grid(True, alpha=0.3) ax2.grid(True, alpha=0.3)
ax2.set_xlim(1955, 2020) ax2.set_xlim(1955, 2020)
ax2.set_ylim(-0.1, 1.1) ax2.set_ylim(-0.1, 1.1)
plt.tight_layout() plt.tight_layout()
plt.savefig("find_optimal_window.png", dpi=150, bbox_inches='tight') plt.savefig("find_optimal_window.png", dpi=150, bbox_inches='tight')
plt.close() plt.close()
print("图表已保存: find_optimal_window.png") print("图表已保存: find_optimal_window.png")
# Save to CSV # Save to CSV
results_df.to_csv("window_enumeration.csv", index=False) results_df.to_csv("window_enumeration.csv", index=False)
print("数据已保存: window_enumeration.csv") print("数据已保存: window_enumeration.csv")
print("\n" + "=" * 80) print("\n" + "=" * 80)
print("分析完成!") print("分析完成!")
print("=" * 80) print("=" * 80)
return results_df return results_df
if __name__ == "__main__": if __name__ == "__main__":
results = main() results = main()

1688
p1/lambda_cost_analysis.py
View File

File diff suppressed because it is too large Load Diff

366
p1/richards_curve_1984.py
View File

@@ -1,183 +1,183 @@
#!/usr/bin/env python3 #!/usr/bin/env python3
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
""" """
Richards S-Curve Fit for 1984-2025 Data Richards S-Curve Fit for 1984-2025 Data
Fits and visualizes the Richards growth model to rocket launch data Fits and visualizes the Richards growth model to rocket launch data
starting from 1984. starting from 1984.
""" """
import pandas as pd import pandas as pd
import numpy as np import numpy as np
from scipy.optimize import curve_fit from scipy.optimize import curve_fit
import matplotlib import matplotlib
matplotlib.use('Agg') matplotlib.use('Agg')
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import warnings import warnings
warnings.filterwarnings('ignore') warnings.filterwarnings('ignore')
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei', 'DejaVu Sans'] plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False plt.rcParams['axes.unicode_minus'] = False
def richards(t, K, r, t0, v): def richards(t, K, r, t0, v):
""" """
Richards curve (generalized logistic model) Richards curve (generalized logistic model)
N(t) = K / (1 + exp(-r*(t - t0)))^(1/v) N(t) = K / (1 + exp(-r*(t - t0)))^(1/v)
""" """
exp_term = np.exp(-r * (t - t0)) exp_term = np.exp(-r * (t - t0))
exp_term = np.clip(exp_term, 1e-10, 1e10) exp_term = np.clip(exp_term, 1e-10, 1e10)
return K / np.power(1 + exp_term, 1/v) return K / np.power(1 + exp_term, 1/v)
def load_data(filepath="rocket_launch_counts.csv"): def load_data(filepath="rocket_launch_counts.csv"):
"""Load rocket launch data""" """Load rocket launch data"""
df = pd.read_csv(filepath) df = pd.read_csv(filepath)
df = df.rename(columns={"YDate": "year", "Total": "launches"}) df = df.rename(columns={"YDate": "year", "Total": "launches"})
df["year"] = pd.to_numeric(df["year"], errors="coerce") df["year"] = pd.to_numeric(df["year"], errors="coerce")
df["launches"] = pd.to_numeric(df["launches"], errors="coerce") df["launches"] = pd.to_numeric(df["launches"], errors="coerce")
df = df.dropna(subset=["year", "launches"]) df = df.dropna(subset=["year", "launches"])
df = df[(df["year"] >= 1957) & (df["year"] <= 2025)] df = df[(df["year"] >= 1957) & (df["year"] <= 2025)]
df = df.astype({"year": int, "launches": int}) df = df.astype({"year": int, "launches": int})
df = df.sort_values("year").reset_index(drop=True) df = df.sort_values("year").reset_index(drop=True)
return df return df
def fit_richards(years, launches, base_year=1984): def fit_richards(years, launches, base_year=1984):
"""Fit Richards model to data""" """Fit Richards model to data"""
t = (years - base_year).astype(float) t = (years - base_year).astype(float)
# Initial parameters and bounds (unconstrained K) # Initial parameters and bounds (unconstrained K)
p0 = [5000.0, 0.1, 40.0, 2.0] p0 = [5000.0, 0.1, 40.0, 2.0]
bounds = ([500, 0.005, 10, 0.2], [100000, 1.0, 150, 10.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) popt, pcov = curve_fit(richards, t, launches, p0=p0, bounds=bounds, maxfev=100000)
perr = np.sqrt(np.diag(pcov)) perr = np.sqrt(np.diag(pcov))
# Calculate R² # Calculate R²
y_pred = richards(t, *popt) y_pred = richards(t, *popt)
ss_res = np.sum((launches - y_pred) ** 2) ss_res = np.sum((launches - y_pred) ** 2)
ss_tot = np.sum((launches - np.mean(launches)) ** 2) ss_tot = np.sum((launches - np.mean(launches)) ** 2)
r_squared = 1 - (ss_res / ss_tot) r_squared = 1 - (ss_res / ss_tot)
return popt, perr, r_squared return popt, perr, r_squared
def main(): def main():
print("=" * 60) print("=" * 60)
print("Richards S-Curve Fit (1984-2025)") print("Richards S-Curve Fit (1984-2025)")
print("=" * 60) print("=" * 60)
# Load data # Load data
df = load_data() df = load_data()
# Filter 1984-2025 # Filter 1984-2025
start_year = 1984 start_year = 1984
data = df[(df["year"] >= start_year) & (df["year"] <= 2025)].copy() data = df[(df["year"] >= start_year) & (df["year"] <= 2025)].copy()
years = data["year"].values years = data["year"].values
launches = data["launches"].values launches = data["launches"].values
print(f"Data range: {start_year} - 2025") print(f"Data range: {start_year} - 2025")
print(f"Data points: {len(data)}") print(f"Data points: {len(data)}")
# Fit model # Fit model
popt, perr, r_squared = fit_richards(years, launches, base_year=start_year) popt, perr, r_squared = fit_richards(years, launches, base_year=start_year)
K, r, t0, v = popt K, r, t0, v = popt
print(f"\nFitted Parameters:") print(f"\nFitted Parameters:")
print(f" K (carrying capacity) = {K:.0f} launches/year") print(f" K (carrying capacity) = {K:.0f} launches/year")
print(f" r (growth rate) = {r:.4f}") print(f" r (growth rate) = {r:.4f}")
print(f" t0 (inflection point) = {start_year + t0:.1f}") print(f" t0 (inflection point) = {start_year + t0:.1f}")
print(f" v (shape parameter) = {v:.3f}") print(f" v (shape parameter) = {v:.3f}")
print(f" R² = {r_squared:.4f}") print(f" R² = {r_squared:.4f}")
# Physical limit # Physical limit
physical_max = 3650 physical_max = 3650
print(f"\nPhysical limit: {physical_max} (10 sites × 365 days)") print(f"\nPhysical limit: {physical_max} (10 sites × 365 days)")
print(f"K / Physical limit = {K/physical_max:.2f}x") print(f"K / Physical limit = {K/physical_max:.2f}x")
# ========== Create Visualization ========== # ========== Create Visualization ==========
fig, ax = plt.subplots(figsize=(12, 7)) fig, ax = plt.subplots(figsize=(12, 7))
# Historical data points # Historical data points
ax.scatter(years, launches, color='#2C3E50', s=80, alpha=0.8, ax.scatter(years, launches, color='#2C3E50', s=80, alpha=0.8,
label='Historical Data (1984-2025)', zorder=3, edgecolor='white', linewidth=0.5) label='Historical Data (1984-2025)', zorder=3, edgecolor='white', linewidth=0.5)
# Generate smooth S-curve # Generate smooth S-curve
years_smooth = np.linspace(start_year, 2100, 500) years_smooth = np.linspace(start_year, 2100, 500)
t_smooth = years_smooth - start_year t_smooth = years_smooth - start_year
pred_smooth = richards(t_smooth, *popt) pred_smooth = richards(t_smooth, *popt)
# S-curve prediction # S-curve prediction
ax.plot(years_smooth, pred_smooth, color='#27AE60', lw=3, ax.plot(years_smooth, pred_smooth, color='#27AE60', lw=3,
label=f'Richards Model (K={K:.0f}, R²={r_squared:.3f})', zorder=2) label=f'Richards Model (K={K:.0f}, R²={r_squared:.3f})', zorder=2)
# Confidence band (approximate using parameter errors) # Confidence band (approximate using parameter errors)
K_low = max(500, K - 2*perr[0]) K_low = max(500, K - 2*perr[0])
K_high = K + 2*perr[0] K_high = K + 2*perr[0]
pred_low = richards(t_smooth, K_low, r, t0, v) pred_low = richards(t_smooth, K_low, r, t0, v)
pred_high = richards(t_smooth, K_high, 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, ax.fill_between(years_smooth, pred_low, pred_high, color='#27AE60', alpha=0.15,
label='95% Confidence Band') label='95% Confidence Band')
# Physical limit line # Physical limit line
ax.axhline(physical_max, color='#E74C3C', ls='--', lw=2.5, ax.axhline(physical_max, color='#E74C3C', ls='--', lw=2.5,
label=f'Physical Limit: {physical_max} (1/site/day)') label=f'Physical Limit: {physical_max} (1/site/day)')
# Mark inflection point # Mark inflection point
inflection_year = start_year + t0 inflection_year = start_year + t0
inflection_value = K / (v + 1) ** (1/v) inflection_value = K / (v + 1) ** (1/v)
ax.scatter([inflection_year], [inflection_value], color='#F39C12', s=150, ax.scatter([inflection_year], [inflection_value], color='#F39C12', s=150,
marker='*', zorder=5, label=f'Inflection Point ({inflection_year:.0f})') marker='*', zorder=5, label=f'Inflection Point ({inflection_year:.0f})')
# Mark key years # Mark key years
ax.axvline(2025, color='gray', ls=':', lw=1.5, alpha=0.7) ax.axvline(2025, color='gray', ls=':', lw=1.5, alpha=0.7)
ax.axvline(2050, color='#3498DB', ls=':', lw=2, alpha=0.8) 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(2026, K*0.95, '2025\n(Now)', fontsize=10, color='gray')
ax.text(2051, K*0.85, '2050\n(Target)', fontsize=10, color='#3498DB') ax.text(2051, K*0.85, '2050\n(Target)', fontsize=10, color='#3498DB')
# Prediction for 2050 # Prediction for 2050
t_2050 = 2050 - start_year t_2050 = 2050 - start_year
pred_2050 = richards(t_2050, *popt) pred_2050 = richards(t_2050, *popt)
ax.scatter([2050], [pred_2050], color='#3498DB', s=100, marker='D', zorder=4) 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), ax.annotate(f'{pred_2050:.0f}', xy=(2050, pred_2050), xytext=(2055, pred_2050-300),
fontsize=10, color='#3498DB', fontweight='bold') fontsize=10, color='#3498DB', fontweight='bold')
# Formatting # Formatting
ax.set_xlabel('Year', fontsize=12) ax.set_xlabel('Year', fontsize=12)
ax.set_ylabel('Annual Launches', fontsize=12) ax.set_ylabel('Annual Launches', fontsize=12)
ax.set_title('Richards S-Curve Model Fit (1984-2025 Data)\nRocket Launch Capacity Projection', ax.set_title('Richards S-Curve Model Fit (1984-2025 Data)\nRocket Launch Capacity Projection',
fontsize=14, fontweight='bold') fontsize=14, fontweight='bold')
ax.legend(loc='upper left', fontsize=10) ax.legend(loc='upper left', fontsize=10)
ax.grid(True, alpha=0.3) ax.grid(True, alpha=0.3)
ax.set_xlim(1982, 2100) ax.set_xlim(1982, 2100)
ax.set_ylim(0, K * 1.15) ax.set_ylim(0, K * 1.15)
# Add text box with model info # Add text box with model info
textstr = f'''Model: N(t) = K / (1 + e^(-r(t-t₀)))^(1/v) textstr = f'''Model: N(t) = K / (1 + e^(-r(t-t₀)))^(1/v)
Data Window: {start_year}-2025 ({len(data)} points) Data Window: {start_year}-2025 ({len(data)} points)
K = {K:.0f} launches/year K = {K:.0f} launches/year
r = {r:.4f} r = {r:.4f}
t₀ = {start_year + t0:.1f} t₀ = {start_year + t0:.1f}
v = {v:.3f} v = {v:.3f}
R² = {r_squared:.4f} R² = {r_squared:.4f}
Note: Physical limit (3,650) shown as Note: Physical limit (3,650) shown as
dashed red line''' dashed red line'''
props = dict(boxstyle='round', facecolor='wheat', alpha=0.8) props = dict(boxstyle='round', facecolor='wheat', alpha=0.8)
ax.text(0.98, 0.35, textstr, transform=ax.transAxes, fontsize=9, ax.text(0.98, 0.35, textstr, transform=ax.transAxes, fontsize=9,
verticalalignment='top', horizontalalignment='right', bbox=props, family='monospace') verticalalignment='top', horizontalalignment='right', bbox=props, family='monospace')
plt.tight_layout() plt.tight_layout()
plt.savefig('richards_curve_1984.png', dpi=150, bbox_inches='tight') plt.savefig('richards_curve_1984.png', dpi=150, bbox_inches='tight')
plt.close() plt.close()
print("\nPlot saved: richards_curve_1984.png") print("\nPlot saved: richards_curve_1984.png")
print("=" * 60) print("=" * 60)
if __name__ == "__main__": if __name__ == "__main__":
main() main()

3000
p3/water_sensitivity_analysis.py
View File

File diff suppressed because it is too large Load Diff

BIN
plane/__pycache__/cost_analysis.cpython-314.pyc Normal file
View File

Binary file not shown.

BIN
plane/__pycache__/cost_analysis.cpython-38.pyc Normal file
View File

Binary file not shown.

137
plane/cost_analysis.py
View File

@@ -556,6 +556,143 @@ def plot_comprehensive_summary(data, save_path='aviation_cost_summary.png'):
return fig 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__": if __name__ == "__main__":

BIN
plane/fuel_share_trend.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 196 KiB

After

Width:  |  Height:  |  Size: 138 KiB

BIN
plane/fuel_total_correlation.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 91 KiB

BIN
prob/p_B.pdf
View File

Binary file not shown.