tidy basic jacobian feature, remove decoupled

Author: ccoffrin

CHANGELOG.md

@@ -2,7 +2,8 @@ PowerModels.jl Change Log
 =========================
 
 ### Staged
-- Improved support in `build_pf` for slack buses with multiple generators
+- Added Jacobian calculation of basic network data (#762)
+- Improved support in `build_pf` for slack buses with multiple generators (#785)
 - Fixed bug in `compute_ac_pf` when reporting slack bus `va` values
 - Fixed bug in `constraint_ohms_yt_from_on_off` for DCP models on branches with negative reactance values
 

docs/src/basic-data-utilities.md

@@ -48,6 +48,7 @@ calc_basic_incidence_matrix
 calc_basic_admittance_matrix
 calc_basic_susceptance_matrix
 calc_basic_branch_susceptance_matrix
+calc_basic_jacobian_matrix
 calc_basic_ptdf_matrix
 calc_basic_ptdf_row
 ```
@@ -104,6 +105,7 @@ basic network data using Julia's native linear equation solver,
 compute_basic_dc_pf
 ```
 
+
 !!! tip
     By default PowerModels uses Julia's SparseArrays to ensure the best
     performance of matrix operations on large power network datasets.

src/core/data_basic.jl

@@ -385,76 +385,6 @@ function calc_basic_ptdf_row(data::Dict{String,<:Any}, branch_index::Int)
     return ptdf_column
 end
 
-"""
-Computes the jacobian submatrices
-H = dP/dδ
-L = dQ/dV
-
-The derivatives are ordered according to the bus number, 
-H[1, 1] is ∂P_1 / ∂δ_1
-H[1, 2] is ∂P_1 / ∂δ_2
-
-It does not exclude PV buses rows and columns
-"""
-function calc_basic_decoupled_jacobian_matrices(data::Dict{String,<:Any})
-    if !get(data, "basic_network", false)
-        Memento.warn(_LOGGER, "calc_basic_decoupled_jacobian_matrices requires basic network data and given data may be incompatible. make_basic_network can be used to transform data into the appropriate form.")
-    end
-    num_bus = length(data["bus"])
-    v = calc_basic_bus_voltage(data)
-    vm, va = abs.(v), angle.(v)
-    Y = calc_basic_admittance_matrix(data)
-    neighbors = [Set{Int}([i]) for i in 1:num_bus]
-    I, J, V = findnz(Y)
-    for nz in eachindex(V)
-        push!(neighbors[I[nz]], J[nz])
-        push!(neighbors[J[nz]], I[nz])
-    end
-    H0_I = Int64[]; H0_J = Int64[]; H0_V = Float64[];
-    L0_I = Int64[]; L0_J = Int64[]; L0_V = Float64[];
-    for i in 1:num_bus
-        for j in neighbors[i]
-            push!(H0_I, i); push!(H0_J, j); push!(H0_V, 0.0)
-            push!(L0_I, i); push!(L0_J, j); push!(L0_V, 0.0)
-        end
-    end
-    H = sparse(H0_I, H0_J, H0_V)
-    L = sparse(L0_I, L0_J, L0_V)
-    for i in 1:num_bus
-        for j in neighbors[i]
-            bus_type = data["bus"]["$(j)"]["bus_type"]
-            if bus_type == 1
-                if i == j
-                    y_ii = Y[i,i]
-                    H[i, j] =                      + vm[i] * sum( -real(Y[i,k]) * vm[k] * sin(va[i] - va[k]) + imag(Y[i,k]) * vm[k] * cos(va[i] - va[k]) for k in neighbors[i] if k != i )
-                    L[i, j] = - 2*imag(y_ii)*vm[i] +         sum(  real(Y[i,k]) * vm[k] * sin(va[i] - va[k]) - imag(Y[i,k]) * vm[k] * cos(va[i] - va[k]) for k in neighbors[i] if k != i )
-                else
-                    y_ij = Y[i,j]
-                    H[i, j] = vm[i] * vm[j] * (  real(y_ij) * sin(va[i] - va[j]) - imag(y_ij) * cos(va[i] - va[j]) )
-                    L[i, j] =         vm[i] * (  real(y_ij) * sin(va[i] - va[j]) - imag(y_ij) * cos(va[i] - va[j]) )
-                end
-            elseif bus_type == 2
-                if i == j
-                    y_ii = Y[i,i]
-                    H[i, j] = vm[i] * sum( -real(Y[i,k]) * vm[k] * sin(va[i] - va[k]) + imag(Y[i,k]) * vm[k] * cos(va[i] - va[k]) for k in neighbors[i] if k != i )
-                    L[i, j] = 1.0
-                else
-                    y_ij = Y[i,j]
-                    H[i, j] = vm[i] * vm[j] * (  real(y_ij) * sin(va[i] - va[j]) - imag(y_ij) * cos(va[i] - va[j]) )
-                    L[i, j] = 0.0
-                end
-            elseif bus_type == 3
-                if i == j
-                    H[i, j] = 1.0
-                    L[i, j] = 1.0
-                end
-            else
-                @assert false
-            end
-        end
-    end
-    return H, L
-end
 
 """
 given a basic network data dict, returns a sparse real valued Jacobian matrix
@@ -465,6 +395,7 @@ function calc_basic_jacobian_matrix(data::Dict{String,<:Any})
     if !get(data, "basic_network", false)
         Memento.warn(_LOGGER, "calc_basic_jacobian_matrix requires basic network data and given data may be incompatible. make_basic_network can be used to transform data into the appropriate form.")
     end
+
     num_bus = length(data["bus"])
     v = calc_basic_bus_voltage(data)
     vm, va = abs.(v), angle.(v)

test/common.jl

@@ -72,7 +72,7 @@ end
 """
 An AC Power Flow Solver from scratch. 
 """
-function compute_basic_ac_pf!(data::Dict{String, Any}; decoupled=false)
+function compute_basic_ac_pf!(data::Dict{String, Any})
     if !get(data, "basic_network", false)
         Memento.warn(_LOGGER, "compute_basic_ac_pf requires basic network data and given data may be incompatible. make_basic_network can be used to transform data into the appropriate form.")
     end
@@ -104,16 +104,11 @@ function compute_basic_ac_pf!(data::Dict{String, Any}; decoupled=false)
         if LinearAlgebra.normInf([delta_P; delta_Q]) < tol
             break
         end
+
         # STEP 2 and 3: Compute the jacobian and update step
-        if !decoupled
-            J = calc_basic_jacobian_matrix(data)
-            x = J \ [delta_P; delta_Q]
-        else
-            H, L = calc_basic_decoupled_jacobian_matrices(data)
-            va = H \ delta_P
-            vm = L \ delta_Q
-            x = [va; vm]
-        end
+        J = calc_basic_jacobian_matrix(data)
+        x = J \ [delta_P; delta_Q]
+
         # STEP 4
         # update voltage variables
         for i in 1:bus_num

test/data-basic.jl

@@ -195,7 +195,7 @@ end
             for (i, bus) in data["bus"] # All buses PQ
                 bus["bus_type"] = 1
             end
-            J = PowerModels.calc_basic_jacobian_matrix(data);
+            J = calc_basic_jacobian_matrix(data)
 
             num_bus = length(data["bus"])
 
@@ -212,7 +212,7 @@ end
 
         @testset "24-bus-case" begin
             data = make_basic_network(PowerModels.parse_file("../test/data/matpower/case24.m"))
-            J = PowerModels.calc_basic_jacobian_matrix(data)
+            J = calc_basic_jacobian_matrix(data)
 
             num_bus = length(data["bus"])
 
@@ -244,15 +244,6 @@ end
             @test isapprox(J[3, num_bus+9], -2.1624; atol=1e-4)          # dP/dvm non-diagonal
             @test isapprox(J[num_bus+3, num_bus+9], -7.9603; atol=1e-4)  # dQ/dvm non-diagonal
         end
-
-        @testset "30-bus-case-decoupled-matrices" begin
-            data = make_basic_network(PowerModels.parse_file("../test/data/matpower/case24.m"))
-            H, L = calc_basic_decoupled_jacobian_matrices(data)
-            J = calc_basic_jacobian_matrix(data)
-            n = length(data["bus"])
-            @test isapprox(J[1:n, 1:n], H; atol=1e-6)
-            @test isapprox(J[n+1:end, n+1:end], L; atol=1e-6)
-        end
     end
 
     @testset "basic ac power flow" begin
@@ -273,23 +264,6 @@ end
             end
         end
 
-        @testset "9-bus-case" begin     
-            data = make_basic_network(PowerModels.parse_file("../test/data/matpower/case9.m"))
-            solution = compute_ac_pf(data)["solution"]
-
-            compute_basic_ac_pf!(data; decoupled=true)
-            
-            for (i, bus) in data["bus"]
-                @test isapprox(solution["bus"][i]["va"], bus["va"]; atol=1e-4)
-                @test isapprox(solution["bus"][i]["vm"], bus["vm"]; atol=1e-4)
-            end
-            
-            for (i, gen) in data["gen"]
-                @test isapprox(solution["gen"][i]["pg"], gen["pg"]; atol=1e-4)
-                @test isapprox(solution["gen"][i]["qg"], gen["qg"]; atol=1e-4)
-            end
-        end
-
         @testset "14-bus-case" begin     
             data = make_basic_network(PowerModels.parse_file("../test/data/matpower/case14.m"))
             solution = compute_ac_pf(data)["solution"]