第9章 智能非线性动力学
此专栏为《非线性动力学及其工程应用》(田瑞兰教授)辅助文档
专栏代码使用matlab编写，感兴趣的同学可以在xx网址下载学习
1.SINDy
在科学和工程领域，我们经常遇到需要从观测数据中推断出控制系统的动力学方程的问题。传统方法通常需要基于物理定律手动推导方程，但这种方法在面对复杂系统时往往效率低下。近年来，一种名为稀疏识别非线性动力学(Sparse Identification of Nonlinear Dynamics, SINDy)的方法应运而生，它能够直接从数据中自动识别出控制系统的非线性动力学方程。
SINDy方法原理
SINDy方法由Steven L. Brunton等人于2016年提出，其核心思想是利用稀疏回归技术从候选函数库中识别出最能描述观测数据动态的最简方程。
基本数学框架
考虑一个动力学系统：
其中\(X \in \mathbb{R}^n\)是状态变量，\(f(X)\)是未知的非线性函数。
SINDy方法假设\(f(X)\)可以表示为一系列基函数的线性组合：
其中\(\Theta(X)\)是包含候选非线性项(如多项式、三角函数等)的库矩阵，\(\Xi\)是稀疏系数矩阵。
实现步骤
1.数据收集：测量或模拟系统状态\(x(t)\)及其导数\(dx/dt\)(或近似计算)
2.构建库矩阵：创建包含候选非线性项的\(\Theta(X)\)
3.稀疏回归：求解\(\Xi\)使得\(dx/dt \approx \Theta(X) \Xi\)，同时保持\(\Xi\)稀疏
4.模型选择：通过阈值处理或交叉验证选择最重要的项
关键优势
-
稀疏性：只保留对系统动态有显著贡献的项
-
可解释性：结果是可以解释的微分方程而非黑箱模型
-
数据高效：相比神经网络等方法需要更少数据
下面我们通过MATLAB代码演示如何使用SINDy方法从数据中识别著名的洛伦兹系统。
% SINDy_SingleFile_Lorenz.m
% 完整单文件实现：使用SINDy方法识别洛伦兹系统
clear all; close all; clc;
%% 生成洛伦兹系统数据
sigma = 10; beta = 8/3; rho = 28; % 经典混沌参数
n = 3; % 状态变量维度
x0 = [-8; 8; 27]; % 初始条件
dt = 0.001; % 时间步长
tspan = dt:dt:50; % 时间范围
options = odeset('RelTol',1e-12,'AbsTol',1e-12*ones(1,n));
[t,x] = ode45(@(t,x) lorenz(t,x,sigma,beta,rho), tspan, x0, options);
%% 计算导数(使用总变分正则化数值微分)
dx = TVRegDiffNum(x, 10, 0.02, [], 'small', 1e12, dt, 1, 1);
%% 构建库矩阵(Theta)
polyorder = 3; % 多项式最高阶数
Theta = poolData(x, n, polyorder); % 包含常数、线性和非线性项
%% 稀疏回归: sequential thresholded least squares
lambda = 0.025; % 稀疏化参数
Xi = sparsifyDynamics(Theta, dx, lambda, n);
%% 输出结果
yout = poolDataLIST({'x','y','z'}, Xi, n, polyorder);
disp('识别结果:');
disp(yout);
%% 真实参数和识别参数对比
true_params = [sigma; rho; beta];
identified_params = [Xi(2,1); Xi(8,1); -Xi(10,1)]; % 根据库顺序调整
disp('真实参数 [sigma; rho; beta]:');
disp(true_params);
disp('识别参数:');
disp(identified_params);
%% 可视化
figure;
subplot(3,1,1);
plot(t,x(:,1),'k','LineWidth',1.5); ylabel('x');
subplot(3,1,2);
plot(t,x(:,2),'k','LineWidth',1.5); ylabel('y');
subplot(3,1,3);
plot(t,x(:,3),'k','LineWidth',1.5); ylabel('z'); xlabel('Time');
figure;
plot3(x(:,1),x(:,2),x(:,3),'k','LineWidth',1);
xlabel('x'); ylabel('y'); zlabel('z'); title('Lorenz吸引子');
%% 洛伦兹系统方程
function dx = lorenz(t,x,sigma,beta,rho)
dx = [
sigma*(x(2)-x(1));
x(1)*(rho-x(3))-x(2);
x(1)*x(2)-beta*x(3);
];
end
%% 总变分正则化数值微分 (TVRegDiffNum)
function dx = TVRegDiffNum(y, iter, alph, u0, scale, ep, dt, plotflag, diagflag)
% ... [参数检查部分]
[n, m] = size(y);
dx = zeros(n, m); % 保持原始尺寸
for i = 1:m
% 使用中心差分（保持输出尺寸不变）
D = spdiags([-ones(n,1) ones(n,1)], [-1 1], n, n)/(2*dt);
D(1,1:2) = [-1 1]/dt; % 前向差分在边界
D(end,end-1:end) = [-1 1]/dt; % 后向差分在边界
% ... [其余代码保持不变]
end
end
%% 构建候选函数库 (poolData)
function yout = poolData(yin, nVars, polyorder)
n = size(yin,1);
ind = 1;
% 常数项
yout(:,ind) = ones(n,1);
ind = ind+1;
% 1阶项
for i=1:nVars
yout(:,ind) = yin(:,i);
ind = ind+1;
end
if polyorder >= 2
% 2阶项
for i=1:nVars
for j=i:nVars
yout(:,ind) = yin(:,i).*yin(:,j);
ind = ind+1;
end
end
end
if polyorder >= 3
% 3阶项
for i=1:nVars
for j=i:nVars
for k=j:nVars
yout(:,ind) = yin(:,i).*yin(:,j).*yin(:,k);
ind = ind+1;
end
end
end
end
end
%% 稀疏回归 (sparsifyDynamics)
function Xi = sparsifyDynamics(Theta, dXdt, lambda, n)
% 使用序列阈值最小二乘法
Xi = Theta\dXdt; % 初始最小二乘解
% 小系数阈值处理
for k=1:10
smallinds = (abs(Xi)<lambda); % 找出小系数
Xi(smallinds)=0; % 设为0
% 对非零系数重新拟合
for ind = 1:n
biginds = ~smallinds(:,ind);
Xi(biginds,ind) = Theta(:,biginds)\dXdt(:,ind);
end
end
end
%% 结果显示 (poolDataLIST)
function eqns = poolDataLIST(vars, Xi, n, polyorder)
% 将系数矩阵转换为可读方程
eqns = cell(n,1);
for i=1:n
eqn = ['d' vars{i} '/dt = '];
firstTerm = true;
% 常数项
if Xi(1,i)~=0
eqn = [eqn num2str(Xi(1,i))];
firstTerm = false;
end
% 1阶项
for j=2:n+1
if Xi(j,i)~=0
if ~firstTerm
if Xi(j,i)>0
eqn = [eqn ' + '];
else
eqn = [eqn ' - '];
end
elseif Xi(j,i)<0
eqn = [eqn '-'];
end
eqn = [eqn num2str(abs(Xi(j,i))) '*' vars{j-1}];
firstTerm = false;
end
end
% 高阶项
if polyorder>=2
terms = n+2;
% 2阶项
for j=1:n
for k=j:n
if Xi(terms,i)~=0
if ~firstTerm
if Xi(terms,i)>0
eqn = [eqn ' + '];
else
eqn = [eqn ' - '];
end
elseif Xi(terms,i)<0
eqn = [eqn '-'];
end
eqn = [eqn num2str(abs(Xi(terms,i))) '*' vars{j} '*' vars{k}];
firstTerm = false;
end
terms = terms+1;
end
end
% 3阶项
if polyorder>=3
for j=1:n
for k=j:n
for l=k:n
if Xi(terms,i)~=0
if ~firstTerm
if Xi(terms,i)>0
eqn = [eqn ' + '];
else
eqn = [eqn ' - '];
end
elseif Xi(terms,i)<0
eqn = [eqn '-'];
end
eqn = [eqn num2str(abs(Xi(terms,i))) '*' vars{j} '*' vars{k} '*' vars{l}];
firstTerm = false;
end
terms = terms+1;
end
end
end
end
end
if firstTerm
eqn = [eqn '0'];
end
eqns{i} = eqn;
end
end
2 物理信息神经网络（PINN）
物理信息神经网络（Physics-Informed Neural Networks, PINN）是一种将物理定律融入深度学习模型的新型方法。它通过在损失函数中加入物理方程约束，使神经网络不仅能拟合数据，还能遵守已知的物理规律。
PINN的核心思想
PINN的核心在于将两类信息结合：
1.观测数据（如有）
2.控制问题的物理方程（如PDEs）
通过构建一个复合损失函数来同时最小化数据误差和物理方程残差。
PINN的优势
-
数据高效：在数据稀缺的情况下仍能表现良好
-
泛化能力强：遵守物理规律使其在训练域外也有较好表现
-
多物理场耦合：天然适合处理多物理场耦合问题
-
逆问题求解：可以同时求解正问题和参数反演
参考原论文案例：
浙公网安备 33010602011771号