☰
分数阶时滞神经网络稳定性分析:Caputo导数、LMI判据与数值验证
2026/10/5 9:23:46 网站建设 项目流程

简介:这份PDF文献《含有离散时滞及分布时滞分数阶神经网络的渐近稳定性分析》面向从事神经网络、深度学习与数据建模方向的研究生、教师及工程技术人员,聚焦分数阶神经网络在同时含离散时滞与分布时滞情形下的渐近稳定性难题。文章在Caputo导数意义下构造Lyapunov函数,结合分数阶Razumikhin定理,给出渐近稳定性的充分条件,并通过四个实例验证定理有效性,为复杂动态系统的建模与预测精度优化提供理论工具。资源包共1个文件,为PDF格式,大小约1011KB,内容完整、便于检索与引用。目前已有104人学习下载,适合希望深入理解分数阶微积分、时滞系统稳定性及机器学习理论基础的读者参考研读。

1. 拆解一份分数阶时滞神经网络稳定性分析:从 Caputo 导数到四个数值验证

如果你正在做神经网络动力学分析,尤其是涉及记忆效应或时滞环节的建模,大概率会遇到一个尴尬:整数阶模型怎么调都跟实验数据对不上。这不是参数没调好,而是模型本身丢掉了系统的遗传和记忆特性。我手里这份《含有离散时滞及分布时滞分数阶神经网络的渐近稳定性分析》,就是专门解决这个问题的。它用 Caputo 导数重新刻画神经元状态演化,同时把离散时滞和分布时滞都塞进模型里,再用 Lyapunov 函数配合分数阶 Razumikhin 定理给出渐近稳定的充分条件。整篇论文最值钱的地方不是公式推导,而是它给了四个数值例子——两个满足条件、两个故意不满足——让你能直接对照着复现稳定与不稳定的边界。适合做控制、复杂网络、非线性动力学方向的从业者,也适合需要给时滞系统做稳定性判据的工程人员。

2. 模型怎么建:从整数阶到分数阶时滞系统的改造逻辑

2.1 为什么非要用 Caputo 导数而不是 Riemann-Liouville

分数阶微积分有好几种定义,Riemann-Liouville 和 Caputo 是最常用的两个。这份论文从头到尾用的是 Caputo 导数,原因很实际:Caputo 导数的初始条件跟整数阶一样,可以直接用函数在零时刻的值,不需要分数阶积分初值。Riemann-Liouville 虽然数学上更“原生”,但它的初始条件物理意义不直观,做数值仿真时初值处理特别容易翻车。

论文里给出的 Caputo 导数定义是:

$$^{C}_{t_0}D^{\alpha}t f(t) = \frac{1}{\Gamma(n-\alpha)}\int{t_0}^{t}\frac{f^{(n)}(\tau)}{(t-\tau)^{\alpha+1-n}}d\tau$$

其中 $0 \le n-1 < \alpha < n$,$n$ 是正整数。当 $0 < \alpha < 1$ 时,$n=1$,公式简化为:

$$^{C}_{t_0}D^{\alpha}t f(t) = \frac{1}{\Gamma(1-\alpha)}\int{t_0}^{t}\frac{f'(\tau)}{(t-\tau)^{\alpha}}d\tau$$

这个形式的好处是:被积函数里出现的是 $f'(\tau)$,也就是普通的一阶导数。你在做数值离散的时候,只需要处理一阶导数的历史积分,不用去碰分数阶导数的初值问题。常见做法是用 Grünwald-Letnikov 格式做离散化,或者用 Adams-Bashforth-Moulton 预测-校正算法,论文里数值例子参考了文献[22]的算法,大概率就是这类方法。

注意:Caputo 导数的阶数 $\alpha$ 必须严格在 $(0,1)$ 区间内,论文里所有定理和例子都默认这个范围。如果你把 $\alpha$ 设成大于 1,整个 Lyapunov 分析框架要重写。

2.2 离散时滞和分布时滞分别怎么进模型

论文考虑了两类系统。第一类是只含离散时滞的:

$$^{C}_{t_0}D^{\alpha}_t x(t) = -Ax(t) + Bf(x(t)) + Cg(x(t-r_1(t)))$$

第二类是同时含离散时滞和分布时滞的:

$$^{C}_{t_0}D^{\alpha}t x(t) = -Ax(t) + Bf(x(t)) + Cg(x(t-r_1(t))) + D\int{t-r_2(t)}^{t}h(x(s))ds$$

这里 $x(t) \in \mathbb{R}^n$ 是状态向量,$A$ 是正定对角矩阵,$A = \text{diag}{a_i}$,$a_i > 0$ 表示第 $i$ 个神经元在没有网络连接和外部输入时电位重置到静止状态的速率。$B$、$C$、$D$ 是权重矩阵,$f$、$g$、$h$ 是激励函数,都满足 $f(0)=g(0)=h(0)=0$。

离散时滞 $r_1(t)$ 和分布时滞 $r_2(t)$ 都是时变的,满足 $0 \le r_1(t) \le r_1$,$0 \le r_2(t) \le r_2$。分布时滞那一项是积分形式,表示过去一段时间内所有历史状态对当前状态的累积影响。这个积分项是整篇论文最麻烦的地方,也是为什么定理 2 比定理 1 多出两个矩阵不等式条件。

激励函数满足 Lipschitz 条件:

$$|f(u)-f(v)| \le l_1|u-v|,|g(u)-g(v)| \le l_2|u-v|,|h(u)-h(v)| \le l_3|u-v|$$

$l_1$、$l_2$、$l_3$ 是非负常数。这个条件在数值例子里的取值很关键,后面会看到。

2.3 平衡点平移和误差系统转化

原系统有外部输入 $I$,平衡点 $y^$ 一般不在原点。论文做了一个标准操作:令 $x(t) = y(t) - y^$,把平衡点平移到原点。变换后的系统变成:

$$^{C}_{t_0}D^{\alpha}_t x(t) = -Ax(t) + Bf(x(t)) + Cg(x(t-r_1(t)))$$

以及带分布时滞的版本。注意变换后的 $f$、$g$、$h$ 定义变成了 $f(x(t)) = f(y(t)) - f(y^*)$,以此类推。这个平移操作在数值仿真时特别重要,如果你直接拿原系统跑仿真,平衡点不在原点,Lyapunov 函数的构造和稳定性判断都会偏。

3. 稳定性判据怎么用:两个定理和 Schur 补的实操细节

3.1 定理 1 的矩阵不等式怎么验证

定理 1 针对只含离散时滞的系统。如果 Lipschitz 条件成立,且存在常数 $k>0$ 和对称正定矩阵 $P$ 满足:

$$\begin{pmatrix} -PA - A^TP + PBB^TP + l_1^2I + PCC^TP + kP & 0 \ 0 & l_2^2I - kP \end{pmatrix} < 0$$

那么系统 (15) 的零解是渐近稳定的。

这个矩阵不等式是一个 $2n \times 2n$ 的对称矩阵,要求负定。直接验证负定性在 $n$ 稍微大一点的时候就很麻烦。论文用了 Schur 补引理来降维。Schur 补说的是:对于分块对称矩阵

$$S = \begin{pmatrix} S_{11} & S_{12} \ S_{12}^T & S_{22} \end{pmatrix}$$

$S < 0$ 等价于 $S_{11} < 0$ 且 $S_{22} - S_{12}^T S_{11}^{-1} S_{12} < 0$,或者 $S_{22} < 0$ 且 $S_{11} - S_{12} S_{22}^{-1} S_{12}^T < 0$。

在定理 1 里,$S_{12} = 0$,所以条件退化成两个独立的块对角条件:左上块负定,右下块负定。右下块是 $l_2^2I - kP$,取 $P=I$、$k=1$、$l_2=1/2$ 时,$l_2^2I - kP = (1/4 - 1)I = -3/4I < 0$,自动满足。左上块需要具体算。

论文例 1 给了具体数值:

$$A = \begin{pmatrix} 2 & 0 \ 0 & 2 \end{pmatrix},B = \begin{pmatrix} 0.2 & 0 \ 0 & 0.3 \end{pmatrix},C = \begin{pmatrix} 0.5 & -0.1 \ -0.1 & 0.3 \end{pmatrix}$$

$f(x) = (\sin x_1, \sin x_2)^T$,$g(x) = \frac{1}{2}(\tanh x_1, \tanh x_2)^T$,所以 $l_1=1$,$l_2=1/2$。取 $k=1$,$P=I_2$,计算左上块:

$$-PA - A^TP + PBB^TP + l_1^2I + PCC^TP + kP$$

逐项算:$-PA - A^TP = -2A = \text{diag}(-4,-4)$。$PBB^TP = BB^T = \text{diag}(0.04, 0.09)$。$l_1^2I = I$。$PCC^TP = CC^T$,算出来是 $\begin{pmatrix} 0.26 & -0.08 \ -0.08 & 0.10 \end{pmatrix}$。$kP = I$。加起来:

$$\begin{pmatrix} -4+0.04+1+0.26+1 & -0.08 \ -0.08 & -4+0.09+1+0.10+1 \end{pmatrix} = \begin{pmatrix} -1.70 & -0.08 \ -0.08 & -1.81 \end{pmatrix}$$

这个矩阵负定,因为对角线都是负的,且行列式 $(-1.70)(-1.81) - (-0.08)^2 = 3.077 - 0.0064 > 0$。所以定理 1 条件满足,系统渐近稳定。

3.2 定理 2 多出来的两个条件怎么处理

定理 2 针对同时含离散时滞和分布时滞的系统。条件变成两个矩阵不等式:

$$\begin{pmatrix} -PA - A^TP + PBB^TP + l_1^2I + PCC^TP + PDD^TP + (k_3r_2 + k_1)P & 0 \ 0 & l_2^2I - k_1P \end{pmatrix} < 0$$

$$\begin{pmatrix} (k_2 - k_3)P & 0 \ 0 & r_2l_3^2I - k_2P \end{pmatrix} < 0$$

其中 $k_1 > 0$,$k_3 > k_2 > 0$。第二个条件比较直观:$(k_2 - k_3)P < 0$ 因为 $k_3 > k_2$ 且 $P$ 正定,自动满足;$r_2l_3^2I - k_2P < 0$ 需要选 $k_2$ 足够大,或者 $r_2$ 和 $l_3$ 足够小。

论文例 3 给了具体数值。$A$、$B$、$C$ 跟例 1 一样,新增 $D = \begin{pmatrix} 0.2 & -0.3 \ 0 & -0.3 \end{pmatrix}$。$h(x) = (\tanh x_1, \tanh x_2)^T$,所以 $l_3=1$。$r_2(t) = 0.5\sin^2 t \le 0.5 = r_2$。取 $k_1=k_2=1$,$k_3=2$,$P=I_2$。

先验证第二个条件:$(k_2-k_3)P = -I_2 < 0$,$r_2l_3^2I - k_2P = 0.5I - I = -0.5I < 0$,都满足。

再算第一个条件的左上块:

$$-PA - A^TP + PBB^TP + l_1^2I + PCC^TP + PDD^TP + (k_3r_2 + k_1)P$$

前面几项跟例 1 一样,新增 $PDD^TP = DD^T$,算出来是 $\begin{pmatrix} 0.13 & -0.09 \ -0.09 & 0.09 \end{pmatrix}$。$(k_3r_2 + k_1)P = (2 \times 0.5 + 1)I = 2I$。加起来:

$$\begin{pmatrix} -4+0.04+1+0.26+0.13+2 & -0.08-0.09 \ -0.08-0.09 & -4+0.09+1+0.10+0.09+2 \end{pmatrix} = \begin{pmatrix} -0.57 & -0.17 \ -0.17 & -0.72 \end{pmatrix}$$

这个矩阵负定。右下块 $l_2^2I - k_1P = 0.25I - I = -0.75I < 0$。所以定理 2 条件满足,系统渐近稳定。

3.3 数值仿真时怎么选离散化格式

论文没有给出具体的离散化代码,但提到了“算法参考了文献[22]”。在分数阶时滞神经网络的数值仿真里,常见做法是 Adams-Bashforth-Moulton 预测-校正算法。下面给一个 Python 实现的骨架,以例 1 的系统为例:

import numpy as np from scipy.special import gamma def caputo_abm_simulate(alpha, A, B, C, f, g, r1_func, x0_func, T, h): """ Adams-Bashforth-Moulton 预测-校正算法求解分数阶时滞神经网络 alpha: 分数阶阶数, 0 < alpha < 1 A, B, C: 系统矩阵 f, g: 激励函数 r1_func: 时滞函数 r1(t) x0_func: 历史函数 x(t), t in [-r1_max, 0] T: 仿真总时长 h: 步长 """ N = int(T / h) n = A.shape[0] x = np.zeros((N+1, n)) # 初始化历史项 for i in range(N+1): t = i * h if t <= 0: x[i] = x0_func(t) # 预计算系数 coeff = h**alpha / gamma(alpha + 2) for k in range(1, N+1): t_k = k * h # 预测步 sum_f = np.zeros(n) for j in range(k): t_j = j * h r1_j = r1_func(t_j) idx_delay = max(0, int((t_j - r1_j) / h)) x_delay = x[idx_delay] rhs = -A @ x[j] + B @ f(x[j]) + C @ g(x_delay) if j == 0: sum_f += (k - j)**alpha - (k - j - alpha) * k**alpha else: sum_f += ((k - j + 1)**alpha - 2*(k - j)**alpha + (k - j - 1)**alpha) * rhs x_pred = x[0] + coeff * sum_f # 校正步 r1_k = r1_func(t_k) idx_delay_k = max(0, int((t_k - r1_k) / h)) rhs_k = -A @ x_pred + B @ f(x_pred) + C @ g(x[idx_delay_k]) sum_corr = np.zeros(n) for j in range(k): t_j = j * h r1_j = r1_func(t_j) idx_delay = max(0, int((t_j - r1_j) / h)) rhs = -A @ x[j] + B @ f(x[j]) + C @ g(x[idx_delay]) if j == 0: sum_corr += (k - j + 1)**alpha - (k - j - alpha) * k**alpha else: sum_corr += ((k - j + 1)**alpha - 2*(k - j)**alpha + (k - j - 1)**alpha) * rhs x[k] = x[0] + coeff * (sum_corr + rhs_k) return np.arange(N+1) * h, x

这段代码的逻辑说明:预测步用历史项加权求和得到 $x_{pred}$,校正步把 $x_{pred}$ 代入右端函数再算一次加权求和。权重系数来自分数阶 Adams 格式的卷积核。时滞项通过查历史数组实现,idx_delay是延迟时刻对应的数组下标。参数alpha控制分数阶阶数,h是步长,一般取 0.01 或更小。r1_func是时滞函数,例 1 里是常数 1,例 3 里 $r_2(t) = 0.5\sin^2 t$。

注意:分布时滞的积分项在离散化时需要额外处理,一般用梯形法或矩形法近似积分。积分区间 $[t-r_2(t), t]$ 内的历史状态都要参与计算,计算量比纯离散时滞大很多。

4. 避坑与排查:四个数值例子背后的血泪经验

4.1 现象:定理条件满足但仿真发散

原因:最常见的是平衡点没有平移到原点。论文里的定理都是针对变换后的误差系统 $x(t) = y(t) - y^$ 给出的。如果你直接拿原系统 $y(t)$ 跑仿真,平衡点 $y^$ 不在原点,Lyapunov 函数 $V = x^TPx$ 的构造就不对,稳定性判据自然失效。

解决:先解方程 $-Ay^* + Bf(y^) + Cg(y^) + I = 0$ 求出平衡点,然后令 $x = y - y^*$ 再做仿真。如果激励函数是 $\sin$ 和 $\tanh$ 这类奇函数且 $I=0$,平衡点就是原点,可以跳过这一步。

4.2 现象:Schur 补条件算出来一个块负定一个块不正定

原因:$l_2^2I - kP$ 这一块对 $k$ 和 $P$ 的取值很敏感。$k$ 太小或者 $P$ 的特征值太大,这一块就可能不正定。论文例 1 里取 $k=1$、$P=I$、$l_2=1/2$,$l_2^2I - kP = -0.75I$ 是负定的。但如果你把 $l_2$ 改成 2,$l_2^2I - kP = 4I - I = 3I$ 正定,条件就不满足了。

解决:先检查 Lipschitz 常数 $l_2$ 是否估计准确。$\tanh$ 的导数最大是 1,所以 $g(x) = \frac{1}{2}\tanh x$ 的 Lipschitz 常数确实是 $1/2$。如果激励函数选得不对,$l_2$ 可能远大于 1。另外可以调大 $k$,让 $kP$ 主导 $l_2^2I$。

4.3 现象:分布时滞的积分项导致仿真速度极慢

原因:分布时滞项 $D\int_{t-r_2(t)}^{t}h(x(s))ds$ 在每个时间步都要对历史区间做积分。如果步长 $h$ 取 0.001,仿真时长 100 秒,每个时间步要算 1000 个历史点的加权和,总计算量是 $10^8$ 量级。

解决:分布时滞的积分可以用递推公式加速。如果 $r_2(t)$ 是常数,积分区间长度固定,可以用滑动窗口的方式更新积分值,每次只减去移出窗口的项、加上新进入窗口的项。如果 $r_2(t)$ 时变,可以用插值法近似。论文例 3 里 $r_2(t) = 0.5\sin^2 t$,最大 0.5,积分区间不长,直接算也能接受。

4.4 现象:例 2 和例 4 的仿真结果跟理论判断对不上

原因:例 2 和例 4 是故意构造的不满足定理条件的例子。例 2 里 $B$ 的第一行第一列从 0.2 改成 2,左上块变成 $\begin{pmatrix} 2.26 & -0.08 \ -0.08 & -1.81 \end{pmatrix}$,这个矩阵不是负定的,因为第一行第一列是正的。例 4 里 $B$ 同样改成 2,左上块变成 $\begin{pmatrix} 3.39 & -0.17 \ -0.17 & -0.72 \end{pmatrix}$,也不是负定的。

解决:这两个例子的仿真结果确实显示系统不收敛,跟理论判断一致。如果你复现时发现例 2 或例 4 收敛了,大概率是数值格式的步长太大或者仿真时间太短。把步长减小到 0.001,仿真时间延长到 50 秒以上,不稳定的系统会慢慢发散出来。

4.5 现象:分数阶阶数 $\alpha$ 接近 1 时结果跟整数阶对不上

原因:$\alpha \to 1$ 时 Caputo 导数退化成普通一阶导数,分数阶系统应该趋近整数阶系统。但数值格式在 $\alpha$ 接近 1 时会有精度问题,卷积核的奇异性虽然减弱了,但权重系数的计算容易出现数值误差。

解决:如果要做 $\alpha \to 1$ 的极限验证,建议用 $\alpha = 0.9, 0.95, 0.99$ 分别跑,看趋势是否收敛到整数阶结果。不要直接取 $\alpha = 1$,因为 Adams 格式在 $\alpha = 1$ 时公式会退化,需要单独处理。

5. 进阶用法:把定理条件写成可复用的 LMI 检查脚本

论文里的定理条件本质上是一组线性矩阵不等式(LMI)。虽然论文没有用 LMI 工具箱,但你可以把定理 1 和定理 2 的条件写成 Python 脚本,用cvxpy或者scipy做可行性检查。下面给一个基于scipy.optimize的简化版本,用于验证定理 1 的条件:

import numpy as np from scipy.optimize import minimize def check_theorem1(A, B, C, l1, l2, n): """ 检查定理1条件是否可行 寻找对称正定矩阵 P 和常数 k > 0 使得 2n x 2n 矩阵负定 """ # 参数化 P = L L^T, L 为下三角矩阵 def unpack(params): L_vals = params[:n*(n+1)//2] k = np.exp(params[-1]) # 保证 k > 0 L = np.zeros((n, n)) idx = 0 for i in range(n): for j in range(i+1): L[i, j] = L_vals[idx] idx += 1 P = L @ L.T + 1e-6 * np.eye(n) # 保证正定 return P, k def objective(params): P, k = unpack(params) # 构造左上块 S11 = -P @ A - A.T @ P + P @ B @ B.T @ P + l1**2 * np.eye(n) + P @ C @ C.T @ P + k * P # 构造右下块 S22 = l2**2 * np.eye(n) - k * P # 要求 S11 < 0 且 S22 < 0 # 用最大特征值作为惩罚 eig11 = np.max(np.linalg.eigvalsh(S11)) eig22 = np.max(np.linalg.eigvalsh(S22)) return max(eig11, eig22) # 初始猜测 x0 = np.random.randn(n*(n+1)//2 + 1) * 0.1 res = minimize(objective, x0, method='Nelder-Mead', options={'maxiter': 10000, 'xatol': 1e-8}) P, k = unpack(res.x) S11 = -P @ A - A.T @ P + P @ B @ B.T @ P + l1**2 * np.eye(n) + P @ C @ C.T @ P + k * P S22 = l2**2 * np.eye(n) - k * P return res.fun < 0, P, k, np.max(np.linalg.eigvalsh(S11)), np.max(np.linalg.eigvalsh(S22)) # 用例1的数据测试 A = np.array([[2, 0], [0, 2]]) B = np.array([[0.2, 0], [0, 0.3]]) C = np.array([[0.5, -0.1], [-0.1, 0.3]]) feasible, P, k, eig11, eig22 = check_theorem1(A, B, C, l1=1.0, l2=0.5, n=2) print(f"可行: {feasible}, k={k:.4f}, max_eig_S11={eig11:.6f}, max_eig_S22={eig22:.6f}")

这段脚本的逻辑说明:把 $P$ 参数化为 $LL^T$ 保证正定,$k$ 用指数变换保证正数。目标函数是左上块和右下块最大特征值的较大者,如果这个值小于 0,说明两个块都负定,定理条件可行。Nelder-Mead是无梯度优化方法,适合这种目标函数不平滑的情况。参数l1、l2是 Lipschitz 常数,n是神经元个数。

用这个脚本跑例 1 的数据,应该输出可行: True,并且max_eig_S11和max_eig_S22都是负数。跑例 2 的数据($B$ 改成 $\begin{pmatrix} 2 & 0 \ 0 & 0.3 \end{pmatrix}$),应该输出可行: False,因为左上块的最大特征值是正的。

这个脚本的局限是只检查了定理 1,定理 2 需要额外加两个矩阵不等式条件。另外Nelder-Mead在高维情况下容易陷入局部最优,神经元个数 $n$ 超过 5 的时候建议换cvxpy或者scipy.optimize.minimize的SLSQP方法。

提示:如果你用cvxpy,可以直接把矩阵不等式写成cp.Variable和cp.constraints,求解器会自动处理正定约束和 LMI 约束,比手写优化目标稳定得多。

从那以后我每次拿到这类稳定性判据的论文,都强制走一遍“先平移平衡点、再验证 LMI、最后跑数值仿真”的流程,少一步都可能翻车。希望帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询