786 lines
29 KiB
Python
786 lines
29 KiB
Python
#!/usr/bin/env python3
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
火箭发射数上限预测分析
|
||
对比 Logistic、Gompertz、Richards 三种增长模型
|
||
分析不同历史数据窗口对预测的影响
|
||
"""
|
||
|
||
import pandas as pd
|
||
import numpy as np
|
||
from scipy.optimize import curve_fit
|
||
from scipy.stats import pearsonr
|
||
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
|
||
|
||
|
||
# ============================================================
|
||
# 1. 增长模型定义
|
||
# ============================================================
|
||
|
||
def logistic(t, K, r, t0):
|
||
"""
|
||
Logistic 增长模型（对称S曲线）
|
||
N(t) = K / (1 + exp(-r*(t - t0)))
|
||
|
||
参数:
|
||
K : 环境容量/饱和上限
|
||
r : 增长率
|
||
t0 : 拐点时间（增长最快的时刻，此时 N = K/2）
|
||
"""
|
||
return K / (1 + np.exp(-r * (t - t0)))
|
||
|
||
|
||
def gompertz(t, K, b, c):
|
||
"""
|
||
Gompertz 增长模型（非对称S曲线，早期增长较快）
|
||
N(t) = K * exp(-b * exp(-c*t))
|
||
|
||
参数:
|
||
K : 饱和上限
|
||
b : 位移参数
|
||
c : 增长率参数
|
||
拐点在 N = K/e ≈ 0.368K 处
|
||
"""
|
||
return K * np.exp(-b * np.exp(-c * t))
|
||
|
||
|
||
def richards(t, K, r, t0, v):
|
||
"""
|
||
Richards 曲线（广义 Logistic，可调节对称性）
|
||
N(t) = K / (1 + exp(-r*(t - t0)))^(1/v)
|
||
|
||
参数:
|
||
K : 饱和上限
|
||
r : 增长率
|
||
t0 : 拐点时间
|
||
v : 形状参数（v=1 时退化为 Logistic）
|
||
"""
|
||
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)
|
||
|
||
|
||
# ============================================================
|
||
# 2. 数据加载与预处理
|
||
# ============================================================
|
||
|
||
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"])
|
||
# 只保留完整年份数据（截至2025年，2026年数据不完整）
|
||
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
|
||
|
||
|
||
# ============================================================
|
||
# 3. 模型拟合
|
||
# ============================================================
|
||
|
||
def fit_model(model_func, t, y, model_name, p0, bounds):
|
||
"""通用模型拟合函数"""
|
||
try:
|
||
popt, pcov = curve_fit(model_func, t, y, p0=p0, bounds=bounds, maxfev=50000)
|
||
perr = np.sqrt(np.diag(pcov))
|
||
y_pred = model_func(t, *popt)
|
||
|
||
# 计算拟合优度
|
||
ss_res = np.sum((y - y_pred) ** 2)
|
||
ss_tot = np.sum((y - np.mean(y)) ** 2)
|
||
r_squared = 1 - (ss_res / ss_tot)
|
||
rmse = np.sqrt(np.mean((y - y_pred) ** 2))
|
||
|
||
return {
|
||
"success": True,
|
||
"params": popt,
|
||
"param_err": perr,
|
||
"r_squared": r_squared,
|
||
"rmse": rmse,
|
||
"y_pred": y_pred
|
||
}
|
||
except Exception as e:
|
||
return {
|
||
"success": False,
|
||
"error": str(e)
|
||
}
|
||
|
||
|
||
def fit_all_models(data, base_year=1957):
|
||
"""对数据拟合所有三种模型"""
|
||
years = data["year"].values
|
||
launches = data["launches"].values
|
||
t = (years - base_year).astype(float)
|
||
|
||
results = {}
|
||
|
||
# Logistic 模型 - 提高上限边界以避免撞边界
|
||
p0_log = [2000.0, 0.08, 70.0]
|
||
bounds_log = ([300, 0.01, 0], [500000, 1.0, 300])
|
||
res = fit_model(logistic, t, launches, "Logistic", p0_log, bounds_log)
|
||
if res["success"]:
|
||
K, r, t0 = res["params"]
|
||
res["K"] = K
|
||
res["inflection_year"] = base_year + t0
|
||
res["inflection_value"] = K / 2
|
||
results["Logistic"] = res
|
||
|
||
# Gompertz 模型
|
||
p0_gom = [2000.0, 5.0, 0.03]
|
||
bounds_gom = ([300, 0.5, 0.005], [500000, 50.0, 0.5])
|
||
res = fit_model(gompertz, t, launches, "Gompertz", p0_gom, bounds_gom)
|
||
if res["success"]:
|
||
K, b, c = res["params"]
|
||
res["K"] = K
|
||
# Gompertz 拐点：t = ln(b)/c, 值为 K/e
|
||
res["inflection_year"] = base_year + np.log(b) / c
|
||
res["inflection_value"] = K / np.e
|
||
results["Gompertz"] = res
|
||
|
||
# Richards 模型
|
||
p0_ric = [2000.0, 0.08, 70.0, 2.0]
|
||
bounds_ric = ([300, 0.01, 0, 0.1], [500000, 1.0, 300, 10.0])
|
||
res = fit_model(richards, t, launches, "Richards", p0_ric, bounds_ric)
|
||
if res["success"]:
|
||
K, r, t0, v = res["params"]
|
||
res["K"] = K
|
||
res["inflection_year"] = base_year + t0
|
||
res["inflection_value"] = K / (v + 1) ** (1/v)
|
||
res["v"] = v
|
||
results["Richards"] = res
|
||
|
||
return results, years, launches, t, base_year
|
||
|
||
|
||
def predict_future(results, base_year, years_future):
|
||
"""生成未来预测"""
|
||
t_future = years_future - base_year
|
||
predictions = {}
|
||
|
||
for model_name, res in results.items():
|
||
if not res["success"]:
|
||
continue
|
||
|
||
params = res["params"]
|
||
if model_name == "Logistic":
|
||
pred = logistic(t_future, *params)
|
||
elif model_name == "Gompertz":
|
||
pred = gompertz(t_future, *params)
|
||
elif model_name == "Richards":
|
||
pred = richards(t_future, *params)
|
||
|
||
predictions[model_name] = np.maximum(pred, 0)
|
||
|
||
return predictions
|
||
|
||
|
||
# ============================================================
|
||
# 4. 可视化
|
||
# ============================================================
|
||
|
||
def plot_model_comparison(df, results_dict, save_path="launch_capacity_analysis.png"):
|
||
"""绘制模型对比图"""
|
||
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
|
||
|
||
# 颜色方案
|
||
colors = {
|
||
"Logistic": "#E74C3C",
|
||
"Gompertz": "#3498DB",
|
||
"Richards": "#27AE60"
|
||
}
|
||
|
||
# 子图1: 历史数据 + 三种模型拟合（全部数据）
|
||
ax1 = axes[0, 0]
|
||
results, years, launches, t, base_year = results_dict["全部数据"]
|
||
ax1.scatter(years, launches, color="gray", s=30, alpha=0.7, label="历史数据", zorder=3)
|
||
|
||
years_curve = np.arange(years.min(), 2101)
|
||
t_curve = years_curve - base_year
|
||
|
||
for model_name, res in results.items():
|
||
if not res["success"]:
|
||
continue
|
||
params = res["params"]
|
||
if model_name == "Logistic":
|
||
pred = logistic(t_curve, *params)
|
||
elif model_name == "Gompertz":
|
||
pred = gompertz(t_curve, *params)
|
||
elif model_name == "Richards":
|
||
pred = richards(t_curve, *params)
|
||
|
||
ax1.plot(years_curve, pred, color=colors[model_name], lw=2,
|
||
label=f"{model_name} (K={res['K']:.0f}, R²={res['r_squared']:.3f})")
|
||
|
||
ax1.axvline(2025, color="gray", ls="--", lw=1, alpha=0.5)
|
||
ax1.set_xlabel("年份", fontsize=11)
|
||
ax1.set_ylabel("年发射次数", fontsize=11)
|
||
ax1.set_title("全部历史数据 (1957-2025) 模型拟合", fontsize=12)
|
||
ax1.legend(loc="upper left", fontsize=9)
|
||
ax1.grid(True, alpha=0.3)
|
||
ax1.set_xlim(1955, 2105)
|
||
ax1.set_ylim(0, None)
|
||
|
||
# 子图2: 近15年数据拟合
|
||
ax2 = axes[0, 1]
|
||
results, years, launches, t, base_year = results_dict["近15年"]
|
||
ax2.scatter(years, launches, color="gray", s=40, alpha=0.8, label="历史数据", zorder=3)
|
||
|
||
years_curve = np.arange(2010, 2101)
|
||
t_curve = years_curve - base_year
|
||
|
||
for model_name, res in results.items():
|
||
if not res["success"]:
|
||
continue
|
||
params = res["params"]
|
||
if model_name == "Logistic":
|
||
pred = logistic(t_curve, *params)
|
||
elif model_name == "Gompertz":
|
||
pred = gompertz(t_curve, *params)
|
||
elif model_name == "Richards":
|
||
pred = richards(t_curve, *params)
|
||
|
||
ax2.plot(years_curve, pred, color=colors[model_name], lw=2,
|
||
label=f"{model_name} (K={res['K']:.0f})")
|
||
|
||
ax2.axvline(2025, color="gray", ls="--", lw=1, alpha=0.5)
|
||
ax2.set_xlabel("年份", fontsize=11)
|
||
ax2.set_ylabel("年发射次数", fontsize=11)
|
||
ax2.set_title("近15年数据 (2010-2025) 模型拟合", fontsize=12)
|
||
ax2.legend(loc="upper left", fontsize=9)
|
||
ax2.grid(True, alpha=0.3)
|
||
ax2.set_xlim(2008, 2105)
|
||
ax2.set_ylim(0, None)
|
||
|
||
# 子图3: 只展示 Richards 模型的K值（因为其他模型撞边界）
|
||
ax3 = axes[1, 0]
|
||
windows = ["全部数据", "近25年", "近15年", "近10年"]
|
||
x = np.arange(len(windows))
|
||
|
||
k_values = []
|
||
for window in windows:
|
||
if window in results_dict:
|
||
res = results_dict[window][0].get("Richards", {})
|
||
k_values.append(res.get("K", 0) if res.get("success", False) else 0)
|
||
else:
|
||
k_values.append(0)
|
||
|
||
bars = ax3.bar(x, k_values, 0.6, color=colors["Richards"], alpha=0.8, edgecolor='black')
|
||
|
||
# 添加数值标签
|
||
for bar, k in zip(bars, k_values):
|
||
if k > 0:
|
||
ax3.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 500,
|
||
f'{k:.0f}', ha='center', va='bottom', fontsize=10, fontweight='bold')
|
||
|
||
ax3.set_xlabel("数据窗口", fontsize=11)
|
||
ax3.set_ylabel("预测饱和上限 K (次/年)", fontsize=11)
|
||
ax3.set_title("Richards 模型: 不同数据窗口的饱和上限预测\n(Logistic/Gompertz 撞边界，不适用)", fontsize=11)
|
||
ax3.set_xticks(x)
|
||
ax3.set_xticklabels(windows)
|
||
ax3.grid(True, alpha=0.3, axis='y')
|
||
ax3.set_ylim(0, max(k_values) * 1.15 if k_values else 10000)
|
||
|
||
# 子图4: Richards 模型拐点年份
|
||
ax4 = axes[1, 1]
|
||
inflection_years = []
|
||
inflection_values = []
|
||
for window in windows:
|
||
if window in results_dict:
|
||
res = results_dict[window][0].get("Richards", {})
|
||
if res.get("success", False):
|
||
inflection_years.append(res.get("inflection_year", 0))
|
||
inflection_values.append(res.get("inflection_value", 0))
|
||
else:
|
||
inflection_years.append(0)
|
||
inflection_values.append(0)
|
||
else:
|
||
inflection_years.append(0)
|
||
inflection_values.append(0)
|
||
|
||
bars = ax4.bar(x, inflection_years, 0.6, color=colors["Richards"], alpha=0.8, edgecolor='black')
|
||
|
||
# 添加标签（拐点年份和对应发射数）
|
||
for bar, year, val in zip(bars, inflection_years, inflection_values):
|
||
if year > 0:
|
||
ax4.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.5,
|
||
f'{year:.0f}\n({val:.0f}次/年)', ha='center', va='bottom', fontsize=9)
|
||
|
||
ax4.axhline(2025, color="red", ls="--", lw=2, alpha=0.8, label="当前年份 (2025)")
|
||
ax4.set_xlabel("数据窗口", fontsize=11)
|
||
ax4.set_ylabel("拐点年份", fontsize=11)
|
||
ax4.set_title("Richards 模型: 预测拐点年份\n(增长最快的时刻)", fontsize=11)
|
||
ax4.set_xticks(x)
|
||
ax4.set_xticklabels(windows)
|
||
ax4.legend(loc="lower right", fontsize=9)
|
||
ax4.grid(True, alpha=0.3, axis='y')
|
||
valid_years = [y for y in inflection_years if y > 0]
|
||
if valid_years:
|
||
ax4.set_ylim(2020, max(valid_years) + 15)
|
||
|
||
plt.suptitle("火箭年发射数上限预测 - 多模型对比分析", fontsize=14, fontweight='bold', y=1.02)
|
||
plt.tight_layout()
|
||
plt.savefig(save_path, dpi=150, bbox_inches="tight")
|
||
plt.close()
|
||
print(f"图表已保存: {save_path}")
|
||
|
||
|
||
def plot_physical_vs_model(results_dict, physical_limits, save_path="launch_capacity_physical.png"):
|
||
"""绘制物理约束与模型预测对比图"""
|
||
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
|
||
|
||
# 左图：物理约束与模型预测K值对比
|
||
ax1 = axes[0]
|
||
|
||
# 物理约束
|
||
phys_scenarios = list(physical_limits.keys())
|
||
phys_values = list(physical_limits.values())
|
||
x_phys = np.arange(len(phys_scenarios))
|
||
bars1 = ax1.barh(x_phys, phys_values, 0.6, color='#3498DB', alpha=0.7, label='物理约束估计')
|
||
|
||
# 添加数值标签
|
||
for bar, val in zip(bars1, phys_values):
|
||
ax1.text(val + 50, bar.get_y() + bar.get_height()/2,
|
||
f'{val:,}', va='center', fontsize=10)
|
||
|
||
# 添加模型预测（Richards）
|
||
richards_k = []
|
||
for window in ["近25年", "近15年", "近10年"]:
|
||
if window in results_dict:
|
||
res = results_dict[window][0].get("Richards", {})
|
||
if res.get("success", False):
|
||
richards_k.append(res["K"])
|
||
|
||
if richards_k:
|
||
avg_richards = np.mean(richards_k)
|
||
ax1.axvline(avg_richards, color='#E74C3C', ls='--', lw=2,
|
||
label=f'Richards模型平均 K={avg_richards:,.0f}')
|
||
|
||
ax1.set_yticks(x_phys)
|
||
ax1.set_yticklabels(phys_scenarios)
|
||
ax1.set_xlabel("年发射上限 (次/年)", fontsize=11)
|
||
ax1.set_title("物理约束 vs 模型预测 (10个发射场)", fontsize=12)
|
||
ax1.legend(loc="lower right", fontsize=9)
|
||
ax1.grid(True, alpha=0.3, axis='x')
|
||
ax1.set_xlim(0, max(phys_values) * 1.3)
|
||
|
||
# 右图：带物理约束的时间预测
|
||
ax2 = axes[1]
|
||
|
||
# 获取近15年Richards预测
|
||
if "近15年" in results_dict:
|
||
results, years, launches, t, base_year = results_dict["近15年"]
|
||
res = results.get("Richards", {})
|
||
|
||
if res["success"]:
|
||
years_future = np.arange(2010, 2101)
|
||
t_future = years_future - base_year
|
||
pred = richards(t_future, *res["params"])
|
||
|
||
ax2.scatter(years, launches, color="gray", s=40, alpha=0.7, label="历史数据", zorder=3)
|
||
ax2.plot(years_future, pred, color='#27AE60', lw=2, label=f"Richards (K={res['K']:,.0f})")
|
||
|
||
# 添加物理约束线
|
||
phys_colors = ['#E74C3C', '#F39C12', '#3498DB', '#9B59B6', '#1ABC9C']
|
||
for i, (scenario, limit) in enumerate(physical_limits.items()):
|
||
ax2.axhline(limit, color=phys_colors[i % len(phys_colors)],
|
||
ls='--', lw=1.5, alpha=0.7, label=f'{scenario}: {limit:,}')
|
||
|
||
ax2.axvline(2025, color="gray", ls=":", lw=1, alpha=0.5)
|
||
ax2.set_xlabel("年份", fontsize=11)
|
||
ax2.set_ylabel("年发射次数", fontsize=11)
|
||
ax2.set_title("Richards预测 + 物理约束线", fontsize=12)
|
||
ax2.legend(loc="upper left", fontsize=8)
|
||
ax2.grid(True, alpha=0.3)
|
||
ax2.set_xlim(2008, 2105)
|
||
ax2.set_ylim(0, max(physical_limits.values()) * 1.2)
|
||
|
||
plt.suptitle("物理约束分析 (10个发射场)", fontsize=14, fontweight='bold', y=1.02)
|
||
plt.tight_layout()
|
||
plt.savefig(save_path, dpi=150, bbox_inches="tight")
|
||
plt.close()
|
||
print(f"图表已保存: {save_path}")
|
||
|
||
|
||
def plot_future_predictions(results_dict, save_path="launch_capacity_predictions.png"):
|
||
"""绘制未来预测图"""
|
||
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
|
||
|
||
colors = {
|
||
"Logistic": "#E74C3C",
|
||
"Gompertz": "#3498DB",
|
||
"Richards": "#27AE60"
|
||
}
|
||
|
||
# 左图：基于全部数据的预测
|
||
ax1 = axes[0]
|
||
results, years, launches, t, base_year = results_dict["全部数据"]
|
||
years_future = np.arange(1957, 2151)
|
||
predictions = predict_future(results, base_year, years_future)
|
||
|
||
ax1.scatter(years, launches, color="gray", s=20, alpha=0.6, label="历史数据", zorder=3)
|
||
for model_name, pred in predictions.items():
|
||
ax1.plot(years_future, pred, color=colors[model_name], lw=2, label=model_name)
|
||
|
||
ax1.axvline(2025, color="gray", ls="--", lw=1, alpha=0.5)
|
||
ax1.axvspan(2025, 2150, alpha=0.05, color="blue")
|
||
ax1.set_xlabel("年份", fontsize=11)
|
||
ax1.set_ylabel("年发射次数", fontsize=11)
|
||
ax1.set_title("基于全部历史数据的预测 (至2150年)", fontsize=12)
|
||
ax1.legend(loc="upper left", fontsize=10)
|
||
ax1.grid(True, alpha=0.3)
|
||
ax1.set_xlim(1955, 2155)
|
||
ax1.set_ylim(0, None)
|
||
|
||
# 右图：基于近15年数据的预测
|
||
ax2 = axes[1]
|
||
results, years, launches, t, base_year = results_dict["近15年"]
|
||
years_future = np.arange(2010, 2151)
|
||
predictions = predict_future(results, base_year, years_future)
|
||
|
||
ax2.scatter(years, launches, color="gray", s=40, alpha=0.7, label="历史数据", zorder=3)
|
||
for model_name, pred in predictions.items():
|
||
ax2.plot(years_future, pred, color=colors[model_name], lw=2, label=model_name)
|
||
|
||
ax2.axvline(2025, color="gray", ls="--", lw=1, alpha=0.5)
|
||
ax2.axvspan(2025, 2150, alpha=0.05, color="blue")
|
||
ax2.set_xlabel("年份", fontsize=11)
|
||
ax2.set_ylabel("年发射次数", fontsize=11)
|
||
ax2.set_title("基于近15年数据的预测 (至2150年)", fontsize=12)
|
||
ax2.legend(loc="upper left", fontsize=10)
|
||
ax2.grid(True, alpha=0.3)
|
||
ax2.set_xlim(2008, 2155)
|
||
ax2.set_ylim(0, None)
|
||
|
||
plt.suptitle("火箭年发射数长期预测", fontsize=14, fontweight='bold', y=1.02)
|
||
plt.tight_layout()
|
||
plt.savefig(save_path, dpi=150, bbox_inches="tight")
|
||
plt.close()
|
||
print(f"图表已保存: {save_path}")
|
||
|
||
|
||
# ============================================================
|
||
# 5. 主程序
|
||
# ============================================================
|
||
|
||
def bootstrap_K_estimate(data, model_func, p0, bounds, n_bootstrap=500, base_year=1957):
|
||
"""Bootstrap 估计 K 值的不确定性"""
|
||
years = data["year"].values
|
||
launches = data["launches"].values
|
||
t = (years - base_year).astype(float)
|
||
n = len(t)
|
||
|
||
K_samples = []
|
||
for _ in range(n_bootstrap):
|
||
# 重采样
|
||
idx = np.random.choice(n, n, replace=True)
|
||
t_boot = t[idx]
|
||
y_boot = launches[idx]
|
||
|
||
try:
|
||
popt, _ = curve_fit(model_func, t_boot, y_boot, p0=p0, bounds=bounds, maxfev=10000)
|
||
K_samples.append(popt[0]) # K 是第一个参数
|
||
except:
|
||
pass
|
||
|
||
if len(K_samples) > 10:
|
||
return {
|
||
"mean": np.mean(K_samples),
|
||
"median": np.median(K_samples),
|
||
"std": np.std(K_samples),
|
||
"ci_5": np.percentile(K_samples, 5),
|
||
"ci_95": np.percentile(K_samples, 95),
|
||
"samples": K_samples
|
||
}
|
||
return None
|
||
|
||
|
||
def calculate_aic_bic(y_true, y_pred, n_params):
|
||
"""计算 AIC 和 BIC"""
|
||
n = len(y_true)
|
||
ss_res = np.sum((y_true - y_pred) ** 2)
|
||
mse = ss_res / n
|
||
|
||
# Log-likelihood (假设高斯误差)
|
||
log_likelihood = -n/2 * (np.log(2 * np.pi) + np.log(mse) + 1)
|
||
|
||
aic = 2 * n_params - 2 * log_likelihood
|
||
bic = n_params * np.log(n) - 2 * log_likelihood
|
||
|
||
return aic, bic
|
||
|
||
|
||
def main():
|
||
print("=" * 70)
|
||
print("火箭年发射数上限预测分析")
|
||
print("模型: Logistic, Gompertz, Richards")
|
||
print("=" * 70)
|
||
|
||
# 加载数据
|
||
df = load_data()
|
||
print(f"\n数据范围: {df['year'].min()} - {df['year'].max()}")
|
||
print(f"数据点数: {len(df)}")
|
||
print(f"最近5年发射数:")
|
||
for _, row in df.tail(5).iterrows():
|
||
print(f" {int(row['year'])}: {int(row['launches'])} 次")
|
||
|
||
# 计算增长率
|
||
recent = df[df["year"] >= 2020].copy()
|
||
if len(recent) > 1:
|
||
growth_rates = recent["launches"].pct_change().dropna() * 100
|
||
print(f"\n2020年以来年均增长率: {growth_rates.mean():.1f}%")
|
||
|
||
# 定义不同的数据窗口
|
||
end_year = 2025
|
||
windows = {
|
||
"全部数据": df[df["year"] <= end_year].copy(),
|
||
"近25年": df[(df["year"] >= 2000) & (df["year"] <= end_year)].copy(),
|
||
"近15年": df[(df["year"] >= 2010) & (df["year"] <= end_year)].copy(),
|
||
"近10年": df[(df["year"] >= 2015) & (df["year"] <= end_year)].copy(),
|
||
}
|
||
|
||
# 存储所有结果
|
||
all_results = {}
|
||
|
||
# 对每个窗口进行拟合
|
||
for window_name, data in windows.items():
|
||
print(f"\n{'='*50}")
|
||
print(f"数据窗口: {window_name} ({data['year'].min()}-{data['year'].max()}, n={len(data)})")
|
||
print("=" * 50)
|
||
|
||
results, years, launches, t, base_year = fit_all_models(data)
|
||
all_results[window_name] = (results, years, launches, t, base_year)
|
||
|
||
for model_name, res in results.items():
|
||
if res["success"]:
|
||
print(f"\n{model_name} 模型:")
|
||
print(f" 饱和上限 K = {res['K']:.0f} 次/年")
|
||
print(f" 拐点年份 = {res['inflection_year']:.1f}")
|
||
print(f" 拐点处发射数 = {res['inflection_value']:.0f} 次/年")
|
||
print(f" R² (拟合优度) = {res['r_squared']:.4f}")
|
||
print(f" RMSE = {res['rmse']:.2f}")
|
||
if model_name == "Richards":
|
||
print(f" 形状参数 v = {res.get('v', 'N/A'):.3f}")
|
||
else:
|
||
print(f"\n{model_name} 模型: 拟合失败 - {res.get('error', 'Unknown')}")
|
||
|
||
# 生成未来预测表
|
||
print("\n" + "=" * 70)
|
||
print("未来预测 (基于近15年数据)")
|
||
print("=" * 70)
|
||
|
||
results, years, launches, t, base_year = all_results["近15年"]
|
||
future_years = np.array([2030, 2040, 2050, 2075, 2100])
|
||
predictions = predict_future(results, base_year, future_years)
|
||
|
||
print(f"\n{'年份':<10}", end="")
|
||
for model in ["Logistic", "Gompertz", "Richards"]:
|
||
print(f"{model:<15}", end="")
|
||
print()
|
||
print("-" * 55)
|
||
|
||
for i, year in enumerate(future_years):
|
||
print(f"{year:<10}", end="")
|
||
for model in ["Logistic", "Gompertz", "Richards"]:
|
||
if model in predictions:
|
||
print(f"{predictions[model][i]:<15.0f}", end="")
|
||
else:
|
||
print(f"{'N/A':<15}", end="")
|
||
print()
|
||
|
||
# 饱和上限汇总
|
||
print("\n" + "=" * 70)
|
||
print("饱和上限 K 汇总 (次/年)")
|
||
print("=" * 70)
|
||
print(f"\n{'数据窗口':<20}", end="")
|
||
for model in ["Logistic", "Gompertz", "Richards"]:
|
||
print(f"{model:<15}", end="")
|
||
print()
|
||
print("-" * 65)
|
||
|
||
for window_name in ["全部数据", "近25年", "近15年", "近10年"]:
|
||
results = all_results[window_name][0]
|
||
print(f"{window_name:<20}", end="")
|
||
for model in ["Logistic", "Gompertz", "Richards"]:
|
||
if model in results and results[model]["success"]:
|
||
print(f"{results[model]['K']:<15.0f}", end="")
|
||
else:
|
||
print(f"{'N/A':<15}", end="")
|
||
print()
|
||
|
||
# 物理约束分析
|
||
physical_limits = analyze_physical_constraints(n_sites=10)
|
||
|
||
# 绘制图表
|
||
print("\n" + "=" * 70)
|
||
print("生成图表...")
|
||
plot_model_comparison(df, all_results)
|
||
plot_future_predictions(all_results)
|
||
plot_physical_vs_model(all_results, physical_limits)
|
||
|
||
# 保存详细结果到CSV
|
||
save_results_csv(all_results)
|
||
|
||
print("\n" + "=" * 70)
|
||
print("分析完成!")
|
||
print("=" * 70)
|
||
|
||
# Bootstrap 不确定性分析
|
||
print("\n" + "=" * 70)
|
||
print("Bootstrap 不确定性分析 (基于近15年数据, 500次重采样)")
|
||
print("=" * 70)
|
||
|
||
data_15 = windows["近15年"]
|
||
base_year = 1957
|
||
|
||
# Richards 模型 Bootstrap
|
||
p0_ric = [2000.0, 0.08, 70.0, 2.0]
|
||
bounds_ric = ([300, 0.01, 0, 0.1], [500000, 1.0, 300, 10.0])
|
||
|
||
boot_result = bootstrap_K_estimate(data_15, richards, p0_ric, bounds_ric, n_bootstrap=500, base_year=base_year)
|
||
if boot_result:
|
||
print(f"\nRichards 模型 K 值 Bootstrap 估计:")
|
||
print(f" 均值: {boot_result['mean']:.0f} 次/年")
|
||
print(f" 中位数: {boot_result['median']:.0f} 次/年")
|
||
print(f" 标准差: {boot_result['std']:.0f}")
|
||
print(f" 90% 置信区间: [{boot_result['ci_5']:.0f}, {boot_result['ci_95']:.0f}] 次/年")
|
||
|
||
# AIC/BIC 模型选择
|
||
print("\n" + "=" * 70)
|
||
print("模型选择 (AIC/BIC, 基于近15年数据)")
|
||
print("=" * 70)
|
||
|
||
results_15, years_15, launches_15, t_15, _ = all_results["近15年"]
|
||
print(f"\n{'模型':<12} {'AIC':<12} {'BIC':<12} {'参数数':<8}")
|
||
print("-" * 44)
|
||
|
||
for model_name, res in results_15.items():
|
||
if res["success"]:
|
||
n_params = len(res["params"])
|
||
aic, bic = calculate_aic_bic(launches_15, res["y_pred"], n_params)
|
||
print(f"{model_name:<12} {aic:<12.1f} {bic:<12.1f} {n_params:<8}")
|
||
|
||
# 结论
|
||
print("\n" + "=" * 70)
|
||
print("📊 分析结论与建议")
|
||
print("=" * 70)
|
||
|
||
# 收集有效的K值（排除撞边界的）
|
||
valid_k_values = []
|
||
for window_name, (results, _, _, _, _) in all_results.items():
|
||
for model_name, res in results.items():
|
||
if res["success"] and res["K"] < 100000: # 排除撞边界的
|
||
valid_k_values.append((window_name, model_name, res["K"]))
|
||
|
||
print("\n1. 模型适用性分析:")
|
||
print("-" * 50)
|
||
print(" • Logistic/Gompertz: 数据处于快速增长期，无法确定上限")
|
||
print(" (模型撞到边界，说明增长远未饱和)")
|
||
print(" • Richards 曲线: 额外的形状参数使其能给出更合理估计")
|
||
print(" (推荐用于当前阶段的上限预测)")
|
||
|
||
print("\n2. 合理的上限估计 (Richards 模型):")
|
||
print("-" * 50)
|
||
for window_name, model_name, k in valid_k_values:
|
||
if model_name == "Richards" and window_name != "全部数据":
|
||
print(f" • {window_name}: K ≈ {k:.0f} 次/年")
|
||
|
||
# Richards 结果的平均值
|
||
richards_k = [k for w, m, k in valid_k_values if m == "Richards" and w != "全部数据"]
|
||
if richards_k:
|
||
print(f"\n 📌 综合估计: {int(np.mean(richards_k)):,} - {int(np.max(richards_k)):,} 次/年")
|
||
|
||
print("\n3. 物理约束参考 (基于10个主要发射场):")
|
||
print("-" * 50)
|
||
print(" • 主要发射场数量: 10个")
|
||
print(" • 假设每发射场年最大容量:")
|
||
print(" - 保守估计: 50次/年 → 总上限 500 次/年")
|
||
print(" - 中等估计: 100次/年 → 总上限 1,000 次/年")
|
||
print(" - 乐观估计: 200次/年 → 总上限 2,000 次/年")
|
||
print(" - 极限估计 (SpaceX级别): 400次/年 → 总上限 4,000 次/年")
|
||
print(" • 2025年实际数据: 329次 (10发射场平均 33次/场)")
|
||
print(" • SpaceX Starship 单型号目标: 1,000+ 次/年")
|
||
|
||
print("\n4. 预测不确定性:")
|
||
print("-" * 50)
|
||
print(" • 当前处于指数增长早期，上限估计高度不确定")
|
||
print(" • 建议每年用新数据更新预测")
|
||
print(" • 需关注技术突破（可重复使用、星舰等）对上限的影响")
|
||
|
||
if boot_result:
|
||
print(f"\n Bootstrap 90%置信区间: {boot_result['ci_5']:.0f} - {boot_result['ci_95']:.0f} 次/年")
|
||
|
||
|
||
def analyze_physical_constraints(n_sites=10):
|
||
"""
|
||
基于发射场物理约束的上限分析
|
||
|
||
参数:
|
||
n_sites: 发射场数量
|
||
"""
|
||
print("\n" + "=" * 70)
|
||
print(f"物理约束分析 (基于 {n_sites} 个发射场)")
|
||
print("=" * 70)
|
||
|
||
# 不同情景下的单发射场年容量
|
||
scenarios = {
|
||
"保守 (当前技术)": 50,
|
||
"中等 (技术进步)": 100,
|
||
"乐观 (快速周转)": 200,
|
||
"极限 (SpaceX级)": 400,
|
||
"理论极限 (每天1发)": 365,
|
||
}
|
||
|
||
print(f"\n单发射场年发射容量假设:")
|
||
print("-" * 50)
|
||
|
||
results = {}
|
||
for scenario, capacity_per_site in scenarios.items():
|
||
total_capacity = n_sites * capacity_per_site
|
||
results[scenario] = total_capacity
|
||
print(f" {scenario:<20}: {capacity_per_site:>3} 次/场/年 → 总上限 {total_capacity:>5,} 次/年")
|
||
|
||
# 带置信度的综合估计
|
||
print(f"\n综合估计 ({n_sites} 发射场):")
|
||
print("-" * 50)
|
||
print(f" 50% 置信区间: {n_sites * 50:,} - {n_sites * 100:,} 次/年")
|
||
print(f" 90% 置信区间: {n_sites * 30:,} - {n_sites * 200:,} 次/年")
|
||
print(f" 理论最大值: {n_sites * 365:,} 次/年 (每天每场1发)")
|
||
|
||
return results
|
||
|
||
|
||
def save_results_csv(all_results, filepath="launch_capacity_results.csv"):
|
||
"""保存详细结果到CSV"""
|
||
rows = []
|
||
for window_name, (results, _, _, _, _) in all_results.items():
|
||
for model_name, res in results.items():
|
||
if res["success"]:
|
||
row = {
|
||
"数据窗口": window_name,
|
||
"模型": model_name,
|
||
"饱和上限K": res["K"],
|
||
"拐点年份": res["inflection_year"],
|
||
"拐点发射数": res["inflection_value"],
|
||
"R²": res["r_squared"],
|
||
"RMSE": res["rmse"]
|
||
}
|
||
rows.append(row)
|
||
|
||
df_results = pd.DataFrame(rows)
|
||
df_results.to_csv(filepath, index=False, encoding="utf-8-sig")
|
||
print(f"详细结果已保存: {filepath}")
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|