Skip to content

Commit 5c96182

Browse files
klamikeklamike
committed
add test with acopf taylor
1 parent 83627df commit 5c96182

2 files changed

Lines changed: 247 additions & 0 deletions

File tree

‎test/Project.toml‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,8 @@ Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
33
Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
44
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
55
CoolPDLP = "9e90bdc5-3073-4e0f-a275-e19f28ac1b53"
6+
Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
7+
ExaModels = "1037b233-b668-4ce9-9b63-f9f681f55dd2"
68
ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
79
GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
810
HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
@@ -13,6 +15,14 @@ KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
1315
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
1416
MathOptBenchmarkInstances = "f7f8d0a1-fd34-491e-a7ac-a4cf52f91fe5"
1517
MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
18+
NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6"
19+
NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
20+
PSLP = "ef7f38f3-b543-4219-892d-e0f472564133"
21+
PSLP_jll = "f21d708b-afb1-4941-8a85-571bbd38ad3d"
22+
PaPILO = "d9c0c21a-f458-4936-9625-08537585e348"
23+
PaPILO_Presolve_jll = "ac3407f7-d818-4657-ac2d-3679289da1ef"
24+
PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655"
25+
QuadraticModels = "f468eda6-eac5-11e8-05a5-ff9e497bcd19"
1626
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
1727
SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
1828
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -21,6 +31,10 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
2131

2232
[sources]
2333
CoolPDLP = {path = ".."}
34+
PSLP = {rev = "main", url = "https://github.com/klamike/PSLP.jl.git"}
35+
PSLP_jll = {rev = "main", subdir = "PSLP_jll", url = "https://github.com/klamike/PSLP.jl.git"}
36+
PaPILO = {rev = "main", url = "https://github.com/klamike/PaPILO.jl.git"}
37+
PaPILO_Presolve_jll = {rev = "main", subdir = "PaPILO_Presolve_jll", url = "https://github.com/klamike/PaPILO.jl.git"}
2438

2539
[compat]
2640
Aqua = "0.8.14"

‎test/extensions/QuadraticModels/power.jl‎

Lines changed: 233 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,233 @@
1+
using ExaModels, KernelAbstractions
2+
using Downloads, Test, SparseArrays
3+
using NLPModels, NLPModelsIpopt
4+
using PowerModels; PowerModels.silence()
5+
6+
using CoolPDLP, QuadraticModels
7+
8+
const TMPDIR = tempdir()
9+
10+
function get_power_case(filename)
11+
if !isfile(filename)
12+
ff = joinpath(TMPDIR, filename)
13+
if !isfile(ff)
14+
@info "Downloading $filename"
15+
Downloads.download(
16+
"https://raw.githubusercontent.com/power-grid-lib/pglib-opf/dc6be4b2f85ca0e776952ec22cbd4c22396ea5a3/$filename",
17+
joinpath(TMPDIR, filename),
18+
)
19+
return joinpath(TMPDIR, filename)
20+
else
21+
return ff
22+
end
23+
else
24+
return filename
25+
end
26+
end
27+
28+
29+
function get_power_data_ref(filename)
30+
case = get_power_case(filename)
31+
data = PowerModels.parse_file(case)
32+
PowerModels.standardize_cost_terms!(data, order = 2)
33+
PowerModels.calc_thermal_limits!(data)
34+
return PowerModels.build_ref(data)[:it][:pm][:nw][0]
35+
end
36+
37+
convert_data(data::N, backend) where {names,N<:NamedTuple{names}} =
38+
NamedTuple{names}(ExaModels.convert_array(d, backend) for d in data)
39+
parse_ac_power_data(filename, backend) =
40+
convert_data(parse_ac_power_data(filename), backend)
41+
42+
43+
function parse_ac_power_data(filename)
44+
ref = get_power_data_ref(filename)
45+
46+
arcdict = Dict(a => k for (k, a) in enumerate(ref[:arcs]))
47+
busdict = Dict(k => i for (i, (k, v)) in enumerate(ref[:bus]))
48+
gendict = Dict(k => i for (i, (k, v)) in enumerate(ref[:gen]))
49+
branchdict = Dict(k => i for (i, (k, v)) in enumerate(ref[:branch]))
50+
51+
return (
52+
bus = [
53+
begin
54+
bus_loads = [ref[:load][l] for l in ref[:bus_loads][k]]
55+
bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][k]]
56+
pd = sum(load["pd"] for load in bus_loads; init = 0.0)
57+
gs = sum(shunt["gs"] for shunt in bus_shunts; init = 0.0)
58+
qd = sum(load["qd"] for load in bus_loads; init = 0.0)
59+
bs = sum(shunt["bs"] for shunt in bus_shunts; init = 0.0)
60+
(i = busdict[k], pd = pd, gs = gs, qd = qd, bs = bs)
61+
end for (k, v) in ref[:bus]
62+
],
63+
gen = [
64+
(
65+
i = gendict[k],
66+
cost1 = v["cost"][1],
67+
cost2 = v["cost"][2],
68+
cost3 = v["cost"][3],
69+
bus = busdict[v["gen_bus"]],
70+
) for (k, v) in ref[:gen]
71+
],
72+
arc = [
73+
(i = k, rate_a = ref[:branch][l]["rate_a"], bus = busdict[i]) for
74+
(k, (l, i, j)) in enumerate(ref[:arcs])
75+
],
76+
branch = [
77+
begin
78+
f_idx = arcdict[i, branch["f_bus"], branch["t_bus"]]
79+
t_idx = arcdict[i, branch["t_bus"], branch["f_bus"]]
80+
g, b = PowerModels.calc_branch_y(branch)
81+
tr, ti = PowerModels.calc_branch_t(branch)
82+
ttm = tr^2 + ti^2
83+
g_fr = branch["g_fr"]
84+
b_fr = branch["b_fr"]
85+
g_to = branch["g_to"]
86+
b_to = branch["b_to"]
87+
c1 = (-g * tr - b * ti) / ttm
88+
c2 = (-b * tr + g * ti) / ttm
89+
c3 = (-g * tr + b * ti) / ttm
90+
c4 = (-b * tr - g * ti) / ttm
91+
c5 = (g + g_fr) / ttm
92+
c6 = (b + b_fr) / ttm
93+
c7 = (g + g_to)
94+
c8 = (b + b_to)
95+
(
96+
i = branchdict[i],
97+
j = 1,
98+
f_idx = f_idx,
99+
t_idx = t_idx,
100+
f_bus = busdict[branch["f_bus"]],
101+
t_bus = busdict[branch["t_bus"]],
102+
c1 = c1,
103+
c2 = c2,
104+
c3 = c3,
105+
c4 = c4,
106+
c5 = c5,
107+
c6 = c6,
108+
c7 = c7,
109+
c8 = c8,
110+
rate_a_sq = branch["rate_a"]^2,
111+
)
112+
end for (i, branch) in ref[:branch]
113+
],
114+
ref_buses = [busdict[i] for (i, k) in ref[:ref_buses]],
115+
vmax = [v["vmax"] for (k, v) in ref[:bus]],
116+
vmin = [v["vmin"] for (k, v) in ref[:bus]],
117+
pmax = [v["pmax"] for (k, v) in ref[:gen]],
118+
pmin = [v["pmin"] for (k, v) in ref[:gen]],
119+
qmax = [v["qmax"] for (k, v) in ref[:gen]],
120+
qmin = [v["qmin"] for (k, v) in ref[:gen]],
121+
rate_a = [ref[:branch][l]["rate_a"] for (k, (l, i, j)) in enumerate(ref[:arcs])],
122+
angmax = [b["angmax"] for (i, b) in ref[:branch]],
123+
angmin = [b["angmin"] for (i, b) in ref[:branch]],
124+
)
125+
end
126+
127+
function exa_ac_power_model(backend, data)
128+
129+
w = ExaModels.ExaCore(backend = backend)
130+
131+
va = ExaModels.variable(w, length(data.bus);)
132+
133+
vm = ExaModels.variable(
134+
w,
135+
length(data.bus);
136+
start = fill!(similar(data.bus, Float64), 1.0),
137+
lvar = data.vmin,
138+
uvar = data.vmax,
139+
)
140+
pg = ExaModels.variable(w, length(data.gen); lvar = data.pmin, uvar = data.pmax)
141+
142+
qg = ExaModels.variable(w, length(data.gen); lvar = data.qmin, uvar = data.qmax)
143+
144+
p = ExaModels.variable(w, length(data.arc); lvar = -data.rate_a, uvar = data.rate_a)
145+
146+
q = ExaModels.variable(w, length(data.arc); lvar = -data.rate_a, uvar = data.rate_a)
147+
148+
o = ExaModels.objective(
149+
w,
150+
g.cost1 * pg[g.i]^2 + g.cost2 * pg[g.i] + g.cost3 for g in data.gen
151+
)
152+
153+
c1 = ExaModels.constraint(w, va[i] for i in data.ref_buses)
154+
155+
c2 = ExaModels.constraint(
156+
w,
157+
p[b.f_idx] - b.c5 * vm[b.f_bus]^2 -
158+
b.c3 * (vm[b.f_bus] * vm[b.t_bus] * cos(va[b.f_bus] - va[b.t_bus])) -
159+
b.c4 * (vm[b.f_bus] * vm[b.t_bus] * sin(va[b.f_bus] - va[b.t_bus])) for
160+
b in data.branch
161+
)
162+
163+
c3 = ExaModels.constraint(
164+
w,
165+
q[b.f_idx] +
166+
b.c6 * vm[b.f_bus]^2 +
167+
b.c4 * (vm[b.f_bus] * vm[b.t_bus] * cos(va[b.f_bus] - va[b.t_bus])) -
168+
b.c3 * (vm[b.f_bus] * vm[b.t_bus] * sin(va[b.f_bus] - va[b.t_bus])) for
169+
b in data.branch
170+
)
171+
172+
c4 = ExaModels.constraint(
173+
w,
174+
p[b.t_idx] - b.c7 * vm[b.t_bus]^2 -
175+
b.c1 * (vm[b.t_bus] * vm[b.f_bus] * cos(va[b.t_bus] - va[b.f_bus])) -
176+
b.c2 * (vm[b.t_bus] * vm[b.f_bus] * sin(va[b.t_bus] - va[b.f_bus])) for
177+
b in data.branch
178+
)
179+
180+
c5 = ExaModels.constraint(
181+
w,
182+
q[b.t_idx] +
183+
b.c8 * vm[b.t_bus]^2 +
184+
b.c2 * (vm[b.t_bus] * vm[b.f_bus] * cos(va[b.t_bus] - va[b.f_bus])) -
185+
b.c1 * (vm[b.t_bus] * vm[b.f_bus] * sin(va[b.t_bus] - va[b.f_bus])) for
186+
b in data.branch
187+
)
188+
189+
c6 = ExaModels.constraint(
190+
w,
191+
va[b.f_bus] - va[b.t_bus] for b in data.branch;
192+
lcon = data.angmin,
193+
ucon = data.angmax,
194+
)
195+
c7 = ExaModels.constraint(
196+
w,
197+
p[b.f_idx]^2 + q[b.f_idx]^2 - b.rate_a_sq for b in data.branch;
198+
lcon = fill!(similar(data.branch, Float64, length(data.branch)), -Inf),
199+
)
200+
c8 = ExaModels.constraint(
201+
w,
202+
p[b.t_idx]^2 + q[b.t_idx]^2 - b.rate_a_sq for b in data.branch;
203+
lcon = fill!(similar(data.branch, Float64, length(data.branch)), -Inf),
204+
)
205+
206+
c9 = ExaModels.constraint(w, b.pd + b.gs * vm[b.i]^2 for b in data.bus)
207+
208+
c10 = ExaModels.constraint(w, b.qd - b.bs * vm[b.i]^2 for b in data.bus)
209+
210+
c11 = ExaModels.constraint!(w, c9, a.bus => p[a.i] for a in data.arc)
211+
c12 = ExaModels.constraint!(w, c10, a.bus => q[a.i] for a in data.arc)
212+
213+
c13 = ExaModels.constraint!(w, c9, g.bus => -pg[g.i] for g in data.gen)
214+
c14 = ExaModels.constraint!(w, c10, g.bus => -qg[g.i] for g in data.gen)
215+
216+
return ExaModels.ExaModel(w; prod = true),
217+
(va, vm, pg, qg, p, q),
218+
(c1, c2, c3, c4, c5, c6, c7, c8, c9, c10)
219+
220+
end
221+
222+
@testset "Taylor of ACOPF via QuadraticModels.jl " begin
223+
backend = CPU()
224+
data = parse_ac_power_data("pglib_opf_case14_ieee.m", backend)
225+
examodel = first(exa_ac_power_model(backend, data))
226+
res = ipopt(examodel, print_level = 0)
227+
milp = CoolPDLP.MILP(QuadraticModel(examodel, res.solution), ignore_islp = true)
228+
sol, cres = CoolPDLP.solve(milp, PDLP(
229+
Float32, Int32, SparseMatrixCSC; backend,
230+
termination_reltol = 1.0f-4, time_limit = 10.0,
231+
))
232+
@test cres.termination_status == CoolPDLP.OPTIMAL
233+
end

0 commit comments

Comments
 (0)