SoS求解器数值异常求助：验证系统RoA时约束违反而返回成功

问题核心 #

验证RoA时，Mosek/CSDP求解器返回成功，但采样点检查显示SOS约束违反，根源是问题尺度失衡（P矩阵元素量级差达3个数量级，rho与V(x)的尺度不匹配）导致的数值病态，求解器精度不足无法识别约束违反。

具体解决建议 #

• 全链路尺度归一化：
  仅缩放P是不够的，需同步调整所有相关变量：

  • 将P缩放1e5后，V(x) = xᵀPx的量级从1e2变为1e7，对应rho需从3e5调整为3e10；
  • eps按比例调整为1e-7 * 1e5 = 1e-2，保证eps*x.dot(x)的量级与其他项匹配；
  • 若nonlinear_system(x)的输出有量级差异，需同步缩放系统动态，确保整个优化问题的所有项量级在同一区间（如1e0~1e3）。

• 严格求解器精度配置：
  调整求解器的可行性检查阈值，以Mosek为例，添加以下配置：

  solver = MosekSolver()
  # 收紧相对间隙和不可行性阈值
  solver.SetSolverOption("MSK_DPAR_INTPNT_TOL_REL_GAP", 1e-10)
  solver.SetSolverOption("MSK_DPAR_INTPNT_TOL_INFEAS", 1e-10)
  

  CSDP可通过调整maxiter或eps参数提升精度。

• 重构SOS约束形式：
  将原约束-V_dot + l*(V - rho) - eps*x.dot(x) ∈ SOS改写为：

  # l为SOS多项式，约束改写为-V_dot - eps*x.dot(x) ≤ l*(rho - V)
  prog.AddSosConstraint(-V_dot - eps*x.dot(x) + l*(rho - V))
  

  避免V - rho为负时的负项叠加，减少数值抵消误差。

• 优化多项式l的配置：

  • 先尝试降低l的度数至2次，减少高次项带来的数值不稳定；
  • 添加l的下界约束，防止系数过小：

    # 强制l在原点处的取值不小于1e-6
    prog.AddConstraint(l.Evaluate({x[0]:0, x[1]:0, x[2]:0}) >= 1e-6)
    

• 采样点合理性验证：
  当前采样点x的绝对值过大（-1000、2000），可能超出Lyapunov函数的有效局部范围。先在V(x) ≤ 1e5的子水平集内采样（如x_sample = [-100, 0, 200]），确认约束是否真的违反，再逐步扩大范围。

附用户提供的原始信息 #

Lyapunov矩阵P #

P = [[ 6.88693913e-05  8.81733881e-09 -1.12596351e-06]
     [ 8.81733881e-09  6.89157664e-05  7.68071960e-05]
     [-1.12596351e-06  7.68071960e-05  5.00915262e-02]]

采样点检查输出 #

x_sample: 
[-1000     0  2000]
V(x_sample): 
200375.96901588634
rho:
300000.0
V - rho: 
-99624.03098411366
l(x_sample):
93125.29963587047
l(x_sample)*(V - rho): 
-9277517736.328829
-V_dot_x: 
-1176876.6274507358
-eps * x.T * x:
-0.5
-V_dot_x + l * (V- rho) - eps ( x.T * x)
-9278694612.95678

验证代码 #

def is_verified(rho):
    prog = MathematicalProgram()
    x = prog.NewIndeterminates(3, "x")
    V = x.dot(P_lyap).dot(x)
    V_dot = 2 * x.dot(P_lyap).dot(nonlinear_system(x))
    l_deg = 4
    assert l_deg % 2 == 0
    l = prog.NewSosPolynomial(Variables(x), l_deg)[0].ToExpression()
    eps = 1e-7
    prog.AddSosConstraint(-V_dot + l * ( V - rho ) - eps * x.dot(x))
    solver = MosekSolver()
    result = solver.Solve(prog)
    x_sample = np.array([-1000, 0, 2000])
    thisdict = {
        x[0]: x_sample[0],
        x[1]: x_sample[1],
        x[2]: x_sample[2]
    }
    l_result = result.GetSolution(l)
    l_sample =  l_result.Evaluate(thisdict)
    V_x_sample = x_sample.dot(P_lyap).dot(x_sample)
    V_dot_sample = 2 * x_sample.dot(P_lyap).dot(nonlinear_system(x_sample))
    print("x_sample: ")
    print(x_sample)
    print("V_x_sample: ")
    print(V_x_sample)
    print("rho")
    print(rho)
    print("V - rho: ") 
    print(V_x_sample - rho)
    print("l:")
    print(l_sample)
    print("l*(V - rho): ")
    print(l_sample*( V_x_sample - rho))
    print("-V_dot_x: ")
    print(-V_dot_sample)
    print("- eps * x.T * x")
    print(- eps * x_sample.dot(x_sample))
    print("-V_dot_x + l * (V- rho) - eps ( x.T * x)")
    print(-V_dot_sample + l_sample*( V_x_sample - rho) - eps * x_sample.dot(x_sample))
    return result.is_success()

内容的提问来源于stack exchange，提问作者Henry Andersen