64. LQ控制的拉格朗日方法#
!pip install quantecon
import numpy as np
from quantecon import LQ
from scipy.linalg import schur
64.1. 概述#
本讲是 线性二次动态规划 的续篇。
它也可以被视为展示了不变子空间（invariant subspace） 技术的讲座，这些技术扩展了我们在此前关于线性理性预期模型中的稳定性中遇到的方法。
我们将介绍一个无限期线性二次无折现动态规划问题的拉格朗日形式。
这类问题有时也被称为最优线性调节器问题。
拉格朗日形式
能够揭示稳定性和最优性之间的联系；
是求解黎卡提方程的快速算法的基础；
为构造那些不直接来源于跨期最优化问题的动态系统解提供了途径。
本讲中的一个关键工具是 \(n \times n\) 辛矩阵的概念。
辛矩阵的特征值以倒数对的形式出现，也就是说，如果 \(\lambda_i \in (-1,1)\) 是一个特征值，那么 \(\lambda_i^{-1}\) 也是特征值。
矩阵特征值成互为倒数配对的这一性质，是一个明显的信号。它表明该矩阵刻画了状态与协态方程系统的联合动态关系。这些方程系统构成了求解无限期线性二次无折现最优控制问题的一阶必要条件。
我们所关心的辛矩阵刻画了最优控制系统的状态与协态向量的一阶动态关系。
在研究这个矩阵的特征值和特征向量时，我们着重分析不变子空间。
这些线性二次动态规划问题的不变子空间表述提供了一个桥梁，连接了递归（即动态规划）公式和线性控制与线性滤波问题的经典公式。后者依赖于矩阵分解方法（参见这节课和这节课）。
虽然本课程大部分内容关注无折现问题，但后面的章节会描述将折现问题转换为无折现问题的便捷方法。
本讲中的这些技术在我们学习这节课中的斯塔克伯格和拉姆齐问题时将会很有用。
64.2. 无折现LQ动态规划问题#
本问题的目标是选择一个控制序列 \(\{u_t\}_{t=0}^\infty\) 来最大化以下准则
约束条件为 \(x_{t+1}=Ax_t+Bu_t\)，其中 \(x_0\) 是给定的初始状态向量。
这里 \(x_t\) 是一个 \((n\times 1)\) 维的状态变量向量，\(u_t\) 是一个 \((k\times 1)\) 维的控制向量，\(R\) 是一个半正定对称矩阵，\(Q\) 是一个正定对称矩阵，\(A\) 是一个 \((n\times n)\) 矩阵，而 \(B\) 是一个 \((n\times k)\) 矩阵。
最优值函数被证明是二次型 \(V(x)= - x'Px\)，其中 \(P\) 是一个半正定对称矩阵。
使用状态转移方程消去下一期的状态变量，贝尔曼方程变为
对于方程 (64.1) 右侧的最大化问题，一阶必要条件是
备注
我们使用以下规则来对二次型和双线性矩阵形式求导： \({\partial x' A x \over \partial x} = (A + A') x; {\partial y' B z \over \partial y} = B z, {\partial y' B z \over \partial z} = B' y\)。
这意味着最优控制规则 \(u\) 可写为
或
其中
将 \(u = - (Q+B'PB)^{-1}B'PAx\) 代入方程 (64.1) 的右侧并重新整理得到
方程 (64.2) 被称为代数矩阵黎卡提方程。
方程 (64.2) 有多个解。
但只有一个解是正定的。
这个正定解对应于本优化问题的最大值。
它将矩阵 \(P\) 表示为矩阵 \(R,Q,A,B\) 的隐函数。
注意，值函数的梯度为
我们稍后将用到公式 (64.3)。
64.3. 拉格朗日量#
对于无折现最优线性调节器问题，构造拉格朗日量
其中 \(2 \mu_{t+1}\) 是时间 \(t\) 转移方程 \(x_{t+1} = A x_t + B u_t\) 的拉格朗日乘子向量。
（我们在 \(\mu_{t+1}\) 前面加上 \(2\) 是为了使其与方程 (64.3) 很好地对应。）
关于 \(\{u_t,x_{t+1}\}_{t=0}^\infty\) 的最大化一阶条件是
定义 \(\mu_0\) 为 \(x_0\) 的影子价格向量，并对 (64.4) 应用包络定理可推导出
这是系统 (64.5) 中第二个方程在时间 \(t=0\) 时的对应形式。
一个重要的事实是
其中 \(P\) 是一个正定矩阵，它是代数黎卡提方程 (64.2) 的解。
因此，根据方程 (64.3) 和 (64.6)，\(- 2 \mu_{t}\) 是值函数对 \(x_t\) 的梯度。
拉格朗日乘子向量 \(\mu_{t}\) 通常被称为与状态向量 \(x_t\) 对应的协态向量。
接下来的推导可按照以下步骤进行：
利用 (64.5) 的第一个方程，将 \(u_t\) 表示为 \(\mu_{t+1}\) 的函数。
将结果代入状态转移方程 \(x_{t+1} = A x_t + B u_t\)。
将得到的方程和 (64.5) 的第二个方程整理成如下形式
其中
当 \(L\) 满秩时（即当 \(A\) 满秩时），我们可以将系统 (64.7) 写作
其中
64.4. 状态-协状态动态#
我们希望求解差分方程系统(64.8)，其解为满足以下条件的序列 \(\{x_t\}_{t=0}^\infty\)：
初始条件为
终端条件为 \(\lim_{t \rightarrow +\infty} x_t =0\)
这个终端条件反映了我们对稳定解的需求，即解不会在 \(t \to \infty\) 时发散。
我们对 \(\{x_t\}\) 序列稳定性的要求，源自于我们希望最大化以下目标函数：
这要求当 \(t \rightarrow + \infty\) 时，\(x_t' R x_t\)收敛于零。
64.5. 倒数对性质#
为了进一步分析，我们研究在(64.9)中定义的 \((2n \times 2n)\) 矩阵 \(M\) 的性质。
为此，我们引入一个 \((2n \times 2n)\) 矩阵
矩阵 \(J\) 的秩为 \(2n\)。
定义： 如果矩阵\(M\)满足
则称其为辛矩阵。
辛矩阵有一些容易验证的的显著性质：
如果 \(M\) 是辛矩阵，那么 \(M^2\) 也是辛矩阵。
如果 \(M\) 是辛矩阵，那么 \(\textrm{det}(M) = 1\)。
可以直接验证，方程(64.9)中的 \(M\) 是辛矩阵。
从方程(64.10)和 \(J^{-1} = J^\prime = -J\) 可以推导出，对任意辛矩阵\(M\)，有
方程(64.11)表明 \(M^\prime\) 与 \(M\) 的逆矩阵通过相似变换相互关联。
对于方阵，我们回顾以下事实：
相似矩阵具有相同的特征值
矩阵的逆的特征值是该矩阵特征值的倒数
一个矩阵和它的转置矩阵具有相同的特征值
从方程(64.11)可以得出，\(M\) 的特征值成倒数对：如果 \(\lambda\) 是 \(M\) 的特征值，那么 \(\lambda^{-1}\) 也是。
将方程(64.8)写作
其中 \(y_t = \begin{pmatrix}x_t\cr \mu_t\cr\end{pmatrix}\)。
考虑 \(M\) 的三角化
其中
右侧的每个块都是 \((n\times n)\) 维的
\(V\) 是非奇异的
\(W_{22}\) 的所有特征值的模都大于1
\(W_{11}\) 的所有特征值的模都小于1
64.6. Schur分解#
Schur分解和特征值分解是两种形如(64.13)的分解。
将方程(64.12)写作
方程(64.14)对任意初始条件\(y_0\)的解显然为
其中，当 \(t=1\) 时 \(W_{12,t} = W_{12}\)；对于 \(t \geq 2\)，\(W_{12,t}\) 遵循递归关系
这里 \(W^t_{ii}\) 是 \(W_{ii}\) 的 \(t\) 次幂。
将方程(64.15)写作
其中 \(y^\ast_t = V^{-1} y_t\)，特别地
这里 \(V^{ij}\) 表示分块矩阵 \(V^{-1}\) 的 \((i,j)\) 部分。
由于 \(W_{22}\) 是一个不稳定矩阵，\(y^\ast_t\) 将发散，除非 \(y^\ast_{20} = 0\)。
令 \(V^{ij}\) 表示分块矩阵 \(V^{-1}\) 的第 \((i,j)\) 块。
为了获得稳定性，我们必须设定 \(y^\ast_{20} =0\)，根据方程 (64.16) 可得
或
该方程在时间上递推成立，也就是说，它意味着
但注意，由于 \((V^{21} \ V^{22})\) 是 \(V\) 的逆矩阵的第二行块，因此
这意味着
因此，
所以我们可以写成
和
然而，我们知道 \(\mu_t = P x_t\)，其中 \(P\) 出现在求解黎卡提方程的矩阵中。
因此，上述推导表明
值得注意的是，公式(64.17)为我们提供了一种高效的方法来计算正定矩阵 \(P\)，该矩阵可以解决动态规划中出现的代数黎卡提方程(64.2)。
即使 \(M\) 的特征值不以倒数对的形式出现，这种方法也可以用来计算任何形如(64.8)的系统的解（如果解存在的话）。
只要 \(M\) 的特征值一半位于单位圆内、一半位于单位圆外，这种方法通常都能奏效。
当系统对应于一个存在扭曲的模型均衡，而这些扭曲使得均衡不再是某个最优问题的解时，经过折现调整后的特征值将不再成互为倒数对。参见[Ljungqvist and Sargent, 2018]第12章。
64.7. 应用#
这里我们通过一个示例来演示计算过程，这个示例是从quantecon讲座中借鉴的确定性版本。
# 模型参数
r = 0.05
c_bar = 2
μ = 1
# 构建为LQ问题
Q = np.array([[1]])
R = np.zeros((2, 2))
A = [[1 + r, -c_bar + μ],
[0, 1]]
B = [[-1],
[0]]
# 构造一个LQ实例
lq = LQ(Q, R, A, B)
给定矩阵 \(A\)、\(B\)、\(Q\)、\(R\)，我们可以计算 \(L\)、\(N\) 和 \(M=L^{-1}N\)。
def construct_LNM(A, B, Q, R):
n, k = lq.n, lq.k
# 构造 L 和 N
L = np.zeros((2*n, 2*n))
L[:n, :n] = np.eye(n)
L[:n, n:] = B @ np.linalg.inv(Q) @ B.T
L[n:, n:] = A.T
N = np.zeros((2*n, 2*n))
N[:n, :n] = A
N[n:, :n] = -R
N[n:, n:] = np.eye(n)
# 计算 M
M = np.linalg.inv(L) @ N
return L, N, M
L, N, M = construct_LNM(lq.A, lq.B, lq.Q, lq.R)
M
array([[ 1.05 , -1. , -0.95238095, 0. ],
[ 0. , 1. , 0. , 0. ],
[ 0. , 0. , 0.95238095, 0. ],
[ 0. , 0. , 0.95238095, 1. ]])
让我们验证 \(M\) 是辛矩阵。
n = lq.n
J = np.zeros((2*n, 2*n))
J[n:, :n] = np.eye(n)
J[:n, n:] = -np.eye(n)
M @ J @ M.T - J
array([[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.]])
我们可以使用np.linalg.eigvals计算矩阵\(M\)的特征值，并按升序排列。
eigvals = sorted(np.linalg.eigvals(M))
eigvals
[np.float64(0.9523809523809523),
np.float64(1.0),
np.float64(1.0),
np.float64(1.05)]
当我们应用Schur分解使得\(M=V W V^{-1}\)时，我们希望
\(W\)的左上块\(W_{11}\)的所有特征值的模小于1，并且
右下块\(W_{22}\)的特征值的模大于1。
为了得到我们想要的结果，让我们定义一个排序函数，告诉scipy.schur将模小于1的对应特征值排序到左上角。
stable_eigvals = eigvals[:n]
def sort_fun(x):
"将模小于1的特征值排序到左上角。"
if x in stable_eigvals:
stable_eigvals.pop(stable_eigvals.index(x))
return True
else:
return False
W, V, _ = schur(M, sort=sort_fun)
W
array([[ 1. , -0.02316402, -1.00085948, -0.95000594],
[ 0. , 0.95238095, -0.00237501, -0.95325452],
[ 0. , 0. , 1.05 , 0.02432222],
[ 0. , 0. , 0. , 1. ]])
V
array([[ 0.99875234, 0.00121459, -0.04992284, 0. ],
[ 0.04993762, -0.02429188, 0.99845688, 0. ],
[ 0. , 0.04992284, 0.00121459, 0.99875234],
[ 0. , -0.99845688, -0.02429188, 0.04993762]])
我们可以检查 \(W_{11}\) 和 \(W_{22}\) 的特征值的模。
由于它们都是三角矩阵，特征值就是对角线元素。
# W11
np.diag(W[:n, :n])
array([1. , 0.95238095])
# W22
np.diag(W[n:, n:])
array([1.05, 1. ])
以下函数封装了 \(M\) 矩阵构建、Schur分解以及稳定性约束下的 \(P\) 计算。
def stable_solution(M, verbose=True):
"""
给定一个线性差分方程系统
y' = |a b| y
x' = |c d| x
该系统可能不稳定，通过施加稳定性约束来求解。
参数
---------
M : np.ndarray(float)
表示线性差分方程系统的矩阵。
"""
n = M.shape[0] // 2
stable_eigvals = list(sorted(np.linalg.eigvals(M))[:n])
def sort_fun(x):
"将模小于1的特征值排序到左上角。"
if x in stable_eigvals:
stable_eigvals.pop(stable_eigvals.index(x))
return True
else:
return False
W, V, _ = schur(M, sort=sort_fun)
if verbose:
print('特征值：\n')
print(' W11: {}'.format(np.diag(W[:n, :n])))
print(' W22: {}'.format(np.diag(W[n:, n:])))
# 计算 V21 V11^{-1}
P = V[n:, :n] @ np.linalg.inv(V[:n, :n])
return W, V, P
def stationary_P(lq, verbose=True):
"""
计算表示值函数的矩阵 :math:`P`
V(x) = x' P x
在无限时域情况下。通过在解路径上施加稳定性约束
并使用舒尔分解来进行计算。
参数
----------
lq : qe.LQ
用于分析无限时域形式的线性二次最优控制问题的
QuantEcon类。
返回值
-------
P : array_like(float)
值函数表示中的P矩阵。
"""
Q = lq.Q
R = lq.R
A = lq.A * lq.beta ** (1/2)
B = lq.B * lq.beta ** (1/2)
n, k = lq.n, lq.k
L, N, M = construct_LNM(A, B, Q, R)
W, V, P = stable_solution(M, verbose=verbose)
return P
# 计算P
stationary_P(lq)
特征值：
W11: [1. 0.95238095]
W22: [1.05 1. ]
array([[ 0.1025, -2.05 ],
[-2.05 , 41. ]])
请注意，以这种方式计算得到的矩阵 \(P\) 与 quantecon 中通过迭代 Riccati 差分方程直至收敛来求解代数 黎卡提方程的程序得到的结果非常接近。
这种微小的差异来自计算误差，可以通过增加最大迭代次数或降低收敛容差来减小。
lq.stationary_values()
(array([[ 0.1025, -2.05 ],
[-2.05 , 41.01 ]]),
array([[-0.09761905, 1.95238095]]),
0)
使用Schur分解效率要高得多。
%%timeit
stationary_P(lq, verbose=False)
90.1 μs ± 337 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
%%timeit
lq.stationary_values()
1.34 ms ± 8.18 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
64.8. 其他应用#
上述对潜在不稳定的线性差分方程系统施加稳定性约束的方法，并不限于线性二次型动态最优化问题
例如，在我们的线性理性预期模型中的稳定性讲座中也使用了相同的方法。
让我们尝试使用本讲中定义的stable_solution函数来求解该讲座中描述的模型。
def construct_H(ρ, λ, δ):
"根据参数构建矩阵H。"
H = np.empty((2, 2))
H[0, :] = ρ,δ
H[1, :] = - (1 - λ) / λ, 1 / λ
return H
H = construct_H(ρ=.9, λ=.5, δ=0)
W, V, P = stable_solution(H)
P
特征值：
W11: [0.9]
W22: [2.]
array([[0.90909091]])
64.9. 折现问题#
64.9.1. 转换状态和控制以消除折现#
有一对十分有用的变量变换，可以将一个含折现因子的最优控制问题转化为一个无折现问题。
假设我们有一个折现问题，其目标函数为
且状态转移方程仍为 \(x_{t +1 }=Ax_t+Bu_t\)。
定义变换后的状态和控制变量
\(\hat x_t = \beta^{\frac{t}{2}} x_t \)
\(\hat u_t = \beta^{\frac{t}{2}} u_t\)
以及变换后的转移方程矩阵
\(\hat A = \beta^{\frac{1}{2}} A\)
\(\hat B = \beta^{\frac{1}{2}} B \)
使得调整后的状态和控制变量满足以下转移规律：
那么由 \(A, B, R, Q, \beta\) 定义的折现最优控制问题，等价于一个由 \(\hat A, \hat B, Q, R\) 定义的无折现问题。变换前的现最优控制问题的最优策略由 \(P, F\) 表征，那么对于变换后的无折现问题，它的解由满足以下方程的 \(\hat F, \hat P\) 表征：
和
从 \(\hat A, \hat B\) 的定义可以直接得出 \(\hat F = F\) 和 \(\hat P = P\)。
通过利用这种变换，我们可以将一个贴现问题转化为一个等价的无贴现问题来求解。
特别地，我们可以先将一个带折现因子的线性二次问题转化为无折现形式，然后利用前文介绍的拉格朗日与不变子空间方法求解相应的无折现最优调节器问题。
例如，当 \(\beta=\frac{1}{1+r}\) 时，我们可以用 \(\hat{A}=\beta^{1/2} A\) 和 \(\hat{B}=\beta^{1/2} B\) 求解 \(P\)。
这些设定在上面定义的 stationary_P 函数中是默认采用的。
β = 1 / (1 + r)
lq.beta = β
stationary_P(lq)
特征值：
W11: [0.97590007 0.97590007]
W22: [1.02469508 1.02469508]
array([[ 0.0525, -1.05 ],
[-1.05 , 21. ]])
我们可以验证该解与使用 quantecon 包中的 LQ.stationary_values 例程得到的解是一致的。
lq.stationary_values()
(array([[ 0.0525, -1.05 ],
[-1.05 , 21. ]]),
array([[-0.05, 1. ]]),
np.float64(0.0))
64.9.2. 折现问题的拉格朗日量#
出于多种目的，有必要简要地明确描述折现问题的拉格朗日量。
因此，对于折现最优线性调节器问题，构建拉格朗日量：
其中 \(2 \mu_{t+1}\) 是状态向量 \(x_{t+1}\) 的拉格朗日乘子向量。
对于 \(\{u_t,x_{t+1}\}_{t=0}^\infty\) 的最大化一阶条件是：
定义 \(2 \mu_0\) 为 \(x_0\) 的影子价格向量，并对(64.18)应用包络定理可得：
这是系统(64.19)中第二个方程在 \(t=0\) 时刻的对应形式。
按照上述无折现系统(64.5)的处理方法，我们可以将一阶条件重新整理为如下系统:
在特殊情况\(\beta = 1\)时，该式与方程(64.5)一致，这是符合预期的。
通过观察系统(64.20)，我们可以推断出一些揭示最优线性调节器问题结构的恒等式。在这篇讲座中，当我们应用和扩展本讲座的方法来研究斯塔克尔伯格和拉姆齐问题时，其中一些恒等式会很有用。
首先,注意方程系统(64.20)的第一个区块表明，当\(\mu_{t+1} = P x_{t+1}\)时，
可以重新整理为
该式所表示的状态最优闭环动态必须与我们通过动态规划推导得到的另一种表达式一致，即：
但使用
可以推导出
因此，我们的两个闭环动态表达式当且仅当满足以下条件时才相等:
矩阵方程(64.22)可以通过应用分块矩阵求逆公式来验证。
备注
选择适当的矩阵 \(a, b, c, d\)， 并使用公式\((a - b d^{-1} c)^{-1} = a^{-1} + a^{-1} b (d - c a^{-1} b)^{-1} c a^{-1}\)即可。
接下来，注意对于任何固定的 \(F\)，只要 \(A-BF\) 的特征值的模小于 \(\frac{1}{\beta}\)，使用这个规则的值函数永远是 \(-x_0 \tilde P x_0\)，其中 \(\tilde P\) 满足以下矩阵方程：
显然，只有当 \(F\) 满足公式(64.21)时，\(\tilde P = P\) 才成立。
接下来，注意系统(64.20)的第二个方程暗示了拉格朗日乘数的”前瞻”方程
其解为
其中
这里我们必须要求 \(F\) 满足方程(64.21)。