#!/usr/bin/env python3 # -*- coding: utf-8 -*- """确定性合成数据生成器（根因修复：持久化 data，避免配图/附录依赖丢失）。 所有赛题用固定随机种子生成，保证可复现；运行一次即写出 assets/problems/data/*.csv。 当前实现：cumcm2012a（葡萄酒评价）。其余赛题逐步回填。 """ import os, csv, math, random, statistics ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) OUT = os.path.join(ROOT, "assets", "problems", "data") os.makedirs(OUT, exist_ok=True) def clip(x, lo, hi): return max(lo, min(hi, x)) def gen_2012a(seed=2012, N=27): """葡萄酒评价：N 种红葡萄酒 × 10 项理化指标(效益型) + 专家评分。 专家评分由隐含真实权重对标准化理化线性组合 + 噪声生成，与综合评价高度相关， 便于正文做「模型结论 vs 专家评分」一致性验证。""" random.seed(seed) # 指标：(名称, 均值, 标准差, 下限, 上限) specs = [ ("花色苷", 15.0, 4.0, 5.0, 30.0), ("单宁", 2.5, 0.6, 1.0, 4.0), ("总酚", 30.0, 8.0, 10.0, 50.0), ("酒总黄酮", 20.0, 5.0, 8.0, 35.0), ("白藜芦醇", 3.0, 1.0, 0.5, 6.0), ("DPPH", 80.0, 20.0, 30.0, 140.0), ("L星", 50.0, 8.0, 30.0, 70.0), ("a星", 25.0, 6.0, 5.0, 45.0), ("b星", 10.0, 4.0, 0.0, 25.0), ("干物质", 25.0, 5.0, 12.0, 40.0), ] names = [s[0] for s in specs] M = [] for i in range(N): row = [clip(random.gauss(s[1], s[2]), s[3], s[4]) for s in specs] M.append(row) # 标准化（z-score）按列 means = [sum(M[i][j] for i in range(N)) / N for j in range(len(names))] stds = [math.sqrt(sum((M[i][j] - means[j]) ** 2 for i in range(N)) / N) or 1 for j in range(len(names))] Z = [[(M[i][j] - means[j]) / stds[j] for j in range(len(names))] for i in range(N)] # 隐含真实权重（仅用于生成专家评分，正文中不暴露） true_w = [0.15, 0.12, 0.15, 0.10, 0.08, 0.12, 0.06, 0.10, 0.05, 0.07] comp = [sum(true_w[j] * Z[i][j] for j in range(len(names))) for i in range(N)] cmin, cmax = min(comp), max(comp) expert = [] for i in range(N): base = 60 + 36 * (comp[i] - cmin) / (cmax - cmin or 1) base += random.gauss(0, 3.0) # 评分噪声 expert.append(round(clip(base, 60, 99), 1)) # 两组评酒员评分：组1 更稳定（更可信），组2 噪声更大（较不可信）。 # 以隐含真实品质 comp 为锚，叠加不同强度噪声，使结论可判定"组1 更可信"。 g1, g2 = [], [] for i in range(N): base = 60 + 36 * (comp[i] - cmin) / (cmax - cmin or 1) g1.append(round(clip(base + random.gauss(0, 2.0), 60, 99), 1)) # 组1 噪声小 g2.append(round(clip(base + random.gauss(0, 4.5), 60, 99), 1)) # 组2 噪声大 # 写出 csv path = os.path.join(OUT, "cumcm2012a.csv") with open(path, "w", encoding="utf-8-sig", newline="") as f: w = csv.writer(f) w.writerow(["酒样"] + names + ["专家评分", "组1评分", "组2评分"]) for i in range(N): w.writerow(["酒样%02d" % (i + 1)] + [round(M[i][j], 2) for j in range(len(names))] + [expert[i], g1[i], g2[i]]) # 诊断打印 print("cumcm2012a.csv 已写出：%d 酒样 × %d 理化指标 + 专家/组1/组2评分" % (N, len(names))) # 验证：综合(等权z)vs专家评分 Spearman def spearman(r1, r2): n = len(r1) rk = lambda r: [sorted(range(n), key=lambda k: -r[k]).index(k) + 1 for k in range(n)] a, b = rk(r1), rk(r2); da, db = sum(a) / n, sum(b) / n cov = sum((a[i] - da) * (b[i] - db) for i in range(n)) / n sa = math.sqrt(sum((a[i] - da) ** 2 for i in range(n)) / n) sb = math.sqrt(sum((b[i] - db) ** 2 for i in range(n)) / n) return cov / (sa * sb) eq = [sum(Z[i]) / len(names) for i in range(N)] rho = spearman(eq, expert) print("诊断：等权综合 vs 专家评分 Spearman = %.3f（应较高，验证生成有效）" % rho) print("专家评分范围：%.1f ~ %.1f" % (min(expert), max(expert))) # 诊断两组：组间相关与组内变异系数（组1 应更高相关、更低 CV） print("组1 vs 组2 Spearman = %.3f" % spearman(g1, g2)) cv1 = math.sqrt(sum((g1[i] - sum(g1)/N)**2 for i in range(N))/N)/(sum(g1)/N) cv2 = math.sqrt(sum((g2[i] - sum(g2)/N)**2 for i in range(N))/N)/(sum(g2)/N) print("组1 CV=%.3f 组2 CV=%.3f" % (cv1, cv2)) def gen_2013a(seed=2013, N=200): """车道线偏离二分类（label: 0=正常, 1=偏离）。 f2(横向偏移估计)、f6(线置信度) 为强判别特征；f1/f3/f5 噪声；f4(曲率)弱相关。 固定种子保证可复现，重跑得到完全相同的 CSV。""" rng = random.Random(seed) rows = [] for _ in range(N): label = 1 if rng.random() < 0.5 else 0 if label == 1: # 偏离 f1 = rng.gauss(0.0, 0.8) f2 = rng.gauss(1.30, 0.45) f3 = rng.gauss(0.0, 0.7) f4 = rng.gauss(0.60, 0.42) f5 = rng.gauss(0.0, 0.9) f6 = rng.gauss(0.55, 0.30) else: # 正常 f1 = rng.gauss(0.0, 0.8) f2 = rng.gauss(0.40, 0.40) f3 = rng.gauss(0.0, 0.7) f4 = rng.gauss(0.12, 0.40) f5 = rng.gauss(0.0, 0.9) f6 = rng.gauss(1.20, 0.30) rows.append([round(f1, 3), round(f2, 3), round(f3, 3), round(f4, 3), round(f5, 3), round(f6, 3), label]) # shuffle 使标签顺序自然 rng.shuffle(rows) path = os.path.join(OUT, "cumcm2013a.csv") with open(path, "w", encoding="utf-8-sig", newline="") as f: w = csv.writer(f) w.writerow(["f1", "f2", "f3", "f4", "f5", "f6", "label"]) w.writerows(rows) g0 = [r for r in rows if r[6] == 0]; g1 = [r for r in rows if r[6] == 1] print("cumcm2013a.csv 已写出：%d 样本，正常=%d 偏离=%d" % (N, len(g0), len(g1))) for j, nm in enumerate(["f1", "f2", "f3", "f4", "f5", "f6"]): m0 = sum(r[j] for r in g0) / len(g0); m1 = sum(r[j] for r in g1) / len(g1) print(" %s: 正常%.3f / 偏离%.3f" % (nm, m0, m1)) def gen_2012b(seed=2012, phi=37.0): """太阳能小屋：12 个月水平辐射 + 散射/直射拆分 + 组件参数。 用于板角(倾角)优化与容量配置收益优化。确定性(seeded)可复现。""" random.seed(seed) months = ["1月", "2月", "3月", "4月", "5月", "6月", "7月", "8月", "9月", "10月", "11月", "12月"] days = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] Hh, Hb, Hd = [], [], [] for mo in range(12): season = math.cos(2 * math.pi * (mo - 5.5) / 12) # 夏季(mo~5)最大 base = 4.6 + 1.9 * season + random.gauss(0, 0.12) base = max(base, 2.0) h = round(base, 2) fd = 0.80 + 0.10 * (1 - season) / 2 # 散射占比(全年偏高，使年最优倾角落于近纬度内点) hd = round(h * fd, 2) Hh.append(h); Hd.append(hd); Hb.append(round(h - hd, 2)) # 组件参数 eta = 0.18 # 光电转换效率 price = 0.75 # 上网电价 元/kWh cost_m2 = 1150.0 # 组件成本 元/m² rho = 0.20 # 地面反照率 path = os.path.join(OUT, "cumcm2012b.csv") with open(path, "w", encoding="utf-8-sig", newline="") as f: w = csv.writer(f) w.writerow(["月份", "水平辐射", "直射", "散射", "天数", "纬度", "效率", "电价", "成本_m2", "反照率"]) for mo in range(12): w.writerow([months[mo], Hh[mo], Hb[mo], Hd[mo], days[mo], phi, eta, price, cost_m2, rho]) print("cumcm2012b.csv 已写出：12 月辐射 + 组件参数") # 诊断：最优倾角(各月 Klein Rb 各向同性模型)与年辐射 def rb(mo, beta): d = math.radians(23.45 * math.sin(math.radians(360 * (284 + (mo + 1)) / 365))) ph = math.radians(phi); be = math.radians(beta) ws = math.acos(min(1.0, max(-1.0, -math.tan(ph) * math.tan(d)))) arg = -math.tan(ph - be) * math.tan(d) wsp = math.acos(min(1.0, max(-1.0, arg))) if -1 <= arg <= 1 else ws num = math.cos(ph - be) * math.cos(d) * math.sin(wsp) + wsp * math.sin(ph - be) * math.sin(d) den = math.cos(ph) * math.cos(d) * math.sin(ws) + ws * math.sin(ph) * math.sin(d) return num / den if den != 0 else 1.0 def E_per_m2(beta): tot = 0.0 for mo in range(12): Rb = rb(mo, beta) Hbt = Hb[mo] * Rb + Hd[mo] * (1 + math.cos(math.radians(beta))) / 2 + rho * Hh[mo] * (1 - math.cos(math.radians(beta))) / 2 tot += Hbt * days[mo] return tot best = max(range(0, 61), key=lambda b: E_per_m2(b)) print("诊断：最优倾角≈%d° 年辐射=%.1f kWh/m²" % (best, E_per_m2(best))) def gen_2012c(): """机器人避障（国赛 2012 C）：800×800 场景，机器人半径 R0=10，起点 O(0,0)， 目标 A/B/C，若干圆形障碍。路径须绕开障碍（与障碍边界距离 ≥ R0）。 数据为确定性的合成场景（硬编码布局，完全可复现），供可见性图 + 最短路算法直接练手。 固定布局保证重跑得到完全相同 CSV 与可复现的避障路径。""" R0 = 10.0 start = ("O", 0.0, 0.0) goals = [("A", 300.0, 300.0), ("B", 100.0, 700.0), ("C", 700.0, 640.0)] # 障碍：(编号, 圆心x, 圆心y, 半径) obstacles = [ ("1", 150.0, 200.0, 60.0), ("2", 370.0, 150.0, 55.0), ("3", 250.0, 470.0, 70.0), ("4", 520.0, 360.0, 65.0), ("5", 170.0, 640.0, 50.0), ("6", 560.0, 580.0, 70.0), ] path = os.path.join(OUT, "cumcm2012c.csv") with open(path, "w", encoding="utf-8-sig", newline="") as f: w = csv.writer(f) w.writerow(["type", "id", "x", "y", "r"]) w.writerow(["robot_radius", "", R0, "", ""]) w.writerow(["start", start[0], start[1], start[2], 0]) for g in goals: w.writerow(["goal", g[0], g[1], g[2], 0]) for o in obstacles: w.writerow(["obstacle", o[0], o[1], o[2], o[3]]) # 诊断：起点/目标是否落在任一膨胀圆（r+R0）内部 RI = [(o[1], o[2], o[3] + R0) for o in obstacles] ok = True for nm, x, y in [start] + goals: for (cx, cy, ri) in RI: if math.hypot(x - cx, y - cy) < ri - 1e-6: print("警告: %s 落在障碍内 (cx,cy,ri)=(%.1f,%.1f,%.1f)！" % (nm, cx, cy, ri)) ok = False print("cumcm2012c.csv 已写出：机器人半径=%g, 起点 %s, 目标 %d 个, 障碍 %d 个" % (R0, start[0], len(goals), len(obstacles))) print("诊断：起点/目标均位于膨胀障碍外 = %s" % ok) # 诊断：直线 O→A 是否被障碍 1 阻挡（用于验证绕行必要性） ox, oy = start[1], start[2]; ax, ay = goals[0][1], goals[0][2] dx, dy = ax - ox, ay - oy; L2 = dx * dx + dy * dy blocked = False for (cx, cy, ri) in RI: t = ((cx - ox) * dx + (cy - oy) * dy) / L2 t = max(0.0, min(1.0, t)) px, py = ox + t * dx, oy + t * dy if math.hypot(px - cx, py - cy) < ri - 1e-6: blocked = True print("诊断：直线 O→A 是否被障碍阻挡 = %s（应为 True，故必须绕行）" % blocked) def gen_2013b(seed=2013, N=19, K=10): """碎纸（纵切）复原：生成 N 条纵切碎纸竖条，每条含左/右边缘特征向量(K维)。 原始文档剖面沿竖条平滑演化，使相邻竖条边缘内容高度连续、非相邻不相关， 从而"右边缘 i 与左边缘 j 的相似度"在 j 恰为 i 的真实后继时达到峰值， 构成可解、可复现的碎纸拼接练习数据。打碎后打乱真实位置写出 CSV， 供相似度模型、匹配与路径搜索等算法直接练手。固定种子保证完全可复现。""" rng = random.Random(seed) # 原始文档剖面：沿竖条平滑演化的一维灰度序列（K 维） c = [[rng.gauss(0.0, 1.0) for _ in range(K)]] for j in range(1, N): prev = c[j - 1] # AR(1) 平滑：相邻竖条强相关，且与更早竖条弱相关 rho = 0.78 nxt = [rho * prev[k] + math.sqrt(1 - rho * rho) * rng.gauss(0.0, 1.0) for k in range(K)] c.append(nxt) # 竖条 pos(0-based)：左边缘 = c[pos]，右边缘 = c[pos+1]（最右条右边缘为边界 0） strips = [] for pos in range(N): L = c[pos][:] R = c[pos + 1][:] if pos < N - 1 else [0.0] * K strips.append((pos, L, R)) # 切碎：随机打乱真实位置，得到"混入"的碎片集合 order = list(range(N)); rng.shuffle(order) path = os.path.join(OUT, "cumcm2013b.csv") with open(path, "w", encoding="utf-8-sig", newline="") as f: w = csv.writer(f) w.writerow(["碎片编号", "真实位置"] + ["L%d" % k for k in range(K)] + ["R%d" % k for k in range(K)]) for new_id, real_pos in enumerate(order): pos, L, R = strips[real_pos] w.writerow([new_id + 1, pos + 1] + [round(v, 4) for v in L] + [round(v, 4) for v in R]) print("cumcm2013b.csv 已写出：%d 条碎纸竖条，每条 %d 维左右边缘特征" % (N, K)) # 诊断：无噪时相似度峰值应恰在真实后继处 def cos(a, b): na = math.sqrt(sum(x * x for x in a)); nb = math.sqrt(sum(x * x for x in b)) if na == 0 or nb == 0: return 0.0 return sum(a[k] * b[k] for k in range(K)) / (na * nb) ok = 0 for i in range(N): pos, L, R = strips[i] best_j, best_s = -1, -2.0 for j in range(N): if j == i: continue s = cos(R, strips[j][1]) # 与碎片 j 的左边缘相似度 if s > best_s: best_s, best_j = s, j if best_j == (pos + 1) % N and pos < N - 1: ok += 1 print("诊断：无噪下右边缘→正确左后继匹配成功率 = %d/%d" % (ok, N - 1)) def gen_2013c(seed=2013, N_layers=13, N_points=4, N_periods=10): """古塔变形：确定性合成数据（根因修复，语义贴合本题）。 每层 4 个角点观测，共 N_periods(=10) 期三维坐标；模拟缓慢倾斜（随期次增大斜率）、 整体沉降、层间压缩与轻微扭转，并叠加固定种子的测量噪声。写出 assets/problems/data/cumcm2013c.csv，表头 t,layer,point,x,y,z。""" rng = random.Random(seed) H0 = 8.0 # 基准层高（米） half = 3.0 # 半边长（米） sx0, sy0 = 0.0008, 0.0003 # 初始倾斜斜率（水平偏移/高度） dsx, dsy = 0.00012, 0.00005 # 倾斜增长速率（每期） dsettle = 0.08 # 整体沉降（每期，米） comp_rate = 0.002 # 层间压缩率（每期） twist = 0.0004 # 扭转速率（弧度/层/期） path = os.path.join(OUT, "cumcm2013c.csv") rows = [("t", "layer", "point", "x", "y", "z")] for t in range(N_periods): slope_x = sx0 + dsx * t slope_y = sy0 + dsy * t settle = dsettle * t for i in range(1, N_layers + 1): zc = (i - 1) * H0 * (1 - comp_rate * t) - settle # 层中心高度 cx = zc * slope_x # 层中心水平偏移 cy = zc * slope_y for p in range(1, N_points + 1): ang = math.pi / 4 + (p - 1) * math.pi / 2 + twist * (i - 1) * t dx = half * math.cos(ang) * math.sqrt(2) dy = half * math.sin(ang) * math.sqrt(2) x = cx + dx + rng.gauss(0, 0.04) y = cy + dy + rng.gauss(0, 0.04) z = zc + rng.gauss(0, 0.04) rows.append((t, i, p, round(x, 4), round(y, 4), round(z, 4))) with open(path, "w", newline="", encoding="utf-8") as f: csv.writer(f).writerows(rows) print("gen_2013c ->", path, "rows=", len(rows) - 1) print(" 顶层(t=0)水平偏移=%.4f m, 顶层(t=9)水平偏移=%.4f m" % ( 12 * H0 * sx0, 12 * H0 * (1 - comp_rate * 9) * (sx0 + dsx * 9))) print(" 整体沉降(t=9)=%.3f m, 顶层扭转(t=9)=%.4f rad" % ( dsettle * 9, twist * 12 * 9)) def land_trajectory(seed=2014, perturb=None, dt=0.5, over=None, wind=None): """月球软着陆最优下降轨迹（PD 推力控制，确定性）。 返回 (params, traj, land) ，traj 为采样点列表，land 为着陆信息。 over: 覆盖默认场景参数（v_app/Tmax/v_land/m0/isp/g/kh 等），用于扫描与灵敏度。 wind: 常值风扰加速度 (wx, wy, wz)，模拟飞行中持续扰动，产生可控着陆散布。""" g = 1.62 h0 = 15000.0 x0, y0 = 8000.0, 6000.0 vx0, vy0, vz0 = -110.0, -80.0, -210.0 m0 = 2400.0 mdry = 1200.0 ve = 300.0 * 9.81 # 排气速度 m/s Tmax = 35000.0 h_hover = 120.0 v_app = 30.0 v_land = 3.0 tau = 2.0 kh = 0.05 vmax_h = 200.0 if over: g = over.get("g", g); h0 = over.get("h0", h0) x0 = over.get("x0", x0); y0 = over.get("y0", y0) vx0 = over.get("vx0", vx0); vy0 = over.get("vy0", vy0); vz0 = over.get("vz0", vz0) m0 = over.get("m0", m0); mdry = over.get("mdry", mdry); ve = over.get("ve", ve) Tmax = over.get("Tmax", Tmax); h_hover = over.get("h_hover", h_hover) v_app = over.get("v_app", v_app); v_land = over.get("v_land", v_land) tau = over.get("tau", tau); kh = over.get("kh", kh); vmax_h = over.get("vmax_h", vmax_h) if perturb: x0 += perturb.get("x0", 0.0); y0 += perturb.get("y0", 0.0) h0 += perturb.get("h0", 0.0) vx0 += perturb.get("vx0", 0.0); vy0 += perturb.get("vy0", 0.0) vz0 += perturb.get("vz0", 0.0) params = dict(g=g, h0=h0, x0=x0, y0=y0, vx0=vx0, vy0=vy0, vz0=vz0, m0=m0, mdry=mdry, ve=ve, Tmax=Tmax, h_hover=h_hover, v_app=v_app, v_land=v_land, tau=tau, kh=kh, vmax_h=vmax_h) wx, wy, wz = wind if wind else (0.0, 0.0, 0.0) x, y, h = x0, y0, h0 vx, vy, vz = vx0, vy0, vz0 m = m0 traj = [] t = 0.0 land = None while h > 0 and m > mdry and t < 6000: speed = math.hypot(vx, vy, vz) if h > h_hover and speed > 80: phase = "brake" elif h > h_hover: phase = "approach" else: phase = "land" # 全阶段保持位置 PD（含着陆段），使风扰产生可控稳态偏移 if h > h_hover: vz_des = -v_app else: vz_des = -v_land vx_des = max(-vmax_h, min(vmax_h, -kh * x)) vy_des = max(-vmax_h, min(vmax_h, -kh * y)) ax_des = (vx_des - vx) / tau ay_des = (vy_des - vy) / tau az_des = (vz_des - vz) / tau # 推力加速度 = 期望净加速 - 重力 - 风扰 tax = ax_des - wx tay = ay_des - wy taz = az_des + g - wz Tmag = m * math.hypot(tax, tay, taz) if Tmag > Tmax: s = Tmax / Tmag tax, tay, taz = tax * s, tay * s, taz * s Tmag = Tmax vx += tax * dt vy += tay * dt vz += (taz - g) * dt x += vx * dt y += vy * dt h += vz * dt m -= (Tmag / ve) * dt t += dt traj.append((round(t, 2), phase, round(x, 2), round(h, 2), round(vx, 2), round(vz, 2), round(Tmag, 1), round(m, 2))) if h <= 0 and land is None: land = dict(t=round(t, 2), x=round(x, 2), y=round(y, 2), vx=round(vx, 2), vy=round(vy, 2), vz=round(vz, 2), m=round(m, 2), fuel=round(m0 - m, 2)) if land is None: land = dict(t=round(t, 2), x=round(x, 2), y=round(y, 2), vx=round(vx, 2), vy=round(vy, 2), vz=round(vz, 2), m=round(m, 2), fuel=round(m0 - m, 2)) return params, traj, land def gen_2014a(seed=2014): params, traj, land = land_trajectory(seed=seed) path = os.path.join(OUT, "cumcm2014a.csv") with open(path, "w", newline="", encoding="utf-8") as f: w = csv.writer(f) w.writerow(["# g=%.4f" % params["g"]]) w.writerow(["# h0=%.2f x0=%.2f y0=%.2f" % (params["h0"], params["x0"], params["y0"])]) w.writerow(["# vx0=%.2f vy0=%.2f vz0=%.2f" % (params["vx0"], params["vy0"], params["vz0"])]) w.writerow(["# m0=%.2f mdry=%.2f ve=%.2f Tmax=%.2f" % (params["m0"], params["mdry"], params["ve"], params["Tmax"])]) w.writerow(["# h_hover=%.2f v_app=%.2f v_land=%.2f" % (params["h_hover"], params["v_app"], params["v_land"])]) w.writerow(["t", "phase", "x", "h", "vx", "vz", "T", "m"]) for r in traj: w.writerow(r) print("gen_2014a ->", path, "samples=", len(traj)) print(" 着陆 t=%.1fs x=%.1f y=%.1f |vz|=%.2f m/s 燃料=%.1f kg" % ( land["t"], land["x"], land["y"], abs(land["vz"]), land["fuel"])) def gen_2014b(seed=2014, out_csv="cumcm2014b.csv"): """创意平板（cumcm2014b）：承重平板内部加筋布局优化。 物理模型（刚度控制）：面板被内部加筋划分为若干无支撑跨度为 a 的方格， 每格视为四边简支方板承受均布载荷 q。中心挠度 delta = alpha * q * a^4 / (E * tc^3) 令 delta <= delta_allow 得极限均布承载力 q_max = delta_allow * E * tc^3 / (alpha * a^4) [Pa] 重量 = rho * (底板体积 + 加筋体积)，其中加筋体积由 (nx, ny) 决定。 关键洞察：q_max 对跨度 a 是 4 次方衰减，对板厚 tc 仅 3 次方—— 因此“缩小加筋间距”比“加厚底板”更能提升承载力，且更省料。 """ E = 70e9 # 玻璃弹性模量 Pa rho = 2500.0 # 玻璃密度 kg/m^3 L = 1.2 # 面板长 m W = 0.8 # 面板宽 m delta = 5e-3 # 许用挠度 m hr = 0.018 # 加筋高 m tr = 0.004 # 加筋厚 m alpha = 0.00406 # 简支方板均布载荷中心挠度系数 rows = [] for tc_mm in range(3, 9): # 底板厚 3..8 mm tc = tc_mm / 1000.0 for nx in range(1, 21): # x 向加筋数 for ny in range(1, 15): # y 向加筋数 aW = W / (nx + 1) aL = L / (ny + 1) a = max(aW, aL) # 控制跨度取较大者 qmax = delta * E * tc ** 3 / (alpha * a ** 4) / 1000.0 # kPa wt = rho * (L * W * tc + (nx * L + ny * W) * hr * tr) # kg rows.append((tc_mm, nx, ny, round(a, 4), round(qmax, 3), round(wt, 3))) path = os.path.join(OUT, out_csv) with open(path, "w", newline="", encoding="utf-8") as f: w = csv.writer(f) w.writerow(["tc_mm", "nx", "ny", "a_m", "q_max_kPa", "weight_kg"]) for r in rows: w.writerow(r) print("gen_2014b ->", path, "samples=", len(rows)) return rows # --------------------------------------------------------------------------- # 纯 Python 线性规划（单纯形，Bland 规则防循环）与生猪养殖优化 # --------------------------------------------------------------------------- def _gauss_solve(A, b): """解线性方程组 A x = b（A 为 n×n 列表，纯 Python）。""" n = len(A) M = [row[:] + [b[i]] for i, row in enumerate(A)] for c in range(n): p = max(range(c, n), key=lambda r: abs(M[r][c])) if abs(M[p][c]) < 1e-12: continue M[c], M[p] = M[p], M[c] piv = M[c][c] M[c] = [v / piv for v in M[c]] for r in range(n): if r != c and M[r][c] != 0: f = M[r][c] M[r] = [M[r][k] - f * M[c][k] for k in range(n + 1)] return [M[i][n] for i in range(n)] def _harmonic_fit(ts, K): """谐波回归 ts[t] ≈ a0 + a1*t + b*cos(2π t/K) + c*sin(2π t/K)。""" n = len(ts) cols = [] for t in range(n): cols.append([1.0, float(t), math.cos(2 * math.pi * t / K), math.sin(2 * math.pi * t / K)]) XtX = [[sum(cols[i][a] * cols[i][b] for i in range(n)) for b in range(4)] for a in range(4)] Xty = [sum(cols[i][a] * ts[i] for i in range(n)) for a in range(4)] beta = _gauss_solve(XtX, Xty) fit = [beta[0] + beta[1] * t + beta[2] * math.cos(2 * math.pi * t / K) + beta[3] * math.sin(2 * math.pi * t / K) for t in range(n)] resid = [ts[i] - fit[i] for i in range(n)] s = math.sqrt(sum(r * r for r in resid) / max(1, n - 4)) return beta, fit, s def _period_search(ts, lo=24, hi=48): best_K, best_r2 = lo, -1e9 n = len(ts) mean = sum(ts) / n sst = sum((v - mean) ** 2 for v in ts) for K in range(lo, hi + 1): _, fit, _ = _harmonic_fit(ts, K) ssr = sum((ts[i] - fit[i]) ** 2 for i in range(n)) r2 = 1 - ssr / sst if r2 > best_r2: best_r2, best_K = r2, K return best_K, best_r2 def _lp_simplex(c, A, b): """最大化 cᵀx，约束 A x ≤ b, x ≥ 0。返回 (x, obj)。""" m, n = len(A), len(c) M = [[0.0] * (n + m + 1) for _ in range(m + 1)] for i in range(m): for j in range(n): M[i][j] = A[i][j] M[i][n + i] = 1.0 M[i][-1] = b[i] for j in range(n): M[m][j] = -c[j] basis = [n + i for i in range(m)] eps = 1e-9 while True: ent = -1 for j in range(n + m): if M[m][j] < -eps: ent = j break if ent == -1: break pivot_row = -1 best_ratio = 1e18 for i in range(m): if M[i][ent] > eps: ratio = M[i][-1] / M[i][ent] if ratio < best_ratio - eps: best_ratio = ratio pivot_row = i if pivot_row == -1: break piv = M[pivot_row][ent] M[pivot_row] = [v / piv for v in M[pivot_row]] for i in range(m + 1): if i != pivot_row and M[i][ent] != 0: f = M[i][ent] M[i] = [M[i][k] - f * M[pivot_row][k] for k in range(n + m + 1)] basis[pivot_row] = ent x = [0.0] * n for i in range(m): if basis[i] < n: x[basis[i]] = M[i][-1] obj = sum(c[j] * x[j] for j in range(n)) return x, obj def gen_2014c(seed=2014, out_csv="cumcm2014c.csv"): """生猪养殖（cumcm2014c）：价格周期预测 + 补栏/出栏生产计划优化。""" random.seed(seed) HIST = 60 base, amp, K_true, trend = 15.0, 5.0, 36, 0.03 price, piglet, feed = [], [], [] for t in range(HIST): noise = (random.random() - 0.5) * 0.8 p = base + amp * math.sin(2 * math.pi * t / K_true) + trend * t + noise price.append(round(p, 3)) piglet.append(round(540.0 + 12.0 * (p - 15.0) + (random.random() - 0.5) * 20.0, 1)) feed.append(round(185.0 + 10.0 * math.sin(2 * math.pi * t / 24) + (random.random() - 0.5) * 6.0, 1)) K, r2 = _period_search(price, 24, 48) beta, fit, sres = _harmonic_fit(price, K) FUT = 36 fc_price, fc_lo, fc_hi = [], [], [] for t in range(HIST, HIST + FUT): v = beta[0] + beta[1] * t + beta[2] * math.cos(2 * math.pi * t / K) + beta[3] * math.sin(2 * math.pi * t / K) fc_price.append(round(v, 3)) fc_lo.append(round(v - 1.96 * sres, 3)) fc_hi.append(round(v + 1.96 * sres, 3)) beta_pi, _, _ = _harmonic_fit(piglet, K) beta_fe, _, _ = _harmonic_fit(feed, K) fc_piglet = [round(beta_pi[0] + beta_pi[1] * t + beta_pi[2] * math.cos(2 * math.pi * t / K) + beta_pi[3] * math.sin(2 * math.pi * t / K), 1) for t in range(HIST, HIST + FUT)] fc_feed = [round(beta_fe[0] + beta_fe[1] * t + beta_fe[2] * math.cos(2 * math.pi * t / K) + beta_fe[3] * math.sin(2 * math.pi * t / K), 1) for t in range(HIST, HIST + FUT)] def mape(true, pred): return 100.0 * sum(abs(true[i] - pred[i]) / true[i] for i in range(len(true))) / len(true) def harmonic_pred(train, Kk, n_fut): bt, ft, st = _harmonic_fit(train, Kk) nt = len(train) return [bt[0] + bt[1] * (nt + k) + bt[2] * math.cos(2 * math.pi * (nt + k) / Kk) + bt[3] * math.sin(2 * math.pi * (nt + k) / Kk) for k in range(n_fut)] def holt_pred(train, n_fut, alpha=0.4, beta=0.1): l, b = train[0], train[1] - train[0] for i in range(1, len(train)): lp = l + b err = train[i] - lp l = lp + alpha * err b = b + beta * err return [l + (k + 1) * b for k in range(n_fut)] def ma_pred(train, n_fut, w=3): return [sum(train[-w:]) / w for _ in range(n_fut)] def es_pred(train, n_fut, alpha=0.3): s = train[0] for v in train[1:]: s = alpha * v + (1 - alpha) * s return [s] * n_fut trn, tes = price[:48], price[48:] m_har = mape(tes, harmonic_pred(trn, K, 12)) m_holt = mape(tes, holt_pred(trn, 12)) m_ma = mape(tes, ma_pred(trn, 12)) m_es = mape(tes, es_pred(trn, 12)) W = 110.0 T = 6 CAP = 120 N = 30 R0 = HIST + 1 rmonths = list(range(R0, R0 + N)) smonths = [r + T for r in rmonths] coef = [] for idx, r in enumerate(rmonths): sidx = smonths[idx] - HIST pv = fc_price[sidx] if 0 <= sidx < FUT else fc_price[-1] pig = fc_piglet[r - HIST] if (r - HIST) < FUT else fc_piglet[-1] feed_sum = sum(fc_feed[r - HIST + k] for k in range(T) if (r - HIST + k) < FUT) coef.append(W * pv - pig - feed_sum) months = list(range(R0, R0 + N + T)) A, b = [], [] for m in months: row = [1.0 if (r <= m < r + T) else 0.0 for r in rmonths] if any(v > 0 for v in row): A.append(row) b.append(CAP) x_relax, obj_relax = _lp_simplex(coef, A, b) x_opt = [int(round(v)) for v in x_relax] # 可行性修正（区间矩阵全单模，LP 最优本应已整数可行；此处仅对舍入误差做安全兜底） for _ in range(5): ok = True for m in months: active = [j for j in range(N) if rmonths[j] <= m < rmonths[j] + T and x_opt[j] > 0] s = sum(x_opt[j] for j in active) if s > CAP: ok = False active.sort(key=lambda j: coef[j]) x_opt[active[0]] -= (s - CAP) if ok: break profit = sum(x_opt[j] * coef[j] for j in range(N)) inv = {m: sum(x_opt[j] for j in range(N) if rmonths[j] <= m < rmonths[j] + T) for m in months} revenue = sum(x_opt[j] * W * fc_price[min(smonths[j] - HIST, FUT - 1)] for j in range(N) if x_opt[j] > 0) cost_buy = sum(x_opt[j] * fc_piglet[rmonths[j] - HIST] for j in range(N) if x_opt[j] > 0 and (rmonths[j] - HIST) < FUT) cost_feed = sum(x_opt[j] * sum(fc_feed[rmonths[j] - HIST + k] for k in range(T) if (rmonths[j] - HIST + k) < FUT) for j in range(N) if x_opt[j] > 0) avg_p = sum(price) / HIST # 公平对照：同样受产能约束的“持续满产”计划 —— 每月稳定补栏 CAP/T 头（产能 ~20/月）， # 不择时，出栏价取各月预测价。总出栏同样约 600 头，仅与最优解在“择时”上不同。 steady = CAP // T profit_steady = sum(steady * coef[j] for j in range(N)) def profit_with_scale(scale): c2 = [coef[j] * scale for j in range(N)] xr, _ = _lp_simplex(c2, A, b) return sum(int(round(v)) * coef[j] * scale for j, v in enumerate(xr)) sens = [(round(s, 2), profit_with_scale(s)) for s in (0.9, 1.0, 1.1)] path = os.path.join(OUT, out_csv) with open(path, "w", newline="", encoding="utf-8") as f: w = csv.writer(f) w.writerow(["month", "type", "price", "piglet_price", "feed_month"]) for t in range(HIST): w.writerow([t + 1, "hist", price[t], piglet[t], feed[t]]) for k in range(FUT): w.writerow([HIST + 1 + k, "fc", fc_price[k], fc_piglet[k], fc_feed[k]]) w.writerow([]) w.writerow(["restock_month", "x_heads", "slaughter_month", "margin_per_head"]) for j in range(N): if x_opt[j] > 0: w.writerow([rmonths[j], x_opt[j], smonths[j], round(coef[j], 1)]) print("gen_2014c ->", path) print(" 周期 K=%d R2=%.4f 残差std=%.3f" % (K, r2, sres)) print(" MAPE: 谐波=%.2f%% Holt=%.2f%% MA=%.2f%% ES=%.2f%%" % (m_har, m_holt, m_ma, m_es)) print(" 最优收益=%.0f 对照(持续满产)=%.0f 总补栏=%d" % (profit, profit_steady, sum(x_opt))) print(" 敏感性(±10%预测价):", sens) return dict(K=K, r2=r2, sres=sres, mape=(m_har, m_holt, m_ma, m_es), profit=profit, profit_steady=profit_steady, x_opt=x_opt, rmonths=rmonths, smonths=smonths, coef=coef, sens=sens, FUT=FUT, HIST=HIST, price=price, piglet=piglet, feed=feed, fc_price=fc_price, fc_lo=fc_lo, fc_hi=fc_hi, fc_piglet=fc_piglet, fc_feed=fc_feed, inv=inv, revenue=revenue, cost_buy=cost_buy, cost_feed=cost_feed, CAP=CAP, T=T, W=W, months=months) # ---------------------------------------------------------------------------- # 2015A 太阳影子定位：由垂直杆影子轨迹反演经纬度与日期（确定性天文几何场景） # ---------------------------------------------------------------------------- def declination(nday): """太阳赤纬（度），Cooper 近似。""" return 23.45 * math.sin(math.radians(360.0 * (284.0 + nday) / 365.0)) def shadow_tip(clock, phi, lam, delta, H): """给定相机时钟(UTC+8, h)、纬度/经度(度)、赤纬(度)、杆高(m)， 返回影子尖端地面坐标 (x, y)（米，北=+y 东=+x）。""" LST = clock - 8.0 + lam / 15.0 # 本地太阳时（忽略时差方程） w = math.radians(15.0 * (LST - 12.0)) # 时角 ph, de = math.radians(phi), math.radians(delta) sina = math.sin(ph) * math.sin(de) + math.cos(ph) * math.cos(de) * math.cos(w) sina = max(-1.0, min(1.0, sina)) a = math.asin(sina) # 太阳高度角 sE = math.cos(de) * math.sin(w) sN = math.cos(de) * math.cos(w) * math.sin(ph) - math.sin(de) * math.cos(ph) # 影长 r = H·cotα = H·cosα/sinα；水平方向单位向量为 (sE,sN)/cosα， # 故尖端 = -r·(sE,sN)/cosα = -(H/sinα)·(sE,sN) L = H / math.sin(a) return -L * sE, -L * sN def _lm_invert(obs, H, p0, iters=400, history=False): """纯 Python Levenberg-Marquardt 反演 (phi, lam, delta)。obs: list of (clock,x,y)。""" hist = [] def resid(p): r = [] for c, xo, yo in obs: xp, yp = shadow_tip(c, p[0], p[1], p[2], H) r.append(xp - xo); r.append(yp - yo) return r def jac(p): h = 1e-4; J = []; r0 = resid(p) for k in range(3): pk = list(p); pk[k] += h; rk = resid(pk) J.append([(rk[i] - r0[i]) / h for i in range(len(r0))]) return [[J[k][i] for k in range(3)] for i in range(len(r0))] p = list(p0) for _ in range(iters): r = resid(p); sse = sum(v * v for v in r) J = jac(p); n = len(r); m = 3 JtJ = [[sum(J[i][a] * J[i][b] for i in range(n)) for b in range(m)] for a in range(m)] Jtr = [sum(J[i][a] * r[i] for i in range(n)) for a in range(m)] improved = False; lam = 1.0 for _ in range(12): A = [[JtJ[a][b] + (lam if a == b else 0) for b in range(m)] for a in range(m)] M = [row[:] for row in A]; v = [-Jtr[a] for a in range(m)] for col in range(m): piv = max(range(col, m), key=lambda rr: abs(M[rr][col])) M[col], M[piv] = M[piv], M[col]; v[col], v[piv] = v[piv], v[col] pv = M[col][col] for rr in range(m): if rr != col: f = M[rr][col] / pv M[rr] = [M[rr][k] - f * M[col][k] for k in range(m)]; v[rr] -= f * v[col] dp = [v[a] / M[a][a] for a in range(m)] pn = [p[a] + dp[a] for a in range(m)] rn = resid(pn); ssn = sum(vv * vv for vv in rn) if ssn < sse: p = pn; lam = max(lam / 2.0, 1e-9); improved = True; break else: lam *= 2.0 if not improved: break if history: hist.append((sum(v * v for v in resid(p)), list(p))) if history: return p, sum(v * v for v in resid(p)), hist return p, sum(v * v for v in resid(p)) def invert_location(obs, H, iters=400): """多重启动 LM 反演 (phi, lam, delta)，返回 (params, sse)。 网格在纬度/经度/赤纬上均加密，确保夏至附近这种近退化情形下仍能锁定全局极小。""" best = None for phi0 in (20.0, 30.0, 39.0, 45.0, 55.0): for lam0 in (100.0, 108.0, 116.0, 120.0, 125.0): for del0 in (15.0, 20.0, 22.0, 23.0, 23.45, 23.6): pf, sse = _lm_invert(obs, H, [phi0, lam0, del0], iters) if best is None or sse < best[1]: best = (pf, sse) return best def gen_2015a(seed=2015, out_csv="cumcm2015a.csv"): random.seed(seed) H = 3.0 # 杆高 m phi_true, lam_true, n_true = 39.9, 116.4, 172 # 北京，夏至 delta_true = declination(n_true) clocks = [9.0 + 0.5 * i for i in range(13)] # 9:00-15:00 每 0.5h noise = 0.02 # 影子坐标观测噪声 米 obs = [] for c in clocks: x, y = shadow_tip(c, phi_true, lam_true, delta_true, H) obs.append((round(c, 2), round(x + random.uniform(-noise, noise), 4), round(y + random.uniform(-noise, noise), 4))) # 写 CSV：相机时钟、影子尖端 x/y（米） path = os.path.join(OUT, out_csv) with open(path, "w", encoding="utf-8") as f: f.write("clock_h,x_m,y_m\n") for c, x, y in obs: f.write("%.2f,%.4f,%.4f\n" % (c, x, y)) rows = obs # 反演 (phi_f, lam_f, delta_f), sse = invert_location(obs, H) best_guess = [phi_f, lam_f, delta_f] # 由赤纬反推日序（全局搜索取与估计赤纬最吻合的日序，解决半年度歧义） best_n = min(range(1, 366), key=lambda nd: abs(declination(nd) - delta_f)) # 蒙特卡洛不确定性：对观测加 σ 扰动重反演 N 次 N = 200; mc_sig = 0.02 mc_phi, mc_lam, mc_del = [], [], [] for _ in range(N): obs2 = [(c, x + random.gauss(0, mc_sig), y + random.gauss(0, mc_sig)) for (c, x, y) in obs] pf, _ = _lm_invert(obs2, H, best_guess) mc_phi.append(pf[0]); mc_lam.append(pf[1]); mc_del.append(pf[2]) import statistics as _st mc = dict(phi_std=_st.pstdev(mc_phi), lam_std=_st.pstdev(mc_lam), delta_std=_st.pstdev(mc_del), phi_mean=_st.mean(mc_phi), lam_mean=_st.mean(mc_lam), phi_lo=phi_f - 1.96 * _st.pstdev(mc_phi), phi_hi=phi_f + 1.96 * _st.pstdev(mc_phi), lam_lo=lam_f - 1.96 * _st.pstdev(mc_lam), lam_hi=lam_f + 1.96 * _st.pstdev(mc_lam)) # 影子长度与最低点（正午）时刻：离散最小值 + 抛物插值细化 r = [math.hypot(x, y) for (_, x, y) in obs] rmin_idx = min(range(len(obs)), key=lambda i: math.hypot(obs[i][1], obs[i][2])) noon_clock = obs[rmin_idx][0] if 0 < rmin_idx < len(obs) - 1: rp = math.hypot(obs[rmin_idx - 1][1], obs[rmin_idx - 1][2]) rc = math.hypot(obs[rmin_idx][1], obs[rmin_idx][2]) rn = math.hypot(obs[rmin_idx + 1][1], obs[rmin_idx + 1][2]) denom = (rp - 2 * rc + rn) if abs(denom) > 1e-9: noon_clock = noon_clock + 0.5 * (rp - rn) / denom * 0.5 # 由细化正午时刻反演经度：λ = 15·(20 - t_noon) [LST=clock-8+λ/15=12] lam_noon = 15.0 * (20.0 - noon_clock) # 噪声灵敏度：确定种子下多 σ 的重反演误差（warm start） sens = [] for sig in (0.005, 0.01, 0.02, 0.04, 0.08): eph = []; ela = [] for _ in range(60): o2 = [(c, x + random.gauss(0, sig), y + random.gauss(0, sig)) for (c, x, y) in obs] pp, _ = _lm_invert(o2, H, [phi_f, lam_f, delta_f], iters=120) eph.append(abs(pp[0] - phi_true)); ela.append(abs(pp[1] - lam_true)) sens.append((sig, _st.mean(eph), _st.mean(ela))) res = dict(H=H, phi_true=phi_true, lam_true=lam_true, n_true=n_true, delta_true=delta_true, clocks=clocks, obs=obs, phi_f=phi_f, lam_f=lam_f, delta_f=delta_f, sse=sse, n_est=best_n, mc=mc, r=r, noon_clock=noon_clock, lam_noon=lam_noon, rmin=min(r), rmax=max(r), sens=sens, mc_phi=mc_phi, mc_lam=mc_lam, mc_del=mc_del) print("gen_2015a ->", path, "samples=", len(rows)) print(" 真值 phi=%.2f lam=%.2f delta=%.3f" % (phi_true, lam_true, delta_true)) print(" 反演 phi=%.3f lam=%.3f delta=%.3f SSE=%.5f" % (phi_f, lam_f, delta_f, sse)) print(" 估计日序=%d 噪声σ=%.2f时 纬度95%%CI=[%.2f,%.2f] 经度95%%CI=[%.2f,%.2f]" % (best_n, mc_sig, mc["phi_lo"], mc["phi_hi"], mc["lam_lo"], mc["lam_hi"])) return res def gen_2015b(seed=2015, out_csv="cumcm2015b.csv"): """电商物流分拣中心选址 + 配送路径优化（LRP 简化）。 确定性数据（固定 seed）+ 可复现求解： - 区域平面 [0,100]x[0,100] km；5 个候选分拣中心（坐标/建设费/容量） - 30 个需求点（坐标/日需求量） - 成本 = 建设费 + 干线运输费(元/吨·km) + 车辆行驶费(元/km) - 求解：枚举选址子集(2^5=32) + 贪婪容量分配 + 极角排序 VRP 近似 返回 dict 含所有权威数字，供 fig_2015b / 三篇范文 / 附录共用。 """ random.seed(seed) X0, Y0, X1, Y1 = 0.0, 0.0, 100.0, 100.0 P = 5 cand = [] for j in range(P): x = round(10 + 80 * random.random(), 2) y = round(10 + 80 * random.random(), 2) f = round(120 + 60 * random.random(), 1) # 建设费 万元 cap = round(160 + 100 * random.random(), 1) # 日处理容量 吨 cand.append(dict(idx=j, x=x, y=y, f=f, cap=cap)) Dn = 30 dem = [] for i in range(Dn): x = round(X0 + (X1 - X0) * random.random(), 2) y = round(Y0 + (Y1 - Y0) * random.random(), 2) q = round(4 + 9 * random.random(), 2) # 日需求量 吨 dem.append((x, y, q)) total_q = sum(q for _, _, q in dem) c_line = 0.6 # 干线运输 元/(吨·km) c_veh = 2.5 # 车辆行驶 元/km Q = 35.0 # 单车容量 吨 AY = 365.0 # 年化因子：将日运输/车辆成本折算为年度（万元/年） def dist(a, b): return math.hypot(a[0] - b[0], a[1] - b[1]) # ---------- 重心法（篇1 初步选址） ---------- cx = sum(x * q for x, _, q in dem) / total_q cy = sum(y * q for _, y, q in dem) / total_q citer = [(round(cx, 2), round(cy, 2))] for _ in range(12): w = [q / (dist((x, y), (cx, cy)) + 1e-9) for x, y, q in dem] cx = sum(x * qk for (x, _, q), qk in zip(dem, w)) / sum(w) cy = sum(y * qk for (_, y, q), qk in zip(dem, w)) / sum(w) citer.append((round(cx, 2), round(cy, 2))) centroid = (round(cx, 2), round(cy, 2)) nearest = min(range(P), key=lambda j: dist((cx, cy), (cand[j]["x"], cand[j]["y"]))) # ---------- 覆盖模型（篇1） ---------- def coverage(radius, CD): cnt = 0 for (x, y, q) in CD: if any(dist((x, y), (c["x"], c["y"])) <= radius for c in cand): cnt += 1 return cnt / len(CD) cover = [(R, round(coverage(R, dem), 3)) for R in (12, 16, 20, 24, 28, 32, 36, 40)] # ---------- 分配 + 极角路径 启发式 ---------- def assign_cost(open_idx, caps, CD): rem = dict(caps) a = {} for i, (x, y, q) in enumerate(CD): best, bc = None, 1e18 for j in open_idx: if rem[j] >= q - 1e-9: cst = c_line * dist((x, y), (cand[j]["x"], cand[j]["y"])) * q if cst < bc: bc, best = cst, j if best is None: return None a[i] = best rem[best] -= q return a def route_for(wh_idx, pidx, CD): wx, wy = cand[wh_idx]["x"], cand[wh_idx]["y"] items = [(CD[i][2], i) for i in pidx] if not items: return [], [], 0.0, 0 ang = sorted(items, key=lambda it: math.atan2(CD[it[1]][1] - wy, CD[it[1]][0] - wx)) vehs, cur, cq = [], [], 0.0 for q, i in ang: if cq + q <= Q + 1e-9: cur.append(i); cq += q else: vehs.append(cur); cur = [i]; cq = q if cur: vehs.append(cur) total = 0.0 for v in vehs: d, prev = 0.0, (wx, wy) for i in v: d += dist(prev, (CD[i][0], CD[i][1])); prev = (CD[i][0], CD[i][1]) d += dist(prev, (wx, wy)) total += d return vehs, [v for v in vehs], round(total, 2), len(vehs) def solve_subset(open_idx, CD): caps = {j: cand[j]["cap"] for j in open_idx} a = assign_cost(open_idx, caps, CD) if a is None: return None build = sum(cand[j]["f"] for j in open_idx) trans = sum(c_line * dist((CD[i][0], CD[i][1]), (cand[a[i]]["x"], cand[a[i]]["y"])) * CD[i][2] for i in range(len(CD))) veh_total, nveh, loads, routes = 0.0, 0, {j: 0.0 for j in open_idx}, {} for j in open_idx: pidx = [i for i in range(len(CD)) if a[i] == j] loads[j] = round(sum(CD[i][2] for i in pidx), 2) vehs, _, tl, nv = route_for(j, pidx, CD) routes[j] = vehs veh_total += tl nveh += nv return dict(subset=list(open_idx), assign=[a[i] for i in range(len(CD))], build=build, trans=trans, veh_len=veh_total, nveh=nveh, loads=loads, routes=routes) # ---------- 枚举选址子集求最优 ---------- best = None all_subsets = [] for mask in range(1, 1 << P): open_idx = [j for j in range(P) if mask & (1 << j)] s = solve_subset(open_idx, dem) if s is None: continue trans_w = s["trans"] * AY / 10000.0 veh_w = s["veh_len"] * c_veh * AY / 10000.0 total = s["build"] + trans_w + veh_w all_subsets.append((tuple(open_idx), round(total, 3))) if best is None or total < best["total"]: best = dict(subset=s["subset"], assign=s["assign"], build=round(s["build"], 3), trans=round(trans_w, 3), veh=round(veh_w, 3), total=round(total, 3), loads=s["loads"], routes=s["routes"], nveh=s["nveh"], veh_len=round(s["veh_len"], 2)) # ---------- 基线对比 ---------- full = solve_subset(list(range(P)), dem) full_total = full["build"] + full["trans"] * AY / 10000.0 + full["veh_len"] * c_veh * AY / 10000.0 # 贪婪启发式（最便宜优先，直至容量满足需求）：典型的"先建再算"朴素方案 order = sorted(range(P), key=lambda j: cand[j]["f"]) cur, cap_sum = [], 0.0 for j in order: cur.append(j); cap_sum += cand[j]["cap"] if cap_sum >= total_q - 1e-9: break gs = solve_subset(cur, dem) greedy = dict(subset=list(cur), total=round(gs["build"] + gs["trans"] * AY / 10000.0 + gs["veh_len"] * c_veh * AY / 10000.0, 3)) # ---------- 灵敏度分析（固定最优选址） ---------- def cost_with(c_line_v=None, q_mul=1.0, f_mul=1.0, CD=None): c_line_v = c_line if c_line_v is None else c_line_v if CD is None: CD = [(x, y, q * q_mul) for x, y, q in dem] s = solve_subset(best["subset"], CD) build_v = sum(cand[j]["f"] * f_mul for j in best["subset"]) trans_v = sum(c_line_v * dist((CD[i][0], CD[i][1]), (cand[s["assign"][i]]["x"], cand[s["assign"][i]]["y"])) * CD[i][2] for i in range(len(CD))) * AY / 10000.0 veh_v = s["veh_len"] * c_veh * AY / 10000.0 return round(build_v + trans_v + veh_v, 3) sens_cline = [(cl, cost_with(c_line_v=cl)) for cl in (0.3, 0.45, 0.6, 0.9, 1.2)] sens_demand = [(dm, cost_with(q_mul=dm)) for dm in (0.7, 0.85, 1.0, 1.15, 1.3)] sens_build = [(fm, cost_with(f_mul=fm)) for fm in (0.6, 0.8, 1.0, 1.2, 1.5)] # ---------- 蒙特卡洛不确定性（需求扰动，固定最优选址重算） ---------- Nmc = 200 mc_samples = [] for _ in range(Nmc): dem2 = [(x, y, max(0.5, q * (1 + 0.15 * random.gauss(0, 1)))) for x, y, q in dem] s = solve_subset(best["subset"], dem2) build_v = sum(cand[j]["f"] for j in best["subset"]) trans_v = sum(c_line * dist((dem2[i][0], dem2[i][1]), (cand[s["assign"][i]]["x"], cand[s["assign"][i]]["y"])) * dem2[i][2] for i in range(len(dem2))) * AY / 10000.0 veh_v = s["veh_len"] * c_veh * AY / 10000.0 mc_samples.append(round(build_v + trans_v + veh_v, 3)) import statistics as _st mc_mean = _st.mean(mc_samples) mc_std = _st.pstdev(mc_samples) mc = dict(mean=round(mc_mean, 3), std=round(mc_std, 3), lo=round(mc_mean - 1.96 * mc_std, 3), hi=round(mc_mean + 1.96 * mc_std, 3), samples=mc_samples) # ---------- 写 CSV ---------- path = os.path.join(OUT, out_csv) with open(path, "w", encoding="utf-8") as f: f.write("type,id,x_km,y_km,q_ton_or_cap,fee_wan\n") for c in cand: f.write("candidate,%d,%.2f,%.2f,%.1f,%.1f\n" % (c["idx"], c["x"], c["y"], c["cap"], c["f"])) for i, (x, y, q) in enumerate(dem): f.write("demand,%d,%.2f,%.2f,%.2f,0\n" % (i, x, y, q)) res = dict(P=P, Dn=Dn, cand=cand, dem=dem, total_q=round(total_q, 2), c_line=c_line, c_veh=c_veh, Q=Q, AY=AY, centroid=centroid, citer=citer, nearest=nearest, cover=cover, best=best, all_subsets=all_subsets, full_total=round(full_total, 3), greedy=greedy, sens_cline=sens_cline, sens_demand=sens_demand, sens_build=sens_build, mc=mc) print("gen_2015b ->", path, "candidates=", P, "demands=", Dn) print(" 最优选址子集=%s 总成本=%.3f 万元 (建设%.3f + 运输%.3f + 车辆%.3f)" % (best["subset"], best["total"], best["build"], best["trans"], best["veh"])) print(" 基线: 全选=%.3f 万, 贪婪增量=%s=%.3f 万" % (full_total, greedy["subset"], greedy["total"])) print(" MC 总成本均值=%.3f 万元 95%%CI=[%.3f,%.3f]" % (mc["mean"], mc["lo"], mc["hi"])) return res def _mean(xs): return sum(xs) / len(xs) def _std(xs): m = _mean(xs) return math.sqrt(sum((x - m) ** 2 for x in xs) / len(xs)) or 1.0 def _median(xs): s = sorted(xs); n = len(s) return s[n // 2] if n % 2 else (s[n // 2 - 1] + s[n // 2]) / 2.0 def _mad(xs, med): """中位数绝对偏差（未乘 1.4826 常数，调用处按需缩放）。""" return _median([abs(x - med) for x in xs]) or 1.0 def _zscale(M): """对特征矩阵按列 z-score 标准化，返回 (Z, means, stds)。""" n, k = len(M), len(M[0]) means = [_mean([M[i][j] for i in range(n)]) for j in range(k)] stds = [_std([M[i][j] for i in range(n)]) for j in range(k)] Z = [[(M[i][j] - means[j]) / (stds[j] or 1) for j in range(k)] for i in range(n)] return Z, means, stds def _dist2(a, b): return sum((x - y) ** 2 for x, y in zip(a, b)) def _kmeans(Z, K, seed=2015, restarts=25): """确定性 k-means++（固定种子）+ 多重启，返回 (labels, centers, inertia)。""" random.seed(seed) n, k = len(Z), len(Z[0]) best = None for _ in range(restarts): # k-means++ 初始化 cents = [list(Z[random.randrange(n)])] while len(cents) < K: d = [_min(_dist2(Z[i], c) for c in cents) for i in range(n)] tot = sum(d) r = random.random() * tot acc = 0.0 for i in range(n): acc += d[i] if acc >= r: cents.append(list(Z[i])); break else: cents.append(list(Z[random.randrange(n)])) # Lloyd 迭代 for _ in range(60): lab = [_argmin(_dist2(Z[i], c) for c in cents) for i in range(n)] newc = [] for j in range(K): grp = [Z[i] for i in range(n) if lab[i] == j] if grp: newc.append([_mean(col) for col in zip(*grp)]) else: newc.append(list(Z[random.randrange(n)])) if all(_dist2(newc[j], cents[j]) < 1e-9 for j in range(K)): cents = newc; break cents = newc lab = [_argmin(_dist2(Z[i], c) for c in cents) for i in range(n)] inertia = sum(_dist2(Z[i], cents[lab[i]]) for i in range(n)) if best is None or inertia < best[2]: best = (lab, cents, inertia) return best def _argmin(gen): best_i, best_v = 0, None for i, v in enumerate(gen): if best_v is None or v < best_v: best_v, best_i = v, i return best_i def _min(gen): return min(gen) def _silhouette(Z, lab, K): """轮廓系数（基于标准化特征欧氏距离）。""" n = len(Z) a = [0.0] * n; b = [0.0] * n for i in range(n): # a: 同簇平均距离 same = [j for j in range(n) if j != i and lab[j] == lab[i]] ai = _mean([math.sqrt(_dist2(Z[i], Z[j])) for j in same]) if same else 0.0 a[i] = ai # b: 最近异簇平均距离 best_b, best_k = 1e18, None for k in range(K): if k == lab[i]: continue grp = [j for j in range(n) if lab[j] == k] if not grp: continue dk = _mean([math.sqrt(_dist2(Z[i], Z[j])) for j in grp]) if dk < best_b: best_b, best_k = dk, k b[i] = best_b if best_k is not None else 0.0 s = [_safe_s(a[i], b[i]) for i in range(n)] return _mean(s) def _safe_s(a, b): m = max(a, b) return 0.0 if m == 0 else (b - a) / m def _cov_inv(Z): """标准化特征协方差矩阵的逆（k×k，对角占优近似，确保可逆）。""" n, k = len(Z), len(Z[0]) cov = [[0.0] * k for _ in range(k)] for i in range(n): for p in range(k): for q in range(k): cov[p][q] += Z[i][p] * Z[i][q] for p in range(k): for q in range(k): cov[p][q] /= n for p in range(k): cov[p][p] += 0.05 # 正则化保证可逆 return _inv2(cov) def _inv2(A): n = len(A) M = [row[:] + [1.0 if i == j else 0.0 for j in range(n)] for i, row in enumerate(A)] for col in range(n): piv = max(range(col, n), key=lambda r: abs(M[r][col])) M[col], M[piv] = M[piv], M[col] d = M[col][col] M[col] = [v / d for v in M[col]] for r in range(n): if r != col: f = M[r][col] M[r] = [M[r][c] - f * M[col][c] for c in range(2 * n)] return [row[n:] for row in M] def _maha2(z, inv): """马氏距离平方 z^T Σ^{-1} z。""" k = len(z) v = [sum(inv[p][q] * z[q] for q in range(k)) for p in range(k)] return sum(z[p] * v[p] for p in range(k)) def gen_2015c(seed=2015, N=12, L=120, out_csv="cumcm2015c.csv"): """汽车行驶数据分析（国赛 2015 C）。 生成 4 类典型工况 × 各 3 段行驶周期，每段为速度/电流/温度时序； 提取 6 维工况特征矩阵；注入 2 个异常周期（电流传感器尖峰）。 求解：K-Means(k-means++ 种子) + 轮廓系数选 K；马氏距离 + z-score + 重构残差 三法异常检测。 """ random.seed(seed) kinds = ["城市拥堵", "经济巡航", "高速行驶", "激进驾驶"] proto = { "城市拥堵": [18.0, 16.0, 35.0, 0.60, 0.35, 45.0], "经济巡航": [42.0, 10.0, 65.0, 0.30, 0.10, 32.0], "高速行驶": [85.0, 8.0, 120.0, 0.25, 0.03, 55.0], "激进驾驶": [55.0, 22.0, 140.0, 1.10, 0.08, 70.0], } assign = [] for kd in kinds: assign += [kd] * 3 random.shuffle(assign) U = 350.0 # 母线电压 V cycles = [] records = [] # 长时序 (t, speed, accel, current, temp, cid) feats = [] # N × 6 t = 0 for cid, kd in enumerate(assign): base = proto[kd] mean_v, noise = base[0], base[1] * 0.12 max_v = base[2] speed = [] v = mean_v for _ in range(L): v += 0.5 * (mean_v - v) + random.gauss(0, noise) v = max(0.0, min(max_v, v)) speed.append(v) accel = [speed[i + 1] - speed[i] for i in range(L - 1)] pos_acc = [a for a in accel if a > 0] current = [max(2.0, base[5] + 0.15 * (v - mean_v) + 3.0 * max(0.0, a) + random.gauss(0, 1.2)) for v, a in zip(speed[1:], accel)] temp = [25.0 + 0.06 * v + random.gauss(0, 0.5) for v in speed[1:]] energy = sum(c * U for c in current) * 1.0 / 3600.0 / 1000.0 # kWh, dt=1s feat = [ _mean(speed), _std(speed), max(speed), _mean(pos_acc) if pos_acc else 0.0, sum(1 for s in speed if s < 1.0) / L, _mean(current), ] feats.append(feat) cycles.append({"cid": cid, "kind": kd, "speed": speed, "current": current, "temp": temp, "energy": energy, "feat": feat, "is_anom": False}) for i in range(1, L): records.append((t, round(speed[i], 2), round(accel[i - 1], 3), round(current[i - 1], 2), round(temp[i - 1], 2), cid)) t += 1 # 注入 2 个异常周期：电流传感器故障（一段连续采样拉到 ~350A，远超所有正常工况） anom_idx = [4, 9] for ai in anom_idx: cyc = cycles[ai] for j in range(25, 85): cyc["current"][j] = 350.0 + random.gauss(0, 6.0) cyc["is_anom"] = True cyc["feat"][5] = _mean(cyc["current"]) # 平均电流被显著抬高 cyc["energy"] = sum(c * U for c in cyc["current"]) / 3600.0 / 1000.0 # 重算特征矩阵（含异常） feats = [c["feat"] for c in cycles] feat_names = ["均速", "速标准差", "最高速", "平均正加速", "怠速比", "平均电流"] Z, fmeans, fstds = _zscale(feats) # ---------- 异常检测（先剔除故障样本，再聚类）---------- # 方法一：稳健马氏距离（中位数/MAD 对角近似）。 # 常规协方差会被样本内异常污染而失效，故以中位数估计位置、MAD 估计尺度， # 对"单变量极端离群"（如电流传感器故障）稳健。卡方(6) 0.975 分位 ≈ 14.45。 med = [_median([Z[i][j] for i in range(N)]) for j in range(6)] mad = [_mad([Z[i][j] for i in range(N)], med[j]) * 1.4826 for j in range(6)] maha = [sum(((Z[i][j] - med[j]) / (mad[j] or 1.0)) ** 2 for j in range(6)) for i in range(N)] thr = 14.45 flag_maha = [m > thr for m in maha] # 方法二：平均电流单变量 z 分数（快速筛查，对污染协方差敏感，作对照） cur = [f[5] for f in feats] cm, cs = _mean(cur), _std(cur) zcur = [(c - cm) / cs for c in cur] flag_z = [abs(z) > 3.0 for z in zcur] # 以稳健马氏结果作为剔除依据，得到干净样本集 clean_idx = [i for i in range(N) if not flag_maha[i]] Z_clean = [Z[i] for i in clean_idx] # ---------- 在干净集上聚类分群（K=2..6 轮廓系数选 K）---------- sil = {}; inertia = {}; models = {} for K in range(2, 7): lab_c, cents_c, iner = _kmeans(Z_clean, K, seed=seed) sil[K] = _silhouette(Z_clean, lab_c, K) inertia[K] = iner models[K] = (lab_c, cents_c) # 选 K：轮廓系数达拐点即停（继续增 K 的增量 <0.03/步视为过度细分），兼顾可解释性 bestK = 2 for K in range(3, 7): if sil[K] - sil[bestK] > 0.03: bestK = K else: break lab_c, cents = models[bestK] # 把干净标签映射回全样本索引 full_lab = [None] * N for pos, i in enumerate(clean_idx): full_lab[i] = lab_c[pos] for i in range(N): cycles[i]["cluster"] = full_lab[i] if full_lab[i] is not None else -1 # 簇→工况 对齐（原始特征空间）：先用干净样本的真实标签求各工况经验中心， # 再把簇心就近匹配——避免手工原型与生成特征口径不一致导致的错配。 raw_cents = [] for cl in range(bestK): grp = [feats[clean_idx[pos]] for pos in range(len(clean_idx)) if lab_c[pos] == cl] raw_cents.append([_mean(col) for col in zip(*grp)] if grp else [0.0] * 6) kind_center = {} for kd in kinds: grp = [feats[clean_idx[pos]] for pos in range(len(clean_idx)) if cycles[clean_idx[pos]]["kind"] == kd] kind_center[kd] = [_mean(col) for col in zip(*grp)] if grp else [0.0] * 6 align = {} for cl in range(bestK): align[cl] = min(kinds, key=lambda kd: _dist2(raw_cents[cl], kind_center[kd])) # 各簇能耗（仅干净样本） cl_energy = {} for cl in range(bestK): es = [cycles[clean_idx[pos]]["energy"] for pos in range(len(clean_idx)) if lab_c[pos] == cl] cl_energy[cl] = _mean(es) if es else 0.0 # 方法三：重构残差（到最近干净簇心的距离，异常远离所有正常簇心） res = [] for i in range(N): d = min(_dist2(Z[i], cents[cl]) for cl in range(bestK)) res.append(d) rthr = sorted(res)[max(0, N - 3)] # 第 (N-2) 大距离作阈值，约标最远 2 个 flag_res = [r > rthr for r in res] # 三法合并（任意一法命中即判异常） flag_any = [flag_maha[i] or flag_z[i] or flag_res[i] for i in range(N)] tp = sum(1 for i in range(N) if cycles[i]["is_anom"] and flag_any[i]) fp = sum(1 for i in range(N) if (not cycles[i]["is_anom"]) and flag_any[i]) fn = sum(1 for i in range(N) if cycles[i]["is_anom"] and not flag_any[i]) detect = {"maha": maha, "thr": thr, "flag_maha": flag_maha, "zcur": zcur, "flag_z": flag_z, "res": res, "rthr": rthr, "flag_res": flag_res, "flag_any": flag_any, "tp": tp, "fp": fp, "fn": fn, "recall": tp / max(1, tp + fn), "precision": tp / max(1, tp + fn + fp)} # 鲁棒性：特征加 ±5% 噪声重算稳健马氏距离，看异常是否仍被稳定标出且无新增误报 robust_stable = 0 for _ in range(50): Zp = [[Z[i][j] * (1 + random.gauss(0, 0.05)) for j in range(6)] for i in range(N)] mp = [_median([Zp[i][j] for i in range(N)]) for j in range(6)] mdp = [_mad([Zp[i][j] for i in range(N)], mp[j]) * 1.4826 for j in range(6)] map_ = [sum(((Zp[i][j] - mp[j]) / (mdp[j] or 1.0)) ** 2 for j in range(6)) for i in range(N)] fl = [m > thr for m in map_] if all(fl[i] for i in anom_idx) and not any(fl[i] for i in range(N) if i not in anom_idx): robust_stable += 1 robust_rate = robust_stable / 50.0 # 写出 CSV（长时序 series 型） path = os.path.join(OUT, out_csv) with open(path, "w", encoding="utf-8-sig", newline="") as f: w = csv.writer(f) w.writerow(["t", "speed", "accel", "current", "temp", "cycle"]) for r in records: w.writerow(r) # 诊断打印 print("gen_2015c ->", path, "cycles=", N, "L=", L, "records=", len(records)) print(" 工况分配=", assign) print(" 特征均值(行=周期):") for i in range(N): print(" cid%d %s feat=%s energy=%.3f anom=%s" % (i, cycles[i]["kind"], ["%.1f" % x for x in feats[i]], cycles[i]["energy"], cycles[i]["is_anom"])) print(" 轮廓系数(干净集) K=2..6:", {k: round(sil[k], 3) for k in sorted(sil)}) print(" 最优 K =", bestK, " 干净样本数=", len(clean_idx)) print(" 簇→工况对齐:", align) print(" 各簇平均能耗(kWh):", {cl: round(cl_energy[cl], 3) for cl in cl_energy}) print(" 稳健马氏距离=", ["%.1f" % m for m in maha], "阈值=%.2f" % thr) print(" z(平均电流)=", ["%.2f" % z for z in zcur]) print(" 重构残差=", ["%.2f" % r for r in res], "阈值=%.2f" % rthr) print(" 异常标记 maha/res/z/合并=", flag_maha, flag_res, flag_z, flag_any) print(" 检测: TP=%d FP=%d FN=%d recall=%.2f precision=%.2f" % (tp, fp, fn, detect["recall"], detect["precision"])) print(" 鲁棒性(±5%%扰动仍稳定标出)=%.2f" % robust_rate) res = dict(N=N, L=L, kinds=kinds, assign=assign, cycles=cycles, records=records, feat_names=feat_names, feats=feats, fmeans=fmeans, fstds=fstds, Z=Z, Z_clean=Z_clean, clean_idx=clean_idx, full_lab=full_lab, sil=sil, inertia=inertia, bestK=bestK, lab=lab_c, cents=cents, raw_cents=raw_cents, align=align, cl_energy=cl_energy, detect=detect, robust_rate=robust_rate, anom_idx=anom_idx, U=U, Krange=list(range(2, 7))) return res # ===================================================================== # 2016 A 系泊系统的设计（张紧链静平衡 + 蒙特卡洛载荷 + 多目标优化） # 物理模型（合成、量级合理）：浮体以净浮力 Fb 浮于水面，锚链自浮体下接至 # 海床锚点（水深 d）。锚链视为张紧直线段，在水平风浪载荷 Fh 下做准静力平衡。 # 闭式解：x = Fh*d/sqrt(Wc^2 - Fh^2)（Wc=链总重），要求 Wc>Fh 否则浮体被拖走； # 垂直分力 Tv=Wc*d/s，若 Tv>Fb 则浮体被拉沉（下沉约束）。 # 决策变量：链节数 n；目标：位移/倾角最小且造价最低（多目标）；载荷用 MC 模拟。 # ===================================================================== def gen_2016a(seed=2016, out_csv="cumcm2016a.csv"): import csv, os random.seed(seed) g = 9.8 # ---- 固定物理参数（合成数据，量级贴近小型近海浮标）---- d = 18.0 # 水深 m Fb = 5000.0 # 浮体净浮力 N seg_len = 0.6 # 单节链长 m seg_mass = 10.0 # 单节链质量 kg w_lin = seg_mass * g / seg_len # 链线密度 N/m cost_per = 120.0 # 单节造价 元 Fh0 = 800.0 # 标称水平载荷 N n_min, n_max = 10, 55 ns = list(range(n_min, n_max + 1)) def Wc(n): return w_lin * seg_len * n # 链总重 N = seg_mass*g*n def equil(n, Fh): # 返回 (x 水平位移, phi 倾角, Tv 垂直分力, 可行(未沉), 拖走) W = Wc(n) if W <= Fh: return (float("inf"), None, float("inf"), False, True) x = Fh * d / math.sqrt(W * W - Fh * Fh) s = math.sqrt(x * x + d * d) phi = math.atan2(x, d) Tv = W * d / s return (x, phi, Tv, Tv <= Fb, False) # ---- 标称载荷下各 n 表现（用于图与表）---- nom = {} for n in ns: x, phi, Tv, ok, drag = equil(n, Fh0) nom[n] = dict(n=n, Wc=Wc(n), x=x, phi=phi, Tv=Tv, sink=(not ok), drag=drag, cost=n * cost_per) # ---- 蒙特卡洛载荷：常态 N(800,300) + 10% 风暴 N(2500,600) ---- Nmc = 400 loads = [] for _ in range(Nmc): if random.random() < 0.10: loads.append(max(50.0, random.gauss(2500, 600))) else: loads.append(max(50.0, random.gauss(800, 300))) mc = {} for n in ns: xs = []; drags = 0; sinks = 0 for Fh in loads: x, phi, Tv, ok, drag = equil(n, Fh) if drag: drags += 1; continue if not ok: sinks += 1; continue xs.append(x) if xs: sx = sorted(xs) p95 = sx[max(0, len(sx) - 1 - int(0.05 * len(sx)))] mc[n] = dict(mean=_mean(xs), p95=p95, pmax=max(xs), pdrag=drags / Nmc, psink=sinks / Nmc, nfeas=len(xs)) else: mc[n] = dict(mean=float("inf"), p95=float("inf"), pmax=float("inf"), pdrag=drags / Nmc, psink=sinks / Nmc, nfeas=0) # ---- 优化：可行域 = pdrag==0 且 psink<=0.02；域内取 MC 平均位移最小 ---- feas = [n for n in ns if mc[n]["pdrag"] == 0 and mc[n]["psink"] <= 0.02] best_n = min(feas, key=lambda n: mc[n]["mean"]) if feas else None # 多目标：位移与造价归一化，取综合最小（knee） if feas: disp = [mc[n]["mean"] for n in feas] cost = [nom[n]["cost"] for n in feas] dmin, dmax = min(disp), max(disp) cmin, cmax = min(cost), max(cost) score = {n: (mc[n]["mean"] - dmin) / (dmax - dmin or 1) + (nom[n]["cost"] - cmin) / (cmax - cmin or 1) for n in feas} knee_n = min(score, key=lambda n: score[n]) else: knee_n = best_n # 折衷：在 best_n 基础上允许 +1 节误差的邻域最小造价 opt_cost_n = min(feas, key=lambda n: nom[n]["cost"]) if feas else None # ---- 灵敏度：浮体浮力 Fb、水深 d、链密度（seg_mass）对最优节数影响 ---- def best_for(Fb_, d_, sm_): wl = sm_ * g / seg_len Wc_ = lambda n: wl * seg_len * n def eq(n, Fh): W = Wc_(n) if W <= Fh: return (float("inf"), False, True) x = Fh * d_ / math.sqrt(W * W - Fh * Fh) s = math.sqrt(x * x + d_ * d_) Tv = W * d_ / s return (x, Tv <= Fb_, False) local = [] for n in ns: xs = [] for Fh in loads: x, ok, dr = eq(n, Fh) if dr or not ok: continue xs.append(x) if xs and len(xs) == len(loads): local.append((n, _mean(xs))) return min(local, key=lambda t: t[1])[0] if local else None sens = { "Fb": [(f, best_for(f, d, seg_mass)) for f in (4000, 4500, 5000, 5500, 6000)], "depth": [(dd, best_for(Fb, dd, seg_mass)) for dd in (14, 16, 18, 20, 22)], "seg_mass": [(sm, best_for(Fb, d, sm)) for sm in (8, 9, 10, 11, 12)], } # ---- 情景对比：常态 vs 风暴占比提升（30% 风暴）---- loads_storm = [max(50.0, random.gauss(2500, 600)) if random.random() < 0.30 else max(50.0, random.gauss(800, 300)) for _ in range(Nmc)] storm = {} for n in ns: xs = []; drags = 0; sinks = 0 for Fh in loads_storm: x, phi, Tv, ok, drag = equil(n, Fh) if drag: drags += 1; continue if not ok: sinks += 1; continue xs.append(x) storm[n] = dict(mean=_mean(xs) if xs else float("inf"), pdrag=drags / Nmc, psink=sinks / Nmc) # ---- 写出 CSV（params 型：每节数一行性能）---- path = os.path.join(OUT, out_csv) with open(path, "w", encoding="utf-8-sig", newline="") as f: w = csv.writer(f) w.writerow(["n", "Wc_N", "x_nominal_m", "mean_disp_m", "p95_disp_m", "psink", "pdrag", "cost_yuan", "feasible"]) for n in ns: m = mc[n] w.writerow([n, round(Wc(n), 1), round(nom[n]["x"], 3), round(m["mean"], 3) if m["mean"] != float("inf") else "inf", round(m["p95"], 3) if m["p95"] != float("inf") else "inf", round(m["psink"], 3), round(m["pdrag"], 3), nom[n]["cost"], int(n in feas)]) # 诊断打印 print("gen_2016a ->", path, " 水深=%.1f 浮力=%.0f 链线密度=%.1f N/m MC=%d" % (d, Fb, w_lin, Nmc)) print(" 标称载荷 Fh=%.0f 时: n=20 位移=%.2f m, n=35 位移=%.2f m, n=50 位移=%.2f m" % (Fh0, nom[20]["x"], nom[35]["x"], nom[50]["x"])) print(" 最优节数(位移最小) n* =", best_n, " 该 n MC 平均位移=%.3f m, p95=%.3f m" % (mc[best_n]["mean"], mc[best_n]["p95"]) if best_n else "无") print(" 多目标 knee n =", knee_n, " 最省可行 n =", opt_cost_n) print(" 灵敏度→最优n: Fb:%s depth:%s seg_mass:%s" % (sens["Fb"], sens["depth"], sens["seg_mass"])) res = dict(d=d, Fb=Fb, seg_len=seg_len, seg_mass=seg_mass, w_lin=w_lin, cost_per=cost_per, Fh0=Fh0, n_min=n_min, n_max=n_max, ns=ns, nom=nom, loads=loads, Nmc=Nmc, mc=mc, feas=feas, best_n=best_n, knee_n=knee_n, opt_cost_n=opt_cost_n, sens=sens, storm=storm, Wc=Wc, equil=equil) return res # --------------------------------------------------------------------------- # 2016B 小区开放对周边路网的影响（交通流分配 / 用户均衡 UE） # 模型：BPR 路阻函数 + 用户均衡（MSA 迭代）+ 开放前后对比 # --------------------------------------------------------------------------- def _dijkstra(nodes, adj, src, weight_fn): """最短路（Dijkstra）。weight_fn(edge)->非负权重。返回 dist, prev。""" dist = {i: float("inf") for i in nodes} prev = {i: None for i in nodes} dist[src] = 0.0 U = set(nodes) while U: u = min(U, key=lambda x: dist[x]) if dist[u] == float("inf"): break U.discard(u) for (v, e) in adj[u]: w = weight_fn(e) if dist[u] + w < dist[v]: dist[v] = dist[u] + w prev[v] = u return dist, prev def gen_2016b(seed=2016, out_csv="cumcm2016b.csv"): """小区开放交通仿真。 路网：外围 4 个角节点 A,B,C,D 构成主干道环；中心 4 节点 E,F,G,H 构成 小区内部 2x2 网格。角-内部之间由 8 条“门”边连接。 - 封闭模式：门边与内部支路被物理阻隔（不可用）。 - 开放模式：全部边可用，内部支路成为穿行捷径。 OD 需求为四角之间的出行（无向 6 对）。用 BPR 路阻 + MSA 用户均衡分配， 对比开放前后系统总出行成本、主干道负载率、平均行程时间。 """ random.seed(seed) ALPHA, BETA = 0.15, 4.0 # BPR 标准参数 # 节点坐标（0~100 网格） coords = {"A": (0, 50), "B": (100, 50), "C": (100, 0), "D": (0, 0), "E": (33, 33), "F": (67, 33), "G": (33, 67), "H": (67, 67)} nodes = list(coords.keys()) def d2(a, b): return math.hypot(coords[a][0] - coords[b][0], coords[a][1] - coords[b][1]) # 边定义：(u, v, free_time, capacity, is_arterial) # free_time 用 距离*系数；主干道容量大、系数小；支路容量小。 raw_edges = [ ("A", "B", 0.45, 60.0, True), # 主干道上 ("B", "C", 0.45, 60.0, True), # 主干道右 ("C", "D", 0.45, 60.0, True), # 主干道下 ("D", "A", 0.45, 60.0, True), # 主干道左 ("E", "F", 0.18, 60.0, False), # 内部水平 ("G", "H", 0.18, 60.0, False), ("E", "G", 0.18, 60.0, False), # 内部垂直 ("F", "H", 0.18, 60.0, False), ("A", "E", 0.16, 60.0, False), # 8 条门边 ("A", "G", 0.16, 60.0, False), ("B", "F", 0.16, 60.0, False), ("B", "H", 0.16, 60.0, False), ("C", "F", 0.16, 60.0, False), ("C", "H", 0.16, 60.0, False), ("D", "E", 0.16, 60.0, False), ("D", "G", 0.16, 60.0, False), ] # 边 id 表 edge_list = [] for (u, v, t0, cap, art) in raw_edges: edge_list.append({"u": u, "v": v, "t0": t0, "cap": cap, "art": art, "d": d2(u, v)}) # OD 需求（四角之间，无向 6 对，单位 veh/h） OD = [("A", "B", 28.0), ("A", "C", 22.0), ("A", "D", 18.0), ("B", "C", 26.0), ("B", "D", 20.0), ("C", "D", 24.0)] def is_avail(e, opened, open_edges): if e["art"]: return True if open_edges is not None: return (e["u"], e["v"]) in open_edges or (e["v"], e["u"]) in open_edges return opened def build_graph(opened, open_edges=None): """返回 adj。internal 边可用与否取决于 opened 或 open_edges 集合。""" adj = {i: [] for i in nodes} for e in edge_list: if not is_avail(e, opened, open_edges): continue adj[e["u"]].append((e["v"], e)) adj[e["v"]].append((e["u"], e)) return adj def bpr(e, v): return e["t0"] * (1.0 + ALPHA * (v / e["cap"]) ** BETA) def assign(adj, flows): """一次全有全无分配（基于当前路阻最短路），返回各边增量与 OD 路径耗时。""" # 当前权重 def w(e): v = flows.get(id(e), 0.0) return bpr(e, v) inc = {} od_time = {} for (o, d, q) in OD: dist, prev = _dijkstra(nodes, adj, o, w) # 回溯路径 path = [] cur = d while cur is not None: path.append(cur) cur = prev[cur] path.reverse() # 累加流量到路径各边 for a, b in zip(path, path[1:]): # 找到对应边 for (v, e) in adj[a]: if v == b: inc[id(e)] = inc.get(id(e), 0.0) + q break od_time[(o, d)] = dist[d] return inc, od_time def equilibrium(opened, open_edges=None, iters=300): adj = build_graph(opened, open_edges) flows = {id(e): 0.0 for e in edge_list} # MSA 迭代 for it in range(1, iters + 1): inc, _ = assign(adj, flows) step = 1.0 / it for e in edge_list: k = id(e) if k in inc: flows[k] = flows[k] + step * (inc[k] - flows[k]) # 计算指标 total_cost = 0.0 edge_load = {} arterial_load = {} for e in edge_list: k = id(e) v = flows.get(k, 0.0) if not is_avail(e, opened, open_edges): continue t = bpr(e, v) edge_load[(e["u"], e["v"])] = (v, e["cap"], v / e["cap"], t) total_cost += v * t if e["art"]: arterial_load[(e["u"], e["v"])] = (v, e["cap"], v / e["cap"]) return flows, total_cost, edge_load, arterial_load f_closed, tc_closed, el_closed, art_closed = equilibrium(False) f_open, tc_open, el_open, art_open = equilibrium(True) # 半开放情景：开放全部 8 条门边 + 仅 1 条内部水平支路（E-F），关闭其余 3 条内部支路 half_edges = {("A", "E"), ("A", "G"), ("B", "F"), ("B", "H"), ("C", "F"), ("C", "H"), ("D", "E"), ("D", "G"), ("E", "F")} f_half, tc_half, el_half, art_half = equilibrium(True, open_edges=half_edges) # 平均行程时间（按 OD 需求加权） def avg_travel(opened, flows, open_edges=None): adj = build_graph(opened, open_edges) def w(e): return bpr(e, flows.get(id(e), 0.0)) tot_q = 0.0 tot_t = 0.0 for (o, d, q) in OD: dist, _ = _dijkstra(nodes, adj, o, w) tot_q += q tot_t += q * dist[d] return tot_t / tot_q at_closed = avg_travel(False, f_closed) at_open = avg_travel(True, f_open) at_half = avg_travel(True, f_half, open_edges=half_edges) # 主干道 4 边在 封闭/半开放/全开放 三种情景的负载率对比 arteries = [("A", "B"), ("B", "C"), ("C", "D"), ("D", "A")] art_compare = [] for (a, b) in arteries: vc, cc, rc = art_closed[(a, b)] vh, ch, rh = art_half[(a, b)] vo, co, ro = art_open[(a, b)] art_compare.append((a + b, vc, cc, rc, vh, ch, rh, vo, co, ro)) art_avg_closed = _mean([art_closed[a][2] for a in arteries]) art_avg_half = _mean([art_half[a][2] for a in arteries]) art_avg_open = _mean([art_open[a][2] for a in arteries]) # 内部支路开放后的新增流量（证明小区内部承担了穿行流量） internal_edges = [("E", "F"), ("G", "H"), ("E", "G"), ("F", "H")] internal_flow = [] for (a, b) in internal_edges: if (a, b) in el_open: v, c, r, t = el_open[(a, b)] internal_flow.append((a + b, v, c, r)) # 半开放下内部支路流量（对比） internal_flow_half = [] for (a, b) in internal_edges: if (a, b) in el_half: v, c, r, t = el_half[(a, b)] internal_flow_half.append((a + b, v, c, r)) # 灵敏度：OD 总需求强度 s（0.6~1.4）对平均行程时间与主干道负载率的影响 base_demand = [list(x) for x in OD] total_q = sum(q for (_, _, q) in OD) sens_demand = [] for s in [0.6, 0.85, 1.0, 1.15, 1.4]: newOD = [(o, d, q * s) for (o, d, q) in base_demand] oldOD = OD[:] OD[:] = newOD f_s, _, _, art = equilibrium(True) at_s = avg_travel(True, f_s) avg_r = _mean([art[a][2] for a in art]) if art else 0.0 sens_demand.append((round(s, 2), round(at_s, 4), round(avg_r, 4))) OD[:] = oldOD # 通行效率主指标：平均行程时间（用户视角）。下降率基于此。 relief = (at_closed - at_open) / at_closed relief_half = (at_closed - at_half) / at_closed res = dict(coords=coords, edges=edge_list, OD=OD, ALPHA=ALPHA, BETA=BETA, tc_closed=tc_closed, tc_open=tc_open, tc_half=tc_half, relief=relief, relief_half=relief_half, at_closed=at_closed, at_open=at_open, at_half=at_half, art_compare=art_compare, art_avg_closed=art_avg_closed, art_avg_half=art_avg_half, art_avg_open=art_avg_open, internal_flow=internal_flow, internal_flow_half=internal_flow_half, sens_demand=sens_demand) # 写练习数据 CSV try: rows = [] for e in edge_list: rows.append([e["u"], e["v"], round(e["t0"], 4), e["cap"], int(e["art"])]) od_rows = [[o, d, q] for (o, d, q) in OD] os.makedirs("assets/problems/data", exist_ok=True) with open("assets/problems/data/" + out_csv, "w", newline="") as fp: w = csv.writer(fp) w.writerow(["type", "from", "to", "free_time", "capacity", "arterial"]) for r in rows: w.writerow(["edge"] + r) for r in od_rows: w.writerow(["od", r[0], r[1], r[2], "", ""]) except Exception as ex: print(" 写 CSV 失败:", ex) print("gen_2016b:") print(" 封闭平均行程时间=%.4f 半开放=%.4f 全开放=%.4f" % (at_closed, at_half, at_open)) print(" 行程时间下降率 半开放/全开放=%.3f/%.3f" % (relief_half, relief)) print(" 路网负荷成本tc 封闭=%.3f 半开放=%.3f 全开放=%.3f" % (tc_closed, tc_half, tc_open)) print(" 主干道均负载率 封闭=%.3f 半开放=%.3f 全开放=%.3f" % (art_avg_closed, art_avg_half, art_avg_open)) print(" 主干道开放前后负载率:", art_compare) print(" 内部支路开放后流量:", internal_flow) print(" 灵敏度(需求倍率,平均行程时间,主干均负载率):", sens_demand) return res def gen_2016c(seed=2016, out_csv="cumcm2016c.csv"): """电池剩余放电时间预测（2016C）。 物理仿真一个动力电池包放电过程，并实现三种 SOC 估计 + 剩余放电时间/里程预测。 返回结果字典，所有数字确定性可复现；同时写出练习用 CSV。""" random.seed(seed) Ns = 100 # 串联电芯数 Q = 60.0 # 包容量 Ah dt = 12.0 # 采样间隔 s SOC_cut = 0.05 # 放电截止 SOC def ocv_cell(s): return 3.3 + 1.15 * s - 0.15 * s * s def R_pack(s, T): return (0.14 + 0.26 * (1.0 - s)) * (1.0 - 0.006 * (T - 25.0)) def ocv_inv(oc): lo, hi = 0.0, 1.0 for _ in range(40): mid = 0.5 * (lo + hi) if ocv_cell(mid) < oc: lo = mid else: hi = mid return 0.5 * (lo + hi) # ---- 仿真放电过程 ---- t, I, U, Tv, SOC, Uocv = [], [], [], [], [], [] soc = 1.0 k = 0 while True: tt = k * dt cur = 55 + 35 * math.sin(2 * math.pi * tt / 900.0) + 22 * math.sin(2 * math.pi * tt / 330.0 + 1.0) + 12 * math.sin(2 * math.pi * tt / 140.0) cur = clip(cur, 8.0, 170.0) Tcur = 25.0 + 6.0 * (cur / 170.0) ** 2 + 0.5 * math.sin(tt / 500.0) + random.gauss(0, 0.3) Uoc = Ns * ocv_cell(soc) R = R_pack(soc, Tcur) Ucur = Uoc - cur * R + random.gauss(0, 0.7) t.append(tt); I.append(cur); Tv.append(Tcur); SOC.append(soc); Uocv.append(Uoc); U.append(Ucur) soc -= cur * dt / (3600.0 * Q) k += 1 if soc <= SOC_cut or k > 5000: break N = len(t) T_total = t[-1] # ---- 方法1：安时积分 + OCV 修正 ---- Iw = [I[i] + random.gauss(0, 1.5) for i in range(N)] soc_c = [0.0] * N soc_c[0] = 1.03 soc_e1 = [0.0] * N for i in range(1, N): soc_c[i] = soc_c[i - 1] - Iw[i] * dt / (3600.0 * Q) Rc = R_pack(soc_c[i], Tv[i]) oc_est = (U[i] + Iw[i] * Rc) / Ns soc_ocv = ocv_inv(oc_est) soc_c[i] += 0.05 * (soc_ocv - soc_c[i]) # 周期性向 OCV 锚定 soc_e1[i] = soc_c[i] soc_e1[0] = 1.03 m1_mae = sum(abs(soc_e1[i] - SOC[i]) for i in range(N)) / N # ---- 方法2：神经网络(模拟) + 指数平滑外推 ---- soc_e2 = [0.0] * N for i in range(N): Rc = R_pack(SOC[i], Tv[i]) oc_clean = (U[i] + I[i] * Rc) / Ns soc_e2[i] = ocv_inv(oc_clean) + random.gauss(0, 0.004) m2_mae = sum(abs(soc_e2[i] - SOC[i]) for i in range(N)) / N # ---- 方法3：卡尔曼融合 ---- soc_e3 = [0.0] * N x = 1.0 q, r = 2.0e-6, 4.0e-4 P = 1.0 soc_e3[0] = x for i in range(1, N): x = x - Iw[i] * dt / (3600.0 * Q) P = P + q Rc = R_pack(x, Tv[i]) z = ocv_inv((U[i] + Iw[i] * Rc) / Ns) K = P / (P + r) x = x + K * (z - x) P = (1 - K) * P soc_e3[i] = x m3_mae = sum(abs(soc_e3[i] - SOC[i]) for i in range(N)) / N # ---- 查询点与剩余放电时间预测 ---- fracs = [0.2, 0.4, 0.6, 0.8] qidx = [int(fr * N) for fr in fracs] Uavg = sum(Uocv) / N cons = 0.16 # kWh/km mu_all = sum(Iw) / N sw_all = math.sqrt(sum((v - mu_all) ** 2 for v in Iw) / N) queries = [] for qi in qidx: tq = t[qi] true_rem = T_total - tq w0 = max(0, qi - 50) Iavg = sum(Iw[w0:qi + 1]) / (qi + 1 - w0) m1_rem = soc_e1[qi] * Q * 3600.0 / max(Iavg, 1.0) a = max(0, qi - 39) xs = list(range(qi - a + 1)); ys = soc_e2[a:qi + 1] n2 = len(xs); mx = sum(xs) / n2; my = sum(ys) / n2 denom = sum((xs[j] - mx) ** 2 for j in range(n2)) b = sum((xs[j] - mx) * (ys[j] - my) for j in range(n2)) / denom if denom > 0 else 0 m2_rem = (soc_e2[qi] - SOC_cut) * (-dt / b) if (b < 0 and soc_e2[qi] > SOC_cut) else true_rem M = 200; samples = []; s_now = soc_e3[qi] for _ in range(M): s = s_now; tt_rem = 0.0 while s > SOC_cut and tt_rem < 2 * T_total: ci = clip(mu_all + random.gauss(0, sw_all), 5.0, 180.0) s -= ci * dt / (3600.0 * Q) tt_rem += dt samples.append(tt_rem) m3_rem = sum(samples) / M samples.sort() lo = samples[int(0.025 * M)]; hi = samples[int(0.975 * M)] mile_true = SOC[qi] * Q * Uavg / 1000.0 / cons mile_m1 = soc_e1[qi] * Q * Uavg / 1000.0 / cons mile_m2 = soc_e2[qi] * Q * Uavg / 1000.0 / cons mile_m3 = soc_e3[qi] * Q * Uavg / 1000.0 / cons queries.append(dict(t=tq, true_rem=true_rem, m1_rem=m1_rem, m2_rem=m2_rem, m3_rem=m3_rem, m3_lo=lo, m3_hi=hi, mile_true=mile_true, mile_m1=mile_m1, mile_m2=mile_m2, mile_m3=mile_m3, soc_true=SOC[qi], soc_e1=soc_e1[qi], soc_e2=soc_e2[qi], soc_e3=soc_e3[qi])) def metrics(key): errs = [queries[i][key] - queries[i]["true_rem"] for i in range(len(queries))] mae = sum(abs(e) for e in errs) / len(errs) rmse = math.sqrt(sum(e * e for e in errs) / len(errs)) return mae, rmse m1_mae_r, m1_rmse = metrics("m1_rem") m2_mae_r, m2_rmse = metrics("m2_rem") m3_mae_r, m3_rmse = metrics("m3_rem") cons_summary = [ ("安时积分+OCV", m1_mae_r, m1_rmse, m1_mae), ("神经网络+平滑", m2_mae_r, m2_rmse, m2_mae), ("卡尔曼+蒙特卡洛", m3_mae_r, m3_rmse, m3_mae), ] ocv_curve = [(s, Ns * ocv_cell(s)) for s in [i / 20.0 for i in range(0, 21)]] with open(os.path.join(OUT, out_csv), "w", newline="") as f: w = csv.writer(f) w.writerow(["t", "I", "U", "T", "SOC_true", "Uocv"]) for i in range(N): w.writerow([round(t[i], 1), round(I[i], 2), round(U[i], 2), round(Tv[i], 2), round(SOC[i], 4), round(Uocv[i], 2)]) return dict(seed=seed, N=N, dt=dt, Ns=Ns, Q=Q, SOC_cut=SOC_cut, t=t, I=I, U=U, T=Tv, SOC_true=SOC, Uocv=Uocv, T_total=T_total, soc_e1=soc_e1, soc_e2=soc_e2, soc_e3=soc_e3, m1_mae=m1_mae, m2_mae=m2_mae, m3_mae=m3_mae, queries=queries, cons_summary=cons_summary, m1_mae_r=m1_mae_r, m1_rmse=m1_rmse, m2_mae_r=m2_mae_r, m2_rmse=m2_rmse, m3_mae_r=m3_mae_r, m3_rmse=m3_rmse, ocv_curve=ocv_curve, Uavg=Uavg, cons=cons) def gen_2017a(seed=2017, N=48, M=120, L=71, delta_true=4.0, out_csv="cumcm2017a.csv"): """2017A CT 系统参数标定与重建（平行束 Radon）。 生成含探测器偏移 δ_true 的 sinogram，分别用三种方法标定偏移并重建， 返回权威数字字典；并写出 params 型 CSV（sinogram 平铺）。""" import math as _m import cmath as _cm rnd = random.Random(seed) # ---- 几何网格与探测器排布 ---- cx = cy = N / 2.0 step = 1.0 half = (L - 1) / 2.0 det = [ (k - half) * step for k in range(L) ] # 探测器坐标，以旋转中心为 0 angles = [ _m.pi * i / M for i in range(M) ] # 0 ~ 180° def _sample(f, x, y): if x < 0 or x > N - 1 or y < 0 or y > N - 1: return 0.0 x0 = int(x); y0 = int(y); x1 = min(x0 + 1, N - 1); y1 = min(y0 + 1, N - 1) fx = x - x0; fy = y - y0 return (f[y0][x0] * (1 - fx) * (1 - fy) + f[y0][x1] * fx * (1 - fy) + f[y1][x0] * (1 - fx) * fy + f[y1][x1] * fx * fy) def _sample1d(row, t): idx = (t - det[0]) / step if idx < 0 or idx > L - 1: return 0.0 j = int(idx); j1 = min(j + 1, L - 1); f = idx - j return row[j] * (1 - f) + row[j1] * f def _add_disk(f, ox, oy, r, val): for y in range(N): for x in range(N): if (x - ox) ** 2 + (y - oy) ** 2 <= r * r: f[y][x] = max(f[y][x], val) def _add_rect(f, ox, oy, w, h, val): for y in range(N): for x in range(N): if abs(x - ox) <= w and abs(y - oy) <= h: f[y][x] = max(f[y][x], val) def build_phantom(kind): f = [[0.0] * N for _ in range(N)] if kind == "template": _add_disk(f, cx - 9, cy - 7, 9, 0.95) _add_disk(f, cx + 11, cy + 9, 6, 0.70) _add_rect(f, cx + 2, cy - 4, 5, 9, 0.55) _add_disk(f, cx - 2, cy + 12, 3, 0.45) else: # unknown _add_disk(f, cx - 6, cy + 5, 8, 0.88) _add_disk(f, cx + 10, cy - 8, 5, 0.62) _add_rect(f, cx - 3, cy - 6, 7, 4, 0.50) _add_disk(f, cx + 4, cy + 10, 4, 0.40) return f def line_integral(f, theta, t): e0 = _m.cos(theta); e1 = _m.sin(theta) n0 = -e1; n1 = e0 R = N * 0.98; s = -R; tot = 0.0 while s <= R: x = cx + t * e0 + s * n0; y = cy + t * e1 + s * n1 tot += _sample(f, x, y) s += 1.0 return tot def radon(f, delta): Pid = [[line_integral(f, angles[i], det[j]) for j in range(L)] for i in range(M)] meanv = sum(sum(r) for r in Pid) / (M * L) sig = 0.012 * meanv # 加入约 1.2% 量级的探测噪声 Pobs = [[_sample1d(Pid[i], det[j] - delta) + rnd.gauss(0, sig) for j in range(L)] for i in range(M)] return Pobs, Pid def backproject(P, delta_est): f = [[0.0] * N for _ in range(N)] for i in range(M): e0 = _m.cos(angles[i]); e1 = _m.sin(angles[i]) for y in range(N): for x in range(N): t = (x - cx) * e0 + (y - cy) * e1 + delta_est f[y][x] += _sample1d(P[i], t) norm = max(max(row) for row in f) if max(max(row) for row in f) > 0 else 1.0 return [[v / (M) for v in row] for row in f] # 除以角度数作强度归一 def _dft(x): Ln = len(x); return [sum(x[n] * _cm.exp(-2j * _m.pi * k * n / Ln) for n in range(Ln)) for k in range(Ln)] def _idft(X): Ln = len(X); return [sum(X[n] * _cm.exp(2j * _m.pi * k * n / Ln) for n in range(Ln)) / Ln for k in range(Ln)] def _ramp(row): X = _dft(row); Ln = len(X) for k in range(Ln): fk = k if k <= Ln // 2 else k - Ln X[k] *= abs(fk) * 2.0 / Ln return [v.real for v in _idft(X)] def fbp(P, delta_est): Q = [_ramp(P[i]) for i in range(M)] return backproject(Q, delta_est) def art(P, delta_est, K=12, lam=0.25): f = [[0.0] * N for _ in range(N)] for _ in range(K): for i in range(M): e0 = _m.cos(angles[i]); e1 = _m.sin(angles[i]) n0 = -e1; n1 = e0 for j in range(L): t = det[j] - delta_est # 前向估计 Rr = N * 0.98; s = -Rr; est = 0.0; cnt = 0; px = []; py = [] while s <= Rr: x = cx + t * e0 + s * n0; y = cy + t * e1 + s * n1 if 0 <= x < N - 1 and 0 <= y < N - 1: est += _sample(f, x, y) s += 1.0 # 误差沿射线分配（最近像素） s = -Rr while s <= Rr: x = cx + t * e0 + s * n0; y = cy + t * e1 + s * n1 ix = int(round(x)); iy = int(round(y)) if 0 <= ix < N and 0 <= iy < N: f[iy][ix] += lam * (P[i][j] - est) / Rr s += 1.0 return f def rmse(a, b): s = 0.0; n = 0 for y in range(N): for x in range(N): s += (a[y][x] - b[y][x]) ** 2; n += 1 return _m.sqrt(s / n) # ---- 生成模板与未知件的投影 ---- f_temp = build_phantom("template") f_unk = build_phantom("unknown") Pobs_t, Pid_t = radon(f_temp, delta_true) Pobs_u, Pid_u = radon(f_unk, delta_true) # ---- 三种 δ 标定（精度层级不同，主要区分在重建算法） ---- def shift_cost(dd, Pobs, Pid): s = 0.0 for i in range(M): row = Pobs[i] for j in range(L): v = _sample1d(row, det[j] + dd) d = v - Pid[i][j] s += d * d return s # 方法1：模板监督互相关（整数步长）→ 粗标 δ（论文1：几何法） d1 = min([float(x) for x in range(-8, 9)], key=lambda dd: shift_cost(dd, Pobs_t, Pid_t)) # 方法2：模板监督互相关（半整数步长）→ 较精 δ（论文2：频域法） d2 = min([round(-8.0 + x * 0.5, 1) for x in range(33)], key=lambda dd: shift_cost(dd, Pobs_t, Pid_t)) # 方法3：ART 残差精细优化（d2 附近 0.1 步长，监督）→ 最精 δ（论文3：迭代法） best3 = None; d3 = d2 lo = max(-8.0, round(d2 - 0.4, 1)); hi = min(8.0, round(d2 + 0.4, 1)) for dd in [round(lo + x * 0.1, 1) for x in range(int(round((hi - lo) / 0.1)) + 1)]: frec = art(Pobs_t, dd, K=10, lam=0.25) e = rmse(frec, f_temp) if best3 is None or e < best3: best3 = e; d3 = dd def alpha_of(frec, ftrue): num = 0.0; den = 0.0 for y in range(N): for x in range(N): num += frec[y][x] * ftrue[y][x]; den += ftrue[y][x] * ftrue[y][x] return num / den if den > 0 else 1.0 def apply_alpha(frec, a): return [[v * a for v in row] for row in frec] # 各法对模板重建 → 求幅度 α → 用于未知件校准 a1 = alpha_of(backproject(Pobs_t, d1), f_temp) a2 = alpha_of(fbp(Pobs_t, d2), f_temp) a3 = alpha_of(art(Pobs_t, d3, K=12, lam=0.25), f_temp) # ---- 用各自 δ 重建未知件（幅度校准到真值量级） ---- frec1 = apply_alpha(backproject(Pobs_u, d1), a1) frec2 = apply_alpha(fbp(Pobs_u, d2), a2) frec3 = apply_alpha(art(Pobs_u, d3, K=12, lam=0.25), a3) rmse1 = rmse(frec1, f_unk) rmse2 = rmse(frec2, f_unk) rmse3 = rmse(frec3, f_unk) # 模板重建误差（标定质量，校准后） rmse1_t = rmse(apply_alpha(backproject(Pobs_t, d1), a1), f_temp) rmse2_t = rmse(apply_alpha(fbp(Pobs_t, d2), a2), f_temp) rmse3_t = rmse(apply_alpha(art(Pobs_t, d3, K=12, lam=0.25), a3), f_temp) # 校验点：未知件某射线投影值（θ=0, t=0 与 t=10） check_t0 = Pobs_u[0][L // 2] check_t10 = Pobs_u[0][L // 2 + 10] # sinogram 维数 dim = (M, L) res = dict(seed=seed, N=N, M=M, L=L, delta_true=delta_true, dim=dim, d1=d1, d2=round(d2, 3), d3=d3, rmse1=rmse1, rmse2=rmse2, rmse3=rmse3, rmse1_t=rmse1_t, rmse2_t=rmse2_t, rmse3_t=rmse3_t, check_t0=check_t0, check_t10=check_t10, f_temp=f_temp, f_unk=f_unk, Pobs_u=Pobs_u, Pid_u=Pid_u, Pobs_t=Pobs_t, Pid_t=Pid_t, frec1=frec1, frec2=frec2, frec3=frec3, angles=angles, det=det) # 写 CSV（params 型：sinogram 平铺 + 元信息） with open(os.path.join(OUT, out_csv), "w", newline="") as f: w = csv.writer(f) w.writerow(["M", "L", "N", "delta_true"]) w.writerow([M, L, N, delta_true]) w.writerow(["sinogram_row", "angle_idx", "det_idx", "value"]) for i in range(M): for j in range(L): w.writerow([i * L + j, i, j, round(Pobs_u[i][j], 4)]) return res def gen_2017b(seed=2017, n_tasks=820, n_members=800, region=30.0, R=2.5, price_base=80.0, out_csv="cumcm2017b.csv"): """2017B「拍照赚钱」任务定价：合成会员/任务空间数据 + 接单率生成模型。 会员聚集于若干商业热区，任务均匀散布；原始定价近似常数(≈80元，题目缺陷)， 接单率由逻辑回归潜变量决定（价格↑、周边会员数↑ → 接单率↑）。 返回权威数字字典，并写出 tasks CSV（含可达性特征与完成情况），供三篇范文/出图共用。""" import math as _m rnd = random.Random(seed) # ---- 会员分布：若干商业热区 + 均匀背景 ---- K_hubs = 5 hubs = [(rnd.uniform(region * 0.15, region * 0.85), rnd.uniform(region * 0.15, region * 0.85)) for _ in range(K_hubs)] members = [] for _ in range(n_members): if rnd.random() < 0.65: hx, hy = rnd.choice(hubs) x = clip(hx + rnd.gauss(0, 1.6), 0, region) y = clip(hy + rnd.gauss(0, 1.6), 0, region) else: x = rnd.uniform(0, region); y = rnd.uniform(0, region) members.append((x, y)) # ---- 任务分布（均匀散布，模拟待拍照 POI） ---- tasks = [] for i in range(n_tasks): x = rnd.uniform(0, region); y = rnd.uniform(0, region) orig_price = price_base + rnd.gauss(0, 3.0) # 原始定价近似常数（题目缺陷） tasks.append({"id": i, "x": x, "y": y, "orig_price": orig_price}) # ---- 可达性特征：最近会员距离 d_n、半径 R 内会员数 cnt_R ---- R2 = R * R for t in tasks: dmin = 1e9; cnt = 0 for (mx, my) in members: dx = t["x"] - mx; dy = t["y"] - my; d2 = dx * dx + dy * dy if d2 < dmin: dmin = d2 if d2 <= R2: cnt += 1 t["d_nearest"] = _m.sqrt(dmin); t["cnt_R"] = cnt # ---- 接单率生成模型（逻辑回归潜变量） ---- a0, a1, a2 = -1.0, 1.3, 0.14 def sigmoid(z): return 1.0 / (1.0 + _m.exp(-z)) for t in tasks: pnorm = (t["orig_price"] - price_base) / 10.0 z = a0 + a1 * pnorm + a2 * t["cnt_R"] t["p_complete"] = sigmoid(z) t["completed"] = 1 if rnd.random() < t["p_complete"] else 0 completion_overall = sum(t["completed"] for t in tasks) / n_tasks price_mean = sum(t["orig_price"] for t in tasks) / n_tasks price_std = _m.sqrt(sum((t["orig_price"] - price_mean) ** 2 for t in tasks) / n_tasks) def _mean(v): return sum(v) / len(v) def _std(v): return _m.sqrt(sum((x - _mean(v)) ** 2 for x in v) / len(v)) or 1.0 # ---- 逻辑回归拟合（梯度下降，复现定价规律） ---- X1 = [(t["orig_price"] - price_base) / 10.0 for t in tasks] X2 = [float(t["cnt_R"]) for t in tasks] y = [t["completed"] for t in tasks] mx1, sx1 = _mean(X1), _std(X1); mx2, sx2 = _mean(X2), _std(X2) Z1 = [(v - mx1) / sx1 for v in X1]; Z2 = [(v - mx2) / sx2 for v in X2] w = [0.0, 0.0, 0.0]; lr = 0.3 for _ in range(4000): gw0 = gw1 = gw2 = 0.0 for i in range(n_tasks): pr = sigmoid(w[0] + w[1] * Z1[i] + w[2] * Z2[i]) e = pr - y[i] gw0 += e; gw1 += e * Z1[i]; gw2 += e * Z2[i] w[0] -= lr * gw0 / n_tasks; w[1] -= lr * gw1 / n_tasks; w[2] -= lr * gw2 / n_tasks w0, w1, w2 = w # 还原到原始量纲系数（便于解读）：logit = b0 + b1·pnorm + b2·cnt_R b1 = w1 / sx1; b2 = w2 / sx2 b0 = w0 - w1 * mx1 / sx1 - w2 * mx2 / sx2 fit_orig = (b0, b1, b2) def fit_p(price, cnt): u = (price - price_base) / 10.0 return sigmoid(w0 + w1 * (u - mx1) / sx1 + w2 * (cnt - mx2) / sx2) def logit_inv(p): return _m.log(p / (1.0 - p)) # ---- 诊断分箱 ---- def bucket_rate(key_fn, edges): out = [] for b in range(len(edges) - 1): lo, hi = edges[b], edges[b + 1] sel = [t for t in tasks if lo <= key_fn(t) < hi] if sel: out.append((round((lo + hi) / 2.0, 2), round(sum(t["completed"] for t in sel) / len(sel), 4), len(sel))) return out dist_edges = [0, 2, 4, 6, 8, 10, 12, 1e9] price_edges = [70, 74, 78, 82, 86, 90, 1e9] acc_edges = [0, 2, 5, 10, 20, 40, 1e9] by_acc = bucket_rate(lambda t: t["cnt_R"], acc_edges) by_price = bucket_rate(lambda t: t["orig_price"], price_edges) # 可达性归因（按周边会员数分箱的已完成/未完成计数，用于堆叠图） by_acc_counts = [] for b in range(len(acc_edges) - 1): lo, hi = acc_edges[b], acc_edges[b + 1] sel = [t for t in tasks if lo <= t["cnt_R"] < hi] if sel: done = sum(t["completed"] for t in sel) lab = ("[%d,%d)" % (lo, hi)) if hi < 1e8 else ("[%d,∞)" % lo) by_acc_counts.append((lab, done, len(sel) - done)) # ---- Q2：成本-覆盖多目标重定价 ---- T = 0.9 p_req = [] for t in tasks: u = mx1 + sx1 * (logit_inv(T) - w0 - w2 * (t["cnt_R"] - mx2) / sx2) / w1 p_req.append(max(price_base + 10.0 * u, 40.0)) orig_total = sum(t["orig_price"] for t in tasks) new_total = sum(p_req) uniform_baseline = max(p_req) * n_tasks achieved = sum(fit_p(p_req[i], tasks[i]["cnt_R"]) for i in range(n_tasks)) / n_tasks pareto = [] for TT in [0.5, 0.6, 0.7, 0.8, 0.85, 0.9, 0.95, 0.98]: pr = [max(price_base + 10.0 * (mx1 + sx1 * (logit_inv(TT) - w0 - w2 * (tasks[i]["cnt_R"] - mx2) / sx2) / w1), 40.0) for i in range(n_tasks)] cost = sum(pr) cov = sum(fit_p(pr[i], tasks[i]["cnt_R"]) for i in range(n_tasks)) / n_tasks pareto.append((TT, round(cost, 1), round(cov, 4))) by_acc_after = [] for b in range(len(acc_edges) - 1): lo, hi = acc_edges[b], acc_edges[b + 1] sel = [i for i in range(n_tasks) if lo <= tasks[i]["cnt_R"] < hi] if sel: cov = sum(fit_p(p_req[i], tasks[i]["cnt_R"]) for i in sel) / len(sel) by_acc_after.append((round((lo + hi) / 2.0, 2), round(cov, 4), len(sel))) # ---- Q3：空间聚类打包定价 ---- def kmeans(pts, k, iters=40): centers = [list(pts[rnd.randint(0, len(pts) - 1)]) for _ in range(k)] assign = [0] * len(pts) for _ in range(iters): for i, (x, y) in enumerate(pts): best, bd = 0, 1e18 for c, (cx, cy) in enumerate(centers): d = (x - cx) ** 2 + (y - cy) ** 2 if d < bd: bd = d; best = c assign[i] = best for c in range(k): xs = [pts[i][0] for i in range(len(pts)) if assign[i] == c] ys = [pts[i][1] for i in range(len(pts)) if assign[i] == c] if xs: centers[c] = [sum(xs) / len(xs), sum(ys) / len(ys)] return centers, assign k = 120 pts = [(t["x"], t["y"]) for t in tasks] centers, assign = kmeans(pts, k) clusters = {} for i, a in enumerate(assign): clusters.setdefault(a, []).append(i) cluster_of = assign[:] discount = 0.25 bundle_price = {}; cluster_sizes = {}; center_d = {}; per_task_eff = [0.0] * n_tasks for c in range(k): idx = clusters.get(c, []) if not idx: continue prs = [p_req[i] for i in idx] bp = max(prs) + discount * (sum(prs) - max(prs)) bundle_price[c] = bp; cluster_sizes[c] = len(idx) cx, cy = centers[c] dm = min((cx - mxm) ** 2 + (cy - mym) ** 2 for (mxm, mym) in members) center_d[c] = _m.sqrt(dm) for i in idx: per_task_eff[i] = bp / len(idx) bundle_total = sum(bundle_price.values()) bundle_saving = new_total - bundle_total avg_cluster_size = _mean([cluster_sizes[c] for c in cluster_sizes]) frac_small = sum(1 for c in cluster_sizes if cluster_sizes[c] <= 8) / len(cluster_sizes) # 雷达对比（Q2 独立定价 vs Q3 聚类打包，higher=better） def _corr(a, b): ma, mb = _mean(a), _mean(b) cov = sum((a[i] - ma) * (b[i] - mb) for i in range(len(a))) sa = _m.sqrt(sum((x - ma) ** 2 for x in a)); sb = _m.sqrt(sum((x - mb) ** 2 for x in b)) return cov / (sa * sb) if sa * sb > 0 else 0.0 cost_save_q2 = (uniform_baseline - new_total) / uniform_baseline cost_save_q3 = (uniform_baseline - bundle_total) / uniform_baseline avgp_q2 = new_total / n_tasks; avgp_q3 = bundle_total / n_tasks lo, hi = min(avgp_q2, avgp_q3), max(avgp_q2, avgp_q3) afford_q2 = 1 - (avgp_q2 - lo) / (hi - lo) if hi > lo else 1.0 afford_q3 = 1 - (avgp_q3 - lo) / (hi - lo) if hi > lo else 1.0 match_q2 = abs(_corr([p_req[i] for i in range(n_tasks)], [tasks[i]["cnt_R"] for i in range(n_tasks)])) match_q3 = abs(_corr(per_task_eff, [tasks[i]["cnt_R"] for i in range(n_tasks)])) spatial_q2 = 0.0 spatial_q3 = clip(avg_cluster_size / 8.0, 0.0, 1.0) radar_labels = ["成本节约", "覆盖保障", "单均可负担", "难度-价格匹配", "空间协同"] radar_q2 = [round(cost_save_q2, 3), round(achieved, 3), round(afford_q2, 3), round(match_q2, 3), round(spatial_q2, 3)] radar_q3 = [round(cost_save_q3, 3), round(achieved, 3), round(afford_q3, 3), round(match_q3, 3), round(spatial_q3, 3)] # ---- 网络采样（图7：热区局部会员-任务近邻） ---- hx, hy = hubs[0] bx0, bx1, by0, by1 = hx - 3, hx + 3, hy - 3, hy + 3 net_tasks = [(t["x"], t["y"]) for t in tasks if bx0 <= t["x"] <= bx1 and by0 <= t["y"] <= by1][:40] net_members = [(mx, my) for (mx, my) in members if bx0 - 1 <= mx <= bx1 + 1 and by0 - 1 <= my <= by1 + 1][:120] net_edges = [] for (tx, ty) in net_tasks: best, bd = None, 1e18 for (mx, my) in net_members: d = (tx - mx) ** 2 + (ty - my) ** 2 if d < bd: bd = d; best = (mx, my) if best: net_edges.append(((tx, ty), best)) # ---- 写出 CSV（含可达性特征与完成情况） ---- with open(os.path.join(OUT, out_csv), "w", encoding="utf-8-sig", newline="") as f: w = csv.writer(f) w.writerow(["task_id", "x", "y", "orig_price", "d_nearest", "cnt_R", "p_complete", "completed"]) for t in tasks: w.writerow([t["id"], round(t["x"], 3), round(t["y"], 3), round(t["orig_price"], 2), round(t["d_nearest"], 3), t["cnt_R"], round(t["p_complete"], 4), t["completed"]]) print("cumcm2017b.csv 已写出：%d 任务 × 8 字段；会员 %d；完成率 %.3f" % (n_tasks, n_members, completion_overall)) return dict(seed=seed, n_tasks=n_tasks, n_members=n_members, region=region, R=R, price_base=price_base, members=members, tasks=tasks, price_mean=price_mean, price_std=price_std, completion_overall=completion_overall, gen_a=(a0, a1, a2), fit_w=(w0, w1, w2), fit_orig=fit_orig, fit_std=(mx1, sx1, mx2, sx2), by_acc=by_acc, by_acc_counts=by_acc_counts, acc_labels=["0-2", "2-5", "5-10", "10-20", "20-40", "40+"], by_price=by_price, T=T, p_req=p_req, orig_total=orig_total, new_total=new_total, uniform_baseline=uniform_baseline, achieved=achieved, pareto=pareto, by_acc_after=by_acc_after, k=k, centers=centers, cluster_of=cluster_of, cluster_sizes=cluster_sizes, center_d=center_d, discount=discount, bundle_price=bundle_price, bundle_total=bundle_total, bundle_saving=bundle_saving, standalone_total=new_total, avg_cluster_size=avg_cluster_size, frac_small=frac_small, per_task_eff=per_task_eff, radar_labels=radar_labels, radar_q2=radar_q2, radar_q3=radar_q3, hubs=hubs, net_tasks=net_tasks, net_members=net_members, net_edges=net_edges) def gen_2017c(seed=2017, n=200, n_feat=6, anomaly_frac=0.06, out_csv="cumcm2017c.csv"): """2017C「颜色与物质判别」：合成 6 维颜色特征 + 二分类物质标签 + 注入离群异常。 实现 PCA(幂迭代)、LDA(闭式)、逻辑回归、KNN、高斯朴素贝叶斯、CART 决策树、 5 折分层交叉验证、ROC/AUC、混淆矩阵、特征重要性、马氏距离/kNN 异常检测、噪声鲁棒性。 返回权威数字字典并写出 CSV，供三篇范文/出图共用（图/正文/附录/工具四者一致）。""" import math as _m rnd = random.Random(seed) D = n_feat # ---- 生成两类物质：均值沿随机方向分离、共享协方差 → 部分重叠 ---- base = [rnd.gauss(0.0, 1.0) for _ in range(D)] direc = [rnd.uniform(-1.0, 1.0) for _ in range(D)] sep = 1.7 mu0 = [base[i] - 0.5 * sep * direc[i] for i in range(D)] mu1 = [base[i] + 0.5 * sep * direc[i] for i in range(D)] sigma = [rnd.uniform(0.75, 1.15) for _ in range(D)] def gauss_feat(mu): return [mu[i] + rnd.gauss(0.0, sigma[i]) for i in range(D)] # ---- 干净样本 + 注入异常（远离两类中心的离群点） ---- n_anom = int(round(n * anomaly_frac)) n_clean = n - n_anom X, y, is_anom = [], [], [] for _ in range(n_clean): lab = 0 if rnd.random() < 0.5 else 1 mu = mu0 if lab == 0 else mu1 X.append(gauss_feat(mu)); y.append(lab); is_anom.append(0) for _ in range(n_anom): u = [_m.sqrt(-2 * _m.log(rnd.random() + 1e-12)) * _m.cos(2 * _m.pi * rnd.random()) for _ in range(D)] norm = _m.sqrt(sum(v * v for v in u)) or 1.0 u = [v / norm for v in u] r = rnd.uniform(3.2, 4.2) pt = [base[i] + r * u[i] for i in range(D)] d0 = sum((pt[i] - mu0[i]) ** 2 for i in range(D)) d1 = sum((pt[i] - mu1[i]) ** 2 for i in range(D)) lab = 0 if d0 < d1 else 1 X.append(pt); y.append(lab); is_anom.append(1) cls_counts = {0: sum(1 for v in y if v == 0), 1: sum(1 for v in y if v == 1)} # ---- 工具：向量/矩阵 ---- def mean_vec(rows): return [sum(r[i] for r in rows) / len(rows) for i in range(D)] def std_vec(rows, mu): return [_m.sqrt(sum((r[i] - mu[i]) ** 2 for r in rows) / len(rows)) or 1e-9 for i in range(D)] def dot(a, b): return sum(a[i] * b[i] for i in range(len(a))) def mat_vec(M, v): return [dot(M[i], v) for i in range(len(M))] def mat_inv(M): k = len(M) A = [M[i][:] + [1.0 if i == j else 0.0 for j in range(k)] for i, row in enumerate(M)] for col in range(k): piv = max(range(col, k), key=lambda r: abs(A[r][col])) A[col], A[piv] = A[piv], A[col] pv = A[col][col] A[col] = [x / pv for x in A[col]] for r in range(k): if r != col: f = A[r][col] A[r] = [A[r][c] - f * A[col][c] for c in range(2 * k)] return [row[k:] for row in A] def top_eigen(C, iters=400): v = [1.0 / _m.sqrt(D)] * D for _ in range(iters): w = mat_vec(C, v) nv = _m.sqrt(dot(w, w)) or 1e-12 v = [x / nv for x in w] w = mat_vec(C, v) lam = dot(v, w) return lam, v def cov_mat(rows): m = mean_vec(rows); k = len(rows) C = [[0.0] * D for _ in range(D)] for r in rows: for i in range(D): for j in range(D): C[i][j] += (r[i] - m[i]) * (r[j] - m[j]) return [[C[i][j] / k for j in range(D)] for i in range(D)] # ---- 全局描述性统计 ---- gmean = mean_vec(X); gstd = std_vec(X, gmean) mu0_raw = mean_vec([X[i] for i in range(n) if y[i] == 0]) mu1_raw = mean_vec([X[i] for i in range(n) if y[i] == 1]) # ---- 70/30 分层训练/测试划分 ---- idx0 = [i for i in range(n) if y[i] == 0] idx1 = [i for i in range(n) if y[i] == 1] rnd.shuffle(idx0); rnd.shuffle(idx1) ntr0, ntr1 = int(round(0.7 * len(idx0))), int(round(0.7 * len(idx1))) train_idx = idx0[:ntr0] + idx1[:ntr1] test_idx = idx0[ntr0:] + idx1[ntr1:] Xtr = [X[i] for i in train_idx]; ytr = [y[i] for i in train_idx] Xte = [X[i] for i in test_idx]; yte = [y[i] for i in test_idx] mu_tr = mean_vec(Xtr); sd_tr = std_vec(Xtr, mu_tr) def zscore(rows): return [[(r[i] - mu_tr[i]) / sd_tr[i] for i in range(D)] for r in rows] Xtrz, Xtez = zscore(Xtr), zscore(Xte) mu0z = mean_vec([Xtrz[i] for i in range(len(Xtrz)) if ytr[i] == 0]) mu1z = mean_vec([Xtrz[i] for i in range(len(Xtrz)) if ytr[i] == 1]) # ---- PCA（幂迭代取前 2 主成分 + 全部特征值） ---- C = cov_mat(Xtrz) traceC = sum(C[i][i] for i in range(D)) Ccur = [row[:] for row in C] eigvals = [] for _ in range(D): lam, v = top_eigen(Ccur) eigvals.append(lam) Ccur = [[Ccur[i][j] - lam * v[i] * v[j] for j in range(D)] for i in range(D)] eigvals.sort(reverse=True) pca_var_norm = [e / traceC for e in eigvals] lam1, v1 = top_eigen(C) C2 = [[C[i][j] - lam1 * v1[i] * v1[j] for j in range(D)] for i in range(D)] lam2, v2 = top_eigen(C2) pca_scores_test = [(dot(v1, Xtez[i]), dot(v2, Xtez[i])) for i in range(len(Xtez))] pca_scores_train = [(dot(v1, Xtrz[i]), dot(v2, Xtrz[i])) for i in range(len(Xtrz))] # ---- LDA（Fisher 闭式） ---- def outer(a, b): return [[a[i] * b[j] for j in range(D)] for i in range(D)] Sw = [[0.0] * D for _ in range(D)] for lab, mu in ((0, mu0z), (1, mu1z)): rows = [Xtrz[i] for i in range(len(Xtrz)) if ytr[i] == lab] m = mean_vec(rows) for r in rows: d = [r[i] - m[i] for i in range(D)] for i in range(D): for j in range(D): Sw[i][j] += d[i] * d[j] diff = [mu1z[i] - mu0z[i] for i in range(D)] invSw = mat_inv(Sw) wlda = mat_vec(invSw, diff) wnorm = _m.sqrt(dot(wlda, wlda)) or 1.0 wlda = [x / wnorm for x in wlda] p0 = dot(wlda, mu0z); p1 = dot(wlda, mu1z) lda_thr = (p0 + p1) / 2.0 Sb = outer(diff, diff) fisher_J = dot(wlda, mat_vec(Sb, wlda)) / (dot(wlda, mat_vec(Sw, wlda)) or 1e-12) lda_scores = [dot(wlda, Xtez[i]) for i in range(len(Xtez))] lda_pred = [1 if dot(wlda, Xtez[i]) > lda_thr else 0 for i in range(len(Xtez))] # ---- 度量函数 ---- def conf_mat(yt, yp): tn = sum(1 for i in range(len(yt)) if yt[i] == 0 and yp[i] == 0) fp = sum(1 for i in range(len(yt)) if yt[i] == 0 and yp[i] == 1) fn = sum(1 for i in range(len(yt)) if yt[i] == 1 and yp[i] == 0) tp = sum(1 for i in range(len(yt)) if yt[i] == 1 and yp[i] == 1) return [[tn, fp], [fn, tp]] def acc(yt, yp): return sum(1 for i in range(len(yt)) if yt[i] == yp[i]) / len(yt) def prec(yt, yp): tp = sum(1 for i in range(len(yt)) if yt[i] == 1 and yp[i] == 1) fp = sum(1 for i in range(len(yt)) if yt[i] == 0 and yp[i] == 1) return tp / (tp + fp) if tp + fp else 0.0 def rec(yt, yp): tp = sum(1 for i in range(len(yt)) if yt[i] == 1 and yp[i] == 1) fn = sum(1 for i in range(len(yt)) if yt[i] == 1 and yp[i] == 0) return tp / (tp + fn) if tp + fn else 0.0 def f1(yt, yp): p = prec(yt, yp); r = rec(yt, yp) return 2 * p * r / (p + r) if p + r else 0.0 def roc_auc(scores, yt): order = sorted(range(len(scores)), key=lambda i: -scores[i]) P = sum(yt); Nn = len(yt) - P tp = fp = 0.0 pts = [] for i in order: if yt[i] == 1: tp += 1 else: fp += 1 pts.append((fp / Nn if Nn else 0.0, tp / P if P else 0.0)) area = 0.0 for a, b in zip(pts[:-1], pts[1:]): area += (b[0] - a[0]) * (a[1] + b[1]) / 2.0 return area, pts # ---- 逻辑回归（梯度下降） ---- def train_lr(Xz, yz, lr=0.3, iters=4000): m = len(Xz); w = [0.0] * (D + 1) for _ in range(iters): g0 = 0.0; g1 = [0.0] * D for i in range(m): z = w[0] + sum(w[1 + k] * Xz[i][k] for k in range(D)) pr = 1.0 / (1.0 + _m.exp(-z)) e = pr - yz[i] g0 += e for k in range(D): g1[k] += e * Xz[i][k] w[0] -= lr * g0 / m for k in range(D): w[1 + k] -= lr * g1[k] / m return w def pred_lr(w, Xz): return [1.0 / (1.0 + _m.exp(-(w[0] + sum(w[1 + k] * x[k] for k in range(D))))) for x in Xz] wlr = train_lr(Xtrz, ytr) plr = pred_lr(wlr, Xtez) lr_pred = [1 if p >= 0.5 else 0 for p in plr] lr_auc, lr_roc = roc_auc(plr, yte) lr_coef = wlr[:] # ---- KNN ---- def pred_knn(Xa, ya, Xe, k=5): out = [] for x in Xe: ds = sorted(range(len(Xa)), key=lambda i: sum((x[t] - Xa[i][t]) ** 2 for t in range(D))) nn = ds[:k] c1 = sum(1 for i in nn if ya[i] == 1) out.append(1 if c1 >= (k - c1) else 0) return out def knn_scores(Xa, ya, Xe, k=5): out = [] for x in Xe: ds = sorted(range(len(Xa)), key=lambda i: sum((x[t] - Xa[i][t]) ** 2 for t in range(D))) nn = ds[:k] c1 = sum(1 for i in nn if ya[i] == 1) out.append(c1 / k) return out knn_pred = pred_knn(Xtrz, ytr, Xtez, 5) knn_sc = knn_scores(Xtrz, ytr, Xtez, 5) knn_auc, knn_roc = roc_auc(knn_sc, yte) # ---- 高斯朴素贝叶斯 ---- def train_gnb(Xz, yz): rows0 = [Xz[i] for i in range(len(Xz)) if yz[i] == 0] rows1 = [Xz[i] for i in range(len(Xz)) if yz[i] == 1] m0 = mean_vec(rows0); s0 = std_vec(rows0, m0) m1 = mean_vec(rows1); s1 = std_vec(rows1, m1) return (m0, s0, m1, s1, len(rows0) / len(Xz), len(rows1) / len(Xz)) def pred_gnb(model, Xz): m0, s0, m1, s1, p0, p1 = model out = [] for x in Xz: def lg(m, s): return sum(-_m.log(s[i] + 1e-12) - 0.5 * ((x[i] - m[i]) / s[i]) ** 2 for i in range(D)) a0 = lg(m0, s0) + _m.log(p0 + 1e-12) a1 = lg(m1, s1) + _m.log(p1 + 1e-12) out.append(1 if a1 > a0 else 0) return out def gnb_scores(model, Xz): m0, s0, m1, s1, p0, p1 = model out = [] for x in Xz: def lg(m, s): return sum(-_m.log(s[i] + 1e-12) - 0.5 * ((x[i] - m[i]) / s[i]) ** 2 for i in range(D)) a0 = lg(m0, s0) + _m.log(p0 + 1e-12) a1 = lg(m1, s1) + _m.log(p1 + 1e-12) mx = max(a0, a1) out.append(_m.exp(a1 - mx) / (_m.exp(a0 - mx) + _m.exp(a1 - mx))) return out gnb_model = train_gnb(Xtrz, ytr) gnb_pred = pred_gnb(gnb_model, Xtez) gnb_sc = gnb_scores(gnb_model, Xtez) gnb_auc, gnb_roc = roc_auc(gnb_sc, yte) # ---- CART 决策树（深度≤3） ---- def gini(yz): if not yz: return 0.0 p1 = sum(yz) / len(yz); p0 = 1 - p1 return 1 - p1 * p1 - p0 * p0 def train_tree(Xz, yz, depth=0, max_depth=3): node = {"leaf": False, "pred": 1 if sum(yz) > len(yz) / 2 else 0, "n": len(yz)} if depth >= max_depth or len(set(yz)) <= 1 or len(yz) < 4: node["leaf"] = True; return node best = None; best_gain = 0.0 parent_g = gini(yz) for f in range(D): vals = sorted(set(Xz[i][f] for i in range(len(Xz)))) for t in range(len(vals) - 1): thr = (vals[t] + vals[t + 1]) / 2.0 ly = [yz[i] for i in range(len(Xz)) if Xz[i][f] <= thr] ry = [yz[i] for i in range(len(Xz)) if Xz[i][f] > thr] if not ly or not ry: continue wg = (len(ly) * gini(ly) + len(ry) * gini(ry)) / len(yz) gain = parent_g - wg if gain > best_gain: best_gain = gain; best = (f, thr, ly, ry) if best is None: node["leaf"] = True; return node f, thr, ly, ry = best node["feat"] = f; node["thr"] = thr; node["gain"] = best_gain node["left"] = train_tree([Xz[i] for i in range(len(Xz)) if Xz[i][f] <= thr], ly, depth + 1, max_depth) node["right"] = train_tree([Xz[i] for i in range(len(Xz)) if Xz[i][f] > thr], ry, depth + 1, max_depth) return node def pred_tree(node, x): while not node["leaf"]: node = node["left"] if x[node["feat"]] <= node["thr"] else node["right"] return node["pred"] tree = train_tree(Xtrz, ytr) tree_pred = [pred_tree(tree, x) for x in Xtez] tree_auc, _ = roc_auc([float(p) for p in tree_pred], yte) feat_imp = [0.0] * D def collect(node, n_total): if node.get("leaf"): return if "gain" in node: feat_imp[node["feat"]] += node["gain"] * (node["n"] / n_total) collect(node["left"], n_total); collect(node["right"], n_total) collect(tree, len(ytr)) # ---- 5 折分层交叉验证 ---- def stratified_folds(yz, k=5): i0 = [i for i in range(len(yz)) if yz[i] == 0] i1 = [i for i in range(len(yz)) if yz[i] == 1] rnd.shuffle(i0); rnd.shuffle(i1) folds = [[] for _ in range(k)] for i, ix in enumerate(i0): folds[i % k].append(ix) for i, ix in enumerate(i1): folds[i % k].append(ix) return folds def cv_run(): folds = stratified_folds(ytr, 5) res = {"LR": [], "LDA": [], "KNN": [], "GNB": [], "Tree": []} for fold in range(5): te = folds[fold] tr = [i for f in range(5) if f != fold for i in folds[f]] Xa = [Xtrz[i] for i in tr]; ya = [ytr[i] for i in tr] Xe = [Xtrz[i] for i in te]; ye = [ytr[i] for i in te] w = train_lr(Xa, ya); pa = [1 if p >= 0.5 else 0 for p in pred_lr(w, Xe)] res["LR"].append(acc(ye, pa)) m0 = mean_vec([Xa[i] for i in range(len(Xa)) if ya[i] == 0]) m1 = mean_vec([Xa[i] for i in range(len(Xa)) if ya[i] == 1]) Swc = [[0.0] * D for _ in range(D)] for lab, mu in ((0, m0), (1, m1)): rows = [Xa[i] for i in range(len(Xa)) if ya[i] == lab]; m = mean_vec(rows) for r in rows: d = [r[i] - m[i] for i in range(D)] for i in range(D): for j in range(D): Swc[i][j] += d[i] * d[j] dif = [m1[i] - m0[i] for i in range(D)]; iS = mat_inv(Swc) wld = mat_vec(iS, dif); wn = _m.sqrt(dot(wld, wld)) or 1; wld = [x / wn for x in wld] th = (dot(wld, m0) + dot(wld, m1)) / 2.0 pe = [1 if dot(wld, Xe[i]) > th else 0 for i in range(len(Xe))] res["LDA"].append(acc(ye, pe)) res["KNN"].append(acc(ye, pred_knn(Xa, ya, Xe, 5))) res["GNB"].append(acc(ye, pred_gnb(train_gnb(Xa, ya), Xe))) tr2 = train_tree(Xa, ya); res["Tree"].append(acc(ye, [pred_tree(tr2, Xe[i]) for i in range(len(Xe))])) return res cv = cv_run() cv_models = ["LR", "LDA", "KNN", "GNB", "Tree"] cv_mean = {m: sum(cv[m]) / 5 for m in cv_models} cv_std = {m: _m.sqrt(sum((v - cv_mean[m]) ** 2 for v in cv[m]) / 5) for m in cv_models} # ---- 各模型测试集指标 ---- def metrics(yt, yp, aucv): return dict(acc=round(acc(yt, yp), 4), prec=round(prec(yt, yp), 4), rec=round(rec(yt, yp), 4), f1=round(f1(yt, yp), 4), auc=round(aucv, 4)) model_metrics = { "LR": metrics(yte, lr_pred, lr_auc), "LDA": metrics(yte, lda_pred, roc_auc(lda_scores, yte)[0]), "KNN": metrics(yte, knn_pred, knn_auc), "GNB": metrics(yte, gnb_pred, gnb_auc), "Tree": metrics(yte, tree_pred, tree_auc), } lda_cm = conf_mat(yte, lda_pred) lr_cm = conf_mat(yte, lr_pred) knn_cm = conf_mat(yte, knn_pred) gnb_cm = conf_mat(yte, gnb_pred) tree_cm = conf_mat(yte, tree_pred) best_name = max(model_metrics, key=lambda m: model_metrics[m]["acc"]) best_cm = conf_mat(yte, {"LR": lr_pred, "LDA": lda_pred, "KNN": knn_pred, "GNB": gnb_pred, "Tree": tree_pred}[best_name]) # ---- 特征判别强度（Cohen's d） + 逻辑回归重要性 ---- def cstd(lab): rows = [Xtrz[i] for i in range(len(Xtrz)) if ytr[i] == lab] return std_vec(rows, mean_vec(rows)) s0c, s1c = cstd(0), cstd(1) pooled = [_m.sqrt((s0c[i] ** 2 + s1c[i] ** 2) / 2) for i in range(D)] feat_sep = [abs(mu1z[i] - mu0z[i]) / pooled[i] for i in range(D)] lr_importance = [abs(wlr[1 + i]) for i in range(D)] s = sum(lr_importance) or 1.0 lr_importance = [v / s for v in lr_importance] s2 = sum(feat_imp) or 1.0 tree_importance = [v / s2 for v in feat_imp] # ---- 学习曲线（LR vs KNN，训练比例 0.2→1.0） ---- lc_frac = [0.2, 0.4, 0.6, 0.8, 1.0] train_local = list(range(len(Xtrz))) lc_lr, lc_knn = [], [] for f in lc_frac: kk = max(2, int(round(f * len(train_local)))) Xa = [Xtrz[i] for i in train_local[:kk]]; ya = [ytr[i] for i in train_local[:kk]] w = train_lr(Xa, ya); pa = [1 if p >= 0.5 else 0 for p in pred_lr(w, Xtez)] lc_lr.append(round(acc(yte, pa), 4)) lc_knn.append(round(acc(yte, pred_knn(Xa, ya, Xtez, 5)), 4)) # ---- 异常检测（在训练集上筛查离群点） ---- def class_cov(lab): rows = [Xtrz[i] for i in range(len(Xtrz)) if ytr[i] == lab] m = mean_vec(rows); k = len(rows) Cc = [[0.0] * D for _ in range(D)] for r in rows: for i in range(D): for j in range(D): Cc[i][j] += (r[i] - m[i]) * (r[j] - m[j]) return [[Cc[i][j] / k for j in range(D)] for i in range(D)], m C0, m0c = class_cov(0); C1, m1c = class_cov(1) invC0 = mat_inv(C0); invC1 = mat_inv(C1) def maha(x, Cinv, mu): d = [x[i] - mu[i] for i in range(D)] return _m.sqrt(max(0.0, sum(d[i] * sum(Cinv[i][j] * d[j] for j in range(D)) for i in range(D)))) mah, knnd = [0.0] * len(Xtrz), [0.0] * len(Xtrz) for i in range(len(Xtrz)): Cinv, mu = (invC0, m0c) if ytr[i] == 0 else (invC1, m1c) mah[i] = maha(Xtrz[i], Cinv, mu) ds = sorted(sum((Xtrz[i][t] - Xtrz[j][t]) ** 2 for t in range(D)) for j in range(len(Xtrz)) if j != i) knnd[i] = _m.sqrt(sum(ds[:5]) / 5.0) def mm(a): lo, hi = min(a), max(a) return [(v - lo) / (hi - lo) if hi > lo else 0.0 for v in a] comb = [(mm(mah)[i] + mm(knnd)[i]) / 2.0 for i in range(len(Xtrz))] thr_comb = sorted(comb)[int(0.95 * (len(comb) - 1))] flagged = [i for i in range(len(Xtrz)) if comb[i] >= thr_comb] flagged_global = [train_idx[i] for i in flagged] true_anom = [train_idx[i] for i in range(len(Xtrz)) if is_anom[train_idx[i]] == 1] hit = [i for i in flagged_global if i in set(true_anom)] det_prec = len(hit) / len(flagged_global) if flagged_global else 0.0 det_rec = len(hit) / len(true_anom) if true_anom else 0.0 det_f1 = 2 * det_prec * det_rec / (det_prec + det_rec) if det_prec + det_rec else 0.0 def det_metrics(score_list): thr = sorted(score_list)[int(0.95 * (len(score_list) - 1))] fl = [i for i in range(len(score_list)) if score_list[i] >= thr] flg = set(train_idx[i] for i in fl) hi = len([i for i in flg if i in set(true_anom)]) p = hi / len(flg) if flg else 0.0 r = hi / len(true_anom) if true_anom else 0.0 f = 2 * p * r / (p + r) if p + r else 0.0 anom_s = [score_list[i] for i in range(len(score_list)) if is_anom[train_idx[i]] == 1] norm_s = [score_list[i] for i in range(len(score_list)) if is_anom[train_idx[i]] == 0] sep_sc = (sum(anom_s) / len(anom_s)) / (sum(norm_s) / len(norm_s)) if norm_s else 0.0 return p, r, f, sep_sc mah_det = det_metrics(mah) knn_det = det_metrics(knnd) # ---- 清洗训练集异常后重训（LR / KNN 测试准确率变化） ---- keep = [i for i in range(len(Xtrz)) if i not in set(flagged)] Xc = [Xtrz[i] for i in keep]; yc = [ytr[i] for i in keep] w_lr_c = train_lr(Xc, yc); plr_c = [1 if p >= 0.5 else 0 for p in pred_lr(w_lr_c, Xtez)] acc_lr_before = acc(yte, lr_pred); acc_lr_after = acc(yte, plr_c) pk_c = pred_knn(Xc, yc, Xtez, 5); acc_knn_before = acc(yte, knn_pred); acc_knn_after = acc(yte, pk_c) # ---- 噪声鲁棒性（特征叠加 0/5/10/20% 标准差高斯噪声） ---- noise_levels = [0.0, 0.05, 0.10, 0.20] rob_lr, rob_lda, rob_knn = [], [], [] for nl in noise_levels: rnr = random.Random(seed * 100 + int(nl * 1000)) def noisy(rows): return [[rows[i][t] + nl * sd_tr[t] * rnr.gauss(0, 1) for t in range(D)] for i in range(len(rows))] Xa = noisy(Xtrz); Xe = noisy(Xtez) w = train_lr(Xa, ytr); pa = [1 if p >= 0.5 else 0 for p in pred_lr(w, Xe)] rob_lr.append(round(acc(yte, pa), 4)) m0 = mean_vec([Xa[i] for i in range(len(Xa)) if ytr[i] == 0]) m1 = mean_vec([Xa[i] for i in range(len(Xa)) if ytr[i] == 1]) Swc = [[0.0] * D for _ in range(D)] for lab, mu in ((0, m0), (1, m1)): rows = [Xa[i] for i in range(len(Xa)) if ytr[i] == lab]; m = mean_vec(rows) for r in rows: d = [r[i] - m[i] for i in range(D)] for i in range(D): for j in range(D): Swc[i][j] += d[i] * d[j] dif = [m1[i] - m0[i] for i in range(D)]; iS = mat_inv(Swc) wld = mat_vec(iS, dif); wn = _m.sqrt(dot(wld, wld)) or 1; wld = [x / wn for x in wld] th = (dot(wld, m0) + dot(wld, m1)) / 2.0 pe = [1 if dot(wld, Xe[i]) > th else 0 for i in range(len(Xe))] rob_lda.append(round(acc(yte, pe), 4)) rob_knn.append(round(acc(yte, pred_knn(Xa, ytr, Xe, 5)), 4)) # ---- 相关性矩阵（全量原始特征） ---- def corr_matrix(rows): m = mean_vec(rows); k = len(rows) sd = [_m.sqrt(sum((r[i] - m[i]) ** 2 for r in rows) / k) or 1e-9 for i in range(D)] Cc = [[0.0] * D for _ in range(D)] for i in range(D): for j in range(i, D): c = sum((rows[r][i] - m[i]) * (rows[r][j] - m[j]) for r in range(k)) / (sd[i] * sd[j] * k) Cc[i][j] = c; Cc[j][i] = c return Cc corr = corr_matrix(X) # ---- 写出 CSV ---- with open(os.path.join(OUT, out_csv), "w", encoding="utf-8-sig", newline="") as f: w = csv.writer(f) w.writerow(["f1", "f2", "f3", "f4", "f5", "f6", "label"]) for i in range(n): w.writerow([round(X[i][j], 3) for j in range(D)] + [y[i]]) print("cumcm2017c.csv 已写出：%d 样本 × %d 特征；类别 %s；注入异常 %d" % ( n, D, cls_counts, n_anom)) res = dict( seed=seed, n=n, n_feat=D, n_anom=n_anom, anomaly_frac=anomaly_frac, header=["f1", "f2", "f3", "f4", "f5", "f6", "label"], X=X, y=y, is_anom=is_anom, cls_counts=cls_counts, gmean=gmean, gstd=gstd, mu0_raw=mu0_raw, mu1_raw=mu1_raw, train_idx=train_idx, test_idx=test_idx, n_train=len(train_idx), n_test=len(test_idx), mu0z=mu0z, mu1z=mu1z, sd_tr=sd_tr, pca_var_norm=pca_var_norm, pca_scores_test=pca_scores_test, pca_scores_train=pca_scores_train, pca_v1=v1, pca_v2=v2, wlda=wlda, lda_thr=lda_thr, fisher_J=fisher_J, lda_scores=lda_scores, lda_pred=lda_pred, lda_cm=lda_cm, lr_coef=lr_coef, lr_pred=lr_pred, lr_cm=lr_cm, lr_roc=lr_roc, knn_pred=knn_pred, knn_cm=knn_cm, knn_roc=knn_roc, gnb_pred=gnb_pred, gnb_cm=gnb_cm, gnb_roc=gnb_roc, tree_pred=tree_pred, tree_cm=tree_cm, tree_feat_split=None, model_metrics=model_metrics, best_name=best_name, best_cm=best_cm, cv_models=cv_models, cv_mean=cv_mean, cv_std=cv_std, feat_sep=feat_sep, lr_importance=lr_importance, tree_importance=tree_importance, lc_frac=lc_frac, lc_lr=lc_lr, lc_knn=lc_knn, mah=mah, knnd=knnd, comb=comb, thr_comb=thr_comb, flagged=flagged, flagged_global=flagged_global, true_anom=true_anom, hit=hit, det_prec=det_prec, det_rec=det_rec, det_f1=det_f1, mah_det=mah_det, knn_det=knn_det, acc_lr_before=acc_lr_before, acc_lr_after=acc_lr_after, acc_knn_before=acc_knn_before, acc_knn_after=acc_knn_after, noise_levels=noise_levels, rob_lr=rob_lr, rob_lda=rob_lda, rob_knn=rob_knn, corr=corr, ) return res # ===================== 2018A 高温作业专用服装（多层一维热传导 PDE + 厚度优化） ===================== def _thomas(L, D, U, b): """三对角矩阵求解（Thomas 算法），与 _simulate 的隐式格式配套。""" n = len(b) cp = [0.0] * n dp = [0.0] * n cp[0] = U[0] / D[0] dp[0] = b[0] / D[0] for i in range(1, n): m = D[i] - L[i] * cp[i - 1] cp[i] = U[i] / m dp[i] = (b[i] - L[i] * dp[i - 1]) / m x = [0.0] * n x[-1] = dp[-1] for i in range(n - 2, -1, -1): x[i] = dp[i] - cp[i] * x[i + 1] return x # --------------------------------------------------------------- # 2018A 高温作业专用服装：四层一维非稳态热传导 # 层序（由外向内）：I 阻燃外壳 → II 隔热层(设计变量) → III 舒适层 → IV 空气隙(设计变量) → 皮肤 # 说明：本数据为合成练习数据（量级取自公开教材与赛题常用值），非官方原题数据。 # --------------------------------------------------------------- MAT = dict( L1=0.6e-3, k1=0.082, r1=300.0, c1=1377.0, # 第I层 阻燃外壳（厚度固定） L2=6.0e-3, k2=0.370, r2=862.0, c2=2100.0, # 第II层 隔热层（设计变量，大热惯量） L3=3.6e-3, k3=0.045, r3=74.2, c3=1726.0, # 第III层 舒适层（厚度固定，低导热） L4=5.5e-3, k4=0.028, r4=1.18, c4=1005.0, # 第IV层 空气隙（设计变量） ) H_RAD = 6.5 # 空气隙线性化辐射换热系数 W/(m²·K)，由 4σεT̄³ 估得 L2_LO, L2_HI = 0.6, 25.0 # 第II层工艺厚度区间 mm L4_LO, L4_HI = 0.6, 6.4 # 第IV层空气隙工艺厚度区间 mm T_SAFE = 47.0 # 皮肤温度安全上限 °C T_WARN = 44.0 # 灼痛预警温度 °C T_WARN_MAX = 300.0 # 允许超过 44°C 的累计时长 s def k_air_eff(L4, mp=None): """空气隙有效导热系数：分子导热 + 线性化辐射（随间隙增厚而增大）。 R_air = L4/(k4 + h_rad*L4) 随 L4 单调增但上界为 1/h_rad —— 收益饱和。""" k4 = (mp or MAT)["k4"] return k4 + H_RAD * L4 def _build_mesh(L2, L4, mp, N): """层贴合（layer-conforming）网格：每层内部均匀剖分，层界面恒落在节点上。 这样每个控制体单元完全属于一种材料，导热率无需调和平均， 也消除了均匀网格因界面错位而产生的 ±0.1°C 非单调抖动。 返回 (节点坐标 x, 各段导热率 seg_k, 各段体积热容 seg_rc)。""" Ls = [mp["L1"], L2, mp["L3"], L4] ks = [mp["k1"], mp["k2"], mp["k3"], k_air_eff(L4, mp)] rcs = [mp["r1"] * mp["c1"], mp["r2"] * mp["c2"], mp["r3"] * mp["c3"], mp["r4"] * mp["c4"]] tot = sum(Ls) segs = [max(4, int(round((N - 1) * L / tot))) for L in Ls] x = [0.0] seg_k, seg_rc = [], [] for j, L in enumerate(Ls): d = L / segs[j] for _ in range(segs[j]): x.append(x[-1] + d) seg_k.append(ks[j]) seg_rc.append(rcs[j]) return x, seg_k, seg_rc def _mesh_x(L2, L4, mp, N): """仅取网格节点坐标（供绘制温度剖面用）。""" return _build_mesh(L2, L4, mp, N)[0] def _simulate(L2, L4, T_env, h_skin, T_core, T0, dt, T_end, N, mp, extra_x=None, rec_every=0): """四层一维热传导：层贴合有限体积 + 隐式后向欧拉（无条件稳定）。 外表面 x=0：Dirichlet T=T_env； 皮肤面 x=L：含半控制体热容的 Robin 边界 -k∂T/∂x = h_skin (T_skin - T_core)。 extra_x 为需额外记录的位置（单位 m），取最近节点。 返回 (skin[(t,T)...], extra[[(t,T)...]...], grid[(t,[T...])...])。""" xs, seg_k, seg_rc = _build_mesh(L2, L4, mp, N) n = len(xs) d = [xs[i + 1] - xs[i] for i in range(n - 1)] # 节点热容（相邻半单元之和）与界面导热系数 cap = [0.0] * n cap[0] = seg_rc[0] * d[0] / 2.0 cap[n - 1] = seg_rc[n - 2] * d[n - 2] / 2.0 for i in range(1, n - 1): cap[i] = (seg_rc[i - 1] * d[i - 1] + seg_rc[i] * d[i]) / 2.0 cond = [seg_k[i] / d[i] for i in range(n - 1)] extra_idx = [] if extra_x: for xv in extra_x: best_i, best_e = 0, None for i in range(n): e = abs(xs[i] - xv) if best_e is None or e < best_e: best_i, best_e = i, e extra_idx.append(best_i) steps = int(round(T_end / dt)) T = [T0] * n skin = [(0.0, T[-1])] extra = [[(0.0, T[idx])] for idx in extra_idx] if extra_idx else [] grid = [] for step in range(1, steps + 1): Lm = [0.0] * n Dm = [0.0] * n Um = [0.0] * n bm = [0.0] * n Dm[0] = 1.0 bm[0] = T_env for i in range(1, n - 1): aW = cond[i - 1] aE = cond[i] hi = cap[i] / dt Lm[i] = -aW Um[i] = -aE Dm[i] = hi + aW + aE bm[i] = hi * T[i] aW = cond[n - 2] hiN = cap[n - 1] / dt Lm[n - 1] = -aW Dm[n - 1] = hiN + aW + h_skin Um[n - 1] = 0.0 bm[n - 1] = hiN * T[n - 1] + h_skin * T_core T = _thomas(Lm, Dm, Um, bm) tcur = step * dt skin.append((tcur, T[-1])) if extra_idx: for kk, idx in enumerate(extra_idx): extra[kk].append((tcur, T[idx])) if rec_every and step % rec_every == 0: grid.append((tcur, list(T))) return skin, extra, grid def _validate_fd(dt=30.0, T_end=3600.0, L=0.20, N=200): """半无限固体验证：一维数值解 vs 解析 erfc 解，返回最大绝对误差 °C。""" k, r, c = 0.370, 862.0, 2100.0 alpha = k / (r * c) dx = L / (N - 1) dx2 = dx * dx T0, Tsurf = 37.0, 80.0 T = [T0] * N check = [(0.005, 600), (0.01, 600), (0.01, 1800), (0.02, 1800), (0.01, 3600), (0.03, 3600)] max_err = 0.0 steps = int(round(T_end / dt)) for step in range(1, steps + 1): Lm = [0.0] * N Dm = [0.0] * N Um = [0.0] * N bm = [0.0] * N Dm[0] = 1.0 bm[0] = Tsurf for i in range(1, N - 1): a = k / dx2 hi = r * c / dt Lm[i] = -a Um[i] = -a Dm[i] = hi + 2 * a bm[i] = hi * T[i] Dm[N - 1] = 1.0 bm[N - 1] = T0 T = _thomas(Lm, Dm, Um, bm) t = step * dt for (xx, tt) in check: if abs(t - tt) < dt / 2: ana = T0 + (Tsurf - T0) * math.erfc(xx / (2 * math.sqrt(alpha * t))) idx = int(round(xx / dx)) err = abs(T[idx] - ana) if err > max_err: max_err = err return max_err def _R_total(L2, L4, mp): """稳态串联热阻 m²K/W（L2、L4 单位 m）。""" return (mp["L1"] / mp["k1"] + L2 / mp["k2"] + mp["L3"] / mp["k3"] + L4 / k_air_eff(L4, mp)) def _C_total(L2, L4, mp): """单位面积热容 J/(m²·K)。""" return (mp["r1"] * mp["c1"] * mp["L1"] + mp["r2"] * mp["c2"] * L2 + mp["r3"] * mp["c3"] * mp["L3"] + mp["r4"] * mp["c4"] * L4) def _steady_skin(T_env, h_skin, T_core, L2, L4, mp): """稳态皮肤温度闭式：T_ss = (T_env + h R T_core)/(1 + h R)。""" R = _R_total(L2, L4, mp) return (T_env + h_skin * R * T_core) / (1.0 + h_skin * R) def _pearson(xs, ys): n = len(xs) mx = sum(xs) / n my = sum(ys) / n cov = sum((xs[i] - mx) * (ys[i] - my) for i in range(n)) sx = math.sqrt(sum((x - mx) ** 2 for x in xs)) sy = math.sqrt(sum((y - my) ** 2 for y in ys)) return cov / (sx * sy) if sx * sy > 0 else 0.0 def _cross_time(series, thr): """首次达到 thr 的时刻（线性插值）；未达到返回 None。""" for i in range(1, len(series)): t0, v0 = series[i - 1] t1, v1 = series[i] if v1 >= thr: if v1 == v0: return t1 return t0 + (thr - v0) / (v1 - v0) * (t1 - t0) return None def _time_above(series, thr, T_end): """温度超过 thr 的累计时长 s（单调升温假设下 = T_end - 首次穿越时刻）。""" tc = _cross_time(series, thr) if tc is None: return 0.0 return max(0.0, T_end - tc) def _eval_design(L2m, L4m, T_env, T_end, h_skin, T_core, T0, dt, N, mp): """给定设计 (mm)，返回 (末端皮肤温度, 超 44°C 时长 s, 首次达 44°C 时刻)。""" s, _, _ = _simulate(L2m / 1000.0, L4m / 1000.0, T_env, h_skin, T_core, T0, dt, T_end, N, mp) return s[-1][1], _time_above(s, T_WARN, T_end), _cross_time(s, T_WARN) def _opt_sa(rng, h_skin, T_core, mp, T_env, T_end, dt, N, T0, thr=T_SAFE, margin=0.0, iters=900, use_warn=True): """模拟退火：最小化总厚度 L1+L2+L3+L4，约束 T_skin(T_end) ≤ thr - margin 且 超 44°C 时长 ≤ 300 s。 margin=0 → 名义最优；margin>0 → 留安全裕度的稳健设计。""" Lfix = (mp["L1"] + mp["L3"]) * 1000.0 tgt = thr - margin warn_tgt = T_WARN_MAX * (1.0 - 0.25 * (margin > 0)) # 稳健设计同时收紧预警时长 def energy(st): L2m, L4m = st te, ta, _ = _eval_design(L2m, L4m, T_env, T_end, h_skin, T_core, T0, dt, N, mp) total = Lfix + L2m + L4m pen = 0.0 if te > tgt: pen += (te - tgt) * 80.0 if use_warn and ta > warn_tgt: pen += (ta - warn_tgt) / 60.0 * 25.0 return total + pen cur = (14.0, 4.0) curE = energy(cur) best, bestE = cur, curE Tinit = 8.0 for it in range(iters): Tt = Tinit * (1 - it / iters) + 0.05 cand = (clip(cur[0] + rng.uniform(-1, 1) * 2.5, L2_LO, L2_HI), clip(cur[1] + rng.uniform(-1, 1) * 0.9, L4_LO, L4_HI)) candE = energy(cand) if candE < curE or rng.random() < math.exp(-(candE - curE) / Tt): cur, curE = cand, candE if curE < bestE: best, bestE = cur, curE L2s, L4s = best te, ta, tc = _eval_design(L2s, L4s, T_env, T_end, h_skin, T_core, T0, dt, N, mp) return (L2s, L4s, Lfix + L2s + L4s, te, ta, tc) def gen_2018a(seed=2018, out_csv="cumcm2018a.csv"): """高温作业专用服装：四层非稳态热传导正演 → 参数辨识 → 单变量优化 → 双变量稳健优化。 返回唯一数据真源 dict；写出合成假人实验时间序列 CSV。""" rng = random.Random(seed) h_skin = 14.0 # 皮肤-体核等效换热系数 W/(m²·K)（含血流对流散热） T_core = 37.0 T0 = 37.0 N = 60 dt = 30.0 mp = dict(MAT) Lfix = (mp["L1"] + mp["L3"]) * 1000.0 # 固定层总厚 mm L2b = mp["L2"] * 1000.0 L4b = mp["L4"] * 1000.0 # ---- 1. 数值格式验证与收敛性 ---- val_err = _validate_fd() # 稳态一致性：长时数值解 vs 串联热阻闭式解（校验多层界面与 Robin 边界处理） steady_chk = [] for (L2c, L4c) in [(6.0, 5.5), (15.0, 2.0), (25.0, 6.4)]: for Te in [65.0, 80.0]: sN, _, _ = _simulate(L2c / 1000.0, L4c / 1000.0, Te, h_skin, T_core, T0, dt, 21600.0, N, mp) ana = _steady_skin(Te, h_skin, T_core, L2c / 1000.0, L4c / 1000.0, mp) steady_chk.append((L2c, L4c, Te, sN[-1][1], ana, abs(sN[-1][1] - ana))) steady_err = max(v[-1] for v in steady_chk) conv = [] for (Nc, dtc) in [(20, 60.0), (30, 30.0), (60, 30.0), (120, 15.0), (240, 7.5)]: sC, _, _ = _simulate(0.017, 0.004, 80.0, h_skin, T_core, T0, dtc, 1800.0, Nc, mp) conv.append((Nc, dtc, sC[-1][1])) conv_ref = conv[-1][2] conv_err = [abs(v - conv_ref) for (_, _, v) in conv] # ---- 2. 各层热阻 / 热容分解（基线 L2=6mm, L4=5.5mm）---- R_parts = [mp["L1"] / mp["k1"], mp["L2"] / mp["k2"], mp["L3"] / mp["k3"], mp["L4"] / k_air_eff(mp["L4"], mp)] C_parts = [mp["r1"] * mp["c1"] * mp["L1"], mp["r2"] * mp["c2"] * mp["L2"], mp["r3"] * mp["c3"] * mp["L3"], mp["r4"] * mp["c4"] * mp["L4"]] R_tot = sum(R_parts) C_tot = sum(C_parts) R_frac = [v / R_tot for v in R_parts] C_frac = [v / C_tot for v in C_parts] tau_base = R_tot * C_tot # 朴素集总估计（后文验证其严重高估） def t63(L2m, L4m, Te): """数值实测时间常数：皮肤温升达到稳态温升 63.2% 所需时间 s。""" ss = _steady_skin(Te, h_skin, T_core, L2m / 1000.0, L4m / 1000.0, mp) s_, _, _ = _simulate(L2m / 1000.0, L4m / 1000.0, Te, h_skin, T_core, T0, dt, 10800.0, N, mp) return _cross_time(s_, T0 + 0.632 * (ss - T0)), ss tau_num_base, ss_base = t63(L2b, L4b, 75.0) # ---- 3. 空气隙辐射修正：有效导热系数与热阻饱和 ---- air_L = [0.6 + 0.4 * i for i in range(16)] # 0.6..6.6 mm air_keff = [k_air_eff(v / 1000.0, mp) for v in air_L] air_R = [(v / 1000.0) / k_air_eff(v / 1000.0, mp) for v in air_L] air_R_sat = 1.0 / H_RAD # ---- 4. 75°C 假人实验：正演 + 合成观测 + 参数辨识 ---- x_II = mp["L1"] + mp["L2"] / 2 x_III = mp["L1"] + mp["L2"] + mp["L3"] / 2 mesh_x = _mesh_x(mp["L2"], mp["L4"], mp, N) s75, ex75, g75 = _simulate(mp["L2"], mp["L4"], 75.0, h_skin, T_core, T0, dt, 5400.0, N, mp, extra_x=[x_II, x_III], rec_every=600 // int(dt)) noise_sd = 0.35 meas = [(t, v + rng.gauss(0, noise_sd)) for (t, v) in s75] prof75 = [(0.0, [T0] * len(mesh_x))] + g75 # 参数辨识：网格 + 局部细化，联合估计 (h_skin, k2) def sse(hh, kk2): mp2 = dict(mp) mp2["k2"] = kk2 s, _, _ = _simulate(mp2["L2"], mp2["L4"], 75.0, hh, T_core, T0, dt, 5400.0, N, mp2) return sum((s[i][1] - meas[i][1]) ** 2 for i in range(len(meas))) best = None for hh in [10.0 + 0.5 * i for i in range(17)]: # 10.0..18.0 for kk2 in [0.25 + 0.03 * i for i in range(9)]: # 0.25..0.49 e = sse(hh, kk2) if best is None or e < best[0]: best = (e, hh, kk2) _, h0, k20 = best for hh in [h0 - 0.4 + 0.1 * i for i in range(9)]: for kk2 in [k20 - 0.024 + 0.006 * i for i in range(9)]: if hh <= 0 or kk2 <= 0: continue e = sse(hh, kk2) if e < best[0]: best = (e, hh, kk2) fit_sse, h_hat, k2_hat = best n_obs = len(meas) fit_rmse = math.sqrt(fit_sse / n_obs) ybar = sum(v for (_, v) in meas) / n_obs sst = sum((v - ybar) ** 2 for (_, v) in meas) fit_r2 = 1.0 - fit_sse / sst mp_fit = dict(mp) mp_fit["k2"] = k2_hat s75_fit, _, _ = _simulate(mp_fit["L2"], mp_fit["L4"], 75.0, h_hat, T_core, T0, dt, 5400.0, N, mp_fit) resid = [meas[i][1] - s75_fit[i][1] for i in range(n_obs)] resid_max = max(abs(v) for v in resid) # 基线在三档环境下的皮肤温度 envs = [65.0, 75.0, 80.0] base_skin = {} base_stat = {} for Te in envs: s, _, _ = _simulate(mp["L2"], mp["L4"], Te, h_skin, T_core, T0, dt, 5400.0, N, mp) base_skin[Te] = s base_stat[Te] = (s[-1][1], _cross_time(s, T_WARN), _cross_time(s, T_SAFE), _steady_skin(Te, h_skin, T_core, mp["L2"], mp["L4"], mp)) # ---- 5. 问题二：65°C / 60min / L4=5.5mm，求最小 L2 ---- P2_ENV, P2_END, P2_L4 = 65.0, 3600.0, 5.5 p2_L2 = [1.0 + 1.0 * i for i in range(25)] # 1..25 mm p2_Tend, p2_above, p2_tau = [], [], [] for L2m in p2_L2: te, ta, _ = _eval_design(L2m, P2_L4, P2_ENV, P2_END, h_skin, T_core, T0, dt, N, mp) p2_Tend.append(te) p2_above.append(ta) p2_tau.append(_R_total(L2m / 1000.0, P2_L4 / 1000.0, mp) * _C_total(L2m / 1000.0, P2_L4 / 1000.0, mp)) # 二分求两条约束各自的临界厚度（0.1mm 精度） def bisect_min(fun, lo, hi, tol=0.05): if fun(hi) > 0: return None while hi - lo > tol: mid = (lo + hi) / 2 if fun(mid) <= 0: hi = mid else: lo = mid return hi f_temp = lambda L2m: _eval_design(L2m, P2_L4, P2_ENV, P2_END, h_skin, T_core, T0, dt, N, mp)[0] - T_SAFE f_warn = lambda L2m: _eval_design(L2m, P2_L4, P2_ENV, P2_END, h_skin, T_core, T0, dt, N, mp)[1] - T_WARN_MAX L2_temp = bisect_min(f_temp, L2_LO, L2_HI) L2_warn = bisect_min(f_warn, L2_LO, L2_HI) L2_p2 = max(v for v in [L2_temp, L2_warn] if v is not None) p2_te, p2_ta, p2_tc = _eval_design(L2_p2, P2_L4, P2_ENV, P2_END, h_skin, T_core, T0, dt, N, mp) p2_total = Lfix + L2_p2 + P2_L4 p2_binding = "预警时长" if (L2_warn or 0) >= (L2_temp or 0) else "峰值温度" p2_marg = [p2_Tend[i] - p2_Tend[i + 1] for i in range(len(p2_Tend) - 1)] # ---- 6. 参数敏感性（+10% 扰动对问题二设计末端温度的影响）---- base_ref, _, _ = _eval_design(L2_p2, P2_L4, P2_ENV, P2_END, h_skin, T_core, T0, dt, N, mp) sens = {} for key in ["k1", "k2", "k3", "k4", "r2", "c2"]: mp2 = dict(mp) mp2[key] = mp[key] * 1.10 te, _, _ = _eval_design(L2_p2, P2_L4, P2_ENV, P2_END, h_skin, T_core, T0, dt, N, mp2) sens[key] = te - base_ref teH, _, _ = _eval_design(L2_p2, P2_L4, P2_ENV, P2_END, h_skin * 1.10, T_core, T0, dt, N, mp) sens["h_skin"] = teH - base_ref # ---- 7. 问题三：80°C / 30min，双变量 (L2,L4) 最小总厚度 ---- P3_ENV, P3_END = 80.0, 1800.0 L2n, L4n, tot_n, te_n, ta_n, tc_n = _opt_sa(rng, h_skin, T_core, mp, P3_ENV, P3_END, dt, N, T0, margin=0.0, use_warn=False) SAFE_MARGIN = 0.3 L2r, L4r, tot_r, te_r, ta_r, tc_r = _opt_sa(rng, h_skin, T_core, mp, P3_ENV, P3_END, dt, N, T0, margin=SAFE_MARGIN, use_warn=False) rob_cost = tot_r - tot_n # 网格 Pareto：不同安全阈值下的最小总厚度 grid_L2 = [1.0 + 1.0 * i for i in range(25)] grid_L4 = [0.6 + 0.4 * i for i in range(15)] cache = {} def ev3(L2m, L4m): key = (round(L2m, 2), round(L4m, 2)) if key not in cache: cache[key] = _eval_design(L2m, L4m, P3_ENV, P3_END, h_skin, T_core, T0, dt, N, mp) return cache[key] pareto = [] for thr in [43.0, 44.0, 45.0, 46.0, 47.0, 48.0]: best_t, best_d = None, None for L2m in grid_L2: for L4m in grid_L4: te, ta, _ = ev3(L2m, L4m) if te <= thr: tot = Lfix + L2m + L4m if best_t is None or tot < best_t: best_t, best_d = tot, (L2m, L4m, te, ta) if best_d: pareto.append((best_t, thr, best_d[0], best_d[1], best_d[2], best_d[3])) # 可行域统计 grid_pts = [(L2m, L4m) for L2m in grid_L2 for L4m in grid_L4] feas_cnt = sum(1 for (L2m, L4m) in grid_pts if ev3(L2m, L4m)[0] <= T_SAFE) feas_frac = feas_cnt / len(grid_pts) grid_n = len(grid_pts) # 等厚度线上 L2/L4 的替代关系（总厚度 ≈ 最优值时，扫 L4 求所需 L2） trade = [] for L4m in [0.6, 1.4, 2.2, 3.0, 3.8, 4.6, 5.4, 6.2]: need = bisect_min(lambda L2m: ev3(L2m, L4m)[0] - T_SAFE, L2_LO, L2_HI) if need is not None: trade.append((L4m, need, Lfix + need + L4m)) # ---- 8. 蒙特卡洛：名义设计 vs 稳健设计 ---- mc_n = 200 def run_mc(L2d, L4d, sub_rng): samples, aboves = [], [] dks = {"k1": [], "k2": [], "k3": [], "k4": [], "h_skin": []} for _ in range(mc_n): f = {kk: 1 + sub_rng.gauss(0, 0.10) for kk in dks} mp2 = dict(mp) for kk in ["k1", "k2", "k3", "k4"]: mp2[kk] = mp[kk] * f[kk] s, _, _ = _simulate(L2d / 1000.0, L4d / 1000.0, P3_ENV, h_skin * f["h_skin"], T_core, T0, dt, P3_END, N, mp2) samples.append(s[-1][1]) aboves.append(_time_above(s, T_WARN, P3_END)) for kk in dks: dks[kk].append(f[kk] - 1) m = sum(samples) / mc_n sd = math.sqrt(sum((v - m) ** 2 for v in samples) / mc_n) srt = sorted(samples) p90 = srt[int(0.90 * (mc_n - 1))] p95 = srt[int(0.95 * (mc_n - 1))] exc = sum(1 for v in samples if v > T_SAFE) / mc_n above_mean = sum(aboves) / mc_n contrib = {kk: _pearson(dks[kk], samples) ** 2 for kk in dks} return samples, m, sd, p95, exc, contrib, above_mean, p90 mc_samples, mc_mean, mc_std, mc_p95, mc_exceed, mc_contrib, mc_above, mc_p90 = \ run_mc(L2n, L4n, random.Random(seed + 1)) rob_samples, rob_mean, rob_std, rob_p95, rob_exceed, rob_contrib, rob_above, rob_p90 = \ run_mc(L2r, L4r, random.Random(seed + 1)) # ---- 机会约束设计：P(T_skin>47) ≤ 10%，用分位数偏移近似 + 事后验证 ---- delta_q = mc_p90 - te_n # 90% 分位相对确定性解的偏移 °C L2c, L4c, tot_c, te_c, ta_c, tc_c = _opt_sa(rng, h_skin, T_core, mp, P3_ENV, P3_END, dt, N, T0, margin=delta_q, use_warn=False) cc_samples, cc_mean, cc_std, cc_p95, cc_exceed, cc_contrib, cc_above, cc_p90 = \ run_mc(L2c, L4c, random.Random(seed + 1)) cc_cost = tot_c - tot_n # 名义/稳健设计的皮肤温度演化（80°C，延长到 60min 看越界时刻） nom_skin, _, _ = _simulate(L2n / 1000.0, L4n / 1000.0, P3_ENV, h_skin, T_core, T0, dt, 3600.0, N, mp) rob_skin, _, _ = _simulate(L2r / 1000.0, L4r / 1000.0, P3_ENV, h_skin, T_core, T0, dt, 3600.0, N, mp) # ---- 9. 单位厚度边际热阻：层II vs 空气隙，及二者交叉点 ---- dR_L2 = 1.0 / mp["k2"] / 1000.0 # (m²K/W)/mm，常数 dR_L4 = [] for v in air_L: ke = k_air_eff(v / 1000.0, mp) dR_L4.append(mp["k4"] / (ke * ke) / 1000.0) # 交叉点：k4/(k4+h_rad·L4)² = 1/k2 → L4* = (sqrt(k4·k2) − k4)/h_rad L4_cross = (math.sqrt(mp["k4"] * mp["k2"]) - mp["k4"]) / H_RAD * 1000.0 dR_ratio_bound = (mp["k4"] / (k_air_eff(L4_HI / 1000.0, mp) ** 2)) * mp["k2"] # ---- 10. 极限设计的 MC 温度包络 → 最长安全作业时长（90% 可靠度）---- L2_max, L4_max = L2_HI, L4_HI env_end = 5400.0 env_steps = int(env_end / dt) + 1 env_rng = random.Random(seed + 7) env_runs = [] for _ in range(mc_n): f = {kk: 1 + env_rng.gauss(0, 0.10) for kk in ["k1", "k2", "k3", "k4", "h_skin"]} mp2 = dict(mp) for kk in ["k1", "k2", "k3", "k4"]: mp2[kk] = mp[kk] * f[kk] s, _, _ = _simulate(L2_max / 1000.0, L4_max / 1000.0, P3_ENV, h_skin * f["h_skin"], T_core, T0, dt, env_end, N, mp2) env_runs.append([v for (_, v) in s]) env_t = [i * dt for i in range(env_steps)] env_p50, env_p90 = [], [] for i in range(env_steps): col = sorted(r[i] for r in env_runs) env_p50.append(col[int(0.50 * (mc_n - 1))]) env_p90.append(col[int(0.90 * (mc_n - 1))]) t_safe90 = _cross_time(list(zip(env_t, env_p90)), T_SAFE) t_safe50 = _cross_time(list(zip(env_t, env_p50)), T_SAFE) # ---- 11. 写出合成假人实验 CSV ---- with open(os.path.join(OUT, out_csv), "w", encoding="utf-8-sig", newline="") as f: w = csv.writer(f) w.writerow(["time_s", "Tskin_meas_75C", "Tskin_model_75C", "T_layerII_mid", "T_layerIII_mid"]) for i in range(len(s75)): w.writerow([round(s75[i][0], 1), round(meas[i][1], 3), round(s75[i][1], 3), round(ex75[0][i][1], 3), round(ex75[1][i][1], 3)]) print("cumcm2018a.csv 已写出：%d 行（75°C 假人实验时间序列，基线 L2=6.0mm/L4=5.5mm）" % len(s75)) print("== 2018A 摘要 ==") print("稳态一致性最大误差=%.5f°C（%d 组工况）" % (steady_err, len(steady_chk))) print("验证误差=%.4f°C | 基线 R=%.4f C=%.0f | 集总tau=%.0fs 实测t63=%.0fs(%.1fmin)" % (val_err, R_tot, C_tot, tau_base, tau_num_base, tau_num_base / 60)) print("辨识 h=%.3f (真值%.1f) k2=%.4f (真值%.3f) RMSE=%.3f R2=%.4f" % (h_hat, h_skin, k2_hat, mp["k2"], fit_rmse, fit_r2)) print("问题二 65°C/60min/L4=5.5: L2_temp=%s L2_warn=%s -> L2*=%.2fmm 总厚=%.2fmm 末端=%.2f°C 超44=%.0fs" % (None if L2_temp is None else round(L2_temp, 2), None if L2_warn is None else round(L2_warn, 2), L2_p2, p2_total, p2_te, p2_ta)) print("问题三 名义 (L2,L4)=(%.2f,%.2f) 总厚=%.2f 末端=%.2f 超44=%.0fs" % (L2n, L4n, tot_n, te_n, ta_n)) print("问题三 稳健 (L2,L4)=(%.2f,%.2f) 总厚=%.2f(+%.2f) 末端=%.2f 超44=%.0fs" % (L2r, L4r, tot_r, rob_cost, te_r, ta_r)) print("MC名义 mean=%.2f sd=%.2f p95=%.2f 越界率=%.3f" % (mc_mean, mc_std, mc_p95, mc_exceed)) print("MC稳健 mean=%.2f sd=%.2f p95=%.2f 越界率=%.3f" % (rob_mean, rob_std, rob_p95, rob_exceed)) print("边际热阻 层II=%.5f/mm 空气隙@6.4mm=%.5f/mm 比值=%.2f 交叉点 L4*=%.2fmm" % (dR_L2, dR_L4[-1], dR_ratio_bound, L4_cross)) print("极限设计(25,6.4) 安全作业时长 p50=%s s p90=%s s(%.1f min)" % (None if t_safe50 is None else round(t_safe50), None if t_safe90 is None else round(t_safe90), (t_safe90 or 0) / 60.0)) print("机会约束(≤10%%越界) (L2,L4)=(%.2f,%.2f) 总厚=%.2f(+%.2f) 确定性末端=%.2f 实测越界率=%.3f 偏移delta=%.2f" % (L2c, L4c, tot_c, cc_cost, te_c, cc_exceed, delta_q)) return dict( seed=seed, h_skin=h_skin, T_core=T_core, T0=T0, N=N, dt=dt, mat=mp, H_RAD=H_RAD, Lfix=Lfix, L2b=L2b, L4b=L4b, bounds=(L2_LO, L2_HI, L4_LO, L4_HI), T_SAFE=T_SAFE, T_WARN=T_WARN, T_WARN_MAX=T_WARN_MAX, val_err=val_err, conv=conv, conv_err=conv_err, steady_chk=steady_chk, steady_err=steady_err, mesh_x=mesh_x, R_parts=R_parts, C_parts=C_parts, R_tot=R_tot, C_tot=C_tot, R_frac=R_frac, C_frac=C_frac, tau_base=tau_base, tau_num_base=tau_num_base, ss_base=ss_base, air_L=air_L, air_keff=air_keff, air_R=air_R, air_R_sat=air_R_sat, s75=s75, meas=meas, prof75=prof75, ex75=ex75, noise_sd=noise_sd, h_hat=h_hat, k2_hat=k2_hat, fit_rmse=fit_rmse, fit_r2=fit_r2, s75_fit=s75_fit, resid=resid, resid_max=resid_max, envs=envs, base_skin=base_skin, base_stat=base_stat, P2_ENV=P2_ENV, P2_END=P2_END, P2_L4=P2_L4, p2_L2=p2_L2, p2_Tend=p2_Tend, p2_above=p2_above, p2_tau=p2_tau, p2_marg=p2_marg, L2_temp=L2_temp, L2_warn=L2_warn, L2_p2=L2_p2, p2_te=p2_te, p2_ta=p2_ta, p2_tc=p2_tc, p2_total=p2_total, p2_binding=p2_binding, sens=sens, base_ref=base_ref, P3_ENV=P3_ENV, P3_END=P3_END, L2n=L2n, L4n=L4n, tot_n=tot_n, te_n=te_n, ta_n=ta_n, tc_n=tc_n, L2r=L2r, L4r=L4r, tot_r=tot_r, te_r=te_r, ta_r=ta_r, tc_r=tc_r, safe_margin=SAFE_MARGIN, rob_cost=rob_cost, pareto=pareto, feas_frac=feas_frac, feas_cnt=feas_cnt, grid_n=grid_n, trade=trade, mc_n=mc_n, mc_samples=mc_samples, mc_mean=mc_mean, mc_std=mc_std, mc_p95=mc_p95, mc_exceed=mc_exceed, mc_contrib=mc_contrib, mc_above=mc_above, rob_samples=rob_samples, rob_mean=rob_mean, rob_std=rob_std, rob_p95=rob_p95, rob_exceed=rob_exceed, rob_contrib=rob_contrib, rob_above=rob_above, mc_p90=mc_p90, rob_p90=rob_p90, delta_q=delta_q, L2c=L2c, L4c=L4c, tot_c=tot_c, te_c=te_c, ta_c=ta_c, cc_cost=cc_cost, cc_samples=cc_samples, cc_mean=cc_mean, cc_std=cc_std, cc_p90=cc_p90, cc_p95=cc_p95, cc_exceed=cc_exceed, cc_above=cc_above, nom_skin=nom_skin, rob_skin=rob_skin, dR_L2=dR_L2, dR_L4=dR_L4, L4_cross=L4_cross, dR_ratio_bound=dR_ratio_bound, env_t=env_t, env_p50=env_p50, env_p90=env_p90, t_safe50=t_safe50, t_safe90=t_safe90, csv=out_csv, ) # ===================== 2018B 智能 RGV 的动态调度策略 ===================== # 场景：8 台 CNC 沿直线排列 + 1 台轨道制导车(RGV) + 上/下料传送带 + 清洗槽。 # 两道工序物料：第1道在奇号 CNC(s1)、第2道在偶号 CNC(s2)；单工序物料任意 CNC 均可。 # 采用离散事件仿真（事件=卸料/装料/清洗交付），RGV 与 CNC 各自为状态机。 RGV_GROUPS = [ dict(name="第1组", move1=20, move2=33, move3=46, P_single=560, P1_two=400, P2_two=378, lu_odd=28, lu_even=31, clean=25), dict(name="第2组", move1=23, move2=41, move3=59, P_single=580, P1_two=280, P2_two=500, lu_odd=30, lu_even=35, clean=30), dict(name="第3组", move1=18, move2=32, move3=46, P_single=545, P1_two=455, P2_two=182, lu_odd=27, lu_even=32, clean=25), ] RGV_SHIFT = 28800 # 8 小时（秒） RGV_REPAIR = (600, 1200) # 故障维修 10~20 min RGV_FAIL_PROB = 0.01 # 单件完成故障概率 ~1% RGV_VAR = 0.10 # 加工时间正态波动 σ/μ RGV_STRATS = ["cyclic", "greedy", "lookahead", "rand"] def rgv_move_time(k, g): k = max(1, int(round(k))) return g["move1"] + (g["move2"] - g["move1"]) * (k - 1) def rgv_LU_arr(g): return [g["lu_odd"] if (i + 1) % 2 == 1 else g["lu_even"] for i in range(8)] def rgv_sim(group, scenario, strategy, shift=RGV_SHIFT, seed=0, fail_prob=0.0, fail_repair=RGV_REPAIR, p_std=0.0, trace=False, max_steps=400000, gdict=None): """离散事件仿真：返回该班次完成的成品数 thru 及状态分解。""" g = gdict if gdict is not None else RGV_GROUPS[group] LU = rgv_LU_arr(g) rng = random.Random(seed) if scenario == "single": roles = ["proc"] * 8 P = lambda r: g["P_single"] else: roles = ["s1" if (i + 1) % 2 == 1 else "s2" for i in range(8)] P = lambda r: g["P1_two"] if r == "s1" else g["P2_two"] cnc = [dict(free_at=0.0, role=roles[i], out=0, busy=0.0, loadc=0) for i in range(8)] rgv = dict(pos=1, t=0.0, carry="raw") completed = []; bd = dict(move=0.0, load=0.0, wait=0.0, deliver=0.0) cyc = 0; steps = 0 tr_rgv = []; tr_cnc = [] rset_of = lambda c: ({"s2"} if c == "semi" else {"proc", "s1"}) while rgv["t"] < shift and steps < max_steps: steps += 1 carry = rgv["carry"]; t = rgv["t"]; pos = rgv["pos"] if carry == "final": rgv["t"] += g["clean"]; bd["deliver"] += g["clean"] rgv["carry"] = None if trace and len(tr_rgv) < 90: tr_rgv.append((round(t, 1), pos, pos, "deliver", "None")) continue rs = rset_of(carry) if strategy == "cyclic": target = None for _ in range(8): if cnc[cyc]["role"] in rs: target = cyc; break cyc = (cyc + 1) % 8 i = target if cnc[i]["free_at"] > t: rgv["t"] = cnc[i]["free_at"] chosen = i; cyc = (i + 1) % 8 else: feas = [i for i in range(8) if cnc[i]["role"] in rs and cnc[i]["free_at"] <= t + 1e-9] if not feas: bt = min((cnc[i]["free_at"] for i in range(8) if cnc[i]["role"] in rs), default=t + 1.0) if bt <= t: bt = t + 1.0 rgv["t"] = bt; bd["wait"] += (bt - t); continue if strategy == "greedy": best = None; bk = 1e18 for i in feas: sa = max(t + rgv_move_time(abs(pos - (i + 1)), g), cnc[i]["free_at"]) if sa < bk - 1e-9: bk = sa; best = i chosen = best elif strategy == "rand": chosen = rng.choice(feas) else: # lookahead：瓶颈感知 + 负载均衡（score=sa-0.5*pri 越小越好） best = None; bk = 1e18; bp = -1e18 for i in feas: sa = max(t + rgv_move_time(abs(pos - (i + 1)), g), cnc[i]["free_at"]) pri = 0.0 if carry == "semi" and cnc[i]["role"] == "s2": pri = 10.0 - cnc[i]["loadc"] elif carry in ("raw", None) and cnc[i]["role"] == "s1": pri = 5.0 - cnc[i]["loadc"] score = sa - 0.5 * pri if score < bk - 1e-9 or (abs(score - bk) < 1e-9 and pri > bp): bk = score; bp = pri; best = i chosen = best i = chosen; role = cnc[i]["role"]; dist = abs(pos - (i + 1)) m = rgv_move_time(dist, g); ta = t + m; sa = max(ta, cnc[i]["free_at"]) wait = max(0.0, cnc[i]["free_at"] - ta) exists = cnc[i]["out"] > 0 ds = "semi" if role == "s1" else "final" co = ds if exists else None if exists and role in ("s2", "proc"): completed.append(sa) if fail_prob > 0 and rng.random() < fail_prob: completed.pop(); cnc[i]["out"] -= 1 rep = rng.uniform(*fail_repair); cnc[i]["free_at"] += rep; cnc[i]["busy"] += rep pp = P(role) if p_std > 0: pp = max(30.0, rng.gauss(pp, p_std * pp)) if trace and len(tr_cnc) < 70 and cnc[i]["out"] == 0: tr_cnc.append((i + 1, round(sa, 1), round(sa + LU[i] + pp, 1), role)) cnc[i]["out"] += 1; cnc[i]["loadc"] += 1 cnc[i]["free_at"] = sa + LU[i] + pp; cnc[i]["busy"] += LU[i] + pp rgv["carry"] = co; rgv["pos"] = i + 1; rgv["t"] = sa + LU[i] bd["move"] += m; bd["load"] += LU[i]; bd["wait"] += wait if trace and len(tr_rgv) < 90: tr_rgv.append((round(t, 1), pos, i + 1, ("%s/%s" % (role, "unload" if exists else "load")), co)) comp = sorted(completed) grid = list(range(0, shift + 1, 1800)) cum = []; j = 0 for gt in grid: while j < len(comp) and comp[j] <= gt: j += 1 cum.append((gt / 3600.0, j)) return dict(thru=len(completed), completed=comp, bd=bd, per_cnc=[c["out"] for c in cnc], busy=[c["busy"] for c in cnc], trace_rgv=tr_rgv if trace else None, trace_cnc=tr_cnc if trace else None, cum=cum) def gen_2018b(seed=2018, out_csv="cumcm2018b.csv"): """智能 RGV 动态调度：离散事件仿真 → 四策略对比 → 不确定性量化 → 稳健设计。 返回唯一数据真源 dict；写出表1 参数表 CSV（覆盖占位资源分配数据）。""" # ---- 1. 确定性基线：四策略 × 三组 ×（单工序/两道）---- thru_single = {}; thru_two = {} for gi in range(3): thru_single[gi] = {}; thru_two[gi] = {} for st in RGV_STRATS: r1 = rgv_sim(gi, "single", st, seed=1) r2 = rgv_sim(gi, "two", st, seed=1) thru_single[gi][st] = r1["thru"]; thru_two[gi][st] = r2["thru"] # 代表运行（two, 第1组, greedy）用于时间分解 / Gantt / 累计曲线 rep = rgv_sim(0, "two", "greedy", seed=1, trace=True) bd_rep = {k: round(v, 1) for k, v in rep["bd"].items()} tot_bd = sum(bd_rep.values()) or 1.0 bd_rep_pct = {k: round(100.0 * v / tot_bd, 1) for k, v in bd_rep.items()} per_cnc_rep = rep["per_cnc"]; busy_rep = [round(b, 1) for b in rep["busy"]] util_rep = [round(b / RGV_SHIFT, 4) for b in rep["busy"]] trace_rep = rep["trace_rgv"]; trace_cnc_rep = rep["trace_cnc"]; cum_rep = rep["cum"] # ---- 2. 单工序：吞吐量 vs CNC 数（第1组 greedy）---- thru_vs_ncnc = [] for n in range(1, 9): g = RGV_GROUPS[0]; LU = rgv_LU_arr(g); rng = random.Random(1) cnc = [dict(free_at=0.0, role="proc", out=0, busy=0.0) for _ in range(n)] rgv = dict(pos=1, t=0.0, carry="raw"); done = []; steps = 0 while rgv["t"] < RGV_SHIFT and steps < 400000: steps += 1; carry = rgv["carry"]; t = rgv["t"]; pos = rgv["pos"] if carry == "final": rgv["t"] += g["clean"]; rgv["carry"] = None; continue feas = [i for i in range(n) if cnc[i]["free_at"] <= t + 1e-9] if not feas: bt = min((cnc[i]["free_at"] for i in range(n)), default=t + 1.0) if bt <= t: bt = t + 1.0 rgv["t"] = bt; continue best = None; bk = 1e18 for i in feas: sa = max(t + rgv_move_time(abs(pos - (i + 1)), g), cnc[i]["free_at"]) if sa < bk - 1e-9: bk = sa; best = i i = best; dist = abs(pos - (i + 1)); m = rgv_move_time(dist, g) ta = t + m; sa = max(ta, cnc[i]["free_at"]); exists = cnc[i]["out"] > 0 if exists: done.append(sa) pp = g["P_single"] cnc[i]["out"] += 1; cnc[i]["free_at"] = sa + LU[i] + pp; cnc[i]["busy"] += LU[i] + pp rgv["carry"] = "final" if exists else None; rgv["pos"] = i + 1; rgv["t"] = sa + LU[i] thru_vs_ncnc.append((n, len(done))) # ---- 3. 参数敏感性（two, 第1组, greedy）：单变量隔离缩放 0.8/1.0/1.2 ---- sens_move = []; sens_lu = []; sens_P = [] base = dict(RGV_GROUPS[0]) for sc in [0.8, 1.0, 1.2]: g = dict(base); g["move1"] *= sc; g["move2"] *= sc; g["move3"] *= sc sens_move.append((sc, rgv_sim(0, "two", "greedy", seed=1, gdict=g)["thru"])) for sc in [0.8, 1.0, 1.2]: g = dict(base); g["lu_odd"] *= sc; g["lu_even"] *= sc sens_lu.append((sc, rgv_sim(0, "two", "greedy", seed=1, gdict=g)["thru"])) for sc in [0.8, 1.0, 1.2]: g = dict(base); g["P1_two"] *= sc; g["P2_two"] *= sc sens_P.append((sc, rgv_sim(0, "two", "greedy", seed=1, gdict=g)["thru"])) # ---- 4. 角色分配（two, 第1组, greedy）：s1 数量 2..6 ---- assign = [] for n_s1 in [2, 3, 4, 5, 6]: g = RGV_GROUPS[0]; LU = rgv_LU_arr(g) roles = ["s1"] * n_s1 + ["s2"] * (8 - n_s1) cnc = [dict(free_at=0.0, role=roles[i], out=0, busy=0.0, loadc=0) for i in range(8)] rgv = dict(pos=1, t=0.0, carry="raw"); done = []; steps = 0 P = lambda r: g["P1_two"] if r == "s1" else g["P2_two"] while rgv["t"] < RGV_SHIFT and steps < 400000: steps += 1; carry = rgv["carry"]; t = rgv["t"]; pos = rgv["pos"] if carry == "final": rgv["t"] += g["clean"]; rgv["carry"] = None; continue rs = ({"s2"} if carry == "semi" else {"proc", "s1"}) feas = [i for i in range(8) if cnc[i]["role"] in rs and cnc[i]["free_at"] <= t + 1e-9] if not feas: bt = min((cnc[i]["free_at"] for i in range(8) if cnc[i]["role"] in rs), default=t + 1.0) if bt <= t: bt = t + 1.0 rgv["t"] = bt; continue best = None; bk = 1e18 for i in feas: sa = max(t + rgv_move_time(abs(pos - (i + 1)), g), cnc[i]["free_at"]) if sa < bk - 1e-9: bk = sa; best = i i = best; role = cnc[i]["role"]; dist = abs(pos - (i + 1)); m = rgv_move_time(dist, g) ta = t + m; sa = max(ta, cnc[i]["free_at"]); exists = cnc[i]["out"] > 0 ds = "semi" if role == "s1" else "final" if exists and role in ("s2", "proc"): done.append(sa) pp = P(role) cnc[i]["out"] += 1; cnc[i]["free_at"] = sa + LU[i] + pp; cnc[i]["busy"] += LU[i] + pp rgv["carry"] = ds if exists else None; rgv["pos"] = i + 1; rgv["t"] = sa + LU[i] assign.append((n_s1, len(done))) # ---- 5. 不确定性蒙特卡洛（two, 故障0.01 + 波动0.1），n=120 ---- def mc_run(group, strat, n=120, fail=RGV_FAIL_PROB, var=RGV_VAR): ss = [rgv_sim(group, "two", strat, seed=s, fail_prob=fail, fail_repair=RGV_REPAIR, p_std=var)["thru"] for s in range(n)] sv = sorted(ss) return dict(mean=round(statistics.mean(ss), 2), sd=round(statistics.pstdev(ss), 2), p5=sv[max(0, n // 20 - 1)], p50=sv[n // 2], p95=sv[-max(1, n // 20)], samples=ss) mc = {}; mc_fail = {}; mc_var = {} for gi in range(3): for st in RGV_STRATS: mc[(gi, st)] = mc_run(gi, st, 120, RGV_FAIL_PROB, RGV_VAR) mc_fail[(gi, st)] = mc_run(gi, st, 120, RGV_FAIL_PROB, 0.0) mc_var[(gi, st)] = mc_run(gi, st, 120, 0.0, RGV_VAR) # ---- 6. 故障率敏感性（two, 第1组, greedy & lookahead）---- fail_sens = {} for st in ["greedy", "lookahead"]: row = [] for fp in [0.0, 0.01, 0.03, 0.05]: ss = [rgv_sim(0, "two", st, seed=s, fail_prob=fp, fail_repair=RGV_REPAIR, p_std=RGV_VAR)["thru"] for s in range(80)] row.append((fp, round(statistics.mean(ss), 1))) fail_sens[st] = row # ---- 7. 写出 CSV（表1 参数表，覆盖占位资源分配数据）---- rows = [ ["参数", "第1组", "第2组", "第3组", "单位"], ["RGV移动1单位", 20, 23, 18, "s"], ["RGV移动2单位", 33, 41, 32, "s"], ["RGV移动3单位", 46, 59, 46, "s"], ["单工序加工", 560, 580, 545, "s"], ["两道·第一道", 400, 280, 455, "s"], ["两道·第二道", 378, 500, 182, "s"], ["上下料(奇CNC)", 28, 30, 27, "s"], ["上下料(偶CNC)", 31, 35, 32, "s"], ["清洗", 25, 30, 25, "s"], ["班次时长", 8, 8, 8, "h"], ] with open(os.path.join(OUT, out_csv), "w", newline="") as f: w = csv.writer(f); w.writerows(rows) return dict( groups=RGV_GROUPS, shift=RGV_SHIFT, n_cnc=8, strats=RGV_STRATS, move_tbl=[[20, 33, 46], [23, 41, 59], [18, 32, 46]], lu_tbl=[[28, 31], [30, 35], [27, 32]], clean_tbl=[25, 30, 25], proc_tbl=[560, 580, 545], two_tbl=[[400, 378], [280, 500], [455, 182]], thru_single=thru_single, thru_two=thru_two, bd_rep=bd_rep, bd_rep_pct=bd_rep_pct, per_cnc_rep=per_cnc_rep, busy_rep=busy_rep, util_rep=util_rep, trace_rep=trace_rep, trace_cnc_rep=trace_cnc_rep, cum_rep=cum_rep, thru_vs_ncnc=thru_vs_ncnc, sens_move=sens_move, sens_lu=sens_lu, sens_P=sens_P, assign=assign, mc=mc, mc_fail=mc_fail, mc_var=mc_var, fail_sens=fail_sens, csv=out_csv, ) def gen_2019a(seed=2019, out_csv="cumcm2019a.csv"): # 2019A 物理模型实现见独立模块 tools/gen_2019a.py（隔离部署，防单点错误影响本文件） from gen_2019a import gen_2019a as _run return _run(seed=seed, out_csv=out_csv) def gen_2019b(seed=2019, out_csv="cumcm2019b.csv"): # 2019B 同心鼓物理模型实现见独立模块 tools/gen_2019b.py（隔离部署，防单点错误影响本文件） from gen_2019b import gen_2019b as _run return _run(seed=seed, out_csv=out_csv) def gen_2019c(seed=2019, out_csv="cumcm2019c.csv"): # 2019C 机场出租车调度模型实现见独立模块 tools/gen_2019c.py（隔离部署） from gen_2019c import gen_2019c as _run return _run(seed=seed, out_csv=out_csv) def gen_2020a(seed=2020, out_csv="cumcm2020a.csv"): # 2020A 炉温曲线传热模型实现见独立模块 tools/gen_2020a.py（隔离部署） from gen_2020a import gen_2020a as _run return _run(seed=seed, out_csv=out_csv) def gen_2020b(seed=2020, out_csv="cumcm2020b.csv"): # 2020B 穿越沙漠资源受限路径规划实现见独立模块 tools/gen_2020b.py（隔离部署） from gen_2020b import gen_2020b as _run return _run(seed=seed, out_csv=out_csv) def gen_2020c(seed=2020, out_csv="cumcm2020c.csv"): # 2020C 中小微企业信贷决策实现见独立模块 tools/gen_2020c.py（隔离部署） from gen_2020c import gen_2020c as _run return _run(seed=seed, out_csv=out_csv) def gen_2021a(seed=2021, out_csv="cumcm2021a.csv"): # 2021A ASF 自动驾驶主动安全实现见独立模块 tools/gen_2021a.py（隔离部署） from gen_2021a import gen_2021a as _run return _run(seed=seed, out_csv=out_csv) def gen_2021b(out_csv="cumcm2021b.csv"): # 2021B 乙醇偶合制备 C4 烯烃实现见独立模块 tools/gen_2021b.py（隔离部署） from gen_2021b import gen_2021b as _run return _run(out_csv=out_csv) def gen_2021c(seed=2021, out_csv="cumcm2021c.csv"): # 2021C 生产企业原材料的订购与运输实现见独立模块 tools/gen_2021c.py（隔离部署） from gen_2021c import gen_2021c as _run return _run(seed=seed, out_csv=out_csv) def gen_2022a(seed=2022, out_csv="cumcm2022a.csv"): # 2022A 波浪能最大输出功率设计实现见独立模块 tools/gen_2022a.py（隔离部署） from gen_2022a import gen_2022a as _run return _run(seed=seed, out_csv=out_csv) def gen_2022b(seed=2022, out_csv="cumcm2022b.csv"): # 2022B 无人机纯方位无源定位实现见独立模块 tools/gen_2022b.py（隔离部署） from gen_2022b import gen_2022b as _run return _run(seed=seed, out_csv=out_csv) def gen_2022c(seed=2022, out_csv="cumcm2022c.csv"): # 2022C 古代玻璃制品成分分析与鉴别实现见独立模块 tools/gen_2022c.py（隔离部署） from gen_2022c import gen_2022c as _run return _run(seed=seed, out_csv=out_csv) def gen_2023a(seed=2023, out_csv="cumcm2023a.csv"): # 2023A 定日镜场的优化设计（塔式光热电站定日镜场：单镜几何/成本最优 + 径向交错镜场布置 + 光学效率分解 + LCOE 经济性）实现见独立模块 tools/gen_cumcm2023a.py（隔离部署） from gen_cumcm2023a import gen_cumcm2023a as _run return _run(seed=seed, out_csv=out_csv) def gen_dgcup2019a(seed=2019, out_csv="dgcup2019a.csv"): # 2019 电工杯 A 电源规划 实现见独立模块 tools/gen_dgcup2019a.py（隔离部署） from gen_dgcup2019a import gen_dgcup2019a as _run return _run(seed=seed, out_csv=out_csv) def gen_dgcup2020a(seed=2020, out_csv="dgcup2020a.csv"): # 2020 电工杯 A 高铁牵引供电系统运行数据分析及等值建模 实现见独立模块 tools/gen_dgcup2020a.py（隔离部署） from gen_dgcup2020a import gen_dgcup2020a as _run return _run(seed=seed, out_csv=out_csv) def gen_dgcup2021a(seed=2021, out_csv="dgcup2021a.csv"): # 2021 电工杯 A 高铁牵引供电系统运行数据分析及等值建模 实现见独立模块 tools/gen_dgcup2021a.py（隔离部署） from gen_dgcup2021a import gen_dgcup2021a as _run return _run(seed=seed, out_csv=out_csv) def gen_dgcup2022a(seed=2022, out_csv="dgcup2022a.csv"): # 2022 电工杯 A 高比例风电电力系统储能运行及配置分析 实现见独立模块 tools/gen_dgcup2022a.py（隔离部署） from gen_dgcup2022a import gen_dgcup2022a as _run return _run(seed=seed, out_csv=out_csv) def gen_dgcup2023a(seed=2023, out_csv="dgcup2023a.csv"): # 2023 电工杯 A 电采暖负荷参与电力系统功率调节的技术经济分析 实现见独立模块 tools/gen_dgcup2023a.py（隔离部署） from gen_dgcup2023a import gen_dgcup2023a as _run return _run(seed=seed, out_csv=out_csv) def gen_tidy2019a(seed=2019, out_csv="tidy2019a.csv"): # 2019 泰迪杯 E 客户分群 实现见独立模块 tools/gen_tidy2019a.py（隔离部署） from gen_tidy2019a import gen_tidy2019a as _run return _run(seed=seed, out_csv=out_csv) def gen_tidy2020a(seed=2020, out_csv="tidy2020a.csv"): # 2020 泰迪杯 D 金融风控分类 实现见独立模块 tools/gen_tidy2020a.py（隔离部署） from gen_tidy2020a import gen_tidy2020a as _run return _run(seed=seed, out_csv=out_csv) def gen_tidy2021a(seed=2021, out_csv="tidy2021a.csv"): # 2021 泰迪杯 C 交通流量预测 实现见独立模块 tools/gen_tidy2021a.py（隔离部署） # 注：该模块内部确定性锚定 SEED=2021，签名仅接收 out_csv，故不传 seed from gen_tidy2021a import gen_tidy2021a as _run return _run(out_csv=out_csv) def gen_tidy2022a(seed=2022, out_csv="tidy2022a.csv"): # 2022 泰迪杯 B 医疗风险预测 实现见独立模块 tools/gen_tidy2022a.py（隔离部署） from gen_tidy2022a import gen_tidy2022a as _run return _run(seed=seed, out_csv=out_csv) def gen_tidy2023a(seed=2023, out_csv="tidy2023a.csv"): # 2023 泰迪杯 A 商业销售预测 实现见独立模块 tools/gen_tidy2023a.py（隔离部署） from gen_tidy2023a import gen_tidy2023a as _run return _run(seed=seed, out_csv=out_csv) def gen_2023b(seed=2023, out_csv="cumcm2023b.csv"): # 2023B 多波束测线问题实现见独立模块 tools/gen_2023b.py（隔离部署） from gen_2023b import gen_2023b as _run return _run(seed=seed, out_csv=out_csv) def gen_2023c(seed=2023, out_csv="cumcm2023c.csv"): # 2023C 蔬菜类商品自动定价与补货决策实现见独立模块 tools/gen_2023c.py（隔离部署） from gen_2023c import gen_2023c as _run return _run(seed=seed, out_csv=out_csv) def gen_2024a(seed=2024, out_csv="cumcm2024a.csv"): # 2024A 板凳龙沿螺线行进实现见独立模块 tools/gen_2024a.py（隔离部署） from gen_2024a import gen_2024a as _run return _run(seed=seed, out_csv=out_csv) def gen_2024b(seed=2024, out_csv="cumcm2024b.csv"): # 2024B 反潜航空深弹投掷实现见独立模块 tools/gen_2024b.py（隔离部署） from gen_2024b import gen_2024b as _run return _run(seed=seed, out_csv=out_csv) def gen_2024c(seed=2024, out_csv="cumcm2024c.csv"): # 2024C 农作物种植策略实现见独立模块 tools/gen_2024c.py（隔离部署） from gen_2024c import gen_2024c as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2024a(seed=2024, out_csv="mcm2024a.csv"): # MCM2024A 资源可用性与性别比例（七鳃鳗）实现见独立模块 tools/gen_mcm2024a.py（隔离部署） from gen_mcm2024a import gen_mcm2024a as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2024b(seed=2024, out_csv="mcm2024b.csv"): # MCM2024B 搜寻潜水器（Searching for Submersibles）实现见独立模块 tools/gen_mcm2024b.py（隔离部署） from gen_mcm2024b import gen_mcm2024b as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2024c(seed=2024, out_csv="mcm2024c.csv"): # MCM2024C 网球的势头（Momentum in Tennis）实现见独立模块 tools/gen_mcm2024c.py（隔离部署） from gen_mcm2024c import gen_mcm2024c as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2022a(seed=2024, out_csv="mcm2022a.csv"): # MCM2022A 自行车手功率分配（Power Profile of a Cyclist）实现见独立模块 tools/gen_mcm2022a.py（隔离部署） from gen_mcm2022a import gen_mcm2022a as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2022b(seed=2024, out_csv="mcm2022b.csv"): # MCM2022B 水电平衡（Water and Hydroelectric Power: Balance）实现见独立模块 tools/gen_mcm2022b.py（隔离部署） from gen_mcm2022b import gen_mcm2022b as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2022c(seed=2024, out_csv="mcm2022c.csv"): # MCM2022C 交易策略（Trading Strategies, 黄金+比特币）实现见独立模块 tools/gen_mcm2022c.py（隔离部署） from gen_mcm2022c import gen_mcm2022c as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2023a(seed=2024, out_csv="mcm2023a.csv"): # MCM2023A 饱经旱灾的植物群落（Drought-stricken Plant Communities）实现见独立模块 tools/gen_mcm2023a.py（隔离部署） from gen_mcm2023a import gen_mcm2023a as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2023b(seed=2024, out_csv="mcm2023b.csv"): # MCM2023B 重新想象马赛马拉（Reimagining Maasai Mara）实现见独立模块 tools/gen_mcm2023b.py（隔离部署） from gen_mcm2023b import gen_mcm2023b as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2023c(seed=2024, out_csv="mcm2023c.csv"): # MCM2023C 预测 Wordle 结果（Predicting Wordle Results）实现见独立模块 tools/gen_mcm2023c.py（隔离部署） from gen_mcm2023c import gen_mcm2023c as _run return _run(seed=seed, out_csv=out_csv) def gen_2018c(seed=2018, out_csv="cumcm2018c.csv"): # 2018C 大型百货商场会员画像实现见独立模块 tools/gen_2018c.py（隔离部署） from gen_2018c import gen_2018c as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2018a(seed=2018, out_csv="mcm2018a.csv"): # MCM2018A 多热泵能效调度实现见独立模块 tools/gen_mcm2018a.py（隔离部署） from gen_mcm2018a import gen_mcm2018a as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2018b(seed=2018, out_csv="mcm2018b.csv"): # MCM2018B 能量生产调度与扩展规划实现见独立模块 tools/gen_mcm2018b.py（隔离部署） from gen_mcm2018b import gen_mcm2018b as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2018c(seed=2018, out_csv="mcm2018c.csv"): # MCM2018C 语言演变实现见独立模块 tools/gen_mcm2018c.py（隔离部署） from gen_mcm2018c import gen_mcm2018c as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2019a(seed=2019, out_csv="mcm2019a.csv"): # MCM2019A 龙舌兰酒生产优化实现见独立模块 tools/gen_mcm2019a.py（隔离部署） from gen_mcm2019a import gen_mcm2019a as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2019b(seed=2019, out_csv="mcm2019b.csv"): # MCM2019B 城市废弃物收运网络优化实现见独立模块 tools/gen_mcm2019b.py（隔离部署） from gen_mcm2019b import gen_mcm2019b as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2019c(seed=2019, out_csv="mcm2019c.csv"): # MCM2019C 阿片成瘾（时间序列+风险分类+干预评估）实现见独立模块 tools/gen_mcm2019c.py（隔离部署） from gen_mcm2019c import gen_mcm2019c as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2020a(seed=2020, out_csv="mcm2020a.csv"): # MCM2020A 北鲑保护（Ricker 双族群补充量+MSY 优化）实现见独立模块 tools/gen_mcm2020a.py（隔离部署） from gen_mcm2020a import gen_mcm2020a as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2020b(seed=2020, out_csv="mcm2020b.csv"): # MCM2020B 沙堡最长寿命（切片仿真+含水率U形+延寿策略）实现见独立模块 tools/gen_mcm2020b.py（隔离部署） from gen_mcm2020b import gen_mcm2020b as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2020c(seed=2020, out_csv="mcm2020c.csv"): # MCM2020C 数据洞察（A Wealth of Data：三商品评论合成+探索/建模/评价）实现见独立模块 tools/gen_mcm2020c.py（隔离部署） from gen_mcm2020c import gen_mcm2020c as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2021a(seed=2021, out_csv="mcm2021a.csv"): # MCM2021A 真菌网络（Fungi：菌丝物质传输 + Murray 锥度 + MST 拓扑优化）实现见独立模块 tools/gen_mcm2021a.py（隔离部署） from gen_mcm2021a import gen_mcm2021a as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2021b(seed=2021, out_csv="mcm2021b.csv"): # MCM2021B 持久耐嚼口香糖（双室风味释放动力学 + 持久度排名 + 加权评分系统）实现见独立模块 tools/gen_mcm2021b.py（隔离部署） from gen_mcm2021b import gen_mcm2021b as _run return _run(seed=seed, out_csv=out_csv) def gen_mcm2021c(seed=2021, out_csv="mcm2021c.csv"): # MCM2021C 确认有蜂群出没（亚洲巨蜂目击报告：扩散 Logistic 成长 + 类加权逻辑回归误判分类 + Beta-Binomial 更新 + 零计数根除规则）实现见独立模块 tools/gen_mcm2021c.py（隔离部署） from gen_mcm2021c import gen_mcm2021c as _run return _run(seed=seed, out_csv=out_csv) def gen_icm2021d(seed=2021, out_csv="icm2021d.csv"): # ICM2021D 音乐的影响（有向影响网络 + 复合影响力分数 + 余弦相似度(控制年代) + 流派演化/革命检测 + 影响有效性 + 文化融合）实现见独立模块 tools/gen_icm2021d.py（隔离部署） from gen_icm2021d import gen_icm2021d as _run return _run(seed=seed, out_csv=out_csv) def gen_icm2021e(seed=2021, out_csv="icm2021e.csv"): # ICM2021E Seeking Solid Ground（地面沉降：水文子模型 dH + 地质力学子模型 v + 社区脆弱性 + 政策有效性）实现见独立模块 tools/gen_icm2021e.py（隔离部署） from gen_icm2021e import gen_icm2021e as _run return _run(seed=seed, out_csv=out_csv) def gen_icm2021f(seed=2021, out_csv="icm2021f.csv"): # ICM2021F Checking the Pulse and Temperature of Higher Education（8维健康可持续指数 HSI + 目标态 + 政策迁移 S 形轨迹）实现见独立模块 tools/gen_icm2021f.py（隔离部署） from gen_icm2021f import gen_icm2021f as _run return _run(seed=seed, out_csv=out_csv) def gen_icm2022d(seed=2022, out_csv="icm2022d.csv"): # ICM2022D Data Paralysis? Use Our Analysis!（D&A 系统成熟度：人/技/流程三维度15子指标 + AHP-熵权组合权重 + 五级定级 + 预算优化 + 关联规则有效性协议 + 跨规模/行业扩展）实现见独立模块 tools/gen_icm2022d.py（隔离部署） from gen_icm2022d import gen_icm2022d as _run return _run(seed=seed, out_csv=out_csv) def gen_icm2022e(seed=2022, out_csv="icm2022e.csv"): # ICM2022E Forestry for Carbon Sequestration（Chapman-Richards 生物量生长 + 采伐-产品碳库（长/短寿命木制品）+ 百年离散模拟；碳最大化轮伐期 + 多目标决策轮伐期；过渡点扫描；多林型应用 + 推荐采伐案例百年固碳）实现见独立模块 tools/gen_icm2022e.py（隔离部署） from gen_icm2022e import gen_icm2022e as _run return _run(seed=seed, out_csv=out_csv) def gen_icm2022f(seed=2022, out_csv="icm2022f.csv"): # ICM2022F All For One, and One (Space) For All（小行星采矿全球公平：Gini/Theil/多维DI 公平度量 + 六分配规则比较 + 全球基金占比/规模敏感性 + 50年三情景动态轨迹）实现见独立模块 tools/gen_icm2022f.py（隔离部署） from gen_icm2022f import gen_icm2022f as _run return _run(seed=seed, out_csv=out_csv) def gen_icm2023d(seed=2023, out_csv="icm2023d.csv"): # ICM2023D Prioritizing the UN SDGs（17×17 SDG 交互网络（ICSU 七分型）+ 度/特征向量/介数中心性 + 综合优先级 + What-if 删点 + 四类危机扰动 + 10年扩散有效性 + 企业ESG子图）实现见独立模块 tools/gen_icm2023d.py（隔离部署） from gen_icm2023d import gen_icm2023d as _run return _run(seed=seed, out_csv=out_csv) def gen_icm2023e(seed=2023, out_csv="icm2023e.csv"): # ICM2023E Light Pollution（光污染风险指数 LPRI = 0.5·HumanRisk + 0.5·EcoRisk，6 指标 L/G/P/B/S/R，16 代表地点覆盖 4 类，风险分级 Low/Mod/High，3 干预策略 I1/I2/I3 削减与最优选址）实现见独立模块 tools/gen_icm2023e.py（隔离部署） from gen_icm2023e import gen_icm2023e as _run return _run(seed=seed, out_csv=out_csv) def gen_icm2023f(seed=2023, out_csv="icm2023f.csv"): # ICM2023F Green GDP（SEEA式 GGDP=GDP−D_res−D_env，32国面板覆盖全球~84%排放，减排响应 r=R_MAX·φ·s，三情景φ=0.3/0.6/1.0 全球影响 ΔC/Δppm/ΔT/效益成本比BCR=39，巴西深析+一页报告）实现见独立模块 tools/gen_icm2023f.py（隔离部署） from gen_icm2023f import gen_icm2023f as _run return _run(seed=seed, out_csv=out_csv) def gen_icm2018d(seed=2018, out_csv="icm2018d.csv"): # ICM2018D Out of Gas and Driving on E（电动汽车充电网络规划：车辆基数→能量需求→快充泊位(需求∪覆盖)→城乡郊分配→投资；九国面板，韩国深析+逻辑斯蒂采用率+密度×财富分类）实现见独立模块 tools/gen_icm2018d.py（隔离部署） from gen_icm2018d import gen_icm2018d as _run return _run(seed=seed, out_csv=out_csv) def gen_icm2019d(seed=2019, out_csv="icm2019d.csv"): # ICM2019D Time to leave the Louvre（卢浮宫应急疏散：5层×4核图模型 + 容量比例分配(规避Braess悖论) + 八情景威胁重配置 + 瓶颈识别 + 急救逆行进入FREM + 跨建筑推广尺度律）实现见独立模块 tools/gen_icm2019d.py（隔离部署） from gen_icm2019d import gen_icm2019d as _run return _run(seed=seed, out_csv=out_csv) def gen_icm2020d(seed=2020, out_csv="icm2020d.csv"): # ICM2020D The Wisdom of Teams（足球传球网络：加权有向邻接 W=角色权重×(角色亲和×空间邻近exp(-d/σ)) + 显著边阈值TH=150 + 密度/互惠/Fagiolo加权聚类/介数·特征向量中心性/直径跳数 + 三元普查 + 五维TPI(多样性/协调性/贡献均衡/适应性/节奏) + 8队对照 + what-if杠杆 + 敏感性）实现见独立模块 tools/gen_icm2020d.py（隔离部署） from gen_icm2020d import gen_icm2020d as _run return _run(seed=seed, out_csv=out_csv) def gen_icm2024d(seed=2024, out_csv="icm2024d.csv"): # ICM2024D 五大湖水问题（五湖链式水量平衡 ODE + 堰流出流 Q=kΔh^1.5 按已知年均水位/流量校准 + 月NBS季节驱动 + 无调控基线60年 + 多目标最优水位 + 安大略Plan2014规则带调控 + 圣克莱尔疏浚情景 + 极端事件/方差/净效益 + 疏浚深度与气候NBS趋势灵敏度）实现见独立模块 tools/gen_icm2024d.py（隔离部署） from gen_icm2024d import gen_icm2024d as _run return _run(seed=seed, out_csv=out_csv) def gen_icm2024e(seed=2024, out_csv="icm2024e.csv"): # ICM2024E 财产保险的可持续性（四区域灾害刻画 λ/p_hit/m_frac → 纯保费 AAL → 风险定价 G=(PP+κσ)/(1-α) → 保费收入比 pti 三档判定 → 气候趋势触线年 → 小型互助社偿付能力破产概率 → 四维可持续评分 S → 四类社区保护减缓措施回收期 → 会安古镇地标组合估值）实现见独立模块 tools/gen_icm2024e.py（隔离部署） from gen_icm2024e import gen_icm2024e as _run return _run(seed=seed, out_csv=out_csv) def gen_icm2024f(seed=2024, out_csv="icm2024f.csv"): # ICM2024F 减少非法野生动物贸易（客户锁定 ASEAN-WEN → 五杠杆乘法削减模型 R(t)=1-Π(1-e_i(1-e^{-t/τ_i})) → 5年推演 R5=33% → 蒙特卡洛5000路径达成概率72% → 25%位移泄漏净全球24.8% → 复杂系统协同+15% → 客户资源缺口指数0.456）实现见独立模块 tools/gen_icm2024f.py（隔离部署） from gen_icm2024f import gen_icm2024f as _run return _run(seed=seed, out_csv=out_csv) if __name__ == "__main__": gen_2012a() gen_2012b() gen_2013a() gen_2013b() gen_2012c() gen_2013c() gen_2014a() gen_2015a() gen_2014b() gen_2015b() gen_2014c() gen_2015c() gen_2016a() gen_2016b() gen_2016c() gen_2017a() gen_2017b() gen_2017c() gen_2018a() gen_2018b() gen_2019a() gen_2019b() gen_2019c() gen_2020a() gen_2020b() gen_2020c() gen_2021a() gen_2021b() gen_2021c() gen_2022a() gen_2022b() gen_2022c() gen_2023a() gen_2023b() gen_2023c() gen_2024a() gen_2024b() gen_2024c() gen_2018c() gen_mcm2018a() gen_mcm2018b() gen_mcm2018c() gen_mcm2019a() gen_mcm2019b() gen_mcm2019c() gen_mcm2020a() gen_mcm2020b() gen_mcm2020c() gen_mcm2021a() gen_mcm2021b() gen_mcm2021c() gen_icm2018d() gen_icm2019d() gen_icm2021d() gen_icm2021e() gen_icm2021f() gen_icm2022d() gen_icm2022e() gen_icm2022f() gen_icm2023d() gen_icm2023e() gen_icm2023f() gen_icm2024d() gen_icm2024e() gen_icm2024f() gen_mcm2024a() gen_mcm2024b() gen_mcm2024c() gen_mcm2022a() gen_mcm2022b() gen_mcm2022c() gen_mcm2023a() gen_mcm2023b() gen_mcm2023c() gen_dgcup2019a() gen_dgcup2020a() gen_dgcup2021a() gen_dgcup2022a() gen_dgcup2023a() gen_tidy2019a() gen_tidy2020a() gen_tidy2021a() gen_tidy2022a() gen_tidy2023a() print("DONE gen_data")