【机器学习】基于稀疏识别方法的洛伦兹混沌系统预测

1. 引言

1.1. DNN模型的来由

从数据中识别非线性动态学意味着什么?
假设我们有时间序列数据,这些数据来自一个(非线性)动态学系统。

  • 识别一个系统意味着基于数据推断该系统的控制方程。换句话说,就是找到动态系统方程 x ˙ = f ( x ) \mathbf{\dot x} = f(\mathbf{x}) x˙=f(x)中的 f f f(其中 x \mathbf{x} x 可能是向量值)。
    在这里插入图片描述
    例如,对于Lorenz系统,我们希望从时间序列数据中学习到方程右边的部分。

1.2. 研究稀疏性的意义

为什么我们需要稀疏性?
在这里,稀疏性意味着控制方程中的项数很少。稀疏性是有益的,因为它更具:

  • 可解释性。在需要理解变量及其相互作用的应用中至关重要,例如在需要安全关键保证的应用中。
  • 泛化能力。如果正确,方程将准确描述训练数据所填充的状态空间区域之外的动态。
    通常,人们可以将SINDy识别的模型视为与物理方程相对的模型,而不是大型、不透明的深度神经网络。

2. SINDy算法

SINDy试图找到适合数据 X ˙ = f ( X ) \mathrm{\dot X} = f(\mathrm{X}) X˙=f(X)的动态系统 f f f。这个函数逼近问题被表述为线性回归 X ˙ = Θ ( X ) Ξ \mathrm{\dot X} = \Theta(\mathrm{X}) \Xi X˙=Θ(X)Ξ,其中系数为 Ξ \Xi Ξ 和回归项库 Θ ( X ) \Theta(X) Θ(X)。算法分为三个步骤:

  • 从动态系统生成数据 X X X 并计算导数 X ˙ \dot X X˙
  • 建立候选项库 Θ ( X ) \Theta(X) Θ(X) 作为 X X X上的函数。
  • 稀疏回归系数 Ξ \Xi Ξ,以最好地描述数据。
  1. SINDy假设是测量了 n n n 维数据点的时间序列 x = ( x 1 , … x n ) \mathbf{x}=(x_1, \ldots x_n) x=(x1,xn) m m m 个时间步 t 1 , … , t m t_1, \ldots, t_m t1,,tm,我们定义数据矩阵 X X X 和导数矩阵 X ˙ \dot X X˙ 为:
    X = [ x 1 ( t 1 ) x 2 ( t 1 ) ⋯ x n ( t 1 ) x 1 ( t 2 ) x 2 ( t 2 ) ⋯ x n ( t 2 ) x 1 ( t 3 ) x 2 ( t 3 ) ⋯ x n ( t 3 ) ⋮ ⋮ ⋱ ⋮ x 1 ( t m ) x 2 ( t m ) ⋯ x n ( t m ) ] X=\begin{bmatrix} x_1(t_1)&x_2(t_1)&\cdots&x_n(t_1)\\ x_1(t_2)&x_2(t_2)&\cdots&x_n(t_2)\\ x_1(t_3)&x_2(t_3)&\cdots&x_n(t_3)\\ \vdots&\vdots&\ddots&\vdots\\ x_1(t_m)&x_2(t_m)&\cdots&x_n(t_m)\\ \end{bmatrix} X= x1(t1)x1(t2)x1(t3)x1(tm)x2(t1)x2(t2)x2(t3)x2(tm)xn(t1)xn(t2)xn(t3)xn(tm)
    X ˙ = [ x ˙ 1 ( t 1 ) x ˙ 2 ( t 1 ) ⋯ x ˙ n ( t 1 ) x ˙ 1 ( t 2 ) x ˙ 2 ( t 2 ) ⋯ x ˙ n ( t 2 ) x ˙ 1 ( t 3 ) x ˙ 2 ( t 3 ) ⋯ x ˙ n ( t 3 ) ⋮ ⋮ ⋱ ⋮ x ˙ 1 ( t m ) x ˙ 2 ( t m ) ⋯ x ˙ n ( t m ) ] \dot{X}=\begin{bmatrix} \dot{x}_1(t_1)&\dot{x}_2(t_1)&\cdots&\dot{x}_n(t_1)\\ \dot{x}_1(t_2)&\dot{x}_2(t_2)&\cdots&\dot{x}_n(t_2)\\ \dot{x}_1(t_3)&\dot{x}_2(t_3)&\cdots&\dot{x}_n(t_3)\\ \vdots&\vdots&\ddots&\vdots\\ \dot{x}_1(t_m)&\dot{x}_2(t_m)&\cdots&\dot{x}_n(t_m)\\ \end{bmatrix} X˙= x˙1(t1)x˙1(t2)x˙1(t3)x˙1(tm)x˙2(t1)x˙2(t2)x˙2(t3)x˙2(tm)x˙n(t1)x˙n(t2)x˙n(t3)x˙n(tm)
  2. 接下来,我们定义库矩阵 Θ ( X ) \Theta(X) Θ(X),其列是应用于数据的一组基函数 { θ l } l = 1 , … , L \{\theta_l\}_{l=1,\ldots, L} {θl}l=1,,L
    Θ ( X ) = [ ∣ ∣ ∣ Θ 1 ( X ) Θ 2 ( X ) ⋯ Θ n ( X ) ∣ ∣ ∣ ] \Theta(X)=\begin{bmatrix} |&|&&|\\ \Theta_1(X)&\Theta_2(X)&\cdots&\Theta_n(X)\\ |&|&&|\\ \end{bmatrix} Θ(X)= Θ1(X)Θ2(X)Θn(X)
    以下是对上述文本的改写:
    一些简单的例子包括多项式的基础,比如在泰勒级数展开中使用的多项式,或者三角函数,如 sin(x_1), cos(x_1), sin(2x_1), 等等,这些在傅里叶级数展开中常见。然而,根据不同的问题,可能需要使用更复杂的基础函数,例如贝塞尔函数。
  • 最后,我们使用稀疏线性回归算法(例如LASSO)找到系数 Ξ \Xi Ξ,使得:
    Ξ = [ ∣ ∣ ∣ ε 1 ( X ) ε 2 ( X ) ⋯ ε n ( X ) ∣ ∣ ∣ ] \Xi=\begin{bmatrix} |&|&&|\\ \varepsilon_1(X)&\varepsilon_2(X)&\cdots&\varepsilon_n(X)\\ |&|&&|\\ \end{bmatrix} Ξ= ε1(X)ε2(X)εn(X)
    因此
    X ˙ = θ ( X ) Ξ \dot{X}=\theta(X)\Xi X˙=θ(X)Ξ
    在这里插入图片描述

3.线性动态系统示例

假设我们测量了一个由以下动力学系统控制的粒子轨迹:

d d t ( x y ) = ( − 2 x y ) = ( − 2 0 0 1 ) ( x y ) \frac{d}{dt} \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} -2x \\ y \end{pmatrix} = \begin{pmatrix} -2 & 0 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} dtd(xy)=(2xy)=(2001)(xy)

3.1.构建数据矩阵

以初始条件 x 0 = 3 x_0=3 x0=3 y 0 = 1 3 y_0=\frac{1}{3} y0=31我们构建数据矩阵 X X X

import numpy as np
import pysindy as pst = np.linspace(0, 1, 100)
x = 3 * np.exp(-2 * t)
y = 0.5 * np.exp(t)
X = np.stack((x, y), axis=-1)

3.2.拟合模型

在构建SINDy模型时,我们有几个关键选项需要确定,包括微分方法、基函数库以及优化器的选择。现在,让我们具体探讨一下以傅里叶基作为基函数的情况。

model = ps.SINDy(differentiation_method=ps.FiniteDifference(order=2),feature_library=ps.FourierLibrary(),optimizer=ps.STLSQ(threshold=0.2),feature_names=["x", "y"]
)
model.fit(X, t=t)

3.3.测试模型

拟合模型后,可以使用 print 成员函数检查模型。

model.print()
(x)' = 0.772 sin(1 x) + 2.097 cos(1 x) + -2.298 sin(1 y) + -3.115 cos(1 y)
(y)' = 1.362 sin(1 y) + -0.222 cos(1 y)

然后使用测试数据验证

import matplotlib.pyplot as pltdef plot_simulation(model, x0, y0):t_test = np.linspace(0, 1, 100)x_test = x0 * np.exp(-2 * t_test)y_test = y0 * np.exp(t_test)sim = model.simulate([x0, y0], t=t_test)plt.figure(figsize=(6, 4))plt.plot(x_test, y_test, label="Ground truth", linewidth=4)plt.plot(sim[:, 0], sim[:, 1], "--", label="SINDy estimate", linewidth=3)plt.plot(x0, y0, "ko", label="Initial condition", markersize=8)plt.xlabel("x")plt.ylabel("y")plt.legend()plt.show()x0 = 6
y0 = -0.1
plot_simulation(model, x0, y0)
x0 = 6
y0 = -0.1
plot_simulation(model, x0,y0)

在这里插入图片描述
正如我们预期的那样,对于这个问题,傅里叶基并不是最佳选择。现在,让我们尝试一个不同的基函数!

3.4.自定义基函数

实验1:研究当我们选择一个不同的基函数时,会有什么样的结果。

为了实现这个目标,我们将编写一个Sindy算法类,并使用一个自定义的基函数库,而不是默认的ps.PolynomialLibrary()

我们将这样初始化ps.SINDy():将feature_library属性设置为我们的自定义傅里叶库(如果可用)。然后,我们将使用.fit(X, t=t)方法来拟合这个实例化的算法。

# 导入所需的库(假设 ps 是 PySINDy 的别名)  
# 注意:PySINDy 库的实际导入可能需要根据你的环境进行调整  
from pysindy import SINDy, FiniteDifference, PolynomialLibrary, STLSQ  
from pysindy.utils import plot_simulation  # 初始化一个 SINDy 模型,使用有限差分法(二阶)进行微分  
# 这里使用了多项式库(一阶)作为特征库  
# 使用了阈值为 0.2 的 STLSQ 优化器  
# 特征名称设置为 'x' 和 'y'  
model_1 = SINDy(  differentiation_method=FiniteDifference(order=2),  # 使用二阶有限差分法进行微分  feature_library=PolynomialLibrary(degree=1),       # 使用一阶多项式库作为特征库  optimizer=STLSQ(threshold=0.2),                    # 使用阈值为 0.2 的 STLSQ 优化器  feature_names=["x", "y"]                           # 特征名称设置为 'x' 和 'y'  
)  # 使用数据 X 和时间 t 拟合模型  
model_1.fit(X, t=t)  # 打印模型结果  
model_1.print()  # 假设 plot_simulation 是一个自定义函数,用于绘制模拟结果  
# 这里设置初始条件 x0=6, y0=-0.1  
x0 = 6  
y0 = -0.1  # 调用 plot_simulation 函数绘制模拟结果  
plot_simulation(model_1, x0, y0)

在这里插入图片描述
在此简单示例中,模型利用多项式基函数准确地识别了背后的数学方程。

我们面对的是一个一阶多项式微分方程,因此,使用与这些模型假设完全一致的回归方法来识别它是合情合理的。

实验1.1:尝试逐步提高多项式的度数,观察在哪个度数时SINDy无法正确识别方程,以及在哪个度数时测试数据集上的预测开始出现偏差。这将帮助我们了解SINDy在不同复杂度下的效能表现。

# 导入 PySINDy 相关的库(假设 ps 是 PySINDy 的别名)  
from pysindy import SINDy, FiniteDifference, PolynomialLibrary, STLSQ  
from pysindy.utils import plot_simulation  # 假设 plot_simulation 函数已经定义或可用  # 初始化一个 SINDy 模型  
# 使用二阶有限差分法进行微分  
# 使用四阶多项式库作为特征库  
# 使用阈值为 0.2 的 STLSQ 优化器  
# 特征名称设置为 'x' 和 'y'  
model_1 = ps.SINDy(  differentiation_method=ps.FiniteDifference(order=2),  # 使用二阶有限差分法进行微分  feature_library=ps.PolynomialLibrary(degree=4),       # 使用四阶多项式库作为特征库  optimizer=ps.STLSQ(threshold=0.2),                    # 使用阈值为 0.2 的 STLSQ 优化器  feature_names=["x", "y"]                              # 特征名称设置为 'x' 和 'y'  
)  # 使用数据 X 和时间 t 拟合模型  
model_1.fit(X, t=t)  # (注意:以下部分是错误的,因为不应该再次初始化 SINDy,并且没有缩进)  
# 正确的方式是上面已经初始化和拟合了 model_1,接下来直接使用它  
# SINDy(differentiation_method=FiniteDifference(),  
#       feature_library=PolynomialLibrary(degree=4), feature_names=['x', 'y'],  
#       optimizer=STLSQ(threshold=0.2))  # 打印模型结果  
model_1.print()  # 设置初始条件 x0=6, y0=-0.1  
x0 = 6  
y0 = -0.1  # 调用 plot_simulation 函数绘制模拟结果  
plot_simulation(model_1, x0, y0)

在这里插入图片描述
实验1.2:探讨阈值设定对模型性能的影响。如果阈值设定过低,将导致哪些问题?相反,如果阈值设定过高,又将引发哪些后果?

STLSQ(序列阈值最小二乘法):这种方法通过迭代进行最小二乘优化,并对于权重向量中小于特定阈值的元素进行置零处理,以此来最小化目标函数 ∣ ∣ y − X w ∣ ∣ 2 2 + α ∣ ∣ w ∣ ∣ 2 2 ||y-X_w||^2_2+\alpha||w||^2_2 ∣∣yXw22+α∣∣w22

阈值:这是权重向量中系数的最小允许幅度。任何幅度小于这个阈值的系数都将被清零,从而有助于实现模型的稀疏性。

# 导入 PySINDy 库(假设 ps 是 PySINDy 的别名)  
from pysindy import SINDy, FiniteDifference, PolynomialLibrary, STLSQ  
# 假设 plot_simulation 函数已经定义或可用  
from pysindy.utils import plot_simulation  # 初始化 SINDy 模型  
# 使用二阶有限差分法进行微分  
# 使用四阶多项式库作为特征库  
# 使用阈值为 0.1 的 STLSQ 优化器  
# 特征名称设置为 'x' 和 'y'  
model_1 = ps.SINDy(  differentiation_method=ps.FiniteDifference(order=2),  # 微分方法:二阶有限差分法  feature_library=ps.PolynomialLibrary(degree=4),       # 特征库:四阶多项式库  optimizer=ps.STLSQ(threshold=0.1),                    # 优化器:STLSQ,阈值设置为 0.1  feature_names=["x", "y"]                              # 特征名称:['x', 'y']  
)  # 使用数据 X 和时间 t 拟合模型  
model_1.fit(X, t=t)  # 注意:以下部分是错误的,它试图再次初始化一个 SINDy 模型,但并未赋值给任何变量  
# 这里我们忽略它,因为 model_1 已经正确初始化和拟合  
# SINDy(differentiation_method=FiniteDifference(),  
#       feature_library=PolynomialLibrary(degree=4), feature_names=['x', 'y'],  
#       optimizer=STLSQ())  # 打印模型结果  
model_1.print()  # 设置初始条件  
x0 = 6  
y0 = -0.1  # 调用 plot_simulation 函数绘制模拟结果  
plot_simulation(model_1, x0, y0)

在这里插入图片描述

4. 洛伦兹吸引子实验

当使用 SINDy 对 洛伦兹吸引子(Lorenz Attractor)进行测试时,我们首先通过其真实的动力学方程来模拟一条轨迹。这样,我们可以使用模拟的数据来验证 SINDy 方法是否能够有效地识别出 Lorenz 系统的动态结构。

import numpy as np  
import matplotlib.pyplot as plt  
from scipy.integrate import odeint  
from mpl_toolkits.mplot3d import Axes3D  
# 假设 ps 是 PySINDy 的别名  
from pysindy import SINDy, STLSQ, PolynomialLibrary  # Lorenz 吸引子的参数  
rho = 28.0  
sigma = 10.0  
beta = 8.0 / 3.0  
dt = 0.01  # Lorenz 系统的动力学方程  
def f(state, t):  x, y, z = state  return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]  # 初始状态  
state0 = [1.0, 1.0, 1.0]  
# 时间步长数组  
time_steps = np.arange(0.0, 40.0, dt)  # 模拟 Lorenz 吸引子的轨迹  
x_train = odeint(f, state0, time_steps)  # 初始化 SINDy 模型  
model = SINDy(  optimizer=STLSQ(threshold=0.05),  # 使用阈值为 0.05 的 STLSQ 优化器  feature_library=PolynomialLibrary(degree=2),  # 使用二阶多项式库  
)  # 使用模拟的数据拟合 SINDy 模型  
model.fit(x_train, t=dt)  # 使用 SINDy 模型进行模拟  
x_sim = model.simulate(x_train[0], time_steps)  # 打印 SINDy 模型的识别结果  
model.print()  # 绘制 Lorenz 吸引子的轨迹  
plt.figure(figsize=(6, 4))  
# 绘制真实轨迹  
plt.plot(x_train[:, 0], x_train[:, 2], label='真实轨迹')  
# 绘制 SINDy 估计的轨迹  
plt.plot(x_sim[:, 0], x_sim[:, 2], '--', label='SINDy 估计轨迹')  
# 绘制初始条件  
plt.plot(x_train[0, 0], x_train[0, 2], "ko", label="初始条件", markersize=8)  
plt.legend()  
plt.xlabel('x')  
plt.ylabel('z')  
plt.title('Lorenz 吸引子轨迹')  
plt.draw()  
plt.show()

在这里插入图片描述

import matplotlib.pyplot as plt# 定义一个函数,用于绘制特定维度的时间序列数据
def plot_dimension(dim, name):# 创建一个图形对象,设置图形大小为9x2英寸fig = plt.figure(figsize=(9, 2))# 获取当前的坐标轴对象ax = fig.gca()# 在坐标轴上绘制训练数据的时间序列ax.plot(time_steps, x_train[:, dim])# 在相同的坐标轴上绘制模拟数据的时间序列,使用虚线表示ax.plot(time_steps, x_sim[:, dim], "--")# 设置x轴标签为"时间"plt.xlabel("时间")# 设置y轴标签为维度名称plt.ylabel(name)# 调用函数绘制x、y、z三个维度的时间序列数据
plot_dimension(0, 'x')  # 绘制x维度
plot_dimension(1, 'y')  # 绘制y维度
plot_dimension(2, 'z')  # 绘制z维度

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

5. 挑战

将SINDy应用于真实世界数据更具挑战性,主要体现在以下几个方面:

  1. 数据:需要多少数据以及数据的质量如何?SINDy默认依赖于干净、快速采样且噪声小的数据。数据需要足够干净,以便计算导数。但这些要求可以相互权衡:例如,如果数据干净,少量即可。反之,对于大量但嘈杂的数据,可以通过积分来平均噪声。
  2. 坐标:应该测量哪些变量?如果可能,使用专家领域知识。或者应用奇异值分解来找到最相关的变量。或者将深度自编码器与SINDy结合训练。
  3. :哪些库项最能捕捉动态?实际上,从简单的开始,例如线性项,然后是二次项等。也可以使用领域知识。注意,不要过度扩展库,以避免列的线性依赖导致 (\Theta(X)) 矩阵病态。
  4. 优化:使用哪种算法?如何实际找到稀疏项?最简单的是LASSO,实践中,顺序阈值最小二乘(迭代硬阈值系数)效果更好。在这里,也可以纳入物理信息约束:例如,可以直接在优化中实施能量守恒或全局稳定性。

在过去的几年中,这方面取得了很多进展,SINDy已成功应用于各种真实世界问题,包括发现由于问题复杂性而无法用解析方法表示的等离子体流体动力学方程。

附录 A:水库计算实例

水库计算是一种简单的方法,用于在没有反向传播通过时间的情况下训练递归神经网络,以及与之相关的臭名昭著的梯度消失和爆炸问题。基本步骤如下:

  1. 随机初始化递归神经网络权重。
  2. 固定隐藏连接权重。
  3. 使用线性回归训练线性输出层。

更正式地说,递归神经网络的状态由其隐藏激活 h ∈ R M \mathbf{h} \in \mathbb{R}^{M} hRM给出,这些激活通过随机初始化并固定的矩阵 W h \mathbf{W}_{\mathrm{h}} Wh 连接。输入序列 X 1 : τ \mathbf{X}_{1:\tau} X1:τ 通过线性映射 W i \mathbf{W}_{\mathrm{i}} Wi嵌入到状态 h τ \mathbf{h}_{\tau} hτ
h 1 = θ ( W h h 0 + W i X 1 ) ⋮ h τ = θ ( W h h τ − 1 + W i X τ ) h_1=\theta(W_hh_0+W_iX_1)\\ \vdots\\ h_{\tau}=\theta(W_hh_{\tau-1}+W_iX_{\tau}) h1=θ(Whh0+WiX1)hτ=θ(Whhτ1+WiXτ)

然后,水库可以基于隐藏状态 h τ \mathbf{h}_{\tau} hτ线性预测下一个输入 x τ + 1 \mathbf{x}_{\tau +1} xτ+1
x ^ τ + 1 = W o h τ \mathbf{\hat {x}}_{\tau +1} = \mathbf{W}_{\mathrm{o}} \mathbf{h}_\tau x^τ+1=Wohτ
其中 W o \mathbf{W}_{\mathrm{o}} Wo 将隐藏激活映射到输出。只有这个输出矩阵 W o \mathbf{W}_{\mathrm{o}} Wo的参数被训练。因此,可以通过岭回归(正则化参数 λ \lambda λ)以解析方式获得最优值:
W o = ( H T H + λ I ) − 1 H T X W_{\mathrm{o}}=(H^TH+\lambda\mathrm{I})^{-1}H^TX Wo=(HTH+λI)1HTX
其中 H \mathbf{H} H 是隐藏状态, I \mathbf{I} I是单位矩阵, X \mathbf{X} X是回归的目标。

A.1. 设置

# 导入scipy库中的sparse模块  
from scipy import sparse  # 定义一些参数  
# 半径  
radius = 0.6  
# 稀疏度  
sparsity = 0.01  
# 输入的维度  
input_dim = 3  
# 水库(或称为储备池)的大小  
reservoir_size = 1000  
# 预热步数  
n_steps_prerun = 10  
# 正则化参数  
regularization = 1e-2   

A.2.随机初始化网络权重

代码是用于初始化水库计算中的网络权重的Python代码。

import numpy as np
from scipy import sparse
from scipy.sparse.linalg import eigs# 定义水库大小和输入维度
reservoir_size = 1000  # 假设的水库大小
input_dim = 3         # 输入数据的维度
sparsity = 0.01      # 稀疏性参数
radius = 0.6          # 权重的半径# 随机初始化隐藏层权重矩阵,大小为reservoir_size x reservoir_size
weights_hidden = sparse.random(reservoir_size, reservoir_size, density=sparsity)# 计算隐藏层权重矩阵的特征值
eigenvalues, _ = eigs(weights_hidden)# 标准化权重以确保最大的特征值的绝对值不超过radius
weights_hidden = weights_hidden / np.max(np.abs(eigenvalues)) * radius# 初始化输入到隐藏层的权重矩阵,大小为reservoir_size x input_dim
weights_input = np.zeros((reservoir_size, input_dim))# 为每个输入维度分配一个子空间
q = int(reservoir_size / input_dim)
for i in range(input_dim):# 为每个输入维度随机生成-1到1之间的值weights_input[i * q:(i + 1) * q, i] = 2 * np.random.rand(q) - 1# 初始化输出层权重矩阵,大小为input_dim x reservoir_size
weights_output = np.zeros((input_dim, reservoir_size))

代码功能说明:

  • 首先,使用scipy.sparse.random函数随机生成一个稀疏的weights_hidden矩阵,该矩阵表示隐藏层神经元之间的连接权重。
  • 然后,使用scipy.sparse.linalg.eigs函数计算weights_hidden矩阵的特征值,以确保权重矩阵的条件数在合理范围内。
  • 接着,根据最大的绝对特征值调整权重矩阵的规模,以满足指定的radius参数。
  • 为输入层到隐藏层的连接初始化一个全零矩阵weights_input,然后为每个输入维度随机分配权重。
  • 最后,初始化输出层的权重矩阵weights_output为一个全零矩阵,这个矩阵将在后续的训练过程中被优化。

这些权重矩阵将用于水库计算模型中,该模型通过固定隐藏层权重并仅训练输出层来简化训练过程。

B.3 嵌入序列

将序列嵌入到网络的隐藏状态中。
代码定义了几个函数,用于初始化和处理水库计算中的隐藏状态。

import numpy as np# 初始化隐藏状态函数
def initialize_hidden(reservoir_size, n_steps_prerun, sequence):# 创建一个全零的隐藏状态矩阵,大小为(reservoir_size, 1)hidden = np.zeros((reservoir_size, 1))# 预热阶段,不更新隐藏状态,只用于填充初始隐藏状态for t in range(n_steps_prerun):# 获取输入序列并重塑形状input = sequence[t].reshape(-1, 1)# 更新隐藏状态,使用tanh激活函数hidden = np.tanh(weights_hidden @ hidden + weights_input @ input)return hidden# 增强隐藏状态的函数,例如通过平方操作
def augment_hidden(hidden):# 复制隐藏状态矩阵h_aug = hidden.copy()# 对每隔一个元素进行平方操作,增强状态表示h_aug[::2] = pow(h_aug[::2], 2.0)return h_aug# 使用序列数据和预热步数调用初始化隐藏状态函数
hidden = initialize_hidden(reservoir_size, n_steps_prerun, sequence)# 初始化存储隐藏状态和目标状态的列表
hidden_states = []
targets = []# 对序列中每个时间步进行循环
for t in range(n_steps_prerun, len(sequence) - 1):# 重塑输入和目标向量形状input = np.reshape(sequence[t], (-1, 1))target = np.reshape(sequence[t + 1], (-1, 1))# 更新隐藏状态hidden = np.tanh(weights_hidden @ hidden + weights_input @ input)# 增强隐藏状态hidden = augment_hidden(hidden)# 将当前隐藏状态添加到列表中hidden_states.append(hidden)# 将目标状态添加到列表中targets.append(target)# 将列表转换为数组,并去除单维度条目
targets = np.squeeze(np.array(targets))
hidden_states = np.squeeze(np.array(hidden_states))

代码功能说明:

  • initialize_hidden 函数用于初始化隐藏状态,通过预热步骤填充初始状态,但不更新状态。
  • augment_hidden 函数用于增强隐藏状态,例如通过平方操作来增加状态的表达能力。
  • 在主循环中,对于序列中的每个时间步,更新隐藏状态,并通过augment_hidden函数增强状态,然后将状态和目标添加到各自的列表中。
  • 最后,将隐藏状态和目标状态的列表转换为NumPy数组,并使用np.squeeze去除单维度条目,以便用于后续的线性回归或其他处理。

B.4.获取线性输出层权重

使用岭回归来获取线性输出层的权重。
代码展示了如何使用Python实现水库计算中的输出权重计算和预测函数,以及如何绘制预测结果与实际数据的对比图。

import numpy as np
import matplotlib.pyplot as plt# 计算输出权重矩阵
weights_output = (np.linalg.inv(hidden_states.T @ hidden_states + regularization * np.eye(reservoir_size)) @hidden_states.T @ targets).T# 定义预测函数
def predict(sequence, n_steps_predict):# 初始化隐藏状态hidden = initialize_hidden(reservoir_size, n_steps_prerun, sequence)# 重塑输入向量input = sequence[n_steps_prerun].reshape((-1, 1))# 初始化输出列表outputs = []for t in range(n_steps_prerun, n_steps_prerun + n_steps_predict):# 更新隐藏状态hidden = np.tanh(weights_hidden @ hidden + weights_input @ input)# 增强隐藏状态hidden = augment_hidden(hidden)# 计算输出output = weights_output @ hidden# 更新输入为当前输出input = output# 将输出添加到列表中outputs.append(output)# 返回预测结果作为数组return np.array(outputs)# 使用predict函数进行预测
x_sim = predict(sequence, 4000)# 绘制真实数据和预测结果的对比图
plt.figure(figsize=(6, 4))
plt.plot(x_train[:4000, 0], x_train[:4000, 2], label="Ground truth")
plt.plot(x_sim[:, 0], x_sim[:, 2], '--', label="Reservoir computing estimate")
plt.plot(x_train[0, 0], x_train[0, 2], "ko", label="Initial condition", markersize=8)
plt.legend()
plt.show()

代码功能说明:

  • weights_output 的计算使用了岭回归,其中 hidden_states 是隐藏状态的集合,targets 是目标值,regularization 是正则化参数,用于防止过拟合。
  • predict 函数接收一个序列和预测步数 n_steps_predict,初始化隐藏状态,并循环进行预测。在每一步,使用当前隐藏状态和输入来计算下一个时间步的输出。
  • initialize_hidden 函数用于初始化隐藏状态,augment_hidden 函数用于增强隐藏状态,例如通过平方操作。
  • 最后,使用 matplotlib 库绘制真实数据和水库计算估计结果的对比图,以及初始条件的标记。

附录B. 洛伦兹吸引子简介

洛伦兹吸引子(Lorenz Attractor)是一个在混沌理论中非常著名的动态系统模型,由气象学家爱德华·洛伦兹(Edward Lorenz)在1963年提出。这个模型描述了一个三变量的自治常微分方程组,通常用于展示混沌现象,特别是在大气对流的研究中。

洛伦兹吸引子的数学模型由以下三个方程组成:

d x d t = σ ( y − x ) \frac{dx}{dt} = \sigma(y - x) dtdx=σ(yx)
d y d t = x ( ρ − z ) − y \frac{dy}{dt} = x(\rho - z) - y dtdy=x(ρz)y
d z d t = x y − β z \frac{dz}{dt} = xy - \beta z dtdz=xyβz

其中,(x), (y), 和 (z) 是状态变量,而 (\sigma), (\rho), 和 (\beta) 是系统的参数。这组方程描述了流体对流过程中的三个关键变量:对流强度((x)),水平温度梯度((y)),以及垂直温度梯度((z))。

当参数 (\sigma), (\rho), 和 (\beta) 被设置为某些值时(例如,(\sigma = 10), (\rho = 28), (\beta = \frac{8}{3})),系统将会展现出混沌行为。此时,洛伦兹吸引子将表现为一个双螺旋状的吸引子,系统状态((x), (y), (z))的轨迹会在这个吸引子上无限循环,但永远不会重复。

洛伦兹吸引子的特性包括:

  1. 敏感性:系统对初始条件非常敏感,即使初始条件只有微小的差异,长期下来也会导致截然不同的结果。这就是所谓的“蝴蝶效应”。
  2. 混沌性:系统具有混沌特性,即其动态行为看起来是随机的,但实际上是由确定性方程产生的。
  3. 吸引性:尽管系统表现出混沌行为,但其轨迹仍然被限制在一个有限的区域内,即洛伦兹吸引子。

洛伦兹吸引子不仅在气象学和流体力学中有着广泛的应用,还成为了研究混沌现象和复杂系统动力学的重要工具。它揭示了自然界中许多复杂系统可能存在的普遍行为,对于理解这些系统的行为和控制它们具有重要的指导意义。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/354947.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

[创业之路-120] :全程图解:软件研发人员如何从企业的顶层看软件产品研发?

目录 一、企业全局 二、供应链 三、团队管理 四、研发流程IPD 五、软件开发流程 六、项目管理 七、研发管理者的自我修炼 一、企业全局 二、供应链 三、团队管理 四、研发流程IPD 五、软件开发流程 六、项目管理 七、研发管理者的自我修炼

系统架构设计师 - 数据库系统(1)

数据库系统 数据库系统数据库模式 ★分布式数据库 ★★★数据库设计阶段 ★★ER模型 ★关系模型 ★ ★结构约束条件完整性约束 关系代数 ★ ★ ★ ★概述自然连接 大家好呀!我是小笙,本章我主要分享系统架构设计师 - 数据库系统(1)知识,希望内…

idea插件开发之在项目右键添加菜单

写在前面 本文看下如何在右键列表中增加菜单。 正戏 首先创建一个Action,要显示的menu选择ProjectViewPopupMenu,如下: action public class CAction extends AnAction {Overridepublic void actionPerformed(AnActionEvent e) { // …

OSPF 动态路由协议(思科、华为)

#交换设备 OSPF 动态路由协议 一、基本概念 1.中文翻译:开放式最短路径优先路由协议(open shortest path first),是一个内部网关路由协议(一个自治系统内)2.也称为:链路状态路由协议&#xf…

火爆全网 LLM大模型教程:从零开始构建大语言模型,git突破18K标星

什么!一本书的Github仓库居然有18.5k的星标!(这含金量不必多说) 对GPT大模型感兴趣的有福了!这本书的名字叫 《Build a Large Language Model (From Scratch)》 也就是 从零开始构建大语言模型! 虽然这是一…

常说的云VR是什么意思?与传统vr的区别

虚拟现实(Virtual Reality,简称VR)是一种利用计算机技术模拟产生一个三维空间的虚拟世界,让用户通过视觉、听觉、触觉等感官,获得与现实世界类似或超越的体验。VR技术发展历程可追溯至上世纪,经历概念提出、…

鸿蒙 Web组件的生命周期(api10、11、12)

概述 开发者可以使用Web组件加载本地或者在线网页。 Web组件提供了丰富的组件生命周期回调接口,通过这些回调接口,开发者可以感知Web组件的生命周期状态变化,进行相关的业务处理。 Web组件的状态主要包括:Controller绑定到Web组…

两行css 实现瀑布流

html <ul ><li><a href"" ><img src"05094532gc6w.jpg" alt"111" /><p>传奇</p></a></li><li><a href"" ><img src"05094532gc6w.jpg" alt"111"…

文件防篡改监控工具 - WGCLOUD全面介绍

WGCLOUD是一款优秀的运维监控软件&#xff0c;免费、轻量、高效&#xff0c;部署容易&#xff0c;上手简单&#xff0c;对新手非常友好 WGCLOUD部署完成后&#xff0c;点击菜单【文件防篡改】&#xff0c;可以看到如下页面 我们点击【添加】按钮&#xff0c;输入监控文件的信息…

深圳比创达EMC|EMC与EMI测试整改:保障电子设备电磁兼容性步骤

随着电子技术的迅猛发展&#xff0c;电子设备在我们的日常生活中扮演着越来越重要的角色。然而&#xff0c;这些设备在运行时产生的电磁干扰&#xff08;EMI&#xff09;以及对外界电磁干扰的敏感性&#xff08;EMC&#xff09;问题&#xff0c;不仅影响着设备本身的性能&#…

Windows电脑部署Jellyfin服务端并进行远程访问配置详细教程

文章目录 前言1. Jellyfin服务网站搭建1.1 Jellyfin下载和安装1.2 Jellyfin网页测试 2.本地网页发布2.1 cpolar的安装和注册2.2 Cpolar云端设置2.3 Cpolar本地设置 3.公网访问测试4. 结语 前言 本文主要分享如何使用Windows电脑本地部署Jellyfin影音服务并结合cpolar内网穿透工…

代理模式(静态代理/动态代理)

代理模式&#xff08;Proxy Pattern&#xff09; 一 定义 为其他对象提供一种代理&#xff0c;以控制对这个对象的访问。 代理对象在客户端和目标对象之间起到了中介作用&#xff0c;起到保护或增强目标对象的作用。 属于结构型设计模式。 代理模式分为静态代理和动态代理。…

[Mysql] 的基础知识和sql 语句.教你速成(上)——逻辑清晰,涵盖完整

目录 前言 上篇的内容概况 下篇的内容概况 数据库的分类 关系型数据库 常见的关系型数据库系统 非关系型数据库 1. 键值对数据库&#xff08;Key-Value Stores&#xff09; 特点&#xff1a; 常见的键值对数据库&#xff1a; 2. 文档数据库&#xff08;Document Store…

好用的linux一键换源脚本

最近发现一个好用的linux一键换源脚本&#xff0c;记录一下 官方链接 大陆使用 bash <(curl -sSL https://linuxmirrors.cn/main.sh)# github地址 bash <(curl -sSL https://raw.githubusercontent.com/SuperManito/LinuxMirrors/main/ChangeMirrors.sh) # gitee地址 …

Redis 网络模型

一、用户空间和内核空间 1.1 linux 简介 服务器大多采用 Linux 系统&#xff0c;这里我们以 Linux 为例来讲解&#xff0c;下面有两个不同的 linux 发行版&#xff0c;分别位 ubuntu 和 centos&#xff0c;其实发行版就是在 Linux 系统上包了一层壳。 任何 Linux 发行版&#…

RS-232协议详解:深入理解与实际应用

RS-232协议详解 RS-232协议&#xff0c;也称为推荐标准232&#xff0c;是一种用于串行通信的标准协议。它在计算机和外围设备之间的通信中广泛应用。本文将详细介绍RS-232协议的各个方面&#xff0c;包括其历史、工作原理、信号类型、连接方式、应用场景等。希望通过这篇文章&a…

如何使用React的lazy和Suspense来实现代码分割?

在React中&#xff0c;使用React.lazy和Suspense可以方便地实现组件的代码分割。代码分割是一种优化技术&#xff0c;它将代码拆分成多个包&#xff0c;然后按需加载这些包&#xff0c;从而加快应用的初始加载时间。下面是如何使用这两个API的基本步骤&#xff1a; 使用React.l…

软考初级网络管理员__网络单选题

1.某人的电子邮箱为 Rjspks163.com,对于Rjspks和163.com的正确理解为(2)&#xff0c;在发送电子邮件时&#xff0c;常用关键词使用中&#xff0c;(3)是错误的&#xff0c;采用的协议是(4)。若电子邮件出现字符乱码现象&#xff0c;以下方法中(5)一定不能解决该问题。 SNMP SM…

【安防天下】模拟视频监控系统——模拟监控系统的构成视频采集设备

文章目录 1 模拟监控系统的构成2 视频采集设备2.1 摄像机相关技术2.1.1 摄像机的工作原理2.1.2 摄像机的分类2.1.3 摄像机的主要参数 2.2 镜头相关介绍2.2.1 镜头的主要分类2.2.2 镜头的主要参数 1 模拟监控系统的构成 模拟视频监控系统又称闭路电视监控系统&#xff0c; 一般…

DB9母头接口定义485

在通信技术中&#xff0c;DB9接口广泛应用于串行通信&#xff0c;尤其是在RS232和RS485标准中。虽然DB9接口最常见于RS232通信&#xff0c;但通过适当的引脚映射&#xff0c;它也可以用于RS485通信。本文将详细介绍如何定义和使用DB9母头接口进行RS485连接。 DB9母头接口简介 …