Commit 13643793 authored by Martin Alonso Moraga's avatar Martin Alonso Moraga
Browse files

new example for state estimation with network.py

parent 561b21ec
__all__ = ["bc_powerflow", "bc_state_estimator", "nv_powerflow", "nv_state_estimator", "network"]
\ No newline at end of file
__all__ = ["bc_powerflow", "bc_state_estimator", "nv_powerflow", "nv_state_estimator", "network", "nv_powerflow_cim"]
\ No newline at end of file
......@@ -77,4 +77,20 @@ def load_python_data(nodes, branches, type):
system.branches.append(Branch(branches.R[branch_idx], branches.X[branch_idx],
system.nodes[branches.start[branch_idx]-1], system.nodes[branches.end[branch_idx]-1]))
return system
\ No newline at end of file
return system
def Ymatrix_calc(system):
nodes_num = len(system.nodes)
branches_num = len(system.branches)
Ymatrix = np.zeros((nodes_num, nodes_num),dtype=np.complex)
Adjacencies = [[] for _ in range(nodes_num)]
for index in range(branches_num):
fr = system.branches[index].start_node.index
to = system.branches[index].end_node.index
Ymatrix[fr][to] -= system.branches[index].y
Ymatrix[to][fr] -= system.branches[index].y
Ymatrix[fr][fr] += system.branches[index].y
Ymatrix[to][to] += system.branches[index].y
Adjacencies[fr].append(to+1)
Adjacencies[to].append(fr+1)
return Ymatrix, Adjacencies
\ No newline at end of file
import numpy as np
import math
from network import BusType
def Ymatrix_calc(system):
nodes_num = len(system.nodes)
branches_num = len(system.branches)
......@@ -17,7 +17,7 @@ def Ymatrix_calc(system):
Adjacencies[fr].append(to+1)
Adjacencies[to].append(fr+1)
return Ymatrix, Adjacencies
def solve(system):
"""It performs Power Flow by using rectangular node voltage state variables."""
......@@ -132,5 +132,4 @@ def solve(system):
S1=np.append(S1, V[system.branches[i].start_node.index]*(I_conj[i]))
S2=np.append(S2, -V[system.branches[i].end_node.index]*(I_conj[i]))
return V, I, Iinj, S1, S2, Sinj, num_iter
\ No newline at end of file
return V, I, Iinj, S1, S2, Sinj, num_iter
\ No newline at end of file
......@@ -465,7 +465,7 @@ def DsseMixed(branch, node, zdata, Ymatrix, Adj):
tbusipmu = zdata.mto[imagpmuidx]
z = zdata.mval
#print(z)
Pinj = z[pidx]
Qinj = z[qidx]
Pbr = z[pfidx]
......@@ -499,7 +499,7 @@ def DsseMixed(branch, node, zdata, Ymatrix, Adj):
idx2 = idx + node.num
H2[i][idx2] = Bmatrix[m][idx]
H3[i][idx2] = - Gmatrix[m][idx]
""" Jacobian for branch Power Measurements (converted to equivalent
rectangualar current measurements)"""
H4 = numpy.zeros((npf,2*node.num))
......@@ -517,7 +517,7 @@ def DsseMixed(branch, node, zdata, Ymatrix, Adj):
n2 = n + node.num
H4[i][n2] = - Bmatrix[m][n]
H5[i][n2] = Gmatrix[m][n]
""" Jacobian for Voltage Pmu Measurements (converted into rectangular) """
H7 = numpy.zeros((nvpmu,2*node.num))
H8 = numpy.zeros((nvpmu,2*node.num))
......@@ -569,7 +569,7 @@ def DsseMixed(branch, node, zdata, Ymatrix, Adj):
H9[i][n2] = - Bmatrix[m][n]
H10[i][m2] = - Gmatrix[m][n]
H10[i][n2] = Gmatrix[m][n]
epsilon = 5
Vr = numpy.ones(node.num)
Vx = numpy.zeros(node.num)
......
This diff is collapsed.
def Zdatatrue_creation(zdata, zdatameas, meas, system, V, I, Iinj, S1, S2, Sinj):
""" It gives the true values at the measurement points."""
import numpy as np
t1 = np.ones(meas.V.num)
t2 = 2*np.ones(meas.Sinj.num)
t3 = 3*np.ones(meas.Sinj.num)
t4 = 4*np.ones(meas.S1.num + meas.S2.num)
t5 = 5*np.ones(meas.S1.num + meas.S2.num)
t6 = 6*np.ones(meas.I.num)
t7 = 7*np.ones(meas.Vpmu_mag.num)
t8 = 8*np.ones(meas.Vpmu_phase.num)
t9 = 9*np.ones(meas.Ipmu_mag.num)
t10 = 10*np.ones(meas.Ipmu_phase.num)
zdata.mtype = np.concatenate((t1,t2,t3,t4,t5,t6,t7,t8,t9,t10),axis=0)
zdatameas.mtype = np.concatenate((t1,t2,t3,t4,t5,t6,t7,t8,t9,t10),axis=0)
z1 = np.abs(V[meas.V.index-1])
z2 = Sinj.real[meas.Sinj.index-1]
z3 = Sinj.imag[meas.Sinj.index-1]
z4_1 = S1.real[meas.S1.index-1]
z4_2 = S2.real[meas.S2.index-1]
z5_1 = S1.imag[meas.S1.index-1]
z5_2 = S2.imag[meas.S2.index-1]
z6 = np.abs(I[meas.I.index-1])
z7 = np.abs(V[meas.Vpmu_mag.index-1])
z8 = np.angle(V[meas.Vpmu_phase.index-1])
z9 = np.abs(I[meas.Ipmu_mag.index-1])
z10 = np.angle(I[meas.Ipmu_phase.index-1])
zdata.mval = np.concatenate((z1,z2,z3,z4_1,z4_2,z5_1,z5_2,z6,z7,z8,z9,z10),axis=0)
"""
b1 = np.zeros(meas.V.num)
b2 = np.zeros(meas.Sinj.num)
b3 = np.zeros(meas.Sinj.num)
b4_1 = branch.code[meas.S1.index-1]
b4_2 = -branch.code[meas.S2.index-1]
b5_1 = branch.code[meas.S1.index-1]
b5_2 = -branch.code[meas.S2.index-1]
b6 = branch.code[meas.I.index-1]
b7 = np.zeros(meas.Vpmu_mag.num)
b8 = np.zeros(meas.Vpmu_phase.num)
b9 = branch.code[meas.Ipmu_mag.index-1]
b10 = branch.code[meas.Ipmu_phase.index-1]
zdata.mbranch = np.concatenate((b1,b2,b3,b4_1,b4_2,b5_1,b5_2,b6,b7,b8,b9,b10),axis=0)
zdata.mbranch = zdata.mbranch.astype(int)
zdatameas.mbranch = np.concatenate((b1,b2,b3,b4_1,b4_2,b5_1,b5_2,b6,b7,b8,b9,b10),axis=0)
zdatameas.mbranch = zdata.mbranch.astype(int)
"""
fr1 = meas.V.index
fr2 = meas.Sinj.index
fr3 = meas.Sinj.index
fr4_1 = np.array([])
fr4_2 = np.array([])
fr5_1 = np.array([])
fr5_2 = np.array([])
for index in range(len(meas.S1.index)):
fr4_1 = np.append(fr4_1, system.branches[meas.S1.index[index]-1].start_node.index)
fr5_1 = np.append(fr5_1, system.branches[meas.S1.index[index]-1].start_node.index)
for index in range(len(meas.S2.index)):
fr4_2 = np.append(fr4_2, system.branches[meas.S2.index[index]-1].end_node.index)
fr5_2 = np.append(fr5_2, system.branches[meas.S2.index[index]-1].end_node.index)
fr6 = np.array([])
for index in range(len(meas.I.index)):
fr6 = np.append(fr6, system.branches[meas.I.index[index]-1].start_node.index)
fr7 = meas.Vpmu_mag.index
fr8 = meas.Vpmu_phase.index
fr9 = np.array([])
fr10 = np.array([])
for index in range(len(meas.Ipmu_mag.index)):
fr9 = np.append(fr9, system.branches[meas.Ipmu_mag.index[index]-1].start_node.index)
for index in range(len(meas.Ipmu_phase.index)):
fr10 = np.append(fr10, system.branches[meas.Ipmu_phase.index[index]-1].start_node.index)
zdata.mfrom = np.concatenate((fr1,fr2,fr3,fr4_1,fr4_2,fr5_1,fr5_2,fr6,fr7,fr8,fr9,fr10),axis=0)
zdata.mfrom = zdata.mfrom.astype(int)
zdatameas.mfrom = np.concatenate((fr1,fr2,fr3,fr4_1,fr4_2,fr5_1,fr5_2,fr6,fr7,fr8,fr9,fr10),axis=0)
zdatameas.mfrom = zdata.mfrom.astype(int)
to1 = np.zeros(meas.V.num)
to2 = np.zeros(meas.Sinj.num)
to3 = np.zeros(meas.Sinj.num)
to4_1 = np.array([])
to4_2 = np.array([])
to5_1 = np.array([])
to5_2 = np.array([])
for index in range(len(meas.S1.index)):
to4_1 = np.append(to4_1, system.branches[meas.S1.index[index]-1].end_node.index)
to5_1 = np.append(to5_1, system.branches[meas.S1.index[index]-1].end_node.index)
for index in range(len(meas.S2.index)):
to4_2 = np.append(to4_2, system.branches[meas.S2.index[index]-1].start_node.index)
to5_2 = np.append(to5_2, system.branches[meas.S2.index[index]-1].start_node.index)
to6 = np.array([])
for index in range(len(meas.I.index)):
to6 = np.append(to6, system.branches[meas.I.index[index]-1].end_node.index)
to7 = np.zeros(meas.Vpmu_mag.num)
to8 = np.zeros(meas.Vpmu_phase.num)
to9 = np.array([])
to10 = np.array([])
for index in range(len(meas.Ipmu_mag.index)):
to9 = np.append(to9, system.branches[meas.Ipmu_mag.index[index]-1].end_node.index)
for index in range(len(meas.Ipmu_phase.index)):
to10 = np.append(to10, system.branches[meas.Ipmu_phase.index[index]-1].end_node.index)
zdata.mto = np.concatenate((to1,to2,to3,to4_1,to4_2,to5_1,to5_2,to6,to7,to8,to9,to10),axis=0)
zdata.mto = zdata.mto.astype(int)
zdatameas.mto = np.concatenate((to1,to2,to3,to4_1,to4_2,to5_1,to5_2,to6,to7,to8,to9,to10),axis=0)
zdatameas.mto = zdata.mto.astype(int)
d1 = z1*(meas.V.unc/300)
d2 = np.absolute(z2*(meas.Sinj.unc/300))
d3 = np.absolute(z3*(meas.Sinj.unc/300))
d4_1 = np.absolute(z4_1*(meas.S1.unc/300))
d4_2 = np.absolute(z4_2*(meas.S2.unc/300))
d5_1 = np.absolute(z5_1*(meas.S1.unc/300))
d5_2 = np.absolute(z5_2*(meas.S2.unc/300))
d6 = z6*(meas.I.unc/300)
d7 = z7*(meas.Vpmu_mag.unc/300)
d8 = (meas.Vpmu_phase.unc/300)*np.ones(meas.Vpmu_phase.num)
d9 = z9*(meas.Ipmu_mag.unc/300)
d10 = (meas.Ipmu_phase.unc/300)*np.ones(meas.Ipmu_phase.num)
zdata.mstddev = np.concatenate((d1,d2,d3,d4_1,d4_2,d5_1,d5_2,d6,d7,d8,d9,d10),axis=0)
zdatameas.mstddev = np.concatenate((d1,d2,d3,d4_1,d4_2,d5_1,d5_2,d6,d7,d8,d9,d10),axis=0)
# idx = np.where(zdata.mstddev<10**(-6))
# zdata.mstddev[idx] = 10**(-6)
return zdata, zdatameas
def Zdatameas_creation(zdata, zdatameas):
""" It gives the measured values (affected by uncertainty) at the measurement points."""
import numpy as np
zmeas = zdata.mval
zdev = zdata.mstddev
err_pu = np.random.normal(0,1,len(zmeas))
zdatameas.mval = zmeas + np.multiply(zdev,err_pu)
return zdatameas
\ No newline at end of file
......@@ -60,8 +60,8 @@ def Zdatatrue_creation(zdata, zdatameas, meas, branch, V, I, Iinj, S1, S2, Sinj)
fr7 = meas.Vpmu_mag.index
fr8 = meas.Vpmu_phase.index
fr9 = branch.start[meas.Ipmu_mag.index-1]
fr10 = branch.start[meas.Ipmu_phase.index-1]
fr10 = branch.start[meas.Ipmu_phase.index-1]
zdata.mfrom = np.concatenate((fr1,fr2,fr3,fr4_1,fr4_2,fr5_1,fr5_2,fr6,fr7,fr8,fr9,fr10),axis=0)
zdata.mfrom = zdata.mfrom.astype(int)
zdatameas.mfrom = np.concatenate((fr1,fr2,fr3,fr4_1,fr4_2,fr5_1,fr5_2,fr6,fr7,fr8,fr9,fr10),axis=0)
......
......@@ -12,11 +12,11 @@ from nv_powerflow_cim import solve as solve_cim
from nv_powerflow import *
class PerUnit:
def __init__(self, S, V):
self.S = S
self.V = V
self.I = S/V
self.Z = S/(V**2)
def __init__(self, S, V):
self.S = S
self.V = V
self.I = S/V
self.Z = S/(V**2)
""" Insert here per unit values of the grid for power and voltage """
......
import numpy as np
import math
import sys
import copy
sys.path.append("../acs/state_estimation")
import py_95bus_network_data
import py_95bus_meas_data
import nv_state_estimator
import nv_powerflow
import network
import nv_powerflow_cim
import cim_py_95bus_meas_data
import nv_state_estimator_cim
class PerUnit:
def __init__(self, S, V):
self.S = S
self.V = V
self.I = S/V
self.Z = S/(V**2)
class Measurements:
def __init__(self, index, unc):
self.index = index.astype(int)
self.unc = unc
self.num = len(index)
class Measurement_set:
def __init__(self, V, I, Sinj, S1, S2, Vpmu_mag, Vpmu_phase, Ipmu_mag, Ipmu_phase):
self.V = V
self.Sinj = Sinj
self.S1 = S1
self.S2 = S2
self.I = I
self.Vpmu_mag = Vpmu_mag
self.Vpmu_phase = Vpmu_phase
self.Ipmu_mag = Ipmu_mag
self.Ipmu_phase = Ipmu_phase
class Zdata_init:
def __init__(self, meas):
nmeas = meas.V.num + meas.I.num + 2*meas.Sinj.num + 2*meas.S1.num + 2*meas.S2.num + meas.Ipmu_mag.num + meas.Ipmu_phase.num + meas.Vpmu_mag.num + meas.Vpmu_phase.num
self.mtype = np.zeros(nmeas)
self.mval = np.zeros(nmeas)
self.mbranch = np.zeros(nmeas)
self.mfrom = np.zeros(nmeas)
self.mto = np.zeros(nmeas)
self.mstddev = np.zeros(nmeas)
""" Insert here per unit values of the grid for power and voltage """
S = 100*(10**6)
V = (11*(10**3))/math.sqrt(3)
slackV = 1.02
Base = PerUnit(S,V)
branch, node = py_95bus_network_data.Network_95_nodes(Base, slackV)
system = network.load_python_data(node, branch, node.type)
Vtrue, Itrue, Iinjtrue, S1true, S2true, Sinjtrue, num_iter = nv_powerflow.solve(branch, node)
Vtrue_cim, Itrue_cim, Iinjtrue_cim, S1true_cim, S2true_cim, Sinjtrue_cim, num_iter_cim = nv_powerflow_cim.solve(system)
""" Write here the indexes of the nodes/branches where the measurements are"""
V_idx = np.array([1,11,55])
I_idx = np.array([13,36])
Sinj_idx = np.linspace(2,node.num,node.num-1)
S1_idx = np.array([1,11,28,55,59])
S2_idx = np.array([10,54])
Ipmu_idx = np.array([1,11,55])
Vpmu_idx = np.array([13,36])
""" Write here the percent uncertainties of the measurements"""
V_unc = 1
I_unc = 2
Sinj_unc = 5
S_unc = 2
Pmu_mag_unc = 0.7
Pmu_phase_unc = 0.7
V = Measurements(V_idx,V_unc)
I = Measurements(I_idx,I_unc)
Sinj = Measurements(Sinj_idx,Sinj_unc)
S1 = Measurements(S1_idx,S_unc)
S2 = Measurements(S2_idx,S_unc)
Ipmu_mag = Measurements(Ipmu_idx,Pmu_mag_unc)
Ipmu_phase = Measurements(Ipmu_idx,Pmu_phase_unc)
Vpmu_mag = Measurements(Vpmu_idx,Pmu_mag_unc)
Vpmu_phase = Measurements(Vpmu_idx,Pmu_phase_unc)
meas = Measurement_set(V, I, Sinj, S1, S2, Ipmu_mag, Ipmu_phase, Vpmu_mag, Vpmu_phase)
zdata = Zdata_init(meas)
zdatameas = Zdata_init(meas)
zdata_cim = copy.deepcopy(zdata)
zdatameas_cim = copy.deepcopy(zdatameas)
zdata, zdatameas = py_95bus_meas_data.Zdatatrue_creation(zdata, zdatameas, meas, branch, Vtrue, Itrue, Iinjtrue, S1true, S2true, Sinjtrue)
zdata_cim, zdatameas_cim = cim_py_95bus_meas_data.Zdatatrue_creation(zdata_cim, zdatameas_cim, meas, system, Vtrue_cim, Itrue_cim, Iinjtrue_cim, S1true_cim, S2true_cim, Sinjtrue_cim)
#compare results:
print("Vtrue==Vtrue_cim?: " + str((Vtrue==Vtrue_cim).all()))
print("Itrue==Itrue_cim?: " + str((Itrue==Itrue_cim).all()))
print("Iinjtrue==Iinjtrue_cim?: " + str((Iinjtrue==Iinjtrue_cim).all()))
print("S1true==S1true_cim?: " + str((S1true==S1true_cim).all()))
print("S2true==S2true_cim?: " + str((S2true==S2true_cim).all()))
print("Sinjtrue==Sinjtrue_cim?: " + str((Sinjtrue==Sinjtrue_cim).all()))
print("zdata.mtype==zdata_cim.mtype?: {}".format(np.array_equal(zdata.mtype, zdata_cim.mtype)))
print("zdata.mval==zdata_cim.mval?: {}".format(np.array_equal(zdata.mval, zdata_cim.mval)))
print("zdatameas.mval==zdatameas_cim.mval?: {}".format(np.array_equal(zdatameas.mval, zdatameas_cim.mval)))
#print("zdata.mbranch==zdata_cim.mbranch?: {}".format(np.array_equal(zdata.mbranch, zdata_cim.mbranch)))
print("zdata.mfrom==zdata_cim.mfrom?: {}".format(np.array_equal(zdata.mfrom, zdata_cim.mfrom)))
#print(zdata.mfrom-zdata_cim.mfrom)
print("zdata.mto==zdata_cim.mto?: {}".format(np.array_equal(zdata.mto, zdata_cim.mto)))
#print(zdata.mto-zdata_cim.mto)
print("zdata.mstddev==zdata_cim.mstddev?: {}".format(np.array_equal(zdata.mstddev, zdata_cim.mstddev)))
MC_trials = 1
iter_counter = np.zeros(MC_trials)
Ymatr, Adj = nv_powerflow.Ymatrix_calc(branch,node)
Ymatr_cim, Adj_cim = network.Ymatrix_calc(system)
for mciter in range(0,MC_trials):
zdatameas = py_95bus_meas_data.Zdatameas_creation(zdata, zdatameas)
zdatameas_cim=copy.deepcopy(zdatameas) #to compare use the same random normal (Gaussian) distribution.
#zdatameas_cim = cim_py_95bus_meas_data.Zdatameas_creation(zdata_cim, zdatameas_cim)
Vest, Iest, Iinjest, S1est, S2est, Sinjest = nv_state_estimator.DsseCall(branch, node, zdatameas, Ymatr, Adj)
Vest_cim, Iest_cim, Iinjest_cim, S1est_cim, S2est_cim, Sinjest_cim = nv_state_estimator_cim.DsseCall(system, zdatameas_cim, Ymatr_cim, Adj_cim)
print("Ymatr==Ymatr_cim?: {}".format(np.array_equal(Ymatr, Ymatr_cim)))
print("Adj==Adj_cim?: {}".format(np.array_equal(Adj, Adj_cim)))
print("Vest==Vest_cim?: {}".format(np.array_equal(Vest.complex, Vest_cim.complex)))
#print(Vest.complex - Vest_cim.complex)
print("Iest==Iest_cim?: {}".format(np.array_equal(Iest.complex, Iest_cim.complex)))
print("Iinjest==Iinjest_cim?: {}".format(np.array_equal(Iinjest.complex, Iinjest_cim.complex)))
print("S1est==S1est_cim?: {}".format(np.array_equal(S1est.complex, S1est_cim.complex)))
print("S2est==S2est_cim?: {}".format(np.array_equal(S2est.complex, S2est_cim.complex)))
print("Sinjest==Sinjest_cim?: {}".format(np.array_equal(Sinjest.complex, Sinjest_cim.complex)))
\ No newline at end of file
Supports Markdown
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