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
