[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  |  
 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)
                 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):
         npzfile = s.load('test_dareSsolves.npz')
         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.

More information about the Scipy-tickets mailing list