Skip to content
GitLab
Commit d01c0b9b authored by martin.moraga's avatar martin.moraga
Browse files

improve nv_state_estimator_cim.py

parent 347be7ee
Hide whitespace changes
Inline Side-by-side
acs/state_estimation/measurement.py
View file @ d01c0b9b
......@@ -2,17 +2,22 @@ from enum import Enum
import numpy as np
class ElemType(Enum):
Node = 1 #Node Voltage
Branch = 2 #Complex Power Injection at node
Node = 1 #Node Voltage
Branch = 2 #Complex Power flow at branch
class MeasType(Enum):
V = 1 #Node Voltage
Sinj= 2 #Complex Power Injection at node
S1 = 4 #Complex Power flow at branch, measured at initial node
S2 = 5 #Complex Power flow at branch, measured at final node
I = 6 #Branch Current
Vpmu = 7 #Node Voltage
Ipmu = 9 #Branch Current
V_mag = 1 #Node Voltage
Sinj_real= 2 #Complex Power Injection at node
Sinj_imag= 3 #Complex Power Injection at node
S1_real = 4 #Active Power flow at branch, measured at initial node (S1.real)
S1_imag = 5 #Reactive Power flow at branch, measured at initial node (S1.imag)
I_mag = 6 #Branch Current
Vpmu_mag = 7 #Node Voltage
Vpmu_phase = 8 #Node Voltage
Ipmu_mag = 9 #Branch Current
Ipmu_phase = 10 #Branch Current
S2_real = 11 #Active Power flow at branch, measured at final node (S2.real)
S2_imag = 12 #Reactive Power flow at branch, measured at final node (S2.imag)
class Measurement():
def __init__(self, element, element_type, meas_type, meas_value, std_dev):
......@@ -35,16 +40,18 @@ class Measurement():
self.element_type = element_type
self.meas_type = meas_type
self.meas_value = meas_value
self.std_dev = self.meas_value*(std_dev/300)
if self.std_dev<10**(-6)
self.std_dev = 10**(-6)
self.mval = 0.0 #measured values (affected by uncertainty)
self.std_dev = std_dev
self.std_dev = std_dev
self.mval = 0.0 #measured values (affected by uncertainty)
class Measurents_set():
def __init__(self):
self.measurements = []
self.measurements = [] #array with all measurements
def create_measurement(self, element, element_type, meas_type, meas_value, std_dev):
"""
to add elements to the measurements array
"""
self.measurements.append(Measurement(element, element_type, meas_type, meas_value, std_dev))
def meas_creation(self):
......@@ -54,79 +61,114 @@ class Measurents_set():
err_pu = np.random.normal(0,1,len(self.measurements))
for index, measurement in enumerate(self.measurements):
measurement.mval = measurement.meas_value + self.std_dev*err_pu[index]
def getNumberOfMeasurements(self)
"""
return number of measurements of each type in the array Measurents_set.measurements
"""
nvi, npi, nqi, npf, nqf, nii, nvpum, nipmu = 0
for elem in self.measurements:
if elem.meas_type is MeasType.V:
nvi = nvi+1
elif elem.meas_type is MeasType.Sinj:
npi = npi+1
nqi = nqi+1
elif elem.meas_type is MeasType.S1 or elem.meas_type is MeasType.S2:
npf = npf+1
nqf = nqf+1
elif elem.meas_type is MeasType.I:
nii = nii+1
elif elem.meas_type is MeasType.Vpmu:
nvpum = nvpum+1
elif elem.meas_type is MeasType.Ipmu:
nipmu = nipmu+1
return nvi, npi, nqi, npf, nqf, nii, nvpum, nipmu
def getMeasuredActiveInjPowers(self):
def meas_creation_test(self, err_pu):
"""
For test purposes.
It calculates the measured values (affected by uncertainty) at the measurement points.
This function takes as paramenter the random gaussian distribution.
"""
for index, measurement in enumerate(self.measurements):
measurement.mval = measurement.meas_value + measurement.std_dev*err_pu[index]
def getMeasurements(self, type):
"""
return an array with the measurements of type Sinj.real
return an array with all measurements of type "type" in the array Measurents_set.measurements.
"""
Pinj = np.array([])
for elem in self.measurements:
if elem.meas_type is MeasType.Sinj:
Pinj = np.append(Pinj, elem.real)
measurements = []
for measurement in self.measurements:
if measurement.meas_type is type:
measurements.append(measurement)
return Pinj
return measurements
def getNumberOfMeasurements(self, type):
"""
return number of measurements of type "type" in the array Measurents_set.measurements
"""
number = 0
for measurement in self.measurements:
if measurement.meas_type is type:
number = number+1
return number
def getMeasuredReactiveInjPowers(self):
def getIndexOfMeasurements(self, type):
"""
return an array with the measurements of type Sinj.imag
return index of all measurements of type "type" in the array Measurents_set.measurements
"""
Qinj = np.array([])
for elem in self.measurements:
if elem.meas_type is MeasType.Sinj:
Qinj = np.append(Qinj, elem.imag)
idx = np.zeros(self.getNumberOfMeasurements(type), dtype=int)
i = 0
for index, measurement in enumerate(self.measurements):
if measurement.meas_type is type:
idx[i] = index
i = i+1
return idx
return Qinj
def getWeightsMatrix(self):
"""
return an array the weights (obtained as standard_deviations^-2)
"""
weights = np.zeros(len(self.measurements))
for index, measurement in enumerate(self.measurements):
#the weight is small and can bring instability during matrix inversion, so we "cut" everything below 10^-6
if measurement.std_dev<10**(-6):
measurement.std_dev = 10**(-6)
weights[index] = measurement.std_dev**(-2)
def getMeasuredActiveBPowers(self):
return weights
def getMeasValues(self):
"""
return an array with the measurements of type S1.real or S2.real
for test purposes
returns an array with all measured values
"""
Pbr = np.array([])
for elem in self.measurements:
if elem.meas_type is MeasType.S1 or elem.meas_type is MeasType.S2:
Pbr = np.append(Pbr, elem.real)
return Pbr
meas_val = np.zeros(len(self.measurements))
for index, measurement in enumerate(self.measurements):
meas_val[index] = measurement.meas_value
return meas_val
def getMeasuredReactiveBPowers(self):
def getmVal(self):
"""
return an array with the measurements of type S1.imag or S2.imag
returns an array with all measured values (affected by uncertainty)
"""
Qbr = np.array([])
for elem in self.measurements:
if elem.meas_type is MeasType.S1 or elem.meas_type is MeasType.S2:
Qbr = np.append(Qbr, elem.imag)
return Qbr
mVal = np.zeros(len(self.measurements))
for index, measurement in enumerate(self.measurements):
mVal[index] = measurement.mval
""" Replace in mVal amplitude and phase of Vpmu by real and imaginary part """
#get all measurements of type MeasType.Vpmu_mag
Vpmu_mag_idx = self.getIndexOfMeasurements(type=MeasType.Vpmu_mag)
#get all measurements of type MeasType.Vpmu_phase
Vpmu_phase_idx = self.getIndexOfMeasurements(type=MeasType.Vpmu_phase)
for vpmu_mag_index, vpmu_phase_index in zip(Vpmu_mag_idx, Vpmu_phase_idx):
vamp = self.measurements[vpmu_mag_index].mval
vtheta = self.measurements[vpmu_phase_index].mval
mVal[vpmu_mag_index] = vamp*np.cos(vtheta)
mVal[vpmu_phase_index] = vamp*np.sin(vtheta)
def getWeightsMatrix(self)
""" Replace in z amplitude and phase of Ipmu by real and imaginary part """
#get all measurements of type MeasType.Ipmu_mag
Ipmu_mag_idx = self.getIndexOfMeasurements(type=MeasType.Ipmu_mag)
#get all measurements of type MeasType.Ipmu_phase
Ipmu_phase_idx = self.getIndexOfMeasurements(type=MeasType.Ipmu_phase)
for ipmu_mag_index, ipmu_phase_index in zip(Ipmu_mag_idx, Ipmu_phase_idx):
iamp = self.measurements[ipmu_mag_index].mval
itheta = self.measurements[ipmu_phase_index].mval
mVal[ipmu_mag_index] = iamp*np.cos(itheta)
mVal[ipmu_phase_index] = iamp*np.sin(itheta)
return mVal
def getStd_Dev(self):
"""
creates the weights matrix (obtained as standard_deviations^-2)
for test purposes
returns an array with all standard deviations
"""
weights = np.zeros(len(self.measurements))
std_dev = np.zeros(len(self.measurements))
for index, measurement in enumerate(self.measurements):
weights[index] = measurement.std_dev**(-2)
std_dev[index] = measurement.std_dev
return np.diag(weights)
\ No newline at end of file
return std_dev
\ No newline at end of file
acs/state_estimation/nv_powerflow_cim.py
View file @ d01c0b9b
......@@ -167,7 +167,7 @@ class PowerflowResults():
for branch in self.branches:
S2 = np.append(S2, branch.power2)
return S2
def solve(system):
"""It performs Power Flow by using rectangular node voltage state variables."""
......@@ -254,59 +254,4 @@ def solve(system):
results.calculateS1()
results.calculateS2()
return results, num_iter
def calculateI(system, V):
"""
To calculate the branch currents
"""
Ymatrix, Adj = Ymatrix_calc(system)
I = np.zeros((len(system.branches)), dtype=np.complex)
for idx in range(len(system.branches)):
fr = system.branches[idx].start_node.index
to = system.branches[idx].end_node.index
I[idx] = - (V[fr] - V[to])*Ymatrix[fr][to]
return I
def calculateInj(system, I):
"""
To calculate current injections at a node
"""
Iinj = np.zeros((len(system.nodes)), dtype=np.complex)
for k in range(0, (len(system.nodes))):
to=[]
fr=[]
for m in range(len(system.branches)):
if k==system.branches[m].start_node.index:
fr.append(m)
if k==system.branches[m].end_node.index:
to.append(m)
Iinj[k] = np.sum(I[to]) - np.sum(I[fr])
return Iinj
def calculateS1(system, V, I):
"""
To calculate powerflow on branches
"""
S1 = np.zeros((len(system.branches)), dtype=np.complex)
for i in range(0, len(system.branches)):
S1[i] = V[system.branches[i].start_node.index]*(np.conj(I[i]))
return S1
def calculateS2(system, V, I):
"""
To calculate powerflow on branches
"""
S2 = np.zeros((len(system.branches)), dtype=np.complex)
for i in range(0, len(system.branches)):
S2[i] = -V[system.branches[i].end_node.index]*(np.conj(I[i]))
return S2
def calculateSinj(V, Iinj):
"""
To calculate power injection at a node
"""
Sinj = np.multiply(V, np.conj(Iinj))
return Sinj
\ No newline at end of file
return results, num_iter
\ No newline at end of file
acs/state_estimation/nv_state_estimator_cim.py
View file @ d01c0b9b
import sys
import numpy
import nv_powerflow_cim
import numpy as np
from results import Results
from measurement import *
sys.path.append("../../../dataprocessing")
from villas.dataprocessing.readtools import read_timeseries_dpsim
def DsseCall(system, zdata, measurements):
""" It identifies the type of measurements present in the measurement set and
calls the appropriate estimator for dealing with them."""
# system: model of the system (nodes, lines, topology)
# zdata: Vector of measurements in Input (voltages, currents, powers)
# 1. Select type of Estimator.
# If at least a PMU is present we launch the combined estimator, otherwise a simple traditional etimator
Vmag_meas = 0
Vpmu_meas = 0
for elem in measurements.measurements_set:
if elem.meas_type=="v":
Vmag_meas = Vmag_meas + 1
elif elem.meas_type=="Vpmu":
Vpmu_meas = Vpmu_meas + 1
trad_code = 1 if Vmag_meas >0 else 0
PMU_code = 2 if Vpmu_meas >0 else 0
est_code = trad_code + PMU_code
# 2. Run Estimator.
if est_code == 1:
Vest = DsseTrad(system, zdata)
elif est_code == 2:
Vest = DssePmu(system, zdata)
else:
Vest = DsseMixed(system, zdata)
# 3. Populate arrays with data in output of estimator
# The Estimator provides:
# Vest: Node Voltage phasor estimated
# Iest: Branch Current phasor estimated
# Iinjest: Injection Current phasor estimated
# S1est: Complex Power flow at branch, measured at initial node
# S2est: Complex Power flow at branch, measured at final node
# Sinjest: Complex Power Injection at node
""" From here on, all the other quantities of the grid are calculated """
results = nv_powerflow_cim.PowerflowResults(system)
results.load_voltages(Vest)
results.calculateI()
results.calculateIinj()
results.calculateSinj()
results.calculateI()
results.calculateS()
return results
def DsseTrad(system, zdata):
""" It performs state estimation using rectangular node voltage state variables
and it is customized to work without PMU measurements"""
# Traditional state estimator
# system: model of the system (nodes, lines, topology)
# zdata: Vector of measurements in Input (voltages, currents, powers)
# number of nodes of the grids, identify also the number of states (2*nodes_num-1)
nodes_num = len(system.nodes)
vidx = numpy.where(zdata.mtype==1) #voltage input measurement
pidx = numpy.where(zdata.mtype==2) #active power input measurement
qidx = numpy.where(zdata.mtype==3) #reactive power input measurement
pfidx = numpy.where(zdata.mtype==4) #active power flow input measurement
qfidx = numpy.where(zdata.mtype==5) #reactive power flow input measurement
iidx = numpy.where(zdata.mtype==6) #current flow input measurement
nvi = len(vidx[0]) # number v.measure
npi = len(pidx[0]) # number power.act.measure
nqi = len(qidx[0]) # number power.react.measure
npf = len(pfidx[0]) # number power.act.flow.measure
nqf = len(qfidx[0]) # number power.react.flow.measure
nii = len(iidx[0]) #number c.flow.measure
busvi = zdata.mfrom[vidx] #node where v.measure is taken
buspi = zdata.mfrom[pidx] #node where power.act.measure is taken
busqi = zdata.mfrom[qidx] #node where power.react.measure is taken
fbuspf = zdata.mfrom[pfidx] #node from where power.act.measure starts
tbuspf = zdata.mto[pfidx] #node to where power.act.measure goes
fbusqf = zdata.mfrom[qfidx] #node from where power.react.measure starts
#tbusqf = zdata.mto[qfidx] #should be added [aan] ? #node to where power.react.measure goes
fbusiamp = zdata.mfrom[iidx] #node from where current.flow.measure starts
tbusiamp = zdata.mto[iidx] #node to where current.flow.measure goes
z = zdata.mval #values of measurements are stored in array z
#extrapolate different types of measurements
Pinj = z[pidx]
Qinj = z[qidx]
Pbr = z[pfidx]
Qbr = z[qfidx]
# zero injection measurements
# the weight is small and can bring instability during matrix inversion, so we "cut" everything below 10^-6
idx = numpy.where(zdata.mstddev<10**(-6))
zdata.mstddev[idx] = 10**(-6)
# weights matrix is obtained as stdandard_deviations^-2
weights = zdata.mstddev**(-2)
W = numpy.diag(weights)
# Admittances of the lines of the network
Gmatrix = system.Ymatrix.real
Bmatrix = system.Ymatrix.imag
Yabs_matrix = numpy.absolute(system.Ymatrix)
Yphase_matrix = numpy.angle(system.Ymatrix)
def DsseCall(system, measurements):
"""
Performs state estimation
It identifies the type of measurements present in the measurement set and
calls the appropriate estimator for dealing with them.
params:
@system: model of the system (nodes, lines, topology)
@measurements: Vector of measurements in Input (voltages, currents, powers)
return: object of class results.Results
"""
# 1. Select type of Estimator.
# If at least a PMU is present we launch the combined estimator, otherwise a simple traditional etimator
Vmag_meas = 0
Vpmu_meas = 0
for elem in measurements.measurements:
if elem.meas_type=="V_mag":
Vmag_meas = Vmag_meas + 1
elif elem.meas_type=="Vpmu_mag":
Vpmu_meas = Vpmu_meas + 1
trad_code = 1 if Vmag_meas >0 else 0
PMU_code = 2 if Vpmu_meas >0 else 0
est_code = trad_code + PMU_code
# 2. Run Estimator.
#number of nodes of the grid
nodes_num = len(system.nodes)
Gmatrix = system.Ymatrix.real
Bmatrix = system.Ymatrix.imag
Yabs_matrix = np.absolute(system.Ymatrix)
Yphase_matrix = np.angle(system.Ymatrix)
Adj = system.Adjacencies
if est_code == 1:
Vest = DsseTrad(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_matrix)
elif est_code == 2:
Vest = DssePmu(nodes_num, measurements, Gmatrix, Bmatrix)
else:
Vest = DsseMixed(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_matrix, Adj)
# 3. Populate arrays with data in output of estimator
# The Estimator provides:
# Vest: Node Voltage phasor estimated
# Iest: Branch Current phasor estimated
# Iinjest: Injection Current phasor estimated
# S1est: Complex Power flow at branch, measured at initial node
# S2est: Complex Power flow at branch, measured at final node
# Sinjest: Complex Power Injection at node
""" From here on, all the other quantities of the grid are calculated """
results = Results(system)
results.load_voltages(Vest)
results.calculateI()
results.calculateIinj()
results.calculateSinj()
results.calculateI()
results.calculateS1()
results.calculateS2()
return results
def DsseTrad(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_matrix):
"""
Traditional state estimator
It performs state estimation using rectangular node voltage state variables
and it is customized to work without PMU measurements
params:
@nodes_num: number of nodes of the grid
@measurements: Vector of measurements in Input (voltages, currents, powers)
@Gmatrix
@Bmatrix
@Yabs_matrix
@Yphase_matrix
return: np.array V - estimated voltages
"""
#calculate weights matrix (obtained as stdandard_deviations^-2)
weights = measurements.getWeightsMatrix()
W = np.diag(weights)
# Jacobian Matrix. Includes the derivatives of the measurements (voltages, currents, powers) with respect to the states (voltages)
#Jacobian for Power Injection Measurements
#H2, H3 = calculateJacobiMatrixSinj(measurements, nodes_num-1, Gmatrix, Bmatrix, Adj)
""" Jacobian for Power Injection Measurements (converted to equivalent
rectangualar current measurements) """
# Derivative of power injection(converted to current injection) with respect to voltage: it is the admittance.
H2 = numpy.zeros((npi,2*nodes_num-1))
H3 = numpy.zeros((nqi,2*nodes_num-1))
H2 = np.zeros((npi,2*nodes_num-1))
H3 = np.zeros((nqi,2*nodes_num-1))
for i in range(npi):
m = buspi[i]-1
m2 = m + nodes_num - 1
......@@ -126,307 +107,197 @@ def DsseTrad(system, zdata):
H2[i][m2] = Bmatrix[m][m]
H3[i][m] = - Bmatrix[m][m]
H3[i][m2] = - Gmatrix[m][m]
idx = numpy.subtract(system.Adj[m],1)
idx = np.subtract(system.Adjacencies[m],1)
H2[i][idx] = - Gmatrix[m][idx]
H3[i][idx] = - Bmatrix[m][idx]
if 0 in idx:
pos = numpy.where(idx==0)
idx = numpy.delete(idx,pos)
pos = np.where(idx==0)
idx = np.delete(idx,pos)
idx2 = idx + nodes_num-1
H2[i][idx2] = Bmatrix[m][idx]
H3[i][idx2] = - Gmatrix[m][idx]
""" Jacobian for branch Power Measurements (converted to equivalent
rectangualar current measurements)"""
# Derivative of branch power flow(converted to current flow) with respect to voltage: it is the admittance.
H4 = numpy.zeros((npf,2*nodes_num-1))
H5 = numpy.zeros((nqf,2*nodes_num-1))
for i in range(npf):
m = fbuspf[i]-1
n = tbuspf[i]-1
H4[i][m] = - Gmatrix[m][n]
H4[i][n] = Gmatrix[m][n]
H5[i][m] = - Bmatrix[m][n]
H5[i][n] = Bmatrix[m][n]
if m > 0:
m2 = m + nodes_num-1
H4[i][m2] = Bmatrix[m][n]
H5[i][m2] = - Gmatrix[m][n]
if n > 0:
n2 = n + nodes_num-1
H4[i][n2] = - Bmatrix[m][n]
H5[i][n2] = Gmatrix[m][n]
epsilon = 5 # treshold to stop Netwon Rapson iterations
Vr = numpy.ones(nodes_num) #initialize voltage real part to 1 per unit
Vx = numpy.zeros(nodes_num) #initialize voltage imaginary part to 0 per unit
V = Vr + 1j*Vx
num_iter = 0 # number of iterations of Newton Rapson
StateVr = numpy.ones(nodes_num) #initialize voltage real part to 1 per unit
StateVx = numpy.zeros(nodes_num-1) #initialize voltage imaginary part to 0 per unit
State = numpy.concatenate((StateVr,StateVx),axis=0)
#Jacobian for branch Power Measurements
H4, H5 = calculateJacobiBranchPower(measurements, nodes_num, Gmatrix, Bmatrix, type=2)
#get array which contains the index of measurements type MeasType.Sinj_real, MeasType.Sinj_imag in the array measurements.measurements
pidx = measurements.getIndexOfMeasurements(type=MeasType.Sinj_real)
qidx = measurements.getIndexOfMeasurements(type=MeasType.Sinj_imag)
#get array which contains the index of measurements type MeasType.S_real, MeasType.S_imag in the array measurements.measurements
p1br = measurements.getIndexOfMeasurements(type=MeasType.S1_real)
p2br = measurements.getIndexOfMeasurements(type=MeasType.S2_real)
q1br = measurements.getIndexOfMeasurements(type=MeasType.S1_imag)
q2br = measurements.getIndexOfMeasurements(type=MeasType.S2_imag)
#get an array with all measured values (affected by uncertainty)
z = measurements.getmVal()
V = np.ones(nodes_num) + 1j*np.zeros(nodes_num)
State = np.concatenate((np.ones(nodes_num), np.zeros(nodes_num)), axis=0)
epsilon = 5
num_iter = 0