AC Optimal Power Flow

JuliaGrid uses the JuMP package to construct optimal power flow models, allowing users to manipulate these models using the standard functions provided by JuMP. As a result, JuliaGrid supports popular solvers mentioned in the JuMP documentation to solve the optimization problem.

Performing AC optimal power flow first requires the PowerSystem type that has been created with the AC model. After that, create the AcOptimalPowerFlow type to establish the AC optimal power flow framework using the function:

• acOptimalPowerFlow.

To solve the AC optimal power flow problem and obtain generator active and reactive power outputs, as well as bus voltage magnitudes and angles, users can use the following function:

• solve!.

After solving the AC optimal power flow, JuliaGrid provides functions for computing powers and currents:

• power!,
• current!.

Alternatively, instead of using functions responsible for solving optimal power flow and computing powers and currents, users can use the wrapper function:

• powerFlow!.

Users can also access specialized functions for computing specific types of powers or currents for individual buses, branches, or generators within the power system.

Optimal Power Flow Model

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

using JuMP, Ipopt

system = powerSystem()

@bus(minMagnitude = 0.95, maxMagnitude = 1.05)
addBus!(system; label = "Bus 1", type = 3, active = 0.1, angle = -0.1)
addBus!(system; label = "Bus 2", reactive = 0.01, magnitude = 1.1)

@branch(minDiffAngle = -pi, maxDiffAngle = pi, reactance = 0.5, type = 1)
addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 2", maxFromBus = 0.15)

@generator(maxActive = 0.5, minReactive = -0.1, maxReactive = 0.1, status = 0)
addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 0.4, reactive = 0.2)
addGenerator!(system; label = "Generator 2", bus = "Bus 2", active = 0.2, reactive = 0.1)

cost!(system; generator = "Generator 1", active = 2, polynomial = [800.0; 200.0; 80.0])
cost!(system; generator = "Generator 2", active = 1, piecewise = [10 12.3; 14.7 16.8; 18 19])

cost!(system; generator = "Generator 1", reactive = 2, polynomial = [2.0])
cost!(system; generator = "Generator 2", reactive = 1, piecewise = [2.0 4.0; 6.0 8.0])

acModel!(system)

Next, the acOptimalPowerFlow function is used to formulate the AC optimal power flow problem:

analysis = acOptimalPowerFlow(system, Ipopt.Optimizer)

analysis = acOptimalPowerFlow(system, Ipopt.Optimizer; interval = false)

Optimization Variables

In the AC optimal power flow model, the active and reactive power outputs of the generators are expressed as nonlinear functions of the bus voltage magnitudes and angles. As a result, the variables in this model include the active and reactive power outputs of the generators, as well as the bus voltage magnitudes and angles:

julia> JuMP.all_variables(analysis.method.jump)8-element Vector{VariableRef}:
 magnitude[1]
 magnitude[2]
 angle[1]
 angle[2]
 active[1]
 active[2]
 reactive[1]
 reactive[2]

It is important to note that this is not a comprehensive set of optimization variables. When the cost function is defined as a piecewise linear function comprising multiple segments, as illustrated in the case of the active power output cost for Generator 2, JuliaGrid automatically generates helper optimization variables named actwise or reactwise, and formulates a set of linear constraints to effectively address these cost functions. For simplicity, this section initially assumes that Generator 2 is out-of-service. Consequently, the helper variable is not included in the set of optimization variables. Later sections activate the generator, introducing the helper variable and additional constraints to the optimization model.

It is worth emphasizing that in instances where a piecewise linear cost function consists of only a single segment, as demonstrated by the reactive power output cost of Generator 2, the function is modeled as a standard linear function, avoiding the need for additional helper optimization variables.

Please be aware that JuliaGrid maintains references to all variables, which are categorized into six fields:

julia> fieldnames(typeof(analysis.method.variable.voltage))(:magnitude, :angle)
julia> fieldnames(typeof(analysis.method.variable.power))(:active, :reactive, :actwise, :reactwise)

Variable Names

Users can define custom variable names for printing equations, which can help present them in a more compact form. For example:

analysis = acOptimalPowerFlow(system, Ipopt.Optimizer; magnitude = "V", angle = "θ")

Add Variables

Once the AcOptimalPowerFlow type is established, users can add new variables representing generator active power outputs by introducing additional generators. For example:

addGenerator!(analysis; label = "Generator 3", bus = "Bus 1", maxActive = 0.2, status = 1)

This command adds both new variables and the corresponding box constraints to the optimization model.

To confirm that the variable has been successfully added, you can use the following function:

julia> JuMP.is_valid(analysis.method.jump, analysis.method.variable.power.active[3])true
julia> JuMP.is_valid(analysis.method.jump, analysis.method.variable.power.reactive[3])true

Constraint Functions

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

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

They fall into two main categories: box constraints and non-box constraints.

Box Constraints

The slack constraint is represented as an equality tied to the fixed voltage angle at the slack bus.

The capability constraints define variable bounds on both active and reactive power generation and are always implemented as two separate constraints: one for the lower bound and one for the upper bound. If the bounds are equal, or if the generator is out-of-service, JuliaGrid models the constraint as an equality instead.

The magnitude field of the voltage constraints defines bounds on bus voltage magnitude values and is also implemented using a lower and an upper bound. If the bounds are equal, JuliaGrid models the constraint as an equality instead.

Non-Box Constraints

The balance constraints correspond to the active and reactive power balance equations defined at each bus and are modeled as equality constraints.

The angle field of the voltage constraints is associated with the minimum and maximum voltage angle difference between the from-bus and to-bus ends of each branch and is modeled as an interval constraint by default. If the bounds are equal, an equality constraint is used instead.

The flow constraints, which refer to branch flow limits at both ends of each branch, are also modeled as interval constraints by default. If the bounds are equal, an equality constraint is used.

If preferred, both the angle field of the voltage constraints and the flow constraints can be represented as two separate one-sided constraints, one for the lower and one for the upper bound, by setting the keyword argument interval = false when calling the acOptimalPowerFlow function.

Finally, the piecewise constraints are introduced when piecewise linear cost functions with multiple segments are defined, and they impose only upper bounds.

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

Users can modify this constraint by changing which bus serves as the slack bus and by adjusting the value of the bus angle. This uses the updateBus! function, for example:

updateBus!(analysis; label = "Bus 1", type = 1)
updateBus!(analysis; label = "Bus 2", type = 3, angle = -0.2)

The updated slack constraint can be inspected:

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

Bus Power Balance Constraints

The balance field contains references to the equality constraints associated with the active and reactive power balance equations defined for each bus. These constraints ensure that the total active and reactive power injected by the generators matches the total active and reactive power demanded at each bus.

The constant term in the active power balance equations is determined by the active keyword within the addBus! function, which defines the active power demanded at the bus. Access the references to the active power balance constraints using the following code snippet:

julia> print(system.bus.label, analysis.method.constraint.balance.active)Bus 1: active[3] + (V[1] * ((V[2] * ((0 * cos(θ[1] - θ[2])) - (2 * sin(θ[1] - θ[2])))))) = 0.1
Bus 2: V[2] * ((V[1] * ((0 * cos(-θ[1] + θ[2])) - (2 * sin(-θ[1] + θ[2]))))) = 0

Similarly, the constant term in the reactive power balance equations is determined by the reactive keyword within the addBus! function, which defines the reactive power demanded at the bus. Access the references to the reactive power balance constraints using the following code snippet:

julia> print(system.bus.label, analysis.method.constraint.balance.reactive)Bus 1: reactive[3] + (V[1] * ((-2 V[1]) + (V[2] * ((0 * sin(θ[1] - θ[2])) + (2 * cos(θ[1] - θ[2])))))) = 0
Bus 2: V[2] * ((-2 V[2]) + (V[1] * ((0 * sin(-θ[1] + θ[2])) + (2 * cos(-θ[1] + θ[2]))))) = 0.01

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!(analysis; label = "Bus 2", active = 0.5)
updateBranch!(analysis; label = "Branch 1", reactance = 0.25)

The updated set of active power balance constraints can be examined:

julia> print(system.bus.label, analysis.method.constraint.balance.active)Bus 1: active[3] + (V[1] * ((V[2] * ((0 * cos(θ[1] - θ[2])) - (4 * sin(θ[1] - θ[2])))))) = 0.1
Bus 2: V[2] * ((V[1] * ((0 * cos(-θ[1] + θ[2])) - (4 * sin(-θ[1] + θ[2]))))) = 0.5

Bus Voltage Constraints

The voltage field contains references to the inequality constraints associated with the voltage magnitude and voltage angle difference limits. These constraints ensure that the bus voltage magnitudes and the angle differences between the from-bus and to-bus ends of each branch are within specified limits.

The minimum and maximum bus voltage magnitude limits are set using the minMagnitude and maxMagnitude keywords within the addBus! function. The constraints associated with these limits can be accessed using:

julia> print(system.bus.label, analysis.method.constraint.voltage.magnitude)Bus 1: V[1] ≤ 1.05, V[1] ≥ 0.95
Bus 2: V[2] ≤ 1.05, V[2] ≥ 0.95

The minimum and maximum voltage angle difference limits between the from-bus and to-bus ends of each branch are set using the minDiffAngle and maxDiffAngle keywords within the addBranch! function. The constraints associated with these limits can be accessed using the following code snippet:

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

Using the updateBus! and updateBranch! functions, users can modify these constraints:

updateBus!(analysis; label = "Bus 1", minMagnitude = 1.0, maxMagnitude = 1.0)
updateBranch!(analysis; label = "Branch 1", minDiffAngle = -1.7, maxDiffAngle = 1.7)

The updated set of constraints can be examined:

julia> print(system.bus.label, analysis.method.constraint.voltage.magnitude)Bus 1: V[1] = 1
Bus 2: V[2] ≤ 1.05, V[2] ≥ 0.95
julia> print(system.branch.label, analysis.method.constraint.voltage.angle)Branch 1: θ[1] - θ[2] ∈ [-1.7, 1.7]

Branch Flow Constraints

The flow field refers to inequality constraints that enforce limits on the apparent power flow, active power flow, or current flow magnitude at the from-bus and to-bus ends of each branch. The type of constraint applied is specified using the type keyword in the addBranch! function:

• type = 1 active power flow,
• type = 2 apparent power flow,
• type = 3 apparent power flow with a squared inequality constraint,
• type = 4 current flow magnitude,
• type = 5 current flow magnitude with a squared inequality constraint.

These limits are specified using the minFromBus, maxFromBus, minToBus, and maxToBus keywords within the addBranch! function. By default, these limit keywords are associated with apparent power (type = 3).

However, the example configures it to use active power flow by setting type = 1. To access the flow constraints of branches at the from-bus end, use the following code snippet:

julia> print(system.branch.label, analysis.method.constraint.flow.from)Branch 1:  - ((V[1]*V[2]) * ((0 * cos(θ[1] - θ[2])) + (-4 * sin(θ[1] - θ[2])))) ∈ [0, 0.15]

Additionally, using the updateBranch! function, modify these specific constraints:

updateBranch!(analysis; label = "Branch 1", minFromBus = -0.15, maxToBus = 0.15)

The updated set of flow constraints can be examined:

julia> print(system.branch.label, analysis.method.constraint.flow.from)Branch 1:  - ((V[1]*V[2]) * ((0 * cos(θ[1] - θ[2])) + (-4 * sin(θ[1] - θ[2])))) ∈ [-0.15, 0.15]
julia> print(system.branch.label, analysis.method.constraint.flow.to)Branch 1:  - ((V[1]*V[2]) * ((0 * cos(θ[1] - θ[2])) - (-4 * sin(θ[1] - θ[2])))) ∈ [0, 0.15]

Generator Power Capability Constraints

The capability field contains references to the inequality constraints associated with the minimum and maximum active and reactive power outputs of the generators.

The constraints associated with the minimum and maximum active power output limits of the generators are defined using the minActive and maxActive keywords within the addGenerator! function. To access the constraints associated with these limits, use the following code snippet:

julia> print(system.generator.label, analysis.method.constraint.capability.active)Generator 1: active[1] = 0
Generator 2: active[2] = 0
Generator 3: active[3] ≤ 0.2, active[3] ≥ 0

Similarly, the constraints associated with the minimum and maximum reactive power output limits of the generators are specified using the minReactive and maxReactive keywords within the addGenerator! function. To access these constraints, use the following code snippet:

julia> print(system.generator.label, analysis.method.constraint.capability.reactive)Generator 1: reactive[1] = 0
Generator 2: reactive[2] = 0
Generator 3: reactive[3] ≤ 0.1, reactive[3] ≥ -0.1

As demonstrated, the active and reactive power outputs of Generator 1 and Generator 2 are currently fixed at zero due to previous actions that set these generators out-of-service. However, modify these specific constraints using the updateGenerator! function, as shown below:

updateGenerator!(analysis; label = "Generator 1", status = 1)
updateGenerator!(analysis; label = "Generator 2", status = 1, minActive = 0.1)

The updated set of constraints can be examined:

julia> print(system.generator.label, analysis.method.constraint.capability.active)Generator 1: active[1] ≤ 0.5, active[1] ≥ 0
Generator 2: active[2] ≤ 0.5, active[2] ≥ 0.1
Generator 3: active[3] ≤ 0.2, active[3] ≥ 0
julia> print(system.generator.label, analysis.method.constraint.capability.reactive)Generator 1: reactive[1] ≤ 0.1, reactive[1] ≥ -0.1
Generator 2: reactive[2] ≤ 0.1, reactive[2] ≥ -0.1
Generator 3: reactive[3] ≤ 0.1, reactive[3] ≥ -0.1

Power Piecewise Constraints

In cost modeling, the piecewise field serves as a reference to the inequality constraints associated with piecewise linear cost functions. These constraints are defined using the cost! function with active = 1 or reactive = 1.

In this example, only the active power cost of Generator 2 is modeled as a piecewise linear function with two segments, and JuliaGrid takes care of setting up the appropriate inequality constraints for each segment:

julia> print(system.generator.label, analysis.method.constraint.piecewise.active)Generator 2: 0.9574468085106385 active[2] - actwise[2] ≤ -2.7255319148936152
Generator 2: 0.6666666666666663 active[2] - actwise[2] ≤ -7.000000000000007

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

As mentioned at the beginning, piecewise linear cost functions with multiple segments will also introduce helper variables that are added to the objective function. In this specific example, the helper variable is:

julia> analysis.method.variable.power.actwise[2]actwise[2]

Add Constraints

Users can introduce additional branch-related constraints into the AC optimal power flow model by adding a new branch with the addBranch! function. For example, a user can include a new branch in an already defined PowerSystem and AcOptimalPowerFlow model:

addBranch!(analysis; label = "Branch 2", from = "Bus 1", to = "Bus 2", reactance = 1)

This will affect all constraints related to branches, but it will also update balance constraints to configure the optimization model to match the current state of the power system. For example, the following updated constraints can be observed:

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

Similarly, the addGenerator! function adds both new variables and their associated box constraints.

Delete Constraints

When a branch or generator is taken out-of-service, JuliaGrid automatically adjusts the optimization problem to reflect that action, which may include removing certain constraints.

In some cases, users may also remove specific constraints manually. This can be done using the remove! function by specifying the constraint type: :slack, :capability, :balance, :voltage, :flow, or :piecewise.

For constraint types such as :capability, :balance, and :piecewise, users must also specify whether the constraint targets :active or :reactive power. Similarly, for :voltage, the options are :magnitude or :angle, and for :flow, the options are :from or :to.

For example, to delete the constraint associated with the voltage angle difference at Branch 2, use:

remove!(analysis, :voltage, :angle; label = "Branch 2")

Alternatively, instead of using a label, constraints can also be deleted by index:

remove!(analysis, :voltage, :angle; index = 4)

After these operations, the remaining voltage angle difference constraints can be displayed:

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

Objective Function

The objective function of the AC optimal power flow is formulated using polynomial and piecewise linear cost functions associated with the generators, defined using the cost! functions.

In the provided example, the objective function minimized to obtain the optimal values is:

julia> JuMP.objective_function(analysis.method.jump)800 active[1]² + 200 active[1] + actwise[2] + reactive[2] + 84

The objective function is stored in analysis.method.objective, where it is organized to separate its quadratic and nonlinear components.

Update Objective Function

With the cost! functions, users can modify the objective function by adjusting polynomial or piecewise linear coefficients or by changing the type of polynomial or piecewise linear function used. For example, consider Generator 1, which uses a quadratic polynomial cost function for active power. Redefine the cost function for this generator as a cubic polynomial and define a nonlinear objective function:

cost!(analysis; generator = "Generator 1", active = 2, polynomial = [63; 25; 4; 0.5])

This leads to an updated objective function, which can be examined:

julia> JuMP.objective_function(analysis.method.jump)(25 active[1]² + actwise[2] + reactive[2] + 4 active[1] + 4.5) + (63 * (active[1] ^ 3))

Setup Initial Values

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

Initial Primal Values

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

julia> generator = analysis.power.generator;
julia> print(system.generator.label, generator.active, generator.reactive)Generator 1: 0.4, 0.2
Generator 2: 0.2, 0.1
Generator 3: 0.0, 0.0
julia> print(system.bus.label, analysis.voltage.magnitude, analysis.voltage.angle)Bus 1: 1.0, -0.1
Bus 2: 1.1, 0.0

Users can adjust these values, which will then be used as the initial primal values when executing the solve! function.

Using AC Power Flow

In this perspective, users can run the AC power flow analysis and use the resulting solution to set initial primal values. For example:

flow = newtonRaphson(system)
powerFlow!(flow; power = true)

After obtaining the solution, use the active and reactive power outputs of the generators, along with bus voltage magnitudes and angles, to set the initial values:

setInitialPoint!(analysis, flow)

Initial 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. By default, dual values are undefined, but users can manually assign them using the addDual! function.

If a constraint is defined as an equality, an interval, or has only a lower or upper bound, it corresponds to a single dual variable. In such cases, an initial value can be set using the dual keyword. For example:

addDual!(analysis, :balance, :active; label = "Bus 1", dual = 1e-3)

For constraints with both lower and upper bounds, users can assign initial dual values separately using the lower and upper keywords. For example:

addDual!(analysis, :capability, :reactive; label = "Generator 1", lower = 500.0, upper = 0.0)

Alternatively, dual variables can be added by specifying the constraint index instead of a label:

addDual!(analysis, :capability, :reactive; index = 1, lower = 500.0, upper = 0.0)

Optimal Power Flow Solution

To establish the AC optimal power flow problem, use the acOptimalPowerFlow function. After setting up the problem, use the solve! function to compute the optimal values for the active and reactive power outputs of the generators and the bus voltage magnitudes and angles:

solve!(analysis)

Executing this function produces the solution with the optimal values for the active and reactive power outputs of the generators, as well as the bus voltage magnitudes and angles.

julia> generator = analysis.power.generator;
julia> print(system.generator.label, generator.active, generator.reactive)Generator 1: -9.373524037461336e-9, 0.05626636529233223
Generator 2: 0.49999999019365665, -0.10000000749409634
Generator 3: 0.10000001917986744, 0.05626636529233222
julia> print(system.bus.label, analysis.voltage.magnitude, analysis.voltage.angle)Bus 1: 1.0, -0.1999999979935736
Bus 2: 0.9774934538830672, -0.2

Objective Value

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

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

Dual Variables

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

julia> print(system.bus.label, analysis.method.dual.balance.active)Bus 1: 9.61257632143584e-15
Bus 2: 0.8096015952603774

Wrapper Function

JuliaGrid provides a wrapper function for AC optimal power flow analysis and also supports the computation of powers and currents using the powerFlow! function:

analysis = acOptimalPowerFlow(system, Ipopt.Optimizer)
powerFlow!(analysis; verbose = 1)

EXIT: The optimal solution was found.

Print Results in the REPL

Users can use 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:

show = Dict("Active Power Balance" => false)
printBusConstraint(analysis; show)

|---------------------------------------------------------------------------|
| Bus Constraint Data                                                       |
|---------------------------------------------------------------------------|
| Label |            Voltage Magnitude             | Reactive Power Balance |
|       |                                          |                        |
|       | Minimum | Solution | Maximum |      Dual |   Solution |      Dual |
|       |    [pu] |     [pu] |    [pu] | [$/pu-hr] |       [pu] | [$/pu-hr] |
|-------|---------|----------|---------|-----------|------------|-----------|
| Bus 1 |  1.0000 |   1.0000 |  1.0000 |   -0.0000 |    -0.0000 |    0.0000 |
| Bus 2 |  0.9500 |   0.9775 |  1.0500 |    0.0000 |     0.0000 |    0.0000 |
|---------------------------------------------------------------------------|

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

printBusConstraint(analysis; label = "Bus 1", header = true)
printBusConstraint(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:

open("bus.txt", "w") do file
    printBusConstraint(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(analysis, io; style = false)
CSV.write("constraint.csv", CSV.File(take!(io); delim = "|"))

Primal and Dual Warm Start

Solving the same AcOptimalPowerFlow type again provides a warm start. Here, the initial primal and dual values for the next solve correspond to the solution obtained from the previous step, including any user-defined data previously integrated in JuliaGrid.

Primal Variables

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

julia> generator = analysis.power.generator;
julia> print(system.generator.label, generator.active, generator.reactive)Generator 1: -9.373593224513171e-9, 0.05626636529232862
Generator 2: 0.49999999019699026, -0.10000000749409634
Generator 3: 0.10000001917660305, 0.05626636529233582
julia> print(system.bus.label, analysis.voltage.magnitude, analysis.voltage.angle)Bus 1: 1.0, -0.1999999979942667
Bus 2: 0.9774934538830672, -0.2

Dual Variables

All dual values are also available. Here, only the dual variables for one type of constraint are listed as an example:

julia> print(system.generator.label, analysis.method.dual.capability.reactive)Generator 1: -4.126311435827629e-8
Generator 1: 0.0
Generator 2: 0.0
Generator 2: 0.9999999449411637
Generator 3: -4.126311435828731e-8
Generator 3: 0.0

Modify Optimal Power Flow

Now introduce changes to the power system from the previous example:

updateGenerator!(analysis; label = "Generator 3", maxActive = 0.05)

Next, solve this modified optimal power flow problem. If solve! is used at this point, the primal and dual initial values will be set to the previously obtained values:

powerFlow!(analysis; verbose = 1)

EXIT: The optimal solution was found.

This produces a new solution:

julia> generator = analysis.power.generator;
julia> print(system.generator.label, generator.active, generator.reactive)Generator 1: 0.04999998074011849, 0.056266365292332234
Generator 2: 0.5000000096192797, -0.10000000749409634
Generator 3: 0.05000000964060179, 0.056266365292332234
julia> print(system.bus.label, analysis.voltage.magnitude, analysis.voltage.angle)Bus 1: 1.0, -0.20000000196815226
Bus 2: 0.9774934538830672, -0.2

Reset Primal and Dual Values

Users can reset initial primal and dual values to their default configurations at any time. This can be done using the active and reactive power outputs of the generators and the initial bus voltage magnitudes and angles extracted from the PowerSystem type, with the setInitialPoint! function:

setInitialPoint!(analysis)

The primal initial values will now be identical to those that would be obtained if the acOptimalPowerFlow function were executed after all the updates have been applied, while all dual variable values will be removed.

Extended Formulation

The JuMP model created by JuliaGrid is stored in the method.jump field of the AcOptimalPowerFlow type. This allows users to modify the model directly using JuMP macros and functions as needed. However, when making such modifications, users become responsible for tasks like setting initial values and extracting solutions, since these changes operate outside the standard JuliaGrid workflow.

Beyond this approach, JuliaGrid also provides a way to extend the standard AC optimal power flow formulation within its own framework. This supports features such as warm start and automatic solution storage, as described below.

Add Variable

User-defined variables can be added to the AC optimal power flow model using the @addVariable macro. It also allows immediate assignment of initial primal and dual values. For example:

@addVariable(analysis, 0.0 <= y <= 0.2, primal = 0.1, lower = 10.0, upper = 0.0)

Collections of variables can also be defined:

@addVariable(analysis, 0.0 <= x[i = 1:2] <= 0.4, primal = [0.1, 0.2], upper = [0.0; -2.5])

Add Constraints

Custom constraints can be added to the AC optimal power flow model using the @addConstraint macro. These constraints are not limited to user-defined variables; any optimization variable defined up to that point can be used. Focus on the voltage angle variables:

θ = analysis.method.variable.voltage.angle

Next, a new constraint can be defined, and at the same time, an initial dual value can be specified:

@addConstraint(analysis, 0.1 <= x[1] + 2 * x[2] + y + θ[2] <= 1.2, dual = 0.0)

Collections of constraints can also be defined:

@addConstraint(analysis, [i = 1:2], x[i] + 2 * θ[i] <= 0.6, dual = [0.0; 0.5])

Delete Constraints

To remove a constraint, use the remove! function with the :constraint symbol. For example, to remove the first added constraint:

remove!(analysis, :constraint; index = 1)

Objective Function

Users can modify the objective function using the set_objective_function function from the JuMP package. In JuliaGrid, the original objective is stored in analysis.method.objective, which can be accessed and customized as needed. This makes it possible to simultaneously remove nonlinear components and adjust the quadratic part of the objective function:

expr = 50 * x[1] - x[2]^2 + y + 123
JuMP.set_objective_function(analysis.method.jump, analysis.method.objective.quadratic - expr)

The updated objective function is:

julia> JuMP.objective_function(analysis.method.jump)25 active[1]² + x[2]² + 4 active[1] + actwise[2] + reactive[2] - 50 x[1] - y - 118.5

Optimal Power Flow Solution

Users can now solve the extended formulation using:

powerFlow!(analysis; verbose = 1)

EXIT: The optimal solution was found.

After solving, users can access the optimal values:

julia> analysis.power.generator.active3-element Vector{Float64}:
 0.04999998029613146
 0.5000000098437287
 0.05000000986013984
julia> analysis.power.generator.reactive3-element Vector{Float64}:
  0.05626636612837133
 -0.10000000909090911
  0.05626636612837133
julia> analysis.voltage.magnitude2-element Vector{Float64}:
 1.0
 0.9774934535486515
julia> analysis.voltage.angle2-element Vector{Float64}:
 -0.20000000201407564
 -0.2
julia> analysis.extended.solution3-element Vector{Float64}:
 0.2000000090909091
 0.40000000998181817
 4.393914459490238e-5

Power and Current Analysis

After obtaining the solution from the AC optimal power flow, calculate various electrical quantities related to buses and branches using the power! and current! functions. For instance, consider the power system used to obtain the AC optimal power flow solution:

using Ipopt

system = powerSystem()

@bus(minMagnitude = 0.9, maxMagnitude = 1.1)
addBus!(system; label = "Bus 1", type = 3, magnitude = 1.05, angle = 0.17)
addBus!(system; label = "Bus 2", active = 0.1, reactive = 0.01, conductance = 0.04)
addBus!(system; label = "Bus 3", active = 0.05, reactive = 0.02)

@branch(resistance = 0.5, reactance = 1.0, conductance = 1e-4, susceptance = 0.01)
addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 2", maxFromBus = 0.15)
addBranch!(system; label = "Branch 2", from = "Bus 1", to = "Bus 3", maxFromBus = 0.10)
addBranch!(system; label = "Branch 3", from = "Bus 2", to = "Bus 3", maxFromBus = 0.25)

@generator(maxActive = 0.5, minReactive = -0.1, maxReactive = 0.1)
addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 3.2, reactive = 0.5)
addGenerator!(system; label = "Generator 2", bus = "Bus 2", active = 0.2, reactive = 0.1)

cost!(system; generator = "Generator 1", active = 2, polynomial = [1100.2; 500; 80])
cost!(system; generator = "Generator 2", active = 1, piecewise = [10 12.3; 14.7 16.8; 18 19])

analysis = acOptimalPowerFlow(system, Ipopt.Optimizer)
powerFlow!(analysis)

Now use the following functions to calculate powers and currents:

power!(analysis)
current!(analysis)

For instance, to show the active power injections and the from-bus current magnitudes, use:

julia> print(system.bus.label, analysis.power.injection.active)Bus 1: -9.583329214392367e-9
Bus 2: 0.08516587136754715
Bus 3: -0.04999999997621973
julia> print(system.branch.label, analysis.current.from.magnitude)Branch 1: 0.026583569097689393
Branch 2: 0.02328034568030229
Branch 3: 0.03814096829046926

Print Results in the REPL

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

@voltage(pu, deg)
@power(MW, MVAr)
show = Dict("Power Generation" => false, "Current Injection" => false)
printBusData(analysis; show)

|-------------------------------------------------------------------------------------------|
| Bus Data                                                                                  |
|-------------------------------------------------------------------------------------------|
| Label |       Voltage       |    Power Demand    |  Power Injection   |    Shunt Power    |
|       |                     |                    |                    |                   |
|   Bus | Magnitude |   Angle |  Active | Reactive |  Active | Reactive | Active | Reactive |
|       |      [pu] |   [deg] |    [MW] |   [MVAr] |    [MW] |   [MVAr] |   [MW] |   [MVAr] |
|-------|-----------|---------|---------|----------|---------|----------|--------|----------|
| Bus 1 |    0.9279 |  9.7403 |  0.0000 |   0.0000 | -0.0000 |   3.0790 | 0.0000 |  -0.0000 |
| Bus 2 |    0.9140 | 11.6476 | 10.0000 |   1.0000 |  8.5166 |  -3.2859 | 3.3417 |  -0.0000 |
| Bus 3 |    0.9000 |  9.1611 |  5.0000 |   2.0000 | -5.0000 |  -2.0000 | 0.0000 |  -0.0000 |
|-------------------------------------------------------------------------------------------|

Active and Reactive Power Injection

To calculate the active and reactive power injection associated with a specific bus, use:

julia> active, reactive = injectionPower(analysis; label = "Bus 1")(-9.583329214392367e-9, 0.030789711880008647)

Active and Reactive Power Injection from Generators

To calculate the active and reactive power injection from generators at a specific bus, use:

julia> active, reactive = supplyPower(analysis; label = "Bus 2")(0.1851658713222783, -0.022858875148501698)

Active and Reactive Power at Shunt Element

To calculate the active and reactive power associated with the shunt element at a specific bus, use:

julia> active, reactive = shuntPower(analysis; label = "Bus 2")(0.033417219109870354, -0.0)

Active and Reactive Power Flow

Similarly, to compute the active and reactive power flow at both the from-bus and to-bus ends of a specific branch, use:

julia> active, reactive = fromPower(analysis; label = "Branch 2")(0.01718399882598644, 0.013091667390527213)
julia> active, reactive = toPower(analysis; label = "Branch 2")(-0.016754093702898584, -0.02075438599122779)

Active and Reactive Power at Charging Admittances

To calculate the total active and reactive power associated with the branch charging admittances of a particular branch, use:

julia> active, reactive = chargingPower(analysis; label = "Branch 1")(8.482535529692716e-5, -0.008482535529692716)

Active powers indicate active losses within the branch's charging admittances. Moreover, charging admittances inject reactive power into the power system due to their capacitive nature, as denoted by a negative sign.

Active and Reactive Power at Series Impedance

To calculate the active and reactive power across the series impedance of the branch, use:

julia> active, reactive = seriesPower(analysis; label = "Branch 2")(0.0003463509187067052, 0.0006927018374134103)

The active power also considers active losses originating from the series resistance of the branch, while the reactive power represents reactive losses resulting from the impedance's inductive characteristics.

Current Injection

To calculate the current injection associated with a specific bus, use:

julia> magnitude, angle = injectionCurrent(analysis; label = "Bus 1")(0.033180637369821964, -1.4007966380459163)

Current Flow

Compute the current flow at both the from-bus and to-bus ends of the specific branch by using:

julia> magnitude, angle = fromCurrent(analysis; label = "Branch 2")(0.02328034568030229, -0.4810434670251617)
julia> magnitude, angle = toCurrent(analysis; label = "Branch 2")(0.029636431574178668, 2.4098367486797905)

Current Through Series Impedance

To calculate the current passing through the series impedance of the branch in the direction from the from-bus end to the to-bus end, use the following function:

julia> magnitude, angle = seriesCurrent(analysis; label = "Branch 2")(0.026319229422865145, -0.6228132099043933)