From fca26b99bbd6d164198481eefda24c80da2ec3f0 Mon Sep 17 00:00:00 2001 From: mcosovic Date: Thu, 21 Dec 2023 14:08:56 +0100 Subject: [PATCH] reusing DC State Estimation --- src/definition/analysis.jl | 20 +- src/measurement/ammeter.jl | 50 +++-- src/measurement/pmu.jl | 139 ++++++++++---- src/measurement/powermeter.jl | 171 ++++++++++++++--- src/measurement/voltmeter.jl | 6 +- src/powerSystem/branch.jl | 117 ++++++++++++ src/powerSystem/bus.jl | 41 +++- src/powerSystem/generator.jl | 41 ++++ src/stateEstimation/badData.jl | 2 +- src/stateEstimation/dcStateEstimation.jl | 192 ++++++++++++------- test/runtests.jl | 3 +- test/stateEstimation/analysis.jl | 65 +++++++ test/stateEstimation/badData.jl | 33 +++- test/stateEstimation/reusing.jl | 231 +++++++++++++++++++++++ 14 files changed, 951 insertions(+), 160 deletions(-) create mode 100644 test/stateEstimation/reusing.jl diff --git a/src/definition/analysis.jl b/src/definition/analysis.jl index bd1a2dc5c..55e66b7d3 100644 --- a/src/definition/analysis.jl +++ b/src/definition/analysis.jl @@ -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 @@ -250,7 +261,6 @@ mutable struct DCStateEstimationMethodLAV jump::JuMP.Model variable::VariableLAV residual::Dict{Int64, JuMP.ConstraintRef} - objective::JuMP.AffExpr layout::DCStateEstimationLayout end diff --git a/src/measurement/ammeter.jl b/src/measurement/ammeter.jl index 8d54d31df..3aac42ed3 100644 --- a/src/measurement/ammeter.jl +++ b/src/measurement/ammeter.jl @@ -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) @@ -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 diff --git a/src/measurement/pmu.jl b/src/measurement/pmu.jl index e868d5a7d..96c3c546d 100644 --- a/src/measurement/pmu.jl +++ b/src/measurement/pmu.jl @@ -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 @@ -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 @@ -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, diff --git a/src/measurement/powermeter.jl b/src/measurement/powermeter.jl index 197fa4b4b..82b3c63e4 100644 --- a/src/measurement/powermeter.jl +++ b/src/measurement/powermeter.jl @@ -396,11 +396,9 @@ function addPowermeter!(system, device, measure, powerBus, powerFrom, powerTo, d measure.status = fill(Int8(0), deviceNumber) basePowerInv = 1 / (system.base.power.value * system.base.power.prefix) - label = collect(keys(sort(system.bus.label; byvalue = true))) - @inbounds for i = 1:system.bus.number + @inbounds for (label, i) in system.bus.label device.number += 1 - labelBus = getLabel(system.bus, label[i], "bus") - setLabel(device, missing, default.label, labelBus) + setLabel(device, missing, default.label, label) device.layout.index[i] = i device.layout.bus[i] = true @@ -410,31 +408,26 @@ function addPowermeter!(system, device, measure, powerBus, powerFrom, powerTo, d measure.status[i] = statusBus end - count = 1 - label = collect(keys(sort(system.branch.label; byvalue = true))) - @inbounds for i = (system.bus.number + 1):2:deviceNumber - labelBranch = getLabel(system.branch, label[count], "branch") - + @inbounds for (label, i) in system.branch.label device.number += 1 - setLabel(device, missing, default.label, labelBranch; prefix = "From ") + setLabel(device, missing, default.label, label; prefix = "From ") - device.number += 1 - setLabel(device, missing, default.label, labelBranch; prefix = "To ") + device.layout.index[device.number] = i + device.layout.from[device.number] = true - device.layout.index[i] = count - device.layout.index[i + 1] = count - device.layout.from[i] = true - device.layout.to[i + 1] = true + measure.variance[device.number] = topu(varianceFrom, default.varianceFrom, prefixPower, basePowerInv) + measure.mean[device.number] = powerFrom[i] + measure.variance[device.number]^(1/2) * randn(1)[1] + measure.status[device.number] = statusFrom - measure.variance[i] = topu(varianceFrom, default.varianceFrom, prefixPower, basePowerInv) - measure.mean[i] = powerFrom[count] + measure.variance[i]^(1/2) * randn(1)[1] - measure.status[i] = statusFrom + device.number += 1 + setLabel(device, missing, default.label, label; prefix = "To ") - measure.variance[i + 1] = topu(varianceTo, default.varianceTo, prefixPower, basePowerInv) - measure.mean[i + 1] = powerTo[count] + measure.variance[i + 1]^(1/2) * randn(1)[1] - measure.status[i + 1] = statusTo + device.layout.index[device.number] = i + device.layout.to[device.number] = true - count += 1 + measure.variance[device.number] = topu(varianceTo, default.varianceTo, prefixPower, basePowerInv) + measure.mean[device.number] = powerTo[i] + measure.variance[device.number]^(1/2) * randn(1)[1] + measure.status[device.number] = statusTo end device.layout.label = device.number @@ -493,6 +486,138 @@ function updateWattmeter!(system::PowerSystem, device::Measurement; label::L, prefix.activePower, basePowerInv) end +function updateWattmeter!(system::PowerSystem, device::Measurement, analysis::DCStateEstimationWLS; + label::L, active::T = missing, variance::T = missing, status::T = missing, + noise::Bool = template.wattmeter.noise) + + dc = system.model.dc + wattmeter = device.wattmeter + method = analysis.method + + indexWattmeter = wattmeter.label[getLabel(wattmeter, label, "wattmeter")] + basePowerInv = 1 / (system.base.power.value * system.base.power.prefix) + + updateMeter(wattmeter.active, indexWattmeter, active, variance, status, noise, + prefix.activePower, basePowerInv) + + constIf = constMeter(wattmeter.active.status[indexWattmeter]) + if isset(status) || isset(active) || isset(variance) + if wattmeter.layout.bus[indexWattmeter] + indexBus = wattmeter.layout.index[indexWattmeter] + if isset(status) + method.done = false + for j in dc.nodalMatrix.colptr[indexBus]:(dc.nodalMatrix.colptr[indexBus + 1] - 1) + method.jacobian[indexWattmeter, dc.nodalMatrix.rowval[j]] = dc.nodalMatrix.nzval[j] * constIf + end + end + if isset(status) || isset(active) + method.mean[indexWattmeter] = (wattmeter.active.mean[indexWattmeter] - dc.shiftPower[indexBus] - system.bus.shunt.conductance[indexBus]) * constIf + end + if isset(status) || isset(variance) + method.weight[indexWattmeter] = constIf / wattmeter.active.variance[indexWattmeter] + end + else + indexBranch = wattmeter.layout.index[indexWattmeter] + if wattmeter.layout.from[indexWattmeter] + addmitance = dc.admittance[indexBranch] * constIf + else + addmitance = -dc.admittance[indexBranch] * constIf + end + if isset(status) + method.done = false + method.jacobian[indexWattmeter, system.branch.layout.from[indexBranch]] = addmitance + method.jacobian[indexWattmeter, system.branch.layout.to[indexBranch]] = -addmitance + end + if isset(status) || isset(active) + method.mean[indexWattmeter] = wattmeter.active.mean[indexWattmeter] * constIf + system.branch.parameter.shiftAngle[indexBranch] * addmitance + end + if isset(status) || isset(variance) + method.weight[indexWattmeter] = constIf / wattmeter.active.variance[indexWattmeter] + end + end + end +end + +function updateWattmeter!(system::PowerSystem, device::Measurement, analysis::DCStateEstimationLAV; + label::L, active::T = missing, variance::T = missing, status::T = missing, + noise::Bool = template.wattmeter.noise) + + dc = system.model.dc + wattmeter = device.wattmeter + method = analysis.method + + indexWattmeter = wattmeter.label[getLabel(wattmeter, label, "wattmeter")] + basePowerInv = 1 / (system.base.power.value * system.base.power.prefix) + + updateMeter(wattmeter.active, indexWattmeter, active, variance, status, noise, + prefix.activePower, basePowerInv) + + if isset(status) || isset(active) + if wattmeter.layout.bus[indexWattmeter] + indexBus = wattmeter.layout.index[indexWattmeter] + else + indexBranch = wattmeter.layout.index[indexWattmeter] + if wattmeter.layout.from[indexWattmeter] + admittance = dc.admittance[indexBranch] + else + admittance = -dc.admittance[indexBranch] + end + end + end + + if isset(status) + if wattmeter.active.status[indexWattmeter] == 1 + if is_fixed(method.variable.residualx[indexWattmeter]) + unfix(method.variable.residualx[indexWattmeter]) + set_lower_bound(method.variable.residualx[indexWattmeter], 0.0) + + unfix(method.variable.residualy[indexWattmeter]) + set_lower_bound(method.variable.residualy[indexWattmeter], 0.0) + + set_objective_coefficient(method.jump, method.variable.residualx[indexWattmeter], 1) + set_objective_coefficient(method.jump, method.variable.residualy[indexWattmeter], 1) + end + + if wattmeter.layout.bus[indexWattmeter] + angleJacobian = @expression(method.jump, AffExpr()) + for j in dc.nodalMatrix.colptr[indexBus]:(dc.nodalMatrix.colptr[indexBus + 1] - 1) + k = dc.nodalMatrix.rowval[j] + add_to_expression!(angleJacobian, dc.nodalMatrix.nzval[j] * (method.variable.angley[k] - method.variable.anglex[k])) + end + + remove!(method.jump, method.residual, indexWattmeter) + method.residual[indexWattmeter] = @constraint(method.jump, angleJacobian + method.variable.residualy[indexWattmeter] - method.variable.residualx[indexWattmeter] == 0.0) + else + from = system.branch.layout.from[indexBranch] + to = system.branch.layout.to[indexBranch] + + angleJacobian = admittance * (method.variable.angley[from] - method.variable.anglex[from] - method.variable.angley[to] + method.variable.anglex[to]) + + remove!(method.jump, method.residual, indexWattmeter) + method.residual[indexWattmeter] = @constraint(method.jump, angleJacobian + method.variable.residualy[indexWattmeter] - method.variable.residualx[indexWattmeter] == 0.0) + end + else + remove!(method.jump, method.residual, indexWattmeter) + + if !is_fixed(method.variable.residualx[indexWattmeter]) + fix(method.variable.residualx[indexWattmeter], 0.0; force = true) + fix(method.variable.residualy[indexWattmeter], 0.0; force = true) + + set_objective_coefficient(method.jump, method.variable.residualx[indexWattmeter], 0) + set_objective_coefficient(method.jump, method.variable.residualy[indexWattmeter], 0) + end + end + end + + if wattmeter.active.status[indexWattmeter] == 1 && (isset(status) || isset(active)) + if wattmeter.layout.bus[indexWattmeter] + JuMP.set_normalized_rhs(method.residual[indexWattmeter], wattmeter.active.mean[indexWattmeter] - dc.shiftPower[indexBus] - system.bus.shunt.conductance[indexBus]) + else + JuMP.set_normalized_rhs(method.residual[indexWattmeter], wattmeter.active.mean[indexWattmeter] + system.branch.parameter.shiftAngle[indexBranch] * admittance) + end + end +end + """ updateVarmeter!(system::PowerSystem, device::Measurement, analysis::Analysis; kwargs...) diff --git a/src/measurement/voltmeter.jl b/src/measurement/voltmeter.jl index 6ffab7da0..6c920583e 100644 --- a/src/measurement/voltmeter.jl +++ b/src/measurement/voltmeter.jl @@ -65,10 +65,10 @@ function addVoltmeter!(system::PowerSystem, device::Measurement; labelBus = getLabel(system.bus, bus, "bus") setLabel(voltmeter, label, default.label, labelBus) - index = system.bus.label[getLabel(system.bus, bus, "bus")] - push!(voltmeter.layout.index, index) + indexBus = system.bus.label[labelBus] + push!(voltmeter.layout.index, indexBus) - baseVoltageInv = 1 / (system.base.voltage.value[index] * system.base.voltage.prefix) + baseVoltageInv = 1 / (system.base.voltage.value[indexBus] * system.base.voltage.prefix) setMeter(voltmeter.magnitude, magnitude, variance, status, noise, default.variance, default.status, prefix.voltageMagnitude, baseVoltageInv) end diff --git a/src/powerSystem/branch.jl b/src/powerSystem/branch.jl index b3e167271..f83c6536f 100644 --- a/src/powerSystem/branch.jl +++ b/src/powerSystem/branch.jl @@ -249,6 +249,60 @@ function addBranch!(system::PowerSystem, analysis::ACOptimalPowerFlow; end end +function addBranch!(system::PowerSystem, analysis::DCStateEstimationWLS; + label::L = missing, from::L, to::L, status::T = missing, + resistance::T = missing, reactance::T = missing, susceptance::T = missing, + conductance::T = missing, turnsRatio::T = missing, shiftAngle::T = missing, + minDiffAngle::T = missing, maxDiffAngle::T = missing, + longTerm::T = missing, shortTerm::T = missing, emergency::T = missing, type::T = missing) + + dc = system.model.dc + branch = system.branch + method = analysis.method + + if isset(shiftAngle) + oldShiftAngleFrom = dc.shiftPower[system.bus.label[getLabel(system.bus, from, "bus")]] + oldShiftAngleTo = dc.shiftPower[system.bus.label[getLabel(system.bus, to, "bus")]] + end + + addBranch!(system; label, from, to, status, resistance, reactance, susceptance, + conductance, turnsRatio, shiftAngle, minDiffAngle, maxDiffAngle, longTerm, shortTerm, + emergency, type) + + if branch.layout.status[end] == 1 + indexBus = Int64[] + if haskey(method.layout.wattmeter.bus, branch.layout.from[end]) + push!(indexBus, branch.layout.from[end]) + if isset(shiftAngle) + for indexWattmeter in method.layout.wattmeter.bus[branch.layout.from[end]] + method.mean[indexWattmeter] += oldShiftAngleFrom + end + end + end + if haskey(method.layout.wattmeter.bus, branch.layout.to[end]) + push!(indexBus, branch.layout.to[end]) + if isset(shiftAngle) + for indexWattmeter in method.layout.wattmeter.bus[branch.layout.to[end]] + method.mean[indexWattmeter] += oldShiftAngleTo + end + end + end + if !isempty(indexBus) + method.done = false + for index in indexBus + for indexWattmeter in method.layout.wattmeter.bus[index] + for j in dc.nodalMatrix.colptr[index]:(dc.nodalMatrix.colptr[index + 1] - 1) + method.jacobian[indexWattmeter, dc.nodalMatrix.rowval[j]] = dc.nodalMatrix.nzval[j] + end + if isset(shiftAngle) + method.mean[indexWattmeter] -= - dc.shiftPower[index] + end + end + end + end + end +end + """ updateBranch!(system::PowerSystem, analysis::Analysis; kwargs...) @@ -539,6 +593,69 @@ function updateBranch!(system::PowerSystem, analysis::ACOptimalPowerFlow; end end +function updateBranch!(system::PowerSystem, analysis::DCStateEstimationWLS; + label::L, status::T = missing, resistance::T = missing, reactance::T = missing, + susceptance::T = missing, conductance::T = missing, turnsRatio::T = missing, + shiftAngle::T = missing, minDiffAngle::T = missing, maxDiffAngle::T = missing, + longTerm::T = missing, shortTerm::T = missing, emergency::T = missing, type::T = missing) + + dc = system.model.dc + branch = system.branch + method = analysis.method + + indexBranch = branch.label[getLabel(branch, label, "branch")] + indexBuses = [branch.layout.from[indexBranch]; branch.layout.to[indexBranch]] + oldShiftPower = [dc.shiftPower[indexBuses[1]]; dc.shiftPower[indexBuses[2]]] + oldShiftAdmittance = branch.parameter.shiftAngle[indexBranch] * dc.admittance[indexBranch] + + updateBranch!(system; label, status, resistance, reactance, susceptance, + conductance, turnsRatio, shiftAngle, minDiffAngle, maxDiffAngle, longTerm, shortTerm, + emergency, type) + + if isset(reactance) || isset(turnsRatio) || isset(shiftAngle) || isset(status) + if haskey(method.layout.wattmeter.from, indexBranch) + for indexWattmeter in method.layout.wattmeter.from[indexBranch] + constStatus = constMeter(method.weight[indexWattmeter]) + + if isset(reactance) || isset(turnsRatio) || isset(status) + method.done = false + method.jacobian[indexWattmeter, branch.layout.from[indexBranch]] = dc.admittance[indexBranch] * constStatus + method.jacobian[indexWattmeter, branch.layout.to[indexBranch]] = -dc.admittance[indexBranch] * constStatus + end + method.mean[indexWattmeter] = ((method.mean[indexWattmeter] - oldShiftAdmittance) + branch.parameter.shiftAngle[indexBranch] * dc.admittance[indexBranch]) * constStatus + end + end + if haskey(method.layout.wattmeter.to, indexBranch) + for indexWattmeter in method.layout.wattmeter.to[indexBranch] + constStatus = constMeter(method.weight[indexWattmeter]) + + if isset(reactance) || isset(turnsRatio) || isset(status) + method.done = false + method.jacobian[indexWattmeter, branch.layout.from[indexBranch]] = -dc.admittance[indexBranch] * constStatus + method.jacobian[indexWattmeter, branch.layout.to[indexBranch]] = dc.admittance[indexBranch] * constStatus + end + method.mean[indexWattmeter] = ((method.mean[indexWattmeter] + oldShiftAdmittance) - branch.parameter.shiftAngle[indexBranch] * dc.admittance[indexBranch]) * constStatus + end + end + + for (k, indexBus) in enumerate(indexBuses) + if haskey(method.layout.wattmeter.bus, indexBus) + for indexWattmeter in method.layout.wattmeter.bus[indexBus] + constStatus = constMeter(method.weight[indexWattmeter]) + + if isset(reactance) || isset(turnsRatio) || isset(status) + method.done = false + for j in dc.nodalMatrix.colptr[indexBus]:(dc.nodalMatrix.colptr[indexBus + 1] - 1) + method.jacobian[indexWattmeter, dc.nodalMatrix.rowval[j]] = dc.nodalMatrix.nzval[j] * constStatus + end + end + method.mean[indexWattmeter] = (method.mean[indexWattmeter] + oldShiftPower[k] - dc.shiftPower[indexBus]) * constStatus + end + end + end + end +end + """ @branch(kwargs...) diff --git a/src/powerSystem/bus.jl b/src/powerSystem/bus.jl index 4e520af78..84ed8badd 100644 --- a/src/powerSystem/bus.jl +++ b/src/powerSystem/bus.jl @@ -129,11 +129,14 @@ function addBus!(system::PowerSystem, analysis::ACPowerFlow; kwargs...) throw(ErrorException("The AC power flow model cannot be reused when adding a new bus.")) end -######### Query About Bus ########## function addBus!(system::PowerSystem, analysis::DCOptimalPowerFlow; kwargs...) throw(ErrorException("The DC optimal power flow model cannot be reused when adding a new bus.")) end +function addBus!(system::PowerSystem, analysis::DCStateEstimation; kwargs...) + throw(ErrorException("The DC state estimation model cannot be reused when adding a new bus.")) +end + """ updateBus!(system::PowerSystem, analysis::Analysis; kwargs...) @@ -443,6 +446,42 @@ function updateBus!(system::PowerSystem, analysis::ACOptimalPowerFlow; end end +function updateBus!(system::PowerSystem, analysis::DCStateEstimationWLS; + label::L, type::T = missing, + active::T = missing, reactive::T = missing, + conductance::T = missing, susceptance::T = missing, + magnitude::T = missing, angle::T = missing, + minMagnitude::T = missing, maxMagnitude::T = missing, + base::T = missing, area::T = missing, lossZone::T = missing) + + bus = system.bus + method = analysis.method + + indexBus = bus.label[getLabel(bus, label, "bus")] + + if isset(type) && bus.layout.slack == indexBus && type != 3 + throw(ErrorException("The DC state estimation model cannot be reused due to required bus type conversion.")) + end + + oldConductance = bus.shunt.conductance[indexBus] + oldSlackAngle = bus.voltage.angle[bus.layout.slack] + + updateBus!(system; label, type, active, reactive, conductance, susceptance, + magnitude, angle, minMagnitude, maxMagnitude, base, area, lossZone) + + if isset(conductance) && haskey(method.layout.wattmeter.bus, indexBus) + for indexWattmeter in method.layout.wattmeter.bus[indexBus] + method.mean[indexWattmeter] = (method.mean[indexWattmeter] + oldConductance) - bus.shunt.conductance[indexBus] + end + end + + if bus.layout.slack == indexBus && isset(angle) + for i = (method.layout.wattmeter.number + 1):method.layout.number + method.mean[i] = (method.mean[i] + oldSlackAngle) - bus.voltage.angle[bus.layout.slack] + end + end +end + """ @bus(kwargs...) diff --git a/src/powerSystem/generator.jl b/src/powerSystem/generator.jl index d751fc867..456ca3ce9 100644 --- a/src/powerSystem/generator.jl +++ b/src/powerSystem/generator.jl @@ -240,6 +240,21 @@ function addGenerator!(system::PowerSystem, analysis::ACOptimalPowerFlow; end end +function addGenerator!(system::PowerSystem, analysis::DCStateEstimationWLS; + label::L = missing, bus::L, area::T = missing, status::T = missing, + active::T = missing, reactive::T = missing, magnitude::T = missing, + minActive::T = missing, maxActive::T = missing, minReactive::T = missing, + maxReactive::T = missing, lowActive::T = missing, minLowReactive::T = missing, + maxLowReactive::T = missing, upActive::T = missing, minUpReactive::T = missing, + maxUpReactive::T = missing, loadFollowing::T = missing, reserve10min::T = missing, + reserve30min::T = missing, reactiveTimescale::T = missing) + + addGenerator!(system; label, bus, area, status, active, reactive, magnitude, + minActive, maxActive, minReactive, maxReactive, lowActive, minLowReactive, + maxLowReactive, upActive, minUpReactive, maxUpReactive, loadFollowing, reserve10min, + reserve30min, reactiveTimescale) +end + """ updateGenerator!(system::PowerSystem, analysis::Analysis; kwargs...) @@ -670,6 +685,32 @@ function updateGenerator!(system::PowerSystem, analysis::ACOptimalPowerFlow; end end +function updateGenerator!(system::PowerSystem, analysis::DCStateEstimationWLS; + label::L, area::T = missing, status::T = missing, + active::T = missing, reactive::T = missing, magnitude::T = missing, + minActive::T = missing, maxActive::T = missing, minReactive::T = missing, + maxReactive::T = missing, lowActive::T = missing, minLowReactive::T = missing, + maxLowReactive::T = missing, upActive::T = missing, minUpReactive::T = missing, + maxUpReactive::T = missing, loadFollowing::T = missing, reserve10min::T = missing, + reserve30min::T = missing, reactiveTimescale::T = missing) + + generator = system.generator + + if isset(status) + checkStatus(status) + index = generator.label[getLabel(generator, label, "generator")] + indexBus = generator.layout.bus[index] + if status == 0 && system.bus.layout.slack == indexBus && length(system.bus.supply.generator[indexBus]) == 1 + throw(ErrorException("The DC state estimation model cannot be reused due to required bus type conversion.")) + end + end + + updateGenerator!(system; label, area, status, active, reactive, magnitude, + minActive, maxActive, minReactive, maxReactive, lowActive, minLowReactive, + maxLowReactive, upActive, minUpReactive, maxUpReactive, loadFollowing, reserve10min, + reserve30min, reactiveTimescale) +end + """ @generator(kwargs...) diff --git a/src/stateEstimation/badData.jl b/src/stateEstimation/badData.jl index fb09046fb..31573c239 100644 --- a/src/stateEstimation/badData.jl +++ b/src/stateEstimation/badData.jl @@ -52,7 +52,7 @@ while analysis.bad.detect end ``` """ -function badData!(system::PowerSystem, device::Measurement, analysis::DCStateEstimation; threshold = 3.0) +function badData!(system::PowerSystem, device::Measurement, analysis::DCStateEstimationWLS; threshold = 3.0) bus = system.bus se = analysis.method bad = analysis.bad diff --git a/src/stateEstimation/dcStateEstimation.jl b/src/stateEstimation/dcStateEstimation.jl index a929ea494..7e1f2ef8a 100644 --- a/src/stateEstimation/dcStateEstimation.jl +++ b/src/stateEstimation/dcStateEstimation.jl @@ -66,10 +66,22 @@ analysis = dcStateEstimation(system, device, Ipopt.Optimizer) ``` """ function dcStateEstimation(system::PowerSystem, device::Measurement, factorization::Type{<:Union{QR, LDLt, LU}} = LU) - layout = DCStateEstimationLayout(Int64[], 0, 0, 0) + layout = DCStateEstimationLayout( + WattmeterLayout( + Dict{Int64,Array{Int64,1}}(), + Dict{Int64,Array{Int64,1}}(), + Dict{Int64,Array{Int64,1}}(), + device.wattmeter.number + ), + PMULayout( + Dict{Int64,Array{Int64,1}}(), + device.pmu.number + ), + device.wattmeter.number + device.pmu.number + ) + jacobian, weight, mean, layout = dcStateEstimationModel(system, device, layout) - numberDevice = lastindex(mean) method = Dict(LU => lu, LDLt => ldlt, QR => qr) return DCStateEstimationWLS( PolarAngle(Float64[]), @@ -94,28 +106,61 @@ end function dcStateEstimation(system::PowerSystem, device::Measurement, (@nospecialize optimizerFactory); bridge::Bool = true, name::Bool = true) - layout = DCStateEstimationLayout(Int64[], 0, 0, 0) + layout = DCStateEstimationLayout( + WattmeterLayout( + Dict{Int64,Array{Int64,1}}(), + Dict{Int64,Array{Int64,1}}(), + Dict{Int64,Array{Int64,1}}(), + device.wattmeter.number + ), + PMULayout( + Dict{Int64,Array{Int64,1}}(), + device.pmu.number + ), + device.wattmeter.number + device.pmu.number + ) jacobian, weight, mean, layout = dcStateEstimationModel(system, device, layout) bus = system.bus + pmu = device.pmu + jump = JuMP.Model(optimizerFactory; add_bridges = bridge) set_string_names_on_creation(jump, name) anglex = @variable(jump, 0 <= anglex[i = 1:bus.number]) angley = @variable(jump, 0 <= angley[i = 1:bus.number]) - residualx = @variable(jump, 0 <= residualx[i = 1:layout.device]) - residualy = @variable(jump, 0 <= residualy[i = 1:layout.device]) + residualx = @variable(jump, 0 <= residualx[i = 1:layout.number]) + residualy = @variable(jump, 0 <= residualy[i = 1:layout.number]) fix(anglex[bus.layout.slack], 0.0; force = true) - fix(angley[bus.layout.slack], bus.voltage.angle[bus.layout.slack]; force = true) + fix(angley[bus.layout.slack], 0.0; force = true) + objective = @expression(jump, AffExpr()) residual = Dict{Int64, JuMP.ConstraintRef}() angleJacobian = jacobian * (angley - anglex) - for i = 1:layout.device - residual[i] = @constraint(jump, angleJacobian[i] + residualy[i] - residualx[i] - mean[i] == 0.0) + for i = 1:layout.wattmeter.number + if device.wattmeter.active.status[i] == 1 + add_to_expression!(objective, residualx[i] + residualy[i]) + residual[i] = @constraint(jump, angleJacobian[i] + residualy[i] - residualx[i] - mean[i] == 0.0) + else + fix(residualx[i], 0.0; force = true) + fix(residualy[i], 0.0; force = true) + end end - objective = @objective(jump, Min, sum(residualx[i] + residualy[i] for i = 1:layout.device)) + for (i, k) in enumerate(layout.wattmeter.number + 1:layout.number) + if pmu.layout.bus[i] + if pmu.angle.status[i] == 1 + add_to_expression!(objective, residualx[k] + residualy[k]) + residual[k] = @constraint(jump, angleJacobian[k] + residualy[k] - residualx[k] - mean[k] == 0.0) + else + fix(residualx[k], 0.0; force = true) + fix(residualy[k], 0.0; force = true) + end + end + end + + @objective(jump, Min, objective) return DCStateEstimationLAV( PolarAngle(Float64[]), @@ -132,7 +177,6 @@ function dcStateEstimation(system::PowerSystem, device::Measurement, (@nospecial jump, VariableLAV(anglex, angley, residualx, residualy), residual, - objective, layout ), BadData(true, 0.0, 0, "") @@ -185,7 +229,7 @@ function solve!(system::PowerSystem, analysis::DCStateEstimationLAV) end for i = 1:system.bus.number - analysis.voltage.angle[i] = value(se.variable.angley[i]::JuMP.VariableRef) - value(se.variable.anglex[i]::JuMP.VariableRef) + analysis.voltage.angle[i] = value(se.variable.angley[i]::JuMP.VariableRef) - value(se.variable.anglex[i]::JuMP.VariableRef) + bus.voltage.angle[bus.layout.slack] end end @@ -240,88 +284,98 @@ function dcStateEstimationModel(system::PowerSystem, device::Measurement, layout end nonZeroElement = 0 - @inbounds for (i, index) in enumerate(wattmeter.layout.index) - if wattmeter.active.status[i] == 1 - if wattmeter.layout.bus[i] - nonZeroElement += (dc.nodalMatrix.colptr[index + 1] - dc.nodalMatrix.colptr[index]) - layout.device += 1 - else - nonZeroElement += 2 - layout.device += 1 - end + for (i, index) in enumerate(wattmeter.layout.index) + if wattmeter.layout.bus[i] + nonZeroElement += (dc.nodalMatrix.colptr[index + 1] - dc.nodalMatrix.colptr[index]) + layout.wattmeter.bus[index] = Int64[] + elseif wattmeter.layout.from[i] + nonZeroElement += 2 + layout.wattmeter.from[index] = Int64[] + else + nonZeroElement += 2 + layout.wattmeter.to[index] = Int64[] end end - layout.wattmeter = copy(layout.device) - @inbounds for i = 1:pmu.number - if pmu.layout.bus[i] && pmu.angle.status[i] == 1 + for i = 1:pmu.number + if pmu.layout.bus[i] nonZeroElement += 1 - layout.device += 1 + layout.pmu.bus[pmu.layout.index[i]] = Int64[] end end - layout.pmu = layout.device - layout.wattmeter row = fill(0, nonZeroElement) col = similar(row) jac = fill(0.0, nonZeroElement) - mean = fill(0.0, layout.device) + mean = fill(0.0, layout.number) weight = similar(mean) - layout.index = fill(0, layout.device) count = 1 rowindex = 1 - @inbounds for (i, index) in enumerate(wattmeter.layout.index) + for (i, index) in enumerate(wattmeter.layout.index) if wattmeter.active.status[i] == 1 - if wattmeter.layout.bus[i] - for j in dc.nodalMatrix.colptr[index]:(dc.nodalMatrix.colptr[index + 1] - 1) - row[count] = rowindex - col[count] = dc.nodalMatrix.rowval[j] - jac[count] = dc.nodalMatrix.nzval[j] - count += 1 - end - mean[rowindex] = wattmeter.active.mean[i] - dc.shiftPower[index] - bus.shunt.conductance[index] - weight[rowindex] = 1 / wattmeter.active.variance[i] - layout.index[rowindex] = i - - rowindex += 1 - else - if wattmeter.layout.from[i] - addmitance = dc.admittance[index] - else - addmitance = -dc.admittance[index] - end + constIf = 1 + else + constIf = 0 + end + if wattmeter.layout.bus[i] + for j in dc.nodalMatrix.colptr[index]:(dc.nodalMatrix.colptr[index + 1] - 1) row[count] = rowindex - col[count] = branch.layout.from[index] - jac[count] = addmitance + col[count] = dc.nodalMatrix.rowval[j] + jac[count] = dc.nodalMatrix.nzval[j] * constIf count += 1 - row[count] = rowindex - col[count] = branch.layout.to[index] - jac[count] = -addmitance - - mean[rowindex] = wattmeter.active.mean[i] + branch.parameter.shiftAngle[index] * addmitance - weight[rowindex] = 1 / wattmeter.active.variance[i] - layout.index[rowindex] = i - - count += 1; rowindex += 1 end + mean[rowindex] = (wattmeter.active.mean[i] - dc.shiftPower[index] - bus.shunt.conductance[index]) * constIf + weight[rowindex] = (1 / wattmeter.active.variance[i]) * constIf + push!(layout.wattmeter.bus[index], i) + + rowindex += 1 + else + if wattmeter.layout.from[i] + addmitance = dc.admittance[index] * constIf + push!(layout.wattmeter.from[index], i) + else + addmitance = -dc.admittance[index] * constIf + push!(layout.wattmeter.to[index], i) + end + + row[count] = rowindex + col[count] = branch.layout.from[index] + jac[count] = addmitance + count += 1 + row[count] = rowindex + col[count] = branch.layout.to[index] + jac[count] = -addmitance + + mean[rowindex] = (wattmeter.active.mean[i] + branch.parameter.shiftAngle[index] * addmitance) * constIf + weight[rowindex] = (1 / wattmeter.active.variance[i]) * constIf + + count += 1; rowindex += 1 end end - @inbounds for i = 1:pmu.number - if pmu.layout.bus[i] && pmu.angle.status[i] == 1 + for i = 1:pmu.number + if pmu.angle.status[i] == 1 + constIf = 1 + else + constIf = 0 + end + + if pmu.layout.bus[i] row[count] = rowindex col[count] = pmu.layout.index[i] - jac[count] = 1.0 + jac[count] = constIf - mean[rowindex] = pmu.angle.mean[i] - bus.voltage.angle[bus.layout.slack] - weight[rowindex] = 1 / pmu.angle.variance[i] - layout.index[rowindex] = i + mean[rowindex] = (pmu.angle.mean[i] - bus.voltage.angle[bus.layout.slack]) * constIf + weight[rowindex] = (1 / pmu.angle.variance[i]) * constIf + + push!(layout.pmu.bus[pmu.layout.index[i]], i) count += 1; rowindex += 1 end end - return sparse(row, col, jac, layout.device, bus.number), weight, mean, layout + return sparse(row, col, jac, layout.number, bus.number), weight, mean, layout end function deleteSlackJacobian(analysis::DCStateEstimation, slack::Int64) @@ -351,4 +405,14 @@ function restoreSlackJacobian(analysis::DCStateEstimation, slackRange::UnitRange @inbounds for (k, i) in enumerate(slackRange) se.jacobian[se.jacobian.rowval[i], slack] = elementsRemove[k] end +end + +function constMeter(flag::Union{Int8, Float64}) + if flag == 0 + constIf = 0.0 + else + constIf = 1.0 + end + + return constIf end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 45269897f..22507eef1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -44,4 +44,5 @@ include("measurement/buildUpdate.jl") ######## State Estimation ########## include("stateEstimation/analysis.jl") -include("stateEstimation/badData.jl") \ No newline at end of file +include("stateEstimation/reusing.jl") +# include("stateEstimation/badData.jl") \ No newline at end of file diff --git a/test/stateEstimation/analysis.jl b/test/stateEstimation/analysis.jl index 9f03210c4..3e90f382c 100644 --- a/test/stateEstimation/analysis.jl +++ b/test/stateEstimation/analysis.jl @@ -166,4 +166,69 @@ system30 = powerSystem(string(pathData, "case30test.m")) analysisQR = dcStateEstimation(system30, device, QR) solve!(system30, analysisQR) @test analysisLU.voltage.angle ≈ analysisQR.voltage.angle +end + +@testset "DC State Estimation: Incomplete Set" begin + @default(template) + @default(unit) + + ############### Modified IEEE 14-bus Test Case ################ + dcModel!(system14) + analysis = dcPowerFlow(system14) + solve!(system14, analysis) + power!(system14, analysis) + + deviceAll = measurement() + devicePart = measurement() + + ################ Wattmeters ################ + for (key, value) in system14.bus.label + if value in [1; 5; 8] + addWattmeter!(system14, deviceAll; bus = key, active = analysis.power.injection.active[value], noise = false, status = 0) + else + addWattmeter!(system14, deviceAll; bus = key, active = analysis.power.injection.active[value], noise = false) + addWattmeter!(system14, devicePart; bus = key, active = analysis.power.injection.active[value], noise = false) + end + end + + for (key, value) in system14.branch.label + if value in [4; 15; 18] + addWattmeter!(system14, deviceAll; from = key, active = analysis.power.from.active[value], noise = false, status = 0) + else + addWattmeter!(system14, deviceAll; from = key, active = analysis.power.from.active[value], noise = false) + addWattmeter!(system14, devicePart; from = key, active = analysis.power.from.active[value], noise = false) + end + if value in [4; 16; 19] + addWattmeter!(system14, deviceAll; to = key, active = analysis.power.to.active[value], noise = false, status = 0) + else + addWattmeter!(system14, deviceAll; to = key, active = analysis.power.to.active[value], noise = false) + addWattmeter!(system14, devicePart; to = key, active = analysis.power.to.active[value], noise = false) + end + end + + analysisAll = dcStateEstimation(system14, deviceAll) + solve!(system14, analysisAll) + + analysisPart = dcStateEstimation(system14, devicePart) + solve!(system14, analysisPart) + + @test analysisAll.voltage.angle ≈ analysisPart.voltage.angle + + ################ PMUs ################ + for (key, value) in system14.bus.label + if value in [1; 6; 10; 14] + addPmu!(system14, deviceAll; bus = key, magnitude = 1.0, angle = analysis.voltage.angle[value], noise = false, statusAngle = 0) + else + addPmu!(system14, deviceAll; bus = key, magnitude = 1.0, angle = analysis.voltage.angle[value], noise = false) + addPmu!(system14, devicePart; bus = key, magnitude = 1.0, angle = analysis.voltage.angle[value], noise = false) + end + end + + analysisAll = dcStateEstimation(system14, deviceAll) + solve!(system14, analysisAll) + + analysisPart = dcStateEstimation(system14, devicePart) + solve!(system14, analysisPart) + + @test analysisAll.voltage.angle ≈ analysisPart.voltage.angle end \ No newline at end of file diff --git a/test/stateEstimation/badData.jl b/test/stateEstimation/badData.jl index 9f56fef4d..48d11443d 100644 --- a/test/stateEstimation/badData.jl +++ b/test/stateEstimation/badData.jl @@ -29,7 +29,7 @@ system30 = powerSystem(string(pathData, "case30test.m")) addPmu!(system14, device; bus = key, magnitude = 1.0, angle = analysis.voltage.angle[value], noise = false) end - ####### WLS: One Outlier ####### + ####### WLS LU: One Outlier ####### updateWattmeter!(system14, device; label = "Wattmeter 2", active = 100, noise = false) analysisSE = dcStateEstimation(system14, device) solve!(system14, analysisSE) @@ -41,7 +41,7 @@ system30 = powerSystem(string(pathData, "case30test.m")) solve!(system14, analysisSE) @test analysis.voltage.angle ≈ analysisSE.voltage.angle - ####### WLS: Two Outliers ####### + ####### WLS LU: Two Outliers ####### updatePmu!(system14, device; label = "PMU 10", angle = 10pi, noise = false) analysisSE = dcStateEstimation(system14, device) solve!(system14, analysisSE) @@ -57,4 +57,33 @@ system30 = powerSystem(string(pathData, "case30test.m")) solve!(system14, analysisSE) @test analysis.voltage.angle ≈ analysisSE.voltage.angle + + ####### WLS QR: One Outlier ####### + updatePmu!(system14, device; label = "PMU 10", angle = analysis.voltage.angle[10], noise = false) + analysisSE = dcStateEstimation(system14, device, QR) + solve!(system14, analysisSE) + + badData!(system14, device, analysisSE; threshold = 3.0) + @test analysisSE.bad.label == "Wattmeter 2" + @test analysisSE.bad.maxNormalizedResidual ≈ 829.9 atol = 1e-1 + + solve!(system14, analysisSE) + @test analysis.voltage.angle ≈ analysisSE.voltage.angle + + ####### WLS QR: Two Outliers ####### + updatePmu!(system14, device; label = "PMU 10", angle = 10pi, noise = false) + analysisSE = dcStateEstimation(system14, device, QR) + solve!(system14, analysisSE) + + badData!(system14, device, analysisSE; threshold = 3.0) + @test analysisSE.bad.label == "PMU 10" + @test analysisSE.bad.maxNormalizedResidual ≈ 5186.3 atol = 1e-1 + + solve!(system14, analysisSE) + badData!(system14, device, analysisSE; threshold = 3.0) + @test analysisSE.bad.label == "Wattmeter 2" + @test analysisSE.bad.maxNormalizedResidual ≈ 829.9 atol = 1e-1 + + solve!(system14, analysisSE) + @test analysis.voltage.angle ≈ analysisSE.voltage.angle end \ No newline at end of file diff --git a/test/stateEstimation/reusing.jl b/test/stateEstimation/reusing.jl new file mode 100644 index 000000000..c200f0071 --- /dev/null +++ b/test/stateEstimation/reusing.jl @@ -0,0 +1,231 @@ +system14 = powerSystem(string(pathData, "case14test.m")) + +@testset "Reusing Meters DC State Estimation" begin + @default(template) + @default(unit) + + ############### Modified IEEE 14-bus Test Case ################ + updateBus!(system14; label = 1, type = 2) + updateBus!(system14; label = 3, type = 3, angle = -0.17) + + dcModel!(system14) + analysis = dcPowerFlow(system14) + solve!(system14, analysis) + power!(system14, analysis) + + ####### Measurements ####### + device = measurement() + + @wattmeter(label = "!") + for (key, value) in system14.bus.label + if value == 1 + addWattmeter!(system14, device; bus = key, active = rand(1)[], noise = false) + elseif value == 3 + addWattmeter!(system14, device; bus = key, active = analysis.power.injection.active[value], noise = false, status = 0) + elseif value == 5 + addWattmeter!(system14, device; bus = key, active = rand(1)[], variance = 1e-5) + elseif value == 9 + addWattmeter!(system14, device; bus = key, active = rand(1)[], noise = false) + else + addWattmeter!(system14, device; bus = key, active = analysis.power.injection.active[value], noise = false) + end + end + + for (key, value) in system14.branch.label + if value == 4 + addWattmeter!(system14, device; from = key, active = rand(1)[], noise = false) + elseif value == 15 + addWattmeter!(system14, device; from = key, active = analysis.power.from.active[value], noise = false, status = 0) + elseif value == 18 + addWattmeter!(system14, device; from = key, active = rand(1)[], variance = 1e-5) + elseif value == 20 + addWattmeter!(system14, device; from = key, active = rand(1)[]) + else + addWattmeter!(system14, device; from = key, active = analysis.power.from.active[value], noise = false) + end + + if value == 5 + addWattmeter!(system14, device; to = key, active = rand(1)[], noise = false) + elseif value == 8 + addWattmeter!(system14, device; to = key, active = analysis.power.to.active[value], noise = false, status = 0) + elseif value == 11 + addWattmeter!(system14, device; to = key, active = rand(1)[], variance = 1e-5) + elseif value == 19 + addWattmeter!(system14, device; to = key, active = rand(1)[]) + else + addWattmeter!(system14, device; to = key, active = analysis.power.to.active[value], noise = false) + end + end + + for (key, value) in system14.bus.label + if value == 2 + addPmu!(system14, device; bus = key, magnitude = 1, angle = rand(1)[], noise = false) + elseif value == 6 + addPmu!(system14, device; bus = key, magnitude = 1, angle = analysis.voltage.angle[value], noise = false, statusAngle = 0) + elseif value == 9 + addPmu!(system14, device; bus = key, magnitude = 1, angle = rand(1)[], varianceAngle = 1e-5) + elseif value == 13 + addPmu!(system14, device; bus = key, magnitude = 1, angle = rand(1)[], noise = false) + else + addPmu!(system14, device; bus = key, magnitude = 1, angle = analysis.voltage.angle[value], noise = false) + end + end + + ####### Original WLS and LAV Models ####### + analysisWLS = dcStateEstimation(system14, device) + analysisLAV = dcStateEstimation(system14, device, Ipopt.Optimizer) + + ####### Update Only Devices ####### + updateWattmeter!(system14, device; label = 1, status = 0) + updateWattmeter!(system14, device; label = 3, status = 1) + updateWattmeter!(system14, device; label = 5, active = analysis.power.injection.active[5], variance = 1e-2, noise = false) + updateWattmeter!(system14, device; label = 9, active = analysis.power.injection.active[9], noise = false) + + updateWattmeter!(system14, device; label = "From 4", status = 0) + updateWattmeter!(system14, device; label = "From 15", status = 1) + updateWattmeter!(system14, device; label = "From 18", active = analysis.power.from.active[18], variance = 1e-2, noise = false) + updateWattmeter!(system14, device; label = "From 20", active = analysis.power.from.active[20], noise = false) + + updateWattmeter!(system14, device; label = "To 5", status = 0) + updateWattmeter!(system14, device; label = "To 8", status = 1) + updateWattmeter!(system14, device; label = "To 11", active = analysis.power.to.active[11], variance = 1e-2, noise = false) + updateWattmeter!(system14, device; label = "To 19", active = analysis.power.to.active[19], noise = false) + + updatePmu!(system14, device; label = 2, statusAngle = 0) + updatePmu!(system14, device; label = 6, statusAngle = 1) + updatePmu!(system14, device; label = 9, angle = analysis.voltage.angle[9], varianceAngle = 1e-5, noise = false) + updatePmu!(system14, device; label = 13, angle = analysis.voltage.angle[13], noise = false) + + ####### Solve Updated WLS and LAV Models ####### + analysisWLSUpdate = dcStateEstimation(system14, device) + solve!(system14, analysisWLSUpdate) + @test analysisWLSUpdate.voltage.angle ≈ analysis.voltage.angle + + analysisLAVUpdate = dcStateEstimation(system14, device, Ipopt.Optimizer) + JuMP.set_silent(analysisLAVUpdate.method.jump) + solve!(system14, analysisLAVUpdate) + @test analysisLAVUpdate.voltage.angle ≈ analysis.voltage.angle + + ###### Update Devices and Original WLS Model ####### + updateWattmeter!(system14, device, analysisWLS; label = 1, status = 0) + updateWattmeter!(system14, device, analysisWLS; label = 3, status = 1) + updateWattmeter!(system14, device, analysisWLS; label = 5, active = analysis.power.injection.active[5], variance = 1e-2, noise = false) + updateWattmeter!(system14, device, analysisWLS; label = 9, active = analysis.power.injection.active[9], noise = false) + + updateWattmeter!(system14, device, analysisWLS; label = "From 4", status = 0) + updateWattmeter!(system14, device, analysisWLS; label = "From 15", status = 1) + updateWattmeter!(system14, device, analysisWLS; label = "From 18", active = analysis.power.from.active[18], variance = 1e-2, noise = false) + updateWattmeter!(system14, device, analysisWLS; label = "From 20", active = analysis.power.from.active[20], noise = false) + + updateWattmeter!(system14, device, analysisWLS; label = "To 5", status = 0) + updateWattmeter!(system14, device, analysisWLS; label = "To 8", status = 1) + updateWattmeter!(system14, device, analysisWLS; label = "To 11", active = analysis.power.to.active[11], variance = 1e-2, noise = false) + updateWattmeter!(system14, device, analysisWLS; label = "To 19", active = analysis.power.to.active[19], noise = false) + + updatePmu!(system14, device, analysisWLS; label = 2, statusAngle = 0) + updatePmu!(system14, device, analysisWLS; label = 6, statusAngle = 1) + updatePmu!(system14, device, analysisWLS; label = 9, angle = analysis.voltage.angle[9], varianceAngle = 1e-5, noise = false) + updatePmu!(system14, device, analysisWLS; label = 13, angle = analysis.voltage.angle[13], noise = false) + + updateWattmeter!(system14, device, analysisWLS; label = 4, status = 0) + updateWattmeter!(system14, device, analysisWLS; label = 4, status = 1) + updateWattmeter!(system14, device, analysisWLS; label = "From 2", status = 0) + updateWattmeter!(system14, device, analysisWLS; label = "From 2", status = 1) + updateWattmeter!(system14, device, analysisWLS; label = "To 12", status = 0) + updateWattmeter!(system14, device, analysisWLS; label = "To 12", status = 1) + updatePmu!(system14, device, analysisWLS; label = 10, statusAngle = 0) + updatePmu!(system14, device, analysisWLS; label = 10, statusAngle = 1) + + solve!(system14, analysisWLS) + @test analysisWLS.voltage.angle ≈ analysis.voltage.angle + + ##### Update Devices and Original LAV Model ####### + updateWattmeter!(system14, device, analysisLAV; label = 1, status = 0) + updateWattmeter!(system14, device, analysisLAV; label = 3, status = 1) + updateWattmeter!(system14, device, analysisLAV; label = 5, active = analysis.power.injection.active[5], variance = 1e-2, noise = false) + updateWattmeter!(system14, device, analysisLAV; label = 9, active = analysis.power.injection.active[9], noise = false) + + updateWattmeter!(system14, device, analysisLAV; label = "From 4", status = 0) + updateWattmeter!(system14, device, analysisLAV; label = "From 15", status = 1) + updateWattmeter!(system14, device, analysisLAV; label = "From 18", active = analysis.power.from.active[18], variance = 1e-2, noise = false) + updateWattmeter!(system14, device, analysisLAV; label = "From 20", active = analysis.power.from.active[20], noise = false) + + updateWattmeter!(system14, device, analysisLAV; label = "To 5", status = 0) + updateWattmeter!(system14, device, analysisLAV; label = "To 8", status = 1) + updateWattmeter!(system14, device, analysisLAV; label = "To 11", active = analysis.power.to.active[11], variance = 1e-2, noise = false) + updateWattmeter!(system14, device, analysisLAV; label = "To 19", active = analysis.power.to.active[19], noise = false) + + updatePmu!(system14, device, analysisLAV; label = 2, statusAngle = 0) + updatePmu!(system14, device, analysisLAV; label = 6, statusAngle = 1) + updatePmu!(system14, device, analysisLAV; label = 9, angle = analysis.voltage.angle[9], varianceAngle = 1e-5, noise = false) + updatePmu!(system14, device, analysisLAV; label = 13, angle = analysis.voltage.angle[13], noise = false) + + updateWattmeter!(system14, device, analysisLAV; label = 4, status = 0) + updateWattmeter!(system14, device, analysisLAV; label = 4, status = 1) + updateWattmeter!(system14, device, analysisLAV; label = "From 2", status = 0) + updateWattmeter!(system14, device, analysisLAV; label = "From 2", status = 1) + updateWattmeter!(system14, device, analysisLAV; label = "To 12", status = 0) + updateWattmeter!(system14, device, analysisLAV; label = "To 12", status = 1) + updatePmu!(system14, device, analysisLAV; label = 10, statusAngle = 0) + updatePmu!(system14, device, analysisLAV; label = 10, statusAngle = 1) + + JuMP.set_silent(analysisLAV.method.jump) + solve!(system14, analysisLAV) + @test analysisLAV.voltage.angle ≈ analysis.voltage.angle +end + + +system14 = powerSystem(string(pathData, "case14test.m")) +@testset "Reusing Power System DC State Estimation" begin + ############### Modified IEEE 14-bus Test Case ################ + dcModel!(system14) + powerFlow = dcPowerFlow(system14) + solve!(system14, powerFlow) + power!(system14, powerFlow) + + ############### Measurements ################ + device = measurement() + + for (key, value) in system14.bus.label + addWattmeter!(system14, device; bus = key, active = powerFlow.power.injection.active[value], noise = false) + addPmu!(system14, device; bus = key, magnitude = 1.0, angle = powerFlow.voltage.angle[value], noise = false) + end + for (key, value) in system14.branch.label + addWattmeter!(system14, device; from = key, active = powerFlow.power.from.active[value], noise = false) + addWattmeter!(system14, device; to = key, active = powerFlow.power.to.active[value], noise = false) + end + addWattmeter!(system14, device; bus = 2, active = rand(1)[]) + addWattmeter!(system14, device; from = 1, active = rand(1)[]) + addWattmeter!(system14, device; from = 1, active = rand(1)[]) + addWattmeter!(system14, device; to = 2, active = rand(1)[]) + addWattmeter!(system14, device; to = 2, active = rand(1)[]) + addPmu!(system14, device; bus = 2, magnitude = 1.0, angle = rand(1)[]) + + ############### Solution ################ + updateBus!(system14; label = 2, conductance = 10) + updateBus!(system14; label = 1, angle = -0.17) + updateBranch!(system14; label = 1, reactance = 15, shiftAngle = 2.1, turnsRatio = 0.8) + updateBranch!(system14; label = 2, status = 0) + updateBranch!(system14; label = 18, status = 1) + analysisCon = dcStateEstimation(system14, device) + solve!(system14, analysisCon) + + + ############### Update Power System ################ + updateBus!(system14; label = 2, conductance = 0) + updateBus!(system14; label = 1, angle = 0.0) + updateBranch!(system14; label = 1, reactance = 1, shiftAngle = 0.0, turnsRatio = 1.0) + updateBranch!(system14; label = 2, status = 1) + updateBranch!(system14; label = 18, status = 0) + analysis = dcStateEstimation(system14, device) + + updateBus!(system14, analysis; label = 2, conductance = 10) + updateBus!(system14, analysis; label = 1, angle = -0.17) + updateBranch!(system14, analysis; label = 1, reactance = 15, shiftAngle = 2.1, turnsRatio = 0.8) + updateBranch!(system14, analysis; label = 2, status = 0) + updateBranch!(system14, analysis; label = 18, status = 1) + + solve!(system14, analysis) + + @test analysis.voltage.angle ≈ analysisCon.voltage.angle +end \ No newline at end of file