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 composite type that has been created with the DC model. After that, create the DCOptimalPowerFlow composite 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, there are specialized functions dedicated to calculating specific types of active powers related to particular buses or branches:

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, longTerm = 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)

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. Here is an example:

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 composite type. These constraints are divided into six fields:

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

Info

We suggest that readers refer to the tutorial on DC Optimal Power Flow for insights into the implementation.

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: angle[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: angle[1] = -0.1

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: active[1] - 120 angle[1] + 20 angle[2] + 100 angle[3] = 0
Bus 2: active[2] + active[3] + 20 angle[1] - 120 angle[2] + 100 angle[3] = 0.14
Bus 3: 100 angle[1] + 100 angle[2] - 200 angle[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: active[1] - 120 angle[1] + 20 angle[2] + 100 angle[3] = 0
Bus 2: active[3] + 20 angle[1] - 120 angle[2] + 100 angle[3] = 0.14
Bus 3: 100 angle[1] + 100 angle[2] - 200 angle[3] = 0.1

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: angle[1] - angle[2] ∈ [-3.1, 3.1]
Branch 2: angle[1] - angle[3] ∈ [-3.1, 3.1]
Branch 3: angle[2] - angle[3] ∈ [-3.1, 3.1]

Info

Please note that if the limit constraints are set to minDiffAngle = -2π and maxDiffAngle = 2π for the corresponding branch, JuliGrid will omit the corresponding inequality constraint.

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: angle[1] - angle[2] ∈ [-1.7, 1.7]
Branch 2: angle[1] - angle[3] ∈ [-3.1, 3.1]
Branch 3: angle[2] - angle[3] ∈ [-3.1, 3.1]

Active Power Flow Constraints

The flow field refers to the inequality constraints linked to active power flow limits at the from-bus and to-bus ends of each branch. These limits are set using the longTerm keyword in the addBranch! function:

julia> print(system.branch.label, analysis.method.constraint.flow.active)Branch 1: angle[1] - angle[2] ∈ [-0.006, 0.006]
Branch 2: angle[1] - angle[3] ∈ [-0.0012, 0.0012]
Branch 3: angle[2] - angle[3] ∈ [-0.0012, 0.0012]

Info

Please note that if the limit constraints are set to longTerm = 0.0 for the corresponding branch, JuliGrid will omit the corresponding inequality constraint.

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, longTerm = 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: angle[1] - angle[3] ∈ [-0.0042, 0.0042]
Branch 3: angle[2] - angle[3] ∈ [-0.0012, 0.0012]

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: active[1] ∈ [0, 0.8]
Generator 2: active[2] = 0
Generator 3: active[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. Consequently, we can adjust these specific constraints using the updateGenerator! function, for example:

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: active[1] ∈ [0, 0.8]
Generator 2: active[2] ∈ [0, 0.5]
Generator 3: active[3] ∈ [0, 0.4]

It is important to note that by bringing back Generator 2 into service, it 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 active[3] - actwise[3] ≤ 3.805263157894739
Generator 3: 1.2 active[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: angle[1] - angle[3] ∈ [-0.0042, 0.0042]
Branch 3: angle[2] - angle[3] ∈ [-0.0012, 0.0012]
Branch 4: angle[1] - angle[2] ∈ [-0.12, 0.12]
julia> print(system.generator.label, analysis.method.constraint.capability.active)Generator 1: active[1] ∈ [0, 0.8]
Generator 2: active[2] ∈ [0, 0.5]
Generator 3: active[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])

Info

In the event that a user deletes a constraint and subsequently executes a function that updates bus, branch, or generator parameters, and if the deleted constraint is affected by these functions, JuliaGrid will automatically reinstate that constraint. Users should exercise caution when deleting constraints, as this action is considered potentially harmful since it operates independently of power system data.

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 active[1]² + 500 active[1] + actwise[3] + active[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 active[1]² + 853.4 active[3]² + 500 active[1] + active[2] + 257 active[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 active[1]² + 853.4 active[3]² + 500 active[1] + active[2] + 257 active[3]

Setup Starting Primal Values

In JuliaGrid, the assignment of starting primal values for optimization variables takes place when the solve! function is executed. 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, value) in system.generator.label
    analysis.power.generator.active[value] = generatorPower(system, flow; label = key)
end

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

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.09539999999999904
Generator 3: 0.0
Generator 4: 0.14460000000000023
julia> print(system.bus.label, analysis.voltage.angle)Bus 1: -0.1
Bus 2: -0.10460000000000001
Bus 3: -0.1042

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.09539999999999904

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 values for the subsequent solving step correspond to the solution obtained from the previous step.

In the previous example, we obtained the following solution:

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

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

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

Next, solve this new power system. During the execution of the solve! function, the primal starting values will first be set, and these values will be defined according to the values given above.

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.08
Generator 3: 0.015399999999999636
Generator 4: 0.14460000000000023
julia> print(system.bus.label, analysis.voltage.angle)Bus 1: -0.1
Bus 2: -0.10460000000000001
Bus 3: -0.1042

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 composite type, employing the startingPrimal! function:

startingPrimal!(system, analysis)

These values are precisely identical to what we would obtain if we executed the dcOptimalPowerFlow function following all the updates we performed:

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

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, longTerm = 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.09000000000000127
Bus 3: -0.05
julia> print(system.branch.label, analysis.power.from.active)Branch 1: -0.007142857142858339
Branch 2: 0.007142857142855563
Branch 3: 0.04285714285714726

Info

To better understand the powers associated with buses and branches that are calculated by the power! function, we suggest referring to the tutorials on DC Optimal Power Flow.

To compute specific quantities for particular components, rather than calculating powers or currents for all components, users can utilize one of the provided functions below.

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.09000000000000127

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.19000000000000128

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")0.007142857142855563
julia> active = toPower(system, analysis; label = "Branch 2")-0.007142857142855563