# 陈茂森论文 ## 随机移动网络系统的稳定性 ### 马尔科夫链与网络平均度推导 #### **1.马尔科夫链的基本概念** 马尔科夫链描述的是这样一种随机过程：系统在若干个可能的状态中变化，**下一时刻所处状态只依赖于当前状态**，而与过去的状态无关，这就是所谓的“无记忆性”或**马尔科夫性**。 **无记忆性**意味着，对于任何 $s, t \ge 0$， $$ P(T > s+t \mid T > s) = P(T > t). $$ 假设你已经等待了 $s$ 分钟，那么再等待至少 $t$ 分钟的概率，和你一开始就等待至少 $t$ 分钟的概率完全相同。 在所有概率分布里，只有指数分布 $$ P(T>t) = e^{-\lambda t} $$ 具有这种“无记忆性”特征： $$ P(T>s+t \mid T>s) = \frac{P(T>s+t)}{P(T>s)} = \frac{e^{-\lambda(s+t)}}{e^{-\lambda s}} = e^{-\lambda t} = P(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)$ 的微分方程写成： $$ \frac{d p_1(t)}{dt} = \lambda_{01} \, (1-p_1(t)) - \lambda_{10} \, p_1(t). $$ 这个方程可以整理为： $$ \frac{d p_1(t)}{dt} + (\lambda_{01}+\lambda_{10}) \, p_1(t) = \lambda_{01}. $$ 这其实是一个一阶线性微分方程，其标准求解方法是求解其齐次解与非齐次解。 #### **3. 求解微分方程** 整个微分方程的通解为： $$ p_1(t)= \frac{\lambda_{01}}{\lambda_{01}+\lambda_{10}} + C\, e^{-(\lambda_{01}+\lambda_{10})t}. $$ 利用初始条件 $p_1(0)=p_1^0$（初始时刻链路连通的概率），我们可以求出 $C$： 即 $$ C = p_1^0 -\frac{\lambda_{01}}{\lambda_{01}+\lambda_{10}}. $$ 所以，链路在任意时刻 $t$ 连通的概率为： $$ p_1(t)= \frac{\lambda_{01}}{\lambda_{01}+\lambda_{10}} + \left( p_1^0 -\frac{\lambda_{01}}{\lambda_{01}+\lambda_{10}} \right)e^{-(\lambda_{01}+\lambda_{10})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)$ 表达式代入，就得到网络平均度随时间变化的表达式： $$ \bar{d}(t)= (N-1)\left[ \frac{\lambda_{01}}{\lambda_{01}+\lambda_{10}} + \left( p_1^0 -\frac{\lambda_{01}}{\lambda_{01}+\lambda_{10}} \right)e^{-(\lambda_{01}+\lambda_{10})t}\right]. $$ 这就是**网络平均度的变化函数**： - 网络开始时每条链路的连通概率为 $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)$，其分布概率密度函数写为 $$ f(x,y). $$ 那么节点在模型子区域 $R_1$ 中出现的概率为 $$ P_{R_1}=\int_{R_1} f(x,y) \,dx\,dy. $$ 在平衡状态下，理论上节点的位置分布 **$f(x,y)$ 保持不变**，即每个区域内节点出现的概率 $P_{R_1}$ 是常数，不随时间变化。证明在平衡状态下节点分布稳定。 #### **扰动后的恢复能力** - 静止节点分布特性满足均匀分布,概率密度函数为 $g(x,y)$； - 运动节点的概率密度函数为 $h(x,y)$； - **在时刻 $t_0$ 时**，网络中共有 $N$ 个节点，其中有 $s$ 个静止，故**静止节点的比例为 $p=\frac{s}{N}$**。 节点整体的分布概率密度函数可写为 $$ f(x,y)=p\, g(x,y)+(1-p)\, h(x,y). $$ 在平衡状态下，$p$ 的理论值为一个常数，所以 $f(x,y)$ 不随时间变化，从而网络连通度稳定。 接下来，**考虑外界扰动**的影响：假设在**时刻 $t_1$** 新加入 $m$ 个**符合均匀分布的节点**， **扰动后的总分布（$t_1$时刻后）** - 新加入的 $m$ 个节点是静止的，其分布为 $g(x,y)$ - 此时网络的总节点数 $N+m$： - **静止节点总数**：$s+m$ - **运动节点总数**：$N-s$（原有运动节点数不变） 因此，扰动后的分布为： $$ f(x,y,t_1) = \frac{s + m}{N + m} g(x,y) + \frac{N - s}{N + m} h(x,y) $$ $$ f(x,y,t_1) = p' \cdot g(x,y) + (1-p') \cdot h(x,y) $$ 其中 $p' = \frac{s + m}{N + m}$。 **近似处理（当 $N, s \gg m$ 时）** $$ p' = \frac{s + m}{N + m} \approx \frac{s}{N} = p $$ 因此，扰动后的分布近似为： $$ f(x,y,t_1) \approx p \cdot g(x,y) + (1-p) \cdot h(x,y) $$ 这与初始平衡态的分布相同，说明网络在扰动后恢复了平衡态。 ### 系统稳定性分析 #### 平衡点及误差坐标 论文第 2.1 节推导出，单条链路连通概率 $p_1(t)$ 满足 $$ \dot p_1(t) = -(\lambda_{01}+\lambda_{10})\,p_1(t) \;+\;\lambda_{01}. \tag{2‑18} $$ 网络有 $N$ 个节点，**平均度** $$ d(t) = (N-1)\,p_1(t). $$ 设平衡连通概率 $p_1^*$ 满足 $\dot p_1=0$，解得 **平衡点** $$ p_1^* = \frac{\lambda_{01}}{\lambda_{01}+\lambda_{10}}, \quad d^* \;=\;(N-1)\,p_1^*. $$ **定义误差（偏离平衡的量）** $$ e(t)=d(t)-d^*. $$ 其中 $$ d(t)=(N-1)\,p_1(t), \qquad d^*=(N-1)\,p_1^*. $$ 将 $d$ 和 $d^*$ 代入 $$ e(t) =d(t)-d^* =(N-1)\,p_1(t)\;-\;(N-1)\,p_1^* =(N-1)\,\bigl[p_1(t)-p_1^*\bigr]. $$ 解得 $$ p_1(t)-p_1^* \;=\;\frac{e(t)}{\,N-1\,} \quad\Longrightarrow\quad p_1(t) =\frac{e(t)}{\,N-1\,}+p_1^*. $$ **误差求导** $$ \dot e(t) = \frac{d}{dt}\bigl[(N-1)(p_1-p_1^*)\bigr] = (N-1)\,\dot p_1(t), $$ 得到 $$ \dot e =(N-1)\Bigl[-(\lambda_{01}+\lambda_{10})\Bigl(\tfrac{e}{N-1}+p_1^*\Bigr) +\lambda_{01}\Bigr].\\\dot e = -(\lambda_{01}+\lambda_{10})\,e. $$ 这就是把原来以 $p_1$ 为自变量的微分方程，转写成以 "偏离平衡量" $e$ 为自变量的形式 记常数 $$ c = \lambda_{01}+\lambda_{10} >0, $$ 则误差模型就是一维线性常微分方程： $$ \dot e = -\,c\,e. $$ #### 构造李雅普诺夫函数 对一维系统 $\dot e=-ce$（$c>0$），自然选取 $$ V(e)=e^2 $$ 作为李雅普诺夫函数，理由是： - $V(e)>0$ 当且仅当 $e\neq0$； - 平衡点 $e=0$ 时，$V(0)=0$。 #### 计算 $V$ 的时间导数 对 $V$ 关于时间求导： $$ \dot V(e) = \frac{d}{dt}\bigl(e^2\bigr) = 2\,e\,\dot e = 2\,e\,\bigl(-c\,e\bigr) = -2c\,e^2. $$ 因为 $c>0$ 且 $e^2\ge0$，所以 $$ \boxed{\dot V(e)\;=\;-2c\,e^2\;\le\;0.} $$ - 当 $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$ ### 瑞利商公式 目的：求对称矩阵A的最大特征值 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矩阵 $$ A = \begin{pmatrix}2 & 1\\1 & 2\end{pmatrix}, \quad x = \begin{pmatrix}2\\1\end{pmatrix}. $$ 1. **集中式计算** $$ y= \frac{x^T A x}{x^T x} = \frac{\begin{pmatrix}2 & 1\end{pmatrix} \begin{pmatrix}5\\4\end{pmatrix}}{2^2 + 1^2} = \frac{10 + 4}{5} = \frac{14}{5} = 2.8. $$ 2. **分布式计算** 各节点分别计算本地观测值 ​ 节点1的计算： $$ b_1 = a_{11}x_1 + a_{12}x_2 = 2 \cdot 2 + 1 \cdot 1 = 5. $$ ​ 节点2的计算： $$ b_2 = a_{21}x_1 + a_{22}x_2 = 1 \cdot 2 + 2 \cdot 1 = 4. $$ ​ 然后通过全网共识计算 $$ y = \frac{2 \cdot 5 + 1 \cdot 4}{2^2 + 1^2} = \frac{14}{5} = 2.8. $$ ### 主要符号表 | 符号 | 类型 | 含义 | 存储/计算位置 | | --------------- | -------- | ----------------------------------------------- | --------------- | | $n$ | 下标 | 当前计算的奇异值序号（从0开始） | 全局共识 | | $K$ | 常量 | 需要计算的前$K$大奇异值总数 | 预设参数 | | $j,k$ | 下标 | 节点编号（$j$表示当前节点） | 本地存储 | | $𝒩_j$ | 集合 | 节点$j$的邻居节点集合 | 本地拓扑信息 | | $a_{jk}$ | 矩阵元素 | 邻接矩阵$A$中节点$j$与$k$的连接权值 | 节点$j$本地存储 | | $v_{n,j}^{(t)}$ | 向量分量 | 第$n$个右奇异向量在节点$j$的分量（第$t$次迭代） | 节点$j$存储 | | $u_{n,j}$ | 向量分量 | 第$n$个左奇异向量在节点$j$的分量 | 节点$j$计算存储 | | $\sigma_n$ | 标量 | 第$n$个奇异值 | 全局共识存储 | | $\delta$ | 标量 | 收敛阈值 | 预设参数 | ### 分布式幂迭代求前$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}$，其奇异值分解为： $$ A = U \Sigma V^T $$ 其中： - $U$ 的列向量 $\{u_n\}$ 是左奇异向量 - $V$ 的列向量 $\{v_n\}$ 是右奇异向量 - $\Sigma$ 是对角矩阵，元素 $\sigma_n$ 为奇异值 **左奇异向量的定义关系**： $$ A v_n = \sigma_n u_n \quad \Rightarrow \quad u_n = \frac{1}{\sigma_n} A v_n $$ 展开为分量形式（对第 $j$ 个分量）： $$ u_{n,j} = \frac{1}{\sigma_n} \sum_{k=1}^N a_{jk} v_{n,k} $$ **输入**：$\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 的滤波估算 | | KF | EKF | 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个呢？ 卡尔曼滤波 观测值怎么来？是否需要拟合历史数据生成观测值？还是根据第三章分布式幂迭代求真实的特征值？