@@ -9,6 +9,20 @@ kernelspec:
99 display_name : Python 3
1010 language : python
1111 name : python3
12+ translation :
13+ title : LQ控制的拉格朗日方法
14+ headings :
15+ Overview : 概述
16+ Undiscounted LQ DP Problem : 无折现LQ动态规划问题
17+ Lagrangian : 拉格朗日量
18+ State-Costate Dynamics : 状态-协状态动态
19+ Reciprocal Pairs Property : 倒数对性质
20+ Schur decomposition : Schur分解
21+ Application : 应用
22+ Other Applications : 其他应用
23+ Discounted Problems : 折现问题
24+ Discounted Problems::Transforming States and Controls to Eliminate Discounting : 转换状态和控制以消除折现
25+ Discounted Problems::Lagrangian for Discounted Problem : 折现问题的拉格朗日量
1226---
1327
1428+++
@@ -186,28 +200,28 @@ $$ (eqn:muPx)
186200* 将得到的方程和 {eq}`lag-lqdp-eq2` 的第二个方程整理成如下形式
187201
188202$$
189- L\ \begin{pmatrix }x_ {t+1}\cr \mu_ {t+1}\cr\end{pmatrix }\ = \ N\ \begin{pmatrix }x_t\cr \mu_t\cr\end{pmatrix }\
203+ L\ \begin{bmatrix }x_ {t+1}\cr \mu_ {t+1}\cr\end{bmatrix }\ = \ N\ \begin{bmatrix }x_t\cr \mu_t\cr\end{bmatrix }\
190204,\ t \geq 0,
191205$$ (eq:systosolve)
192206
193207其中
194208
195209$$
196- L = \ \begin{pmatrix }I & BQ^{-1} B^\prime \cr 0 & A^\prime\cr\end{pmatrix }, \quad N = \
197- \begin{pmatrix }A & 0\cr -R & I\cr\end{pmatrix }.
210+ L = \ \begin{bmatrix }I & BQ^{-1} B^\prime \cr 0 & A^\prime\cr\end{bmatrix }, \quad N = \
211+ \begin{bmatrix }A & 0\cr -R & I\cr\end{bmatrix }.
198212$$
199213
200214当 $L$ 满秩时(即当 $A$ 满秩时),我们可以将系统 {eq}`eq:systosolve` 写作
201215
202216$$
203- \begin{pmatrix }x_ {t+1}\cr \mu_ {t+1}\cr\end{pmatrix }\ = M\ \begin{pmatrix }x_t\cr\mu_t\cr\end{pmatrix }
217+ \begin{bmatrix }x_ {t+1}\cr \mu_ {t+1}\cr\end{bmatrix }\ = M\ \begin{bmatrix }x_t\cr\mu_t\cr\end{bmatrix }
204218$$ (eq4orig)
205219
206220其中
207221
208222$$
209- M\equiv L^{-1} N = \begin{pmatrix }A+B Q^{-1} B^\prime A^{\prime-1}R &
210- -B Q^{-1} B^\prime A^{\prime-1}\cr -A^{\prime -1} R & A^{\prime -1}\cr\end{pmatrix }.
223+ M\equiv L^{-1} N = \begin{bmatrix }A+B Q^{-1} B^\prime A^{\prime-1}R &
224+ -B Q^{-1} B^\prime A^{\prime-1}\cr -A^{\prime -1} R & A^{\prime -1}\cr\end{bmatrix }.
211225$$ (Mdefn)
212226
213227+++
@@ -216,8 +230,8 @@ $$ (Mdefn)
216230
217231我们希望求解差分方程系统{eq}`eq4orig`,其解为满足以下条件的序列 $\{x_t\}_{t=0}^\infty$:
218232
219- * 初始条件为
220- * 终端条件为 $\lim_{t \rightarrow +\infty} x_t =0$
233+ * $x_0$ 的一个初始条件
234+ * 终端条件 $\lim_{t \rightarrow +\infty} x_t =0$
221235
222236这个终端条件反映了我们对**稳定**解的需求,即解不会在 $t \to \infty$ 时发散。
223237
238252为此,我们引入一个 $(2n \times 2n)$ 矩阵
239253
240254$$
241- J = \begin{pmatrix }0 & -I_n\cr I_n & 0\cr\end{pmatrix }.
255+ J = \begin{bmatrix }0 & -I_n\cr I_n & 0\cr\end{bmatrix }.
242256$$
243257
244258矩阵 $J$ 的秩为 $2n$。
280294y_ {t+1} = M y_t
281295$$ (eq658)
282296
283- 其中 $y_t = \begin{pmatrix }x_t\cr \mu_t\cr\end{pmatrix }$。
297+ 其中 $y_t = \begin{bmatrix }x_t\cr \mu_t\cr\end{bmatrix }$。
284298
285299考虑 $M$ 的**三角化**
286300
287301$$
288- V^{-1} M V= \begin{pmatrix }W_ {11} & W_ {12} \cr 0 & W_ {22}\cr\end{pmatrix }
302+ V^{-1} M V= \begin{bmatrix }W_ {11} & W_ {12} \cr 0 & W_ {22}\cr\end{bmatrix }
289303$$ (eqn:triangledecomp)
290304
291305其中
323337将方程{eq}`eq6510`写作
324338
325339$$
326- \begin{pmatrix }y^\ast_ {1t}\cr y^\ast_ {2t}\cr\end{pmatrix }\ =\ \left[ \begin{matrix} W^t_ {11} &
327- W_ {12, t}\cr 0 & W^t_ {22}\cr\end{matrix}\right] \quad \begin{pmatrix }y^\ast_ {10}\cr
328- y^\ast_ {20}\cr\end{pmatrix }
340+ \begin{bmatrix }y^\ast_ {1t}\cr y^\ast_ {2t}\cr\end{bmatrix }\ =\ \left[ \begin{matrix} W^t_ {11} &
341+ W_ {12, t}\cr 0 & W^t_ {22}\cr\end{matrix}\right] \quad \begin{bmatrix }y^\ast_ {10}\cr
342+ y^\ast_ {20}\cr\end{bmatrix }
329343$$
330344
331345其中 $y^\ast_t = V^{-1} y_t$,特别地
358372\mu_t = - (V^{22})^{-1} V^{21} x_t.
359373$$
360374
361- 但注意,由于 $(V^{21} \ V^{22})$ 是 $V$ 的逆矩阵的第二行块,因此
375+ 但注意,由于 $(V^{21}\ V^{22})$ 是 $V$ 的逆矩阵的第二行块,因此
362376
363377$$
364- (V^{21} \ V^{22})\quad \begin{pmatrix }V_ {11}\cr V_ {21}\cr\end{pmatrix } = 0
378+ (V^{21} \ V^{22})\quad \begin{bmatrix }V_ {11}\cr V_ {21}\cr\end{bmatrix } = 0
365379$$
366380
367381这意味着
@@ -406,11 +420,14 @@ $$ (eqn:Pvaughn)
406420
407421## 应用
408422
409- 这里我们通过一个示例来演示计算过程,这个示例是从[quantecon讲座](https://python.quantecon.org/lqcontrol.html)中借鉴的确定性版本。
423+ 这里我们用来自 {doc}`lqcontrol` 的确定性永久收入示例来演示这一计算方法。
424+
425+ 由于该模型是带折现的,我们将不变子空间方法应用于由变换矩阵 $\hat A = \beta^{1/2} A$ 和 $\hat B = \beta^{1/2} B$ 得到的等价*无折现*系统。
410426
411427```{code-cell} ipython3
412428# 模型参数
413429r = 0.05
430+ β = 1 / (1 + r)
414431c_bar = 2
415432μ = 1
416433
@@ -423,15 +440,15 @@ B = [[-1],
423440 [0]]
424441
425442# 构造一个LQ实例
426- lq = LQ(Q, R, A, B)
443+ lq = LQ(Q, R, A, B, beta=β )
427444```
428445
429446给定矩阵 $A$、$B$、$Q$、$R$,我们可以计算 $L$、$N$ 和 $M=L^{-1}N$。
430447
431448```{code-cell} ipython3
432449def construct_LNM(A, B, Q, R):
433450
434- n, k = lq.n, lq.k
451+ n = A.shape[0]
435452
436453 # 构造 L 和 N
437454 L = np.zeros((2*n, 2*n))
@@ -451,7 +468,10 @@ def construct_LNM(A, B, Q, R):
451468```
452469
453470```{code-cell} ipython3
454- L, N, M = construct_LNM(lq.A, lq.B, lq.Q, lq.R)
471+ A_bar = lq.A * lq.beta ** (1/2)
472+ B_bar = lq.B * lq.beta ** (1/2)
473+
474+ L, N, M = construct_LNM(A_bar, B_bar, lq.Q, lq.R)
455475```
456476
457477```{code-cell} ipython3
@@ -472,7 +492,7 @@ M @ J @ M.T - J
472492我们可以使用`np.linalg.eigvals`计算矩阵$M$的特征值,并按升序排列。
473493
474494```{code-cell} ipython3
475- eigvals = sorted(np.linalg.eigvals(M))
495+ eigvals = sorted(np.linalg.eigvals(M), key=lambda z: (abs(z), z.real, z.imag) )
476496eigvals
477497```
478498
@@ -484,18 +504,14 @@ eigvals
484504为了得到我们想要的结果,让我们定义一个排序函数,告诉`scipy.schur`将模小于1的对应特征值排序到左上角。
485505
486506```{code-cell} ipython3
487- stable_eigvals = eigvals[:n]
507+ tol = 1e-10
488508
489509def sort_fun(x):
490510 "将模小于1的特征值排序到左上角。"
511+ return abs(x) < 1 - tol
491512
492- if x in stable_eigvals:
493- stable_eigvals.pop(stable_eigvals.index(x))
494- return True
495- else:
496- return False
497-
498- W, V, _ = schur(M, sort=sort_fun)
513+ W, V, stable_dim = schur(M, sort=sort_fun)
514+ stable_dim
499515```
500516
501517```{code-cell} ipython3
@@ -538,25 +554,24 @@ def stable_solution(M, verbose=True):
538554 表示线性差分方程系统的矩阵。
539555 """
540556 n = M.shape[0] // 2
541- stable_eigvals = list(sorted(np.linalg.eigvals(M))[:n])
557+ tol = 1e-10
542558
543559 def sort_fun(x):
544560 "将模小于1的特征值排序到左上角。"
561+ return abs(x) < 1 - tol
545562
546- if x in stable_eigvals:
547- stable_eigvals.pop(stable_eigvals.index(x))
548- return True
549- else:
550- return False
551-
552- W, V, _ = schur(M, sort=sort_fun)
563+ W, V, stable_dim = schur(M, sort=sort_fun)
564+ if stable_dim != n:
565+ raise ValueError(
566+ f"预期有{n}个模小于1的稳定特征值,但实际找到{stable_dim}个。"
567+ )
553568 if verbose:
554569 print('特征值:\n')
555570 print(' W11: {}'.format(np.diag(W[:n, :n])))
556571 print(' W22: {}'.format(np.diag(W[n:, n:])))
557572
558- # 计算 V21 V11^{-1}
559- P = V[n: , :n] @ np.linalg.inv(V[:n , :n])
573+ # 计算 V21 V11^{-1},不显式构造逆矩阵
574+ P = np.linalg.solve(V[:n , :n].T, V[n: , :n].T).T
560575
561576 return W, V, P
562577
@@ -621,7 +636,7 @@ lq.stationary_values()
621636
622637## 其他应用
623638
624- 上述对潜在不稳定的线性差分方程系统施加稳定性约束的方法,并不限于线性二次型动态最优化问题
639+ 上述对潜在不稳定的线性差分方程系统施加稳定性约束的方法,并不限于线性二次型动态最优化问题。
625640
626641例如,在我们的[线性理性预期模型中的稳定性](https://python.quantecon.org/re_with_feedback.html#another-perspective)讲座中也使用了相同的方法。
627642
677692\hat x_ {t+1} = \hat A \hat x_t + \hat B \hat u_t.
678693$$
679694
680- 那么由 $A, B, R, Q, \beta$ 定义的折现最优控制问题,等价于一个由 $ \hat A, \hat B, Q, R$ 定义的无折现问题。变换前的现最优控制问题的最优策略由 $P, F$ 表征,那么对于变换后的无折现问题,它的解由满足以下方程的 $\hat F, \hat P$ 表征:
695+ 那么由 $A, B, R, Q, \beta$ 定义的折现最优控制问题,其最优策略由 $P, F$ 表征,对应于一个由 $ \hat A, \hat B, Q, R$ 定义的等价无折现问题,其最优策略由满足以下方程的 $\hat F, \hat P$ 表征:
681696
682697$$
683698\hat F=(Q+B'\hat PB)^{-1}\hat B'P \hat A
701716
702717这些设定在上面定义的 `stationary_P` 函数中是默认采用的。
703718
704- ```{code-cell} ipython3
705- β = 1 / (1 + r)
706- lq.beta = β
707- ```
708-
709719```{code-cell} ipython3
710720stationary_P(lq)
711721```
@@ -832,4 +842,3 @@ $$ (eq667)
832842这里我们必须要求 $F$ 满足方程{eq}`eqn:optimalFformula`。
833843
834844方程{eq}`eq666`和{eq}`eq667`为最优值函数提供了不同的视角。
835-
0 commit comments