陈茂森论文

随机移动网络系统的稳定性

马尔科夫链与网络平均度推导

1.马尔科夫链的基本概念

马尔科夫链描述的是这样一种随机过程：系统在若干个可能的状态中变化，下一时刻所处状态只依赖于当前状态，而与过去的状态无关，这就是所谓的“无记忆性”或马尔科夫性。

无记忆性意味着，对于任何 $s, t \ge 0$，

假设你已经等待了 $s$ 分钟，那么再等待至少 $t$ 分钟的概率，和你一开始就等待至少 $t$ 分钟的概率完全相同。

在所有概率分布里，只有指数分布

链路状态的马尔科夫模型

考虑网络中每条链路的动态行为，其状态空间为：

• 状态0：链路断开
• 状态1：链路连通

定义概率函数：

• $p_1(t)$：时刻 $t$ 处于连通状态的概率
• $p_0(t) = 1 - p_1(t)$：断开概率

同时，我们假设链路从一个状态转移到另一个状态需要等待一段时间，这段等待时间通常服从指数分布（论文中通过 KS 检验确认）：

• 从断开（0）到连通（1）的等待时间 $T_{01} \sim \text{Exp}(\lambda_{01})$
• 从连通（1）到断开（0）的等待时间 $T_{10} \sim \text{Exp}(\lambda_{10})$

其中，$\lambda_{01}$ 和 $\lambda_{10}$ 为转移速率，表示单位时间内事件（转移）发生的平均次数

2.推导单条链路的连通概率

根据连续时间马尔科夫链的理论，我们可以写出状态转移的微分方程。对于状态1（连通状态），概率 $p_1(t)$ 的变化率由两个部分组成：

1. 当链路处于状态0时，以速率 $\lambda_{01}$ 变为状态1。这部分概率增加的速率为
   $$ \lambda_{01} , p_0(t)=\lambda_{01} (1-p_1(t)). $$

2. 当链路处于状态1时，以速率 $\lambda_{10}$ 转换为状态0。这部分使 $p_1(t)$ 减少，其速率为
   $$ \lambda_{10} , p_1(t). $$

所以，$p_1(t)$ 的微分方程写成：

这个方程可以整理为：

这其实是一个一阶线性微分方程，其标准求解方法是求解其齐次解与非齐次解。

3. 求解微分方程

整个微分方程的通解为：

利用初始条件 $p_1(0)=p_1^0$（初始时刻链路连通的概率），我们可以求出 $C$：

即

所以，链路在任意时刻 $t$ 连通的概率为：

这就是单条链路的连通概率函数，描述了从任意初始条件出发，经过一段时间后，链路达到平衡状态的过程。

4.推导网络平均度的变化函数

在一个由 $N$ 个节点构成的网络中，每个节点都与其它节点进行通信（不考虑自环），因此每个节点最多有 $N-1$ 个邻居。对于任意一对节点 $i$ 和 $j$，它们之间链路连通的概率 $p_1(t)$（假设所有链路独立且同分布）。

• 某个节点 $i$ 在时刻 $t$ 的度 $d_i(t)$可以写作： $$ d_i(t)= \sum_{\substack{j=1 \ j\neq i}}^N p_{ij}(t), $$ 其中 $p_{ij}(t)=p_1(t)$。

• 因此，每个节点的期望度为： $$ E[d_i(t)]=(N-1)p_1(t). $$

• 网络平均度就是对所有节点的期望度取平均，由于网络中每个节点都遵循相同统计规律，所以网络平均度可表示为： $$ \bar{d}(t) = \frac{1}{N}\sum_{i=1}^{N} E[d_i(t)] = \frac{N \cdot (N-1)p_1(t)}{N} = (N-1)p_1(t) $$

将我们前面得到的 $p_1(t)$ 表达式代入，就得到网络平均度随时间变化的表达式：

这就是网络平均度的变化函数：

• 网络开始时每条链路的连通概率为 $p_1^0$
• 这种调整过程符合指数衰减规律，即偏离平衡值的部分按 $e^{-(\lambda_{01}+\lambda_{10})t}$ 衰减
• 当 $t$ 趋向无穷大时，指数项 $e^{-(\lambda_{01}+\lambda_{10})t}$ 衰减为0，网络平均度趋向于 $(N-1)\frac{\lambda_{01}}{\lambda_{01}+\lambda_{10}}$，这就是网络达到平衡状态后的理论平均度。

特征信号参数的平稳性

证明系统在平衡态下具有统计上的稳定性。

从节点空间分布证明平稳性。

设节点在模型区域的坐标为 $(X,Y)$，其分布概率密度函数写为

那么节点在模型子区域 $R_1$ 中出现的概率为

在平衡状态下，理论上节点的位置分布 $f(x,y)$ 保持不变，即每个区域内节点出现的概率 $P_{R_1}$ 是常数，不随时间变化。证明在平衡状态下节点分布稳定。

扰动后的恢复能力

• 静止节点分布特性满足均匀分布,概率密度函数为 $g(x,y)$；
• 运动节点的概率密度函数为 $h(x,y)$；
• 在时刻 $t_0$ 时，网络中共有 $N$ 个节点，其中有 $s$ 个静止，故静止节点的比例为 $p=\frac{s}{N}$。

节点整体的分布概率密度函数可写为

在平衡状态下，$p$ 的理论值为一个常数，所以 $f(x,y)$ 不随时间变化，从而网络连通度稳定。

接下来，考虑外界扰动的影响：假设在时刻 $t_1$ 新加入 $m$ 个符合均匀分布的节点，

扰动后的总分布（$t_1$时刻后）

• 新加入的 $m$ 个节点是静止的，其分布为 $g(x,y)$
• 此时网络的总节点数 $N+m$：
  • 静止节点总数：$s+m$
  • 运动节点总数：$N-s$（原有运动节点数不变）

因此，扰动后的分布为：

其中 $p' = \frac{s + m}{N + m}$。

近似处理（当 $N, s \gg m$ 时）

系统稳定性分析

平衡点及误差坐标

论文第 2.1 节推导出，单条链路连通概率 $p_1(t)$ 满足

记常数

构造李雅普诺夫函数

对一维系统 $\dot e=-ce$（$c>0$），自然选取

• $V(e)>0$ 当且仅当 $e\neq0$；
• 平衡点 $e=0$ 时，$V(0)=0$。

计算 $V$ 的时间导数

对 $V$ 关于时间求导：

因为 $c>0$ 且 $e^2\ge0$，所以

• 当 $e\neq0$ 时，$\dot V<0$；
• 当 $e=0$ 时，$\dot V=0$。

这正是“半负定”（negative semi-definite）的定义。

结论

李雅普诺夫第二类定理告诉我们：

若存在一个函数 $V(e)$ 在平衡点处为 0、在邻域内正定，且其导数 $\dot V(e)$ 在该邻域内为半负定，则平衡点 $e=0$（即 $d=d^*$）是稳定的。

由于我们已经构造了满足上述条件的 $V(e)=e^2$，并验证了 $\dot V(e)\le0$，故平衡态 $d=d^*$ 是 李雅普诺夫意义下稳定 的。

网络特征谱参数的估算

由于邻接矩阵不能保证半正定性，因此会产生幂迭代估算过程不能收敛的问题。需构造$A^T A$

基于奇异值分解改进幂迭代估算（集中式）

输入：矩阵 $B = A^T A$，目标特征值数量 $k$，收敛阈值 $\delta$
输出：前 $k$ 个特征值 $\lambda_1' \geq \lambda_2' \geq \dots \geq \lambda_k'$ 及对应特征向量 $u_1', u_2', \dots, u_k'$

1. 初始化

1. 随机生成初始非零向量 $v^{(0)}$，归一化： $$ v^{(0)} \gets \frac{v^{(0)}}{|v^{(0)}|_2} $$

2. 设置已求得的特征值数量 $n \gets 0$，剩余矩阵 $B_{\text{res}} \gets B$

2. 迭代求前k个特征值与特征向量

While $n < k$:

1. 幂迭代求当前最大特征值与特征向量

   • 初始化向量 $v^{(0)}$（若 $n=0$，用随机向量；否则用与已求特征向量正交的向量）

   • Repeat: a. 计算 $v^{(t+1)} \gets B_{\text{res}} v^{(t)}$​ b. 归一化： $$ v^{(t+1)}\gets \frac{v^{(t+1)}}{|v^{(t+1)}|2} $$ c. 计算 Rayleigh 商： $$ y^{(t)} = \frac{(v^{(t)})^T B{\text{res}} v^{(t)}}{(v^{(t)})^T v^{(t)}} $$

     d. Until $|y^{(t)} - y^{(t-1)}| < \delta$（收敛）

   • 记录当前特征值与特征向量： $$ \lambda_{n+1}' \gets y^{(t)}, \quad u_{n+1}' \gets v^{(t)} $$

2. 收缩矩阵以移除已求特征分量 每次收缩操作将已求得的特征值从矩阵中“移除”，使得剩余矩阵的谱（特征值集合）中次大特征值“升级”为最大特征值。

   • 更新剩余矩阵： $$ B_{\text{res}} \gets B_{\text{res}} - \lambda_{n+1}' u_{n+1}' (u_{n+1}')^T $$

   • 确保 $B_{\text{res}}$ 的对称性（数值修正）

3. 增量计数

   • $n \gets n + 1$

瑞利商公式

1. 集中式：

   $$ y(k)= \frac{x(k)^T A x(k)}{x(k)^T x(k)} $$

2. 分布式一致性计算：

   $$ y(k) = \frac{\sum_{i=1}^N x_i(k) b_i(k)}{\sum_{i=1}^N x_i^2(k)} $$ 其中 $$ b_i(k) = \sum_j a_{ij} x_j(k) $$

两者是等价的：

考虑一个简单的2×2矩阵

1. 集中式计算

2. 分布式计算
   各节点分别计算本地观测值

​ 节点1的计算：

​ 节点2的计算：

主要符号表

分布式幂迭代求前$K$大奇异对

While $n < K$:

1. 初始化：

   • 若 $n = 0$：

     • 各节点$j$随机初始化 $v_{0,j}^{(0)} \sim \mathcal{N}(0,1)$

   • 若 $n > 0$：

     • 分布式Gram-Schmidt正交化： $$ v_{n,j}^{(0)} \gets v_{n,j}^{(0)} - \sum_{m=0}^{n-1} \underbrace{\text{Consensus}\left(\sum_k v_{m,k} v_{n,k}^{(0)}\right)}{\text{全局内积}\langle v_m, v_n^{(0)} \rangle} v{m,j} $$

     • 分布式归一化： $$ v_{n,j}^{(0)} \gets \frac{v_{n,j}^{(0)}}{\sqrt{\text{Consensus}\left(\sum_k (v_{n,k}^{(0)})^2\right)}} $$

2. 迭代计算：

   • Repeat： a. 第一轮通信（计算$z=Av$）： $$ z_j^{(t)} = \sum_{k \in 𝒩_j} a_{jk} v_{n,k}^{(t)} \quad \text{(邻居交换$v_{n,k}^{(t)}$)} $$ b. 第二轮通信（计算$y=A^T z$）： $$ y_j^{(t+1)} = \sum_{k \in 𝒩_j} a_{kj} z_k^{(t)} \quad \text{(邻居交换$z_k^{(t)}$)} $$ c. 隐式收缩（$n>0$时）： $$ y_j^{(t+1)} \gets y_j^{(t+1)} - \sum_{m=0}^{n-1} \sigma_m^2 v_{m,j} \cdot \underbrace{\text{Consensus}\left(\sum_k v_{m,k} y_k^{(t+1)}\right)}{\text{投影系数计算}} $$ d. 归一化： $$ v{n,j}^{(t+1)} = \frac{y_j^{(t+1)}}{\sqrt{\text{Consensus}\left(\sum_k (y_k^{(t+1)})^2\right)}} $$ e. 计算Rayleigh商： $$ \lambda^{(t)} = \text{Consensus}\left(\sum_k v_{n,k}^{(t)} y_k^{(t+1)}\right) $$ f. 终止条件： $$ \text{If } \frac{|\lambda^{(t)} - \lambda^{(t-1)}|}{|\lambda^{(t)}|} < \delta \text{ then break} $$

3. 保存结果： $$ \sigma_n = \sqrt{\lambda^{(\text{final})}}, \quad v_{n,j} = v_{n,j}^{(\text{final})} $$

   • 所有节点同步 $n \gets n + 1$

分布式计算左奇异向量$u_{n,j}$

对于邻接矩阵 $A \in \mathbb{R}^{N \times N}$，其奇异值分解为：

• $U$ 的列向量 ${u_n}$ 是左奇异向量
• $V$ 的列向量 ${v_n}$ 是右奇异向量
• $\Sigma$ 是对角矩阵，元素 $\sigma_n$ 为奇异值

左奇异向量的定义关系：

输入：$\sigma_n$, $v_{n,j}$（来自幂迭代最终结果）

For $n = 0$ to $K-1$：

1. 本地计算： $$ u_{n,j} = \frac{1}{\sigma_n} \sum_{k \in 𝒩_j} a_{jk} v_{n,k} \quad \text{(需邻居节点发送$v_{n,k}$)} $$

2. 正交归一化：

   • For $m = 0$ to $n-1$： $$ u_{n,j} \gets u_{n,j} - \text{Consensus}\left(\sum_k u_{m,k} u_{n,k}\right) \cdot u_{m,j} $$

   • 归一化： $$ u_{n,j} \gets \frac{u_{n,j}}{\sqrt{\text{Consensus}\left(\sum_k u_{n,k}^2\right)}} $$

分布式重构邻接矩阵$A$

输入：$\sigma_n$, $u_{n,j}$, $v_{n,k}$

For 每个节点$j$并行执行：

1. 对每个邻居$k \in 𝒩_j$：

   • 请求节点$k$发送$v_{n,k}$（$n=0,...,K-1$）

   • 计算： $$ a_{jk} = \sum_{n=0}^{K-1} \sigma_n u_{n,j} v_{n,k} $$

2. 非邻居元素： $$ a_{jk} = 0 \quad \text{for} \quad k \notin 𝒩_j $$

非稳态下动态特征参数的估算

一致性控制策略

1. 异步更新模型

   • 节点仅在离散时刻 $t_k^i$ 接收邻居信息，更新自身状态 $x_i(t)$。
   • 各节点的状态更新时刻是独立的

2. 延时处理

   • 若检测到延时，节点选择最新收到的邻居状态替代旧值（避免使用过期数据）。

3. 一致性协议设计

   • 无时延系统 $$ \dot{x}i = \sum{j \in N(t_k^i, i)} a_{ij}(t_k^i) \left( x_j(t_k^i) - x_i(t) \right) $$

     参数	含义
     $\dot{x}_i$	节点 $i$ 的状态变化率（导数），表示 $x_i$ 随时间的变化速度。
     $x_i(t)$	节点 $i$ 在时刻 $t$ 的本地状态值（如特征估计、传感器数据等）。
     $x_j(t_k^i)$	节点 $i$ 在 $t_k^i$ 时刻收到的邻居节点 $j$ 的状态值。
     $N(t_k^i, i)$	节点 $i$ 在 $t_k^i$ 时刻的邻居集合（可直接通信的节点）。
     $a_{ij}(t_k^i)$	权重因子，控制邻居 $j$ 对节点 $i$ 的影响权重，满足 $\sum_j a_{ij} = 1$。

   • 有时延系统 $$ \dot{x}i = \sum{j \in N(t_k^i, i)} a_{ij}(t_k^i) \left( x_j(t_k^i - \tau_{ij}^k) - x_i(t) \right) $$

     参数	含义
     $\tau_{ij}^k$	节点 $j$ 到 $i$ 在时刻 $t_k^i$ 的信息传输延时。
     $t_k^i - \tau_{ij}^k$	节点 $i$ 实际使用的邻居状态 $x_j$ 的有效时刻（扣除延时）。

   • 权重$\alpha_{ij}$ $$ \text{有有效邻居时 } a_{ij}(t) = \begin{cases} \frac{\alpha_{ij}(t_k^i)}{\sum_{s \in N(t_k^i, i)} \alpha_{is}(t_k^i)}, & \text{若 } j \in N(t_k^i, i) \ 0, & \text{若 } j \notin N(t_k^i, i) \end{cases} $$

     $$ \text{无有效邻居时 } a_{ij}(t) = \begin{cases} 1, & \text{若 } j = i \ 0, & \text{若 } j \neq i \end{cases} $$

4. 通信拓扑定义

   • 引入 $G^0(t)$：实际成功通信的瞬时拓扑（非理想链路 $G(t)$），强调有效信息传递而非物理连通性。

收敛性分析

动态网络的收敛条件：

• 在节点移动导致的异步通信和随机延时下，只要网络拓扑满足有限时间内的联合连通性（即时间窗口内信息能传递到全网），所有节点的状态 $x_i$ 最终会收敛到同一全局值。 （平均代数连通度 > 0 == 动态网络拓扑的平均拉普拉斯矩阵的第二小特征值>0 ）

• 无需时刻连通：允许瞬时断连，但长期需保证信息能通过动态链路传递。

基于 UKF 的滤波估算

本文基于UKF：

• 采用**确定性采样（Sigma点）**直接近似非线性分布
• 完全规避对 f(x) 和 h(x) 的求导需求
• 保持高斯系统假设
• 允许函数不连续/不可微
• 适应拓扑突变等非线性情况

UKF 具体步骤

符号说明

• $i$: 节点索引，$N$ 为总节点数
• $x_i(k)$: 节点 $i$ 在时刻 $k$ 的状态分量 ($x$)
• $b_i(k)$: 节点 $i$ 的本地状态估计值 （相当于$Ax$）
• $a_{ij}$: 邻接矩阵元素（链路权重）
• $Q_k, R_k$: 过程噪声与观测噪声协方差
• $\mathcal{X}_{i,j}$: 节点 $i$ 的第 $j$ 个 Sigma 点
• $W_j^{(m)}, W_j^{(c)}$: Sigma 点权重（均值和协方差）

Step 1: 分布式初始化

1. 节点状态初始化：
   • 每个节点 $i$ 随机生成初始状态分量 $x_i(0)$。
   • 本地状态估计 $b_i(0)$ 初始化为 $x_i(0)$。

Step 2: 生成 Sigma 点（确定性采样）

在每个节点本地执行：

1. 计算 Sigma 点： $$ \begin{aligned} \mathcal{X}{i,0} &= \hat{b}{i,k-1} \ \mathcal{X}{i,j} &= \hat{b}{i,k-1} + \left( \sqrt{(n+\lambda) P_{i,k-1}} \right)j \quad (j=1,\dots,n) \ \mathcal{X}{i,j+n} &= \hat{b}{i,k-1} - \left( \sqrt{(n+\lambda) P{i,k-1}} \right)_j \quad (j=1,\dots,n) \end{aligned} $$

   • $\lambda = \alpha^2 (n + \kappa) - n$（缩放因子，$\alpha$ 控制分布范围，$\kappa$ 通常取 0）
   • $\sqrt{(n+\lambda) P}$ 为协方差矩阵的平方根（如 Cholesky 分解）

2. 计算 Sigma 点权重： $$ \begin{aligned} W_0^{(m)} &= \frac{\lambda}{n + \lambda} \quad &\text{(中心点均值权重)} \ W_0^{(c)} &= \frac{\lambda}{n + \lambda} + (1 - \alpha^2 + \beta) \quad &\text{(中心点协方差权重)} \ W_j^{(m)} = W_j^{(c)} &= \frac{1}{2(n + \lambda)} \quad (j=1,\dots,2n) \quad &\text{(对称点权重)} \end{aligned} $$

   • $\beta$ 为高阶矩调节参数（高斯分布时取 2 最优）

Step 3: 预测步骤（时间更新）

1. 传播 Sigma 点： $$ \mathcal{X}{i,j,k|k-1}^* = f(\mathcal{X}{i,j,k-1}) + q_k \quad (j=0,\dots,2n) $$

   • $f(\cdot)$ 为非线性状态转移函数
   • $q_k$ 为过程噪声 ，反映网络拓扑动态变化（如节点移动导致的链路扰动）。

2. 计算预测均值和协方差： $$ \hat{b}{i,k|k-1} = \sum{j=0}^{2n} W_j^{(m)} \mathcal{X}_{i,j,k|k-1}^* $$

   $$ P_{i,k|k-1} = \sum_{j=0}^{2n} W_j^{(c)} \left( \mathcal{X}{i,j,k|k-1}^* - \hat{b}{i,k|k-1} \right) \left( \mathcal{X}{i,j,k|k-1}^* - \hat{b}{i,k|k-1} \right)^T + Q_k $$

   • $Q_k$ 为过程噪声协方差

Step 4: 分布式观测生成

1. 邻居状态融合：

   • 节点 $ i $ 从邻居 $ j $ 获取其本地观测值 $ b_{j,\text{local}}(k) $：

   $$ b_{j,\text{local}}(k) = \sum_{l=1}^N a_{jl} x_l(k) \quad \text{(节点 $ j $ 对邻居状态的加权融合)} $$

   • 节点 $ i $ 综合邻居信息生成自身观测： $$ b_i^H(k) = \sum_{j=1}^N a_{ji} b_{j,\text{local}}(k) + r_k $$
   • 注：$ r_k $ 为通信噪声，反映信息传输误差（如延时、丢包）。

Step 5: 观测更新（测量更新）

1. 观测 Sigma 点： $$ \mathcal{Z}{i,j,k|k-1} = h(\mathcal{X}{i,j,k|k-1}^*) + r_k \quad (j=0,\dots,2n) $$

   • $h(\cdot)$ 为非线性观测函数

2. 计算观测统计量： $$ \hat{z}{i,k|k-1} = \sum{j=0}^{2n} W_j^{(m)} \mathcal{Z}_{i,j,k|k-1} $$

   $$ P_{i,zz} = \sum_{j=0}^{2n} W_j^{(c)} \left( \mathcal{Z}{i,j,k|k-1} - \hat{z}{i,k|k-1} \right) \left( \mathcal{Z}{i,j,k|k-1} - \hat{z}{i,k|k-1} \right)^T + R_k $$

   $$ P_{i,xz} = \sum_{j=0}^{2n} W_j^{(c)} \left( \mathcal{X}{i,j,k|k-1}^* - \hat{b}{i,k|k-1} \right) \left( \mathcal{Z}{i,j,k|k-1} - \hat{z}{i,k|k-1} \right)^T $$

   • $R_k$ 为观测噪声协方差

3. 计算卡尔曼增益并更新状态： $$ K_{i,k} = P_{i,xz} P_{i,zz}^{-1} $$

   $$ \hat{b}{i,k|k} = \hat{b}{i,k|k-1} + K_{i,k} \left( b_i^H(k) - \hat{z}_{i,k|k-1} \right) $$

   $$ P_{i,k|k} = P_{i,k|k-1} - K_{i,k} P_{i,zz} K_{i,k}^T $$

Step 6: 全局一致性计算

1. 瑞利商计算：

   • 所有节点通过一致性协议交换 $\hat{b}{i,k|k}$，计算全局状态： $$ y(k) = \frac{\sum{i=1}^N x_i(k) \hat{b}{i,k|k}}{\sum{i=1}^N x_i^2(k)} $$

2. 正交化：

   • 更新本地状态分量(相当于幂迭代$x=Ax$再归一化)： $$ x_i(k+1) = \frac{\hat{b}_{i,k|k}}{|\hat{b}(k)|_2} $$

Step 7: 收敛判断

• 若 $y(k)$ 收敛，输出 $\sigma = \sqrt{y(k)}$；否则返回 Step 2。

稳态下动态特征参数的估算

稳态下，网络拓扑变化趋于平稳，奇异值的理论曲线不再随时间变化（实际值因噪声围绕理论值波动）。此时采用集中式多观测值卡尔曼滤波

多观测值滤波算法

• 核心思想：利用相邻奇异值的有序性约束（$\sigma_{n-1} \leq \sigma_n \leq \sigma_{n+1}$），构造双观测值作为上下界，限制估计范围。

• 观测值生成：
  对第$n$大奇异值$\sigma_n$，其观测值$y_n$由相邻奇异值线性组合：
  $$ y_n = C_1 \sigma_{n-1} + C_2 \sigma_{n+1} $$

  • 系数$C_1, C_2$：根据$\sigma_{n-1}$和$\sigma_{n+1}$的权重动态调整（如距离比例）。
  • 物理意义：将$\sigma_{n-1}$和$\sigma_{n+1}$作为$\sigma_n$的下界和上界，避免单观测值因噪声导致的估计偏离。

疑问：

第三章的目的是什么？先分解再重构的意义在？

状态转移函数和观测函数怎么来？UKF每次预测单奇异值，如何同时预测K个呢？

卡尔曼滤波 观测值怎么来？是否需要拟合历史数据生成观测值？还是根据第三章分布式幂迭代求真实的特征值？