Skip to content

Commit

Permalink
reusing DC State Estimation
Browse files Browse the repository at this point in the history
  • Loading branch information
@mcosovic
mcosovic committed Dec 21, 2023
1 parent 5fa7545 commit fca26b9
Show file tree
Hide file tree
Showing 14 changed files with 951 additions and 160 deletions.
20 changes: 15 additions & 5 deletions src/definition/analysis.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,11 +213,22 @@ mutable struct BadData
label::String
end

mutable struct WattmeterLayout
bus::Dict{Int64,Array{Int64,1}}
from::Dict{Int64,Array{Int64,1}}
to::Dict{Int64,Array{Int64,1}}
number::Int64
end

mutable struct PMULayout
bus::Dict{Int64,Array{Int64,1}}
number::Int64
end

mutable struct DCStateEstimationLayout
index::Array{Int64,1}
device::Int64
wattmeter::Int64
pmu::Int64
wattmeter::WattmeterLayout
pmu::PMULayout
number::Int64
end

mutable struct DCStateEstimationWLSMethod
Expand Down Expand Up @@ -250,7 +261,6 @@ mutable struct DCStateEstimationMethodLAV
jump::JuMP.Model
variable::VariableLAV
residual::Dict{Int64, JuMP.ConstraintRef}
objective::JuMP.AffExpr
layout::DCStateEstimationLayout
end

Expand Down
50 changes: 23 additions & 27 deletions src/measurement/ammeter.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,20 +77,20 @@ function addAmmeter!(system::PowerSystem, device::Measurement;
ammeter.number += 1
labelBranch = getLabel(system.branch, location, "branch")

index = system.branch.label[getLabel(system.branch, location, "branch")]
push!(ammeter.layout.index, index)
indexBranch = system.branch.label[getLabel(system.branch, location, "branch")]
push!(ammeter.layout.index, indexBranch)

basePowerInv = 1 / (system.base.power.value * system.base.power.prefix)
if ammeter.layout.from[end]
setLabel(ammeter, label, default.label, labelBranch; prefix = "From ")
defaultVariance = default.varianceFrom
defaultStatus = default.statusFrom
baseVoltage = system.base.voltage.value[system.branch.layout.from[index]] * system.base.voltage.prefix
baseVoltage = system.base.voltage.value[system.branch.layout.from[indexBranch]] * system.base.voltage.prefix
else
setLabel(ammeter, label, default.label, labelBranch; prefix = "To ")
defaultVariance = default.varianceTo
defaultStatus = default.statusTo
baseVoltage = system.base.voltage.value[system.branch.layout.to[index]] * system.base.voltage.prefix
baseVoltage = system.base.voltage.value[system.branch.layout.to[indexBranch]] * system.base.voltage.prefix
end

baseCurrentInv = baseCurrentInverse(basePowerInv, baseVoltage)
Expand Down Expand Up @@ -201,38 +201,34 @@ function addAmmeter!(system::PowerSystem, device::Measurement, analysis::AC;
ammeter.magnitude.variance = similar(ammeter.magnitude.mean)
ammeter.magnitude.status = fill(Int8(0), ammeterNumber)

count = 1
basePowerInv = 1 / (system.base.power.value * system.base.power.prefix)
label = collect(keys(sort(system.branch.label; byvalue = true)))
@inbounds for i = 1:2:ammeterNumber
labelBranch = getLabel(system.branch, label[count], "branch")

ammeter.number += 1
setLabel(ammeter, missing, default.label, labelBranch; prefix = "From ")

@inbounds for (label, i) in system.branch.label
ammeter.number += 1
setLabel(ammeter, missing, default.label, labelBranch; prefix = "To ")
setLabel(ammeter, missing, default.label, label; prefix = "From ")

ammeter.layout.index[i] = count
ammeter.layout.index[i + 1] = count
ammeter.layout.from[i] = true
ammeter.layout.to[i + 1] = true
ammeter.layout.index[ammeter.number] = i
ammeter.layout.from[ammeter.number] = true

baseVoltage = system.base.voltage.value[system.branch.layout.from[count]] * system.base.voltage.prefix
baseVoltage = system.base.voltage.value[system.branch.layout.from[i]] * system.base.voltage.prefix
baseCurrentInv = baseCurrentInverse(basePowerInv, baseVoltage)
ammeter.magnitude.variance[i] = topu(varianceFrom, default.varianceFrom, prefix.currentMagnitude, baseCurrentInv)
ammeter.magnitude.mean[i] = analysis.current.from.magnitude[count] + ammeter.magnitude.variance[i]^(1/2) * randn(1)[1]
ammeter.magnitude.status[i] = statusFrom

baseVoltage = system.base.voltage.value[system.branch.layout.to[count]] * system.base.voltage.prefix
ammeter.magnitude.variance[ammeter.number] = topu(varianceFrom, default.varianceFrom, prefix.currentMagnitude, baseCurrentInv)
ammeter.magnitude.mean[ammeter.number] = analysis.current.from.magnitude[i] + ammeter.magnitude.variance[ammeter.number]^(1/2) * randn(1)[1]
ammeter.magnitude.status[ammeter.number] = statusFrom

ammeter.number += 1
setLabel(ammeter, missing, default.label, label; prefix = "To ")

ammeter.layout.index[ammeter.number] = i
ammeter.layout.to[ammeter.number] = true

baseVoltage = system.base.voltage.value[system.branch.layout.to[i]] * system.base.voltage.prefix
baseCurrentInv = baseCurrentInverse(basePowerInv, baseVoltage)
ammeter.magnitude.variance[i + 1] = topu(varianceTo, default.varianceTo, prefix.currentMagnitude, baseCurrentInv)
ammeter.magnitude.mean[i + 1] = analysis.current.to.magnitude[count] + ammeter.magnitude.variance[i + 1]^(1/2) * randn(1)[1]
ammeter.magnitude.status[i + 1] = statusTo

count += 1
ammeter.magnitude.variance[ammeter.number] = topu(varianceTo, default.varianceTo, prefix.currentMagnitude, baseCurrentInv)
ammeter.magnitude.mean[ammeter.number] = analysis.current.to.magnitude[i] + ammeter.magnitude.variance[ammeter.number]^(1/2) * randn(1)[1]
ammeter.magnitude.status[ammeter.number] = statusTo
end

ammeter.layout.label = ammeter.number
end

Expand Down
139 changes: 106 additions & 33 deletions src/measurement/pmu.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -276,11 +276,9 @@ function addPmu!(system::PowerSystem, device::Measurement, analysis::AC;
pmu.angle.status = similar(pmu.magnitude.status)

prefixInv = 1 / system.base.voltage.prefix
label = collect(keys(sort(system.bus.label; byvalue = true)))
@inbounds for i = 1:system.bus.number
labelBus = getLabel(system.bus, label[i], "bus")
@inbounds for (label, i) in system.bus.label
pmu.number += 1
setLabel(pmu, missing, default.label, labelBus)
setLabel(pmu, missing, default.label, label)

pmu.layout.index[i] = i
pmu.layout.bus[i] = true
Expand All @@ -294,44 +292,41 @@ function addPmu!(system::PowerSystem, device::Measurement, analysis::AC;
pmu.angle.status[i] = statusAngleBus
end

count = 1
basePowerInv = 1 / (system.base.power.value * system.base.power.prefix)
label = collect(keys(sort(system.branch.label; byvalue = true)))
@inbounds for i = (system.bus.number + 1):2:pmuNumber
labelBranch = getLabel(system.branch, label[count], "branch")

@inbounds for (label, i) in system.branch.label
pmu.number += 1
setLabel(pmu, missing, default.label, labelBranch; prefix = "From ")
setLabel(pmu, missing, default.label, label; prefix = "From ")

pmu.number += 1
setLabel(pmu, missing, default.label, labelBranch; prefix = "To ")

pmu.layout.index[i] = count
pmu.layout.index[i + 1] = count
pmu.layout.from[i] = true
pmu.layout.to[i + 1] = true
pmu.magnitude.status[i] = statusMagnitudeFrom
pmu.magnitude.status[i + 1] = statusMagnitudeTo
pmu.angle.status[i] = statusAngleFrom
pmu.angle.status[i + 1] = statusAngleTo

baseVoltage = system.base.voltage.value[system.branch.layout.from[count]] * system.base.voltage.prefix
pmu.layout.index[pmu.number] = i
pmu.layout.from[pmu.number] = true

baseVoltage = system.base.voltage.value[system.branch.layout.from[i]] * system.base.voltage.prefix
baseCurrentInv = baseCurrentInverse(basePowerInv, baseVoltage)
pmu.magnitude.variance[i] = topu(varianceMagnitudeFrom, default.varianceMagnitudeFrom, prefix.currentMagnitude, baseCurrentInv)
pmu.magnitude.mean[i] = analysis.current.from.magnitude[count] + pmu.magnitude.variance[i]^(1/2) * randn(1)[1]

pmu.angle.variance[i] = topu(varianceAngleFrom, default.varianceAngleFrom, prefix.currentAngle, 1.0)
pmu.angle.mean[i] = analysis.current.from.angle[count] + pmu.angle.variance[i]^(1/2) * randn(1)[1]
pmu.magnitude.variance[pmu.number] = topu(varianceMagnitudeFrom, default.varianceMagnitudeFrom, prefix.currentMagnitude, baseCurrentInv)
pmu.magnitude.mean[pmu.number] = analysis.current.from.magnitude[i] + pmu.magnitude.variance[pmu.number]^(1/2) * randn(1)[1]
pmu.magnitude.status[pmu.number] = statusMagnitudeFrom

baseVoltage = system.base.voltage.value[system.branch.layout.to[count]] * system.base.voltage.prefix
pmu.angle.variance[pmu.number] = topu(varianceAngleFrom, default.varianceAngleFrom, prefix.currentAngle, 1.0)
pmu.angle.mean[pmu.number] = analysis.current.from.angle[i] + pmu.angle.variance[pmu.number]^(1/2) * randn(1)[1]
pmu.angle.status[pmu.number] = statusAngleFrom

pmu.number += 1
setLabel(pmu, missing, default.label, label; prefix = "To ")

pmu.layout.index[pmu.number] = i
pmu.layout.to[pmu.number] = true

baseVoltage = system.base.voltage.value[system.branch.layout.to[i]] * system.base.voltage.prefix
baseCurrentInv = baseCurrentInverse(basePowerInv, baseVoltage)
pmu.magnitude.variance[i + 1] = topu(varianceMagnitudeTo, default.varianceMagnitudeTo, prefix.currentMagnitude, baseCurrentInv)
pmu.magnitude.mean[i + 1] = analysis.current.to.magnitude[count] + pmu.magnitude.variance[i + 1]^(1/2) * randn(1)[1]

pmu.angle.variance[i + 1] = topu(varianceAngleTo, default.varianceAngleTo, prefix.currentAngle, 1.0)
pmu.angle.mean[i + 1] = analysis.current.to.angle[count] + pmu.angle.variance[i + 1]^(1/2) * randn(1)[1]
pmu.magnitude.variance[pmu.number] = topu(varianceMagnitudeTo, default.varianceMagnitudeTo, prefix.currentMagnitude, baseCurrentInv)
pmu.magnitude.mean[pmu.number] = analysis.current.to.magnitude[i] + pmu.magnitude.variance[pmu.number]^(1/2) * randn(1)[1]
pmu.magnitude.status[pmu.number] = statusMagnitudeTo

count += 1
pmu.angle.variance[pmu.number] = topu(varianceAngleTo, default.varianceAngleTo, prefix.currentAngle, 1.0)
pmu.angle.mean[pmu.number] = analysis.current.to.angle[i] + pmu.angle.variance[pmu.number]^(1/2) * randn(1)[1]
pmu.angle.status[pmu.number] = statusAngleTo
end

pmu.layout.label = pmu.number
Expand Down Expand Up @@ -408,6 +403,84 @@ function updatePmu!(system::PowerSystem, device::Measurement; label::L,
prefixAngle, 1.0)
end

function updatePmu!(system::PowerSystem, device::Measurement, analysis::DCStateEstimationWLS;
label::L, magnitude::T = missing, angle::T = missing, varianceMagnitude::T = missing,
varianceAngle::T = missing, statusMagnitude::T = missing, statusAngle::T = missing,
noise::Bool = template.pmu.noise)

pmu = device.pmu
method = analysis.method

updatePmu!(system, device; label, magnitude, angle, varianceMagnitude, varianceAngle,
statusMagnitude, statusAngle, noise)

indexPmu = pmu.label[getLabel(pmu, label, "PMU")]
if (isset(statusAngle) || isset(angle) || isset(varianceAngle)) && pmu.layout.bus[indexPmu]
constIf = constMeter(pmu.angle.status[indexPmu])

indexBus = pmu.layout.index[indexPmu]
index = indexPmu + device.wattmeter.number

if isset(statusAngle)
method.jacobian[index, indexBus] = constIf
end
if isset(statusAngle) || isset(angle)
method.mean[index] = (pmu.angle.mean[indexPmu] - system.bus.voltage.angle[system.bus.layout.slack]) * constIf
end
if isset(statusAngle) || isset(varianceAngle)
method.weight[index] = constIf / pmu.angle.variance[indexPmu]
end
end
end

function updatePmu!(system::PowerSystem, device::Measurement, analysis::DCStateEstimationLAV;
label::L, magnitude::T = missing, angle::T = missing, varianceMagnitude::T = missing,
varianceAngle::T = missing, statusMagnitude::T = missing, statusAngle::T = missing,
noise::Bool = template.pmu.noise)

pmu = device.pmu
method = analysis.method

updatePmu!(system, device; label, magnitude, angle, varianceMagnitude, varianceAngle,
statusMagnitude, statusAngle, noise)

indexPmu = pmu.label[getLabel(pmu, label, "PMU")]
index = indexPmu + device.wattmeter.number
if pmu.layout.bus[indexPmu] && isset(statusAngle)
if pmu.angle.status[indexPmu] == 1
indexBus = pmu.layout.index[indexPmu]

if is_fixed(method.variable.residualx[index])
unfix(method.variable.residualx[index])
set_lower_bound(method.variable.residualx[index], 0.0)

unfix(method.variable.residualy[index])
set_lower_bound(method.variable.residualy[index], 0.0)

set_objective_coefficient(method.jump, method.variable.residualx[index], 1)
set_objective_coefficient(method.jump, method.variable.residualy[index], 1)
end

remove!(method.jump, method.residual, index)
method.residual[index] = @constraint(method.jump, method.variable.angley[indexBus] - method.variable.anglex[indexBus] + method.variable.residualy[index] - method.variable.residualx[index] == 0.0)
else
remove!(method.jump, method.residual, index)

if !is_fixed(method.variable.residualx[index])
fix(method.variable.residualx[index], 0.0; force = true)
fix(method.variable.residualy[index], 0.0; force = true)

set_objective_coefficient(method.jump, method.variable.residualx[index], 0)
set_objective_coefficient(method.jump, method.variable.residualy[index], 0)
end
end
end

if pmu.layout.bus[indexPmu] && pmu.angle.status[indexPmu] == 1 && (isset(statusAngle) || isset(angle))
JuMP.set_normalized_rhs(method.residual[index], pmu.angle.mean[indexPmu] - system.bus.voltage.angle[system.bus.layout.slack])
end
end

"""
@pmu(label, varianceMagnitudeBus, statusMagnitudeBus, varianceAngleBus, statusAngleBus,
varianceMagnitudeFrom, statusMagnitudeFrom, varianceAngleFrom, statusAngleFrom,
Expand Down
Loading

0 comments on commit fca26b9

Please sign in to comment.