含多种分布式电源的弱环配电网三相潮流计算_丁明.pdf

第 29 卷 第 13 期 中 国 电 机 工 程 学 报 Vol.29 No.13 May 5, 2009 2009 年 5 月 5 日 Proceedings of the CSEE 2009 Chin.Soc.for Elec.Eng. 35 文章编号：0258-8013 (2009) 13-0035-06 中图分类号：TM 72 文献标志码：A 学科分类号：47040 含多种分布式电源的弱环配电网三相潮流计算 丁明，郭学凤 (合肥工业大学教育部光伏系统工程研究中心，安徽省 合肥市 230009) Three-phase Power Flow for the Weakly Meshed Distribution Network With the Distributed Generation DING Ming, GUO Xue-feng (Photovoltaic System Research Center of Ministry of Education, Hefei University of Technology, Hefei 230009, Anhui Province, China) ABSTRACT: Traditional power flow algorithms are not applicable to the distribution network with distributed generation because of the unbalance in power source, load and line. Based on the analysis of wind turbines, photovoltaic system, fuel cell, storage battery, high frequency micro turbines and power frequency combined heat and power (CHP), their models in power flow calculation were presented. In view of unbalance in three phases and weakly meshed characteristics of the distribution network, the back/forward sweep algorithm in power flow calculation was improved, which was widely used for the radial distribution network. A new three phase back/forward sweep power flow calculation for the weakly meshed distribution network with the distributed generation was presented. The example has proved the efficiency and feasibility of the algorithm. KEY WORDS: distributed generation; three-phase power flow; weakly meshed distribution network 摘要： 由于分布式电力系统中电源、 负荷和线路参数存在不 对称性， 传统的配电网潮流算法已不适用于这类网络。 在对 风力发电机组、光伏发电系统、燃料电池、蓄电池、高频微 型燃气轮机及工频热电联产机组的运行和控制性能详细分 析的基础上， 建立各种分布式电源潮流计算模型， 考虑配电 网三相参数不平衡和环网问题， 对广泛用于辐射型配电网的 前推回代潮流算法进行改进， 提出适用于含多种分布式电源 的弱环配电网的前推回代三相潮流算法。 算例表明， 该文提 出的算法是有效和可行的。 关键词：分布式电源；三相潮流计算；弱环配电网 基金项目：国家自然科学基金项目(50607002)；国家 863 高技术基 金项目(2007AA05Z240)；教育部科学技术研究重大项目(306004)。 Project Supported by National Natural Science Foundation of China(NSFC)(50607002); The National High Technology Research and Development of China 863 Program(2007AA05Z240). 0 引言 近年来电力系统中出现了大量的不对称负荷， 例如电力机车、电弧炉和单相整流设备，配电网采 用单相三相混合供电模式，不对称(如单相)分布式 电源引入配电网也加重了三相潮流的不平衡1。此 时若仍按对称情况来计算配电网潮流分布，必将产 生较大误差。配电网一般是环形设计，开环运行， 但在配电网实际运行和优化计算时，需要计算分段 开关和联络开关闭合时的弱环网潮流。因此，研究 弱环配电网三相不对称潮流的计算模型是非常必 要的。 分布式电源引入配电网，对网损、电压的分布 及有功功率和无功功率的数量及传输方向均有很 大影响2-3。文献4假定通过风电机组和光伏发电 系统中电容器的自动投切，可使功率因数恒定，在 潮流计算中，简单地将它们视为 PQ 节点。文献5 建立了基于异步发电机等值电路的风力发电机组 PQ(V)稳态模型；文献6采用 R-X 迭代模型，将风 电场作为一个以 R-X 表示的阻抗接在母线上； 但都 只分析了传统的异步风机。 文献7建立了异步发电 机、无励磁调节能力的同步发电机和燃料电池在潮 流计算中的模型，并改进直接法，提出了基于灵敏 度补偿的配电网潮流计算方法。 文献8将常见的各 种分布式电源归结为 PV 节点、PI 节点或 PQ(V)节 点，提出了基于牛顿法的配电网三相潮流计算方 法。但以上方法均未考虑环网问题。 含分布式电源的潮流计算的重要特点是分布 式电源模型与传统发电机组模型可能不同。在对风 力发电机组、光伏发电系统、燃料电池、蓄电池、 36 中 国 电 机 工 程 学 报 第 29 卷 高频微型燃气轮机及工频热电联产机组的运行和 控制性能详细分析的基础上，本文建立各种分布式 电源潮流计算模型，考虑配电网三相参数不平衡和 环网问题，对广泛用于辐射型配电网的前推回代潮 流算法进行改进，提出适用于含多种分布式电源的 弱环配电网的前推回代三相潮流算法。给出算例和 结果。 1 改进前推回代法 前推回代法9充分利用配电网辐射型的结构特 点，具有较高的计算效率和优良的收敛特性，但不 能处理环网。传统配电网一般仅包含 2 种节点：平 衡节点和 PQ 节点。变电站出口母线通常视为平衡 节点， 而其它负荷节点和中间节点都视为 PQ 节点。 分布式电源引入配电网的同时也引入了 PV 节点。 鉴于前推回代法不能处理环网和 PV 节点，本文改 进前推回代法， 提出弱环网和 PV 节点的处理方法。 图 1 为一个有 m 个环网和 n 个 PV 节点的配电 网， 若将 PV 节点和环网(一般环网由联络开关和分 段开关闭合时形成， 所以在开关处打开)打开， 出现 m+n 个开环点,存在如下关系： = Z IU (1) 假定所有节点电压标幺值近似等于 1.0，并且 相角很小的情况下，有 * = Z SU (2) 推导可得 = XRQU RXP (3) 式中：Q 和U 的维数为 m+n；P 和 的维数为 m。因为对 PV 节点P 始终为 0，在下次迭代进行 开环网潮流计算前，将Q 和P 叠加到相应的节点 即可。注意叠加时流出功率的节点取正，流入功率 的节点取负。 i1 i2 im j1 j2 jn 图 1 配电网 Fig. 1 A distribution network 2 分布式电源的潮流计算模型 2.1 燃料电池 目前的燃料电池10可用于分布式发电的有 4 种：磷酸燃料电池、固体氧化物燃料电池、熔融碳 酸盐燃料电池和质子交换膜燃料电池。 图 2 为燃料电池的结构图，由 3 部分组成11： 预处理装置、燃料电池堆和并网逆变装置(power conditioning unit, PCU)。单个燃料电池的直流电压 约为 0.7 V， 由一定数量的燃料电池并联成燃料电池 堆产生较高的适于应用的电压， PCU 将燃料电池堆 出口电压转化为交流电压。燃料电池堆的交流输出 电压和输出功率可视为调制度 m 和相角 的函数： ACFC = UmU (4) ACFCssin( )/ PmU UX= (5) 堆 化学能 转化为 电能 空气 天然气 丙烷 甲醇 汽蚀 预处理装置 燃料电池 氧气 氢气 PCU 空气 负荷 + + + IL Us UAC UPC X I 图 2 燃料电池的结构图 Fig. 2 Structure of fuel cell 假设逆变器无能耗，燃料电池堆交流输出电压 的相角和氢气流速之间存在下述关系： 2 H s0 2FvX q mU N = (6) 式中：PAC为逆变器交流侧的输出功率；X为线路 电抗；Us为负荷侧电压；UFC为燃料电池直流输出 电压；qH2为氢气流量；N0为燃料电池堆中所并联 的燃料电池数目；F为常数；v为氢气的利用率。 由式(6)知控制氢气的流量可控制， 由式(5)知控制 可控制燃料电池的有功功率输出；PCU通过调节m 控制交流输出侧的电压幅值，所以在潮流计算中把 燃料电池视为PV节点。 2.2 光伏发电系统和储能系统 光伏发电系统12一般只给电网提供有功功率， 即将太阳能光伏阵列的直流电能转换为与电网同 频同相的交流电能馈送给电网。图3为光伏并网发 电的控制框图。 光伏并网多采用电压源型电流控制逆变器。太 阳电池阵列工作电压决定输出功率，要输出最大功 太阳电 池阵列市网 负载 电网电压同步信号 MPPTAVRACR DC/AC 变换器 IDC UDC If + + DC U p I I 图 3 光伏并网发电的控制框图 Fig. 3 Diagram of PV control and grid interconnection 第 13 期 丁明等： 含多种分布式电源的弱环配电网三相潮流计算 37 率，必须先稳定工作电压为最大功率点电压。最大 功率跟踪(maximum power point tracking, MPPT)用 于确定最大功率点电压 DC U。电压调节控制单元 (automatic voltage regulating, AVR)输出并网电流幅 值 P I给定。 P I与电网电压同步信号合并，与负荷 电流比较，电流调节控制单元(automatic current regulating, ACR)输出DC/AC变换器的控制信号， 控 制DC/AC变换器输出与电网电压同频同相的正弦 电流，并网功率因数为1。 蓄电池的控制方法与光伏发电系统基本相同13。 区别是蓄电池既可作为电源，又可作为电网的负 载。当电路工作在整流状态下，能量从电网侧流向 直流侧的蓄电池，对蓄电池充电；当电路工作在逆 变状态下，将直流侧的能量回馈给电网，蓄电池对 电网放电。 本文所用前推回代法是基于支路电流的， 用节点 注入电流表示各种类型的负荷。 既然光伏发电系统和 蓄电池并网控制的是输入电网的电流， 在潮流计算中 将它们视为向电网输入电能的恒电流负荷(I节点)。 2.3 微型燃气轮机热电联产 微型燃气轮机的功率较小，一般在25100 kW 范围内，对微型燃气轮机热电联产来说，所产生的 功率既有热功率又有电功率，运行模式可是“以热 定电”或“以电定热” 。微型燃气轮机通过对转速 的控制来控制功角和有功功率输出。微型燃气轮机 电力电子设备接口由3部分组成14：AC/DC整流 器、带电容的直流母线和DC/AC变换器。微型燃 气轮机输出三相可变频/可变电压， 通过电力电子设 备将可变电压/可变频率转换为固定电压/固定频 率，向负荷供电。所以微型燃气轮机在潮流计算中 处理为PV节点。 2.4 工频热电联产机组 工频热电联产机组即传统的异步发电机和同 步发电机。同步发电机的励磁控制系统在电力系统 的运行中起着重要作用，控制发电机的端电压、无 功功率和功率因数等参数。励磁调节器的控制方式 有： 自动电压调节、 恒无功调节和恒功率因数调节。 自动电压调节时工频热电联产机组在潮流计算中 处理为PV节点，恒无功调节和恒功率因数调节时 处理为PQ节点。 异步发电机本身没有励磁装置，主要靠电网提 供的无功功率建立磁场，吸收的无功功率与机端电 压有关，在潮流计算中既不能作为PQ节点，也不 能作为PV节点，需要特殊考虑。 图4为异步发电机的近似等效电路， 其中U为 发电机的节点电压幅值，Is为定子电流，IR为转子 电流，Im为励磁电流，R为转子电阻，s为转差率， Re为机械负载等效电阻，xm为励磁电抗，x为发电 机定子电抗与转子电抗之和，Pe为异步发电机有功 功率。有功功率一定的情况下，发电机转差率和无 功功率为 U Is IR Im xm x Re=R/s + 图 4 异步发电机的近似等效电路 Fig. 4 Approximate model of asynchronous generator 2422 e 2 e (4) 2 R UUxP s P x = (7) 22 m e m ()Rxxxs QP Rx s + = (8) 可看出在有功功率一定的情况下，异步发电机 无功功率与电压的变化密切相关，在潮流计算中作 为PQ(V)节点处理，算法在每次迭代时都利用上次 迭代所得电压对功率进行更新。 2.5 风力发电机 风力发电机组大致分为4种15， 其拓扑结构如 图5所示。 图5(a)、(b)的定速型和转差控制型异步风机， 与传统的异步发电机没有本质区别，区别是风力发 电机组通过电容器组的自动投切对无功功率就地 补偿，不再靠电网提供的无功功率建立磁场。所以 与异步发电机的处理方法类似，均处理为PQ(V)节 点，此时参与潮流迭代的节点注入无功Q QQ =， Q 为电容器组输出无功， Q 为风机吸收的电网无 功。算法在每次迭代时都利用上次迭代所得电压对 功率进行更新。 图5(c)的同步直驱型风机与电网完全解耦，采 用全功率变流器实现风机的全范围调速。通过控制 全功率变流器满足电力系统的有功和无功需求，在 潮流计算中视为PQ节点。 图5(d)的双馈型风机(doubly-fed wind power generator，DFIG)16，定子并网，通过与电网相连 双向变流器控制转子励磁电流的频率、幅值和相 角。DFIG并网多采用交流励磁变速恒频控制。此 控制策略下风机的运行包括启动区、风能追踪区、 38 中 国 电 机 工 程 学 报 第 29 卷 齿轮箱 软起 装置 无功 补偿 (a) 定速型 齿轮箱 附加电阻 控制 无功 补偿 (b) 转差控制器 (c) 同步直驱 齿轮箱 变流器 齿轮箱 变流器 (d) 双馈型 图 5 风力发电机的主要机型拓扑结构 Fig. 5 Topology of typical wind turbine 转速恒定区和功率恒定区4大区域。4个区域内 DFIG的控制模式不同： 启动区需调节DFIG输出电 压，以符合并网要求，采用并网控制；其他3个运 行区域内DFIG发电运行，要求DFIG输出电能可 控，能实现有功、无功的解耦控制。并网控制是对 并网之前的DFIG定子输出电压进行调制，使其与 电网电压在幅值、频率和相位上一致，以抑制并网 时的冲击电流。并网控制时DFIG脱网空载运行， 没有功率输出；当DFIG定子电压和电网电压相同 时，DFIG并网空载运行，输出功率仍为0。可认为 并网控制是功率解耦控制的一个特例(P=0，Q=0)， 在潮流计算中将DFIG作为PQ节点处理。 3 含多种分布式电源的弱环配电网前推回 代三相潮流算法的实现 本文采用VC+ + 6.0建立节点、线路、变压器、 联络开关、PV节点、PQ节点、I节点和PQ(V)节 点等的类，编制了含多种分布式电源的弱环配电网 的前推回代三相潮流算法的程序。具体步骤如下： 1）从数据库中读取数据。检测开关的状态判 断环网和PV节点的数目，形成矩阵 XR RX 。 2）初始化网络电压。 3）计算节点注入电流。I节点直接叠加即可。 PV节点、PQ节点和PQ(V)节点转化为恒功率负荷 与普通PQ节点一起进行处理。PV节点的无功初始 化为(Qmax+Qmin)/2，Qmax和Qmin为PV节点的无功 上下限。工频热电联产同步发电机(自动电压调 节)PV节点的无功功率上下限由同步发电机自身决 定，燃料电池和微型燃气轮机无功功率上下限由电 力电子变换器参数决定。 4）电流前推，电压回代，检查两次迭代间节 点电压变化量， 如满足所定收敛条件， 内循环结束； 否则，返回步骤3）继续下一次迭代。 5）计算开环点电压差U 和 、开环点电压差 的变化量。根据式(3)算出修正功率Q 和P，按照 叠加原则叠加到相应的节点上。检查开环点电压差 的变化量，如满足收敛条件，外循环结束；否则， 返回步骤3）继续下一次迭代。 6）计算结束。 4 算例分析 本文对33节点配电系统进行分析，该系统的 三相网络参数见文献17，存在明显的三相负荷和 支路参数不对称的情况。原系统有5个环网，本文 选用了其中的2个(33和35支路)。电网高压侧额 定电压为12.66 kV，低压侧为0.38 kV。本文为对分 布式电源接入弱环三相配电网的性能进行分析，拟 定了11个方案，如表1所示。 图6为接入分布式电源后的网络接线图，其中 的变压器变比为12.66/0.38，归算到低压侧的阻抗 为0.01+j0.04 ，Y0/Y0接线。给出了利用本文所提 出算法得到的计算结果，表2为各方案部分节点电 压(A相)，表3为各方案所需迭代次数。 由表2知：除了方案9比方案1的各节点电压 有所降低外，其余方案的各节点电压都比方案1的 高；方案5的各节点电压高于方案3和4；方案10 的各节点电压高于方案6和8； 方案11的各节点电 压高于方案2、3和10。工频热电联产异步发电机 22 21 20 19 213456789 10 11 12 13 14 15 16 17 18 2324 25 3635 2627 2829 30 32 31 33 39 40 3438 37 41 图 6 含分布式电源的配电网 Fig. 6 Distribution network with distributed generation 第 13 期 丁明等： 含多种分布式电源的弱环配电网三相潮流计算 39 表 1 方案描述 Tab. 1 Description of schemes 方案 接入分布式电源类型 接入分布式电源参数 接入节点 1 无 2 PQ(DFIG) P=300 kW, Q=100 kvar 36 3 PV(燃料电池) P=500 kW, U=12.66 kV 39 4 PQ(V)(定速型异步风机) P=1 000 kW, s=0.033 35 5 方案 3+方案 4 6 I(三相光伏发电系统能或蓄电池) IA=IB = IC =10 A, A=0, B =120, C =240 34 7 I(单相光伏发电系统能或蓄电池) IA =10A, A =0 34 8 PQ(工频热电联产同步发电机(恒功率因数调节) P=400 kW，cos =0.85 37 9 PQ(V)(工频热电联产异步发电机) P=300 kW，s=0.033 38 10 方案 6+方案 8 11 方案 10+方案 2+方案 3 注：电容器QN-Unit=40 kvar，UN=12.66 kV；电力电子变换器PN=600 kW，min=0.75。 表 2 潮流计算结果(节点电压) Tab. 2 Node voltage of power flow kV 节点号 方案 1 方案 2 方案 3 方案 4 方案 5 方案 6 方案 7 方案 8 方案 9 方案 10 方案 11 1 12.66 12.66 12.66 12.66 12.66 12.66 12.66 12.66 12.66 12.66 12.66 2 12.648 3 12.649 2 12.65 12.65 12.651 612.648 312.648 312.649 612.642 4 12.649 7 12.652 2 3 12.593 12.598 6 12.603 6 12.603 512.614 212.593 212.593 212.601 512.555 5 12.601 7 12.617 6 6 12.466 8 12.472 3 12.507 9 12.477 412.518 612.467 412.467 412.499 712.297 7 12.500 3 12.546 5 9 12.412 12.417 6 12.492 5 12.422 612.503 312.413 212.413 312.476 512.056 5 12.477 7 12.562 9 12 12.382 5 12.388 1 12.463 2 12.393 212.474 12.384 212.384 312.470611.920 6 12.472 4 12.557 6 18 12.333 2 12.338 9 12.414 3 12.343 912.425 112.336 612.336 612.513 711.295 2 12.517 12.601 9 22 12.626 7 12.627 5 12.628 3 12.628 312.63 12.626 712.626 712.628 12.620 8 12.628 12.630 5 25 12.536 2 12.564 3 12.546 8 12.577 412.589 112.536 412.536 412.544 712.498 5 12.544 9 12.583 4 33 12.343 4 12.349 12.384 9 12.354 112.395 712.344 12.344 112.376 612.172 5 12.377 3 12.423 9 34 0.370 19 0.370 36 0.372 62 0.370 510.372 950.370 290.370 30.376 160.332 34 0.376 26 0.378 81 35 12.536 2 12.564 3 12.546 8 12.577 412.589 112.536 412.536 412.544 712.498 5 12.544 9 12.583 4 36 12.536 2 12.564 3 12.546 8 12.577 412.589 112.536 412.536 412.544 712.498 5 12.544 9 12.583 4 37 0.370 19 0.370 36 0.372 62 0.370 510.372 950.370 290.370 30.376 160.332 34 0.376 26 0.378 81 38 0.370 19 0.370 36 0.372 62 0.370 510.372 950.370 290.370 30.376 160.332 34 0.376 26 0.378 81 39 12.412 12.417 6 12.492 5 12.422 612.503 312.413 212.413 312.476 512.056 5 12.477 7 12.562 9 40 12.66 12.66 12.66 12.66 12.66 12.66 12.66 12.66 12.66 12.66 12.66 41 0.370 19 0.370 36 0.372 62 0.370 510.372 950.370 290.370 30.376 160.332 34 0.376 26 0.378 81 表 3 潮流计算结果(迭代次数) Tab. 3 Convergence of power flow 循环类别 方案 1 方案 2 方案 3 方案 4方案 5方案 6方案 7方案 8方案 9 方案 10 方案 11 内循环 5 5 5 7 7 5 5 5 10 5 5 外循环 2 2 2 2 2 2 2 2 2 2 2 要靠电网提供的无功建立磁场，所以方案9降低了 各节点的电压水平。其它分布式电源对电网电压都 有很好的支撑作用。 由表3知：方案211与方案1比并没有明显 地增加迭代次数，方案5、10和11(配电网中含多 个分布式电源)迭代次数也没有明显增加。 表明本文 提出的算法对于处理多种分布式电源混合并网是 有效的，迭代次数并未因接入分布式电源的增多而 有明显增加。 对比方案6和7可看出算法可很好地处理电源 的不对称，网络本身三相负荷和支路参数不对称， 每个方案都有环网存在，表明算法可很好地处理配 电网的弱环和三相不对称问题。 5 结论 在对风力发电机组、 光伏发电系统、 燃料电池、 蓄电池、高频微型燃气轮机及工频热电联产机组的 运行和控制性能进行详细分析的基础上，本文建立 了各种分布式电源潮流计算模型，考虑配电网三相 参数不平衡和环网问题，对广泛用于辐射型配电网 40 中 国 电 机 工 程 学 报 第 29 卷 的前推回代潮流算法进行改进，提出了适用于含多 种分布式电源的弱环配电网的前推回代三相潮流 算法。 通过对例系统的大量测试， 得到了如下结论： 1）本文算法可很好地处理配电网的弱环和三 相不对称的问题。 2）本文所述的分布式电源模型和潮流算法， 对于处理多种分布式电源混合并网是有效的，迭代 次数并未因接入分布式电源的增多而有明显增加。 3）工频热电联产异步发电机、定速型异步风 机和转差控制型异步风机有可能降低系统电压水 平，其它的分布式电源对电网电压都有很好的支撑 作用。 参考文献 1 颜伟， 刘方， 王官浩 三相辐射型配电网络的相分量潮流计算J 电 力系统自动化，2002，10(5)：24-27 Yan Wei，Liu Fang，Wang GuanhaoThree-phase power flow using phase component for radial distribution systemJAutomation of Electrical Power Systems，2002，10(5)：24-27(in Chinese) 2 王志群，朱守真，周双喜分布式发电对配电网电压分布的影响 J电力系统自动化，2004，28(16)：56-60 Wang Zhiqun， Zhu Shouzhen， Zhou Shuangxi Impacts of distributed generation on distribution system voltage profileJAutomation of Electrical Power Systems，2004，28(16)：56-60(in Chinese) 3 丁明，王敏分布式发电技术J电力自动化设备，2004，24(7)： 31-36 Ding Ming， Wang Min Distributed generation technologyJ Electrical Power Automation Equipment，2004，24(7)：31-36(in Chinese) 4 王成山，郑海峰，谢莹华，等计及分布式发电的配电系统随机 潮流计算J电力系统自动化，2005，29(24)：39-44 Wang Chengshan，Zheng Haifeng，Xie Yinghua，et alProbabilistic power flow containing distributed generation in distribution systemJAutomation of Electric Power Systems，2005，29(24)： 39-44(in Chinese) 5 王守相， 江兴月， 王成山 含风力发电机组的配电网潮流计算J 电 网技术，2006，30(21)：42-45 Wang Shouxiang，Jiang Xingyue，Wang ChengshanPower flow analysis of distribution network containing wind power generatorsJPower System Technology，2006，30(21)：42-45(in Chinese) 6 Feijo A E，Cidrs JModeling of wind farms in the load flow analysisJIEEE Transactions on Power Systems，2000，15(1)： 110-115 7 陈海焱，陈金富， 段献忠含分布式电源的配电网潮流计算J电 力系统自动化，2006，30(1)：35-40 Chen Haiyan，Chen Jinfu，Duan XianzhongStudy on power flow calculation of distribution system with DGsJ Automation of Electric Power Systems，2006，30(1)：35-40(in Chinese) 8 王守相，黄丽娟，王成山分布式发电系统的不平衡三相潮流计 算J电力自动化设备，2007，27(8)：11-15 Wang Shouxiang，Huang Lijuan，Wang ChengshanUnbalanced three-phase power flow calculation for distributed power generation systemJElectrical Power Automation Equipment，2007，27(8)： 11-15(in Chinese) 9 颜伟，刘方，王官浩，等辐射型网络潮流的分层前推回代算法 J中国电机工程学报，2003，23(8)：76-80 Yan Wei，Liu Fang，Wang Guanhao，et alLayer-by-layer back/forward sweep method for radial distribution load flowJProceedings of the CSEE，2003，23(8)：76-80(in Chinese) 10 Jung J WModeling and control of fuel cell based distributed generation systemsD Columbia： The Ohio State University， 2005 11 El-Sharkh M Y， Rahma