Skip to content

Commit a23176a

Browse files
klamikeklamike
committed
power
1 parent 68435d3 commit a23176a

6 files changed

Lines changed: 284 additions & 20 deletions

File tree

‎Project.toml‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,9 +24,11 @@ DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
2424
FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
2525
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
2626
OpenCL = "08131aa3-fb12-5dee-8b74-c09406e224a2"
27+
PGLib = "07a8691f-3d11-4330-951b-3c50f98338be"
28+
PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655"
2729
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
2830
Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
2931
pocl_jll = "627d6b7a-bbe6-5189-83e7-98cc0a5aeadd"
3032

3133
[targets]
32-
test = ["Test", "LinearAlgebra", "OpenCL", "pocl_jll", "AcceleratedKernels", "DifferentiationInterface", "FiniteDifferences", "Zygote"]
34+
test = ["Test", "LinearAlgebra", "OpenCL", "pocl_jll", "AcceleratedKernels", "DifferentiationInterface", "FiniteDifferences", "Zygote", "PGLib", "PowerModels"]

‎test/api.jl‎

Lines changed: 15 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -165,7 +165,7 @@ function test_batch_model(model::ExaModel, batch_size::Int;
165165
end
166166
end
167167

168-
@testset "API" begin
168+
@testset "API - Luksan" begin
169169
models, names = create_luksan_models()
170170

171171
for (name, model) in zip(names, models)
@@ -178,3 +178,17 @@ end
178178
end
179179
end
180180
end
181+
182+
@testset "API - Power" begin
183+
models_p, names_p = create_power_models()
184+
185+
for (name, model) in zip(names_p, models_p)
186+
@testset "$(name) Model" begin
187+
for batch_size in [1, 2, 4]
188+
@testset "Batch Size $(batch_size)" begin
189+
test_batch_model(model, batch_size; atol=1e-5, rtol=1e-5)
190+
end
191+
end
192+
end
193+
end
194+
end

‎test/power.jl‎

Lines changed: 198 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,198 @@
1+
function _get_case_file(filename::String)
2+
isfile(filename) && return filename
3+
4+
cached = joinpath(PGLib.PGLib_opf, filename)
5+
!isfile(cached) && error("File $filename not found in PGLib/pglib-opf")
6+
return cached
7+
end
8+
function _build_power_ref(filename)
9+
path = _get_case_file(filename)
10+
data = PowerModels.parse_file(path)
11+
PowerModels.standardize_cost_terms!(data, order = 2)
12+
PowerModels.calc_thermal_limits!(data)
13+
return PowerModels.build_ref(data)[:it][:pm][:nw][0]
14+
end
15+
16+
17+
_convert_array(data::N, backend) where {names,N<:NamedTuple{names}} =
18+
NamedTuple{names}(ExaModels.convert_array(d, backend) for d in data)
19+
function _parse_ac_data_raw(filename)
20+
ref = _build_power_ref(filename) # FIXME: only parse once
21+
arcdict = Dict(a => k for (k, a) in enumerate(ref[:arcs]))
22+
busdict = Dict(k => i for (i, (k, _)) in enumerate(ref[:bus]))
23+
gendict = Dict(k => i for (i, (k, _)) in enumerate(ref[:gen]))
24+
branchdict = Dict(k => i for (i, (k, _)) in enumerate(ref[:branch]))
25+
return (
26+
bus = [
27+
begin
28+
loads = [ref[:load][l] for l in ref[:bus_loads][k]]
29+
shunts = [ref[:shunt][s] for s in ref[:bus_shunts][k]]
30+
pd = sum(load["pd"] for load in loads; init = 0.0)
31+
qd = sum(load["qd"] for load in loads; init = 0.0)
32+
gs = sum(shunt["gs"] for shunt in shunts; init = 0.0)
33+
bs = sum(shunt["bs"] for shunt in shunts; init = 0.0)
34+
(i = busdict[k], pd = pd, gs = gs, qd = qd, bs = bs)
35+
end for (k, _) in ref[:bus]
36+
],
37+
gen = [
38+
(
39+
i = gendict[k],
40+
cost1 = v["cost"][1],
41+
cost2 = v["cost"][2],
42+
cost3 = v["cost"][3],
43+
bus = busdict[v["gen_bus"]],
44+
) for (k, v) in ref[:gen]
45+
],
46+
arc = [
47+
(i = k, rate_a = ref[:branch][l]["rate_a"], bus = busdict[i]) for
48+
(k, (l, i, _)) in enumerate(ref[:arcs])
49+
],
50+
branch = [
51+
begin
52+
branch = branch_raw
53+
f_idx = arcdict[i, branch["f_bus"], branch["t_bus"]]
54+
t_idx = arcdict[i, branch["t_bus"], branch["f_bus"]]
55+
56+
g, b = PowerModels.calc_branch_y(branch)
57+
tr, ti = PowerModels.calc_branch_t(branch)
58+
ttm = tr^2 + ti^2
59+
60+
g_fr = branch["g_fr"]; b_fr = branch["b_fr"]
61+
g_to = branch["g_to"]; b_to = branch["b_to"]
62+
63+
(
64+
i = branchdict[i],
65+
j = 1,
66+
f_idx = f_idx,
67+
t_idx = t_idx,
68+
f_bus = busdict[branch["f_bus"]],
69+
t_bus = busdict[branch["t_bus"]],
70+
c1 = (-g * tr - b * ti) / ttm,
71+
c2 = (-b * tr + g * ti) / ttm,
72+
c3 = (-g * tr + b * ti) / ttm,
73+
c4 = (-b * tr - g * ti) / ttm,
74+
c5 = (g + g_fr) / ttm,
75+
c6 = (b + b_fr) / ttm,
76+
c7 = (g + g_to),
77+
c8 = (b + b_to),
78+
rate_a_sq = branch["rate_a"]^2,
79+
)
80+
end for (i, branch_raw) in ref[:branch]
81+
],
82+
ref_buses = [busdict[i] for (i, _) in ref[:ref_buses]],
83+
vmax = [v["vmax"] for (_, v) in ref[:bus]],
84+
vmin = [v["vmin"] for (_, v) in ref[:bus]],
85+
pmax = [v["pmax"] for (_, v) in ref[:gen]],
86+
pmin = [v["pmin"] for (_, v) in ref[:gen]],
87+
qmax = [v["qmax"] for (_, v) in ref[:gen]],
88+
qmin = [v["qmin"] for (_, v) in ref[:gen]],
89+
rate_a = [ref[:branch][l]["rate_a"] for (l, _, _) in ref[:arcs]],
90+
angmax = [b["angmax"] for (_, b) in ref[:branch]],
91+
angmin = [b["angmin"] for (_, b) in ref[:branch]],
92+
)
93+
end
94+
95+
_parse_ac_data(filename) = _parse_ac_data_raw(filename)
96+
function _parse_ac_data(filename, backend)
97+
_convert_array(_parse_ac_data_raw(filename), backend)
98+
end
99+
100+
function create_ac_power_model(
101+
filename::String = "pglib_opf_case14_ieee.m";
102+
prod::Bool = true,
103+
backend = OpenCLBackend(),
104+
)
105+
data = _parse_ac_data(filename, backend)
106+
107+
c = ExaCore(backend = backend)
108+
109+
# Decision variables ------------------------------------------------------
110+
va = variable(c, length(data.bus)) # voltage angles (rad)
111+
vm = variable(c, length(data.bus);
112+
start = fill!(similar(data.bus, Float64), 1.0),
113+
lvar = data.vmin, uvar = data.vmax) # voltage magnitude (p.u.)
114+
115+
pg = variable(c, length(data.gen); lvar = data.pmin, uvar = data.pmax)
116+
qg = variable(c, length(data.gen); lvar = data.qmin, uvar = data.qmax)
117+
118+
p = variable(c, length(data.arc); lvar = -data.rate_a, uvar = data.rate_a)
119+
q = variable(c, length(data.arc); lvar = -data.rate_a, uvar = data.rate_a)
120+
121+
# Objective – quadratic generation cost -----------------------------------
122+
objective(c, g.cost1 * pg[g.i]^2 + g.cost2 * pg[g.i] + g.cost3 for g in data.gen)
123+
124+
# Reference bus angle ------------------------------------------------------
125+
c1 = constraint(c, va[i] for i in data.ref_buses)
126+
127+
# Branch power-flow equations ---------------------------------------------
128+
constraint(
129+
c,
130+
(p[b.f_idx] - b.c5 * vm[b.f_bus]^2 -
131+
b.c3 * (vm[b.f_bus] * vm[b.t_bus] * cos(va[b.f_bus] - va[b.t_bus])) -
132+
b.c4 * (vm[b.f_bus] * vm[b.t_bus] * sin(va[b.f_bus] - va[b.t_bus])) for
133+
b in data.branch),
134+
)
135+
136+
constraint(
137+
c,
138+
(q[b.f_idx] + b.c6 * vm[b.f_bus]^2 +
139+
b.c4 * (vm[b.f_bus] * vm[b.t_bus] * cos(va[b.f_bus] - va[b.t_bus])) -
140+
b.c3 * (vm[b.f_bus] * vm[b.t_bus] * sin(va[b.f_bus] - va[b.t_bus])) for
141+
b in data.branch),
142+
)
143+
144+
constraint(
145+
c,
146+
(p[b.t_idx] - b.c7 * vm[b.t_bus]^2 -
147+
b.c1 * (vm[b.t_bus] * vm[b.f_bus] * cos(va[b.t_bus] - va[b.f_bus])) -
148+
b.c2 * (vm[b.t_bus] * vm[b.f_bus] * sin(va[b.t_bus] - va[b.f_bus])) for
149+
b in data.branch),
150+
)
151+
152+
constraint(
153+
c,
154+
(q[b.t_idx] + b.c8 * vm[b.t_bus]^2 +
155+
b.c2 * (vm[b.t_bus] * vm[b.f_bus] * cos(va[b.t_bus] - va[b.f_bus])) -
156+
b.c1 * (vm[b.t_bus] * vm[b.f_bus] * sin(va[b.t_bus] - va[b.f_bus])) for
157+
b in data.branch),
158+
)
159+
160+
# Angle difference limits --------------------------------------------------
161+
constraint(
162+
c,
163+
va[b.f_bus] - va[b.t_bus] for b in data.branch; lcon = data.angmin, ucon = data.angmax,
164+
)
165+
166+
# Apparent power thermal limits -------------------------------------------
167+
constraint(
168+
c,
169+
p[b.f_idx]^2 + q[b.f_idx]^2 - b.rate_a_sq for b in data.branch;
170+
lcon = fill!(similar(data.branch, Float64, length(data.branch)), -Inf),
171+
)
172+
constraint(
173+
c,
174+
p[b.t_idx]^2 + q[b.t_idx]^2 - b.rate_a_sq for b in data.branch;
175+
lcon = fill!(similar(data.branch, Float64, length(data.branch)), -Inf),
176+
)
177+
178+
# Power balance at each bus -----------------------------------------------
179+
load_balance_p = constraint(c, b.pd + b.gs * vm[b.i]^2 for b in data.bus)
180+
load_balance_q = constraint(c, b.qd - b.bs * vm[b.i]^2 for b in data.bus)
181+
182+
# Map arc & generator variables into the bus balance equations
183+
constraint!(c, load_balance_p, a.bus => p[a.i] for a in data.arc)
184+
constraint!(c, load_balance_q, a.bus => q[a.i] for a in data.arc)
185+
constraint!(c, load_balance_p, g.bus => -pg[g.i] for g in data.gen)
186+
constraint!(c, load_balance_q, g.bus => -qg[g.i] for g in data.gen)
187+
188+
return ExaModel(c; prod = prod)
189+
end
190+
191+
# TODO: parametric version
192+
193+
function create_power_models(backend = OpenCLBackend())
194+
models = ExaModel[]
195+
push!(models, create_ac_power_model("pglib_opf_case14_ieee.m"; backend = backend))
196+
names = ["AC-OPF – IEEE-14"]
197+
return models, names
198+
end

‎test/runtests.jl‎

Lines changed: 12 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,16 @@ using BatchNLPKernels
22
using Test
33
using ExaModels
44
using KernelAbstractions
5+
using DifferentiationInterface
6+
const DI = DifferentiationInterface
57

8+
import Zygote
9+
import FiniteDifferences
10+
11+
using PowerModels
12+
PowerModels.silence()
13+
using PGLib
14+
using LinearAlgebra
615

716
using OpenCL, pocl_jll, AcceleratedKernels
817
ExaModels.convert_array(x, ::OpenCLBackend) = CLArray(x)
@@ -17,7 +26,8 @@ Base.findall(bitarray::CLArray) = Base.findall(identity, bitarray)
1726

1827

1928
include("luksan.jl")
29+
include("power.jl")
30+
include("test_viols.jl")
2031
include("test_diff.jl")
2132
include("api.jl")
22-
include("config.jl")
23-
include("test_viols.jl")
33+
include("config.jl")

‎test/test_diff.jl‎

Lines changed: 21 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,3 @@
1-
using DifferentiationInterface
2-
const DI = DifferentiationInterface
3-
4-
import Zygote
5-
import FiniteDifferences
6-
7-
81
function test_diff_gpu(model::ExaModel, batch_size::Int)
92
bm = BOI.BatchModel(model, batch_size, config=BOI.BatchModelConfig(:full))
103

@@ -111,7 +104,7 @@ function test_diff_cpu(model::ExaModel, batch_size::Int)
111104
end
112105

113106

114-
@testset "AD rules" begin
107+
@testset "AD rules - Luksan" begin
115108
cpu_models, names = create_luksan_models(CPU())
116109
gpu_models, _ = create_luksan_models(OpenCLBackend())
117110

@@ -130,3 +123,23 @@ end
130123
end
131124
end
132125
end
126+
127+
@testset "AD rules - Power" begin
128+
cpu_models_p, names_p = create_power_models(CPU())
129+
gpu_models_p, _ = create_power_models(OpenCLBackend())
130+
131+
for (name, (cpu_model, gpu_model)) in zip(names_p, zip(cpu_models_p, gpu_models_p))
132+
@testset "$(name) Model" begin
133+
for batch_size in [1, 4]
134+
@testset "Batch Size $(batch_size)" begin
135+
@testset "CPU Diff" begin
136+
test_diff_cpu(cpu_model, batch_size)
137+
end
138+
@testset "GPU Diff" begin
139+
test_diff_gpu(gpu_model, batch_size)
140+
end
141+
end
142+
end
143+
end
144+
end
145+
end

‎test/test_viols.jl‎

Lines changed: 35 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,3 @@
1-
using DifferentiationInterface
2-
const DI = DifferentiationInterface
3-
4-
import Zygote
5-
import FiniteDifferences
6-
71
function test_violations_correctness(model::ExaModel, batch_size::Int;
82
atol::Float64=1e-10, rtol::Float64=1e-10)
93
bm = BOI.BatchModel(model, batch_size, config=BOI.BatchModelConfig(:violations))
@@ -36,6 +30,16 @@ function test_violations_correctness(model::ExaModel, batch_size::Int;
3630
@test all(>=(0), Vb)
3731
@test all(isfinite, Vc)
3832
@test all(isfinite, Vb)
33+
OpenCL.@allowscalar begin
34+
if !isempty(model.meta.lvar) && !isempty(model.meta.uvar)
35+
if isfinite(model.meta.lvar[1])
36+
@test Vb[1, :] .≈ 0.1 atol=atol rtol=rtol
37+
end
38+
if isfinite(model.meta.uvar[end])
39+
@test Vb[end, :] .≈ 0.1 atol=atol rtol=rtol
40+
end
41+
end
42+
end
3943
end
4044

4145
Vc, Vb = BOI.all_violations!(bm, X)
@@ -266,7 +270,7 @@ function test_violations_config_errors()
266270
end
267271
end
268272

269-
@testset "Violations API" begin
273+
@testset "Violations API - Luksan" begin
270274
@testset "Config Errors" begin
271275
test_violations_config_errors()
272276
end
@@ -276,7 +280,7 @@ end
276280

277281
for (name, (cpu_model, gpu_model)) in zip(names, zip(cpu_models, gpu_models))
278282
@testset "$name Model" begin
279-
for batch_size in [1, 2, 4]
283+
for batch_size in [1, 4]
280284
@testset "Batch Size $batch_size" begin
281285
@testset "Correctness" begin
282286
test_violations_correctness(gpu_model, batch_size, atol=1e-5, rtol=1e-5)
@@ -291,4 +295,27 @@ end
291295
end
292296
end
293297
end
298+
end
299+
300+
@testset "Violations API - Power" begin
301+
cpu_models_p, names_p = create_power_models(CPU())
302+
gpu_models_p, _ = create_power_models(OpenCLBackend())
303+
304+
for (name, (cpu_model, gpu_model)) in zip(names_p, zip(cpu_models_p, gpu_models_p))
305+
@testset "$(name) Model" begin
306+
for batch_size in [1, 4]
307+
@testset "Batch Size $(batch_size)" begin
308+
@testset "Correctness" begin
309+
test_violations_correctness(gpu_model, batch_size; atol=1e-5, rtol=1e-5)
310+
end
311+
@testset "CPU Differentiability" begin
312+
test_violations_differentiability_cpu(cpu_model, batch_size)
313+
end
314+
@testset "GPU Differentiability" begin
315+
test_violations_differentiability_gpu(gpu_model, batch_size)
316+
end
317+
end
318+
end
319+
end
320+
end
294321
end

0 commit comments

Comments
 (0)