传统工业园区核心企业的能源网络结构简单,主要表现在能源使用结构较单一, 如:钢铁企业[1]、水泥企业[2]及卷烟企业[3]的能源使用结构几乎全部依赖于国家电网及化石能源,可再生能源应用比例较低。此外,园区核心企业还存在能源信息封闭、开放互联程度不高的问题。核心企业能源供给系统大多基于其自身生产过程的冷、热、电的需求状况进行配置[4-5],这样会导致由企业生产波动而引起的冷、热、电盈余与缺口的问题难以及时高效地得到解决,企业单位产值的能源成本高,园区企业之间能源网络系统未能较好地协同与互补,能源系统使用效率低,大部分园区上下游企业之间只有产业链的依存关系,缺少企业能源信息交流与合理交易。
在当今“工业4.0”与ICT信息技术大趋势推动下,使用物联网将能源生产与使用结合起来,辅以能源调度技术是实现智慧用能、智能生产的重要途径。建设智慧工业园区是当今数字化时代必然趋势,其中智慧节能工业园区 (SEEIP) 建设的关键是园区核心企业的能源网络的构建与能源系统的合理配置。文中针对传统工业园区核心企业的能源网络存在的不足,提出了一种改进的核心企业局域智慧能源网络框架,并基于此框架给出能源网络系统动态的优化配置方法。以期为SEEIP核心企业的能源系统优化配置与传统工业园区核心企业的智慧化改造提供一种切实可行的工程方法。
1 核心企业能源网络 1.1 能源网络框架孙彦广[6]对企业能源系统采用模块化分析的方法,将能源系统分为能源供应模块、能源转换模块和能源需求模块,提出的企业能量流网络一般模式,如图 1所示。能量流网络模型中外购能源、自产能源、副产能源为能量流始端节点,基本用户、外销用户为能量流终端节点,储存单元、转换单元构成能量流中间缓冲系统,管网为能量流连结器。
![]() |
图 1 企业能量流网络一般模式 Figure 1 The general mode of enterprise energy flow network |
传统核心企业能源网络不但没有体现ICT信息技术和“可再生能源”的利用而且也没有体现与工业园区其它企业之间的能量互联与互补关系。核心企业能源供给独立于园区其它企业会造成能源系统容量配置过大,引发投资与运行费用偏高,企业能源消费成本高,因而企业生产产品的价格缺乏国际竞争力。
因此,以“ICT技术”和提高“可再生能源”使用比例为标志的核心企业局域智慧能源网络是我国工业园区发展的必然模式。能源的生产、传输、分配、转换、存储、消费等环节通过互联互通的能源网络与信息平台[7]进行协调控制,可大幅度提高能源利用效率,结合不同能源形式的“替代和转化”,充分利用可再生能源,改变园区能源消费结构,实现节能减排目标。
核心企业局域智慧能源网络如图 2所示。
![]() |
图 2 核心企业局域智慧能源网络 Figure 2 The local smart energy network of core enterprise |
图中egin和ecout为进入和离开能源生产模块的能量,egin和ecout为进入和离开能源转换模块的能量,esin和esout为进入和离开能源蓄存模块的能量。
一次能源介质按照其所属能源形式可以划分为国家电网、化石能源 (煤、石油、天然气等)、可再生能源 (风能、潮汐能、太阳能、地热能、生物质能等) 及余热余能 (余热烟气、饱和蒸汽、焦炉煤气、转炉煤气等)4大类。
二次能量介质主要是从产能设备生产出来进入能源转换设备使用的电与热能介质。
三次能量介质主要是从能源转换设备生产出来进入能源蓄存设备的冷、热、电能介质。
末端用户的用能需求包含供暖、供冷、工艺加热、电力、生活热水等需求。
基于ICT技术的能源信息共享平台起到信息反馈与能源调度,以管理企业生产与园区其它企业之间的能源交易,能源交易可发生在企业能源生产环节和转换环节。“虚拟设备”只是为了便于对企业所有能量介质的流向进行分析而提出的一个概念,“虚拟设备”对于进入其中的介质并没有经过任何处理而直接进入了下一环节,介质的物化性质没有任何变化。所以此网络框架中的各能量处理模块并不是传统意义的能量处理环节,是包含该企业所有能量流的一个广义处理环节。
1.2 能源网络数学模型 1.2.1 能源生产模块能源生产模块包含众多能源生产设备,如:燃油发电机、内燃机、燃气轮机、锅炉、风力发电机、太阳能光伏发电、太阳能集热板、余热回收设备等,能源交易出口在此被划入“虚拟设备”的范畴。
在能源生产模块中,各种一次能源介质通过能源生产设备被转化为电能和热能,被二次能量介质所携带。同种一次能源介质的实时消耗总量可由下式计算:
em1,ing, t=∑eq∂eq×em1,ineq, g, t, | (1) |
式中:e为一次能源介质m1逐时总消耗量;g代表能源生产环节;t表示任意时刻;m1为具体的一次能源介质 (包括煤、石油、天然气等形式); in代表进入能源生产环节; eq代表各能源生产设备,包含图 2所示的“虚拟设备”。∂eq为二元变量,其取值为1或者0,取值为1代表选择了这种设备,取值为0代表没有选择。
各种能源生产设备产出的二次能源介质m2(如电、热) 的量能可用下式描述:
em2,outg, t=∂eq×ef, ineq, g, t×ηeq, m2, | (2) |
式中:m2代表产出的二次能源介质的形式,主要包括电能及热能,m2{p, h}, p表示电能,h表示热能;ηeq, m2为一次能源介质转化为二次能源介质的效率,可由设备厂家提供。显然,当忽略输送损失时,经过“虚拟设备”的各种能源介质的转换效率为1。
1.2.2 能源转换模块能源转换模块包含换热器、吸收式制冷机组、电热器、水源热泵等设备,能源交易入口在此被划入“虚拟设备”的范畴。
在能源转换环节二次能源介质 (电能、热能) 通过能源转化设备转换为三次能源介质 (冷能、热能、电能)。同种类型及品位的二次能源介质m2总消耗量可由下式描述:
em2,inc, t=∑eq∂eq×em2,ineq, c,t, | (3) |
其中,c表示能源转换环节。
各种能源转换设备产出的三次能源介质m3(如冷、热、电) 的产量可以下式描述:
em3,outc, t=∂eq×em2,ineq, c, t×ηeq, m3, | (4) |
式中:m3代表产出的三次能源介质主要包括冷能、热能及电能,m3{c, h, p};c表示热能;h表示热能;p表示电能;ηeq, m3为二次能源介质转化为三次能源介质的效率。同样,二次能源介质通过“虚拟设备”的效率也近似为1。
1.2.3 能源蓄存模块能源蓄存模块包含冷、热、电三大类的能源蓄存设备。常见的有冰蓄冷设备、储热罐及蓄电池。
任意时刻进入蓄能模块的能源介质大部分是供给用户直接使用的,这部分能源介质相当于被该模块的“虚拟设备”所处理;任意时刻流入的能源介质中多余用户需求的那部分则进入能源蓄能设备, 而不足用户需求的那部分可由蓄能设备在前面时刻蓄存的能量进行补充。用户的需求波动由能源蓄存模块进行调节,起到削峰填谷的作用,便于应对产能波动引起的能源供需波动变化,提高能源使用效率。
各种类型及品位的三次能源介质m3总消耗量可由下式描述:
em3,ins, t=∑eq∂eq×em3,ineq, s, t, | (5) |
式中:s表示能源环节;es, tm3, in为逐时进入蓄存装置的三次能源产量;设备的能源消耗量同样要考虑“虚拟设备”消耗的三次次能源介质量。
em3,outs, t=∑eq∂eq×em2,outeq, s, t, | (6) |
式中:es, tm3, out为逐时从蓄存装置出来的三次能源产量;设备的能源产出量同样要考虑“虚拟设备”产出的三次能源介质量。
由于蓄能与释能的能力与上一时段和这一时段内蓄能装置内所蓄存的冷 (热) 量有关,是一个动态过程,因此蓄能装置采用动态数学模型[8-9],可由下式表示:
定义
Pm3eq, s, t=βt⋅em3,ineq, s, t−(1−βt)⋅em3,outeq, s, t, | (7) |
Em3eq, s, t+1=Em3eq, s, t(1−μ)Δt+{Δt⋅δ⋅Pm3eq, s, t,Pm3eq, s, t≥0Δtδ⋅Pm3eq, s, t,Pm3eq, s, t≥0 | (8) |
式中:Em3, t为t时刻内蓄能装置所蓄存的冷热量;Δt为t时刻和t+1时刻的时间间隔;Pm3, t为t时刻蓄放功率,规定蓄存为正,释放为负;βt为能量蓄存状态的二元变量,取值为1代表选择了蓄存能量,取值为0代表释放能量;δ为冷 (热) 量传递过程中的有效系数,取0.98;μ为蓄能装置向环境散冷 (热) 损失系数,取0.02。
显然,根据实际情况对于蓄能周期T (一般蓄能周期T取1天) 的时段内任意时刻t,蓄能装置所蓄存的冷 (热) 量上限不大于蓄能装置的设备容量,下限不小于0,同时工程上一般规定每小时内蓄能或释能总量不超过蓄能装置容量的0.3倍; 并且蓄能装置1个工作周期的初始和终止时刻蓄存能量相等。有如下关系:
{0≤Em3,t≤Em3,volume−0.3≤3600⋅Pm3,tEm3,max≤0.3Em3,0=Em3,T=C | (9) |
式中:Em3, volume为蓄能设备的容量;C为1个周期的初始和终止时刻蓄能设备里面蓄存的能量,一般取为蓄能装置容量的70%。
2 优化模型的建立核心企业能源系统配置问题是在满足用户能量种类、品位以及数量需求的前提下,考虑经济性、节能环保及当地的能源资源条件下,能源生产设备、能源转换设备、能源蓄存设备组合配置的优化问题。
核心企业能源系统能源动态调度问题是在已知能源系统组合配置的前提下,在一个相对稳定的生产周期内,考虑经济性、节能环保及当地的能源资源条件,动态的优化能源生产、转化、蓄存设备的启停情况以满足用户能量种类、品位以及数量需求。
2.1 目标函数为了便于分析,将不同能源系统配置对环境的影响转换为经济指标,用碳税来进行量化。
核心企业能源系统配置问题选择年平均费用作为该优化的目标函数。年平均费用主要包括投资费、运行维护费、能源使用费、能源交易费用及碳税这5个部分,费用最低的能源网络系统作为最优配置方案[8]。
核心企业能源系统能源动态调度问题是在一个相对稳定的计划生产周期T0内选择生产周期总费用作为该优化的目标函数。生产周期总费用主要包括运行维护费、能源使用费、能源交易费用及碳税这4个部分,总费用最低的能源设备调度方案作为最优配置方案。
minCtotal=Ccapital+CO & M+Cfuel+Ctrade+CCO2 | (10) |
minCtotal, T0=CO & M+Cfuel+Ctrade+CCO2, | (11) |
其中:
Ccapital=Cfacility+Gpipe=∑eq∂eq×Neq×Feq×Peq×I01−(1+I0)−Teq+∑iεi×Li×pi,pipe×I01−(1+I0)−Ti,pipe, | (12) |
CO & M=∑eq∂eq×PO & M, eq×∑esεes×PO & M, es×Pes+CInternet, | (13) |
Cfuel=∑t(em1,ing, t×Pm1,t+egp, t×Pgp,t), | (14) |
Ctrade=∑t(em1,inin,t×Pm3,t), | (15) |
CCO2=ECO2×PCO2, | (16) |
式中:minCtotal为最低年平均费用;Ccapital为初投资费用;CO&M为操作维护费用;Cfuel为能耗费用;Ctrade为能源交易费用;CCO2为碳税;Cfacility为设备的投资费用;Gpipe为能源管网投资费用;CInternet为能源信息平台的操作维护费用;εi,εes,∂eq均为二元变量其取值为0或者1;Neq为设备台数;Feq为设备单位容量的投资费用;Peq为单台设备的容量;Teq为设备的寿命周期;I0为设备的年折旧率,一般取0.1;Li为管道的长度;Pi, pipe为单位长度的管道的投资费用;PO&M, eq为能量生产与转化设备单位能量操作维护费用;PO&M, es为能量蓄存设备单位容量操作维护费用;Pm1, t、Pgp, t为逐时能源价格和电价; Pm2, t表示能源逐时买卖交易价格;ECO2, t为某时刻t碳总排放量;PCO2为碳税价格。
2.2 约束条件 2.2.1 能源供需约束条件能源的生产量应不小于能源需求量,按照核心企业局域智慧能源网络的架构,必须满足以下3个层次的供需约束:
1) 能源蓄存模块能源的产生量大于等于各用户的使用需求量,满足下式
eh, outs, t≥Dh, t, | (17) |
ec, outs, t≥Dc, t, | (18) |
ep, outs, t≥Dp, t, | (19) |
式中:es, th, out,es, tc, out,es, tp, out表示从蓄存模块供出能源量;Dh, t,Dc, t,Dp, t为用户的需求量。
2) 能源转换模块能源的产生量大于等于能源蓄存模块的能源需求量,满足下式
eh, outc, t≥eh, ins, t, | (20) |
ec, outc, t≥ec, ins, t, | (21) |
ep, outc, t≥ep, ins, t, | (22) |
式中:ec, th, out, ec, tc, out, ec, tp, out为能源转换模块供出的能源量;es, th, in, es, tc, in, es, tp, in为进入蓄能模块的能源量。
3) 能源生产模块能源的产生量大于等于能源转换模块的能源需求量,满足下式
eh, outg, t≥eh, inc, t, | (23) |
ec, outg, t≥ec, inc, t, | (24) |
ep, outg, t≥ep, inc, t, | (25) |
式中:eg, th, out,eg, tc, out,eg, tp, out为能源生产模块供出的能源量;ec, th, in,ec, tc, in,ec, tp, in为进入能源转换模块的能源量。
能源生产和转化的各个环节的能量供入和产出基本满足下面的关系
eoutt=eint×ηeq, | (26) |
式中:etin为供给能源生产和转换环节的能源量;etout为产出能源生产和转换环节的能源量;ηeq为能源生产和转换环节的效率。自末端用户需求的开始逐级约束前一环节技术设备的选择。
2.2.2 能源设备容量约束条件所有设备都要求运行在效率较高的区间运行,满足下式
Pmineq≤eouteq, t≤Pmaxeq, | (27) |
蓄能设备容量取决于蓄能周期内最大蓄能量
Pes≥maxEm3,t。 | (28) |
核心企业局域智慧能源网络优化配置的决策变量是设备类型、设备容量和台数,即:优化的目标函数是关于二元变量∂eq、设备容量Peq及台数Neq的复合函数。能源网络优化模型公式 (10)~式 (28) 中,各参数均是关于∂eq、Peq及Neq的函数,该优化问题可描述为MIN Ctotal=f(∂eq,Peq,Neq) 的混合整数非线性规划 (MINLP) 问题。通过收集现有的各技术设备的样本信息,建立“设备类型”及“设备容量”的数据库,访问数据库选取设备,读取相关技术参数。
核心企业能源系统能源动态调度问题与前者类似,决策变量同样是使用设备的类型、设备容量和台数,只是设备选择约束在企业能源系统已有的设备范围内。
求解这两个问题方法相同,这里以核心企业局域智慧能源网络优化配置问题为例,给出求解MINLP问题的算法。
3.1 MINLP问题的数学描述核心企业局域智慧能源网络优化配置问题的数学函数如下
minf(∂eq,Peq,Neq)。 | (29) |
约束条件如下:
hk(∂eq,Peq,Neq)=0,k=1,2,⋯,me;gl(∂eq,Peq,Neq)≤0,l=1,2,⋯,mc;Neq∈N;Pmineq≤Pwq≤Pmaxeq∂eq=0或1; |
式中, me和mc分别为等式约束和不等式约束个数;∂eq为二元整数变量;Peq,Neq为离散变量。
改进的差异演化算法 (DE) 是在连续空间进行启发式随机搜索的一种优化算法,具有求解效率高、收敛速度快、全局搜索能力强的特点。差异演化算法不宜直接处理有约束的优化问题, 通常采用2种方法对上述规划问题约束条件进行处理[10]:一是限制解在可行域内,二是罚函数法。文中采用拉格朗日二项式罚函数法[11]将目标函数转化为无约束问题。
罚函数参数及朗格朗日乘数对算法的收敛性影响较大,取值过小造成罚函数惩罚过轻,约束条件得不到满足;取值过大造成惩罚过量使算法陷入局部最优解。因此,采用线性自适应策略来调整罚函数各参数[12]在进化初期各参数取较小值,进化到后期加大惩罚力度,强化局部搜索。
3.2 MINLP问题的求解算法差异演化算法 (DE) 要求在连续空间进行优化计算, 是一种求解实数变量优化问题的有效方法。用DE求解混合整数规划问题, 必须对DE进行改进。DE的基本操作包括变异、交叉和选择操作, 依据适应值大小进行操作。根据DE算法的特点, 只要对变异操作进行改进就可以将DE用于整数规划和混合整数规划。对于整数变量对变异后的矢量进行取整运算, 这样使得变异操作可以在实数域进行, 从而扩大了寻优空间, 有利于提高算法的寻优能力。DE在求解混合整数规划问题的基本操作步骤如下:
1) 生成初始种群。
2) 变异操作。差异演化算法最基本的变异成分是父代的差分矢量, 每个矢量对包括父代2个不同的个体。根据变异个体的生成方法不同, 形成了多种不同的差分进化算法方案[13]。混合序列差异演化算法 (HLUDE) 的基本思想是根据各个体的函数值进行升序排列, 排在前面的个体比较接近最优值, 排在后面的个体则相对远离最优值。因此, 对排在前面的个体进行局部搜索, 三角变异操作 (TMO)[14]是一种相当贪婪的局部搜索操作;对排在后面个体实施大范围的搜索, 这样做既可保证快的收敛速度又可增加全局搜索能力;Rainer Storn和Kenneth Price提出的DE/rand/1/bin是一种全局搜索能力最强的变异模式[15],对排在前面个体实施该搜索, 可保证快速的局部搜索能力。
3) 交叉操作。交叉操作可以保持种群的多样性。
4) 选择操作。比较试验个体与原个体的适应值大小,只有试验个体的适应度值较优时,才会被选为子代,否则直接将原个体作为子代。
5) 终止条件。多数的迭代算法都是预先设置一最大的进化代数, 只要迭代过程达到预设的值, 算法就结束。鉴于差异演化的变异操作是基于染色体的差异向量进行的, 随着迭代的进行, 群体差异化变小, 进化到群体间无差异, 也就没必要再继续迭代。文中用欧氏距离来判别群体间差异化的大小。只要进化代数达到最大的进化代数或子代中任意2个个体的欧氏距离均小于指定的欧氏距离算法结束。
3.3 数学模型求解求解MINLP的改进的DE算法流程如图 3所示。
![]() |
图 3 算法流程图 Figure 3 The algorithm flowchart |
求解MINLP的改进的DE算法流程如下:
Step1:初始化种群规模,缩放因子,交叉概率。在每个变量定义域内随机初始化每个个体。确定罚函数参数以及调节因子,设置最大迭代次数及计算精度要求,置当前迭代次数t=0;
Step2:计算每个个体每个约束条件的惩罚量;
Step3:计算每个个体的适应值,求出最优适应值以及最优个体。
Step4:判断是否达到精度要求或者是否达到最大迭代次数,若是则退出,否则执行下一步。
Step5:对每个个体执行6~8步骤生成第t+1代种群;
Step6:在种群中随机选择3个不同的个体与待变异个体,按照变异操作公式生成变异个体;
Step7:进行交叉操作,生成试验个体
Step8:进行选择操作,生成t+1代个体;
Step9:t=t+1, 返回step2。
4 结语在工业“4.0”及能源信息化背景下,针对SEEIP建设及既有工业园区核心企业能源系统的智慧化改造,通过研究得到一些可用于SEEIP核心企业的能源系统优化配置与能源调度的工程解决方法。
1) 提出一种基于互联网平台的局域智慧能源网络框架,通过ICT技术的能源信息共享平台实现能源信息反馈、能源调度与能源交易,以便于实现能源供给需求的品位对口。
2) 给出了能源生产模块、能源转换模块与能源蓄存模块中各种能源介质及品位的能量计算数学模型。
3) 考虑经济成本与环境影响,以投资费、运行维护费、能源使用费、能源交易费用及碳税这5个部分组成的年平均费用最低为目标,给出SEEIP企业能源网络系统配置的优化模型。
4) 采用改进的差异演化算法 (DE) 求解优化模型,并给出了求解步骤与程序流程,通过ICT信息共享平台可动态的给出企业能源系统最佳配置与供给方案,以满足用户冷、热、电动态需求。虽然目前提出了该计算方法,但还需通过应用于实际工业园区实践以验证其实用性与可靠性。
[1] |
蒋育翔, 黄全福, 洪小和.
钢铁企业能源及相互关系研究[J]. 经济问题探索, 2011(1): 125–129.
JIANG Yuxiang, HUANG Quanfu, HONG Xiaohe. Analysis of energy in steel industry and its interrelation[J]. Inquiry into Economic Issues, 2011(1): 125–129. (in Chinese) |
[2] |
中国水泥协会. 中国水泥年鉴[M]. 北京: 中国建材工业出版社, 2015.
China cement association. China cement clmanac[M]. Beijing: China Building Materials Press, 2015. (in Chinese) |
[3] |
鲁端峰, 闫亚明, 刘晓萍, 等. 卷烟企业能源消耗现状及节能措施探讨[C]//中国环境科学学会2013年学术年会. 北京: 中国环境科学出版社, 2013. LU Duanfeng, YAN Yaming, LIU Xiaoping, et al. The research of cigarette enterprises' consumption situation and energy saving measures[C]//The 2013 Academic Conference of China Environmental Science Association. Beijing: China Environmental Science Press, 2013. (in Chinese) |
[4] |
黄文燕. 钢铁企业能源管理系统及其能耗预测的研究[D]. 广州: 华南理工大学, 2013. HUANG Wenyan.The research of the steel enterprises energy management system and its energy consumption prediction[D]. Guangzhou: South China University of Technology, 2013. (in Chinese) |
[5] |
谢骆乐. 基于模块化设计理论的区域能源规划模型研究[D]. 重庆: 重庆大学, 2013. XIE Luole. Research on the model of regional energy planning based on modularization design theory[D]. Chongqing: Chongqing Uninversity, 2013. (in Chinese) |
[6] |
孙彦广. 钢铁企业能量流网络信息模型及多种能源介质动态调控[C]//全国能源环保生产技术会议. 北京: 中国金属学会, 2010. SUN Yanguang. The network information model of steel industry and a variety of energy medium dynamic regulation[C]//The Technical Conference of Energy and Environment Production. Beijing: The Chinese Society for Metals, 2010. (in Chinese) |
[7] |
周海明, 刘广一, 刘超群.
能源互联网技术框架研究[J]. 中国电力, 2014, 47(11): 140–144.
ZHOU Haiming, LIU Guangyi, LIU Chaoqun. Study on the energy Internet technology framework[J]. Electric Power, 2014, 47(11): 140–144. (in Chinese) |
[8] | Zhou Z, Liu P, Li Z, et al. An engineering approach to the optimal design of distributed energy systems in China[J]. Applied Thermal Engineering, 2013, 53(2): 387–396. DOI:10.1016/j.applthermaleng.2012.01.067 |
[9] | Yokoyama R, Ito K. Effect of inlet air cooling by ice storage on unit sizing of a gas turbine cogeneration plant[J]. Journal of Engineering for Gas Turbines and Power, 2002, 126(2): 335–342. |
[10] | Lampinen J. A constraint handling approach for the differential evolution algorithm[C]//Proceedings of the 2002 Congress on Evolutionary Computation. [S.l.]: IEEE, 2002: 1468-1473. |
[11] | Hamida S B, Schoenauer M. ASCHEA: new results using adaptive segregational constraint handling[C]//Proceedings of the 2002 Congress on Evolutionary Computation. [S.l.]: IEEE, 2002: 884-889. |
[12] | Huang H J, Wang F S. Fuzzy decision-making design of chemical plant using mixed-integer hybrid differential evolution[J]. Computers and Chemical Engineering, 2002, 26(12): 1649–1660. DOI:10.1016/S0098-1354(02)00098-4 |
[13] | Babu B V, Jehan M M L. Differential evolution for multi-objective optimization[C]//The 2003 Congress on Evolutionary Computation. [S.l.]: IEEE, 2003:2696-2703. |
[14] | Fan H Y, Lampinen J. A trigonometric mutation operation to differential evolution[J]. Journal of Global Optimization, 2003, 27(1): 105–129. DOI:10.1023/A:1024653025686 |
[15] | Sarimveis H, Nikolakopoulos A. A line up evolutionary algorithm for solving nonlinear constrained optimization problems[J]. Computers and Operations Research, 2005, 32(6): 1499–1514. DOI:10.1016/j.cor.2003.11.015 |