From de1f6fd61205320ca4b7a292f2b1163d2c7793ca Mon Sep 17 00:00:00 2001 From: David M Fobes Date: Wed, 18 Sep 2024 11:57:38 -0600 Subject: [PATCH] UPD: JuMP 1.23.2 --- CHANGELOG.md | 4 + CONTRIBUTING.md | 2 +- Project.toml | 6 +- src/core/constraint.jl | 6 +- src/core/generator.jl | 2 +- src/core/inverter.jl | 322 ++++++++++++++++----------------- src/io/preprocessing.jl | 2 +- test/data/case2_inverters.json | 60 ++---- test/runtests.jl | 4 +- test/test_2bus.jl | 2 +- 10 files changed, 193 insertions(+), 217 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ed4be12..10095f3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,10 @@ - none +## v0.4.0 + +- Updated to JuMP 1.23.2 to use new NL syntax + ## v0.3.3 - Update compat for PowerModelsDistribution v0.15 diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 5bfa202..30e8f34 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -68,7 +68,7 @@ All new dependencies should be carefully considered before being added. It is im All new dependencies are are ultimately approved should also include an entry under `[compat]` indicating the acceptable versions (Julia automerge requirement). This includes test-only dependencies that appear under `[extras]` -The `Manifest.toml` __should not__ be included in the repo. +The `Manifest.toml` **should not** be included in the repo. ## Pull Requests diff --git a/Project.toml b/Project.toml index 03a3349..aa2ee65 100644 --- a/Project.toml +++ b/Project.toml @@ -2,7 +2,7 @@ name = "PowerModelsStability" uuid = "f9e4c324-c3b6-4bca-9c3d-419775f0bd17" authors = ["Haoxiang Yang", "Harsha Nagarajan", "David M Fobes"] repo = "https://github.com/lanl-ansi/PowerModelsStability.jl.git" -version = "0.3.3" +version = "0.4" [deps] InfrastructureModels = "2030c09a-7f63-5d83-885d-db604e0e9cc0" @@ -16,9 +16,9 @@ PowerModelsDistribution = "d7431456-977f-11e9-2de3-97ff7677985e" InfrastructureModels = "0.7" Ipopt = "0.9, 1" JSON = "0.21" -JuMP = "0.22, 0.23, 1" +JuMP = "1.23.2" Memento = "1" -PowerModelsDistribution = "0.14, 0.15" +PowerModelsDistribution = "0.16" julia = "1.6" [extras] diff --git a/src/core/constraint.jl b/src/core/constraint.jl index 85c221c..45561cc 100644 --- a/src/core/constraint.jl +++ b/src/core/constraint.jl @@ -2,7 +2,7 @@ "Given a matrix A with NL expression and an eigenvector v, calculate the multiplication v^\\top A v" function matrixMulti(mp::JuMP.Model, A::Array{Any,2}, v) n = length(v) - multiExpression = JuMP.@NLexpression(mp, sum(sum(v[i] * A[i,j] * v[j] for j in 1:n) for i in 1:n)) + multiExpression = JuMP.@expression(mp, sum(sum(v[i] * A[i,j] * v[j] for j in 1:n) for i in 1:n)) return multiExpression end @@ -36,7 +36,7 @@ function obtainGlobal_var(mpData, pm, rN, omega0) Isub = PMS.obtainI_inverter_global(mpData, rN, omega0, busList, brList, loadList, load_L, load_R, load_X, loadConnections) Atot = PMS.combineSub(busList, brList, inverters, invBusDict, Asub, Bsub, Csub, Dsub, Esub, Fsub, Gsub, Hsub, Isub, 2) - + return Atot end @@ -46,6 +46,6 @@ function constraint_stability(pm::_PMD.AbstractUnbalancedPowerModel, nw::Int, ei for ev in eigenVectorList evMulti_real = PMS.matrixMulti(pm.model, Amg, real(ev)) evMulti_imag = PMS.matrixMulti(pm.model, Amg, imag(ev)) - JuMP.@NLconstraint(pm.model, evMulti_real + evMulti_imag <= 0) + JuMP.@constraint(pm.model, evMulti_real + evMulti_imag <= 0) end end diff --git a/src/core/generator.jl b/src/core/generator.jl index cf6e926..d15939a 100644 --- a/src/core/generator.jl +++ b/src/core/generator.jl @@ -21,7 +21,7 @@ function obtainA_diesel(mpData, opfSol, iBusList, omega0, kD) iQ = Dict() A = Dict() error_list = [] - + for i in iBusList if mpData["bus"][i]["bus_type"] == 3 delta0[i] = opfSol["solution"]["bus"][i]["va"][1] diff --git a/src/core/inverter.jl b/src/core/inverter.jl index 50ca9db..962fda6 100644 --- a/src/core/inverter.jl +++ b/src/core/inverter.jl @@ -4,44 +4,44 @@ import LinearAlgebra: eigvals, eigvecs, dot, inv, norm "given a vector with NL expression and a constant matrix C, calculate the multiplication Cv" function vectorMulti(mp::JuMP.Model, C, v) - + n = length(v) - multiExpression = [JuMP.@NLexpression(mp, sum(C[i,j] * v[j] for j in 1:n)) for i in 1:n] + multiExpression = [JuMP.@expression(mp, sum(C[i, j] * v[j] for j in 1:n)) for i in 1:n] return multiExpression end "invert the matrix" function inverseMat(originMat, connection) - + invM = zeros(3, 3) if LA.norm(originMat) > 0 try invOri = LA.inv(originMat) - invM[connection,connection] = invOri + invM[connection, connection] = invOri catch for i in 1:length(connection) for j in 1:length(connection) - if abs(1 / originMat[i,j]) < Inf - invM[connection[i],connection[j]] = 1 / originMat[i,j] + if abs(1 / originMat[i, j]) < Inf + invM[connection[i], connection[j]] = 1 / originMat[i, j] end end end end end - + return invM end "Obtain the submatrix A, with the opf solution given" function obtainA_inverter_global(mpData, opfSol, rN, omega0, busList, invList, invLine, invConnected, inverters, invBusDict) - + A = Dict() for invItem in inverters # initialize the matrix to be zero matrices, if it is a regular bus - A[invItem,invItem] = zeros(9, 9) + A[invItem, invItem] = zeros(9, 9) end # for every inverter @@ -65,7 +65,7 @@ function obtainA_inverter_global(mpData, opfSol, rN, omega0, busList, invList, i id0 = zeros(3) iq0 = zeros(3) for j in 1:3 - id0[j], iq0[j] = LA.inv([vd0[j] vq0[j];vq0[j] -vd0[j]]) * [P0 / 3;Q0 / 3] + id0[j], iq0[j] = LA.inv([vd0[j] vq0[j]; vq0[j] -vd0[j]]) * [P0 / 3; Q0 / 3] end dvd_delta = -vq0 dvq_delta = vd0 @@ -80,45 +80,45 @@ function obtainA_inverter_global(mpData, opfSol, rN, omega0, busList, invList, i dQ_id = vq0 dQ_iq = -vd0 - A[i,i][1,2] = 1 - A[i,i][2,:] = -1 / mpData["bus"][i]["tau"] * [mpData["bus"][i]["mp"] * dP_delta - 1 - mpData["bus"][i]["mp"] * dP_v - mpData["bus"][i]["mp"] * dP_id[1]; mpData["bus"][i]["mp"] * dP_id[2]; mpData["bus"][i]["mp"] * dP_id[3] - mpData["bus"][i]["mp"] * dP_iq[1]; mpData["bus"][i]["mp"] * dP_iq[2]; mpData["bus"][i]["mp"] * dP_iq[3] - ] - A[i,i][3,:] = -1 / mpData["bus"][i]["tau"] * [mpData["bus"][i]["mq"] * dQ_delta - 0 - 1 + mpData["bus"][i]["mq"] * dQ_v - mpData["bus"][i]["mq"] * dQ_id[1]; mpData["bus"][i]["mq"] * dQ_id[2]; mpData["bus"][i]["mp"] * dQ_id[3] - mpData["bus"][i]["mq"] * dQ_iq[1]; mpData["bus"][i]["mq"] * dQ_iq[2]; mpData["bus"][i]["mp"] * dQ_iq[3] - ] - A[i,i][4:6,1] = LA.inv(L) * dvd_delta - A[i,i][4:6,3] = LA.inv(L) * dvd_v - A[i,i][4:6,4:6] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) - A[i,i][4:6,7:9] = LA.inv(L) * link["br_x"] - A[i,i][7:9,1] = LA.inv(L) * dvq_delta - A[i,i][7:9,3] = LA.inv(L) * dvq_v - A[i,i][7:9,4:6] = -LA.inv(L) * link["br_x"] - A[i,i][7:9,7:9] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) + A[i, i][1, 2] = 1 + A[i, i][2, :] = -1 / mpData["bus"][i]["tau"] * [mpData["bus"][i]["mp"] * dP_delta + 1 + mpData["bus"][i]["mp"] * dP_v + mpData["bus"][i]["mp"] * dP_id[1]; mpData["bus"][i]["mp"] * dP_id[2]; mpData["bus"][i]["mp"] * dP_id[3]; + mpData["bus"][i]["mp"] * dP_iq[1]; mpData["bus"][i]["mp"] * dP_iq[2]; mpData["bus"][i]["mp"] * dP_iq[3] + ] + A[i, i][3, :] = -1 / mpData["bus"][i]["tau"] * [mpData["bus"][i]["mq"] * dQ_delta + 0 + 1 + mpData["bus"][i]["mq"] * dQ_v + mpData["bus"][i]["mq"] * dQ_id[1]; mpData["bus"][i]["mq"] * dQ_id[2]; mpData["bus"][i]["mp"] * dQ_id[3]; + mpData["bus"][i]["mq"] * dQ_iq[1]; mpData["bus"][i]["mq"] * dQ_iq[2]; mpData["bus"][i]["mp"] * dQ_iq[3] + ] + A[i, i][4:6, 1] = LA.inv(L) * dvd_delta + A[i, i][4:6, 3] = LA.inv(L) * dvd_v + A[i, i][4:6, 4:6] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) + A[i, i][4:6, 7:9] = LA.inv(L) * link["br_x"] + A[i, i][7:9, 1] = LA.inv(L) * dvq_delta + A[i, i][7:9, 3] = LA.inv(L) * dvq_v + A[i, i][7:9, 4:6] = -LA.inv(L) * link["br_x"] + A[i, i][7:9, 7:9] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) end # for every pair of inverters in the same bus for iBus in invList if length(invConnected[iBus]) > 1 for i in 1:length(invConnected[iBus]) - for j in (i + 1):length(invConnected[iBus]) - A[i,j] = zeros(9, 9) + for j in (i+1):length(invConnected[iBus]) + A[i, j] = zeros(9, 9) link = mpData["branch"][invLine[i]] L = link["br_x"] ./ omega0 - A[i,j][4:6,4:6] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) - A[i,j][7:9,7:9] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) + A[i, j][4:6, 4:6] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) + A[i, j][7:9, 7:9] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) link = mpData["branch"][invLine[j]] L = link["br_x"] ./ omega0 - A[j,i] = zeros(9, 9) - A[j,i][4:6,4:6] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) - A[j,i][7:9,7:9] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) + A[j, i] = zeros(9, 9) + A[j, i][4:6, 4:6] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) + A[j, i][7:9, 7:9] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) end end end @@ -129,10 +129,10 @@ end "Obtain the submatrix A, with the opf problem variables in" function obtainA_inverter_global_var(mpData, pm, rN, omega0, busList, invList, invLine, invConnected, inverters, invBusDict) - + A = Dict() for invItem in inverters - A[invItem,invItem] = Array{Any,2}(zeros(9, 9)) + A[invItem, invItem] = Array{Any,2}(zeros(9, 9)) end for iBus in inverters @@ -144,10 +144,10 @@ function obtainA_inverter_global_var(mpData, pm, rN, omega0, busList, invList, i delta0 = _PMD.var(pm, 0, :va, i)[1] vmList = _PMD.var(pm, 0, :vm, i) vaList = _PMD.var(pm, 0, :va, i) - vd0 = [JuMP.@NLexpression(pm.model,vmList.data[j] * cos(vaList.data[j])) for j in 1:length(vmList.data)] - vq0 = [JuMP.@NLexpression(pm.model,vmList.data[j] * sin(vaList.data[j])) for j in 1:length(vmList.data)] - P0 = JuMP.@NLexpression(pm.model,0) - Q0 = JuMP.@NLexpression(pm.model,0) + vd0 = [JuMP.@expression(pm.model, vmList.data[j] * cos(vaList.data[j])) for j in 1:length(vmList.data)] + vq0 = [JuMP.@expression(pm.model, vmList.data[j] * sin(vaList.data[j])) for j in 1:length(vmList.data)] + P0 = JuMP.@expression(pm.model, 0) + Q0 = JuMP.@expression(pm.model, 0) for g in keys(mpData["gen"]) if "$(mpData["gen"][g]["gen_bus"])" == i @@ -155,74 +155,74 @@ function obtainA_inverter_global_var(mpData, pm, rN, omega0, busList, invList, i Q0 += sum(_PMD.var(pm, 0, :qg, g)) end end - - id0 = Array{Any,1}([0,0,0]) - iq0 = Array{Any,1}([0,0,0]) + + id0 = Array{Any,1}([0, 0, 0]) + iq0 = Array{Any,1}([0, 0, 0]) for j in 1:3 - id0[j] = JuMP.@NLexpression(pm.model,(vd0[j] * P0 / 3 + vq0[j] * Q0 / 3) / (vmList.data[j]^2)) - iq0[j] = JuMP.@NLexpression(pm.model,(vq0[j] * P0 / 3 - vd0[j] * Q0 / 3) / (vmList.data[j]^2)) + id0[j] = JuMP.@expression(pm.model, (vd0[j] * P0 / 3 + vq0[j] * Q0 / 3) / (vmList.data[j]^2)) + iq0[j] = JuMP.@expression(pm.model, (vq0[j] * P0 / 3 - vd0[j] * Q0 / 3) / (vmList.data[j]^2)) end - dvd_delta = [JuMP.@NLexpression(pm.model,-vq0[j]) for j in 1:length(vq0)] - dvq_delta = [JuMP.@NLexpression(pm.model,vd0[j]) for j in 1:length(vd0)] - dvd_v = [JuMP.@NLexpression(pm.model,cos(vaList.data[j])) for j in 1:length(vaList.data)] - dvq_v = [JuMP.@NLexpression(pm.model,sin(vaList.data[j])) for j in 1:length(vaList.data)] - dP_delta = JuMP.@NLexpression(pm.model,sum(dvd_delta[j] * id0[j] for j in 1:length(dvd_delta)) + sum(dvq_delta[j] * iq0[j] for j in 1:length(dvq_delta))) - dQ_delta = JuMP.@NLexpression(pm.model,-sum(dvd_delta[j] * iq0[j] for j in 1:length(dvd_delta)) + sum(dvq_delta[j] * id0[j] for j in 1:length(dvq_delta))) - dP_v = JuMP.@NLexpression(pm.model,sum(dvd_v[j] * id0[j] for j in 1:length(dvd_v)) + sum(dvq_v[j] * iq0[j] for j in 1:length(dvq_v))) - dQ_v = JuMP.@NLexpression(pm.model,-sum(dvd_v[j] * iq0[j] for j in 1:length(dvd_v)) + sum(dvq_v[j] * id0[j] for j in 1:length(dvq_v))) + dvd_delta = [JuMP.@expression(pm.model, -vq0[j]) for j in 1:length(vq0)] + dvq_delta = [JuMP.@expression(pm.model, vd0[j]) for j in 1:length(vd0)] + dvd_v = [JuMP.@expression(pm.model, cos(vaList.data[j])) for j in 1:length(vaList.data)] + dvq_v = [JuMP.@expression(pm.model, sin(vaList.data[j])) for j in 1:length(vaList.data)] + dP_delta = JuMP.@expression(pm.model, sum(dvd_delta[j] * id0[j] for j in 1:length(dvd_delta)) + sum(dvq_delta[j] * iq0[j] for j in 1:length(dvq_delta))) + dQ_delta = JuMP.@expression(pm.model, -sum(dvd_delta[j] * iq0[j] for j in 1:length(dvd_delta)) + sum(dvq_delta[j] * id0[j] for j in 1:length(dvq_delta))) + dP_v = JuMP.@expression(pm.model, sum(dvd_v[j] * id0[j] for j in 1:length(dvd_v)) + sum(dvq_v[j] * iq0[j] for j in 1:length(dvq_v))) + dQ_v = JuMP.@expression(pm.model, -sum(dvd_v[j] * iq0[j] for j in 1:length(dvd_v)) + sum(dvq_v[j] * id0[j] for j in 1:length(dvq_v))) dP_id = vd0 dP_iq = vq0 dQ_id = vq0 - dQ_iq = [JuMP.@NLexpression(pm.model,-vd0[i]) for i in 1:length(vd0)] - - A[iBus,iBus][1,2] = 1 - A[iBus,iBus][2,:] = [JuMP.@NLexpression(pm.model,-1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mp"] * dP_delta) - -1 / mpData["bus"][iBus]["tau"] - JuMP.@NLexpression(pm.model,-1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mp"] * dP_v) - JuMP.@NLexpression(pm.model,-1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mp"] * dP_id[1]) - JuMP.@NLexpression(pm.model,-1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mp"] * dP_id[2]) - JuMP.@NLexpression(pm.model,-1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mp"] * dP_id[3]) - JuMP.@NLexpression(pm.model,-1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mp"] * dP_iq[1]) - JuMP.@NLexpression(pm.model,-1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mp"] * dP_iq[2]) - JuMP.@NLexpression(pm.model,-1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mp"] * dP_iq[3]) - ] - A[iBus,iBus][3,:] = [JuMP.@NLexpression(pm.model,-1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mq"] * dQ_delta) - 0 - JuMP.@NLexpression(pm.model,-1 / mpData["bus"][iBus]["tau"] * (1 + mpData["bus"][iBus]["mq"] * dQ_v)) - JuMP.@NLexpression(pm.model,-1 / mpData["bus"][iBus]["tau"] * (mpData["bus"][iBus]["mq"] * dQ_id[1])) - JuMP.@NLexpression(pm.model,-1 / mpData["bus"][iBus]["tau"] * (mpData["bus"][iBus]["mq"] * dQ_id[2])) - JuMP.@NLexpression(pm.model,-1 / mpData["bus"][iBus]["tau"] * (mpData["bus"][iBus]["mp"] * dQ_id[3])) - JuMP.@NLexpression(pm.model,-1 / mpData["bus"][iBus]["tau"] * (mpData["bus"][iBus]["mq"] * dQ_iq[1])) - JuMP.@NLexpression(pm.model,-1 / mpData["bus"][iBus]["tau"] * (mpData["bus"][iBus]["mq"] * dQ_iq[2])) - JuMP.@NLexpression(pm.model,-1 / mpData["bus"][iBus]["tau"] * (mpData["bus"][iBus]["mp"] * dQ_iq[3])) - ] - A[iBus,iBus][4:6,1] = PMS.vectorMulti(pm.model, LA.inv(L), dvd_delta) - A[iBus,iBus][4:6,3] = PMS.vectorMulti(pm.model, LA.inv(L), dvd_v) - A[iBus,iBus][4:6,4:6] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) - A[iBus,iBus][4:6,7:9] = LA.inv(L) * link["br_x"] - A[iBus,iBus][7:9,1] = PMS.vectorMulti(pm.model, LA.inv(L), dvq_delta) - A[iBus,iBus][7:9,3] = PMS.vectorMulti(pm.model, LA.inv(L), dvq_v) - A[iBus,iBus][7:9,4:6] = -LA.inv(L) * link["br_x"] - A[iBus,iBus][7:9,7:9] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) + dQ_iq = [JuMP.@expression(pm.model, -vd0[i]) for i in 1:length(vd0)] + + A[iBus, iBus][1, 2] = 1 + A[iBus, iBus][2, :] = [JuMP.@expression(pm.model, -1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mp"] * dP_delta) + -1 / mpData["bus"][iBus]["tau"] + JuMP.@expression(pm.model, -1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mp"] * dP_v) + JuMP.@expression(pm.model, -1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mp"] * dP_id[1]) + JuMP.@expression(pm.model, -1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mp"] * dP_id[2]) + JuMP.@expression(pm.model, -1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mp"] * dP_id[3]) + JuMP.@expression(pm.model, -1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mp"] * dP_iq[1]) + JuMP.@expression(pm.model, -1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mp"] * dP_iq[2]) + JuMP.@expression(pm.model, -1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mp"] * dP_iq[3]) + ] + A[iBus, iBus][3, :] = [JuMP.@expression(pm.model, -1 / mpData["bus"][iBus]["tau"] * mpData["bus"][iBus]["mq"] * dQ_delta) + 0 + JuMP.@expression(pm.model, -1 / mpData["bus"][iBus]["tau"] * (1 + mpData["bus"][iBus]["mq"] * dQ_v)) + JuMP.@expression(pm.model, -1 / mpData["bus"][iBus]["tau"] * (mpData["bus"][iBus]["mq"] * dQ_id[1])) + JuMP.@expression(pm.model, -1 / mpData["bus"][iBus]["tau"] * (mpData["bus"][iBus]["mq"] * dQ_id[2])) + JuMP.@expression(pm.model, -1 / mpData["bus"][iBus]["tau"] * (mpData["bus"][iBus]["mp"] * dQ_id[3])) + JuMP.@expression(pm.model, -1 / mpData["bus"][iBus]["tau"] * (mpData["bus"][iBus]["mq"] * dQ_iq[1])) + JuMP.@expression(pm.model, -1 / mpData["bus"][iBus]["tau"] * (mpData["bus"][iBus]["mq"] * dQ_iq[2])) + JuMP.@expression(pm.model, -1 / mpData["bus"][iBus]["tau"] * (mpData["bus"][iBus]["mp"] * dQ_iq[3])) + ] + A[iBus, iBus][4:6, 1] = PMS.vectorMulti(pm.model, LA.inv(L), dvd_delta) + A[iBus, iBus][4:6, 3] = PMS.vectorMulti(pm.model, LA.inv(L), dvd_v) + A[iBus, iBus][4:6, 4:6] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) + A[iBus, iBus][4:6, 7:9] = LA.inv(L) * link["br_x"] + A[iBus, iBus][7:9, 1] = PMS.vectorMulti(pm.model, LA.inv(L), dvq_delta) + A[iBus, iBus][7:9, 3] = PMS.vectorMulti(pm.model, LA.inv(L), dvq_v) + A[iBus, iBus][7:9, 4:6] = -LA.inv(L) * link["br_x"] + A[iBus, iBus][7:9, 7:9] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) end # for every pair of inverters in the same bus for iBus in invList if length(invConnected[iBus]) > 1 for i in 1:length(invConnected[iBus]) - for j in (i + 1):length(invConnected[iBus]) - A[i,j] = zeros(9, 9) + for j in (i+1):length(invConnected[iBus]) + A[i, j] = zeros(9, 9) link = mpData["branch"][invLine[i]] L = link["br_x"] ./ omega0 - A[i,j][4:6,4:6] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) - A[i,j][7:9,7:9] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) + A[i, j][4:6, 4:6] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) + A[i, j][7:9, 7:9] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) link = mpData["branch"][invLine[j]] L = link["br_x"] ./ omega0 - A[j,i] = zeros(9, 9) - A[j,i][4:6,4:6] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) - A[j,i][7:9,7:9] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) + A[j, i] = zeros(9, 9) + A[j, i][4:6, 4:6] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) + A[j, i][7:9, 7:9] = -LA.inv(L) * link["br_r"] - rN * LA.inv(L) end end end @@ -240,16 +240,16 @@ function obtainB_inverter_global(mpData, rN, omega0, busList, brList, invList, i iBus = invBusDict[i] for k in brList # initialize the matrix to be zero matrices - B[i,k] = zeros(9, 6) + B[i, k] = zeros(9, 6) iInd = parse(Int64, iBus) L = mpData["branch"][invLine[i]]["br_x"] ./ omega0 if mpData["branch"][k]["f_bus"] == iInd - B[i,k][4:6,1:3] = rN * PMS.inverseMat(L, mpData["branch"][invLine[i]]["f_connections"]) - B[i,k][7:9,4:6] = rN * PMS.inverseMat(L, mpData["branch"][invLine[i]]["f_connections"]) + B[i, k][4:6, 1:3] = rN * PMS.inverseMat(L, mpData["branch"][invLine[i]]["f_connections"]) + B[i, k][7:9, 4:6] = rN * PMS.inverseMat(L, mpData["branch"][invLine[i]]["f_connections"]) elseif mpData["branch"][k]["t_bus"] == iInd - B[i,k][4:6,1:3] = -rN * PMS.inverseMat(L, mpData["branch"][invLine[i]]["t_connections"]) - B[i,k][7:9,4:6] = -rN * PMS.inverseMat(L, mpData["branch"][invLine[i]]["t_connections"]) + B[i, k][4:6, 1:3] = -rN * PMS.inverseMat(L, mpData["branch"][invLine[i]]["t_connections"]) + B[i, k][7:9, 4:6] = -rN * PMS.inverseMat(L, mpData["branch"][invLine[i]]["t_connections"]) end end @@ -260,26 +260,26 @@ end "Obtain the submatrix C, with the opf solution given" function obtainC_inverter_global(mpData, rN, omega0, busList, brList, invList, invLine, invConnected, load_L, loadConnections, inverters, invBusDict) - + C = Dict() for i in inverters # initialize the matrix to be zero matrices - C[i,invBusDict[i]] = zeros(9, 6) + C[i, invBusDict[i]] = zeros(9, 6) iInd = parse(Int64, i) if invBusDict[i] in keys(load_L) L = load_L[invBusDict[i]] - C[i,invBusDict[i]][4:6,1:3] = rN * PMS.inverseMat(L, loadConnections[invBusDict[i]]) - C[i,invBusDict[i]][7:9,4:6] = rN * PMS.inverseMat(L, loadConnections[invBusDict[i]]) + C[i, invBusDict[i]][4:6, 1:3] = rN * PMS.inverseMat(L, loadConnections[invBusDict[i]]) + C[i, invBusDict[i]][7:9, 4:6] = rN * PMS.inverseMat(L, loadConnections[invBusDict[i]]) end end return C end - "Obtain the submatrix D, with the opf solution given" +"Obtain the submatrix D, with the opf solution given" function obtainD_inverter_global(mpData, rN, omega0, busList, brList, invList, invLine, inverters, invBusDict) - + D = Dict() for k in brList @@ -287,14 +287,14 @@ function obtainD_inverter_global(mpData, rN, omega0, busList, brList, invList, i for i in inverters iBus = invBusDict[i] # it is a regular bus - D[k,i] = zeros(6, 9) + D[k, i] = zeros(6, 9) if iBus == "$(mpData["branch"][k]["f_bus"])" - D[k,i][1:3,4:6] = rN * PMS.inverseMat(L, mpData["branch"][k]["f_connections"]) - D[k,i][4:6,7:9] = rN * PMS.inverseMat(L, mpData["branch"][k]["f_connections"]) + D[k, i][1:3, 4:6] = rN * PMS.inverseMat(L, mpData["branch"][k]["f_connections"]) + D[k, i][4:6, 7:9] = rN * PMS.inverseMat(L, mpData["branch"][k]["f_connections"]) end if iBus == "$(mpData["branch"][k]["t_bus"])" - D[k,i][1:3,4:6] = -rN * PMS.inverseMat(L, mpData["branch"][k]["t_connections"]) - D[k,i][4:6,7:9] = -rN * PMS.inverseMat(L, mpData["branch"][k]["t_connections"]) + D[k, i][1:3, 4:6] = -rN * PMS.inverseMat(L, mpData["branch"][k]["t_connections"]) + D[k, i][4:6, 7:9] = -rN * PMS.inverseMat(L, mpData["branch"][k]["t_connections"]) end end end @@ -309,30 +309,30 @@ function obtainE_inverter_global(mpData, rN, omega0, brList) for k1 in brList L = mpData["branch"][k1]["br_x"] ./ omega0 for k2 in brList - E[k1,k2] = zeros(6, 6) + E[k1, k2] = zeros(6, 6) if k1 == k2 - E[k1,k2][1:3,1:3] = -2 * rN * PMS.inverseMat(L, mpData["branch"][k1]["f_connections"]) - - PMS.inverseMat(L * mpData["branch"][k1]["br_r"], mpData["branch"][k1]["f_connections"]) - E[k1,k2][1:3,4:6] = PMS.inverseMat(L * mpData["branch"][k1]["br_x"], mpData["branch"][k1]["f_connections"]) - E[k1,k2][4:6,1:3] = -PMS.inverseMat(L * mpData["branch"][k1]["br_x"], mpData["branch"][k1]["f_connections"]) - E[k1,k2][4:6,4:6] = -2 * rN * PMS.inverseMat(L, mpData["branch"][k1]["f_connections"]) - - PMS.inverseMat(L * mpData["branch"][k1]["br_r"], mpData["branch"][k1]["f_connections"]) + E[k1, k2][1:3, 1:3] = -2 * rN * PMS.inverseMat(L, mpData["branch"][k1]["f_connections"]) - + PMS.inverseMat(L * mpData["branch"][k1]["br_r"], mpData["branch"][k1]["f_connections"]) + E[k1, k2][1:3, 4:6] = PMS.inverseMat(L * mpData["branch"][k1]["br_x"], mpData["branch"][k1]["f_connections"]) + E[k1, k2][4:6, 1:3] = -PMS.inverseMat(L * mpData["branch"][k1]["br_x"], mpData["branch"][k1]["f_connections"]) + E[k1, k2][4:6, 4:6] = -2 * rN * PMS.inverseMat(L, mpData["branch"][k1]["f_connections"]) - + PMS.inverseMat(L * mpData["branch"][k1]["br_r"], mpData["branch"][k1]["f_connections"]) elseif (mpData["branch"][k1]["f_bus"] == mpData["branch"][k2]["f_bus"]) - E[k1,k2][1:3,1:3] = -rN * PMS.inverseMat(L, mpData["branch"][k1]["f_connections"]) - E[k1,k2][4:6,4:6] = -rN * PMS.inverseMat(L, mpData["branch"][k1]["f_connections"]) + E[k1, k2][1:3, 1:3] = -rN * PMS.inverseMat(L, mpData["branch"][k1]["f_connections"]) + E[k1, k2][4:6, 4:6] = -rN * PMS.inverseMat(L, mpData["branch"][k1]["f_connections"]) elseif (mpData["branch"][k1]["t_bus"] == mpData["branch"][k2]["t_bus"]) - E[k1,k2][1:3,1:3] = -rN * PMS.inverseMat(L, mpData["branch"][k1]["t_connections"]) - E[k1,k2][4:6,4:6] = -rN * PMS.inverseMat(L, mpData["branch"][k1]["t_connections"]) + E[k1, k2][1:3, 1:3] = -rN * PMS.inverseMat(L, mpData["branch"][k1]["t_connections"]) + E[k1, k2][4:6, 4:6] = -rN * PMS.inverseMat(L, mpData["branch"][k1]["t_connections"]) elseif mpData["branch"][k1]["f_bus"] == mpData["branch"][k2]["t_bus"] if mpData["branch"][k1]["t_bus"] == mpData["branch"][k2]["f_bus"] - E[k1,k2][1:3,1:3] = 2 * rN * PMS.inverseMat(L, mpData["branch"][k1]["f_connections"]) - E[k1,k2][4:6,4:6] = 2 * rN * PMS.inverseMat(L, mpData["branch"][k1]["f_connections"]) + E[k1, k2][1:3, 1:3] = 2 * rN * PMS.inverseMat(L, mpData["branch"][k1]["f_connections"]) + E[k1, k2][4:6, 4:6] = 2 * rN * PMS.inverseMat(L, mpData["branch"][k1]["f_connections"]) else - E[k1,k2][1:3,1:3] = rN * PMS.inverseMat(L, mpData["branch"][k1]["f_connections"]) - E[k1,k2][4:6,4:6] = rN * PMS.inverseMat(L, mpData["branch"][k1]["f_connections"]) + E[k1, k2][1:3, 1:3] = rN * PMS.inverseMat(L, mpData["branch"][k1]["f_connections"]) + E[k1, k2][4:6, 4:6] = rN * PMS.inverseMat(L, mpData["branch"][k1]["f_connections"]) end end end @@ -343,21 +343,21 @@ end "Obtain the submatrix F, with the opf solution given" function obtainF_inverter_global(mpData, rN, omega0, busList, brList, loadList, loadConnections) - + F = Dict() for k in brList L = mpData["branch"][k]["br_x"] ./ omega0 for i in busList - F[k,i] = zeros(6, 6) + F[k, i] = zeros(6, 6) # obtain line L if i in keys(loadList) if mpData["branch"][k]["f_bus"] == parse(Int64, i) - F[k,i][1:3,1:3] = -rN * PMS.inverseMat(L, mpData["branch"][k]["f_connections"]) - F[k,i][4:6,4:6] = -rN * PMS.inverseMat(L, mpData["branch"][k]["f_connections"]) + F[k, i][1:3, 1:3] = -rN * PMS.inverseMat(L, mpData["branch"][k]["f_connections"]) + F[k, i][4:6, 4:6] = -rN * PMS.inverseMat(L, mpData["branch"][k]["f_connections"]) elseif mpData["branch"][k]["t_bus"] == parse(Int64, i) - F[k,i][1:3,1:3] = rN * PMS.inverseMat(L, mpData["branch"][k]["t_connections"]) - F[k,i][4:6,4:6] = rN * PMS.inverseMat(L, mpData["branch"][k]["t_connections"]) + F[k, i][1:3, 1:3] = rN * PMS.inverseMat(L, mpData["branch"][k]["t_connections"]) + F[k, i][4:6, 4:6] = rN * PMS.inverseMat(L, mpData["branch"][k]["t_connections"]) end end end @@ -372,11 +372,11 @@ function obtainG_inverter_global(mpData, rN, omega0, busList, brList, invList, l for i in inverters iBus = invBusDict[i] - G[iBus,i] = zeros(6, 9) + G[iBus, i] = zeros(6, 9) if iBus in keys(loadList) L = load_L[iBus] - G[iBus,i][1:3,4:6] = rN * PMS.inverseMat(L[loadConnections[iBus],loadConnections[iBus]], loadConnections[iBus]) - G[iBus,i][4:6,7:9] = rN * PMS.inverseMat(L[loadConnections[iBus],loadConnections[iBus]], loadConnections[iBus]) + G[iBus, i][1:3, 4:6] = rN * PMS.inverseMat(L[loadConnections[iBus], loadConnections[iBus]], loadConnections[iBus]) + G[iBus, i][4:6, 7:9] = rN * PMS.inverseMat(L[loadConnections[iBus], loadConnections[iBus]], loadConnections[iBus]) end end @@ -389,18 +389,18 @@ function obtainH_inverter_global(mpData, rN, omega0, busList, brList, load_L, lo for i in busList for k in brList - H[i,k] = zeros(6, 6) + H[i, k] = zeros(6, 6) # obtain line L if i in keys(load_L) L = load_L[i] try # if L != 0 if mpData["branch"][k]["f_bus"] == parse(Int64, i) - H[i,k][1:3,1:3] = -rN * PMS.inverseMat(L[loadConnections[i],loadConnections[i]], loadConnections[i]) - H[i,k][4:6,4:6] = -rN * PMS.inverseMat(L[loadConnections[i],loadConnections[i]], loadConnections[i]) + H[i, k][1:3, 1:3] = -rN * PMS.inverseMat(L[loadConnections[i], loadConnections[i]], loadConnections[i]) + H[i, k][4:6, 4:6] = -rN * PMS.inverseMat(L[loadConnections[i], loadConnections[i]], loadConnections[i]) elseif mpData["branch"][k]["t_bus"] == parse(Int64, i) - H[i,k][1:3,1:3] = rN * PMS.inverseMat(L[loadConnections[i],loadConnections[i]], loadConnections[i]) - H[i,k][4:6,4:6] = rN * PMS.inverseMat(L[loadConnections[i],loadConnections[i]], loadConnections[i]) + H[i, k][1:3, 1:3] = rN * PMS.inverseMat(L[loadConnections[i], loadConnections[i]], loadConnections[i]) + H[i, k][4:6, 4:6] = rN * PMS.inverseMat(L[loadConnections[i], loadConnections[i]], loadConnections[i]) end catch Memento.warn(_LOGGER, "L[$(i)] is not invertable! for branch $(k)") @@ -421,12 +421,12 @@ function obtainI_inverter_global(mpData, rN, omega0, busList, brList, loadList, if i in keys(loadList) L = load_L[i] try - I[i][1:3,1:3] = -rN * PMS.inverseMat(L[loadConnections[i],loadConnections[i]], loadConnections[i]) - - PMS.inverseMat(L[loadConnections[i],loadConnections[i]], loadConnections[i]) * load_R[i] - I[i][1:3,4:6] = PMS.inverseMat(L[loadConnections[i],loadConnections[i]], loadConnections[i]) * load_X[i] - I[i][4:6,1:3] = -PMS.inverseMat(L[loadConnections[i],loadConnections[i]], loadConnections[i]) * load_X[i] - I[i][1:3,4:6] = -rN * PMS.inverseMat(L[loadConnections[i],loadConnections[i]], loadConnections[i]) - - PMS.inverseMat(L[loadConnections[i],loadConnections[i]], loadConnections[i]) * load_R[i] + I[i][1:3, 1:3] = -rN * PMS.inverseMat(L[loadConnections[i], loadConnections[i]], loadConnections[i]) - + PMS.inverseMat(L[loadConnections[i], loadConnections[i]], loadConnections[i]) * load_R[i] + I[i][1:3, 4:6] = PMS.inverseMat(L[loadConnections[i], loadConnections[i]], loadConnections[i]) * load_X[i] + I[i][4:6, 1:3] = -PMS.inverseMat(L[loadConnections[i], loadConnections[i]], loadConnections[i]) * load_X[i] + I[i][1:3, 4:6] = -rN * PMS.inverseMat(L[loadConnections[i], loadConnections[i]], loadConnections[i]) - + PMS.inverseMat(L[loadConnections[i], loadConnections[i]], loadConnections[i]) * load_R[i] catch Memento.warn(_LOGGER, "L[$(i)] is not invertable!") end @@ -438,14 +438,14 @@ end "Combine all submatrices to one stability control matrix" function combineSub(busList, brList, inverters, invBusDict, A, B, C, D, E, F, G, H, I, type) - - n = length(busList) + + n = length(busList) ni = length(inverters) - m = length(brList) + m = length(brList) if type == 1 Atot = zeros(9 * ni + 6 * m + 6 * n, 9 * ni + 6 * m + 6 * n) - else + else Atot = Array{Any,2}(zeros(9 * ni + 6 * m + 6 * n, 9 * ni + 6 * m + 6 * n)) end @@ -454,41 +454,41 @@ function combineSub(busList, brList, inverters, invBusDict, A, B, C, D, E, F, G, # fill in A for j in 1:ni if (inverters[i], inverters[j]) in keys(A) - Atot[9 * (i - 1) + 1:9 * i,9 * (j - 1) + 1:9 * j] = A[inverters[i],inverters[j]] + Atot[9*(i-1)+1:9*i, 9*(j-1)+1:9*j] = A[inverters[i], inverters[j]] end end for j in 1:m # fill in B - Atot[9 * (i - 1) + 1:9 * i,9 * ni + 6 * (j - 1) + 1:9 * ni + 6 * j] = B[inverters[i],brList[j]] + Atot[9*(i-1)+1:9*i, 9*ni+6*(j-1)+1:9*ni+6*j] = B[inverters[i], brList[j]] # fill in D - Atot[9 * ni + 6 * (j - 1) + 1:9 * ni + 6 * j,9 * (i - 1) + 1:9 * i] = D[brList[j],inverters[i]] + Atot[9*ni+6*(j-1)+1:9*ni+6*j, 9*(i-1)+1:9*i] = D[brList[j], inverters[i]] end for j in 1:n if busList[j] == invBusDict[inverters[i]] # fill in C - Atot[9 * (i - 1) + 1:9 * i,9 * ni + 6 * m + 6 * (j - 1) + 1:9 * ni + 6 * m + 6 * j] = C[inverters[i],busList[j]] + Atot[9*(i-1)+1:9*i, 9*ni+6*m+6*(j-1)+1:9*ni+6*m+6*j] = C[inverters[i], busList[j]] # fill in G - Atot[9 * ni + 6 * m + 6 * (j - 1) + 1:9 * ni + 6 * m + 6 * j,9 * (i - 1) + 1:9 * i] = G[busList[j],inverters[i]] + Atot[9*ni+6*m+6*(j-1)+1:9*ni+6*m+6*j, 9*(i-1)+1:9*i] = G[busList[j], inverters[i]] end end end for i in 1:n # fill in I - Atot[9 * ni + 6 * m + 6 * (i - 1) + 1:9 * ni + 6 * m + 6 * i,9 * ni + 6 * m + 6 * (i - 1) + 1:9 * ni + 6 * m + 6 * i] = I[busList[i]] + Atot[9*ni+6*m+6*(i-1)+1:9*ni+6*m+6*i, 9*ni+6*m+6*(i-1)+1:9*ni+6*m+6*i] = I[busList[i]] for j in 1:m # fill in F - Atot[9 * ni + 6 * (j - 1) + 1:9 * ni + 6 * j,9 * ni + 6 * m + 6 * (i - 1) + 1:9 * ni + 6 * m + 6 * i] = F[brList[j],busList[i]] + Atot[9*ni+6*(j-1)+1:9*ni+6*j, 9*ni+6*m+6*(i-1)+1:9*ni+6*m+6*i] = F[brList[j], busList[i]] # fill in H - Atot[9 * ni + 6 * m + 6 * (i - 1) + 1:9 * ni + 6 * m + 6 * i,9 * ni + 6 * (j - 1) + 1:9 * ni + 6 * j] = H[busList[i],brList[j]] + Atot[9*ni+6*m+6*(i-1)+1:9*ni+6*m+6*i, 9*ni+6*(j-1)+1:9*ni+6*j] = H[busList[i], brList[j]] end end for j1 in 1:m for j2 in 1:m - # fill in E - Atot[9 * ni + 6 * (j1 - 1) + 1:9 * ni + 6 * j1,9 * ni + 6 * (j2 - 1) + 1:9 * ni + 6 * j2] = E[brList[j1],brList[j2]] + # fill in E + Atot[9*ni+6*(j1-1)+1:9*ni+6*j1, 9*ni+6*(j2-1)+1:9*ni+6*j2] = E[brList[j1], brList[j2]] end end @@ -498,7 +498,7 @@ end "obtain the global matrix of the small-signal stability control" function get_global_stability_matrix(mpData, opfSol, rN, omega0) - if isapprox(omega0, 0, atol=1E-6) + if isapprox(omega0, 0, atol=1E-6) Memento.error(_LOGGER, "Inverter frequency set-point value cannot be zero") end diff --git a/src/io/preprocessing.jl b/src/io/preprocessing.jl index fc1d3ab..a8dc001 100644 --- a/src/io/preprocessing.jl +++ b/src/io/preprocessing.jl @@ -73,7 +73,7 @@ end "Obtain the load parameters from the model data" function get_load_parameters(mpData, loadList, vnomList, omega0, loadConnections) - + θ_list = [0, 2π/3, -2π/3] load_R = Dict() load_X = Dict() diff --git a/test/data/case2_inverters.json b/test/data/case2_inverters.json index 999c06f..97e0129 100644 --- a/test/data/case2_inverters.json +++ b/test/data/case2_inverters.json @@ -10,50 +10,22 @@ "tau": 0.0001, "inverter_bus": true, "r": [ - [0.001,0,0], - [0,0.001,0], - [0,0,0.001] + [0.001, 0, 0], + [0, 0.001, 0], + [0, 0, 0.001] ], "x": [ - [0.001,0,0], - [0,0.001,0], - [0,0,0.001] - ], - "pg": [ - 33.3333, - 33.3333, - 33.3333 - ], - "qg": [ - 33.3333, - 33.3333, - 33.3333 - ], - "pg_lb": [ - 0.0, - 0.0, - 0.0 - ], - "qg_lb": [ - 0.0, - 0.0, - 0.0 - ], - "pg_ub": [ - 100.0, - 100.0, - 100.0 - ], - "qg_ub": [ - 0.0, - 0.0, - 0.0 - ], - "vg": [ - 4.6, - 4.6, - 4.6 - ] - } - ] + [0.001, 0, 0], + [0, 0.001, 0], + [0, 0, 0.001] + ], + "pg": [33.3333, 33.3333, 33.3333], + "qg": [33.3333, 33.3333, 33.3333], + "pg_lb": [0.0, 0.0, 0.0], + "qg_lb": [0.0, 0.0, 0.0], + "pg_ub": [100.0, 100.0, 100.0], + "qg_ub": [0.0, 0.0, 0.0], + "vg": [4.6, 4.6, 4.6] + } + ] } diff --git a/test/runtests.jl b/test/runtests.jl index f1ad6b8..53b9c27 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,8 +10,8 @@ const PMS = PowerModelsStability PMD.silence!() -ipopt_solver = PMD.optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 3, "tol"=>1e-5, "max_iter" => Int(1E4), "sb" => "yes") +ipopt_solver = PMD.optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0, "tol"=>1e-5, "max_iter" => Int(1E4), "sb" => "yes") @testset "PowerModelsStability" begin include("test_2bus.jl") -end \ No newline at end of file +end diff --git a/test/test_2bus.jl b/test/test_2bus.jl index 7532e86..27d600f 100644 --- a/test/test_2bus.jl +++ b/test/test_2bus.jl @@ -1,5 +1,5 @@ @testset "test the two bus system" begin - + println(">>>>>>> Testing 2-bus system <<<<<<<<") # load the data