8000 Fix bug in dare so that it returns a stabilizing solution · python-control/python-control@7cf9630 · GitHub
[go: up one dir, main page]
Skip to content

Commit 7cf9630

Browse files
cwrowleycwrowley
committed
Fix bug in dare so that it returns a stabilizing solution
The new implementation calls the routine scipy.linalg.solve_discrete_are(A, B, Q, R) which is included in scipy versions >= 0.11. The old implementation using slycot did satisfy the Riccati equation, but did not return a stabilizing solution. The scipy implementation apparently works correctly, though. This change fixes #8. Unit tests now make sure closed-loop eigenvalues lie inside the unit circle. Note: the scipy implementation handles only the case S = 0, E = I, the default values. If S and E are specified, the old routine (using slycot) is called. This passes the existing tests, but the tests include only one simple case, so it would be good to test this more extensively.
1 parent 8c1fdeb commit 7cf9630

File tree

2 files changed

+30
-6
lines changed

2 files changed

+30
-6
lines changed

‎control/mateqn.py‎

Lines changed: 15 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,8 @@
4242
"""
4343

4444
from numpy.linalg import inv
45-
from scipy import shape, size, asarray, copy, zeros, eye, dot
45+
from scipy import shape, size, asarray, asmatrix, copy, zeros, eye, dot
46+
from scipy.linalg import eigvals, solve_discrete_are
4647
from .exception import ControlSlycot, ControlArgument
4748

4849
#### Lyapunov equation solvers lyap and dlyap
@@ -667,7 +668,6 @@ def care(A,B,Q,R=None,S=None,E=None):
667668
else:
668669
raise ControlArgument("Invalid set of input parameters.")
669670

670-
671671
def dare(A,B,Q,R,S=None,E=None):
672672
""" (X,L,G) = dare(A,B,Q,R) solves the discrete-time algebraic Riccati
673673
equation
@@ -688,8 +688,19 @@ def dare(A,B,Q,R,S=None,E=None):
688688
where A, Q and E are square matrices of the same dimension. Further, Q and
689689
R are symmetric matrices. The function returns the solution X, the gain
690690
matrix G = (B^T X B + R)^-1 (B^T X A + S^T) and the closed loop
691-
eigenvalues L, i.e., the eigenvalues of A - B G , E. """
692-
691+
eigenvalues L, i.e., the eigenvalues of A - B G , E.
692+
"""
693+
if S is not None or E is not None:
694+
return dare_old(A, B, Q, R, S, E)
695+
else:
696+
Rmat = asmatrix(R)
697+
Qmat = asmatrix(Q)
698+
X = solve_discrete_are(A, B, Qmat, Rmat)
699+
G = inv(B.T.dot(X).dot(B) + Rmat) * B.T.dot(X).dot(A)
700+
L = eigvals(A - B.dot(G))
701+
return X, L, G
702+
703+
def dare_old(A,B,Q,R,S=None,E=None):
693704
# Make sure we can import required slycot routine
694705
try:
695706
from slycot import sb02md

‎control/tests/mateqn_test.py‎

Lines changed: 15 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -44,8 +44,9 @@
4444

4545
import unittest
4646
from numpy import matrix
47-
from numpy.testing import assert_array_almost_equal
48-
from numpy.linalg import inv
47+
from numpy.testing import assert_array_almost_equal, assert_array_less
48+
# need scipy version of eigvals for generalized eigenvalue problem
49+
from scipy.linalg import inv, eigvals
4950
from scipy import zeros,dot
5051
from control.mateqn import lyap,dlyap,care,dare
5152
from control.exception import slycot_check
@@ -178,6 +179,9 @@ def test_dare(self):
178179
A.T * X * A - X -
179180
A.T * X * B * inv(B.T * X * B + R) * B.T * X * A + Q, zeros((2,2)))
180181
assert_array_almost_equal(inv(B.T * X * B + R) * B.T * X * A, G)
182+
# check for stable closed loop
183+
lam = eigvals(A - B * G)
184+
assert_array_less(abs(lam), 1.0)
181185

182186
A = matrix([[1, 0],[-1, 1]])
183187
Q = matrix([[0, 1],[1, 1]])
@@ -190,6 +194,9 @@ def test_dare(self):
190194
A.T * X * A - X -
191195
A.T * X * B * inv(B.T * X * B + R) * B.T * X * A + Q, zeros((2,2)))
192196
assert_array_almost_equal(B.T * X * A / (B.T * X * B + R), G)
197+
# check for stable closed loop
198+
lam = eigvals(A - B * G)
199+
assert_array_less(abs(lam), 1.0)
193200

194201
def test_dare_g(self):
195202
A = matrix([[-0.6, 0],[-0.1, -0.4]])
@@ -206,6 +213,9 @@ def test_dare_g(self):
206213
(A.T * X * B + S) * inv(B.T * X * B + R) * (B.T * X * A + S.T) + Q,
207214
zeros((2,2)) )
208215
assert_array_almost_equal(inv(B.T * X * B + R) * (B.T * X * A + S.T), G)
216+
# check for stable closed loop
217+
lam = eigvals(A - B * G, E)
218+
assert_array_less(abs(lam), 1.0)
209219

210220
A = matrix([[-0.6, 0],[-0.1, -0.4]])
211221
Q = matrix([[2, 1],[1, 3]])
@@ -221,6 +231,9 @@ def test_dare_g(self):
221231
(A.T * X * B + S) * inv(B.T * X * B + R) * (B.T * X * A + S.T) + Q,
222232
zeros((2,2)) )
223233
assert_array_almost_equal((B.T * X * A + S.T) / (B.T * X * B + R), G)
234+
# check for stable closed loop
235+
lam = eigvals(A - B * G, E)
236+
assert_array_less(abs(lam), 1.0)
224237

225238
def suite():
226239
return unittest.TestLoader().loadTestsFromTestCase(TestMatrixEquations)

0 commit comments

Comments
 (0)
0