使用scipy.linalg.solve_discrete_are求解仿射动态系统离散代数Riccati方程时遇有限解不存在错误的问题排查

嘿，我来帮你拆解这个问题！你碰到的LinAlgError: Failed to find a finite solution本质是系统不满足离散Riccati方程求解的核心条件，咱们一步步理清楚：

一、根本原因：系统不可镇定 #

scipy.linalg.solve_discrete_are要求你的扩展系统**(A_new, B_new)是可镇定的**——简单来说，就是所有模大于1的不稳定特征值，对应的系统模态必须是可控的。咱们来验证你的系统：

1. A_new的特征值：原A是2x2全1矩阵，特征值是2（模>1，属于不稳定特征值）和0，加上扩展后新增的特征值1，所以A_new的特征值是2、0、1。
2. 不稳定特征值的可控性检查：对于特征值2，咱们看矩阵(A_new - 2*I)和B_new拼接后的秩：

   A_new_minus_2I = A_new - 2*np.eye(3)
   controllability_check = np.hstack([A_new_minus_2I, B_new])
   print(np.linalg.matrix_rank(controllability_check))  # 输出是2，小于系统状态维度3
   

   这说明特征值2对应的模态是不可控的，系统无法通过控制输入稳定这个模态，所以Riccati方程找不到有限解。

另外你构造的Q矩阵也有问题：np.eye(x_dim+1)会惩罚扩展出来的那个固定为1的状态，但这个状态是咱们人为添加的常数项，根本不需要被代价函数惩罚，这也会干扰求解过程。

二、解决办法 #

1. 先确保原系统可镇定

如果你的实际系统就是当前的A和B，那它本身不可镇定，得调整A或B让不稳定模态可控。比如把B改成：

B = np.array([[1], [0]])  # 此时可控性矩阵秩为2，满秩，系统可控

2. 修正扩展系统的代价函数

对于仿射系统的扩展形式，代价函数应该只惩罚原始状态x_t，不要惩罚新增的常数状态。把Q矩阵修改为：

Q = np.eye(x_dim + 1)
Q[-1, -1] = 0  # 关闭对新增常数状态的惩罚

3. 更合理的思路：直接处理仿射系统的LQR问题

不用扩展成线性系统，而是通过稳态偏差转化来求解：

• 第一步：求系统稳态x_ss和u_ss，满足x_ss = A x_ss + B u_ss + b，解出u_ss（假设B列满秩）：

  # 解稳态方程：(I - A)x_ss - B u_ss = b
  M = np.hstack([(np.eye(x_dim) - A), -B])
  sol = np.linalg.lstsq(M, b, rcond=None)[0]
  x_ss = sol[:x_dim]
  u_ss = sol[x_dim:]
  

• 第二步：定义偏差tilde_x = x_t - x_ss，tilde_u = u_t - u_ss，得到偏差线性系统tilde_x_{t+1} = A tilde_x + B tilde_u
• 第三步：对这个偏差系统求解Riccati方程，得到反馈增益K，最终控制律为u_t = u_ss - K @ (x_t - x_ss)

修正后的完整示例代码 #

import numpy as np
from scipy.linalg import solve_discrete_are, solve

x_dim = 2
u_dim = 1

# 修改B让系统可镇定
A = np.ones((x_dim, x_dim))
B = np.array([[1], [0]])  # 替换原B为满秩可控的矩阵
b = np.ones((x_dim,))

# 构造扩展系统
A_new = np.zeros((x_dim + 1, x_dim + 1))
A_new[:x_dim, :x_dim] = A
A_new[:x_dim, -1] = b
A_new[-1, -1] = 1

B_new = np.zeros((x_dim + 1, u_dim))
B_new[:x_dim, :] = B
B_new[-1, :] = 0

# 修正Q矩阵，不惩罚新增的常数状态
Q = np.eye(x_dim + 1)
Q[-1, -1] = 0
R = np.eye(u_dim)
S = np.zeros((x_dim + 1, u_dim))

# 求解Riccati方程和反馈增益
try:
    X = solve_discrete_are(A_new, B_new, Q, R, e=None, s=S)
    K = solve(B_new.T @ X @ B_new + R, B_new.T @ X @ A_new + S.T)
    print("求解得到的X矩阵：")
    print(X)
    print("\n反馈增益K：")
    print(K)
except np.linalg.LinAlgError as e:
    print(f"求解出错：{e}")

备注：内容来源于stack exchange，提问作者Jabby