Skip to content

Commit

Permalink
Merge branch 'main' into jit-warmup-case
Browse files Browse the repository at this point in the history
  • Loading branch information
ccoffrin authored Jan 27, 2024
2 parents 079eaeb + 2376b01 commit 1ff0e10
Show file tree
Hide file tree
Showing 13 changed files with 945 additions and 432 deletions.
33 changes: 33 additions & 0 deletions .github/workflows/variants.yml
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
name: variants
on:
push:
branches:
- main
pull_request:
types: [opened, synchronize, reopened]
jobs:
test:
name: ${{ matrix.variant }} - ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
runs-on: ${{ matrix.os }}
strategy:
fail-fast: false
matrix:
variant: ['jump-nl', 'jump-symbolic-ad', 'jump-traced', 'optimization-cs-asrd']
version: ['1']
os: [ubuntu-latest]
arch: [x64]
steps:
- uses: actions/checkout@v3
- uses: julia-actions/setup-julia@v1
with:
version: ${{ matrix.version }}
arch: ${{ matrix.arch }}
- name: install
shell: julia --color=yes --project=variants {0}
run: |
using Pkg
Pkg.instantiate()
- name: test
env:
MODELING_FRAMEWORK: ${{ matrix.variant }}
run: julia --color=yes --project=variants variants/runtests.jl
12 changes: 12 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

The AC Optimal Power Flow problem (AC-OPF) is one of the most foundational optimization problems that arises in the design and operations of power networks.
Mathematically the AC-OPF is a large-scale, sparse, non-convex nonlinear continuous optimization problem.
The models make heavy use of scalar indexing in the definition of the nonlinear functions, and most operations are simple floating point operations (`+`, `*`, `cos`, and `sin`) applied to scalar data.
In practice AC-OPF is most often solved to local optimality conditions using interior point methods.
This project proposes AC-OPF as _proxy-application_ for testing the viability of different nonlinear optimization frameworks, as performant solutions to AC-OPF has proven to be a necessary (but not always sufficient) condition for solving a wide range of industrial network optimization tasks.

Expand All @@ -22,6 +23,17 @@ At present this project is focused on comparing NLP modeling layers that are ava

This work is not intended for comparing different nonlinear optimization algorithms, which are often independent of the NLP modeling layer. Consequently, the [Ipopt](https://github.com/jump-dev/Ipopt.jl) solver is used as a standard NLP algorithm whenever it is accessible from the modeling layer.

### Comparison with SciMLBenchmarks

All framework implementations use PowerModels.jl to read input `.m` data files into non-concrete dictionaries, from which the objective and constraint functions are defined.
Reading the data into non-concrete dictionaries simulates a very common non-expert user of NLP tooling, but biases the performance results in favor of symbolic front-ends which perform additional optimizations on the objective and constraint functions before performing the optimization.

Moreover, the `total time` reported by `rosetta-opf` includes Julia's compilation time.
This biases performance against frameworks which generate large amounts of new code that requires compilation.
We do not attempt to remove the compilation overhead, for example, by using a tool like [PackageCompiler.jl](https://github.com/JuliaLang/PackageCompiler.jl).

For a variation of `rosetta-opf` that removes these factors to focus on core implementation differences of the numerical backends, see the [SciMLBenchmarks adaptation](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/OptimizationFrameworks/optimal_powerflow/).

## Mathematical and Data Models
This work adopts the mathematical model and data format that is used in the IEEE PES benchmark library for AC-OPF, [PGLib-OPF](https://github.com/power-grid-lib/pglib-opf). The Julia package [PowerModels](https://github.com/lanl-ansi/PowerModels.jl) is used for parsing the problem data files and making standard data transformations.

Expand Down
3 changes: 1 addition & 2 deletions jump.jl
Original file line number Diff line number Diff line change
Expand Up @@ -110,8 +110,7 @@ function solve_opf(file_name)
model_variables = JuMP.num_variables(model)

# for consistency with other solvers, skip the variable bounds in the constraint count
non_nl_constraints = sum(JuMP.num_constraints(model, ft, st) for (ft, st) in JuMP.list_of_constraint_types(model) if ft != JuMP.VariableRef)
model_constraints = JuMP.num_nonlinear_constraints(model) + non_nl_constraints
model_constraints = JuMP.num_constraints(model; count_variable_in_set_constraints = false)

model_build_time = time() - time_model_start

Expand Down
157 changes: 1 addition & 156 deletions test/runtests.jl
Original file line number Diff line number Diff line change
@@ -1,161 +1,6 @@
using Test

function validate_result(file_name::String, result)
for k in [
"case",
"variables",
"constraints",
"feasible",
"cost",
"time_total",
"time_data",
"time_build",
"time_solve",
"solution",
]
@test haskey(result, k)
end
@test result["feasible"] === true
validate_solution(file_name, result["solution"], result["cost"])
return
end

function validate_solution(
file_name::String,
x::Dict{String,Float64},
cost::Float64;
atol = 1e-6,
rtol = 1e-6,
)
data = PowerModels.parse_file(file_name)
PowerModels.standardize_cost_terms!(data, order=2)
PowerModels.calc_thermal_limits!(data)
ref = PowerModels.build_ref(data)[:it][:pm][:nw][0]
va = Dict{Int,Float64}(i => x["va_$i"] for i in keys(ref[:bus]))
vm = Dict{Int,Float64}()
for i in keys(ref[:bus])
vm[i] = x["vm_$i"]
@test vm[i] <= ref[:bus][i]["vmax"] + atol
@test ref[:bus][i]["vmin"] <= vm[i] + atol
end
pg = Dict{Int,Float64}()
for i in keys(ref[:gen])
pg[i] = x["pg_$i"]
@test pg[i] <= ref[:gen][i]["pmax"] + atol
@test ref[:gen][i]["pmin"] <= pg[i] + atol
end
qg = Dict{Int,Float64}()
for i in keys(ref[:gen])
qg[i] = x["qg_$i"]
@test qg[i] <= ref[:gen][i]["qmax"] + atol
@test ref[:gen][i]["qmin"] <= qg[i] + atol
end
p = Dict{NTuple{3,Int},Float64}()
for (l,i,j) in ref[:arcs]
p[(l, i, j)] = x["p_$(l)_$(i)_$(j)"]
@test p[(l, i, j)] <= ref[:branch][l]["rate_a"] + atol
@test -ref[:branch][l]["rate_a"] <= p[(l, i, j)] + atol
end
q = Dict{NTuple{3,Int},Float64}()
for (l,i,j) in ref[:arcs]
q[(l, i, j)] = x["q_$(l)_$(i)_$(j)"]
@test q[(l, i, j)] <= ref[:branch][l]["rate_a"] + atol
@test -ref[:branch][l]["rate_a"] <= q[(l, i, j)] + atol
end
# Test objective function
@test isapprox(
sum(gen["cost"][1]*pg[i]^2 + gen["cost"][2]*pg[i] + gen["cost"][3] for (i,gen) in ref[:gen]),
cost;
atol,
rtol,
)
for (i,bus) in ref[:ref_buses]
@test isapprox(va[i], 0; atol = atol)
end

for (i,bus) in ref[:bus]
bus_loads = [ref[:load][l] for l in ref[:bus_loads][i]]
bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][i]]
@test isapprox(
sum(p[a] for a in ref[:bus_arcs][i]; init = 0.0),
sum(pg[g] for g in ref[:bus_gens][i]; init = 0.0) -
sum(load["pd"] for load in bus_loads; init = 0.0) -
sum(shunt["gs"] for shunt in bus_shunts; init = 0.0)*vm[i]^2;
atol,
rtol,
)

@test isapprox(
sum(q[a] for a in ref[:bus_arcs][i]; init = 0.0),
sum(qg[g] for g in ref[:bus_gens][i]; init = 0.0) -
sum(load["qd"] for load in bus_loads; init = 0.0) +
sum(shunt["bs"] for shunt in bus_shunts; init = 0.0)*vm[i]^2;
atol,
rtol,
)
end

# Branch power flow physics and limit constraints
for (i,branch) in ref[:branch]
f_idx = (i, branch["f_bus"], branch["t_bus"])
t_idx = (i, branch["t_bus"], branch["f_bus"])

p_fr = p[f_idx]
q_fr = q[f_idx]
p_to = p[t_idx]
q_to = q[t_idx]

vm_fr = vm[branch["f_bus"]]
vm_to = vm[branch["t_bus"]]
va_fr = va[branch["f_bus"]]
va_to = va[branch["t_bus"]]

g, b = PowerModels.calc_branch_y(branch)
tr, ti = PowerModels.calc_branch_t(branch)
ttm = tr^2 + ti^2
g_fr = branch["g_fr"]
b_fr = branch["b_fr"]
g_to = branch["g_to"]
b_to = branch["b_to"]

# From side of the branch flow
@test isapprox(
p_fr,
(g+g_fr)/ttm*vm_fr^2 + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to));
atol,
rtol,
)
@test isapprox(
q_fr,
-(b+b_fr)/ttm*vm_fr^2 - (-b*tr-g*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to));
atol,
rtol,
)

# To side of the branch flow
@test isapprox(
p_to,
(g+g_to)*vm_to^2 + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr));
atol,
rtol,
)
@test isapprox(
q_to,
-(b+b_to)*vm_to^2 - (-b*tr+g*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr));
atol,
rtol,
)

# Voltage angle difference limit
@test branch["angmin"] <= va_fr - va_to + atol
@test va_fr - va_to <= branch["angmax"] + atol

# Apparent power limit, from side and to side
@test p_fr^2 + q_fr^2 <= branch["rate_a"]^2 + atol
@test p_to^2 + q_to^2 <= branch["rate_a"]^2 + atol
end
return
end
include("validator.jl")

@testset "Rosetta OPF" begin
@testset "$framework" for framework in [
Expand Down
159 changes: 159 additions & 0 deletions test/validator.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,159 @@
using Test
import PowerModels

function validate_result(file_name::String, result)
for k in [
"case",
"variables",
"constraints",
"feasible",
"cost",
"time_total",
"time_data",
"time_build",
"time_solve",
"solution",
]
@test haskey(result, k)
end
@test result["feasible"] === true
validate_solution(file_name, result["solution"], result["cost"])
return
end

function validate_solution(
file_name::String,
x::Dict{String,Float64},
cost::Float64;
atol = 1e-6,
rtol = 1e-6,
)
data = PowerModels.parse_file(file_name)
PowerModels.standardize_cost_terms!(data, order=2)
PowerModels.calc_thermal_limits!(data)
ref = PowerModels.build_ref(data)[:it][:pm][:nw][0]
va = Dict{Int,Float64}(i => x["va_$i"] for i in keys(ref[:bus]))
vm = Dict{Int,Float64}()
for i in keys(ref[:bus])
vm[i] = x["vm_$i"]
@test vm[i] <= ref[:bus][i]["vmax"] + atol
@test ref[:bus][i]["vmin"] <= vm[i] + atol
end
pg = Dict{Int,Float64}()
for i in keys(ref[:gen])
pg[i] = x["pg_$i"]
@test pg[i] <= ref[:gen][i]["pmax"] + atol
@test ref[:gen][i]["pmin"] <= pg[i] + atol
end
qg = Dict{Int,Float64}()
for i in keys(ref[:gen])
qg[i] = x["qg_$i"]
@test qg[i] <= ref[:gen][i]["qmax"] + atol
@test ref[:gen][i]["qmin"] <= qg[i] + atol
end
p = Dict{NTuple{3,Int},Float64}()
for (l,i,j) in ref[:arcs]
p[(l, i, j)] = x["p_$(l)_$(i)_$(j)"]
@test p[(l, i, j)] <= ref[:branch][l]["rate_a"] + atol
@test -ref[:branch][l]["rate_a"] <= p[(l, i, j)] + atol
end
q = Dict{NTuple{3,Int},Float64}()
for (l,i,j) in ref[:arcs]
q[(l, i, j)] = x["q_$(l)_$(i)_$(j)"]
@test q[(l, i, j)] <= ref[:branch][l]["rate_a"] + atol
@test -ref[:branch][l]["rate_a"] <= q[(l, i, j)] + atol
end
# Test objective function
@test isapprox(
sum(gen["cost"][1]*pg[i]^2 + gen["cost"][2]*pg[i] + gen["cost"][3] for (i,gen) in ref[:gen]),
cost;
atol,
rtol,
)
for (i,bus) in ref[:ref_buses]
@test isapprox(va[i], 0; atol = atol)
end

for (i,bus) in ref[:bus]
bus_loads = [ref[:load][l] for l in ref[:bus_loads][i]]
bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][i]]
@test isapprox(
sum(p[a] for a in ref[:bus_arcs][i]; init = 0.0),
sum(pg[g] for g in ref[:bus_gens][i]; init = 0.0) -
sum(load["pd"] for load in bus_loads; init = 0.0) -
sum(shunt["gs"] for shunt in bus_shunts; init = 0.0)*vm[i]^2;
atol,
rtol,
)

@test isapprox(
sum(q[a] for a in ref[:bus_arcs][i]; init = 0.0),
sum(qg[g] for g in ref[:bus_gens][i]; init = 0.0) -
sum(load["qd"] for load in bus_loads; init = 0.0) +
sum(shunt["bs"] for shunt in bus_shunts; init = 0.0)*vm[i]^2;
atol,
rtol,
)
end

# Branch power flow physics and limit constraints
for (i,branch) in ref[:branch]
f_idx = (i, branch["f_bus"], branch["t_bus"])
t_idx = (i, branch["t_bus"], branch["f_bus"])

p_fr = p[f_idx]
q_fr = q[f_idx]
p_to = p[t_idx]
q_to = q[t_idx]

vm_fr = vm[branch["f_bus"]]
vm_to = vm[branch["t_bus"]]
va_fr = va[branch["f_bus"]]
va_to = va[branch["t_bus"]]

g, b = PowerModels.calc_branch_y(branch)
tr, ti = PowerModels.calc_branch_t(branch)
ttm = tr^2 + ti^2
g_fr = branch["g_fr"]
b_fr = branch["b_fr"]
g_to = branch["g_to"]
b_to = branch["b_to"]

# From side of the branch flow
@test isapprox(
p_fr,
(g+g_fr)/ttm*vm_fr^2 + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to));
atol,
rtol,
)
@test isapprox(
q_fr,
-(b+b_fr)/ttm*vm_fr^2 - (-b*tr-g*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to));
atol,
rtol,
)

# To side of the branch flow
@test isapprox(
p_to,
(g+g_to)*vm_to^2 + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr));
atol,
rtol,
)
@test isapprox(
q_to,
-(b+b_to)*vm_to^2 - (-b*tr+g*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr));
atol,
rtol,
)

# Voltage angle difference limit
@test branch["angmin"] <= va_fr - va_to + atol
@test va_fr - va_to <= branch["angmax"] + atol

# Apparent power limit, from side and to side
@test p_fr^2 + q_fr^2 <= branch["rate_a"]^2 + atol
@test p_to^2 + q_to^2 <= branch["rate_a"]^2 + atol
end
return
end
Loading

0 comments on commit 1ff0e10

Please sign in to comment.