A Computational Algebraic Geometry Based Global ... - CiteSeerX

A Computational Algebraic Geometry Based Global Optimization Technique to Address Economic Dispatch †Rajesh G. Kavasseri and Parthasarathi Nag ‡ †Department of Electrical and Computer Engineering North Dakota State University, Fargo, ND 58105 - 5285 email: [email protected] ‡Department of Mathematics, Black Hills State University, Spearfish, SD 57799 Abstract— In algebraic geometry, the concept of Gröbner basis allows a systematic study of the solution of a system of polynomial equations. This concept can be applied to find the global (and all local optima) optimum of a nonlinear, not necessarily convex function, the only restriction being that the objective function be polynomial. The method is based on computing a lexicographic (lex) ordered Gröbner basis for the ideal generated by the first order necessary conditions defined by the Lagrangian. Computing the optimal solution is then equivalent to computing the variety corresponding to this ideal. By virtue of the (lex) ordering, the system is transformed in to set of polynomials which can be solved successively to obtain the solutions. Here, we illustrate the application of the method on a non-convex function and identify the global optimum from the set of fifteen stationary points (6 local minima, 2 local maxima and 7 saddles). Then we apply the method to solve the classical economic dispatch problem including a combined cycle heat plant (CCHP) whose piecewise linear cost function is approximated by a smooth tenth order polynomial. Interestingly, the the method yields two possible solutions from which the least cost solution can be picked. While the work reported here is only preliminary, we find the results encouraging and hope that the method will find applicability in identifying the global optimum of non-convex power systems optimization problems.

I. I NTRODUCTION The problem of economic dispatch centered on minimizing the operating costs of a set of generators in a power system subject to constraints has a long and interesting history, a compendium of which is provided in [1], [2]. Depending upon the nature of the objective function and constraints, a variety of approaches including nonlinear, quadratic, linear, mixed-integer programming, Newton based and interior point methods have been proposed to solve the given optimization problem, [1], [2]. A class of problems known as polynomial programming problems (PPP) arises when the objective function and constraints are smooth polynomials. This class covers a wider range of non-convex functions in addition to representing the well known case of quadratic-convex optimization problems. For this class of problems, there are a few interesting approaches to find the global minimum such as the semi-definite programming (SDP) using the sum of

squares technique, [4], [5], homotopy continuation methods, [6], [7] and algebraic techniques, [3], [14]. While the SDP and related approaches can be used in the case of convex approximations to the original problem, the advantage of algebraic methods is that they can identify the exact global minimum at the expense of computational cost. Our motivation to explore the algebraic technique of Gröbner basis to find the global minimum of polynomial functions arises in the wake of recent interest in the economic dispatch involving combined cycle heat plants (CCHP). CCHP plants have nonquadratic cost functions and the resulting loss of convexity has led researchers to explore alternative heuristic methods, [15], hybrid and genetic algorithms [16], [17]. In this paper, we consider the problem of economic dispatch in a simple power system involving a CCHP and approximate its piecewise linear cost curve with a smooth tenth order polynomial. Then the algebraic technique proceeds by computing the Gröbner basis for the ideal generated by the necessary first order Lagrangian conditions. Because of the lex ordering, this set of equations is transformed into a set of polynomial equations, the zeroes of which can be successively computed (one variable at a time) in a fashion similar to Gaussian elimination method. Therefore, the method can find all the solutions to the Lagrangian system, from which the global minimum can be exactly found. This method is first illustrated on a test example with fifteen critical points ( 7 saddles, 6 local minima and 2 local maxima) and then applied to the economic dispatch problem. The rest of this paper is organized as follows. In Sec.II, a brief background on the theory of Gröbner bases is presented. In Sec. III, the proposed method is illustrated on a test system consisting of conventional thermal units (with quadratic cost functions) and a CCHP with a non-smooth piecewise linear cost function. A discussion of the results obtained with this method is also presented. Finally, the conclusions are collected in Sec. IV.

1-4244-1298-6/07/$25.00 ©2007 IEEE. Authorized licensed use limited to: Universidad Pablo de Olavide. Downloaded on May 15, 2009 at 13:45 from IEEE Xplore. Restrictions apply.

¨ II. A BRIEF OVERVIEW OF THE THEORY OF G R OBNER BASES A. Overview The concept of Gröbner bases (GB) was proposed in [10] to address fundamental problems in commutative algebra, especially the theory of polynomial ideals, [11]. In recent years, owing to computational advances, the concept has found several applications in areas including: systems theory (multivariate polynomial matrix factorization, Bezout identities, state space isomorphisms [13]), signal processing (design of multidimensional FIR and IIR filter banks) and control theory (H∞ theory, Lyapunov functions, [13]), [12]. One of its important applications arises in the context of solving an algebraic set of polynomial equations of the form: f1 (x1 , x2 , . . . xn ) = 0 .. .

(2)

fn (x1 , x2 , . . . xn ) = 0

(3)

(1)

where f1 , . . . fn are multivariate polynomials in n indeterminates (x1 , . . . xn ). The GB formalism provides a systematic technique to answer questions such as (i) does the system Eqn.(3) have solutions in the first place ? If so, (ii) what is the cardinality of the set of solutions ? (iii) when the solution set is finite, can we enumerate and determine all the solutions ? Incidentally, it can be noted that questions (i) and (ii) can be answered without having to solve the system. Loosely speaking, the key strategy is to reduce the original set of polynomials to an “equivalent” canonical set of the form: g1 (x1 ) = 0 g2 (x1 , x2 ) = 0 gm (x1 , x2 , . . . xn ) = 0

(4) (5) (6)

where g1 , . . . gm are polynomials in their respective arguments. Note that this canonical set has a structure that resembles the triangular form obtained through Gaussian elimination for linear systems. Because the two systems are equivalent, they possess the same solution set. However, note that it is easier to solve the canonical system because the process involves finding the roots of a polynomial in a single variable. For example, once all the roots of g1 (x1 ) = 0 are found, then these can be substituted in to g2 (x1 , x2 ) = 0 (which is essentially a polynomial in the single variable x2 ) to find the roots x2 and so on, until all the solutions have been determined. The computational process of arriving at this canonical set (called the Gröbner Basis) from the original set is accomplished through Buchberger’s algorithm, [10] which has been implemented in numerous symbolic computation packages such as MAPLE, Singular, Mathematica and CoCoA. In what follows, the relevant background and terminology from computational algebraic geometry are introduced to understand the concept of Gröbner Bases and Buchberger’s algorithm.

B. Background and Terminology Let x denote an n− tuple of indeterminates (x1 , x2 , . . . xn ), α denote an n− tuple of natural numbers (α1 , α2 , . . . αn ) α2 αn 1 and k denote a field. Then x α = xα 1 · x2 · · · xn is said to be a monomial. A polynomial f is then a combination of monomials such that f = rα xα where the coefficients rα α

belong to the field k. The set of all such polynomials denoted by k[x1 , . . . xn ] forms a commutative ring (R) which is referred to as a polynomial ring. A non-empty subset I of R = k[x1 , . . . xn ] which is closed under addition and multiplication is called an ideal. More precisely, a subset I ⊂ R is an ideal if (i) 0 ∈ I (ii) I is an additive subgroup of R. and (iii) r ∈ R and i ∈ I implies ri ∈ I. For example, if f1 and f2 are polynomials, then the set (I(f1 , f2 )) of all linear combinations of the form h1 f1 + h2 f2 where h1 , h2 are polynomials in R forms an ideal which is said to be generated by f1 and f2 and denoted by < f1 , f2 >. Alternatively, one can think of f1 , f2 as basis elements for the ideal that they generate. A formal way of characterizing the common zeroes of a polynomial system is obtained by introducing the idea of a affine algebraic variety. For example, the set of all common zeroes to the system {f1 , f2 } is called the variety corresponding to that system and denoted by V (f1 , f2 ). With the notion of an ideal, one can now note that V (f1 , f2 ) = V (I(f1 , f2 )). Generalizing this, if F = {f1 , . . . fn } is a finite (basis) set of polynomials, then < F > is the set (or the ideal) n generated by these basis elements such that I = {g| g = 1 fi hi , hi ∈ R}. A natural question is to ask if every ideal of the polynomial ring has a finite basis, which the Hilbert Basis Theorem answers in the affirmative and the ring R is said to be Nötherian, (see [8],[9] for details). Thus the problem of finding the zeroes to a system of polynomials F = {f1 , . . . fn } is equivalent to determining the variety for the ideal generated by that set or V (I(f1 , f2 , . . . fn )). Since every ideal of a polynomial ring has a finite basis, and there can be several bases for an ideal, the question now becomes is it possible by changing the basis such that it has a specialized structure (compared to the original system) say the “triangular form” or the canonical form mentioned earlier ? The answer is affirmative and this makes it easier to determine the affine algebraic variety. The basis so constructed is called the Gröbner Basis which can be obtained through an algorithmic computation due to B. Buchberger, [10]. This development hinges on two key ideas namely: (a) the concept of monomial orderings and (b) normal form reduction of a polynomial. In the following subsection, a brief introduction to these concepts is provided. For further expository details, one can refer to [3]. C. Monomial Orderings, Normal Forms and S-polynomials 1) Monomial Orderings: The concept of monomial ordering plays a fundamental role in the computation of a Gröbner Basis. Let Z+ n = {α = (α1 , α2 , . . . , αn ) | αi ≥ 0, αi ∈ Z, ∀i = 1 . . . n} denote the set of n tuples non-negative integers. Then > is an admissible total ordering on Z+ n if (i) α > 0 ∀α ∈ Z+ n and (ii) ∀α, β, γ ∈ Z+ n , α > β ⇒

Authorized licensed use limited to: Universidad Pablo de Olavide. Downloaded on May 15, 2009 at 13:45 from IEEE Xplore. Restrictions apply.

α + γ > β + γ. Using this notion of ordering, monomials can be ordered in different ways, one of the more popular ones being the lexicographic order which is used in this paper. Suppose α, β ∈ Z+ n , then: •

lexicographic order (α >lex β) ⇒ the left most nonzero entry in α − β is positive. Eg: If α = (4, 3, 5), β = (1, 8, 9), then, α >lex β because α − β = (3, −5, −4) and therefore, the monomials are ordered as: x4y 3 z 5 >lex xy 8 z 9 .

Once the monomials within a polynomial are ordered, then the polynomial can be rearranged to reflect that ordering. For example, if f1 (x, y, z) = 7xyz 3 − 3x2 y 3 + 2xy 3 z, then with respect to the lex ordering x > y > z, then the polynomial is reordered in descending order as f1 (x, y, z) = −3x2 y 3 + 2xy 3 z + 7xyz 3 . After such an ordering, one can define the following terms: (i) the multidegree of the polynomial f as the max (α ∈ Z+ n ) that appears in f . For example, multideg(f1 ) = (2,3,0), (ii) the leading monomial LM (f ) of f as xmultideg( f ) . In this case, LM (f1 ) = x2 y 3 , (iii) the leading term (LT (f )) as the monomial rmultideg( f ) xmultideg( f ) . For example, LT (f1 ) = −3x2 y 3 . In Gaussian elimination, variables are eliminated successively in a certain order to arrive at the reduced triangular form. In the case of polynomial systems, monomial ordering serves to identify the sequence in which the leading terms of a polynomial are eliminated. This aspect will become clear when the concept of normal forms and S-polynomials are introduced, which is the second and most important step in computing the GB. 2) Normal Forms and S-polynomials: Let F be a finite collection of polynomials. A polynomial p is said to reduce to another polynomial q modulo F if the leading term LT (p) of polynomial p can be deleted from p by subtracting a certain multiple of a polynomial f from F to obtain q. This can be LT ( p) expressed as q = p − [f LT ( f ) ]. This is possible only if the leading term of p is divisible by the leading term of some polynomial in the set F . In the event this is not possible, i.e. when LT (p) is not divisible by any of LT (f ), f ∈ F , then p is said to be irreducible modulo F . Further, if p reduces to q modulo F and q is irreducible modulo F , then q is said to be the normal form of p, denoted by q = normalf orm(p). Example: Let F = {f1 , f2 } with f1 = x2 y 2 + x2 y + y and f2 = x − y 2 . Suppose p = x6 y 2 + x3 y 4 + y 2 + 7. Choose f = f1 so that LT (f ) = x2 y 2 and LT (p) = x6 y 2 . LT ( p) Then by the computation q = p − [f LT ( f ) ], yields q = −x6 y − x4y + x3 y 4 + y 2 + 7. Note that q can be further reduced with respect to F by choosing f = f2 . This reduction process can be continued until the resultant polynomial q is no longer reducible modulo F in which case we say q is the normal form of p. Also recall that the terms in all the polynomials are ordered with the lex ordering.

With these concepts, it is possible to formulate a definition for a Gröbner Basis. The following is a definition proposed by Buchberger, [10]. Suppose F is a collection of polynomials, then F is a Gröbner Basis iff all normal forms modulo F are unique. That is to say, if g1 and g2 are normal forms of g modulo F , then g1 = g2 . In order to do compute the Gröbner Basis, we introduce the concept of S-polynomials. These S-polynomials are related to computation of syzygies on the elements of the leading terms of the given polynomials. It is to be noted that the “S” in Spolynomials seems to have originated from syzygy. Given two polynomials f and g, the S-polynomial of f and g denoted by spoly(f, g) is defined as: spoly(f, g) = LCM (LM (f ), LM (g)) ×

f g − LT (f ) LT (g)

The least common multiple or LCM of two monomials x α2 αn β 1 = xβ1 1 · xβ2 2 · · · xβnn is given by xα 1 · x2 · · · xn and x LCM (x

α

,x

β

(xam

) = x1

α1 ,β1 )

max(

x2

α2 ,β2 )

· · · x(xamn

αn ,βn )

α

=

(8)

In essence, spoly(f, g) is simply a linear combination of the polynomials f and g chosen such that the leading terms cancel. Example: Let f1 = x4y 2 z 3 + xy and f2 = x3 y 4z − x2 y 2 z 3 . Then LM (f1 ) = x4y 2 z 3 and LM (f2 ) = x3 y 4z and LCM (LM (f1 ), LM (f2 )) = x4x4z 3 . Then from Eqn.(7) yields, spoly(f1 , f2 ) = −x3 y 2 z 5 + xy 2 . D. Buchbergers Algorithm With the concept of S-polynomials, an alternative definition for a GB can be formulated as follows, [10]. Theorem: Suppose G = {g1 , g2 , . . . gn } be a finite set of polynomials and I be the ideal generated by G. Then, the following are equivalent (i) G is a Gröbner Basis, (ii) ∀gi , gj ∈ G, spoly(gi , gj ) reduces to zero modulo G. Criterion (ii) in the above theorem serves as the fulcrum in the algorithm to compute the GB. The basic algorithm to compute the GB as first proposed by Buchberger, [10] is outlined below. F = {f1 , . . . fn } ← input: original set of polynomials that generates an ideal I G = {g1 , . . . gm } → output: set of polynomials that generate the same ideal I with F ⊂ G G ← F (initialize G to F ) R ← {(fi , fj )|fi , fj ∈ F (i = j)} while R = ∅ {(f, g)} ← pick a pair in R R = R\{(f, g)} h = normalf orm(spoly(f, g), G)


(7)

if h = 0 G ← G ∪ {h} R ← R ∪ {(f, h)}, ∀f ∈ G end end return(G)

2 1.5 1

f(x,y)

0.5

The algorithm first initializes the GB to the original set F and forms a set R of all distinct pairs of polynomials in F . For a pair (f, g) from the set R, the spoly(f, g) is computed and reduced to h modulo G. If h is non-zero, it is appended to the set G and the pair set R is updated to include new pairs {h, g} for every g in F and the process continues till h = 0. The algorithm above does not yield a necessarily unique GB and it not optimal in the sense that there could be redundant polynomials in the set GB. In the modified Buchbergers algorithm, this is remedied and a reduced set such that g = normalf orm(g, G\{g}) is obtained. The resultant set is called the reduced GB (RGB), which has been implemented in several symbolic computation packages. While a brief overview of this topic was presented here, the reader is referred to [3], [11], [12] for further details. In the next section, we apply the GB technique to solve the system of polynomials arising in the solution of the economic dispatch problem for a system of generators. The RGB computations in this paper were done in MAPLE-10 by invoking the basis command with the package groebner. III. R ESULTS The results are organized in three subsections. First, we illustrate the application of the technique to determine all the stationary points of a function [4] known to have several local minima. Second we consider the classic economic dispatch problem (lossy case with quadratic cost functions) and finally consider the optimal dispatch with a CCHP and conventional steam generator. A discussion of the results obtained is also presented. A. Illustrative Example We shall consider the following global optimization problem (due to N. Z. Shor [4]) of determining all the local and global minima of the function described below 21 4 1 6 x + x + xy − 4y 2 + 4y 4. (9) 10 3 This is an unconstrained nonlinear non-convex optimization problem which is solved by the following steps: F (x, y) = 4x2 −

Step 1: Compute the critical points by computing ∇F = (

∂F ∂F , )=0 ∂x ∂y

42 3 x + 2x5 + y = 0 5 f2 (x, y) = x − 8y + 16y 3 = 0

−1 −1.5 −2 1 0.5

1 0.5

0 0

−0.5

−0.5 −1

y

Fig. 1.

−1 x

A plot of the function described in Eqn.(2).

respectively. Step 2: We compute the Gröbner basis of the ideal generated by the above equations using MAPLE-10 by choosing a lexicographic order x lex y. The Gröbner basis is as follows: 10485760y 15 + 26214400y 11 − 456704y 5 +3534848y 7 + 22144y 3 − 26214400y 13

g1 (x, y)

=

g2 (x, y)

−325y − 13279232y 9 = x − 8y + 16y 3

(13) (14)

Step 3: We compute the affine algebraic variety V( g1 , g2 ) of the ideal generated by Gröbner basis (3.3a) and (3.3b) as listed in the second column of Table I. TABLE I S UMMARY OF CRITICAL POINTS FOR THE FUNCTION point 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

(x, y) (1.703606716,-.7960835687) (1.109205337,-.7682680925) (0.089842013,-.7126564030) (-1.296070267,-.6050843880) (-1.607104753,-.5686514549) (-1.638067984,-.2286740690) (-1.230229877,-.1623345845) (0,0) (1.230229877,.1623345845) (1.638067984,.2286740690) (1.607104753,.5686514549) (1.296070267,.6050843880) (-0.089842013,.7126564030) (-1.109205337,.7682680925) (-1.703606716,.7960835687)

∆(H) 427.10 -160.9505872 126.7048724 -59.53313938 71.36896398 -68.97094291 47.71854634 -65 47.71854634 -68.97094291 71.36896398 -59.53313938 126.7048724 -160.9505872 427.1022892

Fx 19.09 -7.86720077 7.797247505 -6.11369524 9.62161975 12.38088158 -7.23355214 8.0 -7.23355214 12.38088158 9.62161975 -6.11369524 7.797247505 -7.86720077 19.09469944

(10)

which yields, f1 (x, y) = 8x −

0 −0.5

(11) (12)

Step 4:⎡We compute the ⎤ determinant of the Hessian matrix ∂2 F ∂2 F ⎢ ∂x2 ∂x∂y ⎥ ⎥ at the respective critical points listed H = ⎢ ⎣ ∂2 F ∂2 F ⎦ ∂x∂y

∂y 2


in the set V( g1 , g2 ) so as characterize their properties i.e. maximum, minimum, saddle and degenerate critical points. The determinant of the Hessian matrix ∆(H) and the second second partial derivative with respect to x (Fxx ) are shown in the third and fourth columns of Table I respectively. ∂2 F ∂2 F > 0 or > 0 then ∂x2 ∂y 2 the corresponding critical point is a local minimum. Also, if ∂2 F ∂2 F the determinant is positive and < 0 or < 0 then ∂x2 ∂y 2 the corresponding critical point is a local maximum and if the determinant is negative then the critical point is a saddle point.

If the determinant is positive and

From the above classification scheme, one can note that there are seven saddle points, six local minima and two local maxima. Step 5: Evaluating the function at each of the 15 points yields us the global maximum value of 2.496295351 at the points # 7 and 9 while the global minimum of −1.031628453 is attained at points # 3 and 13. This completes the procedure of computing and identifying the global maxima and minima using Gröbner basis. B. Classic ED problem: lossy case with quadratic cost functions We consider a system of generators with quadratic cost functions, [18] with the total cost given by: n

C(P1 , . . . Pn ) = (β0 i + βi Pi + γi Pi2 ), (15)

The polynomial system (Eqns.(19-20)) are solved by computing the GB for the ideal generated by them. Here, we choose the lex ordering P1 lex P2 lex . . . Pn lex λ to obtain the GB. For example, with n = 3 (three generator case), the GB obtained is as follows: a16 λ6 + a13 λ3 + a12 λ2 + a11 λ − a10 + a14λ4 + a15 λ5

(21)

a2 P3 + a23 λ3 + a22 λ2 + a21 λ − a20 + a24λ4 + a25 λ5 (22) 4 3 2 a3 P2 + a4 + a53 λ5 (23) 3 λ + a3 λ + a23 λ + a13 λ − a03 a4P1 + a45 λ5 + a43 λ3 + a42λ2 a41λ − a40 + a44λ4, (24) where the coefficients, a1 , a10 , a11 , a12 , a13 , a14, a15 , a2 , a20 , a21 , a22 , a23 , a24, a25 , , a45 a3 , a03 , a13 , a23 , a3 , a4 3 , a53 , a4, a40 , a41, a42, a43 , a44 are given in the Appendix I. Notice that the RGB for this system has a structure that essentially decouples the variables. In this case, Eqn.(21) is a polynomial strictly in λ, and Eqns.(22-24) are of the form Pi + g(λ) = 0, i = 1, 3 where g is a polynomial in λ. The optimal multiplier λ is then found by solving for the roots of the polynomial systems Eqn.(21). While Eqn.(21) is a 6th order polynomial with 6 roots, there is only one real root at λ = 7.67, corresponding to a unique global minimum. Therefore iterations are required only once in determining the roots of the first polynomial Eqn.(21). Substituting for λ in Eqns. (22-24) then yields the optimal dispatch values as tabulated in Table II. P1 , P2 , P3 as 35.1, 64.12 and 52.47 MW respectively.

i=1

which is to be minimized subject to the constraint condition n

Pi = PD + P L (16) i=1

where Pi is the real power outputs of generator i and PD is the demand and PL is a quadratic loss function given by PL = P t BP . The loss coefficients and system data are supplied in Appendix II.

The method was also applied to two other test cases with 6 and 9 generators. The test case parameters were adapted from the 26 bus power system network in [18]. A summary of the dispatch results obtained for all cases is presented in Tables II, III and IV. From the results, it is clear that solutions obtained with the proposed method are in close agreement to those obtained by the λ -iteration method. The coefficients of the RGB polynomials for these cases has been omitted from the appendices due to space limitations.

Following the classical Lagrangian approach, we have: n

L = C(P1 , . . . Pn ) + λ(PD + PL − Pi ) (17)

TABLE II E CONOMIC D ISPATCH (n = 3, PD unit 1 2 3 λ

i=1

and the necessary first order conditions for a stationary point is given by ∂L ∂L ∂L )=0 ,... , ∂P1 ∂Pn ∂λ which amounts to solving the polynomial system: ∇L = (

n

βi + 2γi Pi + 2λ( Bij Pj ) + B0 i λ − λ = 0 n

i=1

Pi = 0

Pi(λ-iter) 33.47 64.09 55.10 7.76

(18) C. CCHP - Thermal generator dispatch (19)

j=1

PD + PL −

Pi(GB) 35.1 64.12 52.47 7.67

= 150 MW)

(20)

For simplicity, we consider a two generator system (lossless case) one of which is CCHP whose cost function [17] is a continuous, nondecreasing, piecewise linear function which is approximated by a tenth order polynomial cost C1 = 10 ai P1i . The functional approximation is shown below. i=0


TABLE III E CONOMIC D ISPATCH (n = 6, PD Pi(GB) 454.68 184.49 269.37 138.08 157.69 82.60 13.5754

unit 1 2 3 4 5 6 λ

Following the Lagrangian formulation,

= 1276 MW)

L = C(P1 , P2 ) + λ(PD − (P1 + P2 )).

Pi(λ-iter) 452.59 173.19 267.79 136.27 153.68 80.57 13.5429

the necessary conditions for a stationary point yield, 10

∂L = iai P1i−1 − λ = 0 ∂P1 i=1

Pi(GB) 349.39 96.98 187.23 52.57 75.16 214.80 101.55 132.67 103.84 12.0006

unit 1 2 3 4 5 6 7 8 9 λ

= 1300 MW)

Pi(λ-iter) 350.8 97.91 188.18 53.77 74.20 215.82 103.44 128.17 104.99 12.0200

h1 =

10

9

γi λi = 0

(30)

h2 = 35 + 8P2 − 5λ

(31)

h3 = 5λ + 8P1 − 47

(32)

i=0

The other is a regular generator with the quadratic cost function given by C2 = b0 + b1 P2 + b2 P22 . Thus the resultant cost function is given by C(P1 , P2 ) =

(27)

∂L = b1 + 2b2 P2 − λ = 0 (28) ∂P2 ∂L = PD − (P1 + P2 ) = 0. (29) ∂λ Computing the Gröbner basis for the above set of equations with the lex order P1 lex P2 lex λ. yields,

TABLE IV E CONOMIC D ISPATCH (n = 9, PD

(26)

ai P1i + b0 + b1 P2 + b2 P22 ,

(25)

i=0

subject to the constraint condition P1 + P2 = PD where PD is the demand of the load and is equal to 150M W and P1 , P2 are the generator real power outputs.

where the value of the coefficients γi , i = 0 . . . 9 are given in the Appendix II. Solving Eqn.(30) yields three possible (real) solutions for λ, two of which are positive while the third is negative. Therefore, we get two possible feasible solutions (third solution is infeasible) to the above system of equations, which are {(P1 , P2 , λ) V( h1 , h2 , h3 ) R3 = (0.694, 0.805, 8.28), (1.0774, 0.4225, 7.6761)} Evaluating the cost function at these two points yields the former of the two as the least cost solution. In this case, the least cost solution corresponds to the region when the CCHP is operating in the flat region (63, 82) MW. D. Discussion

Cost

1.6

y = 58*x10 − 7.5e+002*x9 + 4.2e+003*x8 − 1.4e+004*x7 + 2.9e+004*x6 − 4e+004*x5 + 3.8e+004*x4 − 2.4e+004*x3 + 9.4e+003*x2 − 2.1e+003*x + 2.1e+002

1.4

1.2

1 0.5

1

1.5

2

P (p.u) 1

Fig. 2. Polynomial Approximation (dashed line) of the CCHP cost function (solid line). A base of 100 MVA is assumed.

Here, a brief discussion on some of the unique aspects of the GB method and some the challenges associated with it are presented. • The illustrative example (section III. A) conveys two noteworthy features. First, it is evident that the degree (15) of the polynomial in the RGB Eqn.(13) can be much greater than the maximum degree (6) that occurs within the polynomial Eqn.(2). However, what is gained at this expense is a provision to find all critical points where the gradient ∇F vanishes. In this case, it is possible to find all the critical points simply by solving for all real roots of Eqn.(13) and then substituting them in Eqn.(14). In fact, the ability to extract all feasible solutions is indeed a strength of the GB method. Second, notice that the GB based solution process differs from conventional gradient based numerical optimization algorithms. While the former is global, the latter is essentially local and progresses from one candidate solution point to another


iteratively based on the gradient information at that point. •

An inspection of the coefficients obtained for the Gröbner basis polynomials in Appendices I and II reveals that these coefficients can be considerably large numbers. Moreover, in the intermediate step of calculating the S-polynomial, the coefficients tend to be even larger integers. This is an inherent limitation arising from the application of Buchberger’s algorithm which could potentially restrict the application of this technique to large scale problems at present. For the problems considered in this paper, no notable difference in execution time was found between the GB and λ-iteration method. It has been noted in the literature that the choice of ordering has a significant impact on the computation time. Therefore for larger scale problems, a way to reduce execution time could be by choosing an ordering different from lex which requires the most computation time. While there are n ! possible choices of ordering for an n variable problem, the degrevelex ordering has been found to provide some relief in terms of computation time for most problems. However, further exploratory work is required to assess the benefits of different ordering schemes for power system optimization problems.

Another important issue (namely that of scalability) with the proposed method is to assess the growth in computational complexity with the number of variables at hand. Here, the number of generators (all with quadratic cost functions) was varied (from three to nine) and the complexity assessed by noting the highest degree (d) of the polynomial that arises in the GB along with the number of solutions generated. The highest degree of the polynomial d was found to vary with the number of variables n as d = 2(n − 1). With the exception of the three generator case, the number of real solutions to the system was two in all cases with one of the multipliers λ being positive and the other, negative. In these cases, it is easy to pick the feasible solution corresponding to λ > 0. In all cases, we consistently observe a unique structure for the GB in the form Pi + g(λ) = 0 as mentioned earlier. By virtue of this, there is only one optimal solution corresponding to a given λ. Also, the examples considered here involved ten variables utmost, one equality constraint and no inequality constraint, which are atypical from the standpoint of large scale power systems optimization problems. However, the work reported here is only preliminary and an extension of this technique to a realistic power system problem is bound to be challenging and exciting. Currently, the authors are actively exploring methods to circumvent some of the limitations noted here. •

IV. C ONCLUSIONS An algebraic technique based on the concept of Gröbner basis is presented to solve the problem of economic dispatch.

The principal advantage of this method is its ability to find the global minimum of a nonlinear non-convex optimization problems in a non-iterative fashion, provided the objective function and constraints are polynomial. On the other hand, some of the limitations of the method include (i) the large size of the numbers involved in the computations and (ii) the computational burden involved with S polynomial calculations and (iii) the highest degree of the polynomial that needs to be solved in the GB. While improvements to Buchberger’s basic algorithm have been noted to provide some relief in computational time and burden, the degree and coefficients in the final GB polynomials remain unchanged. This issue can pose computational challenges, when applied to a large scale optimization problem with several variables and constraints. Nevertheless, it is anticipated that advances in computational methods and algorithms can help exploit the potency of the GB method as a global optimization tool. Further research is required to illustrate the applicability of the technique in solving power system optimization problems.

A PPENDIX I: C OEFFICIENTS OF POLYNOMIALS IN E QNS .(21) -(24) a16 = 28776237562121 a13 = 23226607885552575000000000 a12 = 39831050641551875000000000000 a11 = 5768329999125000000000000000000 a14 = 5759380162234654625000 a10 = 40618752243750000000000000000000000 a15 = 663448111612650000 a2 = 850865381250000000000000000000, a23 = 20675682795937731677292875000 a22 = 52357844584742850260187500000000 a21 = 23693230617366957466875000000000000 a20 = 59571787558096893562500000000000000000 a24 = 3159686229598841165650000 a25 = 165041477563297494019 a3 a13 a3 a23 a13 a03 a53

= 1347865136718750000000000000 = 3771473965556740275000 = 24616191300306679197875000 = 62142127301584009828125000000 = 27903305423099008671875000000000 = 70579680764575271484375000000000000 = 197376213438587939

a4 = 10435980080625000000000000000000 a41 = 37266062153835996509 a43 = 4807451109574014924359125000 a42 = 12437108362964165009825000000000 a41 = 5941601046547402359375000000000000 a40 = 14334125264693081625000000000000000000 a44= 722428353143394676450000


A PPENDIX II: C OEFFICIENTS OF POLYNOMIALS IN E QN .(30) γ9 γ6 γ2 γ5 γ1

= 283203125, γ8 = 18685546875, γ7 = 546295312500 = 9279508437500, γ4 = 723097388056250 = 10143026055451500, γ0 = 11791046329231715 = 100753369593750, γ3 = 3414103872902500 = 16953467763524957

b0 = 20, b1 = 7, b2 = 10 (in Eqn.25) The loss coefficients B (in pu on a 100 MVA base) for the three generator case (n = 3) [18] ⎡

0.0218 B = ⎣ 0.0093 0.0028

0.0093 0.0228 0.0031

⎤ 0.0028 0.0017 ⎦ 0.0015

A PPENDIX III: G ENERATOR COST FUNCTION COEFFICIENTS TABLE V unit 1 2 3 4 5 6 7 8 9

α 240 200 220 200 190 190 170 200 250

β 7 10 8.5 11 10 12 10 9 10

γ 0.007 0.0095 0.009 0.009 0.009 0.0075 0.009 0.007 0.008

[10] B. Buchberger, “An Algorithm for Finding a Basis for the Residue Class Ring of a Zero-Dimensional Polynomial Ideal (in German), Doctoral Dissertation, Math. Inst. University of Innsbruck, Austria, 1965. [11] B. Buchberger, “Gröbner Bases: A Short Introduction for Systems Theorists”, Proc. EUROCAST 2001 (8th International Conference on Computer Aided Systems Theory - Formal Methods and Tools for Computer Science), , Las Palmas de Gran Canaria, Lecture Notes in Computer Science pp: 2178-2201, Feb. 19-23, 2001. [12] B. Buchberger, “Gröbner Bases and Systems Theory” Multidimensional Systems and Signal Processing, vo.12, no.(3-4), Springer Netherlands, pp: 223-251, July, 2001. [13] Computer Simplification of Formulas in Linear Systems Theory, “J. W. Helton, M. Stankus and J. J. Wavrik, IEEE Trans. Automatic Control, 43(3), pp: 302-314, March 1998. [14] P. A. Parrilo and B. Sturmfels, “Minimizing polynomial functions”, Workshop on Algorithmic and Quantitative aspects of real algebraic geometry in mathematics and computer science, DIMACS, Rutgers, March 2001. [15] J. R. Santos, A. T. Lora, A. G. Exposito and J. L. M. Ramos, “Finding Improved Local Minima of Power System Optimization Problems by interior point methods”, IEEE Trans. on Power Systems, 18(1), pp: 238244, February 2003. [16] F. Gao and G. B. Sheble, “Economic Dispatch Algorithms for thermal unit system involvong combined cycle units”, 15 PSCC, Liege, Belgium, Sec.21, paper 5, August 2005. [17] R Gnanadass, P. Venkatesh, T. G. Palanivelu and K. Manivannan, “Evolutionary Programming Solution of Economic Load Dispatch with Combined Cycle Co-generation Effect”, IE(I) Journal -EL, vol. 85, pp: 124-128, Spetember 2004. [18] H. Saadat, “Power System Analysis”, McGraw Hill, New York, 1999. [19] A. J. Wood and B. F. Wollenberg, “Power Generation, Operation and Control (2nd Ed), 1996.

ACKNOWLEDGMENT The authors would like to thank the reviewers for their constructive suggestions which have helped in substantially improving the quality of this manuscript. R EFERENCES [1] J. A. Mohmoh, M. E. El Hawary and R. Adapa, “A review of selected optimal power flow literature to 1993 Part I: Nonlinear and Quadratic Programming Approaches”, IEEE Trans. on Power Systems, 14(1), pp. 96-104, February 1999. [2] J. A. Mohmoh, M. E. El Hawary and R. Adapa, “A review od selected optimal power flow literature to 1993 Part II: Newton, Linear Programming and Interior Point Methods”, IEEE Trans. on Power Systems, 14(1), pp. 105-111, February 1999. [3] D. Cox, J. Little and D. O’Shea, “Ideals, Varieties, and Algorithms: An Introduction to Algebraic Geometry and Commutative Algebra”, 2nd ed. New York: Springer-Verlag, 1996. [4] N. Z. Shor, “Class of global minimum bounds of polynomial functions, Cybernetics, 23(6), pp: 731-734, 1987. [5] N. Z. Shor and P. I. Stetsyuk, “The use of a modification of the r algorithm for finding the global minimum of polynomial functions”, Cybernetics and Systems Analysis, 33, pp: 482-497, 1997. [6] T. Y. Li, “Numerical solution of multivariate polynomial systems by homotopy continuation methods”, Acta Numerica, 6, pp: 399-436, 1997. [7] J. Verschelde, “Polynomial homotopies for dense, sparse and determinantal systems”, Mathematical Systems Research Institute (MSRI) preprint, # 1999-041. [8] David Eisenbud, “Commutative Algebra with a View toward Algebraic Geometry”, Springer-Verlag, 1999. [9] M. Reid, “Undergraduate Algebraic Geometry”, London Mathematical Society Series, No.12, Dec 1998.
