Commit e5a45c3e authored by admin-mpa's avatar admin-mpa
Browse files

Added DsseAllocation method

parent d0745e12
......@@ -177,7 +177,7 @@ def DssePmu(nodes_num, measurements, Gmatrix, Bmatrix):
# Jacobian for branch Power Measurements
H4, H5 = calculateJacobiBranchPower(measurements, nodes_num, Gmatrix, Bmatrix, type=1)
# Jacobian for branch Power Measurements
# Jacobian for Voltage Pmu Measurements
H7, H8 = calculateJacobiVoltagePmu(measurements, nodes_num, Gmatrix, Bmatrix)
W = update_W_matrix(measurements, weights, W, "Vpmu")
......@@ -264,7 +264,7 @@ def DsseMixed(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_mat
# Jacobian for branch Power Measurements
H4, H5 = calculateJacobiBranchPower(measurements, nodes_num, Gmatrix, Bmatrix, type=1)
# Jacobian for branch Power Measurements
# Jacobian for Voltage Pmu Measurements
H7, H8 = calculateJacobiVoltagePmu(measurements, nodes_num, Gmatrix, Bmatrix)
W = update_W_matrix(measurements, weights, W, "Vpmu")
......@@ -346,7 +346,122 @@ def DsseMixed(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_mat
return V
def calculateJacobiMatrixSinj(measurements, nodes_num, Gmatrix, Bmatrix, type):
def DsseAllocation(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_matrix, meas_code, inj_code):
"""
Traditional state estimator
It performs state estimation using rectangular node voltage state variables
and it is built to work in scenarios where both conventional and PMU measurements
are simultaneously present.
@param nodes_num: number of nodes of the grid
@param measurements: Vector of measurements in Input (voltages, currents, powers)
@param Gmatrix: real part of the admittance matrix
@param Bmatrix: imaginary part of the admittance matrix
@param Yabs_matrix: module of the admittance matrix
@param Yphase_matrix: angles of the admittance matrix
@param meas_code: code to determine if PMU measurements are present or not (1=no PMUs, 2 or 3=PMUs present)
@param inj_code: code to determine if the grid has only loads (1) or both loads and generation units (2) connected
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)
if meas_code == 1:
type=1
elif meas_code == 2 or meas_code == 3:
type=2
# Jacobian for branch Power Measurements
H4, H5 = calculateJacobiBranchPower(measurements, nodes_num, Gmatrix, Bmatrix, type, inj_code)
# Jacobian for Voltage Pmu Measurements
H7, H8 = calculateJacobiVoltagePmu(measurements, nodes_num, Gmatrix, Bmatrix, inj_code)
W = update_W_matrix(measurements, weights, W, "Vpmu")
# Jacobian for Current Pmu Measurements
H9, H10 = calculateJacobiCurrentPmu(measurements, nodes_num, Gmatrix, Bmatrix, inj_code)
W = update_W_matrix(measurements, weights, W, "Ipmu")
# get array which contains the index of measurements type V_mag and I_mag
vidx = measurements.getIndexOfMeasurements(type=MeasType.V_mag)
iidx = measurements.getIndexOfMeasurements(type=MeasType.I_mag)
nvi = len(vidx)
nii = len(iidx)
# 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.getMeasValues()
V = np.ones(nodes_num) + 1j * np.zeros(nodes_num)
Kfactor = np.ones(inj_code)
State = np.concatenate((np.ones(nodes_num), np.zeros(nodes_num), Kfactor), axis=0)
epsilon = 5
num_iter = 0
# Iteration of Netwon Rapson method: needed to solve non-linear system of equation
while epsilon > 10 ** (-6):
""" Computation of equivalent current measurements in place of the power measurements """
z = convertSinjMeasIntoCurrents(measurements, V, z, pidx, qidx)
z = convertSbranchMeasIntoCurrents(measurements, V, z, p1br, q1br, p2br, q2br)
""" Voltage Magnitude Measurements """
h1, H1 = update_h1_vector(measurements, V, vidx, nvi, nodes_num, type, inj_code)
""" Power Injection Measurements """
h2, h3, H2, H3 = update_h2_h3_vector(measurements, nodes_num, V, Gmatrix, Bmatrix, type, inj_code, Kfactor)
""" Power Flow Measurements """
h4 = np.inner(H4, State)
h5 = np.inner(H5, State)
""" Current Magnitude Measurements """
h6, H6 = update_h6_vector(measurements, V, iidx, nii, Yabs_matrix, Yphase_matrix, nodes_num, num_iter, type, inj_code)
""" PMU Voltage Measurements """
h7 = np.inner(H7, State)
h8 = np.inner(H8, State)
""" PMU Current Measurements """
h9 = np.inner(H9, State)
h10 = np.inner(H10, State)
""" WLS computation """
H = np.concatenate((H1, H2, H3, H4, H5, H6, H7, H8, H9, H10), axis=0)
y = np.concatenate((h1, h2, h3, h4, h5, h6, h7, h8, h9, h10), axis=0)
res = np.subtract(z, y)
g = np.inner(H.transpose(), np.inner(W, res))
WH = np.inner(W, H.transpose())
G = np.inner(H.transpose(), WH.transpose())
Ginv = np.linalg.inv(G)
Delta_State = np.inner(Ginv, g)
State = State + Delta_State
epsilon = np.amax(np.absolute(Delta_State))
V.real = State[:nodes_num]
V.imag = State[nodes_num:-inj_code]
Kfactor = State[-inj_code:]
num_iter = num_iter + 1
return V
def calculateJacobiMatrixSinj(measurements, nodes_num, Gmatrix, Bmatrix, type, add_var):
"""
It calculates the Jacobian for Power Injection Measurements
(converted to equivalent rectangualar current measurements)
......@@ -355,7 +470,8 @@ def calculateJacobiMatrixSinj(measurements, nodes_num, Gmatrix, Bmatrix, type):
@param nodes_num: len of system.nodes
@param Gmatrix
@param Bmatrix
@param type: 1 for DssePmu and DssePmu, 2 for DsseTrad
@param type: 1 for DssePmu and DsseMixed, 2 for DsseTrad
@param add_var: additional number of variables to be considered when running the DsseAllocation method
return: 1. H2: Jacobian for Pinj
2. H3: Jacobian for Qinj
"""
......@@ -364,11 +480,11 @@ def calculateJacobiMatrixSinj(measurements, nodes_num, Gmatrix, Bmatrix, type):
# get all measurements of type MeasType.Sinj_real
qinj_meas = measurements.getMeasurements(type=MeasType.Sinj_imag)
if type == 1:
H2 = np.zeros((len(pinj_meas), 2 * nodes_num))
H3 = np.zeros((len(qinj_meas), 2 * nodes_num))
H2 = np.zeros((len(pinj_meas), 2 * nodes_num + add_var))
H3 = np.zeros((len(qinj_meas), 2 * nodes_num + add_var))
elif type == 2:
H2 = np.zeros((len(pinj_meas), 2 * nodes_num - 1))
H3 = np.zeros((len(qinj_meas), 2 * nodes_num - 1))
H2 = np.zeros((len(pinj_meas), 2 * nodes_num - 1 + add_var))
H3 = np.zeros((len(qinj_meas), 2 * nodes_num - 1 + add_var))
for index, measurement in enumerate(pinj_meas):
m = measurement.element.index
......@@ -383,7 +499,7 @@ def calculateJacobiMatrixSinj(measurements, nodes_num, Gmatrix, Bmatrix, type):
return H2, H3
def calculateJacobiBranchPower(measurements, nodes_num, Gmatrix, Bmatrix, type):
def calculateJacobiBranchPower(measurements, nodes_num, Gmatrix, Bmatrix, type, add_var):
"""
It calculates the Jacobian for branch Power Measurements
(converted to equivalent rectangualar current measurements)
......@@ -392,7 +508,8 @@ def calculateJacobiBranchPower(measurements, nodes_num, Gmatrix, Bmatrix, type):
@param nodes_num: len of system.nodes
@param Gmatrix
@param Bmatrix
@param type: 1 for DssePmu and DssePmu, 2 for DsseTrad
@param type: 1 for DssePmu and DsseMixed, 2 for DsseTrad
@param add_var: additional number of variables to be considered when running the DsseAllocation method
return: 1. H4: Jacobian for S_real
2. H5: Jacobian for S_imag
"""
......@@ -404,11 +521,11 @@ def calculateJacobiBranchPower(measurements, nodes_num, Gmatrix, Bmatrix, type):
q2_meas = measurements.getMeasurements(type=MeasType.S2_imag)
if type == 1:
H4 = np.zeros((len(p1_meas) + len(p2_meas), 2 * nodes_num))
H5 = np.zeros((len(q1_meas) + len(q2_meas), 2 * nodes_num))
H4 = np.zeros((len(p1_meas) + len(p2_meas), 2 * nodes_num + add_var))
H5 = np.zeros((len(q1_meas) + len(q2_meas), 2 * nodes_num + add_var))
elif type == 2:
H4 = np.zeros((len(p1_meas) + len(p2_meas), 2 * nodes_num - 1))
H5 = np.zeros((len(q1_meas) + len(q2_meas), 2 * nodes_num - 1))
H4 = np.zeros((len(p1_meas) + len(p2_meas), 2 * nodes_num - 1 + add_var))
H5 = np.zeros((len(q1_meas) + len(q2_meas), 2 * nodes_num - 1 + add_var))
for i, measurement in enumerate(p1_meas):
m = measurement.element.start_node.index
......@@ -461,7 +578,7 @@ def calculateJacobiBranchPower(measurements, nodes_num, Gmatrix, Bmatrix, type):
return H4, H5
def calculateJacobiVoltagePmu(measurements, nodes_num, Gmatrix, Bmatrix):
def calculateJacobiVoltagePmu(measurements, nodes_num, Gmatrix, Bmatrix, add_var):
"""
It calculates the Jacobian for Voltage Pmu Measurements
(converted to equivalent rectangualar current measurements)
......@@ -470,6 +587,7 @@ def calculateJacobiVoltagePmu(measurements, nodes_num, Gmatrix, Bmatrix):
@param nodes_num: len of system.nodes
@param Gmatrix
@param Bmatrix
@param add_var: additional number of variables to be considered when running the DsseAllocation method
return: 1. H7: Jacobian for S_real
2. H8: Jacobian for S_imag
"""
......@@ -478,8 +596,8 @@ def calculateJacobiVoltagePmu(measurements, nodes_num, Gmatrix, Bmatrix):
Vpmu_mag_meas = measurements.getMeasurements(type=MeasType.Vpmu_mag)
# get all measurements of type MeasType.Vpmu_phase
Vpmu_phase_meas = measurements.getMeasurements(type=MeasType.Vpmu_phase)
H7 = np.zeros((len(Vpmu_mag_meas), 2 * nodes_num))
H8 = np.zeros((len(Vpmu_mag_meas), 2 * nodes_num))
H7 = np.zeros((len(Vpmu_mag_meas), 2 * nodes_num + add_var))
H8 = np.zeros((len(Vpmu_mag_meas), 2 * nodes_num + add_var))
# TODO: index of Vmag = index of Vphase???
for index, measurement in enumerate(Vpmu_mag_meas):
......@@ -493,7 +611,7 @@ def calculateJacobiVoltagePmu(measurements, nodes_num, Gmatrix, Bmatrix):
return H7, H8
def calculateJacobiCurrentPmu(measurements, nodes_num, Gmatrix, Bmatrix):
def calculateJacobiCurrentPmu(measurements, nodes_num, Gmatrix, Bmatrix, add_var):
"""
It calculates the Jacobian for Current Pmu Measurements
(converted to equivalent rectangualar current measurements)
......@@ -502,6 +620,7 @@ def calculateJacobiCurrentPmu(measurements, nodes_num, Gmatrix, Bmatrix):
@param nodes_num: len of system.nodes
@param Gmatrix
@param Bmatrix
@param add_var: additional number of variables to be considered when running the DsseAllocation method
return: 1. H9: Jacobian for S_real
2. H10: Jacobian for S_imag
"""
......@@ -510,8 +629,8 @@ def calculateJacobiCurrentPmu(measurements, nodes_num, Gmatrix, Bmatrix):
Ipmu_mag_meas = measurements.getMeasurements(type=MeasType.Ipmu_mag)
# get all measurements of type MeasType.Vpmu_phase
Ipmu_phase_meas = measurements.getMeasurements(type=MeasType.Vpmu_phase)
H9 = np.zeros((len(Ipmu_mag_meas), 2 * nodes_num))
H10 = np.zeros((len(Ipmu_mag_meas), 2 * nodes_num))
H9 = np.zeros((len(Ipmu_mag_meas), 2 * nodes_num + add_var))
H10 = np.zeros((len(Ipmu_mag_meas), 2 * nodes_num + add_var))
for index, measurement in enumerate(Ipmu_mag_meas):
iamp = measurement.meas_value_ideal
......@@ -569,7 +688,7 @@ def update_W_matrix(measurements, weights, W, type):
return W
def update_h1_vector(measurements, V, vidx, nvi, nodes_num, type):
def update_h1_vector(measurements, V, vidx, nvi, nodes_num, type, add_var):
"""
update h1 and H1 vectors
......@@ -578,7 +697,8 @@ def update_h1_vector(measurements, V, vidx, nvi, nodes_num, type):
@param vidx: array which contains the index of measurements type V_mag in measurements.measurements
@param nvi: len of vidx
@param nodes_num: number of nodes of the grid - len(system.nodes)
@param type: 1 for DssePmu and DssePmu, 2 for DsseTrad
@param type: 1 for DssePmu and DsseMixed, 2 for DsseTrad
@param add_var: additional number of variables to be considered when running the DsseAllocation method
return: vector h1 and H1
"""
......@@ -586,9 +706,9 @@ def update_h1_vector(measurements, V, vidx, nvi, nodes_num, type):
h1 = np.zeros(nvi)
# the Jacobian rows where voltage measurements are presents is updated
if type == 1:
H1 = np.zeros((nvi, 2 * nodes_num))
H1 = np.zeros((nvi, 2 * nodes_num + add_var))
elif type == 2:
H1 = np.zeros((nvi, 2 * nodes_num - 1))
H1 = np.zeros((nvi, 2 * nodes_num - 1 + add_var))
for i, index_vmag in enumerate(vidx):
# get index of the node
node_index = measurements.measurements[index_vmag].element.index
......@@ -598,14 +718,14 @@ def update_h1_vector(measurements, V, vidx, nvi, nodes_num, type):
m2 = node_index + nodes_num
H1[i][m2] = np.sin(np.angle(V[node_index]))
elif type == 2:
if m > 0:
if node_index > 0:
m2 = node_index + nodes_num - 1
H1[i][m2] = np.sin(V.phase[m])
H1[i][m2] = np.sin(V.phase[node_index])
return h1, H1
def update_h6_vector(measurements, V, iidx, nii, Yabs_matrix, Yphase_matrix, nodes_num, num_iter, type):
def update_h6_vector(measurements, V, iidx, nii, Yabs_matrix, Yphase_matrix, nodes_num, num_iter, type, add_var):
"""
update h6 and H6 vectors where current flows are present
......@@ -617,7 +737,8 @@ def update_h6_vector(measurements, V, iidx, nii, Yabs_matrix, Yphase_matrix, nod
@param Yphase_matrix:
@param nodes_num: number of nodes of the grid - len(system.nodes)
@param num_iter: number of current iteration
@param type: 1 for DssePmu and DssePmu, 2 for DsseTrad
@param type: 1 for DssePmu and DsseMixed, 2 for DsseTrad
@param add_var: additional number of variables to be considered when running the DsseAllocation method
return: vector h6 and H6
"""
......@@ -627,9 +748,9 @@ def update_h6_vector(measurements, V, iidx, nii, Yabs_matrix, Yphase_matrix, nod
h6complex = np.zeros((nii), dtype=complex)
h6 = np.ones((nii))
if type == 1:
H6 = np.zeros((nii, 2 * nodes_num))
H6 = np.zeros((nii, 2 * nodes_num + add_var))
elif type == 2:
H6 = np.zeros((nii, 2 * nodes_num - 1))
H6 = np.zeros((nii, 2 * nodes_num - 1 + add_var))
for i, index_imag in enumerate(iidx):
# get index of the start node
......@@ -668,6 +789,57 @@ def update_h6_vector(measurements, V, iidx, nii, Yabs_matrix, Yphase_matrix, nod
return h6, H6
def update_h2_h3_vector(measurements, nodes_num, V, Gmatrix, Bmatrix, type, add_var, Kfactor):
"""
It calculates the Jacobian for Power Injection Measurements
(converted to equivalent rectangualar current measurements)
@param measurements: object of class Measurement that contains all measurements (voltages, currents, powers)
@param nodes_num: len of system.nodes
@param V: vector of the estimated voltages
@param Gmatrix
@param Bmatrix
@param type: 1 for DssePmu and DsseMixed, 2 for DsseTrad
@param add_var: additional number of variables to be considered when running the DsseAllocation method
@param Kfactor: temporary value of the scaling factor given by the DsseAllocation method
return: 1. H2: Jacobian for Pinj
2. H3: Jacobian for Qinj
"""
# get all measurements of type MeasType.Sinj_real
pinj_meas = measurements.getMeasurements(type=MeasType.Sinj_real)
# get all measurements of type MeasType.Sinj_real
qinj_meas = measurements.getMeasurements(type=MeasType.Sinj_imag)
h2 = np.zeros((len(pinj_meas)))
h3 = np.zeros((len(qinj_meas)))
if type == 1:
H2 = np.zeros((len(pinj_meas), 2 * nodes_num + add_var))
H3 = np.zeros((len(qinj_meas), 2 * nodes_num + add_var))
elif type == 2:
H2 = np.zeros((len(pinj_meas), 2 * nodes_num - 1 + add_var))
H3 = np.zeros((len(qinj_meas), 2 * nodes_num - 1 + add_var))
for index, measurement in enumerate(pinj_meas):
m = measurement.element.index
if type == 1:
idx = 0
elif type == 2:
idx = 1
#if type(m) == 'load':
# K = Kfactor[0]
# idxK = 0
#elif type(m) == 'generation':
# K = Kfactor[1]
# idxK = 1
H2[index][:nodes_num] = K*Gmatrix[m]
H2[index][nodes_num:-add_var] = -K*Bmatrix[m][idx:]
H2[index][2*nodes_num-idx+idxK] = np.inner(Gmatrix[m],V.real) - np.inner(Bmatrix[m][idx:],V.imag)
H3[index][:nodes_num] = K*Bmatrix[m]
H3[index][nodes_num:-add_var] = K*Gmatrix[m][idx:]
H3[index][2*nodes_num-idx+idxK] = np.inner(Bmatrix[m],V.real) + np.inner(Gmatrix[m][idx:],V.imag)
h2[index] = K*(np.inner(Gmatrix[m],V.real) - np.inner(Bmatrix[m][idx:],V.imag))
h3[index] = K*(np.inner(Bmatrix[m],V.real) + np.inner(Gmatrix[m][idx:],V.imag))
return h2, h3, H2, H3
def convertSinjMeasIntoCurrents(measurements, V, z, pidx, qidx):
"""
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment