On the Implementation of a Primal-Dual Interior Point Filter Line Search Algorithm for Large-Scale

内点法介绍（Interior Point Method）在面对无约束的优化命题时，我们可以采用牛顿法等方法来求解。

而面对有约束的命题时，我们往往需要更高级的算法。

单纯形法（Simplex Method）可以用来求解带约束的线性规划命题（LP），与之类似的有效集法（Active Set Method）可以用来求解带约束的二次规划（QP），而内点法（Interior Point Method）则是另一种用于求解带约束的优化命题的方法。

而且无论是面对LP还是QP，内点法都显示出了相当的极好的性能，例如多项式的算法复杂度。

本文主要介绍两种内点法，障碍函数法（Barrier Method）和原始对偶法（Primal-Dual Method）。

其中障碍函数法的内容主要来源于Stephen Boyd与Lieven Vandenberghe的Convex Optimization一书，原始对偶法的内容主要来源于Jorge Nocedal和Stephen J. Wright的Numerical Optimization一书（第二版）。

为了便于与原书对照理解，后面的命题与公式分别采用了对应书中的记法，并且两者方法针对的是不同的命题。

两种方法中的同一变量可能在不同的方法中有不同的意义，如μ。

在介绍玩两种方法后会有一些比较。

障碍函数法Barrier MethodCentral Path举例原始对偶内点法Primal Dual Interior Point Method Central Path举例几个问题障碍函数法（Barrier Method）对于障碍函数法，我们考虑一个一般性的优化命题：minsubject tof0(x)fi(x)≤0,i=1,...,mAx=b(1) 这里f0,...,fm:Rn→R 是二阶可导的凸函数。

同时我也要求命题是有解的，即最优解x 存在，且其对应的目标函数为p。

此外，我们还假设原命题是可行的（feasible）。

A Stable Primal-Dual Approach for Linear ProgrammingMaria Gonzalez-Lima∗Hua Wei†Henry Wolkowicz‡March10,2007University of WaterlooDepartment of Combinatorics&OptimizationWaterloo,Ontario N2L3G1,CanadaResearch Report CORR2004-26COAP1172-04Key words:Linear Programming,large sparse problems,preconditioned conjugate gradi-ents,stability.AbstractThis paper studies a primal-dual interior/exterior-point path-following approach for linear programming that is motivated on using an iterative solver rather than a direct solver forthe search direction.We begin with the usual perturbed primal-dual optimality equationsFµ(x,y,z)=0.Under nondegeneracy assumptions,this nonlinear system is well-posed,i.e.it has a nonsingular Jacobian at optimality and is not necessarily ill-conditioned asthe iterates approach optimality.We use a simple preprocessing step to eliminate boththe primal and dual feasibility equations.This results in a single bilinear equation thatmaintains the well-posedness property.We then apply both a direct solution techniqueas well as a preconditioned conjugate gradient method(PCG),within an inexact Newtonframework,directly on the linearized equations.This is done without forming the usualnormal equations,NEQ,or augmented system.Sparsity is maintained.The work of aniteration for the PCG approach consists almost entirely in the(approximate)solution of thiswell-posed linearized system.Therefore,improvements depend on efficient preconditioning.。

a rXi v :c s /04744v1[cs.A I]16J ul24Reduced cost-based ranking for generating promising subproblems ano 1and W.J.van Hoeve 21DEIS,University of Bologna,Viale Risorgimento 2,40136Bologna,Italy mmilano@deis.unibo.it http://www-lia.deis.unibo.it/Staff/MichelaMilano/2CWI,P.O.Box 94079,1090GB Amsterdam,The Netherlands w.j.van.hoeve@cwi.nl http://www.cwi.nl/~wjvh/Abstract.In this paper,we propose an effective search procedure that interleaves two steps:subproblem generation and subproblem solution.We mainly focus on the first part.It consists of a variable domain value ranking based on reduced costs.Exploiting the ranking,we generate,in a Limited Discrepancy Search tree,the most promising subproblems first.An interesting result is that reduced costs provide a very precise ranking that allows to almost always find the optimal solution in the first generated subproblem,even if its dimension is significantly smaller than that of the original problem.Concerning the proof of optimality,we exploit a way to increase the lower bound for subproblems at higher discrepancies.We show experimental results on the TSP and its time constrained variant to show the effectiveness of the proposed approach,but the technique could be generalized for other problems.1Introduction In recent years,combinatorial optimization problems have been tackled with hybrid methods and/or hybrid solvers [11,18,13,19].The use of problem relax-ations,decomposition,cutting planes generation techniques in a Constraint Pro-gramming (CP)framework are only some examples.Many hybrid approaches are based on the use of a relaxation R ,i.e.an easier problem derived from the originalone by removing (or relaxing)some constraints.Solving R to optimality provides a bound on the original problem.Moreover,when the relaxation is a linear prob-lem,we can derive reduced costs through dual variables often with no additional computational cost.Reduced costs provide an optimistic esteem (a bound)of each variable-value assignment cost.These results have been successfully used for pruning the search space and for guiding the search toward promising re-gions (see [9])in many applications like TSP [12],TSPTW [10],scheduling with sequence dependent setup times [8]and multimedia applications [4].We propose here a solution method,depicted in Figure 1,based on a two step search procedure that interleaves (i )subproblem generation and (ii )subproblem solution.In detail,we solve a relaxation of the problem at the root node and we2use reduced costs to rank domain values;then we partition the domain of eachvariable X i in two sets,i.e.,the good part D goodi and the bad part D badi.We searchthe tree generated by using a strategy imposing on the left branch the branchingconstraint X i∈D goodi while on the right branch we impose X i∈D badi.At eachleaf of the subproblem generation tree,we have a subproblem which can now be solved(in the subproblem solution tree).Exploring with a Limited Discrepancy Strategy the resulting search space, we obtain that thefirst generated subproblems are supposed to be the most promising and are likely to contain the optimal solution.In fact,if the ranking criterion is effective(as the experimental results will show),thefirst generated subproblem(discrepancy equal to0)P(0),where all variables range on the good domain part,is likely to contain the optimal solution.The following generatedsubproblems(discrepancy equal to1)P(1)i have all variables but the i-th rangingon the good domain and are likely to contain worse solutions with respect to P(0),but still good.Clearly,subproblems at higher discrepancies are supposed to contain the worst solutions.A surprising aspect of this method is that even by using low cardinality good sets,we almost alwaysfind the optimal solution in thefirst gener-ated subproblem.Thus,reduced costs provide extremely useful information indicating for each variable which values are the most promising.Moreover,this property of reduced costs is independent of the tightness of the relaxation.Tight relaxations are essential for the proof of optimality,but not for the quality of reduced costs.Solving only thefirst subproblem,we obtain a very effective in-complete method thatfinds the optimal solution in almost all test instances.To be complete,the method should solve all subproblems for all discrepancies to prove optimality.Clearly,even if each subproblem could be efficiently solved,if all of them should be considered,the proposed approach would not be applicable. The idea is that by generating the optimal solution soon and tightening the lower bound with considerations based on the discrepancies shown in the paper,we do not have to explore all subproblems,but we can prune many of them.In this paper,we have considered as an example the Travelling Salesman Problem and its time constrained variant,but the technique could be applied to a large family of problems.The contribution of this paper is twofold:(i)we show that reduced costs provide an extremely precise indication for generating promising subproblems, and(ii)we show that LDS can be used to effectively order the subproblems. In addition,the use of discrepancies enables to tighten the problem bounds for each subproblem.The paper is organized as follows:in Section2we give preliminaries on Lim-ited Discrepancy Search(LDS),on the TSP and its time constrained variant. In Section3we describe the proposed method in detail.Section4discusses the implementation,focussing mainly on the generation of the subproblems using LDS.The quality of the reduced cost-based ranking is considered in Section5. In this section also the size of subproblems is tuned.Section6presents the computational results.Conclusion and future work follow.3Fig.1.The structure of the search tree.2Preliminaries2.1Limited Discrepancy SearchLimited Discrepancy Search(LDS)wasfirst introduced by Harvey and Ginsberg[15].The idea is that one can oftenfind the optimal solution by exploring onlya small fraction of the space by relying on tuned(often problem dependent) heuristics.However,a perfect heuristic is not always available.LDS addresses the problem of what to do when the heuristic fails.Thus,at each node of the search tree,the heuristic is supposed to provide the good choice(corresponding to the leftmost branch)among possible alternative branches.Any other choice would be bad and is called a discrepancy.In LDS, one tries tofindfirst the solution with as few discrepancies as possible.In fact,a perfect heuristic would provide us the optimal solution immediately.Since this is not often the case,we have to increase the number of discrepancies so as to make it possible tofind the optimal solution after correcting the mistakes made by the heuristic.However,the goal is to use only few discrepancies since in general good solutions are provided soon.LDS builds a search tree in the following way:thefirst solution explored is that suggested by the heuristic.Then solutions that follow the heuristic for every variable but one are explored:these solutions are that of discrepancy equal to one.Then,solutions at discrepancy equal to two are explored and so on.It has been shown that this search strategy achieves a significant cutoffof the total number of nodes with respect to a depthfirst search with chronological backtracking and iterative sampling[20].42.2TSP and TSPTWLet G=(V,A)be a digraph,where V={1,...,n}is the vertex set and A= {(i,j):i,j∈V}the arc set,and let c ij≥0be the cost associated with arc (i,j)∈A(with c ii=+∞for each i∈V).A Hamiltonian Circuit(tour)of G is a partial digraph¯G=(V,¯A)of G such that:|¯A|=n and for each pair of distinct vertices v1,v2∈V,both paths from v1to v2and from v2to v1exist in¯G(i.e. digraph¯G is strongly connected).The Travelling Salesman Problem(TSP)looks for a Hamiltonian circuit G∗= (V,A∗)whose cost (i,j)∈A∗c ij is a minimum.A classic Integer Linear Programming formulation for TSP is as follows:v(T SP)=min i∈V j∈V c ij x ij(1)subject to i∈V x ij=1,j∈V(2)j∈V x ij=1,i∈V(3)i∈Sj∈V\Sx ij≥1,S⊂V,S=∅(4)x ij integer,i,j∈V(5) where x ij=1if and only if arc(i,j)is part of the solution.Constraints(2) and(3)impose in-degree and out-degree of each vertex equal to one,whereas constraints(4)impose strong connectivity.Constraint Programming relies in general on a different model where we have a domain variable Next i(resp.Prev i)that identifies cities visited after(resp. before)node i.Domain variable Cost i identifies the cost to be paid to go from node i to node Next i.Clearly,we need a mapping between the CP model and the ILP model: Next i=j⇔x ij=1.The domain of variable Next i will be denoted as D i. Initially,D i={1,...,n}.The Travelling Salesman Problem with Time Windows(TSPTW)is a time constrained variant of the TSP where the service at a node i should begin within a time window[a i,b i]associated to the node.Early arrivals are allowed,in the sense that the vehicle can arrive before the time window lower bound.However, in this case the vehicle has to wait until the node is ready for the beginning of service.As concerns the CP model for the TSPTW,we add to the TSP model a domain variable Start i which identifies the time at which the service begins at node i.A well known relaxation of the TSP and TSPTW obtained by eliminating from the TSP model constraints(4)and time windows constraints is the Linear Assignment Problem(AP)(see[3]for a survey).AP is the graph theory problem5 offinding a set of disjoint subtours such that all the vertices in V are visited and the overall cost is a minimum.When the digraph is complete,as in our case,AP always has an optimal integer solution,and,if such solution is composed by a single tour,is then optimal for TSP satisfying constraints(4).The information provided by the AP relaxation is a lower bound LB for the original problem and the reduced cost matrix¯c.At each node of the decision tree,each¯c ij estimates the additional cost to pay to put arc(i,j)in the so-lution.More formally,a valid lower bound for the problem where x ij=1isLB|xij =1=LB+¯c ij.It is well-known that when the AP optimal solution isobtained through a primal-dual algorithm,as in our case(we use a C++adap-tation of the AP code described in[2]),the reduced cost values are obtained without extra computational effort during the AP solution.The solution of the AP relaxation at the root node requires in the worst case O(n3),whereas each following AP solution can be efficiently computed in O(n2)time through a single augmenting path step(see[2]for details).However,the AP does not provide a tight bound neither for the TSP nor for the TSPTW.Therefore we will improve the relaxation in Section3.1.3The Proposed MethodIn this section we describe the method proposed in this paper.It is based on two interleaved steps:subproblem generation and subproblem solution.Thefirst step is based on the optimal solution of a(possibly tight)relaxation of the original problem.The relaxation provides a lower bound for the original problem and the reduced cost matrix.Reduced costs are used for ranking(the lower the better)variable domain values.Each domain is now partitioned according to this ranking in two sets called the good set and the bad set.The cardinality of the good set is problem dependent and is experimentally defined.However,it should be significantly lower than the dimension of the original domains.Exploiting this ranking,the search proceeds by choosing at each node the branching constraint that imposes the variable to range on the good domain, while on backtracking we impose the variable to range on the bad domain.By exploring the resulting search tree by using an LDS strategy we generatefirst the most promising problems,i.e.,those where no or few variables range on the bad sets.Each time we generate a subproblem,the second step starts for optimally solving it.Experimental results will show that,surprisingly,even if the sub-problems are small,thefirst generated subproblem almost always contains the optimal solution.The proof of optimality should then proceed by solving the re-maining problems.Therefore,a tight initial lower bound is essential.Moreover, by using some considerations on discrepancies,we can increase the bound and prove optimality fast.The idea of ranking domain values has been previously used in incomplete algorithms,like GRASP[5].The idea is to produce for each variable the so called Restricted Candidate List(RCL),and explore the subproblem generated only by6RCLs for each variable.This method provides in general a good starting point for performing local search.Our ranking method could in principle be applied to GRASP-like algorithms.Another connection can be made with iterative broadening[14],where one can view the breadth cutoffas corresponding to the cardinality of our good sets. Thefirst generated subproblem of both approaches is then the same.However, iterative broadening behaves differently on backtracking(it gradually restarts increasing the breadth cutoff).3.1Linear relaxationIn Section2.2we presented a relaxation,the Linear Assignment Problem(AP), for both the TSP and TSPTW.This relaxation is indeed not very tight and does not provide a good lower bound.We can improve it by adding cutting planes. Many different kinds of cutting planes for these problems have been proposed and the corresponding separation procedure has been defined[22].In this paper, we used the Sub-tour Elimination Cuts(SECs)for the TSP.However,adding linear inequalities to the AP formulation changes the structure of the relaxation which is no longer an AP.On the other hand,we are interested in maintaining this structure since we have a polynomial and incremental algorithm that solves the problem.Therefore,as done in[10],we relax cuts in a Lagrangean way,thus maintaining an AP structure.The resulting relaxation,we call it AP cuts,still has an AP structure,but provides a tighter bound than the initial AP.More precisely,it provides the same objective function value as the linear relaxation where all cuts are added defining the sub-tour polytope.In many cases,in particular for TSPTW instances,the bound is extremely close to the optimal solution.3.2Domain partitioningAs described in Section2.2,the solution of an Assignment Problem provides the reduced cost matrix with no additional computational cost.We recall that the reduced costc ij.Consequently,D badi =D i\D goodi.The ratio defines the size of the good domains,and will be discussed in Section5.Note that the optimal ratio should be experimentally tuned,in order to obtain the optimal solution in thefirst subproblem,and it is strongly problem7 dependent.In particular,it depends on the structure of the problem we are solving.For instance,for the pure TSP instances considered in this paper,a good ratio is0.05or0.075,while for TSPTW the best ratio observed is around 0.15.With this ratio,the optimal solution of the original problem is indeed located in thefirst generated subproblem in almost all test instances.3.3LDS for generating subproblemsIn the previous section,we described how to partition variable domains in a good and a bad set by exploiting information on reduced costs.Now,we show how to explore a search tree on the basis of this domain partitioning.At each node corresponding to the choice of variable X i,whose domain has been parti-tioned in D goodi and D badi,we impose on the left branch the branching constraintX i∈D goodi ,and on the right branch X i∈D badi.Exploring with a Limited Dis-crepancy Search strategy this tree,wefirst explore the subproblem suggested by the heuristic where all variable range on the good set;then subproblems where all variables but one range on the good set,and so on.If the reduced cost-based ranking criterion is accurate,as the experimental results confirm,we are likely tofind the optimal solution in the subproblem P(0) generated by imposing all variables ranging on the good set of values.If this heuristic fails once,we are likely tofind the optimal solution in one of the nsubproblems(P(1)i with i∈{1,...,n})generated by imposing all variables butone(variable i)ranging on the good sets and one ranging on the bad set.Then,we go on generating n(n−1)/2problems P(2)ij all variables but two(namely iand j)ranging on the good set and two ranging on the bad set are considered, and so on.In Section4we will see an implementation of this search strategy that in a sense squeezes the subproblem generation tree shown in Figure1into a con-straint.3.4Proof of optimalityIf we are simply interested in a good solution,without proving optimality,we can stop our method after the solution of thefirst generated subproblem.In this case,the proposed approach is extremely effective since we almost alwaysfind the optimal solution in that subproblem.Otherwise,if we are interested in a provably optimal solution,we have to prove optimality by solving all sub-problems at increasing discrepancies.Clearly, even if all subproblems could be efficiently solved,generating and solving all of them would not be practical.However,if we exploit a tight initial lower bound, as explained in Section3.1,which is successively improved with considerations on the discrepancy,we can stop the generation of subproblems after few trials since we prove optimality fast.An important part of the proof of optimality is the management of lower and upper bounds.The upper bound is decreased as wefind better solutions,and8the lower bound is increased as a consequence of discrepancy increase.The idea is tofind an optimal solution in thefirst subproblem,providing the best possible upper bound.The ideal case is that all subproblems but thefirst can be pruned since they have a lower bound higher than the current upper bound,in which case we only need to consider a single subproblem.The initial lower bound LB0provided by the Assignment Problem AP cuts at the root node can be improved each time we switch to a higher discrepancy k. For i∈{1,...,n},letc∗i=min j∈D badic∗i for all problems at discrepancy greater than or equal to1.We can increase this bound:let L be the nondecreasing ordered list ofc.Now we define a second relaxation of P having cost matrix9 The search for such a set S and the corresponding domain assignments have been‘squeezed’into a constraint,the discrepancy constraint discrgood is the array containing D goodi for all i∈{1,...,n}.The command solvesubproblem is shorthand for solving the subproblem which has been considered in Section3.5.A more traditional implementation of LDS(referred to as‘standard’in Ta-ble1)exploits tree search,where at each node the domain of a variable is split into the good set or the bad set,as described in Section3.3.In Table1,the performance of this traditional approach is compared with the performance of the discrepancy constraint(referred to as discr10standard discr instance time fails0.08440.15570.05140.10690.09420.06300.12710.201520.0660.075csttime failsrbg021.20.1366 rbg021.30.20158 rbg021.40.0940 rbg021.50.16125 rbg021.60.2110 rbg021.70.2270 rbg021.80.1588 rbg021.90.17108 rbg0210.19158 rbg027a0.2153parison of traditional LDS and the discrepancy constraint. goodfirst solution and proving optimality.This is done by tuning the ratio r, which determines the size of thefirst subproblem.We recall from Section3.2that|D goodi |≤rn for i∈{1,...,n}.In Tables2and3we report the quality of the heuristic with respect to the ra-tio.The TSP instances are taken from TSPLIB[23]and the asymmetric TSPTW instances are due to Ascheuer[1].All subproblems are solved to optimality with afixed strategy,as to make a fair comparison.In the tables,‘size’is the actual relative size of thefirst subproblem with respect to the initial problem.The size is calculated by1instance size opt pr fails size opt pr fails0.250120.320170.170110.230160.170120.210120.160110.24013270.14121270.21011k0.131619k0.200165average0.200 1.17680.250175 Table2.Quality of heuristic with respect to ratio for the TSP.The domains in the TSPTW instances are typically much smaller than the number of variables,because of the time window constraints that already remove a number of domain values.Therefore,thefirst subproblem might sometimes be relatively large,since only a few values are left after pruning,and they might be equally promising.The next columns in the tables are‘opt’and‘pr’.Here‘opt’denotes the level of discrepancy at which the optimal solution is found.Typically,we would11 like this to be0.The column‘pr’stands for the level of discrepancy at which optimality is proved.This would preferably be1.The column‘fails’denotes the total number of backtracks during search needed to solve the problem to optimality.Concerning the TSP instances,already for a ratio of0.05all solutions are in thefirst subproblem.For a ratio of0.075,we can also prove optimality at discrepancy1.Taking into account also the number of fails,we can argue that both0.05and0.075are good ratio candidates for the TSP instances.For the TSPTW instances,we notice that a smaller ratio does not necessarily increase the total number of fails,although it is more difficult to prove optimality. Hence,we have a slight preference for the ratio to be0.15,mainly because of its overall(average)performance.An important aspect we are currently investigating is the dynamic tuning of the ratio.6Computational ResultsWe have implemented and tested the proposed method using ILOG Solver and Scheduler[17,16].The algorithm runs on a Pentium1Ghz,256MB RAM,and uses CPLEX6.5as LP solver.The two sets of test instances are taken from the TSPLIB[23]and Ascheuer’s asymmetric TSPTW problem instances[1].Table4shows the results for small TSP instances.Time is measured in seconds,fails again denote the total number of backtracks to prove optimality. The time limit is set to300seconds.Observe that our method(LDS)needs less number of fails than the approach without subproblem generation(No LDS). This comes with a cost,but still our approach is never slower and in some cases considerably faster.The problems were solved both with a ratio of0.05and0.075, the best of which is reported in the table.Observe that in some cases a ratio of 0.05is best,while in other cases0.075is better.For three instances optimality could not be proven directly after solving thefirst subproblem.Nevertheless the optimum was found in this subproblem(indicated by objective‘obj’is‘opt’). Time and the number of backtracks(fails)needed in thefirst subproblem are reported for these instances.In Table5the results for the asymmetric TSPTW instances are shown.Our method(LDS)uses a ratio of0.15to solve all these problems to optimality.It is compared to our code without the subproblem generation(No LDS),and to the results by Focacci,Lodi and Milano(FLM2002)[10].Up to now,FLM2002 has the fastest solution times for this set of instances,to our knowledge.When comparing the time results(measured in seconds),one should take into account that FLM2002uses a Pentium III700MHz.Our method behaves in general quite well.In many cases it is much faster than FLM2002.However,in some cases the subproblem generation does not pay off.This is for instance the case for rbg040a and rbg042a.Although our method finds the optimal solution in thefirst branch(discrepancy0)quite fast,the initial bound LB0is too low to be able to prune the search tree at discrepancy1.In those12ratio=0.05ratio=0.1ratio=0.15ratio=0.2 size opt pr fails size opt pr failsrbg010a0.940150.99017 rbg016a0.8802410.980147 rbg016b0.7203540.910239 rbg017.20.490190.660127 rbg0170.6615700.890249 rbg017a0.7201130.8801122 rbg019a0.98012710130 rbg019b0.7813850.950171 rbg019c0.60131370.760299 rbg019d0.970161016 rbg020a0.840130.94015 rbg021.20.5901480.7401148 rbg021.30.55141850.6803163 rbg021.40.5303330.670287 rbg021.50.55021030.6702206 rbg021.60.48122330.6701124 rbg021.70.4902910.650170 rbg021.80.48125180.630188 rbg021.90.48125740.6301108 rbg0210.60131370.760299 rbg027a0.5902350.7701960.570.673.0517250.740.101.5769instance size opt pr failsrbg010a10181014910149rbg016b0.9901380.7401310.900120rbg0170.9901560.9301143101236rbg019a101300.97017210172rbg019c0.920111210161016rbg020a10150.80011850.9501240rbg021.30.90021420.74021140.940170rbg021.50.88011720.73011320.9501148rbg021.70.8201780.7101890.900195rbg021.90.83011210.82021060.9501119rbg027a0.9001407average0.930 1.0510113 instance time fails ratiogr170.1210.050.0719---gr240.1810.050.1770---bayg290.28280.050.30418opt0.1936dantzig42 1.213660.07512.6115k---gr48∗19.1022k0.075limit limit opt80.3155k∗Optimality not proven directly after solvingfirst subproblem.putational results for the TSP.instance time failsrbg010a0.0660.1210.048 rbg016b0.11270.0170.042 rbg0170.0690.1220.051 rbg019a0.05110.2800.0922 rbg019c0.08170.0320.074 rbg020a0.07110.2440.0815 rbg021.30.11520.31210.0932 rbg021.50.14890.73180.1650instance time failsrbg021.70.19450.62220.1027 rbg021.90.10310.3810.1474 rbg027a0.19452.78410.68119 rbg033a0.627055.213k0.9336 rbg035a.2 5.23 2.6k3.58410.8356 rbg038a0.3742738.1136k1k387k rbg042a29.3611k180.419k 4.21 1.5k rbg055a 4.391634.049325.69128putational results for the asymmetric TSPTW.instance time fails179179014214201481480190202 6.318218201791790opt obj gap(%)rbg021.50.1019 rbg021.60.1549 rbg0210.1044 rbg035a.2 4.75 2.2k rgb040a65.8325k rbg042a33.2611.6kTable putational results for thefirst subproblem of the asymmetric TSPTW14cases we need more time to prove optimality than we would have needed if we did not apply our method(No LDS).In such cases our method can be applied as an effective incomplete method,by only solvingfirst subproblem.Table6shows the results for those instances for which optimality could not be proven directly after solving thefirst subproblem.In almost all cases the optimum is found in thefirst subproblem.7Discussion and ConclusionWe have introduced an effective search procedure that consists of generating promising subproblems and solving them.To generate the subproblems,we split the variable domains in a good part and a bad part on the basis of reduced costs. The domain values corresponding to the lowest reduced costs are more likely to be in the optimal solution,and are put into the good set.The subproblems are generated using a LDS strategy,where the discrepancy is the number of variables ranging on their bad set.Subproblems are considerably smaller than the original problem,and can be solved faster.To prove optimality, we introduced a way of increasing the lower bound using information from the discrepancies.Computational results on TSP and asymmetric TSPTW instances show that the proposed ranking is extremely accurate.In almost all cases the optimal solution is found in thefirst subproblem.When proving optimality is difficult, our method can still be used as an effective incomplete search procedure,by only solving thefirst subproblem.Some interesting points arise from the paper:first,we have seen that reduced costs represent a good ranking criterion for variable domain values.The ranking quality is not affected if reduced costs come from a loose relaxation.Second,the tightness of the lower bound is instead very important for the proof of optimality.Therefore,we have used a tight bound at the root node and increased it thus obtaining a discrepancy-based bound.Third,if we are not interested in a complete algorithm,but we need very good solutions fast,our method turns out to be a very effective choice,since the first generated subproblem almost always contains the optimal solution.Future directions will explore a different way of generating subproblems.Our domain partitioning is statically defined only at the root node and maintained during the search.However,although being static,the method is still very effec-tive.We will explore a dynamic subproblem generation. AcknowledgementsWe would like to thank Filippo Focacci and Andrea Lodi for useful discussion, suggestions and ideas.。

集合覆盖问题的模型与算法王继强【摘要】The set cover problem has favourable applications in areas of network design, but it is NP-hard in computational com-plexity. A 0-1 program model is formulated for the set cover problem. An approximation algorithm deriving from greedy idea is put forward, and is proved from the angle of primal-dual program. A case study of sensor network optimal design based on LINGO software demonstrates correctness of the model and effectiveness of the algorithm.% 集合覆盖问题在网络设计领域中有着良好的应用背景，但它在算法复杂性上却是NP-困难问题。

建立了集合覆盖问题的0-1规划模型，给出了源于贪心思想的近似算法，并从原始-对偶规划的角度进行了证明，基于LINGO软件的传感器网络最优设计案例验证了模型的正确性和算法的有效性。

【期刊名称】《计算机工程与应用》【年(卷),期】2013(000)017【总页数】4页(P15-17,72)【关键词】集合覆盖;近似算法;0-1规划;对偶规划;线性交互式通用优化器(LINGO)【作者】王继强【作者单位】山东财经大学数学与数量经济学院，济南 250014【正文语种】中文【中图分类】Q784;TP301.61 引言集合覆盖问题是组合最优化和理论计算机科学中的一类典型问题，它要求以最小代价将某一集合利用其若干子集加以覆盖。

Unit OneTask 1⑩④⑧③⑥⑦②⑤①⑨Task 2① be consistent with他说，未来的改革必须符合自由贸易和开放投资的原则。

② specialize in启动成本较低，因为每个企业都可以只专门从事一个很窄的领域。

③ d erive from以上这些能力都源自一种叫机器学习的东西，它在许多现代人工智能应用中都处于核心地位。

④ A range of创业公司和成熟品牌推出的一系列穿戴式产品让人们欢欣鼓舞,跃跃欲试。

⑤ date back to置身硅谷的我们时常淹没在各种"新新"方式之中，我们常常忘记了，我们只是在重新发现一些可追溯至涉及商业根本的朴素教训。

Task 3T F F T FTask 4The most common viewThe principle task of engineering: To take into account the customers ‘ needs and to find the appropriate technical means to accommodate these needs.Commonly accepted claims:Technology tries to find appropriate means for given ends or desires;Technology is applied science;Technology is the aggregate of all technological artifacts;Technology is the total of all actions and institutions required to create artefacts or products and the total of all actions which make use of these artefacts or products.The author’s opinion: it is a viewpoint with flaws.Arguments: It must of course be taken for granted that the given simplified view of engineers with regard to technology has taken a turn within the last few decades. Observable changes: In many technical universities, the inter‐disciplinary courses arealready inherent parts of the curriculum.Task 5① 工程师对于自己的职业行为最常见的观点是：他们是通过应用科学结论来计划、开发、设计和推出技术产品的。

On the Smooth Implementation ofComponent-based System Specifications Antonio Albarr´a n†,Francisco Dur´a n‡,and Antonio Vallecillo‡†Junta de Andaluc´ıa,Spain.aalbarran@ceh.junta-andalucia.es‡Universidad de M´a laga,Spain.{duran,av}@lcc.uma.esAbstractMaude is an executable rewriting logic language particularly well suited for the specification of object-oriented open and distributed systems.CORBA is one of theleading object and component platforms in the market,specially well suited for theimplementation of distributed applications in open systems.In this paper we show howto connect both worlds,specification and implementation,so objects in any of themcan interact with objects in the other in a natural way.1IntroductionThere seems to be a growing consensus on the use of ready-to-use components for building software systems in order to reduce development time and effort,and to improve software quality.Whether using components or not,software systems have to be developed following rigorous guidelines to have a chance of being dependable,finished on time,and easy to maintain.Typically,the software development process is seen as the construction of a sequence—up to backtracking and iterations—of more and more detailed descriptions of the software under development,leading to afinal set of elements that contain a set of executable programs and their documentation.Making use of formal specifications is a demanding process and should be suitably tar-geted.From a logical point of view,the correctness of a program is a nonsense without a formal specification against which correctness can be proved.However,there are many cases where formal specifications would be a mere luxury.In most projects were formal specifi-cations are used,formal and informal specifications are mixed,so that some components or refinement steps are formalized and others not.The decision to use formal specifications mainly depends on the criticality of the component,in terms of consequences of a fault (human lives,costs,etc.)and on the complexity of its requirements or its development.The formal or informal global specification of a system is valuable by itself as a way of establishing a description of the system on which the developer and the client agree.In component software,we generally begin with a high-level graphic description of the system, using,for example,a UML-like notation for identifying the different components and their ponent interfaces are then identified and specified.Such specifications can be viewed as contracts between clients of an interface and the provider of an implementation of it,and they are used as descriptions of the components for specifying the interactions among the different components via some coordination or architecture description language. Typically,interface specifications include a description of the functional aspects of the operations,including its syntax and its semantics,given by invariants and pre-and post-conditions.It can be expected that non-functional requirements will also be included into such contracts.Over the last years,there has been a big effort by the research community for increasing the usability of formal methods,developing simpler to use and more powerful tools which can be more easily integrated into the software development process.Theorem proving, prototyping,model checking,static analysis,code generation,and testing are the kinds of activities where formal specifications have shown to play a useful role.For the present paper it is particularly relevant the possibility of prototyping:If specifications can be executed (or at least animated)they may be used for getting rapid prototypes of the systems,which may then be used for completing and clarifying such specifications.On the other hand,one of the benefits of component-based systems is the clear sepa-ration of functionality between their constituent parts,and the‘loose’connection among them.With the use of component platforms such as CORBA,different components can communicate without any assumption on their implementation language,the hardware they run on,or their physical location.This kind of component interaction(only through their interfaces)may help the reasoning process about the systems being built.In this paper we explore the possibility of connecting the worlds of a formal specifi-cation language such as Maude[1]and of a model for open distributed systems such as CORBA[5].Maude is an executable rewriting logic[3]language particularly well suited for the specification of object-oriented open and distributed systems[4].CORBA is one of the leading object and component platforms in the market,specially well suited for the implementation of distributed applications in open systems.By allowing objects of both worlds to seamlessly interoperate,we can obtain several interesting results.In thefirst place,suppose that we have the specification of the system written in Maude, in terms of the specifications of its constituent components.Until now we were able to prove some properties about the system,do some model-checking,or build a prototype by simply executing its Maude specification.However,there is a gap between the system specifications and a running implementation of it.In our proposal,Maude objects can be transparently replaced by CORBA objects(i.e.,their implementations)one by one,so that objects in both worlds coexist,while still being able to reason about the system.Secondly,black-box test suites for CORBA objects can be easily specified in Maude, based on their interfaces only.Those specifications can be later executed directly against the actual implementation of the CORBA object,thus testing its behavior.Third,given a running CORBA system and the specification of a new CORBA object (or a new version of one of its components),we can validate the new system without having to implement such a new object.We can just‘drop’the new specification into the running system,and check its behavior.Finally,and most importantly,in safe-critical systems the vital parts can be formally specified and their behavior formally verified,while non-critical components can be just informally specified and implemented.The system—composed of Maude specifications and some CORBA objects—is now available for execution,without having to wait for all critical components to be implemented.In this paper we discuss how the worlds of Maude and CORBA can be smoothly inte-grated,obtaining all these benefits.Sections2and3briefly describe CORBA and Maude, respectively.Then,Section4presents an example application to illustrate our proposal, and Section5briefly describes how to connect both worlds.Finally,Section6draws some conclusions and describes some future research activities.2CORBACORBA is one of the major distributed object platforms.Proposed by the OMG(http: //),the Object Management Architecture(OMA)attempts to define,at a high level of description,the various facilities required for distributed object-oriented com-puting.The core of the OMA is the Object Request Broker(ORB),a mechanism that provides transparency of object location,activation and communication.The Common Object Request Broker Architecture(CORBA)specification describes the interfaces and services that must be provided by compliant ORBs[5].In the OMA model,objects provide services,and clients issue requests for those services to be performed on their behalf.The purpose of the ORB is to deliver requests to objects and return any output values back to clients,in a transparent way to the client and the server. Clients need to know the object reference of the server object.ORBs use object references to identify and locate objects to redirect requests to them.As long as the referenced object exists,the ORB allows the holder of an object reference to request services from it.Even though an object reference identifies a particular object,it does not necessarily describe anything about the object’s interface.Before an application can make use of an object,it must know what services the object provides.CORBA defines an IDL to describe object interfaces,a textual language with a syntax resembling that of C++.The CORBA IDL provides basic data types(such as short,long,float,...),constructed types(struct, union)and template types(sequence,string).These are used to describe the interface of objects,defined by set of types,attributes and the signature(parameters,return types and exceptions raised)of the object methods,grouped into interface definitions.Finally,the construct module is used to hold type definitions,interfaces,and other modules for name scoping purposes.3Rewriting Logic and MaudeRewriting logic[3]is a logic in which the state space of a distributed system is specified as an algebraic data type in terms of an equational specification(Σ,E),whereΣis a signature of types(sorts)and operations,and E is a set of(conditional)equational axioms.The dynamics of such a system is then specified by rewrite rules of the form t→t ,where t and t areΣ-terms that describe the local,concurrent transitions possible in the system, i.e.,when a part of the system statefits the pattern t then it can change to a new local statefitting pattern t .The guards of conditional rules act as blocking pre-conditions,in the sense that a conditional rule can only befired if the condition is satisfied.Maude[1]is a high-level language and a high-performance interpreter and compiler that supports equational and rewriting logic specification and programming of systems. Thus,Maude integrates an equational style of functional programming with rewriting logic computation.This logic can naturally deal with state and with highly nondeterministic concurrent computations;in particular,it supports very well concurrent object-oriented computation[4].Such a logic allows us to specify functional as well as non-functional properties.In Maude,object-oriented systems are specified by object-oriented modules in which classes and subclasses are declared.Each class is declared using the syntax“class C|a1: S1,...,a n:S n”,where C is the name of the class,a i are attribute identifiers,and S i are the corresponding sorts of the attributes.Objects of a class C are then record-like structures of the form<O:C|a1:v1,...,a n:v n>,where O is the name of the object,and the v i are the current values of its attributes.Objects can interact in a different number of ways, including message passing.In a concurrent object-oriented system the concurrent state,which is called the con-figuration,has the structure of a multiset made up of objects and messages that evolves by concurrent rewriting using rules that describe the effects of the communication events between some objects and messages.The general form of such rewrite rules is crl[r]:M1...M m<O1:C1|atts1>...<O n:C n|atts n>−→<O i1:C i1|atts i1>...<O ik:C ik|atts ik><Q1:D1|atts 1>...<Q p:D p|atts p>M 1...M qif Cond.where r is the rule label,M i are messages,O i and Q j are object identifiers,C i,C j and D h are classes,i1,...,i k is a subset of1...n,and Cond is a boolean condition(the rule ‘guard’).The result of applying such a rule is that:•messages M1...M m disappear,i.e.,they are consumed;•the state,and possibly the classes of objects O i1,...,O ikmay change;•all the other objects O j vanish;•new objects Q1,...,Q p are created;and•new messages M 1...M q are created,i.e.,they are sent.4A Case StudyIn order to illustrate our proposal,let us suppose that we have already designed a system (using UML or any other notation),and that we have decided to implement it in CORBA. Among other entities,the system contains a bank account component,which uses a cal-culator for doing the computations.This is a simple example,but will serve us to model the typical client-server interactions between two objects.The CORBA interfaces of these components follow.interface Account{interface Calculator{void deposit(in int amount);int plus(in int a,in int b);void withdraw(in int amount);int minus(in int a,in int b);int balance();};};Imagine that they are critical components,whose behavior needs to be formally specified in order to reason about ing Maude,their specifications can be written as follows: (omod CALCULATOR isprotecting BASIC-TYPES.class Calculator|.msgs plusRequest minusRequest:int int Oid Oid->Msg.msgs plusResponse minusResponse:int int int Oid Oid->Msg.vars O O’:Oid.vars X Y:int.rl[plus]:<O:Calculator|>plusRequest(X,Y,O,O’)=><O:Calculator|>plusResponse(X,Y,(X+Y),O’,O).rl[minus]:<O:Calculator|>minusRequest(X,Y,O,O’)=><O:Calculator|>minusResponse(X,Y,(X-Y),O’,O).endom)(omod BANK-ACCOUNT isprotecting BASIC-TYPES+DEFAULT[Oid].class Account|bal:int,client:Default[Oid],calc:Calculator.msgs depositRequest withdrawRequest:int Oid Oid->Msg.msg balanceRequest:Oid Oid->Msg.msgs depositResponse withdrawResponse:Oid Oid->Msg.msg balanceResponse:int Oid Oid->Msg.vars C A O:Oid.vars B X Z:int.rl[deposit]:<A:Account|bal:B,client:null,calc:C>depositRequest(X,A,O)=><A:Account|bal:B,client:O,calc:C>plusRequest(B,X,C,A).rl[depositResponse]:<A:Account|bal:B,client:O,calc:C>plusResponse(B,X,Z,A,C)=><A:Account|bal:Z,client:null,calc:C>depositResponse(O,A)....***Rest of the specification omittedendom)As we can see in these specifications,a call to the account’s deposit method causes a call to the calculator’s plus method to do the addition.The bank account is then a client of the calculator,whose reference is kept as part of the account object attributes.Please note the following mechanisms used for modeling the behavior of CORBA objects in Maude:•The natural communication mechanism between objects in Maude is message pass-ing.Thus,we have modeled the RPC mechanism used in CORBA for client-server interactions with two messages,one that carries the request,and one for the reply.•CORBA types(such as“int”)are available in Maude by means of a Maude module BASIC-TYPES that contains their specifications.•Apart from the method arguments,all messages incorporate two special arguments at the end:the object identifiers of the receiving and the calling objects.Suppose that we have proved that the system is safe and that we are ready to implement it.Let us imagine that we have sourced a CORBA implementation of the calculator interfaceFigure1:Replacing Maude objects by their CORBA implementations.from a third party,and that we want to execute it,testing it against the specification of its client,i.e.,the bank account.The next section shows how both worlds—CORBA and Maude—can be integrated,and how their objects interoperate.5Connecting both WorldsOur proposal for the integration is based on the use of proxies.An object in one world is replaced by a proxy,which captures its incoming calls and passes them to another proxy in the second world.This proxy simulates the client calls in the second world,awaits for the server response,and then passes it back to thefirst proxy,which builds the answer to the original request from it.The communication between objects and proxies is accomplished by using the natural communication mechanisms available in each world,and communication between proxies is implemented via TCP,using the new features available in Maude2.01 [2].Figure1shows a gradual‘refinement’of a Maude specification into a running CORBA system.The three Maude components in the left are replaced by their corresponding implementations step by step.All objects are of course unaware of the replacement process, since we use proxies that substitute the Maude components being replaced.This is just an example,since the number of components that can be replaced in each step may vary depending on the implementator’s choice.Let us see how those proxies look like in the particular case of the bank account and the calculator.Suppose that we want to replace the calculator specification by its implemen-tation,making the bank account specification interact with it.In thefirst place we need a proxy of the calculator in the Maude world:(omod CALCULATOR-PROXY isincluding TCPHandler.protecting STRING+DEFAULT[Oid]+BASIC-TYPES.class Calculator|CORBAProxyAddr:String,client:Default[Oid].msgs plusRequest minusRequest:int int Oid Oid->Msg.msgs plusResponse minusResponse:int int int Oid Oid->Msg.1Maude2.0is not available yet.The examples in this paper have been executed using Java objects to pull messages in and out the Maude objects.Our believe is that once Maude2.0is available these wrappers will be removed without affecting the rest of the system.vars O O’:Oid.vars X Y Z:int.rl[plusRequest]:<O:Calculator|CORBAProxyAddr:A,client:null>plusRequest(X,Y,O,O’)=><O:Calculator|client:O’>sendTCP(System,O,A,"<m:plusRequest>"+"<a type=‘"int‘">"+X+"</a>"+"<b type=‘"int‘">"+Y+"</b>"+"</m:plusRequest>").rl[plusResponse]:<O:Calculator|client:O’>receivedTCP(O,System,"<m:plusResponse>"+"<a type=‘"int‘">"+X+"</a>"+"<b type=‘"int‘">"+Y+"</b>"+"<result type=‘"int‘">"+Z"</result>"+"</m:plusResponse>")=><O:Calculator|client:null>plusResponse(X,Y,Z,O’,O).endom)This object is the one that provides now the services originally provided by the Calculator class.However,although it accepts the same methods as the original calculator,it does not carry out the computations any longer:it just passes the information around,packing and unpacking messages into the appropriate formats.Upon receipt of a method call,this proxy sends by TCP the request to the appropriate CORBA proxy(whose IP address keeps as an attribute).And once the response from the CORBA proxy is received,it sends back the answer to the original client.It is important to notice how the automated construction of this kind of proxies is straight forward,directly from the CORBA IDLs of the objects.At the CORBA side,the structure of the CORBA proxy is very similar.The following piece of code shows a possible implementation the CORBA proxy in Java.Calculator calc=....;//obtains the CORBA object ref.for calculatorString addr="150.214.108.1:2020";//IP:socket of MaudeProxyfor(;;){String s=SocketReadLn(addr);if(getMethod(s).equals("<m:plusRequest>")){int a=(new Integer(getArg(s,1))).intValue();int b=(new Integer(getArg(s,2))).intValue();int res=calc.plus(a,b);//call to the calculator CORBA objectSocketWriteLn(addr,"<m:plusResponse>"+"<a type=\"int\">"+a+"</a>"+"<b type=\"int\">"+b+"</b>"+"<result type=\"int\">"+res+"</result>"+"</m:plusResponse>");}else if(getMethod(s).equals("<m:minusRequest>")){...}}}Note that the generation of this proxy can also be automated directly from the CORBA IDL of the original object.6Concluding RemarksIn this paper we have shown how a connection between a formal notation and a commercial component platform can be built.This connection can be very helpful for bridging the gap between system specification and implementation.Our proposal provides a transpar-ent migration between objects in both worlds,so system specifications can be smoothly replaced by their implementations,as well as permitting new objects’specifications to be executed within running CORBA systems,and thus allowing a veryflexible way of sys-tems/components prototyping.As pointed out by one of the anonymous referees,formal specifications are not adequate for describing large real-world systems.However,we believe that their decomposition into suitably sized components in such a way that the critical ones can then be selected and formally specified may help greatly.Our main goal in the piece of work reported in this paper was to study how the con-nection between formal specifications and implementations of them could be implemented. Now,we plan to study the possibility of improving it by using SOAP,SCOAP,or some Web-based RPC protocol,that will extend its capabilities beyond CORBA objects,thus allowing us to connect Maude specifications to other object and component platforms. References[1]Manuel Clavel,Francisco Dur´a n,Steven Eker,Patrick Lincoln,Narciso Mart´ı-Oliet,Jos´e Meseguer,and Jos´e Quesada.Maude:Specification and programming in rewriting logic.Manuscript,SRI International,1999.Available at .[2]Manuel Clavel,Francisco Dur´a n,Steven Eker,Patrick Lincoln,Narciso Mart´ı-Oliet,Jos´e Meseguer,and Jos´e Quesada.Towards Maude2.0.In Kokichi Futatsugi,edi-tor,Proceedings of3rd International Workshop on Rewriting Logic and its Applications (WRLA’00),volume36of Electronic Notes in Theoretical Computer Science.Elsevier, 2000.[3]Jos´e Meseguer.Conditional rewriting logic as a unified model of concurrency.TheoreticalComputer Science,96:73–155,1992.[4]Jos´e Meseguer.A logical theory of concurrent objects and its realization in the Maudelanguage.In Gul Agha,Peter Wegner,and Akinori Yonezawa,editors,Research Direc-tions in Object-Based Concurrency,pages314–390.The MIT Press,1993.[5]OMG.The Common Object Request Broker:Architecture and Specification.ObjectManagement Group,2.4edition,November2000./technology/ documents/formal/corbaiiop.htm.。

12IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 26, NO. 1, FEBRUARY 2011MATPOWER: Steady-State Operations, Planning, and Analysis Tools for Power Systems Research and EducationRay Daniel Zimmerman, Member, IEEE, Carlos Edmundo Murillo-Sánchez, Member, IEEE, and Robert John Thomas, Life Fellow, IEEEAbstract—MATPOWER is an open-source Matlab-based power system simulation package that provides a high-level set of power flow, optimal power flow (OPF), and other tools targeted toward researchers, educators, and students. The OPF architecture is designed to be extensible, making it easy to add user-defined variables, costs, and constraints to the standard OPF problem. This paper presents the details of the network modeling and problem formulations used by MATPOWER, including its extensible OPF architecture. This structure is used internally to implement several extensions to the standard OPF problem, including piece-wise linear cost functions, dispatchable loads, generator capability curves, and branch angle difference limits. Simulation results are presented for a number of test cases comparing the performance of several available OPF solvers and demonstrating MATPOWER’s ability to solve large-scale AC and DC OPF problems. Index Terms—Load flow analysis, optimal power flow, optimization methods, power engineering, power engineering education, power system economics, power system simulation, power systems, simulation software, software tools.I. INTRODUCTIONTHIS paper describes MATPOWER, an open-source Matlab power system simulation package [1]. It is used widely in research and education for AC and DC power flow and optimal power flow (OPF) simulations. It also includes tools for running OPF-based auction markets and co-optimizing reserves and energy. Included in the distribution are numerous example power flow and OPF cases, ranging from a trivial four-bus example to real-world cases with a few thousand buses. MATPOWER consists of a set of Matlab M-files designed to give the best performance possible while keeping the code simple to understand and customize. Matlab has become a popular tool for scientific computing, combining a high-level language ideal for matrix and vector computations, a cross-platform runtime withManuscript received December 22, 2009; revised April 19, 2010; accepted May 17, 2010. Date of publication June 21, 2010; date of current version January 21, 2011. This work was supported in part by the Consortium for Electric Reliability Technology Solutions and the Office of Electricity Delivery and Energy Reliability, Transmission Reliability Program of the U.S. Department of Energy under the National Energy Technology Laboratory Cooperative Agreement No. DE-FC26-09NT43321. Paper no. TPWRS-00995-2009. R. D. Zimmerman and R. J. Thomas are with the Department of Applied Economics and Management and the School of Electrical and Computer Engineering, Cornell University, Ithaca, NY 14853 USA (e-mail: rz10@; rjt1@). C. E. Murillo-Sánchez is with the Universidad Autónoma de Manizales, and with Universidad Nacional de Colombia, both in Manizales, Colombia (e-mail: carlos_murillo@). Digital Object Identifier 10.1109/TPWRS.2010.2051168robust math libraries, an integrated development environment and GUI with excellent visualization capabilities, and an active community of users and developers. As a high-level scientific computing language, it is well suited for the numerical computation typical of steady-state power system simulations. The initial motivation for the development of the Matlabbased power flow and OPF code that would eventually become MATPOWER arose from the computational requirements of the PowerWeb platform [3], [4]. As a web-based market simulation platform used to test electricity markets, PowerWeb requires a “smart market” auction clearing software that uses an OPF to compute the allocations and pricing. Having the clear potential to be useful to other researchers and educators, the software was released in 1997 via the Internet as an open-source power system simulation package, now distributed under the GNU GPL [2]. Even beyond its initial release, much of the ongoing development of MATPOWER continued to be driven in large part by the needs of the PowerWeb project. This at least partially explains the lack of a graphical user interface used by some related tools such as PSAT [5]. While it is often employed as an end-user tool for simply running one-shot simulations defined via an input case file, the package can also be quite valuable as a library of functions for use in custom code developed for one’s own research. At this lower level, MATPOWER provides easy-to-use functions for and matrices, calculating forming standard network power transfer and line outage distribution factors (PTDFs and LODFs), and efficiently computing first and second derivatives of the power flow equations, among other things. At a higher level, the structure of the OPF implementation is explicitly designed to be extensible [6], allowing for the addition of user-defined variables, costs, and linear constraints. The default OPF solver is a high-performance primal-dual interior point solver implemented in pure-Matlab. This solver has application to general nonlinear optimization problems outside of MATPOWER and comes with a convenience wrapper function to make it trivial to set up and solve linear programming (LP) and quadratic programming (QP) problems. To help ensure the quality of the code, MATPOWER includes an extensive suite of automated tests. Some may find the testing framework useful for creating automated tests for their own Matlab programs. A number of Matlab-based software packages related to power system simulation have been developed by others. A nice summary of their features is presented in [5]. The primary distinguishing characteristics of MATPOWER, aside from0885-8950/$26.00 © 2010 IEEEZIMMERMAN et al.: MATPOWER: STEADY-STATE OPERATIONS, PLANNING, AND ANALYSIS TOOLS13With the series admittance element in the model denoted by , the branch admittance matrix can be written (2) If the four elements of this matrix for branch are labeled as follows: (3)Fig. 1. Branch model.being one of the first to be publicly and freely available as open-source, are the extensible architecture of the OPF formulation and its ease of use as a toolbox of functions to incorporate into one’s own programs. It is also compatible with Octave. This paper describes the MATPOWER package as it stands at version 4, detailing the component modeling in Section II, the power flow and optimal power flow formulations in Sections III and IV, and some additional functionality in Section V. Some example results and conclusions are presented in Section VI. II. MODELING MATPOWER employs all of the standard steady-state models typically used for power flow analysis. The AC models are described first, then the simplified DC models. Internally, the magnitudes of all values are expressed in per unit and angles of complex quantities are expressed in radians. Due to the strengths of the Matlab programming language in handling matrices and vectors, the models and equations are presented here in matrix and vector form. A. Data Formats The data files used by MATPOWER are Matlab M-files or MAT-files which define and return a single Matlab struct. The M-file format is plain text that can be edited using any standard text editor. The fields of the struct are baseMVA, bus, branch, gen, and optionally gencost, where baseMVA is a scalar and the rest are matrices. In the matrices, each row corresponds to a single bus, branch, or generator. The columns are similar to the columns in the standard IEEE CDF and PTI formats. The number of rows in bus, branch, and gen are , , and , respectively. B. Branches All transmission lines, transformers, and phase shifters are modeled with a common branch model, consisting of a standard transmission line model, with series impedance and total charging capacitance , in series with an ideal phase shifting transformer. The transformer, whose tap ratio has magnitude and phase shift angle , is located at the from end of the branch, as shown in Fig. 1. The complex current injections and at the from and to ends of the branch, respectively, can be expressed in terms of the 2 2 branch admittance matrix and the respective terminal voltages and (1)then four vectors , , , and can be constructed, where the th element of each comes from the corresponding element of . Furthermore, the sparse connection maand used in building the system admittance matrices th element of and trices can be defined as follows. The the th element of are equal to 1 for each branch , where branch connects from bus to bus . All other elements of and are zero. C. Generators A generator is modeled as a complex power injection at a specific bus. For generator , the injection is (4) be the vector of these generator Let generator connection matrix injections. A sparse can be defined such that its th element is 1 if generator is located at bus and 0 otherwise. The vector of all bus injections from generators can then be expressed as (5) D. Loads Constant power loads are modeled as a specified quantity of real and reactive power consumed at a bus. For bus , the load is (6) and denotes the vector of complex loads at all buses. Constant impedance and constant current loads are not implemented directly, but the constant impedance portions can be modeled as a shunt element described below. Dispatchable loads are modeled as negative generators and appear as negative values in . E. Shunt Elements A shunt connected element such as a capacitor or inductor is modeled as a fixed impedance to ground at a bus. The admittance of the shunt element at bus is given as (7) denotes the and admittances at all buses. F. Network Equations For a network with buses, all constant impedance elements of the model are incorporated into a complex bus advector of shunt14IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 26, NO. 1, FEBRUARY 2011mittance matrix that relates the complex nodal current into the complex node voltages : jections (8) branches, the system Similarly, for a network with and relate the bus voltages to branch admittance matrices the vectors and of branch currents at the from and to ends of all branches, respectively: (9) (10) vector and If is used to denote an operator that takes an creates the corresponding diagonal matrix with the vector elements on the diagonal, these system admittance matrices can be formed as follows: (11) (12) (13) The current injections of (8)–(10) can be used to compute the corresponding complex power injections as functions of the complex bus voltages : (14) (15) (16) The nodal bus injections are then matched to the injections from loads and generators to form the AC nodal power balance equations, expressed as a function of the complex bus voltages and generator injections in complex matrix form as (17) G. DC Modeling The DC formulation [11] (with more detailed derivations in [1]) is based on the same parameters, but with the following three additional simplifying assumptions. • Branches can be considered lossless. In particular, branch resistances and charging capacitances are negligible: (18) • All bus voltage magnitudes are close to 1 p.u. (19) • Voltage angle differences across branches are small enough that (20) By combining (1) and (2) with (18) and (19), the complex current flow in a branch can be approximated as (21)Furthermore, using (19) and this approximate current to compute the complex power flow, then extracting the real part and applying the last of the DC modeling assumptions from (20) yields (22) As expected, given the lossless assumption, a similar derivation . for leads to The relationship between the real power flows and voltage angles for an individual branch can then be summarized as (23)where,, andisdefined in terms of the series reactance and tap ratio for that . branch as With a DC model, the linear network equations relate real power to bus voltage angles, versus complex currents to comvector plex bus voltages in the AC case. Let the be constructed similar to , where the th element is and let be the vector whose th element is equal to . Then the nodal real power injections can be expressed vector of bus voltage angles as a linear function of , the (24) where . Similarly, the branch flows at the from ends of each branch are linear functions of the bus voltage angles (25) and, due to the lossless assumption, the flows at the to ends are given by . The construction of the system matrices is analogous to the system matrices for the AC model: (26) (27) The DC nodal power balance equations for the system can be expressed in matrix form as (28) where approximates the amount of power consumed by the constant impedance shunt elements under the voltage assumption of (19). III. POWER FLOW The standard power flow or loadflow problem involves solving for the set of voltages and flows in a network corresponding to a specified pattern of load and generation. MATPOWER includes solvers for both AC and DC power flow problems, both of which involve solving a set of equations of the form (29)ZIMMERMAN et al.: MATPOWER: STEADY-STATE OPERATIONS, PLANNING, AND ANALYSIS TOOLS15constructed by expressing a subset of the nodal power balance equations as functions of unknown voltage quantities. All of MATPOWER’s solvers exploit the sparsity of the problem and, except for Gauss-Seidel, scale well to very large systems. Currently, none of them include any automatic updating of transformer taps or other techniques to attempt to satisfy typical OPF constraints, such as generator, voltage, or branch flow limits. A. AC Power Flow In MATPOWER, by convention, a single generator bus is typically chosen as a reference bus to serve the roles of both a voltage angle reference and a real power slack. The voltage angle at the reference bus has a known value, but the real power generation at the slack bus is taken as unknown to avoid overspecifying the problem. The remaining generator buses are classified as PV buses, with the values of voltage magnitude and and generator real power injection given. Since the loads are also given, all non-generator buses are PQ buses, with real , , and and reactive injections fully specified. Let denote the sets of bus indices of the reference bus, PV buses, and PQ buses, respectively. In the traditional formulation of the AC power flow problem, the power balance equation in (17) is split into its real and reactive components, expressed as functions of the voltage angles and magnitudes and generator injections and , where the load injections are assumed constant and given: (30) (31) For the AC power flow problem, the function from (29) is formed by taking the left-hand side of the real power balance equations (30) for all non-slack buses and the reactive power balance equations (31) for all PQ buses and plugging in the reference angle, the loads and the known generator injections and voltage magnitudes: (32) The vector consists of the remaining unknown voltage quantities, namely the voltage angles at all non-reference buses and the voltage magnitudes at PQ buses: (33)updated value of by factorizing this Jacobian. This method is described in detail in many textbooks. Also included are solvers based on variations of the fastdecoupled method [8], specifically, the XB and BX methods described in [9]. These solvers greatly reduce the amount of computation per iteration, by updating the voltage magnitudes and angles separately based on constant approximate Jacobians which are factored only once at the beginning of the solution process. These per-iteration savings, however, come at the cost of more iterations. The fourth algorithm is the standard GaussSeidel method from Glimm and Stagg [10]. It has numerous disadvantages relative to the Newton method and is included primarily for academic interest. By default, the AC power flow solvers simply solve the problem described above, ignoring any generator limits, branch flow limits, voltage magnitude limits, etc. However, there is an option that allows for the generator reactive power limits to be respected at the expense of the voltage setpoint. This is done by adding an outer loop around the AC power flow solution. If any generator has a violated reactive power limit, its reactive injection is fixed at the limit, the corresponding bus is converted to a PQ bus, and the power flow is solved again. This procedure is repeated until there are no more violations. B. DC Power Flow For the DC power flow problem [11], the vector the set of voltage angles at non-reference buses consists of (34) and (29) takes the form (35) where is the matrix obtained by simply eliminating from the row and column corresponding to the slack bus and reference angle, respectively. Given that the generator injections are specified at all but the slack bus, can be formed directly from the non-slack rows of the last four terms of (28). The voltage angles in are computed by a direct solution of the set of linear equations. The branch flows and slack bus generator injection are then calculated directly from the bus voltage angles via (25) and the appropriate row in (28), respectively. C. Linear Shift FactorsThis yields a system of nonlinear equations with equations and unknowns, where and are the number of PV and PQ buses, respectively. After solving for , the remaining real power balance equation can be used to compute the generator real power injection at the slack bus. Similarly, the remaining reactive power balance equations yield the generator reactive power injections. MATPOWER includes four different algorithms for solving the AC power flow problem. The default solver is based on a standard Newton’s method [7] using a polar form and a full Jacobian updated at each iteration. Each Newton step involves computing the mismatch , forming the Jacobian based on the sensitivities of these mismatches to changes in and solving for anThe DC power flow model can also be used to compute the sensitivities of branch flows to changes in nodal real power injections, sometimes called injection shift factors (ISF) or generation shift factors [11]. These sensitivity matrices, also called power transfer distribution factors or PTDFs, carry an implicit assumption about the slack distribution. If is used to denote a PTDF matrix, then the element in row and column , , represents the change in the real power flow in branch given a unit increase in the power injected at bus , with the assumption that the additional unit of power is extracted according to some specified slack distribution: (36)16IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 26, NO. 1, FEBRUARY 2011This slack distribution can be expressed as an vector of non-negative weights whose elements sum to 1. Each element specifies the proportion of the slack taken up at each bus. For the special case of a single slack bus , is equal to the vector . The corresponding PTDF matrix can be constructed by first matrix creating the (37) and then inserting a column of zeros at column . Here are obtained from and , respectively, by eliminating their reference bus columns and, in the case of , removing row corresponding to the slack bus. The PTDF matrix , corresponding to a general slack dis, tribution , can be obtained from any other PTDF, such as by subtracting from each column, equivalent to the following simple matrix multiplication: (38) These same linear shift factors may also be used to compute sensitivities of branch flows to branch outages, known as line outage distribution factors or LODFs [12]. Given a PTDF matrix , the corresponding LODF matrix can be constructed as follows, where is the element in row and column , representing the change in flow in branch (as a fraction of its initial flow) for an outage of branch . First, let represent the matrix of sensitivities of branch flows to branch flows, found by multplying the PTDF matrix by the node-branch incidence matrix: (39) If is the sensitivity of flow in branch with respect to flow in branch , then can be expressed as (40) MATPOWER includes functions for computing both the DC PTDF matrix and the corresponding LODF matrix for either a single slack bus or a general slack distribution vector . IV. OPTIMAL POWER FLOW MATPOWER includes code to solve both AC and DC versions of the optimal power flow problem. The standard version of each takes the following form: (41) (42) (43) (44) A. Standard AC OPF The optimization vector for the standard AC OPF problem consists of the vectors of voltage angles and magni-and the tudes power injectionsandvectors of generator real and reactive :(45)The objective function (41) is simply a summation of individual and of real and reactive power polynomial cost functions injections, respectively, for each generator: (46) The equality constraints in (42) are simply the full set of nonlinear real and reactive power balance equations from (30) and (31). The inequality constraints (43) consist of two sets of branch flow limits as nonlinear functions of the bus voltage angles and magnitudes, one for the from end and one for the to end of each branch: (47) (48) The flows are typically apparent power flows expressed in MVA, but can be real power or current flows, yielding the following three possible forms for the flow constraints: (49)where is defined in (9), in (15), , and the vector of flow limits has the appropriate units for the type of constraint. It is likewise for . The variable limits (44) include an equality constraint on any reference bus angle and upper and lower limits on all bus voltage magnitudes and real and reactive generator injections: (50) (51) (52) (53) B. Standard DC OPF When using DC network modeling assumptions and limiting polynomial costs to second order, the standard OPF problem above can be simplified to a quadratic program, with linear constraints and a quadratic cost function. In this case, the voltage magnitudes and reactive powers are eliminated from the problem completely and real power flows are modeled as linear functions of the voltage angles. The optimization variable is (54)ZIMMERMAN et al.: MATPOWER: STEADY-STATE OPERATIONS, PLANNING, AND ANALYSIS TOOLS17and the overall problem reduces to (55)–(60) at the bottom of the page. C. Extended OPF Formulation MATPOWER employs an extensible OPF structure [6] to allow the user to modify or augment the problem formulation without rewriting the portions that are shared with the standard OPF formulation. This is done through optional input parameters, preserving the ability to use pre-compiled solvers. The standard formulation is modified by introducing additional optional user-defined costs , constraints, and variables and can be written in the following form: (61) (62) (63) (64) (65) (66) The user-defined cost function is specified in terms a set of parameters in a pre-defined form described in detail in [6]. This form provides the flexibility to handle a wide range of costs, from simple linear functions of the optimization variables to scaled quadratic penalties on quantities, such as voltages, lying outside a desired range, to functions of linear combinations of variables, inspired by the requirements of price coordination terms found in the decomposition of large loosely coupled problems encountered in our own research. D. Standard Extensions In addition to making this extensible OPF structure available to end users, MATPOWER also takes advantage of it internally to implement several additional capabilities. 1) Piecewise Linear Costs: The standard OPF formulation in (41)–(44) does not directly handle the non-smooth piecewise linear cost functions that typically arise from discrete bids and offers in electricity markets. When such cost functions are convex, however, they can be modeled using a constrained cost variable (CCV) method. The piecewise linear cost function is replaced by a helper variable and a set of linear constraints that form a convex “basin” requiring the cost variable to lie in the epigraph of the function .A convex -segment piecewise linear cost function. . .. . .(67), can be defined by a sequence of points where denotes the slope of the th segment, (68)and and . The “basin” corresponding to this cost function is formed by the following constraints on the helper cost variable : (69) The cost term added to the objective function in place of is simply the variable . For an AC or DC OPF, MATPOWER uses this CCV approach internally to automatically generate the appropriate helper variable, cost term, and corresponding set of constraints for any piecewise linear generator costs. 2) Dispatchable Loads: A simple approach to dispatchable or price-sensitive loads is to model them as negative real power injections with associated negative costs. This is done by specifying a generator with a negative output, ranging from a minimum injection equal to the negative of the largest possible load to a maximum injection of zero. With this model, if the negative cost corresponds to a benefit for consumption, minimizing of generation is equivalent to maximizing social the cost welfare. With an AC network model, there is also the question of reactive dispatch for such loads. In MATPOWER, it is assumed that dispatchable loads maintain a constant power factor and an additional equality constraint is automatically added to enforce this requirement for any “negative generator” being used to model a dispatchable load. 3) Generator Capability Curves: The typical AC OPF formulation includes simple box constraints on a generator’s real and reactive injections. However, the true - capability curves of physical generators usually involve some tradeoff between real and reactive capability. If the user provides the parameters defining this tradeoff for a generator, MATPOWER automatically constructs the corresponding constraints. 4) Branch Angle Difference Limits: The difference between the bus voltage angle at the from end of a branch and the(55) (56) (57) (58) (59) (60)18IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 26, NO. 1, FEBRUARY 2011TABLE I OPF TEST CASESangle at the to end can be bounded above and below to act as a proxy for a transient stability limit, for example. If these limits are provided, MATPOWER creates the corresponding constraints on the voltage angle variables. E. Solvers Early versions of MATPOWER relied on Matlab’s Optimization Toolbox [13] to provide the NLP and QP solvers needed to solve the AC and DC OPF problems, respectively. While they worked reasonably well for very small systems, they did not scale well to larger networks. Eventually, optional packages with additional solvers were added to improve performance, typically relying on Matlab extension (MEX) files implemented in Fortran or C and pre-compiled for each machine architecture. For DC optimal power flow, there is a MEX build [14] of the high performance BPMPD solver [15] for LP/QP problems. For the AC OPF problem, the MINOPF [16] and TSPOPF [17] packages provide solvers suitable for much larger systems. The former is based on MINOS [18] and the latter includes the primal-dual interior point and trust region based augmented Lagrangian methods described in [19]. Beginning with version 4, MATPOWER also includes its own primal-dual interior point solver (MIPS) implemented in pureMatlab code, derived from the MEX implementation of the corresponding algorithms in [19]. If no optional packages are installed, the MIPS solver will be used by default for both the AC OPF and as the QP solver used by the DC OPF. The AC OPF solver also employs a unique technique for efficiently forming the required Hessians via a few simple matrix operations. The MIPS solver has application to general nonlinear optimization problems outside of MATPOWER and comes with a convenience wrapper function to make it trivial to set up and solve LP and QP problems. V. ADDITIONAL FUNCTIONALITY As mentioned earlier, MATPOWER was birthed out of a need for an OPF-based electricity auction clearing mechanism for a “smart market”. In this context, offers to sell and bids to buy power from generators and loads define the “costs” for the OPF that determines the allocations and prices used to clear the auction. MATPOWER includes code that takes bids and offers for real or reactive power, sets up and runs the corresponding OPF, and returns the cleared bids and offers. The standard OPF formulation described above includes no mechanism for completely shutting down generators which are very expensive to operate. Instead they are simply dispatched attheir minimum generation limits. MATPOWER includes the capability to run an OPF combined with a unit de-commitment for a single time period, which allows it to shut down these expensive units and find a least cost commitment and dispatch using an algorithm similar to dynamic programming. In some cases, it may be desirable to further constrain an OPF solution with the requirement to hold a specified level of capacity in reserve to cover contingencies. MATPOWER includes OPF extensions that allow it to co-optimize energy and reserves, subject to a set of fixed zonal reserve requirements. This code also serves as an example of how to customize the standard OPF with additional variables, costs, and constraints. VI. RESULTS AND CONCLUSIONS Several example cases are used to compare the performance of the various OPF solvers on example networks ranging in size from nine buses and three generators to tens of thousands of buses, thousands of generators and tens of thousands of additional user variables and constraints. Table I summarizes the test cases in terms of the order of the cost function (quadratic or linear), numbers of buses, generators and branches ( , , and ), numbers of variables and constraints ( and ) for both AC and DC OPF formulations, and the number of binding lower voltage limits and branch flow limits for for the DC case. the AC problem and flow limits For each case, six different AC OPF solvers and four different DC OPF solvers were used to solve the problem on a laptop with a 2.33-GHz Intel Core 2 Duo processor running Matlab 7.9. Table II gives the run times in seconds for the solvers which were successful, with the fastest time highlighted in bold for each example. The first algorithm listed for each is from Matlab’s Optimization Toolbox, in the case of the AC OPF and or for the DC problem. Next are the standard and step-controlled variants of the pure-Matlab implementation of the primal-dual interior point method, and last are some of the C and Fortran-based MEX solvers distributed as MATPOWER optional packages. For small systems, the clear winners are MINOPF for AC and BPMPD for DC, both Fortran-based MEX files. For larger systems, the primal-dual interior point solvers have the clear advantage, with the pure-Matlab implementation offering respectable performance relative to the C-based MEX versions. MATPOWER provides a high-level set of power flow and optimal power flow tools for researchers, educators, and students. The optimal power flow is extensible, allowing for easy。

内点法迭代原理及工程实例求解应用摘要：内点法是一种求解线性规划和非线性规划问题的多项式算法，其迭代次数与系统规模关系不大。

目前，内点法被扩展运用于求解二次规划模型，其计算速度和处理不等式约束的能力已经超过了求解二次规划模型的经典算法。

本文主要介绍线性规划中内点法的运用以及对工程实例的计算，并且分析了如何运用内点法迭代原理得到最优解。

关键字：线性规划问题；内点法；最优解；二次规划；1 引言1984年，Karmarkar发现了一个关于求解线性规划的方法，这个方法称作内点法。

内点法是罚函数中的一种，与外点法的最大的区别在于该方法利用罚函数生成一系列内点来逼近原约束问题的最优解。

罚函数的作用是对企图脱离可行域的点给予惩罚，相当于在可行域的边界设置了障碍，不让迭代点穿越到可行域之外。

内点法在迭代中总是从可行点出发，并保持在可行域内部进行搜索。

后得出最优解。

对于不等式约束的最优化问题，比较适合用内点法来解决。

经过实际计算结果得出内点法与单纯形法存在着很大的可比性。

在线性规划问题中，内点法比起单纯形法来说迭代次数更少，所以计算速度更快，从求得的结果来看，收敛性也比较好。

内点法中比较常用的方法是最速下降法和牛顿法。

最速下降法在解析法中是属于比较古老的一种，受该方法的启发，渐渐得到了其他不同的解析方法。

最速下降法每次迭代的计算量很小，解法简单。

如果从一个不好的初始点出发，也能收敛到局部极小点。

迭代原理的应用对于解决线性规划和非线性规划问题中具有至关重要的作用。

2 内点法2.1运筹学运筹学[1]到现在都没有一个相对比较统一的定义，这正是因为它使用的复杂性以及使用的广泛性，也凸显出了它另一方面的独特魅力。

以下是我查阅大量书籍后对运筹学所给出的定义：运筹学是一门在现有的技术及理论条件下，对问题现状的分析强调最优化决策的科学方法。

运筹帷幄之中，决胜千里之外这其中的运筹两字是赤壁之战的核心与关键，是整个战争通敌制胜的法宝。