删除或更新信息,请邮件至freekaoyan#163.com(#换成@)

一种求解控制系统可达集的网格投影法

本站小编 哈尔滨工业大学/2019-10-24

一种求解控制系统可达集的网格投影法

邵立珍1,张扬帆1,胡广大2

(1. 北京科技大学 自动化学院,北京100083; 2. 上海大学 理学院,上海200444)



摘要:

控制系统在某一时刻的前向可达集是指从初始状态出发,在该时刻能够达到的状态的集合.为了求解非线性控制系统的前向可达集,首先通过常微分方程数值方法对连续控制系统进行离散,转化为离散系统;其次由于连续控制系统的可达集可用其相应的离散系统的可达集来近似,针对离散控制系统,提出了基于最优化技术的网格投影法近似描述可达集,该方法布置了均匀分布的网格点,并将网格点向可达集边界投影,每个投影问题都对应一个最优化问题,通过求解这些优化问题得到可达集的近似描述;进一步地,理论分析证明了该方法布置的网格点间隔越小得到的可达集近似误差越小;最后通过数值仿真验证了该方法的有效性,并将其与文献中已有的DFOG(Distance fields on grids)方法进行对比.研究表明:网格投影法可以有效地处理具有非凸可达集的控制系统;相较于DFOG方法,该方法能够得到分布均匀的边界点,且求解的优化问题的数量少,计算时间短.

关键词:  控制系统  非凸可达集  数值方法  网格投影  优化

DOI:10.11918/j.issn.0367-6234.201710120

分类号:TP273

文献标识码:A

基金项目:国家自然科学基金 (11371053); 北京市自然科学基金 (4152034); 中央高校基本科研业务费专项资金(FRF-BD-16-005A)



A grid projection method for reachable sets of control systems

SHAO Lizhen1,ZHANG Yangfan1,HU Guangda2

(1. School of Automation,University of Science and Technology Beijing, Beijing 100083, China; 2. School of Science, Shanghai University, Shanghai 200444, China)

Abstract:

The forward reachable set of a control system at a certain time is the set of states that can be reached at this time from the initial state. In order to solve the forward reachable set of a nonlinear control system, numerical methods for ordinary differential equations was used to discretize the continuous control system. Since the reachable set of a continuous system can be approximated by its discrete counterpart, a grid projection method based on optimization technology was proposed to approximate the reachable set. In this method, uniformly distributed grid points were placed and projected to the boundary of reachable set for a discrete control; each projection problem corresponded to an optimization problem; by solving these problems the approximated description of the reachable set was obtained. Theoretical analysis proves that the smaller the grid interval, the smaller the approximation error. Finally, numerical results verified the effectiveness of the proposed method, which were compared with those of the DFOG(distance fields on grids) method in literature. The research shows that the grid projection method can effectively handle control systems with non-convex reachable sets. Compared with the DFOG method, the method gets uniform distributed boundary points, solves less optimization problems and uses less computation time.

Key words:  control system  non-convex reachable set  numerical method  grid projection  optimization


邵立珍, 张扬帆, 胡广大. 一种求解控制系统可达集的网格投影法[J]. 哈尔滨工业大学学报, 2018, 50(10): 88-94. DOI: 10.11918/j.issn.0367-6234.201710120.
SHAO Lizhen, ZHANG Yangfan, HU Guangda. A grid projection method for reachable sets of control systems[J]. Journal of Harbin Institute of Technology, 2018, 50(10): 88-94. DOI: 10.11918/j.issn.0367-6234.201710120.
基金项目 国家自然科学基金(11371053);北京市自然科学基金(4152034);中央高校基本科研业务费专项资金(FRF-BD-16-005A) 作者简介 邵立珍(1976—),女,副教授;
胡广大(1962—),男,教授,博士生导师 通信作者 邵立珍,lshao@ustb.edu.cn 文章历史 收稿日期: 2017-10-24



Contents            -->Abstract            Full text            Figures/Tables            PDF


一种求解控制系统可达集的网格投影法
邵立珍1, 张扬帆1, 胡广大2    
1. 北京科技大学 自动化学院,北京 100083;
2. 上海大学 理学院,上海 200444

收稿日期: 2017-10-24
基金项目: 国家自然科学基金(11371053);北京市自然科学基金(4152034);中央高校基本科研业务费专项资金(FRF-BD-16-005A)
作者简介: 邵立珍(1976—),女,副教授;
胡广大(1962—),男,教授,博士生导师
通信作者: 邵立珍,lshao@ustb.edu.cn


摘要: 控制系统在某一时刻的前向可达集是指从初始状态出发,在该时刻能够达到的状态的集合.为了求解非线性控制系统的前向可达集,首先通过常微分方程数值方法对连续控制系统进行离散,转化为离散系统;其次由于连续控制系统的可达集可用其相应的离散系统的可达集来近似,针对离散控制系统,提出了基于最优化技术的网格投影法近似描述可达集,该方法布置了均匀分布的网格点,并将网格点向可达集边界投影,每个投影问题都对应一个最优化问题,通过求解这些优化问题得到可达集的近似描述;进一步地,理论分析证明了该方法布置的网格点间隔越小得到的可达集近似误差越小;最后通过数值仿真验证了该方法的有效性,并将其与文献中已有的DFOG(Distance fields on grids)方法进行对比.研究表明:网格投影法可以有效地处理具有非凸可达集的控制系统;相较于DFOG方法,该方法能够得到分布均匀的边界点,且求解的优化问题的数量少,计算时间短.
关键词: 控制系统    非凸可达集    数值方法    网格投影    优化    
A grid projection method for reachable sets of control systems
SHAO Lizhen1, ZHANG Yangfan1, HU Guangda2    
1. School of Automation, University of Science and Technology Beijing, Beijing 100083, China;
2. School of Science, Shanghai University, Shanghai 200444, China


Abstract: The forward reachable set of a control system at a certain time is the set of states that can be reached at this time from the initial state. In order to solve the forward reachable set of a nonlinear control system, numerical methods for ordinary differential equations was used to discretize the continuous control system. Since the reachable set of a continuous system can be approximated by its discrete counterpart, a grid projection method based on optimization technology was proposed to approximate the reachable set. In this method, uniformly distributed grid points were placed and projected to the boundary of reachable set for a discrete control; each projection problem corresponded to an optimization problem; by solving these problems the approximated description of the reachable set was obtained. Theoretical analysis proves that the smaller the grid interval, the smaller the approximation error. Finally, numerical results verified the effectiveness of the proposed method, which were compared with those of the DFOG(distance fields on grids) method in literature. The research shows that the grid projection method can effectively handle control systems with non-convex reachable sets. Compared with the DFOG method, the method gets uniform distributed boundary points, solves less optimization problems and uses less computation time.
Keywords: control system    non-convex reachable set    numerical method    grid projection    optimization    
控制系统的可达集有前向可达集和后向可达集之分,其中前向可达集为系统从初始状态出发,到终止时刻的所有状态的集合;后向可达集为能够到达终止时刻目标状态的所有状态的集合[1].本文仅研究前向可达集,文中所提到的可达集均指前向可达集.确定可达集的边界对于分析系统的动态性能、设计系统控制器具有重要的实际意义.对于非线性控制理论和动态系统中不定参数的定量研究而言,可达集的数值计算是关键的问题之一.可达集的研究可广泛地应用在众多的学科领域.例如,在飞行器控制领域,可以对飞行器轨道的可达集进行分析,实现良好的变轨过程[2-3]、确定飞行器着陆区域[4]、为驾驶员完成标准机动动作提供决策支持[5];通过对导弹运动轨迹的可达集研究,可以帮助确定制导导弹的目标命中区域[6];通过研究机器人运动轨迹的可达集,可以使其选取最短的行动路径并躲避障碍物[7].

研究表明, 线性控制系统的可达集是凸集[8],而非线性系统的可达集则可能是凸的或非凸的两种.目前针对线性系统,已经提出了很多种近似可达集的理论及算法.文献[9-10]通过构造若干接触可达集边界的外部或内部椭球,来得到可达集的近似估计;文献[11-12]通过将可达集等价描述为有限个最优控制问题,通过直接求解这些最优控制问题来得到可达集的近似;文献[13]将可达集描述为一个Hamilton-Jacobi偏微分方程,并利用水平集方法近似可达集.但对于非线性系统,尤其是具有非凸可达集的系统,计算其可达集目前受到了以下几方面的阻碍:可达集几何复杂度会随着时间间隔增加而增加,非凸性所带来的计算问题以及非线性系统造成的计算不稳定性.为了求解非线性控制系统的可达集,文献[14]分析了非线性离散系统的可达集为凸的充分和必要条件;文献[15]在二维可达集内部布置了若干点,并迭代更新,使这些点的分段线性轮廓线所围成的区域最大,从而近似可达集;文献[16]指出连续系统的可达集可用其相应的离散系统的可达集来近似,近似精度取决于离散化方法和离散化步长;针对离散非线性系统,文献[17]提出了DFOG(distance fields on grids)算法,通过在包含可达集的区域内布置均匀分布的网格点,除去以外部网格点为圆心、该点到可达集的最短距离为半径的开球集,剩余集合即为可达集的近似集合.然而上述方法在可达集近似精度和计算效率方面都有待进一步的提高.

本文在已有的研究中,提出了基于凸优化技术的外部投影法来求解线性离散控制系统的可达集,但该方法不适用于具有非凸可达集的非线性控制系统,且没有给出近似结果的误差分析.本文将继续研究求解离散控制系统的可达集的方法,重点针对具有非凸可达集的非线性控制系统.文中首先对非线性系统及其可达集进行描述,然后介绍离散化后的系统及其相应的离散可达集,进一步地提出了近似离散系统可达集的网格投影法,并给出了近似误差的理论分析,最后将该方法与文献[17]中的DFOG方法进行数值仿真对比,给出了仿真和分析结果.

1 控制系统及其可达集考虑非线性动态系统为

$\left\{ \begin{array}{l}\mathit{\boldsymbol{\dot x}}\left( t \right) = f\left( {t,\mathit{\boldsymbol{x}}\left( t \right),\mathit{\boldsymbol{u}}\left( t \right)} \right),\\\mathit{\boldsymbol{x}}\left( {{t_0}} \right) = {\mathit{\boldsymbol{x}}_0},\\\mathit{\boldsymbol{u}}\left( t \right) \in \mathit{\boldsymbol{U}}.\end{array} \right.$ (1)

式中:xRn为系统的状态向量;x0Rn为系统的初始状态;U? Rm是一个非空凸紧集;t∈[t0, tf],其中t0tf分别为动态系统的初始时刻和终点时刻,且t0 < tf.

定义1????系统(1)从初始条件(t0, x0)出发,在tf时刻的前向可达集可定义如下

$\Re \left( {{t_{\rm{f}}},{t_0},{\mathit{\boldsymbol{x}}_0}} \right): = \left\{ {\mathit{\boldsymbol{y}} \in {{\bf{R}}^n}:\exists \mathit{\boldsymbol{u}}\left( \cdot \right){\rm{s}}.\;{\rm{t}}.\;\mathit{\boldsymbol{x}}\left( {{t_{\rm{f}}}} \right) \in \mathit{\boldsymbol{y}}} \right\}.$ (2)

对于连续系统,可以用数值方法如欧拉法、龙格-库塔法、线性多步法等[18]对其离散化得到离散系统.以欧拉法为例, 在[t0, tf]时间内,取等间隔时间ti=t0+ihh=(tft0)/Ni=0, …, N.假定控制函数u(t)在[ti, ti+1]间隔内恒定为$\mathit{\boldsymbol{\tilde u}}\left( {{t_i}} \right)$,则采用欧拉法对系统(1)离散化, 可得离散系统为

$\left\{ \begin{array}{l}\mathit{\boldsymbol{\tilde x}}\left( {{t_{i + 1}}} \right) = \mathit{\boldsymbol{\tilde x}}\left( {{t_i}} \right) + hf\left( {{t_i},\mathit{\boldsymbol{\tilde x}}\left( {{t_i}} \right),\mathit{\boldsymbol{\tilde u}}\left( {{t_i}} \right)} \right),\\\mathit{\boldsymbol{\tilde u}}\left( {{t_i}} \right) \in \mathit{\boldsymbol{U}},\\\mathit{\boldsymbol{\tilde x}}\left( {{t_0}} \right) = {\mathit{\boldsymbol{x}}_0},\end{array} \right.$ (3)

式中i=0, …, N-1.

定义2????对于离散系统(3),从初始条件(t0, x0)出发,在tN时刻的可达集可定义如下

${\tilde \Re } \left( {{t_N},{t_0},{\mathit{\boldsymbol{x}}_0}} \right): = \left\{ {\mathit{\boldsymbol{y}} \in {{\bf{R}}^n}:\exists \mathit{\boldsymbol{\tilde u}}\left( \cdot \right){\rm{s}}.\;{\rm{t}}.\;\mathit{\boldsymbol{\tilde x}}\left( {{t_N}} \right) \in \mathit{\boldsymbol{y}}} \right\}.$ (4)

用离散系统的可达集(4)对连续系统的可达集(2)进行近似,其近似精度可以利用豪斯多夫距离(Hausdorff distance)来衡量.

定义3????设集合C, D ? Rn为两个非空紧集,则CD之间的豪斯多夫距离为

${d_{\rm{H}}}\left( {\mathit{\boldsymbol{C}},\mathit{\boldsymbol{D}}} \right): = \max \left\{ {d\left( {\mathit{\boldsymbol{C}},\mathit{\boldsymbol{D}}} \right),d\left( {\mathit{\boldsymbol{D}},\mathit{\boldsymbol{C}}} \right)} \right\},$

式中$d(\mathit{\boldsymbol{C}},\mathit{\boldsymbol{D}}) = \mathop {{\rm{max}}}\limits_{\mathit{\boldsymbol{c}} \in \mathit{\boldsymbol{C}}} \mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{d}} \in \mathit{\boldsymbol{D}}} {\left\| {\mathit{\boldsymbol{c}} - \mathit{\boldsymbol{d}}} \right\|_2}$.

文献[16]指出,连续系统的可达集$\Re $(tf, t0, x0)可以用其离散系统的可达集${\tilde \Re }$ (tN, t0, x0)来近似,近似精度取决于离散步长h.

定理1????设系统(1)中的f满足利普希茨(Lipschitz)条件,取正整数N为离散步数,h=(tft0)/N为离散步长.利用欧拉方法对系统(1)进行离散化处理得到系统(3),则存在一常数M>0使得连续系统可达集$\Re $(tf, t0, x0)与离散系统可达集${\tilde \Re }$(tN, t0, x0)有如下关系:

${d_{\rm{H}}}\left( {\Re \left( {{t_{\rm{f}}},{t_0},{\mathit{\boldsymbol{x}}_0}} \right),\tilde \Re \left( {{t_N},{t_0},{\mathit{\boldsymbol{x}}_0}} \right)} \right) \le Mh.$

本文用离散控制系统可达集近似连续控制系统可达集,重点研究离散控制系统可达集的数值求解方法.根据定义2,离散控制系统的可达集问题可以归结为如下求集合Y的问题:

$\begin{array}{l}{\rm{Find}}\;\mathit{\boldsymbol{Y}}: = \left\{ {\mathit{\boldsymbol{y}} \in {{\bf{R}}^n}:\mathit{\boldsymbol{y}} = \mathit{\boldsymbol{\tilde x}}\left( {{t_N}} \right)} \right\},\\{\rm{s}}.\;{\rm{t}}.\;\;\mathit{\boldsymbol{\tilde x}}\left( {{t_0}} \right) = {\mathit{\boldsymbol{x}}_0},\\\;\;\;\;\;\;\mathit{\boldsymbol{\tilde x}}\left( {{t_{i + 1}}} \right) = \mathit{\boldsymbol{\tilde x}}\left( {{t_i}} \right) + hf\left( {{t_i},\mathit{\boldsymbol{\tilde x}}\left( {{t_i}} \right),\mathit{\boldsymbol{\tilde u}}\left( {{t_i}} \right)} \right),\\\;\;\;\;\;\;\mathit{\boldsymbol{\tilde u}}\left( {{t_i}} \right) \in \mathit{\boldsymbol{U}},\\\;\;\;\;\;\;i = 0, \cdots ,N - 1.\end{array}$ (5)

进一步地,将式(5)中的每个约束都转化为gi(x)≤0的形式,则上述问题描述为

$\begin{array}{l}{\rm{Find}}\;\mathit{\boldsymbol{Y}}: = \left\{ {\mathit{\boldsymbol{Px}}:\mathit{\boldsymbol{x}} \in \mathit{\boldsymbol{X}}} \right\},\\{\rm{s}}.\;{\rm{t}}.\;\;\mathit{\boldsymbol{X}} = \left\{ {\mathit{\boldsymbol{x}} \in {{\bf{R}}^n}:{g_i}\left( \mathit{\boldsymbol{x}} \right) \le 0,i = 1, \cdots ,l} \right\}.\end{array}$ (6)

式中:x=(${\mathit{\boldsymbol{\tilde u}}}$(t0)T, …, ${\mathit{\boldsymbol{\tilde u}}}$(tN-1)T, ${\mathit{\boldsymbol{\tilde x}}}$(t0)T, …, ${\mathit{\boldsymbol{\tilde x}}}$(tN-1)T, ${\mathit{\boldsymbol{\tilde x}}}$(tN)T)T为变量;P为一个最后n列为单位阵而其他元素为0的矩阵.

2 网格投影法网格投影法是一种求解式(6)中Y的方法.该方法从集合某个内部点出发,沿着坐标轴的正负方向,逐渐构造一系列等间距的网格,并将这些网格点沿着各个坐标轴方向向集合边界进行投影以得到边界点.

2.1 网格点生成与投影点计算设集合Y?Rn为一非空紧集, 用intY表示集合Y内部,而bdY表示集合Y边界.取任意一点a∈intY.将点a沿坐标轴i的正、反方向向集合Y的边界进行投影,可以分别通过解决以下最优化问题来求得相应的投影边界点:

$P_{\max }^{i + }\left( \mathit{\boldsymbol{a}} \right)\max {d_{i + }}\;{\rm{s}}.\;{\rm{t}}.\;\mathit{\boldsymbol{a}} + {d_{i + }}{\mathit{\boldsymbol{e}}_i} = \mathit{\boldsymbol{y}},\mathit{\boldsymbol{y}} \in \mathit{\boldsymbol{Y}};$

$P_{\max }^{i - }\left( \mathit{\boldsymbol{a}} \right)\max {d_{i - }}\;{\rm{s}}.\;{\rm{t}}.\;\mathit{\boldsymbol{a}} - {d_{i - }}{\mathit{\boldsymbol{e}}_i} = \mathit{\boldsymbol{y}},\mathit{\boldsymbol{y}} \in \mathit{\boldsymbol{Y}}.$

式中ei为第i个元素值为1,其余元素均为0的n维向量.设di+*d*i分别为求解Pmaxi+(a)、Pmaxi(a)得到的最优目标函数值,则从a点出发沿着坐标轴i正、负方向向Y投影得到的边界点分别为:

${\mathit{\boldsymbol{p}}_{i + }} = \mathit{\boldsymbol{a}} + d_{i + }^ * {\mathit{\boldsymbol{e}}_i},{\mathit{\boldsymbol{p}}_{i - }} = \mathit{\boldsymbol{a}} - d_{i - }^ * {\mathit{\boldsymbol{e}}_i}.$

a点为起点,在pi+pi之间,可沿坐标轴i的正、反方向布置间隔为ds>0的均匀分布的参考点.进一步将这些参考点沿着其它的坐标轴上投影,生成新的参考点,所有的参考点形成了网格点.

图 1展示了投影及网格点的生成过程.图中二维平面内的闭合实线围成的区域为所求集合Y,内部两条带箭头实线的交点为任意选取的初始点a,箭头分别表示4个投影方向,边界线上有4个投影边界点p1-p1+p2-p2+,线上的每个点表示由a生成的等间隔的网格点.

Figure 1
图 1 投影及网格点生成 Figure 1 Projection and generation of grid points


2.2 算法描述网格投影法的算法流程如图 2所示.

Figure 2
图 2 网格投影法的算法流程 Figure 2 Algorithm flow of the grid projection method


算法1的输出结果BY的边界点集,将B中的相邻点依次连接,围成的内部区域记作clB,clB即为Y的近似集.算法1可以有效地解决具有凸可达集的控制系统问题,如线性系统的可达集.但是对于非凸可达集,如图 3所示,二维平面内的实现围成的闭合区域为集合Ya为初始点,闭合区域的边界上p1-, p1+为投影点,两投影点间的圆点为参考点,出现了参考点位于Y外部的情况.这是由于求解Pmaxi+(a)、Pmaxi(a),内部点的投影会跨过可达集非凸部分的边界线,而得到最远投影距离的投影点p1+,致使生成的内部网格点实际落在了集合Y的外部,如图中标注的6个外部参考点.若对外部的参考点继续按照算法1的流程生成网格点,将无法求得正确的边界点,且导致算法无法停止.

Figure 3
图 3 错误生成外部网格点 Figure 3 Erroneous generation of external grid points


鉴于实际问题中,大多数可达集Y具有如下特点:对于任意一点y?Y,将该点沿每个坐标轴的正、负方向Y做投影,至少有一个方向的投影不存在.因此对于这种类型的非凸可达集,可通过在网格投影法中添加额外的判断条件以区分外部投影点.对某个参考点a′沿所有坐标轴求解Pmaxi+(a′)、Pmaxi(a′),若出现沿某个方向的解不存在时,则认为a′?Y.这时,可计算如下的优化问题得到a′沿坐标轴在Y上的投影,进而修正投影边界点.

$P_{\min }^{i + }\left( {\mathit{\boldsymbol{a'}}} \right)\min {d_{i + }}\;{\rm{s}}.\;{\rm{t}}.\;\mathit{\boldsymbol{a'}} + {d_{i + }}{\mathit{\boldsymbol{e}}_i} = \mathit{\boldsymbol{y}},\mathit{\boldsymbol{y}} \in \mathit{\boldsymbol{Y}};$

$P_{\max }^{i - }\left( {\mathit{\boldsymbol{a'}}} \right)\min {d_{i - }}\;{\rm{s}}.\;{\rm{t}}.\;\mathit{\boldsymbol{a'}} - {d_{i - }}{\mathit{\boldsymbol{e}}_i} = \mathit{\boldsymbol{y}},\mathit{\boldsymbol{y}} \in \mathit{\boldsymbol{Y}}.$

具体修正过程如下:当a′沿坐标轴i投影时,若Pmaxi+(a′)的解存在而Pmaxi(a′)的解不存在,则求解Pmini+(a′)得到投影点pi+= a′+di+*ei,并在pi+pi+之间布置网格点;同理,若Pmaxi+(a′)的解不存在而Pmaxi(a′)的解存在,则求解Pmini(a′),得到投影点pi= a′-di*ei,并在pipi之间布置网格点;对该点的其他坐标轴方向,需要求解Pmini+(a′)、Pmini(a′)得到pi+pi来补全边界投影点, 且这两个边界投影点之间的网格点均被判断为外部网格点.修正后的网格投影法算法如图 4所示.

Figure 4
图 4 修正后的网格投影法算法流程 Figure 4 Algorithm flow of the modified grid projection method


算法2在处理外部参考点时如图 5所示,图 5中二维平面内的实线围成的闭合区域为所求集合Ya为初始点,图中有6个外部网格点位于集合Y外部,它们沿纵坐标轴负方向的投影点不存在,其中带有箭头实线的交点为计算投影的网格点g,闭合区域的边界上p1+p1-为对g求解Pmin1+(g)、Pmin1-(g)得到的投影点,p2+p2+为对g求解Pmax2+(g)、Pmin2+(g)得到的投影点,圆点为网格点.

Figure 5
图 5 修正算法处理的外部网格点 Figure 5 Revised algorithm to deal with external grid points


2.3 算法近似误差分析本文把由算法2得到的边界点集B相邻点依次连接,所围成的区域记作clB.clB与真实集合Y间存在如下关系.

定理2????对于一个非空闭集Y?Rn,假定选取参考点间隔为ds>0,且Y与内部参考点集Gint?Y满足dH(Y, Gint),利用算法2近似Y,记所得投影边界点集合为B,则YB满足如下关系:

${d_{\rm{H}}}\left( {\mathit{\boldsymbol{Y}},{\rm{cl}}\mathit{\boldsymbol{B}}} \right) \le \sqrt n {d_{\rm{s}}}.$

证明????根据定义3,有dH(Y, clB)= max{d(Y, clB), d(clB, Y)},其中:

$d\left( {\mathit{\boldsymbol{Y}},{\rm{cl}}\mathit{\boldsymbol{B}}} \right) = \mathop {\max }\limits_{\mathit{\boldsymbol{y}} \in \mathit{\boldsymbol{Y}}} d\left( {\mathit{\boldsymbol{y}},{\rm{cl}}\mathit{\boldsymbol{B}}} \right),$ (7)

$d\left( {{\rm{cl}}\mathit{\boldsymbol{B}},\mathit{\boldsymbol{Y}}} \right) = \mathop {\max }\limits_{\mathit{\boldsymbol{b}} \in {\rm{cl}}\mathit{\boldsymbol{B}}} d\left( {\mathit{\boldsymbol{b}},\mathit{\boldsymbol{Y}}} \right).$ (8)

对于式(7),取任意一点yY,若同时y∈clB,则d(y, clB)=0;若y?clB,由于离散点集Gint?clB和定理2的假定条件,有$d(\mathit{\boldsymbol{y}},{\rm{cl}}\mathit{\boldsymbol{B}}) \le d(\mathit{\boldsymbol{y}},{\mathit{\boldsymbol{G}}_{{\rm{int}}}}) = \mathop {{\rm{max}}}\limits_{\mathit{\boldsymbol{y}} \in \mathit{\boldsymbol{Y}}} \mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{g}} \in {\mathit{\boldsymbol{G}}_{{\rm{int}}}}} d(\mathit{\boldsymbol{y}},\mathit{\boldsymbol{g}}) \le \sqrt n {d_{\rm{s}}}$,式(7)得证.

对于式(8),取任意一点b∈clB,若bY,则d(b, Y)=0;若b?Y, 则有d(b, Y)≤d(b, Gint)≤max{d(F1, Gint), …, d(Fk, Gint)},其中Fi(i=1, …, k)为clB的外围边界面, 且有∪i=1kFi= bd{clB}.n维空间内每个Fk可由相邻的投影边界点bik(i=1, …, n)的凸包来表示.对于每个Fk及其对应的bik而言,假定b1kb2k为相邻的投影点,与它们对应的最近的网格点为g1kg2k,如果生成b1kb2k的两个投影方向相同,则d(g1k, g2k)=ds,否则g1kg2k为同一点.若对每个gik以其自身为起点沿gikbik方向生成新的参考点gik?Y,有d(gik, gik)=ds.将gikgik作为顶点可以构建一个超立方体,该立方体的最长对角线距离为$\sqrt n {d_{\rm{s}}}$.由立方体的构建方式可知bik位于该立方体的边界面上,Fk位于该立方体内.因此有$d({\mathit{\boldsymbol{F}}^k},{\mathit{\boldsymbol{G}}_{{\rm{int}}}}) = \mathop {{\rm{max}}}\limits_{\mathit{\boldsymbol{b}} \in {\mathit{\boldsymbol{F}}^k}} d(\mathit{\boldsymbol{b}},\mathit{\boldsymbol{g}}_i^k) \le \sqrt n {d_{\rm{s}}}$,则式(8)得证,因此定理2成立.

定理2????表明了网格投影法的近似误差与ds有关,ds越小,布置的网格点密度越大,得到的投影边界点数量越多,近似误差越小.

3 数值仿真本文将网格投影法和文献[17]中的DFOG方法进行对比.实验采用MATLAB 2016a仿真了两种算法,选取了同样的可达集内部网格点,利用求解器IPOPT[19]对投影时的非线性优化问题进行求解.将上述程序用于求解以下3个控制系统问题.对于所有的例子,均采用带有分段恒定控制量的欧拉法来离散化.

首先考虑两个二维系统的例子.

例1????考虑如下一个二维双线性系统:

$\left\{ \begin{array}{l}{{\dot x}_1}\left( t \right) = {\rm{ \mathsf{ π} }}{x_2}\left( t \right),\\{{\dot x}_2}\left( t \right) = - {\rm{ \mathsf{ π} }}u\left( t \right){x_1}\left( t \right),\\{x_1}\left( 0 \right) = - 1,\\{x_2}\left( 0 \right) = 0,\\u\left( t \right) \in U = \left[ {0,1} \right],\end{array} \right.$

式中初始时刻t0=0,终点时刻tf=1,取欧拉法的步长h=1/30.

图 6(a)、(b)分别显示了网格投影法和DFOG法的结果.图中圆点为网格点,星点为求得的可达集边界点,图 6(b)中的大圆(不含边界线)展示了DFOG方法计算得到的开球集,DFOG法得到的可达集为圆圈集围成的内部区域.可以看出,两种方法的得到的可达集类似,但网格投影法求得的边界点疏密更为均匀.

Figure 6
图 6 例1的近似可达集 Figure 6 Approximated reachable set for example 1


例2????考虑如下一个二维非线性系统:

$\left\{ \begin{array}{l}{{\dot x}_1}\left( t \right) = x_2^2\left( t \right) + {x_1}\left( t \right) + u\left( t \right),\\{{\dot x}_2}\left( t \right) = {x_2}\left( t \right) + 5\sin \left( {{x_1}\left( t \right)} \right) + u\left( t \right),\\{x_1}\left( 0 \right) = 0,\\{x_2}\left( 0 \right) = 0,\\\left| {u\left( t \right)} \right| \le 1,\end{array} \right.$

式中:初始时刻t0=0,终点时刻tf=1,取欧拉法的步长h=0.01.

图 7(a)、(b)分别显示了网格投影法和DFOG法的结果.由于本系统的可达集在坐标轴上跨度较大,DFOG算法在可达集外部布置的网格点较多;此外也有部分网格点的优化问题没有求得全局最优值,使得相应的开球集半径过大,导致所求可达集比实际可达集小.而网格投影法并不需要求解如此多的优化问题,耗时较短.

Figure 7
图 7 例2的近似可达集 Figure 7 Approximated reachable set for example 2


例3????考虑如下的一个三维线性时不变系统.系统矩阵和控制矩阵分别为:

$\mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}}{ - 1}&1&0\\0&{ - 1}&0\\0&0&{ - 2}\end{array}} \right],\mathit{\boldsymbol{B = }}\left[ {\begin{array}{*{20}{c}}0&0&0\\0&2&0\\0&0&1\end{array}} \right].$

式中:控制约束为U={uR3:‖u2≤1},初始状态x0=(0, 0, 0)T,初始时刻t0=0,终点时刻tf=1,取欧拉步长h=0.02.

图 8(a)、(b)分别显示了网格投影法和DFOG法求得的边界点凸包.图中可以看出DFOG法左侧某些边界点落在了网格投影法得到的可达集内部,表明了DFOG法中某些外部网格点的优化问题没有求得全局最优值,使计算得到的边界点位于真实可达集的内部.

Figure 8
图 8 例3的近似可达集 Figure 8 Approximated reachable set for example 3


表 1进一步总结了用网格投影法和DFOG法计算例1~3得到的结果,表中列出两种方法的计算时间,计算非线性优化问题的个数(NPs),最终得到的边界点的个数(BPs).如对于问题例1,表 1记录了网格投影法的耗时/s、NPs、BPs分别为22.81、54.00、46.00;而DFOG法的相应数据分别为28.97、128、72.这是由于DFOG法需要布置大量外部网格点,且每个网格点都需要求解一个优化问题;而网格投影法中,仅对于非凸集才需要布置外部网格点,且沿每个投影方向只需要计算一次投影.可以说明,在同样内部网格点的情况下,相较于DFOG方法,网格投影法所需的计算时间短,求解优化问题的个数少.进一步地,仿真结果图也显示了网格投影法得到的可达集边界点分布更加均匀.与DFOG法相比,得到的可达集描述更接近真实的可达集.

表 1
表 1 数值仿真结果 Table 1 Numerical simulation results 例子 算法 耗时/s NPs BPs

1 网格投影法 22.81 54 46

DFOG法 28.97 128 72

2 网格投影法 142.79 124 108

DFOG法 841.50 936 760

3 网格投影法 256.13 218 194

DFOG法 868.70 882 791



表 1 数值仿真结果 Table 1 Numerical simulation results


4 结论1) 网格投影法在可达集内部布置均匀分布的网格点,将网格点向可达集边界投影,每个投影点都是可达集的边界点,这些投影点围城的区域即为可达集的近似描述.

2) 网格投影法的近似误差与布置网格点的间隔相关,网格点间隔越小,则近似误差越小,需要求解优化问题个数越多.

3) 相较于DFOG算法,网格投影法具有近似误差小、边界点分布均匀、速度快的优点;当DFOG方法出现优化问题无法求得全局最优时,会使求得的可达集比真实的可达集小,实验结果显示本方法可靠性较高.

4) 提出了网格投影法直接近似离散系统的可达集.该方法简单,易于实施,可以有效地解决可达集为非凸的情况.


参考文献
[1] TRINH H, NAM P T, PATHIRANA P N, et al. On backwards and forwards reachable sets bounding for perturbed time-delay systems[J]. Applied Mathematics and Computation, 2015, 269: 664. DOI:10.1016/j.amc.2015.07.116


[2] 雪丹, 李俊峰, 宝音贺西. 平面脉冲作用下卫星轨道的可达范围研究[J]. 宇航学报, 2009, 30(1): 88.
XUE Dan, LI Junfeng, BAOYIN Hexi. Study on reachable domain for satellite trajectory with coplanar impulse applied[J]. Journal of Astronautics, 2009, 30(1): 88. DOI:10.3873/j.issn.1000-1328.2009.00.015


[3] LI Xuehua, HE Xingsuo, ZHONG Qinfang, et al. Reachable domain for satellite with two kinds of thrust[J]. Acta Astronautica, 2011, 68(11/12): 1860. DOI:10.1016/j.actaastro.2011.01.004


[4] 解永锋, 唐硕. 亚轨道飞行器再入可达域快速计算方法[J]. 飞行力学, 2011, 29(4): 72.
XIE Yongfeng, TANG Shuo. Rapid calculation of entry footprint of suborbital launch vehicles[J]. Flight Dynamics, 2011, 29(4): 72.


[5] 刘瑛, 杜光勋, 全权, 等. 基于Hamilton-Jacobi方程的飞行器机动动作可达集分析[J]. 自动化学报, 2016, 42(3): 347.
LIU Ying, DU Guangxun, QUAN Quan, et al. Reachability calculation for aircraft maneuver using Hamilton-Jacobi function[J]. Acta Automatica Sinica, 2016, 42(3): 347. DOI:10.16383/j.aas.2016.c140888


[6] 寇英信, 陈磊, 李战武, 等. 基于命中概率的制导炸弹可达域定量缩减方法[J]. 电光与控制, 2013, 20(7): 11.
KOU Yingxin, CHEN Lei, LI Zhanwu, et al. Guided bomb accessible region quantitative reduction method based on hit probability[J]. Electronics Optics and Control, 2013, 20(7): 11. DOI:10.3969/j.issn.1671-637X.2013.07.003


[7] GERDTS M, HENRION R, HOMBERG D, et al. Path planning and collision avoidance for robots[J]. Numerical Algebra Control and Optimization, 2012, 2(3): 437. DOI:10.3934/naco.2012.2.437


[8] 邵立珍, 赵方园, 胡广大. 一种求解线性控制系统可达集的数值方法[J]. 控制与决策, 2016, 32(3): 541.
SHAO Lizhen, ZHAO Fangyuan, HU Guangda. A Numercial method for reachable sets of linear control systems[J]. Control and Decision, 2016, 32(2): 541. DOI:10.13195/j.kzyjc.2016.0212


[9] KURZHANSKIY A, VARAIYA P. Ellipsoidal techniques for reachability analysis of discrete-time linear systems[J]. IEEE Transactions on Automatic Control, 2007, 52(1): 26. DOI:10.1109/TAC.2006.887900


[10] KOSTOUSOVA E K. External polyhedral estimates for reachable sets of linear discrete-time systems with integral bounds on controls[J]. International Journal of Pure and Applied Mathematics, 2009, 50(2): 187.


[11] VARAIYA P. Reach set computation using optimal control[M]. INAN M K, KURSHAN R P. Verification of digital and hybrid systems. Berlin, Heidelberg: Springer, 2000: 323. DOI: 10.1007/978-3-642-59615-5_15


[12] BAIER R, BVSKENS C, CHAHMA I A, et al. Approximation of reachable sets by direct solution methods for optimal control problems[J]. Optimization Methods and Software, 2007, 22(3): 433. DOI:10.1080/10556780600604999


[13] MITCHELL I M, BAYEN A M, TOMLIN C J. A time-dependent Hamilton-Jacobi formulation of reachable sets for continuous dynamic games[J]. IEEE Transactions on Automatic Control, 2005, 50(7): 947. DOI:10.1109/TAC.2005.851439


[14] REISSIG G. Convexity of reachable sets of nonlinear discrete-time systems[C]//Proceedings of the 13th IEEE International Conference on Methods and Models in Automation and Robotics. Szczecin, Poland: IEEE, 2007: 199


[15] GORNOV A Y, FINKEL'SHTEIN E A. Algorithm for piecewise-linear approximation of the reachable set boundary[J]. Automation and Remote Control, 2015, 76(3): 385. DOI:10.1134/s0005117915030030


[16] DONTCHEV A L, FARKHI E M. Error estimates for discretized differential inclusions[J]. Computing, 1989, 41(4): 349. DOI:10.1007/bf02241223


[17] BAIER R, GERDTS M, XAUSA I. Approximation of reachable sets using optimal control algorithms[J]. Numerical Algebra, Control and Optimization, 2013, 3(3): 519. DOI:10.3934/naco.2013.3.519


[18] CHAPRA S C, CANALE R P. Numerical methods for engineers[M]. New York: McGraw-Hill, 2010


[19] WACHTER A, BIEGLER L T. On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming[J]. Mathematical Programming, 2006, 106(1): 25. DOI:10.1007/s10107-004-0559-y



相关话题/优化 文献 系统 计算 理学院

  • 领限时大额优惠券,享本站正版考研考试资料!
    大额优惠券
    优惠券领取后72小时内有效,10万种最新考研考试考证类电子打印资料任你选。涵盖全国500余所院校考研专业课、200多种职业资格考试、1100多种经典教材,产品类型包含电子书、题库、全套资料以及视频,无论您是考研复习、考证刷题,还是考前冲刺等,不同类型的产品可满足您学习上的不同需求。 ...
    本站小编 Free壹佰分学习网 2022-09-19
  • 工业纯钛金属织构标准极图的计算及分析
    工业纯钛金属织构标准极图的计算及分析陈亮维,刘状,虞澜,胡劲,易健宏(昆明理工大学材料科学与工程学院,昆明650093)摘要:工业纯钛中的金属织构会引起各向异性,获得织构信息及分析其演变规律对钛材加工与应用非常重要.本文利用单晶钛的晶体结构数据、乌氏网、极图与织构的定义,建立了纯钛的织构与特定晶面极 ...
    本站小编 哈尔滨工业大学 2020-12-05
  • 316H堆焊UNS N10003合金参数优化、组织和硬度的研究
    316H堆焊UNSN10003合金参数优化、组织和硬度的研究杨飞1,2,黎超文2,李志军2,蒋力2,叶祥熙2,刘芳1(1.上海理工大学材料科学与工程学院,上海200093;2.中国科学院上海应用物理研究所,上海201800)摘要:研究异种合金焊接可以降低熔盐堆结构材料的成本并确保其安全性,本文采用钨 ...
    本站小编 哈尔滨工业大学 2020-12-05
  • 响应面分析法优化不锈钢激光切割工艺参数
    响应面分析法优化不锈钢激光切割工艺参数李永亮1,2,3,王敬1,3,梁强1,3(1.重庆工商大学机械工程学院,重庆400069;2.重庆工商大学工程训练中心,重庆400069;3.制造装备机构设计与控制重庆市重点实验室(重庆工商大学),重庆400069)摘要:为了获得良好的不锈钢激光切割质量,确定合 ...
    本站小编 哈尔滨工业大学 2020-12-05
  • 汽车六角球头冷锻工艺优化与数值仿真
    汽车六角球头冷锻工艺优化与数值仿真陈凌翔,李月超(新乡职业技术学院汽车工程系,河南新乡453000)摘要:冷锻成形工艺是一种少无切削的净近成形工艺,以其精度高、生产效率高、低耗节能等优点,大量使用在汽车零配件的生产。六角球头销是汽车转向系统中的关键零件,其六角成形的质量直接影响到产品的使用性能。本文 ...
    本站小编 哈尔滨工业大学 2020-12-05
  • 抗癌药物作用预测计算方法的研究现状与展望
    抗癌药物作用预测计算方法的研究现状与展望顾兆伟1,张立忠2,刘晓峰3,谭先4(1.长春中医药大学附属第三临床医院脑病康复科,长春130000;2.长春市朝阳区清和社区卫生服务中心,长春130000;3.空军杭州特勤疗养中心康复理疗科,杭州310000;4.东北师范大学信息科学与技术学院,长春1300 ...
    本站小编 哈尔滨工业大学 2020-12-05
  • 液压轮毂马达辅助驱动系统控制策略实车验证
    液压轮毂马达辅助驱动系统控制策略实车验证曾小华,崔臣,张轩铭,宋大凤,李立鑫(汽车仿真与控制国家重点实验室(吉林大学),长春130025)摘要:为充分提升重型牵引车辆通过不良路面的能力,对国内某款重型牵引车在传统结构的基础上加装了前轴液压轮毂马达辅助驱动系统,并针对该混合动力系统,开发了工程化的控制 ...
    本站小编 哈尔滨工业大学 2020-12-05
  • 铁锰复合阴极MFC-EF耦合系统产电及降解RhB效能
    铁锰复合阴极MFC-EF耦合系统产电及降解RhB效能史珂1,2,赵庆良1,2,王维业2,王琨1,2(1.城市水资源与水环境国家重点实验室(哈尔滨工业大学),哈尔滨150090;2.哈尔滨工业大学环境学院,哈尔滨150090)摘要:为实现微生物燃料电池(MFC)微电的原位利用,结合电芬顿(EF)技术的 ...
    本站小编 哈尔滨工业大学 2020-12-05
  • 磷虾群优化的改进粒子滤波算法
    磷虾群优化的改进粒子滤波算法朱震曙,蒋长辉,薄煜明,吴盘龙(南京理工大学自动化学院,南京210094)摘要:标准的粒子滤波存在着权值退化问题,重采样可以解决权值退化问题,但也会带来样本贫化现象.为解决样本贫化问题,提出了一种利用磷虾群优化的改进粒子滤波算法.该算法结合粒子滤波的求解过程,以磷虾个体的 ...
    本站小编 哈尔滨工业大学 2020-12-05
  • 桥梁颤振临界风速的概率密度演化计算
    桥梁颤振临界风速的概率密度演化计算姜保宋,周志勇,唐峰(土木工程防灾国家重点实验室(同济大学),上海200092)摘要:针对桥梁结构自身特性以及外部环境的随机性(如刚度、质量、阻尼比、气动导数等因素)所造成的桥梁的颤振临界风速不确定,难以衡量桥梁颤振稳定性问题.将概率密度演化方法与桥梁颤振多模态耦合 ...
    本站小编 哈尔滨工业大学 2020-12-05
  • 城市快速路小间距路段出口预告标志优化设置
    城市快速路小间距路段出口预告标志优化设置黄利华1,2,3,赵晓华1,2,李洋1,4,荣建1,2(1.北京工业大学城市交通学院,北京100124;2.北京市交通工程重点实验室(北京工业大学),北京100124;3.北京城市系统工程研究中心,北京100035;4.北京警察学院,北京102202)摘要:为 ...
    本站小编 哈尔滨工业大学 2020-12-05