Coalition formation feasibility and Nash–Cournot equilibrium problems in electricity markets: A Fuzzy ASA approach

Applied Soft Computing 35 (2015) 1–12

Contents lists available at ScienceDirect

Applied Soft Computing journal homepage: www.elsevier.com/locate/asoc

Coalition formation feasibility and Nash–Cournot equilibrium problems in electricity markets: A Fuzzy ASA approach Hime Aguiar e Oliveira Jr. ∗ National Cinema Agency, Rio de Janeiro, Brazil

a r t i c l e

i n f o

Article history: Received 8 January 2015 Received in revised form 6 June 2015 Accepted 9 June 2015 Available online 22 June 2015 Keywords: Adaptive simulated annealing Coalition formation Generalized Nash equilibrium problems Global optimization Metaheuristics

a b s t r a c t A numerical method based on a stochastic global optimization paradigm is presented and applied to the calculation of Nash–Cournot equilibria in electricity markets. The proposed method uses and solves GNEP’s (generalized Nash equilibrium problems) by means of the Fuzzy Adaptive Simulated Annealing (Fuzzy ASA) algorithm. Concepts of cooperative game theory are used too, such as the bilateral Shapley value. This approach makes it possible to study the feasibility of coalition formation processes in several scenarios, and a case study based on the IEEE 30-bus system is used in order to present and discuss in detail the proposal. The main advantage of the new technique is to simplify and extend current methods, as explained along the text. © 2015 Elsevier B.V. All rights reserved.

1. Introduction During the last decades many countries restructured their electric energy markets, transforming traditional vertically integrated monopolies into deregulated entities, typically oligopolies. Therefore, competition has been incentivized by recently created electricity markets, where buyers and sellers can trade electricity in auctions or through mutual agreements. Cost optimization techniques used by electric utilities in the past are being replaced by efficient bidding methods. Normally, the target of the electric utilities is profit maximization, where prices are determined by suppliers, consumers, transmission line operators, and other agents. Ideal competitive markets are very hard to attain in the real world industry, mainly because of the tiny number of competing players. On the other hand, network constraints affect the competitiveness of the market, since market bidders produce bottlenecks that may induce a substantial increase in unitary prices [1]. In this way, it is sensible to assume that we are dealing with an imperfect market. Oligopolistic market models are being used to study electricity markets since the beginning of the restructuring “wave” that hit the world years ago. The model used here is similar to the ones in [2] and [1], considering that both papers deal with Nash–Cournot equilibria in electricity markets. However, there are some differences between them. While the method in [2] depends on the

∗ Tel.: +55 21 3037 6078. http://dx.doi.org/10.1016/j.asoc.2015.06.013 1568-4946/© 2015 Elsevier B.V. All rights reserved.

existence of a solution to a system of equations and inequalities, which result from mixed complementarity (KKT) conditions, the one in [1] relies on a deterministic function minimization procedure. Another significant difference is in the kind of equilibrium that each study is trying to establish: in [2] the target is to compute a Nash–Cournot equilibrium that would also satisfy a market clearing condition. On the other hand, this work and [1] look for coupled constraint equilibria, which is a relatively new solution concept to game theory problems, in which the strategy space is jointly restricted for all players. Therefore, both works consider the latter kind of equilibrium an adequate solution concept for many relevant electricity market problems. Using the approach of [2] or an analytical solution to a particular concave game of several players with nonlinear profit functions and, possibly, constraints, might be a tough task. The Nikaidô–Isoda function and a relaxation algorithm are combined in [3,4] to create a numerical method for finding solutions of infinite games. Such a method is attractive in that the most complex computational procedure required is minimization of a multivariate function. A sequential improvement of the Nikaidô–Isoda function is obtained through a relaxation algorithm that is shown to converge to a Nash equilibrium for a wide class of problems, including nondifferential payoffs and coupled constraint games [5,4]. The ability of handling games with constrained strategy spaces is of particular importance for electric energy market modeling, taking into account that in typical problems of electricity generation and distribution, the strategy space of competing agents is coupled. This is so mainly because of the capacity constraints and Kirchhoff laws, meaning that, in this

2

H.A. e Oliveira Jr. / Applied Soft Computing 35 (2015) 1–12

Table 1 Limits for power generation and cost coefficients for generating units. g

Company

Generator

Pmax (MW)

1

1

80

2

80

C Pmax (MW)

80

3

50

4 5 6

55 30 40

125

kind of problem, there are joint constraints imposed on the combined strategy space of all agents. Therefore, the set of strategies available to a specific agent depends on the choices of other agents and, if all agents act at the same time, no traditional noncooperative game theory concept can be used to solve such a game. However, the normalized equilibrium of Rosen, designed to solve games subject to coupled constraint sets (called, for short, coupled constraint games) can be applied [5]. Nowadays, this kind of problem is also addressed by the underlying techniques associated to GNEP’s. The main contribution of this paper is to deepen the discussion about coalition formation in certain scenarios – this is done by studying in detail a case study presented in [1], in the context of electricity markets. Besides, the task is handled by solving several GNEP’s, faced as global optimization problems, with the help of the Fuzzy Adaptive Simulated Annealing algorithm.

2

0.035

1.75

0.125

1

0.0166 0.05 0.05

3.25 3 3

(1)

That is, the several generation costs are polynomial functions of delivered power, with constant coefficients (ci , di , ei ), also given in Table 1, except for the ei , that are zero for all generators in this case g C study. The quantities Pmin and Pmin are also null for all generators and companies, respectively. As expected, individual firms want to maximize their profits, although this could mean to cooperate with other market players. Hence, given certain market conditions (elasticity values and price × demand curves), they need to find stable joint configurations (Nash–Cournot equilibria) and compare respective isolated gains, in order to decide what is the best alternative (coalition) to take. Assuming also that the electricity demand is a strictly decreasing function of the price, the demand function in an interval of time during a day of study considered standard can be expressed as 0 (p) + ap, in which P 0 (p) is the total power demand Pload (p) = Pload load level expected for a selected time interval, and a represents the elasticity of the demand with respect to price. In particular, the relationship between price and overall production of the IEEE 30bus system in a selected interval is modeled, in the original work [1], by Pload (p) = 189.2 − 0.5p,

Pload (p) =

ng 

Pgi − Ploss

(4)

i=1

and ng is the total number of generators, Ploss representing the transmission losses throughout the system. Losses are not considered in this study or in [1], but they can be incorporated without interfering with the proposed paradigms, mainly because its effect would be merely to change the form of objective functions. After establishing these premises, the profit j obtained by company j, that owns ngj generating units, can be found as follows ng 

j = (378.4 − 2

S(j)+ng j

Pgi )



S(j)+ng j

Pgk −

k=S(j)+1



2 ((ck /2)Pgk + dk Pgk + ek ),

k=S(j)+1

(5)

In [1], two important scenarios are introduced and analyzed – the first one considers an electricity market that uses the IEEE 30bus system [6] and will be our object of study in the present article. In that case study, it is assumed that there are three generating companies and each of them possesses a number of generating units, as shown in Table 1. Pg and PC are the power generation of a unit and a company, respectively. Using the same notation, it is also assumed that the generation cost of unit i is given by the expression 2 Ci (Pgi ) = (ci /2)Pgi + di Pgi + ei .

in which

i=1

2. The problem under study

(2) 0 (p) Pload

= 189.2. that is to say, a = −0.5 and Putting in a more convenient form, we have p = 378.4 − 2Pload (p),

di ($/MW h)

0.04

130

2

3

ci ($/MW2 h)

(3)

in which the variables corresponding to produced power are subgj gj ject to Pmin ≤ Pgj ≤ Pmax , as shown in Table 1. Also, the terms S(j) are used to make the expressions more compact, and given by S(1) = 0, S(2) = 1, and S(3) = 3. In [1], a specific software package was used to solve the original example, but the same results could have been achieved by applying the traditional Nash–Cournot equilibrium conditions expressed as a system of equations. The results were obtained assuming that the three companies compete against each other. However, some of the generating companies may consider the formation of a cartel, for the sake of increasing their profits. That is to say, it is possible that the final productions of the companies are not necessarily the ones resulting from the shown results, but lower, if they decide to cooperate between themselves. When the number of involved companies is small, it is feasible to study all their possible coalitions. In this fashion, after enumerating all combinations and giving rise to several games in which the players are the coalitions, the proposed algorithm can be applied to each one, computing Nash equilibria corresponding to the several configurations. In this way, cooperation among agents creates several scenarios that have different Nash equilibria. On the other hand, the final assignment of profits (resulting from cooperation) to individual companies is not found by the algorithm. To that end, cooperative game theory concepts, such as the bilateral Shapley value (BSV) [7] and the kernel [8] are useful tools to study how each Nash equilibrium value can be split among the companies participating in each coalition. Both approaches represent a feasible way to handle coalitions among players and to allocate profits after the coalitions are formed. In the real world, seeking profit optimization, companies can join others and become part of a new player, composed of two or more firms. In [1], some scenarios in which the companies can be arranged according to all possible coalition combinations were studied. The coalition values found express the profit obtained by each particular coalition, composed of agents (companies) arranged in different configurations. Such values are obtained as a result of finding Nash

H.A. e Oliveira Jr. / Applied Soft Computing 35 (2015) 1–12

3

Fig. 1. Demand× price curves with several elasticities.

equilibria in a previously cited game. In the previous work [1], the values are calculated by applying a relaxation algorithm to the Nikaidô–Isoda function corresponding to each scenario, and then obtaining the final individual profits per coalition after the convergence of the corresponding iterative algorithm. The adopted approach (also used here) at first finds all possible groupings of individual players, giving rise to several coalition structures, and, after that, faces each coalition as an individual player, in order to make it possible to find stable profit and production levels for each case. Therefore, in order to determine the coalitions that are actually formed, it is needed to obtain the values corresponding to each game, noting that they always correspond to the minimal profits that a coalition can guarantee for itself against any other coalition. Depending on the elasticity values, it may happen that no equilibrium values can be found, not to mention the profit assignments of individual companies. However, when the demand becomes more elastic, bilateral Shapley value solutions do exist. Fig. 1 depicts the changes in elasticity when pivoting around the price intercept. For the different elasticity coefficients, the games can get subadditive, superadditive, or none of them, since sometimes it is better to be alone and sometimes to join other agents, depending on the coalition values and/or joint profits. It is important to highlight that in [1] and here the experiments were done by keeping fixed the intercept with the unitary price axis, varying only the slope of the price–demand line by changing elasticity values. Thus, when the elasticity is smaller (more negative), market prices decrease in a slower pace as a function of increasing overall demand, and vice-versa. Considering that the Shapley value is used in this study, it is worth to mention some of its main characteristics, directing the interested reader to [7] for further details. The Shapley value is a solution concept for transferable utility games that assigns a unique payoff distribution for the grand coalition to every such type of game. It is not strongly based on strategic considerations but, instead, assigns to each player his average marginal contribution to the game. Its calculation is based on obtaining the contributions of each player in all possible “entering” order and averaging all found values. Denoting a transferable

utility game by (N, ), the Shapley value payoff assigned to player i ∈ N may be defined as i (N, ) =

 |S|!(n − |S| − 1)! S⊆N:i/ ∈S

n!

[(S ∪ {i}) − (S)].

(6)

For larger games it may be more convenient to work with the precedent formula than to use the definition based on marginal vectors. The definition of the Shapley value as assigning to each player the average marginal contribution can be viewed as a justification of this solution concept by itself. In the literature there exist a number of axiomatic characterizations for this problem, and this kind of approach is reasonably frequent in cooperative game theory. For example, an axiomatic characterization of the Shapley value can be found in [9]. The aim is to study in detail the elasticity influence on the formation of coalitions, at the same time extending the study in [1] and offering abundant numerical data about this important phenomenon. To that end, many results will be shown to represent final agreements on splitting profits if acting as a cartel. In addition, it is possible to infer that price–demand elasticity is the key to increase the chance of cooperation among oligopolistic firms – if the demand becomes more elastic, higher prices will be paid for a given demand level, and the amount to share will be more significant. This will make a final agreement more likely than when the total profit is smaller, and all players tend to end up in a grand coalition, increasing their profits. In what follows, all topics needed to the full understanding of the proposal will be presented and, at the end, numerical results will be interpreted and discussed. 3. Generalized Nash equilibrium problems In this section it is introduced an approach for solving generalized Nash equilibrium problems (GNEP’s, for short) based on the Fuzzy Adaptive Simulated Annealing algorithm. Taking into account that GNEP’s have been used actively in many fields during the few past decades, and research on this topic is gaining

4

H.A. e Oliveira Jr. / Applied Soft Computing 35 (2015) 1–12

momentum, it is certainly beneficial that innovative methods may be found in order to cope with the complexity of the theme. Here, the method starts with a transformation of the original problem into a constrained optimization one and, by means of a global optimization algorithm, searches for points that represent the desired generalized Nash equilibria for the original task. Hence, the effectiveness of the overall method depends primarily on the efficacy of the chosen optimization algorithm. As GNEP’s are used to model many problems in different disciplines, and researchers in different fields tend to work independently, even today many existing results are not fully known, although some basic concepts were well-established, such as the problem definition and certain approaches considered adequate for solving the problem, as using Nikaidô–Isoda functions, for example. The general issue itself (GNEP) has been referred to as a social equilibrium problem, coupled constraint equilibrium problem, or abstract economy. In this article, the term generalized Nash equilibrium problem will be adopted and defined in the sequel [10]. As is well-known, in GNEP’s, cost functions and strategies of players are real valued and continuous, but more importantly they have the desirable feature that players face constraints dependent on strategies of their opponents. This aspect contrasts with standard Nash equilibrium problems, in which individual actions are not limited by the strategies chosen by others. As expected, the coupled constrained action space in GNEP’s makes them more difficult to handle than conventional Nash games and, naturally, algorithms able to solve them can be applied to standard NEP’s as well. GNEP’s can be situated in an environment consisting of N players, where player i controls ni variables, concentrated into a vector xi ∈ Rni . Denoting by x the vector formed by all these strategic variables, we have

⎛

x1

⎜ x2 ⎜ ⎜ x= ⎜ . ⎜ .. ⎝ 

xN

⎞ ⎟ ⎟ ⎟ ⎟, ⎟ ⎠

(7)

in which x ∈ Rn and n = n1 + n2 + · · · + nN . Also, the symbol x−i represents the vector x without the xi sub-vector,

⎛

x1

⎞

⎜ x2 ⎟ ⎜ ⎟ ⎜ ⎟ ⎜. ⎟ ⎜ .. ⎟ ⎜ ⎟ ⎜ i−1 ⎟  x ⎟ x−i = ⎜ ⎜ ⎟ ⎜ xi+1 ⎟ ⎜ ⎟ ⎜ ⎟ ⎜. ⎟ ⎜ .. ⎟ ⎝ ⎠

Assuming that agents act rationally, that is to say, given the decision vector x−i of other players, they try to choose a decision vector xi that minimizes their cost function, it is natural to expect that a search for the solution of the following optimization problem occurs at each game instance min i (xi , x−i ) s.t. xi ∈ Xi (x−i ),

(9)

xi



in which Xi (x−i )={xi ∈ Rni |(xi , x−i ) ∈ X}, and X is assumed to be a non-empty, closed and convex subset of Rn . X is the overall strategy set. T Therefore, we can define a vector x∗ = (x∗,1 , x∗,2 , . . ., x∗,N ) as being a generalized Nash equilibrium if, for all i = 1, 2, . . ., N, the subvector x*,i solves the following optimization problem min i (xi , x∗,−i ) s.t. xi ∈ Xi (x∗,−i )

(10)

xi

or, equivalently, i (x∗,i , x∗,−i ) ≤ i (xi , x∗,−i ) ∀xi ∈ Xi (x∗,−i ),

(11)

for i = 1, . . ., N. GNEP’s are also known as social equilibrium problems, or Nash equilibrium problems with shared constraints, among other expressions. The GNEP was formally introduced by Debreu [11] in 1952, and since then has been deeply associated to the mathematical economics field. However, only in recent decades interesting applications to other areas started to arise, and researchers recognized its importance outside economics. Although it can be applied to virtually any field, typical applications can be found in telecommunications and electrical engineering, computer science, market regulation, or ecological scenarios, for instance. It is worth mentioning that GNEP’s frequently have multiple solutions. In general, different solutions of specific equilibrium problems (such as the standard Nash equilibrium problems or GNEP’s) may have different implications, differently from an optimization problem in which all optimal solutions can be regarded as having equal value in terms of the objective function. From this viewpoint it may be meaningful to choose a particular generalized equilibrium, featuring a certain additional property. In [5] it is proposed the solution concept named “normalized equilibrium”, that is a special GNE characterized by conditions imposed on Lagrange multipliers associated with the constraints in problems of each player. The uniqueness of a normalized equilibrium can be established under certain conditions [5,12]. 4. Previous approaches and current trends

(8)

xn

Moreover, it is usual to write x = (xi , x−i ). Please, observe that the notation (xi , x−i ) does not mean that the components of x have been reordered, being xi the first block. An additional assumption is that all players are assumed to announce simultaneously their decisions (strategic values), and each one is assigned one loss or cost function i (x) (i : Rn → R). Hence, when overall decision vector is x, player i has to “pay” or “lose” i (x). In this fashion, the cost function of a specific player may depend on decision vectors of other players, the same occurring to the strategy sets, what is not true in standard Nash equilibrium problems.

It is a well-known fact that the standard NEP can be transformed into a variational inequality problem (VIP) [13]. In addition, it is possible to face GNEP’s as quasi-variational inequalities (QVI’s, for short) [14,15]. On the other hand, considering that efficient methods for solving QVI’s are rare, this type of characterization might not be very useful, in practical terms. In another direction, it was stated in [16,17], for example, that certain solutions of the GNEP (the normalized Nash equilibria) can be found by solving a suitable standard VIP associated to the GNEP. However, certain VIP based methods require a higher degree of smoothness of the cost functions than other approaches, based on Nikaidô–Isoda functions [18], that play a central role in the related literature. Relaxation methods using Nikaidô–Isoda functions are studied in [3,4], and a proximal-like method based on them is presented in [17]. In [19], it is reported the use of a regularized version of the Nikaidô–Isoda function in order to get different optimization problems whose global minima are the normalized solutions of the GNEP. In that work, the function and its regularized version

H.A. e Oliveira Jr. / Applied Soft Computing 35 (2015) 1–12

are formally defined, and they obtain a constrained optimization problem equivalent to the GNEP. However, the objective function corresponding to this optimization problem is nonsmooth. Then a smooth optimization problem is derived whose solutions characterize the class of normalized Nash equilibria of the GNEP. Following this, it is shown how the previous techniques can be used in order to get a smooth unconstrained optimization reformulation of the normalized GNEP solutions. The regularized Nikaidô–Isoda function was employed earlier for solving standard NEPs in [20]. The standard NEP is known to be a special case of an equilibrium programming problem [21]. This observation is emphasized in [19] by interpreting the GNEP as a particular instance of such a type of problem. Although GNEP’s represent an important modeling paradigm, the existing approaches aiming at their solution are somewhat dispersed, in terms of the diversity of theoretical tools used when trying to solve them. Rosen studied jointly convex GNEP’s [5] and this context was very frequent in the literature until very recently. He introduced the first algorithm (a projected gradienttype method) for solving jointly convex GNEP’s, and some other algorithms were later developed for this kind of problems, one of the most important is the relaxation method [3,4], based on Nikaidô–Isoda functions. Of course, other approaches may be employed when solving jointly convex problems. Among them we may highlight the recently proposed methods from [12,19,22], which are still based on the Nikaidô–Isoda function and consequently computationally intensive, or one of the variational inequality approaches [16], which are also limited to the jointly convex case. There are many other proposals, and the interested reader can consult [10] and the pointers therein to the related literature. When considering unrestricted GNEP’s the situation becomes even more difficult, because it is known that a GNEP can be reduced to a Quasi Variational Inequality [14,10,15]. On the other hand, the development of globally convergent algorithms for the solution of QVI’s is still a tough task, and this transformation is of little algorithmic use, even if one can use gap functions to transform a given GNEP into an optimization problem. Nikaidô–Isoda functions can be very useful in transforming GNEP’s into optimization problems, but the computational overhead is substantial and the conditions for establishing convergence are somewhat complicated. Unfortunately, we can say that, excepting penalty methods and algorithms for very specific applications, the study of methods for solving general GNEP’s and their convergence in more general conditions is still in its infancy. Penalty algorithms directed to the solution of GNEP’s are based on the common penalization approach, by eliminating the hard inter-related constraints in a GNEP, and reducing it to a somewhat simpler NEP. To this latter problem, we can then apply several methods based on the optimization or variational inequality techniques, for instance. The penalization approach to GNEP’s is relatively recent and was initiated by Fukushima and Pang in [23], in which they propose a sequential penalty approach for solving GNEP’s, and an infinite sequence of differentiable penalized problems is solved. The use of exact penalty approaches in which a single nondifferentiable NEP has to be solved was treated in [24], being this issue also addressed in [25], where among other things, some conditions are stated under which a penalty approach can be employed to find solutions for GNEP’s. Although many problems have to be fully addressed (among them a complete understanding of theoretical conditions under which useful results can be established for a penalty approach), penalty based methods seem to be very promising approaches for solving GNEP’s [26]. In [27] a metaheuristic paradigm is employed to study the Nash normalized equilibrium, which can be obtained by transforming the GNEP into a bi-level program with an optimal value of zero in the upper level. In that paper, it is proposed a Differential Evolution

5

based bi-level programming algorithm that uses Stochastic Ranking to handle constraints and solve the resulting bi-level programming formulation – in this work Nikaidô–Isoda functions are used as well. In [19] it is described a well-known approach, cited previously in this text and detailed here in order to illustrate the computational complexity normally associated to existing methods aimed at solving GNEP’s. In that paper, the main theoretical tool used in obtaining one of the optimization reformulations of the GNEP is the Nikaidô–Isoda function, defined by 

(x, y)=

N 

[i (xi , x−i ) − i (yi , x−i )].

(12)

i=1

From which we can define 

V (x)= sup (x, y),

x ∈ X,

(13)

y∈(x)

in which it is assumed that the supremum always exists for some y ∈ (x). Then it is possible to conclude that V(x) is nonnegative for all x ∈ (x), and that x* is a solution of the GNEP if and only if x* ∈ (x* ) and V(x* ) = 0. Therefore, finding a solution of the GNEP is equivalent to computing a global minimum of the optimization problem min V (x)

(14)

x∈(x)

with 

(x)=X1 (x−1 ) × · · · × XN (x−N ).

(15)

It is worth noting that this optimization problem usually has a nontrivial feasible set, taking into account that (x) depends on the overall vector x. But, in view of a certain result established in [19] (Lemma 2.1), problem (14) happens to be equivalent to the following one min V (x). x∈X

(16)

Note that to each evaluation of V(x) corresponds a full optimization operation, defined by (13). Even though Nikaidô–Isoda functions are very well-established in the literature, they feature certain not so favorable properties. First, given a vector x, the supremum in (13) may not exist unless additional assumptions (like the compactness of X) hold, and, besides, this supremum is usually not attained at an isolated, single point, implying that the mapping V and, therefore, the corresponding optimization reformulation (16), may be nondifferentiable. To overcome these difficulties, a simple regularization of the Nikaidô–Isoda function may be used (please, see [19]) – this idea was employed earlier in other contexts [28,20,29]. In summary, we conclude that the present approaches to the GNEP are computationally expensive, and a new approach that could simplify the search process and reduce the overall processing effort would be welcome. 5. Approach used in the solution of GNEP’s 5.1. Establishing the foundations Our aim is show that it is possible to find Nash equilibria (of games with continuous strategy sets) in GNEP’s after a single optimization run. The general idea can be applied, in principle, by means of any effective, general purpose metaheuristic global optimization method. In this work, the one known as Fuzzy Adaptive Simulated Annealing (Fuzzy ASA) was chosen. The focus is on demonstrating that the Fuzzy ASA method is able to solve GNEP’s, faced as constrained optimization problems, and without using multi-level optimization techniques. The first step is to alter the

6

H.A. e Oliveira Jr. / Applied Soft Computing 35 (2015) 1–12

problem of computing generalized Nash equilibria, that can be faced as a set of coupled and simultaneous constrained optimization problems, into a single one, whose solutions are exactly the optimizers of the original problem. This can be done by means of penalty techniques, and this transformation is shown to be effective, making it possible to succeed in our task. Besides, taking into account the ample scope of Fuzzy ASA, the proposed method tends to be more flexible than some of the already existing ones in many circumstances, taking into account that it is able to deal with nondifferentiable and nonconvex functions defined on nonconvex domains. In order to be more precise, let us consider the following GNEP. T Find x∗ = (x∗,1 , x∗,2 , . . ., x∗,N ) ∈ Rn such that, for all i = 1, 2, . . ., N, the subvector x*,i solves the following optimization problem min i (xi , x∗,−i ) s.t. xi ∈ Xi (x∗,−i ),

(17)

xi

or, equivalently, i (x∗,i , x∗,−i ) ≤ i (xi , x∗,−i ) ∀xi ∈ Xi (x∗,−i ),

(18)

for i = 1, . . ., N. 

in which Xi (x∗,−i )={xi ∈ Rni |(xi , x∗,−i ) ∈ X}, and X, the overall strategy set, is assumed to be a nonempty and closed subset of N Rn (n = n ). i=1 i It is assumed that the functions i are continuous and not necessarily differentiable or convex, being, obviously, bounded on compact domains – in the literature, further conditions are typically imposed on them.

Fig. 2. Fundamental algorithm.

The resulting cost function will have the following general appearance 

5.2. The fundamental algorithm

i=1

• Initialization For each i and Xi (x−i ), determine their specific attributes that could be translated into numerical constraints which, once incorporated in penalized cost functions, could guide the algorithm toward the set of proper solutions. For instance, if the loss functions are differentiable, forcing their partial gradients (player i has control of only ni variables) to zero helps in finding proper equilibria, considering that it is a necessary condition that occurs at extrema of that kind of mappings. The same is valid relatively to Xi (x−i ) sets, becoming easier whenever they are expressed as parametric inequalities, namely, 

Xi (x−i )={xi ∈ Rni | g i (xi , x−i ) ≤ 0}.

N 

C(x)=

(19)

• Step 1 For each attribute found in the previous step, define an additive (penalty) term capable of translating the respective attribute into a numerical value, so that to lower values correspond more favorable regions of the search space and vice-versa. For example, norms of gradients could be proper additive terms when trying to reach local extrema of differentiable functions. • Step 2 Construct the overall objective function by defining weights that will be used in the penalized global optimization step. These weights will be multipliers for each additive term, determining the relative importance of each one during the optimization process, like in typical penalized minimization situations.

N 

+

Fi × i (x) +

N 

FGradi × ||∇ ni i (x)||2 +

i=1

FConstr i × Constr i (x) +

i=1

NAttr 

(20) FAttr i × Attr i (x)

i=1

, in which Fi , FGradi , FConstr i and FAttr i are the mentioned weights,

∇ ni i (x) is the partial gradient of i (x), if it is differentiable, Constri (x) corresponds to the penalty term associated to i (x), Attr i (x) corresponds to i-th auxiliary numerical function capable of driving the optimization process into favorable regions (NAttr instances). • Step 3 Submit C(x) to the chosen global optimization algorithm (in this paper, Fuzzy ASA will be used). • Final decision step Extract the several subvectors from the overall x found in the previous step and assess the outputs. In case of unsatisfactory results, go to step 2, else emit the subvectors as the definitive solution. Fig. 2 depicts a diagram illustrating the proposed algorithm. 5.3. Implementation issues As stated above, any effective global optimization algorithm could, in principle, be used to optimize the global cost function. Our particular choice was motivated by the fact that Fuzzy ASA has shown to be very efficient in many different fields and able to deal with difficult optimization tasks.

H.A. e Oliveira Jr. / Applied Soft Computing 35 (2015) 1–12

7

6. Fuzzy Adaptive Simulated Annealing Fuzzy Adaptive Simulated Annealing is based on the ASA method [30] that uses the simulated annealing concept to optimize numerical functions, but presenting a significant number of improvements. In this section some relevant features of this important implementation will be described. 6.1. Re-annealing It is the dynamical re-scaling of parametric temperatures, adapting generating probability distribution functions to each dimension according to different sensitivities. Therefore, if the cost function does not present significant variations whenever a given parameter is varied in run-time, it is understood that, by extending the search interval amplitude in that particular dimension, the overall performance could be improved, and vice-versa. 6.2. Quenching ASA software architecture offers the possibility of adjusting several structural parameters related to the quenching process, allowing the user or any automatic control mechanism to change the default behavior of the “cooling” process, responsible for driving the evolution of the parametric temperatures. This mechanism is useful in case of stagnation inside suboptimal regions. 6.3. High parametrization degree The ASA implementation allows users to alter virtually any subsystem without expressive programming effort. It is thus possible to change the behavior of generation/acceptance processes, termination criteria, or seed generation, by changing the values of certain easily accessible variables. Consequently, after acquiring some experience on its structure, application development is usually fast. The ASA approach aims to find function minimizers inside a preestablished hyper-rectangle and generates points componentwise, according to the following scheme: xi+1 = xi + xi , with xi = yi (Bi − Ai ), [Ai , Bi ] =interval corresponding to i-th dimension, yi ∈ [−1, 1] is given by yi = sign(ui − 1/2)Ti [(1 + 1/Ti )|2ui −1| − 1], in which ui ∈ [0, 1] is generated by means of the uniform distribution, and Ti is the present temperature relative to dimension i. The compactness of search space is not a serious limitation in practical applications and, in the absence of previous information about optima locations, it suffices to choose sufficiently large hyper-rectangular domains. In addition, the quenching mechanism can benefit, in many cases, the efficiency of the convergence process, although there is always the possibility of reaching prematurely local extrema. In certain settings, however, we may not have alternative ways out of a stagnation situation. In order to overcome this problem, a fuzzy controller was constructed (Fuzzy ASA) [31]. The approach is very simple: the original ASA system is faced as a MISO (Multiple Input Single Output) dynamical system and the additional code simply inserts a feedback loop, by sampling its output (current value of objective function) and acting in its inputs (a subset of run-time adjustable parameters, related to the quenching process) according to a fuzzy control strategy (control algorithm), behaving like human beings whenever subject to similar situations. Hence, active run-time fuzzy control coupled to previously existing ASA mechanisms can regulate temperature evolution, besides being able to take evasive actions in case of premature convergence. In its present version, Fuzzy ASA code rises

Fig. 3. Fuzzy ASA diagram.

the quenching degree after detecting decreasing optimization performance or potential stagnation states, aiming to recover from a possible undesirable convergence to nonglobal optima. It is important to note that this additional module does not try to substitute the many effective devices already present in the original ASA code, just complementing them in atypical situations. Fig. 3 shows a simplified diagram of Fuzzy ASA structure. It is worth to highlight that ASA and its variations, like Fuzzy ASA for example, have been applied to several important fields, such as Signal Processing, Geophysics, or Statistics [30,32–34,31,35]. This can be understood by noting that, in practically all scientific domains, it is possible to transform certain design tasks into parametric optimization ones by finding cost functions capable of conveying all necessary constraints, and whose optimizers correspond to proper solutions for the original problems. 7. Experimental results and their interpretation In this section the results obtained after simulating the adopted model and solving the corresponding GNEP’s are presented. Simulation instances use different elasticity values in the set {−0.1, −0.2, −0.3, −0.4, −0.5, −1.0, −1.2, −1.3, −1.4 }, and some also employ values in { −1.5, −1.6, −1.7, −1.8 }. This was assumed to be sufficient to sweep the many qualitative aspects of the behavior of this model, in terms of rising the interest degree in coalition formation, mainly due to improved profits. As stated in Section 2, the problem under study features three companies (denoted by 1, 2, and 3) that can, at least theoretically, form five different coalition configurations, namely, { {1}, {2}, {3} }, { {1, 2}, {3} }, { {1, 3}, {2} }, { {2, 3}, {1} }, and { {1, 2, 3} }. Apart from the last one, in each case there exists a competitive setting that can be approached by means of the Nash–Cournot equilibrium concept. Besides, after finding the solutions corresponding to each coalition member, there is the problem of dividing the profits among the various companies participating in non-unitary groups. In what follows, results about production levels, profits, and unitary prices relative to each possible coalition structure are presented and illustrated with graphs. Production level measurement unit is MW, profits are in US dollars/h and prices in US dollars/MWh. 7.1. Structure { {1}, {2}, {3} } In this configuration it is possible to see from Figs. 4–6 that with elasticity values lower than −0.5, individual production and profit levels start to behave differently. Basically, there are three distinct regions in the graphs: - On the right and after −0.5, production levels and corresponding profits evolve similarly for the 3 players, and unitary prices rise very slowly in an almost linear way. In this region all players are fully able to cope with additional energy demand. - In the center ([−1.2, −0.5]), players 2 and 3 bifurcate their production levels, and player 1 becomes unable to cope with further

8

H.A. e Oliveira Jr. / Applied Soft Computing 35 (2015) 1–12

140

Prod 1 Prod 2 Prod 3

120

Production

100 Production

Prod {1 + 2 } Prod { 3 }

200

80 60 40

150 100 50

20 0

0 − 1.4

− 1.2

−1

− 0.8 − 0.6 Elasticity

− 0.4

− 0.2

− 1.4

0

−1

− 0.8 − 0.6 Elasticity

− 0.4

− 0.2

0

Fig. 7. Production – coalition structure { {1, 2}, {3} }.

Fig. 4. Production – coalition structure { {1}, {2}, {3} }.

·104

·104

3 Profit 1 Profit 2 Profit 3

Profit { 1 + 2 } Profit { 3 }

2.5 2 Profit

1.5

Profit

− 1.2

1

1.5 1

0.5

0.5 0

0

− 1.4 − 1.4

− 1.2

−1

− 0.8 − 0.6 Elasticity

− 0.4

− 0.2

− 1.2

−1

0

− 0.8 − 0.6 Elasticity

− 0.4

− 0.2

0

Fig. 8. Profit – coalition structure { {1, 2}, {3} }. Fig. 5. Profit – coalition structure { {1}, {2}, {3} }.

7.2. Structure { {1, 2}, {3} } market demand, taking into account that its production limit has been reached. Unitary prices start to “accelerate”, mainly for elasticity values lower than −1.2. - On the left (below −1.2), players 2 and 3 reach their production limits and unitary prices rise at very high rates, reflecting the exhaustion in production ability. As an obvious consequence, all firms experience rising profits.

In summary, the simulations confirm the fact that when a given good gets scarce, its market value tends to rise at expressive rates.

In this configuration it is possible to see that with elasticity values lower than −1.0, individual production and profit levels start to behave differently. Please, refer to Figs. 7–9 to follow the interpretation of numerical results. We can find two distinct regions in the graphs: - On the right and after −1.0, production levels and corresponding profits evolve similarly for both players, and unitary prices rise very slowly in an almost linear way. In this region all players are able to provide additional energy. - On the left (below −1.0), player {3} reaches its production limit and unitary prices rise at very high rates, reflecting the partial

140

140

138 136

120

Price

Price

130

134 132

110

130 100

128 − 1.4

− 1.2

−1

− 0.8 − 0.6 Elasticity

− 0.4

Fig. 6. Prices – coalition structure { {1}, {2}, {3} }.

− 0.2

0

− 1.4

− 1.2

−1

− 0.8 − 0.6 Elasticity

− 0.4

Fig. 9. Prices – coalition structure { {1, 2}, {3} }.

− 0.2

0

H.A. e Oliveira Jr. / Applied Soft Computing 35 (2015) 1–12

Prod { 1 + 3 } Prod { 2 }

200

9

155 150 Price

Production

150 100

145 140 135

50

130 0 − 1.6

− 1.4

− 1.2

−1

− 0.8 − 0.6 Elasticity

− 0.4

− 0.2

− 1.6

0

Fig. 10. Production – coalition structure { {1, 3}, {2} }.

− 1.4

− 1.2

−1

− 0.8 − 0.6 Elasticity

− 0.4

− 0.2

0

Fig. 12. Prices – coalition structure { {1, 3}, {2} }.

exhaustion in production ability. As an obvious consequence, all players get rising profits until reaching the overall production limit.

200 Production

Again, it is possible to conclude that the scarcity of goods is directly linked to rising unitary prices. As a kind of by-product, it is reasonable to digress and say that this type of numerical simulation may shed some light in the inner mechanisms of certain inflationary settings.

Prod { 2 + 3 } Prod { 1 }

250

150 100 50 0 − 1.6 − 1.4 − 1.2

7.3. Structure { {1, 3}, {2} } In this case the qualitative behavior is identical to the previous one, and the numerical results corresponding to this case are shown in Figs. 10–12. Following the same line of reasoning, we can see that with elasticity values lower than −1.0, individual production and profit levels change their behavior, bifurcating and assuming distinct growth rates. There are two distinct regions in the graphs: - On the right and after −1.0, production levels and profits progress similarly for both agents, and unitary prices rise slowly. This region features full production capacity and all players are able to provide more energy. - On the left, player {2} reaches its production limit and unitary prices rise at higher rates, provoked by the partial exhaustion in production ability. Consequently, all players get rising profits until they attain the overall production limit.

− 1 − 0.8 − 0.6 − 0.4 − 0.2 Elasticity

0

Fig. 13. Production – coalition structure { {2, 3}, {1} }.

7.4. Structure { {2, 3}, {1} } According to Figs. 13–15 it is possible to conclude that with elasticity values lower than −0.5, individual production and profit levels start to behave differently. It is also possible to find two distinct regions in the graphs: - On the right and after −0.5, production levels and corresponding profits evolve similarly for the 2 players, and unitary prices rise very slowly. In this region all players are able to provide additional energy. - On the left (below −0.5), player {1} reaches its production limit and unitary prices rise at very high rates, reflecting the partial exhaustion in production ability. Therefore, all players get rising profits until reaching the corresponding overall production limits.

·104

·104 Profit { 1 + 3 } Profit { 2 }

3

Profit { 2 + 3 } Profit { 1 }

4 3

2

Profit

Profit

2.5

1.5 1

2 1

0.5 0 − 1.6

0 − 1.4

− 1.2

−1

− 0.8 − 0.6 Elasticity

− 0.4

Fig. 11. Profit – coalition structure { {1, 3}, {2} }.

− 0.2

0

− 1.6 − 1.4 − 1.2

− 1 − 0.8 − 0.6 − 0.4 − 0.2 Elasticity

Fig. 14. Profit – coalition structure { {2, 3}, {1} }.

0

10

H.A. e Oliveira Jr. / Applied Soft Computing 35 (2015) 1–12

170

191.6 191.4

160 Price

Price

191.2 150

191

140

190.8

130

190.6 − 1.6 − 1.4 − 1.2

− 1 − 0.8 − 0.6 − 0.4 − 0.2 Elasticity

0

− 1.8 − 1.6 − 1.4 − 1.2

Fig. 15. Prices – coalition structure { {2, 3}, {1} }.

− 1 − 0.8 − 0.6 − 0.4 − 0.2 Elasticity

0

Fig. 18. Prices – grand coalition structure { { 1, 2, 3 } }.

·104 300

Isolated participation Shapley value

1.2

Production

250 1

200 150

0.8 Value

100 50 0 − 1.8 − 1.6 − 1.4 − 1.2

− 1 − 0.8 − 0.6 − 0.4 − 0.2 Elasticity

0.6 0.4

0

0.2

Fig. 16. Production – grand coalition structure { {1, 2, 3} }.

0 7.5. Grand coalition {{1, 2, 3 }} In this configuration all variables behave in a practically linear way. This can be associated to the absence of competition, considering that now there is only one agent providing energy to the market. As the demand × price curves are straight lines and there is only one player, it is expected that within the region in which the group of firms is able to satisfy the total demand all three curves get this form. Numerical results are shown in Figs. 16–18. 7.6. Coalition feasibility results

− 1.4

− 1.2

−1

− 0.8 − 0.6 Elasticity

− 0.4

− 0.2

0

Fig. 19. Isolated profit and Shapley value when in grand coalition – Firm 1.

of elasticity values. It is worth to highlight that the comparisons are made based upon very conservative assumptions. That is, on one side we have the profits obtained in a configuration of isolated participation, which is usually a relatively unfavorable circumstance, and these figures are compared to the Shapley values calculated from the grand coalition structure – in this case, the

·104

By analyzing Figs. 19–21, it is possible to draw conclusions about coalition formation benefits for each company as a function

2

Isolated participation Shapley value

·104 1.5

6 Value

5 Profit

4

1

3 0.5

2 1

0

0 − 1.8 − 1.6 − 1.4 − 1.2

− 1 − 0.8 − 0.6 − 0.4 − 0.2 Elasticity

Fig. 17. Profit – grand coalition structure { { 1, 2, 3 } }.

0

− 1.4

− 1.2

−1

− 0.8 − 0.6 Elasticity

− 0.4

− 0.2

Fig. 20. Isolated profit and Shapley value when in grand coalition – Firm 2.

0

H.A. e Oliveira Jr. / Applied Soft Computing 35 (2015) 1–12

2

market – such a tool can be considered an extension of previous initiatives. Themes for future work include the exploration of the proposed method for economic regulation of specific markets, mainly in the detection of already existing oligopolies or indication of their early formation. Yet another relatively easy task may be to adapt the present study for application to other fields, like oil and gas exploration, for example.

·104 Isolated participation Shapley value

1.8 1.6 1.4 1.2 Value

11

1 References

0.8 0.6 0.4 0.2 0 − 1.4

− 1.2

−1

− 0.8 − 0.6 Elasticity

− 0.4

− 0.2

0

Fig. 21. Isolated profit and Shapley value when in grand coalition – Firm 3.

most comfortable situation, taking into account the high level of unitary prices along the spectrum of elasticity values. By observing the graphs it is possible to conclude that, inside the interval [−1.3, −0.3], all 3 firms certainly get more profitable results when acting in an associative form (grand coalition). However, outside that interval, there is a visible trend toward an inversion of that tendency, that is to say, for elasticity values inferior to −1.4 and greater than −0.1, the two settings seem to drive profits to a coincidence state, signaling a potential change in profitability status. 8. Conclusion This article presented a study about the influence of pricedemand elasticity variation in electricity markets, in terms of coalition formation among firms. The method used here is based on facing each equilibrium situation (in a given energy auction, for example) as a GNEP, transforming it into a constrained global optimization problem, and finding at least one global minimizer. Such results correspond to generalized Nash equilibria – solutions that establish (typically) stable configurations in which players cannot improve their positions in a unilateral way. In other words, the proposal makes it possible to estimate how interesting is, for competing companies, to join and cooperate with other players in order to get greater profits than those possible when participating in an isolated way. Although the proposed method is very abrangent, its efficacy was exemplified by means of a very significant problem, previously investigated in [1]. It was shown that it is possible, for a general purpose metaheuristic optimization method, to find solutions for GNEP’s in a single optimization step. So, one major aim was to demonstrate the effectiveness of the Fuzzy ASA method with respect to its ability to find solutions for GNEP’s, without using multi-level optimization techniques. To that end, the well-known optimization formulation of the problem of computing generalized Nash equilibria [12,19,22] was altered, in the sense that the original problem, faced as a set of coupled and simultaneous constrained optimization tasks, is transformed into a single one, by means of a penalty method approach – such a transformation made it possible to succeed in the intended task. In addition, the proposed method can be more abrangent than some of the already existing ones in many circumstances, taking into account that Fuzzy ASA is able to deal with nondifferentiable functions defined on nonconvex domains. In summary, the main contribution of the present work is to furnish an effective tool to evaluate and synthesize economic strategies to cope with significant problems in the electric energy

[1] J. Contreras, M. Klusch, J.B. Krawczyk, Numerical solutions to Nash–Cournot equilibria in coupled constraint electricity markets, IEEE Trans. Power Syst. 19 (2004) 195–206. [2] B.F. Hobbs, Linear complementarity models of Nash–Cournot competition in bilateral and POOLCO power markets, IEEE Trans. Power Syst. 16 (2001) 194–202. [3] J.B. Krawczyk, S. Uryasev, Relaxation algorithms to find Nash equilibria with economic applications, Environ. Model. Assess. 5 (2000) 63–73. [4] S. Uryasev, R.Y. Rubinstein, On relaxation algorithms in computation of noncooperative equilibria, IEEE Trans. Autom. Control 39 (1994) 1263–1267. [5] J.B. Rosen, Existence and uniqueness of equilibrium points for concave Nperson games, Econometrica 33 (1965) 520–534. [6] W. Ongsakul, D.N. Vo, Artificial Intelligence in Power System Optimization, CRC Press, Boca Raton, FL, 2013. [7] H. Peters, Game Theory – A Multi-Leveled Approach, Springer-Verlag, Berlin Heidelberg, 2008. [8] G. Chalkiadakis, E. Elkind, M. Wooldridge, Computational Aspects of Cooperative Game Theory, Morgan & Claypool, San Rafael, CA, 2012. [9] L.S. Shapley, A value for n-person games, in: A.W. Tucker, H.W. Kuhn (Eds.), Contributions to the Theory of Games II, Princeton University Press, Princeton, 1953, pp. 307–317. [10] F. Facchinei, C. Kanzow, Generalized Nash Equilibrium Problems, Institute of Mathematics, University of Würzburg, Würzburg, Germany, 2009, Preprint 290. [11] G. Debreu, A social equilibrium existence theorem, Proc. Natl. Acad. Sci. 38 (1952) 886–893. [12] A. von Heusinger, C. Kanzow, M. Fukushima, Newton’s Method for Computing a Normalized Equilibrium in the Generalized Nash Game Through Fixed Point Formulation, Institute of Mathematics, University of Würzburg, Würzburg, Germany, 2009, Preprint 286. [13] F. Facchinei, J.-S. Pang, Finite-Dimensional Variational Inequalities and Complementarity Problems, vol. I, Springer-Verlag, New York, 2003. [14] A. Bensoussan, Points de Nash dans le cas de fonctionelles quadratiques et jeux differentiels lineaires a N personnes, SIAM J. Control 12 (1974) 460– 499. [15] P.T. Harker, Generalized Nash games and quasivariational inequalities, Eur. J. Oper. Res. 54 (1991) 81–94. [16] F. Facchinei, A. Fischer, V. Piccialli, On generalized Nash games and VIs, Technical Report, Department of Computer and System Sciences “A. Ruberti”, Università di Roma “La Sapienza”, Rome, Italy, 2006. ´ [17] S.D. Flåm, A. Ruszczynski, Noncooperative Convex Games, Computing Equilibria by Partial Regularization. IIASA Working Paper 94-42, Laxenburg, Austria, 1994. [18] H. Nikaidô, K. Isoda, Note on noncooperative convex games, Pac. J. Math. 5 (1955) 807–815. [19] A. von Heusinger, C. Kanzow, Optimization Reformulations of the Generalized Nash Equilibrium Problem Using Nikaido–Isoda Type Functions, Institute of Mathematics, University of Würzburg, Würzburg, Germany, 2006, Preprint 269. [20] G. Gürkan, J.-S. Pang, Approximations of Nash Equilibria, Technical Report, Department of Mathematical Sciences, Rensselaer Polytechnic Institute, Troy, NY, 2006. [21] S.D. Flåm, A.S. Antipin, Equilibrium programming using proximal-like algorithms, Math. Program. 78 (1997) 29–41. [22] A. von Heusinger, C. Kanzow, SC1 optimization reformulations of the generalized Nash equilibrium problem, Optim. Methods Softw. 23 (2008) 953– 973. [23] M. Fukushima, J.-S. Pang, Quasi-variational inequalities, generalized Nash equilibria and multi-leader-follower games, Comput. Manag. Sci. 2 (2005) 21–56. [24] F. Facchinei, J.-S. Pang, Exact penalty functions for generalized Nash problems, in: G. Di Pillo, M. Roma (Eds.), Large Scale Nonlinear Optimization, Springer US, 2006, pp. 115–126. [25] M. Fukushima, Restricted Generalized Nash Equilibria and Controlled Penalty Algorithm, Technical Report 2008-007, Department of Applied Mathematics and Physics, Kyoto University, 2008. [26] F. Facchinei, C. Kanzow, Penalty Methods for the Solution of Generalized Nash Equilibrium Problems (With Complete Test Problems). Technical Report, Institute of Mathematics, University of Würzburg, Würzburg, Germany, 2009. [27] A. Koh, Differential Evolution Based Bi-Level Programming Algorithm for Computing Normalized Nash Equilibrium, in: A. Gaspar-Cunha, et al. (Eds.), Soft

12

[28]

[29] [30] [31]

H.A. e Oliveira Jr. / Applied Soft Computing 35 (2015) 1–12 Computing in Industrial Applications, AISC 96, Springer-Verlag, Berlin Heidelberg, 2011. M. Fukushima, Equivalent differentiable optimization problems and descent methods for asymmetric variational inequality problems, Math. Program. 53 (1992) 99–110. G. Mastroeni, Gap functions for equilibrium problems, J. Global Optim. 27 (2003) 411–426. L. Ingber, Adaptive simulated annealing (ASA): lessons learned, Control Cybern. 25 (1) (1996) 33–54. H.A. e Oliveira Jr., L. Ingber, A. Petraglia, M.R. Petraglia, M.A.S. Machado, Stochastic Global Optimization and Its Applications with Fuzzy Adaptive Simulated Annealing, Springer-Verlag, Berlin-Heidelberg, 2012.

[32] L. Ingber, B. Rosen, Genetic algorithms and very fast simulated reannealing: a comparison, Math. Comput. Model. 16 11 (1992) 87–100. [33] H.A. Oliveira Jr., A. Petraglia, M.R. Petraglia, Frequency domain FIR filter design using fuzzy adaptive simulated annealing, Circuits Syst. Signal Process. 28 (6) (2009) 899–911. [34] H.A. Oliveira Jr., A. Petraglia, Global optimization using space-filling curves and measure-preserving transformations, in: A. Gaspar-Cunha, et al. (Eds.), Soft Computing in Industrial Applications, AISC 96, Springer-Verlag, Berlin Heidelberg, 2011, pp. 121–130. [35] R. Pachter, Z. Wang, Adaptive simulated annealing and its application to protein folding, in: C.A. Floudas, P.M. Pardalos (Eds.), Encyclopedia of Optimization, Springer-Verlag, 2009, pp. 21–26.