当前位置：首页 > news >正文

Hankel alternative view of Koopman (HAVOK) analysis

news 来源：原创 2024/5/18 13:19:54

在这里插入图片描述

文章目录

参考文献
Lorenz 吸引子
Eigen-Time Delay Coordinates
从数据计算导数
对非线性动态系统的稀疏辨识 (SINDY)

参考文献

在这里插入图片描述

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')

在这里插入图片描述

洛伦兹系统的 $x$ 分量


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 H

H = zeros(stackmax,size(xdat,1)-stackmax);
for k=1:stackmax
    H(k,:) = xdat(k:end-stackmax-1+k,1);
end

[U,S,V] = svd(H,'econ');

对 $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确定奇异值的截断值，从而确定 $H$ 的主成分个数 $r$

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)

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

%%// 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)-3
    for k=1:r
        dV(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

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

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

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

文献：https://www.pnas.org/content/pnas/113/15/3932.full.pdf
在这里插入图片描述
用一组人为选定的基函数（多项式函数、正弦函数等）对观测量进行特征扩充， $\to \Theta(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. Kutz

n = 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:nVars
    yout(:,ind) = yin(:,i);
    ind = ind+1;
end

if(polyorder>=2)
    %// poly order 2
    for i=1:nVars
        for j=i:nVars
            yout(:,ind) = yin(:,i).*yin(:,j);
            ind = ind+1;
        end
    end
end

if(polyorder>=3)
    %// poly order 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

if(polyorder>=4)
    %// poly order 4
    for i=1:nVars
        for j=i:nVars
            for k=j:nVars
                for l=k:nVars
                    yout(:,ind) = yin(:,i).*yin(:,j).*yin(:,k).*yin(:,l);
                    ind = ind+1;
                end
            end
        end
    end
end

%// poly order 5, 6, 7 ...

if(usesine)
    for k=1:10;
        yout = [yout sin(k*yin) cos(k*yin)];
    end
end

然后用稀疏回归算法求解：
$\dot{V} = \Theta(V) \Xi$
稀疏线性回归算法可以用 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:10
    smallinds = (abs(Xi)<lambda);   %// find small coefficients
    Xi(smallinds)=0;                %// and threshold
    for ind = 1:n                   %// n is state dimension
        biginds = ~smallinds(:,ind);
        %// Regress dynamics onto remaining terms to find sparse Xi
        Xi(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-1
    Xi(:,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，相当于：
$\dot{V} = [\mathbf{1};V]\Xi$
但在求出 $\Xi$ 之后，并不是建立 $V$ 的所有分量的线性模型，而是将 $V$ 的第 $r$ 项（SVD截断后的最后一项）最为外部扰动项，而构建关于前 $r - 1$ 个分量的线性模型

在这里插入图片描述
即：
$\dot{v} = Av + Bv_r$ 其中 $v=\begin{bmatrix} v_1& \ldots&v_{r-1}\end{bmatrix}^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')