From 9d3b09904fb6eae73332bddce27c554d96e78140 Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Wed, 8 Jan 2025 15:29:14 -0700 Subject: [PATCH] NR PF non-allocating --- src/PowerFlows.jl | 1 + src/newton_ac_powerflow.jl | 517 +++++++++++++++++++----- src/powerflow_types.jl | 3 +- test/Project.toml | 1 + test/runtests.jl | 2 + test/test_newton_ac_powerflow.jl | 226 ++++++++++- test/test_utils/psse_results_compare.jl | 6 +- 7 files changed, 633 insertions(+), 123 deletions(-) diff --git a/src/PowerFlows.jl b/src/PowerFlows.jl index 00907a8b..51289003 100644 --- a/src/PowerFlows.jl +++ b/src/PowerFlows.jl @@ -6,6 +6,7 @@ export PowerFlowData export DCPowerFlow export NLSolveACPowerFlow export KLUACPowerFlow +export LUACPowerFlow export ACPowerFlow export ACPowerFlowSolverType export PTDFDCPowerFlow diff --git a/src/newton_ac_powerflow.jl b/src/newton_ac_powerflow.jl index 014384b3..b47472c8 100644 --- a/src/newton_ac_powerflow.jl +++ b/src/newton_ac_powerflow.jl @@ -167,6 +167,224 @@ function _newton_powerflow( return res.f_converged, res.zero end +function _update_V!(dx::Vector{Float64}, V::Vector{Complex{Float64}}, Vm::Vector{Float64}, + Va::Vector{Float64}, + pv::Vector{Int64}, pq::Vector{Int64}, dx_Va_pv::Vector{Int64}, dx_Va_pq::Vector{Int64}, + dx_Vm_pq::Vector{Int64}) + for (i, j) in zip(pv, dx_Va_pv) + Va[i] -= dx[j] + end + + for (i, j) in zip(pq, dx_Va_pq) + Va[i] -= dx[j] + end + + for (i, j) in zip(pq, dx_Vm_pq) + Vm[i] -= dx[j] + end + + V .= Vm .* exp.(1im .* Va) + + Vm .= abs.(V) + Va .= angle.(V) + return +end + +function _update_F!(F::Vector{Float64}, mis::Vector{Complex{Float64}}, + dx_Va_pv::Vector{Int64}, dx_Va_pq::Vector{Int64}, dx_Vm_pq::Vector{Int64}, + V::Vector{Complex{Float64}}, Ybus::SparseMatrixCSC{Complex{Float64}, Int64}, + Sbus::Vector{Complex{Float64}}, + pv::Vector{Int64}, pq::Vector{Int64}) + + #mis .= V .* conj.(Ybus * V) .- Sbus + LinearAlgebra.mul!(mis, Ybus, V) + mis .= V .* conj.(mis) .- Sbus + + for (i, j) in zip(dx_Va_pv, pv) + F[i] = real(mis[j]) + end + + for (i, j) in zip(dx_Va_pq, pq) + F[i] = real(mis[j]) + end + + for (i, j) in zip(dx_Vm_pq, pq) + F[i] = imag(mis[j]) + end + return +end + +function _update_dSbus_dV!(rows::Vector{Int64}, cols::Vector{Int64}, + V::Vector{Complex{Float64}}, Ybus::SparseMatrixCSC{Complex{Float64}, Int64}, + diagV::LinearAlgebra.Diagonal{Complex{Float64}, Vector{Complex{Float64}}}, + diagVnorm::LinearAlgebra.Diagonal{Complex{Float64}, Vector{Complex{Float64}}}, + diagIbus::LinearAlgebra.Diagonal{Complex{Float64}, Vector{Complex{Float64}}}, + diagIbus_diag::Vector{Complex{Float64}}, + dSbus_dVa::SparseMatrixCSC{Complex{Float64}, Int64}, + dSbus_dVm::SparseMatrixCSC{Complex{Float64}, Int64}, + r_dSbus_dVa::SparseMatrixCSC{Float64, Int64}, + r_dSbus_dVm::SparseMatrixCSC{Float64, Int64}, + i_dSbus_dVa::SparseMatrixCSC{Float64, Int64}, + i_dSbus_dVm::SparseMatrixCSC{Float64, Int64}, + Ybus_diagVnorm::SparseMatrixCSC{Complex{Float64}, Int64}, + conj_Ybus_diagVnorm::SparseMatrixCSC{Complex{Float64}, Int64}, + diagV_conj_Ybus_diagVnorm::SparseMatrixCSC{Complex{Float64}, Int64}, + conj_diagIbus::LinearAlgebra.Diagonal{Complex{Float64}, Vector{Complex{Float64}}}, + conj_diagIbus_diagVnorm::LinearAlgebra.Diagonal{ + Complex{Float64}, + Vector{Complex{Float64}}, + }, + Ybus_diagV::SparseMatrixCSC{Complex{Float64}, Int64}, + conj_Ybus_diagV::SparseMatrixCSC{Complex{Float64}, Int64}) + for i in eachindex(V) + diagV[i, i] = V[i] + diagVnorm[i, i] = V[i] / abs(V[i]) + end + + # manually calculate the diagIbus matrix + LinearAlgebra.mul!(diagIbus_diag, Ybus, V) + for i in eachindex(V) + diagIbus[i, i] = diagIbus_diag[i] + end + + # use the available matrices temporarily to calculate the dSbus_dV matrices + # original formula: + # dSbus_dVm .= diagV * conj.(Ybus * diagVnorm) + conj.(diagIbus) * diagVnorm + # non-allocating version: + + LinearAlgebra.mul!(Ybus_diagVnorm, Ybus, diagVnorm) + conj_Ybus_diagVnorm .= conj.(Ybus_diagVnorm) + LinearAlgebra.mul!(diagV_conj_Ybus_diagVnorm, diagV, conj_Ybus_diagVnorm) + conj_diagIbus .= conj.(diagIbus) + LinearAlgebra.mul!(conj_diagIbus_diagVnorm, conj_diagIbus, diagVnorm) + + dSbus_dVm .= diagV_conj_Ybus_diagVnorm + LinearAlgebra.axpy!(1, conj_diagIbus_diagVnorm, dSbus_dVm) + + # original formula: + # dSbus_dVa .= 1im * diagV * conj.(diagIbus - Ybus * diagV) + # non-allocating version: + LinearAlgebra.mul!(Ybus_diagV, Ybus, diagV) + # Take the conjugate of the result (conj(Ybus * diagV)); conj_diagIbus is already available + conj_Ybus_diagV .= conj.(Ybus_diagV) + + # write the result of conj.(diagIbus - Ybus * diagV) in conj_Ybus_diagV + LinearAlgebra.axpby!(1, conj_diagIbus, -1, conj_Ybus_diagV) + + # Multiply the result by diagV + LinearAlgebra.mul!(dSbus_dVa, diagV, conj_Ybus_diagV) + + # Now multiply by 1im to get the final result + #dSbus_dVa .*= 1im + LinearAlgebra.mul!(dSbus_dVa, dSbus_dVa, 1im) + + for c in cols + for r in rows + r_dSbus_dVa[r, c] = real(dSbus_dVa[r, c]) + i_dSbus_dVa[r, c] = imag(dSbus_dVa[r, c]) + r_dSbus_dVm[r, c] = real(dSbus_dVm[r, c]) + i_dSbus_dVm[r, c] = imag(dSbus_dVm[r, c]) + end + end + + # sometimes can allocate so we have to use the for loop above + # r_dSbus_dVa .= real.(dSbus_dVa) + # r_dSbus_dVm .= real.(dSbus_dVm) + # i_dSbus_dVa .= imag.(dSbus_dVa) + # i_dSbus_dVm .= imag.(dSbus_dVm) + return +end + +# this function is for testing purposes only +function _legacy_dSbus_dV( + V::Vector{Complex{Float64}}, + Ybus::SparseMatrixCSC{Complex{Float64}, Int64}, +) + diagV = LinearAlgebra.Diagonal(V) + diagVnorm = LinearAlgebra.Diagonal(V ./ abs.(V)) + diagIbus = LinearAlgebra.Diagonal(Ybus * V) + dSbus_dVm = diagV * conj.(Ybus * diagVnorm) + conj.(diagIbus) * diagVnorm + dSbus_dVa = 1im * diagV * conj.(diagIbus - Ybus * diagV) + return dSbus_dVa, dSbus_dVm +end + +# this function is for testing purposes only +function _legacy_J( + dSbus_dVa::SparseMatrixCSC{Complex{Float64}, Int64}, + dSbus_dVm::SparseMatrixCSC{Complex{Float64}, Int64}, + pvpq::Vector{Int64}, + pq::Vector{Int64}, +) + j11 = real(dSbus_dVa[pvpq, pvpq]) + j12 = real(dSbus_dVm[pvpq, pq]) + j21 = imag(dSbus_dVa[pq, pvpq]) + j22 = imag(dSbus_dVm[pq, pq]) + J = sparse([j11 j12; j21 j22]) + return J +end + +function _update_submatrix!( + A::SparseMatrixCSC, + B::SparseMatrixCSC, + rows_A::Vector{Int64}, + cols_A::Vector{Int64}, + rows_B::Vector{Int64}, + cols_B::Vector{Int64}, +) + for idj in eachindex(cols_A) + for idi in eachindex(rows_A) + A[rows_A[idi], cols_A[idj]] = B[rows_B[idi], cols_B[idj]] + end + end + return +end + +function _update_J!(J::SparseMatrixCSC, + r_dSbus_dVa::SparseMatrixCSC, + r_dSbus_dVm::SparseMatrixCSC, + i_dSbus_dVa::SparseMatrixCSC, + i_dSbus_dVm::SparseMatrixCSC, + pvpq::Vector{Int64}, + pq::Vector{Int64}, + j_pvpq::Vector{Int64}, + j_pq::Vector{Int64}, +) + _update_submatrix!(J, r_dSbus_dVa, j_pvpq, j_pvpq, pvpq, pvpq) + _update_submatrix!(J, r_dSbus_dVm, j_pvpq, j_pq, pvpq, pq) + _update_submatrix!(J, i_dSbus_dVa, j_pq, j_pvpq, pq, pvpq) + _update_submatrix!(J, i_dSbus_dVm, j_pq, j_pq, pq, pq) + return +end + +function _calc_x( + data::ACPowerFlowData, + V::Vector{Complex{Float64}}, + Va::Vector{Float64}, + Vm::Vector{Float64}, + Ybus::SparseMatrixCSC{Complex{Float64}, Int64}, + n_buses::Int64, +) + # mock the expected x format, where the values depend on the type of the bus: + x = zeros(Float64, 2 * n_buses) + Sbus_result = V .* conj(Ybus * V) # todo preallocate and pass as parameter + for (ix, b) in enumerate(data.bus_type) + if b == PSY.ACBusTypes.REF + # When bustype == REFERENCE PSY.Bus, state variables are Active and Reactive Power Generated + x[2 * ix - 1] = real(Sbus_result[ix]) + data.bus_activepower_withdrawals[ix] + x[2 * ix] = imag(Sbus_result[ix]) + data.bus_reactivepower_withdrawals[ix] + elseif b == PSY.ACBusTypes.PV + # When bustype == PV PSY.Bus, state variables are Reactive Power Generated and Voltage Angle + x[2 * ix - 1] = imag(Sbus_result[ix]) + data.bus_reactivepower_withdrawals[ix] + x[2 * ix] = Va[ix] + elseif b == PSY.ACBusTypes.PQ + # When bustype == PQ PSY.Bus, state variables are Voltage Magnitude and Voltage Angle + x[2 * ix - 1] = Vm[ix] + x[2 * ix] = Va[ix] + end + end + return x +end + function _newton_powerflow( pf::ACPowerFlow{KLUACPowerFlow}, data::ACPowerFlowData; @@ -177,13 +395,15 @@ function _newton_powerflow( tol = get(nlsolve_kwargs, :tol, DEFAULT_NR_TOL) i = 0 + Ybus = data.power_network_matrix.data + # Find indices for each bus type ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type) pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type) pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type) pvpq = [pv; pq] - #nref = length(ref) + # nref = length(ref) npv = length(pv) npq = length(pq) npvpq = npv + npq @@ -191,118 +411,223 @@ function _newton_powerflow( Vm = data.bus_magnitude[:] # prevent unfeasible starting values for Vm; for pv and ref buses we cannot do this: - @. Vm[pq] = clamp.(Vm[pq], 0.9, 1.1) + Vm[pq] .= clamp.(Vm[pq], 0.9, 1.1) Va = data.bus_angles[:] V = zeros(Complex{Float64}, length(Vm)) - @. V = Vm .* exp.(1im * Va) + V .= Vm .* exp.(1im .* Va) - Va_pv = view(Va, pv) - Va_pq = view(Va, pq) - Vm_pq = view(Vm, pq) + # early return if only ref buses present - no need to solve the power flow + converged = npvpq == 0 + if converged + # if only ref buses present, we do not need to enter the power flow loop + x = _calc_x(data, V, Va, Vm, Ybus, n_buses) + return (converged, x) + end - # pre-allocate dx - dx = zeros(Float64, npv + 2 * npq) + # we need to define lookups for mappings of pv, pq buses onto the internal J indexing + pvpq_lookup = zeros(Int64, maximum([ref; pvpq]) + 1) + pvpq_lookup[pvpq] .= 1:npvpq + pq_lookup = zeros(Int64, maximum([ref; pvpq]) + 1) + pq_lookup[pq] .= 1:npq - dx_Va_pv = view(dx, 1:npv) - dx_Va_pq = view(dx, (npv + 1):(npv + npq)) - dx_Vm_pq = view(dx, (npv + npq + 1):(npv + 2 * npq)) + # define the internal J indexing using the lookup arrays + j_pvpq = pvpq_lookup[pvpq] + j_pq = npvpq .+ pq_lookup[pq] - Ybus = data.power_network_matrix.data + # indices for updating of V + dx_Va_pv = Vector{Int64}([1:npv...]) + dx_Va_pq = Vector{Int64}([(npv + 1):(npv + npq)...]) + dx_Vm_pq = Vector{Int64}([(npv + npq + 1):(npv + 2 * npq)...]) Sbus = data.bus_activepower_injection[:] - data.bus_activepower_withdrawals[:] + 1im * (data.bus_reactivepower_injection[:] - data.bus_reactivepower_withdrawals[:]) - + # Pre-allocate mis and F and create views for the respective real and imaginary sections of the arrays: mis = zeros(Complex{Float64}, length(V)) - mis_pvpq = view(mis, pvpq) - mis_pq = view(mis, pq) - F = zeros(Float64, npvpq + npq) - F_real = view(F, 1:npvpq) - F_imag = view(F, npvpq + 1:npvpq + npq) - mis .= V .* conj.(Ybus * V) .- Sbus - @. F_real = real(mis_pvpq) # In-place assignment to the real part, using views - @. F_imag = imag(mis_pq) # In-place assignment to the imaginary part, using views + F .= [real(mis[pvpq]); imag(mis[pq])] + + # preallocate Jacobian matrix and arrays for calculating dSbus_dVa, dSbus_dVm + rows, cols = SparseArrays.findnz(Ybus) + + #diagV = sparse(1:n_buses, 1:n_buses, V) + diagV = LinearAlgebra.Diagonal(V) + diagIbus_diag = zeros(Complex{Float64}, size(V, 1)) + diagIbus = LinearAlgebra.Diagonal(diagIbus_diag) + + diagVnorm = LinearAlgebra.Diagonal(V ./ abs.(V)) + + Ybus_diagVnorm = sparse(rows, cols, Complex{Float64}(0)) + conj_Ybus_diagVnorm = sparse(rows, cols, Complex{Float64}(0)) + diagV_conj_Ybus_diagVnorm = sparse(rows, cols, Complex{Float64}(0)) + conj_diagIbus = conj.(diagIbus) + conj_diagIbus_diagVnorm = conj.(diagIbus) + Ybus_diagV = sparse(rows, cols, Complex{Float64}(0)) + conj_Ybus_diagV = sparse(rows, cols, Complex{Float64}(0)) + + dSbus_dVm = sparse(rows, cols, Complex{Float64}(0)) + dSbus_dVa = sparse(rows, cols, Complex{Float64}(0)) + r_dSbus_dVa = sparse(rows, cols, Float64(0)) + r_dSbus_dVm = sparse(rows, cols, Float64(0)) + i_dSbus_dVa = sparse(rows, cols, Float64(0)) + i_dSbus_dVm = sparse(rows, cols, Float64(0)) + + # maybe use this in the future? + # pvpq_rows = pvpq_lookup[rows][pvpq_lookup[rows] .!= 0] + # pvpq_cols = pvpq_lookup[cols][pvpq_lookup[cols] .!= 0] + # pq_rows = pq_lookup[rows][pq_lookup[rows] .!= 0] + # pq_cols = pq_lookup[cols][pq_lookup[cols] .!= 0] + + J_block = sparse(rows, cols, Float64(0), maximum(rows), maximum(cols)) + J = [J_block[pvpq, pvpq] J_block[pvpq, pq]; J_block[pq, pvpq] J_block[pq, pq]] + + # preallocate the KLU factorization object - symbolic object only + colptr = KLU.decrement(J.colptr) + rowval = KLU.decrement(J.rowval) + n = size(J, 1) + factor_J = KLU.KLUFactorization(n, colptr, rowval, J.nzval) + KLU.klu_analyze!(factor_J) + rf = Ref(factor_J.common) + # factorization for the numeric object does not work here: + # factor_J._numeric = KLU.klu_l_factor(colptr, rowval, J.nzval, factor_J._symbolic, rf) - converged = npvpq == 0 # if only ref buses present, we do not need to enter the loop + while i < maxIter && !converged + i += 1 - # preallocate Jacobian matrix - rows = vcat(1:npvpq, 1:npvpq, npvpq+1:npvpq+npq, npvpq+1:npvpq+npq) - cols = vcat(1:npvpq, npvpq+1:npvpq+npq, 1:npvpq, npvpq+1:npvpq+npq) - J = sparse(rows, cols, Float64(0)) + _update_dSbus_dV!(rows, cols, V, Ybus, diagV, diagVnorm, diagIbus, diagIbus_diag, + dSbus_dVa, dSbus_dVm, r_dSbus_dVa, r_dSbus_dVm, i_dSbus_dVa, i_dSbus_dVm, + Ybus_diagVnorm, conj_Ybus_diagVnorm, diagV_conj_Ybus_diagVnorm, + conj_diagIbus, conj_diagIbus_diagVnorm, Ybus_diagV, conj_Ybus_diagV) + + # todo: improve pvpq, pq, j_pvpq, j_pq (use more specific indices) + _update_J!( + J, + r_dSbus_dVa, + r_dSbus_dVm, + i_dSbus_dVa, + i_dSbus_dVm, + pvpq, + pq, + j_pvpq, + j_pq, + ) - # we need to define lookups for mappings of pv, pq buses onto the internal J indexing - pvpq_lookup = zeros(Int64, maximum([ref; pvpq]) + 1) - pvpq_lookup[pvpq] .= 1:npvpq - pq_lookup = zeros(Int64, maximum([ref; pvpq]) + 1) - pq_lookup[pq] .= 1:npq + # Workaround for the issue with KLU.klu_l_factor + # background: KLU.klu_l_factor does not work properly with the preallocated J matrix with dummy values + # the workaround is to initialize the numeric object here in the loop once and then refactorize the matrix in the loop inplace + if i == 1 + # works when J values are ok: + factor_J._numeric = + KLU.klu_l_factor(colptr, rowval, J.nzval, factor_J._symbolic, rf) + end - # with the pre-allocated J and lookups, we can define views into the sub-matrices of the J matrix for updating the J matrix in the NR loop - j11 = view(J, pvpq_lookup[pvpq], pvpq_lookup[pvpq]) - j12 = view(J, pvpq_lookup[pvpq], npvpq .+ pq_lookup[pq]) - j21 = view(J, npvpq .+ pq_lookup[pq], pvpq_lookup[pvpq]) - j22 = view(J, npvpq .+ pq_lookup[pq], npvpq .+ pq_lookup[pq]) - - # we need views of the diagonals to avoid using LinearAlgebra.Diagonal: - diagV = sparse(1:n_buses, 1:n_buses, Complex{Float64}(1)) - diag_idx = LinearAlgebra.diagind(diagV) - diagV_diag = view(diagV, diag_idx) - - diagIbus = sparse(1:n_buses, 1:n_buses, Complex{Float64}(1)) - diagIbus_diag = view(diagIbus, diag_idx) - - diagVnorm = sparse(1:n_buses, 1:n_buses, Complex{Float64}(1)) - diagVnorm_diag = view(diagVnorm, diag_idx) - - # pre-allocate the dSbus_dVm, dSbus_dVa to have the same structure as Ybus - # they will follow the structure of Ybus except maybe when Ybus has zero values in its diagonal, which we do not expect here - #rows, cols, _ = SparseArrays.findnz(Ybus) - #dSbus_dVm = sparse(rows, cols, Complex{Float64}(0)) - #dSbus_dVa = sparse(rows, cols, Complex{Float64}(0)) - - # preallocate dSbus_dVm, dSbus_dVa with correct structure: - dSbus_dVm = diagV * conj.(Ybus * diagVnorm) + conj.(diagIbus) * diagVnorm - dSbus_dVa = 1im * diagV * conj.(diagIbus - Ybus * diagV) + # factorize the numeric object of KLU inplace, while reusing the symbolic object + KLU.klu_l_refactor( + colptr, + rowval, + J.nzval, + factor_J._symbolic, + factor_J._numeric, + rf, + ) - # create views for the sub-arrays of Sbus_dVa, Sbus_dVm for updating the J: - Sbus_dVa_j11 = view(dSbus_dVa, pvpq, pvpq) - Sbus_dVm_j12 = view(dSbus_dVm, pvpq, pq) - Sbus_dVa_j21 = view(dSbus_dVa, pq, pvpq) - Sbus_dVm_j22 = view(dSbus_dVm, pq, pq) + # solve inplace - the results are written to F, so that we must use F instead of dx for updating V + KLU.klu_l_solve( + factor_J._symbolic, + factor_J._numeric, + size(F, 1), + size(F, 2), + F, + rf, + ) - while i < maxIter && !converged - i += 1 - - ## use the new value of V to update dSbus_dVa, dSbus_dVm: - diagV_diag .= V - diagIbus_diag .= Ybus * V - @. diagVnorm_diag = V ./ abs.(V) - dSbus_dVm .= diagV * conj.(Ybus * diagVnorm) + conj.(diagIbus) * diagVnorm - dSbus_dVa .= 1im * diagV * conj.(diagIbus - Ybus * diagV) + # KLU.solve! overwrites F with the solution instead of returning it as dx, so -F is used here to update V + _update_V!(F, V, Vm, Va, pv, pq, dx_Va_pv, dx_Va_pq, dx_Vm_pq) + + # here F is mismatch again + _update_F!(F, mis, dx_Va_pv, dx_Va_pq, dx_Vm_pq, V, Ybus, Sbus, pv, pq) + + converged = LinearAlgebra.norm(F, Inf) < tol + end + + if !converged + @error("The powerflow solver with KLU did not converge after $i iterations") + else + @info("The powerflow solver with KLU converged after $i iterations") + end + + x = _calc_x(data, V, Va, Vm, Ybus, n_buses) + + return (converged, x) +end + +# legacy NR implementation - here we do not care about allocations, we use this function only for testing purposes +function _newton_powerflow( + pf::ACPowerFlow{LUACPowerFlow}, + data::ACPowerFlowData; + nlsolve_kwargs..., +) + # Fetch maxIter and tol from kwargs, or use defaults if not provided + maxIter = get(nlsolve_kwargs, :maxIter, DEFAULT_NR_MAX_ITER) + tol = get(nlsolve_kwargs, :tol, DEFAULT_NR_TOL) + i = 0 + + Ybus = data.power_network_matrix.data + + # Find indices for each bus type + #ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type) + pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type) + pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type) + pvpq = [pv; pq] - # update the Jacobian by setting values through the pre-defined views for j11, j12, j21, j22 - @. j11 = real(Sbus_dVa_j11) - @. j12 = real(Sbus_dVm_j12) - @. j21 = imag(Sbus_dVa_j21) - @. j22 = imag(Sbus_dVm_j22) + #nref = length(ref) + npv = length(pv) + npq = length(pq) + npvpq = npv + npq + n_buses = length(data.bus_type) - factor_J = KLU.klu(J) - dx .= -(factor_J \ F) + Vm = data.bus_magnitude[:] + # prevent unfeasible starting values for Vm; for pv and ref buses we cannot do this: + Vm[pq] .= clamp.(Vm[pq], 0.9, 1.1) + Va = data.bus_angles[:] + V = zeros(Complex{Float64}, length(Vm)) + V .= Vm .* exp.(1im .* Va) - Va_pv .+= dx_Va_pv - Va_pq .+= dx_Va_pq - Vm_pq .+= dx_Vm_pq + # pre-allocate dx + dx = zeros(Float64, npv + 2 * npq) - @. V = Vm .* exp.(1im * Va) + Ybus = data.power_network_matrix.data + + Sbus = + data.bus_activepower_injection[:] - data.bus_activepower_withdrawals[:] + + 1im * (data.bus_reactivepower_injection[:] - data.bus_reactivepower_withdrawals[:]) + + mis = V .* conj.(Ybus * V) .- Sbus + F = [real(mis[pvpq]); imag(mis[pq])] + + converged = npvpq == 0 + + while i < maxIter && !converged + i += 1 + dSbus_dVa, dSbus_dVm = _legacy_dSbus_dV(V, Ybus) + J = _legacy_J(dSbus_dVa, dSbus_dVm, pvpq, pq) + # using a different factorization that KLU for testing + factor_J = LinearAlgebra.lu(J) + dx .= factor_J \ F + + Va[pv] .-= dx[1:npv] + Va[pq] .-= dx[(npv + 1):(npv + npq)] + Vm[pq] .-= dx[(npv + npq + 1):(npv + 2 * npq)] + V .= Vm .* exp.(1im .* Va) Vm .= abs.(V) Va .= angle.(V) - mis .= V .* conj.(Ybus * V) .- Sbus - @. F_real = real(mis_pvpq) # In-place assignment to the real part - @. F_imag = imag(mis_pq) # In-place assignment to the imaginary part + mis = V .* conj.(Ybus * V) .- Sbus + F .= [real(mis[pvpq]); imag(mis[pq])] + converged = LinearAlgebra.norm(F, Inf) < tol end @@ -312,25 +637,7 @@ function _newton_powerflow( @info("The powerflow solver with KLU converged after $i iterations") end - # mock the expected x format, where the values depend on the type of the bus: - n_buses = length(data.bus_type) - x = zeros(Float64, 2 * n_buses) - Sbus_result = V .* conj(Ybus * V) - for (ix, b) in enumerate(data.bus_type) - if b == PSY.ACBusTypes.REF - # When bustype == REFERENCE PSY.Bus, state variables are Active and Reactive Power Generated - x[2 * ix - 1] = real(Sbus_result[ix]) + data.bus_activepower_withdrawals[ix] - x[2 * ix] = imag(Sbus_result[ix]) + data.bus_reactivepower_withdrawals[ix] - elseif b == PSY.ACBusTypes.PV - # When bustype == PV PSY.Bus, state variables are Reactive Power Generated and Voltage Angle - x[2 * ix - 1] = imag(Sbus_result[ix]) + data.bus_reactivepower_withdrawals[ix] - x[2 * ix] = Va[ix] - elseif b == PSY.ACBusTypes.PQ - # When bustype == PQ PSY.Bus, state variables are Voltage Magnitude and Voltage Angle - x[2 * ix - 1] = Vm[ix] - x[2 * ix] = Va[ix] - end - end + x = _calc_x(data, V, Va, Vm, Ybus, n_buses) - return converged, x + return (converged, x) end diff --git a/src/powerflow_types.jl b/src/powerflow_types.jl index e490e1f5..60d08c18 100644 --- a/src/powerflow_types.jl +++ b/src/powerflow_types.jl @@ -3,13 +3,14 @@ abstract type ACPowerFlowSolverType end struct KLUACPowerFlow <: ACPowerFlowSolverType end struct NLSolveACPowerFlow <: ACPowerFlowSolverType end +struct LUACPowerFlow <: ACPowerFlowSolverType end # Only for testing, a basic implementation using LinearAlgebra.lu, allocates a lot of memory Base.@kwdef struct ACPowerFlow{ACSolver <: ACPowerFlowSolverType} <: PowerFlowEvaluationModel check_reactive_power_limits::Bool = false end -# Create a constructor that defaults to KLUACPowerFlow +# Create a constructor for ACPowerFlow that defaults to KLUACPowerFlow function ACPowerFlow(ACSolver::Type{<:ACPowerFlowSolverType} = KLUACPowerFlow; check_reactive_power_limits::Bool = false, ) diff --git a/test/Project.toml b/test/Project.toml index 7ab27d53..f931fbfa 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -13,6 +13,7 @@ PowerFlows = "94fada2c-fd9a-4e89-8d82-81405f5cb4f6" PowerNetworkMatrices = "bed98974-b02a-5e2f-9fe0-a103f5c450dd" PowerSystemCaseBuilder = "f00506e0-b84f-492a-93c2-c0a9afc4364e" PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/runtests.jl b/test/runtests.jl index 3f8827fc..f798c10e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,6 +10,8 @@ using CSV using DataFrames using JSON3 using DataStructures +import SparseArrays +import SparseArrays: SparseMatrixCSC, sparse const IS = InfrastructureSystems const PSB = PowerSystemCaseBuilder diff --git a/test/test_newton_ac_powerflow.jl b/test/test_newton_ac_powerflow.jl index 7a0b31d8..4564a293 100644 --- a/test/test_newton_ac_powerflow.jl +++ b/test/test_newton_ac_powerflow.jl @@ -1,5 +1,9 @@ @testset "AC Power Flow 14-Bus testing" for ACSolver in - (NLSolveACPowerFlow, KLUACPowerFlow) + ( + NLSolveACPowerFlow, + KLUACPowerFlow, + LUACPowerFlow, +) result_14 = [ 2.3255081760423684 -0.15529254415401786 @@ -48,7 +52,11 @@ end @testset "AC Power Flow 14-Bus Line Configurations" for ACSolver in - (NLSolveACPowerFlow, KLUACPowerFlow) + ( + NLSolveACPowerFlow, + KLUACPowerFlow, + LUACPowerFlow, +) sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false) pf = ACPowerFlow{ACSolver}() base_res = solve_powerflow(pf, sys) @@ -78,7 +86,7 @@ end @testset "AC Power Flow 3-Bus Fixed FixedAdmittance testing" for ACSolver in ( NLSolveACPowerFlow, - KLUACPowerFlow, + KLUACPowerFlow, LUACPowerFlow, ) p_gen_matpower_3bus = [20.3512373930753, 100.0, 100.0] q_gen_matpower_3bus = [45.516916781567232, 10.453799727283879, -31.992561631394636] @@ -93,7 +101,11 @@ end end @testset "AC Power Flow convergence fail testing" for ACSolver in - (NLSolveACPowerFlow, KLUACPowerFlow) + ( + NLSolveACPowerFlow, + KLUACPowerFlow, + LUACPowerFlow, +) pf_sys5_re = PSB.build_system(PSB.PSITestSystems, "c_sys5_re"; add_forecasts = false) remove_component!(Line, pf_sys5_re, "1") remove_component!(Line, pf_sys5_re, "2") @@ -112,7 +124,11 @@ end end @testset "AC Test 240 Case PSS/e results" for ACSolver in - (NLSolveACPowerFlow, KLUACPowerFlow) + ( + NLSolveACPowerFlow, + KLUACPowerFlow, + LUACPowerFlow, +) file = joinpath( TEST_FILES_DIR, "test_data", @@ -148,7 +164,11 @@ end end @testset "AC Multiple sources at ref" for ACSolver in - (NLSolveACPowerFlow, KLUACPowerFlow) + ( + NLSolveACPowerFlow, + KLUACPowerFlow, + LUACPowerFlow, +) sys = System(100.0) b = ACBus(; number = 1, @@ -193,7 +213,11 @@ end end @testset "AC PowerFlow with Multiple sources at PV" for ACSolver in - (NLSolveACPowerFlow, KLUACPowerFlow) + ( + NLSolveACPowerFlow, + KLUACPowerFlow, + LUACPowerFlow, +) sys = System(100.0) b1 = ACBus(; number = 1, @@ -276,7 +300,11 @@ end end @testset "AC PowerFlow Source + non-source at Ref" for ACSolver in - (NLSolveACPowerFlow, KLUACPowerFlow) + ( + NLSolveACPowerFlow, + KLUACPowerFlow, + LUACPowerFlow, +) sys = System(100.0) b = ACBus(; number = 1, @@ -338,7 +366,11 @@ end end @testset "AC PowerFlow Source + non-source at PV" for ACSolver in - (NLSolveACPowerFlow, KLUACPowerFlow) + ( + NLSolveACPowerFlow, + KLUACPowerFlow, + LUACPowerFlow, +) sys = System(100.0) b1 = ACBus(; number = 1, @@ -434,22 +466,184 @@ end ) end - +# in this test, the following aspects are checked: +# 1. The results of the power flow are consistent for the KLU and NLSolve solvers +# 2. The results of the power flow are consistent for the KLU solver and the legacy implementation +# 3. The Jacobian matrix is the same for the KLU solver and the legacy implementation @testset "Compare larger grid results KLU vs NLSolve" begin sys = build_system(MatpowerTestSystems, "matpower_ACTIVSg2000_sys") pf_default = ACPowerFlow() pf_klu = ACPowerFlow(KLUACPowerFlow) pf_nlsolve = ACPowerFlow(NLSolveACPowerFlow) - + + PSY.set_units_base_system!(sys, "SYSTEM_BASE") + data = PowerFlowData( + pf_default, + sys; + check_connectivity = true) + + ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type) + pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type) + pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type) + pvpq = [pv; pq] + + npvpq = length(pvpq) + npq = length(pq) + + # we need to define lookups for mappings of pv, pq buses onto the internal J indexing + pvpq_lookup = zeros(Int64, maximum([ref; pvpq]) + 1) + pvpq_lookup[pvpq] .= 1:npvpq + pq_lookup = zeros(Int64, maximum([ref; pvpq]) + 1) + pq_lookup[pq] .= 1:npq + + # define the internal J indexing using the lookup arrays + j_pvpq = pvpq_lookup[pvpq] + j_pq = npvpq .+ pq_lookup[pq] + + Vm0 = data.bus_magnitude[:] + Va0 = data.bus_angles[:] + V0 = Vm0 .* exp.(1im * Va0) + + Ybus = data.power_network_matrix.data + rows, cols = SparseArrays.findnz(Ybus) + + diagV = LinearAlgebra.Diagonal(V0) + diagIbus_diag = zeros(Complex{Float64}, size(V0, 1)) + diagIbus = LinearAlgebra.Diagonal(diagIbus_diag) + + diagVnorm = LinearAlgebra.Diagonal(V0 ./ abs.(V0)) + + Ybus_diagVnorm = sparse(rows, cols, Complex{Float64}(0)) + conj_Ybus_diagVnorm = sparse(rows, cols, Complex{Float64}(0)) + diagV_conj_Ybus_diagVnorm = sparse(rows, cols, Complex{Float64}(0)) + conj_diagIbus = conj.(diagIbus) + conj_diagIbus_diagVnorm = conj.(diagIbus) + Ybus_diagV = sparse(rows, cols, Complex{Float64}(0)) + conj_Ybus_diagV = sparse(rows, cols, Complex{Float64}(0)) + + dSbus_dVm = sparse(rows, cols, Complex{Float64}(0)) + dSbus_dVa = sparse(rows, cols, Complex{Float64}(0)) + r_dSbus_dVa = sparse(rows, cols, Float64(0)) + r_dSbus_dVm = sparse(rows, cols, Float64(0)) + i_dSbus_dVa = sparse(rows, cols, Float64(0)) + i_dSbus_dVm = sparse(rows, cols, Float64(0)) + + J_block = sparse(rows, cols, Float64(0), maximum(rows), maximum(cols), unique) + J0_KLU = [J_block[pvpq, pvpq] J_block[pvpq, pq]; J_block[pq, pvpq] J_block[pq, pq]] + PF._update_dSbus_dV!(rows, cols, V0, Ybus, diagV, diagVnorm, diagIbus, diagIbus_diag, + dSbus_dVa, dSbus_dVm, r_dSbus_dVa, r_dSbus_dVm, i_dSbus_dVa, i_dSbus_dVm, + Ybus_diagVnorm, conj_Ybus_diagVnorm, diagV_conj_Ybus_diagVnorm, conj_diagIbus, + conj_diagIbus_diagVnorm, Ybus_diagV, conj_Ybus_diagV) + PF._update_J!( + J0_KLU, + r_dSbus_dVa, + r_dSbus_dVm, + i_dSbus_dVa, + i_dSbus_dVm, + pvpq, + pq, + j_pvpq, + j_pq, + ) + + dSbus_dVa0_LU, dSbus_dVm0_LU = PF._legacy_dSbus_dV(V0, Ybus) + J0_LU = PF._legacy_J(dSbus_dVa, dSbus_dVm, pvpq, pq) + + @test all(isapprox.(J0_LU, J0_KLU, rtol = 0, atol = 1e-12)) + @test all(isapprox.(J0_LU.nzval, J0_KLU.nzval, rtol = 0, atol = 1e-12)) + @test J0_KLU.rowval == J0_LU.rowval + @test J0_KLU.colptr == J0_LU.colptr + res_default = solve_powerflow(pf_default, sys) # must be the same as KLU res_klu = solve_powerflow(pf_klu, sys) res_nlsolve = solve_powerflow(pf_nlsolve, sys) + @test all( + isapprox.( + res_klu["bus_results"][!, :Vm], + res_default["bus_results"][!, :Vm], + rtol = 0, + atol = 1e-12, + ), + ) + @test all( + isapprox.( + res_klu["bus_results"][!, :θ], + res_default["bus_results"][!, :θ], + rtol = 0, + atol = 1e-12, + ), + ) + + @test all( + isapprox.( + res_klu["bus_results"][!, :Vm], + res_nlsolve["bus_results"][!, :Vm], + rtol = 0, + atol = 1e-8, + ), + ) + @test all( + isapprox.( + res_klu["bus_results"][!, :θ], + res_nlsolve["bus_results"][!, :θ], + rtol = 0, + atol = 1e-8, + ), + ) + + # test against legacy implementation + pf_legacy = ACPowerFlow(LUACPowerFlow) + res_legacy = solve_powerflow(pf_legacy, sys) + + @test all( + isapprox.( + res_klu["bus_results"][!, :Vm], + res_legacy["bus_results"][!, :Vm], + rtol = 0, + atol = 1e-12, + ), + ) + @test all( + isapprox.( + res_klu["bus_results"][!, :θ], + res_legacy["bus_results"][!, :θ], + rtol = 0, + atol = 1e-12, + ), + ) - @test all(isapprox.(res_klu["bus_results"][!, :Vm], res_default["bus_results"][!, :Vm], rtol=0, atol=1e-12)) - @test all(isapprox.(res_klu["bus_results"][!, :θ], res_default["bus_results"][!, :θ], rtol=0, atol=1e-12)) + Vm1_KLU = res_klu["bus_results"][!, :Vm] + Va1_KLU = res_klu["bus_results"][!, :θ] + V1_KLU = Vm1_KLU .* exp.(1im * Va1_KLU) + + Vm1_LU = res_legacy["bus_results"][!, :Vm] + Va1_LU = res_legacy["bus_results"][!, :θ] + V1_LU = Vm1_LU .* exp.(1im * Va1_LU) + + J1_KLU = [J_block[pvpq, pvpq] J_block[pvpq, pq]; J_block[pq, pvpq] J_block[pq, pq]] + PF._update_dSbus_dV!(rows, cols, V1_KLU, Ybus, diagV, diagVnorm, diagIbus, + diagIbus_diag, dSbus_dVa, dSbus_dVm, r_dSbus_dVa, r_dSbus_dVm, i_dSbus_dVa, + i_dSbus_dVm, Ybus_diagVnorm, conj_Ybus_diagVnorm, diagV_conj_Ybus_diagVnorm, + conj_diagIbus, conj_diagIbus_diagVnorm, Ybus_diagV, conj_Ybus_diagV) + PF._update_J!( + J1_KLU, + r_dSbus_dVa, + r_dSbus_dVm, + i_dSbus_dVa, + i_dSbus_dVm, + pvpq, + pq, + j_pvpq, + j_pq, + ) - @test all(isapprox.(res_klu["bus_results"][!, :Vm], res_nlsolve["bus_results"][!, :Vm], rtol=0, atol=1e-8)) - @test all(isapprox.(res_klu["bus_results"][!, :θ], res_nlsolve["bus_results"][!, :θ], rtol=0, atol=1e-8)) -end \ No newline at end of file + dSbus_dVa1_LU, dSbus_dVm1_LU = PF._legacy_dSbus_dV(V1_LU, Ybus) + J1_LU = PF._legacy_J(dSbus_dVa1_LU, dSbus_dVm1_LU, pvpq, pq) + + @test all(isapprox.(J1_LU, J1_KLU, rtol = 0, atol = 1e-10)) + @test all(isapprox.(J1_LU.nzval, J1_KLU.nzval, rtol = 0, atol = 1e-10)) + @test J1_KLU.rowval == J1_LU.rowval + @test J1_KLU.colptr == J1_LU.colptr +end diff --git a/test/test_utils/psse_results_compare.jl b/test/test_utils/psse_results_compare.jl index 361cda84..3d35e7c8 100644 --- a/test/test_utils/psse_results_compare.jl +++ b/test/test_utils/psse_results_compare.jl @@ -1,4 +1,4 @@ -function psse_bus_results_compare(file_name, pf_results) +function psse_bus_results_compare(file_name::String, pf_results::Dict) pf_result_bus = CSV.read(file_name, DataFrame) v_diff = Float64[] @@ -15,6 +15,10 @@ function psse_bus_results_compare(file_name, pf_results) return v_diff, angle_diff, number end +function psse_bus_results_compare(file_name::String, pf_results::Bool) + throw(ArgumentError("pf_results not available - calculation failed")) +end + function psse_gen_results_compare(file_name, system::PSY.System) base_power = get_base_power(system) pf_result_gen = CSV.read(file_name, DataFrame)