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

本文探讨了SINDy算法及其在非线性动力学系统识别中的应用。SINDy通过稀疏线性回归识别控制方程，强调模型可解释性和泛化能力。通过选择基函数库和优化器，以线性动态系统和洛伦兹吸引子为例验证了SINDy的有效性。尽管面临数据质量、变量选择和特征库构建等挑战，但SINDy仍展现出在复杂系统中应用的潜力。随着研究深入，SINDy有望在更多领域发挥重要作用。

MUKAMO

1759人浏览 · 2024-06-20 18:15:14

MUKAMO · 2024-06-20 18:15:14 发布

1. 引言

1.1. DNN模型的来由

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

识别一个系统意味着基于数据推断该系统的控制方程。换句话说，就是找到动态系统方程 $\mathbf{\dot x} = f(\mathbf{x})$ 中的 $f$ （其中 $\mathbf{x}$ 可能是向量值）。

例如，对于Lorenz系统，我们希望从时间序列数据中学习到方程右边的部分。

1.2. 研究稀疏性的意义

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

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

2. SINDy算法

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

从动态系统生成数据 $X$ 并计算导数 $\dot X$ 。
建立候选项库 $\Theta(X)$ 作为 $X$ 上的函数。
稀疏回归系数 $\Xi$ ，以最好地描述数据。

SINDy假设是测量了 $n$ 维数据点的时间序列 $\mathbf{x}=(x_1, \ldots x_n)$ 在 $m$ 个时间步 $t_1, \ldots, t_m$ ，我们定义数据矩阵 $X$ 和导数矩阵 $\dot X$ 为：
$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}$
$\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}$
接下来，我们定义库矩阵 $\Theta(X)$ ，其列是应用于数据的一组基函数 $\{\theta_l\}_{l=1,\ldots, L}$ ：
$\Theta(X)=\begin{bmatrix} |&|&&|\\ \Theta_1(X)&\Theta_2(X)&\cdots&\Theta_n(X)\\ |&|&&|\\ \end{bmatrix}$
以下是对上述文本的改写：
一些简单的例子包括多项式的基础，比如在泰勒级数展开中使用的多项式，或者三角函数，如 sin(x_1), cos(x_1), sin(2x_1), 等等，这些在傅里叶级数展开中常见。然而，根据不同的问题，可能需要使用更复杂的基础函数，例如贝塞尔函数。

最后，我们使用稀疏线性回归算法（例如LASSO）找到系数 $\Xi$ ，使得：
$\Xi=\begin{bmatrix} |&|&&|\\ \varepsilon_1(X)&\varepsilon_2(X)&\cdots&\varepsilon_n(X)\\ |&|&&|\\ \end{bmatrix}$
因此
$\dot{X}=\theta(X)\Xi$

3.线性动态系统示例

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

$\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}$

3.1.构建数据矩阵

以初始条件 $x_0=3$ 和 $y_0=\frac{1}{3}$ 我们构建数据矩阵 $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)

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

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

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+\alpha||w||^2_2$ 。

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

# 导入 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应用于真实世界数据更具挑战性，主要体现在以下几个方面：

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

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

附录 A：水库计算实例

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

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

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

然后，水库可以基于隐藏状态 $\mathbf{h}_{\tau}$ 线性预测下一个输入 $\mathbf{x}_{\tau +1}$ ：
$\mathbf{\hat {x}}_{\tau +1} = \mathbf{W}_{\mathrm{o}} \mathbf{h}_\tau$
其中 $\mathbf{W}_{\mathrm{o}}$ 将隐藏激活映射到输出。只有这个输出矩阵 $\mathbf{W}_{\mathrm{o}}$ 的参数被训练。因此，可以通过岭回归（正则化参数 $\lambda$ ）以解析方式获得最优值：
$W_{\mathrm{o}}=(H^TH+\lambda\mathrm{I})^{-1}H^TX$
其中 $\mathbf{H}$ 是隐藏状态， $\mathbf{I}$ 是单位矩阵， $\mathbf{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年提出。这个模型描述了一个三变量的自治常微分方程组，通常用于展示混沌现象，特别是在大气对流的研究中。

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

$\frac{dx}{dt} = \sigma(y - x)$
$\frac{dy}{dt} = x(\rho - z) - y$
$\frac{dz}{dt} = xy - \beta z$

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

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