生物医学网络科学Biomedical Network Science
网络可控性与稳定性14 / 17

数学推导与算法

阅读
约 5 分钟

本页给出：A. 最小输入定理的来源与 Hopcroft–Karp 算法；B. 随机失效阈值 fcf_c 的推导与 Newman–Ziff 算法（更系统的渗流推导见网络鲁棒性、渗流与级联失效）；C. PBH 判据与精确可控性的证明要点；D. Gramian、Lyapunov 方程与平均/模态可控性；E. 目标控制的贪心算法；F. 观测节点必要条件的论证；G. 反馈顶点集的求解。Kalman 秩判据与可控性 Gramian 的推导见线性系统的可控性；结构可控性“几乎所有权重”的论证、结构可控与给定参数可控的区别、最大匹配框架与最小支配集框架的比较见结构可控性与最大匹配。

A. 最小输入定理的来源

Lin（1974）的结构可控性定理：(A,B)(A,B) 结构可控，当且仅当（1）无不可达节点：每个状态节点都能从某个输入出发经有向路径到达；（2）无扩张（dilation）：不存在节点子集 SS，使指向 SS 的节点（含输入）个数少于 |S||S|。

Liu、Slotine 与 Barabási（2011）的论证要点：

  1. 由 Hall 定理，条件（2）等价于：在把输入也算进去的二分图中，每个状态节点的“入端”都能被匹配。一个独立输入只能匹配一个节点，所以输入数 ≥N−|M*|\ge N-|M^*|。
  2. 反过来，在每个非匹配节点上加一个独立输入。最大匹配的边把网络分解成若干条以非匹配节点为起点的有向路径（“茎”）和若干个有向环。每个环可以通过一条非匹配边挂到某条茎上；若某个环从任何茎都不可达，就把已有的某个输入信号同时接到环上一个节点——这不增加独立输入的个数。得到的“仙人掌”结构满足 Lin 的两个条件。
  3. 若存在完美匹配，没有非匹配节点，但至少需要 1 个输入，故 ND=max⁡(N−|M*|,1)N_D=\max(N-|M^*|,1)。
算法 1  Hopcroft–Karp 最大匹配与驱动节点
输入：有向图 G=(V,E)；构造二分图：左侧 u⁺，右侧 v⁻，有向边 u→v 对应 u⁺–v⁻
1. 所有 matchL[u⁺] ← NIL，matchR[v⁻] ← NIL
2. while BFS() = true:                 # 每一轮称为一个"阶段"
3.     for 每个未匹配的左顶点 u⁺: DFS(u⁺)
4. D ← { v : matchR[v⁻] = NIL }；若 D 为空，任取一个节点放入 D
5. return 匹配数 |M*|（matchL 中非 NIL 的个数）与驱动节点集 D

BFS()：
  a. 所有未匹配左顶点 dist ← 0 入队，其余 dist ← ∞
  b. 出队 u⁺；对每个邻居 v⁻：w ← matchR[v⁻]
       若 w = NIL：找到最短增广路，记录其长度（不再扩展更深层）
       否则若 dist[w] = ∞：dist[w] ← dist[u⁺]+1，w 入队
  c. 返回是否找到增广路
DFS(u⁺)：
  for 每个邻居 v⁻：w ← matchR[v⁻]
     if w = NIL 或 (dist[w] = dist[u⁺]+1 且 DFS(w))：
         matchL[u⁺] ← v⁻；matchR[v⁻] ← u⁺；return true
  dist[u⁺] ← ∞；return false

复杂度：每个阶段 BFS 与 DFS 合计 O(L)O(L)；可以证明阶段数为 O(N)O(\sqrt N)，总计 O(NL)O(\sqrt N\,L)。注意第 4 步给出的只是一个最小驱动节点集；要得到“所有可能的驱动节点”，需要在终止时沿交替路径继续搜索（见驱动节点的枚举、分类与自适应控制中的输入图框架）。

B. Molloy–Reed 判据与 fcf_c

假设：无度相关的配置模型、N→∞N\to\infty、局部树状（忽略短环）。

第 1 步：分支过程。 沿一条随机边到达一个节点，其度为 kk 的概率为 kP(k)/⟨k⟩kP(k)/\langle k\rangle，它还能继续“向外”延伸 k−1k-1 条边。巨分支存在，当且仅当这一分支过程的平均后代数大于 1： ∑k(k−1)kP(k)⟨k⟩=⟨k2⟩−⟨k⟩⟨k⟩>1⇔κ=⟨k2⟩⟨k⟩>2.\sum_k(k-1)\frac{kP(k)}{\langle k\rangle}=\frac{\langle k^2\rangle-\langle k\rangle}{\langle k\rangle}>1 \iff\kappa=\frac{\langle k^2\rangle}{\langle k\rangle}>2 .

第 2 步：随机删点后的度分布。 删除比例 ff 的节点后，保留下来的节点的每条边以概率 p=1−fp=1-f 独立保留（另一端未被删）。原度为 kk 的节点新度 k′∼Binomial(k,p)k'\sim\text{Binomial}(k,p)，于是 ⟨k′⟩=p⟨k⟩,⟨k′2⟩=p2⟨k2⟩+p(1−p)⟨k⟩.\langle k'\rangle=p\langle k\rangle,\qquad\langle k'^2\rangle=p^2\langle k^2\rangle+p(1-p)\langle k\rangle .

第 3 步：临界点。 新的 κ′=⟨k′2⟩⟨k′⟩=pκ0+1−p\kappa'=\dfrac{\langle k'^2\rangle}{\langle k'\rangle}=p\,\kappa_0+1-p。令 κ′=2\kappa'=2，得 pc(κ0−1)=1p_c(\kappa_0-1)=1，即 fc=1−pc=1−1κ0−1.f_c=1-p_c=1-\frac{1}{\kappa_0-1}. 蓄意攻击（按度删除）同时改变 kmaxk_{\max} 和余下节点的度分布，推导更复杂，见 Cohen 等（2001，PRL）。

算法 2  Newman–Ziff：一次算出整条随机失效曲线
1. 生成节点的随机排列 π（看作删除顺序的逆序）
2. 并查集初始为空；S_max ← 0
3. for t = 1 … N:                     # 按 π 逐个"加回"节点
4.     加入节点 π_t；对其已加入的每个邻居做 union，同步更新集合大小
5.     S_max ← max(S_max, 新集合大小)
6.     记录 P∞(f = 1 − t/N) = S_max / N
7. 对多个随机排列取平均

复杂度：每条曲线 O(Lα(N))O(L\,\alpha(N))，α\alpha 为反阿克曼函数，实际上近似线性（Newman 与 Ziff 2000，PRL）；比“每个 ff 重新删点、重新求连通分支”的 O(N⋅L)O(N\cdot L) 快得多。

C. PBH 判据与精确可控性

PBH 判据的证明要点。（PBH 条件的必要性）若存在左特征向量 𝒒⊤A=λ𝒒⊤\mathbf q^\top A=\lambda\mathbf q^\top 且 𝒒⊤B=𝟎\mathbf q^\top B=\mathbf 0，则 𝒒⊤AkB=λk𝒒⊤B=𝟎\mathbf q^\top A^kB=\lambda^k\mathbf q^\top B=\mathbf 0 对所有 kk 成立，𝒒⊤𝒞=𝟎\mathbf q^\top\mathcal C=\mathbf 0，不可控。（PBH 条件的充分性，逆否形式）若 rank⁡𝒞<N\operatorname{rank}\mathcal C<N，不可控子空间的正交补 {𝒒:𝒒⊤𝒞=0}\{\mathbf q:\mathbf q^\top\mathcal C=0\} 非空，且在 A⊤A^\top 作用下不变（由 Cayley–Hamilton 定理，𝒒⊤𝒞=0⇒𝒒⊤A𝒞=0\mathbf q^\top\mathcal C=0\Rightarrow\mathbf q^\top A\mathcal C=0）。A⊤A^\top 在这个不变子空间中至少有一个特征向量 𝒒\mathbf q，它满足 𝒒⊤A=λ𝒒⊤\mathbf q^\top A=\lambda\mathbf q^\top 且 𝒒⊤B=𝟎\mathbf q^\top B=\mathbf 0。rank⁡[λI−A,B]<N\operatorname{rank}[\lambda I-A,B]<N 与“存在这样的 𝒒\mathbf q”等价，因为 𝒒⊤[λI−A,B]=𝟎\mathbf q^\top[\lambda I-A,B]=\mathbf 0。

精确可控性的计数论证。对特征值 λ\lambda，λI−A\lambda I-A 的左零空间维数为几何重数 μ(λ)\mu(\lambda)。要让 [λI−A,B][\lambda I-A,B] 满秩，BB 的列必须在这 μ(λ)\mu(\lambda) 维左零空间上的投影满秩，因此 rank⁡B≥μ(λ)\operatorname{rank}B\ge\mu(\lambda)；对所有特征值同时成立，ND≥max⁡λμ(λ)N_D\ge\max_\lambda\mu(\lambda)。反过来，可以构造恰好 max⁡λμ(λ)\max_\lambda\mu(\lambda) 列的 BB 同时满足所有特征值的条件（Yuan 等 2013 用初等变换给出构造）。

稀疏网络的近似。若最大重数在 λ=0\lambda=0 处取得，μ(0)=N−rank⁡A\mu(0)=N-\operatorname{rank}A，即 ND=max⁡{1,N−rank⁡A}N_D=\max\{1,N-\operatorname{rank}A\}。对稀疏有向网络，AA 的通用秩等于最大匹配的边数 |M*||M^*|（这是 König 定理与结构秩的经典结论），于是 N−rank⁡A=N−|M*|N-\operatorname{rank}A=N-|M^*|，与最小输入定理一致——这解释了为什么在稀疏网络中两种框架几乎给出相同的 NDN_D。

D. Gramian、Lyapunov 方程与平均/模态可控性

Lyapunov 方程。设 AA 稳定，W∞=∫0∞eAτBB⊤eA⊤τdτW_\infty=\int_0^\infty e^{A\tau}BB^\top e^{A^\top\tau}d\tau。对被积函数求导： ddτ(eAτBB⊤eA⊤τ)=AeAτBB⊤eA⊤τ+eAτBB⊤eA⊤τA⊤.\frac{d}{d\tau}\big(e^{A\tau}BB^\top e^{A^\top\tau}\big)=A\,e^{A\tau}BB^\top e^{A^\top\tau}+e^{A\tau}BB^\top e^{A^\top\tau}A^\top . 从 0 积分到 ∞\infty，左边为 0−BB⊤0-BB^\top，右边为 AW∞+W∞A⊤AW_\infty+W_\infty A^\top，即 AW∞+W∞A⊤+BB⊤=0AW_\infty+W_\infty A^\top+BB^\top=0。离散时间的 W=∑τAτBB⊤(A⊤)τW=\sum_\tau A^\tau BB^\top(A^\top)^\tau 满足 W−AWA⊤=BB⊤W-AWA^\top=BB^\top（把求和的首项单独写出即得）。两者都可用 Bartels–Stewart 等算法在 O(N3)O(N^3) 内求解（SciPy 的 solve_continuous_lyapunov、solve_discrete_lyapunov）。

平均可控性的谱表示。设离散时间对称 A=VΛV⊤A=V\Lambda V^\top，B=𝒆iB=\mathbf e_i，则 tr⁡Wi=∑τ≥0𝒆i⊤A2τ𝒆i=∑jvij2∑τ≥0λj2τ=∑jvij21−λj2.\operatorname{tr}W_i=\sum_{\tau\ge0}\mathbf e_i^\top A^{2\tau}\mathbf e_i=\sum_j v_{ij}^2\sum_{\tau\ge0}\lambda_j^{2\tau}=\sum_j\frac{v_{ij}^2}{1-\lambda_j^2}. 平均可控性对慢模态（|λj||\lambda_j| 接近 1）的分量赋予大权重；模态可控性 ϕi=∑j(1−λj2)vij2\phi_i=\sum_j(1-\lambda_j^2)v_{ij}^2 则对快模态赋予大权重。两者权重互为倒数，这从谱的角度解释了为什么二者与加权度的相关方向相反：度高的节点在主导（慢）特征向量上的分量通常较大。

能量下界。单位方向的最小能量 E=𝒅⊤W−1𝒅E=\mathbf d^\top W^{-1}\mathbf d 满足 Rayleigh 商不等式 1/ηmax≤E≤1/ηmin1/\eta_{\max}\le E\le1/\eta_{\min}；沿随机单位方向的平均能量为 tr⁡(W−1)/N\operatorname{tr}(W^{-1})/N，并且由算术–调和平均不等式，tr⁡(W−1)≥N2/tr⁡W\operatorname{tr}(W^{-1})\ge N^2/\operatorname{tr}W。

E. 目标控制的贪心算法

算法 3  目标控制的贪心算法（Gao 等 2014）
输入：有向图 G=(V,E)，目标集合 T
1. D ← ∅；S ← T；已出现集合 H ← ∅
2. while S 非空 且 S ∉ H:
3.     H ← H ∪ {S}
4.     构造二分图：右侧为 S，左侧为 S 中节点的所有上游邻居 P(S)，边 u→v（u∈P(S), v∈S）
5.     求最大匹配 M
6.     D ← D ∪ {S 中未被 M 匹配的节点}          # 第一轮的未匹配数即下界
7.     S ← {被 M 匹配到的左侧节点}              # 上游节点成为新的"目标"
8. 若 D 为空（目标全部位于环上），任取 T 中一个节点放入 D
9. return D

每轮一次 Hopcroft–Karp，复杂度 O(NL)O(\sqrt N L)；轮数不超过网络中最长的“匹配链”。第 2 步的“S∉HS\notin H”用于在含环网络中终止迭代（环上的匹配会使 SS 重复出现，此时环可由已有输入驱动）。得到的 DD 要用目标可控性秩判据在随机权重下检验（见代码示例）。

F. 观测节点必要条件的论证

设推断图中 RR 是一个根 SCC（没有来自 RR 外的入边）。在推断图中，入边 xi→xjx_i\to x_j 表示“xjx_j 出现在 xix_i 的方程中”；RR 没有来自外部的入边，意味着 RR 之外任何变量的方程都不含 RR 中的变量。因此，若不直接测量 RR 中的节点，所有输出及其各阶导数都与 RR 中变量无关，RR 的状态无法从输出中恢复。所以每个根 SCC 至少需要一个观测节点。对线性系统，还需要消除“扩张”造成的对称性，充分条件由反向网络的最大匹配给出。

G. 反馈顶点集的求解

最小 FVS 是 NP 难问题，但实际网络中可用以下步骤大幅缩小规模：（1）反复删除入度或出度为 0 的节点（它们不在任何环上）；（2）自环节点必须进入 FVS；（3）对只有一条入边或一条出边的节点做“收缩”；（4）对剩余的核心用整数规划、分支定界或模拟退火等启发式求近似最小解（Zañudo 等 2017 采用的就是这类近似算法）。注意 FVS 也不唯一，不同的最小 FVS 对应不同的候选控制集合。

继续阅读