Optimal Power Flow Solution Using Particle Swarm Optimization Technique with Global-Local Best P

合集下载

基于粒子群优化算法和电导增量法的多峰值mppt控制

目前袁传统化石能源造成的环境污染日趋严重袁太阳能由于其普遍尧清洁尧巨大尧长久的特点袁得到了广泛利用[1]遥然而袁光伏发电也存在着受环境影响大尧转换效率低等缺点遥因此袁如何维持复杂环境条件下的光伏最大功率输出成为一项重要研究内容[2]遥传统最大功率点跟踪 MPP境光照均匀一致条
第6期
杨海柱袁等院基于粒子群优化算法和电导增量法的多峰值 MPPT 控制
129
杨海柱袁岳刚伟袁康乐
渊河南理工大学电气工程与自动化学院袁焦作 454000冤
摘要院复杂环境条件下袁光伏阵列由于被遮挡其输出特性呈现多峰值特性袁传统最大功率点跟踪 MPPT 渊maximum power point tracking冤算法不再适用遥为此袁在研究光伏阵列多峰值输出特性的基础上袁提出一种基于粒子群优化 PSO渊particle swarm optimization冤算法和电导增量法 INC渊incremental conductance冤的多峰值 MPPT 算法遥该算法分成 2 步院第 1 步先由 PSO 算法将输入位置调整到最优值附近曰第 2 步再由 INC 算法得到全局最优解袁其中对传统 PSO 算法进行改进袁INC 算法采用变步长扰动遥在 Matlab 中进行仿真袁结果表明该算法可实现复杂环境条件下的最大功率跟踪袁并具备较快的响应速度和稳定的寻优效果遥
粤遭泽贼则葬糟贼院 Under complex environmental conditions, the traditional maximum power point tracking渊MPPT冤 algorithm is no longer applicable due to the multi鄄peak characteristics of output from a PV array caused by shading. To solve this problem, a multi鄄peak MPPT algorithm based on particle swarm optimization渊PSO冤 and incremental conductance渊INC冤 al鄄 gorithms is proposed. This algorithm is divided into two steps: the first step is to adjust the input position to near the op鄄 timal value by using the PSO algorithm, and the second step is to obtain a global optimal solution by using the INC al鄄 gorithm. Specifically, the traditional PSO algorithm is improved, and the INC algorithm uses a variable鄄step size. Simula鄄 tion results in Matlab show that the proposed algorithm can achieve the maximum power tracking under complex envi鄄 ronmental conditions; in addition, it has a faster response speed and stable optimization results.

原-对偶内点法和预测-校正内点法在最优潮流的应用

原-对偶内点法和预测-校正内点法在最优潮流的应用杨利水;杨旭;顾家翠【摘要】Optimal power flow is a nonlinear optimization problem subjected to constraints. The model of OPF includes the objective function, equality constraints and inequality constraints. The paper put Primal-dual IPM and Predictor-corrector IPM in the calculation of Optimal power flow. Primal-dual IPM keeps the primal feasibility and the dual feasibility when searching the optimal solution along the primal-dual path. Predictor-corrector IPM keeps the higher order term when doing the Taylor expansion. Therefore, PCIPM has a better convergence than PDIPM. The Primal-dual IPM and the Predictor-corrector IPM are introduced and the codes written in matlab language are given and examples are presented.%最优潮流问题在数学上是一个带约束条件的优化问题,其模型包括目标函数以及等式约束条件和不等式约束条件.利用原-对偶内点法和预测-校正内点法进行最优潮流的计算,原-对偶内点法是在保持原始可行性和对偶可行性的同时,沿一条原-对偶路径寻找最优解.预测-校正法在进行泰勒展开时保留了高阶项,首先通过修正方程计算仿射方向,在计算得到仿射扰动因子后回代入修正方程得到校正方向,进而得到修正量.预测-校正法具有比原-对偶法更好的收敛性,用Matlab实现了原-对偶内点法和预测-校正内点法进行潮流优化计算,并用算例进行了验证.【期刊名称】《华北电力大学学报(自然科学版)》【年(卷),期】2012(039)006【总页数】6页(P29-34)【关键词】最优潮流;非线性规划;内点法;原-对偶法;预测-校正法【作者】杨利水;杨旭;顾家翠【作者单位】保定电力职业技术学院,河北保定071051;深圳供电局,广东深圳518000;广东电网公司教育培训评价中心,广东广州510520【正文语种】中文【中图分类】TM7440 引言电力系统最优潮流OPF(Optimal Power Flow)是一个复杂的非线性规划问题，要求在满足特定的电力系统运行和安全约束条件下，通过各种控制手段实现预定目标最优的系统稳定运行状态。

matlab潮流计算工具箱使用手册

MATPOWERA M ATLAB™ Power System Simulation PackageVersion 3.2September 21, 2007User’s ManualRay D. Zimmerman Carlos E. Murillo-Sánchez rz10@ carlos_murillo@ © 1997-2007 Power Systems Engineering Research Center (PS ERC)School of Electrical Engineering, Cornell University, Ithaca, NY 14853Table of ContentsTable of Contents (2)1Introduction (3)2Getting Started (4)2.1System Requirements (4)2.2Installation (4)2.3Running a Power Flow (4)2.4Running an Optimal Power Flow (4)2.5Getting Help (4)3Technical Reference (6)3.1Data File Format (6)3.2Modeling (8)3.3Power Flow (11)3.4Optimal Power Flow (12)3.4.1AC OPF Formulation (13)3.4.2DC OPF Formulation (21)3.5Unit Decommitment Algorithm (22)3.6MATPOWER Options (22)3.7Summary of the Files (28)4Acknowledgments (33)5References (33)Appendix A: Notes on LP-Solvers for M ATLAB (34)Appendix B: Additional Notes (34)Appendix C: Auction Code (35)1IntroductionWhat is MATPOWER?MATPOWER is a package of M ATLAB M-files for solving power flow and optimal power flow prob-lems. It is intended as a simulation tool for researchers and educators that is easy to use and modify. MATPOWER is designed to give the best performance possible while keeping the code simple to understand and modify. The MATPOWER home page can be found at:/matpower/Where did it come from?MATPOWER was developed by Ray D. Zimmerman, Carlos E. Murillo-Sánchez and Deqiang Gan of PSERC at Cornell University (/) under the direction of Robert Thomas. The initial need for M ATLAB based power flow and optimal power flow code was born out of the computational requirements of the PowerWeb project (see /powerweb/).Who can use it?•MATPOWER is free. Anyone may use it.•We make no warranties, express or implied. Specifically, we make no guarantees regarding the correctness MATPOWER’s code or its fitness for any particular purpose.•Any publications derived from the use of MATPOWER must cite MATPOWER /matpower/.•Anyone may modify MATPOWER for their own use as long as the original copyright notices remain in place.•MATPOWER may not be redistributed without written permission.•Modified versions of MATPOWER, or works derived from MATPOWER, may not be distrib-uted without written permission.2Getting Started2.1System RequirementsTo use MATPOWER you will need:•M ATLAB version 6 or later1•M ATLAB Optimization Toolbox (required only for some OPF algorithms)Both are available from The MathWorks (see /).2.2InstallationStep 1: Go to the MATPOWER home page (/matpower/) and follow the download instructions.Step 2: Unzip the downloaded file.Step 3: Place the files in a location in your M ATLAB path.2.3Running a Power FlowTo run a simple Newton power flow on the 9-bus system specified in the file case9.m, with the de-fault algorithm options, at the M ATLAB prompt, type:>> runpf('case9')2.4Running an Optimal Power FlowTo run an optimal power flow on the 30-bus system whose data is in case30.m, with the default algo-rithm options, at the M ATLAB prompt, type:>> runopf('case30')To run an optimal power flow on the same system, but with the option for MATPOWER to shut down (decommit) expensive generators, type:>> runuopf('case30')2.5Getting HelpAs with M ATLAB’s built-in functions and toolbox routines, you can type help followed by the name of a command or M-file to get help on that particular function. Nearly all of MATPOWER’s M-files have such documentation. For example, the help for runopf looks like:1Although it is likely that most things work fine in M ATLAB5, this is not supported due to limited testing resources. MATPOWER 3.0 required M ATLAB 5 and MATPOWER 2.0 and earlier required only M ATLAB 4.>> help runopfRUNOPF Runs an optimal power flow.[baseMVA, bus, gen, gencost, branch, f, success, et] = ...runopf(casename, mpopt, fname, solvedcase)Runs an optimal power flow and optionally returns the solved values in the data matrices, the objective function value, a flag which is true if the algorithm was successful in finding a solution, and the elapsed time in seconds. All input arguments are optional. If casename is provided it specifies the name of the input data file or struct (see also 'helpcaseformat' and 'help loadcase') containing the opf data. The defaultvalue is 'case9'. If the mpopt is provided it overrides the defaultMATPOWER options vector and can be used to specify the solutionalgorithm and output options among other things (see 'help mpoption' for details). If the 3rd argument is given the pretty printed output will be appended to the file whose name is given in fname. If solvedcase isspecified the solved case will be written to a case file in MATPOWERformat with the specified name. If solvedcase ends with '.mat' it saves the case as a MAT-file otherwise it saves it as an M-file. MATPOWER also has many options which control the algorithms and the output. Type:>> help mpoptionand see Section 3.6 for more information on MATPOWER's options.3Technical Reference3.1Data File FormatThe data files used by MATPOWER are simply M ATLAB M-files or MAT-files which define and re-turn the variables baseMVA, bus, branch, gen, areas, and gencost. The baseMVA variable is a scalar and the rest are matrices. Each row in the matrix corresponds to a single bus, branch, or generator. The columns are similar to the columns in the standard IEEE and PTI formats. The details of the specifica-tion of the MATPOWER case file can be found in the help for caseformat.m:>> help caseformatCASEFORMAT Defines the MATPOWER case file format.A MATPOWER case file is an M-file or MAT-file which defines the variablesbaseMVA, bus, gen, branch, areas, and gencost. With the exception ofbaseMVA, a scalar, each data variable is a matrix, where a row correspondsto a single bus, branch, gen, etc. The format of the data is similar tothe PTI format described in/research/pstca/formats/pti.txtexcept where noted. An item marked with (+) indicates that it is includedin this data but is not part of the PTI format. An item marked with (-) isone that is in the PTI format but is not included here. Those marked with(2) were added for version 2 of the case file format. The columns foreach data matrix are given below.MATPOWER Case Version Information:A version 1 case file defined the data matrices directly. The last two,areas and gencost, were optional since they were not needed for runninga simple power flow. In version 2, each of the data matrices are storedas fields in a struct. It is this struct, rather than the individualmatrices that is returned by a version 2 M-casefile. Likewise a version 2MAT-casefile stores a struct named 'mpc' (for MATPOWER case). The structalso contains a 'version' field so MATPOWER knows how to interpret thedata. Any case file which does not return a struct, or any struct whichdoes not have a 'version' field is considered to be in version 1 format.See also IDX_BUS, IDX_BRCH, IDX_GEN, IDX_AREA and IDX_COST regardingconstants which can be used as named column indices for the data matrices.Also described in the first three are additional columns that are addedto the bus, branch and gen matrices by the power flow and OPF solvers.Bus Data Format1 bus number (1 to 29997)2 bus typePQ bus = 1PV bus = 2reference bus = 3isolated bus = 43 Pd, real power demand (MW)4 Qd, reactive power demand (MVAr)5 Gs, shunt conductance (MW (demanded) at V = 1.0 p.u.)6 Bs, shunt susceptance (MVAr (injected) at V = 1.0 p.u.)7 area number, 1-1008 Vm, voltage magnitude (p.u.)9 Va, voltage angle (degrees)(-) (bus name)10 baseKV, base voltage (kV)11 zone, loss zone (1-999)(+) 12 maxVm, maximum voltage magnitude (p.u.)(+) 13 minVm, minimum voltage magnitude (p.u.)Generator Data Format1 bus number(-) (machine identifier, 0-9, A-Z)2 Pg, real power output (MW)3 Qg, reactive power output (MVAr)4 Qmax, maximum reactive power output (MVAr)5 Qmin, minimum reactive power output (MVAr)6 Vg, voltage magnitude setpoint (p.u.)(-) (remote controlled bus index)7 mBase, total MVA base of this machine, defaults to baseMVA(-) (machine impedance, p.u. on mBase)(-) (step up transformer impedance, p.u. on mBase)(-) (step up transformer off nominal turns ratio)8 status, > 0 - machine in service<= 0 - machine out of service(-) (% of total VAr's to come from this gen in order to hold V at remote bus controlled by several generators)9 Pmax, maximum real power output (MW)10 Pmin, minimum real power output (MW)(2) 11 Pc1, lower real power output of PQ capability curve (MW)(2) 12 Pc2, upper real power output of PQ capability curve (MW)(2) 13 Qc1min, minimum reactive power output at Pc1 (MVAr)(2) 14 Qc1max, maximum reactive power output at Pc1 (MVAr)(2) 15 Qc2min, minimum reactive power output at Pc2 (MVAr)(2) 16 Qc2max, maximum reactive power output at Pc2 (MVAr)(2) 17 ramp rate for load following/AGC (MW/min)(2) 18 ramp rate for 10 minute reserves (MW)(2) 19 ramp rate for 30 minute reserves (MW)(2) 20 ramp rate for reactive power (2 sec timescale) (MVAr/min)(2) 21 APF, area participation factorBranch Data Format1 f, from bus number2 t, to bus number(-) (circuit identifier)3 r, resistance (p.u.)4 x, reactance (p.u.)5 b, total line charging susceptance (p.u.)6 rateA, MVA rating A (long term rating)7 rateB, MVA rating B (short term rating)8 rateC, MVA rating C (emergency rating)9 ratio, transformer off nominal turns ratio ( = 0 for lines )(taps at 'from' bus, impedance at 'to' bus, i.e. ratio = Vf / Vt)10 angle, transformer phase shift angle (degrees), positive => delay(-) (Gf, shunt conductance at from bus p.u.)(-) (Bf, shunt susceptance at from bus p.u.)(-) (Gt, shunt conductance at to bus p.u.)(-) (Bt, shunt susceptance at to bus p.u.)11 initial branch status, 1 - in service, 0 - out of service(2) 12 minimum angle difference, angle(Vf) - angle(Vt) (degrees)(2) 13 maximum angle difference, angle(Vf) - angle(Vt) (degrees)(+) Area Data Format1 i, area number2 price_ref_bus, reference bus for that area(+) Generator Cost Data FormatNOTE: If gen has n rows, then the first n rows of gencost containthe cost for active power produced by the corresponding generators.If gencost has 2*n rows then rows n+1 to 2*n contain the reactivepower costs in the same format.1 model, 1 - piecewise linear,2 - polynomial2 startup, startup cost in US dollars3 shutdown, shutdown cost in US dollars4 n, number of cost coefficients to follow for polynomialcost function, or number of data points for piecewise linear5 and following, cost data defining total cost functionFor polynomial cost:c2, c1, c0where the polynomial is c0 + c1*P + c2*P^2For piecewise linear cost:x0, y0, x1, y1, x2, y2, ...where x0 < x1 < x2 < ... and the points (x0,y0), (x1,y1),(x2,y2), ... are the end- and break-points of the cost function. Some columns are added to the bus, branch and gen matrices by the solvers. See the help for idx_bus, idx_brch, and idx_gen for more details.3.2ModelingAC FormulationFixed loads are modeled as constant real and reactive power injections,!Pdand!Qdspecified in col-umns 3 and 4, respectively, of the bus matrix. The shunt admittance of any constant impedance shuntelements at a bus are specified by Gsh and Bshin columns 5 and 6, respectively, of the bus matrix!Ysh=Gsh+jBshbaseMVAEach branch, whether transmission line, transformer or phase shifter, is modeled as a standard π-model transmission line, with series resistance R and reactance X and total line charging capacitance!Bc, in series with an ideal transformer and phase shifter, at the from end, with tap ratio!" and phase shift angle!"shift. The parameters R, X,!Bc,!" and!"shift, are found in columns 3, 4, 5, 9 and 10 of the branch matrix, respectively. The branch voltages and currents at the from and to ends of the branch are related by the branch admittance matrix!Ybras follows!IfIt"#$%&'=Y brVfVt"#$%&'(1) where!Ybr=Ys+jBc2"#$%&'1(2)Ys1(e j*shift)Ys1(e)j*shiftYs+jBc2+,---./and!Ys=1R+jX.The elements of the individual branch admittance matrices and the bus shunt admittances are com-bined by MATPOWER to form a complex bus admittance matrix!Ybus, relating the vector of complex bus voltages!Vbuswith the vector of complex bus current injections!Ibus!Ibus=YbusVbusSimilarly, admittance matrices!Yfand!Yt, are formed to compute the vector of complex current injec-tions at the from and to ends of each line, given the bus voltages!Vbus.!If=YfVbusIt=YtVbusThe vectors of complex bus power injections, and branch power injections can be expressed as!Sbus=diag(Vbus)Ibus*Sf=diag(Vf)If*St=diag(Vt)It*where!Vfand!Vtare vectors of the complex bus voltages at the from and to ends, respectively, of all branches, and diag() converts a vector into a diagonal matrix with the specified vector on the diagonal.DC FormulationFor the DC formulation, the same parameters are used, with the exception that the following assump-tions are made:•Branch resistances R and charging capacitances!Bcare negligible (i.e. branches are lossless).•All bus voltage magnitudes are close to 1 p.u.•Voltage angle differences are small enough that!sin"ij#"ij.Combining these assumptions and equation (1) with the fact that!S=VI*, the relationship between thereal power flows and voltage angles for an individual branch can be written as!PfPt"#$%&'=B br(f(t"#$%&'+Pf,shiftPt,shift"#$%&'(2)where!Bbr=1X"1#1#11$%&'()(3)!Pf,shiftPt,shift"#$%&'=(shiftX)1*1"#$%&'.(4)The elements of the individual branch shift injections and!Bbrmatrices are combined by MATPOWERto form a bus!Bbusmatrix and!Pbus,shiftshift injection vector, which can be used to compute bus realpower injections from bus voltage angles!Pbus=Bbus"bus+Pbus,shiftSimilarly, MATPOWER builds the matrix!Bfand the vector!Pf,shiftwhich can be used to compute thevectors !Pfand!Ptof branch real power injections!Pf=Bf"bus+Pf,shiftPt=#Pf3.3Power FlowMATPOWER has five power flow solvers, which can be accessed via the runpf function. In addition to printing output to the screen, which it does by default, runpf optionally returns the solution in out-put arguments:>> [baseMVA, bus, gen, branch, success, et] = runpf(casename);The solution values are stored as follows:bus(:, VM)bus voltage magnitudesbus(:, VA)bus voltage anglesgen(:, PG)generator real power injectionsgen(:, QG)generator reactive power injectionsbranch(:, PF)real power injected into “from” end of branchbranch(:, PT)real power injected into “to” end of branchbranch(:, QF)reactive power injected into “from” end of branchbranch(:, QT)reactive power injected into “to” end of branchsuccess 1 = solved successfully, 0 = unable to solveet computation time required for solutionThe default power flow solver is based on a standard Newton’s method [10] using a full Jacobian, up-dated at each iteration. This method is described in detail in many textbooks. Algorithms 2 and 3 are variations of the fast-decoupled method [9]. MATPOWER implements the XB and BX variations as described in [1]. Algorithm 4 is the standard Gauss-Seidel method from Glimm and Stagg [3], based on code contributed by Alberto Borghetti, from the University of Bologna, Italy. To use one of the power flow solvers other than the default Newton method, the PF_ALG option must be set explicitly. For example, for the XB fast-decoupled method:>> mpopt = mpoption('PF_ALG', 2);>> runpf(casename, mpopt);The last method is a DC power flow [11], which is obtained by executing runpf with the PF_DC op-tion set to 1, or equivalently by executing rundcpf directly. The DC power flow is obtained by a di-rect, non-iterative solution of the bus voltage angles from the specified bus real power injections, based on equations (2), (3) and (4).For the AC power flow solvers, if the ENFORCE_Q_LIMS option is set to 1 (default is 0), then if any generator reactive power limit is violated after running the AC power flow, the corresponding bus is converted to a PQ bus, with the reactive output set to the limit, and the case is re-run. The voltage magnitude at the bus will deviate from the specified value in order to satisfy the reactive power limit. If the generator at the reference bus reaches a reactive power limit and the bus is converted to a PQ bus, the first remaining PV bus will be used as the slack bus for the next iteration. This may result in the real power output at this generator being slightly off from the specified values.Currently, none of MATPOWER’s power flow solvers include any transformer tap changing or han-dling of disconnected or de-energized sections of the network.Performance of the power flow solvers, with the exception of Gauss-Seidel, should be excellent even on very large-scale power systems, since the algorithms and implementation take advantage of M ATLAB’s built-in sparse matrix handling.3.4Optimal Power FlowMATPOWER includes several solvers for the optimal power flow (OPF) problem, which can be ac-cessed via the runopf function. In addition to printing output to the screen, which it does by default, runopf optionally returns the solution in output arguments:>> [baseMVA, bus, gen, gencost, branch, f, success, et] = runopf(casename);In addition to the values listed for the power flow solvers, the OPF solution also includes the follow-ing values:bus(:, LAM_P)Lagrange multiplier on bus real power mismatchbus(:, LAM_Q)Lagrange multiplier on bus reactive power mismatchbus(:, MU_VMAX)Kuhn-Tucker multiplier on upper bus voltage limitbus(:, MU_VMIN)Kuhn-Tucker multiplier on lower bus voltage limitgen(:, MU_PMAX)Kuhn-Tucker multiplier on upper generator real power limitgen(:, MU_PMIN)Kuhn-Tucker multiplier on lower generator real power limitgen(:, MU_QMAX)Kuhn-Tucker multiplier on upper generator reactive power limitgen(:, MU_QMIN)Kuhn-Tucker multiplier on lower generator reactive power limitbranch(:, MU_SF)Kuhn-Tucker multiplier on MVA limit at "from" end of branchbranch(:, MU_ST)Kuhn-Tucker multiplier on MVA limit at "to" end of branchf final objective function valueMATPOWER can make use of a number of different OPF solvers. There are two legacy solvers from early versions of MATPOWER, namely the constr and LP-based solvers, that have been deprecated and will be removed from future versions. The details of the problem formulation and solution algo-rithms used by these solvers can be found in the user's manual included with previous versions of MATPOWER.The current generation of solvers use the generalized AC OPF formulation described below. MATPOWER includes one based on fmincon from M ATLAB’s Optimization Toolbox and there are two optional packages, MINOPF2 and TSPOPF3, that implement higher performance OPF solvers us-ing MEX files. MINOPF, based on the MINOS [7] solver, has been available since mid-2004 and is distributed separately because it has a more restrictive license than MATPOWER. TSPOPF is a collec-tion of three solvers developed by Hongye Wang [11] and is currently distributed separately as well. The performance of MATPOWER’s OPF solvers depends on several factors. First, for problems of this general nature, fmincon does not exploit and preserve sparsity, so it is inherently limited to solving small power systems. The MEX based solvers, on the other hand, do exploit sparsity and are suitable for much larger problems. MINOPF is coded in FORTRAN and evaluates the required Jacobians us-ing an optimized structure that follows the order of evaluation imposed by the compressed-column sparse format which is employed by MINOS. In fact, the new generalized OPF formulation included in MATPOWER 3.0 and later is inspired by the data format used by MINOS. The solvers in the TSPOPF package are implemented in the C language.MATPOWER’s OPF implementation is not currently able to handle unconnected or de-energized sec-tions of the network.2 See /minopf/.3 See /tspopf/.Piecewise linear costs using constrained cost variables (CCV)The OPF formulations in MATPOWER allow for the specification of convex piecewise linear costfunctions for active or reactive generator output. An example of such a cost curve is shown below.This non-differentiable cost is modeled using an extra helper cost variable for each such cost curve and additional constraints on this variable and Pg , one for each segment of the curve. The constraints build a convex “basin” equivalent to requiring the cost variable to lie in the epigraph of the cost curve. When the cost is minimized, the cost variable will be pushed against this basin. If the helper cost vari-able is y , then the contribution of the generator’s cost to the total cost is exactly y . In the above case, the two additional required constraints are1) !y "m 1(P g #x 0)+c 0(y must lie above the first segment) 2) ! y "m 2(P g #x 1)+c 1 (y must lie above the second segment)where m 1 and m 2 are the slopes of the two segments. Also needed, of course, are the box restrictions on P g : P min ≤ P g ≤ P max . The additive part of the cost contributed by this generator is y .This constrained cost variable (CCV) formulation is used by all of the MATPOWER OPF solvers for handling piecewise linear cost functions, with the exception of two that are part of the optional TSPOPF package, namely the step-controlled primal/dual interior point method (SCPDIPM ) and the trust region based augmented Lagrangian method (TRALM ), both of which use a cost smoothing tech-nique instead.3.4.1 AC OPF FormulationThe AC optimal power flow problem solved by MATPOWER is a “smooth” OPF with no discrete variables or controls. The generalized AC OPF formulation, used by the current generation of MATPOWER ’s OPF solvers, offers a number of extra capabilities relative to the traditional formula-tion of minimizing the cost of generation subject to voltage, flow and generator limits, used by the first generation of MATPOWER OPF solvers:• mixed polynomial and piecewise linear costs• dispatchable loads• generator P-Q capability curves• branch angle difference limits • additional user supplied linear constraints• additional user supplied costsNew in MATPOWER 3.2 are the generalized user supplied cost formulation, the generator capability curves, the branch angle difference limits and a simplification of the general linear constraint specifi-cation used in version 3.0.The problem is formulated in terms of two groups of optimization variables, labeled x and z . The x variables are the OPF variables, consisting of the voltage angles ! " and magnitudes V at each bus, and real and reactive generator injections P g and Q g .! x ="V P g Q g # $ % % % % & '( ( ( ( Additional user defined variables are grouped in z .The optimization problem can be expressed as follows: ! min x ,y ,z f 1i (P gi )+f 2i (Q gi )()i"+12w T Hw +C w T wsubject to!g P (x )=P (",V )#P g +P d =0 (active power balance equations)!g Q (x )=Q (",V )#Q g +Q d =0 (reactive power balance equations) ! g S f (x )=S f (",V )#S max $0 (apparent power flow limit of lines, from end) ! g S t (x )=S t (",V )#S max $0 (apparent power flow limit of lines, to end) ! l "A x z # $ % & ' ( "u (general linear constraints)! x min "x "x max (voltage and generation variable limits)! z min "z "z max (limits on user defined variables) Here f 1i and f 2i are the costs of active and reactive power generation, respectively, for generator i at a given dispatch point. Both f 1i and f 2i are assumed to be polynomial or piecewise-linear functions.The most significant additions to the traditional, simple OPF formulation appear in the generalizedcost terms containing w and in the general linear constraints involving the matrix A, described in the next two sections. These two frameworks allow tremendous flexibility in customizing the problem formulation, making MATPOWER even more useful as a research tool.Note: In Optimization Toolbox versions 3.0 and earlier, fmincon seems to be providing inaccurate shadow prices on the constraints. This did not happen with constr and it may be a bug in these ver-sions of the Optimization Toolbox.General Linear ConstraintsIn addition to the standard non-linear equality constraints for nodal power balance and non-linear ine-quality constraints for line flow limits, this formulation includes a framework for additional linear constraints involving the full set of optimization variables.! l "A x z # $ % & '( "u (general linear constraints) Some portions of these linear constraints are supplied directly by the user, while others are generated automatically based on the case data. Automatically generated portions include:• rows for constraints that define generator P-Q capability curves• rows for constant power factor constraints for dispatchable or price-sensitive loads• rows for branch angle difference limits• rows and columns for the helper variables from the CCV implementation of piecewise linear gen-erator costs and their associated constraintsIn addition to these automatically generated constraints, the user can provide a matrix A u and vectors l u and u u to define further linear constraints. These user supplied constraints could be used, for example, to restrict voltage angle differences between specific buses. The matrix A u must have at least n x col-umns where n x is the number of x variables. If A u has more than n x columns, a corresponding z optimi-zation variable is created for each additional column. These z variables also enter into the generalized cost terms described below, so A u and N must have the same number of columns.!l u "A u x z # $ % & '( "u u (user supplied linear constraints) Change from MATPOWER 3.0: The A u matrix supplied by the user no longer includes the (all zero) columns corresponding to the helper variables for piecewise linear generator costs. This should sim-plify significantly the creation of the desired A u matrix. Generalized Cost FunctionThe cost function consists of two parts. The first is the polynomial or piecewise linear cost of genera-tion. A polynomial or piecewise linear cost is specified for each generator’s active output and, option-ally, reactive output in the appropriate row(s) of the gencost matrix. Any piecewise linear costs are implemented using the CCV formulation described above which introduces corresponding helper variables. The general formulation allows generator costs of mixed type (polynomial and piecewise linear) in the same problem.The second part of the cost function provides a general framework for imposing additional costs on the optimization variables, enabling things such as using penalty functions as soft limits on voltages, additional costs on variables involved in constraints handled by Langrangian relaxation, etc.。

法国电网电能质量承诺和电能质量评估_英文_

2009年第3卷第3期南方电网技术特约专稿2009，V ol. 3，No. 3 SOUTHERN POWER SYSTEM TECHNOLOGY FeaturedArticles 文章编号：1674-0629(2009)03-0007-08 中图分类号：TM934.5 文献标志码：A 法国电网电能质量承诺和电能质量评估杨显军（法国电力公司研究与开发部，92141 卡拉马特，法国）摘要：介绍法国电网的电能质量承诺，包括输电运营商（TSO）的质量承诺和配电运营商（DSO）的质量承诺。

TSO和DSO在作出电能质量承诺的同时，还利用入网规则以协议的方式来考虑和限制每个客户的负荷和分布式电源对电能质量的影响。

这些协议的实施表明，只有遵守承诺，加上电网运营商与用户双方共同作出努力，才能获得全面良好的电能质量。

对于电网干扰评估，建议使用基波潮流计算与多相谐波注入相结合的方法。

这个方法利用很短的计算时间就能对非线性负荷建模情形给出可接受的结果。

所提议的频域方法的要点在于，将每个谐波源逐一地同时作为扰动源和受扰者来考察，以便计及扰动源之间的相互作用。

研发了一个即插即用的HVDC组件并将其用在文章最后部分实例研究中。

关键词：电能质量承诺；频域建模；谐波；电力电子设备；高压直流输电Grid Power Quality Commitments and Grid Power Quality Assessmentin FranceYANG Xian-jun (Xavier X. YANG)(R&D, Electricité de France, 92141 Clamart, France)Abstract: Grid Power Quality Commitments in France are introduced, including the Power Quality Commitments of TransmissionSystem Operator (TSO) and the Power Quality Commitments of Distribution System Operator (DSO). Undertaking the commitments,TSO and DSO are using connection rules to take into account and minimize the impact on power quality of each customer (load ordistributed generation) in the way of contracts. The practice of these contracts shows that a good overall power quality can only beachieved by respecting commitments and doing effort from both of grid and customer sides. For grid disturbance assessment, it issuggested to use multi-phase harmonic injection method with fundamental power flow. This method can give an acceptable resultswith very short computing time for non-linear load modelling. The key point of proposed frequency domain method is that eachharmonic source is individually considered as both disturbance source and disturbance receiver in order to take into accountinteractions between sources. A plug & play HVDC component has been developed and used in one of two case-studies in the end ofthis article.Key words: power quality commitment; frequency domain modelling; harmonic; power electronic device; HVDC1 Power Quality Commitments1.1 Power Quality Contracts and CommitmentsStandard EN 50160 gives the main voltage pa-rameters and their permissible deviation ranges at the point of common coupling (PCC) in public low volt-age (LV) and medium voltage (MV) electricity distri-bution systems, under normal operating conditions.Based on EN 50160 and the former French power quality contract Emeraude, two new grid access con-tracts have been put into service in France between system operators (utilities) and customers: called CART for transmission customers and CARD for dis-tribution customers, see Fig. 1. These contracts not only set the power quality commitments for the grid operators but also for the customer disturbance emis-sion level.Fig. 1 Power Quality Commitments in CARD and CART8南方电网技术第3卷The power quality that is handled in grid access contracts includes harmonic, unbalance, flicker, inter-ruption, voltage dip disturbances. Among these dis-turbances, the number of interruption is considered as commitment and others are given as indicative values. 1.2 Power Quality Commitments of Transmis-sion System Operator (TSO)French TSO has the standard and personalized contractual commitments with the customers concern-ing:1） long interruptions (1 per year, 2 per 3 years, 1 per 3 years)；2） short interruptions (1 to 5 per year, 1 to 2 per 3 years)；3） voltage dips (optional from 3 to 5 per year)； 4） RMS voltage magnitude (± 8 %)； 5） frequency (± 1 %)； 6） voltage unbalance (< 2 %)； 7） harmonics ；8） rapid voltage fluctuations, flicker (P lt < 1). For transmission power supply voltage, harmonic limits have been set (Tab. 1). These limits are meas-ured by 10 min mean value at each harmonic voltage (h ≤ 25, the order of harmonics) and total harmonic distortion (THD , h ≤ 40) during 100% of time. For transmission customer, current harmonic emissions are also limited in % compared to contracted rated currenth HRU lim /%h HRU lim /%h HRU lim /%5, 7 2.0 3 2.0 2 1.5 6.0 11, 13 1.5 9 1.0 4 1.0 6.0 171.015, 210.5> 40.56.019 1.0 6.0 23, 25 0.76.01.3 Power Quality Commitments of DistributionSystem Operator (DSO)The Tab. 3 shows the contracted voltage interrup-tion values for the four geographical zones.With personalized contracts signed jointly by system operators and customers, it is possible to in-clude some voltage dips in guaranteed power supply quality by network operators.Tab. 2 Transmission Customer Current Harmonic Commitments Odd Order hHRI lim /% Even Oder h HRI lim /%3 4.0 2 2 5, 7 5.04 1.0 9 2.0 > 4 0.5 11, 133.0>13 2.0Tab. 3 DSO Commitments of Voltage Interruptions per Year Interruption DurationZone 1 Zone 2 Zone 3Zone 4≥ 3 min6 3 3 2 > 1 s and < 3 min301032As grid commitments (Tab. 4), harmonic limits have been set for distribution power supply voltage. These limits are measured by 10 min mean value ateach harmonic voltage (h ≤ 25) and total harmonic distortion (THD , h ≤ 40) during 100% of time. For distribution customer, current harmonic emissions are also limited in % compared to contracted rated current (Tab. 5).Tab. 4 Distribution Voltage Harmonic Commitments Odd Harmonic OrderNot-Triplet Triplet Even HarmonicOrderTHD lim /13 3.0 8.0 17 2.0 8.0 19, 23, 251.5 8.0Tab. 5 Distribution Customer Current Harmonic Commitments Odd Order hHRI lim /% Even Oder h HRI lim /%3 4.0 2 2.0 5, 7 5.04 1.0 9 2.0 > 4 0.5 11, 133.0>13 2.0第3期杨显军：法国电网电能质量承诺和电能质量评估 9For more detail of the power qualify commit-ments of TSO and DSO, please see the following web site: /.1.4 Customer Load ConnectionIn order to meet their power quality commitments, TSO and DSO are using connection rules to take into account and minimize the impact on power quality of each customer (load or distributed generation). DSO has set up a number of technical referential documents and studying tools for load connection. Before per-forming a load or a DG connection, DSO undertakes three level’s studies which are based on standards and site measurements. The aim is to anticipate the possi-ble disturbances which may be provoked by the equipment to be connected and to determine the most adaptive Point of Common Coupling (PCC). Follow-ing these steps, DSO can ensure the optimal perform-ance of the network for all users.The first level study includes identifying the main characteristics of the installation to be connected and effectuating a global analysis. The second level con-sists of evaluating particularly the disturbances. The third level is composed of a detailed study on the iden-tified disturbances. In wind mill connection, this study is based on the standard IEC 61400-21 which evalu-ates the impact in terms of power quality [1].1.5 Power Quality Indices and Evolution ofCommitmentsIn order to measure the evolution of power qual-ity, DSO calculates a number of criterions based on network parameters. These indicators make it possible to evaluate the needs of power supply infrastructure improvement. The following values have been defined as power quality indicators:1）SAIDI or B-Criterion in France, index of mean duration of interruption for LV customers;2）SAIDI for MV customers weighted by the supplied power;3）SAIFI, indices of mean number of interrup-tion of the network, for the long interruption (t id > 3 min), the short interruption (1 s≤ t id≤ 3 min) and the very short interruption (t id < 1 s) . 2 Harmonic Source Modelling and GridDisturbance Assessment2.1 Simulation MethodFrequency domain analysis uses an admittance matrix system model developed from individual com-ponent level models connected according to system topology. The development of admittance matrix mod-els is based on multi-port network theory. Linear grid components are similarly developed from multi-port admittance parameters.In frequency domain simulation, harmonic gener-ating devices are usually modelled by current sources and the whole system admittance matrix is resolved from a start frequency to an end frequency with a fixed or variable frequency step. The deficiencies in the current source method can be partially overcome using a technique that has come to be known as “harmonic power flow”. More often, a fundamental frequency power flow solution is executed using a linear model for all power delivery equipment and loads, and the resultant fundamental frequency load terminal volt-ages are used to “adjust” non-linear load harmonic current vectors automatically without additional user action. However, in frequency approach, initial har-monic current patterns are still required to be known for each nonlinear load.2.2 Harmonic Source ModellingHarmonic model based on a constant current source, see Fig. 2(a), can’t represent real harmonic currents generated by a power electronic device con-nected to an actual network.Fig. 2 General Structure of Three Phase Harmonic Source ModelThe method (a) is too simplified because actual harmonic current of an electric device is always con-10南方电网技术第3卷trolled by network configuration. In order to overcome this drawback, we have developed multiphase har-monic model with a voltage vector controlled current source and a frequency-variable impedance to repre-sent an actual harmonic source, see Fig. 2(b). In pro-posed method (b), J (f ) is a current source controlled by the fundamental voltage crossing frequency-variable impedance Z (f ). Z (f ) represents the behaviour of the power electronic device towards harmonic distur-bances coming from power supply side.The method (b) is detailed in the following para-graphs.Equivalent harmonic current source of a nonlin-ear load is generated by two main steps: the initial cur-rent J 0(f ) and the corrected current source J (f ) by load-flow. In parameter setting before simulation, the initial J 0(f ) is defined by pre-designed waveform pat-tern, rated power and rated voltage of the load. In fre-quency domain, J 0(f ) is represented by the followingformula:00()()cos()k kJ f J k k t ωφ=⋅+.The pre-designed waveform patterns include cur-rent waveforms generated by the most usual distur-bance sources such as 6-pulse and 12-pulse rectifiers, dimmers, arc furnaces, etc. The following procedure summarizes the proposed source modelling method:1） Choosing a pre-designed waveform (database, analytical formulas, wave form generated by Excel data sheet, empirical data table, on-site measurements) which represents the current harmonic distortion to be used;2） Refining the waveform by particular working state of the device. For example, by switching angle of rectifier and by short-circuit power of upstream net-work;3） Spectrum calculation of the refined wave-form by means of Fast Fourier Transform;4） Calculation of the designed fundamental current and phase angle between voltage and current from rated electric values of the device: voltage, pow-ers P and Q ;5） Pre-design of the harmonic spectrum at rated current and creation of a voltage controlled current source model for frequency domain simulation J 0(f ), see example in Fig. 3 for the used current pattern of a 12-pulse rectifier.Fig. 3 Current Waveform and Harmonic Spectrumof a 12-Pulse RectifierDefine the source impedance Z (f ) as a function of frequency. Z (f ) reflects the device behaviour to the disturbances coming from power supply side.For particular case, it is very important to define Z (f ) values within the studied frequencies. For instance, a static converter’s impedance was set an infinite value in the past in performing inter-harmonic studies. This is not true for a number of static converters. Some in-dustrial cases show the ripple control signal at 175 Hz can be partially absorbed by a rectifier with a capacitor filter as in a variable-speed drive [2].The proposed harmonic source model behaviours as follows: when the source J (f ) is activated, imped-ance Z (f ) is naturally set to an infinite value because the device behaviours as a current source. When the source J (f ) is not activated, Z (f ) has to be added into the system admittance matrix. In this case, the device behaviours as a passive component. Value of Z (f ) is defined case by case in accordance with studied phe-nomena:At fundamental frequency, Z (f ) is calculated from rated voltage and powers P , Q of the simulated device.第3期杨显军：法国电网电能质量承诺和电能质量评估 11At studied harmonic or inter harmonic frequen-cies whose source comes from the power supply side, Z(f) can be defined from the real response of the de-vice at these frequencies. Z(f) is depending on device structure, it can be modelled by series (or parallel) RL, RC circuits, by mixed model of CIGRE 36.05[3] or by a data table for irregular cases.If it is difficult to define a general analytical for-mula representing Z(f) for all studied frequencies. For studying a specific frequency, a simple solution is to perform a time domain simulation on the power elec-tronic device in order to obtain the impedance re-sponse. For example, at ripple control frequency of 175 Hz, the equivalent impedance at input of a recti-fier changes with the type of the rectifier structure, DC bus filter, AC filter and rectifier control techniques. A time domain simulation shows that for a 400 V, 250 kV A diode rectifier used in a variable-speed drive, the impedance at 175 Hz (Z175) is mainly determined by the DC bus type: Z175 = (0.1 ~ 0.9) × Z50 for a capaci-tive DC filter, and Z175 = (0.8 ~ 3) × Z50 for an induc-tive DC filter (Z50 is the rated fundament equivalent impedance) [2].2.3 HVDC Link Frequency Domain ModelAn High V oltage DC link (HVDC) model has been studied in frequency domain and a plug & play HVDC component has been created by using studied harmonic source models. A general HVDC link is composed of two power electronic devices (converter stations): one works as rectifier and another as inverter (or generator), see Fig. 4.Fig. 4 HVDC Link Structure between Two GridsThe studied HVDC component is decomposed in Fig. 5. Each side of HVDC is composed of a three-phase controlled harmonic current source with 74 harmonic number, an impedance and a regulator for power regulation. Two current sources are independent and each current source is synchronized with con-nected fundamental voltage after fundamental fre-quency load flow. This model can simulate independ-ently power quality and electromechanical transient phenomena of each side of HVDC link such as har-monic, reactive power, voltage fluctuation, etc. One case study below will show the application of pro-posed model.Fig. 5 Structure of Frequency Domain HVDC Model2.4 Power Quality Meter ModellingFor a power quality software, PQ meter algorithm is also a key component to think about. A general power quality software that gives just the harmonic spectrum is not enough for a site engineer to perform grid disturbance assessment in a short delay. In fact, the software needs to spend a lot of time in post-data processing in order to get necessary PQ values such as inter-harmonic, harmonic grouping, flicker, unbalance, voltage sag, etc. In most commercialised time domain software, there is no integrated power quality meter. Users have to create tool box in post data processing in order to perform grid PQ assessment. Generally, this step is time-consuming.Simulation tools for grid power quality assess-ment should include at least a build-in total harmonic distortion and individual harmonic computation unit. An advanced power quality simulation tool has to in-tegrate other disturbance metering such as in-ter-harmonic, harmonic grouping, IEC flicker meter, etc [4−5]. Tab. 6 gives some electric measurements to be integrated into power quality simulation tools.An advanced power quality meter algorithm has been studied and implanted in our power quality12 南方电网技术第3卷Tab. 6 Electric Measurements Concerning Power QualityBasic Power Quality Values Advanced Power Quality Valuestotal harmonic distortion; individual harmonic magnitude; armonic angle;voltage sag. inter-harmonic level; harmonic grouping; harmonic power; unbalance;flicker level (P st and P lt).software [6−7]. Fig. 6 shows a case study including harmonic and flicker assessments with embedded dis-turbance analyser and IEC flicker meter in a frequency domain software with electromechanical transient analysis. Harmonic distortion and flicker level Pst are displayed directly on the main simulation.Fig. 6 Disturbance Assessment with Embedded PowerQuality Meters in a Simulation Tool3 Case-Studies: Application & ValidationThe developed harmonic source models and power quality meters have been integrated in a power quality simulation tool, which has been used for grid and customer power quality analysis in a number of actual case studies. Two of them are presented here.3.1 Case Study 1: HVDC Link IntegrationA case study has been performed on harmonic level assessment and voltage improvement upon future HVDC integration to two existing grids. The study has been done using software ExpertEC [8−9] (a frequency domain software with electro-mechanical transient analysis) embedded with the developed numeric HVDC model and power quality meters. This software reads directly the substation databases and converts relative grid data into electric circuit models (Fig. 7). Current harmonic patterns at each HVDC converter station have been modelled by both simplified theo-retical data and on-site measurements offered by a manufacturer. The studied HVDC link is inserted be-tween grid 1 and grid 2 modelled by the following parameters:1）Grid 1 — short-circuit power 5 GV A, total number of nodes 1 100;2）Grid 2 — short-circuit power 1 GV A, total number of nodes 320;3）200 MW of HVDC power flow from grid 1 to grid 2;4）Frequency domain simulation: 50 ~ 3 150 Hz;5）Time step in voltage regulation 20 ms.Fig. 7 Simulated Harmonic Levels in Two GridsHVDC Option 1: Thyristor-based HVDC. Con-verter station is composed of 12-pulse rectifier + sim-plified harmonic filters with current harmonic distor-tion THD = 3.2% (Fig. 8).Fig. 8 Current Waveform and Harmonic Spectrum at Input ofThyristor-Based HVDC Converter StationHVDC Option 2: In the same grid conditions, the two converter stations are replaced by IGBT-based con-verters + harmonic filters. In this case, the current har-monic distortion TDH will be lower than 1.2% (Fig. 9).Fig. 9 Current Waveform and Harmonic Spectrum at Input ofIGBT-Based HVDC Converter Station第3期杨显军：法国电网电能质量承诺和电能质量评估 13For these two HVDC options, voltage harmonicassessment has been performed in all nodes of twogrids from 50 to 3 150 Hz. The simulation resultsshow after insertion of the HVDC option 1 (Thyris-tor-based HVDC with simplified AC filters), the totalvoltage harmonic distortion in two grids becomes toohigh (THD > 3.79%) compared to the planning level ofTHD < 3% defined by IEC 61000-3-6. The simplifiedfilters are not enough to reduce current harmonics in this case. The results from HVDC option 2 (IGBT-based HVDC) gives satisfaction: voltage har-monic levels in two grids are 0.35% for grid 1 and 0.71% for grid 2. These results mean that the integra-tion of this HVDC link will not deteriorate substan-tially harmonic levels of the two grids. This HVDC solution could be taken in consideration.Tab. 7 Total Voltage Harmonic Distortion (225 kV Level)Name of Grids Option 1, THD Option2,THD Grid 1 3.79% 0.35%Grid 2 5.91% 0.71%3.2 Case Study 2: Harmonic Current Assess-ment at PCCAnother case study has been performed on a point of common coupling (PCC) of two industrial sites. Site 1 has an important harmonic source (3 MW, power electronic devices). Site 2 has mainly linear load (induction heating by fundamental frequency) no important harmonic source but there is a very impor-tant capacitor bank for var compensation (Fig. 10).Some simulations have been done in order to confirm site engineer’s first assumption: 550 Hz har-monic source is inside site 1 although at PCC the har-monic current measured at feeder 1 is lower than that of feeder 2. Former measurements show that upstream back ground harmonic voltage at PCC is negligible if these two sites are switched off. Simulation results show that overall harmonic power flow at feeder 1 comes from downstream as the harmonic power is negative (Fig. 10). And the overall harmonic power flow at feeder 2 comes from upstream at the feeder 2 (the harmonic power is positive).Fig. 10 Simulation of Harmonic Resonance at 550 HzThis phenomenon is caused by harmonic reso-nance around 550 Hz between capacitor bank and up-stream impedance. Fig. 11 shows harmonic impedance magnitude at PCC.Fig. 11 Harmonic Impedance Magnitude at PCC4 ConclusionThis paper gives an overview of power quality commitments in grid and customer sides. Two con-tracts have been put into service in France between utilities and customers. These power quality contracts set not only the power quality commitments for the grid operator but also on the customer disturbance emission level.Generally, frequency domain software is more adequate in doing power quality assessment than time domain software. Two key points should be taken into account in power quality simulation tools: disturbance load models and power quality meter in order to re-duce computation and post processing time.Power quality meter models and three phase harmonic source models have been developed. By means of advanced build-in non-linear components, it is possible to study grid power quality assessment by14 南方电网技术第3卷dealing with non-linear loads either as disturbance sources or as disturbance receivers. As an extension of introduced general harmonic source model, an HVDC component has been built, which is composed of two independent harmonic sources and two regulators.A number of multi phase harmonic source models have been developed and integrated into a power qual-ity software, which has been used in grid harmonic analysis. Two case-studies have shown applications of harmonic source and power quality meters in actual grid power quality assessment.References:[1]NASLIN C, GONBEAU O, FRAISSE J L, et al. Wind Farm Integra-tion — Power Quality Management for French Distribution Network[C]// The 18th International Conference and Exhibition on ElectricityDistribution (CIRED 2005), June 6−9, 2005, Turin, Italy. London:IEE Conference Publications, 2005, 2: 41.[2]YANG X, DENNETIÈRE S. Modeling of the Behavior of PowerElectronic Equipment to Grid Ripple Control Signal [C]// The 7thInternational Conference on Power Systems Transients (IPST 2007),June 4−7, 2007, Lyon, France. [S.l.]: IPST, 2007.[3]O’NEILL-CARRILLO E, HEYDT G T, KOSTELICH E J, et al.Nonlinear Deterministic Modeling of Highly Varying Loads [J].IEEE Transactions on Power Delivery, 1999, 14 (2): 537−542.[4]IEC 61000-4-30:2008, Electromagnetic Compatibility (EMC): Test-ing and Measurement Techniques — Power quality measurementmethods [S].[5]IEC 61000-4-7:2002, Electromagnetic Compatibility (EMC): Testingand Measurement Techniques — General Guide on Harmonics andInterharmonics Measurements and Instrumentation, for Power Sup-ply Systems and Equipment Connected Thereto [S]. [6]IEC 61000-4-15:2003, Electromagnetic Compatibility (EMC): Test-ing and Measurement Techniques — Flickermeter — Functional andDesign Specifications [S].[7]YANG X. Methodology and Simulation Tool for Flicker PropagationAnalysis [C]// EPRI Power Quality and Applications (PQA) 2006and Advanced Distribution Automation (ADA) Joint Conference,Marriott Marquis, Atlanta, Georgia, USA, July 24−26, 2006.[8]YANG X, BERTHET L. Methodology for Frequency Domain Dis-turbance Source Modelling for Multi Source Harmonic and In-ter-harmonic Studies [C]// The 12th International Conference onHarmonics and Quality of Power (ICHQP), Cascais, Portugal, Octo-ber 1–5, 2006.[9]YANG X. Advanced Methodology and New Tool for MultiphasePower Quality Analysis & Mitigation [C]// 18th International Con-ference on Electricity Distribution (CIRED 2005), Torino, Italy, June 6−9, 2005.Received date：2009-02-28Biography：Xavier X. Yang (1958−), Male. Dr. Yang received his Electrical Engi-neering Diploma B.S. and M.S. degrees in Liaoning Technology University (China), and received the Ph.D. from National Polytechnic Institute of Tou-louse (France). Dr. Yang was employed in French Electricity and Energy Consulting Companies in which he was in charge of site measurement and power quality mitigation. He is now expert researcher in R&D of Electricite de France (EDF). His research interests are in the areas of electric power system modeling, power system harmonics, computer algorithms, fault detection and voltage sag mitigation. He is a member of IEC SC77A Work-ing Group 1 (Harmonic and flicker) and CIGRE Joint Working Groups B4 (HVDC) and C4 (Review of flicker objectives).E-mail: xavier.yang@edf.fr.收稿日期：2009-02-28作者简介：杨显军（1958–)，男。

一种求解最优潮流的改进灰狼优化算法

第１４卷㊀第３期Ｖｏｌ．１４Ｎｏ．３㊀㊀智㊀能㊀计㊀算㊀机㊀与㊀应㊀用ＩｎｔｅｌｌｉｇｅｎｔＣｏｍｐｕｔｅｒａｎｄＡｐｐｌｉｃａｔｉｏｎｓ㊀㊀２０２４年３月㊀Ｍａｒ．２０２４㊀㊀㊀㊀㊀㊀文章编号：２０９５－２１６３（２０２４）０３－００４６－０８中图分类号：ＴＥ３４１文献标志码：Ａ一种求解最优潮流的改进灰狼优化算法王㊀恒，杨㊀婷（铜仁职业技术学院信息工程学院，贵州铜仁５５４３００）摘㊀要：最优潮流是电力系统最关键的问题之一，本文采用一种求解最优潮流的改进灰狼优化算法（ＬＭＧＷＯ）求解最优潮流（ＯＰＦ）问题，该算法引入算术优化算法（ＡｒｉｔｈｍｅｔｉｃＯｐｔｉｍｉｚａｔｉｏｎＡｌｇｏｒｉｔｈｍ，ＡＯＡ）中的乘除算子，利用带透镜成像的反向学习策略增强最优个体的多样性，提高算法跳出局部最优的能力㊂通过与几种常用的算法进行对比实验表明：本文提出的ＬＷＧ⁃ＷＯ算法是有竞争力的，总体上优于对比算法；ＬＭＧＷＯ算法在最小化燃料成本㊁有功输电损耗和改善电压偏差方面更有效地找到了最优潮流（ＯＰＦ）问题的最优解㊂关键词：灰狼优化算法；最优潮流；算术优化算法；燃料成本；有功输电损耗ＡｎｉｍｐｒｏｖｅｄｇｒｅｙｗｏｌｆｏｐｔｉｍｉｚａｔｉｏｎａｌｇｏｒｉｔｈｍｆｏｒｓｏｌｖｉｎｇｏｐｔｉｍａｌｐｏｗｅｒｆｌｏｗＷＡＮＧＨｅｎｇ，ＹＡＮＧＴｉｎｇ（ＳｃｈｏｏｌｏｆＩｎｆｏｒｍａｔｉｏｎＥｎｇｉｎｅｅｒｉｎｇ，ＴｏｎｇｒｅｎＰｏｌｙｔｅｃｈｎｉｃＣｏｌｌｅｇｅ，Ｔｏｎｇｒｅｎ５５４３００，Ｇｕｉｚｈｏｕ，Ｃｈｉｎａ）Ａｂｓｔｒａｃｔ：Ｏｐｔｉｍａｌｐｏｗｅｒｆｌｏｗｉｓｏｎｅｏｆｔｈｅｍｏｓｔｃｒｉｔｉｃａｌｐｒｏｂｌｅｍｓｉｎｐｏｗｅｒｓｙｓｔｅｍ．Ｉｎｔｈｉｓｐａｐｅｒ，ａｎｉｍｐｒｏｖｅｄＧｒｅｙＷｏｌｆＯｐｔｉｍｉｚａｔｉｏｎＡｌｇｏｒｉｔｈｍ（ＬＭＧＷＯ）ｉｓｕｓｅｄｔｏｓｏｌｖｅｔｈｅｏｐｔｉｍａｌｐｏｗｅｒｆｌｏｗ（ＯＰＦ）ｐｒｏｂｌｅｍ．Ｉｎｔｈｉｓａｌｇｏｒｉｔｈｍ，ｍｕｌｔｉｐｌｉｃａｔｉｏｎａｎｄｄｉｖｉｓｉｏｎｏｐｅｒａｔｏｒｓｉｎｔｈｅＡｒｉｔｈｍｅｔｉｃＯｐｔｉｍｉｚａｔｉｏｎＡｌｇｏｒｉｔｈｍ（ＡＯＡ）ａｒｅｉｎｔｒｏｄｕｃｅｄ．Ｔｈｅｒｅｖｅｒｓｅｌｅａｒｎｉｎｇｓｔｒａｔｅｇｙｗｉｔｈｌｅｎｓｉｍａｇｉｎｇｉｓｕｓｅｄｔｏｅｎｈａｎｃｅｔｈｅｄｉｖｅｒｓｉｔｙｏｆｏｐｔｉｍａｌｉｎｄｉｖｉｄｕａｌｓａｎｄｉｍｐｒｏｖｅｔｈｅａｂｉｌｉｔｙｏｆｔｈｅａｌｇｏｒｉｔｈｍｔｏｊｕｍｐｏｕｔｏｆｔｈｅｌｏｃａｌｏｐｔｉｍａｌ．Ｔｈｒｏｕｇｈｃｏｍｐａｒａｔｉｖｅｅｘｐｅｒｉｍｅｎｔａｌａｎａｌｙｓｉｓｏｆｓｅｖｅｒａｌｃｏｍｍｏｎｌｙｕｓｅｄａｌｇｏｒｉｔｈｍｓ，ｔｈｅｐｒｏｐｏｓｅｄＬＷＧＷＯａｌｇｏｒｉｔｈｍｉｓｃｏｍｐｅｔｉｔｉｖｅａｎｄｇｅｎｅｒａｌｌｙｓｕｐｅｒｉｏｒｔｏｒｅｃｅｎｔａｌｇｏｒｉｔｈｍｓ．ＴｈｅｅｘｐｅｒｉｍｅｎｔａｌｒｅｓｕｌｔｓｓｈｏｗｔｈａｔＬＭＧＷＯａｌｇｏｒｉｔｈｍｃａｎｆｉｎｄｔｈｅｏｐｔｉｍａｌｓｏｌｕｔｉｏｎｏｆＯＰＦｐｒｏｂｌｅｍｍｏｒｅｅｆｆｅｃｔｉｖｅｌｙｉｎｔｅｒｍｓｏｆｍｉｎｉｍｉｚｉｎｇｆｕｅｌｃｏｓｔ，ａｃｔｉｖｅｐｏｗｅｒｔｒａｎｓｍｉｓｓｉｏｎｌｏｓｓａｎｄｉｍｐｒｏｖｉｎｇｖｏｌｔａｇｅｄｅｖｉａｔｉｏｎ．Ｋｅｙｗｏｒｄｓ：ｇｒｅｙｗｏｌｆｏｐｔｉｍｉｚａｔｉｏｎａｌｇｏｒｉｔｈｍ；ｏｐｔｉｍａｌｐｏｗｅｒｆｌｏｗ；ａｒｉｔｈｍｅｔｉｃｏｐｔｉｍｉｚａｔｉｏｎａｌｇｏｒｉｔｈｍ；ｆｕｅｌｃｏｓｔ；ａｃｔｉｖｅｐｏｗｅｒｔｒａｎｓｍｉｓｓｉｏｎｌｏｓｓ基金项目：铜仁市科学技术局基础科学研究项目（铜市科研（２０２２）７２号）㊂作者简介：王㊀恒（１９８５－），男，博士研究生，讲师，主要研究方向：智能计算与混合系统㊁人工智能㊁故障诊断研究等㊂Ｅｍａｉｌ：ｗａｎｇｈｅｎｇ＿ｔｒｚｙ＠ｆｏｘｍａｉｌ．ｃｏｍ收稿日期：２０２３－０６－１６０㊀引㊀言最优潮流（ＯＰＦ）问题是电力系统运行过程中备受关注的焦点问题，旨在找到最优的运行方式，使得电力系统的运行成本最低，同时满足安全㊁稳定和环保等约束条件㊂ＯＰＦ问题的求解是在满足一系列物理㊁环境㊁实际和运行的约束条件下，通过优化特定的目标来确定电力系统的运行状态㊂在此之前，许多传统的优化技术的应用已获成功，包括基于梯度的方法㊁牛顿法㊁单纯形法㊁序列线性规划和内点法［１－５］㊂由于ＯＰＦ问题本质上是一个多极㊁多约束㊁非凸的复杂优化问题，使用传统的数值方法来求解，过程复杂㊁耗时且精度较差㊂近年来，元启发式算法的快速发展为解决ＯＰＦ问题提供了更多的选择㊂元启发式算法具有参数少㊁易于操作㊁不需要梯度信息等优点，能够在合理的时间内和高度复杂的约束条件下找到复杂问题的最优解㊂刘自发等学者［６］提出了一种基于混沌粒子群优化方法的电力系统无功最优潮流（ＯＰＦ）问题㊂Ｆａｒｈａｔ等学者［７］提出了一种基于邻域维度学习搜索策略的增强型黏液霉菌算法（ｅｎｈａｎｃｅｄｓｌｉｍｅｍｏｕｌｄａｌｇｏｒｉｔｈｍ，ＥＳＭＡ）用于求解最优潮流（ＯＰＦ）问题等等㊂越来越多的元启发式算法被广泛用于解决电力系统优化相关问题［８－１３］㊂灰狼优化算法（ｇｒｅｙｗｏｌｆｏｐｔｉｍｉｚｅｒ，ＧＷＯ）是由Ｍｉｒｊａｌｉｌｉ等学者［１４］在２０１４年上提出的一种新的元启发式算法㊂灰狼优化算法（ＧＷＯ）原理简单㊁编程容易㊁需要调整的参数少，现已陆续应用于电力系统㊁自动控制㊁能源市场战略招标等领域［１５－１７］㊂然而，与许多元启发式优化算法一样，灰狼优化算法（ＧＷＯ）在求解复杂的非线性问题时容易陷入局部最优且收敛速度慢㊂针对原有灰狼优化算法在求解最优潮流（ＯＰＦ）问题时存在的不足，提出了一种改进的灰狼优化算法（ＬＭＧＷＯ算法）㊂基于镜头成像学习和乘除算子策略对原灰狼优化算法（ＧＷＯ）进行改进，主要有２点改进：（１）为了增强算法的全局探索能力，引入乘除算子策略，提高算法的收敛速度；（２）为增强最优个体的多样性，引入透镜成像修正反向学习策略，提高算法跳出局部最优的能力㊂１㊀最优潮流公式最优潮流（ＯＰＦ）问题是典型的多变量㊁多约束的非线性组合优化问题㊂最优潮流（ＯＰＦ）问题的求解过程是通过寻找最优的控制变量来获得最小的目标函数㊂数学模型定义如下：ｍｉｎＦ（ｕ，ｘ）ｓ．ｔ．ｇ（ｕ，ｘ）＝０ｈ（ｕ，ｘ）ɤ０{㊀㊀其中，Ｆ表示目标函数；ｘ表示控制变量；ｕ表示状态变量；ｇ（ｕ，ｘ）＝０是等式约束；ｈ（ｕ，ｘ）ɤ０是不等式约束㊂１．１㊀控制变量和状态变量最优潮流（ＯＰＦ）问题公式中的控制变量集合为：㊀㊀ｘ＝［ＰＧ２，，ＰＧＮＧ，ＶＧ１，，ＶＧＮＧ，Ｔ１，，ＴＮＴ，ＱＣ１，，ＱＣＮＣ］（１）其中，ＰＧ２，，ＰＧＮＧ为系统除松弛母线外的有功发电量；ＶＧ１，，ＶＧＮＧ为系统的电压幅值；Ｔ１，，ＴＮＴ为变压器分接设定值；ＱＣ１，，ＱＣＮＣ为并联无功补偿；ＮＧ㊁ＮＴ㊁ＮＣ分别为发电机个数㊁调节变压器个数㊁无功补偿器个数㊂最优潮流（ＯＰＦ）问题表述的状态变量集合为：ｕ＝［ＰＧ１，ＶＬ１，，ＶＬＮＬ，ＱＧ１，，ＱＧＮＧ，Ｓｌ１，，Ｓｌｎｌ］（２）其中，ＰＧ为空闲母线输出有功功率；ＶＬ为负载母线电压幅值；ＱＧ为各发电机组输出无功功率；Ｓｌ为输电线路负载㊂１．２㊀目标函数将燃油成本㊁有源输电损耗和电压偏差作为最优潮流（ＯＰＦ）问题的目标函数㊂各目标函数的数学模型定义如下㊂（１）燃料成本（ＦＣ）㊂描述发电成本的目标函数，可得数学建模如下：Ｆ１（ｘ，ｕ）＝ðＮｇｉ＝１（ａｉ＋ｂｉＰＧｉ＋ｃｉＰ２Ｇｉ）（３）㊀㊀其中，Ｎｇ为发电机个数；ａｉ，ｂｉ，ｃｉ为第ｉ台发电机组的燃料成本系数；ＰＧｉ为第ｉ台发电机组的实际发电量㊂（２）有功输电损耗（ＡＰＬ）㊂传输线的ＡＰＬ可表示为：㊀Ｆ２（ｘ，ｕ）＝ðｉ，ｊɪＮｌＧｉｊＶ２ｉ＋Ｖ２ｊ－２ＶｉＶｊｃｏｓ（θｉｊ）()（４）㊀㊀其中，Ｎｌ为输电线路数；Ｇｉｊ为线路ｉｊ的传递电导；Ｖｉ为第ｉ根母线的电压幅值；Ｖｊ为第ｊ根母线的电压幅值；θｉｊ为母线ｉ与ｊ之间的电压相角之差㊂１．３㊀约束条件在最优潮流（ＯＰＦ）问题中，等式约束和不等式约束是电力系统需要满足的约束，通常是每个节点的功率平衡约束，可以通过式（５）和式（６）进行定义：ＰＧｉ－ＰＤｉ＝ＶｉðＮｉ，ｊ＝１Ｖｊ（Ｇｉｊｃｏｓ（δｉ－δｊ）＋Ｂｉｊｓｉｎ（δｉ－δｊ））（５）ＱＧｉ－ＱＤｉ＝ＶｉðＮｉ，ｊ＝１Ｖｊ（Ｇｉｊｓｉｎ（δｉ－δｊ）－Ｂｉｊｃｏｓ（δｉ－δｊ））（６）其中，ＰＤｉ㊁ＱＤｉ分别为第ｉ台母线的有功㊁无功功率；ＰＧｉ和ＱＧｉ为第ｉ台发电机的无功发电量；Ｎ为母线个数；Ｇｉｊ和Ｂｉｊ分别为母线ｉ和ｊ之间的电导和电纳；Ｖｉ和Ｖｊ分别为母线ｉ和ｊ的电压幅值㊂２㊀改进的灰狼优化算法２．１㊀灰狼优化算法灰狼优化算法（ＧＷＯ）是模仿自然界灰狼群体社会等级和捕食行为而衍生的一种元启发式算法［１４］㊂灰狼群体的社会等级为α狼㊁β狼㊁δ狼和ω狼㊂狼的狩猎行为分为跟踪㊁包围和攻击猎物三个步骤㊂狼群包围猎物的数学模型定义为：Ｘ＝Ｘα（ｔ）－Ａ㊃｜Ｃ㊃Ｘα（ｔ）－Ｘ（ｔ）｜（７）㊀㊀其中，Ｘ和Ｘα分别表示狼个体和猎物个体的位置向量，ｔ表示当前迭代次数㊂系数向量Ａ和Ｃ定义为：Ａ＝２ａ㊃ｒ１－ａ（８）Ｃ＝２㊃ｒ２（９）㊀㊀其中，ｒ１和ｒ２是［０，１］之间的随机向量，ａ从２线性递减到０，其数学模型定义为：７４第３期王恒，等：一种求解最优潮流的改进灰狼优化算法ａ＝２－２㊃ｔＴｍａｘ（１０）㊀㊀其中，Ｔｍａｘ为最大迭代次数㊂包围猎物后，β狼和δ狼在α狼的带领下追捕猎物㊂在追捕过程中，狼群的个体位置会随着猎物的逃跑而发生变化㊂因此，灰狼群可以根据α㊁β㊁δ的位置Ｘα，Ｘβ，Ｘδ更新灰狼的位置：Ｘ１＝Ｘα（ｔ）－Ａ１㊃｜Ｃ１㊃Ｘα（ｔ）－Ｘ（ｔ）｜（１１）Ｘ２＝Ｘβ（ｔ）－Ａ２㊃｜Ｃ２㊃Ｘβ（ｔ）－Ｘ（ｔ）｜（１２）Ｘ３＝Ｘδ（ｔ）－Ａ３㊃｜Ｃ３㊃Ｘδ（ｔ）－Ｘ（ｔ）｜（１３）Ｘ（ｔ＋１）＝Ｘ１＋Ｘ２＋Ｘ３３（１４）㊀㊀其中，Ｘ（ｔ＋１）是当前个体的位置㊂２．２㊀改进ＧＷＯ算法的思路和策略２．２．１㊀算术乘除运算符策略２０２１年，Ａｂｕａｌｉｇａｈ等学者［１８］提出的一种新的元启发式算法，即算术优化算法（ＡｒｉｔｈｍｅｔｉｃＯｐｔｉｍｉｚａｔｉｏｎＡｌｇｏｒｉｔｈｍ，ＡＯＡ），主要利用数学中的乘㊁除运算符以及加㊁减运算符四种混合运算㊂ＡＯＡ中的乘除算子具有较强的全局探索能力㊂灰狼种群在更新位置时侧重使用α狼㊁β狼和δ狼作为精英来引导搜索，具有较强的局部开发能力㊂引入算术乘除算子策略，提高ＧＷＯ算法的全局探索能力㊂算术乘除算子策略的数学模型定义为：Ｘｊｉ（ｔ＋１）＝Ｘｊｂｅｓｔː（ＭＯＰ＋ε）㊃［（ｕｂｊ－ｌｂｊ）㊃μ＋ｌｂｊ］，㊀ｒ３ɤ０．５ＸｊｂｅｓｔˑＭＯＰ㊃［（ｕｂｊ－ｌｂｊ）㊃μ＋ｌｂｊ］，㊀㊀㊀㊀ｒ３＞０．５{（１５）㊀㊀其中，Ｘｊｂｅｓｔ表示当前最优解的第ｊ个位置；ｒ３表示介于［０，１］之间的随机数；ε表示防止分母为０的整数；μ表示调节搜索过程的控制参数，μ的值在基本ＡＯＡ中为０．５；ｕｂｊ和ｌｂｊ分别表示第ｉ个位置的上下界㊂ＭＯＰ为概率函数，其数学模型描述为：ＭＯＰ＝１－ｔ１τＴ１τｍａｘ（１６）㊀㊀其中，τ＝５是一个敏感因子，定义了迭代的搜索精度㊂由式（１５）可知，ＡＯＡ可以带来高分布，借助乘除算子实现位置更新，可以大大提高算法的全局探索能力㊂本文设置阈值为０．３㊂２．２．２㊀基于透镜成像的反向学习策略根据灰狼的位置更新公式，由α狼㊁β狼和δ狼带领群体中的其他狼进行位置更新㊂如果α狼㊁β狼和δ狼都处于局部最优，则整个群体会聚集在局部最优区域，导致种群陷入局部最优㊂针对该问题，本文提出一种基于透镜成像原理的反向学习方法，将对立个体与当前最优个体相结合，生成新个体㊂假设在一维空间中，在轴区间［ｌｂ，ｕｂ］上有一个高度为Ｈ的个体Ｐ，其在ｘ轴上的投影为Ｘ（Ｘ为全局最优个体）㊂将焦距为Ｆ的镜头放置在基点位置Ｏ上（本文取基点位置为（ｌｂ＋ｕｂ／２））㊂个体Ｐ通过透镜，以获得高度为Ｈ的倒置图像Ｐ∗，在这点上，第一个倒置的个体ｘ通过透镜成像在Ｘ轴上产生㊂镜头图像的反向学习策略如图１所示㊂㊀㊀在图１中，全局最优个体Ｘ以Ｏ为基点找到其对应的逆个体Ｘ∗㊂因此，可以从透镜成像原理推导出数学模型，推得的公式为：（ｕｂ＋ｌｂ）／２－ＸＸ∗－（ｕｂ＋ｌｂ）／２＝ｈｈ∗（１７）㊀㊀设ｈ／ｈ∗＝ｋ，ｋ表示拉伸因子㊂通过推导式（１７），可以得到反转点Ｘ∗的计算公式：Ｘ∗＝ｕｂ＋ｌｂ２＋ｕｂ＋ｌｂ２ｋ－Ｘｋ（１８）xOh PXl bu b h*X *P*yF图１㊀基于镜头图像的反向学习策略Ｆｉｇ．１㊀Ｒｅｖｅｒｓｅｌｅａｒｎｉｎｇｓｔｒａｔｅｇｙｂａｓｅｄｏｎｌｅｎｓｉｍａｇｅ㊀㊀在算法搜索解时，使用拉伸因子ｋ作为微观调节因子，增强算法的局部开发能力㊂然而，在基本的透镜成像逆学习策略中，拉伸因子一般作为固定值使用，不允许算法探索解空间的全范围㊂为此，本文提出一种基于非线性动态递减的伸缩因子策略，在算法迭代初期可以得到较大的值，有助于算法在不同维度的区域进行更大范围的搜索，以提高种群的多样性㊂非线性动态拉伸因子定义为：㊀ｋ＝ｋｍａｘ－（ｋｍａｘ－ｋｍｉｎ）㊃［１－ｃｏｓ（πｔ２Ｔｍａｘ）］（１９）㊀㊀其中，ｋｍａｘ和ｋｍｉｎ分别表示最大和最小拉伸因子，Ｔｍａｘ表示最大迭代次数㊂可以将式（１８）扩展到Ｄ－维搜索空间，得到数学模型为：８４智㊀能㊀计㊀算㊀机㊀与㊀应㊀用㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀第１４卷㊀Ｘ∗ｊ＝ｕｂｊ＋ｌｂｊ２＋ｕｂｊ＋ｌｂｊ２ｋ－Ｘｊｋ（２０）㊀㊀其中，Ｘｊ和Ｘ∗ｊ分别表示Ｘ和Ｘ∗的的第ｊ维向量，ｕｂｊ和ｌｂｊ分别表示决策变量的第ｊ维向量㊂基于透镜的反向学习策略虽然极大地提高了算法的求解精度，但无法直接判断生成的新反向个体是否优于原始个体㊂因此，本文引入贪心机制来比较新旧个体适应度值，从而筛选出最优个体㊂该方法不断获得更好的解，提高了算法的寻优能力㊂贪婪机制的数学模型描述如下：Ｘｎｅｗ（ｔ）＝Ｘ∗，㊀ｆ（Ｘ）＞ｆ（Ｘ∗）Ｘ，㊀ｆ（Ｘ）ɤｆ（Ｘ∗）{（２１）２．２．３㊀ＬＭＧＷＯ算法实现过程ＬＭＧＷＯ算法实现流程如图２所示㊂计算每只狼的适应度，从狼群中选出α狼、β狼和δ狼开始初始化狼群的位置t =t +1i f t ＜T m a x 结束运行式(19)~（22）执行基于透镜成像的反向学习策略i f r ＜0.3通过式(17）、式(18）执行算术乘除运算符策略通过式(13)~(16)更新狼群的位置计算适应度值更新向量α狼、β狼和δ狼图２㊀ＬＭＧＷＯ算法流程图Ｆｉｇ．２㊀ＦｌｏｗｃｈａｒｔｏｆＬＭＧＷＯａｌｇｏｒｉｔｈｍ３㊀实验３．１㊀实验环境及参数设置在Ｉｎｔｅｌ（Ｒ）Ｃｏｒｅ（ＴＭ）ｉ７－ｉ７－６５００ＵＣＰＵ㊁２．５０ＧＨｚ频率㊁８ＧＢ内存㊁Ｗｉｎｄｏｗｓ１０（６４ｂｉｔ）操作系统上进行仿真实验，编程软件为ＭａｔｌａｂＲ２０１８ａ㊂采用９个基准测试函数，包括５个单峰函数Ｆ１Ｆ５和４个非线性多峰函数Ｆ６Ｆ９，见表１㊂参与对比的灰狼优化算法（ＧＷＯ）［１４］㊁算术优化算法（ＡＯＡ）［１８］㊁正弦余弦算法（ＳＣＡ）［１９］㊁猩猩优化算法（ＣｈＯＡ）［２０］㊁鲸鱼优化算法（ＷＯＡ）［２１］㊁ＬＭＧＷＯ的参数设置见表２㊂表１㊀基准测试函数Ｔａｂｌｅ１Ｂｅｎｃｈｍａｒｋｆｕｎｃｔｉｏｎｓ函数编号名称维度范围最优值Ｆ１Ｓｐｈｅｒｅ３０［－１００，１００］０Ｆ２Ｓｃｈｗｅｆｅｌ．２．２２３０［－１０，１０］０Ｆ３Ｓｃｈｗｅｆｅｌ．１．２３０［－１００，１００］０Ｆ４Ｓｃｈｗｅｆｅｌ．２．２１３０［－１００，１００］０Ｆ５Ｑｕａｒｔｉｃ３０［－１．２８，１．２８］０Ｆ６Ｒａｓｔｒｉｇｉｎ３０［－５．１２，５．１２］０Ｆ７Ａｃｋｌｅｙ３０［－３２，３２］０Ｆ８Ｃｒｉｅｗａｎｋ３０［－６００，６００］０Ｆ９Ａｐｌｉｎｅ３０［－１０，１０］０９４第３期王恒，等：一种求解最优潮流的改进灰狼优化算法表２㊀算法参数设置Ｔａｂｌｅ２㊀Ｐａｒａｍｅｔｅｒｓｅｔｔｉｎｇｓｏｆａｌｇｏｒｉｔｈｍｓ算法名称参数设置ＳＣＡ［１９］Ｍ＝２ＣｈＯＡ［２０］ｆｍａｘ＝２．５，ｆｍｉｎ＝０ＷＯＡ［２１］ａｍａｘ＝２，ａｍｉｎ＝０，ｂ＝１ＡＯＡ［１８］ＭＯＰ＿Ｍａｘ＝１，ＭＯＰ＿Ｍｉｎ＝０．２，α＝５，μ＝０．４９９ＧＷＯ［１４］ａｍａｘ＝２，ａｍｉｎ＝０ＬＭＧＷＯａｍａｘ＝２，ａｍｉｎ＝０３．２㊀算法性能对比分析为了验证了ＬＭＧＷＯ算法的有效性和优越性，将ＬＭＧＷＯ算法与灰狼优化算法（ＧＷＯ）［１４］㊁算术优化算法（ＡＯＡ）［１８］㊁正弦余弦算法（ＳＣＡ）［１９］㊁猩猩优化算法（ＣｈＯＡ）［２０］㊁鲸鱼优化算法（ＷＯＡ）［２１］在９个不同特性的基准测试函数上进行仿真实验㊂在各个算法的测试环境相同的条件下，种群规模Ｎ＝３０，空间维度Ｄｉｍ＝３０，最大迭代次数Ｔｍａｘ＝５００㊂采用均值和标准差作为实验的评价指标，均值和标准差越小，表明算法的性能越好㊂６种算法对９个基准函数的求解结果见表３㊂表３㊀各算法在基准函数上的优化性能比较Ｔａｂｌｅ３㊀Ｏｐｔｉｍｉｚａｔｉｏｎｐｅｒｆｏｒｍａｎｃｅｃｏｍｐａｒｉｓｏｎｏｆｅａｃｈａｌｇｏｒｉｔｈｍｏｎｔｈｅｂｅｎｃｈｍａｒｋｆｕｎｃｔｉｏｎ函数编号指标ＳＣＡＣｈＯＡＷＯＡＡＯＡＧＷＯＬＭＧＷＯＦ１Ｍｅａｎ均值２．８２ˑ１０１５．４５ˑ１０－６２．２０ˑ１０－７２１．５７ˑ１０－７１．８４ˑ１０－２７０Ｓｔｄ标准差７．１５ˑ１０１３．３４ˑ１０－６１．３４ˑ１０－７１４．３６ˑ１０－７２．３５ˑ１０－２８０Ｆ２Ｍｅａｎ均值６．４８ˑ１０－２５．４８ˑ１０－５５．５５ˑ１０－５１４．０８１．０２ˑ１０－１６０Ｓｔｄ标准差３．４５ˑ１０－２５．０２ˑ１０－５９．５４ˑ１０－５１５．１１４．６１ˑ１０－１７０Ｆ３Ｍｅａｎ均值１．２５ˑ１０４６．４５ˑ１０２１．０２ˑ１０４９．６１ˑ１０３５．２１ˑ１０－５０Ｓｔｄ标准差３．１６ˑ１０３８．６４ˑ１０２６．３２ˑ１０４３．２２ˑ１０２１．１７ˑ１０－４０Ｆ４Ｍｅａｎ均值２．７７ˑ１０１９．１５ˑ１０－１４．１１ˑ１０１１．２１１．０４ˑ１０－６０Ｓｔｄ标准差５．６８ˑ１０１５．４７ˑ１０－１２．１９ˑ１０１１．３９１．４７ˑ１０－６０Ｆ５Ｍｅａｎ均值３．２７ˑ１０－２７．６４ˑ１０－３２．４５ˑ１０－３５．１３ˑ１０－１２．３０ˑ１０－３２．４５ˑ１０－５Ｓｔｄ标准差５．９８ˑ１０－２５．１６ˑ１０－３３．０９ˑ１０－３３．１８ˑ１０－２１．７０ˑ１０－３２．０４ˑ１０－５Ｆ６Ｍｅａｎ均值３．０２ˑ１０１８．９９ˑ１０１６．１１ˑ１０－１５４．６７ˑ１０１４．２８０Ｓｔｄ标准差６．４８ˑ１０１１．０２ˑ１０１１．９８ˑ１０－１４２．１３ˑ１０１５．４４０Ｆ７Ｍｅａｎ均值５．５１４．０７ˑ１０１１．１１ˑ１０－１５２．４５ˑ１０－１２．０５ˑ１０－１３８．８８ˑ１０－１６Ｓｔｄ标准差１．８４５．１１ˑ１０－２７．１６ˑ１０－１５４．４１１．１７ˑ１０－１４０Ｆ８Ｍｅａｎ均值３．６５３．４７ˑ１０－２６．３９ˑ１０－２２．５８ˑ１０－２４．６８ˑ１０－３０Ｓｔｄ标准差２．００ˑ１０－１５．１９ˑ１０－２４．７７ˑ１０－２８．１２ˑ１０－２７．５５ˑ１０－３０Ｆ９Ｍｅａｎ均值４．５５ˑ１０－２５．４０ˑ１０－３５．４９ˑ１０－３９４．１１ˑ１０６．７９ˑ１０－４０Ｓｔｄ标准差１．３６ˑ１０－２１．２４ˑ１０－２２．３３ˑ１０－３８２．２８ˑ１０１．１７ˑ１０－４０㊀㊀由表３可以看出，在基准测试中，对于Ｆ１Ｆ４㊁Ｆ６㊁Ｆ８和Ｆ９函数，对比算法均未能找到最优解，而ＬＭＧＷＯ算法达到１００％的求解精度㊂在求解Ｆ５和Ｆ８函数时，ＬＭＧＷＯ的求解精度优于其他５种对比算法，但也与其他算法一样容易陷入局部最优㊂基于以上分析说明ＬＭＧＷＯ算法比其他算法具有更高的求解精度和稳定性，证明了其有效性和优越性㊂３．３㊀ＬＭＧＷＯ算法在高维条件的性能分析为了进一步验证ＬＭＧＷＯ求解高维优化问题的性能，以算法解的均值和平均变化率为评价指标，对９个函数在１００５００维增量下进行测试，将本文提出的ＬＭＧＷＯ算法与原始ＧＷＯ算法独立运行３０次，并记录其均值，实验结果见表４㊂由表４可知，随着维数的增加，ＬＭＧＷＯ的均值基本保持不变，Ｆ１㊁Ｆ２㊁Ｆ３㊁Ｆ４㊁Ｆ６㊁Ｆ９函数的ＬＭＧＷＯ均值保持为０㊂随着维数的增加，ＧＷＯ均值呈现增加趋势㊂在测试函数Ｆ５上，ＬＭＧＷＯ算法的均值基本保持不变，而ＧＷＯ算法的均值变化明显大于ＬＭＧＷＯ算法；在测试函数Ｆ８上，ＬＭＧＷＯ算法的平均变化率均为０，远低于ＧＷＯ算法的平均变化率㊂０５智㊀能㊀计㊀算㊀机㊀与㊀应㊀用㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀第１４卷㊀表４㊀ＬＭＧＷＯ与ＧＷＯ在不同维度下优化函数均值的比较Ｔａｂｌｅ４㊀ＣｏｍｐａｒｉｓｏｎｏｆＬＭＧＷＯａｎｄＧＷＯｏｐｔｉｍｉｚａｔｉｏｎｆｕｎｃｔｉｏｎｍｅａｎｖａｌｕｅｓｉｎｄｉｆｆｅｒｅｎｔｄｉｍｅｎｓｉｏｎｓ函数编号算法名称维数１００２００３００４００５００平均变化率／％Ｆ１ＧＷＯ１．４６ˑ１０－１２１．４３ˑ１０－７５．７９ˑ１０－５８．０８ˑ１０－４１．７９ˑ１０－３４．４８ˑ１０－４ＬＭＧＷＯ００００００Ｆ２ＧＷＯ５．３５ˑ１０－８３．２５ˑ１０－５６．７９ˑ１０－４３．３４ˑ１０－３１．１２ˑ１０－２２．８０ˑ１０－３ＬＭＧＷＯ００００００Ｆ３ＧＷＯ７．３１ˑ１０２２．０２ˑ１０４９．１１ˑ１０４１．９４ˑ１０５３．０９ˑ１０５７．７１ˑ１０４ＬＭＧＷＯ００００００Ｆ４ＧＷＯ８．８２ˑ１０－１２．６１ˑ１０１４．７１ˑ１０１６．０３ˑ１０１６．４８ˑ１０１１．６０ˑ１０１ＬＭＧＷＯ００００００Ｆ５ＧＷＯ７．０３ˑ１０－３１．２６ˑ１０－２３．４９ˑ１０－２６．６３ˑ１０－２９．４６ˑ１０－２２．１９ˑ１０－２ＬＭＧＷＯ３．４１ˑ１０－５３．８７ˑ１０－５４．０５ˑ１０－５４．７２ˑ１０－５６．３９ˑ１０－５７．４５ˑ１０－６Ｆ６ＧＷＯ９．２９２．４２ˑ１０１３．９１ˑ１０１５．０２ˑ１０１７．２０ˑ１０１１．５７ˑ１０１ＬＭＧＷＯ００００００Ｆ７ＧＷＯ６．７７ˑ１０－７２．２２ˑ１０－５５．７４ˑ１０－４９．０９ˑ１０－４２．０２ˑ１０－３５．０５ˑ１０－４ＬＭＧＷＯ８．８８ˑ１０－１６８．８８ˑ１０－１６８．８８ˑ１０－１６８．８８ˑ１０－１６８．８８ˑ１０－１６０Ｆ８ＧＷＯ８．０５ˑ１０－３１．４５ˑ１０－２２．１４ˑ１０－２７．５３ˑ１０－２９．４６ˑ１０－２２．１６ˑ１０－２ＬＭＧＷＯ００００００Ｆ９ＧＷＯ２．８１ˑ１０－３１．１３ˑ１０－２２．５９ˑ１０－２４．５４ˑ１０－２１．６９ˑ１０－１４．１５ˑ１０－２ＬＭＧＷＯ００００００㊀㊀２种算法在不同维度下均值的变化情况如图３所示㊂在９个函数中，ＧＷＯ的均值随着维度变大而显著增加，ＬＭＧＷＯ的均值保持不变㊂这表明维数的不断增加对ＬＭＧＷＯ的寻优能力影响不大，与ＧＷＯ相比寻优性能更加突出，进一步验证了本文所提算法的优越性㊂1.61.41.21.00.80.60.40.20100150200250300350400450500F u n c t i o n d i m e n s i o nA v e r a g e o p t i m i z a t i o n v a l u e /10-3G WO L M G WO(a )F 1变化曲线605040302010100150200250300350400450500F u n c t i o n d i m e n s i o nA v e r a g e o p t i m i z a t i o n v a l u eG WOL M G WO(d )F 4变化曲线2.01.81.61.41.21.00.80.60.40.20100150200250300350400450500F u n c t i o n d i m e n s i o nA v e r a g e o p t i m i z a t i o n v a l u e /10-3G WO L M G WO(g )F 7变化曲线100150200250300350400450500F u n c t i o n d i m e n s i o nA v e r a g e o p t i m i z a t i o n v a l u e G WOL M G WO(h )F 8变化曲线0.090.080.070.060.050.040.030.020.01100150200250300350400450500F u n c t i o n d i m e n s i o nA v e r a g e o p t i m i z a t i o n v a l u eG WOL M G WO(e )F 5变化曲线0.0100.0080.0060.0040.002100150200250300350400450500F u n c t i o n d i m e n s i o nA v e r a g e o p t i m i z a t i o n v a l u eG WOL M G WO(b )F 2变化曲线0.090.080.070.060.050.040.030.020.01100150200250300350400450500F u n c t i o n d i m e n s i o nA v e r a g e o p t i m i z a t i o n v a l u eG WOL M G WO (i )F 9变化曲线0.160.140.120.100.080.060.040.020100150200250300350400450500F u n c t i o n d i m e n s i o nA v e r a g e o p t i m i z a t i o n v a l u eG WOL M G WO (f )F 6变化曲线706050403020100100150200250300350400450500F u n c t i o n d i m e n s i o nA v e r a g e o p t i m i z a t i o n v a l u e /105G WOL M G WO(c )F 3变化曲线3.02.52.01.51.00.5图３㊀基于函数维数变化曲线的函数优化Ｆｉｇ．３㊀Ｆｕｎｃｔｉｏｎｏｐｔｉｍｉｚａｔｉｏｎｂａｓｅｄｏｎｔｈｅｃｕｒｖｅｏｆｆｕｎｃｔｉｏｎｄｉｍｅｎｓｉｏｎｃｈａｎｇｅ１５第３期王恒，等：一种求解最优潮流的改进灰狼优化算法４㊀求解最优潮流（ＯＰＦ）问题为了验证ＬＭＧＷＯ算法的有效性和可行性，在标准ＩＥＥＥ－３０总线测试系统模型上对算法进行了测试㊂该系统包括６台发电机㊁４台变压器㊁９台分流器和４１条支路㊂ＩＥＥＥ３０母线系统单线如图４所示㊂图４中母线１为平衡母线，母线２㊁５㊁８㊁１１㊁１３为电压控制（ＶｏｌｔａｇｅＣｏｎｔｒｏｌ）和无功功率（ＲｅａｃｔｉｖｅＰｏｗｅｒ）母线，其余为有功功率（ＡｃｔｉｖｅＰｏｗｅｒ）和无功功率（ＲｅａｃｔｉｖｅＰｏｗｅｒ）母线㊂本文假设变压器比及无功补偿输出为连续变量，最大迭代次数设置为２００次，种群规模为４０，ＯＰＦ问题维度为２４㊂231314121615181920212210911262524292730286431257817图４㊀ＩＥＥＥ３０总线测试系统单线图Ｆｉｇ．４㊀ＳｉｎｇｌｅｌｉｎｅｄｉａｇｒａｍｏｆＩＥＥＥ３０ｂｕｓｔｅｓｔｓｙｓｔｅｍ４．１㊀案例１：燃料成本（ＦＣ）最小化最小化燃料成本是指通过各种手段和方法，将燃料成本控制在最低水平，以提高经济效益，同时也能够减少对环境的影响㊂将ＬＭＧＷＯ算法与灰狼优化算法（ＧＷＯ）［１４］㊁算术优化算法（ＡＯＡ）［１８］㊁正弦余弦算法（ＳＣＡ）［１９］㊁猩猩优化算法（ＣｈＯＡ）［２０］㊁鲸鱼优化算法（ＷＯＡ）［２１］算法进行对比实验，实验结果见表５㊂由表５可知，优化后的ＬＭＧＷＯ算法燃油成本为７９９．３９４４Ɣ／Ｈ㊂与初始情况相比，燃料成本降低了１１．３７％，具有更加优越的性能㊂表５㊀不同算法在案例１上的比较结果Ｔａｂｌｅ５㊀ＣｏｍｐａｒｉｓｏｎｒｅｓｕｌｔｓｏｆｄｉｆｆｅｒｅｎｔａｌｇｏｒｉｔｈｍｓｉｎＣａｓｅ１算法名称燃油成本／（Ɣ㊃ｈ－１）ＧＷＯ７９９．９６２４ＡＯＡ７９９．９２１７ＳＣＡ８０１．９７００ＣｈＯＡ８００．１８５３ＷＯＡ８００．１０１８ＬＭＧＷＯ７９９．３９４４４．２㊀案例２：有功功率损耗（ＡＰＬ）最小化有功功率损耗（ＡＰＬ）是指电路中有功电流通过负载时所产生的功率损耗㊂有功功率损耗会导致电能转换效率降低，增加能源消耗和运营成本㊂因此，对于电力系统设计和运行来说，减小有功功率损耗是非常重要的㊂将ＬＭＧＷＯ算法与灰狼优化算法（ＧＷＯ）［１４］㊁算术优化算法（ＡＯＡ）［１８］㊁正弦余弦算法（ＳＣＡ）［１９］㊁猩猩优化算法（ＣｈＯＡ）［２０］㊁鲸鱼优化算法（ＷＯＡ）［２１］算法进行对比实验，实验结果见表６㊂根据表６的实验结果，本文提出的ＬＭＧＷＯ算法以有功功率损耗（ＡＰＬ）最小为目标，优于其他用于求解最优潮流（ＯＰＦ）问题的对比算法㊂表６㊀不同算法在案例２上的比较结果Ｔａｂｌｅ６㊀ＣｏｍｐａｒｉｓｏｎｒｅｓｕｌｔｓｏｆｄｉｆｆｅｒｅｎｔａｌｇｏｒｉｔｈｍｓｉｎＣａｓｅ２算法名称有功功率损耗／ＭＷＧＷＯ３．０２６４ＡＯＡ３．１２３２ＳＣＡ３．８２３９ＣｈＯＡ３．１６００ＷＯＡ３．５１６５ＬＭＧＷＯ２．９６９１５㊀结束语本文提出了一种改进的灰狼优化算法（ＬＭＧＷＯ），针对原始ＧＷＯ算法在求解ＯＰＦ问题时的性能进行了２方面的改进㊂将修正反向学习策略与透镜成像学习策略和乘除算子策略相结合，对９个具有不同特性的基准函数进行测试，并与现有元启发式算法进行对比实验㊂实验结果表明，ＬＭＧＷＯ比其他算法具有更好的稳定性和寻优性能㊂在实际应用案例中，将ＬＭＧＷＯ算法和其他对比算法在ＩＥＥＥ３０节点标准测试系统模型上进行对比测试㊂实验结果表明，ＬＭＧＷＯ算法具有较好的性能㊂在未来的工作中，将使用ＬＭＧＷＯ算法解决更困难的最优潮流（ＯＰＦ）问题㊂参考文献［１］ＳＡＬＧＡＤＯＲ，ＢＲＡＭＥＬＬＥＲＡ，ＡＩＴＣＨＩＳＯＮＰ．Ｏｐｔｉｍａｌｐｏｗｅｒｆｌｏｗｓｏｌｕｔｉｏｎｓｕｓｉｎｇｔｈｅｇｒａｄｉｅｎｔｐｒｏｊｅｃｔｉｏｎｍｅｔｈｏｄ．Ｐａｒｔ１：Ｔｈｅｏｒｅｔｉｃａｌｂａｓｉｓ［Ｊ］．ＩＥＴＰｒｏｃｅｅｄｉｎｇｓＣ（Ｇｅｎｅｒａｔｉｏｎ，ＴｒａｎｓｍｉｓｓｉｏｎａｎｄＤｉｓｔｒｉｂｕｔｉｏｎ），１９９０，１３７（６）：４２４－４２８．［２］ＴＩＮＮＥＹＷＦ，ＨＡＲＴＣＥ．ＰｏｗｅｒｆｌｏｗｓｏｌｕｔｉｏｎｂｙＮｅｗｔｏｎᶄｓｍｅｔｈｏｄ［Ｊ］．ＩＥＥＥＴｒａｎｓａｃｔｉｏｎｓｏｎＰｏｗｅｒＡｐｐａｒａｔｕｓａｎｄＳｙｓｔｅｍｓ，１９６７（１１）：１４４９－１４６０．［３］ＬＥＶＩＶＡ，ＮＥＤＩＣＤＰ．Ａｐｐｌｉｃａｔｉｏｎｏｆｔｈｅｏｐｔｉｍａｌｐｏｗｅｒｆｌｏｗｍｏｄｅｌｉｎｐｏｗｅｒｓｙｓｔｅｍｅｄｕｃａｔｉｏｎ［Ｊ］．ＩＥＥＥＴｒａｎｓａｃｔｉｏｎｓｏｎＰｏｗｅｒ２５智㊀能㊀计㊀算㊀机㊀与㊀应㊀用㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀第１４卷㊀Ｓｙｓｔｅｍｓ，２００１，１６（４）：５７２－５８０．［４］ＯＬＯＦＳＳＯＮＭ，ＡＮＤＥＲＳＳＯＮＧ，ＳÖＤＥＲＬ．Ｌｉｎｅａｒｐｒｏｇｒａｍｍｉｎｇｂａｓｅｄｏｐｔｉｍａｌｐｏｗｅｒｆｌｏｗｕｓｉｎｇｓｅｃｏｎｄｏｒｄｅｒｓｅｎｓｉｔｉｖｉｔｉｅｓ［Ｊ］．ＩＥＥＥＴｒａｎｓａｃｔｉｏｎｓｏｎＰｏｗｅｒＳｙｓｔｅｍ，１９９５，１０：１６９１－１６９７．［５］ＤＩＮＧＸｉａｏｙｉｎｇ，ＷＡＮＧＸｉｆａｎ，ＳＯＮＧＹｏｎｇｈｕａ，ｅｔａｌ．Ｔｈｅｉｎｔｅｒｉｏｒｐｏｉｎｔｂｒａｎｃｈａｎｄｃｕｔｍｅｔｈｏｄｆｏｒｏｐｔｉｍａｌｐｏｗｅｒｆｌｏｗ［Ｃ］／／ＰｒｏｃｅｅｄｉｎｇｓｏｆＩｎｔｅｒｎａｔｉｏｎａｌＣｏｎｆｅｒｅｎｃｅｏｎＰｏｗｅｒＳｙｓｔｅｍＴｅｃｈｎｏｌｏｇｙ．Ｋｕｎｍｉｎｇ，Ｃｈｉｎａ：ＩＥＥＥ，２００２，１：６５１－６５５．［６］刘自发，葛少云，余贻鑫．基于混沌粒子群优化方法的电力系统无功最优潮流［Ｊ］．电力系统自动化，２００５，２９（７）：５３－５７．［７］ＦＡＲＨＡＴＭ，ＫＡＭＥＬＳ，ＡＴＡＬＬＡＨＡＭ，ｅｔａｌ．ＥＳＭＡ－ＯＰＦ：Ｅｎｈａｎｃｅｄｓｌｉｍｅｍｏｕｌｄａｌｇｏｒｉｔｈｍｆｏｒｓｏｌｖｉｎｇｏｐｔｉｍａｌｐｏｗｅｒｆｌｏｗｐｒｏｂｌｅｍ［Ｊ］．Ｓｕｓｔａｉｎａｂｉｌｉｔｙ，２０２２，１４（４）：２３０５．［８］ＡｔｔｉａＡＦ，ＥｌＳｅｈｉｅｍｙＲＡ，ＨａｓａｎｉｅｎＨＭ．ＯｐｔｉｍａｌｐｏｗｅｒｆｌｏｗｓｏｌｕｔｉｏｎｉｎｐｏｗｅｒｓｙｓｔｅｍｓｕｓｉｎｇａｎｏｖｅｌＳｉｎｅ－Ｃｏｓｉｎｅａｌｇｏｒｉｔｈｍ［Ｊ］．ＩｎｔｅｒｎａｔｉｏｎａｌＪｏｕｒｎａｌｏｆＥｌｅｃｔｒｉｃａｌＰｏｗｅｒ＆ＥｎｅｒｇｙＳｙｓｔｅｍｓ，２０１８，９９：３３１－３４３．［９］ＷＡＲＩＤＷ．ＯｐｔｉｍａｌｐｏｗｅｒｆｌｏｗｕｓｉｎｇｔｈｅＡＭＴＰＧ－Ｊａｙａａｌｇｏｒｉｔｈｍ［Ｊ］．ＡｐｐｌｉｅｄＳｏｆｔＣｏｍｐｕｔｉｎｇ，２０２０，９１：１０６２５２．［１０］ＷＡＲＩＤＷ，ＨＩＺＡＭＨ，ＭＡＲＩＵＮＮ，ｅｔａｌ．ＯｐｔｉｍａｌｐｏｗｅｒｆｌｏｗｕｓｉｎｇｔｈｅＪａｙａａｌｇｏｒｉｔｈｍ［Ｊ］．Ｅｎｅｒｇｉｅｓ，２０１６，９（９）：６７８．［１１］ＡＢＤＥＳ，ＫＡＭＥＬＳ，ＥＢＥＥＤＭ，ｅｔａｌ．Ａｎｉｍｐｒｏｖｅｄｖｅｒｓｉｏｎｏｆｓａｌｐｓｗａｒｍａｌｇｏｒｉｔｈｍｆｏｒｓｏｌｖｉｎｇｏｐｔｉｍａｌｐｏｗｅｒｆｌｏｗｐｒｏｂｌｅｍ［Ｊ］．ＳｏｆｔＣｏｍｐｕｔｉｎｇ，２０２１，２５：４０２７－４０５２．［１２］ＮＧＵＹＥＮＴＴ．Ａｈｉｇｈｐｅｒｆｏｒｍａｎｃｅｓｏｃｉａｌｓｐｉｄｅｒｏｐｔｉｍｉｚａｔｉｏｎａｌｇｏｒｉｔｈｍｆｏｒｏｐｔｉｍａｌｐｏｗｅｒｆｌｏｗｓｏｌｕｔｉｏｎｗｉｔｈｓｉｎｇｌｅｏｂｊｅｃｔｉｖｅｏｐｔｉｍｉｚａｔｉｏｎ［Ｊ］．Ｅｎｅｒｇｙ，２０１９，１７１：２１８－２４０．［１３］ＡＢＤＥＬ－ＲＡＨＩＭＡＭＭ，ＳＨＡＡＢＡＮＳＡ，ＲＡＧＬＥＮＤＩＪ．Ｏｐｔｉｍａｌｐｏｗｅｒｆｌｏｗｕｓｉｎｇａｔｏｍｓｅａｒｃｈｏｐｔｉｍｉｚａｔｉｏｎ［Ｃ］／／２０１９ＩｎｎｏｖａｔｉｏｎｓｉｎＰｏｗｅｒａｎｄＡｄｖａｎｃｅｄＣｏｍｐｕｔｉｎｇＴｅｃｈｎｏｌｏｇｉｅｓ（ｉ－ＰＡＣＴ）．Ｖｅｌｌｏｒｅ，Ｉｎｄｉａ：ＩＥＥＥ，２０１９，１：１－４．［１４］ＭＩＲＪＡＬＩＬＩＳ，ＭＩＲＪＡＬＩＬＩＳＭ，ＬｅｗｉｓＡ．Ｇｒｅｙｗｏｌｆｏｐｔｉｍｉｚｅｒ［Ｊ］．ＡｄｖａｎｃｅｓｉｎＥｎｇｉｎｅｅｒｉｎｇＳｏｆｔｗａｒｅ，２０１４，６９：４６－６１．［１５］ＮＵＡＥＫＡＥＷＫ，ＡＲＴＲＩＴＰ，ＰＨＯＬＤＥＥＮ，ｅｔａｌ．Ｏｐｔｉｍａｌｒｅａｃｔｉｖｅｐｏｗｅｒｄｉｓｐａｔｃｈｐｒｏｂｌｅｍｕｓｉｎｇａｔｗｏ－ａｒｃｈｉｖｅｍｕｌｔｉ－ｏｂｊｅｃｔｉｖｅｇｒｅｙｗｏｌｆｏｐｔｉｍｉｚｅｒ［Ｊ］．ＥｘｐｅｒｔＳｙｓｔｅｍｓｗｉｔｈＡｐｐｌｉｃａｔｉｏｎｓ，２０１７，８７：７９－８９．［１６］ＰＲＥＣＵＰＲＥ，ＤＡＶＩＤＲＣ，ＰＥＴＲＩＵＥＭ．Ｇｒｅｙｗｏｌｆｏｐｔｉｍｉｚｅｒａｌｇｏｒｉｔｈｍ－ｂａｓｅｄｔｕｎｉｎｇｏｆｆｕｚｚｙｃｏｎｔｒｏｌｓｙｓｔｅｍｓｗｉｔｈｒｅｄｕｃｅｄｐａｒａｍｅｔｒｉｃｓｅｎｓｉｔｉｖｉｔｙ［Ｊ］．ＩＥＥＥＴｒａｎｓａｃｔｉｏｎｓｏｎＩｎｄｕｓｔｒｉａｌＥｌｅｃｔｒｏｎｉｃｓ，２０１７，６４（１）：５２７－５３４．［１７］ＳＡＸＥＮＡＡ，ＫＵＭＡＲＲ，ＤＡＳＳ．β－ｃｈａｏｔｉｃｍａｐｅｎａｂｌｅｄｇｒｅｙｗｏｌｆｏｐｔｉｍｉｚｅｒ［Ｊ］．ＡｐｐｌｉｅｄＳｏｆｔＣｏｍｐｕｔｉｎｇ，２０１９，７５：８４－１０５．［１８］ＡＢＵＡＬＩＧＡＨＬ，ＤＩＡＢＡＴＡ，ＭＩＲＪＡＬＩＬＩＳ，ｅｔａｌ．Ｔｈｅａｒｉｔｈｍｅｔｉｃｏｐｔｉｍｉｚａｔｉｏｎａｌｇｏｒｉｔｈｍ［Ｊ］．ＣｏｍｐｕｔｅｒＭｅｔｈｏｄｓｉｎＡｐｐｌｉｅｄＭｅｃｈａｎｉｃｓａｎｄＥｎｇｉｎｅｅｒｉｎｇ，２０２１，３７６：１１３６０９．［１９］ＭＩＲＪＡＬＩＬＩＳ．ＳＣＡ：Ａｓｉｎｅｃｏｓｉｎｅａｌｇｏｒｉｔｈｍｆｏｒｓｏｌｖｉｎｇｏｐｔｉｍｉｚａｔｉｏｎｐｒｏｂｌｅｍｓ［Ｊ］．Ｋｎｏｗｌｅｄｇｅ－ｂａｓｅｄＳｙｓｔｅｍｓ，２０１６，９６：１２０－１３３．［２０］ＫＨＩＳＨＥＭ，ＭＯＳＡＶＩＭＲ．Ｃｈｉｍｐｏｐｔｉｍｉｚａｔｉｏｎａｌｇｏｒｉｔｈｍ［Ｊ］．ＥｘｐｅｒｔＳｙｓｔｅｍｓｗｉｔｈＡｐｐｌｉｃａｔｉｏｎｓ，２０２０，１４９：１１３３３８．［２１］ＭＩＲＪＡＬＩＬＩＳ，ＬＥＷＩＳＡ．Ｔｈｅｗｈａｌｅｏｐｔｉｍｉｚａｔｉｏｎａｌｇｏｒｉｔｈｍ［Ｊ］．ＡｄｖａｎｃｅｓｉｎＥｎｇｉｎｅｅｒｉｎｇＳｏｆｔｗａｒｅ，２０１６，９５：５１－６７．３５第３期王恒，等：一种求解最优潮流的改进灰狼优化算法。

Optimum

Optimum Gate Design of FreeForm Injection Mould using the Abductive NetworkJ. C. LinDepartment of Mechanical Design Engineering, National Hu-Wei Institute of Technology, Yunlin, TaiwanThis study uses the injection position and size of the gate as the major control parameters for a simulated injection-mould. Once the injection parameters (gate size and gate position) are given, the product performance (deformation) can be accuratelypredicted by the abductive network developed.To avoid the numerous influencing factors, first the part-lineof the parameter equation is created by an abductive networkto limit the range of the gate. The optimal injection parameterscan be searched for by a simulation annealing (SA) optimisation algorithm, with a performance index, to obtain a perfectpart. The major purpose is searching for the optimal gatelocation on the part-surface, and minimising the air-trap anddeformation after part formation. This study also uses a practical example which has been and proved by experiment toachieve a satisfactory result.Keywords: Abductive network; Injection mould; Simulation annealing (SA)1. IntroductionOwing to the rapid development of industry and commerce in recent years, there is a need for rapid and high volume production of goods. The products are manufactured using moulds in order to save the time and cost. Plastic products are the majority. Owing to these products not requiring complicated processes it is possible to cope with market demand speedilyand conveniently.In traditional plastic production, the designs of the portionsof the mould are determined by humans. However, because of the increased performance requirements, the complexity of plastic products has increased. First, the geometric shapes ofthe plastic products are difficult to draw, and the internalshape is often complex which also affects the production ofthe product.Injection processing can be divided into three stages:1. Heat the plastic material to a molten state. Then, by highpressure, force the material to the runner, and then into the mould cavity.2. When the filling of the mould cavity is completed, more molten plastic should be delivered into the cavity at high pressure to compensate for the shrinkage of the plastic. This ensures complete filling of the mould cavity.3. Take out the product after cooling.Though the filling process is only a small proportion of the complete formation cycle, it is very important. If filling in incomplete, there is no pressure holding and cooling is required. Thus, the flow of the plastic fluid should be controlled thoroughly to ensure the quality of the product.The isothermal filling model of a Newtonian fluid is the simplest injection mould-flow filling model. Richardson [1] produced a complete and detailed concept. The major concept is based on the application of lubrication theory, and he simplified the complex 3D flow theory to 2D Hele-Shaw flow. The Hele-Shaw flow was used to simulate the potential flow and was furthermore used in the plasticity flow of the plastic. He assumed the plasticity flow on an extremely thin plate and fully developed the flow by ignoring the speed change through the thickness. Kamal et al. used similar methods to obtain thefilling condition for a rectangular mould cavity, and the analytical result obtained was almost identical to the experimental result.Plastic material follows the Newtonian fluid model for flowin a mould cavity, and Bird et al. [2–4] derived mould-flow theory based on this. When the shape of a mould is complicated and there is variation in thickness, then the equilibrium equations changes to nonlinear and has no analytical solution.Thus, it can be solved only by finite-difference or numerical solutions [2,5]Of course, as the polymer is a visco-elastic fluid, it is bestto solve the flow problem by using visco-elasticity equations.In 1998, Goyal et al. used the White Metzner visco-elasticity model to simulate the disk mould-flow model for central pouring. Metzner, using a finite-difference method to solve the governing equation, fould the visco-elasticity effect would not change the distribution of speed and temperature. However, it affects the stress field very much. If it is a pure visco-elasticflow model, the popular GNF model is generally used toperform numerical simulation.Currently, finite-element methods are mostly used for the solution of mould-flow problems. Other methods are purevisco-elastic models, such as C-FOLW and MOLD-FLOW software. We used this method as well. Some software employs the visco-elastic White–Metzner model, but it is limited to 2D mould-flow analysis. Simple mould flow analysis is limited by CPU time. For the complicated mould shapes, Papthanasionet al. used UCM fluid for filling analysis, using a finitedifference method and BFCC coordination system applicationfor the solution of the more complicated mould-shape andfilling problem, but it was not commercialised [6].Many factors affect plastic material injection. The filling speed, injection pressure and molten temperature, holding pressure [7–12], cooling tube [13,14] and injection gate affect the accuracy of the plastic product, because, when the injection processing is completed, the flow of material in the mouldcavity results in uneven temperature and pressure, and induces residual stress and deformation of the workpiece after cooling.It is difficult to decide on the mould part-surface and gate positions. Generally, the mould part-surface is located at thewidest plane of the mould. Searching for the optimal gate position depends on experience. Minimal modification to the mould is required if you are lucky. However, the time andcost required for the modification of most injection-moulds exceeds the original cost, if the choice of the part-line is poor. For the mould part surface, many workers used various methods to search for the optimal mould part-line, such as geometric shape and feature-based design [15–17]. Some workers used finite-element methods and abductive networks to look for the optimal gate design for a die-casting mould [18].This study used an abductive network to establish the parameter relationship of the mould part-line, and used this formulafor searching for 22 points on the injection mould part-line to serve as the location for an injection gate. Abductive networks are used to match injection pressure and pressure holding time to perform injection formation analysis, and to establish a relationship between these parameters, and the output result of the injection process.It has been shown that prediction accuracy in abductive networks is much higher than that in other networks [19]. Abductive networks based on the abductive modelling techniqueare able to represent complex and uncertain relationships between mould-flow analysis results and injection parameters.It has beeen shown that the injection-strain and injection-stressin a product can be predicted, with reasonable accuracy, basedon the developed networks. The abductive network has been constructed once the relationships of gate location that areinput and simulated have been determined; an appropriate optimisation algorithm with a performance index is then usedto search for the optimal location parameters.In this paper, an optimisation method for simulated annealing [20] is presented. The simulated annealing algorithm is a simulation of the annealing process for minimising the performance index. It has been successfully applied to filtering inimage processing [21], VLSI layout generation [22], discrete tolerance design [23], wire electrical discharge machining [24], deep-draw clearance [25], and casting-die runner design [26],etc. It provides an experimental foundation based on theoryfor the development and application of the technologies.2. Mould-Flow TheoryThe mould flow analysis include four major parts:1. Filling stage.2. Pressure holding stage.3. Cooling and solidification stage.4. Shrinkage and warp, i.e. stress residue stage.Thus, the major mould flow equations are divided into four groups. In the filling stage, the mould cavity is filled with molten plastic fluid at high presssure. Thus, the governing equations include:1. Continuity equation. The plastic deformation or shape change accompany the flow during the filling process (mass conservation):r = plastic density; V = vector velocity2. Momentum equation. Newton’s second law is used to derive the momentum (acceleration condition) or force balance generated by plastic flow:P = flow pressure; f = body force; t = stress tensor.3. Energy equation. The energy conservation of system andlaws of conservation of flow material, if the fluid is incompressible:T = temperature; CP = specific heat of constant pressure; q= heat flux4. Rheology equation=V = deform tensor; (=V)T = transport vector.Holding pressure analysis. The holding pressure process isto maintain the pressure after the mould cavity is filled inorder to inject more plastic, to compensate for the shrinkagein cooling.Cooling analysis. The analysis of the cooling process considers the relationship of the plastic flow distribution and heat transmission. The homogenous mould temperature and the sequence of filling will be affected by the shrinkage of the product formed. If the temperature is distributed non-uniformly,it tends to produce warp. This is mainly due to heat-transferand crystallisation heat of the plastic.r = crystallisation rate; DH = crystallisation heat3. Abductive Network Synthesis and EvaluationMiller [22] observed that human behaviour limits the amountof information considered at a time. The input data are summarised and then the summarised information is passed to a higherreasoning level.In an abductive network, a complex system can be decomposed into smaller, simpler subsystems grouped into severallayers using polynomial function nodes. These nodes evaluate the limited number of inputs by a polynomial function and generate an output to serve as an input to subsequent nodesof the next layer. These polynomial functional nodes are specified as follows:1. NormaliserA normaliser transforms the original input variables into a relatively common region.a1 = q0 + q1x1 (10)Where a1 is the normalised input, q0, q1 are the coefficientsof the normaliser, and x1 is the original input.2. White nodeA white node consists of linear weighted sums of all the outputs of the previous layer.b1 = r0 + r1y1 + r2y2 + r3y3 + % + rnyn (11)Where y1, y2, y3, yn are the input of the previous layer, b1 is the output of the node, and the r0, r1, r2, r3, %, rn are the coefficients of the triple node.3. Single, double, and triple nodesThese names are based on the number of input variables. The algebraic form of each of these nodes is shown in the following:where z1, z2, z3, %, zn, n1, n2, n3, %, nn, o1, o2, o3, %, on are the input of the previous layer, c1, d1, and e1 are the output of the node, and the s0, s1, s2, s3, %, sn, t0, s1, t2, t3, %, tn, u0, u1, u2, u3, %, un are the coefficients of the single, double, and triple nodes.These nodes are third-degree polynomial Eq. and doubles and triples have cross-terms, allowing interaction among the node input variables.4. UnitiserOn the other hand, a unitiser converts the output to a real output.f1 = v0 + v1i1 (15)Where i1 is the output of the network, f1 is the real output,and v0 and v1 are the coefficients of the unitiser.4. Part-Surface ModelThis study uses an actual industrial product as a sample, Fig. 1. The mould part surface is located at the maximum projection area. As shown in Fig. 1, the bottom is the widest plane andis chosen as the mould part surface. However, most important is the searching of gate position on the part surface.This study establishes the parameter equation by using an abductive neuron network, in order to establish the simulated annealing method (SA) to find the optimal gate path position. The parameter equation of a part-surface is expressed byF(Y) = X. First, use a CMM system to measure the XYZ coordinate values (in this study z = 0) of 22 points on the mould part-line on the mould part-surface as illustrated in Table 1, and the gate position is completely on the curve in this space.Prior to developing a space-curve model, a database has tobe trained, and a good relationship msut exist between the control point and abductive network system. A correct andFig. 1. Injection-mould product Table 1. X, Y coordinate.precise curve Eq. is helpful for finding the optimal gate location.To build a complete abductive network, the first requirement is to train the database. The information given by the input and output parameters must be sufficient. A predicted square error (PSE) criterion is then used to determine automatically an optimal structure [23]. The PSE criterion is used to select the least complex but still accurate network.The PSE is composed of two terms:PSE = FSE + KP (16)Where FSE is the average square error of the network forfitting the training data and KP is the complex penalty of the network, shown by the following equation:Where CPM is the complex penalty multiplier, KP is a coefficient of the network; N is the number of training data to beused and s2pis a prior estimate of the model error variance.Based on the development of the database and the predictionof the accuracy of the part-surface, a three-layer abductive network, which comprised design factors (input: various Y coordinate) and output factors (X coordinate) is synthesised automatically. It is capable of predicting accurately the space curve at any point under various control parameters. All polynomial equations used in this network are listed in AppendixA. (PSE = 5.8 ′10-3).Table 2 compares the error predicted by the abductive modeland CMM measurement data. The measurement daa is excluded from the 22 sets of CMM measurement data for establishingthe model. This set of data is used to test the appropriatenessof the model established above. We can see from Table 2 thatthe error is approximately 2.13%, which shows that the modelis suitable for this space curve.Table 2. CMMS-coordinate and neural network predict compared (it is not included in any original 22 sets database).5. Create the Injection-Mould Model Similarly, the relationship is established between input parameters (gate location and gate size) and the output parameter (deformation) during the injection process. To build a complete abductive network, the first requirement is to train the database. The information given by the input and the output data mustbe sufficient. Thus, the training factor (gate location) for abductive network training should be satisfactory for making defect-free products. Figure 2 shows the simulation of FEMmould-flow. Table 3 shows the position of the gate and the maximum deformation of the product obtained from mouldflow analysis.Based on the development of the injection-mould model,three-layer abductive networks, which are comprised of injection- mould conditions and the injection-results (deformation),are synthesised automatically. They are capable of predicting accurately the product strain (the result of injection-moulded product) under various control parameters. All polynomial equations used in this network are listed in Appendix B (PSE= 2.3 ′10-5).Table 4 compares the error predicted by the abductive model and the simulation case. The simulation case is excluded from the 22 sets of simulation cases for establishing the model. This set of data is used to test the appropriateness of the model established above. We can see from Table 4 that the error isFig. 2. The deformation of FEM mould-flow.Table 3. Gate location and the maximum strain relationship.Table 4. Mould-flow simulated and neural network predict compared(it is not included in any original 22 set database).approximately 4.62%, which shows that the model is suitablefor this model requirement6. Simulation Annealing TheoryIn 1983, a theory that was capable of solving the global optimisation problem was developed for the optimised problem. The concept was a powerful optimisation algorithm based onthe annealing of a solid which solved the combinatorial optimisation problem of multiple variables. When the temperature is T and energy E, the thermal equilibrium of the systemis a Boltzman distribution:Z(T) = normalisation factor; KB = Boltzman constant;Exp(-E/KBT) = Boltzman factorMetropolis [24] proposed a criterion for simulating the coolingof a solid to a new state of energy balance. The basiccriterion used by Metropolis is an optimisation algorithm called “simulated annealing”. The algorithm was developed by Kirkpatricket al. [20].In this paper, the simulation-annealing algorithm is used tosearch for the optimal control parameters for gate location.Figure 3 shows the flowchart of the simulated annealing search. First, the algorithm is given an initial temperature Ts and afinal temperature Te, and a set of initial process vectors Ox.The objective function obj is defined, based on the injection parameter performance index. The objective function can be recalculated for all the different perturbed compensation parameters. If the new objective function becomes smaller, thepeturbed process parameters are accepted as the new process parameters and the temperature drops a little in scale. That is:where i is the index for the temperature decrement and the CTis the decay ratio for the temperature (CT , 1).However, if the objective function becomes larger, the probability of acceptance of the perturbed process parameters isgiven as:Where KB is the Boltzman constant and D obj is the differentin the objective function. The procedure is repeated until the temperature Ti approaches zero. It shows the energy droppingto the lowest state.Once the model of the relationship among the functions ofthe gate location, the input parameters and output parametersare established, this model can be used to find the optimal parameter for the gate location. The optimal parameter for the process can be obtained by using the objective function to serveFig. 3. Flowchart of the simulated annealing searching.as a starting point. The objective function obj is formulated as follows:Where w1, w2, and w3 are the weight functions.The control parametric of the X, Y location should complywith the part-surface parameter equation. That means the basic condition of optimisation should fall within a certain range: 1. The X-coordinate value obtained from optimisation should be larger than the minimum X-coordinate value, and smaller than the maximum X-coordinate value.2. The Y-coordinate value obtained from optimisation should be larger than the minimum Y-coordinate value, and smaller than the maximum X-coordinate value. The X, Y-coordinate depends on the Appendix A neural network equation.3. The gate-wide size obtained from optimisation should be larger than the minimum size of gate-width, and smallerthan the maximum gate-width size.4. The gate-height size obtained from optimisation should be larger than the minimum size of gate-height, and smaller than the maximum gate-height size.The inequalities are given as follows:The upper-bound conditions should be kept to an acceptable level so as to find the optimal (acccurate) gate coordinate.7. Results and DiscussionAn example of the simulation is used to illustrate the processof optimising the gate location. In the first case, the weight funtion w1 = w2 = w3 = w4 = 1 (the gate location and size is equal weight). In this case, the parameters used in the simulation annealing algorithm are given as: initial temperature Ts= 100°C, final temperature Te = 0.0001°C, decay ratio CT = 0.98, Boltzman constant Ks = 0.00667, and the upper boundof the X-coordinate is -0.0276 mm, upper bound of Y-coordinate is 32.49 mm (X depends on the Y-coordinate).Simulated annealing is used for finding the optimal gate coordinates and gate size, as shown in Table 5. When gatesize is w = 1.8 mm and h = 0.9 mm, the optimal gate locationis Y = 30.0 mm, and X = -7.824 mm. Predicting the minimum deformation using the optimal parameters is about 0.2563 mm(minimum, Table 5).In Fig. 4, for a fixed gate size w = 1.8 mm and the h =0.9 mm, the optimal value is found by using the Y-coordinate variation. In Fig. 4, the Y-coordinate has the minimum deformation when the optimal parameter of the Y-coordinate is30.0 mm (X = -7.824 mm), and the deformation is 0.2563mm (minimum).Table 5 shows a comparison between the simulation mouldflow error and the optimal value predicted by the model. Theerror is approximately 5.4%. In the foregoing discussion, itTable 5. SA optimal parameters and FEM prediction compared (it is not included in any original 22 set database).Fig. 4. Deformation and Y-coordinate relationship.has been shown clearly that the process parameters for optimum bend performance can be systematically obtained through this approach.8. ConclusionThis paper describes a neural network approach for modelling and optimisation of injection-mould gate parameters.1. Comparing the value of errors using the finite-element method and abductive-network prediction we have designedthe injection-mould model. Based on the best modellingusing an abductive network, the complicated relationships between the injection-gate coordinate parameters and deformationperformance can be obtained.2. A global optimisation algorithm, called simulated annealing (SA), is then applied to the abductive network to obtain the optimal process parameters based on an objective function.3. Lastly, a comparison is made between the FEM simulation mould-flow error values and model predicted values by the optimisation process. It shows that the model fits the FEM simulation mould-flow data and the finite-element and abductive-network predictions. Therefore, the rapidity and efficiency of determining optimal designing prameters for injection-moulding can be used successfully to improve the accuracy of the injection-mould design process.。

Piezoelectric energy harvesting device optimization by synchronous electric charge extraction

of this, the design of vibration-powered electronic devices is made possible by the recent development of ultralow-power integrated circuits, able for instance, to feed sensors, to carry out fast signal processing, and to radio-transmit data (Elvin et al., 2001; Ayers et al., 2003; Roundy et al., 2003a, b; Amirtharajah and Chandrakasan, 2004).
Various techniques for the construction of such piezoelectric microgenerators have been investigated. Mechanical parts constitute an important aspect in design, but previous studies have shown that electrical energy conditioning system has also a determinant effect on the energy conversion effectiveness (Ottman et al., 2002, 2003).
Vibrations take place in most of the structures, but their potential energy often remains very low. In spite

百洁高效清洁牙科附件使用评估指南说明书

POWERUse this assessment guide weekly to ensure top clinical efficiency of your magneto inserts and piezo tips.MAGNETO INSERTS & PIEZO TIPSPERFORMANCE ASSESSMENT GUIDEEVALUATE YOUR INSERTS/TIPS EACH WEEKINSERT/TIP COST PER USE CALCULATORWe’ve provided a calculatorto help you estimate your Calculate the Number of Times an Insert/Tip is Used per YearWhat are signs of o-ring wear?• D ry, brittle or cracked appearance• L eaking at the junction of the handpiece and insert• I nsert doesn’t seem secure in the handpieceCan I use oil, petroleum jelly or other lubricants on the o-ring?Lubricate the o-ring with water only. Oil based products such as petroleum jelly can cause drying/cracking of the o-ring and result in unwanted foreign matter within the water channel.How do I change the o-ring?• U sing a pair ofdressing pliers pull the worn o-ring off over the metal stack.• T o replace the o-ring: Slide a wide drinking straw over the metal stack. Apply a new o-ring to the bottom of the straw and roll the o-ring up until it reaches the o-ring slot.O-RINGSINSERT STACKSWhat should the stack look like?The leaves of the stack should uniformly nest together. Check to make sure stack leaves are not: • Splayed • Bent• Fractured• DetachedIn these cases you may experienceinadequate power — or possibly the tip will not vibrate at all.How do I keep the stack in good condition?• H andle the insert by the grip rather than the stack• U se cassettes during cleaning andsterilization as these protect the insert from the weight of instruments that may be stacked on topWhat should I do if my insert stacks are bent, broken or detached?• I nserts with compromised stacks should be taken out of use and replaced.MAGNETO INSERTS| 4PIEZO TIPSPiezo tips are available in two threadpattern styles.S Series• E Series •Piezo handpieces are compatible with eitherE or S series tips so it is critical to know whichseries tips work with your unit.• E-Series tips are compatible with the Piezon150/250, Piezon® Master 700, AIR-FLOWand the AIR-FLOW Master Piezon• S-Series tips are compatible with SymmetryIQ, Acteon and NSK devices*Differences between E andS Series tips• E-Series tips havefour flatsWRENCHESALL MAGNETO INSERTS AND PIEZO TIPS SHOULD BE CLEANED AND STERILIZED BEFORE EACH USE.BEST PRACTICESOptimal infection control and prevention will be facilitated by using tools that allow for improved efficiency and safety. Instrument cassettes will allow you to move inserts/tips from chairside through cleaning without touching or damaging the instruments. The inserts/tips also receive superior protection from the weight of other instruments while in the protective cassette. Piezo tips should remain within the wrench during the cleaning and care cycle.If cassettes are not used, carefully place the insert/tip in the ultrasonic cleaner either alone or in a location where it will not be subject to the weight of other dental instruments.1. CLEANINGTo clean, rinse thoroughly or fully immerse in a mild ultrasonic cleaning solution (pH 7–10.5) like Hu-Friedy’s IMS® Enzymax® or Enzymax Earth™cleaning solutions. When using the IMS Cassette System, the ultrasonic cleaning time has to be at least 16 minutes, unless a longer exposure time is required by the manufacturer of the detergent. Thoroughly rinse the insert with tap water and dry completely prior to sterilization.2. RINSINGThoroughly rinse with tap water and dry completely prior to sterilization3. STERILIZATIONSterilize inserts in an instrument cassette or paper pouch in a steam autoclave. Steam sterilize for at least 4 minutes at 270°F/132°C or 30 minutes at 250°F/121°C. Do not heat above 275°F/135°C.4. DRYINGThirty minutes of drying time is recommended after steam sterilization, prior to use.CLEANING & CARE | 6What effect will repeated instrument processing have on the lifespan of my inserts/tips?Repeated instrument processing will have little to no effect on the lifespan of inserts/ tips. The lifespan of insert/tips is primarily determined by tip wear.Can I put ultrasonic inserts in an ultrasonic cleaner?Yes, ultrasonic inserts may be cleaned in an ultrasonic cleaning unit with a mild ultrasonic solution. The mild ultrasonic solution willnot harm the insert nor will the cavitation produced by the ultrasonic cleaning unit. Recommended ultrasonic solutions are Enzymax and Enzymax Earth.Can I wipe down or spray the inserts/ tips with a surface disinfectant?Do not use chemical or surface disinfectants as deterioration of the material may occur. Avoid contact with corrosive chemicals such as chlorides, sulfates and surface disinfectants. Do not expose inserts to solutions containing phenols. Chemical vapor sterilization, rapid heat transfer sterilization and cold liquid/chemical sterilants are not recommended and may void the warranty. Can I place inserts/tips in a washer/ disinfector?Please refer to the washer/disinfector instructions and closely follow manufacturer recommendations for magneto inserts and piezo tips.For more information onInstrument Reprocessing, pleasevisit /education.faqs7 |MAGNETO INSERTSISSUE POSSIBLE CAUSE SOLUTIONLow vibration Warped, bent, separated or broken stackleaves Replace insert as stack leaves should be straight and intact for optimal vibrationWorn tips Check tips on wear guide, replace when indicatedWater leak at junction of handpiece and insert Worn o-rings Weekly check o-rings and replace when cracked/brittle Worn handpiece Replace handpieceWater flows excessively after releasing foot pedal Air trapped in water line Remove insert. Flush line for at least 30 seconds. Holdhandpiece vertically and fill to the brim with water whenplacing insertGenerator issues Check generator water line including PSI (refer togenerator instruction manual)Low/no water flow Water line blocked from sediment or highlymineralized water Weekly check external water filter for build-up. Replace when discolored or build-up is visible.Insert clogged from water sediment Replace insertInsert pop-out Worn o-rings Weekly check o-rings and replace when cracked/brittleGenerator water pressure (PSI) too high Refer to generator manual for optimal PSI water settings Handpiece: Heat build-up Generator power setting too high Refer to insert package instructions for appropriatepower setting rangeAir trapped in handpiece Remove insert. Flush line for at least 30 seconds. Holdhandpiece vertically and fill to the brim with water whenplacing insertClogged generator water line filter Weekly check external water filter for build-up. Replacewhen discolored or build-up is visible TROUBLESHOOTING| 8PIEZO TIPSISSUE POSSIBLE CAUSE SOLUTIONLittle to no vibration Tip not firmly attached to handpiece Use wrench to tighten tip to handpieceWorn tips Check tips on wear guide, replace when indicated Tip does not stay in wrench Slot in wrench may be worn Discard wrench; reorder new wrenchTip does not stay in handpiece “Flats” on base of tip may be worn Replace tips Handpiece threads may be worn Replace handpieceTip won’t attach to handpiece Wrong “series” of tip Use E or S series tips according to generatormanufacturer specificationsNo/low water flow Water lines blocked from sediment or highlymineralized water Weekly check external water filter for build-up. Replace when discolored or build-up is visibleTip clogged from water sediment Replace tipWorn tips can significantly diminish your scaling efficiency!• 1 millimeter of tip loss results in a25% loss of efficiency• 2 millimeters of tip loss results in a 50% loss of efficiencyWhat is the clinical result of worn inserts/tips• A pply more pressure on the insert, which can dampen the vibration• T urn up the power on the generator, often to high levels for which the tip may not be designedBoth actions can result in inferior tip performance, more time to remove calculus and stain and more discomfort for you and your patient. Tip extends beyond blue line OPTIMUM EFFICIENCY (Scale Away)Tip touches blue line25%EFFICIENCY LOSS (Reorder)Tip touches red line50%EFFICIENCY LOSS (Discard)WORN OR BENT TIPSUse a Wear Guide to checkyour tips weekly!Maximize the efficiencyof your ultrasonic tipsby regularly using awear guide.• F irst, evaluate the shapeof the tip against the wearguide. If the shape doesnot match the silhouette,discard the insert.• S econd, match the tip tothe corresponding outlineto measure its length anddiscard, reorder or scaleaway depending on the tiplength.INSERT AND TIP WEARGood Tip Worn Tip Good Tip Worn TipGood Tip Worn Tip Good Tip Worn Tip | 10FACTORS AFFECTING THE LIFESPAN OF MAGNETO INSERTS AND PIEZO TIPS include usage frequency and amount/type of hard deposit per patient. An insert/tip used once per day on light debris will likely lastlonger than one used several times per day on heavy debris. In addition, the lifespan of an insert or tip can be shortened by many factors including:• incorrect power settings• applying too much pressure• n ot matching the tip to the application• i mproper sterilization and maintenance ULTRASONIC INSERT & TIPTRADE-IN PROGRAMIT’S EASY AS1-2-3!*Receive one ultrasonic insert or tipfree when you purchase two ultrasonicinserts or tips and trade in any three!MAGNETOINSERTWEARGUIDEPIEZO TIP WEAR GUIDE。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。