非线性动力学工具箱：相空间重构、李雅普诺夫指数与SDE求解的代码整理实践

做数据分析的人多少都会经历一段“手里全是零碎代码”的时间。我印象很深的是，有一年我手头同时跑着几个非线性动力学方向的小任务：要基于观测时间序列做相空间重构，想估计混沌信号里的李雅普诺夫指数；另一个任务是给一组带噪声的观测数据建立随机微分方程模型，参数还得靠智能算法去拟合。这些工作凑到一起，就逼着我必须把散落各处的相空间重构、时序信号分析、随机微分方程求解和智能算法代码做一次系统性整理，做成一个真正能复用的数据驱动工具箱。

整理完之后，我最大的感受是：算法本身其实都有现成文献，真正花时间的反而是“接口设计”和“验证逻辑”。如果你也在做类似的非线性动力学代码整理，或者正准备给手头的脚本建一个像样的代码库，这篇文章应该能帮你少走不少弯路。我会把项目目录结构、关键算法实现思路、随机微分方程求解器的选型对比，以及智能算法做参数辨识时最容易踩的坑，都按实际整理过程的顺序摊开讲。

1. 为什么要收拾成一个“数据驱动工具箱”，而不是继续堆脚本

1.1 整理前让我崩溃的三个问题

整理前的第一反应是“先把活干完再说”。结果活越干越多，代码库成了这个样子：每个实验都有一份.ipynb，里面从头到尾是一大段顺序代码，从读 Excel 开始，然后是互信息计算、虚假最近邻数量判断、嵌入重构、画图，后面紧跟着 SDE 模拟和参数搜索。换个输入文件就得复制一份 notebook，改几个路径再跑一遍。

这里面的问题是系统性的。第一个问题：功能函数和实验参数全都混杂在一起。比如tau = 7这种延迟时间参数直接写在循环里，下个任务想重新用这段代码，必须人工去翻单元格，很容易改错。第二个问题：不同文件的绘图代码相互冲突，有的用matplotlib默认样式，有的用seaborn，放一起跑就报警告。第三个问题，也是最要命的：没有任何一个模块能独立测试。改一个嵌入维数函数，发现得把整个 notebook 从头执行到尾，中间任何一步改了，后面所有数据都得重新算。

这种情况在科研工程里太常见了。本质上不够是因为前面缺少“结构”，而不是因为“代码水平差”。我开始动手整理之后，先做了一件事：把所有算法按职责分块，保证每个模块只回答一个问题。相空间重构模块只负责从一维时间序列生成高维状态空间；不变量估计模块只负责从重构后的轨迹计算李雅普诺夫指数、关联维数；SDE 求解模块只负责给定漂移项和扩散项之后产生一条模拟路径；智能算法模块只负责在参数空间里搜索使目标函数最小的那组参数。

1.2 整理后的目录结构，可以直接抄

这套结构我后来在多个项目里复用过，简化后的形态如下：

nonlinear_toolkit/ ├── dynamics/ │ ├── phase_space.py # 延迟时间、嵌入维数、相空间重构 │ ├── invariants.py # 最大李雅普诺夫指数、关联维数 │ ├── sde_solver.py # Euler-Maruyama / Milstein / 随机RK4 │ └── identifier.py # 基于智能算法的参数辨识 ├── utils/ │ ├── datasets.py # 加载实测数据、生成仿真数据 │ ├── statistics.py # 互信息、自相关、FFT 等基础工具 │ └── validation.py # 收敛阶测试、K折验证、固定随机种子 ├── configs/ │ └── defaults.yaml # 所有算法的默认参数 ├── tests/ │ ├── test_phase_space.py │ ├── test_sde_order.py │ └── test_identifier.py ├── examples/ │ ├── demo_lorenz.py │ └── demo_ornstein_uhlenbeck.py └── run_pipeline.py # 端到端入口脚本

核心思路是：算法代码放dynamics，工具函数放utils，配置全部集中到configs。这样每次拿到新数据，我只需要改defaults.yaml或者写一个新的入口脚本，算法部分基本不用动。

1.3 三条设计原则，是这次整理最大的收获

第一，接口必须统一。相空间重构、不变量计算、SDE 求解，本质上都是“给一段数组，返回另一段数组或标量”的操作。我放弃了大而全的抽象类，只约定函数签名：输入数组、采样间隔、必要参数，输出要么是np.ndarray，要么是浮点数加一个统计字典。这样理解成本最低，调试也方便。

第二，可复现性靠“固定随机种子+固定配置”。智能算法尤其是重灾区。遗传算法、粒子群优化，结果跟随机种子强相关。如果每次跑出来的参数都不一样，下游分析根本没法推进。后来我在所有随机采样入口地方都加了seed参数，并且在utils/validation.py里封装了一个fixed_seed的上下文管理器，统一管理numpy、random和scipy的随机状态。

第三，每个算法都要有对应的验证用例。相空间重构要用已知的 Lorenz 信号验证重构质量，SDE 求解器要测收敛阶，参数辨识要在已知参数的模拟数据上测试能否找回原参数。没有这一步，模块之间一旦配合出问题，很难定位是哪个环节出错了。

2. 相空间重构：延迟时间和嵌入维数，两张图指路

2.1 Takens 定理给我们的现实约束

相空间重构基于 Takens 嵌入定理：一个足够高维的延迟嵌入能够恢复原系统吸引子的拓扑结构。说白了，如果我们只有一维观测x(t)，不能直接看到整个状态空间，但通过X(t) = [x(t), x(t+τ), ..., x(t+(m-1)τ)]这样构造 m 维向量，在合适的 τ 和 m 下，重构出来的轨迹与真实状态空间轨迹是微分同胚的。

这句话很漂亮，但不实用，因为定理没有告诉 τ 和 m 怎么选。这两个参数直接决定重构质量。τ 取得太小，相邻延迟变量几乎线性相关，状态轨迹会挤在一条对角线附近；τ 取太大，相邻延迟变量又几乎独立，重构出来的点云会因为噪声被过度拉伸。m 取得太小，吸引子没有完全展开，轨迹会自己跟自己交叉；m 一旦取得够大，再增加 d 维度也基本不改变几何结构。

2.2 延迟时间 τ：为什么我优先用互信息而不是自相关

很多教材先说自相关法，即找自相关函数第一次降到1/e的时刻。这个办法在线性平稳信号里勉强能用，但在非线性、混沌信号上经常给出误导性的 τ。因为混沌信号自相关可能衰减得非常慢，或者出现周期性回升，1/e准则选出来的延迟并不代表“时间上分得开”。

互信息则会显式地衡量x(t)与x(t+τ)之间共享的信息量。简单说，如果x(t)已经知道，那x(t+τ)的不确定性剩多少。互信息越小，说明两者独立程度越高。实际做法是扫描不同 τ，画出互信息曲线，取第一个局部极小值对应的 τ。

我在代码里用的是最基础的离散化估计：

import numpy as np def mutual_information_curve(x, max_tau=100, bins=16): def entropy_2d(tau): x1 = x[:-tau] x2 = x[tau:] hist, _, _ = np.histogram2d(x1, x2, bins=bins) p = hist / hist.sum() px = p.sum(axis=1, keepdims=True) py = p.sum(axis=0, keepdims=True) p = p[p > 0] px = px[px > 0] py = py[py > 0] hxy = -np.sum(p * np.log(p)) hx = -np.sum(px * np.log(px)) hy = -np.sum(py * np.log(py)) return hx + hy - hxy return [entropy_2d(t) for t in range(1, max_tau + 1)]

正式项目里我通常还会用pyinform或nolds包做二次校验。用模拟 Lorenz 信号测试时，互信息曲线通常会出现明显的第一谷底，谷底处的 τ 大约在 10 到 20 之间，取决于采样频率。这里有个实操心得：如果你的互信息曲线全程平缓，只有缓慢下降，没有明显谷底，多半是信号噪声太大或数据长度太短，这时候不要硬找一个最小点。建议先做带宽滤波或小波去噪，再回来算互信息。

2.3 嵌入维数 m：虚假最近邻法是更可靠的依据

确定 τ 之后就是嵌入维数 m。最常用的方法是虚假最近邻法，英文缩写 FNN。这个方法的直觉是：如果 m 不够大，吸引子还没有展开，那在高维空间看起来是近邻的两个点，很可能是投影造成的假邻居。我们逐渐增加 m，同时统计“假邻居”的比例降到零附近时的最小 m，那就是合适的嵌入维数。

算法简化之后就是这样：

def nearest_neighbor_distance(x, tau, m, idx): n = len(x) - (m - 1) * tau state = np.array([x[i:i + (m - 1) * tau + 1:tau] for i in range(n)]) target = state[idx] ref = state.copy() ref[idx] = np.inf dist2 = np.sum((ref - target) ** 2, axis=1) return np.sqrt(dist2), dist2, ref def false_neighbors_ratio(x, tau, m_max=12, rtol=10.0): ratios = [] for m in range(1, m_max + 1): n = len(x) - (m - 1) * tau state = np.array([x[i:i + (m - 1) * tau + 1:tau] for i in range(n)]) false_count = 0 for i in range(n): d_ref, dist2, _ = nearest_neighbor_distance(x, tau, m, i) d_next = np.sqrt(np.abs(dist2 + (x[i + m * tau] - x[i + (m - 1) * tau]) ** 2)) if d_ref < 1e-10: continue if d_next / d_ref > rtol: false_count += 1 ratios.append(false_count / n) return ratios

严格实现需要考虑不同嵌入维数下距离尺度变化，以及阈值的自适应选择。不过多数情况下，使用d_{m+1} / d_m > 10这一经典阈值就能得到明确信号。我在整理代码时把 FNN 和互信息做进了同一个模块，使用phase_space_parameters(x, fs, max_tau, m_max)返回一个参数字典，而不是各写各的散装函数。

2.4 τ 和 m 的耦合关系，以及仿真数据校验

很多初入非线性动力学领域的人会问“先算 τ 还是先算 m”。经典流程是先算互信息得到 τ，再在固定 τ 下算 FNN。逻辑是：互信息只涉及两维关系，对 m 不敏感，所以可以先确定。FNN 则依赖距离计算，距离空间的几何结构与 τ 有关，所以必须在固定 τ 下做。

整理完这个模块，我用一个简单规则验证：从 Lorenz 系统采样x分量，长度取 8000 点，先做互信息曲线，得到 τ 大约在 13；然后固定 τ=13，跑 FNN，当 m=3 时假邻居比例降到接近 0。这符合 Lorenz 吸引子的真实相空间维数为 3 的预期。如果得到 m=4 或者更大，通常不是算法算错，而是 τ 太小导致延迟坐标间高度冗余，需要更多维度才能展开吸引子。

3. 时序信号分析：最大李雅普诺夫指数和关联维数的实现细节

3.1 为什么要做“信号分析”，而不是直接拟合方程

相空间重构出来的轨迹，本质上是一个高维几何对象。我们希望量化它的几何特性，比如“轨迹发散得有多快”“吸引子的有效维度是多少”。这就是李雅普诺夫指数和关联维数要做的事。

最大李雅普诺夫指数衡量相邻轨迹的平均指数分离速率。指数为正，说明系统对初始条件敏感，即混沌；指数为零，说明系统在临界状态，典型例子是周期轨道；指数为负，说明系统趋于稳定不动点。关联维数则给出吸引子的分形维度下限，帮助我们判断真实动力学大概是低维的，还是高维到无法用延迟嵌入轻易描述。

这个环节的代码整理难度不在于算法本身，而在于参数多，且对噪声敏感。如果没有统一的输入输出接口，很容易出现“图表很好看，但数值随参数抖动”的尴尬。

3.2 最大李雅普诺夫指数：数据量小就选 Rosenstein 方法

计算最大李雅普诺夫指数的经典方法有两类：Wolf 方法直接跟踪状态空间中最邻近的一对轨道，另一种是 Rosenstein、Collins、De Luca 提出的方法，对小数据集更稳健。我在整理时默认选了 Rosenstein 方法，因为它不需要重建整个嵌入向量，只需要对时序信号的延迟坐标做近邻搜索，然后对每个时间点计算相邻轨迹分离度。

简化的实现骨架如下：

def largest_lyapunov_rosenstein(x, fs, tau, m, min_tsep=10, max_iter=100): n = len(x) - (m - 1) * tau state = np.array([x[i:i + (m - 1) * tau + 1:tau] for i in range(n)]) import scipy.spatial tree = scipy.spatial.cKDTree(state) # 为每个参考点找最近邻，同时避开时间太接近的点 # 记录不同时刻 log(平均分离距离) # 返回 (t, log_distance) 曲线，线性段的斜率就是最大李雅普诺夫指数

真实代码比这个骨架复杂，本质上需要处理近邻搜索、间隙排除和线性段选择三个部分。我最开始用暴力法找最近邻，数据长度 20000 点时非常痛苦。换cKDTree之后，同样的计算量从分钟级压到秒级。这段经验写进代码注释之后，每次重构参数都稳定很多。

计算完log(分离距离)对时间步的曲线后，要找一段线性增长的区域做线性拟合。实际操作时，线性段往往只有大约 20 到 60 个时间步，选错了区域会得到错误的指数。我加了一个简单启发式：只取曲线前 60% 的部分，因为混沌信号后期分离饱和，曲线会变平，这时候再拟合就会低估。

3.3 关联维数：把“多少邻居翻倍增长对应多少个盒子”这件事编码化

关联维数的 Grassberger-Procaccia 算法是另一个高价值模块。它统计不同尺度 r 下，有多少点对的距离小于 r，得到关联积分C(r)。在尺度范围内，C(r) ~ r^D，两边取对数就是一条直线，斜率就是关联维数 D。

def correlation_dimension(x, tau, m, r_min=None, r_max=None, n_r=40): n = len(x) - (m - 1) * tau state = np.array([x[i:i + (m - 1) * tau + 1:tau] for i in range(n)]) # 用KDTree或者分块矩阵计算距离矩阵 # 对每个 r 求 C(r)，然后对 log(r)-log(C(r)) 做线性拟合

这里有一个必须注意的问题：对大矩阵直接计算距离会导致内存爆炸，10000 个点的距离矩阵就有上亿个值。我在整理时采用分块策略，每次只处理一个行块，然后累加C(r)。

关联维数对嵌入维数 m 有“饱和效应”：当 m 足够大时，估计出的关联维数会稳定在某一个值附近。这个稳定值就是系统的内在关联维数。如果曲线无法饱和，说明数据长度不够，或者系统本质上是高维/高噪声。我的建议是至少用 3 个不同的 m 做对比，如果 D 随 m 一直线性增长，那大概率不是低维动力学。

3.4 时间序列长度和噪声：信号分析模块里最硬的约束

时序信号分析三个指标的共同弱点是“数据长度不够”。很多 LFP、脑电、金融高频数据在真实项目里只有几千个采样点，这时候最大李雅普诺夫指数的误差会很大，关联维数的估计会偏低。我没有办法变出不存在的数据，但可以做两件事来让结果更可信：一是采用自主采样法，对原始序列做多次重采样，给每个指数一个置信区间；二是先对信号做相位随机化，把真实信号与替代数据的指数分布范围对比，看真实指数是否显著跳出替代数据的范围。

这在代码里也就是utils/validation.py里的一个bootstrap_indicator函数，输入x、指标计算函数和重采样次数，输出平均值和标准差。整理代码的时候顺手把这一层加进去，避免了下游把所有结论搭在偶然的单个数值上。

4. 随机微分方程求解：Euler-Maruyama、Milstein 和随机 Runge-Kutta 哪个更靠谱

4.1 什么时候需要从“时间序列分析”切换到“SDE 建模”

如果只是描述数据，相空间重构和不变量分析已经够用。但数据驱动项目经常会走到下一步：不仅要描述动力学性状，还要建立一个能仿真的模型，生成更多路径来测试控制策略。到这一步，随机微分方程就上场了。

SDE 的通用形式是dX(t) = a(X, t) dt + b(X, t) dW(t)，其中a是漂移项，描述确定性趋势；b是扩散项，描述噪声强度；dW(t)是布朗运动增量。金融里常用几何布朗运动，生物里常用带乘性噪声的 Logistic 型方程，工程里最常见的是 Ornstein-Uhlenbeck 过程。

4.2 三种求解器，以及为什么不能直接套用普通 Runge-Kutta

初学 SDE 最常犯的错误是直接把确定性 ODE 的 Runge-Kutta 方法套到带噪声的方程上。经典 RK4 对布朗运动项的处理并不一致，换句话说，用普通 Runge-Kutta 格式生成的路径并不按 Itô 积分意义上的解收敛到正确 SDE。必须用随机数值格式。

Euler-Maruyama 格式 是最简单的强近似格式：

X_{n+1} = X_n + a(X_n, t_n) Δt + b(X_n, t_n) ΔW_n

其中ΔW_n是均值 0、方差Δt的独立高斯随机变量。这写起来简单，强收敛阶只有 0.5。

Milstein 格式 在标量噪声情况下增加一个二次修正项：

X_{n+1} = X_n + a Δt + b ΔW_n + 0.5 * b * b_x * (ΔW_n^2 - Δt)

其中b_x是扩散项关于状态的偏导。这里多出来的(ΔW_n^2 - Δt)项修正了固有不连续性，把强收敛阶提高到 1.0。

随机 Runge-Kutta 本质上是一类需要满足随机 Taylor 级数匹配条件的格式。单用“路径积分平均”的方式理解它最直观：随机 RK4 不是为了模仿确定性 RK4，而是为了让轨迹的统计矩在更高阶上匹配 Itô 积分。在代码实现里，它比 Milstein 复杂，在噪声是常数或平滑函数时，优势并不显著。

我整理出的三类求解器都遵循同一个接口：

def solve_sde( init, drift_func, diffusion_func, t_end, dt, method="milstein", seed=None, save_every=1 ):

drift_func和diffusion_func都接收(t, x)两个参数，返回与x同形状的数组。这样任意方程式都能通过临时匿名函数接进来，不需要为每个模型重写求解器。

4.3 收敛阶测试：这部分代码写起来比求解器本身更重要

真正让求解器可靠的不是“看起来像那么回事”，而是数值实验验证收敛阶。我在tests/test_sde_order.py里写了一个通用测试流程：

1. 选一个有解析解的 SDE，比如几何布朗运动，漂移项a = μX，扩散项b = σX；
2. 分别用dt = 0.01, 0.005, 0.0025, 0.00125求解；
3. 对每一条路径，记录终点值；
4. 用独立生成的细网格参考解做对比，计算强误差；
5. 对整个路径而非单点做误差比较，会得到更全面的强收敛测试；
6. 拟合强误差对dt的斜率，Euler-Maruyama 应该在 0.5 附近，Milstein 应该在 1.0 附近。

def strong_error_for_dt(solver_func, params, ref_path, dt): paths = [solver_func(**params, dt=dt, seed=s) for s in range(50)] return np.sqrt(np.mean((np.array(paths) - ref_path) ** 2, axis=0)[-1])

这个测试在整理阶段确实抓出过 bug。我一开始把 Milstein 的修正项写成了0.5 * b * b_x * ΔW_n**2，漏了- Δt这一项。单跑一条路径看不出问题，两条路径也看不出来，降低 dt 之后强误差完全不下降。把修正项从ΔW^2改成(ΔW^2 - Δt)后，收敛阶才恢复。如果没有收敛阶测试，这类错误会直接被带到下游参数辨识。

4.4 求解器选型：别迷恋高阶，先看扩散项

做完收敛阶对比，我整理了一张选型表，塞进sde_solver.py的 docstring 里，作为给后续项目使用时的参考：

实际项目中，我经常用 Euler-Maruyama 做探索性仿真，用 Milstein 做正式实验。能上 Milstein 就尽量上 Milstein，尤其是在扩散项不是常数的时候。随机 RK4 在单变量、扩散项平滑的情况下效果不错，但代码复杂度高，调试成本不值得普通项目支付。

5. 智能算法参数辨识：搜参数不是玄学，是目标函数设计问题

5.1 数据驱动参数辨识的本质是做优化

拿到观测时间序列，想反推 SDE 参数，一般是三个步骤：先选择一个候选模型结构，再定义目标函数，最后用优化算法找最优参数。很多人觉得智能算法很玄，其实它的作用很朴素——在目标函数没有解析梯度的背景下，帮你在参数空间里搜索可行解。

整理代码时，我直接把智能算法模块设计成“目标函数 + 约束 + 优化器”三件套。优化器可以用种群类算法，也可以用简单的网格搜索。关键不是选哪个算法，而是选哪个目标函数。这可能是这次代码整理中最重要的领悟。

5.2 目标函数设计：不能简单比较“模拟路径和真实路径逐点相等”

一开始我试过最直接的做法：让模拟器生成一条路径，然后求模拟路径与真实观测路径的逐点均方误差。这个做法有严重问题。SDE 是一条随机路径，两条随机路径之间的逐点差异往往很大，即使模型完全正确，逐点 MSE 也可能很大。真实情况下，观测还有测量噪声，逐点匹配会逼迫优化器去拟合噪声，最后得到错误参数。

我在这里换了一个策略，把目标函数建立在“统计特征”上。常用的特征有三类：平稳分布的均值和方差、自相关函数取几个滞后步长、以及最大李雅普诺夫指数之类的不变量。代码看起来就像这样：

def objective(params, data, features="moments_acf"): sim = solve_sde( init=data[0], drift_func=lambda t, x: params[0] * x, diffusion_func=lambda t, x: params[1], t_end=len(data) / fs, dt=0.01, method="milstein", seed=0, ) # 计算观测数据的均值和ACF # 计算模拟路径的均值和ACF # 返回两者归一化的距离

这样的目标函数，对模拟路径的初始状态不那么敏感，也天然忽略了噪声的逐点随机性。

5.3 为什么我默认用差分进化而不是手写遗传算法

手写遗传算法不难，但在参数维度小于 15 的 SDE 模型里，scipy提供的differential_evolution已经足够强大。差分进化本质上是一种带变异和交叉的种群搜索，它不需要梯度，支持参数边界约束，也天然支持多线程并行。

from scipy.optimize import differential_evolution res = differential_evolution( objective, bounds=[(0.001, 2.0), (0.0001, 1.0)], seed=42, maxiter=200, popsize=15, workers=-1, polish=True, ) print(res.x, res.fun)

polish=True表示在收敛后用局部优化器做一轮精修。差分进化的全局搜索能力强，但最终落点精度一般，精修一下能显著提高参数辨识准确度。如果目标是快筛，则可以设定maxiter=50，并行开很多进程，先跑几轮看参数大概落点，再缩小边界范围做第二轮搜索。

5.4 辨识结果不能只看一组参数，要加“扰动稳定性”检查

智能算法在参数空间里找一个极小值很容易，但找到的极小值是不是稳定，取决于目标函数曲面长什么样。最常见的坑是，真实 SDE 有多个参数组合产生几乎相同的统计特征，也就是参数不可辨识性。比如扩散项几乎为零时，漂移项的某些参数组合在数据里根本不体现，优化器会随便给一个值，但依然把代价函数据降到很低。

我在identifier.py里加了一步“扰动稳定性检查”：把优化得到的最优参数做 ±20% 扰动，重新评估目标函数，如果目标函数变化很小，说明该参数对数据特征不敏感，那这个参数就不应该被过度解读。输出结果除了最优参数，还会附带一列“可辨识性得分”。这步让代码整理项目真正从“能跑”变成“结果可信”。

6. 端到端整理后的工作流，以及我踩过的几个实测坑

6.1 从 CSV 到参数估计：一个命令完成整个流程

模块化整理完成后，我把所有内容封装成了一个run_pipeline.py入口脚本，逻辑非常直白：读取观测数据，进行统一的预处理，计算互信息和 FNN 确定重构参数，然后重构相空间，再估计不变量；如果需要 SDE 建模，就调用求解器和参数辨识模块；最后把主要结果写进一个报告目录。

python run_pipeline.py \ --input data/cog_signal.csv \ --fs 100 \ --tau 12 \ --m 3 \ --sde-model uhlenbeck \ --optimizer differential_evolution \ --seed 7

这个脚本的好处是把散落在 notebook 里的过程固化成一条流水线。参数不再埋在代码里，而是通过命令行或配置文件显式指定。实际使用时，我依然保留 notebook 做探索性可视化，但正式实验都用这份脚本跑，结果全部落到results/目录，方便对比。

6.2 我在整理后最容易踩的两个坑

第一个坑是“重构参数和 SDE 模拟的数据格式不一致”。相空间重构对一维观测数据要求是列向量，但 SDE 模拟器的输出往往是二维数组，行是时间步，列是状态变量。第一次整合时我用错了reshape，导致下游分析出现维度错位，而且错位还很隐蔽，因为一维二维数组在某些 numpy 操作中能自动广播，最终结果只是参数探测不出来，而不是程序直接报错。后来我在datasets.py里统一规定：所有模块输入输出一律是shape = (样本数, 状态变量数)的二维数组，一维信号先补成一个二维列，再进入模块。

第二个坑是“固定随机种子不代表完全固定”。scipy.optimize.differential_evolution的随机状态虽然受seed参数控制，但如果我们在线程并行workers=-1模式下运行，每个进程的随机序列仍然可能不完全一致。只要目标函数内部有用到随机抽样，就必须在用户提供的入口处重新设置全局随机种子。整理后期，我写了一个fixed_seed上下文管理器，才彻底解决 batch 跑时的复现问题。

6.3 后续还能怎么扩展

整理完这个数据结构之后，再往里面加新算法就非常顺手。比如我还想在参数辨识模块中加入基于神经网络的代理模型，把昂贵的目标函数近似掉，从而把差分进化的种群规模提升一个量级。另外，SDE 求解模块未来可以扩展到带跳过程的状态相关噪声，只需要新增一个跳跃扩散求解器，保持现有的solve_sde接口不变就行。模块化的价值就在于：新算法永远是在旁边加一块积木，而不是掀翻整个桌子重来。

整理这批代码花了我大概三周时间，其中大约一半时间几乎都在做旧代码的迁移和测试，而不是在写新算法。但这段整理工作带来的收益是长期稳定的：后续每个新数据集，从拿到 CSV 到完成动力学分析，时间从可能的三天压缩到了三个小时以内。如果你手里也有一堆相空间重构、李雅普诺夫指数或 SDE 拟合的脚本，别急着堆下一个 notebook，先在结构上下点功夫，后面会比自己预想的省力得多。