University of Southern California Dissertations and Theses

Discrete optimization for supply demand matching in smart grids

(USC Thesis Other)

Discrete optimization for supply demand matching in smart grids

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content DISCRETE OPTIMIZATION FOR SUPPLY DEMAND MATCHING IN SMART GRIDS by Sanmukh Rao Kuppannagari A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (COMPUTER ENGINEERING) August 2018 Copyright 2018 Sanmukh Rao Kuppannagari Dedication To my parents for their sacrifices and support. ii Acknowledgments First and foremost, I would like to express my gratitude to my advisor Dr. Viktor Prasanna for his guidance and patience. The path to this dissertation had many turns with more than a few of them leading to dead ends. It is a testament to his mentor ship that I was able to complete my PhD with making enough mistakes to be able to learn from them but not enough to demoralize me. I would also like to thank Dr. Raj- gopal Kannan whose mentor ship greatly enhanced by mathematical maturity and prob- lem solving skills. The zest with which he approaches challenging computer science problems has inspired me greatly and I enjoyed our long discussions on such prob- lems. Moreover, I would like to thank my defense committee members Dr. Mohammed Beshir, Dr. Shaddin Dughmi and Dr. Mihailo Jovanovic for taking an interest in my work and providing with useful feedback which will help improve my dissertation and my research in future. I would also like to thank my qualifier committee members Dr. Paul Bogadan, Dr. Muhammad Naveed and Dr. Insoon Yang for their useful feedback on my qualifier proposal. I would like to thank Prof. Gandhi Puvvada for giving me a chance to teach his course for several semesters thereby improving my teaching skills. I would also like to thank Dr. Hemangee Kapoor and Dr. Krishnamachar Sreenivasan, who had the most impact on me during my under graduation and in my decision of pursuing a PhD. I would like to thank all the EE department staff members and student workers for making my iii life easier by handling all the administrative tasks, especially Kathy Kassar who would help me even after putting the “emergencies only” board, Janice Thompson who helped me get settled in Prof. Prasanna’s group and Tim Boston who helped me get settled in the department. Tim you will be missed. Additionally, I would like to thank my friends and colleagues for all the moral sup- port. I would like to thank Hiteshi whom I turned to during times of emotional distress, Vishnu for his encouragement and support and Sananda for all the long discussions on green buildings. I would like to thank Geoffrey for the countless discussions we had over a bowl of ramen. I would like to thank Manpreet and Surbhi for always being there for me on phone if not in person. I would like to thank all the members of USC Vidushak especially Mallika, Bhamati, Kiran, Satya, Yamini and Albin. I would like to thank people such as Krishna, Sundar etc. who were my roommates at some point or the other during my PhD for the great times I had with them. I would like to thank the members of Prof. Prasanna’s group for all the things I learned from them and all the good times I had with them. I would like to thank my friends Vipul, Megha, Somya, Raghav, Amulya and countless others for their support. And most importantly, I would like to thank Om and Suparna who were like my family away from home. The people I would like to thank the most are my parents: Shridhar and Sailaja. My PhD would not have been possible without all the sacrifices they made to make my life easy, allowing me to solely focus on my research. I would also like to thank all my family members especially my grandparents, my brother Vikrant, my sister Tanushree, my uncles and aunts for all their love and support. iv Contents Dedication ii Acknowledgments iii List of Tables ix List of Figures x Abstract xiii 1 Introduction 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Current Technology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Motivation for Discrete Operational States . . . . . . . . . . . . . . . . 3 1.4 Challenges in Optimizing over Discrete Operational States . . . . . . . 5 1.5 Application of Approximation Algorithms . . . . . . . . . . . . . . . . 6 1.6 Thesis Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.7 Research Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.8 Dissertation Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2 Related Work 10 2.1 Earlier Works on Demand Response . . . . . . . . . . . . . . . . . . . 10 2.2 Tractable Demand Response Solutions for Discrete Curtailment . . . . 11 2.3 Optimal Demand Response using Integer Programming . . . . . . . . . 12 2.4 Solar Curtailment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.5 Applications of Approximation Algorithms for Supply Demand Matching 14 3 Smart Grid Modeling for Discrete Supply Demand Matching Optimizations 15 3.1 Supply Demand Matching via Controllable Curtailment Strategies . . . 16 3.1.1 Smart Grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.1.2 Smart Grid Node . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.1.3 Smart Grid Controller . . . . . . . . . . . . . . . . . . . . . . 17 3.1.4 Demand Curtailment Strategies . . . . . . . . . . . . . . . . . 18 v 3.1.5 Solar Curtailment Strategies . . . . . . . . . . . . . . . . . . . 19 3.1.6 Curtailment Cost . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.1.7 Fairness in Curtailment . . . . . . . . . . . . . . . . . . . . . . 20 3.1.8 Incorporating Strategy Switching Overheads . . . . . . . . . . 20 3.2 Supply Demand Matching Framework . . . . . . . . . . . . . . . . . . 21 3.3 Modeling Smart Grid for Curtailment Selection . . . . . . . . . . . . . 25 3.4 Modeling Smart Grid Network . . . . . . . . . . . . . . . . . . . . . . 26 3.5 Modeling Storage for Supply Demand Matching Algorithms . . . . . . 29 3.6 NP-hardness of Supply Demand Matching Problem . . . . . . . . . . . 31 4 Optimal Curtailment Selection for Demand Response 33 4.1 Problem of Curtailment Strategy Selection for Demand Response . . . . 33 4.2 Traditional Demand Response . . . . . . . . . . . . . . . . . . . . . . 35 4.3 Sustainable Demand Response . . . . . . . . . . . . . . . . . . . . . . 36 4.3.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.3.2 ILP Formulation for Sustainable DR . . . . . . . . . . . . . . . 37 4.4 Sustainable DR with Strategy Overheads . . . . . . . . . . . . . . . . . 38 4.4.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.4.2 ILP formulation for Sustainable DR with Strategy Overheads . . 39 4.5 NO-LESS: Near OptimaLCurtailmEnt Strategy Selection for Supply Demand Matching in Micro Grids . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.5.1 Problem Definition . . . . . . . . . . . . . . . . . . . . . . . . 40 4.5.2 ILP Formulation for NO-LESS . . . . . . . . . . . . . . . . . . 41 4.6 Approximation Algorithms . . . . . . . . . . . . . . . . . . . . . . . . 42 4.6.1 Fast p 2-factor Approximation for Sustainable DR . . . . . . . 42 4.6.2 A PTAS for Sustainable DR . . . . . . . . . . . . . . . . . . . 44 4.6.3 FPTAS for NO-LESS . . . . . . . . . . . . . . . . . . . . . . . 44 4.7 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.7.1 Experimental Setup for TDR and SDR . . . . . . . . . . . . . . 51 4.7.2 Results and Analysis for TDR and SDR . . . . . . . . . . . . . 52 4.7.3 Sustainable DR . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.7.4 Sustainable DR with Strategy Overheads . . . . . . . . . . . . 54 4.7.5 NO-LESS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.8 Automated Dynamic Demand Response Implementation on USC Campus 69 4.8.1 D 2 R Implementation . . . . . . . . . . . . . . . . . . . . . . . 69 4.8.2 Implementation Challenges . . . . . . . . . . . . . . . . . . . 72 4.8.3 Overall DR Evaluation . . . . . . . . . . . . . . . . . . . . . . 72 5 Cost Optimal Supply Demand Matching 77 5.1 Minimum Cost Supply Demand Matching . . . . . . . . . . . . . . . . 77 5.1.1 ILP Formulation . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.1.2 Approximation Algorithm . . . . . . . . . . . . . . . . . . . . 79 vi 5.2 Minimum Cost Supply Demand Matching with Fairness . . . . . . . . 84 5.2.1 ILP Formulation . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.2.2 Approximation Algorithm . . . . . . . . . . . . . . . . . . . . 85 5.3 Online Algorithm for Fair Supply Demand Matching . . . . . . . . . . 88 5.4 Supply Demand Matching with Network Constraints . . . . . . . . . . 90 5.4.1 Transformer Level Supply Demand Matching . . . . . . . . . . 90 5.4.2 Smart Grid Level Supply Demand Matching . . . . . . . . . . 92 5.5 Modeling Storage for Supply Demand Matching Algorithms . . . . . . 99 5.6 Results and Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.6.1 Dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.6.2 Minimum Cost Supply Demand Matching . . . . . . . . . . . . 102 5.6.3 Improvement in PV Penetration . . . . . . . . . . . . . . . . . 108 5.6.4 Minimum Cost Supply Demand Matching with Fairness . . . . 110 5.6.5 Online Algorithm for Minimum Cost Supply Demand Matching with Fairness . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 5.6.6 Supply Demand Matching with Network Constraints . . . . . . 114 5.6.7 Storage Modeling . . . . . . . . . . . . . . . . . . . . . . . . . 116 6 Discrete Supply Demand Matching Under Prediction Uncertainty 118 6.1 Two Stage Stochastic Recourse Model for Discrete Supply Demand Matching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.1.1 Smart Grid Model . . . . . . . . . . . . . . . . . . . . . . . . 118 6.1.2 Supply Demand Matching without Prediction Uncertainty . . . 119 6.1.3 Incorporating Prediction Uncertainty in Supply Demand Matching120 6.2 Approximation Algorithms for Supply Demand Matching Under Pre- diction Uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 6.2.1 Latent State Model . . . . . . . . . . . . . . . . . . . . . . . . 122 6.2.2 Approximation Algorithm for 2-SRL . . . . . . . . . . . . . . 123 6.2.3 Black Box Model . . . . . . . . . . . . . . . . . . . . . . . . . 129 6.2.4 Approximation Algorithm for 2-SRB . . . . . . . . . . . . . . 131 6.2.5 Practical Implementation . . . . . . . . . . . . . . . . . . . . . 132 6.3 Risk Aware Supply Demand Matching in Smart Grids with High DER Penetration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.3.1 Smart Grid Model for Risk Aware Supply Demand Matching . . 134 6.3.2 Input Data Model . . . . . . . . . . . . . . . . . . . . . . . . . 135 6.3.3 Tail End Risk . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 6.3.4 Supply Demand Framework . . . . . . . . . . . . . . . . . . . 136 6.3.5 Risk Aware Supply Demand Matching Framework . . . . . . . 137 6.4 Results and Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 6.4.1 2-SR Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . 141 6.4.2 Risk Aware Supply Demand Matching . . . . . . . . . . . . . . 143 vii 7 Conclusion and Future Work 148 7.1 Broader Impacts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 7.2 Future Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 7.2.1 Discrete Multi Phase Supply Demand Matching . . . . . . . . . 149 7.2.2 Joint Optimization over Multiple Infrastructure for Smart Cities 151 7.3 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 Reference List 153 viii List of Tables 3.1 List of Variables in the Models . . . . . . . . . . . . . . . . . . . . . . 26 3.2 List of Variables in the Storage Model . . . . . . . . . . . . . . . . . . 30 4.1 Error and Scalability (20 nodes) Comparison . . . . . . . . . . . . . . 68 ix List of Figures 1.1 Continuous Operational Values: Node can take any operational value between 0 and 100. Discrete Operational Values: Node can only take the values 0, 20, 40, 60, 80 and 100. . . . . . . . . . . . . . . . . . . . 4 1.2 Supply demand matching in the above figure requires enumerating 3 5 = 243 combinations in the worst case . . . . . . . . . . . . . . . . . . . . 5 3.1 Smart Grid node with 3 + 1 (default) available curtailment strategies. . . 17 3.2 Smart Grid nodes with centralized controller communicating via proto- col such as ANM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.3 Strategy transition diagram for a node with 3 + 1 (default) strategy. A bi-directional arrow signifies that transition can occur in both directions. 21 3.4 Example Supply Demand Curve . . . . . . . . . . . . . . . . . . . . . 22 3.5 High Level Overview of Single Control Mode . . . . . . . . . . . . . . 23 3.6 High Level Overview of Dual Control Mode . . . . . . . . . . . . . . . 24 3.7 Smart Grid Network Model . . . . . . . . . . . . . . . . . . . . . . . . 27 4.1 Error incurred by ILP based Algorithm for various targeted Curtailment Values for Different DR Event . . . . . . . . . . . . . . . . . . . . . . 54 4.2 Error incurred by State-of-the-art Algorithm for various targeted Cur- tailment Values for Different DR Event . . . . . . . . . . . . . . . . . 55 4.3 Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Monday 1-3 pm . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.4 Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Monday 3-5 pm . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.5 Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Tuesday 1-3 pm . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.6 Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Tuesday 3-5 pm . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.7 Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Wednesday 1-3 pm . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.8 Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Wednesday 3-5 pm . . . . . . . . . . . . . . . . . . . . . . . . . 59 x 4.9 Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Thursday 1-3 pm . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.10 Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Thursday 3-5 pm . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.11 Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Friday 1-3 pm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.12 Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Friday 3-5 pm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.13 Absolute errors in kWh incurred for Sustainable DR event on Monday 1-5 pm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.14 Absolute errors in kWh incurred for Sustainable DR event on Tuesday 1-5 pm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.15 Absolute errors in kWh incurred for Sustainable DR event on Wednes- day 1-5 pm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.16 Absolute errors in kWh incurred for Sustainable DR event on Thursday 1-5 pm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.17 Absolute errors in kWh incurred for Sustainable DR event on Friday 1-5 pm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.18 Comparison of demand curtailment values achieved in each interval for SDR and TDR for the targeted curtailment value of 1200 kWh . . . . . 64 4.19 Near Optimality of NO-LESS . . . . . . . . . . . . . . . . . . . . . . 65 4.20 Execution Time vs Number of nodes for different values of epsilon . . . 66 4.21 Execution Time vs epsilon for different number of nodes . . . . . . . . 67 4.22 Execution Time vs Number of Buildings for = 0:01 . . . . . . . . . . 67 4.23 Control and Data Flow For D 2 R Implementation . . . . . . . . . . . . . 70 4.24 Policy Engine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.25 Consumption profile of building THH during DR . . . . . . . . . . . . 73 4.26 Consumption profile of building VKC during DR . . . . . . . . . . . . 74 4.27 Consumption profile of building MRF during DR . . . . . . . . . . . . 75 4.28 Consumption profile of building WPH during DR . . . . . . . . . . . . 75 4.29 Consumption profile of building BHE during DR . . . . . . . . . . . . 76 4.30 Consumption profile of building DRB during DR . . . . . . . . . . . . 76 5.1 Percentage Error of the curtailment target constraints for = 0.5-0.1 (50-10% Error Guarantee) . . . . . . . . . . . . . . . . . . . . . . . . 103 5.2 Percentage Error of the curtailment target constraints for = 0.05-0.02 (5-2% Error Guarantee) . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.3 Ratio of the cost of solution obtained by MinCB and the optimal solution versus the approximation guarantee () . . . . . . . . . . . . . . . . . . 105 5.4 Runtime of MinCB versus: (a) Number of nodes for various values of epsilon . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 xi 5.5 Runtime of MinCB versus: (b) (denoted as percentage error) for vari- ous values of the number of nodes . . . . . . . . . . . . . . . . . . . . 106 5.6 Improvement in PV Penetration . . . . . . . . . . . . . . . . . . . . . 108 5.7 Percentage Error of the Cost of MinCBF w.r.t. the Optimal Cost . . . . 110 5.8 Percentage of Budget Overshoot for 10 Worst Affected Nodes . . . . . 111 5.9 Percentage Undershoot of Interval Curtailment Target for various inter- vals (note: the labels of x axis do not denote the actual interval number) 112 5.10 Gini Coefficient for various values of Alpha () . . . . . . . . . . . . . 113 5.11 Percentage Error of the Cost of Online MinCBF w.r.t the Optimal Cost . 114 5.12 Ratio of Cost of Algorithm 10 to Algorithm 9 versus Feeder Capacity for various values of Distributor capacity. . . . . . . . . . . . . . . . . 115 5.13 Percentage Increase in Runtime For various number of storage nodes for varying epsilon. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 6.1 Supply Demand Matching Framework . . . . . . . . . . . . . . . . . . 137 6.2 Runtime of 2-SRL-Approx for varying values of nodes for fixed = 0:2 142 6.3 Runtime of 2-SRL-Approx for varying values of for fixed number of nodes = 20 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 6.4 Load and Solar Timeseries Data Used for Evaluation . . . . . . . . . . 145 6.5 Relative Percentage Error of Framework w.r.t. Optimal Cost For Vary- ing Initial Storage and Variance . . . . . . . . . . . . . . . . . . . . . . 146 6.6 Time Required for Decision Making w.r.t Storage Capacity for Varying Variance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 xii Abstract Traditional power grids consisted of large grid operator owned generators to supply the demand of the consumers. The inertia of the large generators was able to handle any sudden demand fluctuations followed by an eventual ramp up/down of generation. How- ever, modern smart grids, characterized by the proliferation of small sized Distributed Energy Resources (DER) such as photovoltaics, energy storage etc. lack this ability of handling demand fluctuations due to their lack of inertia. Therefore, they require proac- tive supply demand matching by utilizing the available advanced sensing and remote control capability. A plethora of work has been done assuming that the nodes (supply and demand) in the grid exhibit continuous operational values (generation/consumption values). Fine grained control is required for continuous operational values which might not be eco- nomically available. Thus, supply demand matching often needs to be performed over discrete sets of operational values of the nodes making this an NP-hard problem. While heuristics (with unbounded error) and Integer Program based algorithms (computation- ally expensive) have been studied extensively, the application of approximation algo- rithms to develop fast algorithms with bounded error guarantees has been limited. In this work, we address these limitations by developing a suite of polynomial run- time complexity approximation algorithms to perform supply demand matching in smart grids. We first model the smart grid to capture the available discrete control strategies xiii for each node in the smart grid and the associated costs and uncertainties. We also model practical grid constraints such as the overhead associated with switching strate- gies and the smart grid network capacity constraints. We then develop approximation algorithms with a variety of objectives and constraints such as minimizing curtailment error in Demand Response, minimizing cost of supply demand matching, fairness etc. and theoretically provide the optimality guarantees. In addition to the theoretical analy- sis, we perform practical evaluations of the algorithms and show that our supply demand matching algorithms provide solutions which are close to optimal in a small amount of time. xiv Chapter 1 Introduction 1.1 Motivation Electrical power grids have undergone a drastic transformation since the 1970s in terms of both scale and complexity [69]. Technological advances such as the use of bi- directional AMI meters, allowing real time remote monitoring and control, have trans- formed them into smart grids [48]. Continuous matching of supply and demand is the most critical grid operation. The electricity demand of the consumers fluctuates throughout the day. If the fluctuating demand is not matched using an equal amount of supply, it can lead to an increase (if supply is more than demand) or decrease (if demand is more than supply) in the grid frequency. If the variation in the frequency is not contained quickly, the stability of the grid is severely compromised with potential over-voltages or blackouts. Traditional power grids are characterized by the presence of a centralized grid oper- ator owning large generation units such as steam or water powered turbines to meet the demand. Any fluctuation in frequency due to demand variability is offset by the inertia of the turbines, with a loss in its kinetic energy. The grid operator can then increase or decrease the rotational speed of the turbines to match the demand. The larger the size of the generation units, the larger is the inertia. Larger inertia provides the ability to offset fluctuations for longer times. This allows the grid operators ample amount of time to ramp up/down the generation. 1 Technological advances in modern smart grids have enabled rapid proliferation of small sized consumer owned Distributed Energy Resources (DER) into the grid such as Photovoltaics (PV), Energy storage, fuel cells etc. For example, the adoption of dis- tributed solar energy has increased dramatically due to the falling cost of solar PVs. The installed prices of U.S. residential and commercial PV systems declined 5-7% on average during 1998-2011 [6]. As per the DoE SunShot vision document, solar gener- ated power is expected to grow to 14% of the total power supply in 2030 and 27% by 2050 [67]. As a result, the energy generation in the smart grids has decentralized. The increase in DERs has significantly complicated the operation of supply demand matching due to the following reasons: Generation patterns of renewable energy based DERs are highly influenced by the weather conditions. Hence, in addition to the variability in demand, variability in supply needs to be considered. Grids where a major portion of the demand is met using small sized DERs, have low inertia to handle fluctuations. Thus, less time is available to ramp up/down the supply. The distributed nature of DERs increases the complexity of control as the grid operator needs to consider constraints such as cost optimality, fairness, network constraints etc. for DER scheduling. In addition to supply, the proliferation of controllable loads has enabled consumer participation using techniques such as Demand Response (DR), thereby provid- ing a great opportunity as well as challenge for the operation of supply demand matching. Hence, sophisticated algorithms which pro-actively perform supply demand match- ing in each interval - where an interval is a grid operator determined period of time - 2 by controlling both the demand and supply are required to address the challenges listed above. 1.2 Current Technology A plethora of theoretically sound polynomial runtime complexity optimization algo- rithms have been developed by the smart grid community to perform supply demand matching in a cost effective manner [71]. However, these algorithms assume that the supply and demand nodes exhibit continuous operational values i.e. the generation value from a supply node or a consumption value from a demand node can take any value between its minimum and maximum. Using this assumption, the problem is formulated as linear or convex optimization problem [20] which admit fast polynomial runtime algorithms which provide optimal solutions using methods such as interior point [20]. However, the limitations of these works is that the assumption of continuous operational values requires fine grained control over the nodes which is not economical and might not even be feasible in certain scenarios. 1.3 Motivation for Discrete Operational States To motivate the need for modeling and optimizing on discrete operational states, we consider a toy example. Consider a warehouse which consists only of 100 dimmable LED bulbs each consuming a maximum power of 1 W. We assume that using a single O(1) time complexity operation, the consumption value of an LED can be set to any value between 0 and 1 W. The minimum consumption of the warehouse is 0 W, when all LEDs are set to 0 and the maximum is 100 W, when all LEDs are consuming maximum power (Figure 1.1). Hence, to obtain consumption values which are continuous in the 3 Figure 1.1: Continuous Operational Values: Node can take any operational value between 0 and 100. Discrete Operational Values: Node can only take the values 0, 20, 40, 60, 80 and 100. range 0 to 100 for the warehouse, we need 100 different switches to control LEDs. Assuming identical LEDs, this requires 100 control signals. Now consider the alternate scenario where a set of 20 LEDs are connected to a single switch. Consumption values of 0, 20, 40, 60, 80 and 100 can be obtained under such scenario using only 5 switches. Even if we assume each set to be unique, the number of control signals required are 2 5 = 32 which is much less than 100. In a real world grid, cost of installation is a major factor for a supply or demand node. Hence, even though it might be technically feasible to obtain continuous operational values, it might not be economically viable. Therefore, algorithms for supply demand matching need to operate efficiently over operational values for each node which form a set of discrete values. 4 1.4 Challenges in Optimizing over Discrete Operational States The discrete nature of the operational values transforms the problem of supply demand matching into a NP hard problem (see Section 3.6 for a formal proof). Thus, no optimal polynomial runtime complexity algorithm can exist unless P = NP. In the worst case, a brute force search enumerating all possible solutions, which is exponential in the input size might be required. See Figure 1.2. This approach is followed by current techniques which use Integer Program based formulations. Thus, they provide computationally expensive optimal solutions. This is clearly not possible for grid operations where strict timing constraints need to be adhered. To address the computational complexity, several works have focused on developing fast heuristics. Such heuristics do not provide any guarantee on optimality or on the worst case performance. They can incur unbounded errors in the objective value of the solutions they produce. Power grid operations require high reliability and hence such heuristics cannot be relied upon due to their uncertainty. Figure 1.2: Supply demand matching in the above figure requires enumerating 3 5 = 243 combinations in the worst case 5 1.5 Application of Approximation Algorithms One of the most significant contribution of this work is to model the various supply demand matching problems as variants of packing problems. Several packing problems such as 0-1 knapsack, subset sum problem etc. admit dynamic programming based Fully Polynomial Time Approximation Schemes (FPTAS) or Polynomial Time Approxima- tion Algorithms (PTAS) [72]. Given an accuracy parameter, an FPTAS is an algorithm with a polynomial runtime complexity in the input size and 1 which outputs a solution with an objective value which is within (1) of the optimal value. PTAS provides similar guarantees, however, its runtime complexity is allowed to be exponential in 1 . Modeling the supply demand matching problems as variants of packing problems allows us to use the approximation techniques used for the latter problem. Consider the following simple problem: We are given a set of supply nodes with nodei exhibiting exactly two operational values, 0 ands i with associated costs 0 andc s i . Similarly, we are given demand nodes with operational values, 0 andd i with associated costs 0 and c d i . Now if we assume constant demand and control the supply to meet the demand and ignore the costs, it can be modeled as the subset sum problem. If we consider the costs, the problem can be modeled as the knapsack problem. If we consider costs and control both supply and demand, the problem becomes more complicated, but still admits a dynamic programming based FPTAS. The problems that we consider in this work are higher in complexity as they opti- mize over several time periods, with each node equipped with more than two operational values and they incorporate several constraints and objective values. However, by care- fully analyzing the structure of the problem, we are able to leverage several techniques such as rounding, LP relaxation etc. to develop approximation algorithms and provide guarantees both in terms of the achieved objective values and the constraint violations. 6 1.6 Thesis Statement Rapid proliferation of Distributed Energy Resources (DER) in modern power grids, also known as “smart grids”, enabled by advancement in sensing and mon- itoring technology has introduced unique challenges as well as opportunities for optimal grid operations. Proactive supply demand matching by controlling the DERs and loads is imperative to ensure reliable grid operations. Moreover, as fine grained control is not available economically, supply demand matching algo- rithms need to optimize on discrete sets of operational values. Our research focuses on modeling the modern smart grids consisting of controllable loads, DERs and infrastructure nodes with associated capacity to enable development of discrete supply demand matching algorithms. Our goal is to develop polynomial runtime complexity algorithms for supply demand matching with bounds on the worst case performance both in terms of the achieved objective values and the constraint vio- lations. Our algorithms will enable the development of a smart grid control frame- work which will facilitate DER proliferation and reduce the cost of grid operations while ensuring reliability. 1.7 Research Contributions In this work, we develop fast polynomial runtime complexity approximation algorithms for a variety of objective functions and constraints. We provide sound theoretical anal- ysis of optimality guarantees supplemented with practical evaluations to show that our algorithms provide solutions which are close to optimal in a small amount of time. Specifically, our contributions are as follows: 1. We develop a smart grid model to be utilized by the supply demand matching algorithms to optimize over discrete set of operational values. Additionally, we 7 model the network of the smart grid and capture the capacity constraints of the infrastructure such as transformers and feeders. 2. We develop a suite of supply demand matching approximation algorithms for Demand Response. This is the scenario in which the supply is constant and the demand needs to be controlled (curtailed) such that its value is less than the avail- able supply. We develop algorithms with several objectives such as coarse grained curtailment across the entire horizon of the Demand Response event (Traditional Demand Response) or equitable distribution across all the intervals in the event (Sustainable Demand Response). We also consider several practical constraints such as the overheads incurred by each node in switching strategies and fairness of curtailment across all the nodes. 3. We develop a suite of supply demand matching approximation algorithms under the scenario when both supply and demand are controllable with the objective of minimizing the cost. We again consider several practical constraints such as fairness and network capacity constraints. 4. We model storage to incorporate into the algorithms developed without signifi- cantly affecting the computational complexity. 5. We develop algorithms to perform supply demand matching under prediction uncertainties. 1.8 Dissertation Outline The rest of the dissertation is organized as follows: We discuss the current works that focus on supply demand matching in Chapter 2. We then discuss our smart grid model for discrete supply demand matching optimizations in Chapter 3. This is followed by a 8 discussion of a suite of algorithms developed for performing Demand Response, where supply is constant and load is controllable in Chapter 4 and then the algorithms devel- oped to perform cost optimal supply demand matching with both controllable loads and supply are discussed in Chapter 5. Discrete supply demand matching under prediction uncertainty is discussed in Chapter 6. Finally, we conclude the dissertation with a dis- cussion of future works in Chapter 7. 9 Chapter 2 Related Work As discussed in the previous chapter, several works have focused on developing theo- retically sound algorithms to ensure reliability in smart grids by matching supply and demand. However, the focus of such works has mostly been on the nodes with contin- uous operational values. We refer the readers to the survey paper [71] which performs an exhaustive study of the DR techniques developed in the literature. As our focus in this work is on optimizing over discrete sets, we limit our literature survey to similar works. The works targeting discrete operational states either develop computationally expensive Integer program based solutions or develop fast heuristics with no error guar- antees. In this chapter, we perform a literature review of such works targeting discrete operational states. We also discuss some works which applied approximation algorithm techniques to solve some specific problems related to supply demand matching in smart grids. 2.1 Earlier Works on Demand Response Significant literature exists addressing the challenges, solutions, implementations and estimation methodology for calculating the energy savings for Demand Response (DR) [11, 51]. Early works focused on DR scheduling for individual residential cases [45] or household appliances [34]. These approaches are not scalable to smart grids. 10 Traditionally, DR algorithms have focused on targeting customers based on aggre- gate consumption data, relying on customer selection using billing data or sur- veys [43, 49], employing dynamic programming techniques for load management and minimizing peak load over a period [26], particle swarm optimization based tech- niques [61] and game theoretical solutions constrained by real time pricing [23] and customer comfort levels [17]. However, with data available from smart meters, work such as [64] show that such approaches are very inaccurate. The actual consumption data over a period differs significantly from the data obtained from surveys or billing cycles. Moreover, the selection is done oblivious to the distribution of load through- out the day. Therefore such approaches contribute little to reducing the peak energy consumption and distributing it over other periods. 2.2 Tractable Demand Response Solutions for Discrete Curtailment In order to tractably solve the problem of achieving load curtailment using discrete cur- tailment values, stochastic optimization algorithms have been developed in works such as [38] and [24]. However, such works assume the availability of a large number of nodes in the grid. For example, the authors in [38] propose a stochastic knapsack based algorithm for selecting customers to maximize the probability the desired curtailment value is achieved over the period of the entire DR event while limiting the utility’s cost. The algorithm relies on the central limit theorem to assume the joint customer response is normally distributed and thus is conditioned on the assumption that there are a large number of customers from whom a subset can be selected. For example, in order to achieve high reliability (> 95%) for large curtailment targets, the algorithm requires aroundN = 2000 customers which increases the computation cost. Two heuristics are 11 proposed but the approximation bounds are dependent on the ratio of customer demand response variances which could be arbitrarily large. Although fine grained data is avail- able, the selection algorithm focuses on aggregating the curtailment over a time interval. However, as mentioned before, such an approach can aggressively curtail load in some intervals and concentrate peak in other intervals. Moreover, by maximizing the proba- bility of exceeding the desired curtailment, there might be cases where they overcom- pensate leading to reduced benefits to the utility. Techniques such as [57], [41], [76] and [58] develop fast algorithms which can have arbitrarily large errors in the objective function (utility maximization, cost minimization etc.). Authors in [57] develop a genetic algorithm based heuristic while [76] presents a heuristic based on change making. The algorithm developed in [58] uses Linear Pro- gramming whose solutions need to be rounded to integral values and can have large errors (unbounded integrality gap). The notion of achieving sustainable DR over a peak period divided into subintervals was proposed in [76] using a change making heuris- tic to evenly distribute curtailment over intervals. However, a detailed analysis shows that it achieves consistency between intervals without reference to the target leading to unbounded errors which is also demonstrated by our experimental results. An automated framework employing this heuristic is developed in [77]. 2.3 Optimal Demand Response using Integer Program- ming The other approach is to provide computationally expensive exact solutions, for exam- ple, [28], [46], [18], [75], [25] and [16], where the authors use Integer/Mixed Integer Linear/Non-Linear Programming for their algorithm. The authors in [65] develop a quadratic programming formulation for Demand Response which is then reduced into 12 a distributed algorithm. The paper assumes the availability of continuous curtailment values. However, the practical DR implementation in USC’s smartgrid allows only dis- crete curtailment values which forms the basis of modeling them as such in our work. This constraint makes the problem more complicated as the linear/quadratic program gets converted into an integer linear/quadratic program. 2.4 Solar Curtailment Load Curtailment techniques are ineffective when supply due to solar PVs exceeds the demand. If this is left unmitigated, it causes over-voltages in the system leading to fail- ures. Several works perform reactive solar curtailment in response to rising voltage. VVO [52] increases reactive power to lower the voltage due to real power while iPlug curtails the solar energy input to the grid by redirecting it to charge storage or coordinate with local demand ramp-up resources [56]. The authors in [40, 70, 63] achieve contin- uous curtailment from solar PVs by running them at voltages other than the Maximum Power Point (MPP). This requires fine grained control of the solar panels. For some scenarios, fine grained control might not be available due to limitations of inverter tech- nology. Our work addresses such scenarios through a curtailment model that handles a discrete set of curtailment values and provides bounded polynomial time approxima- tions for achieving discrete curtailment targets. As mentioned in [33], discrete solar curtailment can be performed by simply disconnecting individual PV modules using the micro-inverters installed at PV installations. As opposed to the technique devel- oped in [33], which is reactive to over-voltages, we perform proactive solar curtailment. Works such as [47] develop a scheduling algorithm which determines the schedule of generation from various sources for the entire day. The sources include a mix of renew- able and non-renewable sources and the objective is to minimize cost. The resulting 13 optimization problem is complex as it has to make a decision for each node for each interval of the entire day. This makes it difficult to scale this algorithm for large distri- bution grid or micro-grids with rooftop PVs installed. 2.5 Applications of Approximation Algorithms for Sup- ply Demand Matching As discussed before, the application of approximation algorithms for the problem of supply demand matching is limited. Constant factor approximation algorithms based on strip packing have been developed for peak demand reduction [68, 10, 73, 54]. The problem of strip packing is defined as follows: Given a set of axis aligned rectangles, pack them into fixed width strip while minimizing its height. In smart grid context, the rectangles denote the total energy consumed by each individual appliance with the height of the strip denoting the maximum power (energy) required at any given time. Essentially, these algorithms distribute a fixed amount of energy temporally to minimize the peak consumption. Peak reduction is one special case of the more generalized supply demand matching problem. 14 Chapter 3 Smart Grid Modeling for Discrete Supply Demand Matching Optimizations In this work, we focus on a specific type of smart grid known as discrete microgrids. Discrete microgrids are characterized by the presence of a centralized controller which is responsible for the grid operations. The various supply and demand nodes in such smart grids exhibit discrete operational values i.e. at any given time step, the feasi- ble region of the operational value the node can exhibit forms a discrete set. More- over, due the presence of decentralized generation using small sized Distributed Energy Resources (DER), such grids have very low supply inertia, thereby making proactive supply demand matching imperative. In this chapter, we discuss the smart grid model that we developed to model discrete microgrids. The smart grid model captures important grid characteristics such as the number of supply and demand nodes, the available controllable strategy space and their associated costs, the costs associated with switching strategies, the smart grid network and the capacity constraints of the various infrastructure nodes such as transformers, feeders etc. This model is then utilized by the algorithms developed in the subsequent chapters for supply demand matching. 15 3.1 Supply Demand Matching via Controllable Curtail- ment Strategies 3.1.1 Smart Grid The Smart Grid that we consider in this work consists of several demand nodes: con- sumers of electricity, and supply nodes: electricity producers. The supply nodes that we consider in this work are the customers who have solar PVs installed. A single consumer can act both as a demand node and a supply node. In this case, we model two different nodes for the consumer: one for demand and one for supply. We assume that the smart grid has a high PV penetration i.e. the supply from solar PVs under normal weather conditions meet the demand of the consumers for most of the day. We assume that dur- ing night or during extremely unfavorable weather conditions, conventional sources of electricity are used to meet the demand. 3.1.2 Smart Grid Node As discussed in the previous section, a smart grid node can be a demand node or a supply node. We assume that for each node at any given interval, the operational value when the node is not being controlled is the maximum output for that interval. Any control action is taken in terms of curtailment i.e. a reduction in the operational value. In Sections 3.1.4 and 3.1.5, we discuss how to obtain these output values. We use the term curtailment strategy to denote a curtailment control action. Consider the node shown in Figure 3.1. At a given time t, if no control action is taken, or in other words, default strategy of 0 is followed, the operational value of the node is O(t) with an associated curtailment value of 0. Following strategies 1, 2 or 3 results in an operational value of 1 O(t), 2 O(t) and 3 O(t) respectively, where 16 0 3 2 1 1. The corresponding curtailment values are (1 1 )O(t), (1 2 )O(t) and (1 3 )O(t) respectively. The shaded portion below show the uncon- trollable operational state. For example, for a demand node, this could denotes the baseline demand which is imperative to be met. For some nodes, it is possible that the uncontrollable operational state area is empty. For example, for solar PVs it is possible to curtail the entire generation by disconnecting the PV from the grid. Figure 3.1: Smart Grid node with 3 + 1 (default) available curtailment strategies. 3.1.3 Smart Grid Controller We assume there exists a centralized grid operator/controller with the capability of remotely switching a node into a curtailment strategy. One simple way in which this capability can be achieved is as follows: Each node can be equipped with a controller with the capability to receive and respond to some IP based smart grid specific com- munication protocol such as Active Network Management (ANM). The centralized grid operator/controller can then send signals for switching curtailment strategies using the same protocol. 17 Figure 3.2: Smart Grid nodes with centralized controller communicating via protocol such as ANM 3.1.4 Demand Curtailment Strategies The demand curtailment model considered in this work is based on a real world Demand Response implementation in USC’s Campus Microgrid [36]. In the Smart Grid, each demand node is associated with several demand curtailment strategies. Examples of strategies include Global Zone Temperature Reset (GTR), Variable Frequency Drive Speed Reset (VFD), Equipment Duty Cycling (Duty) and their combinations [36]. Each curtailment strategy for each node in a given time interval exhibits a discrete curtailment value i.e. the reduction in the consumption from the maximum output due to the adop- tion of this curtailment strategy. This complicates the problem as Linear Programming based techniques, which are both fast and optimal can no longer be used. This value can be predicted using algorithms mentioned in [12]. Each node is also associated with a default curtailment strategy of curtailment value 0. Hence, if a node has no curtail- ment strategy available or does not participate in demand curtailment, we assume that it 18 follows the default strategy. Demand curtailment is known in the literature as Demand Response [11]. 3.1.5 Solar Curtailment Strategies Each supply node i.e., a node with PV installation in the Smart Grid consists of several solar panels (each solar panel is called a module). Traditionally, modules are connected in series to an inverter which in turn is connected to the grid. However, this topology affects the efficiency of the PV system as the inverter conditions the output according to the poorest performing module [33]. Therefore newer designs, in which each module is independently connected to a micro-inverter are becoming increasingly popular [27]. Technically, each micro-inverter of a PV installation is an independent grid con- nected generator, turning the PV installation into a segmented generator with discrete generation output. The maximum output of the PV installation will occur when all the solar panels are allowed to feed into the grid. However, at any given time, micro- inverters can be configured such that only a subset of PV modules are connected to the grid. Our objective is to exploit this capability by controlling the micro-inverter con- figuration and enabling discrete curtailment of supply. We refer the reader to [33] for more details on utilizing micro-inverters for solar curtailment. Note that the technique developed in [33] is a reactive technique which reacts to voltage increase and requires high frequency voltage sampling. Our technique is a proactive technique which avoids an increase in voltage by reducing supply in advance. As discussed in [33], the curtailment can be achieved using module tripping by exe- cuting it as an inverter restart. Module tripping is a form of generator tripping where the loss of capacity is less than 100%. Generator tripping is recognized as the most effective way of resolving transient stability issues [33]. Thus, in this work we do not explicitly consider any stability issues related to solar curtailment. 19 3.1.6 Curtailment Cost Each curtailment strategy for each node is associated with a cost value as curtailment leads to a loss in utility. These costs are determined by the grid operator to reflect the loss. Typically, the costs are some function of the curtailment value e.g., if a node, by following a strategy curtails , then the cost of this strategy will bef( ), wheref is some function determined by the grid operator. Linear and quadratic functions are commonly used cost functions in grid operations. The objective of our framework is to minimize cost while performing supply demand matching. 3.1.7 Fairness in Curtailment Curtailment performed by any nodes leads to a loss in utility. Hence, we should ensure that no single node should be disproportionately penalized while achieving a certain global curtailment target for supply demand matching. To ensure that certain nodes of the grid are not unfairly penalized, we assign a budget value B b with each node b. We also assign a lower bound of b B b on the curtailment value. A fairness constraint is added into the algorithms which ensures that no node curtails more than its budget value. 3.1.8 Incorporating Strategy Switching Overheads In order to incorporate strategy switching overheads, for a nodeb, we define a function b to model the cost of switching between allowable strategies. The purpose of the cost is to disallow frequent strategy switching when ramp time is high. If the node is allowed to switch from strategy j to k, then b (j;k) <1 denotes the cost of switching, else b (j;k) =1. Note that b can be represented using a 2D square matrix whosei;j entry represents the cost of switching from strategyi toj. This ensures that a single call of 20 b requiresO(1) amount of time. For a nodeb, we limit the cost incurred in switching strategies by b . For example, consider the state transition diagram shown in Figure 3.3. The follow- ing matrix represents the strategy switching overheads: 0 1 2 3 1 0 1 1 2 1 0 1 3 1 1 0 (3.1) Figure 3.3: Strategy transition diagram for a node with 3 + 1 (default) strategy. A bi- directional arrow signifies that transition can occur in both directions. 3.2 Supply Demand Matching Framework Mitigating supply-demand mismatch within tight timing constraints is critical for smooth operation of a smart grid. As shown in Figure 3.4, during several intervals of the day, such as regions 1 and 3, the demand of the consumers can exceed the solar supply. This can cause blackouts in the grid. Demand curtailment strategies need to be adopted during such intervals to avoid blackouts. The other extreme is shown using region 2 in Figure 3.4. These are the intervals in which the supply due to solar PVs exceeds the demand. This can cause over-voltages in the grid leading to the tripping 21 of fault prevention devices [33]. Under this scenario, supply curtailment strategies are required. Region 1: Demand > Supply Region 3: Demand > Supply Region 2: Supply > Demand 25 27 29 31 33 35 37 6:00 7:00 8:00 9:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 MegaWatts Time of the Day Example Demand Supply Curve Demand Supply Figure 3.4: Example Supply Demand Curve We develop a generalized framework which performs supply demand matching by selecting load or supply curtailment strategies. We define supply demand matching horizon as the time horizon during which the framework is used. Now, each supply node of the smart grid can be associated with a generation prediction model such as ARIMA+ANN ensemble [53] to output the generation value. Supply curtailment values corresponding to the strategies can be obtained as discussed in Section 3.1.5. Simi- larly, each demand node can be associated with a demand prediction model such as ARIMA [14] to output the demand value under default strategy. Each node can also be associated with a curtailment prediction model [12] which determines the curtailment obtained by following curtailment strategies. Determining the best prediction model for each node is a separate research topic and is out of the scope of this work. 22 Smart Grid model Calculate aggregate consumption and prediction for each interval Identify curtailment horizons Identify load curtailment targets and generate curtailment and cost matrices Identify supply curtailment targets and generate curtailment and cost matrices Run Net Load Balancing Algorithm for load curtailment strategy selection Run Net Load Balancing Algorithm for supply curtailment strategy selection Load > Supply For Each Curtailment Horizon Load Curtailment Strategies Supply Curtailment Strategies Figure 3.5: High Level Overview of Single Control Mode Our supply demand framework runs in two different modes differentiated by the nodes that are being controlled in any given time interval of the supply demand matching horizon. The two modes are as follows: 1. Single Control Mode: In this mode, at any given time interval, either demand nodes are controlled or the supply nodes, but not both. A high level overview of single control mode is shown in Figure 3.5. The framework determines the aggre- gate load and supply for each interval in the supply demand matching horizon. It then identifies a list of load curtailment horizons and supply curtailment hori- zons and the respective curtailment targets. A curtailment horizon is defined as a period of time during which either demand is higher than the supply requiring a demand curtailment or vice-versa. For each interval of every curtailment horizon, the supply demand matching framework uses the curtailment prediction models to determine discrete curtailment values for each node. It also determines the cost 23 values associated with them. Then, for each curtailment horizon, it runs one of the supply demand matching algorithms developed in this work. The algorithm to run is pre-determined by the grid operator. The algorithm returns the curtailment strategies to be followed by each node in each time interval of each curtailment horizon. 2. Dual Control Mode: In this mode, at any given time interval, both the demand and supply nodes are controlled. In each interval, the framework determines the cost of achieving all possible supply demand matching values and chooses the one which satisfies all the constraints with minimum cost. See Figure 3.6. The simplicity of single control mode allows us to optimize over several intervals, which is not possible in dual control mode. However, dual control mode can lead to more optimized solutions within a single interval by controlling both the supply side and the demand side. Figure 3.6: High Level Overview of Dual Control Mode 24 3.3 Modeling Smart Grid for Curtailment Selection Now we mathematically define the smart grid model. As per our model, the Smart Grid consists ofM nodes. For each node, there areN curtailment strategies available. The nodes are split into two disjoint sub-setsS: supply nodes andD: demand nodes. LetM s be the number of supply nodes andM d be the number of demand nodes. LetT be the number of time intervals in the curtailment horizon i.e., the intervals during which we schedule the curtailment. Each supply nodeb s 2S is associated with a maximum supply valueS bs (t)8t2f1;:::;Tg. Similarly, each demand nodeb d 2D is associated with a maximum load valueL b d (t)8t2f1;:::;Tg. We are given time varying curtailment matrices s (t)R MsN and d (t)2 R M d N with elements s;bj (t) and d;bj (t) denoting the discrete supply and demand curtailment obtained by node b following curtailment strategy j at time t2f1;:::;Tg respectively. For each time t, we are also given a cost matricesC s (t)2 R MsN andC d (t)2 R M d N wherec s;bj (t) andc d;bj (t) denotes the cost associated with node b following supply and demand curtailment strategy j respectively. Let X s (t) and X d (t) be the decision matrices. Elements x s;bj (t) = 1 if node b follows supply curtailment strategy j and 0 otherwise. Similarly, x d;bj (t) = 1 if node b follows demand curtailment strategy j in interval t and 0 otherwise. If the framework is being used in single control mode, for each interval t, we are given a curtailment target t calculated by taking the difference between the aggregate supply and demand. t represents the desirable curtailment target for each period, however, it might be exceeded. In order to limit wasteful curtailment, we are also given , which denotes the upper bound on the achieved curtailment in the curtailment horizon. The notations used in the following sections are summarized in Table 3.1. 25 Table 3.1: List of Variables in the Models Variable Meaning M s Number of supply nodes M d Number of demand nodes N Number of curtailment strategies T Number of time intervals in curtailment horizon s;bj (t) Curtailment achieved by supply nodeb following curtailment strat- egyj at timet d;bj (t) Curtailment achieved by demand nodeb following curtailment strat- egyj at timet c s;bj (t) Cost of supply node b following curtailment strategy j at time t. Essentially, cost associated with s;bj (t) c d;bj (t) Cost of demand node b following curtailment strategy j at time t. Essentially, cost associated with d;bj (t) x s;bj (t) 0-1 decision variable which denotes whether supply node b should follow (1) strategyj at timet or not (0) x d;bj (t) 0-1 decision variable which denotes whether demand nodeb should follow (1) strategyj at timet or not (0) t Curtailment target for intervalt Upper bound on the curtailment achieved in the curtailment horizon b B b ;B b Lower bound and upper bound on the curtailment budget for nodeb S;D supply and demand node sub-sets F Set of LV feeders T Set of transformers T (tx) Set of all nodes connected to transformertx F(f) Set of all nodes connected to feederf F tx (f) Transformer connected to feederf Cap tx Maximum solar capacity of distributor connected to transformertx Cap f Maximum capacity of power flow from transformer to sub-station via feederf S b (t) Maximum supply for a supply nodeb in time intervalt L b (t) Maximum load (demand) for a demand nodeb in time intervalt 3.4 Modeling Smart Grid Network The network model of the Smart Grid that we consider in this work is adapted from the one used in [40] and is illustrated using Figure 3.7. Power from the high voltage transmission network is fed into a sub-station, where the voltage is stepped down to be 26 delivered to the Smart Grid. Several Low V oltage (LV) Feeders originate from the sub- station and terminate in transformers which further steps down the voltage to deliver to residential consumers using distributors. Figure 3.7: Smart Grid Network Model In a smart grid with a high PV penetration, several consumers with PV installations can exhibit surplus supply which gets redistributed among the loads connected to the same transformer using the shared distributors. We also assume that any remaining surplus can be supplied back through the LV Feeders to the sub-station and eventually to the other transformers. We assume an LV Feeder capacity limit on the surplus that can be supplied back. We also assume that each distributor has an upper limit on the power that can be fed into it from the solar supply. This limit could be artificially applied by the utility to limit the generated power to ensure reliable grid operations as the current grids are not planned for bi-directional power flows. Typically, transformers have a capacity rating of the amount of power that can flow in each direction i.e. from LV feeder to distributor and vice-versa. We assume that this rating is more than the capacities of the LV feeder and the distributor and hence its modeling is redundant. 27 Note that we do not assume any limits on the amount of power that can flow from a sub-station to a transformer using the LV feeders and from a transformer to the loads using the distributors. This is a reasonable assumption as a smart grid infrastructure is typically planned to ensure uninterrupted supply to all the load consumers. The non- zero capacity constraint of LV-feeders and the solar input capacity constraint of the distributors are the additional constraints added to our network model when compared with the one used in [40]. LetF denote the set of the LV-feeders and letT denote the set of all the transformers in the Smart Grid. The capacity of a feederf isCap f and that of a distributor connected to transformertx isCap tx . The set of all the nodes connected to a transformertx2T is denoted usingT (tx). Note thatT (tx)\S andT (tx)\D denote the supply and demand nodes connected to the transformer tx respectively. The set of all the nodes fed by an LV-feederf is denoted usingF(f). This contains all the nodes connected to the transformer which is connected to the sub-station by the LV-feeder.F(f)\S and F(f)\D represent the supply and demand nodes fed by the feederf. We also use the notationF tx (f) to represent the transformer connected to a feederf. Now, the total solar supply input to the distributor connected to a transformer should be within its capacity. Hence, Sup tx (t) = X bs2T (tx)\S (S bs (t) N X j=1 s;bsj (t)x s;bsj (t))8tx2T;8t2f1;:::;Tg (3.2) Sup tx (t)Cap tx 8tx2T;8t2f1;:::;Tg (3.3) Moreover, the total surplus power injected from the transformer into the sub-station via a feeder should be within the capacity of the feeder. Therefore, 28 Sup f (t) = X bs2F(f)\S (S bs (t) N X j=1 s;bsj (t)x s;bsj (t))8f2F;8t2f1;:::;Tg (3.4) Dem f (t) = X b d 2F(f)\D (L b d (t) N X j=1 d;b d j (t)x d;b d j (t))8f2F;8t2f1;:::;Tg (3.5) Sup f (t)Dem f (t)Cap f 8f2F;8t2f1;:::;Tg (3.6) Note that our network model considers only the capacity constraints of the infras- tructure such as feeders and distributors. We do not consider the effects of the structure of the network on local grid operating state variables such as voltage fluctuations or phase imbalances. 3.5 Modeling Storage for Supply Demand Matching Algorithms Energy storage systems are being increasingly adopted in cooperation with renewable energy sources such as solar PVs to provide uninterrupted supply to the consumers [30]. Storage systems can be used as a buffer to act as either supply or load when required, thus reducing the uncertainty of renewable generation. Therefore, modeling storage systems for supply demand matching algorithms will greatly improve the scope of the algorithms. The energy storage model that we consider is developed in [66]. We briefly describe the significant parameters here: Energy Storage Capacity R (kWh): The maximum amount of energy that can be stored. This energy is with reference to a minimum storage capacity i.e. the 29 capacity below which the storage cannot be discharged. We assume the minimum storage capacity to be zero. Storage Discharge EnergyD (kWh): The maximum output/discharging power of the storage multiplied by the interval duration. The storage outputD t for interval t satisfies 0D t D. Storage Charge Energy G (kWh): The maximum input/charging power of the storage multiplied by the interval duration. The storage input G t for interval t satisfies 0G t G. Charge ( G ) and Discharge ( D ) Efficiency: Charge efficiency G 2 (0; 1) is the ratio of charged power to the input power and discharge efficiency D 2 (0; 1) is the ratio of the output power to the discharged power. Using the parameters above, storage capacity at time intervals t and t + 1 satisfy R t+1 =R t + G G t 1 D D t with the following constraints: 0R t+1 R; 0D t D and 0G t G. A summary of the notations used in the energy storage model is listed in Table 3.2. Table 3.2: List of Variables in the Storage Model Variable Meaning R Energy storage capacity R t Energy storage available at timet D Maximum storage discharge energy in a single interval D t storage discharged during intervalt G Maximum storage charge energy in a single interval G t storage charged during intervalt D Discharging efficiency G Charging efficiency The algorithms developed in this work focus on curtailment strategies. To incorpo- rate storage, if the algorithm performs supply curtailment in an interval, storage can be 30 used in charging mode to be interpreted as a supply curtailment. Similarly, if the algo- rithm performs load curtailment, storage can be used in discharging mode. Our storage model can easily be extended to more sophisticated storage scheduling in which a global view of storage capacity across the intervals is considered and where storage can be used both for charging and discharging irrespective of whether the algorithm is performing supply curtailment or load curtailment. 3.6 NP-hardness of Supply Demand Matching Problem In this section, we formally prove the NP hardness of the problem of supply demand matching with the nodes exhibiting discrete operational states. In order to prove that this problem is NP-hard, we will define a simpler version of the problem and reduce the well known knapsack problem, which is an NP-hard problem to it. Adding any additional constraints to this simpler version will only increase the complexity of the problem. The simpler version of the problem is formulated as follows: We are given a set S of node-strategy pairs, where s ij 2 S : i 2 f1;:::;Mg;j 2 f1;:::;Ng denotes the node-i-strategy-j pair, whereM is the number of node andN is the num- ber of strategies. Given a curtailment value , we need to output a S S such that P i;j:s ij 2S (s ij ), where (s ij ) denotes the curtailment obtained by s ij and P N j=1 I(s ij ) 18i2f1;:::;Mg, whereI(s ij ) = 1 ifs ij 2 S and 0 otherwise and P i;j:s ij 2S C(s ij ) is minimized, whereC(s ij ) denotes the cost of curtailment bys ij . Theorem 1. is NP-hard. Proof. A 0-1 Knapsack problem [60] is defined as follows: Given M elements, with elementi having valuev i and sized i , find a sub-set of elements the sum of whose sizes isD and the value is maximized. To reduce this problem into , for each elementi, 31 we adds i1 with (s i1 ) =d i andC(s i1 ) =v i . We set =D. Note thatj2f1g. One can easily observe that 0-1 knapsack problem has a solution if and only if has a solution. 32 Chapter 4 Optimal Curtailment Selection for Demand Response In this chapter, we develop algorithms to perform selection of curtailment strategies for each consumer to perform Demand Response. In this scenario, the supply is fixed and the demand needs to be controlled to ensure that it remains less than supply. Demand curtailment avoids costly supply ramp up thus lowering the cost of grid operations. As supply is not controlled, all the algorithms discussed in this chapter run in single control mode with only demand being controlled. 4.1 Problem of Curtailment Strategy Selection for Demand Response A Smart Grid is typically operated by a utility. The utility is responsible for provid- ing power, controlling and monitoring the SmartGrid. In this chapter, we consider the scenario in which the utility provider has a fixed power generation capacity. Typically, this capacity is sufficient to fulfill the power requirements of the customers. However, when there is a surge in the demand from the customers, the utility needs to ensure that the demand is met by either adding generators or purchasing extra power from the spot market, both of which increase expenditure. Failure to do so compromises the system reliability and leads to blackouts. 33 Power consumption profile of a consumers varies throughout the day with periods of high demand interspersed with periods of low power consumption. Certain periods of the day observe an overlap between the high demands of several customers. We denote such periods with the power requirement of the grid substantially higher than the rest of the day as peak demand period. The demand in a peak period can exceed the power generation capacity. To minimize or avoid the expenditure of purchasing extra power during peak demand periods, utilities adopt the technique of Demand Response. Customers are either incen- tivised to reduce their consumption during a Demand Response Event (DR-Event) or they are penalized by increasing the cost of power during these periods. This reduces the peak power consumption which is now expected to be met by the available genera- tion capacity. Using smart meters, utilities have the power consumption data of each customer. The granularity of the data can be as small as 15 minutes. The power consumption profile of a customer does not change rapidly from day to day, so it is straightforward to predict future pattern. By employing prediction techniques, utilities determine the peak demand periods. They also determine the targeted curtailment required for a DR-Event which should be scheduled during this period. Discussion on the prediction techniques is beyond the scope of this paper. Readers can refer to [13] for further knowledge on this topic. Utilities roll out a program to implement Demand Reponse and enroll customers into it. A customer is provided with a list of strategies to be followed each of which leads to a certain amount of curtailment in power consumption. Strategies can include procedures such as increasing the temperature of the AC systems by 2 degrees or turning off every other light in the hallways, etc. which reduce power consumption. During a DR-Event, the utility signals each customer to follow a particular strategy. A consumer 34 may be penalized if it fails to comply. For instance, the University of Southern California SmartGrid consists of 50,000 sensors across the 170 buildings to monitor electricity usage. Each building can adopt any one of seven available strategies during DR events which occur on Weekdays 1-5 pm [3]. Careful selection of consumers is required to ensure that the targeted curtailment value is met. A good selection algorithm determines the subset of consumers along with the strategies they should follow during the DR-Event such that the achieved curtailment value is as close as possible to the target. The reasons are as follows: 1. Limiting the amount by which the achieved curtailment value overshoots the tar- get ensures that the grid is not underutilized. This avoids any loss of revenues to the utility due to underutilization of grid by aggressive curtailment. 2. Limiting the amount by which the achieved curtailment value undershoots the target ensures that the utility can avoid purchasing power from external sources by bounding the peak demand of the customers. In the following sections, we develop several formulations addressing various use cases to optimally perform Demand Response. 4.2 Traditional Demand Response We formally define the problem of optimal customer selection for Demand Response (Traditional Demand Response) using the parameters defined in Table 3.1. We are given a list ofM d consumers andN strategies. Each consumer can adopt exactly one strategy in the DR event. The decision variablex bj is 1 if consumerb adopts strategyj. We are also given the curtailment in power consumption d;bj obtained by consumerb adopting 35 strategyj. A default strategy with a curtailment value of 0 is also included in the cur- tailment matrix d . A customer adopting a default strategy essentially means that it is not participating in the DR event. A targeted curtailment value for the DR event is provided. The objective it to achieve a curtailment value as close to as possible. The ILP formulation for this problem is as follows: Minimize : j M d X b=1 N X j=1 bj x d;bj j (4.1) Subjectto : N X j=1 x d;bj = 1; b2f1;:::;Mg (4.2) x d;bj 2f0; 1g; 8b;j (4.3) Equation 4.1 minimizes the absolute curtailment error. Equation 4.2 ensures that a consumer cannot adopt more than one strategy in the DR event. Detailed experimental results for consumer selection using the above ILP is shown in Section 4.7. 4.3 Sustainable Demand Response 4.3.1 Motivation The algorithm mentioned in the previous section might aggressively curtail the demand in some intervals while accumulating demands in other intervals. Such assignments have peaks in certain intervals, which can possibly exceed the generation capacity forcing the utility to pay for additional procurement of energy. We define the notion of Sustainable Demand Response (SDR) to address such cases. SDR attempts to evenly smooth the curtailment over the entire period of the DR event. 36 Hence we define SDR as the customer-strategy assignment which minimizes theklk 1 distance between achieved curtailment values and a smoothed target value per interval. As before we are given a setS ofM d consumers,N strategies. The entire DR period is divided into discrete time intervals. Dynamic customer strategies are represented by a time varying curtailment matrix d (t)2 R M d N with element d;bj (t) denoting the discrete curtailment value of consumer b adopting strategy j at time interval t where t2f1;:::;Tg. Let X d (t) be the decision matrix with element x d;bj (t) denoting the corresponding decision variable at time t with denoting the achievable curtailment value across the entire DR event. 4.3.2 ILP Formulation for Sustainable DR We use the following ILP to model a Sustainable DR event. Minimize : T X t=1 t (4.4) Subjectto :j M X b=1 N X j=1 d;bj (t)x d;bj (t) T j t 8t (4.5) N X j=1 x d;bj (t) = 1 8b;t (4.6) 8x d;bj (t)2f0; 1g 8b;j;t The objective is the minimize thekLk 1 norm (Equation 4.4). As before, Equa- tion 4.6 ensures that at any given interval, each consumer adopts exactly one strategy. The various intervals in the ILP above are independent. This makes it trivial to paral- lelize by solving each interval as a separate optimization problem on a single node in a high performance cluster. 37 4.4 Sustainable DR with Strategy Overheads 4.4.1 Motivation Driven by our experience with existing DR implementations on the USC smartgrid, we observe that it is impractical for customers to switch between too many strategies dur- ing the DR event as this leads to additional overhead costs. We model this by associated a cost of 1 for each strategy transition and by limiting the number of transitions for each node by. An additional constraint is added to the ILP to incorporate the strategy switching overheads. Note that under this formulation, a consumer is likely to have con- tiguous strategies across intervals. In our experimental work section, we solve this ILP exactly for reasonable problem sizes (representing the USC microgrid) using the IBM CPLEX Solver. For very large problem sizes, the time required for an exact solution might be large. In such cases, one can use randomized rounding heuristics based on the LP-relaxation of the ILP to obtain approximate solution. 38 4.4.2 ILP formulation for Sustainable DR with Strategy Overheads We use the following ILP to model a Sustainable DR event with strategy overheads. Minimize : T X t=1 t (4.7) Subjectto :j M d X b=1 N X j=1 d;bj (t)x d;bj (t) T j t 8t (4.8) N X j=1 x d;bj (t) = 1 8b;t (4.9) x d;bj (t)2f0; 1g 8b;j S bj (t) =jx d;bj (t)x d;bj (t 1)j 8b;j;t2f2;:::;Tg (4.10) T X t=2 N X j=1 S bj (t) 2 8b (4.11) The new constraints to limit the strategy switching are introduced using 4.10 and 4.11 where 4.10 calculates the number of times customer b switches a particular strategy. Equation 4.11 bounds the total number of times a customer can switch strate- gies. Since the state variableS bj (t) counts both switching into and switching out from strategyj, equation 4.11 uses 2 as the bound. In our experiments, we fix the value of = 2. 39 4.5 NO-LESS: Near OptimaL CurtailmEnt Strategy Selection for Supply Demand Matching in Micro Grids 4.5.1 Problem Definition We are given a set of M d nodes and N strategies. The entire curtailment horizon is divided into discrete time intervals. We are given a time varying curtailment matrix d (t)2R M d N with element d;bj (t) denoting the discrete curtailment value of nodeb adopting strategy j at time interval t where t2f1;:::;Tg. Let B b be the maximum curtailment value for node b. Let X d (t) be the decision matrix with element x d;bj (t) denoting the corresponding decision variable at timet. The achievable curtailment value across the entire curtailment horizon is given by . b denotes the limit on the total cost of strategy switches for a nodeb. We assume that strategy 1 is the default strategy with a curtailment value of 0 i.e. if a node is not included in an interval of the curtailment horizon, it follows strategy 1. Given the model above, the objective is to determine node-strategy pairs for each interval such that: (1) The curtailment target is achieved with minimum curtailment error: difference between the achieved and the targeted curtailment, (2) No node is curtailed more than its maximum allowable value, (3) only allowable strategy switches are performed by each node across each consecutive intervals, and (4) for each node, the cost incurred in switching strategies is within its allowable limit. 40 4.5.2 ILP Formulation for NO-LESS NO-LESS can be formulated using the following ILP: Minimize (4.12) s:t: T X t=1 M X b=1 N X j=1 d;bj (t)x d;bj (t) (4.13) N X j=1 x d;bj (t) = 1 8b; t = 1; 2;:::T (4.14) T X t=1 N X j=1 d;bj (t)x d;bj (t)B b 8b (4.15) S bij (t) =jx d;bi (t 1)x d;bj (t)j 8b; 8t; i;j2f1;:::;Ng (4.16) x d;b1 (0) = 1 8b (4.17) x d;bj (0) = 0 8b; j2f2;:::;Ng (4.18) T X t=2 N X i=1 N X j=1 (x d;bi (t 1) +x d;bj (t)S bij (t)) b (i;j) 2 b ; 8b (4.19) 0 (4.20) x d;bj (t)2f0; 1g8b; 8j; t = 1; 2;:::T (4.21) Equations 4.12, 4.13 ensure that the targeted curtailment is met with minimum error. Equation 4.14 ensures that a node follows only 1 strategy in any interval. Equa- tions 4.15- 4.19 ensure that the curtailment and switching cost for each node is bounded. The ILP above returns the optimal solution, however due to its large execution time we develop an FPTAS for it. 41 4.6 Approximation Algorithms In the previous sections, we formulated problems for various use cases using Integer Linear Programs. ILPs return optimal solution, however, they require large execution time and are not scalable. Hence, in this section we develop approximation algorithms for the ILP algorithms defined above. 4.6.1 Fast p 2-factor Approximation for Sustainable DR We now describe a fast algorithm for computing approximately optimal sustainable DR strategies. Our algorithm provides a p 2-factor approximation to the optimal target dur- ing each curtailment period and therefore for the entire DR event. Theorem 2. Algorithm 1 is a p 2-factor approximation to the optimal sustainable DR solution. Proof. Let Y t denote the curtailment value for period t returned by Algorithm 1. We show that eitherY t 2 [ T p 2 ; p 2 T ] orY t is the optimal curtailment value achievable for that period. If line 2 of the algorithm is satisfied, then this is trivially true. Assume line 2 is not satisfied. From line 8, we must have =(T p 2)Y r = P r1 b=1 d;bk b + d;rkr 2 =(T p 2) = p 2=T . Finally, consider the case when whenr = M andY M =( p 2T ). Sincek b rep- resents the strategy with the largest curtailment value =( p 2T ) for each customerb, by definitionY r is the largest achievable curtailment value< Y p and so eitherY p orY r is the optimal strategy for this period. The following result follows from a straightforward analysis of the algorithm. Theorem 3. Algorithm 1 can be used to compute p 2-approximate sustainable DR solu- tions in O(TM d logN) time when strategies are preprocessed in advance for a given 42 Algorithm 1: Fast p 2-factor Sustainable DR Approximation Preprocessing: Non-decreasing sorted listsf d;b (t)g of consumer-strategies for each consumerb2S and each intervalt 1 for intervalst = 1 toT do 2 if9(k2S)^ (j2 d;k (t)) : d;kj (t)2 [ T p 2 ; p 2 T ] then 3 x d;kj (t) 1 ; // Select consumer k and curtailment strategy d;kj (t) 4 else 5 (p;q p ) argmin b2S;j2 d;b (t) f d;bj (t)j d;bj (t) p 2=Tg ; // p and q p represent the customer and strategy indices of the consumer with the smallest curtailment value p 2=T 6 Y p d;pqp (t) ; 7 For each customerb2S :k b argmax k f d;bk (t) =(T p 2)g ; 8 Letr denote the smallest index such that Y r P r i=1 d;bk b (t);Y r =(T p 2) ; 9 r M if P M i=1 d;bk b (t)< =(T p 2) ; // Select Strategies for Activation as follows 10 ifY p p 2=T =(T p 2)Y r then 11 Setx d;pqp (t) = 1 ; 12 else 13 forb 1 tor do 14 Setx d;bk b (t) = 1 ; 15 end 16 end 17 end 18 end Output: MatrixfX d (t)g of selected consumer-strategies8t. Bounded curtailment valuesb t 2 [ ~ t p 2 ; p 2~ t ], where ~ t = max( t ; =T ); t is the optimal solution to ILP. curtailment target. The one-time preprocessing cost assuming apriori knowledge of curtailment strategies isO(TM d N logN). By precomputing and storing results for a range of target values, we can speed up the retrieval of approximately optimal strategies even further toO(1) time. 43 Note that the same algorithm can be used as an approximation algorithm for Tradi- tional Demand Response. 4.6.2 A PTAS for Sustainable DR While the approximation algorithm above can be used to very quickly compute sus- tainable DR solutions, the error due to the p 2-factor approximation may be unaccept- ably large in some cases. Therefore, using ideas from the subset sum problem [35] we develop a Polynomial Time Approximation Scheme (PTAS) that approximates the opti- mal solution provided by the ILP in Equation 4.4 to within an arbitrarily small-factor in time polynomial in M d N=. Since each customer is restricted to using exactly one strategy, we cannot simply merge all the customer strategies and find the subset that comes closest to the target =T . Theorem 4. Algorithm 2 is a PTAS for the ILP in Equation 4.4. Proof Sketch: The number of intervals l log 1 1 (=T ) is polynomial in ln(=T ) . In line 9, if a curtailment value is already marked feasible (i.e Q (k1) s (t) is 0), then customerk does not contribute to the feasible solution. The total number of iterations isO(lNM d ). From line 7, using induction, we can showV j+1 P d;qr 2B (M) j (t) d;qr (1) M V j+1 (1)V j+1 =V j . and hence P d;qr 2B (M) j (t) d;qr 2 [V j ;V j+1 ]. Thus the algorithm outputs an approximation to the optimal achievable target and is therefore a PTAS. . 4.6.3 FPTAS for NO-LESS We first develop a dynamic programming based FPTAS to determine the set of curtail- ment strategies followed by a single nodeb during the curtailment horizon to achieve a curtailment value of at most (1 +), where is a user determined accuracy parameter. 44 Algorithm 2: PTAS for Sustainable DR during each intervalt Preliminaries: V 0 min i;j d;ij (t);Z min i;j d;ij (t) =T ; Divide [V 0 ;Z] intol intervalsf[V i ;V i+1 ]g,V i+1 = (1 +)V i , 0il 2, V l =Z ; Initialization: 8k2S :Y kj V i if d;kj 2 [V i ;V i+1 ] ; B (0) i (t) ;;Q (0) i (t) 0;i = 0; 1;:::;n 1 ; // B (k) i (t) is a subset of the first k consumers, each with a non-zero strategy selection, that add up to a total curtailment value 2 [V i ;V i+1 ]. 1 for Consumersk = 1 toM d do 2 for all intervalsi, all strategiesr do 3 Z ir Q (k1) i (t) +Y kr ; 4 LetZ ir 2 [V s ;V s+1 ] ; 5 if:Q (k1) s (t) then 6 Q (k) s (t) V s // Making Q (k) s (t) a feasible curtailment value 7 B (k) s (t) B (k1) i (t) S d;kr (t) ; // Adding consumer k strategy r pair to the feasible strategy set 8 end 9 end 10 end Output:B (M d ) j (t): Selection of consumer-strategy pairs, wherej is the closest interval to DR target =T withQ (M) j (t)> 0 The cost of strategy switches will be bounded by b . We then combine the results of all the nodes to achieve the targeted curtailment value from all the nodes. Achieving Curtailment Target for a Single Node Let = T . For each d;bj (t), we defineb d;bj (t) =b d;bj (t) c. We also define b =b c. We define a boolean function which returns true if a curtailment value b t can be 45 achieved using strategyS k at timet incurring a costq t for strategy switches and False otherwise. can be defined using the following dynamic programming formulation: (b t ;t;S k ;q t ) = FALSE ifb t b d;bk (t)< 0jjq t < 0 =jj j (b t b d;bk (t);t 1;S j ;q t b (j;k)) (8j2f1;:::;Ng) otherwise (4.22) If any of ( b ;T;S k ; b )8k2f1;:::;Ng is True, we can determine the required strategies by traversing the recursive formulation and determining the strategies at each time t as described in Algorithm 3. The dynamic program can be solved by creating a table of size b TN b . We will refer to the entire table as to simplify the notations. We can initialize the table using the following equations: ( l ; 1;S k ;q) =TRUE ifb d;bk (1) == l 8k2f1;:::;Ng &&q>= b (1k) 8l FALSE otherwise (4.23) Lemma 1. Algorithm 3 finds the strategies to be followed in each interval by node b to achieve the curtailment value within a cost of b for strategy switches with a time complexity which is polynomial inN;T and 1 . Proof. We omit the proof of correctness as it can easily by argued using the arguments used for dynamic programming algorithms. Now, filling a single entry of the table requiresO(N) time. Hence, the algorithm requiresO( T 2 N 2 b ) which is the dominating term. We assume b = O(T ) i.e. the ratio of the largest to smallest cost is bounded and the number of switches is T , the algorithm is polynomial in N;T and 1 . The total number of entries in the table are b TN b . N b TN b = 46 Algorithm 3: A (1 +) polynomial time approximation algorithm to achieve the curtailment value by a single nodeb 1 Fill entries (b l ;t;S k ;q t )8b l 2f0;:::; b g;8t2f1;:::;Tg;8S k ;k2 f1;:::;Ng;8q t 2f1;:::; b g using equations 4.22 and 4.23 2 X d (t) 8t2f1;:::;Tg 3 cur b 4 q cur b 5 for Timet =T to 1 do 6 if9k2 1;:::;N s.t ( cur ;t;S k ;q cur ) then 7 if !( cur c bk (t);t 1;S k ;q cur ) then 8 q cur q cur 1 9 end 10 cur cur d;bk (t) 11 X d (t) S k 12 end 13 end Output:X d , whereX d (t) denotes the strategy to be followed at timet N T TN b = O( T 2 N 2 b ). If we assume that b = O(T ) the algorithm is polynomial in N;T and 1 . (For loop requires only O(TN) amount of time after the table is filled.) Lemma 2. The total curtailment value obtained by following the strategies output by Algorithm 3 is (1+) where is the total curtailment value obtained by following the strategies output by an optimal algorithm. Proof. Let x be the curtailment value obtained by following the strategies output by Algorithm 3. Let be the optimal curtailment value. Clearly, x P T t=1 d;bX(t) (t). Also, P T t=1 ( d;bX(t) (t) 1). This implies x +T +. Now, since , x (1 +). Using Lemmas 1 and 2, we get the following theorem. 47 Theorem 5. Algorithm 3 is a FPTAS to determine the strategy to be followed in each interval 1;:::;T to achieve a curtailment value by a single node b with the cost of strategy switches is bounded by b . Achieving Total Curtailment by all the Nodes We will run step 1 of Algorithm 3 with curtailment value b B b , the node-specific maximum curtailment value, 0 = M d , and limit on the cost on strategy switch- ing b for each node b 2 f1;:::;M d g to create b . For a node b, define b = [ b l j 9fS k ;qg with b (b l ;T;S k ;q) = TRUE. Assume b is in sorted order. Let L( b ) denote the number of entries in b We will again use a dynamic programming based algorithm for this problem. We define a boolean function which returns True if a curtailment valueb b can be obtained by using nodesf1;:::;bg and False otherwise. can be defined using the following formulation: (b b ;b) = TRUE if9k2f1;:::;L( b )g s:t:(b b b (j);b 1) is TRUE = FALSE otherwise (4.24) can be initialized for all b b 2 f0;:::; b g, where b = b c; = M d T by the following equations: (b b ; 1) = TRUE if9k2f1;:::;L( b )g s:t:b b = b (k) = FALSE otherwise (4.25) 48 Given a total curtailment value , Algorithm 4 can be used to combine the curtailment values that can be achieved from individual nodes to attain the total curtailment value. Algorithm 4: A (1 +) polynomial time approximation algorithm to achieve the curtailment value 1 Run Step 1 of Algorithm 3 for each nodeb with b B b =b b B b c, where = M d T , 0 = M d , and = b and produce b 2 Fill entries (b b ;b)8b b 2f0;:::; b g;8b2f1;:::;M d g using equations 4.24 and 4.25 3 X b (t) 8b 2f1;:::;M d g;8t2f1;:::;Tg 4 cur b 5 forb =M to 2 do 6 if9k2f1;:::;L( b )g s.t. ( cur b (k);b 1) == 1 then 7 Call Algorithm 3 with inputs b (k);T; b to produce outputX 8 X b X 9 cur cur b (k) 10 end 11 end Output:X, whereX b (t) denotes the strategy to be followed by nodeb at timet Theorem 6. Algorithm 4 is a FPTAS for NO-LESS Proof. The correctness can be argued similar to that of Algorithm 3. Let x be the curtailment value obtained by following the strategies output by Algo- rithm 4. Let be the optimal curtailment value. Clearly, x P M b=1 P T t=1 d;bX(t) (t). Also, P M b=1 P T t=1 ( d;bX(t) (t) 1). This implies x +M d T +. Now, since , x (1 +). Step 1 of Algorithm 4 requiresO( M d T 3 N 2 ) time. Step 2 requiresO( b ) time to fill each of the b M d entries hence,O( M 3 d T 2 2 ) time. These being the dominating terms, the total complexity of Algorithm 4 isO( M 3 d T 2 2 + M d T 3 N 2 ) which is polynomial inM d ;N;T and 1 . The For loop containsM d calls to Algorithm 3 and hence requiresO(M d T 2 N 2 ) time. Remark: In a typical micro-grid such as the one used for our experiments, typically 49 the number of strategies isN 10 and the number of time intervals of the DRT 16 15-minute intervals i.e. 4 hours. However, the number of nodes is typically large and keeps on increasing with new nodes being added to DR. Hence, our algorithm with a complexity ofO( M 3 d 2 ) and an optimality guarantee of getting a solution within (1 +) times the optimal for this NP-hard problem is a significant improvement over the state- of-the-art techniques. 50 4.7 Experimental Results 4.7.1 Experimental Setup for TDR and SDR The USC SmartGrid has over 50,000 sensors to monitor electricity usage and equipment status in real time [3]. Demand Response Events occur on weekdays between 1 and 5 pm. We use the data collected by the software developed to support data-driven demand response optimization in USC smartgrid [62]. The software provides us the curtailment values for each strategy that can be adopted by any building (customer) in USC for the queried time interval. For our experiments we use the data from 27 buildings each of which can adopt one of seven strategies. The data is collected for the time intervals 1-3 pm and 3-5 pm for each day from Monday to Friday. We run our experiments for Targeted Curtailment values ranging from 100 kWh to 1400 kWh. We use the Optimization Programming Language [2] to define the Integer Linear Programming formulation developed in this paper. IBM ILOG CPLEX optimization software [1] is used to solve the ILP and produce the set of customers and the strategies they should adopt. We compare our results with the state-of-the-art heuristic [76]. Authors in [76] develop a change making problem based algorithm for customer selection. The change making problem determines how to make a given amount of money using the least amount of coins. The coins in the algorithm are the available customer-strategy pairs and their values the predicted curtailment values. The amount to be made is the tar- geted curtailment value. Customers are grouped into bins differentiated by their values. A greedy algorithm is used to pick customers from the bins with highest values. We choose this heuristic for comparison as it is used in practice by the USC SmartGrid to schedule customers and their strategies for the DR events. 51 4.7.2 Results and Analysis for TDR and SDR The power consumption profile of a building is similar for the same day of a week across different weeks. So by running our experiments on data collected from DR events over a week, we are able to demonstrate our algorithm on a wide range of power consumption profiles. Moreover, a typical DR-event in the USC SmartGrid starts with the selection procedure at 1 pm and then another selection occurs typically around 3 pm. Therefore, we consider them as two separate DR events for our experiments. Figure 4.1 and 4.2 show the errors incurred by our ILP based customer selection algorithm and the State-of-the-art heuristic [76] for every DR event from Monday to Friday for various curtailment target values. As shown in Figure 4.1, the highest error incurred by our ILP based algorithm is around 0.0002 kWh during the DR-events on Tuesday 1-3 pm for a targeted curtailment of 600 kWh and Friday 3-5 pm for a targeted curtailment of 800 kWh. For the State-of-the-art heuristic, the highest error incurred is around 8 kWh during the DR event on Tuesday 3-5 pm for the targeted curtailment of 400 kWh as shown in Figure 4.2. One may note that the errors between the state-of-the-art heuristic and our approach differ by multiple orders of magnitude. so we take the ratio of the error for comparison. A higher value of ratio implies better performance by our approach. In Figures 4.3-4.12 we compare the errors incurred by the two algorithms for various targeted curtailment value for each DR event The customer selection problem can be visualized as a packing problem. We are trying to pack the targeted curtailment with values obtained from the customer-strategy pairs. The ILP produces the best possible packing. Any error is due to the nature of the data. Similarly, the heuristic based approach tries to provide best packing in each of the bins. Error incurred in packing each bin accumulates throughout the algorithm and may lead to very large errors. Since we are using real world data, as seen in the 52 Figures 4.3-4.12 the peaks in the ratio of errors for various DR events varies with the targeted curtailment values with no discernible pattern. The highest ratio observed is around 3 10 7 which occurs during the DR event on Thursday 1-3 pm. Although solving an ILP is computationally intensive, optimal customer selection for each target was obtained in less than 5 seconds on a standard workstation. This time can be significantly reduced by using sophisticated computational platforms. Note that in a typical DR Event, the utility determines the curtailment target well in advance. Thus even a 5 second delay in computing the optimal customer-strategy pairs and signaling the customers is tolerable. 4.7.3 Sustainable DR In Figures 4.13-4.17, we compare the absolute errors in kWh incurred by the Sustainable Demand Response (SDR) ILP formulation and the state-of-the-art heuristic for Targeted curtailment values ranging from 50-1000 kWh. The vertical axis is limited to a maxi- mum error of 3.0 kWh to ensure the readability of the graph. The actual error values are given in the table below the graphs. Our SDR ILP is more restrictive as it ensures that the curtailment is distributed evenly across the intervals which is not a constraint in the state-of-the-art heuristic. Despite this restriction, our SDR ILP perform 4-2000 times better than the heuristic. The only time the heuristic performs better than our ILP is on Wednesday for a Target of 1000 kWh where the error of our SDR ILP is 0.048 kWh while that of the heuristic is 0.001 kWh. To emphasize the significance of Sustainable DR over a DR targeting the entire interval, we compare the demand curtailment values achieved in each interval for the DR event for the targeted curtailment value of 1200 kWh. The DR ILP targeting the entire interval (TDR) incurs an error of 0.002 kWh which is far lower that the 0.031 kWh 53 Figure 4.1: Error incurred by ILP based Algorithm for various targeted Curtailment Values for Different DR Event error incurred by the Sustainable DR ILP. However, most of the curtailment is achieved in intervals 10 and 11 and the rest of the intervals have lower curtailment values. The Sustainable DR achieves a curtailment value of around 75 kWh in each interval. 4.7.4 Sustainable DR with Strategy Overheads The motivation behind using Sustainable DR with strategy overheads is to reduce the number of times each building switches its strategy in the DR interval. We limit the number of times building i can switch strategies to = 2. In Figures 4.13-4.17, we 54 Figure 4.2: Error incurred by State-of-the-art Algorithm for various targeted Curtailment Values for Different DR Event compare the absolute errors in kWh incurred by the Sustainable Demand Response with limits on strategy switching (SDR with switch limit) ILP formulation and the state- of-the-art heuristic for Targeted curtailment values ranging from 50-1000 kWh. The vertical axis is limited to a maximum error of 3.0 kWh to ensure the readability of the graph. The actual error values are given in the table below the graphs. The SDR with strategy overheads is even more restrictive than the SDR ILP. This is reflected in the errors incurred which are an order higher than the SDR ILP. However, SDR with strategy overheads still performs 2-700 times better than the state-of-the-art 55 Figure 4.3: Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Monday 1-3 pm heuristic. The only times the heuristic performs better than our ILP is on Wednesday for a Target of 1000 kWh where the errors 0.083 kWh and 0.001 kWh respectively and on Tuesday for a Target of 400 kWh where the errors are 0.203 kWh and 0.172 kWh respectively for our ILP and the heuristic. The maximum relative error incurred by our ILPs is 1% whereas for the heuristic we observed relative errors as high as 12 %. ILP is a computationally intensive process. To converge to the error rates shown in the results, the IBM CPLEX required 5-10 minutes of processing time. Since DR programs are based on predictive data [13] which are available in advance, the time required can be perceived as reasonable. 56 Figure 4.4: Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Monday 3-5 pm Figure 4.5: Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Tuesday 1-3 pm 57 Figure 4.6: Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Tuesday 3-5 pm Figure 4.7: Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Wednesday 1-3 pm 58 Figure 4.8: Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Wednesday 3-5 pm Figure 4.9: Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Thursday 1-3 pm 59 Figure 4.10: Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Thursday 3-5 pm Figure 4.11: Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Friday 1-3 pm 60 Figure 4.12: Ratio of error of State-of-the-art Heuristic and ILP algorithm for DR event Friday 3-5 pm Figure 4.13: Absolute errors in kWh incurred for Sustainable DR event on Monday 1-5 pm 61 Figure 4.14: Absolute errors in kWh incurred for Sustainable DR event on Tuesday 1-5 pm Figure 4.15: Absolute errors in kWh incurred for Sustainable DR event on Wednesday 1-5 pm 62 Figure 4.16: Absolute errors in kWh incurred for Sustainable DR event on Thursday 1-5 pm Figure 4.17: Absolute errors in kWh incurred for Sustainable DR event on Friday 1-5 pm 63 Figure 4.18: Comparison of demand curtailment values achieved in each interval for SDR and TDR for the targeted curtailment value of 1200 kWh 4.7.5 NO-LESS Experimental Setup for NO-LESS We evaluate NO-LESS using the curtailment dataset obtained by the DR programs con- ducted in the University of Southern California (USC). The USC micro grid consists of more than 150 buildings equipped with smart meters to capture consumption data in 15-minute intervals. Each DR lasted for 4 hours during which each building switched to one out of 6 available strategies [36]. The curtailment data was obtained by using the algorithm developed in [22]. Althoug NO-LESS is developed for load curtailment, we 64 also evaluated its performance for strategy selection for solar curtailment. For solar cur- tailment, we use simulated data generated by using solar irradiance data for Los Angeles County [5] and varying the solar PV are from 10m 2 to 20m 2 and the yield from 5% to 15%. 6 different curtailment strategies: 0, 0:125O; 0:25O; 0:5O; 0:75O;O, with O being the maximum output were generated. The cost of switching between strategies was fixed as 1. And the cost of strategy switching was limited to 10. The allowable strategy switches were 0$ 1$ 5, 0$ 2$ 4$ 5 and 0$ 3$ 5, where 0 is the default strategy with a curtailment value of 0. We implemented the NO-LESS in Java. The experiments were performed on Dell optiplex with 4-cores and 4 GB RAM. 100 200 500 1000 1500 0.00% 1.00% 2.00% 3.00% 4.00% 5.00% 6.00% 7.00% 8.00% є=0.1 (10%) є=0.2 (20%) є=0.02 (2%) є=0.01 (1%) є=0.05 (5%) Targeted Curtailment (kWh) Relative % Error Figure 4.19: Near Optimality of NO-LESS Near Optimality We implemented an Integer Linear Program (ILP) for the problem objective to obtain optimal solutions. We varied the targeted curtailment values from 100 kWh to 2000 kWh. The number of nodes was fixed at 20. Figure 4.19 shows the relative percentage 65 error incurred by NO-LESS algorithm compared against the optimal solution. Relative percentage error is defined as: (PO)100 O , whereP is the solution obtained by NO-LESS andO is the optimal solution obtained by the ILP. As shown in the figure, the relative percentage errors incurred by NO-LESS are much less than the worst case guarantee as determined by. Notice that the relative percentage error increases as the curtailment target increases. This is because our algorithm is essentially a packing problem. For lower curtailment values, the algorithm has enough available node-strategy pairs to choose from to pack. However, as the curtailment value increases, the availability of such pairs becomes sparse. While the optimal algorithm tries all possible combination for optimizing the packing, our algorithm terminates as soon as it finds one which is within the worst case guarantee. Scalability NO-LESS algorithm has two critical parameters which affect scalability: the number of nodesM and the accuracy parameter. Note that our algorithm is independent of the targeted curtailment value. 10 20 30 40 50 60 70 80 90 100 0.000 50.000 100.000 150.000 200.000 250.000 є=0.1 (10%) є=0.2 (20%) є=0.05 (5%) Number of Nodes Time (seconds) Figure 4.20: Execution Time vs Number of nodes for different values of epsilon 66 Figure 4.20 shows the effect on the execution time of the number of nodes for various values of. For each, the execution time increases polynomially with respect to the number of nodes. Even for an accuracy parameter of 5%, the runtime is less than 4 minutes. 0 0.2 0.4 0.6 0.8 1 1.2 0.000 20.000 40.000 60.000 80.000 100.000 120.000 20 60 80 40 10 Accuracy Parameter (є) Time (seconds) Figure 4.21: Execution Time vs epsilon for different number of nodes 10 20 30 40 50 0.000 100.000 200.000 300.000 400.000 500.000 600.000 700.000 e=0.01 (1%) Number of Buildings Time (seconds) Figure 4.22: Execution Time vs Number of Buildings for = 0:01 Figure 4.21 shows the trade-off between the accuracy parameter and the execution time in seconds for various number of nodes. The execution time increases quadratically with respect to 1 . For very small values of such as 0.01 (1%), this leads to very large 67 runtimes such as 600 seconds for 50 nodes. However, our algorithm fills the dynamic programming table in such a way that all entries in a row can be independently populated by looking at the entries of previous columns. This provides us with an opportunity to easily parallelize the algorithm to achieve dramatic improvement in the runtime. We do not perform such analysis in this work. Comparison with State-of-the-art Techniques We also compared our algorithm against a change making scheduler [76] based heuristic and an LP based algorithm [34] whose results are rounded to the nearest integers. The errors for the LP based algorithm range from 20% to 95% while that for the change maker heuristic range from 2% to 10%. Theoretically, both these heuristics can have unbounded errors. In contrast, the errors incurred by our algorithm are less than 1.5% for = 0:05 and less than 0.25% for = 0:01. Table 4.1 compares the theoretical and experimental errors and the computation time of our algorithm for various values of with state-of-the-art techniques. Our algorithm, with a small increase in runtime is able to provide solutions with extremely low errors as shown in the table. Table 4.1: Error and Scalability (20 nodes) Comparison Technique Theoretical Error Experimental Error Time Change Making [76] None 2-10% 1 second LP based algorithm [34] None 20-95% 1 second NO-LESS ( = 0:05) 5% <1.5% 2.2 seconds NO-LESS ( = 0:01) 1% <0.25% 33 seconds 68 4.8 Automated Dynamic Demand Response Implemen- tation on USC Campus We implemented D 2 R on USC campus micro-grid to demonstrate its large scale feasi- bility and identify and resolve the challenges associated with practical deployment. Our D 2 R technique uses learning of occupant energy strategy preferences (at fine grained scales ranging from buildings to floor levels within buildings) to make accurate elec- tricity consumption predictions and individual curtailment recommendations using the algorithms developed in this chapter using only a small subset of consumption data. 4.8.1 D 2 R Implementation Buildings on our campus are instrumented with smart metering and control that can be used to implement a large number of advanced energy curtailment strategies such as resetting temperature set points, reducing air flow, duty cycling. A simplified control and data flow diagram of our D 2 R implementation is shown in Figure 4.23. The micro- grid utility initiates a DR event using OpenADR messages and provides a curtailment target to be achieved over a given DR intervalT , typically 4 hours. The Policy Engine (PE) module provides campus wide curtailment strategy policy recommendations based on the analysis of historical consumption data and curtailment maximization customer (building) selection algorithms. Smart meter data is aggregated over 15 minute intervals into a consumption database. State-of-the-art data-driven models are then used by the PE module to predict energy consumption values over each 15 minute period comprising the entire DR interval T for each building across campus. This information is then provided to the optimization module in the form of a discrete time varying curtailment matrix. The outputs of the PE are sets of buildings-strategy pairs at 15 minute intervals that will achieve the required curtailment target overT . 69 Figure 4.23: Control and Data Flow For D 2 R Implementation Our Policy Engine (PE) is composed of an influence model based demand predic- tion engine that feeds into a building-strategy selection optimization module as shown in Figure 4.24. Historical data from the energy consumption database before the DR day is used in conjuction with time series forecasting techniques such as ARIMA and Lasso-Granger (for learning temporal dependencies among multiple timeseries) to learn occupant energy strategy preferences and load profiles. A significant challenge is to ensure accurate prediction in the absence of high quality consumption data. To ensure quality data we have developed several techniques: interpolation methods for estimating intermittent missing data and sophisticated influence based learning models for estimat- ing systematically unavailable data over larger periods. We have shown that only a small subset (7%) of the meters are required in real-time to make predictions for buildings across the campus micro-grid [13]. The results of our adaptive customer prediction models are fed into our optimiza- tion module which consists of fast approximation algorithms for ILP based fair energy 70 Figure 4.24: Policy Engine curtailment maximization [37]. Given a targeted curtailment value over DR interval T , curtailment matrixM, the optimization module selects sets of building-strategy pairs for each 15 minute interval to achieve the curtailment value overT . Note that optimiz- ing for over intervalT without considering per-period curtailment values might lead to aggressive curtailment in some periods and low curtailment in others. This could be unsustainable to the utility as demands in periods with low curtailment may exceed the generation capacity. Therefore, we have defined and implemented the notion of Sustain- able DR (SDR) in which the targeted curtailment value is distributed proportionally (as f t g) across intervalT [37]. Additionally, we have included the notion of fairness (via 71 building curtailment budgets) and strategy switching overhead as constraints in our opti- mization algorithm. We have shown that this problem is NP-hard and have developed fast polynomial time approximation schemes (PTAS) as well as bounded randomized rounding heuristics with provable error bounds i.e. deviation from the curtailment target . We have shown that our SDR algorithms achieve results with a very low absolute error of 0.001-0.05 kWh range [37]. 4.8.2 Implementation Challenges Real time quality data availability is a key challenge that we have addressed earlier (Sec- tion 4.8.1). Real world challenges such as changing environmental conditions may affect our prediction accuracy, buildings may have to be dropped out of DR due to unrespon- siveness or thermal comfort violation. We address such cases by performing dynamic customer re-selection to offset deviations in the curtailment target. Another challenge is scalability while maintaining curtailment accuracy which we have addressed by devel- oping fast polynomial time approximation algorithms with bounded errors. Finally, extending our implementation to a large city environment will face human behavioral challenges. Determining the right customer incentives to obtain a reasonable compli- ance rate is an open challenge that needs to be addressed. 4.8.3 Overall DR Evaluation In this section, we discuss the performance of a DR event performed on the campus. 28 buildings were available for the 4 hour long event from 1 pm to 5 pm. Each building had anywhere between 1 to 7 strategies available to it. A targeted curtailment of 6400 kWh was set. the Policy Engine provided building-strategy pair recommendations to the FMS using the predicted curtailment values. It estimated that a curtailment of 2332 kWh can be obtained using the available building-strategy pairs. 72 Practical challenges such as increase in temperature above the comfort zone, res- ident complaints etc. lead to dropping out of certain buildings during the DR event. Sometimes buildings fail to respond to the signals for adopting the suggested strategies. Hence, the actual obtained curtailment value was 2103 kWh for the entire DR event as opposed to the estimated 2332 kWh. ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 25 50 75 01:00 03:00 05:00 07:00 09:00 11:00 13:00 15:00 17:00 19:00 21:00 23:00 Time of Day Electricity Consumption (kWh) THH− total drop = 310.8 kWh Figure 4.25: Consumption profile of building THH during DR In Figures 4.25- 4.28, we show the energy consumption profile of a few buildings which achieved the highest curtailment values. The x-axis denotes the time of the day and the y-axis denotes the energy consumption value in kWh. For simplicity, we assume that the consumption in the absence of DR would be the observed value at the start of DR event (1pm). This is denoted by the red line in the figures. The buildings THH, MRF and WPH followed the technique of Variable Frequency Drive Speed Reset (VFD) [50] and obtained a curtailment value of 310, 145 and 98 kWh respectively. The building VKC 73 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 20 40 60 80 01:00 03:00 05:00 07:00 09:00 11:00 13:00 15:00 17:00 19:00 21:00 23:00 Time of Day Electricity Consumption (kWh) VKC− total drop = 169.2 kWh Figure 4.26: Consumption profile of building VKC during DR adopted a combination of Variable Frequency Drive Speed Reset (VFD) and Equipment Duty Cycling (DUTY) [50] and obtained a curtailment value of 169 kWh. However, not every building achieved the expected curtailment value. Figures 4.29- 4.30 show the consumption profile of buildings BHE and DRB which achieved a nega- tive curtailment value. Due to some practical limitations, these buildings failed to adopt the strategies the FMS signaled them to implement. 74 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 10 20 30 40 01:00 03:00 05:00 07:00 09:00 11:00 13:00 15:00 17:00 19:00 21:00 23:00 Time of Day Electricity Consumption (kWh) MRF− total drop = 145.4 kWh Figure 4.27: Consumption profile of building MRF during DR ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 20 30 40 50 60 01:00 03:00 05:00 07:00 09:00 11:00 13:00 15:00 17:00 19:00 21:00 23:00 Time of Day Electricity Consumption (kWh) WPH− total drop = 98.6 kWh Figure 4.28: Consumption profile of building WPH during DR 75 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 12 16 20 24 01:00 03:00 05:00 07:00 09:00 11:00 13:00 15:00 17:00 19:00 21:00 23:00 Time of Day Electricity Consumption (kWh) BHE− total drop = −12.2 kWh Figure 4.29: Consumption profile of building BHE during DR ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 40 50 60 01:00 03:00 05:00 07:00 09:00 11:00 13:00 15:00 17:00 19:00 21:00 23:00 Time of Day Electricity Consumption (kWh) DRB− total drop = −39.8 kWh Figure 4.30: Consumption profile of building DRB during DR 76 Chapter 5 Cost Optimal Supply Demand Matching In this section, we develop algorithms for discrete supply demand matching by con- trolling both supply and demand nodes. Algorithms developed in Sections 5.1, 5.2 and 5.4.1 assume a single control mode while the algorithm developed in Section 5.4.2 can be used in dual control mode. 5.1 Minimum Cost Supply Demand Matching Given the Smart Grid Model above, the objective of the Minimum Cost Supply Demand Matching Algorithm is to determine node-curtailment strategy pairs for each interval of the curtailment horizon such that: (1) The curtailment target t for each interval t is achieved, (2) the cost is minimized for the entire curtailment horizon, and (3) the aggregate curtailment across the entire curtailment horizon is no more than . This problem, as we show in Section 3.6 is NP-hard. We first formulate the problem using an Integer Linear Program (ILP). However, the time complexity for solving ILPs is exponential. Hence, we develop a polynomial time approximation algorithm for the same. More formally, we develop an algorithm with a runtime which is polynomial in the input sizeM;MN andT and 1 , where is an approximation guarantee (accuracy) parameter, which ensures that objective of the problem is minimized and in the worst 77 case the constraints are violated by a maximum factor of (1). Here M = M s or M d based on whether supply curtailment is being performed or demand curtailment. To develop the approximation algorithm, we use a dynamic programming algorithm (Equation 5.6) which for each intervalt, determines the cost of achieving various cur- tailment values, each of which is t . We then use another dynamic programming algorithm (Equation 5.7) to combine the results of each interval to achieve an aggre- gated curtailment value of with minimum cost. The sizes of the tables of both the dynamic programming algorithms are proportional to the maximum possible cost. This leads to very large runtime to solve the problem exactly. Hence, we scale and round the costs. The scaling and rounding causes several curtailment values to be indistinguish- able, this introduces error into the value of our solution. Thus the resulting algorithm is an approximation algorithm (as opposed to exact algorithm) which produces an approx- imate solution in polynomial time which is independent of the maximum cost. We provide the worst case approximation guarantee for the algorithms. In the following sections, to simply notations ;c;x denote the variables correspond- ing to either supply or demand nodes based on whether supply curtailment is being performed or demand curtailment. 5.1.1 ILP Formulation The ILP formulation for the Minimum Cost Supply Demand Matching problem is as follows: 78 Minimize : M X b=1 N X j=1 T X t=1 c bj (t)x bj (t) (5.1) s:t: M X b=1 N X j=1 bj (t)x bj (t) t 8t (5.2) M X b=1 N X j=1 T X t=1 bj (t)x bj (t) (5.3) N X j=1 x bj (t) == 1 8b;t (5.4) x bj (t)2f0; 1g 8b;j;t (5.5) Equation 5.2 ensures that the curtailment target for each time interval is achieved. Equation 5.3 ensures that the aggregate curtailment is less than the maximum limit . Equation 5.4 ensures that each node in each time interval follows exactly one strategy (possibly the default strategy with 0 curtailment value). 5.1.2 Approximation Algorithm Let min = min t f t g be the smallest curtailment target among all the intervals. Define = min M . For each bj (t), defineb bj (t) =d bj (t) e. Similarly, define b and b t 8t. We refer to these values as rounded curtailment values. For each interval t, we define a function t :Z + [f0gf1;:::;Mg. t (b ;b) denotes the minimum cost required to achieve a curtailment value ofb using nodes 1;:::;b whereb2f1;:::;Mg. t can be defined recursively as: 79 t (b ;b) = 8 > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > : min j fc bj (t)g ifb = 1 and 9jjb =b bj (t) 1 ifb = 1 and b ! =b bj (t)8j 1 ifb < 0 min j f t (b b bj (t);b 1) +b bj (t)g otherwise (5.6) Algorithm 5: Determine node strategy pairs given (c;b )2S t Input: (c;b ),t 1 x bj (t)2X(t) 08b;j 2 cur b 3 forb =M to 1 do 4 ifb6= 1 then 5 j argmin j f( cur b bj (t);b 1) +c bj (t)g 6 else 7 j jjb bj (t) == cur 8 end 9 x bj (t) 1 10 cur cur b bj (t) 11 end Output: OutputX(t), the list of curtailment strategies to be followed by each node in intervalt. The dynamic program can be solved by creating a table of sizekM, wherek = b for each interval. For notational simplicity, we refer to table using the same variable t as the recursive function. Once the table is filled, for each interval t, we define a set S t =f( t (b ;M);b )jb b t g. For any element ( t (b ;M);b )2 S t , Algorithm 5 can be used to determine the strategies to be followed by each node to achieveb with cost t (b ;M) in the intervalt. 80 Algorithm 6: Minimum Cost Supply Demand Matching Algorithm Input:C(t); (t); t ;8t; 1 Computeb bj (t); b t 8t; b 2 Fill the table t 8t using equation 5.6 3 ComputeS t =f( t (b ;M);b )jb b t g8t 4 Fill the table using equation 5.7 5 b cur argmin b f(b ;T )jb b g 6 ifb cur == then 7 No curtailment strategies exist 8 Exit Algorithm 9 end 10 fort =T to 1 do 11 ift6= 1 then 12 j argmin jj (c j ;b j )2St f(b cur b j ;t 1) +c j g 13 else 14 j jjb =b j ; (c j ;b j )2S t 15 end 16 Run Algorithm 1 with (c j ;b j ) to getX(t) 17 b cur b cur b j 18 end Output:X(t)8t, the list of curtailment strategies to be followed by each node in each interval Now, givenS t 8t2f1;:::;Tg, we need to select exactly one elemente t = (c t ;b t )2 S t 8t such that P T t=1 b t b and P T t=1 c t is minimized. We define a function :Z + [ f0gf1;:::;Tg. (b ;t) denotes the minimum cost required to achieve the curtailment value ofb and time intervals 1;:::;t. can be defined recursively as: 81 (b ;t) = 8 > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > : min j fc j g ift = 1 and 9jjb =b j ; (c j ;b j )2S t 1 ift = 1 and b ! =b j 8jj (c j ;b j )2S t 1 ifb < 0 min jj (c j ;b j )2St f(b b j ;t 1) +b j g otherwise (5.7) This dynamic program requires a table of sizekT , wherek = b . Again, for nota- tional simplicity, we refer to table using the same variable as the recursive function. Algorithm 6 can be used to determine the curtailment achieved in each interval and the corresponding node strategy pairs. Theorem 7. Algorithm 6 is a polynomial time algorithm for minimum cost supply demand matching which in the worst case violates the maximum curtailment constraint (Equation 5.3) by at most (1 +) factor and violates the per interval curtailment target constraint (Equation 5.2) by at most (1) factor. Proof. For a curtailment value , we say thatb is the rounded curtailment value. Also, we call as the unrounded curtailment value. Correctness: We define Bucket ofb as the rangeb < b i.e. all the curtailment values which get rounded tob . C b denotes the cost assigned to the bucket in table . For each intervalt, using induction on Equation 5.6, it is easy to show thatC t b 8t; 0b b will be the minimum cost required to achieve any curtailment value in the Bucket ofb for table t . Hence, for each elemente t = (c t ;b t ),c t is the minimum possible cost to achieve any curtailment value 82 in the Bucket ofb t corresponding to the table t . Now, if we can show that the range of curtailment values covered by the elements in S t 8t contains the curtailment value chosen by any optimal solution, and that the cumulative curtailment value of any optimal solution at any time t is contained in the entries (:;t) (i.e. entries corresponding to time intervalt), then using induction on table, we can show the aggregate cost of the solutions obtained by Algorithm 6 will be less than or equal to the optimal solution. Now, in each interval, the lowest rounded curtailment valueb considered to create the setS t isb t . Hence, the lowest unrounded curtailment satisfies: t < t i.e. t . Now, the maximum rounded curtailment valueb considered by Algorithm 2 in the table is b . So, the maximum curtailment value considered in the corresponding bucket is b . Hence, the range of curtailment values considered by Algorithm 6 covers the range in which optimal solution can reside and so Algorithm 6 does not miss any possible solution of a lower cost. Runtime: The Algorithm fillsT tables each of which is of sizeM b . We assume that =O(T min ). Each entry of requiresO(N) time. Hence, the total runtime for line 2 and 3 of Algorithm 6 isO( M 2 NT 2 ). The number of entries in table isO(T b ). Each entry requiresO( b ) time. Hence, the total required time to fill isO(T b 2 =O( T 3 M 2 2 )). Once the tables are filled, the for loop requiresO(TMN b ) to output the strategies. Hence, the algorithm is polynomial in the input sizeM;N;T and 1 . Approximation Guarantee: Let b x = P M i=1 P T t=1 b it be the solution from our algo- rithm, whereb it denotes the rounded curtailment by node i in time t. From line 3 of the algorithm, for all t, we know that P M i=1 b it b t . Also, it b it it + 1 by the definition ofb it . This implies that t P M i=1 ( it + 1) P M i=1 it +M. So, t M P M i=1 it =) P M i=1 it t min (1) t as t min . Hence, in 83 each interval, the curtailment target constraint (Equation 5.2) is violated by a maximum factor of (1). Now, from line 5 of the algorithm, b x b . So, x + 1. This means that x + min M (1 + M ) (1 +) as min and M 1. Hence, the aggregate curtailment target constraint (Equation 5.3) is violated by a maximum factor of (1 +). 5.2 Minimum Cost Supply Demand Matching with Fairness Curtailment from a node leads to a loss of utility for the node. Hence, it would be unfair to force some nodes to incur losses due to high curtailment while leaving others with minimal curtailment. The Minimum Cost Supply Demand Matching Algorithm discussed in the previous section does not take fairness into account and can lead to solutions with uneven curtailment values from the nodes. We address the issue of fair- ness in this section by assigning a curtailment budget range (which can be set by the grid operator) to each node. The algorithm, by ensuring that no node incurs a curtailment more or less than its budgeted range over each curtailment horizon, ensures that supply demand is done in a fair manner. Similar to the previous problem, we first develop an ILP formulation for this prob- lem. We then relax the ILP into a Linear Program (which can be solved in polynomial time) and round back the results to integers. This rounding, however, violates certain constraints and increases the objective value. We provide guarantees on the worst case violation of the constraints and the increase in the objective value in the worst case. 84 5.2.1 ILP Formulation Let [ b B b ;B b ] be the curtailment budget for nodeb with b 2 [0; 1]. BothB b and b are determined by the grid operator. The problem of minimum cost supply demand matching with fairness can be formulated using the following ILP: Minimize : M X b=1 N X j=1 T X t=1 c bj (t)x bj (t) (5.8) s:t: M X b=1 N X j=1 bj (t)x bj (t) t 8t (5.9) M X b=1 N X j=1 T X t=1 bj (t)x bj (t) (5.10) b B b N X j=1 T X t=1 bj (t)x bj (t)B b 8b (5.11) N X j=1 x bj (t) == 1 8b;t (5.12) x bj (t)2f0; 1g 8b;j;t (5.13) Equation 5.11 is the additional constraint added that ensures that each node curtails an amount within its budgeted range. 5.2.2 Approximation Algorithm The ILP formulated in Section 5.2.1, when relaxed to a Linear Program will lead to unbounded errors. Hence, we first redefine the ILP. We define 0 b (t) 0 bj (t) = bj (t), where 0 b (t), is the maximum operational value for node b at time t and 0 bj (t) is the 85 reduced operational value of nodeb at timet by following strategyj with curtailment value bj (t). the redefined ILP is as follows: Minimize : M X b=1 N X j=1 T X t=1 c bj (t)x bj (t) (5.14) s:t: M X b=1 N X j=1 0 bj (t)x bj (t) M X b=1 0 b (t) t 8t (5.15) M X b=1 N X j=1 T X t=1 bj (t)x bj (t) (5.16) N X j=1 T X t=1 bj (t)x bj (t)B b 8b (5.17) N X j=1 T X t=1 0 bj (t)x bj (t) T X t=1 0 b (t) b B b 8b (5.18) N X j=1 x bj (t) == 1 8b;t (5.19) x bj (t)2f0; 1g 8b;j;t (5.20) We note that the above ILP is an instance ofk column sparse Packing Integer Prob- lem (PIP) i.e. a packing integer problem in which there in an upper limitk on the number of constraints a variable can occur in [15].k column sparse PIPs have aek+O(k) factor randomized rounding based approximation algorithm for a fixedk [15]. This implies that even if we assume P M b=1 B b i.e. the constraint (5.16) is redundant and that b = 0;8b each variable x ij (t) appears in three constraints: (1) Per interval curtailment target constraint (5.15), (2) Budget constraint (5.17), and (3) constraint (5.19). Hence, using the approximation algorithm developed in [15] we can get a worst case bound of> 8 which is unsatisfactory for reliable grid operations. Hence, to develop an approximation algorithm with tighter theoretical worst case bounds, we first make the following assumption: The costs c bj (t) are a function of 86 bj (t) i.e.,c bj (t) =f( bj (t)). We will derive approximation guarantees when the func- tionf is linear and when it is quadratic. In order to develop a bounded approximation algorithm, we first relax the ILP to a linear program i.e., we replace Equation 5.20 with 0 x bj (t) 18b;j;t and solve the Linear Program. Solving a linear program takes polynomial amount of time using methods such as inter-point or ellipsoid [20]. How- ever, the solution will contain fractional values for the decision variablesx bj (t);8b;j;t which need to be rounded to 0 or 1 for a valid solution. Now, naively rounding the decision variables leads to errors which are unbounded. Hence, we develop Algorithm 7 which is a novel rounding algorithm which guarantees that the constraints are violated by at most a factor of 2 in the worst case. For eachb;t, the algorithm works by com- puting expected curtailment 0 = P N j=1 bj (t)x bj (t) and rounding it to the curtailment value bj (t) nearest to it. We have the following two results for this algorithm. Theorem 8. For a linear cost functionf, Algorithm 7 is a (2,2)-factor Minimum Cost Supply Demand Matching with Fairness Algorithm. The cost of the solution achieved by Algorithm 7 is at most twice the optimal while the budget and targeted curtailment constraints (Eqs 5.16 and 5.17) are violated by at most a factor of two. Theorem 9. For a quadratic cost functionf, Algorithm 7 is (4,2)-factor algorithm. Proof. Let 0 bt = P N j=1 bj (t)x bj (t) be the curtailment value obtained for nodeb in time interval t. Let i bt and i+1 bt be the curtailment values of strategies between which 0 bt falls i.e. i bt 0 bt i+1 bt . Now, if 0 bt is rounded up to i+1 bt , this implies 0 bt ( i bt + i+1 bt )=2 i+1 bt =2. Summing over all values of b;t ensures that Equation 5.16 and the upper bound of Equation 5.17 are violated by at most a factor of 2. Similarly it can be shown that if the objective function is linear, it will be bounded by a factor of 2 and if it is quadratic, it will be bounded by a factor of 4. 87 Now, in order to provide a bound on the constraint violation of Equations 5.15 and 5.18, we note that if 0 bt is rounded up, the bounds are trivially satisfied. How- ever, if they are rounded down, then guarantee no longer holds. Hence, we make an assumption that i+1 bt (2k 1) i bt . Using this assumption, we can ensure that in the worst case i bt 0 bt =k i.e. the constraint is violated at most by a factor ofk. The analysis above can be used by the grid operator to improve the efficiency of supply demand matching. The grid operator can setup curtailment configurations (e.g. load curtailment strategies such as GTR, etc.) such that the difference between two consecutive curtailment values is not vary large. However, it is not always possible to control the curtailment configurations. Hence, the grid operator can develop techniques by noticing that a higher curtailment value can be expected to have less violations for two reasons: (1) very few values of 0 bt will have i bj = 0 and (2) since a large number of non-zero curtailment strategies will be selected, several of them will be rounded up thus reducing the possibility that Equation 5.10 is violated by a large factor. Hence, if the required curtailment value is small for a curtailment horizon, the grid operator can reduce the number of nodes which participate in the curtailment horizon. Note that the above guarantees are worst case guarantees. In practice, the perfor- mance is significantly better as shown using the experimental results. 5.3 Online Algorithm for Fair Supply Demand Match- ing The algorithms discussed in the previous two sections require the availability of the curtailment prediction data for the entire horizon. However, certain scenarios require supply demand matching in an online manner. At the beginning of each interval, the 88 Algorithm 7: Minimum Cost Supply Demand Matching with Fairness Input:C(t); (t); t ;8t; ;B b ; b 8b 1 x bj (t)2X(t) 08b;j;t 2 Relax the ILP to an LP by replacing Equation 5.13 with 0x bj (t) 18b;j;t 3 Solve LP to obtain solutionx bj (t)8b;j;t 4 foreachb;t do 5 0 P N j=1 bj (t)x bj (t) 6 Let bi (t) 0 bi+1 (t) 7 if ( 0 bi (t)) ( bi+1 (t) 0 ) then 8 x bi+1 (t) 1 9 else 10 x bi (t) 1 11 end 12 end Output:X(t)8t, the list of curtailment strategies to be followed by each node in each interval data is made available and supply demand matching needs to be performed in a myopic way. We develop a greedy online heuristic algorithm (Algorithm 8) for this problem. The algorithm finds a minimum cost way to achieve the curtailment target o l = t for the current interval while ensuring that no node curtails more than the budget for the current interval. The upper limit of the budget for the current interval B o b for each node b is determined by multiplying the ratio of B b P T t=1 t with o l . The lower bound is simply b B b . The upper bound on curtailment o u is determined by multiplying the ratio of P T t=1 t with o l . Let o denote the curtailment matrix and C o denote the cost matrix with o bj being the curtailment for node b following strategy j and c o bj being its cost. We remove any values o bj which are outside the curtailment budget range from the curtailment matrix o . Note that the values ; t ;B b 8b are unknown for the current supply demand matching horizon. They can be obtained from some past horizon. 89 Algorithm 8: Online Algorithm for Fair Supply Demand Matching Input: o l ; o ;C o ; ; t 8t;B b ; b 8b 1 o u P T t=1 t o l 2 B o b B b P T t=1 t o l 8b 3 ComputeC o min ;C o max ;b c o bj ; b C o max similar to Algorithm 6 using o in which curtailment values outside curtailment budget range are removed. 4 Fill table using equation 6 5 c min b c fb cj o l (b c;M) o u g 6 Run Algorithm 5 with (c ; (c ;M)) to getX o Output:X o , the list of curtailment strategies 5.4 Supply Demand Matching with Network Con- straints The Supply Demand Matching algorithms discussed in Section 5.1 and 5.2 are agnostic to the underlying network model of the smart grid. This can lead the framework to output solutions which violate network capacity constraints. In this section, we will develop a Supply Demand Matching Framework which explicitly models the Smart Grid network (Section 3.4) and ensures that the solutions output satisfy the constraints imposed by the capacity limitations of the various components of the Smart Grid. 5.4.1 Transformer Level Supply Demand Matching Traditional smart grids are not equipped with bi-directional LV Feeders i.e. the feeders can carry power from the sub-station to the transformers but not in the other direction. This implies Cap f = 08f 2F. Hence, Supply Demand Matching needs to be per- formed at each transformer of the Smart Grid. We develop Algorithm 9 for the same. Algorithm 9 takes as input the time varying cost and curtailment matrices for both supply and demand nodes, the distributor capacities and the per interval maximum sup- ply (demand) for each supply (demand) node. In lines 1 and 2, it initializes the decision 90 Algorithm 9: Transformer Level Supply Demand Matching Algorithm Input:C s (t);C l (t); s (t); l (t)8t;Cap tx 8tx;S bs (t)8b s 2S8t;L b d (t)8b d 2 D;8t 1 X s (t) 0 2 X l (t) 0 3 foreachtx2T do 4 foreacht2f1;:::;Tg do 5 Sup tx (t) = P bs2T (tx)\S (S bs (t)) 6 Dem tx (t) = P b d 2T (tx)\D (L b d (t)) 7 s (t) maxfSup tx (t)Cap tx ; 0g 8 l (t) Dem tx (t)Sup tx (t) + s (t) 9 if l (t) 0 then 10 s (t) l (t) 11 l (t) 0 12 end 13 end 14 s grid operator determined upper bound on solar curtailment. 15 l grid operator determined upper bound on load curtailment. 16 X tx s (t) output of Algorithm 6 with inputsC s (t); s (t); s (t)8t; s 17 X tx l (t) output of Algorithm 6 with inputsC l (t); l (t); l (t)8t; l 18 X s (t) X s (t) +X tx s (t)8t 19 X l (t) X l (t) +X tx l (t)8t 20 end Output:X s (t);X l (t) the list of solar and load curtailment strategies for each node in each interval matrices to 0. Now, for each transformertx, in each time intervalt (foreach blocks at lines 3 and 4), it calculates the total solar supply (line 5) and the total demand (line 6). If the supply is more than the capacity of the distributor, the curtailment target s (t) is set to be the excess supply (line 7). Now, if demand is greater than the total supply (minus the supply curtailment required to satisfy the distributor capacity constraint), a demand curtailment target l (t) is identified (line 8). The If block starting at line 9 han- dles the cases where supply is more than the demand. In such a scenario, the supply should be curtailed for supply demand matching. Algorithm 6 is then used to determine the per interval load and supply curtailment strategies for each node (lines 16 and 17). 91 The upper bounds as required by Algorithm 6 can be determined by the grid operator (lines 14 and 15). Note that in the interest of simplicity, we make an assumption that only the required amount of curtailment is performed and any unnecessary curtailment is avoided. For example, for a transformer, with solar surplus, if a lower cost solution can be found by performing both load and solar curtailment, it is not considered. Our algorithm can easily be extended to handle the cases where the assumption is not true. Theorem 10. Algorithm 9 is a polynomial time algorithm for transformer level mini- mum cost supply demand matching with network constraints which in the worst case violates the distributor constraints (Equation 3.3) by at most (1) factor. Proof. In each transformer, Algorithm 6 is being called withCap tx as the curtailment target. As Algorithm 6 is a polynomial time algorithm which produces a solution of minimum cost while violating the curtailment target by at most (1), so does Algo- rithm 9. 5.4.2 Smart Grid Level Supply Demand Matching The increased distributed generation in future smart grids will lead to a more active participation by the consumers and will required bi-directional power exchanges [29] using LV-feeders. This implies that the maximum power that can be injected from the transformers i.e. Cap f 0. This will allow us to perform more cost efficient supply demand matching by optimally distributing the supply from surplus transformers to the deficit ones. We develop Algorithm 10 to perform smart grid level supply demand matching at any given time interval t. The brief description of Algorithm 10 is as follows: For each transformer, the algorithm first checks whether the solar generation from the nodes attached is within its capacity or not. If not, then the nodes are marked for the required 92 curtailment. Now, in the resulting grid (the grid in which the nodes have performed the required curtailment), the algorithm partitions the feeders into two sets: one with supply surplus (S sup ) and one with load surplus (S dem ). For each set, the minimum costs of performing all possible curtailments are identified. Finally, the algorithm chooses sup- ply and load curtailments such that the cost is minimized and the total supply after solar curtailment equals the total demand after load curtailment. Unlike the transformer level supply demand matching in Section 5.4.1, which distributes curtailment temporally to minimize costs, this algorithm performs spatial curtailment distribution in each interval separately. We define additional notations sf (t); df (t);C sf (t);C df (t);X sf (t);X df (t) denoting the sub-matrices of the original matrices defined in Section 3.4 corresponding to supply (s) or demand (d) nodes attached to feederf for the given time intervalt. Let min = min s;d;f f min sf (t); df (t)g, where min sf (t) and df (t) are defined in lines 8 and 11 of Algorithm 10. Also letF max = max f fjF(f)jg Define = min Fmax . We defineb sf (t) = d sf (t)e,b df (t) =d df (t)e. We refer to these values as rounded curtailment values. Similarly, we defineb min sf (t);b max sf (t). For each feeder f, we define a function f : Z + [f0gf1;:::;jF(f)jg, wherejAj outputs the size of the setA. f (b ;b) denotes the minimum cost required to achieve a curtailment value ofb using nodes 1;:::;b whereb2f1;:::;jF(f)jg. To avoid redundancy in notations, if f is used to denote supply curtailment, any load consuming nodeb is assumed to have 0 supply curtailment strategies. f can be defined recursively as: 93 f (b ;b) = 8 > > > > > > > > > < > > > > > > > > > : min j c bj (t) ifb = 1 and9jjb =b bj (t) 1 ifb = 1 andb ! =b bj (t)8j 1 ifb < 0 min j f f (b b bj (t);b 1) +c bj (t)g otherwise (5.21) Note that this is similar to the dynamic program of Equation 5.6 and hence Algo- rithm 5 can be used to determine the strategies to be followed by each node connected to the feeder f. For each f, if f 2S sup , we define S f =f( f (b ;M);b )j b min sf (t) b b max sf (t)g, where b max sf (t); b min sf (t) are defined in line 11. The minimum curtail- ment requirement ensures that the feeder capacity constraints are met. Iff2S dem , then S f =f( f (b ;M);b )j b df (t)b g. Now, we create two dynamic programming tables sup and dem using the sets S f 8f 2 S sup and S f 8f 2 S dem respectively. The tables are built using functions sup :Z + [f0gf1;:::;jS sup jg and dem :Z + [f0gf1;:::;jS dem jg. The func- tions sup (b ;f) denotes the minimum cost required to achieve the curtailment value ofb using supply nodes from feedersf1;:::;fg and is defined using the recursive function: sup (b ;f) = 8 > > > > > > > > > < > > > > > > > > > : min j c j iff = 1 and9jjb =b j ; (c j ;b j )2S f 1 iff = 1 andb ! =b j 8jj (c j ;b j )2S f 1 ifb < 0 min jj (c j ;b j )2S f f(b b j ;f 1) +c j g otherwise (5.22) 94 dem (b ;f) is defined similarly. Leveraging the tables built using Equations 5.21 and 5.22, Algorithm 10 performs minimum cost smart grid level supply demand match- ing. b S and b D used in equation 19 are the rounded aggregated solar and load surplus in the Smart Grid respectively. A detailed line by line description of Algorithm 10 is provided later in Section 20. We have the following theorem. Theorem 11. Algorithm 10 is a polynomial time algorithm for smart grid level minimum cost supply demand matching with network constraints which in the worst case violates the feeder and distributor constraints (Equations 3.3 and 3.6) by at most (1) factor. Moreover, the curtailment achieved by Algorithm 10 is within (1 2) factor of the optimal curtailment target which can mitigate the supply demand mismatch. Proof. Correctness of the algorithm can be proved using induction, hence we exclude the details. Runtime: The dominating term in the runtime complexity will be the creation of the dem and sup tables. Each of these tables will require creation of tables requiring O(M 2 N b ), by assuming thatjFj = O(M) i.e. the number of nodes and max f fjF(f)jg =O(M) and = max s;d;f f max sf ; df (t)Mg. This is equal toO( M 4 N ) with the assumptions thatF max = O(M) and = O(M min ). The creation of dem and sup tables after creating the tables will requireO(M b 2 ) = O( M 5 2 ). Hence, the algorithm is polynomial in the input size and 1 . The assumption = O(M min ) implicitly implies that no transformer level load/solar curtailment value is ”too small”. If this is not true, we can isolate the trans- former with the smallest curtailment target and solve it independently. We can repeat the process until the condition is satisfied. This will lead to a possibly higher cost solution but will significantly reduce the computation time required. Approximation Guarantee: Let b x = P i b i be the rounded curtailment values selected from the nodes connected to a transformer. Now, we know that b cap b x , 95 where b cap is the minimum curtailment required to satisfy capacity constraints (as deter- mined in Equations 8 and 11). So, cap P i ( i + 1) P i i +F max , where cap and i are the unrounded curtailment values. This implies cap min x =) x (1) cap , where again x is the unrounded curtailment value. Hence, the capacity constraints are violated by a maximum factor of (1). Now, in order to prove the second part of the theorem, without loss of generality, let supply be surplus in the Smart Grid. Let b S be the aggregate supply and b D. Let b sx = P i b si be the solar curtailment values chosen and let b dx = P i b di be the chosen load curtailment values. Let = SD be the required unrounded curtailment value and let x = sx dx be the achieved unrounded curtailment value. Now, we know b Sb sx = b Db dx . By rearranging and unrounding, we can get 1: S + dx D + sx + +M and 2: D + sx S + dx + +M. Using 1, by substitution, we get min (1 + 1 M ) x =) (1 2) x asM 1. Using 2, by substitution, we get x +(1 + 1 M ) min =) (1 + 2) asM 1. So, the curtailment achieved by Algorithm 10 is within (1 2) factor of the optimal curtailment. Algorithm 10 Description Algorithm 10 takes as input the cost C s ;C l and curtailment s ; l matrices for all the supply and demand nodes, the capacity of each distributorCap tx and each feederCap f and the maximum supply S bs and maximum demand L b d for each supply or demand node respectively. It then initializes two empty sets (line 1). S sup denotes the feeders where the maximum possible supply is greater than the maximum possible demand and S dem are the feeders with maximum possible demand greater than supply. The decision variables are all set to 0. For each feeder (for each block starting at line 2), the algo- rithm calculates the maximum possible supply and demand (lines 3 and 4). In line 5, the algorithm identifies the minimum supply curtailment required to meet the distributor 96 Algorithm 10: Smart Grid Level Supply Demand Matching Algorithm for time intervalt Input:C s (t);C l (t); s (t); l (t);Cap tx 8tx;Cap f 8f;S bs (t)8b s 2 S;L b d (t)8b d 2D 1S sup ;S dem ;X sf (t) 0;8s2F(f)\S;8f; X df (t) 0;8d2F(f)\D;8f 2 foreachf2F do 3 Sup f (t) = P bs2F(f)\S (S bs (t)) 4 Dem f (t) = P b d 2F(f)\D (L b d (t)) 5 sf (t) maxfSup f (t)Cap F tx (f) ; 0g 6 ifDem f (t)Sup f (t) sf (t) then 7 S dem f 8 df (t) Dem f (t)Sup f (t) + sf (t) 9 else 10 S sup f 11 max sf (t) Sup f (t)Dem f (t); min sf (t) maxfSup f (t)Cap f ; sf (t); 0g 12 end 13 end 14 foreachf2S dem do 15 X sf (t) output of MinCB with inputsC sf (t); sf (t); sf (t); sf (t) 16 Create dem using Equations 5.21 and 5.22, LetS demf be the sets created after Equation 5.21 17 end 18 foreachf2S sup do 19 Create sup using Equations 5.21 and 5.22, LetS supf be the sets created after Equation 5.21 20 end 21 (b ) sup ( b S b D +b ;jS sup j) + dem (b ;jS dem j)8b 0 22 Call Algorithm 11 to getX sf (t);X df (t) Output:X sf (t);X df (t) the list of curtailment strategies for each node in the intervalt. constraints. Then, in the if-then-else block starting at line 6, if the maximum demand in a feeder is greater than the supply (after curtailment required for satisfying distrib- utor constraints), the feeder is inserted intoS dem and a maximum demand curtailment target is identified to ensure a feasible supply demand matching configuration. Other- wise, the feeder is inserted intoS sup and a minimum and maximum supply curtailment 97 Algorithm 11: Finding node strategy pairs from Input: 1 b cur argmin b f(b )g 2 forf =jS sup j to 1 do 3 iff6= 1 then 4 j argmin fjj (c j ;b j )2S supf g f sup (b cur b j ;f 1) +c j g 5 else 6 j jjb =b j ; (c j ;b j )2S supf 7 end 8 Run Algorithm 5 with (c j ;b j ) over f to getX sf (t) 9 b cur b cur b j 10 end 11 b cur argmin b f(b )g 12 forf =jS dem j to 1 do 13 iff6= 1 then 14 j argmin fjj (c j ;b j )2S demf g f dem (b cur b j ;f 1) +c j g 15 else 16 j jjb =b j ; (c j ;b j )2S demf 17 end 18 Run Algorithm 5 with (c j ;b j ) over f to getX df (t) 19 b cur b cur b j 20 end Output:X sf (t);X df (t) the list of curtailment strategies for each node in the intervalt. is identified. The minimum supply curtailment is required for satisfying the feeder and distributor constraints while the maximum denotes the curtailment after which supply demand matching is not possible due to less supply. In the foreach blocks (lines 14 and 18), setsS demf andS supf are created by filling dynamic programming tables dem and sup respectively. Additionally, for the demand surplus feeders, any supply curtail- ment required to meet the distributor constraints is performed using Algorithm 5 (with time horizon of just 1 interval). Now, in line 21, for each valid rounded demand curtail- mentb , the cost of operation is calculated. Cost of operation is the sum ofb demand curtailment and b S ( b Db ) supply curtailment. The latter curtailment is required to match supply with demand. Here, b S and b D are the rounded aggregate maximum 98 solar and demand surplus in the smartgrid respectively. These are the values which are available to perform smart grid level supply demand matching while ensuring that within each feeder, the distributor and feeder constraints are satisfied and feasible sup- ply demand matching configurations are maintained. Now it calls Algorithm 11 which in lines 1 to 20, using the value ofb identified to ensure minimum cost, the dynamic programming tables are traversed to find suitable strategies for each node in each feeder and the corresponding decision variables are set. As the tables f are similar to the table filled in Algorithm 6, Algorithm 5 is used to traverse the same. 5.5 Modeling Storage for Supply Demand Matching Algorithms Operational states of energy storage systems are continuous in nature as opposed to the discrete operational states of supply and demand nodes considered in the work. It is trivial to incorporate storage in the Integer Linear Program based algorithm of Sec- tion 5.2 by adding continuous variables denoting storage charge/discharge and convert- ing it into a Mixed Integer Linear Program (MILP). However, incorporating storage into algorithms developed in Section 5.1, Section 5.4.1 and Section 5.4.2 requires us to dis- cretize the charge/discharge parameter in such a way that the runtime is not affected. In order to incorporate storage in dynamic programming based approximation algo- rithms, we define b R =d R e; b G =d G e and b D =d D e. For each interval t, we define the charging strategies of storage as b G t =f0; 1;:::; b Gg and discharging strategies as b D t =f0; 1;:::; b Dg. Using these variables, the storage operation can be defined as b R t+1 = b R t + G b G t 1 D b D t (5.23) 99 and the constraints are defined as follows: T 1 D b R t+1 Discharge Mode (5.24) b R t+1 b R G T Charge Mode (5.25) 0 b D t b D (5.26) 0 b G t b G (5.27) Now, if the algorithm chooses b D t as discharge, minf b D t ;Dg is discharged. Sim- ilarly, if the algorithm chooses b C t as charge, minf b C t ;Cg is charged. This ensures that the maximum charging/discharging constraints are met. Equations 5.24 and 5.25 ensure that the storage operation is always within storage capacity. This is because the maximum difference between the charge values obtained by our algorithm: G x and actual charge: G a is bounded by G T and the discharge range is bounded byT 1 D . Moreover, when b R T T 1 D , the storage cannot be used in discharge mode or when b R T b R G T , it cannot be used in charge mode. It is evident that this storage model leads to aT amount of storage being unused, where = ( G + 1 D ). The value min used to define can be modified to provide an upper bound on the unused storage. For example, if the storage capacityR 1 k l M , where l is the aggregate of minimum curtailment target across all the intervals, then modifying min to 1 k min will bound the unused storage byT = min T kM l kM R. The condition R 1 k l M implies that the storage capacity is at least 1 k times the average curtailment produced by a node in an interval to achieve the target l . However, the downside here is that this leads to a (k) 2 factor increase in the runtime. The typical charge/discharge efficiency for Li-ion batteries is 80-90% [4], so 2:15. Now, when the algorithms fill the dynamic programming table in Equations 5.6 and 5.21, instead of iterating throughN strategies to find the minimum, the algorithm 100 iterates through b D + 1 strategies if the storage is used for discharging or b G + 1 strate- gies if it is used for charging. This increases the runtime of filling by a factor of O( b N ). However, the runtime of Algorithm 6 and Algorithm 10 still remain O( T 3 M 2 2 ) andO( M 5 2 ) respectively. 5.6 Results and Analysis In addition to providing theoretical guarantees, we perform practical evaluations of the algorithms developed. We implemented the algorithms using MATLAB [8]. The LP and ILP algorithms were implemented using IBM ILOG Cplex Optimization Studio [1] which is a toolkit for mathematical and constraint programming. The experiments were performed on Dell optiplex with 4-cores and 4 GB RAM. A supply demand matching horizon of 32 15-min intervals was considered with 16 intervals of load curtailment and 16 intervals of solar curtailment. We evaluated the algorithms by varying the (L,U) pair where, for notational simplic- ity, L= P T t=1 t and U= (Section 5.2.1). Note that L and U are the lower and upper bound on the curtailment to be achieved in the curtailment horizon and hence represents the feasible curtailment range. The costs of the strategies were evaluated using the func- tionf( ) = 2 2 , where is the curtailment value of the strategy. Section 5.6.1 describes the input dataset generation. As the dataset was generated from historical data, perfect knowledge of the future was assumed with no prediction errors. 5.6.1 Dataset We obtained the load curtailment data from the demand response implementation on our University Campus. Our campus consists of 150 DR enabled nodes (buildings) each of which can follow 6 load curtailment strategies. The load curtailment values for each 101 node-strategy pair was generated using algorithms mentioned in [12]. We varied the load curtailment target from 500 to 1500 kWh. Unlike load curtailment data for which we had a real world dataset, we had to sim- ulate solar curtailment data. The output of a solar PV is determined mainly by the solar radiance at the PV installation, the PV area and the PV yield [7]. We used the hourly solar radiance data available at [5] for the Los Angeles County. We then used the PV- output calculator available at [7] to calculate the solar generation data by varying the PV area from 10m 2 to 20m 2 . We also varied the solar panel yield from 5% to 15%. Hence, a fixed PV area and yield represents a node in our dataset. To obtain solar curtailment values, if the PV output for a given hour for a node wasO, we generated 6 curtailment values: 0, 0:125O; 0:25O; 0:5O; 0:75O;O. Hence, each curtailment value repre- sents a PV connection/disconnection setting. All the 4 15-min intervals for a given hour were assigned the same curtailment values. For clarity, given an (L;U) pair, we report the error of either solar or load curtailment, whichever one performs worse. We used a precision of 3 decimal places for both load and solar curtailment values. 5.6.2 Minimum Cost Supply Demand Matching We evaluate our minimum cost supply demand matching (MinCB) Algorithm by vary- ing the (L,U) pairs as discussed above. For each (L,U) we study the amount by which the curtailment target constraint (Equation 5.2) is violated. Moreover, we also compare the cost of the solution obtained by our algorithm to the optimal cost obtained by the ILP. We also perform a scalability analysis of our algorithm. Accuracy Analysis Figures 5.1 and 5.2 show the percentage error of the curtailment target constraint (Equa- tion 5.2) violation for various values of the theoretical guarantee . For example, if 102 0.00 5.00 10.00 15.00 20.00 25.00 30.00 35.00 40.00 45.00 Percentage Error Curtailment Range (kWh) 50% 20% 10% Figure 5.1: Percentage Error of the curtailment target constraints for = 0.5-0.1 (50-10% Error Guarantee) = 0:05, the algorithm will incur an error of 5% in the worst case. As we can note from the figures, the errors incurred by our solution are within the theoretical guarantees pro- vided by the number. In practice, the errors are much lower. For = 0:5 (50%), the highest error incurred is 40%. Similarly, for = 0:2 (20%), barring a few cases, the errors are less than 15%. As evident from Theorem 1, the cost of the solutions obtained from MinCB are less than or equal to the optimal solution. This is possible because instead of tightly satisfying the constraint, as the optimal solution does, MinCB tries to find a solution with a lower cost which possibly violates the constraints by a maximum of. Figure 5.3 shows the ratio of the cost of solution obtained by MinCB and the optimal solution 103 0.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00 Percentage Error Curtailment Range (kWh) 5% 2% Figure 5.2: Percentage Error of the curtailment target constraints for = 0.05-0.02 (5-2% Error Guarantee) versus the approximation guarantee (). As we can notice from the figure, for each curtailment range, the ratio increases as the approximation guarantee is tightened i.e. reduced. For lower values of such as 0.02 (2%), the ratio is close to 1. The ratio is never greater than 1 implying that the cost of the solution obtained by MinCB is always less than the optimal cost. We can also note that for a fixed approximation guarantee (), the ratio typically decreases with an increase in the upper bound U of the curtailment range. This trend is more pronounced in higher values of such as 0.5 (50%) and 0.2 (20%). This is because a higher value of U provides a larger error range in which to search for minimum cost 104 50% 20% 10% 5% 2% 250-500 0.71 0.88 0.94 0.97 0.99 375-750 0.67 0.85 0.93 0.96 0.99 500-1000 0.62 0.84 0.92 0.95 0.98 600-1200 0.59 0.84 0.93 0.97 0.99 750-1500 0.60 0.83 0.91 0.96 0.99 0.00 0.20 0.40 0.60 0.80 1.00 1.20 Cost Ratio Approximation Guarantee Parameter 250-500 375-750 500-1000 600-1200 750-1500 Figure 5.3: Ratio of the cost of solution obtained by MinCB and the optimal solution versus the approximation guarantee () solutions. For example, for 50% error guarantee, the error range which is 750 kWh for U = 1500, is three times larger than the error range of 250 kWh for U = 500. Scalability Analysis In order to perform scalability analysis, we fix the values of T : the number of time intervals andN: the number of curtailment strategies per node. We varyM: the number of nodes for various values of. As one can see from Figure 5.4, the algorithm exhibits a near quadratic increase in the runtime with respect to the number of nodes. This is consistent with the runtime complexity analysis. 105 0 100 200 300 400 500 600 700 800 10 20 30 40 50 Runtime (seconds) Number of Nodes 50% 20% 10% Figure 5.4: Runtime of MinCB versus: (a) Number of nodes for various values of epsilon 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 50% 20% 10% 5% 2% Runtime (seconds) Approximation Guarantee (Percentage Error) M = 10 M = 20 M = 30 Figure 5.5: Runtime of MinCB versus: (b) (denoted as percentage error) for various values of the number of nodes 106 We also analyze the scalability with respect to by varying the value of M while keepingT andN fixed. As shown in Figure 5.5, decreasing (increasing accuracy) has a significant impact on runtime. Hence, is a parameter that can be used to trade-off accuracy versus computational complexity. A reader might comment that the runtimes observed, especially for smaller values of are very high. This can be justified as the problem of supply demand matching is NP-hard and hence, high computation capacity is required to increase the accuracy. Our objective in this work is to show that a polynomial time approximation algorithm exists for this NP-hard problem. We did not focus on finding the best optimal solution for the same. Moreover, our experiments are performed on MATLAB. For a real world deployment of this software, using faster programming languages such as C++ will significantly improve the run times (as high as 10-20 times as per the experience of the authors). In the context of real world scenarios, the California ISO’s MRTU applications determine the desired generation changes 5-min ahead of the beginning of the interval and the system needs to start moving towards the set point 2.5 minutes ahead of the interval [44]. Given the typical inverter control delays of the order of milliseconds [33], even our naive MATLAB implementation meets the California ISO constraints for 40 nodes with = 0:2 (20% error) and 25 nodes for = 0:1 (10% error). Similar imple- mentation in C++, assuming a conservative estimate of 10x improvement will meet the constraints for 70 nodes with = 0:1 (10% error) and 40 nodes with = 0:05 (5% error). To further improve the runtime, parallelization techniques can be used as the algo- rithm is amenable to parallelization trivially. TheT t tables can be filled in parallel. Moreover, all the entries corresponding to a nodeb in t and timet in can be filled in parallel leading to dramatic improvement in the performance. 107 We also compare our algorithm against demand curtailment selection techniques such as those developed in [76] and [34]. We observed that these techniques typically incur errors (constraint violations) of around 5-20% and in the worst case can go as high as 95%. We exclude the details of this analysis as comparison against the optimal solutions obtained by solving the ILP makes this analysis redundant. Feb-May July-October No NLB 13.47 13.68 50% 26.94 27.36 10% 134.68 136.78 5% 269.36 273.56 0.00 50.00 100.00 150.00 200.00 250.00 300.00 Percentage PV Penetration Figure 5.6: Improvement in PV Penetration 5.6.3 Improvement in PV Penetration We also analyze the improvement that can be achieved in the percentage of PV pene- tration using our minimum cost supply demand matching (MinCB) algorithm using a simple analysis. We consider a smart grid with a ratio of average load to minimum load 108 as 10 during the day light hours of 10 am to 5 pm. We analyze the solar irradiance time- series data obtained for the Los Angeles County [5] for the periods of February-May 2010 and July-October 2010 to calculate the normalized mean and standard deviation. We define PV Penetration Percentage as PVavg100 Loadavg , wherePV avg is the average supply from solar and Load avg is the average load. To calculate PV Penetration Percentage, we assume that the maximum supply i.e. average supply plus one standard deviation, minus any curtailment should be less than the minimum load. Moreover, we assume that the maximum curtailment that can be performed is 1, where is the accuracy parameter chosen. For example, if = 0:5 i.e. 50% then maximum curtailment that can be performed is 50% of the maximum supply. This is because if the curtailment target is the entire maximum supply, the algorithm can guarantee only 50% curtailment for = 0:5. As shown in Figure 5.6, in the absence of supply demand matching (net load balancing) algorithm (labeled No NLB), the maximum PV penetration is around 13%. This value dramatically increases with the use of supply demand matching framework and can be as high as 270%. The analysis above suggests that by enabling PV curtailment with high accuracy, we can achieve the following objectives: We can ensure that the maximum supply in any time step, assumed to be average supply plus one standard deviation minus supply curtailment, does not exceed the minimum possible load value. A high PV penetration value can be achieved, thereby, reducing the amount of supply from external markets during peak demand periods. For example, a PV penetration value of 270% implies that at any given time step, if the PV is supply- ing its average value and even if the demand is 2.7 times the average load value, the PV supply will be enough to meet the demand and no external supply from the market will be required. 109 5.6.4 Minimum Cost Supply Demand Matching with Fairness In Section 5.6.4, we introduced the notion of fairness by defining curtailment budget ranges for nodes. Here we evaluate the empirical performance of our minimum cost sup- ply demand matching with fairness (MinCBF) algorithm. Note that the budget ranges for each node can be set appropriately by user/grid operator. In our experiments, we set the budgetB b for each nodeb to b max = P B b=1 b max where b max denotes the sum of the maximum curtailment values across all the intervals. We also set b = and vary the value. 250-500 375-750 500-1000 600-1200 750-1500 666.67-1000 833.33-1000 800-1200 1000-1200 -1.00% -0.50% 0.00% 0.50% 1.00% 1.50% 2.00% 2.50% Percentage Error of the Cost of Algorithm 3 w.r.t. the Optimal Cost Curtailment Range (kWh) Percentage Error α = 0 Figure 5.7: Percentage Error of the Cost of MinCBF w.r.t. the Optimal Cost In order to evaluate the accuracy of our algorithm, we compare against the optimal solution obtained by solving the ILP defined in Section 5.2.1. Figure 5.7 shows the percentage error of the cost of the solution produced by MinCBF as compared against the optimal cost. Although, the worst case theoretical guarantee is a factor of 4 for quadratic cost function as provided by Theorem 5.3, in practice the algorithm performs 110 much better with errors varying from -0.79% to 1.88%. The negative percentage error implies that the cost of the solution from our algorithm was less than the cost of the opti- mal solution. Note that this is possible because our solution violates certain constraints which the optimal solution does not. However, the complexity of the ILP prevented CPLEX studio to generate results in a reasonable amount of time. Hence, we compared the objective value of the rounded results obtained from our algorithm against the objective value of the LP with fractional solutions. Since, an LP relaxation of an ILP always produces better objective value (albeit with infeasible integral solutions), the results shown will only improve if com- pared against the optimal ILP solution. 1 2 3 4 5 6 7 8 9 10 0.00% 2.00% 4.00% 6.00% 8.00% 10.00% 12.00% 14.00% sub-title Percentage of Budget Overshot Nodes Percentage of Budget Overshoot Figure 5.8: Percentage of Budget Overshoot for 10 Worst Affected Nodes Figure 5.8 shows the percentage of the budget overshot for the top 10 worst affected nodes across all the values of L and U. The theoretical guarantee of factor 2 (100%) 111 1 2 3 4 5 6 7 8 9 10 0% 1% 2% 3% 4% 5% 6% 7% Percentage Undershoot of Interval Target Time Interval Percentage Error Figure 5.9: Percentage Undershoot of Interval Curtailment Target for various intervals (note: the labels of x axis do not denote the actual interval number) provided by Theorem 5.2 is honored in all the cases. Moreover, in practice the error is less than 13% (1.13 factor) for the worst performing node. Figure 5.9 shows the percentage of the interval curtailment target t undershot for the top 10 worst affected intervals across all the values of L and U. The theoretical guarantee of factor 2 (100%) provided by Theorem 5.2 is honored in all the cases. As we can see from the figure, in practice the error is less than 7% (1.07 factor) for the worst performing interval. This implies that the curtailment target for the worst performing interval could not be met and was deficit by 7%. We also calculated the gini coefficient – which is the most commonly used measure of inequality in economics – of the curtailment achieved by each node as a proportion of its budget to measure the fairness of curtailment. Figure 5.10 shows the results for various values of for a curtailment range of 500-1000 with 20 nodes. The value of 112 0.16 0.17 0.18 0.19 0.2 0.21 0.22 0.05 0.1 0.15 0.2 Gini Alpha Gini Figure 5.10: Gini Coefficient for various values of Alpha () gini decreases with increasing as the dispersion of the curtailment decreases. Above = 0:2, no feasible solution could be found. Mathematically, G = P n i=1 P n j=1 jx i x j j 2n P n i=1 x i (5.28) 5.6.5 Online Algorithm for Minimum Cost Supply Demand Match- ing with Fairness In order to evaluate our online algorithm for minimum cost supply demand matching with fairness (Online MinCBF), we compare the cost of the solutions obtained against the optimal solutions of the ILP defined in Section 5.2.1 as generated in Section 5.6.4. For various values of L and U pairs, we calculate the percentage error in the cost 113 obtained by the online algorithm (Algorithm 8) with respect to the optimal solutions. The budget valuesB b input to Algorithm 8 are same as the ones used in Section 5.6.4. 166.67-250 208.33-250 333.33-500 416.66-500 500-750 625-750 666.66-1000 833.33-1000 800-1200 1000-1200 0% 5% 10% 15% 20% 25% Percentage Error of online algorithm w.r.t. off- line algorithm Curtailment Range (kWh) Percentage Error in Cost Figure 5.11: Percentage Error of the Cost of Online MinCBF w.r.t the Optimal Cost Figure 5.11 shows the results obtained. Even though we do not provide any guar- antee on the worst case bounds for the online algorithm, in practice the error incurred is low. The highest error incurred is around 23% (factor 1.23). This makes the online algorithm a good candidate for supply demand matching when the predictions for the entire horizon are not known in advance. 5.6.6 Supply Demand Matching with Network Constraints We evaluate Algorithms 9 and 10 by considering a smart grid with 10 transformers each connected to 10 nodes. We consider a single time interval in which 5 of the transformers have a load surplus and 5 have solar surplus. For each transformer with solar surplus, we assume that the aggregate load connected is 0.25 times its total solar generation 114 0.00 0.20 0.40 0.60 0.80 1.00 1.20 1 0.75 0.5 0.25 0 Cost Ratio Feeder Capacity 0.75 0.5 0.25 Figure 5.12: Ratio of Cost of Algorithm 10 to Algorithm 9 versus Feeder Capacity for various values of Distributor capacity. S tx . We set the capacity of distributor for each transformerCap tx equal tof 1 S tx , where f 1 2 [0; 1] and vary f 1 . We set the capacity of the corresponding feeder Cap f equal to f 2 f 1 S tx , where f 2 2 [0; 1] and vary f 2 . Figure 5.12 shows the ratio of the cost of solution obtained by Algorithm 10 to that of Algorithm 9 versus the feeder capacity, which is determined by varying f 2 . Each line in the figure corresponds to a single distributor capacity as determined byf 1 . As evident from the figure, for a fixed distributor capacity, a reduction in the feeder capacity leads to a reduction in the cost ratio i.e. cost of Algorithm 10 and Algorithm 9 converge with reduction in the feeder capacity until they are equal for a feeder capacity of 0. Moreover, for a fixed feeder capacity, the gap between the costs reduces with a reduction in the distributor capacity. This is because reduced distributor capacity results in lower solar surplus and hence less surplus redistribution. Note that when the distrib- utor capacity is f 1 = 0:25, the costs are equal. This is because whatever solar power 115 is being produced in a distributor, it is being consumed by the load (which is equal to 0.25 times the solar generation). Hence, the cost incurred by both the algorithms is for curtailing loads in the load surplus transformers. 5.6.7 Storage Modeling We also evaluate the percentage increase in runtime due to incorporation of storage nodes into our algorithm. We fix the curtailment range as (667,1000) kWh and the number of total nodes as 20. We assume an infinitely sized storage with charge/discharge capacity of 1000 kWh and charge/discharge efficiency of 1 in each time interval. We then vary the number of storage nodes for various values of percentage errors as denoted by and calculate the percentage increase in runtime compared against a setup in which no storage nodes are used. Figure 5.13 shows the detailed results. Even for 15 storage nodes and = 0:02 i.e. 2%, the increase in runtime is around 160% i.e. the increased runtime is less than 2.6 times the runtime of the algorithm without any storage nodes. 116 0 20 40 60 80 100 120 140 160 180 1 2 5 10 15 Percentage Increase in Runtime Number of Storage Nodes 10% 5% 2% Figure 5.13: Percentage Increase in Runtime For various number of storage nodes for varying epsilon. 117 Chapter 6 Discrete Supply Demand Matching Under Prediction Uncertainty Previous chapters discussed algorithms for performing discrete supply demand match- ing assuming no errors in prediction data of the operational values of each nodes. How- ever, in practice predictions are uncertain which can lead to errors in the solutions. In this chapter, we explicitly model prediction uncertainty and develop approximation algorithm to perform discrete supply demand matching accounting for the uncertainty. 6.1 Two Stage Stochastic Recourse Model for Discrete Supply Demand Matching We consider a model similar to two stage stochastic recourse model [55] to perform dis- crete supply demand matching while minimizing the expected uncertainty due to errors in the prediction models. We first define a few notations to supplement our smart grid model and then delve into the problem formulation and the solution using an approxi- mation algorithm . 6.1.1 Smart Grid Model As before, our smart grid consists ofM s supply nodes andM d demand nodes each of them withN strategies. The horizon of planning consists ofT intervals. In each interval 118 t2 [T ] where, for notational simplicity, [a] denotes the setf1; 2;:::;ag, the supply node i s (demand nodei d ) can be in exactly one operational statej with a predicted discrete supply (load) value of s isj (t)( d i d j (t)). In this chapter, we consider optimization in a single interval and hence we drop (t) to simplify notations. Moreover, s isj ( d i d j ) denotes the actual generation or consumption for statej instead of the curtailment value as done in the previous chapters. The corresponding utility of the operational state isu s isj (u d i d j ). We assume thatD is the total uncontrollable load in the system. This the minimum load which needs to be matched and hence we assign no utility to it. We use the 0-1 matrix X s (X d ) to represent the decision variables in whichx s isj (x d i d j ) is set to 1 if nodei s (i d ) is in operational statej and 0 otherwise. 6.1.2 Supply Demand Matching without Prediction Uncertainty The grid operator needs to ensure that in any given interval, the total supply in the grid is equal to the total demand while maximizing the utility. This problem can be formulated using the following Integer Program (IP): max X s ;X d U(X s ;X d ) = Ms X is=1 X j2[N] u s isj x s isj + M d X i d =1 X j2O i d u d i d j x d i d j (6.1) s:t: Ms X is=1 X j2O is s isj x s isj = M d X i d =1 X j2O i d d i d j x d i d j +D (6.2) X j2[N] x s isj = 1 8i s (6.3) X j2[N] x d i d j = 1 8i d (6.4) 119 In the above formulation, Equation 6.2 ensures that supply and demand are matched. Equations 6.3 and 6.4 ensure that a node is in exactly one operational state. We assume that the IP is always feasible. 6.1.3 Incorporating Prediction Uncertainty in Supply Demand Matching To model the prediction uncertainty, each discrete supply (load) value is associated with a random variable s isj ( d i d j ), which denotes the error in prediction i.e. the actual real- time supply (load) value is s isj + s isj ( d i d j + d i d j ). The incorporation of the random parameters into our problem coverts the determin- istic problem of Section 6.1.2 into a stochastic one. We employ two state recourse stochastic programming model to make decisions under uncertainty [55]. In a two stage recourse model, a first stage decision is made and a (higher cost) recourse action is taken in second stage after the realizations of the random variables based on the first state decisions. Hence, the first stage decision needs to be made with two objectives: 1) Maximize the utility of implementing first stage decisions, and 2) minimize the cost associated with any recourse action taken in the second stage for all possible realizations of the random variables. We define a scenario as a possible realization of all the random variables. LetS denote the set of all possible realizations. LetK =jSj denote the number of scenarios. For a random variable s isj ( d i d j ), let si isj ( di i d j ) denote the realized value under scenario i2S. Letp i be the probability of the realization of scenarioi. For a first stage deci- sionX s ;X d obtained by solving the optimization problem (6.1), the realization of a scenario i can lead to violation of the constraint 6.2 i.e. a imbalance between supply and demand can occur. The aggregate net load imbalance (load minus supply) is given using: 120 v i (X s ;X d ) = M d X i d =1 X j2O i d di i d j x d i d j Ms X is=1 X j2O is gi igj x g igj (6.5) Hence, in the second stagev i (X s ;X d ) amount of supply needs to injected into the grid (or load increased if negative). We assume this is achieved using battery storage. We assume the cost of the recourse action is given using q(v i (X s ;X d )) or simply q(v i ). We can formulate the two stage recourse model for supply demand matching as the following bi-criteria optimization: 2-Stage Recourse (2-SR): max X s ;X d fU(X s ;X d )g; min X s ;X d f K X i=1 p i q(jv i j)g (6.6) s:t: v i = M d X i d =1 X j2[N] di i d j x d i d j Ms X is=1 X j2[N] si isj x s isj 8i = 1;:::;K (6.7) (6:2) (6:4) The objective of the two stage recourse model is to maximize the utility of the deci- sions taken in stage one and minimize the expected value of the cost incurred due to recourse action. There are two significant challenges in efficiency solving 2-SR: 1) The decision vari- ablesX s ;X d are discrete binary variables. As this problem is an NP-hard problem, no algorithm can exist which has a polynomial complexity runtime in the input parameters and which provides the optimal solution unless P == NP, and 2) the valueK can be very large making the algorithm definition itself vary large. The algorithms developed in the next section address both the challenges. The first challenge is addressed by developing 121 an approximation algorithm which in polynomial time outputs a solution whose value, even in the worst case scenario, is ”close” to the optimal solution. The second challenge is addressed by making certain assumptions on the randomness of the scenarios. 6.2 Approximation Algorithms for Supply Demand Matching Under Prediction Uncertainty We consider two models for the realization of the random variables and develop approx- imation algorithms for the same. 6.2.1 Latent State Model This model assumes the existence of a set of latent statesH. Each latent state h i 2 H is realized with a probability p i and is responsible for the realization of scenario i2S. In other words, all the random load and generation variables are observations corresponding to these latent states. Now, the number of latent states can be infinite. However, using the expectation valueE[ s isj ](E[ d i d j ]) for each random variable s isj ( d i d j ), we can reformulate the 2-SR optimization problem. Assuming a linear cost functionq, Equation 6.6 can be written as: max X s ;X d fU(X s ;X d )g; min X s ;X d f K X i=1 p i (j M d X i d =1 X j2[N] q( di i d j )x d i d j Ns X is=1 X j2[N] q( si isj )x s isj j)g (6.8) 122 switching the summation order and moving the probabilityp inside the cost function we get: max X s ;X d fU(X g ;X d )g; min X s ;X d fj N d X i d =1 X j2[N] q( K X i=1 p i di i d j )x d i d j Ns X is=1 X j2[N] q( K X i=1 p i si isj )x s isj jg (6.9) Replacing the terms P p i i with the expected value and ignoring K, we have the following formulation: 2-Stage Recourse for Latent state model (2-SRL): max X s ;X d fU(X s ;X d )g; min X s ;X d fQ(X s ;X d )g (6.10) s:t:(6:2) (6:4) whereQ(X s ;X d ) =j P M d i d =1 P j2[N] q(E[ d i d j ])x d i d j P Mg is=1 P j2[N] q(E[ s isj ])x s isj j. To simplify notations, the expected cost for node i s (i d ) with operational state j: q(E[ s isj ]) (q(E[ d i d j ])) will be denoted using ~ q s isj (~ q d i d j ). 6.2.2 Approximation Algorithm for 2-SRL In order to develop an approximation algorithm for 2-SRL, 2-SRL-Approx, we will cre- ate two dynamic programs, one for the generation nodes and one for the demand nodes. For each possible generation (demand) value, utility value and uncertainty value triplet, the dynamic program will return true or false denoting whether it is possible to achieve the triplet using all the generation (demand) nodes or not. Using the approximation technique of rounding, we will ensure that the size and the runtime complexity of the dynamic program (and hence the algorithm) is polynomial in the input size and 1 , where 123 is the accuracy parameter. However, this will lead to a loss in accuracy which we will prove that in the worst case is bounded by. We will discuss the dynamic program for the demand nodes. The dynamic program for the generation nodes can be developed analogously. Let t be number such that ~ q 0 d i d j = ~ q d i d j +t > 08i d ;j. We refer to ~ q 0 as translated uncertainty values. As this translation translates the objective function by a constant amount, it does not affect the optimization problem. Let max a bound on the maximum possible demand or supply value, u max be a bound on the maximum possible utility value and ~ q 0 max be a bound on the maximum possible translated uncertainty value that can be achieved using all the demand nodes. Also, let min ;u min ; ~ q 0 min be the minimum possible values greater than zero. We assume the quantities max min ; umax u min and ~ q 0 max ~ q 0 min to be bounded by polynomial in the input size and 1 . LetM = maxfM s ;M d g We create the following partitions: [0]; [ min ; min (1 +) 1 M ); [ min (1 +) 1 M ; min (1 +) 2 M );:::; [ min (1 +) k d 1 M ; min (1 +) k d M ) [0]; [u min ;u min (1 +) 1 M ); [u min (1 +) 1 M ;u min (1 +) 2 M );:::; [u min (1 +) k d 1 M ;u min (1 +) ku M ) [0]; [ ~ q 0 min ; ~ q 0 min (1 +) 1 M ); [ ~ q 0 min (1 +) 1 M ; ~ q 0 min (1 +) 2 M );:::; [ ~ q 0 min (1 +) k d 1 M ; ~ q 0 min (1 +) kq M ) where k d = dM log 1+ ( max min )e, k u = dM log 1+ ( umax u min )e and k q = dM log 1+ ( ~ q 0 max ~ q 0 min )e. The intervals in the partitions are indexed using b 2 f0; min ;:::; min (1 + ) k d 1 M g , b u 2 f0;u min ;:::;u min (1 + ) ku1 M g and 0; ~ q 0 min ;:::; ~ q 0 min (1 +) kq1 M respectively, i.e. the lowest value of the intervals. Without loss of generality, we assume that the demand nodes are ordered as 1; 2;:::;M d . We now define a dynamic programming function d (b ;b u; b q 0 ;i d ) which returns true (1) if there exists a triplet ;u;q 0 such that each element ,u andq 0 of the triplet belongs to 124 the intervals indexed byb ,b u and b q 0 , respectively and the triplet values can be achieved using demand nodes 1;:::;i d . The function can be defined recursively as follows: d (b ;b u;b q 0 ;i d ) = 8 > > > > > > > > > < > > > > > > > > > : 1 ifi d = 1 and9j2O i d j i d j 2b ;u i d j 2b u; ~ q i d j 2 b q 0 1 if9b ;b u ;b q 0 ;j2O i d j d (b ;b u ;b q 0 ;i d 1) = 1 and i d j +b 2b ; u i d j +b u 2b u; ~ q i d j +b q 0 2b q 0 0 otherwise (6.11) Now, all the 1 entries in d (b ;b u;b q 0 ;M d ) denote the feasibleb ;b u;b q 0 triplet for demand nodes. We can obtain similar triplet for the generation nodes. To distinguish them, we use the notationsb d ;b u d ;b q 0 d andb g ;b u g ;b q 0 g for demand and generation nodes respectively. We define a function (b ) which outputs a list of the utility and the asso- ciated uncertainty for each value ofb . can be populated using Algorithm 12 . is used to determine the dominating solutions using Algorithm 13. A simple dynamic pro- gramming table traversal algorithm can then be used to determine the operational states to be followed by each node for any given dominating solution from the pareto frontier. We have the following theorem for 2-SRL-Approx. Theorem 12. 2-SRL-Approx is an algorithm with runtime complexity polynomial in the input size and 1 , where is an accuracy parameter, that produces a pareto frontier for discrete supply demand matching with uncertainty such that for each pareto solution, the utility is within (1 ) of the corresponding solution, and the maximum supply demand mismatch is within factor of both actual supply and demand. Additionally, the uncertainty associated with either supply or demand are within (1 +) of the optimal supply or demand uncertainty. 125 Algorithm 12: Populating using d and s Input: d ; s ;b 1 foreachb do 2 L d allb u d ;b q 0 d s.t. d (b ;b u d ;b q 0 d ;M d ) = 1 3 L s allb u s ;b q 0 s s.t. s (b ;b u s ;b q 0 s ;M s ) = 1 4 L b 5 foreachb u d ;b q 0 d 2L d do 6 foreachb u s ;b q 0 s 2L s do 7 L b L b [ (b u d +b u s ;jb q 0 d b q 0 s j) 8 end 9 end 10 (b ) L b 11 end Output: Proof. Correctness: We will prove the correctness of the algorithm in two steps. In step 1, we will prove that the dynamic program runs correctly to fill ( d and s ) and step 2, we will prove that if the dynamic program runs on actual values instead of the partition indices and an actual value can be achieved, then this information is not missed in dynamic program run on partition indices. We will prove step 1 using induction. Essentially, we need to prove that (b ;b u;b q 0 ;i) = 1 iff it is possible to achieveb ;b u;b q 0 using first i nodes. For i = 1, it is true trivially, as we will set (b ;b u;b q 0 ; 1) = 1 iff there exists a j in [N] such that i d j 2 b ;u i d j 2 b u; ~ q i d j 2 b q 0 . Let the inductive hypothesis be true for i = k. For k + 1, (b ;b u;b q 0 ;k + 1) = 1 if and only if (b ;b u;b q 0 ;k) = 1 and ij +b 2b ; u ij +b u 2b u; ~ q ij +b q 0 2b q 0 . So, if (b ;b u;b q 0 ;k + 1) = 1, we can achieveb ;b u;b q 0 by achievingb ;b u ;b q 0 usingk nodes and ij ;u ij ; ~ q ij by followingj th operational state of (k + 1) st node. If no suchj exists or (b ;b u;b q 0 ;k) = 0, it cannot be achieved. To proves step 2, we show that if usingk nodes, an actual value of a = P k i=1 i (proof is analogous for u and ~ q) can be achieved, thenb P k i=1 i < (1 +) k M b 126 Algorithm 13: Algorithm to generate dominating solutions using Input: 1 L allb u;b q 0 s.t. (b )6= 8b 2 SortL in descending order usingb u as key, if multiple entries haveb u, keep the one with lowestb q 0 and remove all other entries 3 LetL 0 be the entry with largestb u andL K be the entry with smallestb u 4 D L 0 5 fori = 1 toK do 6 b u;b q 0 L i 7 foreachb u l ;b q 0 l l2D do 8 ifb q 0 b q 0 l then 9 break from for each loop 10 end 11 end 12 D D[L i 13 end Output:D: set of dominating solutions and some entries (b : (1 +) k M b ; :; :;k) will be set to 1. We again use induction. For k = 1, if 2 b can be achieved then some entries of (b ; :; :; 1) are true and b < (1 +) 1 M b . Let the inductive hypothesis be true for k 0 i.e. if a value of 0 k = P k 0 i=1 i can be achieved by node i exhibiting an operational value of i , then some entries of (b 0 k : (1 +) k 0 M b 0 k ; :; :;k 0 ) are true andb 0 k P k 0 i=1 i < (1 +) k 0 M b 0 k . Fork 0 + 1, as per the definition of the dynamic program, the entry whereb 0 k + k 0 +1 falls will be true. So, we just need to show thatb k 0 +1 k 0 +1 = P k 0 +1 i=1 i < (1 +) k 0 M b k 0 +1 . We know the following: (1)b k 0 +1 b k + k 0 +1 < (1 +) 1 M b k 0 +1 , and (2) P k 0 i=1 i < (1+) k 0 M b 0 k . So, P k 0 +1 i=1 i < (1+) k 0 M b 0 k + k 0 +1 < (1+) k 0 M (b 0 k + k 0 +1 )< (1+) k 0 M (1+ ) 1 M b k 0 +1 < (1 +) k 0 +1 M b k 0 +1 as (1 +) k 0 M > 1, when;k 0 > 0. Runtime: The dynamic programming table haskk q k u M entries. Each entry requireskk q k u N time to fill. Hence, the total runtime isk 2 k 2 q k 2 u MN. Since, k d =dM log 1+ ( max min )e, k u =dM log 1+ ( umax u min )e andk q =dM log 1+ ( ~ q 0 max ~ q 0 min )e, the runtime is polynomial in the input size and 1 . 127 Approximation Guarantee: Consider a dominating solutionb ;b u;b q 0 . Let u be actual reward value that can be achieved and belongs to the partitionb u. It is evident from the correctness proof above that such au will always exist. Now, as shown in step 2 of the correctness proof,b u u < (1 +) k M b u fork nodes. So, forM s ;M d M nodes this implies,b u > 1 1+ u = 1 1 2 u. Hence,b u (1)u. Similarly, for uncertainty we can prove thatb q 0 (1 +)q 0 . The above proof establishes that d and s are correct, can be filled in polynomial time and the solution values are within (1 +) of the optimal. d and s are used to fill . Since utilities are a simple summation and uncertainty is simple absolute different, correctness of follows. Algorithm 12 has polynomial time complexity. Now, b u d (1)u d andb u s (1)u s impliesb u d +b u s (1)(u d +u s ). As the final uncertainty is obtained using subtraction, we cannot provide a guarantee for the same. Now, consider the supply and demand partition indicesb s =b d =b . Let s and d be the actual generation and demand values.b s ; d < (1 +)b . So, s d < (1 +)b b b s ; s . Hence, the mismatch in generation and demand is within factor of both actual generation and consumption. Discussion on translation of uncertainty and upper bounds max ,u max and ~ q 0 max The uncertainty value ~ q ij associated with a node i and operational state j can have a negative value. Thus it cannot be directly used for indexing the dynamic programming table . Hence, we add a value t =j minf0; min ij f~ q ij ggj to ~ q ij for each i;j to pro- duce ~ q 0 ij which are non-negative. The following relationship holds: P M i=1 P j2[N] ~ q 0 ij = P M i=1 P j2[N] ~ q ij +Mt, where M d is the number of generation nodes M s or demand nodesM d . 128 The values max ,u max and ~ q 0 max denote the maximum possible value of supply (or demand) , utility u and translated uncertainty ~ q 0 respectively. These values have an impact on the number of entries in the dynamic programming table . In this work, we set the values as follows: max = P M i=1 max j2[N] ij , u max = P M i=1 max j2[N] u ij and ~ q 0 max = P M i=1 max j2[N] ~ q 0 ij i.e. we take the sum of maximum possible value that each node can exhibit. We further make the following assumptions: max = poly(M; 1 ; min ), u max = poly(M; 1 ;u min ) and ~ q 0 max =poly(M; 1 ; ~ q 0 min ). This assumption implies the the min- imum possible non zero value exhibited by any node is not ”too small” compared to the maximum achievable value from all the nodes. For the nodes in a smart grid, this is reasonable assumption to make. 6.2.3 Black Box Model We now relax the assumption that all the random load and generation variables are observations corresponding to a set of latent states. The implication here is that each random variable potentially has its own state governing its realization and hence the number of possible scenarios in the system explodes combinatorially. In this case, it is not even possible to define the problem in polynomial time (due to exponential number of constraints of the form (6.7)). Therefore, we consider a black box model to study such a problem. The analysis here is a minor modification of the analysis performed in [21]. Under this model, we are given a black box, which can be used to sample scenarios, where each scenario corresponds to a particular realization of all the random variables. Sampling of a single scenario is assumed to be anO(1) time complexity operation. We show that by sampling only polynomial number of scenarios, we can obtain a solution which is within 1+O() of the optimal solution. This technique is also known as Sample Average Approximation 129 (SAA) [59]. The techniques simply drawsK number of scenarios and assumes that each scenario occurs with a probability of 1 K to populate the expected values for defining the problem. We first define some notations and assumptions. Let the optimal 2-SR formula- tion, i.e. the 2-SR formulation which considers the exponential number of scenarios be defined as: maxfU(x)g; minfE[v(x)]g s:t:x2X (6.12) whereX denotes the feasible range of decision variablesX s ;X d ,U(x) denotes the utility for a feasible solutionx andE[v(x)] denotes the expected value of the uncertainty forx. Now, let the 2-SR model under blackbox model i.e. the formulation obtained by samplingK scenarios, whereK is a polynomial in the input size be defined as: maxf b U(x)g; minf 1 K K X i=1 v i (x)g s:t:x2X (6.13) We make the following assumptions: 1. X d =0 i.e. load is uncontrollable. 2. The feasible range ofx2 X is such that the predicted supply meets the demand in the first stage i.e. P Ms is=1 P j2O is s isj x s isj =D 3. Let Z be the minimum value offE[v(x)]g8x2 X. Let x be the minimizer. We assume: 0 v i (x) Z 8x;i, where ; > 0. This assumption implies 130 that the minimum (non zero) amount of storage required in second stage under any first stage decision with any possible scenario realization is not ”too small” compared to the maximum amount of storage required in the second stage. Note that 0v i (x) can be achieved by translating all the uncertainty values by a fixed amount. Also note thatZ > 0 as it is a value corresponding to expectation value of non-negative terms. Now, we want to study the error introduced due the use of sample mean instead the expected value:jU(x) b U(x)j andjE[v(x)] 1 K P K i=1 v i (x)j to compare the solutions (6.12) and (6.13). AsU(x) is independent of the realized scenario,U(x) = b U(x). So, we just focus onjE[v(x)] 1 K P K i=1 v i (x)j8x2X. Forx2X, we defineY i = v i (x) Z 8i. It follows from assumption (3) thatY i 2 [0; 1]. For notational simplicity we dropx and write the terms asE[v] and 1 K P K i=1 v i . 6.2.4 Approximation Algorithm for 2-SRB 2-SRB-Approx: In order to develop an approximation algorithm for 2-SRB, we will sam- pleK = ( 2 4 logjXj log 1 ) scenarios and calculate the expected values 1 K P K i=1 for all the random variables. We will then use 2-SRL-Approx to generate the dominating solutions. Now, as per Chernoff’s Bound [21], Given independent random variablesY 1 ;:::;Y K withY i 2 [0; 1]8i, andY = P K i=1 Y i for every> 0: PrfjYE[Y ]j>Kg 2e 2 K (6.14) 131 Now, note that E[Y ] = E[ P K i=1 Y i ] = K Z E[ 1 K P K i=1 v i ] = K Z E[v]. So, PrfjE[v] 1 K P K i=1 v i j>Z g =Prfj Z K E[Y ] Z K K X i=1 Y i j>Z g (6.15) =PrfjE[ K X i=1 Y i ] K X i=1 Y i j> 2 K g (6.16) By equation 6.14 and substituting the value of K = ( 2 4 logjXj log 1 ), this probability is less than jXj . Taking union bound over allx2X, with probability 1, 8x2X: jE[v(x)] 1 K K X i=1 v i (x)jZ (6.17) Now, letE 0 [v(x)] = 1 K P K i=1 v i (x). Letx 0 be the minimizer ofE 0 [v(x)]. Note that x 0 is not an exact minimizer and the approximation algorithm 2-SRL is used to obtain x 0 . Using Equation 6.17 withx , we getE 0 [v(x )]Z Z . Similarly, withx 0 , we getE[v(x 0 )] Z +E 0 [v(x 0 )]. This implies, E[v(x 0 )] Z + (1 +)E 0 [v(x )] Z + (1 +) 2 Z . SoE[v(x 0 )] (1 +)(2 +)Z . As each sampling is assumed to be anO(1) operation, we have the following theorem Theorem 13. 2-SRB-Approx, with probability (1), gives a (1 +)(2 +)-factor bound on uncertainty for 2-SRB with runtime polynomial in, 1 , 1 and the input size. 6.2.5 Practical Implementation 2-SRL-Approx and 2-SRB-Approx algorithms, while providing both worst case error and polynomial runtime guarantee might not scale very well for large problem sizes due to the large dimensionality of the dynamic program table of the orderkk q k u M. 132 In this section, we discuss a few techniques to improve the runtime performance of the algorithm. Note that the filling of the dynamic program table is highly parallelizable, as for a node, each entry can be filled independently just using the values of the previous node. However, the large size of the table can quickly fill up the available memory on the system. This can be addressed using the following techniques: 1. When the entries corresponding to node i is being filled, only the entries corre- sponding to node i 1 are needed. Entries for any prior node can be deleted thereby reducing the size of the table by a factor of M 2 . 2. The algorithm can be used in single mode similar to algorithms developed in Section 5.1. This will allow us to use the following memory optimizations: (a) The dimension corresponding to uncertainty can be removed. Instead, each entry of the table (b ;b u;i)can now store the minimum possible uncertainty value that can be achieved by nodes 1;:::;i at utilityb u with operational valueb . The subtraction of demand and supply uncertainties to produce the final uncertainty value does not allow this optimization in dual mode. (b) The dimension corresponding to operational values can be removed if we assume that the constraints governing the aggregate operational value from all the nodes are of one of the following two forms (but not both): (i) the aggregate operational value is greater than a threshold, (ii) the aggregate operational value is below some threshold. More sophisticated techniques for runtime performance improvement is beyond the scope of this work as our objective is to prove the existence of polynomial runtime com- plexity approximation algorithm for discrete supply demand matching with prediction uncertainty and not necessarily the development of a runtime optimized framework. 133 6.3 Risk Aware Supply Demand Matching in Smart Grids with High DER Penetration We now consider an alternative approach to incorporate prediction uncertainty into deci- sion making. In this approach, a storage schedule is developed to minimize the expected values of the cost of the decisions over a horizon. Moreover, the loads and local genera- tion using DERs such as PVs are assumed to be uncontrollable and storage is assumed to be the only controllable DER. Moreover, this approach explicitly considers the presence of an external electricity market which was not considered in the previous models. 6.3.1 Smart Grid Model for Risk Aware Supply Demand Matching The smart grid that we consider in this approach is defined as a micro grid that best represents a University or Industrial campus. It is characterized by the presence of sev- eral load consuming nodes such as buildings, and several producers of non-commercial scale (cheap) electricity. In this work, we assume that the producers supply electric- ity to the micro grid using PV installations over the buildings. The micro grid con- sists of distributed storage with an aggregate capacityR. The storage has a cumulative charge/discharge capacity of c. We do not model storage specific parameters such as Depth of Discharge (DoD), efficiency as we assume the presence of a number of dis- tributed storage systems rather than a single battery system. A centralized controller is responsible for ensuring smooth grid operations over a decision making horizon defined over a set of daytime intervals t2f1;:::;Tg. Note that the horizon is defined over daytime since supply within the micro grid using PVs is only available during that time. We assume that during the night time, all demand needs to be met using the storage accumulated during the day. In our framework, the controller responds to per interval 134 supply demand mismatch by either procuring/selling electricity from the external mar- ket or by charging/discharging the storage. For each intervalt, a buy actionb t needs to be performed, whereb t is the amount of electricity bought from the market (in kWh). Moreover, a discharge actionsr t is performed, wheresr t is the amount of the electricity discharged from the storage (in kWh). If b t and sr t are negative, opposite actions of selling the electricity or charging the storage, respectively are performed. The action performed att is represented usinga t =< b t ;sr t >. The cost of buying (selling) elec- tricity at intervalt isC b t . We do not associate any cost with storage charging/discharging. We also assume that the set of costsC b t ; 0tT is known a priori. 6.3.2 Input Data Model The inputs to our model are the load and generation predictions. As the predictions are prone to errors, they need to be modeled carefully so that the framework can make decision accounting for the errors. We assume that at the beginning of each interval t of the decision making horizon, the controller receives an input data vector Y t =< D t ;D t+1 ;:::;D T >. Here D t 0 = L t 0P t 0; t 0 2ft;:::;Tg denotes the per interval supply demand mismatch; L t 0 and P t 0 denote the per interval aggregate consumption and PV generation respectively. At timet, D t denotes the actual supply demand mismatch (net load imbalance) as observed in real time whereas theD t 0; t 0 2ft+1;:::;Tg are random variables denoting the predicted imbalances. The relation between actiona t andD t isD t =b t +sr t . Thus, if we want to buy excess storage for the future at timet, we can chooseb t >D t . 6.3.3 Tail End Risk At the end of the decision making horizon, we assume that the minimum amount of storage required for night time operations is R 0 . Hence, the tail end event whose risk 135 needs to be minimized is the unavailability of storage capacity of more thanR 0 by the end of the decision making horizon. Depending upon the severity of the implications of failing to avoid the risk, the grid operator can associate it with a cost function as follows: C(R T ;R 0 ) = 8 > > < > > : 0 R T R 0 C risk otherwise (6.18) where R T denotes the storage available after the interval T and C risk denotes the cost of not avoiding the risk. C risk needs to be calculated appropriately to provide the operator with a choice to meet or not to meetR 0 for a given input data. In this work, we assume that the risk needs to be avoided at all cost and henceC risk =1. We will extend our framework to handle partial risk tolerance in our future works. 6.3.4 Supply Demand Framework A high level overview of our supply demand matching framework is shown in Fig- ure 6.1. A storage capacity ofR 0 is available at the beginning of the decision making horizon. In each time intervalt, using the most updated input data Y t , the framework makes a decision and outputs action a t for buying (selling) from the external market and discharging (charging) the storage. The actiona t for each interval is determined by formulating and solving a Markov Decision Process (MDP) [19] as described in the next section. MDP is a widely used mathematical framework for solving sequential decision making problems. At the end of the horizon, the amount of storage availableR T should beR 0 . 136 0 ′ - - - 1 1 … 1 2 2 2 ≤ Figure 6.1: Supply Demand Matching Framework 6.3.5 Risk Aware Supply Demand Matching Framework In this section, we discuss the details of the Risk Aware Supply Demand Matching Framework that we developed to perform minimum cost per interval supply demand matching using storage and electricity market while ensuring that the tail end risk of the shortfall of storage is avoided. Objective Given the supply demand matching model defined in Section 6.3.1, the objective of our framework is to minimize the cost of grid operations while ensuring that in each interval, the power ingress is equal to the power egress and at the end of the horizon, the minimum storage requirement is met. 137 Solution MDP Formulation For each intervalt2f1;:::;Tg, with information Y t , we define the MDP to determine actiona t using the following parameters: Decision Epoch The decision epochs – the time intervals during which the MDP makes decisions – areft;:::;Tg withT <1. States MDP uses the information regarding the current state to make a decision. For each intervalt 0 2ft;:::;Tg, we define states t 0 =<d t 0;R t 0 >, whered t 0vD t 0 denotes the predicted net load imbalances. R t 0; 0 R t 0 R denotes the available storage capacity and is defined asR t 0 =R t 0 1 sr t 0. Initial and Terminal States The initial state is denoted using a single states t1 =< 0;R t1 >. The terminal states are denoted usings T =< 0;R T > whereR T R. Note that the net load imbalance for initial and terminal states are zero. This is because the initial state imbalance is addressed before invoking the current MDP and the terminal state is outside the decision epoch. Actions MDP outputs an action from all available actions in the action space for each time interval after making a decision. The action ensures that the load imbalance for the interval is mitigated while making sure that the storage charging/discharging con- straints are satisfied. For a state s t 0 =< d t 0;R t 0 >, action space contains actions a t 0 =< b t 0;sr t 0 >;t 0 2ft;:::;Tg such that d t 0 = b t 0 +sr t 0, maxfR t 0R;cg sr t 0 minfR t 0;cg wherec denotes the per interval charge/discharge capacity of the storage system. Only the actiona 0 t fort 0 =t i.e., corresponding to the current interval is executed. The remaining actions are part of the computation framework while solving the objective function. 138 State Transition Probabilities Given states t 0 =<d t 0;R t 0 > at timet 0 , the transition probability fors t 0 +1 =<d t 0 +1 ;R t 0 +1 > is as follows: p(s t 0 +1 js t 0;a t 0) = 8 > > < > > : 0 ifR t 0 +1 6=R t 0sr t 0 P [d t 0 +1 j Y t ] otherwise whereP [d t 0 +1 j Y t ] denotes the probability distribution ofd t 0 +1 v D t 0 +1 under input data Y t . Objective The MDP solves the following objective problem: minE[ T X t 0 =t C b t 0:b t 0 +C(R T ;R 0 )] (6.19) subject to the conditions defined above. MDP Solution When input data Y t is available at the beginning of intervalt,the MDP is used to deter- mine the actiona t =<b t ;sr t > which are executed immediately. If the storage capacity at the beginning of the interval was R t1 , then the storage capacity at the end of the intervalR t =R t1 sr t . The MDP is solved as detailed below. 139 Solving MDP We define the cost to go function at timet 0 2ft;:::;Tg as follows: J t 0(R t 0;d t 0) =C b t 0:(d t 0sr t 0)+ E[C(R T ;R 0 ) + T X i=t 0 +1 C b i (D i sr i )j Y t ] (6.20) IfJ t 0 is the optimal value of the function, then it can be written recursively as: J t 0(R t 0;d t 0) = inf sr2SR fC b t 0:(d t 0sr t 0)+ X dvD t 0 +1 P [dj Y t ]J t 0 +1 (R t 0sr t 0;d)g (6.21) whereSR =fsr : maxfR t R;cgsr t minfR t ;cgg. The recurrence relation above implies that at timet 0 the best course of action is the one which minimizes the sum of the cost of operation in the current interval and the expected value of the future cost of operations. The above recurrence can be solved using standard dynamic programming techniques. The initial condition for the recurrence relation is as follows: J T +1 (R T ; 0) = 8 > > < > > : 0 ifR T R 0 1 otherwise (6.22) The optimal action at timet can be determined by solvingJ t (R t1 ; 0) after setting up the dynamic program for the MDP. The variables need to be discretized in order to solve the dynamic program. Runtime Let [a;b] denote the number of discrete entries when the range [a;b] is discretized using units i.e. [a;b] = d ba e. The total number of entries to 140 be filled are: T [0;R] max t f[d min t ;d max t ] g, where d min=max t denote the min- imum/maximum value attained by the random variable D t . Each entry requires O(jSR j max t f[d min t ;d max t ] g) time. Hence, the total time complexity is O(T [0;R] jSR j max t f[d min t ;d max t ] g 2 ). 6.4 Results and Analysis 6.4.1 2-SR Algorithms In order to perform practical evaluations, we implemented the algorithms using python. The experiments were performed on Dell optiplex with 4-cores and 4 GB RAM. Dataset We obtained the load curtailment data from the demand response implementation on our University Campus. We also used the solar dataset as generated in Section 5.6.1. The costs of the strategies were evaluated using the function f( ) = 2 2 , where is the curtailment value of the strategy. The uncertainty values were evaluated using the functionf( ) = 2 . Scalability In order to perform scalability analysis, we calculate the runtime of the algorithm by varying the number of nodes and the value of (denoted as Epsilon in the figures). We first fixed the value of = 0:2 and varied the number of nodes. As shown in Figure 6.2, the increase in runtime is polynomial (cubic) with an increase in the number of nodes. 141 Figure 6.2: Runtime of 2-SRL-Approx for varying values of nodes for fixed = 0:2 We then fix the number of nodes to 20 and vary . As shown in Figure 6.3, the runtime increases significantly with the decrease in . However, the increase is still polynomial in 1 . The runtimes of the algorithms might seem large at the first look. However, note that these are calculated using a naive python implementation. For real world scenar- ios, an implementation in faster programming languages such as C/C++ can lead to significantly better runtimes. Moreover, several trivial parallelization techniques can be adopted to get drastic improvement in the runtime. 142 Figure 6.3: Runtime of 2-SRL-Approx for varying values of for fixed number of nodes = 20 6.4.2 Risk Aware Supply Demand Matching We evaluated the framework developed in Section 6.3.5 to qualitatively assess its perfor- mance in terms of cost minimization and scalability. The framework was implemented using C++ on Dell Optiplex with 4 cores running at 2 Ghz clock frequency. The frame- work was run over a decision making horizon of 20 15-intervals (11 am - 4 pm). Dataset We used electricity consumption (load) time series data available from USC’s micro- grid [36]. We selected a subset of buildings such that the average load during peak periods (1-4 pm) was 90 kWh and during off peak periods was 70 kWh. We assumed 143 a storage capacity of 1500 kWh with a charge/discharge rate of 140 kW. The storage requirement at the end of the horizon was set as 1400 kWh, which is sufficient to supply off peak average load for 19 hours i.e. outside the decision making horizon. We used solar irradiance data available at [5] for Los Angeles USC Downtown area to calculate the solar generation data. We varied the area of solar panels and used PV output cal- culator [7] to generate three time series solar generation data which were then used to create three load imbalance time series. The three time series differed in the fraction of times either load or solar generation was surplus during the decision making horizon. Figure 6.4 shows the time series data of load and solar generation used to create the load imbalance data. The three solar time series are as follows: S1) Both load and solar surplus almost equal number of times, S2) load is surplus most of the time, and S3) solar generation is surplus most of the time. We assumed additive errors using gaussian distri- bution with mean = 0 and variance 2 = 5; 10; 20 i.e.,Pr[e t 0j Y t ]vN (0; (t 0 t) 2 ), such thatD t 0 = ^ d t 0 +e t 0, with ^ d t 0 denoting the actual value of load imbalance. Optimality We compare the cost of operation across the entire decision making horizon, i.e. 11 am - 4pm, of our framework against an optimal framework which has the correct data for the entire horizon in advance. The optimal framework solves a single MDP for decision epochsf1;:::;Tg using the available correct data. The cost of buying was set at $0.75/kWh for 11am - 1pm and $0.95/kWh for the remaining horizon. We evaluated our framework using 6 timeseries data generated i.e. three load imbalance profiles with two different variance of 10 and 20. We set the initial value of storage to 0 and 750 kWh. Figure 6.5 shows the relative percentage error of the cost of solutions obtained from our framework with respect to the optimal algorithm. As we can see from the figure, 144 11:15:00 AM 11:30:00 AM 11:45:00 AM 12:00:00 PM 12:15:00 PM 12:30:00 PM 12:45:00 PM 01:00:00 PM 01:15:00 PM 01:30:00 PM 01:45:00 PM 02:00:00 PM 02:15:00 PM 02:30:00 PM 02:45:00 PM 03:00:00 PM 03:15:00 PM 03:30:00 PM 03:45:00 PM 04:00:00 PM 0 20 40 60 80 100 120 140 160 Load S1 S2 S3 Time Interval Energy (kWh) Figure 6.4: Load and Solar Timeseries Data Used for Evaluation for each pair of initial storage-load imbalance curve, a higher variance leads to slightly higher errors. This is expected as increasing uncertainty in the input data increases the errors in decision making. However, note that the error difference is not significant implying that the framework is able to compensate for the increased uncertainty. Also, note that the errors increase with the change in the initial storage value. This is due to the fact that when the initial storage is 0 kWh, the major focus of the frame- work is to just charge the storage to its final value (1400 kWh) and fulfill the deficit. However, when the initial storage is 750 kWh, charging requirements are reduced and the framework tries to focus on participating in the external market operations by selling the excess energy. The additional actions open up new avenues for the uncertainty in the input data to introduce errors as manifested in the higher relative percentage errors. 145 0-S1 0-S2 0-S3 750-S1 750-S2 750-S3 0.00% 5.00% 10.00% 15.00% 20.00% 25.00% 30.00% Variance = 10 Variance = 20 (Initial Storage-Load Imbalance Curve) Pairs Relative Percentage Error Figure 6.5: Relative Percentage Error of Framework w.r.t. Optimal Cost For Varying Initial Storage and Variance In our future work, we will focus on reducing the errors furthers by using reinforce- ment learning based techniques which learn from the past experience to improve the accuracy for the future operations. Scalability We perform scalability analysis of the framework by varying the storage capacity and the variance of the errors. We calculate the time required to make decision at time period 0 for the entire horizon. For each storage capacity value R, we fix the charge discharge rate to 0:1R. Figure 6.6 shows the time required to make decisions with respect to the storage capacityR for various values of variance. The increase in runtime is quadratic inR. This is in agreement with the runtime which isO(RjSRj), when other variables are fixed andjSRj = O(R). The increase in runtime with variance is 146 100 200 300 400 500 600 700 800 900 1000 0 50 100 150 200 250 300 350 Var = 5.0 Var = 10.0 Var = 20.0 Storage Capacity (kWh) Time (seconds) Figure 6.6: Time Required for Decision Making w.r.t Storage Capacity for Varying Variance significant. However, even for a storage capacity of 1000 kWh and variance of 20, the runtime is just around 340 seconds. A few techniques can be followed to further improve the scalability. The entries of the dynamic table have low interdependencies. Hence, they are easy to parallelize. Moreover, as the runtime is significantly affected by the storage capacity, the smart grid can be partitioned into several smaller grids with smaller storage capacity. This improvement in runtime would be quadratic in nature. As the focus of this paper is to develop a new supply demand matching framework with reasonable performance, a detailed analysis of such techniques is omitted. 147 Chapter 7 Conclusion and Future Work In this chapter, we will provide concluding remarks to the dissertation with a discussion on the broader impacts of the work and the future directions it can lead into. 7.1 Broader Impacts The rapid increase in the penetration of Distributed Energy Resources such as PVs and batteries in the grid has decentralized energy generation [31]. This has enabled the development of microgrids wherein a localized group of electricity sources and loads constitute a subsystem which can be islanded i.e. operated without any import or export of electricity from outside the subsystem [39]. We envision future smart grids consist- ing of several microgrids operating as separate monolithic systems which interact with the larger distribution system as a single entity. A holistic microgrid controller which consists of an integrated modeling, prediction and optimization framework is impera- tive for reliable, secure and economic operation of such microgrids [9]. Optimization frameworks along the line developed in this work will be a significant component of such microgrid controllers. Moreover, in the absence of fine grained control, several components of the grid exhibit discrete operational values. As evident from Table V in the survey paper [71], while a plethora of optimization algorithms have been developed for grids with contin- uous operational values (for example, see column titled ”convex optimization problem” in the table), similar algorithms and analysis techniques are lacking for the discrete 148 case. The research community has focused only on Integer/Mixed Integer Programming based algorithms which are computationally intractable or heuristics which have no the- oretical optimality guarantees. The fast and optimal algorithms developed in this work and accompanying theoretical analyses will enable further development of sophisticated optimization frameworks for controlling discrete smart grids. 7.2 Future Directions Approximation algorithms are powerful tools to develop provably optimal polynomial runtime algorithms for a wide range of NP hard problems. In this section, we discuss a couple of open problems in the domain of smart grid and the larger domain of smart cities/infrastructure. 7.2.1 Discrete Multi Phase Supply Demand Matching The algorithms developed in this work focus on supply demand matching of active sin- gle phase power assuming no power losses due to resistance and reactance of the dis- tributors and feeders. The assumption here is that any fluctuations in reactive power that occurs by modifying the active power can be addressed by reactive power control techniques such as V olt-Var Optimization [52]. In order to create a holistic smart grid 149 control framework, this assumption needs to be addressed. Consider the relaxed branch flow model and the optimal power flow (OPF) problem defined over it [32]: min p;q;P;Q;l;v f (7.1) subject to p j = X k:j!k P jk X i:i!j (P ij r ij l ij ) +g j v j ; 8j (7.2) q j = X k:j!k Q jk X i:i!j (Q ij x ij l ij ) +b j v j ; 8j (7.3) v j =v i 2(r ij P ij +x ij Q ij ) + (r 2 ij +x 2 ij )l ij 8(i;j)2E (7.4) l ij P 2 ij +Q 2 ij v i 8(i;j)2E (7.5) The above OPF is defined on a smart grid modeled using a directed graphG = (V;E) with edges denoted using (i;j) ori!j, andi andj being the vertices.p j andq j denote the active and reactive power injected at vertexj.v j denotes the voltage magnitude.P ij , Q ij andl ij denote the active power, reactive power and the magnitude of current flowing over edgei!j respectively. Assuming all the variables as continuous entities, the above OPF is a convex relax- ation. It is an exact relaxation for radial network under certain assumptions with respect to the objective function and the maximum value ofp j andq j . Now, in order to solve the discrete version of the above OPF, the variables p j and q j will belong to a discrete set instead of being continuous variables. Sophisticated approximation algorithms to solve this formulation is an open research problem. 150 7.2.2 Joint Optimization over Multiple Infrastructure for Smart Cities The proliferation of Internet-of-Things (IoT) based technologies is able to incorporate seamlessly a large number of different and heterogeneous systems such as smart grids, smart transportation [74] thereby evolving the current cities into smart cities. Sev- eral components of such smart cities will exhibit discrete operational values and thus sophisticated approximation algorithms will be required to jointly optimize over multi- ple infrastructures. For example, the charging rates of Electric vehicles form a discrete set of values. Similarly, the problem of smart parking needs to optimize over binary variables denoting the occupancy or reservation status of parking spots [42]. Consider the following problem. A transportation company has a fleet of Electric Vehicles (trucks or cars). From a starting location, each EV needs to visit a number of pre-defined locations (for example to deliver goods) and return to the original location. Given a graphical model of the transportation network, and assuming the costs are in metric space, this problem is the classical traveling salesman problem (TSP) with an approximation factor of 1.5. In the context of smart cities, the problem can be extended as follows. In addition to visiting the locations for delivery, the EVs can also visit charging station. In the charging station, the EVs can charge themselves at some cost or provide the grid with a portion of the energy stored in their batteries, thereby receiving a monetary compensation. Hence, the problem that needs to solved is to minimize the sum of the charging costs, the travel- ing costs minus the discharging compensation while ensuring all the locations are visited within a fixed deadline. 151 7.3 Concluding Remarks Approximation algorithms are strong techniques to solve several intractable real world problems in a reasonable amount of time while providing theoretical worst case guar- antees. However, their application in the domain of smart grids and the more general domain of smart cities has been limited. In this work, we show that using these tech- niques for solving the problem of discrete supply demand matching in smart grids can lead to significant improvements over the state of the art. The work done in this dissertation are first steps towards wider applications of approximation algorithms in these domains. We hope that this dissertation can spark the interest of research community in the applications of approximation algorithms for smart grids and smart cities. This will lead to the development of scalable and reliable solutions for various challenging problems in these domains. Moreover, new theoretical insights developed as a result will broaden the understanding of several core problems of interest to the theoretical computer science community. 152 Reference List [1] Ilog cplex optimization studio welcome page. Online: http://www- 01.ibm.com/support/knowledgecenter/SSSA5P/welcome. [2] Opl language reference manual. [3] Usc smart grid. Online. [4] Car battery efficiencies. Online: http://large.stanford.edu/courses/2010/ph240/sun1/, 2010. [5] National solar radiation data base, 2010. [6] Photovoltaic (pv) pricing trends: Historical, recent, and near-term projections, 2013. [7] How to calculate the annual solar energy output of a photovoltaic system?, 2013- 2016. [8] Matlab-mathworks. https://www.mathworks.com/products/matlab.html, 2017. [9] Deep solar. Online: http://deepsolar.usc.edu/, 2018. [10] S. Alamdari, T. Biedl, T. M. Chan, E. Grant, K. R. Jampani, S. Keshav, A. Lubiw, and V . Pathak. Smart-grid electricity allocation via strip packing with slicing. In Workshop on Algorithms and Data Structures, pages 25–36. Springer, 2013. [11] M. H. Albadi and E. El-Saadany. Demand response in electricity markets: An overview. In IEEE power engineering society general meeting, volume 2007, pages 1–5, 2007. [12] S. Aman, C. Chelmis, and V . Prasanna. Learning to reduce: A reduced electricity consumption prediction ensemble. In Workshops at the Thirtieth AAAI Conference on Artificial Intelligence, 2016. 153 [13] S. Aman, C. Chelmis, and V . K. Prasanna. Influence-driven model for time series prediction from partial observations. In Twenty-Ninth AAAI Conference on Artifi- cial Intelligence, 2015. [14] S. Aman, M. Frincu, C. Chelmis, M. Noor, Y . Simmhan, and V . K. Prasanna. Prediction models for dynamic demand response. In IEEE International Confer- ence on Smart Grid Communications (SmartGridComm): Data Management, Grid Analytics, and Dynamic Pricing, 2015. [15] N. Bansal, N. Korula, V . Nagarajan, and A. Srinivasan. Solving packing inte- ger programs via randomized rounding with alterations. Theory of Computing, 8(1):533–565, 2012. [16] A. Barbato, C. Bolchini, M. Delfanti, A. Geronazzo, G. Accetta, A. Dede, G. Massa, and M. Trioni. An energy management framework for optimal demand response in a smart campus. ICGREEN, page 11, 2015. [17] A. Barbato, A. Capone, L. Chen, F. Martignon, and S. Paris. A power scheduling game for reducing the peak demand of residential users. In Online Conference on Green Communications (GreenCom), 2013 IEEE, pages 137–142. IEEE, 2013. [18] M. Behrangrad, H. Sugihara, and T. Funaki. Analyzing the system effects of opti- mal demand response utilization for reserve procurement and peak clipping. In Power and Energy Society General Meeting, 2010 IEEE, pages 1–7. IEEE, 2010. [19] R. Bellman. A markovian decision process. Indiana Univ. Math. J., 6:679–684, 1957. [20] S. Boyd and L. Vandenberghe. Convex optimization. Cambridge university press, 2004. [21] M. Charikar, C. Chekuri, and M. P´ al. Sampling bounds for stochastic optimization. In Approximation, Randomization and Combinatorial Optimization. Algorithms and Techniques, pages 257–269. Springer, 2005. [22] C. Chelmis, S. Aman, M. R. Saeed, M. Frincu, and V . K. Prasanna. Estimating reduced consumption for dynamic demand response. In Proceedings of the Twenty- Ninght AAAI Conference on Artificial Intelligence. AAAI Press, 2015. [23] J. Chen, B. Yang, and X. Guan. Optimal demand response scheduling with stack- elberg game approach under load uncertainty for smart grid. In Smart Grid Com- munications (SmartGridComm), 2012 IEEE Third International Conference on, pages 546–551. IEEE, 2012. 154 [24] Z. Chen, L. Wu, and Y . Fu. Real-time price-based demand response manage- ment for residential appliances via stochastic optimization and robust optimization. IEEE Transactions on Smart Grid, 3(4):1822–1831, 2012. [25] S. Choi, S. Park, D.-J. Kang, S.-j. Han, and H.-M. Kim. A microgrid energy management system for inducing optimal demand response. In Smart Grid Com- munications (SmartGridComm), 2011 IEEE International Conference on, pages 19–24. IEEE, 2011. [26] A. I. Cohen and C. C. Wang. An optimization method for load management scheduling. IEEE Trans. Power Syst.;(United States), 3(2), 1988. [27] C. Deline. Partially shaded operation of multi-string photovoltaic systems. In Pho- tovoltaic Specialists Conference (PVSC), 2010 35th IEEE, pages 000394–000399. IEEE, 2010. [28] Y . M. Ding, S. H. Hong, and X. H. Li. A demand response energy management scheme for industrial facilities in smart grid. IEEE Transactions on Industrial Informatics, 10(4):2257–2269, 2014. [29] J. Driesen and F. Katiraei. Design for distributed energy resources. IEEE Power and Energy Magazine, 6(3), 2008. [30] A. T. Elsayed, C. R. Lashway, and O. A. Mohammed. Advanced battery manage- ment and diagnostic system for smart grid infrastructure. IEEE Transactions on Smart Grid, 7(2):897–905, 2016. [31] H. Farhangi. The path of the smart grid. IEEE power and energy magazine, 8(1), 2010. [32] M. Farivar and S. H. Low. Branch flow model: Relaxations and convexificationpart i. IEEE Transactions on Power Systems, 28(3):2554–2564, 2013. [33] O. Gagrica, P. H. Nguyen, W. L. Kling, and T. Uhl. Microinverter curtailment strategy for increasing photovoltaic penetration in low-voltage networks. IEEE Transactions on Sustainable Energy, 6(2):369–379, 2015. [34] N. Gatsis and G. B. Giannakis. Cooperative multi-residence demand response scheduling. In Information Sciences and Systems (CISS), 2011 45th Annual Con- ference on, pages 1–6. IEEE, 2011. [35] O. H. Ibarra and C. E. Kim. Fast approximation algorithms for the knapsack and sum of subset problems. Journal of the ACM (JACM), 22(4):463–468, 1975. 155 [36] S. R. Kuppannagari, R. Kannan, C. Chelmis, and V . K. Prasanna. Implementation of learning-based dynamic demand response on a campus micro-grid. In The 25th International Joint Conference on Artificial Intelligence, 2016. [37] S. R. Kuppannagari, R. Kannan, C. Chelmis, A. S. Tehrani, and V . K. Prasanna. Optimal customer targeting for sustainable demand response in smart grids. Pro- cedia Computer Science, 80:324–334, 2016. [38] J. Kwac and R. Rajagopal. Demand response targeting using big data analytics. In Big Data, 2013 IEEE International Conference on, pages 683–690. IEEE, 2013. [39] R. H. Lasseter and P. Paigi. Microgrid: A conceptual solution. In Power Electron- ics Specialists Conference, 2004. PESC 04. 2004 IEEE 35th Annual, volume 6, pages 4285–4290. IEEE, 2004. [40] S. Lee, S. Iyengar, D. Irwin, and P. Shenoy. Distributed rate control for smart solar arrays. In Proceedings of the Eighth International Conference on Future Energy Systems, pages 34–44. ACM, 2017. [41] T. Logenthiran, D. Srinivasan, and T. Z. Shun. Demand side management in smart grid using heuristic optimization. Smart Grid, IEEE Transactions on, 3(3):1244– 1252, 2012. [42] R. Lu, X. Lin, H. Zhu, and X. Shen. Spark: a new vanet-based smart parking scheme for large parking lots. In INFOCOM 2009, IEEE, pages 1413–1421. IEEE, 2009. [43] L. Lutzenhiser, L. Cesafsky, H. Chappells, M. Gossard, M. Moezzi, D. Moran, J. Peters, M. Spahic, P. Stern, E. Simmons, et al. Behavioral assumptions under- lying california residential sector energy efficiency programs. Prepared for the California Institute for Energy and Environment Behavior and Energy Program, 2009. [44] Y . V . Makarov, C. Loutan, J. Ma, and P. De Mello. Operational impacts of wind generation on california power systems. IEEE Transactions on Power Systems, 24(2):1039–1050, 2009. [45] F. Mangiatordi, E. Pallotti, P. D. Vecchio, and F. Leccese. Power consumption scheduling for residential buildings. In Environment and Electrical Engineering (EEEIC), 2012 11th International Conference on, pages 926–930. IEEE, 2012. [46] A. Molderink, V . Bakker, M. G. Bosman, J. L. Hurink, and G. J. Smit. Domes- tic energy management methodology for optimizing efficiency in smart grids. In PowerTech, 2009 IEEE Bucharest, pages 1–7. IEEE, 2009. 156 [47] H. Moradi, D. Groff, and A. Abtahi. Optimal energy schadualing of a stand-alone multi-sourced microgid considering environmental aspects. In 2017 IEEE Innova- tive Smart Grid Technologies Conference (ISGT), 2017. [48] K. Moslehi and R. Kumar. A reliability perspective of the smart grid. Smart Grid, IEEE Transactions on, 1(1):57–64, 2010. [49] S. J. Moss, M. Cubed, and K. Fleisher. Market segmentation and energy efficiency program design. Berkeley, California Institute for Energy and Environment, 2008. [50] N. Motegi, M. A. Piette, D. S. Watson, S. Kiliccote, and P. Xu. Introduction to commercial building control strategies and techniques for demand response. Lawrence Berkeley National Laboratory LBNL-59975, 2007. [51] F. Rahimi and A. Ipakchi. Demand response as a market resource under the smart grid paradigm. Smart Grid, IEEE Transactions on, 1(1):82–88, 2010. [52] S. Rahimi, M. Marinelli, and F. Silvestro. Evaluation of requirements for volt/var control and optimization function in distribution management systems. In Energy Conference and Exhibition (ENERGYCON), 2012 IEEE International, pages 331– 336. IEEE, 2012. [53] P. Ramsami and V . Oree. A hybrid method for forecasting the energy output of photovoltaic systems. Energy Conversion and Management, 95:406–413, 2015. [54] A. Ranjan, P. Khargonekar, and S. Sahni. Offline preemptive scheduling of power demands to minimize peak power in smart grids. In Computers and Communica- tion (ISCC), 2014 IEEE Symposium on, pages 1–6. IEEE, 2014. [55] W. Romeijnders, L. Stougie, and M. H. van der Vlerk. Approximation in two-stage stochastic integer programming. Surveys in Operations Research and Management Science, 19(1):17–33, 2014. [56] S. Rongali, T. Ganuy, M. Padmanabhan, V . Arya, S. Kalyanaraman, and M. I. Petra. iplug: Decentralised dispatch of distributed generation. In Communica- tion Systems and Networks (COMSNETS), 2016 8th International Conference on, pages 1–8. IEEE, 2016. [57] T. Roy, A. Das, and Z. Ni. Optimization in load scheduling of a residential com- munity using dynamic pricing. In 2017 IEEE Innovative Smart Grid Technologies Conference (ISGT), 2017. [58] N. Ruiz, I. Cobelo, and J. Oyarzabal. A direct load control model for virtual power plant management. IEEE Transactions on Power Systems, 24(2):959–966, 2009. 157 [59] A. Ruszczy´ nski and A. Shapiro. Stochastic programming models. Handbooks in operations research and management science, 10:1–64, 2003. [60] S. Sahni. Approximate algorithms for the 0/1 knapsack problem. Journal of the ACM (JACM), 22(1):115–124, 1975. [61] A. Sepulveda, L. Paull, W. G. Morsi, H. Li, C. Diduch, and L. Chang. A novel demand side management program using water heaters and particle swarm opti- mization. In Electric Power and Energy Conference (EPEC), 2010 IEEE, pages 1–5. IEEE, 2010. [62] Y . Simmhan, V . Prasanna, S. Aman, S. Natarajan, W. Yin, and Q. Zhou. Toward data-driven demand-response optimization in a campus microgrid. In Proceedings of the Third ACM Workshop on Embedded Sensing Systems for Energy-Efficiency in Buildings, pages 41–42. ACM, 2011. [63] A. Singh, S. Lee, D. Irwin, and P. Shenoy. Sunshade: enabling software-defined solar-powered systems. In Proceedings of the 8th International Conference on Cyber-Physical Systems, pages 61–70. ACM, 2017. [64] B. A. Smith, J. Wong, and R. Rajagopal. A simple way to use interval data to segment residential customers for energy efficiency and demand response program targeting. In ACEEE Summer Study on Energy Efficiency in Buildings, 2012. [65] P. Srikantha and D. Kundur. Distributed demand curtailment via water-filling. In IEEE International Conference on Smart Grid Communications, 2015. [66] H.-I. Su and A. El Gamal. Modeling and analysis of the role of fast-response energy storage in the smart grid. In Communication, Control, and Computing (Allerton), 2011 49th Annual Allerton Conference on, pages 719–726. IEEE, 2011. [67] SunShot. Sunshot vision study a comprehensive analysis of the potential for u.s. solar electricity generation, 2012. [68] S. Tang, Q. Huang, X.-Y . Li, and D. Wu. Smoothing the energy consumption: Peak demand reduction in smart grid. In INFOCOM, 2013 Proceedings IEEE, pages 1133–1141. IEEE, 2013. [69] C.-W. Ten, C.-C. Liu, and G. Manimaran. Vulnerability assessment of cybersecu- rity for scada systems. IEEE Transactions on Power Systems, 23(4):1836–1846, 2008. [70] R. Tonkoski, L. A. Lopes, and T. H. El-Fouly. Coordinated active power curtail- ment of grid connected pv inverters for overvoltage prevention. IEEE Transactions on Sustainable Energy, 2(2):139–147, 2011. 158 [71] J. S. Vardakas, N. Zorba, and C. V . Verikoukis. A survey on demand response programs in smart grids: Pricing methods and optimization algorithms. IEEE Communications Surveys & Tutorials, 17(1):152–178, 2015. [72] V . V . Vazirani. Approximation algorithms. Springer Science & Business Media, 2013. [73] S. Yaw, B. Mumey, E. McDonald, and J. Lemke. Peak demand scheduling in the smart grid. In Smart Grid Communications (SmartGridComm), 2014 IEEE International Conference on, pages 770–775. IEEE, 2014. [74] A. Zanella, N. Bui, A. Castellani, L. Vangelista, and M. Zorzi. Internet of things for smart cities. IEEE Internet of Things journal, 1(1):22–32, 2014. [75] Z. Zhu, J. Tang, S. Lambotharan, W. H. Chin, and Z. Fan. An integer linear pro- gramming based optimization for home demand-side management in smart grid. In Innovative Smart Grid Technologies (ISGT), 2012 IEEE PES, pages 1–5. IEEE, 2012. [76] V . Zois, M. Frincu, C. Chelmis, M. R. Saeed, and V . Prasanna. Efficient customer selection for sustainable demand response in smart grids. In IGCC, pages 1–6. IEEE, 2014. [77] V . Zois, M. Frincu, and V . Prasanna. Integrated platform for automated sustainable demand response in smart grids. In Intelligent Energy Systems (IWIES), 2014 IEEE International Workshop on, pages 64–69. IEEE, 2014. 159

Abstract (if available)

Abstract Traditional power grids consisted of large grid operator owned generators to supply the demand of the consumers. The inertia of the large generators was able to handle any sudden demand fluctuations followed by an eventual ramp up/down of generation. However, modern smart grids, characterized by the proliferation of small sized Distributed Energy Resources (DER) such as photovoltaics, energy storage etc. lack this ability of handling demand fluctuations due to their lack of inertia. Therefore, they require proactive supply demand matching by utilizing the available advanced sensing and remote control capability. ❧ A plethora of work has been done assuming that the nodes (supply and demand) in the grid exhibit continuous operational values (generation/consumption values). Fine grained control is required for continuous operational values which might not be economically available. Thus, supply demand matching often needs to be performed over discrete sets of operational values of the nodes making this an NP-hard problem. While heuristics (with unbounded error) and Integer Program based algorithms (computationally expensive) have been studied extensively, the application of approximation algorithms to develop fast algorithms with bounded error guarantees has been limited. ❧ In this work, we address these limitations by developing a suite of polynomial runtime complexity approximation algorithms to perform supply demand matching in smart grids. We first model the smart grid to capture the available discrete control strategies for each node in the smart grid and the associated costs and uncertainties. We also model practical grid constraints such as the overhead associated with switching strategies and the smart grid network capacity constraints. We then develop approximation algorithms with a variety of objectives and constraints such as minimizing curtailment error in Demand Response, minimizing cost of supply demand matching, fairness etc. and theoretically provide the optimality guarantees. In addition to the theoretical analysis, we perform practical evaluations of the algorithms and show that our supply demand matching algorithms provide solutions which are close to optimal in a small amount of time.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

A function-based methodology for evaluating resilience in smart grids

PDF

Dynamic graph analytics for cyber systems security applications

PDF

Data-driven methods for increasing real-time observability in smart distribution grids

PDF

Prediction models for dynamic decision making in smart grid

PDF

Demand response management in smart grid from distributed optimization perspective

PDF

Exploiting variable task granularities for scalable and efficient parallel graph analytics

PDF

The smart grid network: pricing, markets and incentives

PDF

Scaling up temporal graph learning: powerful models, efficient algorithms, and optimized systems

PDF

Distribution system reliability analysis for smart grid applications

PDF

A complex event processing framework for fast data management

PDF

Defending industrial control systems: an end-to-end approach for managing cyber-physical risk

PDF

Landscape analysis and algorithms for large scale non-convex optimization

PDF

I. Asynchronous optimization over weakly coupled renewal systems

PDF

Adaptive and resilient stream processing on cloud infrastructure

PDF

Optimizing task assignment for collaborative computing over heterogeneous network devices

PDF

Model-driven situational awareness in large-scale, complex systems

PDF

Scalable exact inference in probabilistic graphical models on multi-core platforms

PDF

User-centric smart sensing for non-intrusive electricity consumption disaggregation in buildings

PDF

Scaling up deep graph learning: efficient algorithms, expressive models and fast acceleration

PDF

From matching to querying: A unified framework for ontology integration

Asset Metadata

Creator Kuppannagari, Sanmukh Rao (author)

Core Title Discrete optimization for supply demand matching in smart grids

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Computer Engineering

Publication Date 07/26/2018

Defense Date 06/11/2018

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag approximation algorithms,discrete supply demand matching,OAI-PMH Harvest,smart grids

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Prasanna, Viktor K. (committee chair), Beshir, Mohammed (committee member), Dughmi, Shaddin (committee member), Jovanovic, Mihailo R. (committee member), Kannan, Rajgopal (committee member)

Creator Email i.sanmukh@gmail.com,kuppanna@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-34412

Unique identifier UC11671739

Identifier etd-Kuppannaga-6517.pdf (filename),usctheses-c89-34412 (legacy record id)

Legacy Identifier etd-Kuppannaga-6517.pdf

Dmrecord 34412

Document Type Dissertation

Format application/pdf (imt)

Rights Kuppannagari, Sanmukh Rao

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA