Commit ddb8c423 authored by Marco Pau's avatar Marco Pau
Browse files

Modifications for the DsseAllocation method

parent 2bd321eb
......@@ -40,14 +40,14 @@ def solve(system):
H[m][i] = 1
H[m + 1][i2] = 1
elif node_type is BusType.PQ:
H[m][:node.num] = np.real(system.Ymatrix[i])
H[m][node.num:] = - np.imag(system.Ymatrix[i])
H[m+1][:node.num] = np.imag(system.Ymatrix[i])
H[m+1][node.num:] = np.real(system.Ymatrix[i])
H[m][:nodes_num] = np.real(system.Ymatrix[i])
H[m][nodes_num:] = - np.imag(system.Ymatrix[i])
H[m+1][:nodes_num] = np.imag(system.Ymatrix[i])
H[m+1][nodes_num:] = np.real(system.Ymatrix[i])
elif node_type is BusType.PV:
z[m + 1] = np.real(node.power)
H[m][:node.num] = np.real(system.Ymatrix[i])
H[m][node.num:] = - np.imag(system.Ymatrix[i])
H[m][:nodes_num] = np.real(system.Ymatrix[i])
H[m][nodes_num:] = - np.imag(system.Ymatrix[i])
epsilon = 10 ** (-10)
#epsilon = 0.01
......
......@@ -3,7 +3,7 @@ from acs.state_estimation.results import Results
from acs.state_estimation.measurement import *
def DsseCall(system, measurements):
def DsseCall(system, measurements, solver_type="conventional"):
"""
Performs state estimation
It identifies the type of measurements present in the measurement set and
......@@ -37,13 +37,20 @@ def DsseCall(system, measurements):
Yabs_matrix = np.absolute(system.Ymatrix)
Yphase_matrix = np.angle(system.Ymatrix)
# run Estimator.
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)
if solver_type == "conventional":
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)
elif solver_type == "advanced":
# TODO: derive from system inj_code analyzing whether load and gens connected to all nodes
inj_code = 1
if inj_code == 1 or inj_code == 2:
Vest = DsseAllocation(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_matrix, est_code, inj_code)
# calculate all the other quantities of the grid
results = Results(system)
......@@ -410,7 +417,10 @@ def DsseAllocation(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphas
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)
if type == 1:
State = np.concatenate((np.ones(nodes_num), np.zeros(nodes_num-1), Kfactor), axis=0)
elif type == 2:
State = np.concatenate((np.ones(nodes_num), np.zeros(nodes_num), Kfactor), axis=0)
epsilon = 5
num_iter = 0
......@@ -443,6 +453,16 @@ def DsseAllocation(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphas
""" WLS computation """
H = np.concatenate((H1, H2, H3, H4, H5, H6, H7, H8, H9, H10), axis=0)
if num_iter == 0:
if meas_code == 1:
H = numpy.delete(H,2*node.num-1,1)
else:
H = numpy.delete(H,2*node.num,1)
if inj_code == 2:
if meas_code == 1:
H = numpy.delete(H,2*node.num-1,1)
else:
H = numpy.delete(H,2*node.num,1)
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))
......@@ -451,6 +471,8 @@ def DsseAllocation(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphas
Ginv = np.linalg.inv(G)
Delta_State = np.inner(Ginv, g)
if num_iter == 0:
Delta_State = numpy.concatenate((Delta_State,numpy.zeros((inj_code))),axis=0)
State = State + Delta_State
epsilon = np.amax(np.absolute(Delta_State))
......@@ -827,6 +849,9 @@ def update_h2_h3_vector(measurements, nodes_num, V, Gmatrix, Bmatrix, inj_code,
idx = 0
elif type == 2:
idx = 1
K = Kfactor[0]
idxK = 0
# TODO:
#if type(m) == 'load':
# K = Kfactor[0]
# idxK = 0
......@@ -835,12 +860,12 @@ def update_h2_h3_vector(measurements, nodes_num, V, Gmatrix, Bmatrix, inj_code,
# idxK = 1
H2[index][:nodes_num] = K*Gmatrix[m]
H2[index][nodes_num:-inj_code] = -K*Bmatrix[m][idx:]
H2[index][2*nodes_num-idx+idxK] = np.inner(Gmatrix[m],V.real) - np.inner(Bmatrix[m][idx:],V.imag)
H2[index][2*nodes_num-idx+idxK] = np.inner(Gmatrix[m],V.real) - np.inner(Bmatrix[m][idx:],V.imag[idx:])
H3[index][:nodes_num] = K*Bmatrix[m]
H3[index][nodes_num:-inj_code] = 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))
H3[index][2*nodes_num-idx+idxK] = np.inner(Bmatrix[m],V.real) + np.inner(Gmatrix[m][idx:],V.imag[idx:])
h2[index] = K*(np.inner(Gmatrix[m],V.real) - np.inner(Bmatrix[m][idx:],V.imag[idx:]))
h3[index] = K*(np.inner(Bmatrix[m],V.real) + np.inner(Gmatrix[m][idx:],V.imag[idx:]))
return h2, h3, H2, H3
......
......@@ -45,7 +45,6 @@ class Results():
self.branches = []
self.Ymatrix = system.Ymatrix
self.Bmatrix = system.Bmatrix
self.Adjacencies = system.Adjacencies
for node in system.nodes:
if node.ideal_connected_with == '':
self.nodes.append(ResultsNode(topo_node=node))
......
......@@ -26,7 +26,7 @@ for file in xml_files:
xml_files_abs.append(os.path.abspath(file))
# read cim files and create new network.Systen object
res, _ = cimpy.cim_import(xml_files_abs, "cgmes_v2_4_15")
res, _, _ = cimpy.cim_import(xml_files_abs, "cgmes_v2_4_15")
system = network.System()
base_apparent_power = 25 # MW
system.load_cim_data(res, base_apparent_power)
......@@ -40,7 +40,7 @@ V_unc = 0
I_unc = 0
Sinj_unc = 0
S_unc = 0
Pmu_mag_unc = 1
Pmu_mag_unc = 0
Pmu_phase_unc = 0
# Create measurements data structures
......
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