DC Optimal Power Flow

Similar to AC Optimal Power Flow, JuliaGrid utilizes the JuMP package to construct optimal power flow models, enabling users to manipulate these models using the standard functions provided by JuMP. JuliaGrid supports popular solvers mentioned in the JuMP documentation to solve the optimization problem.

To perform the DC optimal power flow, we first need to have the PowerSystem type that has been created with the DC model. After that, create the DCOptimalPowerFlow type to establish the DC optimal power flow framework using the function:

• dcOptimalPowerFlow.

To solve the DC optimal power flow problem and acquire generator active power outputs and bus voltage angles, make use of the following function:

• solve!.

After obtaining the solution for DC optimal power flow, JuliaGrid offers a post-processing analysis function to compute powers associated with buses and branches:

• power!.

Additionally, specialized functions are available for calculating specific types of powers for individual buses, branches, or generators.

Optimal Power Flow Model

To set up the DC optimal power flow, we begin by creating the model. To illustrate this, consider the following:

using JuMP, HiGHS

system = powerSystem()

addBus!(system; label = "Bus 1", type = 3, angle = 0.17)
addBus!(system; label = "Bus 2", active = 0.1, conductance = 0.04)
addBus!(system; label = "Bus 3", active = 0.05)

@branch(minDiffAngle = -3.1, maxDiffAngle = 3.1, minFromBus = -0.12, maxFromBus = 0.12)
addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 2", reactance = 0.05)
addBranch!(system; label = "Branch 2", from = "Bus 1", to = "Bus 3", reactance = 0.01)
addBranch!(system; label = "Branch 3", from = "Bus 2", to = "Bus 3", reactance = 0.01)

@generator(minActive = 0.0)
addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 0.6, maxActive = 0.8)
addGenerator!(system; label = "Generator 2", bus = "Bus 2", active = 0.1, maxActive = 0.3)
addGenerator!(system; label = "Generator 3", bus = "Bus 2", active = 0.2, maxActive = 0.4)

cost!(system; label = "Generator 1", active = 2, polynomial = [1100.2; 500; 80])
cost!(system; label = "Generator 2", active = 1, piecewise = [8.0 11.0; 14.0 17.0])
cost!(system; label = "Generator 3", active = 1, piecewise = [6.8 12.3; 8.7 16.8; 11.2 19.8])

dcModel!(system)

Next, the dcOptimalPowerFlow function is utilized to formulate the DC optimal power flow problem:

analysis = dcOptimalPowerFlow(system, HiGHS.Optimizer)

Optimization Variables

In DC optimal power flow, generator active power outputs are linear functions of bus voltage angles. Thus, the model's variables include generator active power outputs and bus voltage angles:

julia> JuMP.all_variables(analysis.method.jump)7-element Vector{VariableRef}:
 active[1]
 active[2]
 active[3]
 angle[1]
 angle[2]
 angle[3]
 actwise[3]

It is important to highlight that when dealing with linear piecewise cost functions comprising multiple segments, as exemplified in the case of Generator 3, JuliaGrid automatically generates helper optimization variables, such as actwise[3], and formulates a set of linear constraints to appropriately handle these cost functions.

However, in instances where a linear piecewise cost function consists of only a single segment, as demonstrated by Generator 2, the function is modelled as a standard linear function, eliminating the necessity for additional helper optimization variables.

Please note that JuliaGrid keeps references to all variables categorized into three fields:

julia> fieldnames(typeof(analysis.method.variable))(:active, :angle, :actwise)

Variable Names

Users have the option to define custom variable names for printing and writing equations, which can help present them in a more compact form. For example:

analysis = dcOptimalPowerFlow(system, HiGHS.Optimizer; active = "P", angle = "θ")

Add Variables

The user has the ability to easily add new variables to the defined DC optimal power flow model by using the @variable macro from the JuMP package:

JuMP.@variable(analysis.method.jump, newVariable)

We can verify that the new variable is included in the defined model by using the function:

julia> JuMP.is_valid(analysis.method.jump, newVariable)true

Delete Variables

To delete a variable, the delete function from the JuMP package can be used:

JuMP.delete(analysis.method.jump, newVariable)

After deletion, the variable is no longer part of the model:

julia> JuMP.is_valid(analysis.method.jump, newVariable)false

Constraint Functions

JuliGrid keeps track of all the references to internally formed constraints in the constraint field of the DCOptimalPowerFlow type. These constraints are divided into six fields:

julia> fieldnames(typeof(analysis.method.constraint))(:slack, :balance, :voltage, :flow, :capability, :piecewise)

Slack Bus Constraint

The slack field contains a reference to the equality constraint associated with the fixed bus voltage angle value of the slack bus. This constraint is set within the addBus! function using the angle keyword:

julia> print(system.bus.label, analysis.method.constraint.slack.angle)Bus 1: θ[1] = 0.17

Users have the flexibility to modify this constraint by changing which bus serves as the slack bus and by adjusting the value of the bus angle. This can be achieved using the updateBus! function, for example:

updateBus!(system, analysis; label = "Bus 1", angle = -0.1)

Subsequently, the updated slack constraint can be inspected as follows:

julia> print(system.bus.label, analysis.method.constraint.slack.angle)Bus 1: θ[1] = -0.1

Bus Active Power Balance Constraints

The balance field contains references to the equality constraints associated with the active power balance equations defined for each bus. The constant terms in these equations are determined by the active and conductance keywords within the addBus! function. Additionally, if there are phase shift transformers in the system, the constant terms can also be affected by the shiftAngle keyword within the addBranch! function:

julia> print(system.bus.label, analysis.method.constraint.balance.active)Bus 1: P[1] - 120 θ[1] + 20 θ[2] + 100 θ[3] = 0
Bus 2: P[2] + P[3] + 20 θ[1] - 120 θ[2] + 100 θ[3] = 0.14
Bus 3: 100 θ[1] + 100 θ[2] - 200 θ[3] = 0.05

During the execution of functions that add or update power system components, these constraints are automatically adjusted to reflect the current configuration of the power system, for example:

updateBus!(system, analysis; label = "Bus 3", active = 0.1)
updateGenerator!(system, analysis; label = "Generator 2", status = 0)

Subsequently, the updated set of active power balance constraints can be examined as follows:

julia> print(system.bus.label, analysis.method.constraint.balance.active)Bus 1: P[1] - 120 θ[1] + 20 θ[2] + 100 θ[3] = 0
Bus 2: P[3] + 20 θ[1] - 120 θ[2] + 100 θ[3] = 0.14
Bus 3: 100 θ[1] + 100 θ[2] - 200 θ[3] = 0.1

Bus Voltage Angle Difference Constraints

The voltage field contains references to the inequality constraints associated with the minimum and maximum bus voltage angle difference between the from-bus and to-bus ends of each branch. These values are specified using the minDiffAngle and maxDiffAngle keywords within the addBranch! function:

julia> print(system.branch.label, analysis.method.constraint.voltage.angle)Branch 1: θ[1] - θ[2] ∈ [-3.1, 3.1]
Branch 2: θ[1] - θ[3] ∈ [-3.1, 3.1]
Branch 3: θ[2] - θ[3] ∈ [-3.1, 3.1]

Additionally, by employing the updateBranch! function, we have the ability to modify these constraints as follows:

updateBranch!(system, analysis; label = "Branch 1", minDiffAngle = -1.7, maxDiffAngle = 1.7)

Subsequently, the updated set of voltage angle difference constraints can be examined as follows:

julia> print(system.branch.label, analysis.method.constraint.voltage.angle)Branch 1: θ[1] - θ[2] ∈ [-1.7, 1.7]
Branch 2: θ[1] - θ[3] ∈ [-3.1, 3.1]
Branch 3: θ[2] - θ[3] ∈ [-3.1, 3.1]

Branch Active Power Flow Constraints

The flow field refers to the inequality constraints associated with active power flow limits at the from-bus end of each branch. These limits are set using the minFromBus and maxFromBus keywords in the addBranch! function:

julia> print(system.branch.label, analysis.method.constraint.flow.active)Branch 1: 20 θ[1] - 20 θ[2] ∈ [-0.12, 0.12]
Branch 2: 100 θ[1] - 100 θ[3] ∈ [-0.12, 0.12]
Branch 3: 100 θ[2] - 100 θ[3] ∈ [-0.12, 0.12]

By employing the updateBranch! function, we have the ability to modify these specific constraints, for example:

updateBranch!(system, analysis; label = "Branch 1", status = 0)
updateBranch!(system, analysis; label = "Branch 2", reactance = 0.03, maxFromBus = 0.14)

Subsequently, the updated set of active power flow constraints can be examined as follows:

julia> print(system.branch.label, analysis.method.constraint.flow.active)Branch 2: 33.333333333333336 θ[1] - 33.333333333333336 θ[3] ∈ [-0.12, 0.14]
Branch 3: 100 θ[2] - 100 θ[3] ∈ [-0.12, 0.12]

Generator Active Power Capability Constraints

The capability field contains references to the inequality constraints associated with the minimum and maximum active power outputs of the generators. These limits are specified using the minActive and maxActive keywords within the addGenerator! function:

julia> print(system.generator.label, analysis.method.constraint.capability.active)Generator 1: P[1] ∈ [0, 0.8]
Generator 2: P[2] = 0
Generator 3: P[3] ∈ [0, 0.4]

As demonstrated, the active power output of Generator 2 is currently fixed at zero due to the earlier action of setting this generator out-of-service. Let us adjust this specific constraint using the updateGenerator! function:

updateGenerator!(system, analysis; label = "Generator 2", status = 1, maxActive = 0.5)

Subsequently, the updated set of active power capability constraints can be examined as follows:

julia> print(system.generator.label, analysis.method.constraint.capability.active)Generator 1: P[1] ∈ [0, 0.8]
Generator 2: P[2] ∈ [0, 0.5]
Generator 3: P[3] ∈ [0, 0.4]

It is important to note that bringing back Generator 2 into service will also have an impact on the balance constraint, which will once again be influenced by the generator's output.

Active Power Piecewise Constraints

In the context of active power modelling, the piecewise field serves as a reference to the inequality constraints related to linear piecewise cost functions. These constraints are created using the cost! function with active = 1 specified when dealing with linear piecewise cost functions comprising multiple segments. JuliaGrid takes care of establishing the appropriate inequality constraints for each segment of the linear piecewise cost:

julia> print(system.generator.label, analysis.method.constraint.piecewise.active)Generator 3: 2.3684210526315796 P[3] - actwise[3] ≤ 3.805263157894739
Generator 3: 1.2 P[3] - actwise[3] ≤ -6.360000000000001

It is worth noting that these constraints can also be automatically updated using the cost! function, and readers can find more details in the section about the objective function.

Add Constraints

Users can effortlessly introduce additional constraints into the defined DC optimal power flow model by utilizing the addBranch! or addGenerator! functions. Specifically, if a user wishes to include a new branch or generator in an already defined PowerSystem and DCOptimalPowerFlow type:

addBranch!(system, analysis; label = "Branch 4", from = "Bus 1", to = "Bus 2", reactance = 1)
addGenerator!(system, analysis; label = "Generator 4", bus = "Bus 1", maxActive = 0.2)

As a result, the flow and capability constraints will be adjusted as follows:

julia> print(system.branch.label, analysis.method.constraint.flow.active)Branch 2: 33.333333333333336 θ[1] - 33.333333333333336 θ[3] ∈ [-0.12, 0.14]
Branch 3: 100 θ[2] - 100 θ[3] ∈ [-0.12, 0.12]
Branch 4: θ[1] - θ[2] ∈ [-0.12, 0.12]
julia> print(system.generator.label, analysis.method.constraint.capability.active)Generator 1: P[1] ∈ [0, 0.8]
Generator 2: P[2] ∈ [0, 0.5]
Generator 3: P[3] ∈ [0, 0.4]
Generator 4: active[4] ∈ [0, 0.2]

Add User-Defined Constraints

Users also have the option to include their custom constraints within the established DC optimal power flow model by employing the @constraint macro. For example, the addition of a new constraint can be achieved as follows:

JuMP.@constraint(analysis.method.jump, 0.0 <= analysis.method.variable.active[4] <= 0.3)

Delete Constraints

To delete a constraint, users can make use of the delete function from the JuMP package. When handling constraints that have been internally created, users can refer to the constraint references stored in the constraint field of the DCOptimalPowerFlow type.

For example, if the intention is to eliminate constraints related to the capability of Generator 4, the following code snippet can be employed:

JuMP.delete(analysis.method.jump, analysis.method.constraint.capability.active[4])

Objective Function

The objective function of the DC optimal power flow is constructed using polynomial and linear piecewise cost functions of the generators, which are defined using the cost! functions. It is important to note that only polynomial cost functions up to the second degree are included in the objective. If there are polynomials of higher degrees, JuliaGrid will exclude them from the objective function.

In the provided example, the objective function that needs to be minimized to obtain the optimal values of the active power outputs of the generators and the bus voltage angles is as follows:

julia> JuMP.objective_function(analysis.method.jump)1100.2 P[1]² + 500 P[1] + actwise[3] + P[2] + 83

Additionally, JuliaGrid stores the objective function in a separate variable, allowing users to access it by referencing the variable analysis.objective.

Update Objective Function

By utilizing the cost! functions, users have the flexibility to modify the objective function by adjusting polynomial or linear piecewise cost coefficients or by changing the type of polynomial or linear piecewise function employed. For instance, consider Generator 3, which incorporates a piecewise cost structure with two segments. Now, we can define a polynomial function for this generator and activate it by specifying the keyword active = 2 as shown:

cost!(system, analysis; label = "Generator 3", active = 2, polynomial = [853.4; 257; 40])

This results in the updated objective function, which can be observed as follows:

julia> analysis.method.objective1100.2 P[1]² + 853.4 P[3]² + 500 P[1] + P[2] + 257 P[3] + 123

User-Defined Objective Function

Users can modify the objective function using the set_objective_function function from the JuMP package. This operation is considered destructive because it is independent of power system data; however, in certain scenarios, it may be more straightforward than using the cost! function for updates. Moreover, using this methodology, users can combine a defined function with a newly defined expression. Here is an example of how it can be achieved:

expr = 100.2 * analysis.method.variable.active[1] * analysis.method.variable.active[1] + 123
JuMP.set_objective_function(analysis.method.jump, analysis.method.objective - expr)

We can now observe the updated objective function as follows:

julia> JuMP.objective_function(analysis.method.jump)1000 P[1]² + 853.4 P[3]² + 500 P[1] + P[2] + 257 P[3]

Setup Starting Values

In JuliaGrid, the assignment of starting primal and dual values for optimization variables takes place when the solve! function is executed.

Starting Primal Values

Starting primal values are determined based on the generator and voltage fields within the DCOptimalPowerFlow type. By default, these values are initially established using the active power outputs of the generators and the initial bus voltage angles:

julia> print(system.generator.label, analysis.power.generator.active)Generator 1: 0.6
Generator 2: 0.1
Generator 3: 0.2
Generator 4: 0.0
julia> print(system.bus.label, analysis.voltage.angle)Bus 1: -0.1
Bus 2: 0.0
Bus 3: 0.0

Users have the flexibility to adjust these values according to their specifications, which will then be used as the starting primal values when executing the solve! function.

Using DC Power Flow

In this perspective, users have the capability to conduct the DC power flow analysis and leverage the resulting solution to configure starting primal values. Here is an illustration of how this can be achieved:

flow = dcPowerFlow(system)
solve!(system, flow)

After obtaining the solution, we can calculate the active power outputs of the generators and utilize the bus voltage angles to set the starting values. In this case, the generator and voltage fields of the DCOptimalPowerFlow type can be employed to store the new starting values:

for (key, idx) in system.generator.label
    analysis.power.generator.active[idx] = generatorPower(system, flow; label = key)
end

for i = 1:system.bus.number
    analysis.voltage.angle[i] = flow.voltage.angle[i]
end

Starting Dual Values

Dual variables, often referred to as Lagrange multipliers or Kuhn-Tucker multipliers, represent the shadow prices or marginal costs associated with constraints. The assignment of initial dual values occurs when the solve! function is executed. Initially, the starting dual values are unknown, but users can access and manually set them. For example:

analysis.method.dual.balance.active[1] = 0.4

Optimal Power Flow Solution

To establish the DC optimal power flow problem, we can utilize the dcOptimalPowerFlow function. After setting up the problem, we can use the solve! function to compute the optimal values for the active power outputs of the generators and the bus voltage angles. Also, to turn off the solver output within the REPL, we use the set_silent function before calling solve! function. Here is an example:

JuMP.set_silent(analysis.method.jump)
solve!(system, analysis)

By executing this function, we will obtain the solution with the optimal values for the active power outputs of the generators and the bus voltage angles:

julia> print(system.generator.label, analysis.power.generator.active)Generator 1: 0.0
Generator 2: 0.09539999999999843
Generator 3: 0.0
Generator 4: 0.14460000000000084
julia> print(system.bus.label, analysis.voltage.angle)Bus 1: -0.1
Bus 2: -0.10460000000000004
Bus 3: -0.10420000000000001

Objective Value

To obtain the objective value of the optimal power flow solution, we can use the objective_value function:

julia> JuMP.objective_value(analysis.method.jump)0.09539999999999843

Dual Variables

The values of the dual variables are stored in the dual field of the DCOptimalPowerFlow type. For example:

julia> analysis.method.dual.balance.active[1]1.4459979675179438e-8

Print Results in the REPL

Users can utilize the functions printBusData and printGeneratorData to display results. Additionally, the functions listed in the Print Constraint Data section allow users to print constraint data related to buses, branches, or generators in the desired units. For example:

@power(MW, MVAr, pu)
printBusConstraint(system, analysis)

|------------------------------|
| Bus Constraint Data          |
|------------------------------|
| Label | Active Power Balance |
|       |                      |
|       | Solution |      Dual |
|       |     [MW] | [$/MW-hr] |
|-------|----------|-----------|
| Bus 1 |   0.0000 |    0.0000 |
| Bus 2 |  14.0000 |    0.0100 |
| Bus 3 |  10.0000 |    0.0101 |
|------------------------------|

Next, users can easily customize the print results for specific constraint, for example:

printBusConstraint(system, analysis; label = "Bus 1", header = true)
printBusConstraint(system, analysis; label = "Bus 2", footer = true)

Save Results to a File

Users can also redirect print output to a file. For example, data can be saved in a text file as follows:

open("bus.txt", "w") do file
    printBusConstraint(system, analysis, file)
end

Save Results to a CSV File

For CSV output, users should first generate a simple table with style = false, and then save it to a CSV file:

using CSV

io = IOBuffer()
printBusConstraint(system, analysis, io; style = false)
CSV.write("constraint.csv", CSV.File(take!(io); delim = "|"))

Primal and Dual Warm Start

Utilizing the DCOptimalPowerFlow type and proceeding directly to the solver offers the advantage of a "warm start". In this scenario, the starting primal and dual values for the subsequent solving step correspond to the solution obtained from the previous step.

Primal Variables

In the previous example, the following solution was obtained, representing the values of the primal variables:

julia> print(system.generator.label, analysis.power.generator.active)Generator 1: 0.0
Generator 2: 0.09539999999999843
Generator 3: 0.0
Generator 4: 0.14460000000000084
julia> print(system.bus.label, analysis.voltage.angle)Bus 1: -0.1
Bus 2: -0.10460000000000004
Bus 3: -0.10420000000000001

Dual Variables

We also obtained all dual values. Here, we list only the dual variables for one type of constraint as an example:

julia> print(system.branch.label, analysis.method.dual.flow.active)Branch 2: -1.039999994152338
Branch 3: 0.0
Branch 4: 0.0

Modify Optimal Power Flow

Now, let us introduce changes to the power system from the previous example:

updateGenerator!(system, analysis; label = "Generator 2", maxActive = 0.08)

Next, we want to solve this modified optimal power flow problem. If we use solve! at this point, the primal and dual starting values will be set to the previously obtained values:

solve!(system, analysis)

As a result, we obtain a new solution:

julia> print(system.generator.label, analysis.power.generator.active)Generator 1: 0.0
Generator 2: 0.0008
Generator 3: 0.09459999999999802
Generator 4: 0.1446000000000005
julia> print(system.bus.label, analysis.voltage.angle)Bus 1: -0.1
Bus 2: -0.10460000000000004
Bus 3: -0.10420000000000001

Reset Primal and Dual Values

Users retain the flexibility to reset these initial primal values to their default configurations at any juncture. This can be accomplished by utilizing the active power outputs of the generators and the initial bus voltage angles extracted from the PowerSystem type, employing the startingPrimal! function:

startingPrimal!(system, analysis)

The primal starting values will now be identical to those that would be obtained if the dcOptimalPowerFlow function were executed after all the updates have been applied.

Using the startingDual! function, users can clear all dual variable values, resetting them to their default state:

startingDual!(system, analysis)

Power Analysis

After obtaining the solution from the DC optimal power flow, we can calculate powers related to buses and branches using the power! function. For instance, let us consider the power system for which we obtained the DC optimal power flow solution:

using HiGHS

system = powerSystem()

addBus!(system; label = "Bus 1", type = 3, angle = 0.17)
addBus!(system; label = "Bus 2", active = 0.1, conductance = 0.04)
addBus!(system; label = "Bus 3", active = 0.05)

@branch(minDiffAngle = -pi, maxDiffAngle = pi, minFromBus = -0.12, maxFromBus = 0.12)
addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 2", reactance = 0.05)
addBranch!(system; label = "Branch 2", from = "Bus 1", to = "Bus 3", reactance = 0.01)
addBranch!(system; label = "Branch 3", from = "Bus 2", to = "Bus 3", reactance = 0.01)

@generator(minActive = 0.0)
addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 3.2, maxActive = 0.5)
addGenerator!(system; label = "Generator 2", bus = "Bus 2", active = 0.2, maxActive = 0.2)

cost!(system; label = "Generator 1", active = 2, polynomial = [1100.2; 500; 80])
cost!(system; label = "Generator 2", active = 1, piecewise = [10.8 12.3; 14.7 16.8])

analysis = dcOptimalPowerFlow(system, HiGHS.Optimizer)
solve!(system, analysis)

Now we can calculate the active powers using the following function:

power!(system, analysis)

Finally, to display the active power injections and from-bus active power flows, we can use the following code:

julia> print(system.bus.label, analysis.power.injection.active)Bus 1: 0.0
Bus 2: 0.0008999999999991246
Bus 3: -0.0005
julia> print(system.branch.label, analysis.power.from.active)Branch 1: -7.142857142816705e-5
Branch 2: 7.142857142983239e-5
Branch 3: 0.00042857142857066766

Print Results in the REPL

Users can utilize any of the print functions outlined in the Print Power System Data or Print Power System Summary. For example, to create a bus data with the desired units, users can use the following function:

@voltage(pu, deg, V)
@power(MW, MVAr, pu)
printBusData(system, analysis)

|---------------------------------------------------------------------|
| Bus Data                                                            |
|---------------------------------------------------------------------|
| Label | Voltage | Power Generation | Power Demand | Power Injection |
|       |         |                  |              |                 |
|   Bus |   Angle |           Active |       Active |          Active |
|       |   [deg] |             [MW] |         [MW] |            [MW] |
|-------|---------|------------------|--------------|-----------------|
| Bus 1 |  9.7403 |           0.0000 |       0.0000 |          0.0000 |
| Bus 2 |  9.7405 |           0.1900 |       0.1000 |          0.0900 |
| Bus 3 |  9.7402 |           0.0000 |       0.0500 |         -0.0500 |
|---------------------------------------------------------------------|

Active Power Injection

To calculate active power injection associated with a specific bus, the function can be used:

julia> active = injectionPower(system, analysis; label = "Bus 2")0.0008999999999991246

Active Power Injection from Generators

To calculate active power injection from the generators at a specific bus, the function can be used:

julia> active = supplyPower(system, analysis; label = "Bus 2")0.0018999999999991246

Active Power Flow

Similarly, we can compute the active power flow at both the from-bus and to-bus ends of the specific branch by utilizing the provided functions below:

julia> active = fromPower(system, analysis; label = "Branch 2")7.142857142983239e-5
julia> active = toPower(system, analysis; label = "Branch 2")-7.142857142983239e-5