Skip to content
New issue

Have a question about this project? # for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “#”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? # to your account

BUG: solve_discrete_riccati returns non-stabilizing solution #356

Closed
oyamad opened this issue Oct 22, 2017 · 5 comments
Closed

BUG: solve_discrete_riccati returns non-stabilizing solution #356

oyamad opened this issue Oct 22, 2017 · 5 comments
Labels

Comments

@oyamad
Copy link
Member

oyamad commented Oct 22, 2017

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.

@oyamad
Copy link
Member Author

oyamad commented Oct 23, 2017

The same result obtains with QuantEcon.jl, by the way.

@oyamad oyamad changed the title Unstable? solution by solve_discrete_riccati BUG: solve_discrete_riccati returns non-stabilizing solution Oct 25, 2017
@oyamad
Copy link
Member Author

oyamad commented Oct 25, 2017

It turned out that qe.solve_discrete_riccati returns a non-stabilizing solution for the Riccati equation given by the inputs A0, B0, R, Q, N above, as opposed to the claim in the reference Chiang, Fan, and Lin.

A solution P of the Riccati equation with A, B, Q, R, N is stabilizing if S = R + B.T @ P @ B is invertible and the absolute values of all the eigenvalues of A_F = A - B @ F are smaller than one, where F is the matrix defined in the LQ control lecture (so is different from the one in Chiang, Fan, and Lin in its sign).

def compute_F_S_A_F(P, A, B, R, N):
    S1 = R + B.T @ P @ B
    S2 = B.T @ P @ A + N
    F = np.linalg.solve(S1, S2)
    A_F = A - B @ F
    return F, S1, A_F

For (A, B, Q, R, N) = (A0, B0, R, Q, N) above (note the reversion of Q and R):

P_qe = qe.solve_discrete_riccati(A0, B0, R, Q, N)
F_qe, S_qe, A_F_qe = compute_F_S_A_F(P_qe, A0, B0, Q, N)
F_qe, S_qe, A_F_qe
(array([[ 0.49499965, -1.        ]]),
 array([[ 8.16243]]),
 array([[ 0.9486833 ,  0.        ],
        [-0.30832118,  1.0540915 ]]))
w, v = np.linalg.eig(A_F_qe)
w
array([ 1.0540915,  0.9486833])

For the solution by scipy.linalg.solve_discrete_are:

P_sp = scipy.linalg.solve_discrete_are(A0, B0, R, Q, s=N.reshape(2, 1))
F_sp, S_sp, A_F_sp = compute_F_S_A_F(P_sp, A0, B0, Q, N)
F_sp, S_sp, A_F_sp
(array([[ 0.20250284, -0.9000018 ]]),
 array([[ 9.06934853]]),
 array([[  9.48683298e-01,   0.00000000e+00],
        [ -2.77492114e-06,   9.48684247e-01]]))
w, v = np.linalg.eig(A_F_sp)
w
array([ 0.94868425,  0.9486833 ])

@oyamad oyamad added the bug label Oct 25, 2017
@oyamad
Copy link
Member Author

oyamad commented Oct 25, 2017

@jstac How did you choose the values of gamma in candidates (in the initialization step in solve_discrete_riccati)?

@oyamad
Copy link
Member Author

oyamad commented Oct 25, 2017

Removing gamma=0.0 from candidates, I got the same output as scipy.linalg.solve_discrete_are (i.e., the stabilizing solution).

The paper Chiang, Fan, and Lin says gamma > 0. Can we just remove 0.0 from candidates?

@jstac
Copy link
Contributor

jstac commented Oct 25, 2017

Well done @oyamad , thanks for spotting this. I can't remember how I chose candidates but it seems that 0.0 should not be in there. Let's remove it.

mmcky pushed a commit that referenced this issue Oct 25, 2017
* FIX: Remove 0.0 from `candidates` in `solve_discrete_riccati`

Fix #356

* TEST: Set atol in `test_robust_rule_vs_simple`
# for free to join this conversation on GitHub. Already have an account? # to comment
Labels
Projects
None yet
Development

No branches or pull requests

2 participants