10
10
2. Fix warnings from checking conditioning of matrices
11
11
"""
12
12
import numpy as np
13
- from numpy import dot
14
13
from numpy .linalg import solve
15
14
from scipy .linalg import solve_discrete_lyapunov as sp_solve_discrete_lyapunov
16
15
from scipy .linalg import solve_discrete_are as sp_solve_discrete_are
@@ -172,19 +171,19 @@ def solve_discrete_riccati(A, B, Q, R, N=None, tolerance=1e-10, max_iter=500,
172
171
# == Choose optimal value of gamma in R_hat = R + gamma B'B == #
173
172
current_min = np .inf
174
173
candidates = (0.01 , 0.1 , 0.25 , 0.5 , 1.0 , 2.0 , 10.0 , 100.0 , 10e5 )
175
- BB = dot (B .T , B )
176
- BTA = dot (B .T , A )
174
+ BB = np . dot (B .T , B )
175
+ BTA = np . dot (B .T , A )
177
176
for gamma in candidates :
178
177
Z = R + gamma * BB
179
178
cn = np .linalg .cond (Z )
180
179
if cn * EPS < 1 :
181
- Q_tilde = - Q + dot (N .T , solve (Z , N + gamma * BTA )) + gamma * I
182
- G0 = dot (B , solve (Z , B .T ))
183
- A0 = dot (I - gamma * G0 , A ) - dot (B , solve (Z , N ))
184
- H0 = gamma * dot (A .T , A0 ) - Q_tilde
180
+ Q_tilde = - Q + np . dot (N .T , solve (Z , N + gamma * BTA )) + gamma * I
181
+ G0 = np . dot (B , solve (Z , B .T ))
182
+ A0 = np . dot (I - gamma * G0 , A ) - np . dot (B , solve (Z , N ))
183
+ H0 = gamma * np . dot (A .T , A0 ) - Q_tilde
185
184
f1 = np .linalg .cond (Z , np .inf )
186
185
f2 = gamma * f1
187
- f3 = np .linalg .cond (I + dot (G0 , H0 ))
186
+ f3 = np .linalg .cond (I + np . dot (G0 , H0 ))
188
187
f_gamma = max (f1 , f2 , f3 )
189
188
if f_gamma < current_min :
190
189
best_gamma = gamma
@@ -199,10 +198,10 @@ def solve_discrete_riccati(A, B, Q, R, N=None, tolerance=1e-10, max_iter=500,
199
198
R_hat = R + gamma * BB
200
199
201
200
# == Initial conditions == #
202
- Q_tilde = - Q + dot (N .T , solve (R_hat , N + gamma * BTA )) + gamma * I
203
- G0 = dot (B , solve (R_hat , B .T ))
204
- A0 = dot (I - gamma * G0 , A ) - dot (B , solve (R_hat , N ))
205
- H0 = gamma * dot (A .T , A0 ) - Q_tilde
201
+ Q_tilde = - Q + np . dot (N .T , solve (R_hat , N + gamma * BTA )) + gamma * I
202
+ G0 = np . dot (B , solve (R_hat , B .T ))
203
+ A0 = np . dot (I - gamma * G0 , A ) - np . dot (B , solve (R_hat , N ))
204
+ H0 = gamma * np . dot (A .T , A0 ) - Q_tilde
206
205
i = 1
207
206
208
207
# == Main loop == #
@@ -212,9 +211,9 @@ def solve_discrete_riccati(A, B, Q, R, N=None, tolerance=1e-10, max_iter=500,
212
211
raise ValueError (fail_msg .format (i ))
213
212
214
213
else :
215
- A1 = dot (A0 , solve (I + dot (G0 , H0 ), A0 ))
216
- G1 = G0 + dot (dot (A0 , G0 ), solve (I + dot (H0 , G0 ), A0 .T ))
217
- H1 = H0 + dot (A0 .T , solve (I + dot (H0 , G0 ), dot (H0 , A0 )))
214
+ A1 = np . dot (A0 , solve (I + np . dot (G0 , H0 ), A0 ))
215
+ G1 = G0 + np . dot (np . dot (A0 , G0 ), solve (I + np . dot (H0 , G0 ), A0 .T ))
216
+ H1 = H0 + np . dot (A0 .T , solve (I + np . dot (H0 , G0 ), np . dot (H0 , A0 )))
218
217
219
218
error = np .max (np .abs (H1 - H0 ))
220
219
A0 = A1
0 commit comments