一种新的基于数据驱动的神经动态规划方法

A New Data-Driven Neural Dynamic Programming Algorithm

摘要: 为了实现无模型离散时间非线性动态系统的最优控制，提出了一种新的基于数据驱动的神经动态规划方法。该方法利用Q函数的残差与基函数的内积为零，同时控制策略的残差与基函数的内积也为零，从而得到控制方程。接着使用离线数据集与在线数据来迭代更新神经网络的系数，从而得到近似最优的控制策略，本文还证明了该算法是收敛的。

Abstract: A new data-driven neural dynamic programming method for model-free discrete-time nonlinear dynamic system is proposed in this paper. The residual of the Q-function and the control strategy are operated to be zero with the basis function through the inner product. Then the coefficients of the neural network are updated by the offline trained data and the online data. Finally the optimal control strategy is obtained and the convergence of this algorithm is proved.

1. 引言

无论是在控制理论还是工程领域，最优化控制问题都是一个重难点。最优控制问题的求解依赖于求解HJB方程，由于系统一般是非线性的，导致对于求解HJB方程很困难，难以得到解析解 [1] 。

Bellman在1957年首次提出了动态规划方法。动态规划对于求解最优控制问题是一个行之有效的方法。该方法适应广泛：离散系统，连续系统，非线性系统，以及随机系统等等。但是该方法是一种逆向计算方法，随着解的时间维度增加，且由于控制变量维数大容易造成“维数灾难”问题。于是近年来，针对未知的复杂非线性系统的最优控制问题，自适应动态规划应运而生。与传统的动态规划相比，由于其采用函数近似结构来逼近系统模型，评价指标 [2] ，和控制策略，避免了因时间逆向计算造成的“维数灾难”问题。自适应动态规划可用于大规模非线性系统 [3] ，其应用也相当广泛；例如：机器人运动规划 [4] ，无人驾驶汽车 [5] ，污水处理系统 [6] 等等。

LIU等采用值迭代的自适应动态规划，并证明其收敛条件是迭代性能指标函数初始化为任意半正定函数 [7] 。根据此收敛条件，提出了一种基于自适应动态规划的协同优化算法 [8] 。该协同优化算法令迭代的残差快速减小，大幅提高了自适应动态规划的收敛速度。LIU等对于离散的非线性系统采用基于策略迭代的自适应动态规划方法 [9] ，证明了该迭代控制律可以使系统稳定，并说明其是收敛的。LUO等则分别对于连续时间非线性系统采用基于策略迭代的自适应动态规划方法 [10] ，该策略迭代的实现使用基于神经网络的最小二乘方法并证明该方法是收敛的。文献 [11] 针对离散时间非线性系统，提出了策略梯度自适应动态规划算法，并证明了该算法中Q函数能收敛到最优值。

由于神经网络良好的自适应与自学习等特点 [12] ，本文采用神经网络与动态规划结合的自适应动态规划对无模型的最优控制问题进行求解。对于神经网络的权重系数采用离线数据集与在线数据结合进行更新。

2. 问题描述

本文所描述的离散时间非线性系统为：

$x_{k+1}=f\left(x_k+u_k\right),k=0,1,2,\cdots$ (1)

其中 $x_k\in R^n$ 为状态向量， $u_k\in R^m$ 为控制向量。此外X和U为一个完备集合， $D=\left\{\left(x,u\right)|x\in X,u\in U\right\}$ 。系统(1)在X上稳定，函数 $f\left(x,u\right)$ 在集合D上连续，且 $f\left(0,0\right)=0$ ， $x=0$ 是系统(1)的一个稳定状态。本文考虑的是无模型的最优控制问题，即函数 $f\left(x,u\right)$ 的具体解析式是不知道的。优化目标是找到一个稳定的反馈控制 $u_k=u\left(x_k\right)$ ，使得性能指标函数(2)最小。

$V_u\left(x_0\right)={\displaystyle {\sum }_{l=0}^{\infty }P\left(x_l,u_l\right)}.$ (2)

其中 $P\left(x,u\right)=S\left(x\right)+W\left(u\right)$ ，和 $W\left(u\right)$ 为正定函数。

3. 策略梯度自适应动态规划

对于上述最优控制问题，由于该系统是非线性的且系统模型 $f\left(x,u\right)$ 的解析式并不知道。为了克服这些困难，引入了策略梯度自适应动态规划来求解该最优控制问题。

定义1对于系统(1)，如果 $u\left(x\right)$ 在X上连续，且 $\forall x\in X$ ， $V_u\left(x\right)<\infty$ ，则称控制 $u\left(x\right)$ 是可控的，记 $u\left(x\right)\in U\left(X\right)$ 。基于自适应动态规划理论，定义值函数 $V_u\left(x_k\right)$ 为：

$V_u\left(x_k\right)={\displaystyle {\sum }_{l=k}^{\infty }P\left(x_l,u_l\right)}.$ (3)

易知：

$V_u\left(x_k\right)=P\left(x_k,u_k\right)+V_u\left(x_{k\text{+1}}\right).$ (4)

记最优值函数为：

$V^{*}\left(x\right)=\underset{u}{\min }{V}_u\left(x\right).$ (5)

对应的最优控制为：

$u^{\text{*}}\left(x\right)=\arg \underset{u}{\min }{V}_u\left(x_k\right).$ (6)

根据自适应动态规划，为了更好求解该最优控制问题，下面引入Q函数：

$Q_u\left(x_k,a\right)=P\left(x_k,a\right)+{\displaystyle {\sum }_{l=k+\text{1}}^{\infty }P\left(x_l,u_l\right)}.$ (7)

且易知 $Q_u\left(x_k,u\right)=V_u\left(x_k\right)$ ，Q函数也可以表示为：

$Q_u\left(x_k,a\right)=P\left(x_k,a\right)+Q_u\left(x_{k+1},u\right)=P\left(x_k,a\right)+V_u\left(x_{k+1}\right).$ (8)

对应最优控制 $u^{*}\left(x\right)$ 的最优Q函数为：

$Q^{*}\left(x_k,a\right)=Q_{u^{*}}\left(x_k,a\right)=P\left(x_k,a\right)+V^{*}\left(x_{k+1}\right).$ (9)

相应的最优控制为：

$u^{*}\left(x\right)=\arg \underset{u}{\min }{V}_u\left(x\right)=\arg \underset{a}{\min }{Q}^{*}\left(x,a\right).$ (10)

为了求解该最优控制问题，通过迭代的思想，希望利用Q函数的梯度信息来更新控制策略u，然后把控制策略u带入Q函数。策略梯度自适应动态规划算法的具体步骤如下：

步骤1：给定初始控制策略 $u^{\left(0\right)}\in U\left(X\right)$ 和允许误差 $\epsilon$ ，且令 $i=0$ 。

步骤2：估计Q函数：

$Q^{\left(i\right)}\left(x_k,a\right)=P\left(x_k,a\right)+Q^{\left(i\right)}\left(x_{k\text{+1}},u^{\left(i\right)}\right).$ (11)

步骤3：更新控制策略u：

$u^{\left(i+1\right)}\left(x\right)=u^{\left(i\right)}\left(x\right)-{\alpha {\nabla }_a{Q}^{\left(i\right)}\left(x,a\right)|}_{a=u^{\left(i\right)}}.$ (12)

其中 $\alpha$ 为常数。

步骤4：若 $\Vert Q^{\left(i\right)}\left(x_k,a\right)-Q^{\left(i+1\right)}\left(x_k,a\right)\Vert \le \epsilon$ ，则输出Q函数和控制策略，否则令 $i=i+1$ ，返回步骤2。

对于系统(1)，首先定义其Hamiltonian函数为：

$\begin{array}{c}H\left(x_k,u,V\right)=V\left(x_{k+1}\right)-V\left(x_k\right)+P\left(x_k,u\right)\\ =V\left(f\left(x_k,u\right)\right)-V\left(x_k\right)+P\left(x_k,u\right).\end{array}$ (13)

且 $V^{\left(i\right)}\left(x_k\right)=V_{u^{\left(i\right)}}\left(x_k\right)={\displaystyle {\sum }_{l=k}^{\infty }P\left(x_l,u^{\left(i\right)}\left(x_l\right)\right)}$ 。

引理1 [11] 给定 $u^{\left(0\right)}\left(x\right)\in U\left(X\right)$ ，根据策略梯度自适应动态规划算法得到控制序列 $\left\{u^{\left(i\right)}\left(x\right)\right\}$ 。假设 ${\nabla }_a{V}^{\left(i\right)}\left(f\left(x_k,a\right)\right)$ ， ${\nabla }_{aa}{V}^{\left(i\right)}\left(f\left(x_k,a\right)\right)$ ， ${\nabla }_{a}W\left(a\right)$ 和 ${\nabla }_{aa}W\left(a\right)$ 存在，且 $A={\nabla }_u^{\text{T}}H{\nabla }_{uu}H{\nabla }_{u}H$ ， $B={\nabla }_u^{\text{T}}H{\nabla }_{u}H$ ， $C=h$ 。 $\underset{\_}{\alpha }=\frac{B-\sqrt{B^2-AC}}{A}$ ， $\bar{\alpha }=\frac{B+\sqrt{B^2-AC}}{A}$ 。如果对于 $\forall i,x,a$ ， $B^2-AC\ge 0$ 且 $\underset{\_}{\alpha }\le \alpha \le \bar{\alpha }$ ，则：

1)，

2) $u^{\left(i\right)}\in U\left(X\right)$ 。

证明：1) 将Hamiltonian函数 $H\left(x_k,u,V^{\left(i\right)}\right)$ 在 $u^{\left(i\right)}$ 处进行二阶泰勒展开：

$\begin{array}{l}H\left(x_k,u,V^{\left(i\right)}\right)=H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right)+{\left(u-u^{\left(i\right)}\right)}^{\text{T}}{\nabla }_{u}H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right)\\ +\frac{1}{2}{\left(u-u^{\left(i\right)}\right)}^{\text{T}}{\nabla }_{uu}H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right)\left(u-u^{\left(i\right)}\right)+h\left(u\right).\end{array}$ (14)

根据式(8)和(13)易知：

${\nabla }_a{Q}^{\left(i\right)}\left(x_k,a\right)={\nabla }_a{V}^{\left(i\right)}\left(f\left(x_k,a\right)\right)+{\nabla }_{a}W\left(a\right)={\nabla }_{u}H\left(x_k,u,V^{\left(i\right)}\right).$ (15)

${\nabla }_{aa}{Q}^{\left(i\right)}\left(x_k,a\right)={\nabla }_{aa}{V}^{\left(i\right)}\left(f\left(x_k,a\right)\right)+{\nabla }_{aa}W\left(a\right)={\nabla }_{uu}H\left(x_k,u,V^{\left(i\right)}\right).$ (16)

由于 ${\nabla }_a{V}^{\left(i\right)}\left(f\left(x_k,a\right)\right),{\nabla }_{aa}{V}^{\left(i\right)}\left(f\left(x_k,a\right)\right),{\nabla }_{a}W\left(a\right),{\nabla }_{aa}W\left(a\right)$ ，对于 $\forall i$ 都存在，因此 ${\nabla }_{u}H,{\nabla }_{uu}H$ 都存在。由式(12)和(14)~(16)易知：

$\begin{array}{l}H\left(x_k,u^{\left(i+1\right)},V^{\left(i\right)}\right)\\ =H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right)+{\left(u^{\left(i+1\right)}-u^{\left(i\right)}\right)}^{\text{T}}{\nabla }_{u}H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right)\\ +\frac{1}{2}{\left(u^{\left(i+1\right)}-u^{\left(i\right)}\right)}^{\text{T}}{\nabla }_{uu}H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right)\left(u^{\left(i+1\right)}-u^{\left(i\right)}\right)+h\left(u^{\left(i+1\right)}\right)\end{array}$

$\begin{array}{l}=H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right)-\alpha {\nabla }_a^T{Q}^{\left(i\right)}\left(x_k,u^{\left(i\right)}\right){\nabla }_{u}H(x_k,u^{\left(i\right)},V^{\left(i\right)})\\ +\frac{1}{2}{\alpha }^2{\nabla }_a^{\text{T}}{Q}^{\left(i\right)}\left(x_k,u^{\left(i\right)}\right){\nabla }_{uu}H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right){\nabla }_a{Q}^{\left(i\right)}\left(x_k,u^{\left(i\right)}\right)+h\left(u^{\left(i+1\right)}\right)\\ =H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right)-\alpha {\nabla }_a^{\text{T}}H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right){\nabla }_{u}H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right)\\ +\frac{1}{2}{\alpha }^2{\nabla }_a^{T}H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right){\nabla }_{uu}H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right)H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right)+h\left(u^{\left(i+1\right)}\right).\end{array}$

令 $A={\nabla }_u^{\text{T}}H{\nabla }_{uu}H{\nabla }_{u}H,B={\nabla }_u^{\text{T}}H{\nabla }_{u}H,C=h$ ，且 $A\ge 0$ ，则：

$H\left(x_k,u^{\left(i+1\right)},V^{\left(i\right)}\right)=H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right)+\frac{1}{2}A{\alpha }^2-B\alpha +C.$ (17)

如果 $B^2-AC\ge 0$ ，当 $\underset{\_}{\alpha }\le \alpha \le \bar{\alpha }$ 时，其中 $\underset{\_}{\alpha }=\frac{B-\sqrt{B^2-AC}}{A}$ ， $\bar{\alpha }=\frac{B+\sqrt{B^2-AC}}{A}$ ，则：

$\frac{\text{1}}{\text{2}}A{\alpha }^2-B\alpha +C\le 0.$ (18)

即 $H\left(x_k,u^{\left(i+1\right)},V^{\left(i\right)}\right)\le H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right)$ 。根据Hamiltonian函数定义知： $H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right)=0$ 。

故： $H\left(x_k,u^{\left(i+1\right)},V^{\left(i\right)}\right)\le 0$ 。

2) 假设 $u^{\left(i\right)}\left(x\right)\in U\left(X\right)$ ，对于系统 $x_{k+1}=f\left(x_k,u^{\left(i+1\right)}\right)$ ，则对于Lyapunov函数 $V^{\left(i\right)}\left(x_k\right)$ 有：

$\begin{array}{c}\Delta V^{\left(i\right)}\left(x_k\right)=V^{\left(i\right)}\left(f\left(x_k,u^{\left(i+1\right)}\right)\right)-V^{\left(i\right)}\left(x_k\right)\\ =V^{\left(i\right)}\left(f\left(x_k,u^{\left(i+1\right)}\right)\right)-V^{\left(i\right)}\left(x_k\right)+P\left(x_k,u^{\left(i+1\right)}\right)-P\left(x_k,u^{\left(i+1\right)}\right)\\ =H\left(x_k,u^{\left(i+1\right)},V^{\left(i\right)}\right)-P\left(x_k,u^{\left(i+1\right)}\right).\end{array}$ (19)

由 $H\left(x_k,u^{\left(i+1\right)},V^{\left(i\right)}\right)\le 0$ ，知： $\Delta V^{\left(i\right)}\left(x_k\right)\le -P\left(x_k,u^{\left(i+1\right)}\right)\le 0$ 。即： $V^{\left(i\right)}\left(f\left(x_k,u^{\left(i+1\right)}\right)\right)\le V^{\left(i\right)}\left(x_k\right)<\infty$ 。

对于所有 $x,u\ne 0$ ，根据定义1知： $u^{\left(i+1\right)}\left(x\right)\in U\left(X\right)$ 。由数学归纳法知，当 $u^{\left(0\right)}\left(x\right)\in U\left(X\right)$ ，有 $u^{\left(i\right)}\left(x\right)\in U\left(X\right)$ 。

引理2 [11] 对于所有 $\left(x,a\right)\in D$ ，根据策略梯度自适应动态规划算法得到序列 $\left\{Q^{\left(i\right)}\left(x,a\right)\right\}$ 和 $$

满足：1) $Q^{\left(i\right)}\left(x,a\right)\ge Q^{\left(i+1\right)}\left(x,a\right)\ge Q^{*}\left(x,a\right)$ ，

2) $\underset{i\to \infty }{\lim }{Q}^{\left(i\right)}\left(x,a\right)=Q^{*}\left(x,a\right)$ 。

证明 1)由式(4)和(8)知：

$\begin{array}{c}{Q}^{\left(i\right)}\left(x_k,u^{\left(i+1\right)}\right)=P\left(x_k,u^{\left(i+1\right)}\right)+V^{\left(i\right)}\left(f\left(x_k,u^{\left(i+1\right)}\right)\right)\\ =P\left(x_k,u^{\left(i+1\right)}\right)+V^{\left(i\right)}\left(f\left(x_k,u^{\left(i+1\right)}\right)\right)-V^{\left(i\right)}\left(x_k\right)+V^{\left(i\right)}\left(x_k\right)\\ =H\left(x_k,u^{\left(i+1\right)},V^{\left(i\right)}\right)+V^{\left(i\right)}\left(x_k\right)\\ \le H\left(x_k,u^{\left(i\right)},V^{\left(i\right)}\right)+V^{\left(i\right)}\left(x_k\right)\\ =V^{\left(i\right)}\left(x_k\right).\end{array}$ (20)

且：

$\begin{array}{c}{V}^{\left(i+1\right)}\left(x_k\right)=P\left(x_k,u^{\left(i+1\right)}\right)+V^{\left(i\right)}\left(x_{k+1}\right)-V^{\left(i\right)}\left(x_{k+1}\right)+V^{\left(i+1\right)}(x_{k+1})\\ =Q^{\left(i\right)}\left(x_k,u^{\left(i+1\right)}\right)-V^{\left(i\right)}\left(x_{k+1}\right)+V^{\left(i+1\right)}\left(x_{k+1}\right)\\ \le V^{\left(i\right)}(x_k)-V^{\left(i\right)}\left(x_{k+1}\right)+V^{\left(i+1\right)}\left(x_{k+1}\right)\\ =P\left(x_k,u^{\left(i\right)}\right)+V^{\left(i+1\right)}\left(x_{k+1}\right)\\ \le P\left(x_k,u^{\left(i\right)}\right)+P\left(x_{k+1},u^{\left(i\right)}\right)+V^{\left(i+1\right)}\left(x_{k+2}\right)\\ \le {\displaystyle {\sum }_{l=k}^{\infty }P\left(x_l,u^{\left(i\right)}\right)}=V^{\left(i\right)}\left(x_k\right).\end{array}$ (21)

则对于所有 $\left(x_k,a\right)\in D$ ，有：

$\begin{array}{c}{Q}^{\left(i+1\right)}\left(x_k,a\right)=P\left(x_k,a\right)+V^{\left(i+1\right)}\left(x_{k+1}\right)\\ \le P\left(x_k,a\right)+V^{\left(i\right)}\left(x_{k+1}\right)\\ =Q^{\left(i\right)}\left(x_k,a\right).\end{array}$ (22)

由于式(9)知：

$\begin{array}{c}{Q}^{\left(i\right)}\left(x_k,a\right)=P\left(x_k,a\right)+V^{\left(i\right)}\left(x_{k+1}\right)\\ \ge P\left(x_k,a\right)+V^{*}\left(x_{k+1}\right)\\ =Q^{\text{*}}\left(x_k,a\right).\end{array}$ (23)

故可得： $Q^{\left(i\right)}\left(x,a\right)\ge Q^{\left(i+1\right)}\left(x,a\right)\ge Q^{*}\left(x,a\right)$ 。

2) 由式(12)知，当i趋于 $\infty$ 时有：

$u^{\infty }=u^{\infty }-{\alpha {\nabla }_a{Q}^{\infty }\left(x,a\right)|}_{a=u^{\infty }}.$ (24)

即：

${\alpha {\nabla }_a{Q}^{\infty }\left(x,a\right)|}_{a=u^{\infty }}=0.$ (25)

由式(8)知：

$$ (26)

则有：

${{\nabla }_{a}W\left(a\right)|}_{a=u^{\infty }}+{{\nabla }_a{V}^{\infty }\left(f\left(x_k,a\right)\right)|}_{a=u^{\infty }}=0.$ (27)

又由式(4)知：

(28)

上式对u求导得：

${{\nabla }_{u}W\left(u\right)|}_{u=u^{*}}+{{\nabla }_u{V}^{\text{*}}\left(f\left(x_k,u\right)\right)|}_{u=u^{*}}=0.$ (29)

易知式(27)和(29)是一样的，由唯一性得：

$V^{\infty }\left(x\right)=V^{\text{*}}\left(x\right).$ (30)

由式(26)得：

$\begin{array}{c}{Q}^{\infty }\left(x_k,a\right)=P\left(x_k,a\right)+V^{\infty }\left(x_{k+1}\right)\\ =P\left(x_k,a\right)+V^{*}\left(x_{k+1}\right)\\ =Q^{*}\left(x_k,a\right).\end{array}$ (31)

由结论1)知，序列是单减序列，记 $\underset{i\to \infty }{\lim }{Q}^{\left(i\right)}\left(x,a\right)=Q^{\infty }\left(x,a\right)$ ，则有：

$\underset{i\to \infty }{\lim }{Q}^{\left(i\right)}\left(x,a\right)=Q^{*}\left(x,a\right)$ 。

4. 基于数据驱动的神经动态规划及其实现

首先定义一个数据 $\left(x,a,\tilde{x}\right)$ 其中 $$ 为当前状态x执行控制动作a后得到状态向量。 $\tilde{x}$ 的获得是通过实际系统输入控制动作a后得到，而不是使用系统模型f的数学解析式。如图1，是策略梯度自适应动态规划算法的结构。其中包括两个部分：离线数据集和在线数据。离线数据集。M为数据的数量，其中离线数据集 $S_M$ 可以通过实际系统随机采样获得，其结构如图2。在线数据 $s_k=\left(x_{k-1},u_{k-1},x_k\right)$ ，其分别是在时刻 $k-1$ 和k的在线状态和控制信息，其获得的结构图如图3。图1的大概流程是：首先，给定一个初始控制 $u^{\left(0\right)}$ ，根据离线数据集 $S_M$ 通过策略梯度自适应动态规划算法步骤2，得到 $Q^{\left(0\right)}\left(x,a\right)$ 。其次，将初始控制 $u^{\left(0\right)}$ 通过实际系统作用于当前状态 $x_0$ ，得到下一状态 $x_1$ ，即得到在线数据 $s_1$ 。然后，将 $Q^{\left(0\right)}\left(x,a\right)$ 和 $u^{\left(0\right)}$ 通过策略梯度自适应动态规划算法步骤3，得到 $u^{\left(1\right)}$ 。最后，将 $s_1$ 加入离线数据集 $S_M$ 作为新的离线数据集，以此重复循环，则可得到在线数据 $s_2,s_3,\cdots$ ，控制序列 $u^{\left(2\right)},u^{\left(3\right)},\cdots$ 和Q函数序列 $Q^{\left(2\right)}\left(x,a\right),Q^{\left(3\right)}\left(x,a\right),\cdots$ 。

Figure 1. Policy gradient adaptive dynamic programming algorithm

图1. 策略梯度自适应动态规划算法

Figure 2. Collect offline data sets

图2. 收集离线数据集

Figure 3. Collect online data

图3. 收集在线数据

4.1. 神经动态规划的实现设计

由于根据策略梯度自适应动态规划算法，要求去计算未知的Q函数 $$ 和控制策略 $u^{\left(i\right)}\left(x\right)$ ，为了实现该过程，下面引入执行―评价神经网络。其中用执行网络来逼近控制策略 $u\left(x\right)$ ，用评价网络来逼近Q函数 $Q\left(x,a\right)$ 。首先介绍俩组线性无关的基函数： $\Phi \left(x\right)={\left\{{\varphi }_j\left(x\right)\right\}}_{j=1}^{\infty }$ ， $\Psi \left(x,a\right)={\left\{{\psi }_j\left(x,a\right)\right\}}_{j=1}^{\infty }$ 。其中 ${\varphi }_j\left(0\right)=0,{\psi }_j\left(0,0\right)=0$ 。根据代数理论，控制策略 $u\left(x\right)$ 和Q函数 $Q\left(x,a\right)$ 可以如下线性表示：

$u\left(x\right)={\displaystyle {\sum }_{j=1}^{\infty }{\upsilon }_j{\varphi }_j\left(x\right)}.$ (32)

$Q\left(x,a\right)={\displaystyle {\sum }_{j=1}^{\infty }{\theta }_j{\psi }_j\left(x,a\right)}.$ (33)

根据函数逼近理论，Q函数 $Q^{\left(i\right)}\left(x,a\right)$ 和控制策略 $u^{\left(i\right)}\left(x\right)$ 可以被有限维基函数近似表示：

$$ (34)

${\bar{Q}}^{\left(i\right)}\left(x,a\right)={\displaystyle {\sum }_{j=1}^{L_2}{\theta }_j^{\left(i\right)}{\psi }_j\left(x,a\right)}={\Psi }_L^{\text{T}}\left(x,a\right){\theta }^{\left(i\right)}.$ (35)

其中 ${\upsilon }^{\left(i\right)}={\left[{\upsilon }_1^{\left(i\right)}\cdots {\upsilon }_{L_1}^{\left(i\right)}\right]}^{\text{T}}$ 和 ${\theta }^{\left(i\right)}={\left[{\theta }_1^{\left(i\right)}\cdots {\theta }_{L_2}^{\left(i\right)}\right]}^{\text{T}}$ 分别为执行网络和评价网络的权重系数向量，但其是未知的。 ${\Phi }_L\left(x\right)={\left[{\varphi }_1\left(x\right)\cdots {\varphi }_{L_1}\left(x\right)\right]}^{\text{T}}$ 和 ${\Psi }_L\left(x,a\right)={\left[{\psi }_1\left(x,a\right)\cdots {\psi }_{L_2}\left(x,a\right)\right]}^{\text{T}}$ 分别为执行网络和评价网络的激活函数向量。则该神经网络的输出可以表示为：

${\hat{u}}^{\left(i\right)}\left(x\right)={\displaystyle {\sum }_{j=1}^{L_1}{\hat{\upsilon }}_j^{\left(i\right)}{\varphi }_j\left(x\right)}={\Phi }_L^{\text{T}}\left(x\right){\hat{\upsilon }}^{\left(i\right)}.$ (36)

${\hat{Q}}^{\left(i\right)}\left(x,a\right)={\displaystyle {\sum }_{j=1}^{L_2}{\hat{\theta }}_j^{\left(i\right)}{\psi }_j\left(x,a\right)}={\Psi }_L^{\text{T}}\left(x,a\right){\hat{\theta }}^{\left(i\right)}.$ (37)

其中 $$ 和 ${\hat{\theta }}^{\left(i\right)}={\left[{\hat{\theta }}_1^{\left(i\right)}\cdots {\hat{\theta }}_{L_2}^{\left(i\right)}\right]}^{\text{T}}$ 分别为 ${\upsilon }^{\left(i\right)}$ 和 ${\theta }^{\left(i\right)}$ 的近似估计。由于神经网络有误差，在用 ${\hat{u}}^{\left(i\right)}\left(x\right)$ 和 ${\hat{Q}}^{\left(i\right)}\left(x,a\right)$ 估计 $u^{\left(i\right)}\left(x\right)$ 和 $Q^{\left(i\right)}\left(x,a\right)$ 时会产生残差，定义Q函数残差为：

$\begin{array}{c}{\sigma }_Q^{\left(i\right)}\left(x,a,\tilde{x}\right)={\hat{Q}}^{\left(i\right)}\left(x,a\right)-{\hat{Q}}^{\left(i\right)}\left(\tilde{x},{\hat{u}}^{\left(i\right)}\right)-P\left(x,a\right)\\ ={\Psi }_L^{\text{T}}\left(x,a\right){\hat{\theta }}^{\left(i\right)}-{\Psi }_L^{\text{T}}\left(\tilde{x},{\hat{u}}^{\left(i\right)}\right){\hat{\theta }}^{\left(i\right)}-P\left(x,a\right)\\ ={\Psi }_L^{\text{T}}\left(x,a\right){\hat{\theta }}^{\left(i\right)}-{\Psi }_L^{\text{T}}\left(\tilde{x},{\Phi }_L^{\text{T}}\left(\tilde{x}\right){\hat{\upsilon }}^{\left(i\right)}\right){\hat{\theta }}^{\left(i\right)}-P\left(x,a\right).\end{array}$ (38)

下面的目标是，在满足残差趋于0的条件下，基于数据驱动来计算 ${\hat{\upsilon }}^{\left(i\right)}$ ， ${\hat{\theta }}^{\left(i\right)}$ 。任意 $s_1\left(x,u\right)\in U$ 和 $$ ，定义内积 ${\langle s_{\text{1}}\left(x,u\right),s_{\text{2}}\left(x,u\right)\rangle }_D={\displaystyle {\int }_D{s}_{\text{1}}\left(x,u\right),s_{\text{2}}\left(x,u\right)\text{d}\left(x,u\right)}$ 。令：

(39)

将(38)式带入(39)式得：

$\begin{array}{l}\left[{\langle {\Psi }_L\left(x,a\right),{\Psi }_L^{\text{T}}\left(x,a\right)\rangle }_D-{\langle {\Psi }_L\left(x,a\right),{\Psi }_L^{\text{T}}\left(\tilde{x},{\Phi }_L^{\text{T}}\left(\tilde{x}\right){\hat{\upsilon }}^{\left(i\right)}\right)\rangle }_D\right]{\hat{\theta }}^{\left(i\right)}\\ -{\langle {\Psi }_L\left(x,a\right),P\left(x,a\right)\rangle }_D=0.\end{array}$ (40)

则可得：

(41)

在计算 ${\hat{\theta }}^{\left(i\right)}$ 时，其中涉及许多积分，根据蒙特卡洛积分方法，令： $I_D={\displaystyle {\int }_D\text{d}\left(x,a\right)}$ 。首先基于离线数据集 $S_M$ 来计算 ${\hat{\theta }}^{\left(0\right)}$ ：

$\begin{array}{c}{\langle {\Psi }_L\left(x,a\right),{\Psi }_L^{\text{T}}\left(x,a\right)\rangle }_D={\displaystyle {\int }_D{\Psi }_L\left(x,a\right){\Psi }_L^{\text{T}}\left(x,a\right)\text{d}\left(x,a\right)}\\ =\frac{I_D}{M}{\displaystyle {\sum }_{l=1}^M{\Psi }_L\left(x_l,a_l\right){\Psi }_L^{\text{T}}\left(x_l,a_l\right)}\\ =\frac{I_D}{M}{\eta }_0.\end{array}$ (42)

其中 ${\eta }_0={\displaystyle {\sum }_{l=1}^M{\Psi }_L\left(x_l,a_l\right){\Psi }_L^{\text{T}}\left(x_l,a_l\right)}$ 。同理可得：

${\langle {\Psi }_L\left(x,a\right),{\Psi }_L^{\text{T}}\left(\tilde{x},{\Phi }_L^{\text{T}}\left(\tilde{x}\right){\hat{\upsilon }}^{\left(0\right)}\right)\rangle }_D=\frac{I_D}{M}{\rho }_0.$ (43)

$$ (44)

且 $$ ， ${\xi }_0={\displaystyle {\sum }_{l=1}^M{\Psi }_L\left(x_l,a_l\right)P\left(x_l,a_l\right)}$ 。则可得：

${\hat{\theta }}^{\left(0\right)}={\left({\eta }_0-{\rho }_0\right)}^{-1}{\xi }_0.$ (45)

如图4，基于离线数据集 $S_M$ 可以计算出 ${\hat{\theta }}^{\left(0\right)}$ ，当 $i=k$ 时，对于在线数据 $s_k$ ，此时把在线数据 $s_k$ 加入到离线数据集 $S_M$ 作为新的离线数据集 $S_M+s_k$ ，且用其来计算 ${\hat{\theta }}^{\left(k\right)}$ ：

${\langle {\Psi }_L\left(x,a\right),{\Psi }_L^{\text{T}}\left(x,a\right)\rangle }_D=\frac{I_D}{M+1}{\eta }_k.$ (46)

${\langle {\Psi }_L\left(x,a\right),{\Psi }_L^{\text{T}}\left(\tilde{x},{\Phi }_L^{\text{T}}\left(\tilde{x}\right){\hat{\upsilon }}^{\left(k\right)}\right)\rangle }_D=\frac{I_D}{M+1}{\rho }_k.$ (47)

${\langle {\Psi }_L\left(x,a\right),P\left(x,a\right)\rangle }_D=\frac{I_D}{M\text{+1}}{\xi }_k.$ (48)

Figure 4. Neural dynamic programming algorithm

图4. 神经动态规划算法

其中 ${\eta }_k={\eta }_0+{\Psi }_L\left(x_{k-1},a_{k-1}\right){\Psi }_L^{\text{T}}\left(x_{k-1},a_{k-1}\right)$ ，

${\rho }_k={\displaystyle {\sum }_{l=1}^M{\Psi }_L\left(x_l,a_l\right){\Psi }_L^{\text{T}}\left({\tilde{x}}_l,{\Phi }_L^{\text{T}}\left({\tilde{x}}_l\right){\hat{\upsilon }}^{\left(k\right)}\right)}+{\Psi }_L\left(x_{k-1},a_{k-1}\right){\Psi }_L^{\text{T}}\left(x_k,u_k\right).$

则由(41)式知，当 $i=k$ 时：

${\hat{\theta }}^{\left(k\right)}={\left({\eta }_k-{\rho }_k\right)}^{-1}{\xi }_k.$ (49)

接着需要计算执行网络的权重系数 ${\hat{\upsilon }}^{\left(i\right)}$ ，由神经动态规划算法步骤3更新控制策略时，用 ${\hat{u}}^{\left(i\right)}\left(x\right)$ 替换 $u^{\left(i\right)}\left(x\right)$ 会产生误差，定义控制策略残差为：

$\begin{array}{c}{\sigma }_u^{\left(i\right)}\left(x\right)=u^{\left(i+1\right)}\left(x\right)-u^{\left(i\right)}\left(x\right)+{\alpha {\nabla }_a{Q}^{\left(i\right)}\left(x,a\right)|}_{a=u^{\left(i\right)}}\\ ={\Phi }_L^{\text{T}}\left(x\right){\hat{\upsilon }}^{\left(i+1\right)}-{\Phi }_L^{\text{T}}\left(x\right){\hat{\upsilon }}^{\left(i\right)}+\alpha {\nabla }_a{\Phi }_L^{\text{T}}\left(x,{\Phi }_L^{\text{T}}\left(x\right){\hat{\upsilon }}^{\left(i\right)}\right){\hat{\theta }}^{\left(i\right)}.\end{array}$ (50)

同理在计算 ${\hat{\upsilon }}^{\left(i\right)}$ 时要满足控制策略残差趋于0：

$$ (51)

即：

$\begin{array}{l}{\langle {\Phi }_L\left(x\right),{\Phi }_L^{\text{T}}\left(x\right)\rangle }_X{\hat{\upsilon }}^{\left(i+1\right)}-{\langle {\Phi }_L\left(x\right),{\Phi }_L^{\text{T}}\left(x\right)\rangle }_X{\hat{\upsilon }}^{\left(i\right)}\\ +\alpha {\langle {\Phi }_L\left(x\right),{\nabla }_a{\Phi }_L^{\text{T}}\left(x,{\Phi }_L^{\text{T}}\left(x\right){\hat{\upsilon }}^{\left(i\right)}\right)\rangle }_X{\hat{\theta }}^{\left(i\right)}=0.\end{array}$ (52)

则可得：

${\hat{\upsilon }}^{\left(i+1\right)}={\hat{\upsilon }}^{\left(i\right)}-\alpha {\langle {\Phi }_L\left(x\right),{\Phi }_L^{\text{T}}\left(x\right)\rangle }_X^{-1}\times {\langle {\Phi }_L\left(x\right),{\nabla }_a{\Phi }_L^{\text{T}}\left(x,{\Phi }_L^{\text{T}}\left(x\right){\hat{\upsilon }}^{\left(i\right)}\right)\rangle }_X{\hat{\theta }}^{\left(i\right)}.$ (53)

根据蒙特卡洛积分方法，令： $I_X={\displaystyle {\int }_X\text{d}x}$ 。基于离线数据集 $$ ，当 $i=k$ 时，对于在线数据 $s_k$ ，此时把在线数据 $s_k$ 与离线数据集 $S_M$ 结合，且用其来计算 ${\hat{\upsilon }}^{\left(k\right)}$ ：

$$ (54)

${\langle {\Phi }_L\left(x\right),{\nabla }_a{\Psi }_L^{\text{T}}\left(x,{\Phi }_L^{\text{T}}\left(x\right){\hat{\upsilon }}^{\left(k\right)}\right)\rangle }_X=\frac{I_X}{M+1}{F}_k.$ (55)

其中 ${\Gamma }_k={\Gamma }_0+{\Phi }_L\left(x_k\right){\Phi }_L^{\text{T}}\left(x_k\right)$ ，且 ${\Gamma }_0={\displaystyle {\sum }_{l=1}^M{\Phi }_L\left(x_l\right){\Phi }_L^{\text{T}}\left(x_l\right)}$ ，

$F_k={\displaystyle {\sum }_{l=1}^M{\Phi }_L\left(x_l\right){\nabla }_a{\Psi }_L^{\text{T}}\left(x_l,{\Phi }_L^{\text{T}}\left(x_l\right){\hat{\upsilon }}^{\left(k\right)}\right)}+{\Phi }_L\left(x_k\right){\nabla }_a{\Psi }_L^{\text{T}}\left(x_k,{\Phi }_L^{\text{T}}\left(x_k\right){\hat{\upsilon }}^{\left(k\right)}\right)$ 。则由(53)可得：

${\hat{\upsilon }}^{\left(k+1\right)}={\hat{\upsilon }}^{\left(k\right)}-\alpha {\Gamma }_k^{-1}{F}_k{\hat{\theta }}^{\left(k\right)}.$ (56)

4.2. 神经动态规划算法

步骤1：收集离线数据集 $S_M$ ，计算 $$ 。

步骤2：给定一个初始允许误差 $\epsilon$ 和初始控制策略 ${\hat{u}}^{\left(i\right)}\left(x\right)\in U\left(X\right)$ ，并令 $i=0$ 。

步骤3： 使用 $S_M$ 和 $s_{i+1}$ 计算 ${\eta }_i,{\rho }_i,{\xi }_i$ 。 并根据(49)式计算 ${\hat{\theta }}^{\left(i\right)}$ 。

步骤4：使用离线数据集 $$ 和在线状态数据 $x_i$ ，计算 ${\Gamma }_i,F_i$ ，并根据(56)式计算 ${\hat{\upsilon }}^{\left(i+1\right)}$ 。

步骤5：若 $\Vert {\hat{\theta }}^{\left(i\right)}-{\hat{\theta }}^{\left(i+1\right)}\Vert \le \epsilon$ ，则输出 ${\hat{\theta }}^{\left(i+1\right)}$ 和 ${\hat{\upsilon }}^{\left(i+1\right)}$ ；否则令 $i=i+1$ ，返回步骤3，继续循环。

5. 结束语

本文提出了一种基于数据驱动的神经动态规划方法。该方法不依赖于系统的数学解析式，采用神经网络与动态规划结合的方式对最优控制问题进行求解。其分别利用Q函数的残差和Q函数的基函数做内积为零，控制策略的残差与控制策略的基函数做内积为零；并使用离线数据集与在线数据来迭代更新神经网络的系数，最后得到所需的控制策略。该方法能将离线数据与在线数据有效结合，使得系数更新更加完善。并且证明该算法是收敛的；且收敛到最优值。

基金项目

广东省自然科学基金项目(No.2018A030313505)，广东省科技计划项目(No.2017B010124003, No.2017 B090909001)。