Skip to content

Commit

Permalink
new version
Browse files Browse the repository at this point in the history
  • Loading branch information
@mcosovic
mcosovic committed Sep 5, 2024
1 parent 03a7151 commit c20ff6a
Show file tree
Hide file tree
Showing 4 changed files with 158 additions and 108 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "JuliaGrid"
uuid = "040ba96c-7d3d-4d4a-931c-37903f53a565"
authors = ["Mirsad Cosovic <[email protected]>"]
version = "0.1.6"
version = "0.1.7"

[deps]
HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
Expand Down
146 changes: 80 additions & 66 deletions src/stateEstimation/badData.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ outlier = residualTest!(system, device, analysis; threshold = 4.0)
solve!(system, analysis)
```
"""
function residualTest!(system::PowerSystem, device::Measurement, analysis::DCStateEstimation{LinearWLS{T}}; threshold = 3.0) where T <: Union{Normal, Orthogonal}
function residualTest!(system::PowerSystem, device::Measurement, analysis::DCStateEstimation{LinearWLS{T}}; threshold::Float64 = 3.0) where T <: Union{Normal, Orthogonal}
errorVoltage(analysis.voltage.angle)

bad = BadData(false, 0.0, "", 0)
Expand Down Expand Up @@ -77,10 +77,13 @@ function residualTest!(system::PowerSystem, device::Measurement, analysis::DCSta
bad.maxNormalizedResidual = 0.0
bad.index = 0
@inbounds for i = 1:se.number
normResidual = abs(se.mean[i] - h[i]) / sqrt(abs((1 / se.precision.nzval[i]) - c[i]))
if normResidual > bad.maxNormalizedResidual
bad.maxNormalizedResidual = normResidual
bad.index = i
residual = se.mean[i] - h[i]
if residual != 0.0
normResidual = abs(residual) / sqrt(abs((1 / se.precision.nzval[i]) - c[i]))
if normResidual > bad.maxNormalizedResidual
bad.maxNormalizedResidual = normResidual
bad.index = i
end
end
end

Expand All @@ -97,22 +100,24 @@ function residualTest!(system::PowerSystem, device::Measurement, analysis::DCSta
se.mean[bad.index] = 0.0
end

if bad.index <= device.wattmeter.number
(bad.label, index),_ = iterate(device.wattmeter.label, bad.index)
if bad.detect
device.wattmeter.active.status[index] = 0
end
else
(bad.label, index),_ = iterate(device.pmu.label, bad.index - device.wattmeter.number)
if bad.detect
device.pmu.angle.status[index] = 0
if bad.index != 0
if bad.index <= device.wattmeter.number
(bad.label, index),_ = iterate(device.wattmeter.label, bad.index)
if bad.detect
device.wattmeter.active.status[index] = 0
end
else
(bad.label, index),_ = iterate(device.pmu.label, bad.index - device.wattmeter.number)
if bad.detect
device.pmu.angle.status[index] = 0
end
end
end

return bad
end

function residualTest!(system::PowerSystem, device::Measurement, analysis::PMUStateEstimation{LinearWLS{T}}; threshold::Union{Float64, Int64} = 3.0) where T <: Union{Normal, Orthogonal}
function residualTest!(system::PowerSystem, device::Measurement, analysis::PMUStateEstimation{LinearWLS{T}}; threshold::Float64 = 3.0) where T <: Union{Normal, Orthogonal}
errorVoltage(analysis.voltage.angle)

bad = BadData(false, 0.0, "", 0)
Expand Down Expand Up @@ -144,10 +149,13 @@ function residualTest!(system::PowerSystem, device::Measurement, analysis::PMUSt
bad.maxNormalizedResidual = 0.0
bad.index = 0
@inbounds for i = 1:se.number
normResidual = abs(se.mean[i] - h[i]) / sqrt(abs((1 / se.precision[i, i]) - c[i]))
if normResidual > bad.maxNormalizedResidual
bad.maxNormalizedResidual = normResidual
bad.index = i
residual = se.mean[i] - h[i]
if residual != 0.0
normResidual = abs(residual) / sqrt(abs((1 / se.precision[i, i]) - c[i]))
if normResidual > bad.maxNormalizedResidual
bad.maxNormalizedResidual = normResidual
bad.index = i
end
end
end

Expand Down Expand Up @@ -179,7 +187,9 @@ function residualTest!(system::PowerSystem, device::Measurement, analysis::PMUSt
pmuIndex = trunc(Int, (bad.index + 1) / 2)
end

(bad.label, ),_ = iterate(device.pmu.label, pmuIndex)
if bad.index != 0
(bad.label, ),_ = iterate(device.pmu.label, pmuIndex)
end
if bad.detect
device.pmu.magnitude.status[pmuIndex] = 0
device.pmu.angle.status[pmuIndex] = 0
Expand All @@ -188,7 +198,7 @@ function residualTest!(system::PowerSystem, device::Measurement, analysis::PMUSt
return bad
end

function residualTest!(system::PowerSystem, device::Measurement, analysis::ACStateEstimation{NonlinearWLS{T}}; threshold::Union{Float64, Int64} = 3.0) where T <: Union{Normal, Orthogonal}
function residualTest!(system::PowerSystem, device::Measurement, analysis::ACStateEstimation{NonlinearWLS{T}}; threshold::Float64 = 3.0) where T <: Union{Normal, Orthogonal}
errorVoltage(analysis.voltage.angle)

bad = BadData(false, 0.0, "", 0)
Expand Down Expand Up @@ -224,62 +234,66 @@ function residualTest!(system::PowerSystem, device::Measurement, analysis::ACSta
bad.maxNormalizedResidual = 0.0
bad.index = 0
@inbounds for i = 1:length(se.mean)
normResidual = abs(se.residual[i]) / sqrt(abs((1 / se.precision[i, i]) - c[i]))
if normResidual > bad.maxNormalizedResidual
bad.maxNormalizedResidual = normResidual
bad.index = i
if se.residual[i] != 0.0
normResidual = abs(se.residual[i]) / sqrt(abs((1 / se.precision[i, i]) - c[i]))
if normResidual > bad.maxNormalizedResidual
bad.maxNormalizedResidual = normResidual
bad.index = i
end
end
end

if bad.maxNormalizedResidual > threshold
bad.detect = true
end

if se.range[1] <= bad.index < se.range[2]
(bad.label, index),_ = iterate(device.voltmeter.label, bad.index)
if bad.detect
device.voltmeter.magnitude.status[index] = 0
end
elseif se.range[2] <= bad.index < se.range[3]
(bad.label, index),_ = iterate(device.ammeter.label, bad.index - device.voltmeter.number)
if bad.detect
device.ammeter.magnitude.status[index] = 0
end
elseif se.range[3] <= bad.index < se.range[4]
(bad.label, index),_ = iterate(device.wattmeter.label, bad.index - device.voltmeter.number - device.ammeter.number)
if bad.detect
device.wattmeter.active.status[index] = 0
end
elseif se.range[4] <= bad.index < se.range[5]
(bad.label, index),_ = iterate(device.varmeter.label, bad.index - device.voltmeter.number - device.ammeter.number - device.wattmeter.number)
if bad.detect
device.varmeter.reactive.status[index] = 0
end
elseif se.range[5] <= bad.index < se.range[6]
badIndex = bad.index - device.voltmeter.number - device.ammeter.number - device.wattmeter.number - device.varmeter.number
if badIndex % 2 == 0
pmuIndex = trunc(Int, badIndex / 2)
alsoBad = bad.index - 1
else
pmuIndex = trunc(Int, (badIndex + 1) / 2)
alsoBad = bad.index + 1
end
if bad.index != 0
if se.range[1] <= bad.index < se.range[2]
(bad.label, index),_ = iterate(device.voltmeter.label, bad.index)
if bad.detect
device.voltmeter.magnitude.status[index] = 0
end
elseif se.range[2] <= bad.index < se.range[3]
(bad.label, index),_ = iterate(device.ammeter.label, bad.index - device.voltmeter.number)
if bad.detect
device.ammeter.magnitude.status[index] = 0
end
elseif se.range[3] <= bad.index < se.range[4]
(bad.label, index),_ = iterate(device.wattmeter.label, bad.index - device.voltmeter.number - device.ammeter.number)
if bad.detect
device.wattmeter.active.status[index] = 0
end
elseif se.range[4] <= bad.index < se.range[5]
(bad.label, index),_ = iterate(device.varmeter.label, bad.index - device.voltmeter.number - device.ammeter.number - device.wattmeter.number)
if bad.detect
device.varmeter.reactive.status[index] = 0
end
elseif se.range[5] <= bad.index < se.range[6]
badIndex = bad.index - device.voltmeter.number - device.ammeter.number - device.wattmeter.number - device.varmeter.number
if badIndex % 2 == 0
pmuIndex = trunc(Int, badIndex / 2)
alsoBad = bad.index - 1
else
pmuIndex = trunc(Int, (badIndex + 1) / 2)
alsoBad = bad.index + 1
end

(bad.label, index),_ = iterate(device.pmu.label, pmuIndex)
if bad.detect
if device.pmu.layout.polar[index]
if se.type[bad.index] in [2; 3; 10]
device.pmu.magnitude.status[index] = 0
(bad.label, index),_ = iterate(device.pmu.label, pmuIndex)
if bad.detect
if device.pmu.layout.polar[index]
if se.type[bad.index] in [2; 3; 10]
device.pmu.magnitude.status[index] = 0
else
device.pmu.angle.status[index] = 0
end
else
device.pmu.magnitude.status[index] = 0
device.pmu.angle.status[index] = 0
end
else
device.pmu.magnitude.status[index] = 0
device.pmu.angle.status[index] = 0

se.mean[alsoBad] = 0.0
se.residual[alsoBad] = 0.0
se.type[alsoBad] = 0
se.mean[alsoBad] = 0.0
se.residual[alsoBad] = 0.0
se.type[alsoBad] = 0
end
end
end
end
Expand Down
114 changes: 75 additions & 39 deletions src/utility/precompile.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@ import PrecompileTools
PrecompileTools.@setup_workload begin
system = powerSystem()
device = measurement()
pseudo = measurement()

addBus!(system; label = 1, type = 3, active = 0.1)
addBus!(system; label = 2, type = 1, reactive = 0.05)
Expand All @@ -11,74 +12,109 @@ PrecompileTools.@setup_workload begin

addPmu!(system, device; bus = 1, magnitude = 1.0, angle = 0.0)
addPmu!(system, device; bus = 2, magnitude = 1.0, angle = 0.0)
addWattmeter!(system, pseudo; bus = 1, active = -0.1)

accompile = Dict(
:injectionPower => injectionPower,
:supplyPower => supplyPower,
:shuntPower => shuntPower,
:fromPower => fromPower,
:toPower => toPower,
:chargingPower => chargingPower,
:seriesPower => seriesPower,
:generatorPower => generatorPower,
:injectionCurrent => injectionCurrent,
:fromCurrent => fromCurrent,
:toCurrent => toCurrent,
:seriesCurrent => seriesCurrent,
)

dccompile = Dict(
:injectionPower => injectionPower,
:supplyPower => supplyPower,
:fromPower => fromPower,
:toPower => toPower,
:generatorPower => generatorPower,
)

PrecompileTools.@compile_workload begin
########## DC Power Flow ###########
analysis = dcPowerFlow(system)
solve!(system, analysis)
power!(system, analysis)

analysis = dcPowerFlow(system, QR)
solve!(system, analysis)

analysis = dcPowerFlow(system, LDLt)
solve!(system, analysis)
########## HDF5 Files ###########
powerSystem("case14.h5")
measurement("measurement14.h5")

########## AC Power Flow ###########
analysis = newtonRaphson(system)
startingVoltage!(system, analysis)
solve!(system, analysis)
power!(system, analysis)
current!(system, analysis)

analysis = newtonRaphson(system, QR)
solve!(system, analysis)
for (name, func) in accompile
func(system, analysis; label = 1)
end

analysis = fastNewtonRaphsonBX(system)
solve!(system, analysis)
power!(system, analysis)
current!(system, analysis)
for (name, func) in accompile
func(system, analysis; label = 1)
end

analysis = gaussSeidel(system)
solve!(system, analysis)

########## DC State Estimation ###########
analysis = dcWlsStateEstimation(system, device)
solve!(system, analysis)
power!(system, analysis)
current!(system, analysis)
for (name, func) in accompile
func(system, analysis; label = 1)
end

analysis = dcWlsStateEstimation(system, device, QR)
solve!(system, analysis)

analysis = dcWlsStateEstimation(system, device, LDLt)
solve!(system, analysis)

analysis = dcWlsStateEstimation(system, device, Orthogonal)
solve!(system, analysis)

########### PMU State Estimation ###########
analysis = pmuWlsStateEstimation(system, device)
########## DC Power Flow ###########
analysis = dcPowerFlow(system)
solve!(system, analysis)
power!(system, analysis)
for (name, func) in dccompile
func(system, analysis; label = 1)
end

analysis = pmuWlsStateEstimation(system, device, QR)
analysis = dcPowerFlow(system, QR)
solve!(system, analysis)

analysis = pmuWlsStateEstimation(system, device, LDLt)
analysis = dcPowerFlow(system, LDLt)
solve!(system, analysis)

analysis = pmuWlsStateEstimation(system, device, Orthogonal)
solve!(system, analysis)
########### Observability Analysis ###########
islands = islandTopologicalFlow(system, device)
restorationGram!(system, device, pseudo, islands)
delete!(accompile, :generatorPower)
delete!(dccompile, :generatorPower)

########### AC State Estimation ###########
########## AC State Estimation ###########
analysis = gaussNewton(system, device)
solve!(system, analysis)
solve!(system, analysis)
residualTest!(system, device, analysis)
power!(system, analysis)
current!(system, analysis)
for (name, func) in accompile
func(system, analysis; label = 1)
end

analysis = gaussNewton(system, device, QR)
solve!(system, analysis)

analysis = gaussNewton(system, device, LDLt)
########## PMU State Estimation ###########
analysis = pmuWlsStateEstimation(system, device)
solve!(system, analysis)
residualTest!(system, device, analysis)
power!(system, analysis)
current!(system, analysis)
for (name, func) in accompile
func(system, analysis; label = 1)
end

analysis = gaussNewton(system, device, Orthogonal)
########### DC State Estimation ###########
analysis = dcWlsStateEstimation(system, device)
solve!(system, analysis)
residualTest!(system, device, analysis)
power!(system, analysis)
for (name, func) in dccompile
func(system, analysis; label = 1)
end
end
end
4 changes: 2 additions & 2 deletions src/utility/routine.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,9 +42,9 @@ end
end

######### Error Voltage ##########
@inline function errorVoltage(voltage::Array{Float64,1})
function errorVoltage(voltage::Array{Float64,1})
if isempty(voltage)
error("The voltage values are missing.")
throw(ErrorException(("The voltage values are missing.")))
end
end

Expand Down

0 comments on commit c20ff6a

Please sign in to comment.