潘公宇 刘一
摘要:为了实现无人驾驶汽车的路径跟踪和避障功能,基于模型预测轮廓控制提出一种路径跟踪控制策略。首先将二自由度单轨车辆模型和Pacejka轮胎模型结合在一起搭建整车模型,遵循模型预测的思路将路径规划和跟踪问题结合为一个非线性优化问题,利用线性化得到线性时变模型,在每个采样时间内以凸二次规划的形式建立非线性优化问题的局部逼近。其次,为了使车辆在跟踪路径的时候维持稳定,结合车辆稳定包络线,施加横向稳定性约束,并通过施加车道线边界约束,实现自动驾驶车辆的避障功能。最后将Matlab/Simulink和Carsim结合进行联合仿真。仿真结果表明,所设计的控制器在较宽路面能自主规划出一条最短路径,实现了跟踪和避障功能,同时维持了车辆在行驶中的稳定性和安全性,在赛道竞速下具有较明显的优势。研究结果可为结构化道路下的汽车自动驾驶提供技术参考,具有一定的应用价值。
关键词:车辆工程;模型预测;轮廓控制;车辆避障;路径跟踪
中图分类号:U463.06;TP273文献标识码:ADOI: 10.7535/hbgykj.2021yx04003
Automatic vehicle path tracking control based on model
predictive contouring control
PAN Gongyu, LIU Yi
(School of Automobile and Traffic Engineering, Jiangsu University, Zhenjiang, Jiangsu 212013,China)
Abstract:In order to realize the function of tracking and obstacle avoidance of unmanned vehicles, a path tracking control strategy based on model predictive contour control was proposed. Firstly, the 2-DOF single vehicle model and the Pacejka tire model were combined to build the vehicle model, the path planning and tracking problems were combined into a nonlinear optimization problem by following the idea of model prediction, and the linear time-varying model was obtained by linearization.The local approximation of the nonlinear optimization problem was established in the form of convex quadratic programming in each sampling time. Then, in order to maintain the stability of the vehicle when tracking the path, the lateral stability constraint was applied in combination with the vehicle stability envelope. By applying lane boundary constraint, the obstacle avoidance function of the autonomous driving vehicle was realized. Finally, Matlab/Simulink and CarSim were combined to conduct co-simulation. The simulation results show that the designed controller can independently plan a shortest path on a wide road surface, realize the function of tracking and obstacle avoidance, and maintain the stability and safety of the vehicle in the running, which has a good advantage in the race track.The results provide technical reference for vehicle autonomous driving on structured roads and has a certain application value.
Keywords:vehicle engimeering; model predictive; contouring control; obstacle avoidance; path tracking
目前智能化是汽車发展不可逆转的趋势,智能化的无人驾驶汽车是一项重大的技术创新。无人驾驶汽车的发展可以有效降低交通事故和缓解交通拥堵,并且随着智能化程度的不断提高,无人驾驶汽车必将产生更大的社会效益。
无人驾驶汽车通过布置在车上的感知系统获取周围环境,由控制器规划出路径,控制车辆沿着规划出的路径行驶,实现无人驾驶的目的。目前无人驾驶控制器是两级控制,一级为路径规划层,另一级为跟踪控制层。文献[1]采用改进粒子群算法,提高了轨迹规划后期的收敛速度和效率,并解决了陷入局部最优解的问题。CHEN等[2]用贝塞尔曲线参数的约束问题规划无人驾驶车辆的路径规划问题,并采用全局搜索能力的遗传算法对其进行求解,从而得到满足约束条件的避让可行路径。LI等[3]在路径规划中采用Sigmoid函数,并且将Sigmoid函数参数与汽车侧向加速度、侧向加速度变化率和避障距离约束条件建立映射关系,使无人驾驶车辆获得了较好的安全性和舒适性。在跟踪控制中,吴艳等[4]使用滑膜控制设计跟踪策略,但滑膜控制的抖振问题降低了跟踪的平顺性。WNAG等[5]采用μ综合方法的跟踪策略,能在极端情况下保持稳定,但是跟踪性能较差。在使用模型预测控制(MPC)方法,对外界干扰具有较强的鲁棒性,MATRAJI等[6]采用超扭算法的二阶滑模控制来消除机器人轨迹跟踪中的抖振。SAHA等[7]利用前馈控制和PID控制来跟踪所需的偏航率和所需的航向角。由于两层的控制策略设计是先得到规划路径再根据规划路径得出控制量,跟踪性能依赖规划的路径,因此跟踪控制无法对路径规划的结果产生影响,从而无法得到最优的控制量。
模型预测控制方法能系统地考虑预测信息和处理多约束能力[8],因此被认为是处理多点预瞄跟踪控制的有效方法。TCHAMNA等[9]结合电子控制系统和悬架控制系统,考虑悬架系统的非线性问题建立模型,通过横摆角速度和侧偏角的约束,实现在稳定跟踪的同时减小侧倾角,改善了车辆的平顺性和舒适性。但是由于模型过于复杂,增加了MPC控制器的求解负担,求解效率低下。LINIGER等[10]对比了避障和规划两种MPC控制架构,分层式架构通过上层规划模块获取一条避障路径,通过下层MPC控制器进行轨迹跟踪,提高了MPC的求解速度;单层式架构直接在MPC控制器的目标函数中加入障碍物信息,由于使用一个预测模型,跟踪效果更优秀,但是计算负担较大。BROWN等[11]用前轮侧向力作为MPC控制器的输入,在目标函数中增加车辆稳定包络线约束和环境包络线约束,实现局部的避障功能,但是跟踪的参考路径要依赖上层路径规划层的输入。
根据以上研究现状,提出了基于模型预测轮廓控制[12]的无人驾驶汽车跟踪控制策略。第4期潘公宇,等:基于模型预测轮廓控制的自动驾驶车辆路径跟踪控制河北工业科技第38卷
1构建车辆模型
车辆二自由度单轨模型被广泛应用于轨迹跟踪控制器的设计,因此该控制器采用具有非线性轮胎力的单轨模型,单轨模型如图1所示。
该单轨模型具有以下假设:1)车辆在平整路面上行驶;2)忽略载荷转移;3)忽略滑移;4)纵向驱动力作用在重心。运动学方程如式(1)所示:
X˙
Y˙
φ˙
v˙x
v˙y
r˙=vxcos φ-vysin φ
vxsin φ+vycos φ
r
1mFFysin δ+mvyr
1mFRy+FFycos δ-mvxr
1IzFFylFcos δ-FRylR。(1)
式中:m和Iz分别表示车辆的质量和转动惯量;lR和lF分别表示质心至后轮和前轮的距离;FRy和FFy分別表示后轮和前轮的侧向力;Fx是作用在质心的驱动力。模型状态量为x=X,Y,φ,vx,vy,rT。式中(X,Y)和φ表示全局坐标系中的质心位置和航向角;vx和vy表示纵向和侧向速度;r表示横摆角速度。控制输入u=[δ],其中δ是前轮转向角。
为了满足计算精度和计算要求之间的平衡,选择简化的Pacejka轮胎模型[13],公式如下:
αR=arctan(vy-lRrvx),
αF=arctan(vy+lFrvx)-δ,
FRy=DRsin(CRarctan(BRαR)),
FFy=DFsin(CFarctan(BFαF)),(2)
式中:BR,CR,DR,BF,CF,DF为通过实验确定的模型参数;αR,αF分别表示后轮和前轮的侧偏角。
2构建轮廓误差
为了获得高性能的自动驾驶控制器,采用模型预测轮廓控制算法,因为其能在预测范围内的参考路径上行驶最远距离。该方法的优点是可以将路径规划和路径跟踪结合成一个非线性问题,在每个采样时间用局部凸二次规划逼近非线性优化问题进行实时求解[14]。在构建模型预测轮廓控制算法之前,首先要进行参考路径参数化和轮廓误差计算。
2.1建立约束
参考路径是由左右车道线计算出车道中心线上的点经过三阶样条多项式拟合得到弧长θ∈[0,L],其中L表示总长度。通过求解三阶多项式,可以得到参考路径上的任意点Xrefθ,Yref(θ)。参考路径上点的切线与大地坐标系X轴的夹角可以用式(3)得到:
ΦθarctanYrefθXrefθ,(3)
这种参数化在参考路径已知点内的插值更加精确。
2.2测量轮廓误差
为了建立模型预测轮廓问题,需要误差测量来定义车辆当前位置X,Y和期望参考点Xrefθ,Yref(θ)的偏差[14]。P: 瘙 綆 2→[0,L]为参考路径上的投影算子,用来确定车辆已经行驶过的参考路径长度,表示为
P(X,Y)
argminθ((X-Xref(θ))2+(Y-Yref(θ))2)。(4)
为了简洁,定义θPP(X,Y)。车辆位置(X,Y)至参考轨迹的正交距离可以得到式(5):
ecX,Y,θP
sinΦθPX-XrefθP-
cosΦθPY-YrefθP,(5)
式中Φ·在式3中定义,ec为轮廓误差,描述如图2所示。
投影算子不适合在在线优化算法中使用,因为它本身就是一个优化问题,投影算子是参考路径上距车辆最近的点,得到投影算子需要从预测范围内参考路径起点开始搜索,增加了计算的负担。因此,引入一个投影算子θP的逼近投影算子θA,这是一个由控制器决定的自变量,为了使θA有用,通过一个滞后误差el连接θA到θP。滞后误差el如图3所示。
elX,Y,θAθA-θP。(6)
为了确保θA是真实弧长的一个很好的近似值,引入了一个成本函数,使沿着路径的误差(滞后误差e︿l)和垂直参考路径的误差(轮廓误差e︿c)最小化,如图3所示。根据参考路径的几何关系,e︿l和e︿c可以表示为式(7)和式(8):
e︿cX,Y,θAsinΦθAX-XrefθA-
cosΦθAY-YrefθA,(7)
e︿lX,Y,θA-cosΦθAX-XrefθA-
sinΦθAY-YrefθA,(8)
式中ΦθA是与参考路径相切的角度。对e︿l进行约束,在求解时使e︿l趋近于0,由此θA也将趋近于θP,得到e︿c≈ec,e︿l≈el,因此只需将e︿l抑制在一个较小的范围内,就不需要再对投影算子进行搜索。
2.3约束车辆稳定包络线
当车辆发生横向运动的时候,非常容易发生甩尾和侧滑等危险工况,所以在设计控制器时需要对车辆的状态进行实时监控,防止车辆失稳。但是实时监控会增大控制器的计算量,增加计算负担。车辆状态的质心侧偏角β和横摆角速度r对车辆的稳定性有着至关重要的作用。BOBIER等[15]提出的车辆稳定包络线如图4所示。只要车辆质心侧偏角和横摆角速度控制在包络线范围内,就可以保证其行驶的稳定性。β-r边界线如式(9)—式(13)所示。
CD:β=a0r-a1,
AB:β=a0r+a1,
BC:r=rmax,
AD:r=rmin。(9)
a0=lRvx,(10)
a1=tanαr,peak,(11)
rmax=μgvx,(12)
rmin=-μgvx。(13)
由于后轮的载荷大,后轮相比于前轮更容易饱和,因此,稳定性边界由后轮力峰值决定。参考Pacejka轮胎模型[11],后轮最大侧偏角公式如式(14)所示:
αr,peak=arctan3mgμCRlFlF+lR,(14)
式中μ是路面摩擦系数。将r代入CD和AB段,最大质心侧偏角βmax和最小质心侧偏角βmin可以表示为
βmax=lRvxr+tanαr,peak,(15)
βmin=lRvxr-tanαr,peak,(16)
式中β=vyvx,将vx设定为一个常数,可以将式(15)和式(16)改写为关于vy和r的不等式:
-αpeak≤vy-lRrvx≤αpeak。(17)
当车速增加或者轮胎磨损与路面的摩擦减弱,在执行转向指令时后轮容易达到非线性区域,在这种情况下轮胎力容易饱和,进而发生甩尾或侧滑,影响驾驶的安全性。为了使车辆状态一直维持在稳定包络线中,可以使用一个线性不等式约束来实现:
Hvxk|≤Gv,k=1,2,…,N,(18)
式中
Hv=000001
00001vx-lRvx,(19)
Gv=μgvx
αpeak。(20)
2.4避障
除了规划稳定的轨迹外,控制器还利用改变车道线边界约束的方法,在每个预测步长时间内,将预测的车辆轨迹约束在改变后的边界内,确保车辆在行驶中不会与障碍物发生碰撞。同时在没有障碍物时,由于车道线边界的存在可以确保车辆不会驶出道路边界。调整车道线边界如图5所示。
图5说明了调整车道线边界的方法。通过对预测步长路径上的离散点采样来确定每个k中道路边界和无障碍物之间的安全区域,如图5 b)所示。将环境中的目标长度扩展到该离散采样长度,识别出障碍物之间大于车辆宽度的可行间隙,将预测范围内连续的安全区域结合在一起,形成无障碍物的通过走廊,如图5 c)所示。对于每个通过走廊,側向上的边界定义了安全区域,如下:
Ymink+12w+fb 式中:w是汽车宽度;fb是障碍物或者道路边界距车辆的最小距离。将车辆的横向位置Y约束在调整后的边界内,达到控制车辆行驶在车道线内和实现避障的目的。这个不等式可以用式(22)来表示: Hexk≤Ge,(22) 式中He=[010000],Ge为通过走廊上边界和下边界的值。 3构建目标函数 结合车辆的动力学模型(1)、参数化参考路径、轮廓误差,构建模型预测轮廓控制算法,目的是在路径跟踪的质量(低轮廓误差)和N个采样时间的有限范围内取得的预测量之间进行权衡。为了减轻计算复杂度,在上述轮廓误差中提到的用θA代替θP。由于在预测范围内的每个时间步长形成滞后误差,有必要引入动态积分状态器θA,k+1=θA,k+vkTs,其中vk为投影速度,θA,k为在K时刻的状态,Ts为采样时间。构建最优控制形式的非线性优化问题,如式(23)所示: min∑Nk=1‖e︿c,kXk,Yk,θA,k‖2qc+‖e︿l,kXk,Yk,θA,k‖2ql-γvkTs+‖Δuk‖2Ru+‖Δvk‖2Rv-γv0Ts, s.t.θ0=θ, θA,k+1=θA,k+vkTs,k=0,…,N-1, 0≤θk≤L,k=1,2,…,N, 0≤vk≤vmax,k=0,1,…,N,23 式中:Δukuk-uk-1和Δvkvk-vk-1。 ∑N-1k=0vkTs=θP,N,θP,N是车辆位置Xk,Yk在参考路径上正交投影点Xref(θP,k),Yref(θP,k)的弧长。uk是控制输入量。θk和vk的边界值是为了避免出现错误的非线性问题的解。θ0为初始值。γ∈ 瘙 綆 ++和qc∈ 瘙 綆 ++是权重值。为了确保精确的近似,从而使成本函数和汽车模型之间有强耦合ql∈ 瘙 綆 ++,选择较高的权重值。为了惩罚快速的变化,加入一个关于输入变化的成本项,有助于获得平滑的控制输入,用来防止未建模的动力学问题放大。 为了实时求解非线性最优控制问题(23),通过对非线性项进行线性化,在每个采样时间建立具有以下二次规划形式的局部凸逼近,如式(24)所示: min∑Nkxk θA,kTΓkxk θA,k+cTkxk θA,k-γvkTs+Δuk ΔvkTRΔuk Δvk-γv0Ts+qvσv+qeσe, s.t.xo=x,θA,0=θ, xk+1=Akxk+Bkuk+gk,k=0,…,N-1, θA,k+1=θA,k+vkTs,k=0,1,…,N-1, 0≤θA,k≤L,k=1,2,…,N, xmin≤xk≤xmax,k=1,2,…,N, umin≤uk≤umax,k=0,1,…,N, 0≤vk≤vmax,k=0,1,…,N, Hvxk≤Gv+σv,k=0,1,…,N, Hexk≤Ge+σe,k=0,1,…,N, (24) (a) (b) 式中Γk∈S7+表示线性化轮廓线和滞后误差成本函数的二次项部分,ck∈ 瘙 綆 7为线性部分。 为了保证线性化误差足够小,使用线性时变估计的车辆动力学(1)以及轮廓和滞后误差。测量值x=x0为第一个线性化点。对一个时间步长的非线性化模型进行仿真,计算线性化点的终端状态xN。将控制輸入uk限制在物理约束范围内并且将状态变量xk限制在合理范围内,以避免求解器的收敛问题。最后,添加车辆稳定包络线和车道线约束,式(a)和式(b)分别为车辆稳定包络线约束和车道边界线约束,使用非负的松弛变量σv和σe,以确保优化问题总是可行的。权重qv和qe分别设置为5和3×104,这样即使松弛变量较小,松弛变量项依然能影响成本函数。 4仿真和结果分析 对弯道转弯进行仿真,在Matlab/Simulink平台上搭建所设计的控制器,在Carsim中搭建B级车模型,结合在一起进行联合仿真,验证所提出的弯道路径跟踪系统的有效性。实验中使用50 ms的采样时间和40个时间步长的预测范围,这相当于超前2 s观测。所使用的车辆参数如表1所示。 4.1弯道路径跟踪仿真 为了体现本文控制器的自主规划能力,对设计的控制器进行仿真,通过输入车道线参数规划出最短路径。为有足够的空间进行规划,采用较宽的道路,因此选用有连续弯道并且路面较宽的赛道地图,宽为16 m,路面的摩擦系数设为0.85,车辆速度为10 m/s。在Matlab中导入左右车道线的信息,将左右车道线的信息输入给所设计的控制算法,以未施加稳定包络线的控制器做对比进行仿真。 路径跟踪的仿真结果如图6所示,红色为施加稳定包络约束的控制器,蓝色为对比控制器,黑色为车道边界线,虚线为车道中心线,①为第1个弯,②为第2个弯。明显可以看出红色轨迹线在连续弯道中的转弯半径更大,从而能使车辆在经过连续弯道时获得更好的自动驾驶性能,红色轨迹路程长度为85.04 m,蓝色轨迹路程长度为99.44 m,2个连续弯道行驶路程缩短了14.48%,也更加符合驾驶员入弯和出弯的驾驶轨迹。由于施加了车道线边界约束,红色轨迹在获得更优轨迹的同时和车道边界线保持了安全距离,防止车辆驶出道路。 横摆角速度结果如图7所示,图中①和②分别为2个控制器在第1个弯和第2个弯时横摆角速度达到最大。由于施加稳定包络约束的车辆入弯更早,达到最大横摆角速度时间也会更早。通过第1个弯时红色曲线到达峰值的时间比黑色曲线提前1.7 s,在4.8 s至5.4 s时横摆角速度超过包络线值,控制器控制转向角变化,将横摆角速度推回包络线内,转向角的变化如图8中的①所示。通过第2个弯是红色曲线到达峰值比黑色的曲线提前了1.88 s,在6.28 s至7.65 s时,红色曲线的峰值超过包络线,控制器控制转向角变化将横摆角速度推回包络线内,转向角的变化如图8中的②所示。横摆角的峰值减小了62.5%,使车辆的稳定性大幅提高。 β-r包络线如图9所示,可以明显看到上下边界为车辆通过弯道时横摆角速度被控制在包络线之内。质心侧偏角和横摆角速度被约束在β-r包络线中,保持了车辆在连续弯道中的稳定性。 4.2避障仿真 本节对设计的控制器进行避障仿真,选用有障碍物的直线赛道路面,赛道宽为16 m,路面的摩擦系数设为0.85,车辆速度为10 m/s。 车辆避障轨迹如图10所示,图中有3个障碍车辆,蓝色曲线为车辆轨迹,很好地避开了障碍物,在横坐标为15 m处车辆躲避障碍车辆1,在5 m处车辆躲避障碍车辆3。 横摆角速度如图11所示,在22 s车辆躲避障碍车辆1,横摆角速度达到波谷,但是没有超过边界值。23 s至24 s车辆躲避第2个障碍物,控制器将横摆角速度控制在边界值内,避免了车辆发生失稳情况。前轮转向角如图12所示,在23 s至24 s时控制器控制转向角减小,将横摆角速度推回包络线中。β-r包络线如图13所示,躲避障碍车辆1时,横摆角速度被控制器控制在包络线内。质心侧偏角和横摆角速度被约束在β-r包络线中,使车辆在避障中保持了稳定。 5结语 1)采用模型预测轮廓控制将路径规划和跟踪控制结合在一起的一层控制策略,通过对控制器输入车道线和车辆状态,实现跟踪功能,消除了跟踪层性能对规划层结果的依赖,能在行驶中根据实际情况实时调整跟踪轨迹,提高了跟踪效果;模型预测结构允许约束处理和抑制潜在干扰,模型预测成本函数解决了跟踪速度和准确性之间的权衡,以获得最佳的行驶路径;采用线性时变公式来降低计算复杂度;考虑到车辆在行驶中的横向稳定性,施加车辆稳定包络线约束。 2)为了在单层控制器中实现避障功能,将车道线边界施加在约束中,通过改变车道线边界实现了无人驾驶汽车的避障功能。 仿真结果表明,控制器具有良好的自主规划能力,在较宽路面的连续弯道能沿最短路径行驶,为了获取更大的转弯半径,在靠近车道线行驶时,由于车道线边界的约束,没有超出车道线,并且很好地保持了行驶稳定性。避障仿真表明,本文控制器通过改变车道线边界约束,能较好地躲过障碍物,并且没有发生失稳情况。 控制器在较宽路面能自主规划出一条最短路径,并实现跟踪和避障功能,在赛道竞速下具有很好的优势,实际道路宽度较窄,虽然也能实现跟踪功能,但体现不出本文控制器的优势,是以后需要改进的地方。 參考文献/References: [1]杜云, 刘冰, 邵士凯,等. 基于改进粒子群算法的无人机航迹规划[J]. 河北工业科技, 2019, 36(5):335-340. DU Yun, LIU Bing, SHAO Shikai, et al. UAV flight path planning based on improved particle swarm optimization[J]. Hebei Journal of Industrial Science and Technology, 2019, 36(5):335-340. [2]CHEN Long,QIN Dongfang,XU Xing,et al.A path and velocity planning method for lane changing collision avoidance of intelligent vehicle based on cubic 3-D Bezier curve[J].Advances in Engineering Software,2019,132:65-73. [3]LI Shaosong,LI Zheng,YU Zhixin,et al.Dynamic trajectory planning and tracking for autonomous vehicle with obstacle avoidance based on model predictive control[J].IEEE Access,2019,7:132074-132086. [4]吴艳,王丽芳,李芳.基于滑模自抗扰的智能车路径跟踪控制[J].控制与决策,2019,34(10):2150-2156. WU Yan,WANG Lifang,LI Fang. Intelligent vehicle path following control based on sliding mode active disturbance rejection control[J].Control and Decision,2019,34(10):2150-2156. [5]WNAG Chunyan,ZHAO Wanzhong,XU Zhijiang,et al.Path planning and stability control of collision avoidance system based on active front steering[J].Science China Technological Sciences,2017,60:1231-1243. [6]MATRAJI I,AL-DURRA A,HARYONO A,et al.Trajectory tracking control of Skid-Steered Mobile Robot based on adaptive Second Order Sliding Mode Control[J].Control Engineering Practice,2018,72:167-176. [7]SAHA J,BEST M,BENMIMOUN A,et al.Autonomous rear-end collision avoidance using an electric power steering system[J].Proceedings of the Institution of Mechanical Engineers,Part D:Journal of Automobile Engineering,2015,229(12):1638-1655. [8]BORRELLI F,FALCONE P,KEVICZKY T,et al.MPC-based approach to active steering for autonomous vehicle systems[J].International Journal of Vehicle Autonomous Systems,2005,3(2/3/4):265-291. [9]TCHAMNA R,YOUN E,YOUN I.Combined control effects of brake and active suspension control on the global safety of a full-car nonlinear model[J].Vehicle System Dynamics,2014,52(sup1):69-91. [10]LINIGER A,DOMAHIDI A,MORARI M.Optimization-based autonomous racing of 1:43 scale RC cars[J].Optimal Control Applications and Methods,2015,36(5):628-647. [11]BROWN M,FUNKE J,ERLIEN S,et al.Safe driving envelopes for path tracking in autonomous vehicles[J].Control Engineering Practice,2016,61:307-316. [12]LAM D,MANZIE C,GOOD M.Model predictive contouring control[C]//49th IEEE Conference on Decision and Control(CDC).[S.l]:[s.n.],2010:6137-6142. [13]PACEJKA H B,BAKKER E.The magic formula tyre model[J].Vehicle System Dynamics,2007,21(sup1):1-18. [14]DIEHL M,BOCK H G,SCHLDER J P,et al.Real-time optimization and nonlinear model predictive control of processes governed by differential-algebraic equations[J].Journal of Process Control,2002,12(4):577-585. [15]BOBIER C G,GERDES J C.Staying within the nullcline boundary for vehicle envelope control using a sliding surface[J].Vehicle System Dynamics,2013,51(2):199-217.