求解Dynamic Riccati Equation得到零矩阵解的问题排查求助

核心错误：矩阵连续乘法的写法错误 #

你的代码中，np.matmul(sigma, Md, sigma)是错误用法——np.matmul仅接受两个矩阵作为输入，第三个参数是用于指定轴的可选参数，并非第三个待相乘的矩阵。这会导致Riccati方程中的二次项计算完全错误，进而得到不符合预期的零矩阵解。

正确的三个矩阵连续乘法写法应该是嵌套调用np.matmul，或者更简洁地使用Python矩阵乘法运算符@：

# 嵌套matmul写法
2 * np.matmul(np.matmul(sigma, Md), sigma)
# 或者用@运算符（更直观）
2 * (sigma @ Md @ sigma)

其他优化点 #

1. 避免全局变量依赖：deriv函数直接使用全局变量s不够鲁棒，建议在函数内部根据输入的一维sigma计算矩阵形状：

   def deriv(sigma, t, A, B, Md):
       n = A.shape[0]
       sigma_mat = sigma.reshape(n, n)
       return (-A @ sigma_mat - sigma_mat @ A.T + B - 2 * (sigma_mat @ Md @ sigma_mat)).flatten()
   

2. 初始条件与收敛性：全1矩阵作为初始条件时，若积分时间不足可能无法收敛到稳态。可以尝试延长积分时间，或者检查矩阵A的稳定性——若A是不稳定矩阵，正向积分可能无法收敛，此时需要考虑反向积分（将时间取反，从终值条件往回积分）。

修正后的完整代码 #

import numpy as np
from scipy.integrate import odeint
import networkx as nx
import matplotlib.pyplot as plt

def deriv(sigma, t, A, B, Md):
    n = A.shape[0]
    sigma_mat = sigma.reshape(n, n)
    # 修正矩阵乘法，并使用@运算符简化写法
    dsigma = -A @ sigma_mat - sigma_mat @ A.T + B - 2 * (sigma_mat @ Md @ sigma_mat)
    return dsigma.flatten()

n = 100
alpha = 5
beta = 5

G = nx.barabasi_albert_graph(n=n, m=2, seed=10374196, initial_graph=None)
A = nx.adjacency_matrix(G).toarray()  # Convert to dense matrix
B = alpha * np.identity(n)
Md = np.identity(n)
Md[1, 1] = beta

sigma0 = np.ones(n**2)  # 全1初始矩阵
t = np.linspace(0, 200, 201)  # 延长积分时间，确保收敛

sol = odeint(deriv, sigma0, t, args=(A, B, Md)) 
sol = np.reshape(sol, (len(t), n, n))  

# 获取最终稳态解的特征值并排序
Y = np.linalg.eigvals(sol[-1])  
Y.sort()
X = np.arange(n)

plt.plot(X, Y)
plt.xlabel("eigenvalue rank")
plt.ylabel("eigenvalue")
plt.show()

验证稳态解 #

修正后，你可以通过检查最终解是否满足稳态方程来验证：

sigma_ss = sol[-1]
residual = A @ sigma_ss + sigma_ss @ A.T + 2 * (sigma_ss @ Md @ sigma_ss) - B
print("稳态方程残差的范数：", np.linalg.norm(residual))

若残差范数接近0，说明解已收敛到正确的稳态。

内容的提问来源于stack exchange，提问作者Destroctive Blade