add tests

mcosovic · May 10, 2024 · 16acb0a · 16acb0a
1 parent 10c2c26
commit 16acb0a
Show file tree

Hide file tree

Showing 10 changed files with 201 additions and 64 deletions.
diff --git a/src/powerSystem/branch.jl b/src/powerSystem/branch.jl
@@ -179,10 +179,7 @@ function addBranch!(system::PowerSystem, analysis::ACPowerFlow{FastNewtonRaphson
     minDiffAngle::A = missing, maxDiffAngle::A = missing,
     longTerm::A = missing, shortTerm::A = missing, emergency::A = missing, type::A = missing)
 
-    bus = system.bus
     branch = system.branch
-    active = analysis.method.active
-    reactive = analysis.method.reactive
 
     addBranch!(system; label, from, to, status, resistance, reactance, susceptance,
     conductance, turnsRatio, shiftAngle, minDiffAngle, maxDiffAngle, longTerm, shortTerm,
@@ -602,7 +599,7 @@ macro branch(kwargs...)
         parameter::Symbol = kwarg.args[1]
 
         if hasfield(BranchTemplate, parameter)
-            if !(parameter in [:status; :type; :label])
+            if !(parameter in [:status; :type; :label; :turnsRatio])
                 container::ContainerTemplate = getfield(template.branch, parameter)
                 if parameter in [:resistance; :reactance]
                     prefixLive = prefix.impedance
@@ -627,6 +624,8 @@ macro branch(kwargs...)
             else
                 if parameter == :status
                     setfield!(template.branch, parameter, Int8(eval(kwarg.args[2])))
+                elseif parameter == :turnsRatio
+                    setfield!(template.branch, parameter, Float64(eval(kwarg.args[2])))
                 elseif parameter == :label
                     label = string(kwarg.args[2])
                     if contains(label, "?")

diff --git a/src/powerSystem/bus.jl b/src/powerSystem/bus.jl
@@ -71,7 +71,7 @@ function addBus!(system::PowerSystem;
 
     push!(bus.layout.type, unitless(type, default.type))
     if !(bus.layout.type[end] in [1, 2, 3])
-        throw(ErrorException("The value $type of the type keyword is illegal."))
+        throw(ErrorException("The value $type of the bus type keyword is illegal."))
     end
     if bus.layout.type[end] == 3
         if bus.layout.slack != 0

diff --git a/src/powerSystem/generator.jl b/src/powerSystem/generator.jl
@@ -450,7 +450,7 @@ function updateGenerator!(system::PowerSystem, analysis::Union{ACPowerFlow{Newto
     if system.bus.layout.type[indexBus] in [2, 3]
         index = system.bus.supply.generator[indexBus][1]
         analysis.voltage.magnitude[indexBus] = generator.voltage.magnitude[index]
-     end
+    end
 end
 
 function updateGenerator!(system::PowerSystem, analysis::ACPowerFlow{GaussSeidel};

diff --git a/src/stateEstimation/acStateEstimation.jl b/src/stateEstimation/acStateEstimation.jl
@@ -127,22 +127,22 @@ function acStateEstimationWLS(system::PowerSystem, device::Measurement)
     nonZeroJacobian = voltmeter.number + 4 * ammeter.number
     nonZeroPrecision = copy(measureNumber)
 
-    for (i, index) in enumerate(wattmeter.layout.index)
+    @inbounds for (i, index) in enumerate(wattmeter.layout.index)
         if wattmeter.layout.bus[i]
             nonZeroJacobian += 2 * (ac.nodalMatrix.colptr[index + 1] - ac.nodalMatrix.colptr[index])
         else
             nonZeroJacobian += 4
         end
     end
-    for (i, index) in enumerate(varmeter.layout.index)
+    @inbounds for (i, index) in enumerate(varmeter.layout.index)
         if varmeter.layout.bus[i]
             nonZeroJacobian += 2 * (ac.nodalMatrix.colptr[index + 1] - ac.nodalMatrix.colptr[index])
         else
             nonZeroJacobian += 4
         end
     end
 
-    for i = 1:pmu.number
+    @inbounds for i = 1:pmu.number
         if pmu.layout.bus[i]
             if pmu.layout.polar[i]
                 nonZeroJacobian += 2
@@ -166,7 +166,7 @@ function acStateEstimationWLS(system::PowerSystem, device::Measurement)
     index = fill(0, measureNumber)
     range = fill(1, 6)
 
-    for (i, k) in enumerate(voltmeter.layout.index)
+    @inbounds for (i, k) in enumerate(voltmeter.layout.index)
         mean[i] = voltmeter.magnitude.status[i] * voltmeter.magnitude.mean[i]
         prec = precisionDiagonal(prec, voltmeter.magnitude.variance[i])
 
@@ -175,7 +175,7 @@ function acStateEstimationWLS(system::PowerSystem, device::Measurement)
     end
     range[2] = jac.idx
 
-    for (i, k) in enumerate(ammeter.layout.index)
+    @inbounds for (i, k) in enumerate(ammeter.layout.index)
         mean[jac.idx] = ammeter.magnitude.status[i] * ammeter.magnitude.mean[i]
         prec = precisionDiagonal(prec, ammeter.magnitude.variance[i])
 
@@ -184,7 +184,7 @@ function acStateEstimationWLS(system::PowerSystem, device::Measurement)
     end
     range[3] = jac.idx
 
-    for (i, k) in enumerate(wattmeter.layout.index)
+    @inbounds for (i, k) in enumerate(wattmeter.layout.index)
         mean[jac.idx] = wattmeter.active.status[i] * wattmeter.active.mean[i]
         prec = precisionDiagonal(prec, wattmeter.active.variance[i])
 
@@ -198,7 +198,7 @@ function acStateEstimationWLS(system::PowerSystem, device::Measurement)
     end
     range[4] = jac.idx
 
-    for (i, k) in enumerate(varmeter.layout.index)
+    @inbounds for (i, k) in enumerate(varmeter.layout.index)
         mean[jac.idx] = varmeter.reactive.status[i] * varmeter.reactive.mean[i]
         prec = precisionDiagonal(prec, varmeter.reactive.variance[i])
 
@@ -212,7 +212,7 @@ function acStateEstimationWLS(system::PowerSystem, device::Measurement)
     end
     range[5] = jac.idx
 
-    for (i, k) in enumerate(pmu.layout.index)
+    @inbounds for (i, k) in enumerate(pmu.layout.index)
         if pmu.layout.polar[i]
             mean[jac.idx] = pmu.magnitude.status[i] * pmu.magnitude.mean[i]
             mean[jac.idx + 1] = pmu.angle.status[i] * pmu.angle.mean[i]
@@ -359,7 +359,7 @@ function acLavStateEstimation(system::PowerSystem, device::Measurement, (@nospec
 
     objective = @expression(se.jump, AffExpr())
 
-    for (k, index) in enumerate(voltmeter.layout.index)
+    @inbounds for (k, index) in enumerate(voltmeter.layout.index)
         index += bus.number
         if voltmeter.magnitude.status[k] == 1
             add_to_expression!(objective, se.residualx[k] + se.residualy[k])
@@ -370,7 +370,7 @@ function acLavStateEstimation(system::PowerSystem, device::Measurement, (@nospec
     end
 
     idx = voltmeter.number + 1
-    for (k, index) in enumerate(ammeter.layout.index)
+    @inbounds for (k, index) in enumerate(ammeter.layout.index)
         if ammeter.magnitude.status[k] == 1
             add_to_expression!(objective, residualx[idx] + residualy[idx])
             addAmmeterResidual!(system, ammeter, se, index, idx, k)
@@ -380,7 +380,7 @@ function acLavStateEstimation(system::PowerSystem, device::Measurement, (@nospec
         idx += 1
     end
 
-    for (k, index) in enumerate(wattmeter.layout.index)
+    @inbounds for (k, index) in enumerate(wattmeter.layout.index)
         if wattmeter.active.status[k] == 1
             add_to_expression!(objective, residualx[idx] + residualy[idx])
             addWattmeterResidual!(system, wattmeter, se, index, idx, k)
@@ -390,7 +390,7 @@ function acLavStateEstimation(system::PowerSystem, device::Measurement, (@nospec
         idx += 1
     end
 
-    for (k, index) in enumerate(varmeter.layout.index)
+    @inbounds for (k, index) in enumerate(varmeter.layout.index)
         if varmeter.reactive.status[k] == 1
             add_to_expression!(objective, residualx[idx] + residualy[idx])
             addVarmeterResidual!(system, varmeter, se, index, idx, k)
@@ -400,7 +400,7 @@ function acLavStateEstimation(system::PowerSystem, device::Measurement, (@nospec
         idx += 1
     end
 
-    for (k, index) in enumerate(pmu.layout.index)
+    @inbounds for (k, index) in enumerate(pmu.layout.index)
         if pmu.layout.polar[k]
             if pmu.layout.bus[k]
                 if pmu.magnitude.status[k] == 1
@@ -476,7 +476,7 @@ function addWattmeterResidual!(system::PowerSystem, wattmeter::Wattmeter, se::LA
         Vi = se.statex[system.bus.number + index] - se.statey[system.bus.number + index]
         expr = @expression(se.jump, Vi * real(system.model.ac.nodalMatrixTranspose[index, index]))
 
-        for ptr in system.model.ac.nodalMatrix.colptr[index]:(system.model.ac.nodalMatrix.colptr[index + 1] - 1)
+        @inbounds for ptr in system.model.ac.nodalMatrix.colptr[index]:(system.model.ac.nodalMatrix.colptr[index + 1] - 1)
             j = system.model.ac.nodalMatrix.rowval[ptr]
             if index != j
                 Gij = real(system.model.ac.nodalMatrixTranspose.nzval[ptr])
@@ -508,7 +508,7 @@ function addVarmeterResidual!(system::PowerSystem, varmeter::Varmeter, se::LAV,
         Vi = se.statex[system.bus.number + index] - se.statey[system.bus.number + index]
         expr = @expression(se.jump, -Vi * imag(system.model.ac.nodalMatrixTranspose[index, index]))
 
-        for ptr in system.model.ac.nodalMatrix.colptr[index]:(system.model.ac.nodalMatrix.colptr[index + 1] - 1)
+        @inbounds for ptr in system.model.ac.nodalMatrix.colptr[index]:(system.model.ac.nodalMatrix.colptr[index + 1] - 1)
             j = system.model.ac.nodalMatrix.rowval[ptr]
             if index != j
                 Gij = real(system.model.ac.nodalMatrixTranspose.nzval[ptr])
@@ -770,7 +770,7 @@ function normalEquation!(system::PowerSystem, analysis::ACStateEstimation)
     voltage = analysis.voltage
     jacobian = se.jacobian
 
-    for col = 1:bus.number
+    @inbounds for col = 1:bus.number
         for lin in jacobian.colptr[col]:(jacobian.colptr[col + 1] - 1)
             row = jacobian.rowval[lin]
             index = se.index[row]
@@ -1015,13 +1015,13 @@ function normalEquation!(system::PowerSystem, analysis::ACStateEstimation)
         end
     end
 
-    for row = se.range[1]:(se.range[2] - 1)
+    @inbounds for row = se.range[1]:(se.range[2] - 1)
         if se.type[row] == 1
             se.residual[row] = se.mean[row] - voltage.magnitude[se.index[row]]
         end
     end
 
-    for row = se.range[5]:(se.range[6] - 1)
+    @inbounds for row = se.range[5]:(se.range[6] - 1)
         if se.type[row] == 10
             se.residual[row] = se.mean[row] - voltage.magnitude[se.index[row]]
         elseif se.type[row] == 11
@@ -1061,7 +1061,7 @@ function jacobianInitialize(jac::SparseModel, val::Int8, col::Int64)
 end
 
 function jacobianInitialize(jac::SparseModel, nodalMatrix::SparseMatrixCSC{ComplexF64,Int64}, bus::Int64, busNumber::Int64)
-    for j in nodalMatrix.colptr[bus]:(nodalMatrix.colptr[bus + 1] - 1)
+    @inbounds for j in nodalMatrix.colptr[bus]:(nodalMatrix.colptr[bus + 1] - 1)
         jac.row[jac.cnt] = jac.idx
         jac.col[jac.cnt] = nodalMatrix.rowval[j]
         jac.row[jac.cnt + 1] = jac.idx

diff --git a/src/stateEstimation/badData.jl b/src/stateEstimation/badData.jl
@@ -68,7 +68,7 @@ function residualTest!(system::PowerSystem, device::Measurement, analysis::DCSta
     JGi = se.coefficient * gainInverse
     idx = findall(!iszero, se.coefficient)
     c = fill(0.0, size(se.coefficient, 1))
-    for i in idx
+    @inbounds for i in idx
         c[i[1]] += JGi[i] * se.coefficient[i]
     end
 
@@ -91,7 +91,7 @@ function residualTest!(system::PowerSystem, device::Measurement, analysis::DCSta
         bad.detect = true
 
         colIndecies = findall(!iszero, se.coefficient[bad.index, :])
-        for col in colIndecies
+        @inbounds for col in colIndecies
             se.coefficient[bad.index, col] = 0.0
         end
         se.mean[bad.index] = 0.0
@@ -133,7 +133,7 @@ function residualTest!(system::PowerSystem, device::Measurement, analysis::PMUSt
     JGi = se.coefficient * gainInverse
     idx = findall(!iszero, se.coefficient)
     c = fill(0.0, size(se.coefficient, 1))
-    for i in idx
+    @inbounds for i in idx
         c[i[1]] += JGi[i] * se.coefficient[i]
     end
 
@@ -161,7 +161,7 @@ function residualTest!(system::PowerSystem, device::Measurement, analysis::PMUSt
         end
 
         colIndecies = findall(!iszero, se.coefficient[bad.index, :])
-        for col in colIndecies
+        @inbounds for col in colIndecies
             se.coefficient[bad.index, col] = 0.0
         end
         se.mean[bad.index] = 0.0
@@ -217,7 +217,7 @@ function residualTest!(system::PowerSystem, device::Measurement, analysis::ACSta
     JGi = se.jacobian * gainInverse
     idx = findall(!iszero, se.jacobian)
     c = fill(0.0, size(se.jacobian, 1))
-    for i in idx
+    @inbounds for i in idx
         c[i[1]] += JGi[i] * se.jacobian[i]
     end
 
@@ -335,7 +335,7 @@ function etree(A)
     n = size(A, 2)
     parent = fill(0, n)
     ancestor = fill(0, n)
-    for k in 1:n, p in A.colptr[k]:(A.colptr[k + 1] - 1)
+    @inbounds for k in 1:n, p in A.colptr[k]:(A.colptr[k + 1] - 1)
         i = A.rowval[p]
         while i != 0 && i < k
             inext = ancestor[i]
@@ -349,15 +349,15 @@ function etree(A)
 
     head = fill(0, n)
     next = fill(0, n)
-    for j in n:-1:1
+    @inbounds for j in n:-1:1
         if parent[j] == 0
             continue
         end
         next[j] = head[parent[j]]
         head[parent[j]] = j
     end
     stack = Int64[]
-    for j in 1:n
+    @inbounds for j in 1:n
         if parent[j] != 0
             continue
         end
@@ -389,7 +389,7 @@ function symbfact(A, parent)
     row = Int64[]; sizehint!(row, n)
 
     visited = falses(n)
-    for k = 1:m
+    @inbounds for k = 1:m
         visited = falses(n)
         visited[k] = true
         for p in Ap[k]:(Ap[k + 1] - 1)
@@ -413,7 +413,7 @@ end
 ### The sparse inverse subset of a real sparse square matrix
 # Copyright 2011, Timothy A. Davis
 # http://www.suitesparse.com
-@inbounds function sparseinv(L, U, d, p, q, Rs, R)
+function sparseinv(L, U, d, p, q, Rs, R)
     Zpattern = R + R'
     n = size(Zpattern, 1)
 
@@ -428,7 +428,7 @@ end
     Zcolptr = Zpattern.colptr
     Zrowval = Zpattern.rowval
     flag = true
-    for j = 1:n
+    @inbounds for j = 1:n
         pdiag = -1
         for p = Zcolptr[j]:(Zcolptr[j + 1] - 1)
             if pdiag == -1 && Zrowval[p] == j
@@ -443,7 +443,7 @@ end
         Zdiagp[j] = pdiag
     end
 
-    if flag
+    @inbounds if flag
         for k = 1:n
             Lmunch[k] = L.colptr[k + 1] - 1
         end
@@ -486,7 +486,7 @@ end
     end
 
     idx = findall(!iszero, Zpattern)
-    for (k, i) in enumerate(idx)
+    @inbounds for (k, i) in enumerate(idx)
         Zpattern[i] = Zx[k]
     end
 

diff --git a/src/stateEstimation/dcStateEstimation.jl b/src/stateEstimation/dcStateEstimation.jl
@@ -104,15 +104,15 @@ function dcStateEstimationWLS(system::PowerSystem, device::Measurement)
     end
 
     nonZeroElement = 0
-    for (i, index) in enumerate(wattmeter.layout.index)
+    @inbounds for (i, index) in enumerate(wattmeter.layout.index)
         if wattmeter.layout.bus[i]
             nonZeroElement += (dc.nodalMatrix.colptr[index + 1] - dc.nodalMatrix.colptr[index])
         else
             nonZeroElement += 2
         end
     end
 
-    for i = 1:pmu.number
+    @inbounds for i = 1:pmu.number
         if pmu.layout.bus[i]
             nonZeroElement += 1
         end
@@ -126,7 +126,7 @@ function dcStateEstimationWLS(system::PowerSystem, device::Measurement)
     precision = spdiagm(0 => mean)
 
     count = 1
-    for (i, k) in enumerate(wattmeter.layout.index)
+    @inbounds for (i, k) in enumerate(wattmeter.layout.index)
         precision.nzval[i] = (1 / wattmeter.active.variance[i])
 
         status = wattmeter.active.status[i]
@@ -163,7 +163,7 @@ function dcStateEstimationWLS(system::PowerSystem, device::Measurement)
 
     rowindex = wattmeter.number + 1
     slackAngle = bus.voltage.angle[bus.layout.slack]
-    for i = 1:pmu.number
+    @inbounds for i = 1:pmu.number
         if pmu.layout.bus[i]
             status = pmu.angle.status[i]
 
@@ -256,7 +256,7 @@ function dcLavStateEstimation(system::PowerSystem, device::Measurement, (@nospec
 
     objective = @expression(jump, AffExpr())
     residual = Dict{Int64, JuMP.ConstraintRef}()
-    for (i, k) in enumerate(wattmeter.layout.index)
+    @inbounds for (i, k) in enumerate(wattmeter.layout.index)
         if device.wattmeter.active.status[i] == 1
             if wattmeter.layout.bus[i]
                 angleCoeff = @expression(jump, AffExpr())
@@ -286,7 +286,7 @@ function dcLavStateEstimation(system::PowerSystem, device::Measurement, (@nospec
     end
 
     slackAngle = bus.voltage.angle[bus.layout.slack]
-    for (i, k) in enumerate(wattmeter.number + 1:deviceNumber)
+    @inbounds for (i, k) in enumerate(wattmeter.number + 1:deviceNumber)
         if pmu.layout.bus[i]
             if pmu.angle.status[i] == 1
                 busIndex = pmu.layout.index[i]
@@ -427,7 +427,7 @@ function solve!(system::PowerSystem, analysis::DCStateEstimation{LAV})
 
     JuMP.optimize!(se.jump)
 
-    for i = 1:system.bus.number
+    @inbounds for i = 1:system.bus.number
         analysis.voltage.angle[i] = value(se.statex[i]::JuMP.VariableRef) - value(se.statey[i]::JuMP.VariableRef) + slackAngle
     end
 end