用 Kalman 秩判据检验可控性

# 系统 dx/dt = A x + B u，A[i, j] != 0 表示有向边 j -> i；B 的每一列是一个独立输入，作用在一个驱动节点上。
# 可控性矩阵 C = [B, AB, A^2 B, ..., A^{N-1} B]，满秩（rank = N）则可控。
# 检查三种情形：(a) 链 0->1->2->3 只驱动节点 0；(b) 星 0->1, 0->2, 0->3 只驱动节点 0；
# (c) 星形再加驱动节点 2、3。解读：(a) 可控，(b) 不可控（叶子只能收到同一信号的不同倍数），(c) 可控。
import numpy as np

def system_matrix(edges, n, seed=0):
    rng = np.random.default_rng(seed)
    A = np.zeros((n, n))
    for u, v in edges:
        A[v, u] = rng.uniform(0.5, 1.5)        # 随机非零权重
    return A

def input_matrix(drivers, n):
    B = np.zeros((n, len(drivers)))
    for j, d in enumerate(drivers):
        B[d, j] = 1.0
    return B

def ctrb_rank(A, B):
    n = A.shape[0]
    C = np.hstack([np.linalg.matrix_power(A, k) @ B for k in range(n)])
    return np.linalg.matrix_rank(C)

chain = [(0, 1), (1, 2), (2, 3)]
star = [(0, 1), (0, 2), (0, 3)]
cases = []
for name, edges, drv in [("链, 驱动{0}", chain, [0]), ("星, 驱动{0}", star, [0]), ("星, 驱动{0,2,3}", star, [0, 2, 3])]:
    r = ctrb_rank(system_matrix(edges, 4), input_matrix(drv, 4))
    cases.append((edges, drv, r))
    print("%-16s rank = %d / 4  -> %s" % (name, r, "可控" if r == 4 else "不可控"))

import networkx as nx
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 3, figsize=(10, 3.4))
titles = ["chain, drive {0}", "star, drive {0}", "star, drive {0,2,3}"]
for ax, (edges, drv, r), title in zip(axes, cases, titles):
    G = nx.DiGraph(edges)
    pos = {0: (0, 0), 1: (1, 0), 2: (2, 0), 3: (3, 0)} if edges is chain else {0: (0, 0), 1: (1, 1), 2: (1, 0), 3: (1, -1)}
    nx.draw_networkx_edges(G, pos, ax=ax, edge_color="0.5", arrows=True, arrowsize=14, node_size=400)
    nx.draw_networkx_nodes(G, pos, ax=ax, node_size=400, node_color=["tab:red" if v in drv else "0.85" for v in G])
    nx.draw_networkx_labels(G, pos, ax=ax, font_size=9)
    ax.set_title("%s\nrank = %d / 4 (%s)" % (title, r, "controllable" if r == 4 else "not controllable"), fontsize=10)
    ax.axis("off"); ax.margins(0.2)
plt.tight_layout(); plt.show()