# [Scipy-tickets] [SciPy] #1732: solve_discrete_are in scipy.linalg does (sometimes) not solve correctly

SciPy Trac scipy-tickets@scipy....
Wed Sep 19 23:24:06 CDT 2012

```#1732: solve_discrete_are in scipy.linalg does (sometimes) not solve correctly
-----------------------------------------+----------------------------------
Reporter:  inde                         |       Owner:  somebody
Type:  defect                       |      Status:  new
Priority:  normal                       |   Milestone:  Unscheduled
Component:  scipy.linalg                 |     Version:  devel
Keywords:  riccati, solve_discrete_are  |
-----------------------------------------+----------------------------------
Hi,
this is for scipy 0.11.0rc2. It seemed to be a normal release on
sourceforge, but there is no 0.11 here, so I am filing under devel.

So the problem I am observing when solving the discrete algebraic Riccati
Equation for the linear quadratic regulator problem is that sometimes the
solution returned by scipy.linalg.solve_discrete_are is not positive
definite or even correct.

I am testing with 3x3 random matrices (code below). The attached npz file
contains one particular case that did not solve for me.

def dare_zero_compare(self,A,B,Q,R,S):
AT_S_A = s.dot(A.T,s.dot(S,A))
AT_S_B = s.dot(A.T,s.dot(S,B))
BT_S_B_plusR_inv = s.linalg.inv(s.dot(B.T,s.dot(S,B)) + R)
BT_S_A = s.dot(B.T,s.dot(S,A))

zeroS = AT_S_A - S - s.dot(AT_S_B, s.dot(BT_S_B_plusR_inv,
BT_S_A)) + Q
return zeroS

def test_dareSsolves(self):
for i in range(100):
n = 3
A = s.random.rand(n,n)
B = s.random.rand(n,1)
R = s.eye(B.shape[1])
Q = s.eye(A.shape[0])

S = s.linalg.solve_discrete_are(A,B,Q,R)

zeroS = self.dare_zero_compare(A,B,Q,R,S)
try:
for ai in zeroS:
for aij in ai:
self.assertAlmostEqual(aij, 0.0)
except AssertionError:
with open ('test_dareSsolves.m','w') as f:
f.write(string_dare_test(A,B,Q,R,S) )
s.savez('test_dareSsolves.npz', A=A,B=B,Q=Q,R=R,S=S)

def test_dareSsolves_repeat(self):
A = npzfile['A']
B = npzfile['B']
Q = npzfile['Q']
R = npzfile['R']
S = npzfile['S']

S_new = s.linalg.solve_discrete_are(A,B,Q,R)
print A,B, Q,R,S,S_new

zeroS = self.dare_zero_compare(A,B,Q,R,S)
zeroS_new = self.dare_zero_compare(A,B,Q,R,S_new)

for ai in zeroS:
for aij in ai:
self.assertAlmostEqual(aij, 0.0)
raise AssertionError

--
Ticket URL: <http://projects.scipy.org/scipy/ticket/1732>
SciPy <http://www.scipy.org>
SciPy is open-source software for mathematics, science, and engineering.
```