## Description

This assignment will provide hands-on practice for optimization with application to the economic dispatch

problem in power systems. The assignment is organized in a tutorial fashion, thereby allowing you to practice

optimization theory on a relevant real-world energy system example.

## Background

In this assignment you will practice optimization with the DistFlow equations of Baran & Wu [1] (1)-(9).

The DistFlow equations model power flow through radial distribution power networks. The beauty of the

DistFlow equations is that they are relatively simple (i.e. they yield convex programs). Yet, they model

active & reactive (AC) power, branch power flow limits, node voltage limits, etc. Despite their simplicity, the

DistFlow equations are much more sophisticated that the “supply = demand” single equation used by many

energy policy researchers. In this assignment, your objective is to optimally schedule distributed energy

generators in a distribution feeder to minimize economic costs, while maintaining safe operating constraints.

0

5 4 1 2 3

6

7

11 10 8 9

12

Figure 1: IEEE 13-Node Test Feeder

Consider the IEEE 13-node Test Feeder shown in

Fig. 1, adopted from [2]. We mathematically represent this test feeder by a radial undirected graph

G = (N ,L). Symbol N = {0, 1, 2, · · · , 12} represents the set of nodes (a.k.a. “buses” in power systems jargon). Symbol L ⊂ N × N represents the

set of edges (a.k.a. “lines” in power systems jargon). For notational convenience, we introduce the

concepts of parent nodes and the adjacency matrix.

The parent of j, ρ(j), is the adjacent node in the

direction toward node 0. The adjacency matrix, A,

encodes the network topology. Mathematically, it is

a 13 × 13 matrix of zeros and ones. Aij = 1 if node

i is the parent node of node j. Otherwise, Aij = 0.

At each node j ∈ N , we have l

P

j

, l

Q

j

active and reactive power consumed, respectively.

Additionally,

pj , qj active and reactive power are generated, respectively. Note that total power is complex number (due to AC power flow), where the real part is

given by active power pj and the imaginary part is given by reactive power qj . The magnitude of complex

power is given by sj in equation (6).

Voltage at each node is a complex number. We mathematically model this by the squared voltage magnitude

Page 1 of 6

CE 295: Energy Systems and Control

Vj . In general, we must schedule generators to regulate nodal voltages around a nominal value (see (8)).

Each edge (i, j) (or “line”) has a characteristic impedance, decomposed into resistance rij and reactance xij .

The active and reactive power flowing along line (i, j) is given by Pij and Qij , respectively. Finally, the

squared magnitude of complex current on line (i, j) is given by Lij . Lines have a maximum current capacity,

given by constraint (9).

Power is supplied in two ways. First, power can be imported from the transmission grid connected to node

0. Second, power can be generated from distributed generators within the feeder network. We have a gas

generator on node 3, and solar generator on node 9.

Your objective is to supply all electricity demand at minimum cost. Simultaneously, you must ensure all

generator output powers, nodal voltages, and line currents are within safe operating bounds.

Table 1: Nomenclature

Symbol Description Units

N Set of nodes (a.k.a. buses) [-]

E Set of edges (a.k.a. lines) [-]

ρ(i) Parent of node i, e.g. ρ(6) = {1}, ρ(12) = {10} [-]

A Adjacency Matrix (13×13) encoding network structure [-]

pi active power generated at node i [MW]

qi reactive power generated at node i [MVAr]

si apparent power generated at node i [MVA]

si,max apparent power generation capacity at node i [MVA]

l

P

i

, lQ

i

active, reactive power consumed at node i [MW], [MVAr]

Vi Squared voltage magnitude at node i [p.u.]

vmin, vmax Minimum, maximum nodal voltage [p.u.]

ci Marginal cost of apparent power generation at node i [USD/MVA]

rij resistance of line (i, j) [p.u.]

xij reactance of line (i, j) [p.u.]

Pij active power flowing on line (i, j) [MW]

Qij reactive power flowing on line (i, j) [MVAr]

Lij Squared magnitude of complex current on line (i, j) [p.u.]

Iij Maximum magnitude of complex current on line (i, j) [p.u.]

WA, WB Uncertain & weather-dependent power capacity of PV panels A,B [MVA]

σA, σB Percentage power output of PV panels A,B [%]

Note: The acronym “p.u.” stands for per-unit. It’s power systems jargon for “normalized to a unitless

quantity”. For your convenience, all the data has been normalized to simplify your analysis. For interested

readers, the base parameters for normalization in this HW are Sbase = 1 MW, Vbase = 4.17 kV.

Page 2 of 6

CE 295: Energy Systems and Control

Pij = (l

P

j − pj ) + rijLij +

X

k∈N

AjkPjk ∀ j ∈ N , i = ρ(j) (1)

Qij = (l

Q

j − qj ) + xijLij +

X

k∈N

AjkQjk ∀ j ∈ N , i = ρ(j) (2)

Vj = Vi + (r

2

ij + x

2

ij )Lij − 2(rijPij + xijQij ) ∀ j ∈ N , i = ρ(j) (3)

Lij =

Pij

2 + Qij

2

Vj

∀ j ∈ N , i = ρ(j) (4)

pj ≥ 0, qj ≥ 0 ∀ j ∈ N (5)

q

pj

2 + qj

2 = k[pj , qj ]k2 = sj ∀ j ∈ N (6)

sj ≤ sj,max ∀ j ∈ N (7)

v

2

min ≤ Vj ≤ v

2

max ∀ j ∈ N (8)

Lij ≤ I

2

ij,max ∀ j ∈ N , i = ρ(j) (9)

## Problem 1: Network Parameters

Download and open the data file HW_Data.xls. Copy over the test network parameters from the xls file

into the Matlab/Python skeleton code file.

(a) Create a bar plot of the active and reactive power consumption. Make the x-axis the node index

number, while the y-axis is the power consumption. Place the active & reactive powers side-by-side

for each node. Add a legend.

(b) Fill out the adjacency matrix with zeros and ones. Include in your report.

## Problem 2: Balancing Supply & Demand without a Network

Start simple: In this problem, we will optimally dispatch our generators to minimize cost, while disregarding

the network completely. That is, we seek to balance active & reactive power supply & demand, while

minimizing generation cost and completely ignoring line losses and constraints.

(a) What are the optimization variables?

(b) Write down the objective function, using the notation in Table 1.

(c) Write down ALL the constraints, using the notation from Table 1. Label the physical meaning of each

constraint. For this problem, ignore voltages and all constraints associated with line flows.

(d) Is this a linear program (LP), quadratic program (QP), or convex program (CP)? Why or why not?

What happens when we relax (6) from an equality constraint to an inequality (≤) constraint?

(e) Code and numerically solve this problem in Matlab or Python using cvx or cvxpy, respectively. Use

the relaxed version of (6), so your optimization program is convex. In your report provide (i) the

optimal active & reactive generator powers, and (ii) the minimum generating cost.

Page 3 of 6

CE 295: Energy Systems and Control

## Problem 3: Add Line Power Flows

Next, we add line power flows Pij , Qij but still neglect the nodal voltage Vj and Lij terms.

(a) What are the optimization variables?

(b) Write down ALL the constraints, using the notation from Table 1. Label the physical meaning of each

constraint. For this problem, ignore (3)-(4) and drop the Lij terms from (1)-(2).

(c) Code and numerically solve this problem in Matlab or Python using cvx or cvxpy, respectively. In your

report provide (i) the optimal active & reactive generator powers, and (ii) the minimum generating

cost. Should the minimum and minimizers be different or the same as Problem 2? Why?

(d) In your code, declare a dual variable µs corresponding to the inequality (7). Re-compute the optimal

solution and dual variable. If we could increase the solar power capacity by 1MW, then how much

money would we save?

## Problem 4: The Complete Optimal Economic Dispatch with DistFlow Equations

Now we add the nodal voltages Vj , squared current magnitudes Lij , and their bounds (8)-(9). This incorporates impedance (i.e. losses) across the network, along with nodal voltage and line transmission limits.

(a) What are the optimization variables?

(b) Write down ALL the constraints, using the notation from Table 1. Label the physical meaning of each

constraint.

(c) Is this a convex program (CP)? Why or why not? What happens when we relax (4) from an equality

constraint to an inequality (≥) constraint?

(d) Code and numerically solve this problem in Matlab or Python using cvx or cvxpy, respectively. In your

report provide (i) the optimal active & reactive generator powers, and (ii) the minimum generating

cost. Use vmin = 0.95, vmax = 1.05 as your voltage limits. Is the solution equivalent to the solution in

Problem 3? Why or why not?

(e) Use the dual variables to determine which constraints are active. Specifically, identify at which nodes

the voltage constraint (8) and line current constraint (9) are active.

(f) Now re-solve the problem with vmin = 0.98, vmax = 1.02 as your voltage limits. In your report provide

(i) the optimal active & reactive generator powers, and (ii) the minimum generating cost. Why did

the solution change?

## Problem 5: [OPTIONAL] Robust Economic Dispatch with Renewables.

Next, we explore robust

economic dispatch in the face of uncertain renewable generation. Imagine the solar generator at node 9 is

comprised of two solar panels, A and B. The output of each panel is uncertain, and weather dependent.

Mathematically, we write:

s9 ≤ WA · σa + WB · σb (10)

where σa, σb ∈ [0, 1] represent the percentage output of each panel. Symbols WA, WB are random variables

that represent the uncertain power capacity of each panel, due to weather. We hypothesize that WA, WB

Page 4 of 6

CE 295: Energy Systems and Control

vary between 1 MW and 1.5 MW. Our goal, in this problem, is to optimally schedule the generators in the

face of uncertain solar capacity. Note, this problem is similar to Example 3.4 in the CH3 notes.

(a) Re-arrange (10) into the form a

T y ≤ b. What are a, y, and b? Hint: a ∈ R

3

, y ∈ R

3

, b ∈ R.

(b) We now hypothesize that vector a is uncertain, but lies within an ellipsoid

a ∈ E = {a¯ + Eu | kuk2 ≤ 1} (11)

Provide the values for a¯ and E. Hint #1: a¯ ∈ R

3

represents the center of the ellipsoid. Hint #2:

E 0 ∈ R

3×3

encodes the lengths of the semi-axes. If E is diagonal, then the diagonal elements

represent the semi-axis lengths along each coordinate of a.

(c) The robust version of (10) is

a

T

y ≤ b, ∀ a ∈ E (12)

Convert this robust linear inequality into a second-order cone constraint. In your report, it is only

necessary to provide the final written form of the second order cone constraint.

(d) In your report, write down the new robust optimization problem. What are the optimization variables?

Write down ALL the constraints, using the notation from Table 1. Label the physical meaning of each

constraint.

(e) Now solve the optimization program with vmin = 0.95, vmax = 1.05 as your voltage limits. In your

report provide (i) the optimal active & reactive generator powers, and (ii) the minimum generating

cost. How did the solution change?

## Interesting Remarks

• Power system operators utilize day-ahead load predictions to solve the economic dispatch problem and

procure generation on an hourly basis. Mismatch between predicted and actual load is compensated

by “spinning reserves”.

• This homework is not trivial. However, by completing it, you have quickly acquired sophisticated

knowledge and skills in power systems and optimization. DistFlow equations, convex optimization,

and robust optimization are skills one normally finds in power systems / optimization PhD students.

Well done!

## Deliverables

Submit the following on bCourses. Be sure that the function files are named exactly as specified (including

spelling and case). Please do NOT zip files.

LASTNAME_FIRSTNAME_HW3.PDF

LASTNAME_FIRSTNAME_HW3.xyz which contains your respective Matlab or Python files, i.e. xyz ∈ {m, ipynb}.

Page 5 of 6

CE 295: Energy Systems and Control

## How to Download and Install CVX

Matlab Instructions

• Go to http://cvxr.com/cvx/download/

• In the “Download Matrix”, focus your attention on the “Standard bundles, including Gurobi and/or

MOSEK”.

• Click the package corresponding to your system. For example, I have a MacBook Pro Early 2015, so

I’ll select mexmaci64. You can download a .zip or .tar.gz file.

• Download and unpack anywhere you like. The Downloads folder is a reasonable option.

• Start Matlab.

• Change directories to the top of the CVX distribution. Hint: Use command » cd /Users/scottmoura/

Downloads/cvx

• Run Matlab command » cvx_setup

• The cvx_setup command runs a variety of checks to verify your installation is correct.

• After successfully installing CVX, please go to http://cvxr.com/news/2014/02/cvx-demo-video/ to

watch the Getting Started Demo from Prof. Stephen Boyd.

• To confirm a successful understanding, try coding and solving the example shown at http://cvxr.com/cvx/

## Python Instructions

• Ensure you have installed Anaconda, as recommended in this class

• Go to https://www.cvxpy.org/install/index.html

• Follow the instructions corresponding to your system

• After successfully installing CVXPY, fire-up the iPython notebook. You can try the Portfolio optimization example linked here.

## References

[1] M. Baran and F. Wu, “Network reconfiguration in distribution systems for loss reduction and load

balancing,” IEEE Transactions on Power Delivery, vol. 4, no. 2, pp. 1401–1407, apr 1989. [Online]. Available:

http://ieeexplore.ieee.org/document/25627/

[2] J. Fuller, Y. Xu, B. Kersting, R. Dugan, and S. Carneiro Jr., “IEEE PES Distribution Test Feeders,” 2010.

[Online]. Available: https://ewh.ieee.org/soc/pes/dsacom/testfeeders/

Page 6 of 6