参考文献

Lorenz 吸引子

使用 ode45 生成序列

%// generate Data
sigma = 10;  %// Lorenz's parameters (chaotic)
beta = 8/3;
rho = 28;
n = 3;
x0=[-8; 8; 27];  %// Initial condition%// Integrate
dt = 0.001;
tspan=[dt:dt:200];
N = length(tspan);
options = odeset('RelTol',1e-12,'AbsTol',1e-12*ones(1,n));
[t,xdat]=ode45(@(t,x) lorenz(t,x,sigma,beta,rho),tspan,x0,options);

%// Plot
figure
L = 1:200000;
plot3(xdat(L,1),xdat(L,2),xdat(L,3),'Color',[.1 .1 .1],'LineWidth',1.5)
axis on
view(-5,12)
axis tight
xlabel('x'), ylabel('y'), zlabel('z')
set(gca,'FontSize',14)
set(gcf,'Position',[100 100 600 400])
set(gcf,'PaperPositionMode','auto')

洛伦兹系统的 xxx 分量


figure
plot(tspan,xdat(:,1),'k','LineWidth',2)
xlabel('t'), ylabel('x')
set(gca,'XTick',[0 10 20 30 40 50 60 70 80 90 100],'YTick',[-20 -10 0 10 20])
set(gcf,'Position',[100 100 550 300])
xlim([0 100])
set(gcf,'PaperPositionMode','auto')

Eigen-Time Delay Coordinates

将lorenz系统的第一个分量为分析对象，通过堆叠时滞序列建立 Hankel 矩阵，然后做奇异值分解：

stackmax = 100;  %// the number of shift-stacked rows%%// EIGEN-TIME DELAY COORDINATES
clear V, clear dV, clear HH = zeros(stackmax,size(xdat,1)-stackmax);
for k=1:stackmaxH(k,:) = xdat(k:end-stackmax-1+k,1);
end[U,S,V] = svd(H,'econ');

对 V∗V^*V∗ 的前三行（奇异值最大的三个动态模式）作图

figure
L = 1:170000;
plot3(V(L,1),V(L,2),V(L,3),'Color',[.1 .1 .1],'LineWidth',1.5)
axis tight
xlabel('v_1'), ylabel('v_2'), zlabel('v_3')
set(gca,'FontSize',14)
view(34,22)
set(gcf,'Position',[100 100 600 400])
set(gcf,'PaperPositionMode','auto')

从数据计算导数

利用最优SVHT确定奇异值的截断值，从而确定 HHH 的主成分个数 rrr

sigs = diag(S);
beta = size(H,1)/size(H,2);
thresh = optimal_SVHT_coef(beta,0) * median(sigs);
r = length(sigs(sigs>thresh))
r=min(rmax,r)

利用四阶中心差分法计算 VVV 的导数 V˙\dot{V}V˙
V˙(t−2)=−V(t+2)+8∗V(t+1)−8∗V(t−1)+V(t−2)12dt\dot{V}(t-2) = \frac{-V(t+2)+8*V(t+1)-8*V(t-1)+V(t-2)}{12dt} V˙(t−2)=12dt−V(t+2)+8∗V(t+1)−8∗V(t−1)+V(t−2)

%%// COMPUTE DERIVATIVES
%// compute derivative using fourth order central difference
%// use TVRegDiff if more error
dV = zeros(length(V)-5,r);
for i=3:length(V)-3for k=1:rdV(i-2,k) = (1/(12*dt))*(-V(i+2,k)+8*V(i+1,k)-8*V(i-1,k)+V(i-2,k));end
end

完了之后令 x=V,x˙=V˙x = V, \dot{x}=\dot{V}x=V,x˙=V˙

%// concatenate
x = V(3:end-3,1:r);
dx = dV;

对非线性动态系统的稀疏辨识 (SINDY)

文献：https://www.pnas.org/content/pnas/113/15/3932.full.pdf

用一组人为选定的基函数（多项式函数、正弦函数等）对观测量进行特征扩充，X→Θ(X)X \to \Theta(X)X→Θ(X)，用如下的 poolData 函数实现：

function yout = poolData(yin,nVars,polyorder,usesine)
%// Copyright 2015, All Rights Reserved
%// Code by Steven L. Brunton
%// For Paper, "Discovering Governing Equations from Data:
%//        Sparse Identification of Nonlinear Dynamical Systems"
%// by S. L. Brunton, J. L. Proctor, and J. N. Kutzn = size(yin,1);
%// yout = zeros(n,1+nVars+(nVars*(nVars+1)/2)+(nVars*(nVars+1)*(nVars+2)/(2*3))+11);ind = 1;
%// poly order 0
yout(:,ind) = ones(n,1);
ind = ind+1;%// poly order 1
for i=1:nVarsyout(:,ind) = yin(:,i);ind = ind+1;
endif(polyorder>=2)%// poly order 2for i=1:nVarsfor j=i:nVarsyout(:,ind) = yin(:,i).*yin(:,j);ind = ind+1;endend
endif(polyorder>=3)%// poly order 3for i=1:nVarsfor j=i:nVarsfor k=j:nVarsyout(:,ind) = yin(:,i).*yin(:,j).*yin(:,k);ind = ind+1;endendend
endif(polyorder>=4)%// poly order 4for i=1:nVarsfor j=i:nVarsfor k=j:nVarsfor l=k:nVarsyout(:,ind) = yin(:,i).*yin(:,j).*yin(:,k).*yin(:,l);ind = ind+1;endendendend
end%// poly order 5, 6, 7 ...if(usesine)for k=1:10;yout = [yout sin(k*yin) cos(k*yin)];end
end

然后用稀疏回归算法求解：
V˙=Θ(V)Ξ\dot{V} = \Theta(V) \Xi V˙=Θ(V)Ξ
稀疏线性回归算法可以用 LASSO，也可以用序贯阈值最小二乘法 (sequential thresholded least-squares algorithm)，即每次最小二乘求出权重后，将低于阈值的权重强制设为0，然后用剩下的特征再做最小二乘，迭代若干次

function Xi = sparsifyDynamics(Theta,dXdt,lambda,n)
%// Copyright 2015, All Rights Reserved
%// Code by Steven L. Brunton
%// For Paper, "Discovering Governing Equations from Data:
%//        Sparse Identification of Nonlinear Dynamical Systems"
%// by S. L. Brunton, J. L. Proctor, and J. N. Kutz%// compute Sparse regression: sequential least squares
Xi = Theta\dXdt;  %// initial guess: Least-squares%// lambda is our sparsification knob.
for k=1:10smallinds = (abs(Xi)<lambda);   %// find small coefficientsXi(smallinds)=0;                %// and thresholdfor ind = 1:n                   %// n is state dimensionbiginds = ~smallinds(:,ind);%// Regress dynamics onto remaining terms to find sparse XiXi(biginds,ind) = Theta(:,biginds)\dXdt(:,ind); end
end

回归代码，计算出 Ξ\XiΞ

%%//  BUILD HAVOK REGRESSION MODEL ON TIME DELAY COORDINATES
%// This implementation uses the SINDY code, but least-squares works too
%// Build library of nonlinear time series
polyorder = 1;
Theta = poolData(x,r,1,0);
%// normalize columns of Theta (required in new time-delay coords)
for k=1:size(Theta,2)normTheta(k) = norm(Theta(:,k));Theta(:,k) = Theta(:,k)/normTheta(k);
end
m = size(Theta,2);
%// compute Sparse regression: sequential least squares
%// requires different lambda parameters for each column
clear Xi
for k=1:r-1Xi(:,k) = sparsifyDynamics(Theta,dx(:,k),lambda*k,1);  %// lambda = 0 gives better results
end
Theta = poolData(x,r,1,0);
for k=1:length(Xi)Xi(k,:) = Xi(k,:)/normTheta(k);
end

由于设置 polyorder = 1，相当于：
V˙=[1;V]Ξ\dot{V} = [\mathbf{1};V]\Xi V˙=[1;V]Ξ
但在求出 Ξ\XiΞ 之后，并不是建立 VVV 的所有分量的线性模型，而是将 VVV 的第 rrr 项（SVD截断后的最后一项）最为外部扰动项，而构建关于前 r−1r-1r−1 个分量的线性模型

即：
v˙=Av+Bvr\dot{v} = Av + Bv_r v˙=Av+Bvr其中 v=[v1…vr−1]Tv=\begin{bmatrix} v_1& \ldots&v_{r-1}\end{bmatrix}^Tv=[v1…vr−1]T

A = Xi(2:r+1,1:r-1)';  %// 第一项为常数项对应的稀疏，所以从2开始
B = A(:,r);
A = A(:,1:r-1);
%
L = 1:50000;
sys = ss(A,B,eye(r-1),0*B);
[y,t] = lsim(sys,x(L,r),dt*(L-1),x(1,1:r-1));

%%//  Part 4:  Model Time Series
L = 300:25000;
L2 = 300:50:25000;
figure
subplot(2,1,1)
plot(tspan(L),x(L,1),'Color',[.4 .4 .4],'LineWidth',2.5)
hold on
plot(tspan(L2),y(L2,1),'.','Color',[0 0 .5],'LineWidth',5,'MarkerSize',15)
xlim([0 max(tspan(L))])
ylim([-.0051 .005])
ylabel('v_1')
box on
subplot(2,1,2)
plot(tspan(L),x(L,r),'Color',[.5 0 0],'LineWidth',1.5)
xlim([0 max(tspan(L))])
ylim([-.025 .024])
xlabel('t'), ylabel('v_{15}')
box on
set(gcf,'Position',[100 100 550 350])
set(gcf,'PaperPositionMode','auto')

Hankel alternative view of Koopman (HAVOK) analysis相关推荐

A Collection of 100+ Writing Task 2 Essays for IELTS
EDITION 2019 A Collection of 100+ Writing Task 2 Essays IELTS ESSAYS FROM EXAMINERS VERSION 3.0 OREM ...
pca 主成分分析_超越普通PCA：非线性主成分分析
pca 主成分分析 TL;DR: PCA cannot handle categorical variables because it makes linear assumptions about t ...
移动端分步注册_移动应用程序的可用性测试：分步指南
移动端分步注册 Written by Justin Mifsud 由贾斯汀·米夫苏德 ( Justin Mifsud)撰写 The mobile market is huge and growing ...
大数据多维分析常用操作图解 OLAP Operations
多维数据模型中的 OLAP 操作 OLAP Operations in the Multidimensional Data Model 在多维模型中,记录被组织成不同的维度,每个维度包括由概念层次结构 ...
史上R语言最强--资源（免费课程、书籍、教程和各种高级图表）
史上R语言最强–资源(免费课程.书籍.教程和各种高级图表) 下面是一些免费R的书籍.教程.软件包.图表2单和其他材料,可以在网上学习编程和改进工作流程.有单独的概述 Python resources, ...
ICML 2018 paper(oral)
参考链接 icml 2018 oral Paperlist Optimal Tuning for Divide-and-conquer Kernel Ridge Regression with Mas ...
ai人工智能市场客户_投资管理中的人工智能可提升客户关系和回报
ai人工智能市场客户 Let's be honest. An investment manager's clients probably won't care about the fancy AI t ...
pytorch 查看模型参数，查看模型特定层输入输出，模型结构图绘制总结
1 参考链接大杂烩 https://zhuanlan.zhihu.com/p/33992733 绘制图形时候记得安装graphviz插件 https://stackoverflow.com/qu ...
Kinect+OpenNI学习笔记之12(简单手势所表示的数字的识别)
引自:http://www.cnblogs.com/tornadomeet/archive/2012/11/04/2753185.html 前言这篇文章是本人玩kinect时做的一个小实验,即不采用 ...

Hankel alternative view of Koopman (HAVOK) analysis

文章目录