精华内容
下载资源
问答
  • 利用ADP(自适应动态规划)中的HDP,实现对非线性离散时间系统的稳定控制。代码利用python实现。构建评价网络(Critic Network)、模型网络(Model Network)和执行网络(Action Network)三个网络。需要安装Pytorch;由于...
  • 自适应动态规划综述

    2017-11-09 13:40:56
    这是一篇介绍自适应动态规划的综述,里面详细介绍了离散和连续系统的动态规划方法,同时讲解了PI与VI,并给出具体的算法,还是很不错的
  • 自适应动态规划matalab简单代码实现,适合初学者,代码可运行
  • Adaptive Dynamic Programming 自适应动态规划的入门介绍。令初学者简明扼要的了解ADP的核心思想。
  • 本代码运用自适应动态规划理论,结合BP神经网络,设计实现多智能体系统的一致控制。其中的控制率是由HDP框架的BP神经网络基于智能体的实时状态数据自适应产生的。
  • 使用自适应动态规划实现单极倒立摆的控制,可供学习参考
  • 利用ADP中的HDP(with two Critic Network)求解离散非线性系统,代码利用python实现。构建评价网络(Critic Network)、模型网络(Model Network),Model Network直接用原系统代替。(需要安装Pytorch;...
  • 针对一类典型的带有控制约束的非线性离散时间系统, 提出了一种基于自适应动态规划(adaptive dynamic program,ADP) 算法的多设定值跟踪控制方法,并对其收敛性和稳定性做了严格分析.在ADP迭代跟踪控制的基础上, 根据多...
  • 为了求解有限时域最优控制问题, 自适应动态规划(ADP) 算法要求受控系统能一步控制到零. 针对不能一步控制到零的非线性系统, 提出一种改进的ADP 算法, 其初始代价函数由任意的有限时间容许序列构造. 推导了算法的迭代...
  • 基于自适应动态规划的多模型自适应跟踪控制
  • 自适应动态规划(一)

    千次阅读 多人点赞 2020-10-08 09:42:52
    自适应动态规划(一) 先立一个flag,这个算法我一定要研究透彻,连续更新。 动态规划 参考书籍《最优控制理论与系统》第四章 动态规划 递推方程 JN(x)=min⁡SN(x){d[x,sN(x)]+JN−1[SN(x)]}J1(x)=d(x,F) J_N(x)=\...

    自适应动态规划(一)

    先立一个flag,这个算法我一定要研究透彻,连续更新。

    动态规划

    参考书籍《最优控制理论与系统》第四章 动态规划

    递推方程
    J N ( x ) = min ⁡ S N ( x ) { d [ x , s N ( x ) ] + J N − 1 [ S N ( x ) ] } J 1 ( x ) = d ( x , F ) J_N(x)=\min_{S_N(x)}\{d[x,s_N(x)]+J_{N-1}[S_N(x)]\} \\ J_1(x)=d(x,F) JN(x)=SN(x)min{d[x,sN(x)]+JN1[SN(x)]}J1(x)=d(x,F)
    最优性原理

    贝尔曼指出,多级决策过程的最优策略具有这样的性质;不论初始状态和初始决策如何,当把其中的任何一级和状态再作为初始级和初始状态时,其余的决策对此必定也是一个最优策略。

    动态规划本质上告诉人们:整体最优闭为局部最优。

    动态规划要求按照时间阶段逆向计算,而同时动态系统的状态又要求根据系统函数按照时间正向顺序计算,这使得传统动态规划的直接应用变得困难。

    离散系统例题

    已知离散系统方程
    x ( k + 1 ) = 2 x ( k ) + u ( k ) ,         x ( 0 ) = 1 x(k+1)=2x(k)+u(k),~~~~~~~x(0)=1 x(k+1)=2x(k)+u(k),       x(0)=1
    代价函数
    J = x 2 ( 3 ) + ∑ k = 0 2 [ x 2 ( k ) + u 2 ( k ) ] J=x^2(3)+\sum_{k=0}^2[x^2(k)+u^2(k)] J=x2(3)+k=02[x2(k)+u2(k)]
    式中的状态x(k)及控制u(k)均不受约束。试求最优控制序列 u ∗ ( 0 ) , u ∗ ( 1 ) , u ∗ ( 2 ) {u^*(0),u^*(1),u^*(2)} u(0),u(1),u(2),使代价函数极小。

    解:

    N = 3 , k = 2 N=3,k=2 N=3,k=2
    J 1 ∗ [ x ( 2 ) ] = min ⁡ u ( 2 ) { [ x 2 ( 2 ) + u 2 ( 2 ) ] + J 0 ∗ [ x ( 3 ) ] } J 0 ∗ [ x ( 3 ) ] = x 2 ( 3 ) = [ 2 x ( 2 ) + u ( 2 ) ] 2 \begin{aligned} J_{1^*}[x(2)]&=\min_{u(2)}\{[x^2(2)+u^2(2)]+J_{0^*}[x(3)]\} \\ J_{0^*}[x(3)]&=x^2(3)=[2x(2)+u(2)]^2 \end{aligned} J1[x(2)]J0[x(3)]=u(2)min{[x2(2)+u2(2)]+J0[x(3)]}=x2(3)=[2x(2)+u(2)]2
    所以
    J 1 ∗ [ x ( 2 ) ] = min ⁡ u ( 2 ) { x 2 ( 2 ) + u 2 ( 2 ) + [ 2 x ( 2 ) + u ( 2 ) ] 2 } J_{1^*}[x(2)]=\min_{u(2)}\{x^2(2)+u^2(2)+[2x(2)+u(2)]^2\} J1[x(2)]=u(2)min{x2(2)+u2(2)+[2x(2)+u(2)]2}
    由于 u ( k ) u(k) u(k)无约束,故令
    ∂ { ⋅ } ∂ u ( 2 ) = 2 u ( 2 ) + 2 [ 2 x ( 2 ) + u ( 2 ) ] = 0 \frac{\partial \{\cdot\}}{\partial u(2)}=2u(2)+2[2x(2)+u(2)]=0 u(2){}=2u(2)+2[2x(2)+u(2)]=0
    求得
    u ∗ ( 2 ) = − x ( 2 ) u^*(2)=-x(2) u(2)=x(2)
    将上述结果带入 J 1 ∗ [ x ( 2 ) ] J_{1^*}[x(2)] J1[x(2)]方程,易得
    J 1 ∗ [ x ( 2 ) ] = 3 x 2 ( 2 ) J_{1^*}[x(2)]=3x^2(2) J1[x(2)]=3x2(2)
    N = 2 , k = 1 N=2,k=1 N=2,k=1
    J 2 ∗ [ x ( 1 ) ] = min ⁡ u ( 1 ) { [ x 2 ( 1 ) + u 2 ( 1 ) ] + J 1 ∗ [ x ( 2 ) ] } = min ⁡ u ( 1 ) { [ x 2 ( 1 ) + u 2 ( 1 ) ] + 3 [ 2 x ( 1 ) + u ( 1 ) ] 2 } \begin{aligned} J_{2^*}[x(1)]&=\min_{u(1)}\{[x^2(1)+u^2(1)]+J_{1^*}[x(2)]\} \\ &=\min_{u(1)}\{[x^2(1)+u^2(1)]+3[2x(1)+u(1)]^2\} \end{aligned} J2[x(1)]=u(1)min{[x2(1)+u2(1)]+J1[x(2)]}=u(1)min{[x2(1)+u2(1)]+3[2x(1)+u(1)]2}
    同理取偏导数为0解得
    u ∗ ( 1 ) = − 3 2 x ( 1 ) J 2 ∗ [ x ( 1 ) ] = 4 x 2 ( 1 ) u^*(1)=-\frac{3}{2}x(1) \\ J_{2^*}[x(1)]=4x^2(1) u(1)=23x(1)J2[x(1)]=4x2(1)
    N = 1 , k = 0 N=1,k=0 N=1,k=0
    J 3 ∗ [ x ( 0 ) ] = min ⁡ u ( 0 ) { [ x 2 ( 0 ) + u 2 ( 0 ) ] + J 1 ∗ [ x ( 1 ) ] } = min ⁡ u ( 0 ) { [ x 2 ( 0 ) + u 2 ( 0 ) ] + 4 [ 2 x ( 0 ) + u ( 0 ) ] 2 } \begin{aligned} J_{3^*}[x(0)]&=\min_{u(0)}\{[x^2(0)+u^2(0)]+J_{1^*}[x(1)]\} \\ &=\min_{u(0)}\{[x^2(0)+u^2(0)]+4[2x(0)+u(0)]^2\} \end{aligned} J3[x(0)]=u(0)min{[x2(0)+u2(0)]+J1[x(1)]}=u(0)min{[x2(0)+u2(0)]+4[2x(0)+u(0)]2}
    同理取偏导数为0解得
    u ∗ ( 0 ) = − 8 5 x ( 0 ) J 3 ∗ [ x ( 0 ) ] = 21 5 x 2 ( 0 ) u^*(0)=-\frac{8}{5}x(0) \\ J_{3^*}[x(0)]=\frac{21}{5}x^2(0) u(0)=58x(0)J3[x(0)]=521x2(0)
    带入已知的 x ( 0 ) = 1 x(0)=1 x(0)=1,按正向顺序算出
    u ∗ ( 0 ) = − 8 5 ,       x ∗ ( 1 ) = 2 ( 0 ) + u ∗ ( 0 ) = 2 5 ,       J 3 ∗ [ x ( 0 ) ] = 21 5 u ∗ ( 1 ) = − 3 5 ,       x ∗ ( 1 ) = 2 ( 1 ) + u ∗ ( 1 ) = 1 5 ,       J 2 ∗ [ x ( 1 ) ] = 16 25 u ∗ ( 2 ) = − 1 5 ,       x ∗ ( 3 ) = 2 ( 2 ) + u ∗ ( 2 ) = 1 5 ,       J 1 ∗ [ x ( 2 ) ] = 3 25 u^*(0)=-\frac{8}{5},~~~~~x^*(1)=2(0)+u^*(0)=\frac{2}{5},~~~~~J_{3^*}[x(0)]=\frac{21}{5} \\ u^*(1)=-\frac{3}{5},~~~~~x^*(1)=2(1)+u^*(1)=\frac{1}{5},~~~~~J_{2^*}[x(1)]=\frac{16}{25} \\ u^*(2)=-\frac{1}{5},~~~~~x^*(3)=2(2)+u^*(2)=\frac{1}{5},~~~~~J_{1^*}[x(2)]=\frac{3}{25} \\ u(0)=58,     x(1)=2(0)+u(0)=52,     J3[x(0)]=521u(1)=53,     x(1)=2(1)+u(1)=51,     J2[x(1)]=2516u(2)=51,     x(3)=2(2)+u(2)=51,     J1[x(2)]=253
    于是,最优轨迹及最优代价分别为
    u ∗ = − 8 5 , − 3 5 , − 1 5 x ∗ = 1 , 2 5 , 1 5 , 1 5 J ∗ = J 3 ∗ [ x ( 0 ) ] = 21 5 u^*={-\frac{8}{5},-\frac{3}{5},-\frac{1}{5}} \\ x^*={1,\frac{2}{5},\frac{1}{5},\frac{1}{5}} \\ J^*=J_3^*[x(0)]=\frac{21}{5} u=58,53,51x=1,52,51,51J=J3[x(0)]=521
    这种计算是比较复杂的,可以利用动态规划数值计算法来求解离散最优控制问题。但是计算量和存储量会急剧增长。

    维数灾难

    对于状态向量为 n n n维、控制向量为 m m m维、时间离散段为N的离散系统,在状态向量的每个元取 p p p个值、控制向量的每个元取 q q q个值的情况下,计算性能指标的求值次数为 N p n q m Np^nq^m Npnqm次,要求存储容量为 2 p n 2p^n 2pn个字。取 N = 10 , p = q = 20 , n = 3 , m = 1 N=10,p=q=20,n=3,m=1 N=10,p=q=20,n=3,m=1,则需要存储容量为1.6万个字,编制表格的时间约需16s。此时还是一个低维的问题,因此计算量限制了动态规划去解决高维问题,自适应动态规划可以解决这个问题。

    自适应动态规划

    自适应动态规划,利用前向计算更加符合系统的计算模型。

    ADP

    • HDP启发式动态规划
      • 最优控制要通过评判网络权值 ω \omega ω和输入输出关系式得出
    • DHP二次启发式动态规划
      • 最优控制可以通过协状态( ∂ J ( x ( k ) ) ∂ x ( k ) \frac{\partial J(x(k))}{\partial x(k)} x(k)J(x(k)))直接获得
    • ADHDP执行依赖启发式动态规划也被称为Q学习
      • 和HDP相似,仅仅是评判网络多了控制变量做为输入
    • ADDHP执行依赖二次启发式规划(同理)

    DHP相对于HDP具有更高的控制精度,但需要更高的计算量。

    启发式动态规划(HDP)

    HDP包含三个神经网络,分别为”执行网络“,”模型网络“,”评判网络“。

    • 执行网络:输入系统当前系统状态,输出控制变量。代表了控制策略。
    • 模型网络:用神经网络逼近实际的模型,利用神经网络的强大能力可以不必对实际的对象进行建模,降低算法部署的复杂性。
    • 评判网络:输入系统状态变量,输出性能指标。利用历史数据,用神经网络逼近代价函数,使系统可以通过这个网络进行前向计算。

    在HDP中,性能指标函数可以写成如下表达式
    J ( x ( k ) ) = l ( x ( k ) , u ( x ( k ) ) ) + J ( x ( k + 1 ) ) J(x(k))=l(x(k),u(x(k)))+J(x(k+1)) J(x(k))=l(x(k),u(x(k)))+J(x(k+1))
    其实评判网络的作用就是去拟合这个等式

    评判网络的loss一般都是平方误差,可以写为(用 J ω J_{\omega} Jω代表评价网络)
    l o s s = 1 2 [ J ω ( x ( k ) ) − l ( x ( k ) , u ( x ( k ) ) ) + J ( x ( k + 1 ) ) ] 2 loss = \frac{1}{2}[J_{\omega}(x(k))-l(x(k),u(x(k)))+J(x(k+1))]^2 loss=21[Jω(x(k))l(x(k),u(x(k)))+J(x(k+1))]2
    然后利用梯度下降的方法对网络进行迭代更新。用数学表示为:
    ω ∗ = a r g min ⁡ ω { ∣ J ( x ( k ) , ω ) − d ( x ( k ) , ω ) ∣ 2 } \omega^*=arg \min_{\omega}\{|J(x(k),\omega)-d(x(k),\omega)|^2\} ω=argωmin{J(x(k),ω)d(x(k),ω)2}
    其中
    d ( x ( k ) , ω ) = l ( x ( k ) , u ( x ( k ) ) ) + J ( x ( k + 1 ) , ω ) d(x(k),\omega)=l(x(k),u(x(k)))+J(x(k+1),\omega) d(x(k),ω)=l(x(k),u(x(k)))+J(x(k+1),ω)
    对控制网络的loss
    ∂ J ∗ ( x ( k ) ) ∂ u ( k ) = ∂ l ( x ( ( k ) , u ( k ) ) ) ∂ u ( k ) + ∂ J ∗ ( x ( k + 1 ) ) ∂ u ( k ) = ∂ l ( x ( ( k ) , u ( k ) ) ) ∂ u ( k ) + ∂ J ∗ ( x ( k + 1 ) ) ∂ x ( k + 1 ) ∂ x ( k + 1 ) ∂ u ( k ) \begin{aligned} \frac{\partial J^*(x(k))}{\partial u(k)}&=\frac{\partial l(x((k),u(k)))}{\partial u(k)}+\frac{\partial J^*(x(k+1))}{\partial u(k)} \\ &=\frac{\partial l(x((k),u(k)))}{\partial u(k)}+\frac{\partial J^*(x(k+1))}{\partial x(k+1)}\frac{\partial x(k+1)}{\partial u(k)} \end{aligned} u(k)J(x(k))=u(k)l(x((k),u(k)))+u(k)J(x(k+1))=u(k)l(x((k),u(k)))+x(k+1)J(x(k+1))u(k)x(k+1)

    ω ∗ = arg ⁡ min ⁡ ω { l ( x ( k ) , u ( x ( k ) , ω ) ) + J ∗ ( x ( k + 1 ) ) } \omega^*=\arg \min_{\omega}\{l(x(k),u(x(k),\omega))+J^*(x(k+1))\} ω=argωmin{l(x(k),u(x(k),ω))+J(x(k+1))}

    loss就可以表示为
    l o s s = 1 2 [ u ( x ( k ) , ω ) − arg ⁡ min ⁡ u { l ( x ( k ) , u ( x ( k ) ) ) + J ∗ ( x ( k + 1 ) ) } ] 2 loss=\frac{1}{2}[u(x(k),\omega)-\arg \min_{u}\{l(x(k),u(x(k)))+J*(x(k+1))\}]^2 loss=21[u(x(k),ω)argumin{l(x(k),u(x(k)))+J(x(k+1))}]2
    从上面可以看到,HDP主要是利用神经网络模拟了评价函数和控制函数,通过采样产生大量的数据来提前拟合原本需要实时在线进行反向递推的代价函数。

    在这里插入图片描述
    在这里插入图片描述

    在这里插入图片描述在这里插入图片描述

    HDP的参考代码

    代码是复制粘贴别人的ADP(自适应动态规划)-HDP,在此进行学习。

    import torch
    from torch.autograd import Variable
    import matplotlib.pyplot as plt
    import numpy as np
    
    torch.manual_seed(1)                                                           # 设置随机种子,使得每次生成的随机数是确定的
    state_dim = 2                                                                  # 状态维度
    v_dim = 1                                                                      # 价值维度
    action_dim = 1                                                                 # 动作维度
    learing_rate = 0.005                                                           # 学习率
    learing_num = 500                                                              # 学习次数
    sim_num = 20                                                                   # 仿真步长
    x0 = np.array([2,-1])                                                          # 初始状态
    epislon = 0.0001                                                               # 阈值
    Fre_V1_Paras = 5                                                               # 更新V1的频率
    
    ########################################################################################################################
    # 定义神经网络类
    ########################################################################################################################
    import torch
    from torch.autograd import Variable
    import matplotlib.pyplot as plt
    import numpy as np
    
    state_dim = 2                                                                  # 状态维度
    v_dim = 1                                                                      # 价值维度
    action_dim = 1                                                                 # 动作维度
    A_learing_rate = 0.01                                                          # A网络学习率
    V_learing_rate = 0.01                                                          # V网络学习率
    learing_num = 1000                                                              # 学习次数
    sim_num = 20                                                                   # 仿真步长
    x0 = np.array([2,-1])                                                          # 初始状态
    epislon = 0.01                                                                 # 阈值
    model_net_train_num = 50                                                       # 模型网络训练测次数
    
    torch.manual_seed(1)                                                           # 设置随机种子,使得每次生成的随机数是确定的
    
    # 采样状态  将状态定义在x1 [-2,2]   x2 [-1,1]
    x = np.arange(-2, 2, 0.1)
    y = np.arange(-1, 1, 0.1)
    xx, yy = np.meshgrid(x, y)  # 为一维的矩阵
    state = np.transpose(np.array([xx.ravel(), yy.ravel()]))  # 所有状态
    state_num = state.shape[0]  # 状态个数
    
    ####################################################################################################################
    # 定义原始模型函数
    ####################################################################################################################
    def model(current_state, u):
        next_state = np.zeros([current_state.shape[0], current_state.shape[1]])  # 初始化下一个状态
        for index in range(current_state.shape[0]):  # 对每个样本计算下一个状态 根据输入的u
            next_state[index, 0] = 0.2 * current_state[index, 0] * np.exp(current_state[index, 1] ** 2)
            next_state[index, 1] = 0.3 * current_state[index, 1] ** 3 - 0.2 * u[index]
            pass
        return next_state
    
    # 动作采样  将输入定在[-10 10] 内
    action = np.arange(-10, 10, 0.05)
    action_num = action.shape[0]
    
    # 计算状态-动作对
    a = np.random.rand(4,2)
    b = np.random.rand(3)
    state1 = np.repeat(state,action.shape[0],axis=0)
    action1 = np.tile(action,(state.shape[0],1)).reshape((-1,1))
    model_input = np.zeros([state.shape[0]*action.shape[0],2+1])
    model_input[:,0:2] = state1
    model_input[:,2] = action1.reshape(state.shape[0]*action.shape[0])
    
    # 计算模型网络的label
    model_target =  model(state1,action1)
    
    
    ########################################################################################################################
    # 定义神经网络类
    ########################################################################################################################
    class Net1(torch.nn.Module):
        # 初始化
        def __init__(self):
            super(Net1, self).__init__()
            self.lay1 = torch.nn.Linear(state_dim, 10, bias = False)               # 线性层
            self.lay1.weight.data.normal_(0,0.5)                                   # 权重初始化
            self.lay2 = torch.nn.Linear(10, action_dim, bias = False)                       # 线性层
            self.lay2.weight.data.normal_(0, 0.5)                                  # 权重初始化
    
        def forward(self, x):
            layer1 = self.lay1(x)                                                  # 第一隐层
            layer1 = torch.nn.functional.relu(layer1,inplace= False)               # relu激活函数
            output = self.lay2(layer1)                                             # 输出层
            return output
    
    class Net2(torch.nn.Module):
        # 初始化
        def __init__(self):
            super(Net2, self).__init__()
            self.lay1 = torch.nn.Linear(state_dim+action_dim, 10, bias = False)    # 线性层
            self.lay1.weight.data.normal_(0,0.5)                                   # 权重初始化
            self.lay2 = torch.nn.Linear(10, state_dim, bias = False)               # 线性层
            self.lay2.weight.data.normal_(0, 0.5)                                  # 权重初始化
    
        def forward(self, x):
            layer1 = self.lay1(x)                                                  # 第一隐层
            layer1 = torch.nn.functional.relu(layer1,inplace = False)              # relu激活函数
            output = self.lay2(layer1)                                             # 输出层
            return output
    
    ########################################################################################################################
    # 定义价值迭代类
    ########################################################################################################################
    class HDP():
        def __init__(self):
            self.V_model = Net1()                                                  # 定义V1网络
            self.A_model = Net1()                                                  # 定义A网络
            self.modelnet = Net2()                                                 # 定义模型网络
    
            self.criterion = torch.nn.MSELoss(reduction='mean')                    # 平方误差损失
            # 训练一定次数,更新Critic Net的参数
            # 这里只需要定义A网络和V2网络的优化器
            self.optimizerV = torch.optim.SGD(filter(lambda p: p.requires_grad, self.V_model.parameters()), lr = V_learing_rate)    # 利用梯度下降算法优化model.parameters
            self.optimizerA = torch.optim.SGD(filter(lambda p: p.requires_grad, self.A_model.parameters()), lr = A_learing_rate)    # 利用梯度下降算法优化model.parameters
            self.optimizer_model = torch.optim.SGD(self.modelnet.parameters(), lr=0.01)      # 利用梯度下降算法优化model.parameters
    
            self.learn_model()                                                               # 训练模型网络
    
            self.cost = []                                                                   # 初始化误差矩阵
            # for p in self.V_model.parameters():
            #     p.requires_grad = True  # 固定模型网络参数
            #     pass
            # print('Model Net Training has finished!')
            # for p in self.A_model.parameters():
            #     p.requires_grad = True  # 固定模型网络参数
            #     pass
            # print('Model Net Training has finished!')
            # pass
    
        def learn_model(self):
                print('Model Net Training has begun!')
                self.model_loss = []
                for learn_index in range(model_net_train_num):                               # 训练模型
                    model_predict = self.modelnet(Variable(torch.Tensor(model_input)))       # 预测值
                    loss = self.criterion(model_predict, Variable(torch.Tensor(model_target)))  # 计算损失
                    self.model_loss.append(loss)
                    #print('loss:   ', np.array(loss.data))
                    if np.abs(loss.data) < 0.001:
                        break
                    self.optimizer_model.zero_grad()                                         # 对模型参数做一个优化,并且将梯度清0
                    loss.backward(retain_graph=True)                                         # 计算梯度
                    self.optimizer_model.step()                                              # 权重更新
                    pass
                print('Model Net Training has finished!')
    
                pass
    
        ####################################################################################################################
        # J_loss函数
        ####################################################################################################################
        def J_loss1(self,sk,uk):
            Vk = np.zeros(sk.shape[0])                                                       # x_k 的V值
            for index in range(uk.shape[0]):                                                 # 对每个样本计算下一个状态 根据输入的u
                Vk[index] = sk[index,0] ** 2 + sk[index,1] ** 2 + uk[index] ** 2
                pass
            return Vk
            pass
    
        def J_loss2(self,sk,uk,Vk_1):
            Vk = np.zeros(uk.shape[0])                                                      # x_k 的V值
            for index in range(uk.shape[0]):                                                # 对每个样本计算下一个状态 根据输入的u
                Vk[index] = sk[index,0] ** 2 + sk[index,1] ** 2 + uk[index] ** 2 + Vk_1[index,0]
                pass
            return Vk
            pass
    
        def fitparas(self,p1,is_required):
            if is_required==1:
                for p in p1.parameters():
                  p.requires_grad = True  # 固定模型网络参数
                  pass
            else:
                for p in p1.parameters():
                    p.requires_grad = False  # 固定模型网络参数
                    pass
            pass
    
        ####################################################################################################################
        # 定义学习函数
        ####################################################################################################################
        def learning(self):
           self.loss = []
           for train_index in range(learing_num):
               print('the ' , train_index + 1 , ' --th  learing start')
    
               u_star = np.zeros([state_num, 1])
               for index in range(state_num):                                              # 循环计算所有状态的U*
                   st = model_input[(index)*action_num:(index+1)*action_num,:]
                   #print(st)
                   new_xk_1 = self.modelnet(Variable(torch.Tensor(st)))
                   next_V1 = self.V_model(Variable(torch.Tensor(new_xk_1)))
    
                   A1 = self.J_loss2(st, action, np.array(next_V1.data))
                   u_star_index = np.argmin(A1)
                   u_star[index] = action[u_star_index]
                   pass
    
               # 计算target
               Vk = self.V_model(Variable(torch.Tensor(state)))                            # 计算Vk
               Gd = self.J_loss1(state , u_star)                                           # 计算gD
               target = np.array(Vk.data).reshape(Vk.shape[0],1)-np.array(Gd).reshape(Vk.shape[0],1)                                 #计算标签
    
               # 计算预测值  # 这里需要对两个网络A和V网络计算loss 两个loss是一样的
               # 但是因为A网络需要经过Model网络计算梯度,不能直接中间结果保存到某个变量中,再带入另一个net 这样会出现 梯度为 None
               # 具体的操作如 predict1
               # 开始有用predict1 直接对V网络计算梯度,但是好像报了一个关于Pytorch版本的错误 然后 就利用中间变量保存
               #model_input_with_u = Variable(torch.Tensor(np.zeros([state_num,state_dim+action_dim])))
               model_input_with_u = torch.randn(state_num, state_dim + action_dim)
               #model_input_with_u[:,0:state_dim] = Variable(torch.Tensor(state))
               model_input_with_u[:,0:state_dim] = torch.tensor(state).type(torch.FloatTensor)
               model_input_with_u[:,state_dim] = self.A_model(Variable(torch.Tensor(state))).view(state_num).type(torch.FloatTensor)
               next_xk_1 = self.modelnet(Variable(torch.Tensor(model_input_with_u))).type(torch.FloatTensor)
               predict2 = self.V_model(Variable(torch.Tensor(next_xk_1)))                   # 计算Vk+1
    
               predict1 = self.V_model(\
                   self.modelnet(\
                       torch.cat((Variable(torch.tensor(state)).type(torch.FloatTensor),\
                                  self.A_model(Variable(torch.Tensor(state))).view(state_num,1).type(torch.FloatTensor)\
                                  ),1)
                   )
               )
    
               # predict2 = self.V_model( \
               #     self.modelnet( \
               #         torch.cat((Variable(torch.tensor(state)).type(torch.FloatTensor), \
               #                    self.A_model(Variable(torch.Tensor(state))).view(state_num, 1).type(torch.FloatTensor) \
               #                    ), 1)
               #     )
               # )
               # for tt in range(next_xk_1.shape[0]):
               #     for ttt in range(next_xk_1.shape[1]):
               #         next_xk_1[tt,ttt].requires_grad = True
               #         pass
               #     pass
               #
               # for t in range(predict.shape[0]):
               #     predict[t].requires_grad = True
    
               #print('是否求梯度', next_xk_1[5,1].requires_grad)
    
               #############################################################################################################
               # 更新Crictic Actor网络
               #############################################################################################################
    
               # for p in self.A_model.named_parameters():
               #     p.requires_grad
               #self.fitparas(self.A_model, 1)
               #self.fitparas(self.V_model, 0)
               #self.fitparas(self.modelnet, 0)                                             # Medel Net 不更新
               # A_loss = self.criterion(self.A_model(Variable(torch.Tensor(state))).view(state_num,1),Variable(torch.Tensor(target)))             # 计算损失
    
               A_loss = self.criterion(predict1,torch.tensor(target).type(torch.FloatTensor))             # 计算损失
               #A_loss.requires_grad = True
               self.optimizerA.zero_grad()                                                 # 对模型参数做一个优化,并且将梯度清0
               A_loss.backward(retain_graph=True)                                          # 计算梯度
               self.optimizerA.step()                                                      # 权重更新
               print('        the ', train_index + 1, ' Action Net have updated')
    
               # for name, parms in self.A_model.named_parameters():
               #     print('A_model-->name:', name, '-->grad_requirs:', parms.requires_grad, \
               #           ' -->grad_value:', parms.grad)
    
    
               # self.fitparas(self.V_model, 1)                                              # 只对V求梯度
               # self.fitparas(self.A_model, 0)
               #self.fitparas(self.modelnet, 0)
               V_loss = self.criterion(predict2, Variable(torch.Tensor(target)))            # 计算损失
               self.optimizerV.zero_grad()  # 对模型参数做一个优化,并且将梯度清0
               V_loss.backward(retain_graph=True)                                          # 计算梯度
               self.optimizerV.step()                                                      # 权重更新
               print('        the ' , train_index+1 , ' Critic Net have updated')
    
               # for name, parms in self.V_model.named_parameters():
               #     print('V_model-->name:', name, '-->grad_requirs:', parms.requires_grad, \
               #           ' -->grad_value:', parms.grad)
    
               # print('A paras:\n', list(self.A_model.named_parameters()))
               # print('V paras:\n', list(self.V_model.named_parameters()))
               # print('model paras:\n', list(self.modelnet.named_parameters()))
    
               print("        AC net loss: ",A_loss.data)
               self.loss.append(np.array(A_loss.data))
    
           # 保存和加载整个模型
           # 每次训练完可以保存模型,仿真时候 直接load训练好的模型 或者 继续训练可以接着上一次训练的结果继续训练
           #torch.save(model_object, 'model.pth')
           #model = torch.load('model.pth')
           pass
        #######################################################################################################
        # 定义仿真函数
        # 通过得到的Actor选择动作
        # 同时利用Critic计算V
        #######################################################################################################
        def simulator(self):
            print('the simulation is start')
            #self.restore(self.path)
            State_traject = np.zeros([sim_num + 1, state_dim])
            State_traject[0, :] = x0
            u_traject = np.zeros([sim_num, 1])
            for index in range(sim_num):
                print('the ', index, ' --th  time start')
                print('当前状态:', Variable(torch.Tensor(State_traject[index,:])).data)
                sim_actor = self.A_model(Variable(torch.Tensor(State_traject[index,:])))
                print('当前输入:',sim_actor)
                u_traject[index] = sim_actor.data
                sim_nexstate = model(State_traject[index, :].reshape(1, 2), sim_actor.data)  #带入实际系统中
                print('下一时刻状态:', sim_nexstate)
                State_traject[index + 1, :] = sim_nexstate
                pass
            pass
            V_traject = self.V_model(Variable(torch.Tensor(State_traject))).data
            print('the simulation is over')
            self.plot_curve(State_traject , u_traject , V_traject , self.loss,self.model_loss)
            pass
    
        #######################################################################################################
        # 绘图函数
        # 分别绘制状态轨迹 控制输入u轨迹 值函数V轨迹
        # 并将结果保存!
        #######################################################################################################
        def plot_curve(self, s, u, V,cost,modelloss):
            # print('\nstate\n',s)
            # print('\nu\n', u)
            # print('\nV\n', V)
            # 绘制状态轨迹
            plt.figure(1)
            plt.plot(s[:, 0], 'r', label='State_1')
            plt.plot(s[:, 1], 'b', label='State_2')
            plt.title('State_Trajecteory')
            plt.xlabel('iter')
            plt.ylabel('State')
            plt.legend()
            plt.grid()
            plt.savefig(r'ADPresultfig\HDP_state.png')
            plt.show()
    
            # 绘制控制输入u轨迹
            plt.figure(2)
            plt.plot(u, )
            plt.title('U_Trajecteory')
            plt.xlabel('iter')
            plt.ylabel('u')
            plt.grid()
            plt.savefig(r'ADPresultfig\HDP_u.png')
            plt.show()
    
            # 绘制值函数V的轨迹
            plt.figure(3)
            plt.plot(V, 'r')
            plt.title('V_Trajecteory')
            plt.xlabel('iter')
            plt.ylabel('V')
            plt.grid()
            plt.savefig(r'ADPresultfig\HDP_V.png')
            plt.show()
    
            print(cost)
            # 绘制值函数V的轨迹
            plt.figure(4)
            plt.plot(cost, 'r')
            plt.title('Train_loss_Trajecteory')
            plt.xlabel('iter')
            plt.ylabel('Train_loss')
            plt.grid()
            plt.savefig(r'ADPresultfig\HDP_loss.png')
            plt.show()
    
            # 绘制模型网络loss的轨迹
            plt.figure(5)
            plt.plot(modelloss, 'r')
            plt.title('Model_loss_Trajecteory')
            plt.xlabel('iter')
            plt.ylabel('Model_loss')
            plt.grid()
            plt.savefig(r'ADPresultfig\HDP_Model_loss.png')
            plt.show()
            pass
    
    ########################################################################################################################
    # 函数起始运行
    # 在仿真时候,直接调用最优的模型进行仿真
    # 最优的模型根据损失函数进行判断
    ########################################################################################################################
    if __name__ == '__main__':
        Agent = HDP()                                                            # 值迭代类实例化
        Agent.learning()                                                         # 学习
        Agent.simulator()                                                        # 仿真
    

    参考文献


    1. Discrete-time nonlinear HJB solution using Approximate dynamic programming: Convergence Proof
    2. ADP(自适应动态规划)-HDP
    3. 基于近似动态规划的非线性系统最优控制研究
    展开全文
  • (自适应动态规划综述)

    千次阅读 2020-01-08 09:36:49
    (自适应动态规划综述) 摘要:自适应动态规划(Adaptive/Approximate Dynamic Programming,ADP)是最优控制领域新兴起的一种近似最优方法,它在人工智能领域、强化学习、人工神经网络、模糊系统、演化计算等方面蓬勃...

    (自适应动态规划综述)
    摘要:自适应动态规划(Adaptive/Approximate Dynamic Programming,ADP)是最优控制领域新兴起的一种近似最优方法,它在人工智能领域、强化学习、人工神经网络、模糊系统、演化计算等方面蓬勃发展,为求解非线性系统优化问题提供了很多解决思路和具体技术方法,是当前国际最优化领域的研究热点。本文将按照自适应动态规划的研究背景意义、国内外研究现状、理论发展及应用四个方面对其进行介绍及总结。
    关键词:自适应动态规划,非线性系统,稳定性
    1 ADP的研究背景及意义
    动态运动、动态系统在自然界中普遍存在,对我们来说,要认识、理解并改善一个系统,对动态系统稳定性的研究必不可少。本世纪50∼60年代, 在空间技术发展和数字计算机实用化的推动下, 动态系统的优化理论得到了迅速的发展,形成了一个重要的学科分支:最优控制。它在空间技术、系统工程、经济管理与决策、人口控制、多级工艺设备的优化等许多领域都有越来越广泛的应用。1
    20世纪50年代初美国数学家Bellman等人在研究多阶段决策过程(multistep decision process)的优化问题时,提出了著名的最优化原理(principle of optimality),即把多阶段过程转化为一系列单阶段问题,利用各阶段之间的关系,逐个求解,创立了解决这类过程优化问题的新方法——动态规划(Dynamic Programming,DP)。动态规划,从本质上讲是一种非线性规划方法,其核心是贝尔曼最优性原理。这个原理可以归结为一个基本递推公式,从而使决策过程连续递推,并将一个多步(

    展开全文
  • 通过自适应动态规划求解相应的Hamilton-Jacobi-Bellman(HJB)方程,该方程由最小二乘技术,神经网络逼近器和策略迭代(PI)算法组成。 我们方法的主要思想是通过最小二乘技术对状态,状态导数和输入信息进行采样以...
  • 自适应动态规划(三) 值迭代稳定性证明 自适应动态规划的核心就是去求解除下面的序列,但是这个序列一定是收敛的吗?论文中给出了证明。 V(x(k))=min⁡u(k){U(x(k),u(k))+V(x(k+1))} V(x(k))=\min_{u(k)}\{U(x(k),u...

    自适应动态规划(三)

    值迭代稳定性证明

    自适应动态规划的核心就是去求解除下面的序列,但是这个序列一定是收敛的吗?论文中给出了证明。
    V ( x ( k ) ) = min ⁡ u ( k ) { U ( x ( k ) , u ( k ) ) + V ( x ( k + 1 ) ) } V(x(k))=\min_{u(k)}\{U(x(k),u(k))+V(x(k+1))\} V(x(k))=u(k)min{U(x(k),u(k))+V(x(k+1))}
    在这个证明中,首先确定的是 V 0 ( x ) = 0 V_0(x)=0 V0(x)=0 的初始条件,HDP的迭代公式如下:
    u 0 ( x ( k ) ) = arg ⁡ min ⁡ u ( x T ( k ) Q x ( k ) + u T R u + V 0 ( x ( k + 1 ) ) u_0(x(k))=\arg \min_{u}(x^T(k)Qx(k)+u^TRu+V_0(x(k+1)) u0(x(k))=argumin(xT(k)Qx(k)+uTRu+V0(x(k+1))
    然后再求代价
    V 1 = x T ( k ) Q x ( k ) + u 0 T ( x ( k ) ) R u 0 ( x ( k ) ) + V 0 ( x ( k + 1 ) ) V_1=x^T(k)Qx(k)+u_0^T(x(k))Ru_0(x(k))+V_0(x(k+1)) V1=xT(k)Qx(k)+u0T(x(k))Ru0(x(k))+V0(x(k+1))
    迭代过程
    u i ( x ( k ) ) = arg ⁡ min ⁡ u ( x T ( k ) Q x ( k ) + u T R u + V i ( x ( k + 1 ) ) V i + 1 = min ⁡ ( x T ( k ) Q x ( k ) + u T R u + V i ( x ( k + 1 ) ) = x T ( k ) Q x ( k ) + u i T ( x ( k ) ) R u i ( x ( k ) ) + V i ( x ( k + 1 ) ) \begin{aligned} u_i(x(k))&=\arg \min_{u}(x^T(k)Qx(k)+u^TRu+V_i(x(k+1)) \\ V_{i+1}&=\min(x^T(k)Qx(k)+u^TRu+V_i(x(k+1)) \\ &=x^T(k)Qx(k)+u_i^T(x(k))Ru_i(x(k))+V_i(x(k+1)) \end{aligned} ui(x(k))Vi+1=argumin(xT(k)Qx(k)+uTRu+Vi(x(k+1))=min(xT(k)Qx(k)+uTRu+Vi(x(k+1))=xT(k)Qx(k)+uiT(x(k))Rui(x(k))+Vi(x(k+1))
    定义
    Λ i + 1 ( x ( k ) ) = x ( k ) Q x ( k ) + u T ( i ) R u ( i ) + Λ i ( x ( k + 1 ) ) \Lambda_{i+1}(x(k))=x(k)Qx(k)+u^T(i)Ru(i)+\Lambda_i(x(k+1)) Λi+1(x(k))=x(k)Qx(k)+uT(i)Ru(i)+Λi(x(k+1))
    V 0 = Λ 0 = 0 V_0=\Lambda_0=0 V0=Λ0=0,则 V i ≤ Λ i V_i\leq\Lambda_i ViΛi。因为 V i + 1 V_{i+1} Vi+1每次取得都是最小值,因此显然上面的结论成立。

    假设 η ( x ( k ) ) \eta (x(k)) η(x(k)) 是稳定的控制, V 0 = Z 0 = 0 V_0=Z_0=0 V0=Z0=0
    Z i + 1 ( x ( k ) ) = x ( k ) Q x ( k ) + η T ( x ( k ) ) R η ( x ( k ) ) + Z i ( x ( k + 1 ) ) Z i ( x ( k ) ) = x ( k ) Q x ( k ) + η T ( x ( k ) ) R η ( x ( k ) ) + Z i − 1 ( x ( k + 1 ) ) Z_{i+1}(x(k))=x(k)Qx(k)+\eta^T (x(k))R\eta (x(k))+Z_i(x(k+1)) \\ Z_{i}(x(k))=x(k)Qx(k)+\eta^T (x(k))R\eta (x(k))+Z_{i-1}(x(k+1)) Zi+1(x(k))=x(k)Qx(k)+ηT(x(k))Rη(x(k))+Zi(x(k+1))Zi(x(k))=x(k)Qx(k)+ηT(x(k))Rη(x(k))+Zi1(x(k+1))
    对上式求差值
    Z i + 1 ( x ( k ) ) − Z i ( x ( k ) ) = Z i ( x ( k + 1 ) ) − Z i − 1 ( x ( k + 1 ) ) = Z i − 1 ( x ( k + 1 ) ) − Z i − 2 ( x ( k + 1 ) ) = Z i − 2 ( x ( k + 1 ) ) − Z i − 3 ( x ( k + 1 ) ) ⋅ ⋅ = Z 1 ( x ( k + i ) ) − Z 0 ( x ( k + i ) ) \begin{aligned} Z_{i+1}(x(k))-Z_i(x(k))&=Z_i(x(k+1))-Z_{i-1}(x(k+1)) \\ &=Z_{i-1}(x(k+1))-Z_{i-2}(x(k+1)) \\ &=Z_{i-2}(x(k+1))-Z_{i-3}(x(k+1)) \\ &\cdot \\ &\cdot \\ &=Z_1(x(k+i))-Z_0(x(k+i)) \end{aligned} Zi+1(x(k))Zi(x(k))=Zi(x(k+1))Zi1(x(k+1))=Zi1(x(k+1))Zi2(x(k+1))=Zi2(x(k+1))Zi3(x(k+1))=Z1(x(k+i))Z0(x(k+i))
    Z 0 ( x ( k ) ) = 0 Z_0(x(k))=0 Z0(x(k))=0 Z i + 1 ( x ( k ) ) − Z i ( x ( k ) ) = Z 1 ( x ( k + i ) ) − Z 0 ( x ( k + i ) ) Z_{i+1}(x(k))-Z_i(x(k))=Z_1(x(k+i))-Z_0(x(k+i)) Zi+1(x(k))Zi(x(k))=Z1(x(k+i))Z0(x(k+i))
    Z i + 1 ( x ( k ) ) = Z 1 ( k + i ) + Z i ( x ( k ) ) = Z 1 ( x ( k + i ) ) + Z 1 ( x ( k + i − 1 ) ) + Z i − 1 ( x ( k ) ) = Z 1 ( x ( k + i ) ) + Z 1 ( x ( k + i − 1 ) ) + Z 1 ( x ( k + i − 2 ) ) + Z i − 2 ( x ( k ) ) = Z 1 ( x ( k + i ) ) + Z 1 ( x ( k + i − 1 ) ) + Z 1 ( x ( k + i − 2 ) ) + ⋯ + Z 1 ( x ( k ) ) \begin{aligned} Z_{i+1}(x(k))&=Z_1(k+i)+Z_i(x(k)) \\ &=Z_1(x(k+i))+Z_1(x(k+i-1))+Z_{i-1}(x(k)) \\ &=Z_1(x(k+i))+Z_1(x(k+i-1))+Z_1(x(k+i-2))+Z_{i-2}(x(k)) \\ &=Z_1(x(k+i))+Z_1(x(k+i-1))+Z_1(x(k+i-2))+\dots +Z_1(x(k)) \end{aligned} Zi+1(x(k))=Z1(k+i)+Zi(x(k))=Z1(x(k+i))+Z1(x(k+i1))+Zi1(x(k))=Z1(x(k+i))+Z1(x(k+i1))+Z1(x(k+i2))+Zi2(x(k))=Z1(x(k+i))+Z1(x(k+i1))+Z1(x(k+i2))++Z1(x(k))
    因此上式可以写为
    Z i + 1 ( x ( k ) ) = ∑ j = 0 i Z 1 ( x ( k + j ) ) = ∑ j = 0 i ( x T ( k + j ) Q x ( k + j ) + η T ( x ( k + j ) ) R η ( x ( k + j ) ) ) ≤ ∑ j = 0 ∞ ( x T ( k + j ) Q x ( k + j ) + η T ( x ( k + j ) ) R η ( x ( k + j ) ) ) \begin{aligned} Z_{i+1}(x(k))&=\sum_{j=0}^iZ_1(x(k+j)) \\ &=\sum_{j=0}^i(x^T(k+j)Qx(k+j)+\eta^T (x(k+j))R\eta (x(k+j))) \\ &\leq\sum_{j=0}^{\infty}(x^T(k+j)Qx(k+j)+\eta^T (x(k+j))R\eta (x(k+j))) \end{aligned} Zi+1(x(k))=j=0iZ1(x(k+j))=j=0i(xT(k+j)Qx(k+j)+ηT(x(k+j))Rη(x(k+j)))j=0(xT(k+j)Qx(k+j)+ηT(x(k+j))Rη(x(k+j)))
    默认系统是稳定的,即 x ( k ) → 0 x(k)\rightarrow0 x(k)0,当 k → ∞ k\rightarrow \infty k

    因此可得对任意的 i i i
    Z i + 1 ( x ( k ) ) ≤ ∑ i = 0 ∞ Z 1 ( x ( k + i ) ) ≤ Y Z_{i+1}(x(k))\leq\sum_{i=0}^{\infty}Z_1(x(k+i))\leq Y Zi+1(x(k))i=0Z1(x(k+i))Y
    有上面的可知
    V i + 1 ( x ( k ) ) ≤ Z i + 1 ( x ( k ) ) ≤ Y      ( 1 ) V_{i+1}(x(k))\leq Z_{i+1}(x(k))\leq Y ~~~~ (1) Vi+1(x(k))Zi+1(x(k))Y    (1)
    定义
    Φ i + 1 ( x ( k ) ) = x ( k ) Q x ( k ) + u T ( i + 1 ) R u ( i + 1 ) + Φ i ( x ( k + 1 ) ) \Phi_{i+1}(x(k))=x(k)Qx(k)+u^T(i+1)Ru(i+1)+\Phi_i(x(k+1)) Φi+1(x(k))=x(k)Qx(k)+uT(i+1)Ru(i+1)+Φi(x(k+1))
    使用归纳法进行证明 V i + 1 ( x ( k ) ) ≥ Φ i ( x ( k ) ) V_{i+1}(x(k))\geq \Phi_i(x(k)) Vi+1(x(k))Φi(x(k))

    V 0 = Φ 0 = 0 V_0=\Phi_0=0 V0=Φ0=0
    V 1 ( x ( k ) ) = x T ( k ) Q x ( k ) + u 0 T ( x ( k ) ) R u 0 ( x ( k ) ) V 1 ( x ( k ) ) − Φ 0 ( x ( k ) ) = x T ( k ) Q x ( k ) + u 0 T ( x ( k ) ) R u 0 ( x ( k ) ) ≥ 0 V 1 ( x ( k ) ) ≥ Φ 0 ( x ( k ) ) V_1(x(k))=x^T(k)Qx(k)+u_0^T(x(k))Ru_0(x(k))\\ V_1(x(k))-\Phi_0(x(k))=x^T(k)Qx(k)+u_0^T(x(k))Ru_0(x(k))\geq 0 \\ V_1(x(k))\geq\Phi_0(x(k)) V1(x(k))=xT(k)Qx(k)+u0T(x(k))Ru0(x(k))V1(x(k))Φ0(x(k))=xT(k)Qx(k)+u0T(x(k))Ru0(x(k))0V1(x(k))Φ0(x(k))
    假设 V i ( x ( k ) ) ≥ Φ i − 1 ( x ( k ) ) V_{i}(x(k))\geq \Phi_{i-1}(x(k)) Vi(x(k))Φi1(x(k))
    Φ i ( x ( k ) ) = x ( k ) Q x ( k ) + u T ( i ) R u ( i ) + Φ i − 1 ( x ( k + 1 ) ) V i + 1 ( x ( k ) ) = x T ( k ) Q x ( k ) + u i T ( x ( k ) ) R u i ( x ( k ) ) + V i ( x ( k + 1 ) ) V i + 1 ( x ( k ) ) − Φ i ( x ( k ) ) = V i ( x ( k + 1 ) ) − Φ i − 1 ( x ( k + 1 ) ) ≥ 0 \Phi_{i}(x(k))=x(k)Qx(k)+u^T(i)Ru(i)+\Phi_{i-1}(x(k+1)) \\ V_{i+1}(x(k))=x^T(k)Qx(k)+u_i^T(x(k))Ru_i(x(k))+V_i(x(k+1)) \\ V_{i+1}(x(k))-\Phi_{i}(x(k))=V_i(x(k+1))-\Phi_{i-1}(x(k+1))\geq 0 Φi(x(k))=x(k)Qx(k)+uT(i)Ru(i)+Φi1(x(k+1))Vi+1(x(k))=xT(k)Qx(k)+uiT(x(k))Rui(x(k))+Vi(x(k+1))Vi+1(x(k))Φi(x(k))=Vi(x(k+1))Φi1(x(k+1))0
    因此可以得出
    V i + 1 ( x ( k ) ) ≥ Φ i ( x ( k ) )           ( 2 ) V_{i+1}(x(k))\geq \Phi_i(x(k))~~~~~~~~~(2) Vi+1(x(k))Φi(x(k))         (2)
    有(1)和(2)可以得出
    V i ( x ( k ) ) ≤ Φ i ( x ( k ) ) ≤ V i + 1 ( x ( k ) ) ⇒ V i ( x ( k ) ) ≤ V i + 1 ( x ( k ) ) V_{i}(x(k))\leq \Phi_i(x(k))\leq V_{i+1}(x(k)) \\ \Rightarrow V_{i}(x(k))\leq V_{i+1}(x(k)) Vi(x(k))Φi(x(k))Vi+1(x(k))Vi(x(k))Vi+1(x(k))
    可以看出 V i V_i Vi 是非增的并且有上界,即 V i → V ∗ V_i\rightarrow V^* ViV,当 i → ∞ i\rightarrow \infty i

    证闭

    上面的证明说明数列 V i + 1 ( x ( k ) ) = arg ⁡ min ⁡ { x T ( k ) Q x ( k ) + u T ( x ( k ) ) R u ( x ( k ) ) + V i ( x ( k + 1 ) ) } V_{i+1}(x(k))=\arg \min \{x^T(k)Qx(k)+u^T(x(k))Ru(x(k))+V_i(x(k+1))\} Vi+1(x(k))=argmin{xT(k)Qx(k)+uT(x(k))Ru(x(k))+Vi(x(k+1))}一定收敛,并且收敛到贝尔曼方程,这是值迭代收敛严格的证明,但是值函数的初值必须为零。

    参考文献

    • Discrete-time nonlinear HJB solution using Approximate dynamic programming: Convergence Proof
    展开全文
  • 基于数据的自适应动态规划——总体最小二乘法
  • @TOC 自适应动态规划学习笔记(3) 第三天 ADP的三个部分 &emnp;书接上回,上图展示了ADP的三个基本的组成,其中Critic Network输出对函数$J$的估计值

    @TOC 自适应动态规划学习笔记(3)

    第三天(图全是偷的)

    在这里插入图片描述

    图1 ADP的三个部分

    Model Network

     书接上回,图(1)中所示的Model Network就是对于系统公式(1) x k + 1 = F ( x k , u k ) x_{k+1}=F(x_k,u_k) xk+1=F(xk,uk) F F F的拟合,可以提前离线训练好,也可以和Critic Network、Action Network一起训练。这部分就是利用神经网络系统辨识部分,等后面有空会再记录系统辨识的知识。(推荐《系统辨识理论及MATLAB仿真》刘金琨著)
    p:如果 F F F有确定的表达,依旧需要建立Model Network,这是必须的。至于为什么,我母鸡啊。

    Critic Network

     图(1)中展示展示的ADP的三个基本组成中,Critic Network输出对函数 J J J的估计值 J ^ \hat{J} J^,也就是对下式的估计(详见第一天)

    J ( x k , u k ‾ ) = ∑ i = k ∞ γ i − k U ( x i , u i ) (2) J(x_k,\underline{u_k})=\sum_{i=k}^{\infty}\gamma^{i-k}U(x_i,u_i)\tag{2} J(xk,uk)=i=kγikU(xi,ui)(2)

    Critic Network的损失函数或者说平方误差由下式确定
    ∣ ∣ E h ∣ ∣ = 1 2 ∑ k E k 2 = 1 2 ∑ k ( J ^ k − U k − γ J ^ k + 1 ) 2 (6) ||E_h||=\frac{1}{2}\sum_{k}E_k^2=\frac{1}{2}\sum_k(\hat{J}_k-U_k-\gamma\hat{J}_{k+1})^2\tag{6} Eh=21kEk2=21k(J^kUkγJ^k+1)2(6)其中 J ^ k = J ^ ( x k , W c ) \hat{J}_k=\hat{J}(x_k,W_c) J^k=J^(xk,Wc), W c W_c Wc表示Critic Network的参数(权值),函数 U k U_k Uk是与式(2)中相同的效用函数(utility function)用来反映系统的性能。在实际中,函数 U k U_k Uk通常是 x k x_k xk u k u_k uk的函数,可以表示为 U k = U ( x k , u k ) U_k=U(x_k,u_k) Uk=U(xk,uk),当对于任意的 k k k,都有 E k = 0 E_k=0 Ek=0时,式(6)变为
    J ^ k = U k + γ J ^ k + 1 = U k + γ ( U k + 1 + γ J ^ k + 2 ) = ⋯ = ∑ i = k ∞ γ i − k U i , (7) \begin{aligned} \hat{J}_k &= U_k +\gamma\hat{J}_{k+1}\\ &= U_k+\gamma(U_{k+1}+\gamma\hat{J}_{k+2})\\ &=\cdots\\ &=\sum_{i=k}^\infty\gamma^{i-k}U_i,\tag{7} \end{aligned} J^k=Uk+γJ^k+1=Uk+γ(Uk+1+γJ^k+2)==i=kγikUi,(7)
     这与式(2)完全相同,因此可以将(6)作为损失函数来训练Critic Network,从而使其输出 J ^ \hat{J} J^成为公式(2)定义的成本函数 J J J的估计。
     得到Model Network后,对Critic Network进行训练,Critic Network会给出代价函数的估计,训练过程中很多标准的算法都是可以使用的,需要注意式(6)中的 J ^ k + 1 \hat{J}_{k+1} J^k+1是Critic Network在 k + 1 k+1 k+1时刻对代价函数 J J J的估计, J ^ k + 1 = J ^ ( x ^ k + 1 , W c ) \hat{J}_{k+1}=\hat{J}(\hat{x}_{k+1},W_c) J^k+1=J^(x^k+1,Wc),注意这个神经网络的参数 W c W_c Wc与对 J ^ \hat{J} J^进行估计的神经网络参数完全相同,而 x ^ k + 1 \hat{x}_{k+1} x^k+1则是Model Network对于 k + 1 k+1 k+1时刻状态 x k + 1 x_{k+1} xk+1的估计。
     Critic Network有两种训练方法:前向时间法(a forward-in-time approach)和后向时间法(a backward-in-time approach),图(2)表示了前向时间方法示意图,在这种方法中,待训练的Critic Network的输出是 J ^ k \hat{J}_k J^k,训练标签是 U k + γ J ^ k + 1 U_k+\gamma \hat{J}_{k+1} Uk+γJ^k+1,注意图中的对 J ^ k \hat{J}_k J^k J ^ k + 1 \hat{J}_{k+1} J^k+1的估计所用的神经网络参数完全一样,但是其输入却完全不同。图(3)表示了后向时间方法的示意图,这种方法中待训练网络输出 J ^ k + 1 \hat{J}_{k+1} J^k+1,标签为 ( J ^ k − U k ) / γ (\hat{J}_k-U_k)/\gamma (J^kUk)/γ,这两种方法的目的都是为了在满足式(7)的情况下最小化式(6). x ^ k + 1 \hat{x}_k+1 x^k+1是Model Network的输出。
    前向时间


    图2 前向时间方法

    后向时间


    图3 后向时间方法

     参考强化学习中的TD算法(temporal difference method ),其通过使用 ∣ r t + 1 + γ V ( s t + 1 ) ∣ |r_{t+1}+\gamma V(s_t+1)| rt+1+γV(st+1)作为标签来使得学习目标是 ∣ r t + 1 + γ V ( s t + 1 ) − V ( s t ) ∣ |r_{t+1}+\gamma V(s_t+1)-V(s_t)| rt+1+γV(st+1)V(st)最小化。这与图(2)所示的前向时间方法的标签 ∣ U k + γ J ^ k + 1 ∣ |U_k+\gamma \hat{J}_{k+1}| Uk+γJ^k+1相似。(所以说这俩同源?要学的东西又增多了(微笑) )

    Action Network

     在Critic Network训练完成后,就可以以最小化 U k + γ J ^ k + 1 U_k+\gamma\hat{J}_{k+1} Uk+γJ^k+1来训练Action Network,从而输出控制信号 u k = u ( x k , W a ) u_k=u(x_k,W_a) uk=u(xk,Wa),其中 W a W_a Wa是Action Network 的参数(权值),如果通过这种方式训练完成Action Network,我们就可以得到一个可以根据Critic Network的性能产生一个最优的,至少是次优的控制信号作为其输出。回想下动态规划的目标:获得式(5)(见第二天)中的最优控制序列,这个控制序列会使得(2)的函数 J J J最小。

    u k ∗ = arg min ⁡ u k { U ( x k , u k ) + γ J ∗ ( x k + 1 ) } (5) u_k^*=\argmin_{u_k}\lbrace U(x_k,u_k)+\gamma J^*(x_{k+1})\rbrace\tag{5} uk=ukargmin{U(xk,uk)+γJ(xk+1)}(5)

     这里的关键就是通过对代价函数 J J J的估计,从而建立了当前行动(action)和未来的结果(consequences)之间的一种可以互相“交流”的联系。

    The key here is to interactively build a link between present actions and future consequences via an estimate of the cost function

     在Action Network训练完成后,就可以检查系统性能是否达标,不达标的话可以通过调整Critic Network来继续训练,重复这一过程直到性能达标。这三个神经网络是通过图(1)的方式连接在一起的,在三个网络同时训练的时候,作为过程的一部分,控制信号 u k u_{k} uk将会作用在外部环境得到 x k + 1 x_{k+1} xk+1同时作用在Model Network上得到 x ^ k + 1 \hat{x}_{k+1} x^k+1,与此同时,通过最小化 ∣ ∣ x k + 1 − x ^ k + 1 ∣ ∣ ||x_{k+1}-\hat{x}_{k+1}|| xk+1x^k+1来更新Model Network。
     Action Network的训练通过最小化 U k + γ J ^ k + 1 U_k+\gamma\hat{J}_{k+1} Uk+γJ^k+1来完成,同时要保持Critic Network和Model Network参数固定不变,梯度信息通过Critic Network和Model Network反向传播(backward propagation)到Action Network,这三个神经网络组成了一个大的前馈网络(feedforward network)。这就意味着即使 F F F有着明确的表达式,Model Network是实现ADP所必需的,便于通过反向传播算法训练Action Network。(等我搞懂我偏不用神经网络

    ADP分类(都是我瞎翻译的)

    根据Werbos的论文,ADP可以分为几种主要的方案:启发式动态规划(heuristic dynamic programming ,HDP),动作依赖的HDP(action-dependent HDP ,ADHDP)、双启发HDP(dual HDP ,DHP)、动作依赖的DHP(ADDHP)、全球化DHP(globalized DHP,GDHP)以及动作依赖的GDHP(ADGDHP)。HDP是最基本的版本,就是图(1)所示的内容。Werbos论证了ADP和强化学习的关系(TD与HDP,ADHDP和Q-learning)。
    (Werbos yyds,具体论证有空再看,先学会HDP再说)

    展开全文
  • 自适应动态规划(二)

    千次阅读 热门讨论 2020-10-10 08:31:44
    自适应动态规划(二) 贝尔曼公式和离散LQR 一个离散系统 x(k+1)=Ax(k)+Bu(k) x(k+1)=Ax(k)+Bu(k) x(k+1)=Ax(k)+Bu(k) 性能指标函数 J(k)=12∑i=k∞(xT(i)Qx(i)+uT(i)Ru(i)) J(k)=\frac{1}{2}\sum_{i=k}^{\infty}(x^...
  • @[toc]自适应动态规划学习笔记(2) 第二天 动态规划的基本原则是贝尔曼的最优性原则,简单描述为:  多级决策过程的最优策略,不论其初始状态和初始决策如何,当把其中任何一级和状态作为初始级和初始状态时,其...
  • The main purpose of this book is to introduce the recently developed framework, known as robust adaptive dynamic programming (RADP), for datadriven, nonmodel-based adaptive optimal control design for ...
  • 在线自适应动态规划之持续激励条件 在在线自适应动态规划算法中,需要通过满足持续激励条件(Persistent Excitation Condition)保证critic网络权值的收敛性,一般的做法都是在算法初期在输入端引入探测噪声...
  • ADP(自适应动态规划)-HDP

    千次阅读 多人点赞 2020-07-10 16:43:17
    ADP   自适应动态规划(Adaptive/Approximate Dynamic Programming,ADP)是在动态规划的基础上发展起来的,最早由 Paul J. Werbos 提出。ADP 的模型框架有:启发式动态规划(Heuristic Dynamic Programming,HDP),...
  • 强化学习和自适应动态规划 本文主要记录一下控制领域强化学习和自适应动态规划的发展,主要分为如下几个方向展开: 以早期Werbos提出Actor-Critic结构的Adaptive Dynamic Programming,并大致分成四类结构,包括...
  • ADP(自适应动态规划)-扩展HDP

    千次阅读 2020-07-13 00:33:08
    ADP   自适应动态规划(Adaptive/Approximate Dynamic Programming,ADP)是在动态规划的基础上发展起来的,最早由 Paul J. Werbos 提出。ADP 的模型框架有:启发式动态规划(Heuristic Dynamic Programming,HDP),...
  • 策略迭代 策略迭代稳定性证明 单调不增的证明 迭代过程 取一个随机容许初始控制律v0(xk)v_0(x_k)v0​(xk​) V0(xk)=U(xk,v0(xk))+V0(xk+1) V_0(x_k)=U(x_k,v_0(x_k))+V_0(x_{k+1}) V0​(xk​)=U(xk​,v0​(xk​))+V0...
  • ADP(自适应动态规划)-值迭代

    千次阅读 多人点赞 2020-07-10 16:55:44
        看网上的ADP的代码挺少的,最近写了一个ADP值迭代的代码,分享一下,接下来也准备写Actor-Critic框架。 1、ADP值迭代原理1     ADP值迭代和强化学习的值迭代很类似,在ADP中的值迭代分为传统的值迭代和广义...
  • 为了获得非线性离散时间系统的最优控制策略, 基于自适应动态规划的原理, 提出了一种带误差限的自适 应动态规划方法. 对于一个任意的状态, 用一个有限长度的控制序列近似最优控制序列, 使性能指标与最优性能指标...
  • 自适应动态规划国内研究人员

    千次阅读 2018-01-27 21:02:36
    ADP与制导律结合 导弹,作为一种强威慑与大威力打击武器,是军事强国军事变革...为了提升导弹制导全方位智能化和自主化要求,将自适应动态规划(ADP)技术与导弹制导方法相结合,从而将自主智能思想引入制导律设计过
  • 对此,微软 AI 认知服务团队提出了动态卷积,与传统的静态卷积(每层单个卷积核)相比,根据注意力动态叠加多个卷积核不仅显著提升了表达能力,额外的计算成本也很小,因而对高效的 CNN 更加友好,同时可以容易地整合...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 8,207
精华内容 3,282
关键字:

自适应动态规划