Ref ravens schema new tests

src/core/constraint_template.jl

@@ -718,6 +718,7 @@ function constraint_mc_transformer_power(pm::AbstractUnbalancedPowerModel, i::In
     tm_scale = calculate_tm_scale(transformer, ref(pm, nw, :bus, f_bus), ref(pm, nw, :bus, t_bus))
 
     #TODO change data model
+    #look at less redundant way to specify polarity
     # there is redundancy in specifying polarity seperately on from and to side
     #TODO change this once migrated to new data model
     pol = transformer["polarity"]

src/core/constraint_template_en.jl

@@ -118,8 +118,10 @@ function constraint_mc_transformer_voltage(pm::ExplicitNeutralModels, i::Int; nw
     tm_scale = calculate_tm_scale(transformer, ref(pm, nw, :bus, f_bus), ref(pm, nw, :bus, t_bus))
 
     #TODO change data model
+    #note SR 07/01: change which data model? make what change? more info needed
     # there is redundancy in specifying polarity seperately on from and to side
     #TODO change this once migrated to new data model
+    #SR note 7/7/26: the above comment was 4 years ago. ready to be migrated?
     pol = transformer["polarity"]
 
     if configuration == WYE

src/core/data.jl

@@ -213,8 +213,12 @@ end
 
 
 "wraps angles in radians to pi"
-function _wrap_to_pi(radians)
-    return radians - 2*pi*floor.((radians .+ pi)/(2*pi))
+function _wrap_to_pi(x)
+    return mod(x + pi, 2*pi) - pi
+end
+
+function _wrap_to_pi(x::AbstractVector{<:Real})
+    return mod.(x .+ pi, 2*pi) .- pi
 end
 
 
@@ -434,22 +438,28 @@ end
 
 
 "Calculates the tap scale factor for the non-dimensionalized equations."
-function calculate_tm_scale(trans::Dict{String,Any}, bus_fr::Dict{String,Any}, bus_to::Dict{String,Any})
+function calculate_tm_scale(
+    trans::Dict{String,Any},
+    bus_fr::Dict{String,Any},
+    bus_to::Dict{String,Any},
+)
     tm_nom = trans["tm_nom"]
 
-    f_vbase = haskey(bus_fr, "vbase") ? bus_fr["vbase"] : bus_fr["base_kv"]
-    t_vbase = haskey(bus_to, "vbase") ? bus_to["vbase"] : bus_to["base_kv"]
+    f_vbase = haskey(bus_fr, "vbase") ? bus_fr["vbase"] : error("from bus is missing vbase")
+    t_vbase = haskey(bus_to, "vbase") ? bus_to["vbase"] : error("to bus is missing vbase")
     config = trans["configuration"]
 
-    tm_scale = tm_nom*(t_vbase/f_vbase)
+    tm_scale = tm_nom * (t_vbase / f_vbase)
+
     if config == DELTA
-        #TODO is this still needed?
+        # TODO: confirm whether delta-side voltage bases are line-to-line or line-to-neutral
+        #this is for kron reduce, line to line
         tm_scale *= sqrt(3)
     elseif config == "zig-zag"
         error("Zig-zag not yet supported.")
     end
 
-    return tm_nom
+    return tm_scale #not returning tm_nom
 end
 
 
@@ -800,91 +810,81 @@ function _mat2ltrivec!(m::Union{Matrix{T}, LinearAlgebra.Symmetric{T}}) where T
 end
 
 
-"makes a hermitian matrix variable from diagonal, and lower and upper triangular vectors"
-function _make_hermitian_matrix_variable(diag, lowertrianglereal, lowertriangleimag)
-    #TODO clean up
-    matrixreal = []
-    if length(diag) == 3
-        matrixreal = [
-        diag[1]                 lowertrianglereal[1]    lowertrianglereal[2];
-        lowertrianglereal[1]    diag[2]                 lowertrianglereal[3];
-        lowertrianglereal[2]    lowertrianglereal[3]    diag[3]
-        ]
-    elseif length(diag) == 4
-        matrixreal = [
-        diag[1]                 lowertrianglereal[1]    lowertrianglereal[2]    lowertrianglereal[4];
-        lowertrianglereal[1]    diag[2]                 lowertrianglereal[3]    lowertrianglereal[5];
-        lowertrianglereal[2]    lowertrianglereal[3]    diag[3]                 lowertrianglereal[6];
-        lowertrianglereal[4]    lowertrianglereal[5]    lowertrianglereal[6]    diag[4]
-        ]
-    elseif length(diag) == 5
-        matrixreal = [
-        diag[1]                 lowertrianglereal[1]    lowertrianglereal[2]    lowertrianglereal[4]    lowertrianglereal[7];
-        lowertrianglereal[1]    diag[2]                 lowertrianglereal[3]    lowertrianglereal[5]    lowertrianglereal[8];
-        lowertrianglereal[2]    lowertrianglereal[3]    diag[3]                 lowertrianglereal[6]    lowertrianglereal[9];
-        lowertrianglereal[4]    lowertrianglereal[5]    lowertrianglereal[6]    diag[4]                 lowertrianglereal[10];
-        lowertrianglereal[7]    lowertrianglereal[8]    lowertrianglereal[9]    lowertrianglereal[10]    diag[5]
-        ]
-    end
-
-    matriximag = []
-    if length(diag) == 3
-        matriximag = [
-        0                       -lowertriangleimag[1]   -lowertriangleimag[2];
-        lowertriangleimag[1]    0                       -lowertriangleimag[3];
-        lowertriangleimag[2]    lowertriangleimag[3]    0
-        ]
-    elseif length(diag) == 4
-        matriximag = [
-        0                       -lowertriangleimag[1]   -lowertriangleimag[2]   -lowertriangleimag[4];
-        lowertriangleimag[1]    0                       -lowertriangleimag[3]   -lowertriangleimag[5];
-        lowertriangleimag[2]    lowertriangleimag[3]    0                       -lowertriangleimag[6];
-        lowertriangleimag[4]    lowertriangleimag[5]    lowertriangleimag[6]    0
-        ]
-    elseif length(diag) == 5
-        matriximag = [
-        0                       -lowertriangleimag[1]   -lowertriangleimag[2]   -lowertriangleimag[4]   -lowertriangleimag[7];
-        lowertriangleimag[1]    0                       -lowertriangleimag[3]   -lowertriangleimag[5]   -lowertriangleimag[8];
-        lowertriangleimag[2]    lowertriangleimag[3]    0                       -lowertriangleimag[6]   -lowertriangleimag[9];
-        lowertriangleimag[4]    lowertriangleimag[5]    lowertriangleimag[6]    0                       -lowertriangleimag[10];
-        lowertriangleimag[7]    lowertriangleimag[8]    lowertriangleimag[9]    lowertriangleimag[10]    0
-        ]
+"""
+Makes a Hermitian matrix variable from a diagonal and lower-triangular real/imaginary vectors.
+
+Vector ordering is assumed to be:
+(2,1), (3,1), (3,2), (4,1), (4,2), (4,3), ...
+"""
+function _make_hermitian_matrix_variable(
+    diag::AbstractVector{Tdiag},
+    lowertrianglereal::AbstractVector{Treal},
+    lowertriangleimag::AbstractVector{Timag},
+) where {Tdiag, Treal, Timag}
+
+    n = length(diag)
+    expected = n * (n - 1) ÷ 2
+
+    @assert length(lowertrianglereal) == expected
+    @assert length(lowertriangleimag) == expected
+
+    Trealout = promote_type(Tdiag, Treal)
+
+    matrixreal = LinearAlgebra.diagm(Trealout.(diag))
+    matriximag = fill(zero(Timag), n, n)
+
+    k = 1
+    for i in 2:n
+        for j in 1:i-1
+            matrixreal[i, j] = lowertrianglereal[k]
+            matrixreal[j, i] = lowertrianglereal[k]
+
+            matriximag[i, j] = lowertriangleimag[k]
+            matriximag[j, i] = -lowertriangleimag[k]
+
+            k += 1
+        end
     end
+
     return matrixreal, matriximag
 end
 
+"""
+Makes a full matrix variable from a diagonal, lower-triangular vector,
+and upper-triangular vector.
+
+Vector ordering is assumed to be:
+(2,1), (3,1), (3,2), (4,1), (4,2), (4,3), ...
+"""
+function _make_full_matrix_variable(
+    diag::AbstractVector{Tdiag},
+    lowertriangle::AbstractVector{Tlower},
+    uppertriangle::AbstractVector{Tupper},
+) where {Tdiag, Tlower, Tupper}
+
+    n = length(diag)
+    expected = n * (n - 1) ÷ 2
+
+    @assert length(lowertriangle) == expected
+    @assert length(uppertriangle) == expected
+
+    Tout = promote_type(Tdiag, Tlower, Tupper) #general type that can hold these elements
+    matrix = Matrix{Tout}(undef, n, n)
+
+    matrix = LinearAlgebra.diagm(Tout.(diag)) #cast "diag" to this type, use as diagonal
+
+    k = 1
+    for i in 2:n
+        for j in 1:i-1
+            matrix[i, j] = lowertriangle[k]
+            matrix[j, i] = uppertriangle[k]
+            k += 1
+        end
+    end
 
-"makes a full matrix variable from a diagonal, and lower and upper triangular vectors"
-function _make_full_matrix_variable(diag::Vector{T}, lowertriangle::Vector{T}, uppertriangle::Vector{T}) where T
-    #TODO clean up
-    matrix = []
-    if length(diag) == 3
-        matrix = [
-        diag[1]             uppertriangle[1]    uppertriangle[2];
-        lowertriangle[1]    diag[2]             uppertriangle[3];
-        lowertriangle[2]    lowertriangle[3]    diag[3]
-        ]
-    elseif length(diag) == 4
-        matrix = [
-        diag[1]             uppertriangle[1]    uppertriangle[2]    uppertriangle[4];
-        lowertriangle[1]    diag[2]             uppertriangle[3]    uppertriangle[5];
-        lowertriangle[2]    lowertriangle[3]    diag[3]             uppertriangle[6];
-        lowertriangle[4]    lowertriangle[5]    lowertriangle[6]    diag[4]
-        ]
-    elseif length(diag) == 5
-        matrix = [
-        diag[1]             uppertriangle[1]    uppertriangle[2]    uppertriangle[4]    uppertriangle[7];
-        lowertriangle[1]    diag[2]             uppertriangle[3]    uppertriangle[5]    uppertriangle[8];
-        lowertriangle[2]    lowertriangle[3]    diag[3]             uppertriangle[6]    uppertriangle[9];
-        lowertriangle[4]    lowertriangle[5]    lowertriangle[6]    diag[4]             uppertriangle[10]
-        lowertriangle[7]    lowertriangle[8]    lowertriangle[9]    lowertriangle[10]    diag[5]
-        ]
-    end
-    # matrix = LinearAlgebra.diagm(0 => diag) + _vec2ltri!(lowertriangle) + _vec2utri!(uppertriangle)
     return matrix
 end
 
-
 """
     correct_mc_voltage_angle_differences!(data::Dict{String,<:Any}, default_pad::Real=deg2rad(10.0))
 
@@ -983,6 +983,7 @@ end
 
 "helper function to build bus shunt matrices for power balance constraints"
 function _build_bus_shunt_matrices(pm::AbstractUnbalancedPowerModel, nw::Int, terminals::Vector{Int}, bus_shunts::Vector{<:Tuple{Int,Vector{Int}}})::Tuple{Matrix{<:Real},Matrix{<:Real}}
+    #secondary place to put off-diag shunt matrix element code
     ncnds = length(terminals)
     Gs = fill(0.0, ncnds, ncnds)
     Bs = fill(0.0, ncnds, ncnds)
@@ -1358,7 +1359,7 @@ end
 function _check_branch_loops(data_math::MathematicalModel{NetworkModel})
     for (i, branch) in data_math["branch"]
         if branch["f_bus"] == branch["t_bus"]
-            error("both sides of branch $(i) connect to bus $(branch["f_bus"])")
+            error("both sides of branch $(i) connect to bus $(branch["f_bus"])\n include `transformations = [transform_loops!]` in the `parse_file` call")
         end
     end
 end
@@ -1915,7 +1916,7 @@ end
 
 
 "infer the internal dimension of a winding, load or generator based on the connections and the configuration"
-function _infer_int_dim(connections::Vector, configuration::ConnConfig, kron_reduced)
+function _infer_int_dim(connections::Vector, configuration::ConnConfig, kron_reduced::Bool=false)
     if configuration==WYE
         if kron_reduced
             return length(connections)

src/core/solution.jl

@@ -34,7 +34,34 @@ function _IM.build_solution_values(var::LinearAlgebra.Symmetric{JuMP.VariableRef
     return JuMP.value.(var.data)
 end
 
+function _order_phase_angles_abc(va::AbstractVector{<:Real})::Vector{Float64}
+    @assert length(va) == 3
+
+    # a near 0, b near -120 degrees, c near +120 degrees
+    refs = [0.0, -2*pi/3, 2*pi/3]
+
+    perms = (
+        (1, 2, 3),
+        (1, 3, 2),
+        (2, 1, 3),
+        (2, 3, 1),
+        (3, 1, 2),
+        (3, 2, 1),
+    )
+
+    best_perm = perms[1]
+    best_cost = Inf
+
+    for p in perms
+        cost = sum(abs2(_wrap_to_pi(va[p[k]] - refs[k])) for k in 1:3)
+        if cost < best_cost
+            best_cost = cost
+            best_perm = p
+        end
+    end
 
+    return [_wrap_to_pi(va[best_perm[k]]) for k in 1:3]
+end
 "converts w models voltages to standard voltage magnitude (sqrt)"
 function _sol_data_model_w!(solution::Dict{String,<:Any})
     if haskey(solution, "nw")
@@ -67,9 +94,8 @@ function _sol_data_model_w!(solution::Dict{String,<:Any})
                             va = LinearAlgebra.pinv(t)*[atan(bus["Wi"][2,1], bus["Wr"][2,1]);
                                                         atan(bus["Wi"][3,1], bus["Wr"][3,1]);
                                                         atan(bus["Wi"][3,2], bus["Wr"][3,2])]
-                            bus["va"] = [va[findmin(abs.(va .- 0))[2]],
-                                         va[findmin(abs.(va .+ 2*pi/3))[2]],
-                                         va[findmin(abs.(va .- 2*pi/3))[2]]] # TODO: better way to get angles in order
+
+                            bus["va"] = _order_phase_angles_abc(va)
                         end
                     end
                     delete!(bus, "Wr")

src/data_model/checks.jl

@@ -514,10 +514,25 @@ function _check_transformer(data_eng::Dict{String,<:Any}, name::Any)
         conf = transformer["configuration"][w]
         conns = transformer["connections"][w]
         nph = conf==WYE ? length(conns)-1 : length(conns)
+        context = "transformer $name winding $w"
         @assert all(nph.==nphs) "transformer $name: winding $w has a different number of phases than the previous ones."
 
+        @assert length(transformer["rw"][w]) == nph "$context: rw length must match number of phases"
+        @assert length(transformer["tm_fix"][w]) == nph "$context: tm_fix length must match number of phases"
+        @assert length(transformer["tm_set"][w]) == nph "$context: tm_set length must match number of phases"
+        @assert length(transformer["tm_lb"][w]) == nph "$context: tm_lb length must match number of phases"
+        @assert length(transformer["tm_ub"][w]) == nph "$context: tm_ub length must match number of phases"
+        @assert length(transformer["tm_step"][w]) == nph "$context: tm_step length must match number of phases"
+
+        @assert all(transformer["rw"][w] .>= 0.0) "$context: rw must be nonnegative"
+        @assert all(transformer["tm_lb"][w] .<= transformer["tm_set"][w]) "$context: tm_lb must be <= tm_set"
+        @assert all(transformer["tm_set"][w] .<= transformer["tm_ub"][w]) "$context: tm_set must be <= tm_ub"
+        @assert all(transformer["tm_step"][w] .> 0.0) "$context: tm_step must be positive"
+
+        @assert transformer["vm_nom"][w] > 0.0 "$context: vm_nom must be positive"
+        @assert transformer["sm_nom"][w] > 0.0 "$context: sm_nom must be positive"
+
         push!(nphs, nph)
-        #TODO check length other properties
     end
 
     _check_connectivity(data_eng, transformer; context="transformer_nw $name")

src/data_model/dss/edge_constructors.jl

@@ -81,6 +81,7 @@ function create_dss_object(::Type{T}, property_pairs::Vector{Pair{String,String}
     line.cmatrix = :cmatrix ∈ raw_fields ? line.cmatrix : imag(Yc) / (2 * pi * line.basefreq)
 
     # TODO: support length mismatch between line and linecode?
+    # DSS _does_ allow this, here would likely be good place to throw error
     # Currently this does not change the values of rmatrix and xmatrix due to
     # lenmult=1, code is only in place for future use.
     lenmult = 1.0
@@ -153,6 +154,7 @@ function create_dss_object(::Type{T}, property_pairs::Vector{Pair{String,String}
     end
 
     # TODO: handle `parallel`
+    # parallel is a dss reactor type, we don't support that here
     if (:kv ∈ raw_fields && :kvar ∈ raw_fields) || :x ∈ raw_fields || :lmh ∈ raw_fields || :z ∈ raw_fields
         if :kvar ∈ raw_fields && :kv ∈ raw_fields
             kvarperphase = reactor.kvar / reactor.phases
@@ -195,6 +197,9 @@ function create_dss_object(::Type{T}, property_pairs::Vector{Pair{String,String}
         reactor.x = reactor.xmatrix[1, 1]
 
         # TODO: account for off-diagonal and single phase
+        # likely not a big issue
+        # trying to collapse full matrix into dss-style scalars
+        # assumes all matrices are diagonal, all phases are identical
         reactor.z = reactor.z1 = reactor.z2 = reactor.z0 = [reactor.r, reactor.x]
         reactor.lmh = reactor.x / (2 * pi) / reactor.basefreq * 1e3
     elseif :z1 ∈ raw_fields

src/data_model/transformations/dss2eng.jl

@@ -75,6 +75,7 @@ function _dss2eng_load!(data_eng::EngineeringModel{NetworkModel}, data_dss::Open
         nphases = dss_obj["phases"]
         bus = _parse_bus_id(dss_obj["bus1"])[1]
         conf = nphases==1 && dss_obj["kv"]==0.24 ? DELTA : dss_obj["conn"] # check if load is connected between split-phase terminals of triplex node (nominal line-line voltage=240V), TODO: better generalization
+        #^logic for center tap xfmr?
 
         if conf==DELTA
             @assert(nphases in [1, 3], "$id: only 1 and 3-phase delta loads are supported!")
@@ -544,6 +545,7 @@ function _dss2eng_xfmrcode!(data_eng::EngineeringModel{NetworkModel}, data_dss::
         nphases = dss_obj["phases"]
         nrw = dss_obj["windings"]
 
+        #maybe assert n windings not greater than 3?
         eng_obj = Dict{String,Any}(
             "tm_set" => Vector{Vector{Float64}}([fill(tap, nphases) for tap in dss_obj["taps"]]),
             "tm_lb" => Vector{Vector{Float64}}(fill(fill(dss_obj["mintap"], nphases), nrw)),
@@ -557,7 +559,8 @@ function _dss2eng_xfmrcode!(data_eng::EngineeringModel{NetworkModel}, data_dss::
             "rw" => Vector{Float64}(dss_obj["%rs"] ./ 100),
             "noloadloss" => dss_obj["%noloadloss"] / 100,
             "cmag" => dss_obj["%imag"] / 100,
-            "xsc" => nrw == 2 ? [dss_obj["xhl"] / 100] : [dss_obj["xhl"], dss_obj["xht"], dss_obj["xlt"]] ./ 100,
+            #what if there are more than 3 windings?
+            "xsc" => nrw == 2 ? [dss_obj["xhl"] / 100] : [dss_obj["xhl"], dss_obj["xht"], dss_obj["xlt"]] ./ 100, #units will be SI
             "source_id" => "xfmrcode.$id",
         )
 

src/data_model/transformations/eng2math.jl

@@ -657,16 +657,16 @@ end
 
 "converts engineering generic shunt components into mathematical shunt components"
 function _map_eng2math_shunt!(data_math::MathematicalModel{NetworkModel}, data_eng::EngineeringModel{NetworkModel}; pass_props::Vector{String}=String[])
+    #start here for allowing off diag matrix shunts
     for (name, eng_obj) in get(data_eng, "shunt", Dict{Any,Dict{String,Any}}())
         math_obj = _init_math_obj("shunt", name, eng_obj, length(data_math["shunt"])+1; pass_props=pass_props)
 
-        # TODO change to new capacitor shunt calc logic
         math_obj["shunt_bus"] = data_math["bus_lookup"][eng_obj["bus"]]
 
-        math_obj["gs"] = get(eng_obj, "gs", zeros(size(eng_obj["bs"])))
+        math_obj["gs"] = get(eng_obj, "gs", zeros(size(eng_obj["bs"]))) #fill missing g's with zero array shaped like bs
 
         data_math["shunt"]["$(math_obj["index"])"] = math_obj
-
+        
         # add capcontrol items to math model
         if haskey(eng_obj,"controls")
             math_obj["controls"] = deepcopy(eng_obj["controls"])