DC Optimal Power Flow
This example uses the power system shown in Figure 1. As in the AC optimal power flow example, we adjust constraints and modify the topology to show how JuliaGrid efficiently handles these scenarios.
Figure 1: The 4-bus power system.
Users can download a Julia script containing the scenarios from this section using the following link.
We start by defining the unit system. Since DC optimal power flow considers only active power and voltage angles, we specify the relevant units:
@power(MW, pu)
@voltage(pu, deg)Next, we define the bus parameters, including the slack bus (type = 3), where the voltage angle is fixed at zero, active power loads, and shunt conductance values. With these definitions, we construct the power system model:
system = powerSystem()
addBus!(system; label = "Bus 1", type = 3)
addBus!(system; label = "Bus 2", active = 20.2)
addBus!(system; label = "Bus 3", conductance = 0.1)
addBus!(system; label = "Bus 4", active = 50.8)We then define the transmission line parameters by specifying reactance values. For the phase-shifting transformer, we set the shift angle with shiftAngle. We also impose bus voltage angle difference constraints between the from-bus and to-bus ends of each branch with minDiffAngle and maxDiffAngle:
@branch(reactance = 0.2, minDiffAngle = -4.1, maxDiffAngle = 4.1)
addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 3")
addBranch!(system; label = "Branch 2", from = "Bus 1", to = "Bus 2")
addBranch!(system; label = "Branch 3", from = "Bus 2", to = "Bus 3")
addBranch!(system; label = "Branch 4", from = "Bus 3", to = "Bus 4", shiftAngle = -2.3)At this stage, no active power flow constraints are imposed; they are introduced later in the example.
Next, we define the generator active power outputs, which serve as initial values for the optimization variables. The outputs are constrained with minActive and maxActive:
@generator(label = "Generator ?")
addGenerator!(system; bus = "Bus 1", active = 63.1, minActive = 10.0, maxActive = 65.5)
addGenerator!(system; bus = "Bus 2", active = 3.0, minActive = 7.0, maxActive = 20.5)
addGenerator!(system; bus = "Bus 2", active = 4.1, minActive = 7.0, maxActive = 22.4)Finally, we define the generator active power supply costs in polynomial form by setting active = 2. The quadratic coefficients are specified with polynomial:
cost!(system; generator = "Generator 1", active = 2, polynomial = [0.04; 20.0; 0.0])
cost!(system; generator = "Generator 2", active = 2, polynomial = [1.00; 20.0; 0.0])
cost!(system; generator = "Generator 3", active = 2, polynomial = [1.00; 20.0; 0.0])Once the power system data are defined, we generate the DC model, including the key matrices and vectors used in the analysis:
dcModel!(system)Display Data Settings
Before running simulations, we configure the numeric format for selected data, including branch active power flows and generator outputs:
fmt = Dict("From-Bus Power" => "%.2f", "To-Bus Power" => "%.2f", "Power Output" => "%.2f")Base Case Analysis
First, we create the DC optimal power flow model with the Ipopt solver. The optimization variables are bus voltage angles and generator active power outputs, denoted as θ and Pg:
analysis = dcOptimalPowerFlow(system, Ipopt.Optimizer; angle = "θ", active = "Pg")We can then print the optimization problem:
julia> print(analysis.method.jump)Min 400 Pg[1]² + 10000 Pg[2]² + 10000 Pg[3]² + 2000 Pg[1] + 2000 Pg[2] + 2000 Pg[3] Subject to -10 θ[1] + 5 θ[2] + 5 θ[3] + Pg[1] = 0 5 θ[1] - 10 θ[2] + 5 θ[3] + Pg[2] + Pg[3] = 0.202 5 θ[1] + 5 θ[2] - 15 θ[3] + 5 θ[4] = 0.2017128639793479 5 θ[3] - 5 θ[4] = 0.3072871360206521 θ[1] - θ[3] ∈ [-0.0715584993317675, 0.0715584993317675] θ[1] - θ[2] ∈ [-0.0715584993317675, 0.0715584993317675] θ[2] - θ[3] ∈ [-0.0715584993317675, 0.0715584993317675] θ[3] - θ[4] ∈ [-0.0715584993317675, 0.0715584993317675] θ[1] = 0 Pg[1] ≥ 0.1 Pg[2] ≥ 0.07 Pg[3] ≥ 0.07 Pg[1] ≤ 0.655 Pg[2] ≤ 0.20500000000000002 Pg[3] ≤ 0.224
We solve the DC optimal power flow model to obtain the bus voltage angles and generator active power outputs, then compute the remaining bus and branch active powers:
powerFlow!(analysis; power = true, verbose = 1)EXIT: The optimal solution was found.After obtaining the solution, we inspect the bus results, including the optimal voltage angles:
printBusData(analysis)|---------------------------------------------------------------------|
| Bus Data |
|---------------------------------------------------------------------|
| Label | Voltage | Power Generation | Power Demand | Power Injection |
| | | | | |
| Bus | Angle | Active | Active | Active |
| | [deg] | [MW] | [MW] | [MW] |
|-------|---------|------------------|--------------|-----------------|
| Bus 1 | 0.0000 | 56.4378 | 0.0000 | 56.4378 |
| Bus 2 | -2.3673 | 14.6622 | 20.2000 | -5.5378 |
| Bus 3 | -4.1000 | 0.0000 | 0.0000 | -0.0000 |
| Bus 4 | -7.6213 | 0.0000 | 50.8000 | -50.8000 |
|---------------------------------------------------------------------|The voltage angle difference constraint on Branch 1 reaches its upper limit. The branch constraint data confirm this with a nonzero associated dual variable:
printBranchConstraint(analysis; label = "Branch 1", header = true, footer = true)|------------------------------------------------------|
| Label | Voltage Angle Difference |
| | |
| | Minimum | Solution | Maximum | Dual |
| | [deg] | [deg] | [deg] | [$/deg-hr] |
|----------|---------|----------|---------|------------|
| Branch 1 | -4.1000 | 4.1000 | 4.1000 | -265.6535 |
|------------------------------------------------------|The optimal generator active power outputs are:
printGeneratorData(analysis; fmt)|---------------------------------------------|
| Generator Data |
|---------------------------------------------|
| Label | Power Output | Status |
| | | |
| Generator | Bus | Active | |
| | | [MW] | |
|-------------|-------|--------------|--------|
| Generator 1 | Bus 1 | 56.44 | 1 |
| Generator 2 | Bus 2 | 7.33 | 1 |
| Generator 3 | Bus 2 | 7.33 | 1 |
|---------------------------------------------|The generator constraint data confirm that all generators operate within their active power limits, with all dual variables equal to zero:
printGeneratorConstraint(analysis)|--------------------------------------------------------|
| Generator Constraint Data |
|--------------------------------------------------------|
| Label | Active Power Capability |
| | |
| | Minimum | Solution | Maximum | Dual |
| | [MW] | [MW] | [MW] | [$/MW-hr] |
|-------------|---------|----------|---------|-----------|
| Generator 1 | 10.0000 | 56.4378 | 65.5000 | -0.0000 |
| Generator 2 | 7.0000 | 7.3311 | 20.5000 | 0.0000 |
| Generator 3 | 7.0000 | 7.3311 | 22.4000 | 0.0000 |
|--------------------------------------------------------|Because Generator 1 has the lowest cost, it supplies most of the power. Generator 2 and Generator 3 produce equal active power because they have identical cost functions.
Finally, we inspect the branch flow results:
printBranchData(analysis; fmt)|-----------------------------------------------------------------------|
| Branch Data |
|-----------------------------------------------------------------------|
| Label | From-Bus Power | To-Bus Power | Status |
| | | | |
| Branch | From-Bus | To-Bus | Active | Active | |
| | | | [MW] | [MW] | |
|----------|----------|--------|----------------|--------------|--------|
| Branch 1 | Bus 1 | Bus 3 | 35.78 | -35.78 | 1 |
| Branch 2 | Bus 1 | Bus 2 | 20.66 | -20.66 | 1 |
| Branch 3 | Bus 2 | Bus 3 | 15.12 | -15.12 | 1 |
| Branch 4 | Bus 3 | Bus 4 | 50.80 | -50.80 | 1 |
|-----------------------------------------------------------------------|The resulting active power flows are shown in Figure 2.
Figure 2: Active power flows in the 4-bus power system for the base case scenario.
Modifying Demands
Next, we update the active power demands. These changes modify both the power system model and the DC optimal power flow model:
updateBus!(analysis; label = "Bus 2", active = 25.2)
updateBus!(analysis; label = "Bus 4", active = 43.3)We then solve the DC optimal power flow again without recreating the model. This enables a warm start because the initial primal and dual values come from the base case:
powerFlow!(analysis; power = true, verbose = 1)EXIT: The optimal solution was found.We can now inspect the generator power outputs:
printGeneratorData(analysis; fmt)|---------------------------------------------|
| Generator Data |
|---------------------------------------------|
| Label | Power Output | Status |
| | | |
| Generator | Bus | Active | |
| | | [MW] | |
|-------------|-------|--------------|--------|
| Generator 1 | Bus 1 | 54.60 | 1 |
| Generator 2 | Bus 2 | 7.00 | 1 |
| Generator 3 | Bus 2 | 7.00 | 1 |
|---------------------------------------------|Compared with the base case, Generator 1 increases its output, while Generator 2 and Generator 3 reduce production to their minimum limits. All voltage angle difference constraints remain within their limits:
printBranchConstraint(analysis)|------------------------------------------------------|
| Branch Constraint Data |
|------------------------------------------------------|
| Label | Voltage Angle Difference |
| | |
| | Minimum | Solution | Maximum | Dual |
| | [deg] | [deg] | [deg] | [$/deg-hr] |
|----------|---------|----------|---------|------------|
| Branch 1 | -4.1000 | 3.7433 | 4.1000 | -0.0000 |
| Branch 2 | -4.1000 | 2.5134 | 4.1000 | -0.0000 |
| Branch 3 | -4.1000 | 1.2299 | 4.1000 | -0.0000 |
| Branch 4 | -4.1000 | 2.6618 | 4.1000 | -0.0000 |
|------------------------------------------------------|We then inspect the branch results for additional insight into power flows:
printBranchData(analysis; fmt)|-----------------------------------------------------------------------|
| Branch Data |
|-----------------------------------------------------------------------|
| Label | From-Bus Power | To-Bus Power | Status |
| | | | |
| Branch | From-Bus | To-Bus | Active | Active | |
| | | | [MW] | [MW] | |
|----------|----------|--------|----------------|--------------|--------|
| Branch 1 | Bus 1 | Bus 3 | 32.67 | -32.67 | 1 |
| Branch 2 | Bus 1 | Bus 2 | 21.93 | -21.93 | 1 |
| Branch 3 | Bus 2 | Bus 3 | 10.73 | -10.73 | 1 |
| Branch 4 | Bus 3 | Bus 4 | 43.30 | -43.30 | 1 |
|-----------------------------------------------------------------------|The resulting active power flows are shown in Figure 3.
Figure 3: Active power flows in the 4-bus power system with modified demands.
Modifying Generator Costs
We adjust the cost functions for Generator 1 and Generator 3, making Generator 1 the highest-cost generator and Generator 3 the lowest-cost generator in the system. Updating both models enables a warm start for this scenario:
cost!(analysis; generator = "Generator 1", active = 2, polynomial = [2.0; 40.0; 0.0])
cost!(analysis; generator = "Generator 3", active = 2, polynomial = [0.5; 10.0; 0.0])Next, we solve the updated problem and compute the resulting powers:
powerFlow!(analysis; power = true, verbose = 1)EXIT: The optimal solution was found.The optimal generator active power outputs are:
printGeneratorData(analysis; fmt)|---------------------------------------------|
| Generator Data |
|---------------------------------------------|
| Label | Power Output | Status |
| | | |
| Generator | Bus | Active | |
| | | [MW] | |
|-------------|-------|--------------|--------|
| Generator 1 | Bus 1 | 25.70 | 1 |
| Generator 2 | Bus 2 | 20.50 | 1 |
| Generator 3 | Bus 2 | 22.40 | 1 |
|---------------------------------------------|In this scenario, the higher cost of Generator 1 decreases its output, while Generator 2 and Generator 3 increase their output. Because Generator 3 has the lowest cost, it produces more active power than Generator 2. Although Generator 1 might be expected to reduce output further and Generator 3 to increase output more, the solution must also satisfy constraints such as active power balance at each bus.
We can also inspect the branch results for this scenario:
printBranchData(analysis; fmt)|-----------------------------------------------------------------------|
| Branch Data |
|-----------------------------------------------------------------------|
| Label | From-Bus Power | To-Bus Power | Status |
| | | | |
| Branch | From-Bus | To-Bus | Active | Active | |
| | | | [MW] | [MW] | |
|----------|----------|--------|----------------|--------------|--------|
| Branch 1 | Bus 1 | Bus 3 | 23.03 | -23.03 | 1 |
| Branch 2 | Bus 1 | Bus 2 | 2.67 | -2.67 | 1 |
| Branch 3 | Bus 2 | Bus 3 | 20.37 | -20.37 | 1 |
| Branch 4 | Bus 3 | Bus 4 | 43.30 | -43.30 | 1 |
|-----------------------------------------------------------------------|Figure 4 shows the power flows for this scenario. Compared with the previous scenario, Branch 2 has significantly lower active power flow, while Branch 3 becomes more heavily loaded.
Figure 4: Active power flows in the 4-bus power system with modified generator costs.
Adding Branch Flow Constraints
To limit active power flow, we add constraints to Branch 2 and Branch 3 by setting type = 1 and specifying the from-bus limit with maxFromBus:
updateBranch!(analysis; label = "Branch 2", type = 1, maxFromBus = 15.0)
updateBranch!(analysis; label = "Branch 3", type = 1, maxFromBus = 15.0)Next, we solve the updated DC optimal power flow:
powerFlow!(analysis; power = true, verbose = 1)EXIT: The optimal solution was found.We can now inspect the generator outputs:
printGeneratorData(analysis; fmt)|---------------------------------------------|
| Generator Data |
|---------------------------------------------|
| Label | Power Output | Status |
| | | |
| Generator | Bus | Active | |
| | | [MW] | |
|-------------|-------|--------------|--------|
| Generator 1 | Bus 1 | 41.80 | 1 |
| Generator 2 | Bus 2 | 7.00 | 1 |
| Generator 3 | Bus 2 | 19.80 | 1 |
|---------------------------------------------|The power flow limit on Branch 3 forces Generator 1 to increase its active power output despite its higher cost than Generator 2 and Generator 3.
The branch data show that active power at the from-bus end of Branch 3 reaches the defined limit, while the power flow on Branch 2 remains within its specified limit:
printBranchData(analysis; fmt)|-----------------------------------------------------------------------|
| Branch Data |
|-----------------------------------------------------------------------|
| Label | From-Bus Power | To-Bus Power | Status |
| | | | |
| Branch | From-Bus | To-Bus | Active | Active | |
| | | | [MW] | [MW] | |
|----------|----------|--------|----------------|--------------|--------|
| Branch 1 | Bus 1 | Bus 3 | 28.40 | -28.40 | 1 |
| Branch 2 | Bus 1 | Bus 2 | 13.40 | -13.40 | 1 |
| Branch 3 | Bus 2 | Bus 3 | 15.00 | -15.00 | 1 |
| Branch 4 | Bus 3 | Bus 4 | 43.30 | -43.30 | 1 |
|-----------------------------------------------------------------------|The resulting active power flows are shown in Figure 5.
Figure 5: Active power flows in the 4-bus power system with added branch flow constraints.
Modifying Network Topology
Finally, we set Branch 2 out-of-service:
updateBranch!(analysis; label = "Branch 2", status = 0)We then solve the updated DC optimal power flow:
powerFlow!(analysis; power = true, verbose = 1)EXIT: The optimal solution was found.We can now inspect the updated generator outputs:
printGeneratorData(analysis; fmt)|---------------------------------------------|
| Generator Data |
|---------------------------------------------|
| Label | Power Output | Status |
| | | |
| Generator | Bus | Active | |
| | | [MW] | |
|-------------|-------|--------------|--------|
| Generator 1 | Bus 1 | 28.40 | 1 |
| Generator 2 | Bus 2 | 17.80 | 1 |
| Generator 3 | Bus 2 | 22.40 | 1 |
|---------------------------------------------|Because Branch 2 is out-of-service and Branch 3 is flow-limited, Generator 1 has less ability to supply the load at Bus 2, so its output decreases. As a result, Generator 2 and Generator 3 increase their output.
The branch data show that active power flows in the remaining in-service branches remain largely unchanged. After the outage of Branch 2, Generator 2 and Generator 3 supply the load at Bus 2, effectively displacing Generator 1:
printBranchData(analysis; fmt)|-----------------------------------------------------------------------|
| Branch Data |
|-----------------------------------------------------------------------|
| Label | From-Bus Power | To-Bus Power | Status |
| | | | |
| Branch | From-Bus | To-Bus | Active | Active | |
| | | | [MW] | [MW] | |
|----------|----------|--------|----------------|--------------|--------|
| Branch 1 | Bus 1 | Bus 3 | 28.40 | -28.40 | 1 |
| Branch 2 | Bus 1 | Bus 2 | 0.00 | -0.00 | 0 |
| Branch 3 | Bus 2 | Bus 3 | 15.00 | -15.00 | 1 |
| Branch 4 | Bus 3 | Bus 4 | 43.30 | -43.30 | 1 |
|-----------------------------------------------------------------------|Figure 6 shows the resulting active power flows with Branch 2 out-of-service.
Figure 6: Active power flows in the 4-bus power system with modified network topology.