Skip to content

Commit

Permalink
Merge pull request #437 from hinayuki64/tanahashi2
Browse files Browse the repository at this point in the history 
 Add method option in robustlq.py
  • Loading branch information
@mmcky
mmcky committed Oct 3, 2018
2 parents ba191fb + ca7feea commit f303822
Show file tree
Hide file tree
Showing 2 changed files with 49 additions and 43 deletions.
26 changes: 19 additions & 7 deletions quantecon/robustlq.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,7 +152,7 @@ def b_operator(self, P):

return F, new_P

def robust_rule(self):
def robust_rule(self, method='doubling'):
"""
This method solves the robust control problem by tricking it
into a stacked LQ problem, as described in chapter 2 of Hansen-
Expand All @@ -165,6 +165,12 @@ def robust_rule(self):
And the value function is :math:`-x'Px`
Parameters
----------
method : str, optional(default='doubling')
Solution method used in solving the associated Riccati
equation, str in {'doubling', 'qz'}.
Returns
-------
F : array_like(float, ndim=2)
Expand All @@ -189,7 +195,7 @@ def robust_rule(self):
lq = LQ(-beta*I*theta, R, A, C, beta=beta)

# == Solve and convert back to robust problem == #
P, f, d = lq.stationary_values()
P, f, d = lq.stationary_values(method=method)
F = np.zeros((self.k, self.n))
K = -f[:k, :]

Expand All @@ -199,7 +205,7 @@ def robust_rule(self):
lq = LQ(Qa, R, A, Ba, beta=beta)

# == Solve and convert back to robust problem == #
P, f, d = lq.stationary_values()
P, f, d = lq.stationary_values(method=method)
F = f[:k, :]
K = -f[k:f.shape[0], :]

Expand Down Expand Up @@ -254,14 +260,17 @@ def robust_rule_simple(self, P_init=None, max_iter=80, tol=1e-8):

return F, K, P

def F_to_K(self, F):
def F_to_K(self, F, method='doubling'):
"""
Compute agent 2's best cost-minimizing response K, given F.
Parameters
----------
F : array_like(float, ndim=2)
A k x n array
method : str, optional(default='doubling')
Solution method used in solving the associated Riccati
equation, str in {'doubling', 'qz'}.
Returns
-------
Expand All @@ -276,18 +285,21 @@ def F_to_K(self, F):
A2 = self.A - dot(self.B, F)
B2 = self.C
lq = LQ(Q2, R2, A2, B2, beta=self.beta)
neg_P, neg_K, d = lq.stationary_values()
neg_P, neg_K, d = lq.stationary_values(method=method)

return -neg_K, -neg_P

def K_to_F(self, K):
def K_to_F(self, K, method='doubling'):
"""
Compute agent 1's best value-maximizing response F, given K.
Parameters
----------
K : array_like(float, ndim=2)
A j x n array
method : str, optional(default='doubling')
Solution method used in solving the associated Riccati
equation, str in {'doubling', 'qz'}.
Returns
-------
Expand All @@ -302,7 +314,7 @@ def K_to_F(self, K):
Q1 = self.Q
R1 = self.R - self.beta * self.theta * dot(K.T, K)
lq = LQ(Q1, R1, A1, B1, beta=self.beta)
P, F, d = lq.stationary_values()
P, F, d = lq.stationary_values(method=method)

return F, P

Expand Down
66 changes: 30 additions & 36 deletions quantecon/tests/test_robustlq.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,9 +51,7 @@ def setUp(self):
self.rblq_test = RBLQ(Q, R, A, B, C, beta, theta)
self.rblq_test_pf = RBLQ(Q_pf, R, A, B_pf, C, beta, theta)
self.lq_test = LQ(Q, R, A, B, C, beta)

self.Fr, self.Kr, self.Pr = self.rblq_test.robust_rule()
self.Fr_pf, self.Kr_pf, self.Pr_pf = self.rblq_test_pf.robust_rule()
self.methods = ['doubling', 'qz']

def tearDown(self):
del self.rblq_test
Expand All @@ -65,52 +63,48 @@ def test_pure_forecasting(self):
def test_robust_rule_vs_simple(self):
rblq = self.rblq_test
rblq_pf = self.rblq_test_pf
Fr, Kr, Pr = self.Fr, self.Kr, self.Pr
Fr_pf, Kr_pf, Pr_pf = self.Fr_pf, self.Kr_pf, self.Pr_pf

Fs, Ks, Ps = rblq.robust_rule_simple(P_init=Pr, tol=1e-12)
Fs_pf, Ks_pf, Ps_pf = rblq_pf.robust_rule_simple(P_init=Pr_pf, tol=1e-12)
for method in self.methods:
Fr, Kr, Pr = self.rblq_test.robust_rule(method=method)
Fr_pf, Kr_pf, Pr_pf = self.rblq_test_pf.robust_rule(method=method)

Fs, Ks, Ps = rblq.robust_rule_simple(P_init=Pr, tol=1e-12)
Fs_pf, Ks_pf, Ps_pf = rblq_pf.robust_rule_simple(P_init=Pr_pf, tol=1e-12)

assert_allclose(Fr, Fs, rtol=1e-4)
assert_allclose(Kr, Ks, rtol=1e-4)
assert_allclose(Pr, Ps, rtol=1e-4)
assert_allclose(Fr, Fs, rtol=1e-4)
assert_allclose(Kr, Ks, rtol=1e-4)
assert_allclose(Pr, Ps, rtol=1e-4)

atol = 1e-10
assert_allclose(Fr_pf, Fs_pf, rtol=1e-4)
assert_allclose(Kr_pf, Ks_pf, rtol=1e-4, atol=atol)
assert_allclose(Pr_pf, Ps_pf, rtol=1e-4, atol=atol)
atol = 1e-10
assert_allclose(Fr_pf, Fs_pf, rtol=1e-4)
assert_allclose(Kr_pf, Ks_pf, rtol=1e-4, atol=atol)
assert_allclose(Pr_pf, Ps_pf, rtol=1e-4, atol=atol)


def test_f2k_and_k2f(self):
rblq = self.rblq_test
Fr, Kr, Pr = self.Fr, self.Kr, self.Pr

K_f2k, P_f2k = rblq.F_to_K(Fr)
F_k2f, P_k2f = rblq.K_to_F(Kr)

assert_allclose(K_f2k, Kr, rtol=1e-4)
assert_allclose(F_k2f, Fr, rtol=1e-4)
assert_allclose(P_f2k, P_k2f, rtol=1e-4)
for method in self.methods:
Fr, Kr, Pr = self.rblq_test.robust_rule(method=method)
K_f2k, P_f2k = rblq.F_to_K(Fr, method=method)
F_k2f, P_k2f = rblq.K_to_F(Kr, method=method)
assert_allclose(K_f2k, Kr, rtol=1e-4)
assert_allclose(F_k2f, Fr, rtol=1e-4)
assert_allclose(P_f2k, P_k2f, rtol=1e-4)

def test_evaluate_F(self):
rblq = self.rblq_test
Fr, Kr, Pr = self.Fr, self.Kr, self.Pr

Kf, Pf, df, Of, of = rblq.evaluate_F(Fr)

# In the future if we wanted, we could check more things, but I
# think the other pieces are basically just plugging these into
# equations so if these hold then the others should be correct
# as well.
assert_allclose(Pf, Pr)
assert_allclose(Kf, Kr)





for method in self.methods:
Fr, Kr, Pr = self.rblq_test.robust_rule(method=method)

Kf, Pf, df, Of, of = rblq.evaluate_F(Fr)

# In the future if we wanted, we could check more things, but I
# think the other pieces are basically just plugging these into
# equations so if these hold then the others should be correct
# as well.
assert_allclose(Pf, Pr)
assert_allclose(Kf, Kr)

if __name__ == '__main__':
suite = unittest.TestLoader().loadTestsFromTestCase(TestRBLQControl)
Expand Down

0 comments on commit f303822

Please sign in to comment.