PBH 判据与精确可控性——同一拓扑、不同权重需要的输入数不同

# PBH 判据：(A,B) 可控 <=> 对 A 的每个特征值 lambda，rank[lambda I - A, B] = N。
# Yuan 等（2013）的精确可控性：给定 A 时最少独立输入数 N_D = max_lambda [N - rank(lambda I - A)]（最大几何重数）。
# 比较 (a) 全 1 权重的完全有向图（N_D = N-1）、(b) 同一拓扑的随机权重（N_D = 1，与最大匹配一致）、
# (c) 稀疏 ER 网络：用 N - rank(A) 近似（Yuan 2013 式 (5)），并与最大匹配给出的 N_D 对比。
# 解读：结构可控性给的是“几乎所有权重”下的答案；权重相等这类约束会把系统推进例外集合。
import numpy as np
import networkx as nx

def exact_ND(A, tol=1e-8):
    """最大几何重数：对每个（数值上不同的）特征值计算 N - rank(lambda I - A)。"""
    N = A.shape[0]
    lams = np.linalg.eigvals(A)
    best, seen = 1, []
    for lam in lams:
        if any(abs(lam - s) < 1e-6 for s in seen):
            continue
        seen.append(lam)
        r = np.linalg.matrix_rank(lam * np.eye(N) - A, tol=tol * max(1.0, np.abs(A).max()) * N)
        best = max(best, N - r)
    return best

def matching_ND(G):
    Bg = nx.Graph()
    top = [("o", u) for u in G]
    Bg.add_nodes_from(top); Bg.add_nodes_from([("i", v) for v in G])
    Bg.add_edges_from((("o", u), ("i", v)) for u, v in G.edges())
    M = nx.bipartite.hopcroft_karp_matching(Bg, top_nodes=top)
    return max(1, sum(1 for v in G if ("i", v) not in M))

N = 6
A_unit = np.ones((N, N)) - np.eye(N)                      # 完全图，权重全为 1
rng = np.random.default_rng(1)
A_rand = A_unit * rng.uniform(0.5, 1.5, size=(N, N))      # 同一拓扑，独立随机权重
print("完全图 N=%d：全 1 权重 N_D = %d；随机权重 N_D = %d；最大匹配 N_D = %d"
      % (N, exact_ND(A_unit), exact_ND(A_rand), matching_ND(nx.complete_graph(N, create_using=nx.DiGraph))))

# PBH 检查：全 1 权重时，只驱动 1 个节点在 lambda = -1 处秩不足
B1 = np.zeros((N, 1)); B1[0, 0] = 1
print("全 1 权重、驱动 1 个节点：rank[-I-A, B] = %d / %d"
      % (np.linalg.matrix_rank(np.hstack([-np.eye(N) - A_unit, B1])), N))

G = nx.gnm_random_graph(300, 450, directed=True, seed=3)
A = nx.to_numpy_array(G).T                                 # A[v,u] = 1 表示 u -> v
print("稀疏有向 ER（N=300, L=450，单位权重）：N - rank(A) = %d；最大几何重数 = %d；最大匹配 N_D = %d"
      % (300 - np.linalg.matrix_rank(A), exact_ND(A), matching_ND(G)))

import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
ev_u, ev_r = np.linalg.eigvals(A_unit), np.linalg.eigvals(A_rand)
ax1.plot(ev_u.real, ev_u.imag, "o", ms=12, mfc="none", mew=2, label="unit weights (N_D = %d)" % exact_ND(A_unit))
ax1.plot(ev_r.real, ev_r.imag, "x", ms=8, mew=2, label="random weights (N_D = %d)" % exact_ND(A_rand))
ax1.annotate("lambda = -1, multiplicity %d" % (N - 1), (-1, 0), xytext=(0.2, 0.3), arrowprops=dict(arrowstyle="->"))
ax1.set_ylim(-0.6, 0.6)
ax1.set_xlabel("Re lambda"); ax1.set_ylabel("Im lambda"); ax1.legend(fontsize=8)
ax1.set_title("Complete digraph N = %d: degenerate spectrum" % N)
ev = np.linalg.eigvals(A)
ax2.plot(ev.real, ev.imag, ".", c="tab:blue", alpha=0.6)
ax2.set_xlabel("Re lambda"); ax2.set_ylabel("Im lambda")
ax2.set_title("Sparse ER (unit weights): %d of %d eigenvalues near 0" % ((np.abs(ev) < 1e-6).sum(), len(ev)), fontsize=10)
plt.tight_layout(); plt.show()