# MATPOWER Manual

### MATPOWER-manual

### MATPOWER-manual

User Manual:

Open the PDF directly: View PDF .

Page Count: 246 [warning: Documents this large are best viewed by clicking the View PDF Link!]

- Introduction
- Getting Started
- Modeling
- Power Flow
- Continuation Power Flow
- Optimal Power Flow
- Extending the OPF
- Unit De-commitment Algorithm
- Miscellaneous Matpower Functions
- Acknowledgments
- Appendix MIPS – Matpower Interior Point Solver
- Appendix Data File Format
- Appendix Matpower Options
- Appendix Matpower Files and Functions
- Appendix Extras Directory
- Appendix ``Smart Market'' Code
- Appendix Optional Packages
- BPMPD_MEX – MEX interface for BPMPD
- CLP – COIN-OR Linear Programming
- CPLEX – High-performance LP and QP Solvers
- GLPK – GNU Linear Programming Kit
- Gurobi – High-performance LP and QP Solvers
- Ipopt – Interior Point Optimizer
- KNITRO – Non-Linear Programming Solver
- MINOPF – AC OPF Solver Based on MINOS
- MOSEK – High-performance LP and QP Solvers
- Optimization Toolbox – LP, QP, NLP and MILP Solvers
- PARDISO – Parallel Sparse Direct and Multi-Recursive Iterative Linear Solvers
- SDP_PF – Applications of a Semidefinite Programming Relaxation of the Power Flow Equations
- TSPOPF – Three AC OPF Solvers by H. Wang

- Appendix Release History
- Pre 1.0 – released Jun 25, 1997
- Version 1.0 – released Sep 17, 1997
- Version 1.0.1 – released Sep 19, 1997
- Version 2.0 – released Dec 24, 1997
- Version 3.0 – released Feb 14, 2005
- Version 3.2 – released Sep 21, 2007
- Version 4.0 – released Feb 7, 2011
- Version 4.1 – released Dec 14, 2011
- Version 5.0 – released Dec 17, 2014
- Version 5.1 – released Mar 20, 2015
- Version 6.0 – released Dec 16, 2016
- Version 7.0 – beta 1 released released Oct 31, 2018

- References

MATP WER

User’s Manual

Version 7.0b1

Ray D. Zimmerman Carlos E. Murillo-S´anchez

October 31, 2018

©2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 Power Systems Engineering Research Center (PSerc)

All Rights Reserved

Contents

1 Introduction 10

1.1 Background ................................ 10

1.2 License and Terms of Use ........................ 11

1.3 Citing Matpower ............................ 12

1.4 Matpower Development ........................ 12

2 Getting Started 13

2.1 System Requirements ........................... 13

2.2 Getting Matpower ........................... 14

2.2.1 Versioned Releases ........................ 14

2.2.2 Current Development Version .................. 14

2.3 Installation ................................ 16

2.4 Running a Simulation ........................... 18

2.4.1 Preparing Case Input Data .................... 18

2.4.2 Solving the Case ......................... 18

2.4.3 Accessing the Results ....................... 19

2.4.4 Setting Options .......................... 20

2.5 Documentation .............................. 20

3 Modeling 23

3.1 Data Formats ............................... 23

3.2 Branches .................................. 23

3.3 Generators ................................. 25

3.4 Loads ................................... 25

3.5 Shunt Elements .............................. 26

3.6 Network Equations ............................ 26

3.7 DC Modeling ............................... 27

4 Power Flow 30

4.1 AC Power Flow .............................. 30

4.2 DC Power Flow .............................. 32

4.3 Distribution Power Flow ......................... 32

4.3.1 Radial Power Flow ........................ 33

4.3.2 Current Summation Method ................... 34

4.3.3 Power Summation Method .................... 35

4.3.4 Admittance Summation Method ................. 35

4.3.5 Handling PV Buses ........................ 37

4.4 runpf ................................... 38

4.5 Linear Shift Factors ............................ 39

5 Continuation Power Flow 43

5.1 Parameterization ............................. 44

5.2 Predictor .................................. 44

5.3 Corrector ................................. 45

5.4 Step Length Control ........................... 45

5.5 Event Detection and Location ...................... 46

5.6 runcpf ................................... 47

5.6.1 CPF Callback Functions ..................... 50

5.6.2 CPF Example ........................... 53

6 Optimal Power Flow 56

6.1 Standard AC OPF ............................ 56

6.1.1 Cartesian vs. Polar Coordinates for Voltage .......... 57

6.1.2 Current vs. Power for Nodal Balance Constraints ....... 58

6.2 Standard DC OPF ............................ 59

6.3 Extended OPF Formulation ....................... 60

6.3.1 User-deﬁned Variables ...................... 61

6.3.2 User-deﬁned Constraints ..................... 61

6.3.3 User-deﬁned Costs ........................ 62

6.4 Standard Extensions ........................... 64

6.4.1 Piecewise Linear Costs ...................... 64

6.4.2 Dispatchable Loads ........................ 66

6.4.3 Generator Capability Curves ................... 68

6.4.4 Branch Angle Diﬀerence Limits ................. 69

6.5 Solvers ................................... 69

6.6 runopf ................................... 70

7 Extending the OPF 76

7.1 Direct Speciﬁcation ............................ 76

7.1.1 User-deﬁned Variables ...................... 76

7.1.2 User-deﬁned Constraints ..................... 76

7.1.3 User-deﬁned Costs ........................ 78

7.1.4 Additional Comments ...................... 79

7.2 Callback Functions ............................ 79

7.2.1 User-deﬁned Variables ...................... 80

7.2.2 User-deﬁned Costs ........................ 80

2

7.2.3 User-deﬁned Constraints ..................... 81

7.3 Callback Stages and Example ...................... 82

7.3.1 ext2int Callback ......................... 83

7.3.2 formulation Callback ...................... 84

7.3.3 int2ext Callback ......................... 88

7.3.4 printpf Callback ......................... 91

7.3.5 savecase Callback ........................ 93

7.4 Registering the Callbacks ......................... 95

7.5 Summary ................................. 97

7.6 Example Extensions ........................... 97

7.6.1 Fixed Zonal Reserves ....................... 97

7.6.2 Interface Flow Limits ....................... 99

7.6.3 DC Transmission Lines ...................... 100

7.6.4 OPF Soft Limits ......................... 103

8 Unit De-commitment Algorithm 110

9 Miscellaneous Matpower Functions 112

9.1 Input/Output Functions ......................... 112

9.1.1 loadcase ............................. 112

9.1.2 savecase ............................. 112

9.1.3 cdf2mpc .............................. 113

9.1.4 psse2mpc ............................. 113

9.1.5 save2psse ............................. 114

9.2 System Information ............................ 114

9.2.1 case info ............................. 114

9.2.2 compare case ........................... 114

9.2.3 find islands ........................... 115

9.2.4 get losses ............................ 115

9.2.5 margcost ............................. 116

9.2.6 isload ............................... 116

9.2.7 loadshed ............................. 116

9.2.8 printpf .............................. 116

9.2.9 total load ............................ 117

9.2.10 totcost .............................. 117

9.3 Modifying a Case ............................. 117

9.3.1 extract islands ......................... 117

9.3.2 load2disp ............................. 118

3

9.3.3 modcost .............................. 118

9.3.4 scale load ............................ 118

9.3.5 apply changes .......................... 119

9.3.6 savechgtab ............................ 122

9.4 Conversion between External and Internal Numbering ......... 123

9.4.1 ext2int,int2ext ......................... 123

9.4.2 e2i data,i2e data ........................ 123

9.4.3 e2i field,i2e field ...................... 124

9.5 Forming Standard Power Systems Matrices ............... 125

9.5.1 makeB ............................... 125

9.5.2 makeBdc .............................. 125

9.5.3 makeJac .............................. 125

9.5.4 makeLODF ............................. 126

9.5.5 makePTDF ............................. 126

9.5.6 makeYbus ............................. 126

9.6 Miscellaneous ............................... 127

9.6.1 define constants ........................ 127

9.6.2 feval w path ........................... 127

9.6.3 have fcn ............................. 127

9.6.4 mpopt2qpopt ........................... 128

9.6.5 mpver ............................... 129

9.6.6 nested struct copy ....................... 129

10 Acknowledgments 130

Appendix A MIPS – Matpower Interior Point Solver 131

A.1 Example 1 ................................. 133

A.2 Example 2 ................................. 135

A.3 Quadratic Programming Solver ..................... 137

A.4 Primal-Dual Interior Point Algorithm .................. 138

A.4.1 Notation .............................. 138

A.4.2 Problem Formulation and Lagrangian .............. 139

A.4.3 First Order Optimality Conditions ............... 140

A.4.4 Newton Step ........................... 141

Appendix B Data File Format 144

Appendix C Matpower Options 150

C.1 Mapping of Old-Style Options to New-Style Options .......... 166

4

Appendix D Matpower Files and Functions 170

D.1 Directory Layout and Documentation Files ............... 170

D.2 Matpower Functions .......................... 172

D.3 Example Matpower Cases ....................... 184

D.4 Automated Test Suite .......................... 187

Appendix E Extras Directory 191

Appendix F “Smart Market” Code 193

F.1 Handling Supply Shortfall ........................ 195

F.2 Example .................................. 195

F.3 Smartmarket Files and Functions .................... 200

Appendix G Optional Packages 201

G.1 BPMPD MEX – MEX interface for BPMPD .............. 201

G.2 CLP – COIN-OR Linear Programming ................. 201

G.3 CPLEX – High-performance LP and QP Solvers ............ 202

G.4 GLPK – GNU Linear Programming Kit ................ 203

G.5 Gurobi – High-performance LP and QP Solvers ............ 203

G.6 Ipopt – Interior Point Optimizer .................... 204

G.7 KNITRO – Non-Linear Programming Solver .............. 205

G.8 MINOPF – AC OPF Solver Based on MINOS ............. 206

G.9 MOSEK – High-performance LP and QP Solvers ........... 206

G.10 Optimization Toolbox – LP, QP, NLP and MILP Solvers . . . . . . . 207

G.11 PARDISO – Parallel Sparse Direct and Multi-Recursive Iterative Lin-

ear Solvers ................................. 207

G.12 SDP PF – Applications of a Semideﬁnite Programming Relaxation of

the Power Flow Equations ........................ 208

G.13 TSPOPF – Three AC OPF Solvers by H. Wang ............ 208

Appendix H Release History 210

H.1 Pre 1.0 – released Jun 25, 1997 ..................... 210

H.2 Version 1.0 – released Sep 17, 1997 ................... 210

H.3 Version 1.0.1 – released Sep 19, 1997 .................. 210

H.4 Version 2.0 – released Dec 24, 1997 ................... 211

H.5 Version 3.0 – released Feb 14, 2005 ................... 212

H.6 Version 3.2 – released Sep 21, 2007 ................... 213

H.7 Version 4.0 – released Feb 7, 2011 .................... 215

H.8 Version 4.1 – released Dec 14, 2011 ................... 218

5

List of Figures

3-1 Branch Model ............................... 24

4-1 Oriented Ordering ............................ 33

4-2 Branch Representation: branch kbetween buses i(sending) and k

(receiving) and load demand and shunt admittances at both buses . . 34

5-1 Nose Curve of Voltage Magnitude at Bus 9 ............... 53

6-1 Relationship of wito rifor di= 1 (linear option) ............ 63

6-2 Relationship of wito rifor di= 2 (quadratic option) ......... 64

6-3 Constrained Cost Variable ........................ 65

6-4 Marginal Beneﬁt or Bid Function .................... 66

6-5 Total Cost Function for Negative Injection ............... 67

6-6 Generator P-QCapability Curve .................... 68

7-1 Adding Constraints Across Subsets of Variables ............ 87

7-2 DC Line Model .............................. 101

7-3 Equivalent “Dummy” Generators .................... 101

7-4 Feasible Region for Branch Flow Violation Constraints ........ 104

List of Tables

4-1 Power Flow Results ............................ 39

4-2 Power Flow Options ........................... 40

4-3 Power Flow Output Options ....................... 41

5-1 Continuation Power Flow Results .................... 48

5-2 Continuation Power Flow Options .................... 49

5-3 Continuation Power Flow Callback Input Arguments ......... 51

5-4 Continuation Power Flow Callback Output Arguments ........ 52

5-5 Continuation Power Flow State ..................... 52

6-1 Optimal Power Flow Results ....................... 71

6-2 Optimal Power Flow Solver Options ................... 73

6-3 Other OPF Options ............................ 74

6-4 OPF Output Options ........................... 75

7-1 User-deﬁned Nonlinear Constraint Speciﬁcation ............ 78

7-2 Names Used by Implementation of OPF with Reserves ........ 83

7-3 Results for User-Deﬁned Variables, Constraints and Costs ....... 89

7-4 Callback Functions ............................ 97

7-5 Input Data Structures for Fixed Zonal Reserves ............ 98

7-6 Output Data Structures for Fixed Zonal Reserves ........... 98

7

7-7 Input Data Structures for Interface Flow Limits ............ 99

7-8 Output Data Structures for Interface Flow Limits ........... 100

7-9 Soft Limit Formulation .......................... 105

7-10 Input Data Structures for OPF Soft Limits ............... 106

7-11 Default Soft Limit Values ........................ 107

7-12 Possible Hard-Limit Modiﬁcations .................... 107

7-13 Output Data Structures for OPF Soft Limits .............. 108

9-1 Columns of chgtab ............................ 120

9-2 Values for CT TABLE Column ....................... 120

9-3 Values for CT CHGTYPE Column ...................... 121

9-4 Values for CT COL Column ........................ 121

A-1 Input Arguments for mips ........................ 132

A-2 Output Arguments for mips ....................... 133

A-3 Options for mips ............................. 134

B-1 Bus Data (mpc.bus)............................ 145

B-2 Generator Data (mpc.gen)........................ 146

B-3 Branch Data (mpc.branch)........................ 147

B-4 Generator Cost Data (mpc.gencost)................... 148

B-5 DC Line Data (mpc.dcline)....................... 149

C-1 Top-Level Options ............................ 152

C-2 Power Flow Options ........................... 153

C-3 Continuation Power Flow Options .................... 154

C-4 OPF Solver Options ........................... 155

C-5 General OPF Options .......................... 156

C-6 Power Flow and OPF Output Options ................. 157

C-7 OPF Options for MIPS .......................... 158

C-8 OPF Options for CLP .......................... 158

C-9 OPF Options for CPLEX ........................ 159

C-10 OPF Options for fmincon ........................ 160

C-11 OPF Options for GLPK ......................... 160

C-12 OPF Options for Gurobi ......................... 161

C-13 OPF Options for Ipopt ......................... 161

C-14 OPF Options for KNITRO ........................ 162

C-15 OPF Options for MINOPF ........................ 163

C-16 OPF Options for MOSEK ........................ 164

C-17 OPF Options for PDIPM ......................... 165

C-18 OPF Options for TRALM ........................ 165

C-19 Old-Style to New-Style Option Mapping ................ 166

8

D-1 Matpower Directory Layout and Documentation Files . . . . . . . 171

D-2 Top-Level Simulation Functions ..................... 172

D-3 Input/Output Functions ......................... 172

D-4 Data Conversion Functions ........................ 173

D-5 Power Flow Functions .......................... 173

D-6 Continuation Power Flow Functions ................... 174

D-7 OPF and Wrapper Functions ...................... 175

D-8 OPF Model Objects ........................... 176

D-9 Deprecated @opt model Methods ..................... 177

D-10 OPF Solver Functions .......................... 177

D-11 Other OPF Functions ........................... 178

D-12 OPF User Callback Functions ...................... 179

D-13 Power Flow Derivative Functions .................... 179

D-14 NLP, LP & QP Solver Functions .................... 180

D-15 Matrix Building Functions ........................ 181

D-16 Utility Functions ............................. 182

D-17 Other Functions .............................. 183

D-18 Small Transmission System Test Cases ................. 184

D-19 Small Radial Distribution System Test Cases .............. 184

D-20 ACTIV Synthetic Grid Test Cases .................... 185

D-21 Polish System Test Cases ......................... 185

D-22 PEGASE European System Test Cases ................. 186

D-23 RTE French System Test Cases ..................... 186

D-24 Automated Test Functions from MP-Test ................ 187

D-25 MIPS Tests ................................ 187

D-26 Test Data ................................. 188

D-27 Miscellaneous Matpower Tests ..................... 189

D-28 Matpower OPF Tests ......................... 190

F-1 Auction Types .............................. 194

F-2 Generator Oﬀers ............................. 196

F-3 Load Bids ................................. 196

F-4 Generator Sales .............................. 199

F-5 Load Purchases .............................. 199

F-6 Smartmarket Files and Functions .................... 200

9

MATP WER

1 Introduction

1.1 Background

Matpower [1] is a package of Matlab®M-ﬁles for solving power ﬂow and optimal

power ﬂow problems. It is intended as a simulation tool for researchers and educators

that is easy to use and modify. Matpower is designed to give the best performance

possible while keeping the code simple to understand and modify. The Matpower

website can be found at:

http://www.pserc.cornell.edu/matpower/

Matpower was initially developed by Ray D. Zimmerman, Carlos E. Murillo-

S´anchez and Deqiang Gan of PSerc1at Cornell University under the direction of

Robert J. Thomas. The initial need for Matlab-based power ﬂow and optimal power

ﬂow code was born out of the computational requirements of the PowerWeb project2.

Many others have contributed to Matpower over the years and it continues to be

developed and maintained under the direction of Ray Zimmerman.

Beginning with version 6, Matpower includes a framework for solving general-

ized steady-state electric power scheduling problems. This framework is known as

MOST, for Matpower Optimal Scheduling Tool [2,3].

MOST can be used to solve problems as simple as a deterministic, single pe-

riod economic dispatch problem with no transmission constraints or as complex as

a stochastic, security-constrained, combined unit-commitment and multiperiod op-

timal power ﬂow problem with locational contingency and load-following reserves,

ramping costs and constraints, deferrable demands, lossy storage resources and un-

certain renewable generation.

MOST is documented separately from the main Matpower package in its own

manual, the MOST User’s Manual.

1http://pserc.org/

2http://www.pserc.cornell.edu/powerweb/

10

1.2 License and Terms of Use

Beginning with version 5.1, the code in Matpower is distributed under the 3-clause

BSD license3[5]. The full text of the license can be found in the LICENSE ﬁle at

the top level of the distribution or at http://www.pserc.cornell.edu/matpower/

LICENSE.txt and reads as follows.

Copyright (c) 1996-2016, Power Systems Engineering Research Center

(PSERC) and individual contributors (see AUTHORS file for details).

All rights reserved.

Redistribution and use in source and binary forms, with or without

modification, are permitted provided that the following conditions

are met:

1. Redistributions of source code must retain the above copyright

notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright

notice, this list of conditions and the following disclaimer in the

documentation and/or other materials provided with the distribution.

3. Neither the name of the copyright holder nor the names of its

contributors may be used to endorse or promote products derived from

this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS

"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT

LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS

FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE

COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,

INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,

BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;

LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER

CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT

LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN

ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE

POSSIBILITY OF SUCH DAMAGE.

3Versions 4.0 through 5.0 of Matpower were distributed under version 3.0 of the GNU General

Public License (GPL) [6] with an exception added to clarify our intention to allow Matpower to

interface with Matlab as well as any other Matlab code or MEX-ﬁles a user may have installed,

regardless of their licensing terms. The full text of the GPL can be found at http://www.gnu.

org/licenses/gpl-3.0.txt.Matpower versions prior to version 4 had their own license.

11

Please note that the Matpower case ﬁles distributed with Matpower are not

covered by the BSD license. In most cases, the data has either been included with

permission or has been converted from data available from a public source.

1.3 Citing Matpower

While not required by the terms of the license, we do request that publications derived

from the use of Matpower explicitly acknowledge that fact by citing reference [1].

R. D. Zimmerman, C. E. Murillo-S´anchez, and R. J. Thomas, “Matpower: Steady-

State Operations, Planning and Analysis Tools for Power Systems Research and Ed-

ucation,” Power Systems, IEEE Transactions on, vol. 26, no. 1, pp. 12–19, Feb. 2011.

DOI: 10.1109/TPWRS.2010.2051168

Similarly, we request that publications derived from the use of MOST explicitly

acknowledge that fact by citing both the main Matpower reference above and

reference [2]:

C. E. Murillo-S´anchez, R. D. Zimmerman, C. L. Anderson, and R. J. Thomas, “Secure

Planning and Operations of Systems with Stochastic Sources, Energy Storage and

Active Demand,” Smart Grid, IEEE Transactions on, vol. 4, no. 4, pp. 2220–2229,

Dec. 2013.

DOI: 10.1109/TSG.2013.2281001

1.4 Matpower Development

Following the release of Matpower 6.0, the Matpower project moved to an open

development paradigm, hosted on the Matpower GitHub project page:

https://github.com/MATPOWER/matpower

The Matpower GitHub project hosts the public Git code repository as well

as a public issue tracker for handling bug reports, patches, and other issues and

contributions. There are separate GitHub hosted repositories and issue trackers for

Matpower, MOST, MIPS and the testing framework used by all of them, MP-Test,

all available from https://github.com/MATPOWER/.

12

2 Getting Started

2.1 System Requirements

To use Matpower 7.0b1 you will need:

•Matlab®version 7.3 (R2006b) or later4, or

•GNU Octave version 4 or later5

See Section 2.1 in the MOST User’s Manual for any requirements speciﬁc to

MOST.

For the hardware requirements, please refer to the system requirements for the

version of Matlab6or Octave that you are using. If the Matlab Optimization

Toolbox is installed as well, Matpower enables an option to use it to solve optimal

power ﬂow problems.

In this manual, references to Matlab usually apply to Octave as well. At the

time of writing, none of the optional MEX-based Matpower packages have been

built for Octave, but Octave does typically include GLPK.

4Matlab is available from The MathWorks, Inc. (http://www.mathworks.com/). Mat-

power 5 and Matpower 6 required Matlab 7 (R14), Matpower 4 required Matlab 6.5 (R13),

Matpower 3.2 required Matlab 6 (R12), Matpower 3.0 required Matlab 5 and Matpower 2.0

and earlier required only Matlab 4. Matlab is a registered trademark of The MathWorks, Inc.

5GNU Octave [4] is free software, available online at http://www.gnu.org/software/octave/.

Matpower 5 and Matpower 6 required Octave 3.4, Matpower 4 required Octave 3.2.

6http://www.mathworks.com/support/sysreq/previous_releases.html

13

2.2 Getting Matpower

You can either download an oﬃcial versioned release or you can obtain the current

development version, which we also attempt to keep stable enough for everyday use.

The development version includes new features and bug ﬁxes added since the last

versioned release.

2.2.1 Versioned Releases

•Download the ZIP ﬁle of the latest oﬃcial versioned release from the Mat-

power website7.

–Go to the Matpower website.

–Click the Download Now button.

2.2.2 Current Development Version

There are also two options for obtaining the most recent development version of

Matpower from the master branch on GitHub.

1. Clone the Matpower repository from GitHub.

•From the command line:

–git clone https://github.com/MATPOWER/matpower.git

•Or, from the Matpower GitHub repository page8:

–Click the green Clone or download button, then Open in Desk-

top.

Use this option if you want to be able to easily update to the current development

release, with the latest bug ﬁxes and new features, using a simple git pull

command, or if you want to help with testing or development. This requires

that you have a Git client9(GUI or command-line) installed.

2. Download a zip ﬁle of the Matpower repository from GitHub.

•Go to the Matpower GitHub repository page.

7http://www.pserc.cornell.edu/matpower/

8https://github.com/MATPOWER/matpower

9See https://git-scm.com/downloads for information on downloading Git clients.

14

•Click the green Clone or download button, then Download ZIP.

Use this option if you need features or ﬁxes introduced since the latest versioned

release, but you do not have access to or are not ready to begin using Git (but

don’t be afraid to give Git a try).10

See CONTRIBUTING.md for information on how to get a local copy of your own

Matpower fork, if you are interesting in contributing your own code or modiﬁca-

tions.

10A good place to start is https://git-scm.com.

15

2.3 Installation

Installation and use of Matpower requires familiarity with the basic operation of

Matlab or Octave. Make sure you follow the installation instructions for the version

of Matpower you are installing. The process was simpliﬁed with an install script

following version 6.0.

Step 1: Get a copy of Matpower as described above.

Clone the repository or download and extract the ZIP ﬁle of the Mat-

power distribution and place the resulting directory in the location of

your choice and call it anything you like.11 We will use <MATPOWER>

as a placeholder to denote the path to this directory (the one containing

install matpower.m). The ﬁles in <MATPOWER>should not need to be

modiﬁed, so it is recommended that they be kept separate from your own

code.

Step 2: Run the installer.

Open Matlab or Octave and change to the <MATPOWER>directory.

Run the installer and follow the directions to add the required directories

to your Matlab or Octave path, by typing:

install matpower

Step 3: That’s it. There is no step 3.

But, if you chose not to have the installer run the test suite for you in step 2,

you can run it now to verify that Matpower is installed and functioning

properly, by typing:12

test matpower

The result should resemble the following, possibly including extra tests,

depending on the availablility of optional packages, solvers and extras.

11Do not place Matpower’s ﬁles in a directory named 'matlab'or 'optim'(both case-

insensitive), as these can cause Matlab’s built-in ver command to behave strangely in ways that

aﬀect Matpower.

12The MOST test suite is run separately by typing test most. See the MOST User’s Manual for

details.

16

>> test_matpower

t_test_fcns.............ok

t_nested_struct_copy....ok

t_feval_w_path..........ok

t_mpoption..............ok

t_loadcase..............ok

t_ext2int2ext...........ok

t_jacobian..............ok

t_hessian...............ok

t_margcost..............ok

t_totcost...............ok

t_modcost...............ok

t_hasPQcap..............ok

t_mplinsolve............ok (6 of 44 skipped)

t_mips..................ok

t_mips_pardiso..........ok (60 of 60 skipped)

t_qps_matpower..........ok (288 of 360 skipped)

t_miqps_matpower........ok (240 of 240 skipped)

t_pf....................ok

t_pf_radial.............ok

t_cpf...................ok

t_islands...............ok

t_opf_model.............ok

t_opf_model_legacy......ok

t_opf_default...........ok

t_opf_mips..............ok (278 of 1296 skipped)

t_opf_dc_mips...........ok

t_opf_dc_mips_sc........ok

t_opf_userfcns..........ok

t_opf_softlims..........ok

t_runopf_w_res..........ok

t_dcline................ok

t_get_losses............ok

t_load2disp.............ok

t_makePTDF..............ok

t_makeLODF..............ok

t_printpf...............ok

t_vdep_load.............ok

t_total_load............ok

t_scale_load............ok

t_apply_changes.........ok

t_psse..................ok

t_off2case..............ok

t_auction_mips..........ok

t_runmarket.............ok

All tests successful (7961 passed, 872 skipped of 8833)

Elapsed time 29.39 seconds. 17

2.4 Running a Simulation

The primary functionality of Matpower is to solve power ﬂow and optimal power

ﬂow (OPF) problems. This involves (1) preparing the input data deﬁning the all of

the relevant power system parameters, (2) invoking the function to run the simulation

and (3) viewing and accessing the results that are printed to the screen and/or saved

in output data structures or ﬁles.

2.4.1 Preparing Case Input Data

The input data for the case to be simulated are speciﬁed in a set of data matrices

packaged as the ﬁelds of a Matlab struct, referred to as a “Matpower case” struct

and conventionally denoted by the variable mpc. This struct is typically deﬁned in a

case ﬁle, either a function M-ﬁle whose return value is the mpc struct or a MAT-ﬁle

that deﬁnes a variable named mpc when loaded13. The main simulation routines,

whose names begin with run (e.g. runpf,runopf), accept either a ﬁle name or a

Matpower case struct as an input.

Use loadcase to load the data from a case ﬁle into a struct if you want to make

modiﬁcations to the data before passing it to the simulation.

>> mpc = loadcase(casefilename);

See also savecase for writing a Matpower case struct to a case ﬁle.

The structure of the Matpower case data is described a bit further in Section 3.1

and the full details are documented in Appendix Band can be accessed at any time

via the command help caseformat. The Matpower distribution also includes many

example case ﬁles listed in Table D-18.

2.4.2 Solving the Case

The solver is invoked by calling one of the main simulation functions, such as runpf

or runopf, passing in a case ﬁle name or a case struct as the ﬁrst argument. For

example, to run a simple Newton power ﬂow with default options on the 9-bus system

deﬁned in case9.m, at the Matlab prompt, type:

>> runpf('case9');

13This describes version 2 of the Matpower case format, which is used internally and is the

default. The version 1 format, now deprecated, but still accessible via the loadcase and savecase

functions, deﬁnes the data matrices as individual variables rather than ﬁelds of a struct, and some

do not include all of the columns deﬁned in version 2.

18

If, on the other hand, you wanted to load the 30-bus system data from case30.m,

increase its real power demand at bus 2 to 30 MW, then run an AC optimal power

ﬂow with default options, this could be accomplished as follows:

>> define_constants;

>> mpc = loadcase('case30');

>> mpc.bus(2, PD) = 30;

>> runopf(mpc);

The define constants in the ﬁrst line is simply a convenience script that deﬁnes a

number of variables to serve as named column indices for the data matrices. In this

example, it allows us to access the “real power demand” column of the bus matrix

using the name PD without having to remember that it is the 3rd column.

Other top-level simulation functions are available for running DC versions of

power ﬂow and OPF, for running an OPF with the option for Matpower to shut

down (decommit) expensive generators, etc. These functions are listed in Table D-2

in Appendix D.

2.4.3 Accessing the Results

By default, the results of the simulation are pretty-printed to the screen, displaying

a system summary, bus data, branch data and, for the OPF, binding constraint

information. The bus data includes the voltage, angle and total generation and load

at each bus. It also includes nodal prices in the case of the OPF. The branch data

shows the ﬂows and losses in each branch. These pretty-printed results can be saved

to a ﬁle by providing a ﬁlename as the optional 3rd argument to the simulation

function.

The solution is also stored in a results struct available as an optional return value

from the simulation functions. This results struct is a superset of the Matpower

case struct mpc, with additional columns added to some of the existing data ﬁelds

and additional ﬁelds. The following example shows how simple it is, after running a

DC OPF on the 118-bus system in case118.m, to access the ﬁnal objective function

value, the real power output of generator 6 and the power ﬂow in branch 51.

>> define_constants;

>> results = rundcopf('case118');

>> final_objective = results.f;

>> gen6_output = results.gen(6, PG);

>> branch51_flow = results.branch(51, PF);

19

Full documentation for the content of the results struct can be found in Sec-

tions 4.4 and 6.6.

2.4.4 Setting Options

Matpower has many options for selecting among the available solution algorithms,

controlling the behavior of the algorithms and determining the details of the pretty-

printed output. These options are passed to the simulation routines as a Matpower

options struct. The ﬁelds of the struct have names that can be used to set the

corresponding value via the mpoption function. Calling mpoption with no arguments

returns the default options struct, the struct used if none is explicitly supplied.

Calling it with a set of name and value pairs modiﬁes the default vector.

For example, the following code runs a power ﬂow on the 300-bus example in

case300.m using the fast-decoupled (XB version) algorithm, with verbose printing of

the algorithm progress, but suppressing all of the pretty-printed output.

>> mpopt = mpoption('pf.alg', 'FDXB', 'verbose', 2, 'out.all', 0);

>> results = runpf('case300', mpopt);

To modify an existing options struct, for example, to turn the verbose option oﬀ

and re-run with the remaining options unchanged, simply pass the existing options

as the ﬁrst argument to mpoption.

>> mpopt = mpoption(mpopt, 'verbose', 0);

>> results = runpf('case300', mpopt);

See Appendix Cor type:

>> help mpoption

for more information on Matpower’s options.

2.5 Documentation

There are four primary sources of documentation for Matpower. The ﬁrst is

this manual, which gives an overview of Matpower’s capabilities and structure

and describes the modeling and formulations behind the code. It can be found in

your Matpower distribution at <MATPOWER>/docs/MATPOWER-manual.pdf and the

latest version is always available at: http://www.pserc.cornell.edu/matpower/

MATPOWER-manual.pdf.

20

Secondly, the MOST User’s Manual describes MOST and its problem formulation,

features, data formats and options. It is located at <MATPOWER>/docs/MOST-manual.pdf

in your Matpower distribution and the latest version is also available online at:

http://www.pserc.cornell.edu/matpower/MOST-manual.pdf.

The Matpower Online Function Reference is the third source of documen-

tation, allowing you to view not only the help text, but also the code itself for

each Matpower function directly from a web browser. It is available online at:

http://www.pserc.cornell.edu/matpower/docs/ref

Last, but certainly not least, is the built-in help command. As with Matlab’s

built-in functions and toolbox routines, you can type help followed by the name of

a command or M-ﬁle to get help on that particular function. Nearly all of Mat-

power’s M-ﬁles have such documentation and this should be considered the main

reference for the calling options for each individual function. See Appendix Dfor a

list of Matpower functions.

As an example, the help for runopf looks like:

21

>> help runopf

RUNOPF Runs an optimal power flow.

[RESULTS, SUCCESS] = RUNOPF(CASEDATA, MPOPT, FNAME, SOLVEDCASE)

Runs an optimal power flow (AC OPF by default), optionally returning

a RESULTS struct and SUCCESS flag.

Inputs (all are optional):

CASEDATA : either a MATPOWER case struct or a string containing

the name of the file with the case data (default is 'case9')

(see also CASEFORMAT and LOADCASE)

MPOPT : MATPOWER options struct to override default options

can be used to specify the solution algorithm, output options

termination tolerances, and more (see also MPOPTION).

FNAME : name of a file to which the pretty-printed output will

be appended

SOLVEDCASE : name of file to which the solved case will be saved

in MATPOWER case format (M-file will be assumed unless the

specified name ends with '.mat')

Outputs (all are optional):

RESULTS : results struct, with the following fields:

(all fields from the input MATPOWER case, i.e. bus, branch,

gen, etc., but with solved voltages, power flows, etc.)

order - info used in external <-> internal data conversion

et - elapsed time in seconds

success - success flag, 1 = succeeded, 0 = failed

(additional OPF fields, see OPF for details)

SUCCESS : the success flag can additionally be returned as

a second output argument

Calling syntax options:

results = runopf;

results = runopf(casedata);

results = runopf(casedata, mpopt);

results = runopf(casedata, mpopt, fname);

results = runopf(casedata, mpopt, fname, solvedcase);

[results, success] = runopf(...);

Alternatively, for compatibility with previous versions of MATPOWER,

some of the results can be returned as individual output arguments:

[baseMVA, bus, gen, gencost, branch, f, success, et] = runopf(...);

Example:

results = runopf('case30');

See also RUNDCOPF, RUNUOPF. 22

3 Modeling

Matpower employs all of the standard steady-state models typically used for power

ﬂow analysis. The AC models are described ﬁrst, then the simpliﬁed DC models. In-

ternally, the magnitudes of all values are expressed in per unit and angles of complex

quantities are expressed in radians. Internally, all oﬀ-line generators and branches

are removed before forming the models used to solve the power ﬂow or optimal power

ﬂow problem. All buses are numbered consecutively, beginning at 1, and generators

are reordered by bus number. Conversions to and from this internal indexing is done

by the functions ext2int and int2ext. The notation in this section, as well as Sec-

tions 4and 6, is based on this internal numbering, with all generators and branches

assumed to be in-service. Due to the strengths of the Matlab programming lan-

guage in handling matrices and vectors, the models and equations are presented here

in matrix and vector form.

3.1 Data Formats

The data ﬁles used by Matpower are Matlab M-ﬁles or MAT-ﬁles which deﬁne

and return a single Matlab struct. The M-ﬁle format is plain text that can be edited

using any standard text editor. The ﬁelds of the struct are baseMVA,bus,branch,gen

and optionally gencost, where baseMVA is a scalar and the rest are matrices. In the

matrices, each row corresponds to a single bus, branch, or generator. The columns

are similar to the columns in the standard IEEE CDF and PTI formats. The number

of rows in bus,branch and gen are nb,nland ng, respectively. If present, gencost

has either ngor 2ngrows, depending on whether it includes costs for reactive power

or just real power. Full details of the Matpower case format are documented

in Appendix Band can be accessed from the Matlab command line by typing

help caseformat.

3.2 Branches

All transmission lines14, transformers and phase shifters are modeled with a com-

mon branch model, consisting of a standard πtransmission line model, with series

impedance zs=rs+jxsand total charging susceptance bc, in series with an ideal

phase shifting transformer. The transformer, whose tap ratio has magnitude τand

14This does not include DC transmission lines. For more information the handling of DC trans-

mission lines in Matpower, see Section 7.6.3.

23

phase shift angle θshift, is located at the from end of the branch, as shown in Fig-

ure 3-1. The parameters rs,xs,bc,τand θshift are speciﬁed directly in columns BR R

(3), BR X (4), BR B (5), TAP (9) and SHIFT (10), respectively, of the corresponding row

of the branch matrix.15

The complex current injections ifand itat the from and to ends of the branch,

respectively, can be expressed in terms of the 2 ×2 branch admittance matrix Ybr

and the respective terminal voltages vfand vt

if

it=Ybr vf

vt.(3.1)

With the series admittance element in the πmodel denoted by ys= 1/zs, the branch

admittance matrix can be written

Ybr =ys+jbc

21

τ2−ys1

τe−jθshift

−ys1

τejθshift ys+jbc

2.(3.2)

Figure 3-1: Branch Model

If the four elements of this matrix for branch iare labeled as follows:

Yi

br =yi

ff yi

ft

yi

tf yi

tt (3.3)

then four nl×1 vectors Yff ,Yft,Ytf and Ytt can be constructed, where the i-th element

of each comes from the corresponding element of Yi

br. Furthermore, the nl×nbsparse

15A value of zero in the TAP column indicates that the branch is a transmission line and not a

transformer, i.e. mathematically equivalent to a transformer with tap ratio set to 1.

24

connection matrices Cfand Ctused in building the system admittance matrices can

be deﬁned as follows. The (i, j)th element of Cfand the (i, k)th element of Ctare

equal to 1 for each branch i, where branch iconnects from bus jto bus k. All other

elements of Cfand Ctare zero.

3.3 Generators

A generator is modeled as a complex power injection at a speciﬁc bus. For generator i,

the injection is

si

g=pi

g+jqi

g.(3.4)

Let Sg=Pg+jQgbe the ng×1 vector of these generator injections. The MW and

MVAr equivalents (before conversion to p.u.) of pi

gand qi

gare speciﬁed in columns

PG (2) and QG (3), respectively of row iof the gen matrix. A sparse nb×nggenerator

connection matrix Cgcan be deﬁned such that its (i, j)th element is 1 if generator j

is located at bus iand 0 otherwise. The nb×1 vector of all bus injections from

generators can then be expressed as

Sg,bus =Cg·Sg.(3.5)

A generator with a negative injection can also be used to model a dispatchable load.

3.4 Loads

Constant power loads are modeled as a speciﬁed quantity of real and reactive power

consumed at a bus. For bus i, the load is

si

d=pi

d+jqi

d(3.6)

and Sd=Pd+jQddenotes the nb×1 vector of complex loads at all buses. The

MW and MVAr equivalents (before conversion to p.u.) of pi

dand qi

dare speciﬁed in

columns PD (3) and QD (4), respectively of row iof the bus matrix. These ﬁelds can

also take on negative quantities to represent ﬁxed (e.g. distributed) generation.

Constant impedance and constant current loads are not implemented directly,

but the constant impedance portions can be modeled as a shunt element described

below. Dispatchable loads are modeled as negative generators and appear as negative

values in Sg.

25

3.5 Shunt Elements

A shunt connected element such as a capacitor or inductor is modeled as a ﬁxed

impedance to ground at a bus. The admittance of the shunt element at bus iis given

as

yi

sh =gi

sh +jbi

sh (3.7)

and Ysh =Gsh +jBsh denotes the nb×1 vector of shunt admittances at all buses.

The parameters gi

sh and bi

sh are speciﬁed in columns GS (5) and BS (6), respectively,

of row iof the bus matrix as equivalent MW (consumed) and MVAr (injected) at a

nominal voltage magnitude of 1.0 p.u and angle of zero.

A shunt element can also be used to model a constant impedance load and, though

correctly, Matpower does not currently report these quantities as “load”.

3.6 Network Equations

For a network with nbbuses, all constant impedance elements of the model are

incorporated into a complex nb×nbbus admittance matrix Ybus that relates the

complex nodal current injections Ibus to the complex node voltages V:

Ibus =YbusV. (3.8)

Similarly, for a network with nlbranches, the nl×nbsystem branch admittance

matrices Yfand Ytrelate the bus voltages to the nl×1 vectors Ifand Itof branch

currents at the from and to ends of all branches, respectively:

If=YfV(3.9)

It=YtV. (3.10)

If [ ·] is used to denote an operator that takes an n×1 vector and creates the

corresponding n×ndiagonal matrix with the vector elements on the diagonal, these

system admittance matrices can be formed as follows:

Yf= [Yff ]Cf+ [Yf t]Ct(3.11)

Yt= [Ytf ]Cf+ [Ytt]Ct(3.12)

Ybus =Cf

TYf+Ct

TYt+ [Ysh].(3.13)

The current injections of (3.8)–(3.10) can be used to compute the corresponding

26

complex power injections as functions of the complex bus voltages V:

Sbus(V) = [V]I∗

bus = [V]Y∗

busV∗(3.14)

Sf(V) = [CfV]I∗

f= [CfV]Y∗

fV∗(3.15)

St(V) = [CtV]I∗

t= [CtV]Y∗

tV∗.(3.16)

The nodal bus injections are then matched to the injections from loads and generators

to form the AC nodal power balance equations, expressed as a function of the complex

bus voltages and generator injections in complex matrix form as

gS(V, Sg) = Sbus(V) + Sd−CgSg= 0.(3.17)

3.7 DC Modeling

The DC formulation [17] is based on the same parameters, but with the following

three additional simplifying assumptions.

•Branches can be considered lossless. In particular, branch resistances rsand

charging capacitances bcare negligible:

ys=1

rs+jxs≈1

jxs

, bc≈0.(3.18)

•All bus voltage magnitudes are close to 1 p.u.

vi≈ejθi.(3.19)

•Voltage angle diﬀerences across branches are small enough that

sin(θf−θt−θshift)≈θf−θt−θshift.(3.20)

Substituting the ﬁrst set of assumptions regarding branch parameters from (3.18),

the branch admittance matrix in (3.2) approximates to

Ybr ≈1

jxs1

τ2−1

τe−jθshift

−1

τejθshift 1.(3.21)

Combining this and the second assumption with (3.1) yields the following approxi-

mation for if:

if≈1

jxs

(1

τ2ejθf−1

τe−jθshift ejθt)

=1

jxsτ(1

τejθf−ej(θt+θshift)).(3.22)

27

The approximate real power ﬂow is then derived as follows, ﬁrst applying (3.19) and

(3.22), then extracting the real part and applying (3.20).

pf=<{sf}

=<vf·i∗

f

≈ <ejθf·j

xsτ(1

τe−jθf−e−j(θt+θshift))

=<j

xsτ1

τ−ej(θf−θt−θshift)

=<1

xsτsin(θf−θt−θshift) + j1

τ−cos(θf−θt−θshift)

≈1

xsτ(θf−θt−θshift) (3.23)

As expected, given the lossless assumption, a similar derivation for the power injec-

tion at the to end of the line leads to leads to pt=−pf.

The relationship between the real power ﬂows and voltage angles for an individual

branch ican then be summarized as

pf

pt=Bi

br θf

θt+Pi

shift (3.24)

where

Bi

br =bi1−1

−1 1 ,

Pi

shift =θi

shiftbi−1

1

and biis deﬁned in terms of the series reactance xi

sand tap ratio τifor branch ias

bi=1

xi

sτi.

For a shunt element at bus i, the amount of complex power consumed is

si

sh =vi(yi

shvi)∗

≈ejθi(gi

sh −jbi

sh)e−jθi

=gi

sh −jbi

sh.(3.25)

28

So the vector of real power consumed by shunt elements at all buses can be approx-

imated by

Psh ≈Gsh.(3.26)

With a DC model, the linear network equations relate real power to bus voltage

angles, versus complex currents to complex bus voltages in the AC case. Let the

nl×1 vector Bff be constructed similar to Yff , where the i-th element is biand let

Pf,shift be the nl×1 vector whose i-th element is equal to −θi

shiftbi. Then the nodal

real power injections can be expressed as a linear function of Θ, the nb×1 vector of

bus voltage angles

Pbus(Θ) = BbusΘ + Pbus,shift (3.27)

where

Pbus,shift = (Cf−Ct)TPf,shift.(3.28)

Similarly, the branch ﬂows at the from ends of each branch are linear functions of

the bus voltage angles

Pf(Θ) = BfΘ + Pf,shift (3.29)

and, due to the lossless assumption, the ﬂows at the to ends are given by Pt=−Pf.

The construction of the system Bmatrices is analogous to the system Ymatrices

for the AC model:

Bf= [Bff ] (Cf−Ct) (3.30)

Bbus = (Cf−Ct)TBf.(3.31)

The DC nodal power balance equations for the system can be expressed in matrix

form as

gP(Θ, Pg) = BbusΘ + Pbus,shift +Pd+Gsh −CgPg= 0 (3.32)

29

4 Power Flow

The standard power ﬂow or loadﬂow problem involves solving for the set of voltages

and ﬂows in a network corresponding to a speciﬁed pattern of load and generation.

Matpower includes solvers for both AC and DC power ﬂow problems, both of

which involve solving a set of equations of the form

g(x) = 0,(4.1)

constructed by expressing a subset of the nodal power balance equations as functions

of unknown voltage quantities.

All of Matpower’s solvers exploit the sparsity of the problem and, except for

Gauss-Seidel, scale well to very large systems. Currently, none of them include any

automatic updating of transformer taps or other techniques to attempt to satisfy

typical optimal power ﬂow constraints, such as generator, voltage or branch ﬂow

limits.

4.1 AC Power Flow

In Matpower, by convention, a single generator bus is typically chosen as a refer-

ence bus to serve the roles of both a voltage angle reference and a real power slack.

The voltage angle at the reference bus has a known value, but the real power gen-

eration at the slack bus is taken as unknown to avoid overspecifying the problem.

The remaining generator buses are typically classiﬁed as PV buses, with the values

of voltage magnitude and generator real power injection given. These are speciﬁed

in the VG (6) and PG (3) columns of the gen matrix, respectively. Since the loads Pd

and Qdare also given, all non-generator buses are classiﬁed as PQ buses, with real

and reactive injections fully speciﬁed, taken from the PD (3) and QD (4) columns of

the bus matrix. Let Iref ,IPV and IPQ denote the sets of bus indices of the reference

bus, PV buses and PQ buses, respectively. The bus type classiﬁcation is speciﬁed in

the Matpower case ﬁle in the BUS TYPE column (2) of the bus matrix. Any isolated

buses must be identiﬁed as such in this column as well.

In the traditional formulation of the AC power ﬂow problem, the power balance

equation in (3.17) is split into its real and reactive components, expressed as functions

of the voltage angles Θ and magnitudes Vmand generator injections Pgand Qg, where

the load injections are assumed constant and given:

gP(Θ, Vm, Pg) = Pbus(Θ, Vm) + Pd−CgPg= 0 (4.2)

gQ(Θ, Vm, Qg) = Qbus(Θ, Vm) + Qd−CgQg= 0.(4.3)

30

For the AC power ﬂow problem, the function g(x) from (4.1) is formed by taking

the left-hand side of the real power balance equations (4.2) for all non-slack buses

and the reactive power balance equations (4.3) for all PQ buses and plugging in the

reference angle, the loads and the known generator injections and voltage magnitudes:

g(x) = "g{i}

P(Θ, Vm, Pg)

g{j}

Q(Θ, Vm, Qg)#∀i∈ IPV ∪ IPQ

∀j∈ IPQ.(4.4)

The vector xconsists of the remaining unknown voltage quantities, namely the volt-

age angles at all non-reference buses and the voltage magnitudes at PQ buses:

x=θ{i}

v{j}

m∀i /∈ Iref

∀j∈ IPQ.(4.5)

This yields a system of nonlinear equations with npv + 2npq equations and un-

knowns, where npv and npq are the number of PV and PQ buses, respectively. After

solving for x, the remaining real power balance equation can be used to compute

the generator real power injection at the slack bus. Similarly, the remaining npv + 1

reactive power balance equations yield the generator reactive power injections.

Matpower includes four diﬀerent algorithms for solving the general AC power

ﬂow problem.16 The default solver is based on a standard Newton’s method [8] using a

polar form and a full Jacobian updated at each iteration. Each Newton step involves

computing the mismatch g(x), forming the Jacobian based on the sensitivities of

these mismatches to changes in xand solving for an updated value of xby factorizing

this Jacobian. This method is described in detail in many textbooks.

Also included are solvers based on variations of the fast-decoupled method [9],

speciﬁcally, the XB and BX methods described in [10]. These solvers greatly reduce

the amount of computation per iteration, by updating the voltage magnitudes and

angles separately based on constant approximate Jacobians which are factored only

once at the beginning of the solution process. These per-iteration savings, however,

come at the cost of more iterations.

The fourth algorithm is the standard Gauss-Seidel method from Glimm and

Stagg [11]. It has numerous disadvantages relative to the Newton method and is

included primarily for academic interest.

By default, the AC power ﬂow solvers simply solve the problem described above,

ignoring any generator limits, branch ﬂow limits, voltage magnitude limits, etc. How-

ever, there is an option (pf.enforce q lims) that allows for the generator reactive

16Three more that are speciﬁc to radial networks typical of distribution systems are described in

Section 4.3.

31

power limits to be respected at the expense of the voltage setpoint. This is done

in a rather brute force fashion by adding an outer loop around the AC power ﬂow

solution. If any generator has a violated reactive power limit, its reactive injection is

ﬁxed at the limit, the corresponding bus is converted to a PQ bus and the power ﬂow

is solved again. This procedure is repeated until there are no more violations. Note

that this option is based solely on the QMAX and QMIN parameters for the generator,

from columns 4 and 5 of the gen matrix, and does not take into account the trape-

zoidal generator capability curves described in Section 6.4.3 and specifed in columns

PC1–QC2MAX (11–16). Note also that this option aﬀects generators even if the bus

they are attached to is already of type PQ.

4.2 DC Power Flow

For the DC power ﬂow problem [17], the vector xconsists of the set of voltage angles

at non-reference buses

x=θ{i},∀i /∈ Iref (4.6)

and (4.1) takes the form

Bdcx−Pdc = 0 (4.7)

where Bdc is the (nb−1) ×(nb−1) matrix obtained by simply eliminating from Bbus

the row and column corresponding to the slack bus and reference angle, respectively.

Given that the generator injections Pgare speciﬁed at all but the slack bus, Pdc can

be formed directly from the non-slack rows of the last four terms of (3.32).

The voltage angles in xare computed by a direct solution of the set of linear

equations. The branch ﬂows and slack bus generator injection are then calculated

directly from the bus voltage angles via (3.29) and the appropriate row in (3.32),

respectively.

4.3 Distribution Power Flow

Distribution systems are diﬀerent from transmission systems in a number of respects,

such as the xs/rsbranch ratio, magnitudes of xsand rsand most importantly the

typically radial structure. Due to these diﬀerences, a number of power ﬂow solu-

tion methods have been developed to account for the speciﬁc nature of distribution

systems and most widely used are the backward/forward sweep methods [12,13].

Matpower includes an additional three AC power ﬂow methods that are speciﬁc

to radial networks.

32

4.3.1 Radial Power Flow

When solving radial distribution networks it is practical to number branches with

numbers that are equal to the receiving bus numbers. An example is given in Fig-

ure 4-1, where branches are drawn black. Furthermore, the oriented branch ordering

[14] oﬀers a possibility for fast and eﬃcient backward/forward sweeps. All branches

are always oriented from the sending bus to the receiving bus and the only re-

quirement is that the sending bus number should be smaller than the receiving bus

number. This means that i<kfor branch i–k. The indices of the sending nodes of

branches are stored in vector Fsuch that i=fk.

Loop 1

Loop 2

011223

3

4477PV2

6

6

5PV1

5

Figure 4-1: Oriented Ordering

As usual, the supply bus (slack bus) is given index 1, meaning that branch indices

should go from 2 to nbwhich is the number of buses in the network. Introducing a

ﬁctitious branch with index 1 and zero impedance, given with dashed black line in

Figure 4-1, the number of branches nlbecomes equal to the number of buses nb.

For the example of Figure 4-1 vector Fis the following

F=0123324T,

and it oﬀers an easy way to follow the path between any bus and bus 0. If we

consider bus 4, the path to the slack bus consists of following branches: branch 4,

since considered bus is their receiving bus; branch 3, since f4= 3; branch 2, since

f3= 2; and branch 1, since f2= 1.

The representation of branch k, connecting buses iand k, is given in Figure 4-2,

where it is modeled with its serial impedance zk

s. At both ends there are load

demands si

dand sk

d, and shunt admittances yi

dand yk

dcomprised of admittances due

capacitance of all lines and shunt elements connected to buses iand k

yk

d=jXbk

lines +Xbk

shunt.

33

isk

fzk

ssk

tk

si

dyi

d

+

−

vi

sk

d

yk

d

+

−

vk

Figure 4-2: Branch Representation: branch kbetween buses i(sending) and k(re-

ceiving) and load demand and shunt admittances at both buses

4.3.2 Current Summation Method

The voltage calculation procedure with the Current Summation Method is performed

in 5 steps as follows [12,13].

1. Set all voltages to 1 p.u. (ﬂat start). Set iteration count ν= 1.

2. Set branch current ﬂow equal to the sum of current of the demand at receiving

end (sk

d) and the current drawn in the admittance (yk

d) connected to bus k

jk

b=sk

d

vk∗

+yk

d·vk, k = 1,2, . . . , nb.(4.8)

3. Backward sweep: Perform current summation, starting from the branch with

the biggest index and heading towards the branch whose index is equal to 1.

The current of branch kis added to the current of the branch whose index is

equal to i=fk.

ji

b,new =ji

b+jk

b, k =nl, nl−1,...,2 (4.9)

4. Forward sweep: The receiving end bus voltages are calculated with known

branch currents and sending bus voltages.

vk=vi−zk

s·jk

b, k = 2,3, . . . , nl(4.10)

5. Compare voltages in iteration νwith the corresponding ones from iteration

ν−1. If the maximum diﬀerence in magnitude is less than the speciﬁed toler-

ance

max

i=1...nbvν

i−vν−1

i< ε (4.11)

the procedure is ﬁnished. Otherwise go to step 2.

34

4.3.3 Power Summation Method

The voltage calculation procedure with the Power Summation Method is performed

in 5 steps as follows [14].

1. Set all voltages to 1 p.u. (ﬂat start). Set iteration count ν= 1.

2. Set receiving end branch ﬂow equal to the sum of the demand at receiving end

(sk

d) and the power drawn in the admittance (yk

d) connected to bus k

sk

t=sk

d+yk

d∗

v2

k

, k = 1,2, . . . , nb.(4.12)

3. Backward sweep: Calculate sending end branch power ﬂows as a sum of re-

ceiving end branch power ﬂows and branch losses via (4.13). Perform power

summation, starting from the branch with the biggest index and heading to-

wards the branch whose index is equal to 1. The sending power of branch kis

added to the receiving power of the branch whose index is equal to i=fkas

in (4.14).

sk

f=sk

t+zk

s·

sk

t

vk

2

k=nl, nl−1,...,2 (4.13)

si

t,new =si

t+sk

fk=nl, nl−1,...,2 (4.14)

4. Forward sweep: The receiving end bus voltages are calculated with known

sending powers and voltages.

vk=vi−zk

s· sk

f

vi!∗

k= 2,3, . . . , nl(4.15)

5. Compare voltages in iteration νwith the corresponding ones from iteration

ν−1, using (4.11). If the maximum diﬀerence in magnitude is less than the

speciﬁed tolerance the procedure is ﬁnished. Otherwise go to step 2.

4.3.4 Admittance Summation Method

For each node, besides the known admittance yk

d, we deﬁne the admittance yk

eas the

driving point admittance of the part of the network fed by node k, including shunt

35

admittance yk

d. We also deﬁne an equivalent current generator jk

efor the part of the

network fed by node k. The current of this generator consists of all load currents fed

by node k. The process of calculation of bus voltages with the admittance summation

method consists of the following 5 steps [15].

1. Set all voltages to 1 p.u. (ﬂat start). Set iteration count ν= 1. Set initial

values

yk

e=yk

d, k = 1,2, . . . , nb.(4.16)

2. For each node k, calculate equivalent admittance yk

e. Perform admittance sum-

mation, starting from the branch with the biggest index and heading towards

the branch whose index is equal to 1. The driving point admittance of branch k

is added to the driving point admittance of the branch whose index is equal to

i=fkas in (4.18).

dk

b=1

1 + zk

s·yk

e

k=nl, nl−1,...,2 (4.17)

yi

e,new =yi

e+dk

b·yk

ek=nl, nl−1,...,2 (4.18)

3. Backward sweep: For each node kcalculate equivalent current generator jk

e, ﬁrst

set it equal to load current jk

dand perform current summation over equivalent

admittances using factor dk

bas in (4.17).

jk

e=jk

d=sk

d

vk∗

k=nl, nl−1,...,2 (4.19)

ji

e,new =ji

e+dk

b·jk

ek=nl, nl−1,...,2 (4.20)

4. Forward sweep: The receiving end bus voltages are calculated with known

equivalent current generators and sending bus voltages.

vk=dk

b·(vi−zk

s·jk

e)k= 2,3, . . . , nl(4.21)

5. Compare voltages in iteration νwith the corresponding ones from iteration

ν−1, using (4.11). If the maximum diﬀerence in magnitude is less than the

speciﬁed tolerance the procedure is ﬁnished. Otherwise go to step 3.

36

4.3.5 Handling PV Buses

The methods explained in the previous three subsections are applicable to radial

networks without loops and PV buses. These methods can be used to solve the

power ﬂow problem in weakly meshed networks if a compensation procedure based on

Thevenin equivalent impedance matrix is added [12,13]. In [14] a voltage correction

procedure is added to the process.

The list of branches is expanded by a set of npv ﬁctitious links, corresponding to

the PV nodes. Each of these links starts at the slack bus and ends at a corresponding

PV bus, thus forming a loop in the network. A ﬁctitious link going to bus kis

represented by a voltage generator with a voltage magnitude equal to the speciﬁed

voltage magnitude at bus k. Its phase angle is equal to the calculated phase angle

at bus k.

A loop impedance matrix Zlis formed for the loops made by the ﬁctitious links

and it has the following properties

•Element zmm

lis equal to the sum of branch impedances oﬀ all branches related

to loop m,

•Element zmk

lis equal to the sum of branch impedances of mutual branches of

loops mand k.

As an illustration, in Figure 4-1 there a two PV generators at buses 5 and 7. The

ﬁctitious links and loops orientation are drawn in red. The Thevenin matrix for this

case is

Zl=z2

b+z3

b+z5

bz2

b+z3

b

z2

b+z3

bz2

b+z3

b+z4

b+z7

b.

First column elements are equal to PV bus voltages when the current injection

at bus 5 is 1 p.u. and v1= 0. Bus voltages can be calculated with the current

summation method in a single iteration. By repeating the procedure for bus 7 one

can calculate the elements of second column.

By breaking all links the network becomes radial [13] and the three backward/forward

sweep methods are applicable. Since all link are ﬁctitious, only the injected reactive

power at their receiving bus mis determined [14] by the following equation

∆qm

pv = ∆dm

pv

v2

m

<{vm},(4.22)

which is practically an increment in reactive power injection of the corresponding

PV generator for the current iteration.

37

The incremental changes of the imaginary part of PV generator current ∆dm

pv can

be obtained by solving the matrix equation

={Zl} · ∆Dpv = ∆Epv,(4.23)

where

∆em

pv =vm

g

|vm|−1· <{vm}, m = 1,2, . . . , npv.(4.24)

In order to ensure 90◦phase diﬀerence between voltage and current at PV gen-

erators in [16] it was suggested to calculate the real part of PV generator current

as

∆cm

pv = ∆dm

pv ={vm}

<{vm}.(4.25)

In such a way the PV generator will inject purely reactive power, as it is supposed

to do. Its active power is added before as a negative load.

Before proceeding with the next iteration, the bus voltage corrections are cal-

culated. In order to do that, the radial network is solved by applying incremental

current changes ∆Ipv = ∆Cpv +j∆Dpv at the PV buses as excitations and setting

v1= 0. After the backward/forward sweep is performed with the current summation

method, the voltage corrections at all buses are known. They are added to the latest

voltages in order to obtain the new bus voltages, which are used in the next iteration

[14].

4.4 runpf

In Matpower, a power ﬂow is executed by calling runpf with a case struct or case

ﬁle name as the ﬁrst argument (casedata). In addition to printing output to the

screen, which it does by default, runpf optionally returns the solution in a results

struct.

>> results = runpf(casedata);

The results struct is a superset of the input Matpower case struct mpc, with some

additional ﬁelds as well as additional columns in some of the existing data ﬁelds.

The solution values are stored as shown in Table 4-1.

Additional optional input arguments can be used to set options (mpopt) and

provide ﬁle names for saving the pretty printed output (fname) or the solved case

data (solvedcase).

38

Table 4-1: Power Flow Results

name description

results.success success ﬂag, 1 = succeeded, 0 = failed

results.et computation time required for solution

results.iterations number of iterations required for solution

results.order see ext2int help for details on this ﬁeld

results.bus(:, VM)†bus voltage magnitudes

results.bus(:, VA) bus voltage angles

results.gen(:, PG) generator real power injections

results.gen(:, QG)†generator reactive power injections

results.branch(:, PF) real power injected into “from” end of branch

results.branch(:, PT) real power injected into “to” end of branch

results.branch(:, QF)†reactive power injected into “from” end of branch

results.branch(:, QT)†reactive power injected into “to” end of branch

†AC power ﬂow only.

>> results = runpf(casedata, mpopt, fname, solvedcase);

The options that control the power ﬂow simulation are listed in Table 4-2 and those

controlling the output printed to the screen in Table 4-3.

By default, runpf solves an AC power ﬂow problem using a standard Newton’s

method solver. To run a DC power ﬂow, the model option must be set to 'DC'. For

convenience, Matpower provides a function rundcpf which is simply a wrapper

that sets the model option to 'DC'before calling runpf.

Internally, the runpf function does a number of conversions to the problem data

before calling the appropriate solver routine for the selected power ﬂow algorithm.

This external-to-internal format conversion is performed by the ext2int function,

described in more detail in Section 7.3.1, and includes the elimination of out-of-service

equipment, the consecutive renumbering of buses and the reordering of generators

by increasing bus number. All computations are done using this internal indexing.

When the simulation has completed, the data is converted back to external format

by int2ext before the results are printed and returned.

4.5 Linear Shift Factors

The DC power ﬂow model can also be used to compute the sensitivities of branch

ﬂows to changes in nodal real power injections, sometimes called injection shift factors

(ISF) or generation shift factors [17]. These nl×nbsensitivity matrices, also called

39

Table 4-2: Power Flow Options

name default description

model 'AC'AC vs. DC modeling for power ﬂow and OPF formulation

'AC'– use AC formulation and corresponding alg options

'DC'– use DC formulation and corresponding alg options

pf.alg 'NR'AC power ﬂow algorithm:

'NR'– Newtons’s method

'FDXB'– Fast-Decoupled (XB version)

'FDBX'– Fast-Decouple (BX version)

'GS'– Gauss-Seidel

'PQSUM'– Power Summation (radial networks only)

'ISUM'– Current Summation (radial networks only)

'YSUM'– Admittance Summation (radial networks only)

pf.tol 10−8termination tolerance on per unit P and Q dispatch

pf.nr.max it 10 maximum number of iterations for Newton’s method

pf.nr.lin solver '' linear solver option for mplinsolve for computing Newton update step

(see mplinsolve for complete list of all options)

'' – default to '\'for small systems, 'LU3'for larger ones

'\'– built-in backslash operator

'LU'– explicit default LU decomposition and back substitution

'LU3'– 3 output arg form of lu, Gilbert-Peierls algorithm with

approximate minimum degree (AMD) reordering

'LU4'– 4 output arg form of lu, UMFPACK solver (same as

'LU')

'LU5'– 5 output arg form of lu, UMFPACK solver w/row scaling

pf.fd.max it 30 maximum number of iterations for fast-decoupled method

pf.gs.max it 1000 maximum number of iterations for Gauss-Seidel method

pf.radial.max it 20 maximum number of iterations for radial power ﬂow methods

pf.radial.vcorr 0 perform voltage correction procedure in distribution power ﬂow

0 – do not perform voltage correction

1 – perform voltage correction

pf.enforce q lims 0 enforce gen reactive power limits at expense of |Vm|

0 – do not enforce limits

1 – enforce limits, simultaneous bus type conversion

2 – enforce limits, one-at-a-time bus type conversion

power transfer distribution factors or PTDFs, carry an implicit assumption about

the slack distribution. If His used to denote a PTDF matrix, then the element in

row iand column j,hij , represents the change in the real power ﬂow in branch i

given a unit increase in the power injected at bus j,with the assumption that the

additional unit of power is extracted according to some speciﬁed slack distribution:

∆Pf=H∆Pbus.(4.26)

40

Table 4-3: Power Flow Output Options

name default description

verbose 1 amount of progress info to be printed

0 – print no progress info

1 – print a little progress info

2 – print a lot of progress info

3 – print all progress info

out.all -1 controls pretty-printing of results

-1 – individual ﬂags control what is printed

0 – do not print anything†

1 – print everything†

out.sys sum 1 print system summary (0 or 1)

out.area sum 0 print area summaries (0 or 1)

out.bus 1 print bus detail, includes per bus gen info (0 or 1)

out.branch 1 print branch detail (0 or 1)

out.gen 0 print generator detail (0 or 1)

out.force 0 print results even if success ﬂag = 0 (0 or 1)

out.suppress detail -1 suppress all output but system summary

-1 – suppress details for large systems (>500 buses)

0 – do not suppress any output speciﬁed by other ﬂags

1 – suppress all output except system summary section†

†Overrides individual ﬂags, but (in the case of out.suppress detail) not out.all = 1.

This slack distribution can be expressed as an nb×1 vector wof non-negative

weights whose elements sum to 1. Each element speciﬁes the proportion of the slack

taken up at each bus. For the special case of a single slack bus k,wis equal to the

vector ek. The corresponding PTDF matrix Hkcan be constructed by ﬁrst creating

the nl×(nb−1) matrix

e

Hk=e

Bf·B−1

dc (4.27)

then inserting a column of zeros at column k. Here e

Bfand Bdc are obtained from Bf

and Bbus, respectively, by eliminating their reference bus columns and, in the case

of Bdc, removing row kcorresponding to the slack bus.

The PTDF matrix Hw, corresponding to a general slack distribution w, can be

obtained from any other PTDF, such as Hk, by subtracting Hk·wfrom each column,

equivalent to the following simple matrix multiplication:

Hw=Hk(I−w·1T).(4.28)

These same linear shift factors may also be used to compute sensitivities of branch

ﬂows to branch outages, known as line outage distribution factors or LODFs [18].

41

Given a PTDF matrix Hw, the corresponding nl×nlLODF matrix Lcan be con-

structed as follows, where lij is the element in row iand column j, representing the

change in ﬂow in branch i(as a fraction of the initial ﬂow in branch j) for an outage

of branch j.

First, let Hrepresent the matrix of sensitivities of branch ﬂows to branch endpoint

injections, found by multplying the PTDF matrix by the node-branch incidence

matrix:

H=Hw(Cf−Ct)T.(4.29)

Here the individual elements hij represent the sensitivity of ﬂow in branch iwith

respect to injections at branch jendpoints, corresponding to a simulated increase in

ﬂow in branch j. Then lij can be expressed as

lij =

hij

1−hjj

i6=j

−1i=j.

(4.30)

Matpower includes functions for computing both the DC PTDF matrix and

the corresponding LODF matrix for either a single slack bus kor a general slack

distribution vector w. See the help for makePTDF and makeLODF and Sections 9.5.5

and 9.5.4, respectively, for details.

42

5 Continuation Power Flow

Continuation methods or branch tracing methods are used to trace a curve given an

initial point on the curve. These are also called predictor-corrector methods since

they involve the prediction of the next solution point and correcting the prediction

to get the next point on the curve.

Consider a system of nnonlinear equations g(x) = 0, x ∈Rn. By adding a

continuation parameter λand one more equation to the system, xcan be traced by

varying λ. The resulting system f(x, λ) = 0 has n+ 1 dimensions. The additional

equation is a parameterized equation which identiﬁes the location of the current

solution with respect to the previous or next solution.

The continuation process can be diagrammatically shown by (5.1).

xj, λjP redictor

−−−−−→ (ˆxj+1,ˆ

λj+1)Corrector

−−−−−→ xj+1, λj+1(5.1)

where, (xj, λj) represents the current solution at step j, (ˆxj+1,ˆ

λj+1) is the predicted

solution for the next step, and (xj+1, λj+1) is the next solution on the curve.

Continuation methods are employed in power systems to determine steady state

stability limits [19] in what is called a continuation power ﬂow17. The limit is de-

termined from a nose curve where the nose represents the maximum power transfer

that the system can handle given a power transfer schedule. To determine the steady

state loading limit, the basic power ﬂow equations

g(x) = P(x)−Pinj

Q(x)−Qinj = 0,(5.2)

are restructured as

f(x, λ) = g(x)−λb = 0 (5.3)

where x≡(Θ, Vm) and bis a vector of power transfer given by

b=Pinj

target −Pinj

base

Qinj

target −Qinj

base .(5.4)

The eﬀects of the variation of loading or generation can be investigated using the

continuation method by composing the bvector appropriately.

17Thanks to Shrirang Abhyankar, Rui Bo, and Alexander Flueck for contributions to Mat-

power’s continuation power ﬂow feature.

43

5.1 Parameterization

The values of (x, λ) along the solution curve can parameterized in a number of ways

[20,21]. Parameterization is a mathematical way of identifying each solution so that

the next solution or previous solution can be quantiﬁed. Matpower includes three

parameterization scheme options to quantify this relationship, detailed below, where

σis the continuation step size parameter.

•Natural parameterization simply uses λdirectly as the parameter, so the

new λis simply the previous value plus the step size.

pj(x, λ) = λ−λj−σj= 0 (5.5)

•Arc length parameterization results in the following relationship, where

the step size is equal to the 2-norm of the distance from one solution to the

next.

pj(x, λ) = X

i

(xi−xj

i)2+ (λ−λj)2−(σj)2= 0 (5.6)

•Pseudo arc length parameterization [23] is Matpower’s default pa-

rameterization scheme, where the next point (x, λ) on the solution curve is

constrained to lie in the hyperplane running through the predicted solution

(ˆxj+1,ˆ

λj+1) orthogonal to the tangent line from the previous corrected solution

(xj, λj). This relationship can be quantiﬁed by the function

pj(x, λ) = x

λ−xj

λjT

¯zj−σj= 0,(5.7)

where ¯zjis the normalized tangent vector at (xj, λj) and σjis the continuation

step size parameter.

5.2 Predictor

The predictor is used to produce an estimate for the next solution. The better the

prediction, the faster is the convergence to the solution point. Matpower uses a

tangent predictor for estimating the curve to the next solution. The tangent vector

zj=dx dλ j

Tat the current solution (xj, λj) is found by solving the linear

system fxfλ

pj−1

xpj−1

λzj=0

1.(5.8)

44

The matrix on the left-hand side is simply the standard power ﬂow Jacobian with an

additional column and row added. The extra column fλis simply the negative of the

power transfer vector band the extra row, required to make the system non-singular

and deﬁne the magnitude of zj, is the derivative of the the parameterization function

at the previous solution point pj−1.

The resulting tangent vector is then normalized

¯zj=zj

||zj||2

(5.9)

and used to compute the predicted approximation (ˆxj+1,ˆ

λj+1) to the next solution

(xj+1, λj+1) using ˆxj+1

ˆ

λj+1 =xj

λj+σj¯zj,(5.10)

where σjis the continuation step size.

5.3 Corrector

The corrector stage ﬁnds the next solution (xj+1, λj+1) by correcting the approxi-

mation estimated by the predictor (ˆxj+1,ˆ

λj+1). Newton’s method is used to ﬁnd

the next solution by solving the n+ 1 dimensional system in (5.11), where one of

(5.5)–(5.7) has been added as an additional constraint to the parameterized power

ﬂow equations of (5.3). f(x, λ)

pj(x, λ)= 0 (5.11)

5.4 Step Length Control

Step length control is a key element aﬀecting the computational eﬃciency of a contin-

uation method. It aﬀects the continuation method with two issues: (1) speed – how

fast the corrector converges to a speciﬁed accuracy, and (2) robustness – whether

the corrector converges to a true solution given a predicted point. Matpower’s

continuation power ﬂow can optionally use adaptive steps, where the step size σis

adjusted by a scaling factor αwithin the limits speciﬁed.

σj+1 =αjσj, σmin ≤σj+1 ≤σmax (5.12)

This scaling factor αjfor step jis limited to a maximum of 2 and is calculated from

an error estimation between the predicted and corrected solutions γjas follows,

αj= 1 + βcpf cpf

γj−1, αj≤2,(5.13)

45

where βcpf is a damping factor, cpf is a speciﬁed tolerance, and γjis given by

γj=

xj+1, λj+1−ˆxj+1,ˆ

λj+1

∞.(5.14)

5.5 Event Detection and Location

A continuation power ﬂow event is triggered when the value of one of the elements

of an event function changes sign from one continuation step to the next. The

event occurs at the point where the corresponding value of the event function passes

through zero. Matpower provides event functions to detect the location at which

the continuation power ﬂow reaches the following:

•a speciﬁed target λvalue

•the nose point

•the end of a full trace

•a generator reactive power limit

•a generator active power limit

•a bus voltage magnitude limit

•a branch ﬂow limit

Each event function is registered with an event name, a ﬂag indicating whether

or not the location of the event should be pinpointed, and if so, to within what

tolerance. For events that are to be located, when an event interval is detected, that

is, when an element of the event function value changes sign, Matpower adjusts

the continuation step size via a False Position or Regula Falsi method until it locates

the point of the zero-crossing to within the speciﬁed tolerance.

The detection of an event zero, or even an event interval, can be used to trigger

further actions. The CPF callback function capability was extended in Matpower

6.x to include the ability to handle events by including information about any event

intervals or zeros detected. For example, Matpower includes a callback that ﬁxes

the reactive output of a generator and converts its bus from PV to PQ when the

corresponding event function indicates that its reactive power limit has been reached.

Another responds to the detection of the nose point by signaling the termination of

the continuation. In fact, continuation power ﬂow termination for nose point, target

46

lambda or full trace modes are all based on CPF callback functions in conjunction

with event detection.

While Matpower does include a mechanism for supplying user deﬁned callback

functions, it does not yet have a corresponding mechanism for user speciﬁed event

functions.

5.6 runcpf

In Matpower, a continuation power ﬂow is executed by calling runcpf with two

Matpower cases (case structs or case ﬁle names) as the ﬁrst two arguments,

basecasedata and targetcasedata, respectively. The ﬁrst contains the base load-

ing/generation proﬁle while the second contains the target loading/generation pro-

ﬁle. In addition to printing output to the screen, which it does by default, runcpf

optionally returns the solution in a results struct.

>> results = runcpf(basecasedata, targetcasedata);

Additional optional input arguments can be used to set options (mpopt) and

provide ﬁle names for saving the pretty printed output (fname) or the solved case

data (solvedcase).

>> results = runcpf(basecasedata, targetcasedata, mpopt, fname, solvedcase);

The results struct is a superset of the input Matpower case struct mpc, with

some additional ﬁelds as well as additional columns in some of the existing data

ﬁelds. In addition to the solution values included in the results for a simple power

ﬂow, shown in Table 4-1 in Section 4.4, the following additional continuation power

ﬂow solution values are stored in the cpf ﬁeld as shown in Table 5-1.

The options that control the continuation power ﬂow simulation are listed in

Table 5-2. All the power ﬂow options for Newton’s method (tolerance, maximum

iterations) and for controlling the output on the screen (see Tables 4-2 and 4-3) are

also available with the continuation power ﬂow.

47

Table 5-1: Continuation Power Flow Results

name description

results.cpf CPF output struct whose content depends on any user callback functions,

where default contains ﬁelds:

done msg string with message describing cause of continuation termination

events(eidx) a structure array of size ne, where neis the number of events located, with

ﬁelds:

kcontinuation step number at which event was located

name name of event

idx index(es) of critical elements in corresponding event function, e.g. index

of generator reaching VAr limit

msg descriptive text detailing the event

iterations nsteps, number of continuation steps performed

lam 1×nvector of λvalues from correction steps†

lam hat 1×nvector of λvalues from prediction steps†

max lam maximum value of λfound in results.cpf.lam

steps 1×nvector of step sizes for each continuation step performed†

Vnb×nmatrix of complex bus voltages from correction steps†

V hat nb×nmatrix of complex bus voltages from prediction steps†

†nis one more than the number of continuation steps, i.e. nsteps + 1, so the ﬁrst element corresponds to the

starting point.

48

Table 5-2: Continuation Power Flow Options

name default description

cpf.parameterization 3 choice of parameterization

1 — natural

2 — arc length

3 — pseudo arc length

cpf.stop at 'NOSE'determines stopping criterion

'NOSE'— stop when nose point is reached

'FULL'— trace full nose curve

λstop — stop upon reaching target λvalue λstop

cpf.enforce p lims 0 enforce gen active power limits

0 — do not enforce limits

1 — enforce limits

cpf.enforce q lims 0 enforce gen reactive power limits at expense of Vm

0 — do not enforce limits

1 — enforce limits

cpf.enforce v lims 0 enforce bus voltage magnitude limits

0 — do not enforce limits

1 — enforce limits

cpf.enforce flow lims 0 enforce branch MVA ﬂow limits

0 — do not enforce limits

1 — enforce limits

cpf.step 0.05 default value for continuation power ﬂow step size σ

cpf.step min 10−4minimum allowed step size, σmin

cpf.step max 0.2 maximum allowed step size, σmax

cpf.adapt step 0 toggle adaptive step size feature

0 — adaptive step size disabled

1 — adaptive step size enabled

cpf.adapt step damping 0.7 damping factor βcpf from (5.13) for adaptive step sizing

cpf.adapt step tol 10−3tolerance cpf from (5.13) for adaptive step sizing

cpf.target lam tol 10−5tolerance for target lambda detection

cpf.nose tol 10−5tolerance for nose point detection (p.u.)

cpf.p lims tol 10−2tolerance for generator active power limit detection

(MW)

cpf.q lims tol 10−2tolerance for generator reactive power limit detection

(MVAr)

cpf.v lims tol 10−4tolerance for bus voltage magnitude limit detection (p.u)

cpf.flow lims tol 0.01 tolerance for branch ﬂow limit detection (MVA)

cpf.plot.level 0 control plotting of nose curve

0 — do not plot nose curve

1 — plot when completed

2 — plot incrementally at each iteration

3 — same as 2, with pause at each iteration

cpf.plot.bus empty index of bus whose voltage is to be plotted

cpf.user callback empty string or cell array of strings with names of user callback

functions†

†See help cpf default callback for details. 49

5.6.1 CPF Callback Functions

Matpower’s continuation power ﬂow provides a callback mechanism to give the

user access to the iteration process for executing custom code at each iteration,

for example, to implement custom incremental plotting of a PV nose curve or to

handle a detected event. This callback mechanism is used internally to handle de-

fault plotting functionality as well as to handle the standard CPF events. The

cpf default callback function, for example, is used to collect the λand Vresults

from each predictor and corrector iteration and optionally plot the PV nose curve.

The prototype for a CPF callback function is

function [nx, cx, done, rollback, evnts, cb_data, results] = ...

cpf_user_callback(k, nx, cx, px, done, rollback, evnts, ...

cb_data, cb_args, results)

and the input and output arguments are described in Tables 5-3 through 5-5 and in

the help for cpf default callback. Each registered CPF callback function is called in

three diﬀerent contexts, distinguished by the value of the ﬁrst argument kas follows:

1. initial – called with k= 0, without results input/output arguments, after base

power ﬂow, before ﬁrst CPF step.

2. iterations – called with k>0, without results input/output arguments, at

each iteration, after predictor-corrector step

3. ﬁnal – called with k<0, with results input/output arguments, after exiting

predictor-corrector loop, inputs identical to last iteration call, except knegated

The user can deﬁne their own callback functions which take the same form and

are called in the same contexts as cpf default callback. User callback functions are

speciﬁed via the Matpower option cpf.user callback. This option can be a string

containing the name of the callback function, or a struct with the following ﬁelds,

where all but the ﬁrst are optional:

•fcn - string with name of callback function

•priority - numerical value specifying callback priority,18 default = 20

•args - arbitrary value (any type) passed to the callback as cb args each time

it is invoked

Multiple user callbacks can be registered by assigning a cell array of such strings

and/or structs to the cpf.user callback option.

18See cpf register callback for details.

50

Table 5-3: Continuation Power Flow Callback Input Arguments

name description

kcontinuation step iteration count

nx next CPF state*

, corresponding to proposed next step

cx current CPF state*

, corresponding to most recent successful step

px previous CPF state*

, corresponding to last successful step prior to cx

done struct, with ﬂag to indicate CPF termination and reason, with ﬁelds:

.flag termination ﬂag, 1 →terminate, 0 →continue

.msg string containing reason for termination

rollback scalar ﬂag to indicate that the current step should be rolled back and

retried with a diﬀerent step size, etc.

evnts struct array listing any events detected for this step‡

cb data struct containing static data§

, with the following ﬁelds (all based on internal

indexing):

.mpc base Matpower case struct of base case

.mpc target Matpower case struct of target case

.Sbusb handle of function returning nb×1 vector of complex base case bus injec-

tions in p.u. and derivatives w.r.t. |V|

.Sbust handle of function returning nb×1 vector of complex target case bus

injections in p.u. and derivatives w.r.t. |V|

.Ybus bus admittance matrix

.Yf branch admittance matrix, “from” end of branches

.Yt branch admittance matrix, “to” end of branches

.pv list of indices of PV buses

.pq list of indices of PQ buses

.ref list of indices of reference buses

.idx pmax vector of gen indices of gens at PMAX

.mpopt Matpower options struct

cb args arbitrary data structure containing user-provided callback arguments

results initial value of output struct to be assigned to cpf ﬁeld of results struct

returned by runcpf

*See Table 5-5 for details of the CPF state.

‡See cpf detect events for details of the evnts ﬁeld.

§Please note that any callback that modiﬁes the underlying problem is responsible to update the contents of

cb data accordingly. E.g. converting a bus from PV to PQ requires updates to mpc base,mpc target,Sbusb,

Sbust,pv,pq, and possibly ref. So, cb data should only be thought of as static for a ﬁxed base and target

case pair.

51

Table 5-4: Continuation Power Flow Callback Output Arguments

name description

All are updated versions of the corresponding input arguments, see Table 5-3 for more details.

nx next CPF state*

, user state (cb ﬁeld) should be updated here if rollback

is false

cx current CPF state*

, may contain updated this step or this parm values

to be used if rollback is true

done callback may have requested termination and set the msg ﬁeld

rollback callback can request a rollback step, even if it was not indicated by an

event function†

evnts msg ﬁeld for a given event may be updated‡

cb data this data should only be modiﬁed if the underlying problem has been

changed (e.g. generator limit reached), in which case it should always

be followed by a step of zero length, i.e. set nx.this step to 0§

results updated version of results input argument to be assigned to cpf ﬁeld of

results struct returned by runcpf

*See Table 5-5 for details of the CPF state.

†In this case, the callback should also modify the step size or parameterization to be used for the re-try, by

setting the this step or this parm ﬁelds in cx.

‡See cpf detect events for details of the evnts ﬁeld.

§It is the job of any callback modifying cb data to ensure that all data in cb data is kept consistent.

Table 5-5: Continuation Power Flow State

name description

lam hat value of λafter predictor step

V hat vector of complex bus voltages after predictor step

lam value of λafter corrector step

Vvector of complex bus voltages after corrector step

znormalized tangent predictor, ¯z

default step default step size

default parm default parameterization*

this step step size for this step only

this parm parameterization*for this step only

step current step size

parm current parameterization*

events event log

cb user state for callbacks†, a user-deﬁned struct containing any information

a callback function would like to pass from one invokation to the next

ef cell array of event function values

*Corresponding to the cpf.parameterization option in Table 5-2.

†Replaces cb state from Matpower 5.

52

5.6.2 CPF Example

The following is an example of running a continuation power ﬂow using a version of

the 9-bus system, looking at increasing all loads by a factor of 2.5. This example

plots the nose curve shown in Figure 5-1.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Voltage at Bus 9

λ

Voltage Magnitude

Figure 5-1: Nose Curve of Voltage Magnitude at Bus 9

53

define_constants;

mpopt = mpoption('out.all', 0, 'verbose', 2);

mpopt = mpoption(mpopt, 'cpf.stop_at', 'FULL', 'cpf.step', 0.2);

mpopt = mpoption(mpopt, 'cpf.plot.level', 2);

mpcb = loadcase(t_case9_pfv2); % load base case

mpct = mpcb; % set up target case with

mpct.gen(:, [PG QG]) = mpcb.gen(:, [PG QG]) * 2.5; % increased generation

mpct.bus(:, [PD QD]) = mpcb.bus(:, [PD QD]) * 2.5; % and increased load

results = runcpf(mpcb, mpct, mpopt);

This should result in something like the following output to the screen.

MATPOWER Version 7.0b1, 31-Oct-2018 -- AC Continuation Power Flow

step 0 : lambda = 0.000, 4 Newton steps

step 1 : stepsize = 0.2 lambda = 0.181 2 corrector Newton steps

step 2 : stepsize = 0.2 lambda = 0.359 2 corrector Newton steps

step 3 : stepsize = 0.2 lambda = 0.530 2 corrector Newton steps

step 4 : stepsize = 0.2 lambda = 0.693 3 corrector Newton steps

step 5 : stepsize = 0.2 lambda = 0.839 3 corrector Newton steps

step 6 : stepsize = 0.2 lambda = 0.952 3 corrector Newton steps

step 7 : stepsize = 0.2 lambda = 0.988 3 corrector Newton steps

step 8 : stepsize = 0.2 lambda = 0.899 3 corrector Newton steps

step 9 : stepsize = 0.2 lambda = 0.776 3 corrector Newton steps

step 10 : stepsize = 0.2 lambda = 0.654 3 corrector Newton steps

step 11 : stepsize = 0.2 lambda = 0.533 2 corrector Newton steps

step 12 : stepsize = 0.2 lambda = 0.413 2 corrector Newton steps

step 13 : stepsize = 0.2 lambda = 0.294 2 corrector Newton steps

step 14 : stepsize = 0.2 lambda = 0.176 2 corrector Newton steps

step 15 : stepsize = 0.2 lambda = 0.060 2 corrector Newton steps

step 16a : stepsize = 0.2 lambda = -0.054 2 corrector Newton steps ^ ROLLBACK

step 16 : stepsize = 0.0604 lambda = 0.000 3 corrector Newton steps

CPF TERMINATION: Traced full continuation curve in 16 continuation steps

The results of the continuation power ﬂow are then found in the cpf ﬁeld of the

returned results struct.

54

>> results.cpf

ans =

V_hat: [9x17 double]

lam_hat: [1x17 double]

V: [9x17 double]

lam: [1x17 double]

steps: [1x17 double]

iterations: 16

max_lam: 0.9876

events: [1x1 struct]

done_msg: 'Traced full continuation curve in 16 continuation steps'

55

6 Optimal Power Flow

Matpower includes code to solve both AC and DC versions of the optimal power

ﬂow problem. The standard version of each takes the following form:

min

xf(x) (6.1)

subject to

g(x) = 0 (6.2)

h(x)≤0 (6.3)

xmin ≤x≤xmax .(6.4)

In both cases, the objective function f(x) consists of the polynomial cost of gener-

ator injections, the equality constraints g(x) are the power balance equations, the

inequality constraints h(x) are the branch ﬂow limits, and the xmin and xmax bounds

include reference bus angles, voltage magnitudes (for AC) and generator injections.

6.1 Standard AC OPF

The optimization vector xfor the standard AC OPF problem consists of the nb×1

vectors of voltage angles Θ and magnitudes Vmand the ng×1 vectors of generator

real and reactive power injections Pgand Qg.

x=

Θ

Vm

Pg

Qg

(6.5)

The objective function f(x) in (6.1) is simply a summation of individual polynomial

cost functions fi

Pand fi

Qof real and reactive power injections, respectively, for each

generator:

f(Pg, Qg) =

ng

X

i=1

fi

P(pi

g) + fi

Q(qi

g).(6.6)

The equality constraints in (6.2) are simply the full set of 2 ·nbnonlinear real and

reactive power balance equations from (4.2) and (4.3).

gP(Θ, Vm, Pg) = Pbus(Θ, Vm) + Pd−CgPg= 0 (6.7)

gQ(Θ, Vm, Qg) = Qbus(Θ, Vm) + Qd−CgQg= 0 (6.8)

56

The inequality constraints (6.3) consist of two sets of nlbranch ﬂow limits as non-

linear functions of the bus voltage angles and magnitudes, one for the from end and

one for the to end of each branch:

hf(Θ, Vm) = |Ff(Θ, Vm)| − Fmax ≤0 (6.9)

ht(Θ, Vm) = |Ft(Θ, Vm)| − Fmax ≤0.(6.10)

The ﬂows are typically apparent power ﬂows expressed in MVA, but can be real power

or current ﬂows, yielding the following three possible forms for the ﬂow constraints:

Ff(Θ, Vm) =

Sf(Θ, Vm),apparent power

Pf(Θ, Vm),real power

If(Θ, Vm),current

(6.11)

where Ifis deﬁned in (3.9), Sfin (3.15), Pf=<{Sf}and the vector of ﬂow limits

Fmax has the appropriate units for the type of constraint. It is likewise for Ft(Θ, Vm).

The values used by Matpower’s OPF for the ﬂow limits Fmax are speciﬁed in the

RATE A column (6) of the branch matrix,19 and the selection of ﬂow constraint type

in (6.11) is determined by the opf.flow lim option.

The variable limits (6.4) include an equality constraint on any reference bus angle

and upper and lower limits on all bus voltage magnitudes and real and reactive

generator injections:

θref

i≤θi≤θref

i, i ∈ Iref (6.12)

vi,min

m≤vi

m≤vi,max

m, i = 1 . . . nb(6.13)

pi,min

g≤pi

g≤pi,max

g, i = 1 . . . ng(6.14)

qi,min

g≤qi

g≤qi,max

g, i = 1 . . . ng.(6.15)

The voltage reference angle θref

iand voltage magnitude bounds vi,max

mand vi,min

mare

speciﬁed in columns VA (9), VMAX (12) and VMIN (13), respectively, of row iof the bus

matrix. Similarly, the generator bounds qi,max

g,qi,min

g,pi,max

gand pi,min

gare specﬁed in

columns QMAX (4), QMIN (5), PMAX (9) and PMIN (10), respectively, of row iof the gen

matrix.

6.1.1 Cartesian vs. Polar Coordinates for Voltage

Another variation of the standard AC OPF problem represents the bus voltages

in cartesian, rather than polar, coordinates. That is, instead of Θ and Vm, the

19Setting the RATE A column (6) of branch to zero is the preferred way to indicate a completely

unconstrained line.

57

optimization vector xincludes the real and imaginary parts of the complex voltage,

denoted respectively by Uand W, where V=U+jW .

x=

U

W

Pg

Qg

(6.16)

The objective function remains unchanged, but the nodal power balance con-

straints (6.7) and (6.8) and branch ﬂow constraints (6.9) and (6.10) are implemented

as functions of Uand W.

gP(U, W, Pg) = Pbus(U, W ) + Pd−CgPg= 0 (6.17)

gQ(U, W, Qg) = Qbus(U, W ) + Qd−CgQg= 0 (6.18)

hf(U, W ) = |Ff(U, W )| − Fmax ≤0 (6.19)

ht(U, W ) = |Ft(U, W )| − Fmax ≤0.(6.20)

In this formulation, the voltage angle reference constraint (6.12) and voltage mag-

nitude limits (6.13) cannot be simply applied as bounds on optimization variables.

These constrained quantities also become functions of Uand W.

θref

i≤θi(ui, wi)≤θref

i, i ∈ Iref (6.21)

vi,min

m≤vi

m(ui, wi)≤vi,max

m, i = 1 . . . nb(6.22)

In Matpower setting the opf.v cartesian option to 1 (0 by default) selects the

cartesian representation for voltages when running an AC OPF.20

6.1.2 Current vs. Power for Nodal Balance Constraints

Another variation of the standard AC OPF problem uses current balance constraints

in place of the power balance constraints (6.7)–(6.8) or (6.17)–(6.18). If we let M

and Nrepresent the real and imaginary parts, respectively, of the current, we can

express the current balance functions for the polar form as

gM(Θ, Vm, Pg, Qg) = <{Ibus(Θ, Vm)+[V∗]−1(Sd−CgSg)∗}= 0 (6.23)

gN(Θ, Vm, Pg, Qg) = ={Ibus(Θ, Vm)+[V∗]−1(Sd−CgSg)∗}= 0 (6.24)

20This option only applies to solvers based on MIPS,fmincon,Ipopt and KNITRO.

58

and for the cartesian form as

gM(U, W, Pg, Qg) = <{Ibus(U, W )+[V∗]−1(Sd−CgSg)∗}= 0 (6.25)

gN(U, W, Pg, Qg) = ={Ibus(U, W )+[V∗]−1(Sd−CgSg)∗}= 0 (6.26)

where Sd=Pd+jQd,Sg=Pg+jQgand [V∗]−1is a diagonal matrix whose i-th

diagonal entry is 1/v∗

i, that is 1

vi

mejθior 1/(ui−jwi).

In this formulation, which can be selected by setting the opf.current balance

option to 1,21 the objective function and other constraints are not aﬀected. This

option can be used in conjuntion with either the polar or cartesian representation of

bus voltages.

6.2 Standard DC OPF

When using DC network modeling assumptions and limiting polynomial costs to

second order, the standard OPF problem above can be simpliﬁed to a quadratic

program, with linear constraints and a quadratic cost function. In this case, the

voltage magnitudes and reactive powers are eliminated from the problem completely

and real power ﬂows are modeled as linear functions of the voltage angles. The

optimization variable is

x=Θ

Pg(6.27)

and the overall problem reduces to the following form.

min

Θ,Pg

ng

X

i=1

fi

P(pi

g) (6.28)

subject to

gP(Θ, Pg) = BbusΘ + Pbus,shift +Pd+Gsh −CgPg= 0 (6.29)

hf(Θ) = BfΘ + Pf,shift −Fmax ≤0 (6.30)

ht(Θ) = −BfΘ−Pf,shift −Fmax ≤0 (6.31)

θref

i≤θi≤θref

i, i ∈ Iref (6.32)

pi,min

g≤pi

g≤pi,max

g, i = 1 . . . ng(6.33)

21This option only applies to solvers based on MIPS,fmincon,Ipopt and KNITRO.

59

6.3 Extended OPF Formulation

Matpower employs an extensible OPF structure to allow the user to modify or

augment the problem formulation without rewriting the portions that are shared

with the standard OPF formulation described above. The standard formulation is

modiﬁed by introducing additional variables, user-deﬁned costs, and/or user-deﬁned

constraints. The full extended formulation can be written as follows:

min

ˆxf(x) + fu(ˆx) (6.34)

subject to

ˆg(ˆx) = 0 (6.35)

ˆ

h(ˆx)≤0 (6.36)

ˆxmin ≤ˆx≤ˆxmax (6.37)

l≤Aˆx≤u(6.38)

The ﬁrst diﬀerence to note is that the optimization variable xfrom the standard

OPF formulation has been augmented with additional variables zto form a new

optimization variable ˆx, and likewise with the lower and upper bounds.

ˆx=x

zˆxmin =xmin

zmin ˆxmax =xmax

zmax (6.39)

Second, there is an additional user-deﬁned cost term fu(ˆx) in the objective function.

This cost consists of three pieces that will be described in more detail below.

fu(ˆx) = fq(ˆx) + fnln(ˆx) + flegacy(ˆx) (6.40)

Third, the nonlinear constraints gand hare augmented with user deﬁned additions

guand huto give ˆgand ˆ

h.

ˆg(ˆx) = g(x)

gu(ˆx),ˆ

h(ˆx) = h(x)

hu(ˆx)(6.41)

And ﬁnally, a new set of linear constraints are included in (6.38).

Up through version 6.0 of Matpower, the OPF extensions were handled via op-

tional input parameters that deﬁne any additional variables,22 linear constraints,23

22Parameters zmin and zmax in (6.39).

23Parameters A,l,uin (6.38).

60

and costs of the pre-speciﬁed form deﬁned by flegacy.24 This preserved the ability to

use solvers that employ pre-compiled MEX ﬁles to compute all of the costs and con-

straints. This is referred to as Matpower’s legacy extended OPF formulation [24].

For the AC OPF, subsequent versions also include the general nonlinear con-

straints guand hu, and the quadratic and general nonlinear costs fqand fnln. The

new quadratic cost terms can be handled by all of Matpower’s AC OPF solvers,

but the general nonlinear costs and constraints require a solver that uses Matlab

code to implement the function, gradient and Hessian evaluations.25

Section 7describes the mechanisms available to the user for taking advantage of

the extensible formulation described here.

6.3.1 User-deﬁned Variables

The creation of additional user-deﬁned zvariables can be done explicitly or implicitly

based on the diﬀerence between the number of columns in Aand the dimension of

the standard OPF optimization variable x. The optional vectors zmin and zmax are

available to impose lower and upper bounds on z, respectively.

6.3.2 User-deﬁned Constraints

•Linear Constraints – The user-deﬁned linear constraints (6.38) are general

linear restrictions involving all of the optimization variables and are speciﬁed

via matrix Aand lower and upper bound vectors land u. These parameters

can be used to create equality constraints (li=ui) or inequality constraints

that are bounded below (ui=∞), bounded above (li=∞) or bounded on

both sides.

•Nonlinear Constraints – The user-deﬁned general nonlinear constraints take

the form

gj

u(x)=0 ∀j∈ Gu(6.42)

hj

u(x)≤0∀j∈ Hu,(6.43)

where Guand Huare sets of indices for user-deﬁned equality and inequality

constraint sets, respectively.

Each of these constraint sets is deﬁned by two M-ﬁle functions, similar to those

required by MIPS, one that computes the constraint values and their gradients

(Jacobian), and the other that computes Hessian values.

24Parameters H,C,N, ˆr,k,d,min (6.48)–(6.51).

25At the time of this writing, this includes solvers based on MIPS,fmincon,Ipopt and KNITRO.

61

6.3.3 User-deﬁned Costs

The user-deﬁned cost function fuconsists of three terms for three diﬀerent types of

costs: quadratic, general nonlinear, and legacy. Each term is a simple summation

over all of the cost sets of that type.

fq(ˆx) = Pjfj

q(ˆx) (6.44)

fnln(ˆx) = Pjfj

nln(ˆx) (6.45)

flegacy(ˆx) = Pjfj

legacy(ˆx) (6.46)

•Quadratic Costs for a cost set jare speciﬁed by parameters Qj,cjand kj

that deﬁne a quadratic function of the optimization variable ˆx.

fj

q(ˆx) = ˆxTQjˆx+cj

Tˆx+kj(6.47)

•General Nonlinear Costs for a cost set jconsist of a cost function fj

nln(ˆx)

provided in the form of a function handle to an M-ﬁle function that evaluates

the cost and its gradients and Hessian for a given value of ˆx.

•Legacy Costs fj

legacy(ˆx) are speciﬁed in terms of parameters Hj,Cj,Nj, ˆrj,

kj,djand mj. For simplicity of presentation, we will drop the jsubscript for

the rest of this discussion. All of the parameters are nw×1 vectors except the

symmetric nw×nwmatrix Hand the nw×nˆxmatrix N. The legacy user cost

function takes the form

flegacy(ˆx) = 1

2wTHw +CTw(6.48)

where wis deﬁned in several steps as follows. First, a new vector uis created

by applying a linear transformation Nand shift ˆrto the full set of optimization

variables

u=Nˆx−ˆr, (6.49)

then a scaled function with a “dead zone” is applied to each element of uto

produce the corresponding element of w.

wi=

mifdi(ui+ki), ui<−ki

0,−ki≤ui≤ki

mifdi(ui−ki), ui> ki

(6.50)

62

Here kispeciﬁes the size of the “dead zone”, miis a simple scale factor and fdiis

a pre-deﬁned scalar function selected by the value of di. Currently, Matpower

implements only linear and quadratic options:

fdi(α) = α, if di= 1

α2,if di= 2 (6.51)

as illustrated in Figure 6-1 and Figure 6-2, respectively.

wi

mi

ri

ˆri

ki

ki

Figure 6-1: Relationship of wito rifor di= 1 (linear option)

This form for flegacy provides the ﬂexibility to handle a wide range of costs, from

simple linear functions of the optimization variables to scaled quadratic penal-

ties on quantities, such as voltages, lying outside a desired range, to functions

of linear combinations of variables, inspired by the requirements of price coor-

dination terms found in the decomposition of large loosely coupled problems

encountered in our own research.

Some limitations are imposed on the parameters in the case of the DC OPF

since Matpower uses a generic quadratic programming (QP) solver for the

optimization. In particular, ki= 0 and di= 1 for all i, so the “dead zone” is

not considered and only the linear option is available for fdi. As a result, for

the DC case (6.50) simpliﬁes to wi=miui.

63

wi

ri

ˆri

ki

ki

Figure 6-2: Relationship of wito rifor di= 2 (quadratic option)

6.4 Standard Extensions

In addition to making this extensible OPF structure available to end users, Mat-

power also takes advantage of it internally to implement several additional capa-

bilities.

6.4.1 Piecewise Linear Costs

The standard OPF formulation in (6.1)–(6.4) does not directly handle the non-

smooth piecewise linear cost functions that typically arise from discrete bids and

oﬀers in electricity markets. When such cost functions are convex, however, they

can be modeled using a constrained cost variable (CCV) method. The piecewise lin-

ear cost function c(x) is replaced by a helper variable yand a set of linear constraints

that form a convex “basin” requiring the cost variable yto lie in the epigraph of the

function c(x).

Figure 6-3 illustrates a convex n-segment piecewise linear cost function

c(x) =

m1(x−x1) + c1, x ≤x1

m2(x−x2) + c2, x1< x ≤x2

.

.

..

.

.

mn(x−xn) + cn, xn−1< x

(6.52)

deﬁned by a sequence of points (xj, cj), j= 0 . . . n, where mjdenotes the slope of

64

x

x0

x1

x2

c

c0

c1

c2

y

cn

xn

Figure 6-3: Constrained Cost Variable

the j-th segment

mj=cj−cj−1

xj−xj−1

, j = 1 . . . n (6.53)

and x0< x1<··· < xnand m1≤m2≤ ··· < mn.

The “basin” corresponding to this cost function is formed by the following n

constraints on the helper cost variable y:

y≥mj(x−xj) + cj, j = 1 . . . n. (6.54)

The cost term added to the objective function in place of c(x) is simply the variable y.

Matpower uses this CCV approach internally to automatically generate the

appropriate helper variable, cost term and corresponding set of constraints for any

piecewise linear costs on real or reactive generation. All of Matpower’s OPF

solvers, for both AC and DC OPF problems, use the CCV approach with the ex-

ception of two that are part of the optional TSPOPF package [25], namely the

step-controlled primal/dual interior point method (SCPDIPM) and the trust region

based augmented Lagrangian method (TRALM), both of which use a cost smoothing

technique instead [26].

Note that Matpower (in opf setup) automatically converts any single-segment

piecewise linear costs into polynomial (linear) form to avoid unnecessarily creating

extra variables and constraints. It is this modiﬁed cost, rather than the original

piecewise linear equivalent, that is returned in the gencost ﬁeld of the results struct.

65

6.4.2 Dispatchable Loads

A simple approach to dispatchable or price-sensitive loads is to model them as nega-

tive real power injections with associated negative costs. This is done by specifying

a generator with a negative output, ranging from a minimum injection equal to the

negative of the largest possible load to a maximum injection of zero.

Consider the example of a price-sensitive load whose marginal beneﬁt function is

shown in Figure 6-4. The demand pdof this load will be zero for prices above λ1,p1

for prices between λ1and λ2, and p1+p2for prices below λ2.

λ(marginal beneﬁt)

$/MW

MW

λ1

λ2

p1

p2

p(load)

Figure 6-4: Marginal Beneﬁt or Bid Function

This corresponds to a negative generator with the piecewise linear cost curve

shown in Figure 6-5. Note that this approach assumes that the demand blocks can

be partially dispatched or “split”. Requiring blocks to be accepted or rejected in

their entirety would pose a mixed-integer problem that is beyond the scope of the

current Matpower implementation.

It should be noted that, with this deﬁnition of dispatchable loads as negative

generators, if the negative cost corresponds to a beneﬁt for consumption, minimizing

the cost f(x) of generation is equivalent to maximizing social welfare.

With an AC network model, there is also the question of reactive dispatch for

such loads. Typically the reactive injection for a generator is allowed to take on any

value within its deﬁned limits. Since this is not normal load behavior, the model used

66

MW

2

1

p1

p2

$

p(injection)

c(total cost)

1p1

2p2

Figure 6-5: Total Cost Function for Negative Injection

in Matpower assumes that dispatchable loads maintain a constant power factor.

When formulating the AC OPF problem, Matpower will automatically generate

an additional equality constraint to enforce a constant power factor for any “negative

generator” being used to model a dispatchable load.

The power factor, which can be lagging or leading, is determined by the ratio of

reactive to active power for the load and is speciﬁed by the active and reactive limits

deﬁning the nominal load in the gen matrix. For all dispatchable loads, the value

in the PMIN column is negative and deﬁnes the nominal active power of the load.

And PMAX is zero, allowing the load to be fully curtailed depending on the price.

The reactive limits, QMIN and QMAX, depend on whether the power ﬂow is lagging

or leading. One deﬁnes the nominal reactive load and the other must be zero to

allow the load to be fully curtailed. The values of PG and QG must be deﬁned to be

consistent with the nominal power factor.

•Lagging Power Factor – The reactive injection is negative, meaning that

reactive power is consumed by the load. Hence, QMIN is negative, QMAX is zero,

and PG and QG must be set so that QG is equal to PG * QMIN/PMIN.

•Leading Power Factor – The reactive injection is positive, that is, reactive

power is produced by the load. Hence, QMAX is positive, QMIN is zero, and PG

and QG must be set so that QG is equal to PG * QMAX/PMIN.

67

6.4.3 Generator Capability Curves

The typical AC OPF formulation includes box constraints on a generator’s real and

reactive injections, speciﬁed as simple lower and upper bounds on p(pmin and pmax)

and q(qmin and qmax). On the other hand, the true P-Qcapability curves of phys-

ical generators usually involve some tradeoﬀ between real and reactive capability,

so that it is not possible to produce the maximum real output and the maximum

(or minimum) reactive output simultaneously. To approximate this tradeoﬀ, Mat-

power includes the ability to add an upper and lower sloped portion to the standard

box constraints as illustrated in Figure 6-6, where the shaded portion represents the

feasible operating region for the unit.

qmax

1

qmax

2

qmin

2

qmin

1

qmin

qmax

pmax

pmin

p1

p2

q

p

Figure 6-6: Generator P-QCapability Curve

The two sloped portions are constructed from the lines passing through the two

pairs of points deﬁned by the six parameters p1,qmin

1,qmax

1,p2,qmin

2, and qmax

2.

If these six parameters are speciﬁed for a given generator in columns PC1–QC2MAX

(11–16), Matpower automatically constructs the corresponding additional linear

inequality constraints on pand qfor that unit.

68

If one of the sloped portions of the capability constraints is binding for genera-

tor k, the corresponding shadow price is decomposed into the corresponding µPmax

and µQmin or µQmax components and added to the respective column (MU PMAX,MU QMIN

or MU QMAX) in the kth row of gen.

6.4.4 Branch Angle Diﬀerence Limits

The diﬀerence between the bus voltage angle θfat the from end of a branch and

the angle θtat the to end can be bounded above and below to act as a proxy for a

transient stability limit, for example. If these limits are provided in columns ANGMIN

(12) and ANGMAX (13) of the branch matrix, Matpower creates the corresponding

constraints on the voltage angle variables.26

6.5 Solvers

Early versions of Matpower relied on Matlab’s Optimization Toolbox [27] to

provide the NLP and QP solvers needed to solve the AC and DC OPF problems,

respectively. While they worked reasonably well for very small systems, they did not

scale well to larger networks. Eventually, optional packages with additional solvers

were added to improve performance, typically relying on Matlab extension (MEX)

ﬁles implemented in Fortran or C and pre-compiled for each machine architecture.

Some of these MEX ﬁles are distributed as optional packages due to diﬀerences in

terms of use. For DC optimal power ﬂow, there is a MEX build [28] of the high

performance interior point BPMPD solver [29] for LP/QP problems. For the AC

OPF problem, the MINOPF [30] and TSPOPF [25] packages provide solvers suitable

for much larger systems. The former is based on MINOS [31] and the latter includes

the primal-dual interior point and trust region based augmented Lagrangian methods

described in [26]. Matpower version 4 and later also includes the option to use

the open-source Ipopt solver27 for solving both AC and DC OPFs, based on the

Matlab MEX interface to Ipopt28. It also includes the option to use CPLEX29 or

MOSEK30 for DC OPFs. Matpower 4.1 added the option to use KNITRO [32]31

26The voltage angle diﬀerence for branch kis taken to be unbounded below if branch(k, ANGMIN)

is less than or equal to −360 and unbounded above if branch(k, ANGMAX) is greater than or equal

to 360. If both parameters are zero, the voltage angle diﬀerence is unconstrained.

27Available from http://www.coin-or.org/projects/Ipopt.xml.

28See https://projects.coin-or.org/Ipopt/wiki/MatlabInterface.

29See http://www.ibm.com/software/integration/optimization/cplex-optimizer/.

30See http://www.mosek.com/.

31See http://www.ziena.com/.

69

for AC OPFs and the Gurobi Optimizer [33]32 for DC OPFs and Matpower 5

added GLPK [34] and 5.1 added CLP [35]. See Appendix Gfor more details on

these optional packages.

Beginnning with version 4, Matpower also includes its own primal-dual interior

point method implemented in pure-Matlab code, derived from the MEX implemen-

tation of the algorithms described in [26]. This solver is called MIPS (Matpower

Interior Point Solver) and is described in more detail in Appendix A. If no optional

packages are installed, MIPS will be used by default for both the AC OPF and

as the QP solver used by the DC OPF. The AC OPF solver also employs a unique

technique for eﬃciently forming the required Hessians via a few simple matrix opera-

tions [36–38]. This solver has application to general nonlinear optimization problems

outside of Matpower and can be called directly as mips. There is also a conve-

nience wrapper function called qps mips making it trivial to set up and solve LP and

QP problems, with an interface similar to quadprog from the Matlab Optimization

Toolbox.

6.6 runopf

In Matpower, an optimal power ﬂow is executed by calling runopf with a case

struct or case ﬁle name as the ﬁrst argument (casedata). In addition to printing

output to the screen, which it does by default, runopf optionally returns the solution

in a results struct.

>> results = runopf(casedata);

The results struct is a superset of the input Matpower case struct mpc, with

some additional ﬁelds as well as additional columns in some of the existing data

ﬁelds. In addition to the solution values included in the results for a simple power

ﬂow, shown in Table 4-1 in Section 4.4, the following additional optimal power ﬂow

solution values are stored as shown in Table 6-1.

32See http://www.gurobi.com/.

70

Table 6-1: Optimal Power Flow Results

name description

results.f ﬁnal objective function value

results.x ﬁnal value of optimization variables (internal order)

results.om OPF model object†

results.bus(:, LAM P) Lagrange multiplier on real power mismatch

results.bus(:, LAM Q) Lagrange multiplier on reactive power mismatch

results.bus(:, MU VMAX) Kuhn-Tucker multiplier on upper voltage limit

results.bus(:, MU VMIN) Kuhn-Tucker multiplier on lower voltage limit

results.gen(:, MU PMAX) Kuhn-Tucker multiplier on upper Pglimit

results.gen(:, MU PMIN) Kuhn-Tucker multiplier on lower Pglimit

results.gen(:, MU QMAX) Kuhn-Tucker multiplier on upper Qglimit

results.gen(:, MU QMIN) Kuhn-Tucker multiplier on lower Qglimit

results.branch(:, MU SF) Kuhn-Tucker multiplier on ﬂow limit at “from” bus

results.branch(:, MU ST) Kuhn-Tucker multiplier on ﬂow limit at “to” bus

results.mu shadow prices of constraints‡

results.g (optional) constraint values

results.dg (optional) constraint 1st derivatives

results.raw raw solver output in form returned by MINOS, and more‡

results.var.val ﬁnal value of optimization variables, by named subset‡

results.var.mu shadow prices on variable bounds, by named subset‡

results.nle shadow prices on nonlinear equality constraints, by named subset‡

results.nli shadow prices on nonlinear inequality constraints, by named subset‡

results.lin shadow prices on linear constraints, by named subset‡

results.cost ﬁnal value of user-deﬁned costs, by named subset‡

†See help for opf model and opt model for more details.

‡See help for opf for more details.

71

Additional optional input arguments can be used to set options (mpopt) and

provide ﬁle names for saving the pretty printed output (fname) or the solved case

data (solvedcase).

>> results = runopf(casedata, mpopt, fname, solvedcase);

Some of the main options that control the optimal power ﬂow simulation are listed

in Tables 6-2 and 6-3. There are many other options that can be used to control the

termination criteria and other behavior of the individual solvers. See Appendix C

or the mpoption help for details. As with runpf the output printed to the screen can

be controlled by the options in Table 4-3, but there are additional output options

for the OPF, related to the display of binding constraints that are listed Table 6-4,

along with an option that can be used to force the AC OPF to return information

about the constraint values and Jacobian and the objective function gradient and

Hessian.

For OPF problems, the preferred way to eliminate the ﬂow limit on a branch is

to set the RATE A column (6) of the corresponding row in the branch matrix to zero.

This indicates a completely unconstrained ﬂow (as opposed to zero ﬂow).

By default, runopf solves an AC optimal power ﬂow problem using a primal dual

interior point method. To run a DC OPF, the model option must be set to 'DC'. For

convenience, Matpower provides a function rundcopf which is simply a wrapper

that sets the model to 'DC'before calling runopf.

Internally, the runopf function does a number of conversions to the problem

data before calling the appropriate solver routine for the selected OPF algorithm.

This external-to-internal format conversion is performed by the ext2int function,

described in more detail in Section 7.3.1, and includes the elimination of out-of-service

equipment, the consecutive renumbering of buses and the reordering of generators

by increasing bus number. All computations are done using this internal indexing.

When the simulation has completed, the data is converted back to external format

by int2ext before the results are printed and returned. In addition, both ext2int

and int2ext can be customized via user-supplied callback routines to convert data

needed by user-supplied variables, constraints or costs into internal indexing.

72

Table 6-2: Optimal Power Flow Solver Options

name default description

opf.ac.solver 'DEFAULT'AC optimal power ﬂow solver:

'DEFAULT'– choose default solver, i.e. 'MIPS'

'MIPS'– MIPS, Matpower Interior Point Solver,

primal/dual interior point method†

'FMINCON'–Matlab Optimization Toolbox, fmincon

'IPOPT'–Ipopt*

'KNITRO'– KNITRO*

'MINOPF'– MINOPF*

, MINOS-based solver

'PDIPM'– PDIPM*

, primal/dual interior point method‡

'SDPOPF'– SDPOPF*

, solver based on semideﬁnite relax-

ation

'TRALM'– TRALM*

, trust region based augmented Lan-

grangian method

opf.dc.solver 'DEFAULT'DC optimal power ﬂow solver:

'DEFAULT'– choose default solver based on availability

in the following order: 'GUROBI','CPLEX',

'MOSEK','OT','GLPK'(linear costs only),

'BPMPD','MIPS'

'MIPS'– MIPS, Matpower Interior Point Solver,

primal/dual interior point method†

'BPMPD'– BPMPD*

'CLP'– CLP*

'CPLEX'– CPLEX*

'GLPK'– GLPK*(no quadratic costs)

'GUROBI'– Gurobi*

'IPOPT'–Ipopt*

'MOSEK'– MOSEK*

'OT'–Matlab Opt Toolbox, quadprog,linprog

*Requires the installation of an optional package. See Appendix Gfor details on the corresponding package.

†For MIPS-sc, the step-controlled version of this solver, the mips.step control option must be set to 1.

‡For SC-PDIPM, the step-controlled version of this solver, the pdipm.step control option must be set to 1.

73

Table 6-3: Other OPF Options

name default description

opf.current balance 0 use current, as opposed to power, balance formulation for

AC OPF, 0 or 1

opf.v cartesian 0 use cartesian, as opposed to polar, representation for volt-

ages for AC OPF, 0 or 1

opf.violation 5×10−6constraint violation tolerance

opf.use vg 0 respect generator voltage setpoint, 0 or 1†

0 – use voltage magnitude limits speciﬁed in bus, ig-

nore VG in gen

1 – replace voltage magnitude limits speciﬁed in bus

by VG in corresponding gen

opf.flow lim 'S'quantity to limit for branch ﬂow constraints

'S'– apparent power ﬂow (limit in MVA)

'P'– active power ﬂow (limit in MW)

'I'– current magnitude (limit in MVA at 1 p.u.

voltage)

'2'– same as 'P', but implemented using square

of active ﬂow, rather than simple max

opf.ignore angle lim 0 ignore angle diﬀerence limits for branches

0 – include angle diﬀerence limits, if speciﬁed

1 – ignore angle diﬀerence limits even if speciﬁed

opf.softlims.default 1 behavior of OPF soft limits for which parameters are not

explicitly provided

0 – do not include softlims if not explicitly speciﬁed

1 – include softlims with default values if not ex-

plicitly speciﬁed

opf.init from mpc*-1 specify whether to use the current state in Matpower

case to initialize OPF

-1 – Matpower decides based on solver/algorithm

0 – ignore current state in Matpower case‡

1 – use current state in Matpower case

opf.start 0 strategy for initializing OPF starting point

0 – default, Matpower decides based on solver,

(currently identical to 1)

1 – ignore current state in Matpower case‡

2 – use current state in Matpower case

3 – solve power ﬂow and use resulting state

opf.return raw der 0 for AC OPF, return constraint and derivative info in

results.raw (in ﬁelds g,dg,df,d2f)

*Deprecated. Use opf.start instead.

†Using a value between 0 and 1 results in the limits being determined by the corresponding weighted average of

the 2 options.

‡Only applies to fmincon,Ipopt, KNITRO and MIPS solvers, which use an interior point estimate; others use

current state in Matpower case, as with opf.start = 2.

74

Table 6-4: OPF Output Options

name default description

out.lim.all -1 controls constraint info output

-1 – individual ﬂags control what is printed

0 – do not print any constraint info†

1 – print only binding constraint info†

2 – print all constraint info†

out.lim.v 1 control output of voltage limit info

0 – do not print

1 – print binding constraints only

2 – print all constraints

out.lim.line 1 control output of line ﬂow limit info‡

out.lim.pg 1 control output of gen active power limit info‡

out.lim.qg 1 control output of gen reactive power limit info‡

†Overrides individual ﬂags.

‡Takes values of 0, 1 or 2 as for out.lim.v.

75

7 Extending the OPF

The extended OPF formulation described in Section 6.3 allows the user to modify the

standard OPF formulation to include additional variables, costs and/or constraints.

There are two primary mechanisms available for the user to accomplish this. The ﬁrst

is by directly constructing the full parameters for the addional costs or constraints

and supplying them either as ﬁelds in the case struct or directly as arguments to

the opf function. The second, and more powerful, method is via a set of callback

functions that customize the OPF at various stages of the execution. Matpower

includes several examples of using the latter method, for example to add a ﬁxed zonal

reserve requirement, to implement interface ﬂow limits, dispatchable DC transmission

lines, or branch ﬂow soft limits.

7.1 Direct Speciﬁcation

This section describes the additional ﬁelds that can be added to mpc (the Matpower

case struct) or, alternatively, the input parameters that can be passed directly to the

opf function. The former is the preferred approach, as it allows the normal use of

runopf and related functions.

7.1.1 User-deﬁned Variables

In the case of direct speciﬁcation, additional variables are created implicitly based on

the diﬀerence between the number of columns in Aand the number nxof standard

OPF variables. If Ahas more columns than xhas elements, the extra columns are

assumed to correspond to a new zvariable. The initial value and lower and upper

bounds for zcan also be speciﬁed in the optional ﬁelds or arguments, z0,zl and zu,

respectively.

7.1.2 User-deﬁned Constraints

•Linear Constraints can be added by specifying the A,land uparameters of

(6.38) directly as ﬁelds or arguments of the same names, A,land u, respectively,

where Ais sparse.

•Nonlinear Constraints for a constraint set j, that is, gj

u(x) or hj

u(x) from

(6.42) or (6.43), are implemented by deﬁning two M-ﬁle functions, similar to

76

those required by MIPS. The ﬁrst is a function to compute the constraint values

and their gradients (Jacobian),33 with calling syntax:

[g, dg] = my constraint fcn(x, <p1>, <p2>, ...)

Here <p1> and <p2> represent arbitrary optional input parameters that remain

constant for all calls to the function. The second is a function to compute the

Hessian of the corresponding term in the Lagrangian function, that is the term

corresponding to λTgj

u(x) or µThj

u(x). This function has the following calling

syntax:

d2G = my constraint hess(x, lambda, <p1>, <p2>, ...)

Once again this is similar to the form of the Hessian evaluation function ex-

pected by MIPS.34

In order to add a set of user-deﬁned nonlinear constraints to the formula-

tion through direct speciﬁcation, an entry must be added to the optional

user constraints.nle (nonlinear equality) or user constraints.nli (nonlin-

ear inequality) ﬁelds of the Matpower case struct. These ﬁelds are cell ar-

rays, in which each element (also a cell array) speciﬁes a constraint set of the

corresponding kind. The format is as follows, where the details are described

in Table 7-1.

mpc.user_constraints.nle = {

{name, N, g_fcn, hess_fcn, varsets, params},

...

};

There is an example in the OPF testing code (e.g. t opf mips) near the end,

that implements a nonlinear relationship between three diﬀerent generator out-

puts by adding the following line.

33The diﬀerences between this function and the one required by MIPS is that (1) this one repre-

sents a single set of equality constraints or inequality constraints, not both, (2) it returns in dg the

m×nJacobian (for mconstraints and nvariables), whereas MIPS expects the n×mtranspose,

(3) this function allows arbitrary additional parameters, and (4) xcan be a cell array of sub-vectors

of optimization variable ˆx(see varsets in Table 7-1).

34The diﬀerences between this function and the one required by MIPS is that (1) this one returns

the Hessian of a single term of the Lagrangian, not of the full Lagrangian, (2) the lambda argument

is a simple vector of multipliers (λor µ) corresponding just to this set of constraints, not a struct

containing the full λand µvectors, in ﬁelds eqnonlin and ineqnonlin, respectively, and (3) this

function allows arbitrary additional parameters, and (4) xcan be a cell array of sub-vectors of

optimization variable ˆx(see varsets in Table 7-1).

77

Table 7-1: User-deﬁned Nonlinear Constraint Speciﬁcation

name description

name string with name of constraint set, used to label multipliers in results struct

Nnumber of constraints, i.e. dimension of gj

u(x) (or hj

u(x) as the case may be)

g fcn string containing name of function to evaluate constraint and gradients (Jacobian)

hess fcn string containing name of function to evaluate Hessian of corresponding term of

Lagrangian

varsets cell array of variable set names, specifying sub-vectors of optimization vector ˆxto

be passed as inputs in x†

params cell array of optional, arbitrary parameters to pass to each call to the constraint and

Hessian evaluation functions

†If varsets is empty, xwill be the full optimization vector ˆx, otherwise it will be a cell array of sub-vectors of ˆxfor

the speciﬁed variable sets. Valid names include 'Va','Vm','Pg', and 'Qg'. It can include others depending on

the OPF extensions in use. See the variable names displayed by results.om for a complete list for your problem.

mpc.user_constraints.nle = {

{'Pg_usr', 1, 'opf_nle_fcn1', 'opf_nle_hess1', {'Pg'}, {}}

};

This adds a single constraint named 'Pg usr'as a function of the vector of

generator active injections (xwill be {Pg}, where Pg is a sub-vector of the

optimization vector ˆxcontaining only the generator active injections). The

constraints and gradients are evaluated by a function named 'opf nle fcn1'

and the Hessian of the corresponding term of the Lagrangian by a function

named 'opf nle hess1', neither of which expect any additional parameters.

7.1.3 User-deﬁned Costs

•Quadratic Costs by direct speciﬁcation are not currently supported.

•General Nonlinear Costs by direct speciﬁcation are not currently supported.

•Legacy Costs – To add legacy costs directly, the parameters H,C,N, ˆr,k,

dand mof (6.48)–(6.51) described in Section 6.3.3 are speciﬁed as ﬁelds or

arguments H,Cw,Nand fparm, respectively, where fparm is the nw×4 matrix

fparm =dˆr k m .(7.1)

When specifying additional costs, Nand Cw are required, while Hand fparm are

optional. The default value for His a zero matrix, and the default for fparm is

78

such that dand mare all ones and ˆrand kare all zeros, resulting in simple

linear cost, with no shift or “dead-zone”. Nand Hshould be speciﬁed as sparse

matrices.

7.1.4 Additional Comments

For a simple formulation extension to be used for a small number of OPF cases, this

method has the advantage of being direct and straightforward. While Matpower

does include code to eliminate the columns of Aand Ncorresponding to Vmand Qg

when running a DC OPF35, as well as code to reorder and eliminate columns appro-

priately when converting from external to internal data formats, this mechanism still

requires the user to take special care in preparing the Aand Nmatrices to ensure

that the columns match the ordering of the elements of the optimization vectors x

and z. All extra constraints and variables must be incorporated into a single set of

parameters that are constructed before calling the OPF. The bookkeeping needed to

access the resulting variables and shadow prices on constraints and variable bounds

must be handled manually by the user outside of the OPF, along with any pro-

cessing of additional input data and processing, printing or saving of the additional

result data. Making further modiﬁcations to a formulation that already includes

user-supplied costs, constraints or variables, requires that both sets be incorporated

into a new single consistent set of parameters.

7.2 Callback Functions

The second method, based on deﬁning a set of callback functions, oﬀers several

distinct advantages, especially for more complex scenarios or for adding a feature

for others to use, such as the zonal reserve requirement or the interface ﬂow limits

mentioned previously. This approach makes it possible to:

•deﬁne and access variable/constraint/cost sets as individual named blocks

•deﬁne constraints, costs only in terms of variables directly involved

•pre-process input data and/or post-process result data

•print and save new result data

•simultaneously use multiple, independently developed extensions (e.g. zonal

reserve requirements and interface ﬂow limits)

35Only if they contain all zeros.

79

With this approach the OPF formulation is modiﬁed in the formulation callback,

which is described and illustrated below in Section 7.3.2. The modiﬁcations to the

formulation are handled by adding variables, costs and/or constraints to the OPF

Model object (om) using one of the add * methods. Please see the documentation for

each of these methods for more details.

7.2.1 User-deﬁned Variables

Additional variables are added using the add var method.

om.add_var('myV', N, myV0, myV_min, myV_max);

Here myV is the name of the variable set, Nthe number of variables being added (i.e.

this is an N×1 vector being appended to the current optimization variable ˆx), and

myV0,myV min,myV max are vectors of initial values, lower bounds and upper bounds,

respectively.

7.2.2 User-deﬁned Costs

•Quadratic Costs are added using the add quad cost method.

om.add_quad_cost('myQcost', Q, c, k, varsets);

Here myQcost is the name of the cost set, Q,cand kare the Qj,cjand kj

parameters from (6.47), respectively, and var sets is an optional cell array of

variable set names. If var sets is provided, the dimensions of Qand cwill

correspond, not to the full optimization vector ˆx, but to the subvector formed

by stacking only the variable sets speciﬁed in var sets.

•General Nonlinear Costs are added using the add nln cost method.

om.add_nln_cost('myNonlinCost', 1, fcn, varsets);

Here myNonlinCost is the name of the cost set, fcn is a function handle and

var sets is an optional cell array of variable set names. The function handle

fcn must point to a function that implements the fj

nln(ˆx) function from (6.46),

returning the value of the function, along with its gradient and Hessian. This

function has a format similar to that required by MIPS for its cost function.

If var sets is provided, instead of the full optimization vector ˆx, the xpassed

to the function fcn will be a cell array of sub-vectors of ˆxcorresponding to the

speciﬁed variable sets.

80

•Legacy Costs are added using the add legacy cost method.

om.add_legacy_cost('myLegacyCost', cp, varsets);

Here myLegacyCost is the name of the cost set, cp is a struct containing the cost

parameters H,C,N, ˆr,k,dand mfrom (6.48)–(6.51) in ﬁelds H,Cw,N,rh,kk,

dd and mm, respectively, and var sets is an optional cell array of variable set

names. If var sets is provided, the number of columns in Nwill correspond,

not to the full optimization vector ˆx, but to the subvector formed by stacking

only the variable sets speciﬁed in var sets.

7.2.3 User-deﬁned Constraints

•Linear Constraints are added using the add lin constraint method.

om.lin_constraint('myLinCons', A, l, u, varsets);

Here myLinCons is the name of the constraint set, A,land ucorrespond to

additional rows to add to the A,land uparameters in (6.38), respectively, and

var sets is an optional cell array of variable set names. If var sets is provided,

the number of columns in Awill correspond, not to the full optimization vector

ˆx, but to the subvector formed by stacking only the variable sets speciﬁed in

var sets.

•Nonlinear Constraints are added using the add nln constraint method.

om.nln_constraint('myNonlinCons', N, iseq, fcn, hess, varsets);

Here myNonlinCons is the name of the constraint set, Nis the number of con-

straints being added, iseq is a true or false value indicating whether this is a set

of equality (or else inequality) constraints, fcn and hess are function handles

and var sets is an optional cell array of variable set names. The function han-

dle fcn must point to a function that implements the corresponding constraint

function gj

u(x) or hj

u(x) and its gradients, and hess to one that evaluates the

corresponding term in the Hessian of the Lagrangian function. This function

has a format similar to that required by MIPS for its cost function. See more

details in Section 7.1.2.

If var sets is provided, instead of the full optimization vector ˆx, the xpassed

to the function fcn will be a cell array of sub-vectors of ˆxcorresponding to the

speciﬁed variable sets.

81

7.3 Callback Stages and Example

Matpower deﬁnes ﬁve stages in the execution of a simulation where custom code

can be inserted to alter the behavior or data before proceeding to the next stage.

This custom code is deﬁned as a set of “callback” functions that are registered via

add userfcn for Matpower to call automatically at one of the ﬁve stages. Each

stage has a name and, by convention, the name of a user-deﬁned callback func-

tion ends with the name of the corresponding stage. For example, a callback for

the formulation stage that modiﬁes the OPF problem formulation to add reserve

requirements could be registered with the following line of code.

mpc = add_userfcn(mpc, 'formulation', @userfcn_reserves_formulation);

The sections below will describe each stage and the input and output arguments

for the corresponding callback function, which vary depending on the stage. An

example that employs additional variables, constraints and costs will be used for

illustration.

Consider the problem of jointly optimizing the allocation of both energy and

reserves, where the reserve requirements are deﬁned as a set of nrz ﬁxed zonal MW

quantities. Let Zkbe the set of generators in zone kand Rkbe the MW reserve

requirement for zone k. A new set of variables rare introduced representing the

reserves provided by each generator. The value ri, for generator i, must be non-

negative and is limited above by a user-provided upper bound rmax

i(e.g. a reserve

oﬀer quantity) as well as the physical ramp rate ∆i.

0≤ri≤min(rmax

i,∆i), i = 1 . . . ng(7.2)

If the vector ccontains the marginal cost of reserves for each generator, the user

deﬁned cost term from (6.34) is simply

fu(ˆx) = cTr. (7.3)

There are two additional sets of constraints needed. The ﬁrst ensures that, for

each generator, the total amount of energy plus reserve provided does not exceed the

capacity of the unit.

pi

g+ri≤pi,max

g, i = 1 . . . ng(7.4)

The second requires that the sum of the reserve allocated within each zone kmeets

the stated requirements.

X

i∈Zk

ri≥Rk, k = 1 . . . nrz (7.5)

82

Table 7-2: Names Used by Implementation of OPF with Reserves

name description

mpc Matpower case struct

reserves additional ﬁeld in mpc containing input parameters for zonal reserves in

the following sub-ﬁelds:

cost ng×1 vector of reserve costs, cfrom (7.3)

qty ng×1 vector of reserve quantity upper bounds, ith element is rmax

i

zones nrz ×ngmatrix of reserve zone deﬁnitions

zones(k,j) =1 if gen jbelongs to reserve zone k(j∈Zk)

0 otherwise (j /∈Zk)

req nrz ×1 vector of zonal reserve requirements, kth element is Rkfrom (7.5)

om OPF model object, already includes standard OPF setup

results OPF results struct, superset of mpc with additional ﬁelds for output data

ng ng, number of generators

Rname for new reserve variable block, ith element is ri

Pg plus R name for new capacity limit constraint set (7.4)

Rreq name for new reserve requirement constraint set (7.5)

Table 7-2 describes some of the variables and names that are used in the example

callback function listings in the sections below.

7.3.1 ext2int Callback

Before doing any simulation of a case, Matpower performs some data conversion

on the case struct in order to achieve a consistent internal structure, by calling the

following.

mpc = ext2int(mpc, mpopt);

All isolated buses, out-of-service generators and branches are removed, along with

any generators or branches connected to isolated buses. The buses are renumbered

consecutively, beginning at 1, and the in-service generators are sorted by increasing

bus number. All of the related indexing information and the original data matrices

are stored in an order ﬁeld in the case struct to be used later by int2ext to perform

the reverse conversions when the simulation is complete.

The ﬁrst stage callback is invoked from within the ext2int function immediately

after the case data has been converted. Inputs are a Matpower case struct (mpc)

freshly converted to internal indexing, a Matpower options struct mpopt,36 and

36The mpopt may be empty if ext2int() is called manually without providing it.

83

any (optional) args value supplied when the callback was registered via add userfcn.

Output is the (presumably updated) mpc. This is typically used to reorder any input

arguments that may be needed in internal ordering by the formulation stage. The

example shows how e2i field can also be used, with a case struct that has already

been converted to internal indexing, to convert other data structures by passing in

2 or 3 extra parameters in addition to the case struct. In this case, it automatically

converts the input data in the qty,cost and zones ﬁelds of mpc.reserves to be

consistent with the internal generator ordering, where oﬀ-line generators have been

eliminated and the on-line generators are sorted in order of increasing bus number.

Notice that it is the second dimension (columns) of mpc.reserves.zones that is being

re-ordered. See the on-line help for e2i field and e2i data for more details on what

all they can do.

function mpc = userfcn_reserves_ext2int(mpc, mpopt, args)

mpc = e2i_field(mpc, {'reserves', 'qty'}, 'gen');

mpc = e2i_field(mpc, {'reserves', 'cost'}, 'gen');

mpc = e2i_field(mpc, {'reserves', 'zones'}, 'gen', 2);

This stage is also a good place to check the consistency of any additional input

data required by the extension and throw an error if something is missing or not as

expected.

7.3.2 formulation Callback

This stage is called at the end of opf setup after the OPF Model (om) object has been

initialized with the standard OPF formulation, but before calling the solver. This

is the ideal place for modifying the problem formulation with additional variables,

constraints and costs, using the add var,add lin constraint,add nln constraint,

add quad cost,add nln cost and add legacy cost methods of the OPF Model ob-

ject.37 Inputs are the om object, the Matpower options struct mpopt and any

(optional) args supplied when the callback was registered via add userfcn. Output

is the updated om object.

The om object contains both the original Matpower case data as well as all of

the indexing data for the variables and constraints of the standard OPF formulation.

37It is perfectly legitimate to register more than one callback per stage, such as when enabling

multiple independent OPF extensions. In this case, the callbacks are executed in the order they

were registered with add userfcn. E.g. when the second and subsequent formulation callbacks

are invoked, the om object will reﬂect any modiﬁcations performed by earlier formulation callbacks.

84

See the on-line help for opf model and opt model for more details on the OPF model

object and the methods available for manipulating and accessing it.

In the example code, a new variable block named Rwith ngelements and the

limits from (7.2) is added to the model via the add var method. Similarly, two linear

constraint blocks named Pg plus R and Rreq, implementing (7.4) and (7.5), respec-

tively, are added via the add lin constraint method. And ﬁnally, the add quad cost

method is used to add to the model a quadratic (actually linear) cost block corre-

sponding to (7.3).

Notice that the last argument to add lin constraint and add quad cost allows

the constraints and costs to be deﬁned only in terms of the relevant parts of the

optimization variable ˆx. For example, the Amatrix for the Pg plus R constraint

contains only columns corresponding to real power generation (Pg) and reserves (R)

and need not bother with voltages, reactive power injections, etc. As illustrated in

Figure 7-1, this allows the same code to be used with both the AC OPF, where

ˆxincludes Vmand Qg, and the DC OPF where it does not. This code is also

independent of any additional variables that may have been added by Matpower

(e.g. yvariables from Matpower’s CCV handling of piece-wise linear costs) or by

the user via previous formulation callbacks. Matpower will place the constraint

and cost matrix blocks in the appropriate place when it constructs the aggregated

constraint and cost matrices at run-time. This is an important feature that enables

independently developed Matpower OPF extensions to work together.

85

function om = userfcn_reserves_formulation(om, mpopt, args)

%% initialize some things

define_constants;

mpc = om.get_mpc();

r = mpc.reserves;

ng = size(mpc.gen, 1); %% number of on-line gens

%% variable bounds

Rmin = zeros(ng, 1); %% bound below by 0

Rmax = r.qty; %% bound above by stated max reserve qty ...

k = find(mpc.gen(:, RAMP_10) > 0 & mpc.gen(:, RAMP_10) < Rmax);

Rmax(k) = mpc.gen(k, RAMP_10); %% ... and ramp rate

Rmax = Rmax / mpc.baseMVA;

%% constraints

I = speye(ng); %% identity matrix

Ar = [I I];

Pmax = mpc.gen(:, PMAX) / mpc.baseMVA;

lreq = r.req / mpc.baseMVA;

%% cost

Cw = r.cost * mpc.baseMVA; %% per unit cost coefficients

%% add them to the model

om.add_var('R', ng, [], Rmin, Rmax);

om.add_lin_constraint('Pg_plus_R', Ar, [], Pmax, {'Pg', 'R'});

om.add_lin_constraint('Rreq', r.zones, lreq, [], {'R'});

om.add_quad_cost('Rcost', [], Cw, 0, {'R'});

86

Va

A = 0A10 00 A2

A = 0A1 0A2

Ar = A1 A2

Va PgVm yQg R

Pg y R

AC OPF

DC OPF

Figure 7-1: Adding Constraints Across Subsets of Variables

87

7.3.3 int2ext Callback

After the simulation is complete and before the results are printed or saved, Mat-

power converts the case data in the results struct back to external indexing by

calling the following.

results = int2ext(results, mpopt);

This conversion essentially undoes everything that was done by ext2int. Generators

are restored to their original ordering, buses to their original numbering and all

out-of-service or isolated generators, branches and buses are restored.

This callback is invoked from int2ext immediately before the resulting case is

converted from internal back to external indexing. At this point, the simulation has

been completed and the results struct, a superset of the original Matpower case

struct passed to the OPF, contains all of the results. This results struct is passed to

the callback, along with the Matpower options struct mpopt,38 and any (optional)

args supplied when the callback was registered via add userfcn. The output of the

callback is the updated results struct. This is typically used to convert any results

to external indexing and populate any corresponding ﬁelds in the results struct.

The results struct contains, in addition to the standard OPF results, solution

information related to all of the user-deﬁned variables, constraints and costs. Ta-

ble 7-3 summarizes where the various data is found. Each of the ﬁelds listed in

the table is actually a struct whose ﬁelds correspond to the named sets created by

add var,add lin constraint,add nln constraint,add quad cost,add nln cost and

add legacy cost.

In the example code below, the callback function begins by converting the reserves

input data in the resulting case (qty,cost and zones ﬁelds of results.reserves) back

to external indexing via calls to i2e field. See the help for i2e field and i2e data

for more details on how they can be used.

38The mpopt may be empty if int2ext() is called manually without providing it.

88

Table 7-3: Results for User-Deﬁned Variables, Constraints and Costs

name description

results.var.val ﬁnal value of user-deﬁned variables

results.var.mu.l shadow price on lower limit of user-deﬁned variables

results.var.mu.u shadow price on upper limit of user-deﬁned variables

results.lin.mu.l shadow price on lower (left-hand) limit of linear constraints

results.lin.mu.u shadow price on upper (right-hand) limit of linear constraints

results.nle.lambda shadow price on nonlinear equality constraints

results.nli.mu shadow price on nonlinear inequality constraints

results.cost ﬁnal value of legacy user costs

results.nlc ﬁnal value of general nonlinear costs

results.qdc ﬁnal value of quadratic costs

89

Then the reserves results of interest are extracted from the appropriate sub-ﬁelds

of results.var,results.lin and results.cost, converted from per unit to per MW

where necessary, and stored with external indexing for the end user in the chosen

ﬁelds of the results struct.

function results = userfcn_reserves_int2ext(results, mpopt, args)

%%----- convert stuff back to external indexing -----

%% convert all reserve parameters (zones, costs, qty, rgens)

results = i2e_field(results, {'reserves', 'qty'}, 'gen');

results = i2e_field(results, {'reserves', 'cost'}, 'gen');

results = i2e_field(results, {'reserves', 'zones'}, 'gen', 2);

r = results.reserves;

ng = size(results.gen, 1); %% number of on-line gens (internal)

ng0 = size(results.order.ext.gen, 1); %% number of gens (external)

%%----- results post-processing -----

%% get the results (per gen reserves, multipliers) with internal gen indexing

%% and convert from p.u. to per MW units

[R0, Rl, Ru] = results.om.params_var('R');

R = results.var.val.R * results.baseMVA;

Rmin = Rl * results.baseMVA;

Rmax = Ru * results.baseMVA;

mu_l = results.var.mu.l.R / results.baseMVA;

mu_u = results.var.mu.u.R / results.baseMVA;

mu_Pmax = results.lin.mu.u.Pg_plus_R / results.baseMVA;

%% store in results in results struct

z = zeros(ng0, 1);

results.reserves.R = i2e_data(results, R, z, 'gen');

results.reserves.Rmin = i2e_data(results, Rmin, z, 'gen');

results.reserves.Rmax = i2e_data(results, Rmax, z, 'gen');

results.reserves.mu.l = i2e_data(results, mu_l, z, 'gen');

results.reserves.mu.u = i2e_data(results, mu_u, z, 'gen');

results.reserves.mu.Pmax = i2e_data(results, mu_Pmax, z, 'gen');

results.reserves.prc = z;

for k = 1:ng0

iz = find(r.zones(:, k));

results.reserves.prc(k) = sum(results.lin.mu.l.Rreq(iz)) / results.baseMVA;

end

results.reserves.totalcost = results.cost.Rcost;

90

7.3.4 printpf Callback

The pretty-printing of the standard OPF output is done via a call to printpf after

the case has been converted back to external indexing. This callback is invoked from

within printpf after the pretty-printing of the standard OPF output. Inputs are

the results struct, the ﬁle descriptor to write to, a Matpower options struct, and

any (optional) args supplied via add userfcn. Output is the results struct. This is

typically used for any additional pretty-printing of results.

In this example, the out.all ﬂag in the options struct is checked before printing

anything. If it is non-zero, the reserve quantities and prices for each unit are printed

ﬁrst, followed by the per-zone summaries. An additional table with reserve limit

shadow prices might also be included.

91

function results = userfcn_reserves_printpf(results, fd, mpopt, args)

%% define named indices into data matrices

[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...

MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...

QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;

%%----- print results -----

r = results.reserves;

ng = length(r.R);

nrz = size(r.req, 1);

if mpopt.out.all ~= 0

fprintf(fd, '\n=======================================================');

fprintf(fd, '\n| Reserves |');

fprintf(fd, '\n=======================================================');

fprintf(fd, '\n Gen Bus Status Reserves Price');

fprintf(fd, '\n # # (MW) ($/MW)');

fprintf(fd, '\n---- ----- ------ -------- --------');

for k = 1:ng

fprintf(fd, '\n%3d %6d %2d ', k, results.gen(k, GEN_BUS), ...

results.gen(k, GEN_STATUS));

if results.gen(k, GEN_STATUS) > 0 && abs(results.reserves.R(k)) > 1e-6

fprintf(fd, '%10.2f', results.reserves.R(k));

else

fprintf(fd, ' - ');

end

fprintf(fd, '%10.2f ', results.reserves.prc(k));

end

fprintf(fd, '\n --------');

fprintf(fd, '\n Total:%10.2f Total Cost: $%.2f', ...

sum(results.reserves.R(r.igr)), results.reserves.totalcost);

fprintf(fd, '\n');

fprintf(fd, '\nZone Reserves Price ');

fprintf(fd, '\n # (MW) ($/MW) ');

fprintf(fd, '\n---- -------- --------');

for k = 1:nrz

iz = find(r.zones(k, :)); %% gens in zone k

fprintf(fd, '\n%3d%10.2f%10.2f', k, sum(results.reserves.R(iz)), ...

results.lin.mu.l.Rreq(k) / results.baseMVA);

end

fprintf(fd, '\n');

%% print binding reserve limit multipliers ...

end

92

7.3.5 savecase Callback

The savecase is used to save a Matpower case struct to an M-ﬁle, for example,

to save the results of an OPF run. The savecase callback is invoked from savecase

after printing all of the other data to the ﬁle. Inputs are the case struct, the ﬁle

descriptor to write to, the variable preﬁx (typically 'mpc.') and any (optional) args

supplied via add userfcn. Output is the case struct. The purpose of this callback is

to write any non-standard case struct ﬁelds to the case ﬁle.

In this example, the zones,req,cost and qty ﬁelds of mpc.reserves are written

to the M-ﬁle. This ensures that a case with reserve data, if it is loaded via loadcase,

possibly run, then saved via savecase, will not lose the data in the reserves ﬁeld.

This callback could also include the saving of the output ﬁelds if present. The

contributed serialize function39 can be very useful for this purpose.

39http://www.mathworks.com/matlabcentral/fileexchange/1206

93

function mpc = userfcn_reserves_savecase(mpc, fd, prefix, args)

%

% mpc = userfcn_reserves_savecase(mpc, fd, mpopt, args)

%

% This is the 'savecase' stage userfcn callback that prints the M-file

% code to save the 'reserves' field in the case file. It expects a

% MATPOWER case struct (mpc), a file descriptor and variable prefix

% (usually 'mpc.'). The optional args are not currently used.

r = mpc.reserves;

fprintf(fd, '\n%%%%----- Reserve Data -----%%%%\n');

fprintf(fd, '%%%% reserve zones, element i, j is 1 iff gen j is in zone i\n');

fprintf(fd, '%sreserves.zones = [\n', prefix);

template = '';

for i = 1:size(r.zones, 2)

template = [template, '\t%d'];

end

template = [template, ';\n'];

fprintf(fd, template, r.zones.');

fprintf(fd, '];\n');

fprintf(fd, '\n%%%% reserve requirements for each zone in MW\n');

fprintf(fd, '%sreserves.req = [\t%g', prefix, r.req(1));

if length(r.req) > 1

fprintf(fd, ';\t%g', r.req(2:end));

end

fprintf(fd, '\t];\n');

fprintf(fd, '\n%%%% reserve costs in $/MW for each gen\n');

fprintf(fd, '%sreserves.cost = [\t%g', prefix, r.cost(1));

if length(r.cost) > 1

fprintf(fd, ';\t%g', r.cost(2:end));

end

fprintf(fd, '\t];\n');

if isfield(r, 'qty')

fprintf(fd, '\n%%%% max reserve quantities for each gen\n');

fprintf(fd, '%sreserves.qty = [\t%g', prefix, r.qty(1));

if length(r.qty) > 1

fprintf(fd, ';\t%g', r.qty(2:end));

end

fprintf(fd, '\t];\n');

end

%% save output fields for solved case ...

94

7.4 Registering the Callbacks

As seen in the ﬁxed zonal reserve example, adding a single extension to the standard

OPF formulation is often best accomplished by a set of callback functions. A typical

use case might be to run a given case with and without the reserve requirements

active, so a simple method for enabling and disabling the whole set of callbacks as a

single unit is needed.

The recommended method is to deﬁne all of the callbacks in a single ﬁle containing

a “toggle” function that registers or removes all of the callbacks depending on whether

the value of the second argument is 'on'or 'off'. The state of the registration of any

callbacks is stored directly in the mpc struct. In our example, the toggle reserves.m

ﬁle contains the toggle reserves function as well as the ﬁve callback functions.

95

function mpc = toggle_reserves(mpc, on_off)

%TOGGLE_RESERVES Enable, disable or check status of fixed reserve requirements.

% MPC = TOGGLE_RESERVES(MPC, 'on')

% MPC = TOGGLE_RESERVES(MPC, 'off')

% T_F = TOGGLE_RESERVES(MPC, 'status')

if strcmp(upper(on_off), 'ON')

% <code to check for required 'reserves' fields in mpc>

%% add callback functions

mpc = add_userfcn(mpc, 'ext2int', @userfcn_reserves_ext2int);

mpc = add_userfcn(mpc, 'formulation', @userfcn_reserves_formulation);

mpc = add_userfcn(mpc, 'int2ext', @userfcn_reserves_int2ext);

mpc = add_userfcn(mpc, 'printpf', @userfcn_reserves_printpf);

mpc = add_userfcn(mpc, 'savecase', @userfcn_reserves_savecase);

mpc.userfcn.status.dcline = 1;

elseif strcmp(upper(on_off), 'OFF')

mpc = remove_userfcn(mpc, 'savecase', @userfcn_reserves_savecase);

mpc = remove_userfcn(mpc, 'printpf', @userfcn_reserves_printpf);

mpc = remove_userfcn(mpc, 'int2ext', @userfcn_reserves_int2ext);

mpc = remove_userfcn(mpc, 'formulation', @userfcn_reserves_formulation);

mpc = remove_userfcn(mpc, 'ext2int', @userfcn_reserves_ext2int);

mpc.userfcn.status.dcline = 0;

elseif strcmp(upper(on_off), 'STATUS')

if isfield(mpc, 'userfcn') && isfield(mpc.userfcn, 'status') && ...

isfield(mpc.userfcn.status, 'dcline')

mpc = mpc.userfcn.status.dcline;

else

mpc = 0;

end

else

error('toggle_dcline: 2nd argument must be ''on'', ''off'' or ''status''');

end

Running a case that includes the ﬁxed reserves requirements is as simple as

loading the case, turning on reserves and running it.

mpc = loadcase('t_case30_userfcns');

mpc = toggle_reserves(mpc, 'on');

results = runopf(mpc);

96

7.5 Summary

The ﬁve callback stages currently deﬁned by Matpower are summarized in Ta-

ble 7-4.

Table 7-4: Callback Functions

name invoked . . . typical use

ext2int . . . from ext2int immediately after

case data is converted from external

to internal indexing.

Check consistency of input data, con-

vert to internal indexing.

formulation . . . from opf after OPF Model (om)

object is initialized with standard

OPF formulation.

Modify OPF formulation, by adding

user-deﬁned variables, constraints,

costs.

int2ext . . . from int2ext immediately before

case data is converted from internal

back to external indexing.

Convert data back to external index-

ing, populate any additional ﬁelds in

the results struct.

printpf . . . from printpf after pretty-printing

the standard OPF output.

Pretty-print any results not included

in standard OPF.

savecase . . . from savecase after printing all of

the other case data to the ﬁle.

Write non-standard case struct ﬁelds

to the case ﬁle.

7.6 Example Extensions

Matpower includes three OPF extensions implementing via callbacks, respectively,

the co-optimization of energy and reserves, interface ﬂow limits and dispatchable DC

transmission lines.

7.6.1 Fixed Zonal Reserves

This extension is a more complete version of the example of ﬁxed zonal reserve

requirements used for illustration above in Sections 7.3 and 7.4. The details of

the extensions to the standard OPF problem are given in equations (7.2)–(7.5) and

a description of the relevant input and output data structures is summarized in

Tables 7-5 and 7-6, respectively.

The code for implementing the callbacks can be found in toggle reserves. A

wrapper around runopf that turns on this extension before running the OPF is

97

Table 7-5: Input Data Structures for Fixed Zonal Reserves

name description

mpc Matpower case struct

reserves additional ﬁeld in mpc containing input parameters for zonal reserves in the

following sub-ﬁelds:

cost (ngor ngr†)×1 vector of reserve costs in $/MW, cfrom (7.3)

qty (ngor ngr†)×1 vector of reserve quantity upper bounds in MW, ith element is

rmax

i

zones nrz ×ngmatrix of reserve zone deﬁnitions

zones(k,j) =1 if gen jbelongs to reserve zone k(j∈Zk)

0 otherwise (j /∈Zk)

req nrz ×1 vector of zonal reserve requirements in MW, kth element is Rkfrom (7.5)

†Here ngr is the number of generators belonging to at least one reserve zone.

Table 7-6: Output Data Structures for Fixed Zonal Reserves

name description

results OPF results struct, superset of mpc with additional ﬁelds for output data

reserves zonal reserve data, with results in the following sub-ﬁelds:‡

Rreserve allocation ri, in MW

Rmin lower bound of riin MW from (7.2), i.e. all zeros

Rmax upper bound of riin MW from (7.2), i.e. min(rmax

i,∆i)

mu shadow prices on constraints in $/MW

lshadow price on lower bound of riin (7.2)

ushadow price on upper bound of riin (7.2)

Pmax shadow price on capacity constraint (7.4)

prc reserve price in $/MW, computed for unit ias the sum of the shadow prices of

constraints kin (7.5) in which unit iparticipates (k|i∈Zk)

‡All zonal reserve result sub-ﬁelds are ng×1 vectors, where the values are set to zero for units not participating

in the provision of reserves.

provided in runopf w res, allowing you to run a case with an appropriate reserves

ﬁeld, such as t case30 userfcns, as follows.

results = runopf_w_res('t_case30_userfcns');

See help runopf w res and help toggle reserves for more information. Exam-

ples of using this extension and a case ﬁle deﬁning the necessary input data can be

found in t opf userfcns and t case30 userfcns, respectively. Additional tests for

runopf w res are included in t runopf w res.

98