非线性系统自学习最优控制：自适应动态规划方法(英文版)思维导图

非线性系统自适应最优切换控制方法毛艳岭 1富 月1摘 要 针对具有未知动态和M 个平衡点的连续时间非线性系统, 将线性自适应最优切换控制器和未建模动态补偿器相结合, 基于嵌入转换技术和近似动态规划思想, 提出一种自适应最优切换控制方法. 首先在非线性系统的M 个平衡点建立M 个线性化模型, 当模型参数已知时, 提出由线性最优切换控制器、切换准则、未建模动态补偿器以及非线性系统组成的控制系统结构; 当模型参数未知时, 在每个平衡点附近采集输入和状态数据, 利用黎卡提方程的迭代求解公式、最小二乘方法、极小值原理以及二次规划技术得到非线性系统的自适应最优切换控制器和最优切换序列; 最后进行仿真实验, 实验结果验证了所提方法的有效性、优越性和实际可应用性.关键词 非线性系统, 切换控制, 自适应最优控制, 嵌入转换引用格式 毛艳岭, 富月. 非线性系统自适应最优切换控制方法. 自动化学报, 2023, 49(10): 2122−2135DOI 10.16383/j.aas.c220180Adaptive Optimal Switching Control of Nonlinear SystemsMAO Yan-Ling 1 FU Yue 1Abstract In this paper, for continuous-time nonlinear systems with unknown dynamics and M equilibrium points,based on embedding-transformation and approximate dynamic programming, an adaptive optimal switching control method is proposed by combining a linear adaptive optimal switching controller and an unmodeled dynamic com-pensator. Firstly, M linearized models are established at M equilibrium points of the nonlinear system. When the model parameters are known, a control system structure consisting of a linear optimal switching controller, a switching mechanism, an unmodeled dynamic compensator, and the nonlinear system is proposed. When the model parameters are unknown, the input and state data are collected at the neighborhood of each equilibrium point.Then the adaptive optimal switching controller and optimal switching sequence are obtained by using the iterative Riccati equation, least square method, minimum principle, and quadratic programming. Finally, simulations are conducted, and the results verify the effectiveness, superiority and applicability of the proposed method.Key words Nonlinear systems, switching control, adaptive optimal control, embedding-transformationCitation Mao Yan-Ling, Fu Yue. Adaptive optimal switching control of nonlinear systems. Acta Automatica Sin-ica , 2023, 49(10): 2122−2135实际工业过程的被控对象大多是非线性的, 比如电镕镁砂熔炼过程的电极、钢球磨煤机制粉过程的磨机等等. 非线性系统结构复杂, 往往难以得到精确的数学模型, 其控制问题一直是控制领域相关学者和工程师的研究难点和热点之一.经典的非线性控制方法, 如反馈线性化方法[1−2],由于需要已知精确的数学模型, 无法应用到实际的工业过程中. 为了解决这个问题, 文献[3]针对具有全状态约束的高阶非线性随机系统, 利用模糊逻辑系统逼近未知非线性函数, 提出了一种新的模糊自适应反步控制方法. 文献[4]在文献[3]的基础上,针对具有指数型性能函数的高阶非线性随机系统,提出了基于模糊逻辑系统和反步法的模糊自适应有限时间跟踪控制方法. 当被控对象的非线性较弱或在某一平衡点附近运行时, 通常采用近似线性模型进行描述, 并针对该模型设计控制器. 例如, 文献[5]利用递归近似理论, 将非线性系统看作线性时变序列系统的极限, 针对线性时变序列系统设计线性二次最优序列控制器, 从而实现原非线性系统的二次最优控制. 文献[6]利用泰勒公式将非线性系统在某一平衡点附近表示为线性模型与高阶非线性项的组合, 将开环解耦补偿器、非线性神经网络补偿器和一步超前最优加权自适应控制器结合, 提出了非线性系统基于神经网络的自适应动态解耦控制方法. 文献[7]考虑到模型阶次的不匹配问题, 通过引入降阶模型, 采用带死区的归一化投影算法对线性收稿日期 2022-03-16 录用日期 2023-02-24Manuscript received March 16, 2022; accepted February 24,2023国家自然科学基金(62333004, 61991403, 61991400, 61873052)资助Supported by National Natural Science Foundation of China (62333004, 61991403, 61991400, 61873052)本文责任编委 赵旭东Recommended by Associate Editor ZHAO Xu-Dong1. 东北大学流程工业综合自动化国家重点实验室 沈阳 1108191. State Key Laboratory of Synthetical Automation for Pro-cess Industries, Northeastern University, Shenyang 110819第 49 卷 第 10 期自 动 化 学 报Vol. 49, No. 102023 年 10 月ACTA AUTOMATICA SINICAOctober, 2023模型参数进行辨识, 利用高阶神经网络估计高阶非线性项, 将带有滤波器的极点配置自适应比例积分微分(Proportional integral derivative)控制器与神经网络补偿器相结合, 提出了非线性系统基于神经网络的自适应PID 控制方法. 神经网络收敛速度较慢且容易陷入局部极小点, 高阶非线性项的估计精确度较低. 为了解决这一问题, 文献[8]首次引入了控制器驱动模型和虚拟未建模动态的概念, 基于线性控制器驱动模型构造一步超前最优自适应控制器, 结合虚拟未建模动态补偿器, 提出了非线性系统自适应切换控制方法. 文献[9]针对复杂的热交换过程, 设计了具有虚拟未建模动态补偿的一步最优比例积分(Proportional integral)控制器, 并提出了数据驱动的双速率控制方法. 上述控制方法虽然能够取得良好的控制效果, 但是当系统的非线性较强或平衡点发生变化时, 这种只考虑单一平衡点的控制方法往往会使控制性能下降甚至导致整个系统失稳.M k -很多实际工业过程的平衡点都会随着工况的不同而发生变化, 比如电熔镁砂熔炼过程的平衡点随着原料成分和加料阶段的不同会发生变化; 钢球磨煤机制粉系统中磨机的平衡点随着原煤成分和湿度的不同而发生变化. 本文针对一类具有 个平衡点的非线性系统, 研究基于多模型切换的自适应控制方法. 多模型自适应控制方法一般用于改善系统的暂态性能或解决参数跳变系统的控制问题, 如文献[10]针对一类连续时间线性系统, 为改善系统的暂态性能, 提出了基于直接模型参考自适应控制的多模型切换控制方法. 文献[11]针对一类参数跳变离散时间线性系统, 提出了基于间接自校正控制的多模型切换控制方法. 文献[12]针对一类参数跳变离散时间非线性系统, 通过引入 差分算子, 分别设计了线性自适应控制器和基于神经网络的非线性自适应控制器, 通过两个控制器之间的切换, 可以提高系统的性能和稳定性. 为了避免不良切换行为,文献[13]采用滞后切换逻辑消除了参数估计器对初始条件的依赖, 通过利用鲁棒线性时不变工具实现高性能的控制目标, 结合控制器混合策略, 提出了多模型自适应混合控制方法. 针对文献[13]所提方法需要模型数量大的问题, 文献[14]采用分离处理原则, 充分利用所有辨识模型信息, 采用二级自适应方法建立自适应控制器. 为了消除系统非线性项对控制输入应严格线性的限制, 文献[15]针对离散时间非线性系统, 采用极点配置控制方法, 提出了由线性间接自校正控制器、基于神经网络的非线性间接自校正控制器和切换机制组成的多模型自适应控制器. 很多研究将多模型自适应控制方法应用到实际系统中, 并且取得了较好的控制效果. 文献[16]将多模型自适应切换控制方法应用于电力系统低频振荡中, 建立了不同工况下的线性小信号模型, 采用递归贝叶斯方法计算每个模型代表实际电力系统的概率, 根据这个概率得到每个控制器输出的占比权重, 最终的控制输出即为每个控制器输出的概率加权平均值. 文献[17]针对动态特性随不同负载状态而变化的柔性传送系统, 分别在不同负载状态处建立线性模型, 提出了基于闭环输出误差最小化的参数估计算法和基于极点配置的多模型自适应切换控制方法. 文献[18]以钢球磨煤机制粉系统为例, 针对一类具有多变量强耦合强非线性且动态特性随不同运行条件而变化的复杂工业过程, 将其在不同平衡点处用不同的线性模型和非线性未建模动态项组成的估计模型来描述, 提出了由非线性解耦控制器、线性解耦控制器和多模型切换机制组成的智能解耦控制方法. 文献[19]针对串联电容补偿输电线路的风力系统次同步谐振问题, 采用传统线性控制方法设计控制器, 根据系统条件设计该控制器的监控控制器, 该方法之后被拓展到了双馈异步发电机在串联补偿输电系统中的次同步振荡问题[20].上述多模型控制方法中, 用于切换的控制器是针对单一时刻的性能指标设计的, 具有次优性, 无法保证切换序列和控制系统的最优性.M M M M M 在实际工业生产过程中, 保证控制系统性能最优对实现工业过程整体优化控制是至关重要的. 本文针对具有未知动态和 个平衡点的连续时间非线性系统, 将嵌入转换法和近似动态规划技术相结合, 提出了一种自适应最优切换控制方法, 一方面能够保证切换序列的最优性, 另一方面可以实现控制系统的最优性能, 改善控制系统的动态品质. 首先在非线性系统的 个平衡点附近采集 组输入和状态数据, 利用黎卡提方程的迭代求解公式和最小二乘方法得到针对每个线性模型的最优控制器增益的估计, 利用极小值原理得到 个近似线性化模型. 然后利用嵌入转换法将 个近似线性化模型嵌入到一个连续时间大系统中, 通过二次规划技术得到非线性系统的线性自适应最优切换控制器和最优切换序列. 最后, 将线性自适应最优切换控制器和未建模动态补偿器相结合, 实现了控制目标. 仿真实验验证了本文所提方法的有效性、优越性和实际可应用性.M 本文针对具有未知动态和 个平衡点的连续时间非线性系统, 提出了自适应最优切换控制方法.主要创新点如下:1) 提出了由线性最优切换控制器、切换准则和未建模动态补偿器组成的控制器结构;10 期毛艳岭等: 非线性系统自适应最优切换控制方法2123M M 2) 模型参数已知时, 基于嵌入转换技术提出了由 个模型、 个最优控制器和切换准则组成的线性最优切换控制器;M M 3) 模型参数未知时, 基于嵌入转换技术和近似动态规划思想提出了由 个近似线性化模型、 个自适应最优控制器和切换准则组成的线性自适应最优切换控制器.1 问题描述M 考虑由如下模型描述的具有 个平衡点的连续时间非线性非仿射系统:x (t )=[x 1(t ),x 2(t ),···,x n (t )]T n u (t )=[u 1(t ),u 2(t ),···,u m (t )]T m f (x (t ),u (t ))=[f 1(·,·),f 2(·,·),···,f n (·,·)]T :R n ×R m →R n 其中 是 维状态向量, 是 维控制输入向量, 表示连续可微的未知非线性向量函数.M u (t )本文的目标是针对具有 个平衡点的未知非线性系统(1), 寻找最优切换序列和自适应最优切换控制律 , 使得闭环系统渐近稳定.M i ∈{1,2,···,M }(x i ,u i )非线性非仿射系统结构复杂, 很难直接根据它的模型设计控制器. 通常的做法是将非线性系统在某一平衡点附近线性化, 针对等价的近似线性模型设计控制器, 从而实现对原非线性系统的有效控制,如文献[4−5]等. 为此本文将非线性系统(1)在 个平衡点附近泰勒展开, 得到第 个平衡点 附近的等价近似线性模型:˙x (t )=A i x (t )+B i u (t )i i A i =∂f ∂xu =u i x =x iB i=∂f ∂uu =u i x =x i(A i ,B i )v i (t )i M 等价模型(2)包括两部分, 第一部分 表示第 个平衡点附近的线性化模型 , 其中 和 为适当维数的未知常值矩阵且 可控; 第二部分 为第 个平衡点附近的未建模动态. 为建立非线性系统(1)在 个平衡点附近的控制器设计模型, 引入如下单位脉冲序列记号b 其中 是整数, 则系统(1)可表示为σ(t )∈{1,2,···,M }其中 表示切换信号. 与此同时,本文所提出的控制器结构也包括两部分, 第一部分根据基于线性化模型建立的如下控制器设计模型进行设计:第二部分根据线性化产生的建模误差来设计, 用于消除未建模动态影响, 实现闭环系统渐近稳定.∑M i =1δ(σ(t )−i )∑i δ(σ(t )−i )在不引起混淆的情况下, 接下来我们将 简化为 .2 自适应最优切换控制器设计2.1 参数已知时的最优切换控制器A iB i i =1,···,M 当 和 ( )已知时, 我们提出了如图1所示的由线性最优切换控制器、切换准则、未建模动态补偿器以及非线性系统组成的控制系统结构, 其中线性最优切换控制器和切换准则根据控制器设计模型(4), 利用嵌入转换法[21]、极小值原理和二次规划方法获得; 未建模动态补偿器根据非线性系统状态和最优模型状态之间的误差设计.δ(σ(t )−i )首先令 在区间[0, 1]内连续变化, 利用嵌入转换法将式(4)嵌入到一个连续时间大系统中. 然后根据该嵌入式连续时间大系统的最优控制问题:δ(σ(t )−i )∈[0,1]Q 、R (A σ(t ),√Q )其中 ,为适当维数的参数矩阵且 可观, 采用极小值原理和二次规划方法得到切换准则函数:P σ(t )其中 根据如下黎卡提方程求解:J σ(t )J σ(t )每一时刻,比较 , 选择与最小的 对应的线性最优切换控制律:σ(t )K σ(t )其中 为最优切换序列, 表示线性最优切换控制器的增益, 通过下式求解:接下来, 为消除未建模动态对控制系统性能的影响, 我们设计了如下未建模动态补偿器:2124自 动 化 学 报49 卷a 1∈R m ×n a 2e m =x −x ∗x ∗σ(t )其中 为可调参数矩阵, 为可调参数, 为建模误差, 为最优线性化模型 的状态.A iB i i =1,···,M 综上, 和 ( )已知时最优切换控制律为:注 1. 线性最优切换控制律和最优切换序列推导过程见附录A.δ(σ(t )−i )δ(σ(t )−i )注 2. 针对控制器设计模型(4), 通过嵌入扩大 的取值范围, 令 在区间[0, 1]内连续变化, 将由多个近似线性模型组成的式(4)嵌入到一个连续时间大系统中; 通过转换将针对控制器设计模型(4)的最优切换控制问题转化为针对该嵌入式连续时间大系统的最优切换控制问题.2.2 参数未知时的自适应最优切换控制器A i B i (i =1,···,M )P σ(t )M M σ(t )ˆKσ(t )ˆPσ(t )P σ(t )A σ(t )M M 当 和 未知时, 无法通过式(7)得到 , 无法得到如式(6)所示的切换准则函数和式(8)所示的线性最优切换控制律. 为解决这一问题, 本文提出了一种自适应最优切换控制方法. 首先在非线性系统的 个平衡点附近采集 组输入、状态数据, 利用黎卡提方程的迭代求解公式和最小二乘算法得到针对线性化模型 的自适应最优控制器增益 以及黎卡提方程近似解, 并根据贝尔曼方程得到 的估计, 从而得到 个平衡点附近的 个线性化模型; 然后M 将 个线性化模型嵌入到一个连续时间大系统中,针对该嵌入式连续时间大系统基于极小值原理和二次规划技术设计线性二次型最优控制律, 进而得到最优切换序列和线性自适应最优切换控制律; 最后将线性自适应最优切换控制律和未建模动态补偿器相结合应用到非线性系统中, 实现对未知动态非线性系统的自适应最优切换控制.A iB i i =1,···,M 针对控制器设计模型(4), 当 和 ( )已知时, 根据Kleinman 定理[22], 很容易得到如下推论:K σ(t ),0∈R m ×n σ(t )P σ(t ),k 推论 1. 令 为针对线性化模型 的稳定反馈控制器增益矩阵, 为下面李雅普诺夫方程的对称正定解:δ(σ(t )−i )∈{0,1}i δ(σ(t )−i )=1k =1,2,···K σ(t ),k 其中 且 , 表示迭代次数, 满足K σ(t ),k P σ(t ),k σ(t )K σ(t )P σ(t )则 和 分别收敛于针对线性化模型 的最优控制器增益 和黎卡提方程解 , 即A iB i 图 1 和 已知时的控制系统结构A iB i Fig. 1 Control system structure when and are known10 期毛艳岭等: 非线性系统自适应最优切换控制方法2125∑∑A iB i(i =1,···,M )定理 1. 针对控制器设计模型(4), 当 和未知时, 使性能指标最小的切换准则函数为:ˆKσ(t )σ(t )K σ(t )N σ(t )P σ(t )A σ(t )其中 是针对线性化模型 的最优控制器增益的估计, 根据式(17)求解; 是矩阵 的估计, 根据式(18)求解:Θσ(t )ˆ¯Pσ(t )Ξσ(t )vec (C )m ×n C mn ⊗其中 , 和 的定义见后文, 是把维矩阵 按列的顺序一列接一列地组成的 维向量, 代表克罗内克积,线性自适应最优切换控制律为:σ(t )J σ(t )其中 为与最小的 对应的最优切换序列.M σ(t )ˆKσ(t )ˆP σ(t )证明. 首先根据离线采集的 组输入、状态数据, 计算针对线性化模型 的自适应最优控制器增益 以及黎卡提方程近似解 . 受文献[23]启发, 将式(4)等价表示为:A σ(t ),k =A σ(t )−B σ(t )K σ(t ),k 其中 . 根据式(12)和式(13), 沿着式(20)的解, 可以得到∫Q σ(t ),k =Q +i δ(σ(t )−i )K Tσ(t ),k RK σ(t ),k A σ(t )B σ(t )∑i δ(σ(t )−i )[x T (A T σ(t ),k P σ(t ),k +P σ(t ),k A σ(t ),k )x ]−x T×Q σ(t ),k x B σ(t )∑i δ(σ(t )−i )B Tσ(t )P σ(t ),k ∑i δ(σ(t )−i )RK σ(t ),k +1其中 . 由此, 可以将包含未知矩阵 和 的 项用 代替. 同理, 可以将包含未知矩阵 的用 代替.由克罗内克积的定义, 可知I n n 其中, 表示 维单位矩阵. 定义如下运算l 对于正整数 , 定义矩阵2126自 动 化 学 报49 卷[∫∫0≤t 0<t 1<···<t l 其中 . 由式(22)和式(23)可知, 式(21)可等价表示为:Θσ(t),k 当 为列满秩矩阵时,σ(t )k K σ(t ),k +1由此, 可以得到线性化模型 第 次迭代的自适应最优控制器增益 .ˆKσ(t )K σ(t )Θσ(t )Ξσ(t )ˆKσ(t )ˆ¯Pσ(t )¯P σ(t )ˆPσ(t )令 为迭代终止时的自适应最优控制器增益并作为 的估计, 和 为迭代终止时的数据向量, 可以得到 的计算公式如式(17)所示. 令 为迭代终止时的 的估计, 由此可以得到黎卡提方程的近似解 .σ(t )P σ(t )B σ(t )P σ(t )A σ(t )M M A i B i 接下来, 针对线性化模型 求解矩阵 和 的估计, 从而得到 个平衡点附近的 个线性化模型. 当 和 已知时, 易知A iB i L σ(t )P σ(t )B σ(t ) 计, 则Pσ(t )A σ(t )σ(t ) 的估计可根据线性化模型 的贝尔曼方程得到, 易知σ(t )P σ(t )A σ(t )N σ(t )D σ(t )N σ(t )ˆPσ(t )M 将式(27)代入上式, 利用离线采集的第 组输入、状态数据, 通过求取最小二乘解可以得到如式(18)所示的矩阵 的估计 . 根据 , 以及 , 可以很容易得到 个平衡点附近的近似控制器设计模型:δ(σ(t )−i )[0,1]δ(t )=[δ(σ(t )−1),···,δ(σ(t )−M )]T W ={δ∈R M :∑i δ(σ(t )−i )=1,δ(σ(t )−i )≥0}最后求取最优切换序列和线性自适应最优切换控制律. 针对模型(29), 应用嵌入变换法, 使 在 内连续变化, 为此令 并记. 定义哈密顿函数:易知, 针对嵌入式近似控制器设计模型的最优控制律为:将式(31)代入式(30), 化简可得δ(σ(t )−i )H (x,δ)下面将 作为决策变量, 通过最小化, 可以得到最优切换序列.δ(σ(t )−i )H (x,δ)实际上, 选择 使 最小等价为使式(33)最小W ¯Hδ(σ(t )−i )∈{0,1}σ(t )这是一个二次规划问题, 由于 是凸集, 是凹函数, 该问题的全局极小值一定在 取得[21], 且该全局极小值对应的 即为最优切换序列. 由此可以得到如式(16)的切换准则函数和式(19)的线性自适应最优切换控制律. □L σ(t )ˆKσ(t )注3. 由式(27)可知 的估计精度由 10 期毛艳岭等: 非线性系统自适应最优切换控制方法2127ˆK σ(t )K σ(t )L σ(t )P σ(t )B σ(t )N σ(t )ˆKσ(t )的估计精度决定. 由文献[23]易知, 收敛于参数已知时的最优控制器增益 , 因此 收敛于 . 由式(28)可知 的估计精度由最小二乘估计算法的精度和 的估计精度共同决定.l 0l ≥l 0rank ([I xx ,I xu ])=n (n +1)2+mn Θσ(t ),k {P σ(t ),k }∞k =0{K σ(t ),k }∞k =0P σ(t )K σ(t )注4. 在每个平衡点附近, 如果存在正整数 ,使得对于任意 , 都有 , 即矩阵 是满秩的, 那么序列 和 分别收敛到黎卡提方程的解 和最优控制器增益 [23].未建模动态补偿器的设计与线性模型参数已知时的情况类似, 即a 1∈R m ×n a 2ˆe m =x −ˆx ∗ˆx ∗σ(t )其中 为可调参数矩阵, 为可调参数,为建模误差, 为最优线性化模型 的状态.A iB i (i =1,···,M )综上, 和 未知时自适应最优切换控制律为:自适应最优切换控制器设计流程如图2所示.3 仿真实验为了验证本文所提方法的有效性, 我们分别进行了模型参数已知时最优切换控制和模型参数未知时自适应最优切换控制的数值仿真实验, 并分别与单一的针对一个模型的最优控制器和自适应最优控制器进行了对比. 除此之外, 为了验证本文所提方法的实际可应用性, 我们进行了模型参数未知时双容水箱液位系统的自适应最优切换控制仿真实验.3.1 参数已知时最优切换控制数值仿真实验考虑如下连续时间非线性系统:x =[x 1,x 2]T ∈R 2u =[u 1,u 2]T ∈R 2其中 是状态向量, 是输入向量.u =[u 1,u 2]T =[−3,10]T ,[−2,10]T [−1,10]T ˙x=[˙x 1,˙x 2]T =[0,0]T [u r 1,u r 2,x r 1,x r 2]T [−3,10,−4.4685,0.5592]T [−2,10,−4.2642,0.7565]T [−1,10,我们的目标是针对已知的非线性系统(36), 寻找最优切换序列和最优切换控制律, 使得闭环系统渐近稳定. 为此, 首先分别将 和 施加到非线性系统(36)上, 并令 得到非线性系统(36)的三个平衡点, 即 = , 和 −4.0264,1.1119]T δ(σ(t )−i )∈{0,1}∑3i =1δ(σ(t )−i )=13. 将式(36)分别在上述三个平衡点处泰勒展开, 并令 且 , 可以得到非线性系统(36)在个平衡点附近的控制器设计模型:其中图 2 自适应最优切换控制器设计算法流程Fig. 2 Flow chart of adaptive optimal switchingcontrol algorithm2128自 动 化 学 报49 卷[]x (0)=[x 1(0),x 2(0)]T =[−4.4685,0.5592]T 接下来给定随机初始状态 , 并选择控制器参数矩阵和未建模动态补偿器参数t =50s t =100s 最后将最优切换控制器(6) ~ (11)加入到系统(36), 得到如图3所示的状态曲线, 如图4所示的控制输入曲线和如图5所示的最优切换序列. 结合图3和图4, 在 和 时, 虽然系统的平衡点发生变化, 但是采用本文提出的最优切换控制方法仍能够将状态很快调节到平衡点附近并保持不变.[u r 1,u r 2,x r 1,x r 2]T [−1,10,−4.0264,1.1119]T x (0)=[x 1(0),x 2(0)]T =[−4.4685,0.5592]T 为了验证本文所提最优切换控制方法的优越性, 我们与单一的针对一个模型的最优控制方法进行了对比实验. 以针对平衡点 = 处的线性化模型为例,给定初始状态 , 选择控制器参数矩阵如式(38)所示, 未建模动态补偿器参数如式(39)所示.[u r 1,u r 2,x r 1,x r 2]T [−1,10,−4.0264,1.1119]T 图6和图7分别为所得到的状态曲线和控制输入曲线. 根据图6和图7可以看出, 针对平衡点 = 处的线性化模型设计的控制器只能将状态调节到对应的平衡点附近. 当平衡点发生变化时, 系统的状态存在稳态误差. 但是由于平衡点的变化引起的建模误差可近似为常数, 因此状态曲线虽然偏离平衡点但恒定不变.3.2 参数未知时自适应最优切换控制数值仿真实验本节的目标是针对未知非线性系统(36), 寻找[u r 1,u r 2,x r 1,x r 2]T =[−2,10,−4.2642,0.7565]T [u r 1,u r 2,x r 1,x r 2]T =[2,10,−2.5517,3.6570]T [u 1,u 2]T =[sin (0.1t ),sin (0.5t )]T t =0s t =2s δxx ,I xx ,I xu 最优切换序列和自适应最优切换控制律, 使得闭环系统渐近稳定. 不失一般性, 这里我们以两个平衡点为例进行仿真实验. 结合图2, 首先分别在平衡点 和 附近施加激励输入信号, 即 ,从 到 , 以0.01 s 为采样周期, 分别采集201组输入和状态数据, 计算 . 选择−3.8−4.24−4.01−4.081.151.08100101100.8−4.290.760.7550.05050.6−4.2−4.61.4050100150050100150x 1, x r 1x 2, x r 20.90.4Time /sx 2x r 2x 1x r 1图 3 采用最优切换控制器时系统的状态Fig. 3 State curves of the system when usingthe optimal switching controller−u 1u 2Time /s图 4 采用最优切换控制器时系统的控制输入Fig. 4 Input curves of the system when usingthe optimal switching controller01234s (t )50100150Time /s图 5 采用最优切换控制器时系统的最优切换序列Fig. 5 Optimal switching sequence of the system whenusing the optimal switching controller10 期毛艳岭等: 非线性系统自适应最优切换控制方法2129控制器参数矩阵||P σ(t ),k −P σ(t ),k −1||≤10−3σ(t )=1,2k ˆPσ(t )ˆK σ(t )终止循环的条件为 , 其中 ; 代表迭代次数. 根据式(17)分别得到针对两个模型的 和 , 即:N σ(t )然后利用所采集的输入和状态数据求解式(18),分别得到针对两个模型的 , 即:[]最后, 根据式(29)可以得到两个线性化模型如下式所示:x (0)=[x 1(0),x 2(0)]T =[−4,0]T t 0=0s t max =100s 将两个线性化模型嵌入到一个连续时间大系统中, 结合图2, 给定初始状态 和初始时间 , 设置 , 选择未建模动态补偿器参数t ≥t max t =50s 将自适应最优切换控制器(35)加入到非线性系统, 当满足 时, 可以得到如图8所示的状态曲线, 如图9所示的控制输入曲线和如图10所示的切换序列. 在 , 由于平衡点突变, 切换序列发生改变, 导致系统的状态震荡, 经过1.8 s 的调节时间, 系统的状态被调节到平衡点附近并保持不变.[u r 1,u r 2,x r 1,x r 2]T =[−2,10,−4.2642,0.7565]T ˆP1ˆK 1N 1[u r 1,u r 2,x r 1,x r 2]T =[−2,10,−4.2642,0.7565]T t =50s 为了验证本文所提自适应最优切换控制方法的优越性, 我们以平衡点 为例, 与单一的针对一个模型的自适应最优控制方法进行了对比实验. 选择控制器参数矩阵如式(40), 根据式(17)和式(18)可以得到 , 和 分别如式(41)和式(42)所示, 根据式(29)可以得到线性化模型如式(43)所示, 选择未建模动态补偿器参数如式(44)所示. 所得到的状态曲线和控制输入曲线如图11和图12所示. 从图11和图12可以看出, 针对平衡点 设计的自适应最优控制器只能将状态调节到对应的平衡点附近. 与模型参数已知时情况相同, 当 时, 平衡点发生变化, 系统状态存在稳态误差. 但是由于平衡点的变化引起的建模误差可近似为常数, 因此状态曲−4.6−4.2−3.80.40.91.4x 1, x r 1x 2, x r 2x 1x r 1x 2x r 250100150050100150Time /s图 6 采用最优控制器时系统的状态Fig. 6 State curves of the system when usingthe optimal controller−Time /su 1u 2图 7 采用最优控制器时系统的控制输入Fig. 7 Input curves of the system when usingthe optimal controller2130自 动 化 学 报49 卷。

⾃适应动态规划（ADP）基础1 基础概念动态规划是利⽤最优性原理来解决最优和最优控制问题的⼀个⾮常有⽤的⼯具。

最优性原则可以表⽰为：“最优策略具有这样的性质:⽆论初始状态和初始决策是什么，其余决策都必须构成与第⼀个决策产⽣的状态相关的最优策略。

”动态规划有⼏个⽅⾯。

⼈们可以考虑离散时间系统或连续时间系统，线性系统或⾮线性系统，时不变系统或时变系统，确定性系统或随机系统，等等。

1.1 举例说明A.⾮线性离散时间(时变)动态(确定性)系统(1)系统定义x(k+1)=F[x(k),u(k),k],k=0,1,⋯（1）其中x∈R n代表系统的状态向量，u∈R n表⽰控制动作，F是系统函数。

(2)代价函数的定义（⽤于衡量控制控制系统的性能好坏，越⼩越好）J(x(i),i)=∞∑k=iγk−i U[x(k),u(k),k](2)其中U表⽰效⽤函数，γ为折扣因⼦，其取值范围为，0<γ≤1。

这⾥代价函数依赖于初始时间i和初始状态x(i)。

(3)动态规划的⽬标选择⼀个控制序列u(k),k=i,i+1,.....，使得代价函数J最⼩化(4)贝尔曼最优⽅程根据贝尔曼最优性原理，可以得出k时刻的最优代价等于：J⋆(x(k))=minu(k)(U(x(k),u(k))+γJ⋆(x(k+1)))(3)在k时刻的最优控制u⋆(k)是达到上述最⼩值的u(k)u⋆(k)=arg min u(k){U(x(k),u(k))+γJ∗(x(k+1))}(4)⽅程(3)是离散时间系统的最优性原则。

它的重要性在于，它允许通过回溯时间来⼀次优化⼀个控制向量。

B.在⾮线性连续时间的情况下(1)系统定义˙x(t)=F[x(t),u(t),t],t≥t0(5)(2)代价函数J(x(t))=∫∞t U(x(τ),u(τ))dτ(6)(3)最优代价满⾜哈密顿-雅可⽐-贝尔曼⽅程(Hamilton-Jacobi-Bellman Equation)−∂J⋆(x(t))∂t=minu∈U(U(x(t),u(t),t)+(∂J⋆(x(t))∂x(t))T×F(x(t),u(t),t)) =U(x(t),u⋆(t),t)+(∂J⋆(x(t))∂x(t))T×F(x(t),u⋆(t),t))(7)⽅程(3)和(7)称为动态规划的最优性⽅程，是实现动态规划的基础。

第四十一章最优控制和自适应控制及其智能控制第一节最优控制［!、"］一、基于变分法的最优控制问题求解最优控制是经典控制理论发展到现代控制理论的重要标志之一。

这里“最优”一词指的是相对于某一给定性能指标最优，如使控制过程的时间最短，燃料消耗最少，或者误差最小，而不是任何性能指标下都是最优的。

给定受控系统的状态方程!"##（"，$，%）寻求不受约束的控制向量$，使系统从初始状态"（%$）#"$在时间间隔［%$，%%］内转移到"（%%）且满足等式约束&［"（%%），%%］#$这里&为’维向量函数；并使指标&取极值(#)［"（%%），%%］’!%%%$*（"，$，%）(%利用变分法求解最优控制时，首先构造哈密尔顿函数+和增广泛函(,。

+#*（"，$，%）’!)#（"，$，%）#+（"，$，!，%）(,#)［"（%%），%%］’")&［"（%%），%%］’!%%%$［+（"，$，!，%）*!)!"］(%式中，"为-维、.为’维拉格郎日乘子向量。

由变分#(,#$导出的极值必要条件为：伴随方程!!#*"+"!状态方程!"##（"，$，%）#"+"!控制方程"+"$#$终端约束&［"（%%），%%］#$横截条件!（%%）#")""（%%）’"&)""（%%）.用计算机联立求解上面五个方程，可得到最优控制问题的数值解。

二、极小值原理与动态规划用变分法求解最优控制问题时，均假定控制!不受约束，并且存在惟一的偏导数!"#!!。

然而任何实际的控制量均限制在允许范围内变化，即!"!或!!$!#%&，$"#，…，’有些问题中!"#!!不存在，在这些情况下，可利用极小值原理求解。

运筹学B（双语）复习纲要I 概念汇总1）运筹学模型的三要素2）LP标准形式3）（非）基变量（入基变量、出基变量）4）基解、基可行解、解基逆矩阵5）退化6）人工变量法、两阶段法7）LP解的四种情况8）LP对偶问题的形式及最优解9）LP对偶问题的经济解释10）运输表格及表上作业法（（非）基变量、入基变量、出基变量）11）网络的基本概念（点、边、权、有向边、链、道路、圈、回路、树、生成树、最小生成树、连通图、割）12）目标规划的含义（模型、偏差变量）13）整数规划模型（分支定界法、割平面法的思路）II 方法汇总1）LP问题图上求解法2）单纯形法3）对偶单纯形法4）运输问题表上作业法（三大步骤）5）求最小生成树6）求最短路问题7）求最大流问题8）目标规划的图解法9）中国邮路问题III题型1）多选题：20分左右2）判断题：10分左右3）简答题：30分左右4）计算题：50分左右《运筹学B》双语课程词汇表Chapter 1 What is Operations Research?Operations Research 运筹学Mathematic model 数学模型decision alternative 决策选择decision variable 决策变量restriction，constraint 约束条件objective criterion 目标准则objective function 目标函数linear programming 线性规划integer programming 整数规划dynamic programming 动态规划network programming 网络规划nonlinear programming 非线性规划algorithm 算法iteration 迭代Chapter 2 Introduction to Linear Programming Graphical solution 图解法Graphical sensitivity analysis 图上灵敏度分析nonnegativity restrictions 非负约束条件feasible solution 可行解optimal feasible solution 最优可行解coefficient 系数denominator 分母infeasible 不可行unit worth 单位价值Chapter 3 The Simplex Methodsolution space 解空间algebraic solution 代数解graphical solution 几何解optimal solution 最优解equation 方程corner point 顶点basic variable 基变量nonbasic variable 非基变量basic solution 基解The Simplex Method 单纯形法iterative 迭代的origin 原点leaving variable 出基变量entering variable 入基变量ratio 比率Gauss-Jordan row operation 高斯-约当行变换pivot column 主列pivot row 主行pivot element 主元素artificial variable 人工变量M-Method 大M方法Two-Phase Method 两阶段方法penalty 罚数degeneracy 退化degenerate 退化的alternative optima 多重最优解infinity 无穷unbounded 无界的pseudo-optimal solution 伪解Chapter 4 Duality and Sensitivity Analysis dual problem 对偶问题primal problem 原问题matrix 矩阵vector 向量identity matrix 单位矩阵verify 证明dual simplex method 对偶单纯形法generalized simplex method 广义单纯形法Chapter 5 Transportation Model and Its Variants Transportation Model 运输模型nontraditional Transportation Model 非典型运输模型The Transportation Algorithm 运输算法source 出发地destination 目的地node 节点arc 边，弧transportation tableau 运输表格balanced 平衡的Northwest-Corner Method 西北角法Least-Cost Method 最小费用法V ogel Approximation Method 沃格尔法The Assignment Model 指派模型Hungarian Method 匈牙利方法Chapter 6 Network Modelsnetwork 网络Minimal Spanning Algorithm 最小生成树算法Shortest-Route Algorithm 最短路算法path 链connected network 连通网络cycle 回路spanning tree 生成树maximal flow 最大流residue network 剩余网络breakthrough path 关键路线Chapter 8 Goal Programminggoal programming 目标规划deviational variable 偏差变量Chapter 9 Integer Linear ProgrammingInteger Linear Programming 整数线性规划integer variables 整数变量Cutting-Plane Algorithm 割平面法B&B Algorithm 分支定界法。

综述非线性系统最优控制理论近年来，最优控制理论[1,2]的研究，无论在深度和广度上，都有了很大的发展，已成为系统与控制领域最热门的研究课题之一，取得了许多研究成果。

同时，也在与其他控制理论相互渗透，出现了许多新的最优控制方式，形成了更为实用的学科分支。

例如鲁棒最优控制[3]、随机最优控制[4]、分布参数系统的最优控制[5]、大系统的次优控制[6]、离散系统的最优控制及最优滑模变结构控制[7,8]等。

而对于非线性系统，其最优控制求解相当困难，需要求解非线性HJB方程或非线性两点边值问题，除简单情况外[9]，这两个问题都无法得到解析解。

因此，许多学者都致力于寻求近似的求解方法[10~13]，通过近似解得到近似的最优控，即次优控制。

1、非线性最优控制理论研究成果分类目前，较为流行的近似最优控制求解方法主要有以下几类[6][13]。

1）幂级数展开法：幂级数展开方法通过一个幂级数来构造控制律，得到序列形式的近似最优解，或者将系统中的非线性项以幂级数形式分解，或者通过引进一个临时变量并围绕它展开。

将上式代入HJB方程求得级数近似解，也可利用Adomian分解将非线性项进行分解，由此寻求非线性HJB方程级数的近似解。

2）Galerkin逐次逼近方法：由动态规划得到的一般性偏微分HJB方程，引入一个迭代过程来求解一般非线性HJB方程的一个近似解序列。

3）广义正交多项式级数展开法：其主要思想是将最优控制问题中的状态变量，控制输入，性能指标和各个参数分别用广义正交多项式展开，利用广义正交多项式的积分、乘积运算阵将描述系统的微分方程转化为一系列的代数方程。

然后，得到，T非奇异时由得到的控制律是一个多项式级数解。

该方法将最优控制问题转化为代数极值问题，从而避免了求解时变非线性Riccati方程。

4）有限差分和有限元方法：经典的有限差分和有限元方法可以用来近似求解非线性HJB方程。

近年来，这类方法用来近似求取非线性HJB方程的粘性解。

Vol.36,No.1ACTA AUTOMATICA SINICA January,2010An Optimal Control Scheme for a Class of Discrete-time Nonlinear Systems with Time Delays Using AdaptiveDynamic ProgrammingWEI Qing-Lai1ZHANG Hua-Guang2LIU De-Rong1ZHAO Yan3Abstract In this paper,an optimal control scheme for a class of nonlinear systems with time delays in both state and control variables with respect to a quadratic performance index function is proposed using a new iterative adaptive dynamic programming (ADP)algorithm.By introducing a delay matrix function,the explicit expression of the optimal control is obtained using the dynamic programming theory and the optimal control can iteratively be obtained using the adaptive critic technique.Convergence analysis is presented to prove that the performance index function can reach the optimum by the proposed method.Neural networks are used to approximate the performance index function,compute the optimal control policy,solve delay matrix function,and model the nonlinear system,respectively,for facilitating the implementation of the iterative ADP algorithm.Two examples are given to demonstrate the validity of the proposed optimal control scheme.Key words Adaptive dynamic programming(ADP),approximate dynamic programming,time delay,optimal control,nonlinear system,neural networksDOI10.3724/SP.J.1004.2010.00121The optimal control problem of nonlinear systems has always been a key focus in the controlfield in the last several decades.Coupled with this is the fact that noth-ing can happen instantaneously,as is so often presumed in many mathematical models.So strictly speaking,time delays exist in the most practical control systems.Time delays may result in degradation in the control efficiency even instability of the control systems.So there have been many studies on the control systems with time delay in various researchfields,such as electrical,chemical engineer-ing,and networked control[1−2].The optimal control prob-lem for the time-delay systems always attracts considerable attention of the researchers and many results have been obtained[3−5].In general,the optimal control for the time-delay systems is an infinite-dimensional control problem[3], which is very difficult to solve.So many analysis and appli-cations are limited to a very simple case:the linear systems with only state delays[6].For nonlinear case with state de-lays,the traditional method is to adopt fuzzy method and robust method,which transforms the nonlinear time-delay systems to linear systems[7].For systems with time delays both in states and controls,it is still an open problem[4−5]. The main difficulty lies in the formulation of the optimal controller which must use the information of the delayed control term so as to obtain an efficient control.This makes the analysis of the system much more difficult,and there is no method strictly facing this problem even in the linear cases.This motivates our research.Adaptive dynamic programming(ADP)is a powerful tool in solving optimal control problems[8−9]and has at-tached considerable attention from many researchers in re-cent years,such as[10−16].However,most of the results focus on the optimal control problems without delays.To Manuscript received September5,2008;accepted March3,2009Supported by National High Technology Research and Development Program of China(863Program)(2006AA04Z183),National Nat-ural Science Foundation of China(60621001,60534010,60572070, 60774048,60728307),and the Program for Changjiang Scholars and Innovative Research Groups of China(60728307,4031002)1.Key Laboratory of Complex Systems and Intelligence Sci-ence,Institute of Automation,Chinese Academy of Sciences,Beijing 100190,P.R.China2.School of Information Science and Engi-neering,Northeastern University,Shenyang110004,P.R.China3. Department of Automatic Control Engineering,Shenyang Institute of Engineering,Shenyang110136,P.R.China the best of our knowledge,there are no results discussing how to use ADP to solve the time-delay optimal control problems.In this paper,the time-delay optimal control problem is solved by the iterative ADP algorithm for the first time.By introducing a delay matrix function,the explicit expression of the optimal control function is ob-tained.The optimal control can iteratively be obtained us-ing the proposed iterative ADP algorithm which avoids the infinite-dimensional computation.Also,it is proved that the performance index function converges to the optimum using the proposed iterative ADP algorithm.This paper is organized as follows.Section1presents the preliminaries.In Section2,the time-delay optimal control scheme is proposed based on iterative ADP algorithm.In Section3,the neural network implementation for the con-trol scheme is discussed.In Section4,two examples are given to demonstrate the effectiveness of the proposed con-trol scheme.The conclusion is drawn in Section5.1PreliminariesBasically,we consider the following discrete-time affine nonlinear system with time delays in state and control vari-ables as follows:x(k+1)=f(x(k),x(k−σ))+g0(x(k),x(k−σ))u(k)+ g1(x(k),x(k−σ))u(k−τ)(1)with the initial condition given by x(s)=φ(s),s=−σ,−σ+1,···,0,where x(k)∈R n is the state vector, f:R n×R n→R n and g0,g1:R n×R n→R n×m are dif-ferentiable functions and the control u(k)∈R m.The state and control delaysσandτare both nonnegative integral numbers.Assume that f(x(k),x(k−σ))+g0(x(k),x(k−σ))u(k)+g1(x(k),x(k−σ))u(k−τ)is Lipschitz continuous on a setΩin R n containing the origin,and that system(1) is controllable in the sense that there exists a bounded con-trol onΩthat asymptotically stabilizes the system.In this paper,how to design an optimal state feedback controller for this class of delayed discrete-time systems is mainly dis-cussed.Therefore,it is desired tofind the optimal control u(x)satisfying u(x(k))=u(k)to minimize the generalized performance functional as follows:122ACTA AUTOMATICA SINICA Vol.36V(x(0),u)=∞k=0x T(k)Q0x(k)+2x T(k)Q1x(k−σ)+x T(k−σ)Q2x(k−σ)+u T(k)R0u(k)+2u T(k)R1u(k−τ)+u T(k−τ)R2u(k−τ)(2)whereQ0Q1Q T1Q2≥0,R0R1R T1R2>0,and l(x(k),x(k−σ),u(k),u(k−τ))=x T(k)Q0x(k)+2x T(k)Q1x(k−σ)+ x T(k−σ)Q2x(k−σ)+u T(k)R0u(i)+2u T(k)R1u(k−τ)+ u T(k−τ)R2u(k−τ)is the utility function.Let V∗(x)de-note the optimal performance index function which satisfiesV∗(x)=minuV(x,u)(3)According to the Bellman s optimal principle,we can get the following Hamilton-Jacobi-Bellman(HJB)equationV∗(x(k))=minu(k)x T(k)Q0x(k)+2x T(k)Q1x(k−σ)+x T(k−σ)Q2x(k−σ)+u T(k)R0u(k)+2u T(k)R1u(k−τ)+u T(k−τ)R2u(k−τ)+V∗(x(k+1))}(4) For the optimal control problem,the state feedback con-trol u(x)must not only stabilize the system onΩbut also guarantee that(2)isfinite,i.e.,u(x)must be admissible[17].Definition1.A control u(x)is defined to be admis-sible with respect to(3)onΩif u(x)is continuous on Ω,u(0)=0,u(x)stabilizes(1)onΩ,and∀x(0)∈Ω, V(x(0))isfinite.2Properties of the iterative ADP ap-proachSince the nonlinear delayed system(1)is infinite-dimensional[3]and the control variable u(k)couples with u(k−τ),it is nearly impossible to obtain the expression of the optimal control by solving the HJB equation(1).To overcome the difficulty,a new iterative algorithm is pro-posed in this paper.The following lemma is necessary to apply the algorithm.Lemma1.For the delayed nonlinear system(1)with respect to the performance index function(2),if there exists a control u(k)=0at time point k,then there exists a bounded matrix function M(k)that makesu(k−τ)=M(k)u(k)(5) hold for j=0,1,···,n.Proof.As u(k)and u(k−τj),j=0,1,···,n are bounded real vectors,can construct a function that sat-isfiesu(k−τ)=h(u(k))(6) where j=0,1,···,n.Then,using the method of undeter-mined coefficients,let M(u(k))satisfyh(u(k))=M(u(k))u(k)(7) Then,we can obtain M(u(k))expressed asM(u(k))=h(u(k))u T(k)u(k)u T(k)−1(8)whereu(k)u T(k)−1means the generalized inverse ma-trix ofu(k)u T(k).On the other side,u(k)and u(k−τ)are both bounded real vectors,then we have h(u(k))andu(k)u T(k)−1are bounded.So M(k)=M(u(k))is thesolution.According to Lemma1,the HJB equation becomesV∗(x(k))=x T(k)Q0x(k)+2x T(k)Q1x(k−σ)+x T(k−σ)Q2x(k−σ)+u∗T(k)R0u∗(k)+2u∗T(k)R1M∗(k)u∗(k)+u∗T(k)M∗T(k)R2M∗(k)u∗(k)+V∗(x(k+1))(9)where u∗(k)is the optimal control and u∗(k−τ)=M∗(k)u∗(k).2.1Derivation of the iterative ADP algorithmAccording to the Bellman s principle of optimality,wecan obtain the optimal control by differentiating the HJBequation(9)with respect to control u.Then,we can obtain the optimal control u∗(k)formu-lated asu∗(k)=−1R0+2R1M∗(k)+M∗T(k)R2M∗(k)−1×g0(x(k),x(k−σ))+g1(x(k),x(k−σ))M∗(k)T×∂V∗(x(k+1))∂x x(k+1)(10)In(10),the inverse of the termR0+2R1M∗(k)+M∗T(k)R2M∗(k)should existand a proof is presented in Appendix to guarantee theexistence of the inverse.From(10),the explicit optimal control expression u∗isobtained by solving the HJB equation(9).We can see thatthe optimal control u∗depends on M∗and V∗(x),whereV∗(x)is a solution to the HJB equation(9).While how tosolve the HJB equation is still open,there is currently nomethod for rigorously seeking for this performance indexfunction of this delayed optimal control problem.Further-more,the optimal delay matrix function M∗is also un-known which makes the optimal control u∗more difficultto obtain.So an iterative index i is introduced into theADP approach to obtain the optimal control iteratively.Firstly,for i=0,1,···,letu(i+1)(k−τ)=M(i)(k)u(i+1)(k)(11)where M(0)(k)=I and u(0)(k−τ)=M(0)(k)u(0)(k).Westart with initial performance index V(0)(x(k))=0,andthe control u(0)(k)can be computed as followsu(0)(x(k))=arg minuΓ0+V(0)(x(k+1))(12)whereΓ0=x T(k)Q0x(k)+2x T(k)Q1x(k−σ)+x T(k−σ)Q2x(k−σ)+u(0)T(k)R0u(0)(k)+2u(0)T(k)R1M(0)(k)u(0)(k)+u(0)T(k)M(0)T(k)R2M(0)(k)u(0)(k)Then,the performance index function is updated asV(1)(x(k))=Γ0+V(0)(x(k+1))(13)No.1WEI Qing-Lai et al.:An Optimal Comtrol Scheme for a Class of (123)Thus,for i =1,2,···,the iterative ADP can be used to implement the iteration betweenu (i )(x (k ))=arg min u Γ(i )+V (i )(x (k +1))=−12R 0+2R 1M (i −1)(k )+M (i −1)T (k )R 2M (i −1)(k ) −1 g 0(x (k ),x (k −σ))+g 1(x (k ),x (k −σ))M(i −1)(k )T ∂V (i )(x (k +1))∂xx (k +1)(14)whereΓ(i )=x T (k )Q 0x (k )+2x T (k )Q 1x (k −σ)+x T (k −σ)Q 2x (k −σ)+u (i )T (k )R 0u (i )(k )+2u (i )T (k )R 1M (i −1)(k )u (i )(k )+u (i )T (k )M (i −1)T (k )R 2M (i −1)(k )u (i )(k )andV (i +1)(x (k ))=Γ(i )+V (i )(x (k +1))(15)Then,the optimal control can be obtained iteratively.From (14)and (15),it can be seen that during the iteration pro-cess,the control actions for different control steps obey different control laws.After the iteration number of i ,the obtained control laws sequence is (u (0),u (1),···,u (i )).For the infinite-horizon problem,both the optimal performance index function and the optimal control law are unique.Therefore,it is necessary to show that the iterative per-formance index function V (i )(x (k ))will converge when the iteration number i →∞under the iterative control u (i )(k )and this will be proved in the following subsection.2.2Properties of the iterative ADP algorithm In this subsection,we focus on the proof of convergence of the iteration between (14)and (15),with the perfor-mance index V (i )(x (k ))→V ∗(x (k )),∀k .Lemma 2[17].Let ˜u(i )(k ),k =0,1,···be any se-quence of control,and u (i )(k )be expressed as (14).Define V (i +1)(x (k ))as (15)and Λ(i +1)(x (k ))asΛ(i +1)(x (k ))=x T (k )Q 0x (k )+2x T (k )Q 1x (k −σ)+x T (k −σ)Q 2x (k −σ)+˜u (i )T (k )R 0˜u (i )(k )+2˜u(i )T (k )R 1M (i −1)(k )˜u (i )(k )+˜u(i )T (k )M (i −1)T (k )R 2M (i −1)(k )˜u (i )(k )+Λ(i )(x (k +1))(16)If V (0)(x (k ))=Λ(0)(x (k ))=0,then V (i )(x (k ))≤Λ(i )(x (k )),∀i .In order to prove the convergence of the performance index function,the following theorem is also necessary.Theorem 1.Let the performance index function V (i )(x (k ))be defined by (15).If x (k )for system (1)is controllable,then there exists an upper bound Y such that 0≤V (i )(x (k ))≤Y ,∀i .Proof.As system (1)is Lipschitz,M (i )(k )is a bounded matrix for i =0,1,···.Define a delay matrix function¯M(k )which makes χTR 0+2R 1¯M(k )+¯M T (k )R 2¯M (k ) χ−χT (R 0+2R 1M (i )(k )+M (i )T (k )R 2M (i )(k )χ≥0(17)hold for ∀i ,where χis any nonzero m -dimensional vector.Let ¯u(k ),k =0,1···be any admissible control input.De-fine a new sequence P (i )(x (k ))as follows:P (i +1)(x (k ))=x T (k )Q 0x (k )+2x T (k )Q 1x (k −σ)+x T (k −σ)Q 2x (k −σ)+¯uT (k )R 0¯u (k )+2¯uT (k )R 1¯M (k )¯u (k )+¯uT (k )¯M T (k )R 2¯M (k )¯u (k )+P (i )(x (k +1))(18)where P (0)(x (k ))=V (0)(x (k ))=0and ¯u(k −τ)=¯M (k )¯u(k ).V (i )(x (k ))is updated by (15).Thus,we can obtainP (i +1)(x (k ))−P (i )(x (k ))=P (i )(x (k +1))−P (i −1)(x (k +1))...=P (1)(x (k +i ))−P (0)(x (k +i ))(19)Because P (0)(x (k +i ))=0,we haveP (i +1)(x (k ))=P (1)(x (k +i ))+P (i )(x (k ))+i j =0P (1)(x (k +j ))(20)According to (18),(20)can be rewritten asP(i +1)(x (k ))=i j =0Ξ(k +j )≤∞ j =0Ξ(k +j )(21)whereΞ(k +j )=x T (k +j )Q 0x (k +j )+2x T (k +j )Q 1x (k +j −σ)+x T (k +j −σ)Q 2x (k +j −σ)+¯uT (k +j )R 0¯u (k +j )+2¯uT (k +j )R 1¯M (k +j )¯u (k +j )+¯uT (k +j )¯M T (k +j )R 2¯M (k +j )u (k +j )Noting that the control input ¯u(k ),k =0,1,···is an ad-missible control,we can obtainP (i +1)(x (k ))≤∞ j =0P (1)(x (k +j ))≤Y,∀i (22)From Lemma 1,we haveV (i +1)(x (k ))≤P (i +1)(x (k ))≤Y,∀i(23)With Lemma 1and Theorem 1,the following main the-orem can be derived.Theorem 2.Define the performance index function V (i )(x (k ))as (15),with V (0)(x (k ))=0.If x (k )for system124ACTA AUTOMATICA SINICA Vol.36 (1)is controllable,then V(i)(x(k))is a nondecreasing se-quence that is V(i)(x(k))≤V(i+1)(x(k))and V(i)(x(k))isconvergent as i→∞.Proof.For the convenience of analysis,define a newsequenceΦ(i)(x(k))as follows:Φ(i+1)(x(k))=x T(k)Q0x(k)+2x T(k)Q1x(k−σ)+x T(k−σ)Q2x(k−σ)+u(i+1)T(k)R0u(i+1)(k)+2u(i+1)T(k)R1M(i)(k)u(i+1)(k)+u(i+1)T(k)M(i)T(k)R2M(i)(k)u(i+1)(k)+Φ(i)(x(k+1))(24)with u(i)(k)obtained by(14)andΦ0(x(k))=V0(x(k))=0.V(i)(x(k))is updated by(15).In the following part,we proveΦ(i)(x(k))≤V(i+1)(x(k))by mathematical induction.First,we prove it holds for i=0.Note thatV(1)(x(k))−Φ(0)(x(k))=x T(k)Q0x(k)+2x T(k)Q1x(k−σ)+x T(k−σ)Q2x(k−σ)≥0(25)Thus for i=0,we can getV(1)(x(k))≥Φ(0)(x(k))(26)Second,we assume it holds for i−1,i.e.,V(i)(x(k))−Φ(i−1)(x(k))≥0,∀x(k).Then,for i,from(15)and(24),we can obtainV(i+1)(x(k))−Φ(i)(x(k))=V(i)(x(k+1))−Φ(i−1)(x(k+1))≥0(27)i.e.,Φ(i)(x(k))≤V(i+1)(x(k))(28)Therefore,the mathematical induction proof is com-pleted.Moreover,from Lemma1,we know that V(i)(x(k))≤Φ(i)(x(k))and therefore we can obtainV(i)(x(k))≤Φ(i)(x(k))≤V(i+1)(x(k))(29)which proves that V(i)(x(k))is a nondecreasing sequencebounded by(23).Hence,we conclude that V(i)(x(k))is anondecreasing convergent sequence as i→∞.We note the obvious corollary.Corollary1.If Theorem2holds,then the delay matrixfunction M(i)(k)is a convergent sequence,as i→∞.According to Corollary1,we defineM(∞)(k)=limi→∞M(i)(k)(30)Next,we will prove that the performance index functionsequence V(i)(x(k))converges to V∗(x(k))as i→∞.AsV(i)(x(k))is a convergent sequence as i→∞,we defineV(∞)(x(k))=limi→∞V(i)(x(k))(31)Let¯u l be the l-th admissible control.Similar to theproof of Theorem1,we can construct the performance in-dex function sequence P(i)l(x)as follows:P(i+1)l(x(k))=x T(k)Q0x(k)+2x T(k)Q1x(k−σ)+x T(k−σ)Q2x(k−σ)+¯u T l(k)R0¯u l(k)+2¯u l(k)R1M(∞)(k)¯u l(k)+¯u l(k)M(∞)T(k)R2M(∞)(k)¯u l(k)+P(i)l(x(k+1))(32)with P(0)l(·)=0and¯u l(k)=M(∞)(k)¯u l(k−τ).Accordingto Theorem1,we haveP(i+1)l(x(k))=ij=0x T(k+j)Q0x(k+j)+2x T(k+j)Q1x(k+j−σ)+x T(k+j−σ)Q2x(k+j−σ)+¯u T l(k+j)R0¯u l(k+j)+2¯u T l(k+j)R1M(∞)(k+j)¯u l(k+j)+¯u T l(k+j)M(∞)T(k+j)R2×M(∞)(k+j)¯u l(k+j)(33)LetP(∞)l(x(k))=limi→∞P(i+1)l(x(k))(34)So,we haveP(i)l(x(k))≤P(∞)l(x(k))(35)Theorem3.Define P(∞)l(x(k))as in(34),and definethe performance index function V(i)(x(k))as in(15)withV(0)(·)=0.For any state vector x(k),define V∗(x(k))=min lP(∞)l(x(k))starting from x(k)for all admissiblecontrol sequences.Then,we can conclude that V∗(x(k))is the limit of the performance index function V(i)(x(k))asi→∞.Proof.For any l,there exists an upper bound Y l,suchthatP(i+1)l(x(k))≤P(∞)l(x(k))≤Y l(36)According to(23),for∀l,we haveV(∞)(x(k))≤P(∞)l(x(k))≤Y l(37)Since V∗(x(k))=min lP(∞)l(x(k)),for any >0,there exists an admissible control¯u K,where K is a nonneg-ative number such that the associated performance indexfunction satisfies P(∞)K(x(k))≤V∗(x(k))+ .According to(23),we have V(∞)(x(k))≤P(∞)l(x(k))for any l.Thus,we can obtain V(∞)(x(k))≤P(∞)K(x(k))≤V∗(x(k))+ .Noting that is chosen arbitrarily,we haveV(∞)(x(k))≤V∗(x(k))(38)On the other hand,since V(i)(x(k))is bounded for∀i,according to the definition of admissible control,the con-trol sequence associated with the performance index func-tion V(∞)(x(k))must be an admissible control,i.e.,thereNo.1WEI Qing-Lai et al.:An Optimal Comtrol Scheme for a Class of ···125exists an admissible control ¯u (i )N such that V(∞)(x (k ))=P (∞)N (x (k )).Combining with the definition V ∗(x (k ))=min l P (∞)l(x (k )),we can obtainV (∞)(x (k ))≥V ∗(x (k ))(39)Therefore,combining (38)and (39),we can conclude thatV (∞)(x (k ))=lim i →∞V (i )(x (k ))=V ∗(x (k ))(40)namely,V ∗(x (k ))is the limit of the performance indexfunction V (i )(x (k )),as i →∞. Based on Theorem 3,we will prove that the performance index function V ∗(x (k ))satisfies the principle of optimality,which shows that V (i )(x (k ))can reach the optimum as i →∞.Theorem 4.For any state vector x (k ),the “op-timal”performance index function V ∗(x (k ))satisfies V ∗(x (k ))=min u (k ){x T (k )Q 0x (k )+2x T (k )Q 1x (k −σ)+x T (k −σ)Q 2x (k −σ)+u T (k )R 0u (k )+2u T (k )R 1M (k )u (k )+u T (k )M (k )R 2M (k )u (k )+V ∗(x (k +1))},where u (k −τ)=M (k )u (k ).Proof.For any u (k )and i ,based on Bellman s optimal-ity principle,we haveV (i )(x (k ))≤Υ(i −1)+V (i −1)(x (k +1))(41)whereΥ(i −1)=x T (k )Q 0x (k )+2x T(k )Q 1x (k −σ)+x T(k −σ)Q 2x (k −σ)+u T(k )R 0u (k )+2u T(k )R 1M (i −1)(k )u (k )+u T(k )M(i −1)T(k )R 2M(i −1)(k )u (k )As V (i )(x (k ))≤V (i +1)(x (k ))≤V (∞)(x (k ))andV (∞)(x (k ))=V ∗(x (k )),we can obtainV (i )(x (k ))≤Υ(i −1)+V ∗(x (k +1))(42)If i →∞,then we haveV ∗(x (k ))≤Υ(∞)+V ∗(x (k +1))(43)Since u (k )in the above equation is chosen arbitrarily,the following equation holdsV ∗(x (k ))≤min u (k )Υ(∞)+V ∗(x (k +1)) (44)On the other hand,for any i ,the performance index func-tion satisfiesV (i )(x (k ))=Ω(i −1)+V (i −1)(x (k +1))(45)whereΩ(i −1)=x T (k )Q 0x (k )+2x T (k )Q 1x (k −σ)+x T (k −σ)Q 2x (k −σ)+u (i −1)T (k )R 0u (i −1)(k )+2u (i )T (k )R 1M (i −2)(k )u (i −1)T (k )+u(i −1)T(k )M(i −2)T(k )R 2M(i −2)(k )u(i −1)T(k )Combining with V (i )(x (k ))≤V ∗(x (k )),∀i ,we haveV ∗(x (k ))≥Ω(i −1)+V (i −1)(x (k +1))(46)Let i →∞,thenV ∗(x (k ))≥limi →∞Ω(i −1)+V(i −1)(x (k +1))≥min u (k )Ω(∞)+V ∗(x (k +1)) (47)Combining (44)with (47),we have V ∗(x (k ))=min u (k ){Ω(∞)+V ∗(x (k +1))}=x T (k )Q 0x (k )+2x T (k )Q 1x (k −σ)+x T (k −σ)Q 2x (k −σ)+u ∗T (k )R 0u ∗(k )+2u ∗T (k )R 1M (∞)(k )u ∗(k )+u ∗T (k )M (∞)T (k )R 2M (∞)(k )u ∗(k )+V ∗(x (k +1))(48)Thus,we have that u (i )(k )→u ∗(k )as i →∞so does u (i )(k −τ).On the other hand,we also have M (i )(k )→M (∞)(k )and u (i )(k −τ)=M (i −1)(k )u (i )(k ).Letting i →∞,we getu ∗(k −τ)=M (∞)(k )u ∗(k )(49)Therefore,we have M (∞)(k )=M ∗(k )and (48)can be written asV ∗(x (k ))=x T (k )Q 0x (k )+2x T (k )Q 1x (k −σ)+x T (k −σ)Q 2x (k −σ)+u ∗T (k )R 0u ∗(k )+2u ∗T (k )R 1M ∗(k )u ∗(k )+u ∗T (k )M ∗T (k )R 2M ∗(k )u ∗(k )+V ∗(x (k +1))(50)where u ∗(k −τ)=M ∗(k )u ∗(k ). Therefore,we can conclude that the performance index function V (i )(x (k ))converges to the optimum V ∗(x (k ))as i →∞.2.3Implementation of iterative ADP algorithm Given the above preparation,we may formulate the de-sired iterative ADP approach to nonlinear systems with delays as follows.Step 1.Give initial state x (s )=φ(s ),s =−σ,−σ+1,···,0,initial control u (ρ),ρ=0,1,···,k −1;give i max and computation accuracy ε.Step 2.Set the iterative step i =0,M (0)(k )=I ,and V (0)(·)=0.Step pute u (0)(k )by (12)and the performance index function V (1)(x (k ))by (13).Step 4.For the iterative step i ≥1,compute u (i )(k )by (14).Step pute the performance index function V (i )(x (k ))by (15).Step 6.IfV (i )(x (k ))−V (i −1)(x (k )) 2<ε(51)go to Step 9;otherwise,go to Step 7.Step 7.If i >i max ,go to Step 9;otherwise,compute M (i )(k )byM (i )(k )=u (i )(k −τ)u (i )T (k ) u (i )(k )u (i )T (k ) −1(52)126ACTA AUTOMATICA SINICA Vol.36Step 8.Set i =i +1and go to Step 4.Step 9.Stop.In (52)of the above algorithm,the term u (i )(k )u (i )T (k ) −1can be obtained by the Moore-Penrose pseudoinverse technique to compute the delay matrix function M (i )(k ).There are two other methods to compute M (i )(k ).One choice is to introduce a small zero-mean Gaussian noise with variances γ2denoted by δ(0,γ2)into the control u (k −τ)[18].The other choice is to use a neural network to approximate delay matrix function M (i )(k ).In this paper,we use the neural network approximation method and the details will be shown in the next section.3Neural network implementationIn the case of linear systems,the performance index func-tion is quadratic and the control policy is linear.In the nonlinear case,this is not necessarily true and therefore we use neural networks to approximate u (i )(k )and V (i )(x (k )).Assume the number of hidden layer neurons is denoted by l ,that the weight matrix between the input layer and hidden layer is denoted by V ,and that the weight matrix between the hidden layer and output layer is denoted by W .Then the output of three-layer neural network (NN)is represented byˆF(X ,V,W )=W T σ(V T X )(53)where σ(V T X )∈R l ,[σ(z )]i =e zi −e−z ie zi +e −z i,i =1,···,l,are the activation functions.The NN estimation error can be expressed byF (X )=F (X ,V ∗,W ∗)+ε(X )(54)where V ∗and W ∗are the ideal weight parameters,andε(X )is the reconstruction error.Here,there are four neural networks,which are critic network,model network,action network,and delay matrix function network (M network),respectively.All the neu-ral networks are chosen as three-layer feedforward network.The whole structure diagram is shown in Fig.1.The utility term in the figure denotes x T (k )Q 0x (k )+2x T (k )Q 1x (k −σ)+x T (k −σ)Q 2x (k −σ)+u T (k )R 0u (k )+2u T (k )R 1u (k −τ)+u T (k −τ)R 2u (k −τ).Fig.1The structure diagram of the algorithm3.1The model networkThe model network is to approximate the system dy-namic and it should be trained before the implementation of the iterative ADP algorithm.The update rule of the model network is adopted as the gradient decent method.The training process is simple and general.The details can be seen in [13,19]and it is omitted here.After the model network is trained,its weights are kept unchanged.3.2The M networkThe M network is to approximate the delay matrix func-tion M (k ).The output of the M network is denoted asˆu (k −τ)=W T M σ(V TM u (k ))(55)We define the error function of the model network ase M (k )=ˆu(k −τ)−u (k −τ)(56)Define the performance error measure asE M (k )=12e TM (k )e M (k )(57)Then,the gradient-based weight update rule for thecritic network can be described byw M (k +1)=w M (k )+∆w M (k )(58)∆w M (k )=αM−∂E M (k )∂w M (k )(59)where αM is the learning rate of the M network.3.3The critic networkThe critic network is used to approximate the perfor-mance index function V (i )(x (k )).The output of the critic network is denoted asˆV (i )(x (k ))=W T ci σ(V T ci z (k ))(60)The target function can be written asV (i +1)(x (k ))=Γ(i )+ˆV(i )(x (k +1))(61)Then,we define the error function for the critic networkase ci (k )=ˆV(i +1)(x (k ))−V (i +1)(x (k ))(62)And the objective function to be minimized in the criticnetwork isE ci (k )=12e 2ci (k )(63)So the gradient-based weight update rule for the critic network is given byw c (i +1)(k )=w ci (k )+∆w ci (k )(64)∆w ci (k )=αc−∂E ci (k )∂w ci (k )(65)∂E ci (k )∂w ci (k )=∂E ci (k )∂ˆV(i )(x (k ))∂ˆV(i )(x (k ))∂w ci (k )(66)where αc >0is the learning rate of critic network andw c (k )is the weight vector in the critic network.No.1WEI Qing-Lai et al.:An Optimal Comtrol Scheme for a Class of ···1273.4The action networkIn the action network,the state x (k )is used as input to create the optimal control as the output of the network.The output can be formulated asˆu(i )(k )=W T ai σ(V Tai x (k ))(67)The target of the output of the action network is givenby (14).So,we can define the output error of the action network ase ai (k )=ˆu (i )(k )−u (i )(k )(68)where u (i )(k )is the target function which can be describedbyu (i )(k )=−12R 0+2R 1M (i −1)(k )+M (i −1)T (k )R 2M (i −1)(k ) −1×g 0(x (k ),x (k −σ))+g 1(x (k ),x (k −σ))M(i −1)(k )T ∂ˆV (i )(x (k +1))∂xx (k +1)As u (i )(k −τ)=M (i −1)(k )u (i )(k ),we have ∂uu (i )(k −τ)∂uu (i )(k )=M (i −1)(k ).Then,according to (55),M (i −1)(k )can be ex-pressed asM (i −1)ij (k )=V T Mi 1− σ(V TM u (k ))2iW Mj (69)for i,j =1,2,···,m .M (i −1)ij(k )denotes the element at row i ,column j of matrix M (i −1)(k );V Mi and W Mj mean the column i and column j of the weight matrices V M and W M ,respectively; σ(V TMu (k )) i is the i th element of the vector σ(V TMu (k )).The weighs in the action network are updated to mini-mize the following performance error measure:E ai (k )=12e Tai (k )e ai (k )(70)The weights updating algorithm is similar to the one for the critic network.By the gradient descent rule,we can obtainw a (i +1)(k )=w ai (k )+∆w ai (k )(71)∆w ai (k )=βa−∂E ai (k )ai(72)∂E ai (k )∂w ai (k )=∂E ai (k )∂e ai (k )∂e ai (k )∂uu (i )(k )∂uu (i )(k )∂w ai (k )(73)where βa >0is the learning rate of then action network.4SimulationIn this section,two examples are provided to demon-strate the effectiveness of the control scheme proposed inthis paper.4.1Optimal control for state delayed system For the first example,the nonlinear system is a modifi-cation of Example 1in [13],which introduces state delays into the system.Consider the following affine nonlinear system:x (k +1)=f (x (k ),x (k −σ))+g (x (k ),x (k −σ))u (k )(74)where x (k )=x 1(k )x 2(k )T ,u (k )=u 1(k )u 2(k ) T ,and f (x (k ),x (k −σ))= x 1(k )exp(x 32(k ))x 2(k −2)x 32(k )x 1(k −2) ,g (x (k ),x (k −σ))= −0.200−0.2.The time delay in thestate is σ=2and the initial condition is x (k )=[1−1]T for −2≤k ≤0.The performance index function is defined as (2),where Q 0=Q 2=R 0=I and Q 1=R 1=R 2=0.We implement the algorithm at the time instant k =5.We choose three-layer neural networks as the critic network,the action network,and the model network with the struc-tures 4-10-2,2-10-1,and 6-10-2,respectively.The initial weights of the action network,critic network,and model network are all set to be random in [−0.5,0.5].It should be mentioned that the model network should be trained first.For the given initial state,we train the model net-work for 3000steps under the learning rate αm =0.05.After the training of the model network,the weights keep unchanged.Then,the critic network and the action net-work are trained for 3000steps so that the given accuracy ε=10−6is reached.In the training process,the learning rate βa =αc =0.05.The convergence curve of the perfor-mance index function is shown in Fig.2.Then,we apply the optimal control to the system for T f =30time steps and obtain the following results.The state trajectories are given as Fig.3and the corresponding control curves are given as Fig.4.Fig.2The convergence of performance indexfunctionFig.3The state variable trajectories。