CS229 Lesson 19 微分动态规划

课程视频地址：http://open.163.com/special/opencourse/machinelearning.html

课程主页：http://cs229.stanford.edu/

更具体的资料链接：https://www.jianshu.com/p/0a6ef31ff77a

笔记参考自中文翻译版：https://github.com/Kivy-CN/Stanford-CS-229-CN

这一讲介绍了LQR以及Ricatti方程。

3.从非线性动力系统到LQR

事实证明，即使动力系统是非线性的，很多问题也可以简化到LQR。 虽然LQR是一个很好的形式，因为我们能够得到一个很好的精确解决方案，但它不够一般化。 让我们以倒立摆的情况为例。 状态之间的转换看起来像

$$\left( \begin{matrix} x_{t+1} \\ \dot x_{t+1} \\ \dot \theta_{t+1}\\ \dot \theta_{t+1} \end{matrix} \right)= F \left( \begin{matrix} \left( \begin{matrix} x_{t+1} \\ \dot x_{t+1} \\ \dot \theta_{t+1}\\ \dot \theta_{t+1} \end{matrix} \right),a_t \end{matrix} \right)$$

其中函数$F$取决于角度的余弦等等。现在，我们可能会问的问题是，我们能否线性化这个系统？

3.1 线性化动力系统

让我们假设在时间$t$，系统将大部分时间花在某个状态$\bar s_t$上并且执行的操作是$\bar a_t$。 对于倒立摆，如果我们达到某种最佳状态，这是正确的：我们的行动很小并且我们不会偏离竖直方向。

​ 我们将使用泰勒展开来线性化动力系统。 在状态为一维且转换函数$F$不依赖于动作的简单情况下，我们将得到

$$s_{t+1}=F(s_t)\approx F(\bar s_t) + F'(\bar s_t).(s_t -\bar s_t)$$

​ 在更一般的设定下，公式看起来相同，只不过是使用了梯度而不是导数

$$\begin{eqnarray*} s_{t+1}\approx F(\bar s_t,\bar a_t) +\nabla_s F(\bar s_t,\bar a_t).(s_t -\bar s_t) +\nabla_a F(\bar s_t,\bar a_t).(a_t -\bar a_t) \tag{3} \end{eqnarray*}$$

​ 现在，$s_{t+1}​$关于$s_t​$和$a_t​$是线性的，因为我们可以将等式（3）重写为

$$s_{t+1}\approx As_t + Bs_t + \kappa$$

​ 其中$\kappa$是某个常数，$A,B$是矩阵。 现在，这种写法看起来非常类似于LQR的假设。 我们只需要抛弃不变的常数$\kappa$！ 事实证明，常数项可以通过人为增加一个维度将其合并到$s_t$。 这与我们在课程开始时用于线性回归的技巧相同……

3.2 差分动态规划（DDP）

之前的方法适用于目标是留在某些状态的情况（想想倒立摆，或者车辆必须停留在车道中间）。 但是，在某些情况下，目标可能更复杂。

​ 我们将介绍一种适用于我们的系统必须遵循某些轨迹（例如考虑火箭）的方法。 这种方法将轨迹离散化为离散的时间步长，并创建中间目标，我们将能够使用以前的技术！ 该方法称为差分动态规划。 主要步骤是

• 步骤1：使用控制器产生一个虚拟轨迹，它接近我们想要遵循的轨迹。 换句话说，我们的控制器能够逼近完美轨迹按照如下方式

  $$s_0' ,a_0' \to s_1' ,a_1' \to ...$$

• 步骤2：线性化每个轨迹点$s_t’$周围的动态，换句话说

  $$s_{t+1}\approx F(s_t',a_t') +\nabla_s F( s_t',a_t').(s_t -s_t') +\nabla_a F( s_t',a_t').(a_t -a_t')$$

  其中$s_t,a_t$是我们目前的状态和行动。 现在我们对每个点都有一个线性近似，我们可以使用上一节并将上式重写为

  $$s_{t+1}\approx A_ts_t + B_ts_t$$

  （注意，在这种情况下，我们使用我们在这部分讲义开头提到的非平稳动力学设定）

  注意我们可以对奖励$R(t)$应用类似的推导，并使用二阶泰勒展开：

  $$\begin {aligned} R(s_t, a_t) \approx& R(s_t', a_t') +\nabla_s R(s_t', a_t')(s_t -s_t') +\nabla_a R(s_t', a_t')(a_t -a_t')\\ +& \frac 12 (s_t -s_t')^T H_{ss} (s_t -s_t') + (s_t -s_t')^T H_{sa} (a_t -a_t') \\ +& \frac 12 (a_t -a_t')^T H_{aa} (a_t -a_t') \end {aligned}$$

  其中$H_{xy}$是$R$在$(s_t’,a_t’)$处关于$x,y$的Hessian的矩阵。 通过添加额外维度的技巧，上式可以重写为

  $$R(s_t, a_t)=-s_t^T U_t s_t -a_t^T W_t a_t$$

对于某些矩阵$U_t,W_t$。这是因为如下简单事实：

$$\left(\begin{matrix} 1 & x \end{matrix} \right) . \left(\begin{matrix} a & b\\ b & c \end{matrix} \right). \left(\begin{matrix} 1 \\ x \end{matrix} \right) =a+2bx + cx^2$$

• 步骤3：现在，你可以说服自己，我们的问题在LQR框架中已经严格重写。 让我们使用LQR找出最优策略$\pi_t$。 因此，我们的新控制器（希望）会更好！

  注意：如果LQR轨迹与轨迹的线性化近似偏差太大，可能会出现一些问题，但可以通过奖励形变来解决…

• 步骤4：现在我们得到了一个新的控制器（我们的新策略$\pi_t$），我们用它来产生一个新的轨迹

  $$s_0' ,\pi_0(s_0') \to s_1' ,\pi_1(s_1') \to \ \to s_T'$$

  请注意，当我们生成这个新轨迹时，我们使用真实的$F$而不是它的线性近似来计算转换，这意味着

  $$s_{t+1}'= F(s_t' ,a_t')$$

  然后，回到步骤2并重复直到某些停止标准。

4.线性二次高斯(LQG)

通常，在真实的世界中，我们无法观察到完整的状态。 例如，自动驾驶汽车可以从摄像机接收图像，这仅仅是观察，而不是世界的完整状态。 到目前为止，我们假设状态已知。 由于这可能不适用于大多数现实问题，我们需要一种新工具来模拟这种情况：部分可观察的MDP(Partially Observable MDP)。

​ POMDP是具有额外观察层的MDP。 换句话说，我们引入一个新的变量$o_t$，给定当前状态$s_t$，$o_t$服从条件分布

$$o_t | s_t \sim O(o|s)$$

​ 形式上，有一个有限范围的POMDP由如下元组给出

$$(\mathcal S,\mathcal O, P_{sa},T, R)$$

​ 在这个框架内，一般策略是在观测结果$o_1,…,o_t$上维持信念状态（在状态上的分布）。 然后，POMDP中的策略将信念状态映射到动作。

​ 在本节中，我们将LQR扩展到此设定下。假设我们观察到$y_t \in \mathbb R^m$，其中$m<n$，并且

$$\begin{cases} y_t &=C.s_t +v_t\\ s_{t+1}& =A.s_t + B.a_t +w_t \end{cases}$$

其中$C\in\mathbb R^{m\times n}$是压缩矩阵，$v_t$是传感器噪声（像$w_t$一样也服从高斯分布）。 注意，奖励函数$R^{(t)}$保持不变，仍然作为状态（不是观测值）和动作的函数。 此外，由于噪声服从高斯分布的，因此信念状态也服从是高斯分布。 在这个新框架中，让我们概述一下我们将采用哪种方法来找到最优策略：

• 步骤1：首先，根据我们的观察结果计算可能状态（信念状态）的分布。 换句话说，我们想要计算如下分布的平均值$s_{t|t}​$和协方差$\Sigma_{t|t}​$

  $$s_t | y_1,...,y_t \sim \mathcal N(s_{t|t},\Sigma_{t|t} )$$

  为了有效地执行计算，我们将使用卡尔曼滤波算法（这个算法在阿波罗登月计划中也被使用）。

• 步骤2：现在我们有了分布，我们将使用平均值$s_{t|t}​$作为$s_t​$的最佳近似值

• 步骤3：然后将动作设置为$a_t:= L_t s_{t|t}$，其中$L_t$来自常规的LQR算法。

​ 直观地，要理解为什么这是有效的，请注意$s_{t|t}$是$s_t$的噪声近似（相当于向LQR添加更多噪声），但我们证明了LQR与噪声无关！

​ 步骤1需要详细说明。 我们将介绍一个简单的情况，即我们的动力系统中不依赖动作（但一般情形遵循相同的想法）。 假设

$$\begin{cases} s_{t+1}& =A.s_t +w_t, w_t \sim \mathcal N(0,\Sigma_s)\\ y_t &=C.s_t +v_t,v_t \sim \mathcal N(0,\Sigma_y) \end{cases}$$

​ 由于噪声是服从高斯分布，我们可以很容易地证明联合分布也服从高斯分布（利用高斯分布的线性变换也服从高斯分布）

$$\left( \begin{matrix} s_1 \\ \vdots \\ s_t \\ y_1 \\ \vdots \\ y_t \\ \end{matrix} \right) \sim \mathcal N (\mu, \Sigma) \\ 对某些\mu,\Sigma$$

​ 然后，使用高斯分布的边际分布公式（参见因子分析讲义），我们会得到

$$s_t | y_1,...,y_t \sim \mathcal N(s_{t|t}, \Sigma_{t|t})$$

​ 但是，使用该公式计算边际分布的参数计算成本很高！ 它需要对$t\times t$矩阵做操作。 回想一下，求逆矩阵可以在$O(t^3)$中完成，然后该步骤需要执行$t$次，从而产生$O(t^4)$的时间成本！

​ 卡尔曼滤波器算法提供了一种更好的计算均值和方差的方法，该算法基于两个基本步骤。 假设我们已知$s_t | y_1,…,y_t​$，然后执行如下两个步骤

• 预测步骤计算$s_{t+1} | y_1,…,y_t​$

• 更新步骤计算$s_{t+1} | y_1,…,y_t,y_{t+1}$

​ 预测和更新步骤的组合更新了我们的信念状态。 换句话说，这个过程看起来像

$$(s_t | y_1,...,y_t)\overset{\text{predict}}\longrightarrow (s_{t+1} | y_1,...,y_t)\overset{\text{update}}\longrightarrow (s_{t+1} | y_1,...,y_t,y_{t+1})\overset{\text{predict}}\longrightarrow...$$

下面分别讨论这两个步骤。

预测步骤：

假设我们已知分布

$$s_t | y_1,...,y_t \sim \mathcal N(s_{t|t}, \Sigma_{t|t})$$

因为

$$s_{t+1} =A.s_t +w_t, w_t \sim \mathcal N(0,\Sigma_s)$$

所以$s_{t+1} | y_1,…,y_t​$也服从高斯分布，所以我们只要计算其均值方差即可

$$\begin{aligned} s_{t+1|t}&=\mathbb E[s_{t+1} | y_1,...,y_t] \\ &=\mathbb E[A.s_t +w_t| y_1,...,y_t] \\ &=A. \mathbb E[s_t| y_1,...,y_t] + \mathbb E[w_t]\\ &=A. s_{t|t}\\ \Sigma_{t+1|t} &=\mathbb E[(s_{t+1}-s_{t+1|t})(s_{t+1}-s_{t+1|t})^T ] \\ &=\mathbb E[(A.s_t +w_t-A. s_{t|t})(A.s_t +w_t-A. s_{t|t})^T ] \\ &=A.\mathbb E[(s_t- s_{t|t})(s_t- s_{t|t})^T ]. A^T + \mathbb E[w_t w_t^T] \\ &=A. \Sigma_{t|t}. A^T + \Sigma_s \end{aligned}$$

从而

$$s_{t+1} | y_1,...,y_t \sim \mathcal N(s_{t+1|t}, \Sigma_{t+1|t})$$

其中

$$\begin{cases} s_{t+1|t} &= A. s_{t|t} \\ \Sigma_{t+1| t} &= A. \Sigma_{t|t}. A^T + \Sigma_s \end{cases}$$

更新步骤：

给定$s_{t+1|t},\Sigma_{t+1| t}​$使得

$$s_{t+1} | y_1,...,y_t \sim \mathcal N(s_{t+1|t}, \Sigma_{t+1|t})$$

记$b_t= s_{t+1} | y_1,…,y_t$，考虑向量$ \left( \begin{matrix} b_t\\ y_{t+1}\end{matrix} \right) $的分布，因为

$$y_{t+1} =C.s_{t+1} +v_{t+1},v_{t+1} \sim \mathcal N(0,\Sigma_y)$$

所以不难看出上述向量服从高斯分布，下面计算其均值以及协方差。

均值：

$$\begin{aligned} \mathbb E[b_t]&= s_{t+1|t}\\ \mathbb E[y_{t+1}] &= \mathbb E[\mathbb E[y_{t+1} |y_1,...,y_t]]\\ &=\mathbb E[\mathbb E[C.s_{t+1} +v_{t+1} |y_1,...,y_t]] \\ &=C.\mathbb E[\mathbb E[s_{t+1} |y_1,...,y_t]] \\ &= C. s_{t+1|t} \end{aligned}$$

协方差：

$$\begin{aligned} \text{Cov}(b_t, b_t)&= \Sigma_{t+1|t} \\ \text{Cov}(b_t, y_{t+1}) &= \mathbb E[(b_t-s_{t+1|t})(y_{t+1}-C. s_{t+1|t})^T]\\ &=\mathbb E[(b_t-s_{t+1|t})(C.s_{t+1} +v_{t+1}-C. s_{t+1|t})^T] \\ &=\mathbb E[\mathbb E[(b_t-s_{t+1|t})(C.s_{t+1} +v_{t+1}-C. s_{t+1|t})^T|y_1,...,y_t]] \\ &=\mathbb E[\mathbb E[(b_t-s_{t+1|t})(C.b_t +v_{t+1}-C. s_{t+1|t})^T]] \\ &=\mathbb E[\mathbb E[(b_t-s_{t+1|t})(C(b_t-s_{t+1|t}))^T]]+ \mathbb E[\mathbb E[(b_t-s_{t+1|t})v_{t+1}]]\\ &=\mathbb E[\mathbb E[(b_t-s_{t+1|t})(b_t -s_{t+1|t})^T]].C^T (由独立性可得第二项为0) \\ &=\Sigma_{t+1|t}.C^T\\ \text{Cov}(y_{t+1}, y_{t+1}) &= \mathbb E[(y_{t+1}-C. s_{t+1|t})(y_{t+1}-C. s_{t+1|t})^T]\\ &=\mathbb E[(C.s_{t+1} +v_{t+1}-C. s_{t+1|t})(C.s_{t+1} +v_{t+1}-C. s_{t+1|t})^T] \\ &=\mathbb E[\mathbb E[(C.s_{t+1} +v_{t+1}-C. s_{t+1|t})(C.s_{t+1} +v_{t+1}-C. s_{t+1|t})^T|y_1,...,y_t]] \\ &=\mathbb E[C.\mathbb E[(b_t - s_{t+1|t})(b_t -s_{t+1|t})^T]C^T + \mathbb E[v_{t+1}v_{t+1}^T]] \\ &=C.\Sigma_{t+1|t}.C^T + \Sigma_y \end{aligned}$$

所以

$$\left( \begin{matrix} b_t\\ y_{t+1}\end{matrix} \right) \sim \mathcal N( \mu_t,\Sigma_t )\\ \mu_t= \left( \begin{matrix} s_{t+1|t} \\ C. s_{t+1|t} \end{matrix} \right), \Sigma_t = \left( \begin{matrix} \Sigma_{t+1|t}& \Sigma_{t+1|t}.C^T\\ C.\Sigma_{t+1|t} &C.\Sigma_{t+1|t}.C^T + \Sigma_y \end{matrix} \right)$$

在因子分析章节中，我们有如下结论：

如果

$$x=\left( \begin{matrix} x_1 \\ x_2 \end{matrix} \right),x\sim \mathcal N(\mu, \Sigma)\\\mu=\left( \begin{matrix} \mu_1 \\ \mu_2 \end{matrix} \right) , \Sigma= \left( \begin{matrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} &\Sigma_{22} \end{matrix} \right)$$

那么

$$x_1 |x_2 \sim \mathcal N(\mu_{1|2},\Sigma_{1|2}) \\ 其中 \mu_{1|2} = \mu_1 + \Sigma_{12}\Sigma_{22}^{-1}(x_2 -\mu_2) \\ \Sigma_{1|2} = \Sigma_{11}- \Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}$$

注意到此处

$$b_t | y_{t+1}= (s_{t+1} | y_1,...,y_t)| y_{t+1}= s_{t+1} | y_1,...,y_t,y_{t+1}$$

所以

$$s_{t+1} | y_1,...,y_{t+1} \sim \mathcal N(s_{t+1|t+1},\Sigma_{t+1|t+1})$$

其中

$$\begin{cases} s_{t+1|t+1} & = s_{t+1|t} +K_t(y_{t+1} -C s_{t+1|t})\\ \Sigma_{t+1|t+1} & = \Sigma_{t+1|t} -K_t .C.\Sigma_{t+1|t} \end{cases}$$

且

$$K_t : = \Sigma_{t+1|t} C^T(C\Sigma_{t+1|t} C^T +\Sigma_y)^{-1}$$

矩阵$K_t$被称为卡尔曼收益。

​ 现在，如果我们仔细研究一下公式，我们会注意到在第$t​$个时间戳之前我们不需要观测值！更新步骤仅取决于先前的分布。 将所有这些放在一起，算法首先运行前向传递来计算$K_t​$，$\Sigma_{t|t}​$和$s_{t|t}​$（在文献中有时称为$\hat s​$）。 然后，它运行一个反向传递（LQR更新）来计算$\Psi_t,\Phi_t​$和$L_t​$。最后，我们用$a_t’= L_t s_{t|t}​$恢复最优策略。