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