非完整井(精选7篇)
非完整井 篇1
随着时代的进步, 地下水数值模拟也开发出自己的研究软件, 用MATLAB这个数学软件来研究地下水数值模拟也是可行的。
承压含水层完整井非稳定流可归纳为非线性偏微分方程的求解问题, 用MATLAB编写的有限元或有限差分等算法的程序来计算偏微分方程需要很大的工作量。所以对于一些偏微分方程的模型, MATLAB提供了两种方法来解决这些问题。一种方法是用PDE工具箱来解决特殊的偏微分方程求解问题, 但它有较大的局限性, 只能求解部分的椭圆型、抛物型、双曲型、特征值型的二阶偏微分方程问题。另一种方法是采用MATLAB中的函数来解决偏微分方程求解问题, MATLAB已经编写好了这些函数的程序, 在脚本文件中直接调用就可以, 这种方法可以解决不同类型的偏微分方程, 应用范围很广泛。
1 pdepe函数用法
根据承压含水层完整井非稳定流的方程模型, 可采用MAT-LAB中的函数pdepe来求解。pdepe是用来求解一维偏微分方程或方程组定解问题的命令函数, 编写程序时, 先布置结点, 再调用函数。在解决定解问题时需先化成标准形式:
其调用格式为:
m表示标准方程中的指数, xmesh表示标准方程中的, tspan表示时间, 其中xmesh和tmesh应取长度大于等于3的单调递增一维数组。
pdefun是表示偏微分方程的函数, 即生成,
函数。编写格式为:[c, f, s]=pdefun (x, t, u, du) , 其中, x, t是对应标准方程中的与, u对应所求解, du对应的一阶导数, 给定几个输入变量即可表示出c, f, s这三个函数。
icfun是表示偏微分方程初值条件的函数, 定义初值u0 (x) 。编写格式为:u0=icfun (x) 。
bcfun是表示偏微分方程边界条件的函数, 即生成p (x, t, u) 和q (x, t) 函数。编写格式为:[pl, ql, pr, qr]=bcfun (xl, ul, xr, ur, t) , 其中xl、ul和xr、ur分别对应左边界和右边界方程中的与所求解;pl、ql和pr、qr分别对应标准边界条件中左边界和右边界的函数。
对于偏微分方程定解问题来说, pdepe函数的输出参数sol是一个二维数组, 表示所求解u, 即u (i, j) 对应的x (i) 和t (j) 的解为sol (i, j) ;对于元偏微方程组定解问题来说, u= (u1, u2…, un) , 输出参数sol是一个三维数组, sol (:, :, n) 表示方程组的解, 即uk对应x (i) 和t (j) 时的解为sol (i, j, k) 。
2 承压含水层完整井非稳定流模型求解
2.1 定降深模型求解
考虑一均质、等向的承压含水层定流量完整井非稳定流问题, 井里的水位保持定降深。在这个模型中, 用降深作为待求函数, 则模型的定解问题可以写成:
M为厚度, K为渗透系数, Sm为贮水系数, SW为抽水前降深, R为抽水井的影响半径, rW为滤水管半径, r取值范围是:[rw, R]。
为用pdepe函数来进行编程求解, 现将方程化为标准形式为:
根据方程的标准形式可知m=0, 布置空间节点r和时间节点t。在M文件或函数文件中, 直接输入pdepe函数的调用格式即可求出偏微分方程定解S。
程序如下:
编写pdefun函数。根据方程的标准形式, 输出函数:, 在编写格式中, d Sdr表示所求解S的一阶导数。函数文件如下:
编写icfun函数。在本问题中, 初始条件为0, 所以函数文件如下:
编写bcfun函数。根据边界条件标准形式, 输出函数pl=S-SW, ql=0, pr=S, qr=0。函数文件如下:
最后再观察水位在不同时刻的图像, 用H=H0-S来求出水位, H为水头, H0为抽水前含水层测管水头。程序如下:
H0=5;H=H0-sol;
2.2 定流模型求解
在承压含水层完整井中, 抽水流量保持固定的, 而降深是可以变化的。在分析定流量非稳定流问题时, 和分析定降深非稳定流一样, 设均质、等向、等厚、无界的, 不透水底层和原始水头线都是水平线, 无越流补给, 符合裘布依假设。因此, 定流量问题的方程、初始条件和定降深问题一样, 唯边界条件稍有区别, 即定流量问题在井壁处的边界条件是已知抽水定流量。现以定流量抽水, 则方程如下:
用pdepe函数解决这个问题时, 只需修改边界条件即可。
编写bcfun函数。根据边界条件标准形式, 输出函数, ql=, pr=, qr=0。函数文件如下:
定流量不同时刻的水位图像
2.3 阶梯状定流量模型求解
承压含水层阶梯状定流量完整井非稳定流应用于抽水流量呈阶梯状变化的情况。将抽水流量从固定不变推广到呈阶梯状变化即可。设t0到t1时刻, 抽水流量为Q1, t1到t时刻, 为Q2。Q1作用到t时刻引起降深为S1, (Q2-Q1) 从开始作用到t时刻, 引起降深为S2, 因此, t时刻的总降深为S1+S2。设在时刻前, 井的抽水流量Q1, Q2, …, Qm共有n段, 这一变化过程可以看作是n个流动作用过程的叠加。
用pdepe函数编程求解时, 只需将bcfun函数文件中Q这个固定数值变成分段函数Q1, Q2, …, Qn即可。
当Qi>Qi-1时产生的降深水位图像
当Qi>Qi-1时产生的恢复水位图像
3结论
采用MATLAB中求解偏微分方程的函数来解决承压含水层完整井非稳定流问题是有效的, 在计算过程中, 绝对误差不大于0.2, 并且它不需要复杂的数据运算过程, 只需求解问题满足所给函数的标准公式, 输入参数即可。这个方法使用起来方便简单, 也适用于其他地下水模型问题, 在此基础上也可以解决复杂的地下水问题, 对于解决地下水数值模拟问题有很大的应用前景。
李瑞 (1988-) , 女, 辽宁大连人, 辽宁师范大学数学学院在读硕士研究生, 研究方向:偏微分数值解法。
《资源节约与环保》2013年第11期
参考文献
[1]王君连.工程地下水计算[M].北京:中国水利水电出版社, 2004.
[2]苏金明, 阮沈勇, 王永利.MATLAB工程数学[M].北京:电子工业出版社, 2005, 8.
[3]陈崇希, 林敏.地下水动力学[M].武汉:中国地质大学出版社, 1999, 10.
非完整井 篇2
一、下井带班领导必须和职工同上同下,要把保证矿井安全生产作为第一责任,切实掌握当班井下的安全生产状况,加强对重点部位、关键环节的检查巡视,及时发现和组织消除隐患,及时制止违章违纪行为,严禁违章指挥、违章作业、当发现有危及职工生命安全的重大隐患和严重问题时,带班领导必须立即组织采取停产、撤人、排除隐患等紧急处置措施,坚决做到不安全不生产。
二、建立健全“矿井负责人下井带班交接班制度”、矿井负责人下井带班档案管理制度、矿井负责人下井带班考评制度、矿井负责人当月下井带班计划表和矿井负责人下井带班公示栏。
三、矿井负责人当月下井带班计划表由矿山生产部门根据工作情况每月进行调整,并在月初上墙公示。
四、下井带班作业种类包括采掘工作面、辅助运输线路、主运输线路、主提升、通风排水及其它零星工程。带班人员根据各自岗位和分工职责以及重点进行带班。当井下遇到巷道贯通、初采、末采、采掘工作面通过老空区、小窑巷道、处理火区、探放水、边远地点、过地质构造区、检修工程、反风、烧焊、密闭采空区、放炮、突发事故等重点工程和特殊时期时,带班下井人员要到重点现场检查、把关、指导。
五、带班下井人员范围:矿级领导、副总工程师和矿长助理;生产技术室、安监科、调度室、地测科、机电科等科室的副科级及以上领导干部。
六、每月带班下井的次数规定:正矿级领导带班次数不少于5个/月,副矿级领导带班次数不少于6个/月,副总工程师、矿长助理带班次数不少于7个/月,各单位行政正职带班次数不少于3个/月,其他领导干部带班次数由调度室根据实际情况进行合理安排。
七、根据公司实际出发,以每班都安排有矿领导带班下井为原则,与工人同时下井、同时升井,保证矿井井下24小时有矿领导带班。
非完整井 篇3
深水井作业受地层压力难以预测、低温环境、浅层水气流动、储层流体不易控制等复杂因素的影响[1],井完整性极难保证,且事故危害较大,严重时可导致井眼报废、船毁人亡等恶性后果。2010年4月20日,“深水地平线”钻井平台爆炸事故造成11人死亡,经济损失达140亿美元的,大面积海域遭到污染。调查表明该事故发生的直接原因是在替海水过程中井完整性失效,油气从储层大量泄漏所致[2]。保证油气井完整性的核心是强化井屏障系统对储层流体的有效隔挡,防止储层流体失控流动。目前相应的研究主要有油气井完整性管理[3]、井完整性失效分类[4]以及多属性决策进行最佳井屏障选择[5],但缺乏对多因素协同作用下井完整性失效的深入研究。本文运用SBFT方法,逻辑演绎了深水井完整性失效的安全屏障及影响因素之间的关系,建立了安全屏障模型和事故树模型,通过对最小割集、最小径集和重要度的计算,讨论预防井完整性失效的主要途径,为有效防止因地层流体无控制流动对人员、环境和财产造成损失,保证深水井作业高效进行提供了理论依据。
1 SBFT模型机理
本文所提出的SBFT分析方法是基于安全屏障模型(SBE)和事故树分析(FTA)发展而来的,是一种定性和定量相结合的方法,可用于深水井作业中井的完整性的失效分析。SBFT分析流程如图1。
SBE模型最早来源于“瑞士奶酪模型(Swiss cheese model)”,该模型是英国曼彻斯特大学教授James Reason[6]于1990年在“Human Error”提出,并获得核能、宇航、海运、医学、石化、采矿等行业广泛认可和应用。Swiss cheese model由4片具有许多孔洞的奶酪组成,每一片奶酪代表一种防御体系;每一片奶酪上的孔洞代表该防御体系的漏洞或缺陷,主要包括组织因素、不安全的监督、不安全行为的前兆和不安全的操作行为等。本文运用SBE模型分析导致井完整性失效的安全屏障以及其逻辑关系,其功能原理如图2所示。每层安全屏障代表保持井完整性的井屏障,由于深水环境因素、操作因素、设计因素、组织管理等因素的影响,每层井屏障中存在着许多漏洞或缺陷(失效因子),这些失效因子随着时间在不断变化。不同层面的井屏障对于缺陷或漏洞可以互相拦截,井完整性就不会因为单一的失效因子而失效。当某一瞬间每层井屏障上的这些缺陷排列在一条直线上,井的完整性就失效了。
FTA方法是由美国贝尔实验室于20世纪60年代基于图论提出的一种安全评价方法[7]。基本原理是:以系统失效为目标,运用逻辑树进行演绎分析失效过程,寻求系统失效的原因。由于事故树刻画事故的因果关系直观、清晰及逻辑性强,因此被用于复杂系统的风险分析[8]。事故树定性和定量分析的技术指标有最小割集、最小径集、重要度和顶事件失效概率分析[9,10]。
2 深水井完整性失效危险源分析及井屏障
2.1 井的完整性失效危险源分析
2.1.1 井完整性定义
挪威NORSOK标准D-010[11]《钻井和作业的油气井完整性》定义油气井完整性为“在井眼整个寿命期间,应用技术、操作和组织措施最大限度地减少地层流体不可控制地流动的风险,以确保油气井始终处于安全可靠的状态。” 所谓深水井完整性是指深水井处于地层流体被有效控制的安全运行状态。一旦地层流体发生不可控制地层间流动或流向地表,深水井在功能上不具有完整性,即深水井完整性失效。
2.1.2 井完整性失效危险源辨识
许多研究机构均曾对井完整性所面临的挑战开展过研究工作[12,13],研究结果表明导致井完整性失效主要危险源是井屏障失效,井屏障主要失效因子有以下几方面:①设备设计缺陷:油套管接头螺纹、管柱、胶筒和井口阀体设计缺陷等。②操作因素:操作压力不足、上扣扭矩不足和下放速度过快等。③环境因素:海水侵蚀、酸性气体腐蚀、高温高压和高压流体冲刷井口阀体等。④设备部件失效:密封元件失效、封隔器胶筒损坏和井下安全阀控制管线故障等。⑤其他因素。
2.2 井屏障(WBEs)
有效的WBEs系统是保证深水井完整性的关键。国际标准化组织(ISO)和美国石油学会(API)标准提出确保WBEs系统的可靠性、可用性和可维护性来保障油气井的长期完整性,确保钻井及相关作业的顺利进行。挪威大陆架(NCS)关于石油作业的HSE设计和装备等相关法规设施条例第47条规定有关井的屏障条款,在井的寿命周期内确保WBEs有效,WBEs设计应确保该井的完整性[14]。保证油气井完整性的核心在于利用有效的WBEs系统控制地层流体不可控制的流动。
NORSOK D-010定义WBEs系统为一个或几个相互依附的WBEs组件的集合。该标准定义了从根本上防止地层流体流出的初级WBEs和二级WBEs,不同作业过程WBEs规定有所不同。本文对起出BOP,下入水下采油树这一过程的WBEs元件、WBEs功能、监测参数及失效形式进行分析如表1。
3 SBFT建模及失效分析
3.1 SBFT模型建立
3.1.1 SBE模型建立
深水井完整性失效分析面临着一系列的挑战。主要挑战是储层流体不可控制地流向地表存在多种泄漏方式和泄漏路径。每一种泄漏路径均会导致井完整性失效。为了便于分析失效因子、WBEs之间的逻辑因果关系,将流体从储层泄漏到地表归纳为3种泄漏方式和6种泄漏途径:
(1)管柱(CS)泄漏。包括3种泄漏途径: CS和TH失效,TP、SCSSV和TH失效,TP、SCSSV、CS和TH失效。
(2)井眼(WB)泄漏。包括2种泄漏途径:CP和TH失效,CC、CA和TH失效。
(3)套管与套管环空(AN)泄漏。考虑1种泄漏途径:CC和WH失效。
储层流体泄漏途径如图3所示。在深水井作业中,如果井眼中两个地层带之间的压差使流体产生无控制流动时,至少应有一个可利用的有效WBEs进行控制;如果该压差可能会导致流体从储层无控制地流到外部环境中时,至少应该有两个可利用的WBEs进行控制。初级WBEs失效,二级WBEs如果能及时控制,井完整性就不会失效;但初级和二级WBEs失效之后,井完整性则会失效。本文以起出BOP下入采油树这一过程为例,建立SBE模型(图4)。
3.1.2 事故树建模
根据SBE模型,本文针对起出BOP下入采油树这一过程建立了以井完整性失效为顶事件,以39个WBEs失效因子(表3)为底事件的事故树模型如图 5-图6所示。
3.2 井完整性失效分析
3.2.1 最小割集和最小径集计算及分析
最小割集代表井完整性失效的基本事件组合,每一个最小割集即一种失效途径。运用布尔代数法对井完整性失效事故树简化求解,得到该事故树的80个最小割集:{ (X5X11 ),(X5X12),…}12个,{ (X12 X16X20),(X12X16X19),…}46个,{ (X9X10X26X27),(X9X10X30X34),…}20个,{ (X9X10X28X29X34),(X9X10X28X29X35),…}2个。说明该过程存在80种可能的井完整性失效模式,也反应了影响井完整性失效的原因复杂。
最小径集表征了井完整性失效的最佳预防措施,借助“成功树”的最小割集间接得到事故树的最小径集:{ (X10X11X12X37X39 ), (X9X11 X12X13X14X15X16X29X30X31),…}共 33个。其中(X10X11X12X37X39 )是防止井完整性失效的最佳方案。只要保证卡瓦功能正常、TH密封件功能有效、TH没有被腐蚀失效、密封垫环紧密、井口阀体功能有效,即保证井屏障TH、WH防御功能正常,这样不同井屏障对于缺陷或漏洞可以互相拦截,井完整性就不会失效。其他方案的基本事件均较多,能做到逐一控制不是一件容易的事情,所以在某种意义上不理想。
3.2.2 顶事件概率计算及分析
计算顶事件概率首先要计算各基本事件的概率,基本事件发生概率是基于OREDA(Offshore Reliability Data Handbook, 2002)[15]和世界海洋事故数据库(WOAD)[16]等失效数据统计,利用基本事件发生的频率近似代替其概率[17]所得。基本事件的概率如图7所示。根据油气井完整性论坛对井完整性分类[4],对由于初级屏障和二级屏障共同失效导致井完整性失效的井数进行统计分析,其所占比值为0.01。
顶上事件概率Q近似计算式如下:
式中:j、s—最小割集序数; i—基本事件序号;k—最小割集数;Xi∈Ej—属于第j个最小割集的第i个基本事件。
可得顶上事件发生的概率Q=0.0103,经比较,顶上事件概率的计算结果与上述统计值(0.01)的相对误差为2.912%,在工程允许的范围之内,可知该事故树的结构合理,同时也表明运用SBFT方法分析深水井完整性失效是可行的。若将基本事件概率提高1倍,则顶事件概率提高至原来的5倍。若将基本事件概率降低1倍,则顶事件概率降低为原来的1/5。因此要严格控制基本事件的发生概率,进而降低井完整性失效概率。
3.2.3 重要度计算与主要失效因子分析
(1)结构重要度I(i)利用最小割集求解结构重要度系数,近似计算式:
式中:ni—基本事件Xi所属最小割集包含的基本事件数。
(2)概率重要度Ig(i)是顶上事件发生概率函数Q对自变量qi求偏导,计算式为:
(3)临界重要度I
通过式(2)-(4)计算各基本事件3个重要度,分析其对深水井完整性失效的影响程度。基本事件的各个重要度值及顺序比较如图8。依据图8可找出39个基本事件中,3个重要度值均是最大的基本事件是TH受到CO2腐蚀(X12)及密封部件损坏(X11)。从结构重要度分析, X12、X11对顶事件的发生影响最大;从概率重要度分析知,降低X12、X11的发生概率,能迅速有效地降低顶事件的发生概率;从临界重要度分析知, X12、X11不仅敏感性强,而且本身发生概率较大,所以其重要程度仍然最高。可知 X12、X11是井完整性失效的主要失效因子,在制定失效预防措施时应首先考虑该类基本事件。
其次概率重要度和临界重要度均较大,结构重要度较低的基本事件为封隔器密封部件损坏(X5)、油管柱粘扣(X25)、管柱蠕动(X8)、胶皮损坏(X7)、上扣扭矩不足(X24)、油管接头设计缺陷(X23)、层间压差大(X6)。降低该类事件的发生概率对顶事件概率影响显著,要避免其发生即要确保CP和CS功能正常。再者基本事件X33、X15、X13……X14结构重要度较其他二者比较大,要尽量从结构上避免。其他基本事件重要度值均较小,改变该类基本事件对顶事件影响不显著,一般不予考虑。
4 结论
(1)本文以起出BOP下入水下采油树作业过程为研究对象,建立了井完整性失效的定性和定量分析模型。经计算得到井完整性失效的概率与统计值相符,证明运用SBFT风险分析方法分析深水井完整性失效是可行的。
(2)确立了深水井完整性失效的3种泄漏方式和6条泄漏途径,并对WBEs之间的逻辑关系进行了梳理,以此构建了防止储层流体泄漏到外部环境的SBE模型。
(3)在SBE模型的基础上,建立了井完整性失效事故树模型,对最小割集和最小径集分析表明井完整性失效有80种可能的途径和33种预防措施组合;结构重要度、概率重要度和临界重要度的计算结果表明该过程的最薄弱环节是油管悬挂器,主要失效因子为CO2冲刷腐蚀以及密封部件损坏,从而使得预防措施具有较强的理论依据。
非完整井 篇4
该文的基本思路是在运动学方程的基础上, 主要考虑动力学方面的影响, 建立动力学模型, 对机器人的轨迹跟踪问题进行研究, 并运用模型预测控制的方法寻找最优控制函数, 得出最优输入。最后, 利用MATLAB软件进行仿真, 验证预测控制方法的可行性。
1 轮式移动机器人建模
该文主要采用的是具有一个前轮驱动和两个后轮的非完整移动机器人作为被控对象。
假设移动机器人的初始位姿为I=[x, y, φ]。其运动学方程为:
因为非完整移动机器人受到纯滚动无滑动这个约束条件, 所以上述运动学方程可以表示为:
接下来将建立在运动学方程基础上的动力学模型。一般都将非完整移动机器人的动力学模型表达成一下的方程式:
对运动学方程 (2.2) 两边求导可得:
整理得:
最后, 结合运动学模型, 就可以得出非完整运动机器人的动力学模型, 如下:
2 终端状态控制器以及控制函数的设计
模型预测控制基本理论就是根据历史信息与未来的输入预测出未来输入的一种计算机算法, 它最主要的就是求解出系统的最优输入, 获得最优跟踪轨迹。而最优输入是通过最优目标函数获得的, 建立如下二次目标函数:
但是, 在 (3.1) 中发现的控制器是不能保证稳定的, 所以通过添加终端状态惩罚到值函数, 同时添加一个终端状态约束到移动机器人控制器来保证稳定。于是, 跟踪问题的值函数变化如下:
在时间t, 在滚动时域控制框架中, 在线解决的开环优化问题可用如下公式表达:
如下是控制算法内容的描述:
(4) 上述过程将继续下去, 直到控制实现了令人满意的性能。
3 仿真
该文是选择Matlab软件在个人计算机上进行仿真测试的, 主要跟踪了正弦轨迹, 优化步骤的最大数字选为80, 优化步骤中最大数字越大, 则跟踪性能越好, 然而也就需要越多的计算时间。
跟踪正弦曲线
设定正弦的参考轨迹为:
进行正弦曲线测试, 仿真轨迹如图 (a) 、 (c) 所表。
由仿真结果可以看出, 虽然在开始跟踪的时候速度和状态有误差, 但是最终机器人的轨迹将收敛到参考曲线, 控制信号收敛到参考控制信号, 以及误差状态收敛于零。
4 结语
该文提出了一种稳定的终端状态控制器, 用来跟踪控制非完整移动机器人。其能够很好地解决局部稳定性的问题, 实现机器人对于轨迹的跟踪。但是该方法也存在一定的缺陷, 即计算量是使用该控制器实时系统中的一个重要问题。如何提高计算效率仍在进一步研究中。提出的终端状态控制器需要一个初始可行的解决方案。目前, 正在使用试验和误差的方法。初始解决方案的可行性分析是人们未来的工作。
参考文献
[1]韩光信, 陈虹, 马苗苗, 等.约束非完整移动机器人轨迹跟踪的非线性预测控制[J].吉林大学学报 (工学版) , 2009 (1) :177-181.
[2]曾志文, 卢惠民, 张辉, 等.基于模型预测控制的移动机器人轨迹跟踪[J].控制工程, 2011 (S1) :80-85.
[3]谷东兵, 胡豁生, Michael Brady.移动机器人的运动预测控制[J].仪器仪表学报, 2000 (2) :155-158.
[4]马海涛.非完整轮式移动机器人的运动控制[D]合肥:中国科学技术大学, 2009.
非完整井 篇5
近年来, 自动引导车 (automatic guided vehicle, AGV) 作为非完整约束轮式移动机器人 (wheeled mobile robot, WMR) 的一种, 由于实际应用的广泛性, 其控制问题也得到了国内外学者的关注。机器人的运动控制大致可以分为镇定控制与轨迹跟踪控制两类。作为一种具有非完整约束的典型系统, WMR不满足反馈稳定的Brockett定理的必要条件, 因此不能通过连续时不变状态反馈控制方法实现机器人系统的镇定, 故从理论而言, 设计镇定控制器更具有挑战性。但从实践意义上来说, 轨迹跟踪更具有实际应用价值, 因为即使对于镇定控制任务, 机器人也要沿某条可行轨迹运动到目标点。为了解决WMR的运动控制问题, 人们已经提出各种控制方法, 包括滑模控制方法[1]、反演控制方法[2]等在内的非线性控制方法, 自适应控制方法[3], 神经网络控制方法[4], 时变反馈控制法[5], 非连续反馈控制法[6]等, 取得了较好的控制效果。然而, 上述传统方法的可调参数与控制性能之间的关系不是很直观, 参数整定很费力, 且在考虑实际存在控制输入约束或状态约束时, 比较难处理。而模型预测控制 (model predictive control, MPC) 可以克服以上缺陷, 模型预测控制利用求惩罚函数的最小化来确定最优控制输入, 参数调整容易, 且状态约束及控制输入约束可以在优化过程中直接考虑在内, 简单明了。模型预测控制亦称为滚动优化控制 (receding horizon control, RHC) , 是近年来发展起来的一类新型计算机控制算法, 由于它采用模型预测、滚动优化和反馈校正等控制策略, 且其处理约束能力强, 因而控制效果好, 已成功应用于石油、化工、机械等比较复杂且模型不容易精确建立的工业生产过程。但是, 由于源于工程控制的预测控制在设计时没有在理论上考虑系统的闭环镇定性, 预测控制是在线反复求解开环优化问题, 而开环最优不能保证闭环镇定, 因此, 如何保证预测控制的镇定性成了一个关键问题。
本文在研究AGV系统的非线性模型预测控制问题时, 在基本的预测控制原理基础上着重考虑了预测控制的镇定性, 设计了终端控制器和终端域, 提出一种跟踪与镇定统一的控制策略。同时, 考虑到运动过程中总会存在噪声干扰, 使得AGV自定位不够准确, 影响轨迹跟踪, 本文用一个简单的状态观测器对AGV状态进行估计, 从而进一步提高控制精度。最后本文还增加了一个避障功能模块, 以此来进一步完善整个控制器功能。
1 AGV模型和预测控制模型
以具有非完整约束的三轮移动机器人小车作为控制对象, 其前轮既是转向轮又是驱动轮, 结构如图1所示。
假设小车当前实际位姿为X= (x, y, θ) T, 任意给定的参考位姿为Xr= (xr, yr, θr) T, 位姿误差Xe= (xe, ye, θe) T;当前实际速度U= (v, ω) T, 参考速度Ur= (vr, ωr) T, 其中, vr为参考线速度, ωr为参考角速度。图1中, D为后轮轴与中心轴之间的交点, M为小车质心。非完整约束意味着小车的轮子只有纯滚动而无滑动, 即
其中, 系统状态矢量和系统控制输入矢量分别表示为X= (x, y, θ) T和U= (v, ω) T。
移动机器人小车在坐标内描述的位姿误差Xe为
利用式 (1) 、式 (2) , 可以得到如下位姿误差微分方程:
为了方便下文的表达, 笔者重新定义控制输入为
式 (3) 可改写成如下形式:
2 轨迹跟踪的模型预测控制算法设计
模型预测控制是一种基于模型的控制, 其机理可以描述为:在每个采样时刻, 预测模型根据获得的当前系统实际状态信息, 在线求解一个有限时域开环最优控制问题, 将所得到的控制序列的第一项应用于被控系统, 剩余的控制输入丢弃, 直至下一采样周期。在下一采样时刻测得系统的实际状态, 重复求解优化过程, 即所谓的滚动优化。
然而, 上述开环的过程并不能保证被控系统的镇定。为了解决这个问题, 可以在求解最优问题时在惩罚函数中加入一个终端惩罚函数项以及终端约束, 并引入一个终端域的概念。本文的控制思想是把整个控制过程分成两个控制阶段:首先, 在有限时域内求解最优控制问题, 将获得的最优控制输入到系统, 驱动系统使其进入到终端域;然后, 利用一个终端线性反馈控制器将系统控制到平衡点。预测控制过程如图2所示, 图中, t为当前时刻, Tp为预测时域。
此外, 由于观测误差直接影响控制效果, 为了能更加准确地估计AGV的当前位姿, 本文为其设计了状态观测器
因此, 本文的非线性模型预测跟踪控制思路是:系统给定AGV的参考速度Ur和参考位姿Xr, 它们可以由路径规划算法求得;期望位姿Xr与由状态观测器观测得到的AGV当前位姿
2.1 MPC优化控制器设计
考虑由以下微分方程来描述的非线性系统:
式中, Xe (t) ∈Rn为系统状态;Ue (t) ∈Rm为系统控制输入。
系统状态及控制输入约束如下:
Xe (t) ∈χ, Ue (t) ∈ζ ∀t≥0
式中, χ、ζ为包含原点在内的密集。
选取惩罚函数为预测状态及相应控制输入的二次函数:
式中,
为了保证闭环镇定, 根据文献[7,8], 可以引入一个终端惩罚项V (·) 。改进后的惩罚函数和优化问题为:
式中,
本文选如下一个Lyapunov函数作为终端惩罚项:
式中, P为对应的加权矩阵。
2.2 MPC终端控制器设计
根据文献[9]的定理, 如果存在一个终端线性反馈控制器ULe (t) = (u1e (t) , u2e (t) ) T, 使满足值函数
则系统能渐进镇定。
本文定义如下一个终端线性状态反馈控制器:
式中, α、β为正常数。
为了表示区别, 我们用XeT= (xeT, yeT, θeT) T代表终端状态下的位姿误差, 其中, xeT为终端状态下的x方向位置误差; yeT为终端状态下的y方向位置误差;θeT为终端状态下的θ姿态角误差。
假设函数
把式 (8) 代入式 (9) , 得
为了满足镇定条件
同时终端域Ω应被定义为
此外, 终端域Ω还受到控制输入约束的限制。由于受到来自电机输出扭矩的上下限约束的限制, 所以设本文控制输入约束为
结合式 (4) 及式 (8) 可得到如下终端域Ω的限制:
综上, 结合控制框图3, 控制算法步骤可以概括如下:
(1) 在某一时刻t, 基于通过反馈得到的系统状态信息
(2) 设δ为采样时间, 在下一采样周期t+δ时刻, 根据反馈得到的此刻的系统实际状态信息
(3) 重复上面的步骤直至达到满意的控制效果。
2.3 状态观测器设计
上述模型预测算法需要在每个采样时刻由状态观测器观测系统状态
设已经求得最优问题的解为U*e, 根据Ue与U的关系式 (4) 容易得到最优控制U*= (v*, ω*) T。
根据移动机器人小车的运动学方程式 (1) , 设
对于系统
令φk=g′ (Uk, k) 即
Xk+1≜AXk+Bφk (15)
系统输出为
Yk=CXk (16)
式中A 、B和C为具有相同维数的常系数矩阵。
设计如下观测器:
式中, L为观测器增益矩阵。
定义观测器误差为
则系统为
式中, φek为计算误差, φek=g′ (U*k, k) -φ*k;N为预测步数, N=Tp/δ。
2.4 避障模块设计
假设AGV小车上部环绕安装有若干个红外传感器, 它们能随时随地检测到前方的障碍物信息并反馈给小车控制器, 以便能及时调整轨迹跟踪控制方案使之避过障碍物而不发生碰撞。对于移动障碍物, 可以利用预测控制的方法来预测下一个Tp的障碍物位姿信息以及障碍物移动速度, 并通过计算类似上面的优化问题来获得系统的最优控制输入。设小车共环绕安装有Ns个红外传感器, 第j个的传感器测得的障碍物信息为[θob, jdob, j] (j=1, 2, …, N) , 其中, θob, j为第j个传感器与前方障碍物的夹角;dob, j为第j个传感器距离前方障碍物的直线距离。
避障惩罚项为
式中, kob 、c1、c2为正常数。
改进之后的惩罚函数J为如下形式:
3 仿真研究
为了验证本文所提出的控制算法的正确性和有效性, 在MATLAB环境下编程, 对AGV进行仿真计算和性能分析比较。
AGV期望轨迹取以 (0, 0) 为圆心, 以1为半径的圆, 参考速度为vr=0.2m/s, ωr=0.2m/s。AGV小车实际初始值取x (0) =1.2m, yr (0) =1.2m, yr (0) =-0.3m, θ (0) =2π/3, v (0) =0.4m/s, ω (0) =0.3rad/s。
(1) 在仅使用基本的模型预测跟踪控制率而不使用所设计的终端控制器的情况下, 对AGV小车进行跟踪控制仿真, 取基本模型预测控制器性能加权矩阵Q=diag (0.5, 0.5, 0.5) , R=diag (0.1, 0.1) , 预测时域与控制时域取为Tp=Tc=1s, 采样周期取为δ=0.1s, 即向前预测步数为N=10。跟踪效果如图4所示。
(2) 结合所设计的预测模型终端控制器和状态观测器对AGV进行仿真试验, 取模型预测控制器性能加权矩阵P=diag (1, 1, 1) ;性能系数取为α=2, β=1。预测时域Tp、控制时域Tc和采样周期δ参数设置同上。图5为AGV小车的跟踪效果图。
此外, 为了验证所设计避障模块的有效性, 对期望轨迹为简单的直线运动的AGV进行避障仿真试验。设期望轨迹上有四个障碍, 避障控制器的系数取为kob=1.5, c1=0.1, c2=10, 图6是小车躲避静态障碍物的效果图。
比较图4和图5仿真结果可见, 由于本文所设计的控制律在设计时就考虑控制约束和终端镇定, 而且在相应控制输入作用下, 跟踪轨迹具有最优意义, 同时状态观测器也进一步提高了反馈精度, 因此在满足输入约束的前提下, 系统能够在较短时间段内跟踪指定轨迹, 取得良好的跟踪性能, 表现为跟踪路径平滑、跟踪快、速度平稳。从图4可以看出, 如果在设计控制器时不考虑实际存在的控制约束和终端镇定, 闭环系统就很有可能会失去了稳定性保证, 也就很难获得令人满意的跟踪效果。此外, 图6避障仿真结果同时也验证了所设计的避障功能模块的有效性。
4 结论
以三轮AGV为研究对象, 研究存在非完整约束和控制输入约束的AGV轨迹跟踪的非线性模型预测控制问题, 设计了一种终端镇定控制器, 并结合AGV的运动模型, 使用状态观测器对含噪声的动态系统进行状态估计, 进一步增强轨迹跟踪控制的效果, 实现了系统全局渐进稳定。此外还考虑AGV的避障功能, 设计了一种避障控制算法。计算机仿真结果证实了所设计轨迹跟踪算法的有效性和实时性。
摘要:针对存在非完整约束和控制输入约束的三轮AGV模型, 研究其非线性模型预测控制策略, 提出了一种跟踪与镇定统一控制算法。首先由模型预测控制原理产生一个优化控制器, 接着设计终端控制器及不变终端域来保证系统的镇定, 然后设计状态观测器以进一步提高跟踪精度, 最后设计避障控制模块来完善整个控制器功能。该控制算法具有一般性, 对所有AGV及其他轮式移动机器人的运动学模型都适用。计算机仿真结果验证了该控制算法的正确性和有效性。
关键词:AGV,模型预测控制,镇定,观测器,轨迹跟踪
参考文献
[1]闫茂德, 吴青云, 贺昱曜.非完整移动机器人的自适应滑模轨迹跟踪控制[J].系统仿真学报, 2007, 19 (3) :579-581.
[2]Dierks T, Jagannathan S.Control of NonholonomicMobile Robot Formations:Backstepping Kinematicsinto Dynamics[C]//16th IEEE International Con-ference on Control Applications, Part of IEEE Multi-conference on Systems and Control.Singapore, 2007:94-99.
[3]Farzad P, Mattias P K.Adaptive Control of DynamicMobile Robots with Nonholonomic Constraints[J].Computers and Electrical Engineering, 2002, 28:241-253.
[4]Lin Sheng, Goldenberg A A.Neural-network Con-trol of Mobile Manipulators[J].IEEE Transactionson Neural Networks Manipulators, 2001, 12 (5) :1121-1133.
[5]Samson C.Control of Chained Systems:Applicationto Path Following and Time-varying Point-stabi-lization of Mobile Robots[J].IEEE Transactions onAutomatic Control, 1995, 40:64-77.
[6]Fontes F A C C, Magni L.Min-max Model Predic-tive Control of Nonlinear Systems Using Discontinu-ous Feedbacks[J].IEEE Transactions on ControlSystem Technology, 2003, 48 (10) :1750-1755.
[7]Chen H, Allgower F A.Quasi-infinite HorizonNonlinear Modelpredictive Control Scheme withGuaranteed Stability[J].Automatica, 1998, 14 (10) :1205-1217.
[8]Mayne D, Rawling J, Rao C, et al.Const RainedModel Predictive Control:Stability and Optimality[J].Automatica, 2000, 36:789-814.
非完整井 篇6
1 套管-水泥环-地层系统力学模型分解
对处于非均匀应力状态下的水泥环受力分析可由图1表征,图组中由内到外三个圆环分别表示套管,水泥环和地层,套管内外壁半径分别是R1和R2(mm);水泥环外壁内径是R3(mm);为了方便引入无限远地层半径R4(R4>>R1)(mm)。套管受内压为Pw(MPa);地层受非均匀水平地应力σH和σh(MPa)。图1中(a)最终可分解为(b)+(c),(b)表示套管-水泥环-地层系统受内压和受平均远场应力的外压,(c)表示水泥环系统受偏差远场应力的外压,其中(b)是个轴对称问题,(c)是个非对称问题。
如果引入记号σ和s:
远场(R=R4)的正应力和切应力分布表达式可表达成[6]。
式中θ为计算方位与最小水平地应力方位的夹角,求解套管-水泥环-地层应力分布,可分解为平均应力分量和偏差应力分量两部分,平均应力分量由套管壁受均匀内压Pw和受远场均匀外压σ两部分组成,为简单起见,平均应力分量部分可合写成式(1)形式:
偏差分量写成式(2)形式
式(1)对应于图1(b)中受均匀外压σ和内压Pw,可由拉梅经典解法[7]解出来,拉梅问题在文中称为问题A,式(2)对应于图1(c),受偏差分量外压,称为问题B。这样原问题可分解为问题A和问题B的叠加。
2 套管-水泥环-地层系统应力分布求解
2.1 问题A的解
针对问题A,套管-水泥环-地层系统的应力场和位移场可简化为以下形式
式(4)中
式中,i=1,2,3分别代表套管,水泥环和地层;σri、σθi和Uri是各部分对应的径向应力,周向应力和径向位移;Gi和υi是各部分对应的剪切模量与泊松比;r是距井眼中心的距离。
2.2 问题B的解
考虑到远场偏差应力边界如式(2)所示,可设应力函数[8]
针对问题B,得出在以上给定应力函数条件下的水泥环系统对应的应力和位移通解形式
2.3 水泥环系统应力和位移分布总形式
求出水泥环系统的径向应力,周向应力和剪应力之后便可以得到其轴向应力分布,相比于地层,由于水泥环和套管处于失重状态,所以其轴向应力分布表达式不同,水泥环系统应力和位移分布总形式如式(7)和式(8)。
式(7)中,α是BIOT系数,水泥环如果当作考虑渗透性的孔隙介质,应该考虑水泥环的有效应力,所以BIOT系数不应该取为0;Pp是地层孔隙压力;σV是地层上覆岩层压力。
方程中总共有18个未知系数,对应18个线性方程,可通过编写程序求解。其中,问题A对应6个待定系数,共有2个应力边界条件,2个应力接触条件和2个位移连续条件,问题B对应12个待定系数,共有4个应力边界条件,4个应力接触条件和4个位移连续条件。具体可参考文献[9]。
3 水泥环水力胶结失效的判别准则
对水泥环来说,受力后可发生剪切,拉伸和滑移破坏,发生剪切破坏的判定准则为
式(9)中,σ1,σ3分别为水泥环最大和最小主应力,C、φ则为水泥环的黏聚力和内摩擦角。
水泥环发生拉伸破坏的强度准则为
式(10)中,σt为水泥环抗拉强度。
水泥环发生拉伸破坏的强度准则为
式(11)中,j=1,2分别代表套管-水泥环和水泥环-地层界面;τj',σj分别为对应界面的剪应力和正应力;τj为对应界面抗剪强度;fj为对应界面摩擦系数。
4 算例分析
以东海宁波区块的某井(C1井)的储层段为研究对象,该区块储层段井底深度为4 725 m,由于水泥环的渗透性很差,本研究把水泥环BIOT系数取为0,孔隙压力为1.545 g/cm3,水平最大地应力为92.27 MPa,水平最小地应力为81.82 MPa,地层岩石弹性模量为16.4 GPa,泊松比为0.182,黏聚力13.8 MPa,内摩擦角37.2°;套管弹性模量210 GPa,泊松比为0.21;水泥环弹性模量6.4 GPa,泊松比0.169,黏聚力10.8 MPa,内摩擦角26.1°;井眼直径为215.9 mm,套管直径为177.8 mm,套管壁厚为9.19 mm。
图2和图3给出了水平最大和最小地应力方位水泥环-地层系统应力分布水平规律,可以看出,径向应力是连续的,而周向应力和轴向应力呈阶跃特征,剪应力为0,在地层远处,径向应力和周向应力趋近于原始地应力。
图4和图5分别给出了水泥环内壁和外壁应力分布随方位角的变化规律,可以看出,在水平最小地应力方位,水泥环内壁主差应力最大,水泥环最容易进入塑性,这可根据式(9)判断;另外,在水平最大地应力方位,水泥环内壁周向应力和轴向应力相比于径向应力要小得多,但是依然为压应力,所以水泥环产生拉伸破坏的可能比较小,这与国外一些文献中的研究结论不同[1,2],因为这些研究结论中并未真实地模拟地层地应力对水泥环造成的挤压作用,这样高估了套管内压对水泥环膨胀造成的拉应力,水泥环是否产生拉伸破坏可根据式(10)判断;水泥环内壁处的剪应力在45°方位最大,水泥环内壁有产生剪切滑移的可能性,水泥环是否产生滑移可根据式(11)判断,无论哪种破坏形式可看出水泥环内壁是危险面。
由于水泥环在最小水平地应力方位下,内壁处最容易进入塑性,本研究关于水泥环是否进入塑性的结论都从最小地应力方位研究。
图6、图7和图8分别给出了水平最小地应力方位分别改变试压压力,弹性模量和泊松比时水泥环区域应力分布规律,可以看出,当增大试压压力时,水泥环区域径向应力增大,周向应力减小,而轴向应力几乎不变,由于主差应力变大,所以水泥环进入塑性的风险增大;当增大弹性模量时,水泥环区域三种应力都增大,但总体表现为主差应力变小,所以水泥环进入塑性的风险增大;当增大泊松比时;水泥环区域径向应力几乎不变,而周向和轴向应力增大,总体表现为主差应力变小,所以水泥环进入塑性的风险减小。
如果把水泥环当作孔隙介质,水泥环应该具有一定的有效应力系数,图10给出了水平最大地应力方位改变BIOT系数时水泥环区域周向和轴向应力分布规律,可以看出,当增大水泥环BIOT系数时,水泥环的应力急剧降低,当两者中的较小应力小于水泥环抗拉强度时会产生拉伸破坏,所以,在评价水泥环完整性的时候有必要准确地预测其有效应力系数。
5 结果与讨论
(1)可由本文中的数学模型求解水泥环系统的应力和位移分布,其中由于套管和水泥环失重,其轴向应力分布表达式不同于地层;
(2)在最小地应力方位下,水泥环内壁处更容易进入塑性,在水泥环内壁45°方位上有产生滑移的可能性,水泥环由于应力状态为压应力,产生拉伸破坏的可能性较小;
(3)如果把水泥环当作孔隙介质,水泥环产生拉伸破坏的可能性增加,所以在评价水泥环完整性时有必要准确地预测其有效应力系数;
(4)为了防止水泥环产生拉伸破坏或者进入塑性,水泥环应具有低模量,高抗拉和抗压强度,高泊松比的特性,另外为了防止界面产生滑移,应提高第一和第二界面的摩擦系数,为此要保证良好的固井质量;
(5)综合可以看出:水泥环弹性模量对应力分布的影响相比于泊松比要大,并且合理的水泥环力学参数选取应该综合考虑套管-水泥环-地层系统的力学参数匹配性和以及该系统所承受的内压和原始地应力相对关系。
参考文献
[1] Zinkham R E,Goodwin R J.Burst resistance of pipe cemented into the earth.SPE 291,1962
[2] Goodwin K J,Crook R J.Cement sheath stress failure.SPE20453,1992
[3] Thiercelin M J,Dargaud B,Baret J F,et al..Cement desing based on cement mechanical response.SPE 52890,1998
[4] 李军,陈勉,张辉,等.水泥环弹性模量对套管外挤载荷的影响分析.石油大学学报(自然科学版),2005;29(6):41—44Li Jun,Chen Mian,Zhang Hui,et al.Effects of cement sheath elastic modulus on casing external collapse load.Journal of the University of Petroleum,China,2005;29(6):41—44
[5] 徐守余,李茂华,牛卫东.水泥环性质对套管抗挤强度影响的有限元分析.石油钻探技术,2007;35(3):5—8Xu Shouyu,Li Maohua,Niu Weidong.Finite element analysis of effect of cement sheath property on casing collapsing strength.Petroleum Drilling Techniques,2007;35(3):5—8
[6] 杨桂通.弹塑性力学引论.北京:清华大学出版社,2004:44—49Yang Guitong.Elastic-plastic mechanics introduction.Beijing:Tsinghua University Press,2004:44—49
[7] 王敏中,王玮,武际可.弹性力学教程.北京:北京大学出版社,2002Wang Minzhong,Wang Wei,Wu Jike.Elastic mechanics.Beijing:Peking University Press,2002
[8] 徐芝纶.弹性力学.第4版.北京:高等教育出版社,2006:73 —76Xu Zhilun.Elasticity.4th edition.Beijing:Higher Education Press,2006:73—76
非完整井 篇7
关键词:移动机器人,非完整约束,运动规划,人工免疫算法
1、引言
伴随着机器人技术的进步与发展,移动式机器人的应用越来越广泛。移动式机器人的移动基通常由三轮或四轮式结构组成。设轮子与地面之间只有纯滚动接触,因而存在非完整束。由于非完整约束是不可积的,系统的控制及运动规划问题将变的困难和复杂。Murray和Sastry[1]采用向量场李代数的分析方法给出了规划问题的解。Mukerjee和Anderson[2]利用Green公式得到了运动规划问题的解。Fernandes等[3]运用变分方法研究了非完整运动规划问题,其方法的思想来源于泛函分析中的Ritz近似理论。文献[4]利用遗传算法讨论了带有非完整约束的航天器系统姿态运动规划的最优控制问题。文献[5]在移动机器人的运动规划中引入了遗传算法。人工免疫系统[6]以自然界中的生物系统的免疫功能为机理,从生物界中得到启示,并将其中的有用成分应用到人工智能中。它作为一种新型的、模拟免疫系统的随机化搜索和优化方法,近几年来在优化领域得到广泛的研究和应用,并已在解决诸多优化问题中显示了良好的性能和效果。
本文利用人工免疫算法,讨论轮式移动机器人非完整运动规划问题。由人工免疫算法取代传统的牛顿迭代方法,解决轮式移动机器人非完整运动规划最优控制中的寻优与搜索问题。首先将系统的非完整约束方程转化为非线性控制系统的状态方程,利用人工免疫算法搜索和寻求最优控制输入,从而得到系统的非完整运动状态转换的优化轨迹。通过数值仿真表明,该方法能有效地求解轮式移动机器人的运动规划问题。
2、系统模型与运动规划问题
2.1 数学模型
设移动机器人的移动基由三个轮子组成,其中前轮为导向轮,后二轮为驱动轮,俯视图如图1所示。以O为原点建立惯性坐标系xoy,在移动基驱动轮轴线上选择一点P为原点建立连体坐标系x′Py′,则移动基的位形可以用广义坐标向量q=(x,y,θ)T描述,其中(x,y)和θ分别表示移动基在惯性坐标系中的位置和方位。为了保证P点的速度方向与x′轴保持一致,设导向轮的转动轴始终固定平行y′轴,驱动轮的转动轴为自由的。当轮子在地面上作纯滚动时,P点的速度与移动基的方位角θ之间存在下列约束关系:
显然,此约束条件是不可积的,因此式(1)为轮式机器人移动基的非完整约束方程,因此,该机器人系统为非完整约束系统。从物理意义上说,机器人非完整约束的含义为:在纯滚动无滑动的约束条件下,机器人在任何时刻都不能沿着轴向作侧向移动。
非完整系统的运动规划问题可以归结为设计适当的控制输入使系统沿某一轨线从初始位置到另一目标位置。因此本质上是两点边值问题,由于其线性化系统是不可控的,所以要得到边值问题的解是相当困难的。本文利用最优控制方法和免疫进化算法,寻找满足边值条件的解。
结合机器人的运动特性,由式(1)得到机器人的数学模型如下式:
其中,VL,VR分别表示左右轮的速度,L为左右轮之间的距离。
记于是将上述规划问题归结为选取适当的1r和2r使得:
该系统为非线性系统,将式(3)中的ri(i=,1)2取作控制输入变量,记作u,u=(r1,r2)T。定义移动基系统的位形q=(x,y,θ)T为状态变量,系统的状态方程由式(2)可表示为:
2.2 运动规划的最优控制描述
根据最小能量控制原理,选择性能指标函数为:
其中,u(t)为Hilbert空间L2([.0T])的可测向量函数,实际计算时,只需考虑有限维数情况,利用Fourier基向量作为正交基向量{ai}Ni=1,则u(t)可以表示为正交基向量的线性组合:
其中为函数u在基上的投影,并且:
将α视为新的控制变量,引入惩罚因子λ,并考虑Fourier基向量的正交性,指标函数(5)可以写为:
其中,q(T)是方程(4)在给定控制输入u时系统在t=T时的状态,显然q(T)是α的函数。设q(T)=F(α),当N和λ给定时,式(8)变为:
这样寻找控制输入u使式(5)为极小值问题就转化为寻找α使式(9)为极小值问题。
3、最优控制的人工免疫算法
3.1 抗体的编码形式
在求解最优控制问题时,采用二进制编码的人工免疫算法需要将连续的空间离散化,若要提高解的精度,就必须增加个体长度,而且在算法的计算过程中还需要不停地进行编码和解码操作,计算时间比较长。而基于实数编码的算法不存在编码和解码过程,能够提高解的精度和运算速度。因此本文算法采用实数编码。
对式(6)中函数u在Fourier基上的投影α采用实数编码,抗体编码α为αi(k)(i=1,2,…,N)组成的N维向量。有:
其中,αi(k)为实数,变量k为代数。
3.2 亲合力函数的建立
由于人工免疫算法中的每一步免疫操作都是基于亲和力函数来进行路径抗体群体更新的,所以亲和力函数的选取有着重要的地位。人工免疫算法的最终目的是找出亲合力高的抗体作为最优解,所以亲合力函数构造如下:
其中,α(k)为抗体,J[α(k)]为式(9)的目标函数。
3.3 算法的具体实现步骤
算法的操作步骤如下:
(1)抗原识别;输入目标函数和各种约束条件作为免疫算法的抗原。
(2)产生初始抗体群;利用式(10)产生n个实数编码的抗体α作为初始群体A(k),在第一次迭代时,抗体通常在解空间中用随机的方法产生。
(3)抗体选择操作;根据式(11)计算群体A(k)中每个抗体的亲合力函数值,并且计算抗体的浓度,基于抗体浓度的概率公式选择其中m个抗体组成群体B(k)(其中m
并且令0q=0。在[01,]区间上产生一个均匀分布的随机数w,若qk-1≤w≤qk,则抗体αk被选中,产生m个这样的随机数,则有m个抗体被选中到群体B中。
(4)抗体在趋同半径内进行扩增操作;群体B(k)中m抗体在趋同半径QR局部范围内进行扩增操作,产生n个抗体组成群体C(k)。群体中的任一抗体αi的局部范围QutongR(αi)构造为:
其中||·||为欧几里德范数,Ω为解空间,QutongR(αi)是由与αi的欧氏距离不大于常数QR的所有可行解构成,解空间是以αi为中心,QR为半径的局部球形区域,定义QR为趋同半径。
(5)抗体在异化半径内进行突变操作;群体C(k)中n-m个抗体在异化半径YR全局范围内进行突变操作,保持群体C(k)的规模不变并组成群体D(k)。群体中的任一抗体αi的全局范围Yihua R(αi)构造为:
其中YihuaR(αi)是以αi为中心,YR为半径的全局球形区域,定义YR为异化半径。
(6)抗体更新;群体D(k)中亲合力最低的m1个抗体替换成步骤3中的免疫记忆抗体,组成新一代群体A(k+1)。
(7)终止条件;如果新一代群体A(k+1)中有满足要求的抗体,则输出最优抗体,算法结束;否则转到第3步重复执行。
4、仿真实验
在仿真实验中,设定抗体群规模n=30,m=10,m1=5,抗体长度N=10。选取10个Fourier正交基组成2×10的Φ矩阵,其中{ai(t)}5i=1为:
由上式各个基矢量行轮换得到,式中k=2π/T。移动基完成规定时间设为T=5s。
设系统初始状态为q0=(-,4-,4π/)4T,要求移动基到达目标位形为qf=(.2,5.3,5π/)3T。仿真结果如图2、3、4所示,为移动基运动的优化轨线。通过100代的人工免疫算法计算得到目标函数最优值J(u)=10.59445112,其误差为ER=.0001014。
5、结语
将人工免疫算法其引入非线性系统最优控制中,这是一种有益的尝试。通过仿真计算表明,人工免疫算法对对于解决轮式移动机器人非完整运动规划的最优控制问题是有效和可行的,使用实数编码尽管在理论上没有二进制编码成熟,但是对于多参数优化问题,如果从计算精度考虑,它还是一种合适的选择。与传统的迭代算法相比,人工免疫算法不需要导数信息,使用范围广泛,而且是一种全局优化算法。本文的人工免疫算法也为其它最优控制问题提供了新的思路与途径。
参考文献
[1]Murray R,Sastry S.Steering nonholonomic systems systems us-ing sinusoids[A].In:Proc.of 29th Conf.on Decision and Control[C],1990,2097-2109.
[2]戈新生,张奇志.航天器太阳阵展开过程最优控制的遗传算法[J].力学季刊,2000,21(1):134-138.
[3]厉虹,胡兵.轮式移动机器人非完整运动规划的遗传算法[J].自动化技术与应用,2005,24(2):13-15.