ELECTRIC POWER BIDDING MODELS FOR COMPETITIVE MARKETS Except where reference in made to the work of others, the work described in this dissertation is my own or was done in collaboration wit my advisory committee. This dissertation does not include proprietary or classified information. ___________________________ Ahmet D. Yucekaya Certificate of Approval: ____________________________ ____________________________ Chan S. Park Jorge Valenzuela, Chair Professor Associate Professor Industrial and Systems Engineering Industrial and Systems Engineering ____________________________ ____________________________ Gerry V. Dozier George T. Flowers Professor Dean Computer Science Graduate School ELECTRIC POWER BIDDING MODELS FOR COMPETITIVE MARKETS Ahmet D. Yucekaya A Dissertation Submitted to the Graduate Faculty of Auburn University in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy Auburn, Alabama December 19, 2008 iii ELECTRIC POWER BIDDING MODELS FOR COMPETIVE MARKETS Ahmet D. Yucekaya Permission is granted to Auburn University to make copies of this dissertation at its discretion, upon request of individuals or institutions and at their expense. The author reserves all publication rights. ____________________________ Signature of Author ____________________________ Date of Graduation iv VITA Ahmet D. Yucekaya, son of Husne and Mehmet Yucekaya, was born on December 2, 1978, in Kahramanmaras, Turkey. He attended Istanbul Technical University, Istanbul where he earned his Bachelor of Science Degree in Industrial Engineering in 2000, and Master of Science Degree in Strategy Development in 2004. He has worked as a consultant and a Planning manager in Industry during 2000 and 2004. He enrolled in Graduate School at University of Pittsburgh in August 2004 where he earned his Master of Science Degree in Industrial Engineering in 2005. He then enrolled in Graduate School at Auburn University in August 2005. v DISSERTATION ABSTRACT ELECTRIC POWER BIDDING MODELS FOR COMPETIVE MARKETS Ahmet D. Yucekaya Doctor of Philosophy, December 19, 2008 (M.S., University of Pittsburgh, 2005) (M.S., Istanbul Technical University, Turkey, 2004) (B.S., Istanbul Technical University, Turkey, 2000) 121 Typed Pages Directed by Jorge Valenzuela Profit maximization for power companies is highly related to the bidding strategies used. In order to sell electricity at high prices and maximize profit, power companies need suitable bidding models that consider power operating constraints and price uncertainty within the market. Therefore, models that include price uncertainty are needed to analyze the market and to make better bidding decisions. vi In this dissertation, the main objective is to develop bidding models for electric power generators and wholesale power suppliers under price uncertainty. A quadratic programming model and a nonlinear programming model were developed to find a solution to the bidding problem. However, in these models the computational time increases exponentially as the size of the problem increases. To overcome this limitation, two particle swarm optimization models are developed. The first method uses a conventional particle swarm optimization technique to find bid prices and quantities under the rules of a competitive power market. The second method uses a decomposition technique in conjunction with the particle swarm optimization approach. In addition, a spreadsheet based simulation algorithm is developed to evaluate a bid offer under given price samples. It is shown that for nonlinear cost functions particle swarm optimization solutions provide higher expected profits than marginal cost based bidding. A model to find an equilibrium solution in competitive power markets for power suppliers bidding into day-ahead market under forecasted demand is also developed. vii ACKNOWLEDGEMENTS First, I would like to thank Dr. Jorge Valenzuela, my advisor and mentor, for his invaluable guidance and support throughout this dissertation. I also would like to thank Dr. Chan S. Park and Dr. Gerry Dozier for their invaluable suggestions. Thanks also to Dr. Latif Kalin for being my out-side reader. Finally, I would like to express my endless appreciation to my parents Husne and Mehmet Yucekaya and to my other family members, for their continuous support and motivation all along my life. viii Style manual or journal used Bibliography conforms to those of Institute of Electrical and Electronics Engineers (IEEE) Transactions Computer software used AMPL CPLEX 7.0 MATLAB MINITAB Microsoft Visual Studio C/C++ Microsoft Office Excel Microsoft Office Word ix TABLE OF CONTENTS LIST OF TABLES xii LIST OF FIGURES xv NOMENCLATURE xvii CHAPTER 1. INTRODUCTION 1 CHAPTER 2. REVIEW OF LITERATURE AND BACKGROUND 5 2.1 Literature Review 5 2.2 Background 8 2.2.1 Auction Theory 8 2.2.2 Market Design 10 2.2.3 Locational Marginal Pricing 12 2.2.4 Cost of Generation 13 2.2.5 PJM Power Market 14 2.2.6 Bidding into Market 15 2.2.7 Power Market Equilibrium 19 CHAPTER 3. SPREADSHEET BASED SIMULATION FOR 22 OFFER EVALUATION 3.1 Bid Simulator 22 3.2 Simulation Functionality 23 3.3 Numerical Example 28 x 3.4 Conclusion 31 CHAPTER 4. ELECTRIC POWER BIDDING UNDER PRICE UNCERTAINITY 32 4.1 Price Uncertainty 32 4.2 The bidding model 34 4.3 Quadratic Programming Model 37 4.4 Nonlinear Programming Model 37 4.5 Marginal Cost Biding 39 4.6 Numerical Example and Analysis 40 4.6.1 Quadratic Programming 40 4.6.2 Nonlinear Programming 40 4.6.3 Marginal Cost Bidding 42 4.7 Conclusion 42 CHAPTER 5. PARTICLE SWARM OPTIMIZATION FOR BIDDING INTO 44 MARKET UNDER PRICE UNCERTAINITY 5.1 Rationale for Heuristics Approach 44 5.2 Particle Swarm Optimization 44 5.3 Conventional Particle Swarm Optimization 47 5.4 Decomposition Based Particle Swarm Optimization 49 5.5 dBPSO Genetic Representation 50 5.6 Numerical Example and Analysis 52 5.6.1 Comparison Analysis of cPSO and Marginal 54 Cost based Bidding 5.6.2. Comparison Analysis of cPSO and dBPSO 56 xi 5.6.3. Converging Process of dBPSOB and dBPSOQ 59 5.6.4 Impact of the order of dBPSOB and dBPSOQ 61 5.6 Conclusion 63 CHAPTER 6. AGENT BASED PARTICLE SWARM OPTIMIZATION FOR 65 SUPPLY FUNCTION EQUILIBRIUM 6.1 Introduction 65 6.2 Supply Function Equilibrium Model 65 6.3 Agent Based Modeling and Simulation 69 6.4 Agent based particle swarm optimization model 72 6.5 Numerical Example and Analysis 78 6.5.1 ABPSO Experiment with Duopoly 78 6.5.2 ABPSO Experiment with m=5 firms 83 6.5.3 ABPSO Experiment with m=10 firms 89 6.6 Conclusion 93 CHAPTER 7. CONCLUSION AND FUTURE RESEARCH 94 REFERENCES 96 xii LIST OF TABLES Table 2.1. Comparison of offer types 17 Table 2.2. Bid components in offer types 18 Table 3.1. An example bidding solution to the problem 29 Table 3.2. Descriptive statistics for daily profits 29 Table 4.1. Price sample for the example problem 40 Table 4.2. Optimum solution to the problem 40 Table 4.3. Price samples used in the problem 41 Table 4.4. Optimum solution to the problem 41 Table 4.5. Bid prices and quantities for marginal cost bidding 42 Table 5.1. Day-ahead market price samples for PSO 53 Table 5.2. Quantity-price solutions for GEN-1 and GEN-2 55 Table 5.3. Parameter set for the experiment 57 Table 5.4. Statistical summary of the results of the cPSO solutions 57 Table 5.5. Bid prices and quantities for best solutions of cPSO 57 Table 5.6. ANOVA test results for cPSO 58 Table 5.7. Statistical summary of dBPSO solutions 58 Table 5.8. Bid prices and quantities for best solutions of dBPSO 58 Table 5.9. ANOVA test results for dBPSO 59 Table 5.10. ANOVA test results for method comparisons 59 Table 5.11. Statistical summary of dBPSO solutions 63 xiii Table 6.1. Day-ahead demand for next day 78 Table 6.2. Cost functions and capacities of the firms (m=2) 78 Table 6.3. Equilibrium prices for the first stopping rule (m=2) 79 Table 6.4. Results found in the first stopping rule (m=2) 79 Table 6.5. Equilibrium prices for the second stopping rule (m=2) 80 Table 6.6. Results found in the second stopping rule (m=2) 81 Table 6.7. Equilibrium prices for the third stopping rule (m=2) 81 Table 6.8. Results found in the third stopping rule (m=2) 82 Table 6.9.Overview of the results found in each method (m=2) 82 Table 6.10. Profit increases in each rule (m=2) 83 Table 6.11. Cost functions and market capacities of the firms for (m=5) 83 Table 6.12. Equilibrium prices for first stopping rule (m=5) 84 Table 6.13. Results found for firm 1, 2 and 3 in first stopping rule (m=5) 84 Table 6.14. Results found for firm 4 and 5 in first stopping rule (m=5) 85 Table 6.15. Equilibrium prices for second stopping rule (m=5) 85 Table 6.16. Results found for firm 1, 2 and 3 in second stopping rule (m=5) 86 Table 6.17. Results found for firm 4 and 5 in second stopping rule (m=5) 86 Table 6.18. Equilibrium prices for third stopping rule (m=5) 87 Table 6.19. Results found for firm 1, 2 and 3 in third stopping condition (m=5) 87 Table 6.20. Results found for firm 4 and 5 in third stopping condition (m=5) 88 Table 6.21. Overview of the results found in each method (m=5) 88 Table 6.22. Profit increases in each rule (m=5) 89 Table 6.23. Cost functions and market capacities of the firms (m=10) 89 xiv Table 6.24. Equilibrium prices for second stopping rule (m=10) 90 Table 6.25. Results found for firm 1, 2 and 3 (m=10) 90 Table 6.26. Results found for firm 4, 5 and 6 (m=10) 91 Table 6.27. Results found for firm 7 and 8 (m=10) 91 Table 6.28. Results found for firm 9 and 10 (m=10) 92 Table 6.29. Unit?s profits (m=10) 92 Table 6.30. Profit increases in rule 2 (m=10) 92 xv LIST OF FIGURES Figure 2.1. Services in a typical SMD 10 Figure 2.2. Electric power system 11 Figure 2.3. A generator?s offer curve in the PJM day-ahead market 16 Figure 3.1. Pseudo code of the simulation model 23 Figure 3.2. Bid Simulator main screen 24 Figure 3.3. Bid Simulator menu 24 Figure 3.4. Bid Simulator simulation and price generation 25 Figure 3.5. Bid Simulator hourly profits 26 Figure 3.6. Bid Simulator daily profits 26 Figure 3.7. Bid Simulator output analysis 28 Figure 3.8. Expected hourly price graph from generated prices 29 Figure 3.9. Expected hourly profits based on generated prices 30 Figure 3.10. Histogram of daily profits 30 Figure 3.11. Cumulative histogram of daily profits 31 Figure 4.1. Day-ahead and real time prices for February 29, 2008 33 Figure 4.2. Bid prices and power quantity intervals 35 Figure 5.1. Pseudo code for the PSO method 45 Figure 5.2. Pseudo code for the cPSO method 48 Figure 5.3. Pseudo code for the dBPSO method 51 Figure 5.4. Hourly day-ahead market price samples 53 xvi Figure 5.5a. Cost function of GEN-1 54 Figure 5.5b. Cost function of GEN-2 54 Figure 5.6. Percentage of increase of cPSO solution in relation to MC 56 Figure 5.7. Converging process of dBPSOB and dBPSOQ to the solution (C=200) 60 Figure 5.8. Converging process of dBPSOB and dBPSOQ to the solution (C=100) 60 Figure 5.9. Converging process of dBPSOB and dBPSOQ to the solution (C=50) 61 Figure 5.10. Converging process of dBPSOB and dBPSOQ to the solution (C=200) 62 Figure 5.11. Converging process of dBPSOB and dBPSOQ to the solution (C=100) 62 Figure 5.12. Converging process of dBPSOB and dBPSOQ to the solution (C=50) 63 Figure 6.1. Equilibrium process in day-ahead power market 66 Figure 6.2. Agent based modeling and simulation 69 Figure 6.3. General flow of a typical ABMS 70 Figure 6.4. Components of ABMS method 73 Figure 6.5. The flow of the ABPSO 74 xvii NOMENCLATURE Parameters i Block number of bidding offer j Firms in the market t Time of bidding period m Number of firms r Iteration number N Number of blocks in bidding offer T Time period that the bid is valid for R Number of rounds or iterations Bmax Price Cap or Maximum allowable bid price in the market Qmax Maximum generation capacity of the generator (MWh) Qjmax Maximum market capacity of the firm j (MWh) a1 No-load cost a2 Linear cost coefficient a3 Quadratic cost coefficient aj1 No-load cost for firm j aj2 Linear cost coefficient for firm j aj3 Quadratic cost coefficient for firm j Dt Demand at time t ?j Profit function of firm j xviii Pt Market price at time t k tP Market clearing price at time t of scenario k ($/MWh) r tP Market clearing price at time t of iteration r ($/MWh) C(q) Cost of generating q MWh ($/MWh) ?1 Stopping condition value for first rule ?2 Stopping condition value for second rule ?3 Stopping condition value for third rule C Cycle value for particle swarm optimization analysis R Number of replication for particle swarm optimization analysis Decision variables qjt(Pt) Dispatched bid quantity of firm j at time t at market price Pt q-jt(Pt) Cumulative Bid quantity of competitor?s of firm j at time t at market price Pt ?qi Bid energy amount increase of block i (MWh) bi Bid price of block i ($/MWh) qt Total energy generated at hour t (MWh) ?qji Bid energy amount increase of block i of firm j (MWh) bji Bid price of block i of firm j ($/MWh) ?t Stopping condition criterion value for first rule t ? ? Stopping condition criterion value for second rule tw? Stopping condition criterion value for third rule 1 CHAPTER 1 INTRODUCTION After the 90s, many countries changed the economics of their electricity markets from monopolies to oligopolies in an effort to increase competition. One of the main market competition structures used in the new deregulated environments is the poolco [1]. A poolco market is a central auction that brings regional buyers and sellers together. All competitive power generators (supply) and buyers (demand) are required to submit blocks of energy amounts and corresponding prices they are willing to receive from or pay to the pool. The prices and quantities submitted by the market participants are binding obligations as they require financial commitments to the market. Once all the supply and demand bids have been submitted and the bidding period ends, an Independent System Operator (ISO) ranks these quantity-price offers based on the least- cost for selling bids and the highest price for buying bids. The ISO then matches the selling bids with buying offers such that the highest offers are matched with the lowest selling bids. The point that supply meets demand determines the market clearing price (MCP). Perfect competition and oligopoly are two models of interest in the deregulation of the electricity market. Under a perfect competition model, power suppliers are expected to bid their marginal costs. For generating units with a nonlinear cost function (such as coal and gas fired units), the marginal cost depends on the quantity of electricity 2 produced by the unit. This implies the necessity of knowing the amount of energy that will be offered to the market before the marginal production cost can be computed. However, if the electricity market is not perfectly competitive, a power supplier may increase benefits by bidding a price higher than its marginal production costs [2]. This behavior is called strategic bidding [3]. This strategy imposes the risk of producing no profit at all if the bid price is too high. There is a risk that the supplier?s bid might be placed in jeopardy. Thus, the strategic bidding problem (SBP) is to determine proper sizes and bid prices such as to maximize expected profits. In both situations described above, the MCP plays an important role in the SBP since it determines what blocks will be selected by the market clearing mechanism. The MCP is the result of a complex interaction among producers and consumers. When power suppliers have the ability to affect the MCP by altering its power output, oligopoly models such as Cournot [4] or supply function equilibrium [5] are usually adopted. These models can incorporate detailed economic information about the system, but they are difficult to solve and they present theoretical problems related to the lack of equilibrium or the existence of multiple equilibriums. In addition, these models require suppliers? cost data, which may not be openly available. An alternative approach, suggested in [6], is to assume that the future values of the MCP are actually unknown by the market participants since the interaction of processes that governs the MCP is too complex to model. Thus, the MCP can be modeled as a random variable to represent the complexity and uncertainty involved in current electricity markets. The advantage of this modeling approach is that it allows the inclusion of the MCP as an exogenous variable to the SBP. In [6], it has been shown that when the MCP is assumed to be exogenous each 3 generating unit can be considered separately. Nevertheless, this modeling approach can be considered as an approximation to the game theoretic method. Its accuracy will depend in part on whether there is a chance that the generating unit that is bidding into the market will be the unit that will set the MCP. In this dissertation, bidding models for two power producer behaviors are developed. The first set of bidding models assumes that power producers cannot influence the price, i.e. they behave as price takers; while the second set assumes that power producers influence the prices with their bids. One contribution of this dissertation is to model the SBP and find a solution using quadratic and nonlinear programming given that the power producer has imperfect price estimations. However, an optimal solution can only be obtained within a reasonable computational time for a limited number of price scenarios. A second contribution is the development of a spreadsheet based simulation algorithm to evaluate bids and help power companies with their decision analysis. A third contribution is the application of heuristic approaches to solve the SBP. Two models are developed to demonstrate the effectiveness of particle swarm optimization method to solve the SBP. A fourth contribution is the application of agent- based simulation method for finding bids when power producers compete for fixed demand. Numerical results show that the agent based simulation model is capable of finding an equilibrium solution for each power supplier. The remainder of this dissertation is organized as follows: Chapter 2 provides a review of related research currently available about auction designs, bidding processes, pricing and market equilibrium. Chapter 3 provides details about the developed spreadsheet based simulation method and its application on bid evaluation. Chapter 4 4 describes the model formulation of the SBP and implemented solution techniques. Chapter 5 provides a description of the conventional particle swarm optimization (cPSO) and decomposition based particle swarm optimization (dBPSO) methods. Chapter 6 describes the formulation and model of the proposed market equilibrium model and developed solution technique. The conclusion and implications for future research may be found in Chapter 7. 5 CHAPTER 2 REVIEW OF LITERATURE AND BACKGROUND 2.1 Literature Review Since the 1980s, much effort has been made to restructure the traditional monopoly of the power industry with the objectives of introducing fair competition and improving economic efficiency [7]. The electricity supply industry has unique features such as a limited number of producers, large investment size (which poses a barrier to entering the market), transmission constraints (which are obstacles for consumers to effectively reach many generators), and transmission losses (which discourage consumers from purchasing power from distant suppliers) [7].These features force market players to be more aggressive on their bidding strategies; it also makes them construct models that carefully consider their constraints and the uncertainty of the market price. The deregulated electricity market usually has a few generators (market suppliers) that usually dominate the market. This makes the market seem more similar to an oligopoly. In such oligopolistic markets, an individual generator can exercise market power and manipulate the market price via its strategic bidding behavior [8]. Companies have to determine bidding strategies so that they can profit even if they are price takers and do not dominate the market. 6 There are several approaches to analyze and develop bidding models. The bidding strategies used in the market are discussed in [9] which contains a literature survey of the current approaches to the bidding problem. The performance of a power market is measured by the common term ?social welfare?. Social welfare is the benefit of a commodity to society: to both customers and suppliers. The game theory approach is commonly used in the literature to model market participant?s behaviors [10]. The approach assumes that each market player tries to maximize its profit. The behavior of the market player is affected by other players? behavior. Several methods used in modeling bidding strategies are explained in [11] as they compare the game theory approach with the conjectural variation based method. In both approaches, each firm in the market rationally tries to maximize its profit while considering the reactions of its competitors. They show that firms can increase their profit by using conjectural variation based method and the equilibrium found corresponds to Nash equilibrium. Several solution methods have been used to solve the bidding problem. In [12], a genetic algorithm is used to solve the bidding problem. Although the solution obtained is a heuristic one, it could be used by a company in its daily bidding process. The authors explain the process of bidding and how the equilibrium price is determined. They construct their model based on the assumption of exogenous prices. In [13], an optimization tool to determine a bidding curve for the Ontario Power Market was developed. The authors used different scenarios of market prices and load. The decision process for the generator is based on the probability distribution of forecasted prices. The model chooses the block of the curve (the price and the corresponding quantity) to be 7 submitted in order to maximize expected profits. The model assumes exogenous prices and includes operational constraints such as ramp-up limits, start-up limits and minimum up-down times. In [14], bids are represented as quadratic functions of power levels. The model optimizes the parameters of each function during a two-phase process. In the first phase, the ISO minimizes the total system cost in which the parameters for other generators are known. In the second phase, solutions are plugged into the generator?s model. The Lagrangian relaxation procedure is used to solve the expected cost minimization problem. In [15], the authors assume that suppliers bid linear supply functions; the coefficients of the functions are chosen for each supplier in such a way that the expected profit is maximized subject to the behavior of one?s rivals. They formulate a stochastic optimization model and use a Monte Carlo based method to tune the parameters of the function. They also include the level of information known for each generator in a symmetric and asymmetric market. In [16], [17] and [18] bidding strategies are developed for price taker generating units. The papers mentioned above deal with bidding problems that consist of bids of one or two blocks. The models developed in this dissertation consider that companies submit up to ten blocks in their bids and include multiple price scenarios. These additional features increase the computational time required to solve these problems. In addition, in some papers the bids are modeled as a price-quantity function. Thus, the optimization consists of finding the coefficients of this function. In contrast, the price- quantity is function free in this dissertation. 8 2.2 Background In order to better analyze the SBP, one should understand the fundamentals of auction theory, power market, electricity pricing, cost of electricity production and bid structures. The fundamentals might not be exactly the same for all power markets, but they are similar. Auction based markets are also used in some other markets such as the stock exchange, agricultural wholesale, and goods wholesale on the Internet. The solution approaches developed in this dissertation could also be applied to such markets. 2.2.1 Auction Theory Auction theory deals with how participants of an endeavor behave in auction markets where game-theoretic behaviors are involved. The objective of an auction is generally determined by an operator. It can be the maximization of the outcome (revenue) like government licenses, or it can be the minimization of the cost like public service by giving equal market opportunity to each competing player. The players, rules, outcomes and payoffs of the auction along with their mission might change depending on the objective, but the same essence of competition remains [19]. An auction design requires careful research and experiments on efficiency, optimum and equilibrium bidding strategies and revenue so that an effective auction can be created and manipulations can be eliminated. It is important that each player does not have perfect information about their competitors? bid. The operator is supposed to provide a confidential environment to provide equal opportunity for players [19]. 9 Typical auctions are classified by designs (rules). Some examples are first-price sealed-bid auctions, second-price sealed-bid auctions, open ascending-bid auctions (English auctions) and open descending-bid auctions (Dutch auctions). In first-price sealed-bid auctions, bids are submitted simultaneously in sealed envelopes by all bidders to the operator. The individual with the highest bid wins and pays the proposed amount. In second-price sealed-bid auctions, bids are submitted in sealed envelopes simultaneously. The individual with the highest bid wins and pays the amount equal to the second highest bid. In open ascending-bid auctions, the price is steadily raised by the operator. Some players drop out as the price becomes too high and the last player wins the auction at the current price. This is more common in revenue maximizing auctions. In open descending-bid auctions, the price starts at a relatively high level and is steadily decreased by the operator until one player is willing to accept the offer. This approach is more common in cost minimization auctions which often are an issue for public service providers. Many auctions are hybrids of these four types. Since Electric power auctions aims to minimize the cost of electricity provided to the market, they are a hybrid of open descending-bid auctions and first-price sealed-bid auctions [20]. The success of a competetive market is related to the design of the auction mechanism used [19]. Auctions can also be classified as single-round auctions and multi-round auctions. In single-round auctions, sell bids are matched with buy bids to reach equilibrium at once. On the other hand, in multi-round auctions bidders are asked to update their offers at each iteration/round so that a more effective equilibrium can be reached [21]. 10 2.2.2 Market Design Federal Energy Regulatory Commission proposed a design on July 31, 2002 titled ?Standard Market Design? (SMD) for the standardization of electric power markets in USA [22], [23]. The major power markets in the US such as PJM, New England (NEEPOL), New York (NYPX), and Midwest ISO (MISO) are variants of the SMD [19]. The Electric Reliability Council of Texas (ERCOT) and California power market (CAISO) are in the process of implementing SMD rules as well. The objective of a typical SMD is to develop a market structure that brings together the physical system and the economic financial operations. This is achieved by defining the roles and the interaction of system components. SMD also deals with the system governance, market operations, risk management, market monitoring and conflict resolutions that might occur among the members [23], [24]. Figure 2.1 illustrates general services involved in SMD [24]. ENERGY Settlements, Auctions, Dispatches MARKET MONITORING Price Caps, Market mitigation INVESTMENT Transmission, Generation Expansion TRANMISSION Congestion, Transmission, Hedging ANCILLARY SERVICES Auction Structure RETAIL SERVICES Demand, Customers Figure 2.1. Services in a typical SMD 11 An electric power system consists of four main parts: generation, transmission, distribution and customers. Figure 2.2 shows the flow of electricity from generation to customers [24]. SMD governs the processes through scheduling coordinator (SC), power exchange (PX), independent system operator (ISO), transmission owner (TO) and load serving entity (LSE). The role of the ISO might slightly change in each SMD, but the major role remains the same. Competition Generation Transmission Distribution Demand Federally Regulated (FERC) State Regulated Customers Figure 2.2. Electric power system The ISO is a neutral entity responsible for maintaining the instantaneous balance of the grid system [24], [25]. It performs its function by controlling the dispatch of flexible plants to ensure that loads match resources available to the system. It also coordinates the day-ahead market and real-time balancing market and monitors compliance with all regional operating and reliability standards. The members of the ISO 12 are power suppliers, wholesale power customers and, transmission line owners. The day- ahead market and real-time market are used to equate supply and demand based upon the sell and buy bids submitted by the members [25]. In a typical SMD, pricing is handled by a nodal mechanism [24]. 2.2.3 Locational Marginal Pricing Locational Marginal Pricing (LMP) is a market-pricing method that is used to manage the efficient use of the transmission system when congestion occurs between source and sink. In a SMD, congestion occurs when restrictions such as capacity of the line and losses prevent the economic or least expensive supply of energy from serving the demand. It also means that, if the system was entirely unconstrained and there were no losses, all the LMPs would be same and it would reflect only the energy price [26], [27]. In an SMD, after offers and bids are submitted and the market is settled, the LMPs are usually calculated at three types of locations, at the node, the load zone and the hub. Nodes are the places on the system where generators inject power into the system or demand (load) withdraws from the system. Each node is connected to one or more buses, which are specific components of the power system at which generators, loads or the transmission system are connected. Prices are made up of three components energy, congestion and losses. The energy component is the cost to serve the next increment of demand at a specific location or node which should be produced by the least expensive generating unit in the system that still has available capacity [26], [28]. If the transmission network is congested, the next increment of energy cannot be delivered from the least expensive unit on the system because it would violate transmission operating 13 criteria or cause overloading. The extra cost which is called congestion cost is calculated at a node as the difference between the energy component of the price and the cost of the providing the additional more expensive energy that can be delivered at that location [24],[25]. Losses occur during the transmission of power from one location to another and incurred costs are also included to the calculation. Generators are paid nodal LMPs and assured by market rules to recover their costs. Load zones consist of aggregations of nodes in a region. SMD requires a load to pay the price calculated for that particular zone. The prices calculated for each zone are a load-weighted average of the nodal prices located within each zone. They are less volatile than nodal prices since they are aggregated from nodes that reflect the true cost for delivering power by different location [25]. 2.2.4 Cost of Generation The cost of the energy produced by the generating unit depends on the amount of fuel consumed and is typically approximated by a quadratic cost function 2321)( qaqaaqC ++= ($/MWh), where q is the amount of energy generated in one hour [27],[29],[[30]. The coefficient a1 represents the fixed cost or no-load cost for each hour. This cost includes the labor and the cost of non-direct goods necessary to produce power for that hour. The value a2 represents the linear cost which is proportional to the amount of power produced. The parameter a3 is the quadratic cost coefficient and it is related with the amount of fuel used to produce electricity [29]. 14 Generators use a single cost function when they bid into market but it is also possible to combine number of cost functions and get an approximated cost function. Most of the time this approach is used by firms which have several generators and prefer to bid into the market using portfolio-based cost functions. Generators also have start-up costs, minimum-load operating costs and minimum up-down constraints. These are also used by the ISO when the bids are evaluated [31]. 2.2.5 PJM Power Market The PJM interconnection is a federally regulated and nonprofit organization that manages the transmission of wholesale electricity in Pennsylvania, New Jersey, Maryland (PJM), Delaware, Illinois, Indiana, Kentucky, Michigan, North Carolina, Ohio, Tennessee, Virginia, West Virginia and the District of Columbia, involving more than 51 million people. Its dispatching capacity is more than 164,000 MW [32]. PJM?s members, totaling more than 450, include power generators, transmission owners, electricity distributors, power marketers and large consumers. PJM?s role as a federally regulated RTO means that it acts independently and impartially in managing the regional transmission system and the wholesale electricity market. The PJM energy market includes two markets ? day-ahead and real-time markets. In addition to these markets, there is daily capacity market, Monthly capacity market, fixed transmission rights (FTRs) auction market, regulation market and spinning reserve market. In the day-ahead market, bilateral transaction schedules, generator offers, and consumer demands are submitted twelve hours before the actual delivery of electricity. 15 Each installed capacity in the day-ahead market has an obligation to submit an offer to the market even if it is unavailable or in outage [24], [25]. All participants must submit bids offers until 12:00 p.m. for the next operating day. The ISO evaluates the bid offers between 12:00-4:00 p.m. No offer is accepted during this time. PJM announces the accepted bids at 4:00 p.m. Non-winning participants have the chance to modify their bids until 6:00 p.m. Demand bids also follow the same process for the day-ahead market [32]. Based on these offers and demands, market clearing prices are determined for each hour of the next operating day. The day-ahead market is considered a forward market because the formation of the generation and consumption is determined the day before the operating day [25]. On the other hand, the real-time energy market balances the deviations occurred in the day-ahead market and the actual generation. Unlike the day-ahead market, market clearing prices in the real-time market are calculated every 5 minutes based on the actual system operations. The methods developed in this dissertation assume that the bids are submitted to the PJM power market. 2.2.6 Bidding into the PJM Market A generator offer for the PJM market is composed of two components, the price and quantity of electricity that a supplier is willing to generate. Offers are submitted in blocks of price quantity pairs. PJM allows submitting at the most ten blocks for a generator offer [32]. Figure 2.3 illustrates a valid offer curve in PJM power market. 16 0 50 100 150 200 250 300 0 50 100 150 200 250 300 Quantity (MWh) Pr ice ($ /M W h) Figure 2.3. A generator?s offer curve in the PJM day-ahead market Each generating unit also submits its minimum run time, minimum down time, no-load costs and start-up costs to the PJM market [32]. PJM runs the ?two-settlement? software to determine the hourly commitment schedules and the LMPs. Generating units that have minimum run times that exceed 24 hours are asked by PJM to submit binding offer prices for the next seven hours. PJM supports mainly three offer types: Cost-capped offers, historic LMP based offers and market-based (or price-based) offers. Each unit submits its operating constraints in its profiles however its usage differs for each offer type [32]. Cost-based offer in the PJM market consists of the incremental operating cost of the generation resource, plus a 10% and, plus variable operations and maintenance costs. The production cost method is the method of capping. Generators use no-load cost and start-up costs as classified hot, intermediate and cold in their offers. Generator offers consist of three parts, offers for energy ($/MWh), start-up ($/day) and no-load ($/hr). 17 Historic LMP capped offers are determined by calculating the average LMP at the generation bus during all hours over the past six months in which the resource was dispatched above minimum. Generators also use start-up cost and no-load cost in their offers. Table 2.1 shows the summary of the offer characteristics [32]. Table 2.1. Comparison of offer types Offer Component Cost-Capped Offers LMP-capped offers Market-based offers Energy offer ($/MWh) Production cost, Average historic LMP No cap, support up to 10 blocks Start-up cost ($/day) Production cost, plus 10% Production cost, plus 10% No cap No-load cost ($/hr) Production cost, plus 10% Production cost, plus 10% No cap Price-capped offers are offers that are not necessarily capped with costs. Firms bid on the market price and get paid the price determined with the respective auction procedure. The start-up cost and no-load cost are still submitted to the PJM but their usage is not same as the former two offer types. However, cost components for cost- capped offers and historic LMP-capped offers directly impact their selection chance since if cost recovery is not possible for those generator, ISO chooses cheaper offers. Figure 2.7 summarizes the characteristics of three types of offers including differences [32]. 18 Table 2.2. Bid components in offer types Component Cost & LMP-capped offers Market-based offers Start-up Daily- 3 types Optional- every 6 months No-load Yes Optional- every 6 months Maintain Minimum Yes No Cooling Requirement Basic Optional- every 6 months Incremental Cost ? Energy from min to max ? Energy from 0 to max Maximum Offer Based on cap methodology $1000/MWh Notice that in Table 2.2 there is a cap for start-up and no-load costs for cost- capped offers, while there is no-cap for market-based offers. This is because market- based offers bid on the price while cost based offers bid on their cost. PJM audits the costs if it is necessary. The unit?s offer type is initially set by the generator owner to indicate whether the unit is to be scheduled as a market-based offer or a cost-based offer. PJM also imposes the rule if the generator once chooses market-based offer and it cannot switch to another offer type [32]. SMDs usually use uniform-price auctions and pay-as-bid auctions to govern the market mechanism. After bids are submitted and the market is settled by the ISO, all dispatched generators in the uniform price auction are paid the market price where as they got paid their bid price in pay-as-bid auction. The selection process for winning generators and the equilibrium price are the same for both designs with the difference that the generators would make different revenues. 19 The pay-as-bid auction method motivates the lower cost generators to bid too high to increase the market price. However, in uniform price auctions generators need to bid lower than market price in order to be selected. Thus, it is accepted that uniform price auctions generally end up with lower market prices because of price pressure [31]. 2.2.7 Power Market Equilibrium SMD aims to increase competition and hence it is a good place where suppliers and consumers meet under the supervision of an ISO and economic fundamentals. The balancing of supply and demand is always crucial in an economic market. However it is vital for an electricity market since the lack of electricity when needed can cause very costly consequences. After wholesale power suppliers bid sell offers and wholesale power customers bid demand offers into the market, the next step is to find an equilibrium point where supply and demand meet. The price at this point not only sets the market clearing price but also sets the market clearing quantity that makes social welfare an optimum. A power market is said to be in equilibrium if i) all suppliers maximize their return at the given market clearing price and market determined dispatch schedule ii) all consumers maximize their utility at the given price and schedule iii) the total supply equals total demand [33]. Equilibrium approaches for power markets are constructed on the behavior of system participants, which involve the game theory approach. In [34], authors classify the modeling approaches of power market problems. They use exogenous price and demand- price function methods to model optimization problems for a single firm. 20 Simulation models that consider all firms behaviors? are classified into equilibrium models and agent-based models. There have also been efforts to model power market equilibrium as Cournot, Bertrand and Supply Function Equilibrium (SFE). The SFE has been of most use in power market modeling [31], [35]. Cournot models have not been found satisfactory for the power markets since quantity produced by each player, the decision variable, is not responsive to the effects of price sensitivity. The cournot model also has the assumption that the residual demand is elastic. But it is not considered as an issue for electricity markets [31], [35]. Cournot models also expect that player?s output will not change the rival?s output, but the offered quantity along with its price by the market participant actually affect the rival?s expected revenue in the market. On the other hand, Bertrand model which has the market price as a decision variable requires each player to build its strategy on expected market price. Bertrand leads to perfect competition if there were player with unlimited capacity. But it is known that as the market converges to equilibrium the players with small capacities will be out of the game. This strategy also does not cover all market?s issues [31], [35]. Klemperer and Meyer first developed SFE [36] and showed that each firm can express its decisions in terms of a quantity and a price in the absence of certainty and having an idea about competitors? strategic variables. It is after that SFE was applied to power markets by Green and Newbery showing that if a firm expose its decision tool in a form of supply function indication prices at which it is willing to offer various quantities to the market for a given demand curve, it can expect in general a greater profit in return [37]. SFE is more accurate comparing to former mentioned approaches since it reflects the bidding rules in SMD where players submit price-quantity offers as decision 21 variables. In [38], authors analyze SFE applications by assuming that the functions provided are linear and there are strategic players in the market with price and capacity cap. They show that SFE will be more effective than the other approaches in terms of representation and reaching to the equilibrium. In [39], authors work on SFE by modeling the market players as non-degreasing supply function providers and competing in a game. In [40], authors analyze the market power by modeling the equilibrium for large scale power systems. They briefly explain the equilibrium approaches and show that SFE can fit best to the power market. In [41], it is showed that if there are multiple players with identical marginal costs and asymmetric capacities, a unique piece-wise symmetric supply function exists. The authors show that small players will be eliminated at some point and larger players will use this advantage in their supply functions. In [31], authors analyze the equilibrium models in terms of transmission network, generator cost function and operating characteristics, bidding, demand and uncertainty. These are evaluated under the umbrella of economic, physical and commercial modeling. 22 CHAPTER 3 SPREADSHEET BASED SIMULATION FOR ELECTRIC POWER OFFER EVALUATION 3.1 Bid Simulator In order to evaluate a bidding strategy for given market price scenarios, a simulation method needs to be developed. The simulation method should include different price samples and should be able to work for different cost functions. We develop a simulation model called Bid Simulator in Excel that includes all parameters needed to evaluate a bidding curve. Bid Simulator is a spreadsheet based simulation model to assess the value of a given bid using market price samples. It provides supportive statistical outputs to the decision maker. The simulation model includes market price scenarios and calculates hourly profits according to the market prices. The pseudo code for the simulation is given in Figure 3.1. If the market price at a particular hour is larger or equal to any given price bid, the supplier would sale power. Otherwise it would not sell power at that hour. 23 Figure 3.1. Pseudo code of the simulation model 3.2 Simulation Functionality The model is implemented in Excel. The simulation spreadsheet consists of five modules: main screen, simulation, hourly profits, daily profits and output analysis. Figure 3.2 and Figure 3.3 show the main screen and menu for the model respectively. Determine bid prices and quantities Generate N Price Samples each for 24 hours For each N For each hour For each block in bidding curve If bidding price <= Market price Calculate hourly profit Else Hourly profit=0 Next block Next hour Daily Profit = Sum (Hourly profit) Next sample Average Profit = Sum (Daily Profit )/N 24 Figure 3.2. Bid Simulator main screen Figure 3.3. Bid Simulator menu 25 The simulation model generates up to 1000 market price samples for 24 hour, i.e. a total of 24000 market prices. The price sampling can be defined in terms of a probabilistic distribution of interest before the simulation. Figure 3.4 shows the simulation model. The results shown are for a generator whose cost function is 20.004245)( qqqc += and maximum capacity 300 MW. The market prices are generated using a normal distribution with mean equal to day-ahead prices of April 20th , 2008 of the PJM market and standard deviation 2 $/MWh. Figure 3.4. Bid Simulator simulation and price generation For the hourly profits module, a given bid is evaluated for each hour using the sampled prices. Figure 3.5 shows the hourly profits. 26 Figure 3.5. Bid Simulator hourly profits Daily profits are calculated for each day using hourly profits in the daily profits module of the package, i.e. 1000 profits. Figure 3.6 shows the hourly profits. Figure 3.6. Bid Simulator daily profits 27 In the output analysis module, the descriptive statistics are calculated for the hourly generated prices, hourly profits and daily profits. Figure 3.7 shows the output analysis. The output analysis includes the following components. Statistics on hourly generated prices: -Expected market price for each hour - Expected hourly market price graph Statistics on hourly profits: - Expected hourly profit - Expected hourly profit graph Statistics on daily profits: - Minimum daily profit - Maximum daily profit - Expected daily profit - Standard deviation of daily profits - Variance of daily profits - 5% and 95% percentile of daily profits - 5% confidence interval on mean - Probabilistic distribution of profits - Histogram of probability density function (PDF) - Histogram of cdf (cumulative distribution function) 28 Figure 3.7. Bid Simulator output analysis 3.3 Numerical Example We evaluate an offer using the Bid Simulator and we compute the probabilistic distribution of profits. The descriptive statistics and graphs mentioned above indicate to the decision maker how good the bidding solution is. We evaluate the bid given in Table 3.1. Figure 3.8 gives expected hourly prices and Figure 3.9 gives expected hourly profits. Table 3.2 gives the descriptive statistics. A 90% confidence interval mean profit lies between $55,934 and the $56,207. Figure 3.10 gives the plot of the histogram and Figure 3.11 gives the cumulative distribution. 29 Table 3.1. An example bidding solution to the problem Block 1 2 3 4 5 6 7 8 9 10 bi 45.12 45.50 45.75 46.00 46.26 46.51 46.76 47.01 47.26 47.52 qi 30 60 90 120 150 180 210 240 270 300 0 10 20 30 40 50 60 70 80 90 100 1 3 5 7 9 11 13 15 17 19 21 23 Hour Price ($/MWh) Figure 3.8. Expected hourly price graph from generated prices Table 3.2. Descriptive Statistics for Daily Profits Maximum Profit ($) 64,066.00 Minimum Profit ($) 49,248.06 Expected Profit ($) 56,070.85 Standard Deviation 2,204.29 Variation 4,858,884.50 95% Percentile ($) 59,652.06 5% Percentile ($) 52,335.05 30 0 2000 4000 6000 8000 10000 12000 14000 16000 1 3 5 7 9 11 13 15 17 19 21 23 Hour Profit($) Figure 3.9. Expected hourly profits based on generated prices 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 47,248 54,657 62,066 Probability Profits ($) Figure 3.10. Histogram of daily profits 31 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 47,248 54,657 62,066 Profits($) Cumulative Probability Figure 3.11. Cumulative histogram of daily profits 3.4 Conclusion In this chapter, the use of simulation to evaluate a bidding curve was explained. Its effectiveness was showed using a numerical example. A decision maker can use the simulator to estimate expected profit and variation under the market price uncertainty. The tool also gives expected profit and price for each hour which can be used as an input in unit scheduling problems. Especially if an efficient market price forecasting method is available, Bid Simulator method can be used to help decision maker to reach more accurate results. Different offers can be compared to further analyze the sensitivity. The Bid Simulator will be used in next chapters to verify the results found in the analysis. 32 CHAPTER 4 ELECTRIC POWER BIDDING UNDER PRICE UNCERTAINITY 4.1 Price Uncertainty Electricity is generally accepted as different from other commodities. It is not storable and its demand is instantaneous so it must be produced and used in real time. These unique characteristics of electricity and necessity of real time balance create a need for coordinated markets. As explained in previous chapters, LMP create diversified market prices by location. The price is strongly load-dependent, highly volatile, seasonal and consumption dependent [16], [17], [42]. The parameters are stochastic which gives a stochastic behavior to the electricity price [43]. Energy consumption, fuel costs, availability of fuels, equipment capacity and market participants? behavior are stochastic [31], [44],[49]. Figure 4.1 shows an example of day-ahead market prices and real time market prices for one day in the PJM power market. 33 0 20 40 60 80 100 120 140 160 180 200 0 2 4 6 8 10 12 14 16 18 20 22 24 Hour Ma rk et Pr ice ($ /M wh ) Day Ahead Real Time Figure 4.1. Day-ahead and real time prices for February 29, 2008 In the literature there are many models that are used to forecast market clearing prices. In [45], the authors develop wavelets and multivariate time series based price modeling. They analyze the market price data statistically to determine the model parameters. In [46], the authors use artificial neural networks to model the market price changes. Neural networks are useful to reflect nonlinear changes that are difficult to predict otherwise. Automatic dynamic harmonic regression model is used in [47] to handle the regression among market prices. In [48], the authors develop an algorithm to calculate the mean and variance of the electricity market price. They also give a stochastic method for load estimation. The bidding models that consider the market price exogenous usually include the electricity price as an input. Market price forecasting methods can be used to determine the prices used in the model. As an alternative, electricity prices can be generated based on a probability distribution function. 34 The market price scenarios then can be included into the model and the developed bidding strategy is the offer that maximizes the expected profit across all scenarios. 4.2 The bidding model We consider an electric power producer with a set of generating units and assume that the producer wishes to submit an offer curve to the day-ahead market for each of its units. We assume that the producer is a price-taker market player that obtains its revenues by selling power at the market clearing prices of the PJM pool. That is, if the power supplier produces electricity with its generating units at a particular hour, it is then willing to take the price prevailing in the market at that hour. We also assume a lack of market power, i.e. the power supplier does not perceive its decisions as affecting market prices. The SBP is formulated under the following additional assumptions: i) An offer or bid, which consists of N price-quantity blocks at the most, needs to be determined for each generator separately. ii) The market clearing prices are considered exogenous to the model, i.e. they are not affected by the bidding decisions of the unit for which the model is being solved. iii) The market clearing price at each hour is considered a random variable whose probability distribution has known parameters. iv) The offer is determined before the market closes at 12 noon and is valid for the next twenty four hours, starting at 12 midnight the same day. 35 In finding an optimal offer curve for a particular generator (see Figure 4.2), there are N pairs of decision variables ib and iq? (i=1,..,N) that need to be determined. Figure 4.2 shows the relations of ib and iq? (i=1,..,N) in the PJM market. The variable iq? denotes the amount of energy increase to get the bid price ib for delivery at any hour of the next day. These values are represented by the vectors?q and b, respectively. If the market clearing price at hour t is equal to or higher than the offered price ib , then all energy blocks offered at this price or lower are accepted by the market operator. Thus, the total energy to be produced at time t and sold to the market at a price Pt is given by: ? = ?= )( 1 tPI i it qq , where tjt PbjPI ?= such that Max )( for t=1..T; i=1..I(Pt) (4.1) Figure 4.2. Bid prices and power quantity intervals 36 The objective function is total profit (revenue minus generation cost) over a 24- hour period. The revenue during hour t is obtained from selling the quantity stipulated under the offer (only if the generator is dispatched) at the market price, Pt ($/MWh). The cost includes those of producing the energy. As mentioned above, Pt is a random variable, and therefore the total profit over a period of T hours is also a random variable. We assume that K samples of the hourly prices are available and they have equal probability of occurrence. We denote the price at hour t of sample k as Ptk. Thus, the objective is to maximize the expected profit over the time period T (usually 24 hours). The bidding problem, which is called ),P( b?q , is stated as follows: ])([1E[Profit] Max),P( 1 1, ?? = = ?== K k T t k t k t k t qCqPKb?qb?q for k=1..K; t=1..T; (4.2) Subject to the following constraints: max 1 Qq N i i ??? = (4.3) max0 Bb i ?? for i=1..N. (4.4) max0 Qq i ??? for i=1..N. (4.5) ? = ?= )( 1 ktPI i i k t qq , ktikt PbiPI ?= such that Max )( for k=1..K; t=1..T; i=1..I(Pt k). (4.6) 2 321 )()( k t k t k t qaqaaqC ++= for k=1..K; t=1..T. (4.7) 37 4.3 Quadratic Programming Model Mathematical Programming is one way to find an optimal solution to the bidding problem. However, it can solve relatively small sized problems. By setting the number of samples to the one and the number of maximum bidding blocks equal to the number of hours of the time horizon, a quadratic programming model can be formulated and solved using a commercial software package such as Cplex. Notice that when the market price consists of one sample and the number of blocks is equal to the number of hours, the optimal bidding price of a block of power is equal to one of the market prices. Therefore, the bidding problem ),P( b?q reduces to the following mathematical representation: Max Z = ? = ?? T t tttt qaqaqP 1 2 32 ][ for t=1,..,N (4.8) Subject to ? = ?= t i it qq 1 for t=1,..,N (4.9) max 1 Qq N i i ??? = for i=1,..,N (4.10) 0?? iq for i=1,..,N (4.11) 4.4 Nonlinear Programming Model Nonlinear Programming (NLP) is the process of solving a problem that includes equalities, inequalities, constraints and an objective function some of which is nonlinear. 38 The process finds a set of unknown real variables that makes the objective function maximized or minimized. The bidding model is formulated as follows: ])([1E[Profit] Max),P( 1 1, ?? = = ?== K k T t k t k t k t qCqPKb?qb?q (4.2) Subject to the following constraints: max 1 Qq N i i ??? = (4.3) i k t k t bPiMz ?? 001.1)( for i = 1..N ; t = 1..T ; k = 1..K. (4.12) i k t k t bPizM ??? )1)(( for i = 1..N ; t = 1..T ; k = 1..K. (4.13) 0)1()( ?+? iziz ktkt for i = 1..N ; t = 1..T ; k = 1..K. (4.14) 1+? ii bb for i=1..N. (4.15) max0 Bb i ?? for i=1..N. (4.4) max0 Qq i ??? for i=1..N. (4.5) ii qMr ???1 for i=1..N. (4.16) )1(max ii rQq ??? for i=1..N. (4.17) 1+? ii rr for i=1..N. (4.18) for i = 1..N ; t = 1..T ; k = 1..K. (4.19) and (4.7). ? = ?= N i i k t k t qizq 1 )( 39 Note that M is a large number here, and z and r are binary variables. When we model the problem in AMPL and solve it using the NEOS solver MINLP, an optimal solution can be found for a limited number of price samples. 4.5 Marginal Cost Bidding A power producer also can submit its marginal cost of production as its bid offer. As a matter of fact, in a perfectly competitive market it is expected that each player submits its marginal costs. To do so, a power producer could offer, for each generating unit, one energy block consisting of the maximum capacity and price equal to the marginal cost of producing this amount. Alternatively, the power supplier could split the maximum capacity into N blocks of identical size and offer them at prices equal to the marginal costs of producing each block. PJM accepts a maximum of ten energy blocks in its daily bidding process, so maximum capacity can be split into 10 blocks and marginal cost of these quantities can be offered to the market. ])([1E[Profit] ),P( 1 1 ?? = = ?== K k T t k t k t k t qCqPKb?q for k=1..K; t=1..T; (4.2) Subject to the following constraints: N Qq i max =? for i=1..N. (4.20) k ti qaab 32 2+= for k=1..K; t=1..T; i=1..N. (4.21) (4.4), (4.6) and (4.7). 40 4.6 Numerical Example and Analysis We now present results obtained with the developed models. We use CPLEX for solving the quadratic programming model, NEOS MINLP for solving the nonlinear programming problem. 4.6.1 Quadratic Programming Model In order to solve the quadratic programming problem we use the generator used in Chapter 3 whose cost function is 20.004245)( qqqc += and maximum capacity 300 MW. The time horizon is set to ten hours. The ten hourly market prices are given in Table 4.1. Table 4.1. Price sample for the example problem Hour 1 2 3 4 5 6 7 8 9 10 Price($/MWh) 35.20 36.80 39.30 45.00 45.50 45.90 46.10 46.80 47.10 50.20 After solving the above model using Cplex 7.0, the optimal profit is found to be $1772.48. The optimal solution to the problem is given in Table 4.2. Table 4.2. Optimum solution to the problem Block 1 2 3 4 5 6 bi ($/Mwh) 45.50 45.90 46.10 46.80 47.10 50.20 qi (Mwh) 59.52 107.14 130.95 214.28 250.00 300.00 4.6.2 Nonlinear Programming Model In order to solve the bidding problem under market price uncertainty, we use the same generator used in section 4.6.1 but the problem is solved for three day price 41 samples. The problem was coded in AMPL and submitted to the one of the NEOS Servers MINLP to get a solution. Table 4.3 gives the price samples used in the model and Table 4.4 provides the optimum solution. Table 4.3. Price samples used in the problem Hour 1 2 3 1 49.21 49.58 47.14 2 47.56 46.78 46.52 3 50.44 40.68 46.1 4 43.48 44.59 46.59 5 41.97 46.87 46.16 6 43.78 48.8 49.83 7 47.25 44.27 45.56 8 57.28 54.37 55.42 9 57.56 53.37 53.29 10 58.59 52.87 55.02 11 51.66 48.85 46.39 12 45.54 48.36 49.36 13 40.14 38.55 42.55 14 37.39 40.44 38.28 15 37.52 42.54 38.82 16 39.45 34.4 37.2 17 40.64 44.53 40.38 18 53.45 53.41 53.71 19 79.44 77.53 75.92 20 72.7 75.02 74.5 21 71.83 70.09 67.66 22 60.39 62.78 64.41 23 48.94 50.19 50.53 24 40.8 44.97 45.43 Table 4.4. Optimum solution to the problem Block 1 2 3 4 5 6 7 8 9 bi ($/MWh) 45.01 45.54 46.10 46.16 46.39 46.52 46.56 46.87 47.56 qi (MWh) 51.19 65.47 130.94 138.08 165.46 180.93 211.88 222.59 261.28 42 The objective function of the optimum solution was $44,779.83. However, it takes about 5 hours to solve the problem with 3 price samples and 300 Mwh capacity. If we increase the capacity to 1500 Mwh and try to solve the problem with same price samples, we could not find an optimal solution after 24-hour run. Results show that it is not likely to solve the problem with more than 3 price samples. 4.6.3 Marginal Cost Bidding The Marginal cost bidding model requires splitting the maximum capacity into equal block sizes. We solve the problem using the same generator with the price samples given in Table 4-3. Table 4.5 gives the bid prices and quantities for marginal cost bidding. Table 4.5. Bid prices and quantities for marginal cost bidding Block 1 2 3 4 5 6 7 8 9 10 bi ($/MWh) 45.25 45.50 45.75 46.00 46.26 46.51 46.76 47.01 47.26 47.52 qi (MWh) 30 60 90 120 150 180 210 240 270 300 Marginal cost model is evaluated for both price samples given in Table 4.1 and 4.3. The profit found for 10-hour price sample is $1766.58 where the optimum solution is $1,772.48. The profit found in for 3 day price samples is $44,750.86 where the optimum solution is $44,779.83. 4.7 Conclusion In this chapter the strategic bidding model was defined and possible solution approaches and their limitations were explained. One can solve a problem with 10 hours of market prices using quadratic programming. We also showed that it is possible to find 43 an optimum solution using 3 days price samples for a generator with 300 Mw capacity. It is clear that a more effective method is needed to solve the problems with more than 3 days price samples and more generating capacity. The solution method should require low computational time since the bidding process is done daily. It might be that in same cases there are not much profit difference between bidding the marginal cost and optimum solution. We will show in Chapter 5 that this is not always the case. 44 CHAPTER 5 PARTICLE SWARM OPTIMIZATION FOR BIDDING INTO MARKET UNDER PRICE UNCERTAINITY 5.1 Rationale for Heuristics Approach Heuristic approaches are commonly used when it is not possible to find optimum solutions to the problems, or when it takes much computational time to find optimum solutions. The approaches are generally accepted as easy to implement, easy to apply to the problems and require less computational time. They do not guarantee an optimum solution though. But in cases where an optimum solution and a good solution don?t make much difference, heuristics approaches are preferred because of the less effort that they require. It is possible to find a solution using more than 3 days price samples using heuristics approaches whereas nonlinear programming could not find an optimum solution. We will show in this chapter that it is possible to find a good solution for every generator regardless of generator capacity. It will also be showed that the required computational time is dramatically less than the nonlinear programming method. 5.2 Particle Swarm Optimization Particle Swarm Optimization (PSO) is a computation technique, introduced by Kennedy and Eberhart in 1995 [50], which was inspired by social behavior of bird 45 flocking or fish schooling. Like genetic and evolutionary algorithms, PSO is a population-based search method, i.e. it moves from a set of points (particles? positions) to another set of points. The particles move through a D-dimensional space and each particle has a velocity that acts as an operator to obtain a new set of individuals. The particles adjust their movements depending on both their own experience and the population?s experience. The following pseudo code describes the PSO approach: Figure 5.1. Pseudo code for the PSO method (The symbol ? denotes the multiplication of two vectors component by component.) 46 At each iteration a particle moves in a direction computed from its best visited position and the best visited position of all the particles in its neighborhood. Among the several variants of PSO, the global variant considers the neighborhood as the whole population, called the swarm, which enables the global sharing of information. The basic elements of the PSO technique are particle, population, velocity, inertia weight, individual best, global, learning coefficients, and stopping criteria best [51]. These are briefly discussed below. I. Particle, Xj(t): a particle j represents an m-dimensional vector candidate solution. The value of m is determined by the number of decision variables. At time t the particle j can be described as Xj(t) = [x1,j,?, xm,j] where the x components are the decision variables. A value of xi,j denotes the position of particle j in the ith coordinate in the search space, i.e. the value of the ith decision variable in the candidate solution j. II. Population, POP(t): The population is a set of n particles at any given time t and it can be represented as POP(t)=[ X1(t) ?, Xn(t)]. III. Velocity, Vj(t): The velocity of moving particles at time t represented by an m-dimensional vector Vj(t) = [v1,j,?, vm,j]. IV. Inertia weight, w: The parameter that controls and directs the impact of the previous velocities on the current velocity. If the inertia weight is large, the search becomes more global, while for smaller w the search becomes more local. V. Individual best, )(* tjX : When a particle flies through the search space it compares its fitness value at the current position to the best fitness value so far. The best position visited by the particle i at time t is denoted by ],...,[)( *,*,1* jmjj xxt =X . 47 VI. Global best, G(t): Represents the best position that gives the best fitness among all individual best positions achieved so far. It is defined by G(t) = [g1,?, gm]. VII. Learning factors, c1 and c2: these coefficients help particles to accelerate towards better areas of the solution space. VIII. Stopping Criteria: represent the conditions for which the search process will terminate and lead to a result. 5.3 Conventional Particle Swarm Optimization The conventional PSO approach (cPSO) is used for solving the whole problem. The population size is set to ? particles. Decision variables bi and iq? are evaluated using the price scenarios in fitness function below. ])([1]Profit[EMax 1 1 ?? = = ?= K k T t k t k t k t qCqPK for k=1..K; t=1..T; (5.1) Subject to: max0 Bb i ?? for i=1..N. (5.2) maxQqk t ? for k=1..K; t=1..T (5.3) ? = ?= )( 1 ktPI i i k t qq , ktikt PbiPI ?= such that Max )( for k=1..K; t=1..T; i=1..I(Pt k). (5.4) 2 321 )()( k t k t k t qaqaaqC ++= for k=1..K; t=1..T. (5.5) 1+? ii bb for i=1..N (5.6) 48 )(penalty]))(()([1))(),((P of Fitness 1 1 jjqCjqPKjj K k T t k t k t k t +?= ?? = = ?qb (5.7) Where ?? ??? >?>= ? = Otherwise 0 )(or )( if )(penalty 1 maxmax N i ii QjqBjbMj (5.8) Notice that the function penalty (j) is used to penalize the objective function due to the violation of the maximum bid price and/or the maximum available capacity of the generating unit. Penalty functions are generally used to eliminate solutions that violate constraints. The following pseudo code describes the main procedure: Figure 5.2. Pseudo code for the cPSO method Main Procedure Randomly generate q? , b Set q?=*?q and b=*b While Run < NRuns Run PSO and obtain solution?q, b If fitness of P(?q,b) > fitness of P( *?q , *b ) then Set ?q?q =* and bb =* Endif EndWhile Output *?q and *b as the best solution of P 49 5.4 Decomposition Based Particle Swarm Optimization The decomposition-based PSO (dBPSO) consists of separating P into two sub- problems. One of the sub-problems assumes that the values of the decision variables bi are known for each i. This sub-problem is called PQ. Using the given values of bi?s and price scenarios, the following variable is computed as in (5.4): k ti k t PbiPI ?= such that Max )( for k=1..K; t=1..T. (5.9) Thus, the formulation of PQ is as follows: ])([1]Profit[EMax 1 1 ?? = = ?= K k T t k t k t k t qCqPK?q for k=1..K; t=1..T; (5.10) Subject to: (5.3), (5.4), (5.5) and 0?? iq for i=1..N (5.11) The second sub-problem assumes that the values of the decision variables iq? are known for each i. The objective of this optimization problem is to find the values of bi that maximize the firm?s profits. This sub-problem is called PB. ])([1]Profit[EMax 1 1 ?? = = ?= K k T t k t k t k t qCqPKb for k=1..K; t=1..T; (5.12) Subject to: (5.2), (5.3), (5.4), (5.5), and (5.6) constraints. 50 5.5 dBPSO Genetic Representation The dBPSO is used for solving both sub-problems, PQ and PB. The sub-problem PB is solved first using the algorithm dBPSOB and its solution is then used to solve PQ using the algorithm dBPSOQ. Then the new solution of PQ is used to re-solve PB. This process is applied successively until no improvement is observed after two iterations. For both problems, the population size is set to ? particles. To initialize the population of dBPSOQ, we proceed as follows: we first sample ?q1 as the amount at which the marginal production cost is equal to the given bid price b1 plus a random uniformly distributed quantity in the interval [-Q, Q], where Q is a user defined parameter. Thus, ),(2 )( 3 21 1 QQUa abq ?+?=? (5.13) and all other values ?qi (for i=2,..,N) are sampled as )1,0(2 )( 3 1 U a bbq ii i ??=? . for i=2,..,N (5.14) Similarly, we use the marginal cost function to initialize the population of dBPSOB. ),(2 1321 BBUqaab ?+?+= (5.15) )1,0(2 3 Uqab ii ?= for i=2,..,N (5.16) Since the objective functions of P, PQ, and PB are identical, we define the same fitness function for each problem. When b(j) and )(j?q denote a particular solution j (again a particle is not necessarily feasible), the fitness function is calculated as in (5.4), (5.5), (5.7) and (5.8). 51 Notice that the same penalty function penalty (j) is used to penalize the objective function due to the violation of the maximum bid price and/or the maximum available capacity of the generating unit. The following pseudo code describes the main procedure: Figure 5.3. Pseudo code for the dBPSO method In this procedure, the initial values of ?qare randomly selected from a uniform distribution as follows: ?q1 = U(0, Qmax/N) (5.17) ?qi = U(qi-1, Qmax) ? q i-1, where ? = ?= i k ki qq 1 for i=2,..,N (5.18) Main Procedure Set 0?q =* and 0b =* Randomly generate ?q NoImprovement = 0 While NoImprovement < 2 Run dBPSOB and obtain solution b Using this b as input, run dBPSOQ and obtain solution ?q If fitness of P(?q,b) > fitness of P( *?q , *b ) then Set ?q?q =* and bb =* NoImprovement = 0 Else NoImprovement = NoImprovement + 1 Endif EndWhile Output *?q and *b as the best solution of P 52 5.6 Numerical Example and Analysis The cPSO and dBPSO methods were coded in C. We run the methods using c1=2, c2=2, w=0.5 [52], and ?=30 [51]. We compute the quantity-price offers for two different generators say GEN-1 and GEN-2. The maximum capacity of GEN-1 is 400 MW and its cost function is equal to 56.52q + 0.0139q2 ($/MWh). This cost function is obtained by multiplying the heat-input function of the generator given in [29] by the current fuel cost of oil-fired units (7.2 $/MBtu) [52]. Similarly, the maximum capacity of GEN-2 is 600 MW and its cost function is equal to 43.2q + 0.108q2 ($/MWh). This cost function is obtained by multiplying the function of the generator given in [53] by 7.2 as well. In both cases, the aim is to have the marginal cost of the unit within the range of the sampled prices. Otherwise, the solution of the SBP would be trivial. In market price generation procedure, we first choose the date of May 17th, 2007, arbitrarily, and use the PJM day-ahead market prices to produce twelve 24-hour price scenarios. At each hour twelve prices are generated by sampling values from a normal distribution with mean equal to the price at that hour of May 17th and a standard deviation equal to 4 ($/MWh). For example, for hour 1 AM twelve prices are sampled from a normal distribution with mean 27.14 ($/MWh) and standard deviation 4. The price scenarios for all twenty four hours are given in Table 5.1. To show the variability among samples and the hourly fluctuations in the market price, we plot the prices in Figure 5.4. 53 Table 5.1. Day-ahead market price samples for PSO Sample Hour 1 2 3 4 5 6 7 8 9 10 11 12 1 29.36 26.54 29.47 25.30 25.09 31.33 24.26 27.05 23.79 25.17 28.28 26.59 2 24.62 24.04 24.26 24.74 26.44 31.88 28.00 26.45 21.09 24.37 24.28 24.87 3 26.09 21.12 25.66 20.57 20.29 21.28 18.79 20.80 21.50 23.20 21.96 22.31 4 21.04 16.80 19.94 21.71 23.34 21.25 22.54 21.38 19.46 21.54 18.82 20.58 5 24.16 24.30 21.24 18.89 22.24 20.51 18.20 21.63 18.25 22.61 25.04 24.29 6 24.89 23.56 25.21 26.88 25.16 25.70 25.15 27.25 25.60 23.45 27.76 23.77 7 36.87 33.02 33.47 34.77 33.24 36.83 31.67 34.86 32.51 31.01 35.13 35.82 8 43.27 46.09 44.46 44.99 44.43 44.75 46.92 43.09 47.58 46.53 47.04 48.78 9 50.33 48.31 45.03 47.86 47.85 45.83 48.64 48.90 48.45 48.01 47.12 47.92 10 50.09 51.61 54.90 54.44 52.08 55.02 52.34 51.52 54.34 50.53 50.23 50.53 11 58.65 54.42 56.90 58.72 57.59 52.62 56.27 56.17 58.57 56.40 53.20 56.75 12 58.97 57.09 57.96 57.88 56.81 57.02 58.60 56.35 54.99 63.38 56.13 59.27 13 62.51 55.20 59.87 56.92 59.85 59.02 54.83 57.97 56.97 58.22 57.82 59.38 14 60.64 63.48 61.66 55.30 60.85 58.96 58.68 60.22 58.82 64.52 58.21 61.70 15 59.63 59.17 62.39 60.45 59.99 60.60 59.86 60.68 63.71 62.35 59.85 59.16 16 56.70 58.07 59.67 58.83 61.24 61.37 56.58 61.36 60.53 58.92 63.00 61.66 17 65.77 61.57 65.12 63.23 64.56 67.54 65.22 64.65 62.40 61.77 62.87 61.79 18 58.31 59.13 59.22 57.99 64.11 56.78 56.24 57.55 59.78 55.91 57.63 58.52 19 49.60 48.38 50.14 50.82 48.86 49.97 52.34 52.97 52.01 50.92 53.74 53.00 20 47.06 48.48 49.13 48.45 47.35 49.46 49.05 48.82 48.19 48.52 48.74 48.76 21 59.66 61.85 55.34 58.47 60.59 58.58 58.35 57.67 59.93 61.16 59.85 62.18 22 61.14 55.29 59.71 56.65 61.09 56.54 58.43 58.51 63.87 60.14 54.70 62.94 23 34.77 36.46 33.50 34.29 37.06 35.50 29.80 32.79 36.34 35.31 35.58 36.44 24 19.95 28.45 24.52 26.64 23.69 26.28 29.70 30.15 26.51 26.80 29.55 27.38 Figure 5.4. Hourly day-ahead market price samples 54 5.6.1 Comparison Analysis of cPSO and Marginal Cost based Bidding We first compare the results provided by the cPSO and the marginal cost method (MC). The MC method consists of dividing the maximum capacity into ten blocks of equal size and then offering them at prices equal to their marginal cost. We plotted GEN- 1 and GEN-2 in Figure 5.5a and Figure 5.5b respectively. We compute the quantity-price offer for two different generators say GEN-1 and GEN-2. Notice that although the cost function of GEN-1 is quadratic, the curve is approximately linear in the range 0 to 400 MWh; whereas, the curve of GEN-2 is clearly quadratic. In addition, we assume that each generator bids into the PJM market (Bmax =999) where the maximum number of blocks is ten (N=10). Figure 4.5.a. Cost function of GEN-1 Figure 4.5.b. Cost function of GEN-2 We run the cPSO algorithm using 10 replications and 400 cycles, and choose the solution with the highest fitness value. The MC and cPSO solutions are given in Table 5.2. For GEN-1, the expected profit of the MC solution is $2,766. The cPSO provides a 55 solution with a slightly higher expected profit of $2,819. For GEN-2, the expected profit of the MC solution is $5,525, whereas cPSO provides a solution with a higher expected profit of $6,619. Table 5.2. Quantity-price solutions for GEN-1 and GEN-2 Block GEN-1- MC GEN-1- cPSO GEN-2 ?MC GEN-2- cPSO i bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) 1 57.64 40 57.37 68.43 56.16 60 40.10 24.55 2 58.75 80 59.20 124.17 69.12 120 52.75 64.29 3 59.87 120 60.83 177.19 82.08 180 59.55 86.30 4 60.99 160 62.02 220.51 95.04 240 98.79 179.08 5 62.11 200 63.05 258.47 108.00 300 98.89 479.77 6 63.22 240 64.41 294.94 120.96 360 99.20 480.27 7 64.34 280 65.69 312.82 133.92 420 99.30 480.77 8 65.46 320 65.79 332.54 146.88 480 99.49 481.27 9 66.58 360 66.68 396.18 159.84 540 155.02 481.77 10 67.64 400 67.64 400.00 172.80 600 172.80 600.00 The cPSO approach is more effective for GEN-2 than for GEN-1. The small difference in profit for GEN-1 seems to be due to the small value of the quadratic coefficient of its cost function. To show that, we change the quadratic coefficient a3 of GEN-1 between 0.0139 and 0.1139 in steps of 0.01 and compute the quantity-price solutions for each of these values using the MC and cPSO methods. We plot the difference of profits in relation to the MC solution in Figure 5.6 below. 56 0 2 4 6 8 10 12 14 16 18 0 0.02 0.04 0.06 0.08 0.1 0.12 Quadratic Cost Coefficient (a3) Pr of it di ffe re nc e ( % ) Figure 5.6. Percentage of increase of cPSO solution in relation to MC 5.6.2. Comparison Analysis of cPSO and dBPSO For the comparison analysis of cPSO and dBPSO, we use GEN-2 and the same price samples and PSO parameters. To compare the performance of both methods in terms of the fitness value, we use three different combinations of replications and cycles (see Table 5.3). Notice that the dBPSO algorithm evaluates two fitness functions at each iteration, which we refer to as a ?cycle?, while cPSO evaluates one function at each cycle. We compare both methods using different numbers of cycles and replications. A replication is a complete run of the algorithm using a different starting seed for generating the random numbers. We select three combinations of numbers of cycles and replications to reach a total of 4,000 function evaluations in each PSO. These values are given in Table 5.3. To assure that these combinations are indeed different settings, we conduct an ANOVA test for cPSO and dBPSO. This procedure allows us to properly compare both approaches under similar parameters and diverse setting. 57 Table 5.3: Parameter set for the experiment Parameter cPSO dBPSO Number of Replication (R) 10 20 40 10 20 40 Number of Cycles (C ) 400 200 100 200 100 50 Total Function Evaluations 4000 4000 4000 4000 4000 4000 First, we run the cPSO for each combination. In Table 5.4 I give the averages, standard deviations, and the best fitness values. The bid prices and quantities that provide the best fitness are given in Table 5.5. Second, we conduct the one-way ANOVA test on means to determine whether there is a statistical difference among these results. The values of the ANOVA table are given in Table 5.6. The results of the test show that with a p-value=0.15 and level of significance ?=0.05 there is no statistical evidence to show that the setting conditions are different. This result indicates that cPSO is robust with regard to changes to R and C. Table 5.4. Statistical summary of the results of the cPSO solutions Parameter set Average fitness($) Standard deviation ($) Best fitness ($) R=10, C=400 6,346.84 143.08 6,530.02 R=20, C=200 6,304.85 125.24 6,529.35 R=40, C=100 6,315.40 127.89 6,524.05 Table 5.5. Bid prices and quantities for best solutions of cPSO R=10, C=400 R=20, C=200 R=40, C=100 Block i bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) 1 46.72 31.57 46.89 31.60 45.86 27.66 2 54.75 75.58 54.80 75.58 54.70 75.69 3 114.76 109.02 107.84 348.03 77.03 290.36 4 136.21 298.12 131.37 430.19 77.13 386.62 5 136.31 432.95 166.70 448.01 80.01 417.34 6 136.41 477.58 166.80 448.51 80.11 481.06 7 136.51 478.08 166.90 449.01 80.39 503.23 8 136.61 480.47 167.00 449.51 80.49 503.73 9 136.71 480.97 167.10 450.01 80.59 504.23 10 172.80 600.00 172.80 600.00 172.80 600.00 58 Table 5.6. ANOVA test results for cPSO Source of Variation Sum of Square DF Mean Square F0 P-Value Factor 11959 2 5980 0.36 0.701 Error 1120190 67 16719 Total 1132149 69 Similarly, we run dBPSO for each of the setting conditions. In Table 5.7, I show the averages, standard deviations, and the best fitness values. The bid prices and quantities that provide the best fitness are given in Table 5.8. The same ANOVA test on means was applied to the experimental data of dBPSO. The results of the statistical test are given in Table 5.9. This ANOVA table also shows no evidence of statistical difference among the setting conditions with a level of significance ?=0.05 and p- value=0.488. dBPSO is also robust with regard to changes to its parameters. Table 5.7. Statistical summary of dBPSO solutions Parameter Set Average fitness ($) Standard Deviation ($) Best fitness ($) R=10, C=200 6,607.18 64.91 6,691.43 R=20, C=100 6,630.04 52.10 6,701.11 R=40, C= 50 6,612.69 61.29 6,689.76 Table 5.8: Bid prices and quantities for best solutions of dBPSO R=10, C=200 R=20, C=100 R=40, C=50 Block i bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) 1 46.45 27.71 46.17 27.72 46.20 27.02 2 52.48 57.18 52.48 57.18 52.29 56.36 3 57.54 74.03 57.18 71.46 57.46 73.54 4 57.64 74.53 59.84 81.74 61.08 89.17 5 61.57 91.34 59.94 82.24 61.18 89.67 6 61.67 91.84 60.04 82.74 61.28 90.17 7 61.77 92.34 62.85 96.58 61.38 90.67 8 61.87 92.84 62.95 97.08 61.48 91.17 9 61.97 93.34 63.05 97.58 61.58 91.67 10 172.80 600.00 172.80 600.00 172.80 600.00 59 Table 5.9. ANOVA test results for dBPSO Source of Variation Sum of Square Df Mean Square F0 P-Value Factor 5,106 2 2,553 0.72 0.488 Error 236,012 67 3,523 Total 241,118 69 Using an analysis of means and a t-test, we test whether the mean of the fitness values of dBPSO is greater than that of cPSO. The results are given in Table 5.10. It shows strong evidence with p=0.00 and confidence level ?=0.05 that dBPSO provides better solutions than cPSO on the average. Table 5.10. ANOVA test results for dBPSO for method comparisons Sample Sample Size Mean St. Deviation T-Value P-Value dBPSO 70 6616.9 59.10 16.57 0.00 cPSO 70 6316.9 128.10 Difference 70 300.00 151.50 5.6.3. Converging Process of dBPSOB and dBPSOQ To illustrate the converging process of the two sub problems, dBPSOB and dBPSOQ, of the dBPSO approach, we plot the evolution of the fitness function of each sub problem with respect to the number of iterations. In Figures 5.7 through 5.9, we show the evolution of the fitness value for the number of cycles equal to 200, 100, and 50, respectively. Notice that in the three plots the fitness value of the dBPSOB starts at a low value, whereas the dBPSOQ starts at a higher value. This is because we begin by solving dBPSOB using randomly selected bi values, while dBPSOQ uses ?qi values that are determined by dBPSOB. Also notice that for the three cases the two sub-problems converge before the 25th iteration. 60 -1000 1000 3000 5000 7000 0 5 10 15 20 25 Cycle (C) Pr ofi t($ ) dBPSOQ dBPSOB Figure 5.7: Converging process of dBPSOB and dBPSOQ to the solution (C=200) -1000 1000 3000 5000 7000 0 5 10 15 20 25 Cycle (C) Pr ofi t($ ) dBPSOQ dBPSOB Figure 5.8: Converging process of dBPSOB and dBPSOQ to the solution (C=100) 61 -1000 1000 3000 5000 7000 0 5 10 15 20 25 Cycle (C) Pr ofi t($ ) dBPSOQ dBPSOB Figure 5.9: Converging process of dBPSOB and dBPSOQ to the solution (C=50) 5.6.4 Impact of the order of dBPSOB and dBPSOQ In order to analyze the impact of the order of solving dBPSOB and dBPSOQ, we solve the problem using dBPSOB first and then dBPSOQ. In this procedure, the initial values of b are randomly selected from a uniform distribution as follows: b1 = U(0, Bmax/N) (5.19) bi = U(bi-1, Bmax) ? b i-1, for i=2,..,N (5.20) In Figures 13 through 15, we show the evolution of the fitness value for number of cycles equal to 200, 100, and 50, respectively. 62 -1000 1000 3000 5000 7000 0 5 10 15 20 25 Cycle (C) Pr of it ($ ) dBPSOQ dBPSOB Figure 5.10: Converging process of dBPSOB and dBPSOQ to the solution (C=200) -1000 1000 3000 5000 7000 0 5 10 15 20 25 Cycle (C) Pr of it ($ ) dBPSOQ dBPSOB Figure 5.11: Converging process of dBPSOB and dBPSOQ to the solution (C=100) 63 -1000 1000 3000 5000 7000 0 5 10 15 20 25 Cycle (C) Pr of it ($ ) dBPSOQ dBPSOB Figure 5.12: Converging process of dBPSOB and dBPSOQ to the solution (C=50) Similarly to the former case of dBPSO, in the three plots the fitness value of the dBPSOQ starts at a low value, whereas the dBPSOB starts at a higher value. This is because we begin by solving dBPSOQ using randomly selected ?qi values, while dBPSOB uses bi values that are determined by dBPSOQ. Notice that for the three cases the two sub-problems converge again before the 25th iteration. Table 5.11 gives a summary of the results. After comparing the plots and Table 5.11 with Table 5.7, the order of the sub- problems is not relevant for solving the problem. Table 5.11: Statistical summary of dBPSO solutions Parameter Set Average fitness ($) Standard Deviation ($) Best fitness ($) R=10, C=200 6,615.90 73.83 6,688.56 R=20, C=100 6,635.57 55.67 6,688.72 R=40, C= 50 6,614.85 69.28 6,689.48 5.7 Conclusion Results showed that PSO outperforms the MC method on determining price- quantity pairs that will be submitted to the day-ahead market. The percentage of 64 improvement is accentuated when the quadratic coefficient of the generator?s cost function is significant. In terms of time and experimentation burden, both models took three minutes on average to find a solution. All three experiments of dBPSO showed that we need less than 25 iterations to obtain a good solution. The results also showed that dBPSO give much better solutions than cPSO. This tells us that the decomposition technique can be applied to problems that have two or more decision variable sets. The problems can be decomposed into two or more parts, and one or more decision sets can be used as components of the solution of the other parts of the problem. This process continues until no improvement is observed. The model discussed here can be further improved by including a forecasting technique for the market prices. If the characteristics of the generating unit and constraints change, the model and solution approach still may remain valid. 65 CHAPTER 6 AGENT BASED PARTICLE SWARM OPTIMIZATION FOR SUPPLY FUNCTION EQUILIBRIUM 6.1 Introduction In the previous chapters, we developed models for an individual generator which submits a bidding curve to the PJM Day-ahead market. We developed two PSO methods to find a heuristic solution to the problem. The models did not include the behavior of all market participants. In this chapter, we include the strategic change in competitors? behavior for a particular generator. The model assumes that the strategy employed by one player is affected by the others? behavior. Game theory and agent based models are two ways to represent this market interaction. We develop an agent based simulation method to simulate the behaviors of all firms in the market. We combine the dBPSO approach and agent based model to compute an equilibrium solution. 6.2 Supply Function Equilibrium Model This model assumes that there are m participants in a power pool who may be referred to as power suppliers. These power suppliers may either have individual generators or even a portfolio of generators. All these participants bid into the day-ahead market and aims to maximize their profit by using bidding strategies that represent their expectations best. The ISO collects the buy bids simultaneously and it starts the security 66 constrained optimum dispatch algorithm to set the equilibrium for the market. The ISO sorts the sell offers starting from minimum price offers to more expensive ones and sorts the buy offers starting from the maximum price offers to less expensive ones. The ISO sets the equilibrium market price where the aggregated supply and demand meet. Figure 6.1 shows the equilibrium process in the day-ahead market [54]. In this figure, supply represents the power quantities offered and demand represents the respective buy offers; bij and qij represent the offered bid price and the power quantity respectively. Notice that the day-ahead equilibrium price is used for uniform bid auctions, and the winners will be paid the MCP. bij Equilibrium Pt Market Price ($/Mwh) Power Allocation (Bid Quantities Mwh) Demand Supply qij (Bid Prices) The model is formulated under the following additional assumptions: Figure 6.1. Equilibrium process in day-ahead power market 67 i. An offer or bid, which consists of N price-quantity blocks at the most, needs to be determined for each firm separately. ii. Each firm can build its strategy based on separate resources or on its portfolio of power resources. iii. The equilibrium for each day is determined before the market closes at noon and is valid for the next twenty four hours, starting at midnight the same day. iv. Demand can be forecasted and is known for the analysis v. The transmission constraints are not included in the model. vi. Equilibrium of interest occurs in a single-round auction market. In finding an equilibrium, there are Nxm pairs of decision variables bji and ?qji (j=1,?M) (i=1,..,N) for each firm j that need to be determined.. The variable ?qji denotes the amount of energy increase in firm j in block i, to get the bid price bji for delivery at any hour of the next day. Total energy to be produced at time t and sold to the market at a price Pt by the firm j is given by: ? = ?= )( 1 )( tPI i jitjt qPq , tjit PbiPI ?= such that Max )( for j=1...m; t=1...T; i=1..I(Pt) (6.1) In equilibrium, it is accepted that firms cannot make more profit by bidding other than their current bid. Also the cost of dispatching is minimized for the system operator at this state [33]. 68 The maximum profit for the firm j in equilibrium can be expressed as: ? = ?=?= T t tjttjttjj PqCPqP jj 1 j, ))(()( Max),P( b?qb?q for j=1...m; t=1..T; (6.2) Subject to the following constraints [55]: ))(),(())(),(( *** tjttjtjtjttjtj PqPqPqPq ?? ??? for j=1...m; t=1...T; (6.3) max 1 j N i ji Qq ??? = for i=1..N; j=1...m; (6.4) max0 Bb ji ?? for j=1...m; i=1...N; (6.5) max0 jji Qq ??? for j=1...m; i=1...N; (6.6) ? = ?= )( 1 )( tPI i jitjt qPq , (6.7) Where tjit PbiPI ?= such that Max )( for j=1..m; t=1..T; i=1...N; 2 321 ))(()())(( tjtjtjtjjtjt PqaPqaaPqC ++= for j=1..m; t=1..T; (6.8) The equilibrium in an economic system requires supply and demand to be equal. As equilibrium constraint, the total amount of power generated is equal to total supply as given below: t m j tjt DPq =? =1 )( t=1..T; (6.9) 69 6.3 Agent Based Modeling and Simulation Agent-Based Modeling and Simulation (ABMS) is a computational approach to model economic systems which have interacting components or dynamic agents. Agents usually interact among themselves and between environments by updating themselves sequentially rather than simultaneously. In building an ABMS, the definition of agents and their interaction environment are crucial. Figure 6.2 shows the ABMS building process [56]. Figure 6.2. Agent based modeling and simulation 70 ABMS is applied to very complex systems where interdependencies are difficult to capture and traditional models are hard to apply. ABMS replaces its framework with an individual agent?s behavioral rules that are updated over time [57]. ABMS is a descriptive method which aims to model the behavior of agents rather than optimality. ABMS models are useful in economics models. In a micro-economic point of view, agents assume that 1) they behave in a rational manner that aims to optimize their well-defined objectives 2) they have identical characteristics that make them alike 3) they will have decreasing marginal utility as the number of agents increase 4) long-run equilibrium of the system is of primary interest to the model [58], [59]. Figure 6.3 shows the process of a typical ABMS [57], [58], [59]. Figure 6.3. General flow of a typical ABMS 71 ABMS, like many other heuristic methods, searches for the best feasible solution by an updating process. The process is evolutionary, in which updating, learning and convergence are involved. ABMS has been used in flow management of evacuation, traffic, stock market, strategic simulation of market, organizational design, and in other areas where players dynamically move [56], [58], [60]. In [61], an agent-based model is developed for a supply chain in which the flow of a commodity finds a way between factories, distributors, wholesalers, retailers and customers. The goal of supply chain agents is to minimize their cost on their way. These are the research assumptions that are made for the agent based modeling to work. In [57], authors develop an agent based model for modeling the human immune system. Authors in [62] address the agent-based modeling approach in financial markets by explaining trading behaviors of firms. In [54], [59] authors develop an agent-based simulation approach for modeling the day-ahead power trading in the US wholesale power market. The model is developed for wholesale power suppliers, individual power generators and wholesale power consumers that are bidding into the day-ahead market. Their hypothesis is that agent- based approach is as good as other approaches like neural networks. They show the effectiveness of the model using data provided in PJM west. Agents bid into the market and they update their bidding strategies in each run of the simulation based on a learning factor until the equilibrium is reached. In [63], an agent-based model was developed for the wholesale electricity market, operating in a short-term environment under capacity conditions and double auction rules. A simulation model that uses a multi-adaptive agent model for generators bidding in the UK power market was developed in [64]. It shows that agents learn bidding 72 strategies in a manner similar to their behavior in real world. In [65], authors use an agent based simulation model to show how companies use market power during the electricity market bidding process. They compare the production cost bidding with the bidding strategies based on physical and economic withholding, including congestion management. 6.4 Agent based particle swarm optimization model It was demonstrated in Chapter 5 that the dBPSO method is an effective way to find a solution for an individual generator. Now let?s assume that m firms compete in the day-ahead market using dBPSO. The objective of this chapter is to show that the model can reach an equilibrium point at which each player?s payoff is maximized and the equilibrium conditions are met using ABMS and dBPSO. Each firm participating in the day-ahead market is modeled as an agent. Figure 6.4 shows the details for the ABMS method with dBPSO applied. Each agent has a unique cost function, a capacity and pairs of quantity-price bids as attributes. The agents? interaction occurs in the pool where offered quantities and corresponding prices are submitted. Agents aim to allocate their price-quantity offers in a way that their profit is maximized. In other words, their interaction occurs based on power quantities and offered prices. 73 Figure 6.4: Components of ABMS method The model assumes that the day-ahead market demand is known and the final objective is to reach equilibrium for the next day. The model starts with an initial price scenario and given demand for 24 hours. It runs until the equilibrium is reached. Each agent?s objective is to maximize its profit by bidding into the environment using dBPSO. Figure 6.5 shows a flow chart of the Agent Based Particle Swarm Optimization (ABPSO). 74 Figure 6.5: The flow of the ABPSO 75 Notice that each agent in the model is a price-taker while bidding into market. The cleared prices that are set at each iteration are the result of bidding strategies submitted by each agent. Since the analysis is for a single-round auction market, an agent can use the simulation to test the behavior of its rivals. Agents update their price and quantity bids at each iteration until the equilibrium is reached. The particular agents simulating the model would be able to observe the behavior of their competitors during this process. This iterative process ensures that each generator minimizes its risk of cost recovery and infeasibility of its offer. The definition of stopping criteria is important in ABSM since it shows whether the equilibrium is reached or not. In [21], authors propose two kinds of stopping criteria. First criterion stops the equilibrium process when the prices are too low. This is because the demand is not sufficiently covered at this price since some units will be eliminated due to low prices. At this point, iteration goes 2 steps back and defines the point as equilibrium point. The second criterion calculates demand-weighted average price and total financial loss for all players and it compares the values with those of the previous ones. When the current average prices are higher than those of the previous two iterations the process stops. The submission of the same or similar strategies at each iteration could lead to similar market prices. It might indicate that market price is converging and strategies submitted are reaching equilibrium. However, we need to verify that this state actually satisfies equilibrium conditions. We evaluate three new stopping rules to end the agent- based process in order to find the best stopping rule that satisfies equilibrium conditions most of the time. In the first stopping rule, we calculate the percentage of differences of 76 resulting market prices and previous iterations? market prices for each hour. When the absolute value percentage of market price differences for the 24 hours are less than a value ?1 the process stops. It can be represented mathematically as: 1 1|| ???=? ? r t r t r tr t P PP t=1..T. (6.10) In the second stopping rule, we calculate the average of differences of resulting market prices and previous iterations? market prices for each hour. At each iteration, current market prices are compared with previous market prices. When the absolute average value of market price differences for 24 hours are less than a small value ?2 the process stops. It can also be represented mathematically as: 2 1 1 ||1 ???=? ? = ? T t r t r t r PP T (6.11) In the third stopping rule, we calculate the weight of each hour?s demand with regard to the total demand for 24 hours. At each iteration, we multiply this percentage for each hour with the price difference percentage found in the first stopping condition. When the average value of this calculation is less than a small value ?3 the process stops. The mathematical representation is: 77 3 1 1 1 ||1 ???=? ? ? = ? = T t t r t t r t r t T t r DP DPP Tw (6.12) Notice that at each iteration, agents could have available strategies from which they choose the best strategy for themselves. This strategy is selected based on interactions with other agents and interaction with the environment. The interaction with the environment occurs in a way that total supply of all agents should be equal to total demand. The amount of power allocated to an agent affects other agents. The interaction between agents occurs based on the offer prices and offer quantities. If the offer of an agent is selected, the agent can either maintain this strategy or update it in the next iteration in order to get better results. This process continues until the price difference in the two iterations is so low that it might lead to convergence. In order to test whether the results of the experiments actually satisfy the equilibrium conditions, we define conditions for equilibrium. I. Supply should be equal to demand as an essential rule of equilibrium. II. In equilibrium, all suppliers maximize their return at the given prices. Also according to the Nash equilibrium, if a player deviates from its strategy it will loose in the long run [38], [39], [40]. We will let each firm behave separately and run the dBPSO to get a separate bidding strategy and profit. The results found in dBPSO will be compared with the ABPSO results for each firm. It is expected that firms find a close strategy in dBPSO to the strategy found in equilibrium. If the difference is less than 1%, we will accept this as a satisfied condition. 78 6.5 Numerical Example and Analysis We coded the model in C and used the same dBPSO parameters we used in Chapter 5. We are more interested in the supply side bidding rather than demand side bidding. The demand for next day is known and is given in Table 6.1. Table 6.1. Day-Ahead demand for next day Hour Demand(Mwh) Hour Demand(Mwh) 1 3115 13 4510 2 3711 14 5142 3 3346 15 3424 4 3771 16 3287 5 3298 17 4501 6 4266 18 5236 7 4117 19 5790 8 5176 20 6084 9 5751 21 6561 10 6513 22 6411 11 6280 23 4411 12 4472 24 4664 6.5.1 ABPSO Experiment with Duopoly We start with a duopoly SMD that has m=2 firms competing. The cost functions of the firms and their capacities are given in Table 6.2. Table 6.2. Cost functions and market capacities of the firms for m=2 Firm(m) a1 a2 a3 Unit 1 0 41.73 0.0063 8085 2 0 47.25 0.0057 6281 79 We start with the first stopping condition set to the value ?1 = 1%, in other words, the percentage price difference for each of the 24 hours is less than 1%. The equilibrium prices found are given in Table 6.3. Table 6.3. Equilibrium prices for first stopping rule (m=2) Hour Price($/Mwh) Hour Price($/Mwh) 1 77.15 13 77.15 2 77.15 14 77.15 3 77.15 15 77.15 4 77.15 16 77.15 5 77.15 17 77.15 6 77.15 18 77.15 7 77.15 19 80.29 8 77.15 20 80.29 9 80.29 21 83.98 10 83.68 22 83.68 11 83.68 23 77.15 12 77.15 24 77.15 Table 6.4 shows the bids found for each firm at the equilibrium and their profits. Table 6.4. Results found in first stopping rule (m=2) Firm-1 Firm-2 Block i bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) 1 70.79 2856.23 77.15 2642.64 2 70.89 2861.75 77.42 2887.47 3 70.99 2862.25 83.68 3200.27 4 80.29 3254.57 83.78 3200.77 5 80.47 3338.97 83.88 3218.77 6 84.57 3450.30 83.98 3247.07 7 86.50 3503.35 84.08 3254.22 8 86.60 3507.39 84.18 3254.72 9 86.70 3507.89 84.28 3255.22 10 143.60 8085.00 118.85 6281.00 80 ?1 = 1% was chosen because values smaller than 1% increase the computational time and in some cases the stopping condition could not be reached in a reasonable time. In the second stopping rule we set ?2 = $0.25. In other words, the average price difference for the 24 hours should be less than $0.25. The values less than $0.25 also increase computational time and many times return no results. The equilibrium prices found are given in Table 6.5. Table 6.5. Equilibrium prices for second stopping rule (m=2) Hour Price($/Mwh) Hour Price($/Mwh) 1 77.15 13 77.15 2 77.15 14 77.15 3 77.15 15 77.15 4 77.15 16 77.15 5 77.15 17 77.15 6 77.15 18 77.15 7 77.15 19 80.29 8 77.15 20 80.29 9 80.29 21 83.98 10 83.68 22 83.68 11 83.68 23 77.15 12 77.15 24 77.15 Table 6.6 shows the bids submitted by each firm at the equilibrium and firms? profits. 81 Table 6.6. Results found in second stopping rule (m=2) Firm-1 Firm-2 Block i bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) 1 70.79 2856.23 77.15 2642.64 2 70.89 2861.75 77.42 2887.47 3 70.99 2862.25 83.68 3200.27 4 80.29 3254.57 83.78 3200.77 5 80.47 3338.97 83.88 3218.77 6 84.57 3450.30 83.98 3247.07 7 86.50 3503.35 84.08 3254.22 8 86.60 3507.39 84.18 3254.72 9 86.70 3507.89 84.28 3255.22 10 143.60 8085.00 118.85 6281.00 In the third stopping rule we set ?3 =1%. In other words, the load weighted average for the 24 hours should be less than 1%.The equilibrium prices found are given in Table 6.7. Table 6.7. Equilibrium prices for the third stopping rule (m=2) Hour Price($/Mwh) Hour Price($/Mwh) 1 78.45 13 78.45 2 78.45 14 78.45 3 78.45 15 78.45 4 78.45 16 78.45 5 78.45 17 78.45 6 78.45 18 78.45 7 78.45 19 79.77 8 78.45 20 79.77 9 79.77 21 84.17 10 84.17 22 82.6 11 82.6 23 78.45 12 78.45 24 78.45 Table 6.8 gives the bids submitted by each firm at the equilibrium and profits. 82 Table 6.8: Results found in third stopping condition (m=2) Firm-1 Firm-2 Block i bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) 1 58.01 1343.58 78.45 2744.24 2 78.32 2955.50 78.59 2795.03 3 79.77 3318.26 82.60 3168.58 4 93.90 3832.65 84.17 3247.16 5 94.00 3840.36 84.27 3275.26 6 94.38 3840.86 84.37 3275.76 7 94.48 3841.36 84.47 3276.26 8 94.58 3841.86 84.57 3276.76 9 94.68 3863.55 84.67 3277.26 10 143.60 8085.00 118.85 6281.00 Using the equilibrium prices found, we run dBPSO to find a solution for each firm separately. The summary of the results are given in Table 6.9. Table 6.9: Overview of the results found in each method (m=2) Method Firm-1 Profit ($) Firm-2 Profit ($) First stopping rule 1,302,689 1,043,974 Second stopping rule 1,302,689 1,043,974 Third stopping rule 1,355,414 1,094,031 All three stopping rules satisfy the balance condition, which is supply equals demand. The second and third stopping rules generally require close computational times which is around 15 minutes. Notice that rule-1 and rule-2 return the same values, i.e., they both satisfy each rule. We let each firm develop a separate strategy. To do so, we use the equilibrium prices found in each rule and run dBPSO. Table 6.10 shows the percentage of profit increase in dBPSO comparing with the equilibrium solution. 83 Results show that three rules give similar profit increases which are smaller than acceptable level 1%. However, the profits found in rule-1 are almost same with the ABPSO solution. Table 6.10: Profit increases in each rule (m=2) Rule 1 Rule 2 Rule 3 Firm 1 0.08% 0.08% 0.15% Firm 2 0.00% 0.00% 0.01% 6.5.2 ABPSO experiment with m=5 We now suppose that there is a SMD that has m=5 firms computing. The cost functions of the firms and their capacities are given in Table 6.11. Table 6.11. Cost functions and market capacities of the firms for (m=5) Firm(m) a1 a2 a3 Unit 1 0 47.273 0.0074 3450 2 0 45.18 0.0048 1600 3 0 44.76 0.0066 2935 4 0 45.35 0.0087 4950 5 0 46.72 0.0061 2281 We set the ?1 = 1% and the find the results with firs stopping condition. The equilibrium prices found are given in Table 6.12 and the results for each firm are given in Table 6.13 and Table 6.14. 84 Table 6.12. Equilibrium prices for first stopping rule (m=5) Hour Price($/Mwh) Hour Price($/Mwh) 1 56.24 13 57.83 2 56.30 14 59.18 3 56.24 15 56.24 4 56.30 16 56.24 5 56.24 17 57.83 6 57.17 18 60.46 7 56.30 19 63.11 8 60.46 20 63.21 9 63.11 21 63.90 10 63.51 22 63.24 11 63.24 23 57.83 12 57.83 24 58.36 Table 6.13. Results found for firm 1, 2 and 3 in first stopping rule (m=5) Firm-1 Firm-2 Firm-3 Block i bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) 1 54.76 636.21 52.77 1090.48 53.57 922.97 2 57.83 829.31 54.67 1090.98 60.46 1365.98 3 63.11 1061.94 55.88 1199.05 64.78 1610.22 4 63.21 1082.94 58.36 1397.38 64.88 1617.56 5 63.31 1083.44 58.86 1500.82 64.98 1618.06 6 63.41 1088.85 58.96 1501.32 65.08 1627.13 7 63.51 1101.42 59.06 1576.97 65.18 1627.63 8 63.61 1105.74 59.16 1577.47 65.28 1648.61 9 63.71 1128.88 59.26 1595.29 65.38 1649.11 10 98.33 3450.00 60.54 1600.00 83.50 2935.00 85 Table 6.14. Results found for firm 4 and 5 in first stopping rule (m=5) Firm-4 Firm-5 Block i bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) 1 56.24 690.41 56.30 788.31 2 63.24 1031.96 56.40 814.21 3 63.34 1061.27 57.17 908.10 4 64.42 1112.65 59.18 1128.05 5 64.68 1181.62 63.16 1359.23 6 65.25 1182.12 63.33 1383.06 7 65.35 1182.62 63.43 1386.21 8 65.45 1185.39 63.53 1386.75 9 65.55 1185.89 63.90 1414.27 10 131.48 4950.00 74.55 2281.00 In the second stopping rule we keep the same value and set the ?2 = $0.25. The equilibrium prices found are given in Table 6.15. Table 6.16 and Table 6.17 below show the bids submitted by each firm at the equilibrium. Table 6.15. Equilibrium prices for second stopping rule (m=5) Hour Price($/Mwh) Hour Price($/Mwh) 1 54.61 13 57.69 2 54.95 14 59.65 3 54.61 15 54.61 4 56.16 16 54.61 5 54.61 17 57.69 6 57.29 18 60.54 7 57.29 19 61.86 8 59.65 20 62.18 9 61.86 21 63.95 10 63.75 22 63.75 11 62.67 23 57.69 12 57.69 24 57.89 86 Table 6.16. Results found for firm 1, 2 and 3 in second stopping rule (m=5) Table 6.17. Results found for firm 4 and 5 in second stopping rule (m=5) Firm-4 Firm-5 Block i bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) 1 46.40 532.01 54.61 663.56 2 50.27 604.27 54.71 685.44 3 56.16 764.99 57.29 928.55 4 61.86 967.41 58.93 1078.69 5 63.72 967.91 61.41 1287.72 6 63.82 994.97 63.75 1396.49 7 63.93 1009.70 63.85 1396.99 8 64.03 1057.48 63.95 1400.27 9 64.13 1057.98 64.10 1414.97 10 131.48 4950.00 74.55 2281.00 For the third stopping rule we also go with the same stopping value and set ?3 = 1%. The equilibrium prices found are given in Table 6.18. Firm-1 Firm-2 Firm-3 Block i bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) 1 53.21 527.86 53.12 835.03 50.71 781.89 2 56.61 769.26 53.28 877.46 57.69 1052.93 3 61.25 1001.89 54.07 982.34 62.18 1338.63 4 61.35 1006.83 54.95 1113.52 62.67 1433.25 5 61.87 1007.33 57.89 1330.21 62.77 1439.68 6 61.97 1007.83 57.99 1332.43 67.37 1463.38 7 62.60 1040.61 58.13 1390.31 67.47 1474.14 8 62.70 1064.30 59.65 1511.66 67.57 1475.78 9 63.17 1114.08 59.75 1569.86 67.67 1479.80 10 98.33 3450.00 60.54 1600.00 83.50 2935.00 87 Table 6.18. Equilibrium prices for third stopping rule (m=5) Hour Price($/Mwh) Hour Price($/Mwh) 1 54.07 13 58.75 2 54.07 14 60.85 3 54.07 15 54.07 4 54.10 16 54.07 5 54.07 17 58.75 6 57.36 18 60.85 7 57.13 19 62.23 8 60.85 20 62.43 9 62.23 21 63.88 10 62.76 22 62.43 11 62.43 23 57.36 12 58.75 24 60.54 Table 6.19 and Table 6.20 below show the bids submitted by each firm at the equilibrium with third stopping rule. The summary of the results are given in Table 6.21. Table 6.19. Results found for firm 1, 2 and 3 in third stopping condition (m=5) Firm-1 Firm-2 Firm-3 Block i bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) 1 54.12 623.77 55.47 1075.35 55.24 814.10 2 59.82 934.01 55.67 1224.07 55.88 1010.89 3 63.23 1084.28 58.32 1377.30 60.30 1237.79 4 63.33 1088.76 58.42 1377.80 63.27 1404.40 5 63.43 1095.11 58.52 1378.30 63.37 1404.90 6 63.56 1095.61 58.62 1378.80 63.47 1412.75 7 63.66 1115.29 58.72 1379.30 63.59 1418.85 8 63.76 1115.79 58.82 1460.73 63.69 1432.33 9 63.86 1116.29 58.92 1461.23 63.79 1440.80 10 98.33 3450.00 60.54 1600.00 83.50 2935.00 88 Table 6.20. Results found for firm 4 and 5 in third stopping condition (m=5) Firm-4 Firm-5 Block i bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) 1 54.38 641.24 46.25 0.50 2 59.40 894.78 54.62 802.07 3 62.08 965.90 60.51 1134.85 4 62.33 998.53 60.61 1159.17 5 62.67 999.03 60.71 1170.19 6 62.77 999.53 61.58 1263.43 7 63.30 1034.49 63.29 1367.16 8 63.40 1042.20 63.39 1367.66 9 63.50 1059.84 63.49 1398.41 10 131.48 4950.00 74.55 2281.00 Table 6.21. Overview of the results found in each method (m=5) Firm-1 Profit($) Firm-2 Profit ($) Firm-3 Profit($) Firm-4 Profit($) Firm-5 Profit($) First rule 121,337 252,427 195,799 133,962 141,900 Second rule 112,390 235,851 182,632 127,743 149,483 Third rule 110,166 240,275 182,493 126,768 149,629 In terms of computational time, again first stopping rule is the most time consuming. The second and third stopping rules took around 35 minutes. Again we let each firm develop a separate strategy. Table 6.22 shows the percentage of profit increase in dBPSO comparing with the equilibrium solution. 89 Table 6.22. Profit increases in each rule (m=5) 6.5.3 ABPSO Experiment with m=10 We now suppose that there is a SMD that has m=10 firms competing. Based on the results found for m=2 and m=5 it is better to use rule-2 as stopping rule. The cost functions of the firms and their capacities are given in Table 6.23 below. Table 6.23. Cost functions and market capacities of the firms (m=10) Rule 1 Rule 2 Rule 3 Firm 1 0.46% 0.14% 4.28% Firm 2 0.17% 0.40% 0.11% Firm 3 0.51% 0.85% 2.20% Firm 4 2.87% 0.72% 2.30% Firm 5 13.51% 0.06% 1.54% Firm(m) a1 a2 a3 Capacity 1 0 46.18000 0.00477 500 2 0 32.95070 0.002357 600 3 0 42.40000 0.004664 250 4 0 40.12000 0.004364 1100 5 0 41.75679 0.003896 585 6 0 46.26748 0.007919 3000 7 0 42.71000 0.016201 1528 8 0 44.68000 0.017737 2000 9 0 42.45727 0.006044 403 10 0 43.19774 0.006153 4400 90 For given ?2 = $0.45, the equilibrium prices found are given in Table 6.24. Table 6.25 through Table 6.28 below show the bids submitted by each firm at the equilibrium. The summary of the results and profit increase percentages are given in Table 6.29 and table 6.30 respectively. Table 6.24:.Equilibrium prices for second stopping rule (m=10) Hour Price($/MWh) Hour Price($/MWh) 1 47.18 13 50.28 2 48.56 14 53.26 3 47.18 15 47.18 4 48.56 16 47.18 5 47.18 17 50.28 6 49.04 18 53.26 7 49.03 19 57.32 8 53.26 20 57.33 9 57.32 21 59.15 10 59.08 22 58.94 11 57.42 23 49.92 12 50.18 24 51.57 Table 6.25. Results found for firm 1, 2 and 3 (m=10 Firm-1 Firm-2 Firm-3 Block i bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) 1 46.90 117.39 20.91 0.50 14.03 1.96 2 48.56 283.36 22.02 85.45 14.13 3.18 3 49.92 394.32 22.27 99.43 14.23 3.68 4 50.18 425.57 22.37 114.02 14.33 12.13 5 50.28 455.75 22.55 122.16 14.43 23.87 6 50.38 456.25 22.65 124.67 14.53 33.95 7 50.48 478.10 22.75 145.18 14.63 41.91 8 50.58 478.78 22.85 149.52 14.73 51.55 9 50.68 479.28 22.95 152.46 14.83 56.99 10 50.95 500.00 35.78 600.00 44.73 250.00 91 Table 6.26. Results found for Firm 4, 5 and 6 (m=10) Table 6.27. Results found for firm 7 and 8 (m=10) Firm-7 Firm-8 Block i bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) 1 46.71 187.75 44.57 116.03 2 52.99 347.34 51.57 261.76 3 56.30 451.20 54.82 390.04 4 58.18 479.35 59.08 452.94 5 58.28 479.85 59.19 453.44 6 58.38 480.35 59.29 453.94 7 58.48 480.85 59.39 456.96 8 58.58 502.14 59.49 457.46 9 58.68 502.64 59.59 457.96 10 92.22 1528.00 115.63 2000.00 Firm-4 Firm-5 Firm-6 Block i bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) 1 47.18 809.78 24.65 0.50 49.03 363.19 2 47.28 825.15 27.70 343.13 57.33 784.18 3 48.32 973.01 27.95 368.73 58.94 866.18 4 49.04 1028.94 28.07 383.73 59.04 878.95 5 49.14 1040.43 28.17 385.00 59.14 889.98 6 49.24 1042.02 28.27 396.61 59.24 913.65 7 49.34 1047.96 28.37 400.68 59.34 914.15 8 49.44 1057.33 28.47 411.19 59.44 921.36 9 49.54 1057.83 28.57 423.03 59.54 921.86 10 49.72 1100.00 46.32 585.00 90.37 3000.00 92 Table 6.28. Results found for Firm 9 and 10 (m=10) Firm-9 Firm-10 Block i bi ($/MWh) qi (MWh) bi ($/MWh) qi (MWh) 1 46.41 343.69 45.34 454.83 2 46.51 347.94 53.26 875.22 3 46.61 351.35 57.32 1149.35 4 46.71 354.55 57.42 1216.91 5 46.85 368.00 57.52 1222.98 6 46.95 373.18 58.85 1255.57 7 47.05 383.41 58.95 1256.07 8 47.15 392.15 59.05 1258.21 9 47.25 399.01 59.15 1284.19 10 47.33 403.00 97.34 4400.00 Table 6.29. Unit?s profits (m=10) Unit Profit ($) Unit Profit ($) 1 46,867 6 38,681 2 254,940 7 38,759 3 51,019 8 24,115 4 190,774 9 69,412 5 112,786 10 93,351 Table 6.30. Profit increases in rule 2 (m=10) Rule 2 Firm 1 0.05% Firm 2 0.00% Firm 3 0.00% Firm 4 0.01% Firm 5 0.00% Firm 6 0.96% Firm 7 0.30% Firm 8 0.97% Firm 9 0.00% Firm 10 0.37% 93 6.6 Conclusion In this chapter, we showed that ABMS with dBPSO applied can be used to simulate the bidding process and to find a nash equilibrium solution. The defined stopping conditions and their results confirm the equilibrium conditions. In terms of computational time, relatively not much time is required to reach to the results. The model can further be applied to the real market if anyhow real operational cost and demand data is provided. 94 CHAPTER 7 CONCLUSION AND FUTURE RESEARCH In this research, models of electric power bidding in power markets were introduced. The fundamentals of a market design were described and related current literature was discussed. In order to evaluate a given bid, a spreadsheet based simulation algorithm was developed. The results found in the numerical examples were verified using the Bid Simulator. It was shown that the nonlinear and quadratic programming models were able to give an optimal solution for a limited number of price samples. However, due to the stochastic nature of market prices more price samples needed to be included. This limitation was overcome using two particle swarm optimization models. From experimental results, both approaches did not carry much computational time burden. Based on the statistical analysis, the decomposition based approach gave better results than the conventional particle swarm optimization. The bids obtained were compared with the marginal cost bidding method. The comparison showed that the quadratic term of the cost function plays an important role in determining the bid strategy. An additional method that models the bidding behaviors of power suppliers for a fixed demand was developed. The method uses the decomposition based particle swarm optimization and agent based simulation. Three stopping conditions to find the Nash 95 equilibrium were tested. The results found were analyzed using the equilibrium conditions needed for a competitive power market. Although the models described cover the fundamental process of bidding, they can be improved in future research. One improvement to the models is to include transmission constraints and congestion. Thus, the new models could be used to analyze the effect of transmission constraints on market prices as well as on bidding behaviors. Another improvement is to include operational constraints such as minimum up-down times, start-up costs and ramp-up limits of generating units. The models could provide a more realistic competition environment. Another extension could be considering operating reserve and contracts in the power market. Firms might choose to sell their power with fixed contracts or can be part of the operating reserve market. Thus, the bidding model could include these additional markets. The models developed in this dissertation consider uniform price auctions. However, there are markets that use the pay-as bid price auction mechanism. Therefore, new bidding models can be developed that include the pay-as bid auction type. Also, each power market has their own rules for the bidding process including number of blocks and the time period that the bid is valid for. The developed models can further be applied to those markets and perform economic comparisons. 96 REFERENCES [1] Ilic, M., Galiana, F., and Fink, L. (1998). Power Systems Restructuring: Engineering and Economics. Boston: Kluwer. [2] Cain, B.M, Alvarado, L.A., (2004) ?Implications of Cost and Bid Format on Electricity Market Studies: Linear Versus Quadratic Costs?, Large Engineering Systems Conference on Power Engineering, July, Halifax, Canada, pp. 1-5. [3] David, A.K, Wen, Fushuan. (2000). ?Strategic Bidding in Competitive Electricity Markets: a Literature Survey?, IEEE Conference Proceedings, 0-7803-6420. [4] Kian,A., Cruz, J.B.,Thomas, R.J., (2005). ?Bidding Strategies in Oligopolistic Dynamic Electricity Double-Sided Auctions?, IEEE Transactions on Power Systems, Vol. 20, No.1, pp. 50- 58. [5] Rajaraman, R., Alvarado, F., (2003) ?Optimal Bidding Strategy in Electricity Markets Under Uncertain Energy and Reserve Prices?, Power Systems Engineering Research Center Report, WI. [6] Valenzuela, J. and Mazumdar, M. (2003). "Commitment of Electric Power Generators under Stochastic Market Prices," Operations Research, Vol. 51, Issue 6, pp. 880-893. [7] Wen, F., David, A.K., (2001), ?Optimal bidding strategies for competitive generators and large consumers?, Electrical Power and Energy Systems 23 (2001) 37?43 97 [8] Liu, Y., Wu, F.F., (2006) ?Generator Bidding in Oligopolistic Electricity Markets Using Optimal Control: Fundamentals and Application?, IEEE Transactions on Power Systems, Vol. 21, No.3, pp. 1050- 1061. [9] David, A.K, Wen, F. (2000). ?Strategic Bidding in Competitive Electricity Markets: a Literature Survey?, IEEE Conference Proceedings, pp. 2168-2173. [10] Kian,A., Cruz, J.B.(2005). ?Bidding Strategies in dynamic electricity markets?, Decision Support Systems, Vol. 40, No:3-4, pp. 543? 551. [11] Song, Y., Ni, Y., Wen, F., Hou, Z., Wu, F. (2003). ?Conjectural variation based bidding strategy in spot markets: fundamentals and comparison with classical game theoretical bidding strategies?, Electric Power Systems Research, Vol. 67, Issue 1, pp. 45- 51. [12] Richter, C. W., Sheble, G..B., (1998), ?Genetic Algorithm Evolution of Utility Bidding strategies for the competitive marketplace?, IEEE Transactions on Power Systems, Vol. 13, No. 1, pp.256-261. [13] Rodriguez, C. P. , Anders, G.. J., (2004),?Bidding Strategy Design for Different Types of Electric Power Market Participants? IEEE Transactions on Power Systems, Vol. 19, No. 2, pp.964-971. [14] Zhang, D., Wang, Y., Luh, P.B.,(2000), ?Optimization Based Bidding Strategies in the Deregulated Market?, IEEE Transactions on Power Systems, Vol. 15, No. 3, pp.981- 986. [15] Wen, F., David, A.K.,(2001),? Optimal Bidding Strategies and Modeling of Imperfect Information Among Competitive Generators?, IEEE Transactions on Power Systems, Vol. 16, No. 1, pp.15-21. 98 [16] Conejo, A.J., Nogales, F.J., Arroyo, J.M., (2002), ?Price-Taker Bidding Strategy Under Price Uncertainity?, IEEE Transactions On Power Systems, Vol. 17, No. 4, pp. 1081-1088. [17] Ladurantaye, D., Gendreau, M., Potvin J., (2007), ?Strategic Bidding for Price- Taker Hydroelectricity Producers?, IEEE Transactions on Power Systems, Vol. 22, No. 4, pp. 2187-2203. [18] Fleten, S.E., Pettersen, E., (2005), ?Constructing Biding Curves for a Price- Taking Retailer in the Norwegian Electricity Market?, IEEE Transactions On Power Systems, Vol. 20, No. 2, pp. 701-708. [19] Oren, S., Sioshansi, R. (2005), ?Market Engineering and Applications?, pp. 115, available online at: http://www.ieor.berkeley.edu/~ramteen/ieor190d/handouts/book.pdf [20] Klemperer, P., (2004), ?Auctions: Theory and Practice?, Princeton University Press, pp. 9-26 available online at: http://www.nuff.ox.ac.uk/users/klemperer/VirtualBook/VirtualBookCoverSheet.asp [21] Contreras, J., Candiles, O., Fuente, J., Gomez, T., (2001), ?Auction Design in Day-Ahead Electricity Markets?, IEEE Transactions On Power Systems, Vol. 16, No. 3, pp. 409-417. [22] FERC, Notice of Proposed Rulemaking, (2002) ?Promoting Transmission Investment through Pricing Reform?, pp. 1-56, available online at: http://www.ferc.gov/industries/electric/indus-act/trans-invest.asp 99 [23] Hogan, W., (2002), ?Electricity Market Design and Structure: Working Paper on Standardized Transmission Service and Wholesale Electric Market Design?, Comments submitted to FERC, Docket No. RM01-12-000, pp.1-64, available online at: http://ksghome.harvard.edu/~WHogan/Hogan_FERC_041002.pdf [24] Oren, S., (2005), ?Market Design for Competitive Electricity?, pp.1-29, available online at: http://www.ieor.berkeley.edu/~ramteen/ieor190d/handouts/ elect_mkt_design.pdf [25] Ott, A.L. (2003).?Experience with PJM Market Operation, System Design, and Implementation?, IEEE Transactions on Power Systems, Vol. 18, No. 2, pp. 528- 534 [26] Alvarado, F.,(2002), ?Is System Control Entirely by Price Feasible??,Proceedings of the 36th Hawaii International Conference on System Sciences, pp.6 [27] Baldick, R., (2002), ?Electricity market equilibrium models: The effect of parametrization", IEEE Transactions on Power Systems, Vol. 17, No: 4, pp. 1170-1176 [28] Conejo, A., Castillo, E., Minguez, R., Milano, F., (2005), ?Locational Marginal Price Sensitivities?, IEEE Transactions on Power Systems, Vol. 20, No: 4, pp. 2026-2033 [29] Wood A., J. and Wollenberg B.F. (1996). Power Generation, Operation and Control, 2nd Ed. ? New York: John Wiley & Sons, Inc., pp. 31-32. [30] Maloney, M.T., (2001) ?Economies and Diseconomies: Estimating Electricity Cost Functions?, Review of Industrial Organization, Vol.19, No.2, pp. 165-180 [31] Baldick, R., (2006), ?Computing the Electricity Market Equilibrium: Uses of market equilibrium models?, Power Systems Conference and Exposition, pp. 66-73 [32] PJM Interconnections, (2007), http://www.pjm.com. 100 [33] Bouffard, F., Galiaba, F., (2006), ?Power System Security and Short-Term Electricity Market Equilibrium?, Power Systems Conference and Exposition, pp. 90-95 [34] Ventosa, M., Baillo, A., Ramos, A., Rivier, M., (2005), ?Electricity market modeling trends?, Energy Policy, Vol. 33, No:7, pp. 897-913. [35] Li, T., Shahidehpour, M., Keyhani, A., (2004), ?Market Power Analysis in Electricity Markets Using Supply Function Equilibrium Model?, IMA Journal of Management Mathematics, Vol. 15, no:4, pp.339-354. [36] Klemperer, P., Meyer, M., (1989), ?Supply Function Equilibrium?, Econometrica, Vol. 57, no:6, pp. 1243-1277. [37] Green, R.J., Newbery, D.M.,(1992), ?Competition in the British Electric Spot Market?, Journal of Political Economy, vol 100, no:5, pp. 929-953. [38] Baldick, R., Grant, R., Kahn, D., (2004), "Theory and Application of Linear Supply Function Equilibrium in Electricity Markets", Journal of Regulatory Economics, vol:25, No:2, pp.143-167 [39] Rudkevich, A., (2003). ?Supply function equilibrium: theory and applications?, Proceedings of the 36th Annual Hawaii International Conference on System Sciences, pp.10. [40] Hobbs, B.F., Helman, U., Pang, J., (2001), ?Equilibrium Market Power Modeling for Large Scale Power Systems?, Power Engineering Society Summer Meeting, Vol. 1, Issue: 2001 pp. 558 ? 563. [41] Holmberg, P., (2005), ?Asymmetric Supply Function Equilibrium with Constant Marginal Costs?, Working Paper : 16, available online at: http://www.nek.uu.se 101 [42] Nogales, F.J., Contreras,J., Conejo, A.J., Espinola, R., (2002), ?Forecasting next- Day Electricity Prices by Time Series Models?, IEEE Transactions on Power Systems, Vol: 17, no:2, pp.342-350 [43] Valenzuela, J. and Mazumdar, M. (2001). "Stochastic Monte Carlo Computation of Power Generation Production Costs under Operating Constraints," IEEE Transactions on Power Systems, Vol. 16, no. 4, pp. 671-677. [44] Kian, A. R., Keyhani,A.,(2001) "Stochastic price modeling of electricity in deregulated energy markets", presented and also published in the proceedings of HICSS- 34. [45] Xu, H.,Niimura T.,(2004). ?Short-Term Electricity Price Modeling and Forecasting Using Wavelets and Multivariate Time Series?, IEEE Transactions on Power Systems, vol.1, pp. 208- 212 [46] Pao, H.T.,(2006) ?Forecasting electricity market pricing using artificial neural networks?, Energy Conversion and Management, Vol. 48, no.3, pp. 907?912 [47] Pedregal, D.J., Trapero, J.R., (2007), ?Electricity prices forecasting by automatic dynamic harmonic regression models?, Energy Conversion and Management, Vol. 48, no.5, pp. 1710?1719 [48] Ruibal, C.M., Mazumdar, M., (2008), ?Forecasting the Mean and the Variance of Electricity Prices in Deregulated Markets?, IEEE Transactions on Power Systems, Vol. 23, No: 1, pp. 25-32. [49] Geman, H., Roncoroni, A. (2006). ?Understanding the Fine Structure of Electricity Prices?, Journal of Business, Vol. 79, No. 6, pp. 1225-1262. 102 [50] Kennedy, J. and Eberhart, R. C. (1995). ?Particle Swarm Optimization?. Proc. IEEE Int'l. Conf. on Neural Networks, IV, Piscataway, NJ: IEEE Service Center, pp.1942?1948. [51] Carlisle, A. and Dozier, G. (2001). ?An Off-The-Shelf PSO?, Proceedings of the 2001 Workshop on Particle Swarm Optimization, Indianapolis, IN, pp. 1-6. [52] Kennedy, J. and Eberhart, R. C. (2001). Swarm Intelligence, Proc. Of the Genetic and Evolutionary Computation Conference, Academic Press, 2001, pp. 75-100. [53] Energy Information Administration, (2008), http://www.eia.doe.gov/ [53] OMalley, M. J. and Liu, C. -C. (2001). Competitive wholesale electricity markets in Power System Restructuring and Deregulation: Trading, Performance, and Information Technology, Book Chapter, pp. 76-109. [54] Sueyoshi, T., Tadiparthi, G.R.,(2007), ?Agent-Based Approach to Handle Business Complexity in U.S. Wholesale Power Trading?, IEEE Transactions on Power Systems, Vol. 22, no. 2, pp. 532-543. [55] Niu, H., Baldick, R., (2005), ?Supply Function Equilibrium Bidding Strategies with Fixed Forward Contracts?, IEEE Transactions on Power Systems, Vol. 20, no. 4, pp. 1859-1867. [56] Macal, C.M., North, M.J.,(2006), ?Tutorial on Agent-Based Modeling and Simulation Part 2: How to model with Agents?, Proceedings of the 2006 Winter Simulation Conference, Orlando,FL, pp.2-15. [57] Bonabeau, E., (2002), ?Agent-based modeling: Methods and techniques for simulating human systems?, PNAS, Vol.99, no:3, pp.7280-7287. 103 [58] Axelrod, R., (2005), ?Agent-Based Modeling as a Bridge Between Disciplines?, Handbook of Computational Economics, Vol. 2: Agent-Based Computational Economics, Handbooks in Economics Series, North-Holland, forthcoming [59] Sueyoshi, T., Tadiparthi, G.R.,(2005), ?A Wholesale Power Trading Simulator With Learning Capabilities?, IEEE Transactions on Power Systems, Vol. 20, no. 3, pp. 1330-1340. [60] Tesfatsion, L., (2001), ?Agent-Based Modeling of Evolutionary Economic Systems?, IEEE Transactions on Evolutionary Computation, Vol. 5, no, 5, pp. 437-441. [61] Lau, J.S., Huang, G.Q., Mak, K.L., Liang, L., (2006), ?Agent-Based Modeling of Supply Chain Distribution Scheduling?, IEEE Transactions on Systems, Man, and Cybernetics-Part A:Systems and Human Systems, Vol:36, no.5, pp. 847-861 [62] LeBaron, B., (2002), ?Short-memory traders and their impact on learning group learning in financial markets?, PNAS, Vol.99, no.3, pp.7201-7206. [63] Nicolaisen, J., Petrov, V., Tesfatsion, L., (2000), ?Market Power and Efficiency in a Computational Electricity Market with Discriminatory Double-Auction Pricing?, IEEE Transactions on Evolutionary Computation, Vol. 5, no. 5, pp. 504-523. [64] Bagnall, A. J. ,(2000), ?A multi-adaptive agent model of generator bidding in the UK market in electricity,? in Proc. Genetic and Evolutionary Computation Conf., Las Vegas, NV, 2000, pp. 605?612. [65] Botterud, A., Thimmapuram, P.R., Yamakado, M., (2006), Argonne National Laboratory, Chicago, IL, available online at:http://www.dis.anl.gov/publications/ articles/ceeesa_EMCAS_IAEE72005NorwayPaper.pdf