赞
踩
从数据中识别非线性动态学意味着什么?
假设我们有时间序列数据,这些数据来自一个(非线性)动态学系统。
为什么我们需要稀疏性?
在这里,稀疏性意味着控制方程中的项数很少。稀疏性是有益的,因为它更具:
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)。算法分为三个步骤:
sin(x_1)
, cos(x_1)
, sin(2x_1)
, 等等,这些在傅里叶级数展开中常见。然而,根据不同的问题,可能需要使用更复杂的基础函数,例如贝塞尔函数。假设我们测量了一个由以下动力学系统控制的粒子轨迹:
d d t ( x y ) = ( − 2 x y ) = ( − 2 0 0 1 ) ( x y ) \frac{d}{dt} (xy) = (−2xy) = (−2001) (xy) dtd(xy)=(−2xy)=(−2001)(xy)
以初始条件 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 ps
t = np.linspace(0, 1, 100)
x = 3 * np.exp(-2 * t)
y = 0.5 * np.exp(t)
X = np.stack((x, y), axis=-1)
在构建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)
拟合模型后,可以使用 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 plt def 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)
正如我们预期的那样,对于这个问题,傅里叶基并不是最佳选择。现在,让我们尝试一个不同的基函数!
实验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 ∣∣y−Xw∣∣22+α∣∣w∣∣22。
阈值:这是权重向量中系数的最小允许幅度。任何幅度小于这个阈值的系数都将被清零,从而有助于实现模型的稀疏性。
# 导入 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)
当使用 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维度
将SINDy应用于真实世界数据更具挑战性,主要体现在以下几个方面:
在过去的几年中,这方面取得了很多进展,SINDy已成功应用于各种真实世界问题,包括发现由于问题复杂性而无法用解析方法表示的等离子体流体动力学方程。
水库计算是一种简单的方法,用于在没有反向传播通过时间的情况下训练递归神经网络,以及与之相关的臭名昭著的梯度消失和爆炸问题。基本步骤如下:
更正式地说,递归神经网络的状态由其隐藏激活
h
∈
R
M
\mathbf{h} \in \mathbb{R}^{M}
h∈RM给出,这些激活通过随机初始化并固定的矩阵
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是回归的目标。
# 导入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
代码是用于初始化水库计算中的网络权重的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
为一个全零矩阵,这个矩阵将在后续的训练过程中被优化。这些权重矩阵将用于水库计算模型中,该模型通过固定隐藏层权重并仅训练输出层来简化训练过程。
将序列嵌入到网络的隐藏状态中。
代码定义了几个函数,用于初始化和处理水库计算中的隐藏状态。
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
函数增强状态,然后将状态和目标添加到各自的列表中。np.squeeze
去除单维度条目,以便用于后续的线性回归或其他处理。使用岭回归来获取线性输出层的权重。
代码展示了如何使用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
库绘制真实数据和水库计算估计结果的对比图,以及初始条件的标记。洛伦兹吸引子(Lorenz Attractor)是一个在混沌理论中非常著名的动态系统模型,由气象学家爱德华·洛伦兹(Edward Lorenz)在1963年提出。这个模型描述了一个三变量的自治常微分方程组,通常用于展示混沌现象,特别是在大气对流的研究中。
洛伦兹吸引子的数学模型由以下三个方程组成:
d
x
d
t
=
σ
(
y
−
x
)
\frac{dx}{dt} = \sigma(y - x)
dtdx=σ(y−x)
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 x x, y y y, 和 z z z 是状态变量,而 σ \sigma σ, ρ \rho ρ, 和 β \beta β 是系统的参数。这组方程描述了流体对流过程中的三个关键变量:对流强度 x x x,水平温度梯度 y y y,以及垂直温度梯度 z z z。
当参数 σ \sigma σ, ρ \rho ρ 和 β \beta β 被设置为某些值时(例如, σ = 10 \sigma = 10 σ=10, ρ = 28 \rho = 28 ρ=28, β = 8 3 \beta = \frac{8}{3} β=38),系统将会展现出混沌行为。此时,洛伦兹吸引子将表现为一个双螺旋状的吸引子,系统状态( x x x, y y y, z z z)的轨迹会在这个吸引子上无限循环,但永远不会重复。
洛伦兹吸引子的特性包括:
洛伦兹吸引子不仅在气象学和流体力学中有着广泛的应用,还成为了研究混沌现象和复杂系统动力学的重要工具。它揭示了自然界中许多复杂系统可能存在的普遍行为,对于理解这些系统的行为和控制它们具有重要的指导意义。
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。