Logistic混沌MATLAB实战：SINDy稀疏辨识与慢时间尺度发现

简介这是一份面向非线性动力学、混沌分析与数据驱动建模学习者的 MATLAB 开源代码包围绕 SINDy 方法实现慢时间尺度动态发现可用于复现 Bramburger、Dylewsky 与 Kutz 在 Physical Review E 2020 年论文中的玩具模型实验与图示。包内共 10 个文件以 7 个 .m 脚本为主配合 2 个 .mat 数据文件和 1 个 README 说明压缩包约 24.17MB脚本覆盖玩具模型数值积分、连续时间 SINDy 建模、粗粒度演化离散映射发现以及木星、土星、Logistic 混沌等示例便于按模块理解算法流程与数据组织。已有 381 人学习下载适合具备一定 MATLAB 基础、希望把 SINDy、滑动窗口 DMD 与慢发现方法落到代码实践的研究生或科研人员参考读者可获得可运行脚本、示例数据与复现思路用于对照论文结果并搭建自己的动力学发现流程。1. 先把混沌序列变成可辨识对象这份 Logistic SINDy 代码解决什么问题这份 Logistic 混沌 MATLAB 代码把 SINDy 稀疏辨识用在慢时间尺度动态发现上不是跑一条混沌序列看波形而是把快变量里的慢结构拎出来再用稀疏回归写成可读的方程。它解决的是“数据看起来乱、但慢流有低维结构”的场景先做延迟嵌入再估导数、建多项式库、用 STLSQ 迭代阈值最后用一步预测和长期轨迹验证。适合已经会一点 MATLAB、想做非线性系统辨识的从业者也适合把系统开源代码拆开复现的学生。新手先跑数据生成和库构建熟手重点看慢窗口与阈值扫描的边界。2. 从 Logistic 映射到延迟嵌入数据生成、采样与导数估计2.1 为什么先做延迟嵌入而不是直接回归Logistic 映射本身是一维离散迭代形式不复杂但直接把 (x_{n1}-x_n) 当成导数、只拿 (x_n) 一列做回归得到的只是一维映射的近似很难把相空间里的结构展开。SINDy 需要的是状态矩阵和对应的导数矩阵它假设 (\dot{X}\Theta(X)\Xi)其中 (\Theta(X)) 是候选函数库。如果状态只有一维库里的交叉项、二次项就缺少可解释的输入慢时间尺度也更难分离。延迟嵌入的作用是把一维序列按延迟 (\tau) 和嵌入维数 (d) 构造成多列状态让后续稀疏回归有更大的候选空间。常见做法是先看自相关或互信息选 (\tau)再用虚假最近邻选 (d)。在这份 Logistic 混沌 MATLAB 代码里我一般先用 (\tau5)、(d3) 跑通链路再回头扫参数。原因很实际Logistic 映射的混沌区对初值敏感嵌入参数稍微变化导数估计和库矩阵的共线性都会变。(\tau) 太小相邻列高度相关最小二乘解不稳定(\tau) 太大列之间关系被拉散局部动态被抹平。(d) 太小慢结构展不开(d) 太大库列数迅速膨胀阈值稍微一大就全零。这里不是追求理论最优而是先让 SINDy 的系数路径可读再谈物理解释。还要注意Logistic 映射是离散迭代不存在天然连续时间导数。代码里把迭代步长设为 (dt1)得到的是“每步变化量”不是物理时间导数。如果换成连续系统采样数据(dt) 必须跟采样周期一致否则系数尺度会整体偏移后面阈值扫描会误判。这个细节在复现时最容易翻车同一套 (\lambda) 在 (dt1) 下能用换成 (dt0.01) 后系数放大一百倍稀疏阈值就显得过小模型不稀疏长期轨迹还会漂。2.2 生成混沌序列与控制采样参数先跑数据生成不要一上来就调 SINDy。Logistic 映射 (x_{k1}r x_k(1-x_k)) 在 (r3.8) 附近已经进入混沌区但初值、暂态长度和序列长度都会影响后续嵌入。下面这段代码把生成、截断暂态和后续嵌入放在一个脚本里参数集中放在开头方便一次改完再跑。r 3.8; % 混沌参数常见起点 3.7~4.0 N 20000; % 总迭代点数太短会导致验证段不足 x zeros(N,1); x(1) 0.413; % 非零且不等于 1避免落入不动点 for k 1:N-1 x(k1) r * x(k) * (1 - x(k)); end x x(2001:end); % 丢弃前 2000 点暂态这段代码的逻辑很直白先给一个非零初值迭代出整条序列再砍掉前 2000 点。参数 (r) 决定系统在周期、混沌还是满混沌区间(r3.8) 适合做混沌动态演示(N) 不能太小因为延迟嵌入会消耗 ((d-1)\tau) 个点训练/验证切分还要留 gap。初值不要取 0、1 或 0.5 附近的特殊点否则容易短时间落进不动点或周期窗口后面求导全是噪声。丢弃暂态的长度也不是玄学2000 点对 20000 点总量来说够用若 (r) 接近分岔点建议加到 5000 点再截断。接下来做延迟嵌入。这里把 (x) 变成 (M\times d) 的状态矩阵每一列是不同延迟的序列。代码里用下标向量直接拼列避免循环里反复拼接大矩阵。tau 5; % 延迟步数 d 3; % 嵌入维数 M numel(x) - (d-1)*tau; X zeros(M,d); for j 1:d X(:,j) x((1:M) (j-1)*tau); endM是嵌入后可用样本数X(:,1)是当前时刻X(:,2)是往前推 (\tau) 步X(:,3)再往前推 (2\tau) 步。参数 (\tau5) 和 (d3) 不是固定答案只是让库矩阵列数可控的起点。若后面发现系数路径很不稳定先把 (\tau) 在 1 到 10 之间扫一遍看自相关第一次过零的位置再把 (d) 从 2 加到 4观察验证误差是否明显下降。不要同时扫 (\tau)、(d) 和 (\lambda)否则出问题时分不清是嵌入还是阈值导致的。2.3 用 Savitzky-Golay 求导和训练/验证切分SINDy 对导数估计很敏感。直接对原始混沌序列做一阶差分噪声会被放大先平滑再求导更稳。常见做法是用 Savitzky-Golay 滤波做局部多项式平滑再对平滑后的状态求梯度。下面这段代码把平滑、求导、按时间块切分和 gap 留白一次做完。dt 1; % 迭代步长连续采样数据要改成实际采样周期 win 11; % 平滑窗口必须为奇数 order 3; % 局部多项式阶数要小于 win if exist(sgolayfilt,file) ~ 2 error(缺少 sgolayfilt请安装 Signal Processing Toolbox 或换平滑函数); end Xs sgolayfilt(X, order, win); dX gradient(Xs, dt); gap (d-1)*tau win; % 防止延迟嵌入重叠造成验证段泄漏 nTrain round(0.7*size(X,1)); trainIdx 1:nTrain; valIdx (nTraingap1):size(X,1); if isempty(valIdx) error(验证段为空减小 gap 或增加数据长度); end Xtrain X(trainIdx,:); dXtrain dX(trainIdx,:); Xval X(valIdx,:); dXval dX(valIdx,:);这里的win11、order3是常用起点。窗口太小平滑不够导数里全是高频抖动SINDy 会挑出一堆虚假项窗口太大真实动态被抹平慢时间尺度还没分离就先被滤波吃掉。order一般取 2 或 3必须小于win否则局部多项式退化成插值。gradient(Xs, dt)对每一列求中心差分边界用单侧差分精度比手写前向差分好一些。切分方式要重点说不能随机打乱样本。延迟嵌入后的相邻样本共享原始点随机切分会让训练集和验证集“见过同一段轨迹”验证误差看起来很好长期轨迹却发散。这里按时间先后切 70/30并在训练段和验证段之间留gap把延迟嵌入重叠和滤波窗口的影响隔开。gap的计算里加了win是因为平滑窗口也会把边界信息带过去。这个细节不做后面调阈值会得到很虚的乐观结果。提示如果手头没有 Signal Processing Toolbox可以把sgolayfilt换成smoothdata(X,gaussian,win)先跑通再手写局部多项式回归替代。不要用未平滑的diff直接喂给 STLSQ除非数据本身几乎没有噪声。3. SINDy 稀疏回归落地多项式库、STLSQ 与阈值扫描3.1 库函数怎么选常数、线性、二次项够不够SINDy 的核心是把候选函数库 (\Theta(X)) 建出来再用稀疏回归挑项。Logistic 映射本身含 (x(1-x))展开后是线性项加二次项所以这份代码用“常数 线性 二次”作为默认库是合理的。库太小真实动态写不进去库太大列之间共线性严重阈值一调就全零或全满。下面这个poly_library函数按列拼库同时返回项名方便后面看系数时知道对应哪一项。function [Theta, names] poly_library(X, polyorder) if isvector(X) X X(:); end [n,m] size(X); cols {ones(n,1)}; names {1}; for i 1:m cols{end1} X(:,i); names{end1} sprintf(x%d,i); end if polyorder 2 for i 1:m for j i:m cols{end1} X(:,i).*X(:,j); names{end1} sprintf(x%d*x%d,i,j); end end end Theta cell2mat(cols); end这段函数的输入X可以是 (n\times m) 矩阵也可以是单行向量如果是向量会先转成行向量保证长期轨迹积分时也能调用。polyorder2时(m3) 的库列数是 (13610)分别是常数、三个线性项、三个平方项和三个交叉项。若后面换成耦合系统可以把polyorder提到 3但列数会涨到 20阈值扫描和共线性检查必须跟上。常见做法是先固定二次库把系数路径跑出来再决定是否加三次项不要一上来堆三角函数和指数项除非有明确的物理先验。库建好后先做一次普通最小二乘看量级。如果某些列的量级差了几个数量级STLSQ 的硬阈值会偏向大尺度列小尺度项还没轮到就被压成零。Logistic 序列在 ([0,1]) 内二次项和线性项量级接近问题不大换成多变量工程数据时我一般会先对Theta的每一列做标准化回归完再把系数换回原尺度。这个步骤不做阈值扫描曲线会看起来很奇怪训练误差降不下去非零项却很少。3.2 STLSQ 迭代阈值实现STLSQ 的思路不复杂先最小二乘求初始系数把绝对值小于 (\lambda) 的系数置零再在剩余活跃项上重新最小二乘重复若干轮。下面这个函数按输出维度逐个回归适合多变量状态。代码里用硬阈值不搞软阈值因为工程上更容易解释“这一项被删了还是留了”。function Xi stlsq(Theta, dX, lambda, maxIter) p size(Theta,2); n size(dX,2); Xi zeros(p,n); for target 1:n coef Theta \ dX(:,target); for it 1:maxIter small abs(coef) lambda; coef(small) 0; active ~small; if ~any(active) break; end coef(active) Theta(:,active) \ dX(:,target); end Xi(:,target) coef; end endTheta是库矩阵dX是导数矩阵lambda是稀疏阈值maxIter防止死循环。每一轮先按当前系数做硬阈值再只保留活跃列做最小二乘。注意small是基于上一轮系数算的活跃列重新回归后可能出现新的小系数下一轮会继续压。这个迭代通常 5 到 10 轮就稳定maxIter10够用如果跑到 10 轮还在变多半是库列共线性太强或阈值落在临界点先检查库而不是继续加迭代。调用时先建训练库和验证库再跑一次固定阈值看看量级polyorder 2; [Theta, names] poly_library(Xtrain, polyorder); [ThetaVal, ~] poly_library(Xval, polyorder); lambda 0.02; maxIter 10; Xi stlsq(Theta, dXtrain, lambda, maxIter); predTrain Theta * Xi; predVal ThetaVal * Xi; relTrain norm(predTrain - dXtrain, fro) / norm(dXtrain, fro); relVal norm(predVal - dXval, fro) / norm(dXval, fro);这里用 Frobenius 范数算相对误差能同时看所有状态维度的整体拟合。lambda0.02只是起点不是推荐值。若nnz(Xi)很少、相对误差很大先把 (\lambda) 降到 (10^{-4}) 看最小二乘解若nnz(Xi)接近库列数、验证误差也不低说明阈值没起作用库可能太大或导数估计有问题。不要把训练误差最小当成目标SINDy 要的是稀疏且验证稳定的系数。3.3 阈值与系数稳定性怎么读阈值扫描是这套流程里最值得花时间的一步。做法是让 (\lambda) 在对数空间从 (10^{-4}) 扫到 (10^{-1})每个 (\lambda) 都记录训练误差、验证误差和非零系数个数然后看曲线拐点。下面代码直接给出扫描和画图。lambdas logspace(-4, -1, 25); errTrain zeros(size(lambdas)); errVal zeros(size(lambdas)); nnzCoef zeros(size(lambdas)); for ii 1:numel(lambdas) XiTmp stlsq(Theta, dXtrain, lambdas(ii), maxIter); errTrain(ii) norm(Theta*XiTmp - dXtrain,fro) / norm(dXtrain,fro); errVal(ii) norm(ThetaVal*XiTmp - dXval,fro) / norm(dXval,fro); nnzCoef(ii) nnz(XiTmp); end figure; semilogx(lambdas, errTrain, o-, lambdas, errVal, s-); xlabel(\lambda); ylabel(relative error); legend(train,validation); grid on;读图时不要只挑验证误差最低点。更实用的判断是验证误差开始明显上升之前非零系数个数已经降到一个可解释的数量再增大 (\lambda)验证误差突然跳高说明删掉了关键项。这个拐点通常比最低点更稳。若曲线很平说明库列之间高度相关或者导数估计太差阈值怎么调都差不多。另一个检查是系数路径把每个 (\lambda) 下的非零系数按项名打印出来看关键项是不是连续存在还是在不同阈值下来回跳。跳变严重时先回退到嵌入参数和求导窗口不要继续调 (\lambda)。下面这张表给出常用参数的起点和调整方向实际项目里按验证误差和系数稳定性微调不要照抄。参数典型起点调整方向观察指标tau51~10 扫自相关过零库列相关性、系数稳定性d32~4逐步加验证误差是否明显下降win117~21 奇数导数平滑度、长期轨迹order32 或 3边界振荡、过平滑polyorder2先 2 后 3非零项数量、共线性lambda(10^{-2})(10^{-4})~(10^{-1}) 对数扫验证误差拐点、非零项数maxIter105~15系数是否收敛注意lambda的量级和dt强相关。dt改成 0.01 后导数整体放大 100 倍同样的lambda相对效果会变小。换采样周期后必须重扫阈值不能直接套旧值。4. 慢时间尺度分离滑动平均、降采样与快慢耦合验证4.1 慢变量不是低通滤波的同义词慢时间尺度发现最容易做成“随便平滑一下再回归”。平滑确实能压掉快变量但慢变量必须保留可辨识的动态否则 SINDy 拿到一条几乎常数的序列导数接近零系数全被阈值压掉。判断慢变量是否合格我一般看两个量自相关衰减到 (1/e) 的滞后长度以及频谱低端是否还有明显峰值。滑动窗口 (W) 至少覆盖几个快时间尺度但不能大到把慢周期也抹平。常见做法是先算自相关再把 (W) 设为自相关特征 lag 的 2 到 5 倍并且保持奇数。下面代码手写自相关避免依赖额外工具箱。maxLag不要超过序列长度的一半否则尾部样本太少自相关估计不稳。maxLag 500; xc x - mean(x); ac zeros(maxLag1,1); for L 0:maxLag ac(L1) (xc(1:end-L) * xc(1L:end)) / (xc * xc); end idx1e find(ac exp(-1), 1, first); if isempty(idx1e) tauCross maxLag; else tauCross idx1e - 1; end W round(3 * max(tauCross, 1)); W max(11, 2*floor(W/2)1);ac是归一化自相关tauCross是第一次低于 (1/e) 的滞后。W取 3 倍并强制奇数是为了让movmean窗口对称。若tauCross很小说明序列去相关很快慢结构可能不明显这时硬做慢变量分离意义不大若tauCross接近maxLag说明maxLag选短了或者系统本身有长记忆需要加长数据再判断。这个检查比拍脑袋定 (W101) 靠谱。4.2 多尺度窗口选择与重构确定 (W) 后用滑动平均得到慢变量 (s)再对 (s) 做延迟嵌入、平滑、求导、建库和 STLSQ。下面代码复用前面的poly_library和stlsq把慢变量链路串起来。s movmean(x, W, Endpoints,discard); tau_s 5; d_s 3; M_s numel(s) - (d_s-1)*tau_s; S zeros(M_s, d_s); for j 1:d_s S(:,j) s((1:M_s) (j-1)*tau_s); end Ss sgolayfilt(S, 3, 11); dS gradient(Ss, dt); [Theta_s, names_s] poly_library(S, 2); nTrain_s round(0.7*size(S,1)); Xi_s stlsq(Theta_s(1:nTrain_s,:), dS(1:nTrain_s,:), 0.02, 10);movmean的Endpointsdiscard会丢掉两端不完整窗口避免边界用缩窗平均造成假动态。W太大时s的动态幅度变小dS也变小lambda0.02可能直接把所有项压零W太小时s还带着快变量慢时间尺度没分离出来。实际调参时我会把 (W) 设成[21 51 101 201]扫一遍记录验证误差和非零项数量再选一个系数路径最稳定的窗口。不要只看一条平滑曲线“看起来顺眼”要看回归结果是否可复现。慢变量嵌入参数tau_s和d_s也不要完全照搬快变量。慢变量采样率相对变低tau_s5可能已经跨过较长一段慢动态若W很大tau_s要相应减小否则嵌入矩阵列之间关系太远。一个实用检查是看慢变量自相关tau_s不要超过慢变量自相关第一次过零的滞后。这个检查花不了几分钟但能省掉后面反复调阈值的麻烦。4.3 用一步预测和长期轨迹双重验证慢变量模型不能只看训练拟合。先做一步预测用验证段的库矩阵乘系数和验证段导数比较。再做长期轨迹从验证段某个初值出发用 SINDy 方程迭代若干步看轨迹是否发散或塌缩。下面代码给出最小验证片段。valStart nTrain_s 20; idx valStart:min(valStart199, size(S,1)-1); pred_dS Theta_s(idx,:) * Xi_s; relErr norm(pred_dS - dS(idx,:), fro) / norm(dS(idx,:), fro); z S(valStart,:); traj zeros(300, d_s); traj(1,:) z; for k 1:299 th poly_library(z, 2); z z dt * (th * Xi_s); traj(k1,:) z; endrelErr是一步预测相对误差能反映局部导数拟合。长期轨迹用显式欧拉迭代dt1与数据生成一致若换成连续采样dt要改小否则长期轨迹会数值发散。poly_library(z,2)这里传入单行向量函数内部会转成行向量并返回 (1\times p) 的库行乘Xi_s得到 (1\times d_s) 的导数。长期轨迹发散不一定是模型错也可能是初值落在训练数据覆盖不到的区域或者慢变量窗口把动态抹得太平。判断时要同时看一步预测和轨迹一步预测好、长期轨迹很快发散多半是库缺项或阈值删掉了稳定项一步预测就差先回去查导数和嵌入。5. 避坑与排查五个让系数全错的细节5.1 现象系数稀疏到全零现象很直观阈值扫描一开始nnz(Xi)很快掉到 0训练误差和验证误差都很大画出来的系数路径空空的。原因通常不是 SINDy 本身而是 (\lambda) 相对导数尺度太大或者库列量级差异太大硬阈值把本来有用的项一起压掉。Logistic 序列在 ([0,1]) 内导数幅度不大lambda0.1就可能全零换成标准化后的多变量数据lambda0.02又可能显得太小。另一个常见原因是平滑窗口太大导数被抹得接近零所有项都显得不重要。解决顺序我一般这样走先把lambda设成 0跑一次普通最小二乘看系数最大绝对值是多少再把初始 (\lambda) 设成最大系数的 (1%) 到 (5%)然后检查Theta每列的范数差异超过两个数量级就先标准化列回归完再换回原尺度最后把win本文还有配套的精品资源点击获取