Skip to content

BUG: solve_discrete_riccati returns non-stabilizing solution #356

Closed
@oyamad

Description

@oyamad

I have got an LQ problem for which LQ.stationary_values returns some suspicious solution, where solve_discrete_riccati is a suspect.

Different solutions: qe.solve_discrete_riccati versus scipy.linalg.solve_discrete_are

Consider the following inputs:

A = np.array([[1.0, 0.0], [0.225, 0.0]])
B = np.array([[0.0], [1.11111]])
R = np.array([[2.89634, -4.0404], [-4.0404, 8.16243]])
Q = np.array([[8.16243]])
N = np.array([4.0404, -8.16243])
beta = 0.9
A0, B0 = np.sqrt(beta) * A, np.sqrt(beta) * B

qe.solve_discrete_riccati returns a solution different from one from scipy.linalg.solve_discrete_are.

Solution by qe.solve_discrete_riccati:

P_qe = qe.solve_discrete_riccati(A0, B0, R, Q, N)
P_qe
array([[ 8.96343411,  0.        ],
       [ 0.        ,  0.        ]])

Solution by scipy.linalg.solve_discrete_are:

P_sp = scipy.linalg.solve_discrete_are(A0, B0, R, Q, s=N.reshape(2, 1))
P_sp
array([[ 15.94687678,  -2.38748478],
       [ -2.38748478,   0.81622831]])

Both are indeed solutions:

def riccati_rhs(X, A, B, Q, R, N):
    out = A.T @ X @ A - \
        (N + B.T @ X @ A).T @ np.linalg.inv(B.T @ X @ B + R) @ \
        (N + B.T @ X @ A) + Q
    return out
riccati_rhs(P_qe, A0, B0, R, Q, N)

array([[ 8.96343411,  0.        ],
       [ 0.        ,  0.        ]])
riccati_rhs(P_sp, A0, B0, R, Q, N)

array([[ 15.94687678,  -2.38748478],
       [ -2.38748478,   0.81622831]])

Stability of A - B @ F

Consider the LQ control problem defined by the inputs above, and suppose that F is the matrix that gives the optimal control. Then G = A - B @ F is the matrix that generates state transition:

def compute_F(P, A, B, Q, N, beta):
    S1 = Q + beta * np.dot(B.T, np.dot(P, B))
    S2 = beta * np.dot(B.T, np.dot(P, A)) + N
    F = np.linalg.solve(S1, S2)
    return F
F_qe = compute_F(P_qe, A, B, Q, N, beta)
F_qe

array([[ 0.49499965, -1.        ]])
F_sp = compute_F(P_sp, A, B, Q, N, beta)
F_sp

array([[ 0.20250284, -0.9000018 ]])

Eigenvalues of G_qe:

G_qe = A - B @ F_qe
w, v = np.linalg.eig(G_qe)
w

array([ 1.11111,  1.     ])

"Discounted eigenvalues":

w * beta

array([ 0.999999,  0.9     ])

Eigenvalues of G_sp:

G_sp = A - B @ F_sp
w, v = np.linalg.eig(G_sp)
w

array([ 1.000001,  1.      ])

Original approximating LQ control problem

The inputs came from an LQ approximation of a model growth model. See this notebook for details.

Metadata

Metadata

Assignees

No one assigned

    Labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions