DC Optimal Power Flow
This example utilizes the power system shown in Figure 1. Similar to the AC optimal power flow, we adjust constraints and modify the topology to highlight JuliaGrid’s ability to efficiently handle such 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 powers and voltage angles, we specify the relevant units:
@power(MW, pu)
@voltage(pu, deg)Next, we define bus parameters for the analysis. This includes setting the slack bus (type = 3), where the voltage angle is fixed at zero, and specifying active power loads and shunt elements with 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 transmission line parameters by specifying reactance values. For phase-shifting transformer, we include the shift angle using the shiftAngle keyword. Additionally, we set bus voltage angle difference constraints between the from-bus and to-bus ends of each branch using minDiffAngle and maxDiffAngle keywords:
@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, but they will be introduced later in the example.
Next, we define the active power outputs of the generators, which serve as initial values for the optimization variables. Generator outputs are constrained using minActive and maxActive keywords:
@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 active power supply costs of the generators in polynomial form by setting active = 2. Then, we express the polynomial as a quadratic using the polynomial keyword:
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 is defined, we generate a DC model that includes key matrices and vectors for analysis, such as the nodal admittance matrix. This model is automatically updated when data changes and can be shared across multiple analyses:
dcModel!(system)Display Data Settings
Before running simulations, we configure the numeric format for specific data type of interest including active power flow at branches and generator outputs:
fmt = Dict("From-Bus Power" => "%.2f", "To-Bus Power" => "%.2f", "Power Output" => "%.2f")Base Case Analysis
The process starts by formulating the DC optimal power flow model and selecting the Ipopt solver. The optimization variables include bus voltage angles and generator active power outputs, denoted as θ and Pg for data printing:
analysis = dcOptimalPowerFlow(system, Ipopt.Optimizer; angle = "θ", active = "Pg")The optimization problem being solved can then be printed:
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 Pg[1] - 10 θ[1] + 5 θ[2] + 5 θ[3] = 0 Pg[2] + Pg[3] + 5 θ[1] - 10 θ[2] + 5 θ[3] = 0.202 5 θ[1] + 5 θ[2] - 15 θ[3] + 5 θ[4] = 0.2017128639793479 5 θ[3] - 5 θ[4] = 0.3072871360206521 Pg[1] ∈ [0.1, 0.655] Pg[2] ∈ [0.07, 0.20500000000000002] Pg[3] ∈ [0.07, 0.224] θ[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
Solving the DC optimal power flow model provides the bus voltage angles and generator active power outputs. The next step involves computing active power flows across buses and branches:
solve!(system, analysis)
power!(system, analysis)After obtaining the solution, bus-related results, including the optimal bus voltage angles, are examined:
printBusData(system, 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 |
|---------------------------------------------------------------------|We observe that the voltage angle difference constraint at Branch 1 reaches its upper limit. This can also be confirmed by examining the branch constraint data, where the associated dual variable takes a nonzero value:
printBranchConstraint(system, 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 active power outputs of the generators are:
printGeneratorData(system, 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 |
|---------------------------------------------|All generators operate within their active power limits, as confirmed by the generator constraint data, where all dual variables remain zero:
printGeneratorConstraint(system, 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 |
|--------------------------------------------------------|Furthermore, Generator 1, with the lowest cost, supplies most of the power, while Generator 2 and Generator 3 produce equal power amounts due to identical cost functions.
Finally, we review the results related to branch flows:
printBranchData(system, 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 |
|-----------------------------------------------------------------------|Thus, we obtained the active power flows, as illustrated in Figure 2.
Figure 2: Active power flows in the 4-bus power system for the base case scenario.
Modifying Demands
Let us now introduce a new state by updating the active power demands of consumers. These updates modify both the power system model and the DC optimal power flow model simultaneously:
updateBus!(system, analysis; label = "Bus 2", active = 25.2)
updateBus!(system, analysis; label = "Bus 4", active = 43.3)Next, we solve the DC optimal power flow again to compute the new solution without recreating the model. This step enables a warm start, as the initial primal and dual values correspond to those obtained in the base case:
solve!(system, analysis)
power!(system, analysis)Now, we observe the power output of the generators:
printGeneratorData(system, 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 to the base case, Generator 1 increases output, while Generator 2 and Generator 3 reduce production to their minimum limits. Meanwhile, all voltage angle difference constraints remain within limits:
printBranchConstraint(system, 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 |
|------------------------------------------------------|At the end of this scenario, we can review branch-related results for a more comprehensive insight into power flows:
printBranchData(system, 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 obtained results allow us to illustrate the active power flows 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 one in the system. By updating both the power system model and the DC optimal power flow model simultaneously, we enable a warm start for solving this new scenario:
cost!(system, analysis; generator = "Generator 1", active = 2, polynomial = [2.0; 40.0; 0.0])
cost!(system, analysis; generator = "Generator 3", active = 2, polynomial = [0.5; 10.0; 0.0])Next, we solve the updated problem and calculate the resulting powers:
solve!(system, analysis)
power!(system, analysis)The optimal active power outputs of the generators are as follows:
printGeneratorData(system, 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, we observe that, due to the higher cost of Generator 1, its output decreases, while the outputs of Generator 2 and Generator 3 increase. Notably, Generator 3 now produces a higher amount of power compared to Generator 2 due to its lower cost. While one might expect Generator 1 to decrease supplies more drastically and Generator 3 to increase more dramatically, this is not the case due to the need to satisfy other constraints, such as active power balance at each bus.
We can also review the results related to branches for this scenario:
printBranchData(system, 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 illustrates the power flows for this scenario. Compared to the previous scenario, Figure 4 shows that Branch 2 has significantly lower active power flow, while Branch 3 has become considerably more 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 introduce constraints on Branch 2 and Branch 3 by setting type = 1, where the active power flow at the from-bus end of these branches is limited using the maxFromBus keyword:
updateBranch!(system, analysis; label = "Branch 2", type = 1, maxFromBus = 15.0)
updateBranch!(system, analysis; label = "Branch 3", type = 1, maxFromBus = 15.0)Next, we recalculate the DC optimal power flow:
solve!(system, analysis)
power!(system, analysis)Now, let us observe the generator outputs:
printGeneratorData(system, 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 at Branch 3 forces Generator 1 to increase its active power output despite its higher cost compared to Generator 2 and Generator 3, due to the need to satisfy all constraints.
We can review the branch data and observe that the active power at the from-bus end of Branch 3 reaches the defined limit, while the power flow at Branch 2 stays within the specified limits:
printBranchData(system, 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 |
|-----------------------------------------------------------------------|Based on the obtained results, we can illustrate the power flows in Figure 5.
Figure 5: Active power flows in the 4-bus power system with added branch flow constraints.
Modifying Network Topology
At the end, we set Branch 2 out-of-service:
updateBranch!(system, analysis; label = "Branch 2", status = 0)We then recalculate the DC optimal power flow:
solve!(system, analysis)
power!(system, analysis)We can now observe the updated generator outputs:
printGeneratorData(system, 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 |
|---------------------------------------------|Due to the outage of Branch 2 and the flow limit at Branch 3, Generator 1 faces difficulties supplying load at Bus 2, reducing its output. Consequently, the only solution is to increase the output of Generator 2 and Generator 3.
Upon reviewing the branch data, we observe that the active power flows in the remaining in-service branches remain largely unchanged. This is because, following the outage of Branch 2, Generator 2 and Generator 3 have taken over the responsibility of supplying the load at Bus 2, effectively displacing Generator 1:
printBranchData(system, 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 illustrates these results under the outage of Branch 2.
Figure 6: Active power flows in the 4-bus power system with modified network topology.