Commit 76dc9bd6 authored by Markus Mirz's avatar Markus Mirz
Browse files

added package to acs namespace

parent de27f84e
__all__ = ["bc_powerflow", "bc_state_estimator", "nv_powerflow", "nv_state_estimator"]
\ No newline at end of file
......@@ -125,7 +125,7 @@ def BC_power_flow(branches, nodes):
Sinj_rx = np.multiply(V, np.conj(Iinj))
Sinj = np.real(Sinj_rx) + 1j * np.imag(Sinj_rx)
S1 = np.multiply(V[branchs.start-1], np.conj(I))
S2 = np.multiply(V[branchs.end-1], np.conj(I))
S1 = np.multiply(V[branches.start-1], np.conj(I))
S2 = np.multiply(V[branches.end-1], np.conj(I))
return V, I, Iinj, S1, S2, Sinj, num_iter
\ No newline at end of file
......@@ -58,10 +58,10 @@ def NV_power_flow(branch, node):
epsilon = 10**(-10)
diff = 5
V = np.ones(nodes.num) + 1j* np.zeros(nodes.num)
V = np.ones(node.num) + 1j* np.zeros(node.num)
num_iter = 0
State = np.ones(2*branches.num)
State = np.ones(2*branch.num)
State = np.concatenate((np.array([1,0]),State),axis=0)
while diff > epsilon:
......@@ -69,20 +69,20 @@ def NV_power_flow(branch, node):
i = k-1
m = 2*i
i2 = i + node.num
if nodes.type[i] == 'slack':
if node.type[i] == 'slack':
h[m] = np.inner(H[m],State)
h[m+1] = np.inner(H[m+1],State)
elif nodes.type[i] == 'PQ':
z[m] = (node.pwr_flow_values[1][i]*V.real[i] + node.pwr_flow_values[2][i]*V.imag[i])/(V.mag[i]**2)
z[m+1] = (node.pwr_flow_values[1][i]*V.imag[i] - node.pwr_flow_values[2][i]*V.real[i])/(V.mag[i]**2)
elif node.type[i] == 'PQ':
z[m] = (node.pwr_flow_values[1][i]*V[i].real + node.pwr_flow_values[2][i]*V[i].imag)/(np.abs(V[i])**2)
z[m+1] = (node.pwr_flow_values[1][i]*V[i].imag - node.pwr_flow_values[2][i]*V[i].real)/(np.abs(V[i])**2)
h[m] = np.inner(H[m],State)
h[m+1] = np.inner(H[m+1],State)
elif nodes.type[i] == 'PV':
z[m] = (node.pwr_flow_values[1][i]*V.real[i] + node.pwr_flow_values[2][i]*V.imag[i])/(V.mag[i]**2)
elif node.type[i] == 'PV':
z[m] = (node.pwr_flow_values[1][i]*V[i].real + node.pwr_flow_values[2][i]*V[i].imag)/(np.abs(V[i])**2)
h[m] = np.inner(H[m],State)
h[m+1] = V.mag[i]
H[m+1][i] = np.cos(V.phase[i])
H[m+1][i2] = np.sin(V.phase[i])
h[m+1] = np.abs(V[i])
H[m+1][i] = np.cos(np.angle(V[i]))
H[m+1][i2] = np.sin(np.angle(V[i]))
r = np.subtract(z,h)
Hinv = np.linalg.inv(H)
......@@ -115,8 +115,8 @@ def NV_power_flow(branch, node):
Sinj_rx = np.multiply(V, np.conj(Iinj))
Sinj = np.real(Sinj_rx) + 1j * np.imag(Sinj_rx)
S1 = np.multiply(V[branchs.start-1], np.conj(I))
S2 = - np.multiply(V[branchs.end-1], np.conj(I))
S1 = np.multiply(V[branch.start-1], np.conj(I))
S2 = - np.multiply(V[branch.end-1], np.conj(I))
return V, I, Iinj, S1, S2, Sinj, num_iter
\ No newline at end of file
import numpy
import math
import matplotlib as plt
import Network_data
import powerflow
from acs.state_estimation import *
import py_95bus_network_data
class PerUnit:
def __init__(self, S, V):
......@@ -17,6 +18,6 @@ V = (11*(10**3))/math.sqrt(3)
slackV = 1.02
Base = PerUnit(S,V)
branch, node = Network_data.Network_95_nodes(Base, slackV)
branch, node = py_95bus_network_data.Network_95_nodes(Base, slackV)
Vtrue, Itrue, Iinjtrue, S1true, S2true, Sinjtrue, num_iter = powerflow.BC_power_flow(branch, node)
\ No newline at end of file
Vtrue, Itrue, Iinjtrue, S1true, S2true, Sinjtrue, num_iter = bc_powerflow.BC_power_flow(branch, node)
\ No newline at end of file
def Zdatatrue_creation(zdata, zdatameas, meas, branch, V, I, Iinj, S1, S2, Sinj):
""" It gives the true values at the measurement points."""
import numpy
import numpy as np
t1 = numpy.ones(meas.V.num)
t2 = 2*numpy.ones(meas.Sinj.num)
t3 = 3*numpy.ones(meas.Sinj.num)
t4 = 4*numpy.ones(meas.S1.num + meas.S2.num)
t5 = 5*numpy.ones(meas.S1.num + meas.S2.num)
t6 = 6*numpy.ones(meas.I.num)
t7 = 7*numpy.ones(meas.Vpmu_mag.num)
t8 = 8*numpy.ones(meas.Vpmu_phase.num)
t9 = 9*numpy.ones(meas.Ipmu_mag.num)
t10 = 10*numpy.ones(meas.Ipmu_phase.num)
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 = numpy.concatenate((t1,t2,t3,t4,t5,t6,t7,t8,t9,t10),axis=0)
zdatameas.mtype = numpy.concatenate((t1,t2,t3,t4,t5,t6,t7,t8,t9,t10),axis=0)
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 = V.mag[meas.V.index-1]
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 = I.mag[meas.I.index-1]
z7 = V.mag[meas.Vpmu_mag.index-1]
z8 = V.phase[meas.Vpmu_phase.index-1]
z9 = I.mag[meas.Ipmu_mag.index-1]
z10 = I.phase[meas.Ipmu_phase.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 = numpy.concatenate((z1,z2,z3,z4_1,z4_2,z5_1,z5_2,z6,z7,z8,z9,z10),axis=0)
zdata.mval = np.concatenate((z1,z2,z3,z4_1,z4_2,z5_1,z5_2,z6,z7,z8,z9,z10),axis=0)
b1 = numpy.zeros(meas.V.num)
b2 = numpy.zeros(meas.Sinj.num)
b3 = numpy.zeros(meas.Sinj.num)
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 = numpy.zeros(meas.Vpmu_mag.num)
b8 = numpy.zeros(meas.Vpmu_phase.num)
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 = numpy.concatenate((b1,b2,b3,b4_1,b4_2,b5_1,b5_2,b6,b7,b8,b9,b10),axis=0)
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 = numpy.concatenate((b1,b2,b3,b4_1,b4_2,b5_1,b5_2,b6,b7,b8,b9,b10),axis=0)
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
......@@ -62,56 +62,56 @@ def Zdatatrue_creation(zdata, zdatameas, meas, branch, V, I, Iinj, S1, S2, Sinj)
fr9 = branch.start[meas.Ipmu_mag.index-1]
fr10 = branch.start[meas.Ipmu_phase.index-1]
zdata.mfrom = numpy.concatenate((fr1,fr2,fr3,fr4_1,fr4_2,fr5_1,fr5_2,fr6,fr7,fr8,fr9,fr10),axis=0)
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 = numpy.concatenate((fr1,fr2,fr3,fr4_1,fr4_2,fr5_1,fr5_2,fr6,fr7,fr8,fr9,fr10),axis=0)
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 = numpy.zeros(meas.V.num)
to2 = numpy.zeros(meas.Sinj.num)
to3 = numpy.zeros(meas.Sinj.num)
to1 = np.zeros(meas.V.num)
to2 = np.zeros(meas.Sinj.num)
to3 = np.zeros(meas.Sinj.num)
to4_1 = branch.end[meas.S1.index-1]
to4_2 = branch.start[meas.S2.index-1]
to5_1 = branch.end[meas.S1.index-1]
to5_2 = branch.start[meas.S2.index-1]
to6 = branch.end[meas.I.index-1]
to7 = numpy.zeros(meas.Vpmu_mag.num)
to8 = numpy.zeros(meas.Vpmu_phase.num)
to7 = np.zeros(meas.Vpmu_mag.num)
to8 = np.zeros(meas.Vpmu_phase.num)
to9 = branch.end[meas.Ipmu_mag.index-1]
to10 = branch.end[meas.Ipmu_phase.index-1]
zdata.mto = numpy.concatenate((to1,to2,to3,to4_1,to4_2,to5_1,to5_2,to6,to7,to8,to9,to10),axis=0)
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 = numpy.concatenate((to1,to2,to3,to4_1,to4_2,to5_1,to5_2,to6,to7,to8,to9,to10),axis=0)
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 = numpy.absolute(z2*(meas.Sinj.unc/300))
d3 = numpy.absolute(z3*(meas.Sinj.unc/300))
d4_1 = numpy.absolute(z4_1*(meas.S1.unc/300))
d4_2 = numpy.absolute(z4_2*(meas.S2.unc/300))
d5_1 = numpy.absolute(z5_1*(meas.S1.unc/300))
d5_2 = numpy.absolute(z5_2*(meas.S2.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)*numpy.ones(meas.Vpmu_phase.num)
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)*numpy.ones(meas.Ipmu_phase.num)
d10 = (meas.Ipmu_phase.unc/300)*np.ones(meas.Ipmu_phase.num)
zdata.mstddev = numpy.concatenate((d1,d2,d3,d4_1,d4_2,d5_1,d5_2,d6,d7,d8,d9,d10),axis=0)
zdatameas.mstddev = numpy.concatenate((d1,d2,d3,d4_1,d4_2,d5_1,d5_2,d6,d7,d8,d9,d10),axis=0)
# idx = numpy.where(zdata.mstddev<10**(-6))
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
import numpy as np
zmeas = zdata.mval
zdev = zdata.mstddev
err_pu = numpy.random.normal(0,1,len(zmeas))
zdatameas.mval = zmeas + numpy.multiply(zdev,err_pu)
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
import numpy
import numpy as np
import math
import matplotlib as plt
import Network_data
import Power_flow
import Meas_data
import NvStateEstimator
#import BcStateEstimator
import matplotlib.pyplot as plt
from acs.state_estimation import *
from acs.state_estimation import Power_flow
# scenario data
import py_95bus_network_data
import py_95bus_meas_data
class PerUnit:
def __init__(self, S, V):
......@@ -35,12 +36,12 @@ class Measurement_set:
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 = numpy.zeros(nmeas)
self.mval = numpy.zeros(nmeas)
self.mbranch = numpy.zeros(nmeas)
self.mfrom = numpy.zeros(nmeas)
self.mto = numpy.zeros(nmeas)
self.mstddev = numpy.zeros(nmeas)
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)
......@@ -48,18 +49,18 @@ V = (11*(10**3))/math.sqrt(3)
slackV = 1.02
Base = PerUnit(S,V)
branch, node = Network_data.Network_95_nodes(Base, slackV)
branch, node = py_95bus_network_data.Network_95_nodes(Base, slackV)
Vtrue, Itrue, Iinjtrue, S1true, S2true, Sinjtrue, num_iter = Power_flow.NV_power_flow(branch, node)
Vtrue, Itrue, Iinjtrue, S1true, S2true, Sinjtrue, num_iter = bc_powerflow.BC_power_flow(branch, node)
""" Write here the indexes of the nodes/branches where the measurements are"""
V_idx = numpy.array([1,11,55])
I_idx = numpy.array([13,36])
Sinj_idx = numpy.linspace(2,node.num,node.num-1)
S1_idx = numpy.array([1,11,28,55,59])
S2_idx = numpy.array([10,54])
Ipmu_idx = numpy.array([1,11,55])
Vpmu_idx = numpy.array([13,36])
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
......@@ -82,26 +83,26 @@ meas = Measurement_set(V, I, Sinj, S1, S2, Ipmu_mag, Ipmu_phase, Vpmu_mag, Vpmu_
zdata = Zdata_init(meas)
zdatameas = Zdata_init(meas)
zdata, zdatameas = Meas_data.Zdatatrue_creation(zdata, zdatameas, meas, branch, Vtrue, Itrue, Iinjtrue, S1true, S2true, Sinjtrue)
zdata, zdatameas = py_95bus_meas_data.Zdatatrue_creation(zdata, zdatameas, meas, branch, Vtrue, Itrue, Iinjtrue, S1true, S2true, Sinjtrue)
MC_trials = 1
#Iest_matrix = numpy.zeros((MC_trials,branch.num),dtype=numpy.complex_)
#Vest_matrix = numpy.zeros((MC_trials,node.num),dtype=numpy.complex_)
#S1est_matrix = numpy.zeros((MC_trials,branch.num),dtype=numpy.complex_)
#S2est_matrix = numpy.zeros((MC_trials,branch.num),dtype=numpy.complex_)
#Sinjest_matrix = numpy.zeros((MC_trials,node.num),dtype=numpy.complex_)
iter_counter = numpy.zeros(MC_trials)
#Iest_matrix = np.zeros((MC_trials,branch.num),dtype=np.complex_)
#Vest_matrix = np.zeros((MC_trials,node.num),dtype=np.complex_)
#S1est_matrix = np.zeros((MC_trials,branch.num),dtype=np.complex_)
#S2est_matrix = np.zeros((MC_trials,branch.num),dtype=np.complex_)
#Sinjest_matrix = np.zeros((MC_trials,node.num),dtype=np.complex_)
iter_counter = np.zeros(MC_trials)
Ymatr, Adj = Power_flow.Ymatrix_calc(branch,node)
Ymatr, Adj = nv_powerflow.Ymatrix_calc(branch,node)
for mciter in range(0,MC_trials):
if mciter%10 == 0:
print(mciter)
zdatameas = Meas_data.Zdatameas_creation(zdata, zdatameas)
zdatameas = py_95bus_meas_data.Zdatameas_creation(zdata, zdatameas)
Vest, Iest, Iinjest, S1est, S2est, Sinjest = NvStateEstimator.DsseCall(branch, node, zdatameas, Ymatr, Adj)
Vest, Iest, Iinjest, S1est, S2est, Sinjest = nv_state_estimator.DsseCall(branch, node, zdatameas, Ymatr, Adj)
# Iest_matrix[mciter] = Iest.complex
# Vest_matrix[mciter] = Vest.complex
......
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