Commit dd495ff5 authored by Steffen Vogel's avatar Steffen Vogel 🎅🏼

initial data from Alberto

parents
# Copyright © 2016-2017, RWTH Aachen University
# Authors: Andrea Angioni
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# A copy of the GNU General Public License is in the LICENSE file
# in the top level directory of this source tree.
# coding=utf8
import psycopg2,time,datetime
import random
import numpy
import cmath
import math
from datetime import datetime,timedelta
import sys, getopt
def testDSSE():
nameDB = "Cloud_DB"
userDB = "aan-iko"
passDB = "password"
IPDB = "134.130.169.25"
portDB = "5432"
groupgrid = "1"
try:
opts, args = getopt.getopt(sys.argv[1:], "hn:u:p:i:g:", ["nameDB=", "userDB=", "passDB=", "IPDB=", "groupgrid="])
except getopt.GetoptError:
print("not recognized input, please type: \n function.py -n <nameDB> -u <userDB> -p <passDB> -i <IP DB> -g <groupgrid>")
sys.exit(2)
for opt, arg in opts:
if opt == "-h":
print("please type: \n function.py -n <nameDB> -u <userDB> -p <passDB> -i <IP DB>")
sys.exit()
elif opt in ("-n", "--nameDB"):
nameDB = arg
print("name of DB is "+nameDB+"")
elif opt in ("-u", "--userDB"):
userDB = arg
print("user name for DB is "+userDB+"")
elif opt in ("-p", "--passDB"):
passDB = arg
print("pasword for DB is *not printed*")
elif opt in ("-i", "--IPDB"):
print("IP address for DB is "+IPDB+"")
IPDB = arg
elif opt in ("-g", "--groupgrid"):
print("Group grid is "+groupgrid+"")
groupgrid = arg
#print("The input parameters are: name of DB "+nameDB+", user name "+userDB+", pasword "+passDB+", IP address "+IPDB+", group_grid "+groupgrid+"")
databaseconn1 = databaseconn(nameDB,userDB,passDB,IPDB,portDB) #class for DB connection
conn = connDB(databaseconn1) #function for connection to DB
#griddata1 = readGridData(conn,groupgrid) #function to read grid data from DB table and rearrange in classes for python
griddata1 = readGridDataFLISR(conn) #function to read grid data, slightly adapted for FLISR
setdsse1=initDSSE(conn)
while True:
ttot0=0;
measdata1 = readMeasData(conn,setdsse1,griddata1,groupgrid) #function to retrieve measurements from DB
dssedata1 = DSSEVRI1 (griddata1,measdata1) #function with the actual state estimator, requires grid data and measurements as input
ttot3 = writeDB(conn,dssedata1)
print("time to read measurements "+str(measdata1.time)+" s;\ntime to run DSSE "+str(dssedata1.time)+" s;\ntime to write into DB "+str(ttot3)+" s")
if (setdsse1.treal - dssedata1.time - measdata1.time -ttot3)>0:
time.sleep(setdsse1.treal - dssedata1.time - measdata1.time -ttot3) # it waits a time equal to the difference between the real time object and the time required b the DSSE, this should allow to run in real time every second
if setdsse1.treal > setdsse1.trealset and( dssedata1.time + measdata1.time + ttot3) < float(setdsse1.treal/2): #if the real time target was changed but we have again very short computation time we bring it back to the original
setdsse1.treal=math.ceil( dssedata1.time + measdata1.time + ttot3)
print("elapsed time is "+str(dssedata1.time + measdata1.time + ttot3)+", real time target is again "+ str(setdsse1.treal) +".\n")
SQLtext0=""
SQLtext0+="INSERT INTO management.algorithm_log(algorithm_id,timestamp,log_message) "
SQLtext0 += " VALUES('1','"+str(datetime.utcnow())+"','DSSE - CHANGE REAL TIME OBJECT TO "+ str(setdsse1.treal) +"');\n"
SQLtext0 += "UPDATE management.algorithm_parameter SET parameter_value ="+ str(setdsse1.treal) +" WHERE "
SQLtext0 += "'parameter_name' ='DSSE_T_REAL_TIME' "
else:
setdsse1.treal=math.ceil( dssedata1.time + measdata1.time + ttot3)
print("elapsed time is "+str(dssedata1.time + measdata1.time + ttot3)+", real time target readapted to "+ str(setdsse1.treal) +".\n")
SQLtext0=""
SQLtext0+="INSERT INTO management.algorithm_log(algorithm_id,timestamp,log_message) "
SQLtext0 += " VALUES('1','"+str(datetime.utcnow())+"','DSSE - CHANGE REAL TIME OBJECT TO "+ str(setdsse1.treal) +"');\n"
SQLtext0 += "UPDATE management.algorithm_parameter SET parameter_value ="+ str(setdsse1.treal) +" WHERE "
SQLtext0 += "'parameter_name' ='DSSE_T_REAL_TIME' "
cur = conn.cursor()
cur.execute(SQLtext0)
conn.commit()
#print("elapsed time is "+str(dssedata1.time + measdata1.time + ttot3)+"")
def readPowerData(conn,setdsse,griddata,groupgrid): #function to read power measruements, particularly needed by FLISR
tlarge=60; #if there are few measurements we extend the search in the last minute
start2=datetime.utcnow()
print("current time: "+str(start2)+"")
cur = conn.cursor()
measdata1 = measdata()
#1st part. we read the measurements that were updated before the last execution
Nmeas = 0;
ref_angle = 0;
mtype1=['9','10','11','12']
for index1,row1 in enumerate(mtype1):
mtype=mtype1[index1]
SQLtext = selectSQL(mtype,str(datetime.utcnow() - timedelta(seconds=setdsse.treal) ),str(datetime.utcnow()),groupgrid)
#print(SQLtext)
try:
cur.execute(SQLtext)
conn.commit()
rows = cur.fetchall()
for index,row in enumerate(rows):
acc1 = float(row[3]);
if row[1]=='AN':
ph1=1
elif row[1]=='BN':
ph1=2
elif row[1]=='CN':
ph1=3
if int(mtype)==9 or int(mtype)==10 or int(mtype)==11 or int(mtype)==12 or int(mtype)==13 :
measdata1.addmeasurementinstance(Nmeas,row[4],ph1,int(row[5]),acc1,row[2]/griddata.pbase) #index represents the itemid for both nodes and voltages, both of them starts from zero
Nmeas=Nmeas+1;
except psycopg2.OperationalError as e:
print('Unable to execute query!\n{0}').format(e)
sys.exit(1)
measdata1.nmeas = Nmeas #this indicates the total number of measurements available
#2nd part. we verify if we have enough measurements, if not we look more in the past
if (griddata.nnode/(measdata1.nmeas+1)) > setdsse.xmeasmin: #if there are few measurements we extend the search in the last minute
print("Not enough measurements in real time, looking back at "+str(setdsse.treal*setdsse.xmeaspast)+" seconds")
SQLtext0=""
SQLtext0+="INSERT INTO management.flag(ftype_id,flag_timestamp,flag_message) "
SQLtext0 += " VALUES((SELECT ftype_id FROM management.flag_index WHERE ftype_name = 'SE_FM'),'"+str(datetime.utcnow())+"','FEW measurements available for DSSE');\n"
cur.execute(SQLtext0)
conn.commit()
measdata1 = measdata()
Nmeas=0;
for index1,row1 in enumerate(mtype1):
mtype=mtype1[index1]
SQLtext = selectSQL(mtype,str(datetime.utcnow() - timedelta(seconds=setdsse.xmeaspast) ),str(datetime.utcnow()),groupgrid)
try:
cur.execute(SQLtext)
conn.commit()
rows = cur.fetchall()
for index,row in enumerate(rows):
acc1=float(row[3]);
if row[1]=='AN':
ph1=1
elif row[1]=='BN':
ph1=2
elif row[1]=='CN':
ph1=3
if int(mtype)==9 or int(mtype)==10 or int(mtype)==11 or int(mtype)==12 or int(mtype)==13 :
measdata1.addmeasurementinstance(Nmeas,row[4],ph1,int(row[5]),acc1,row[2]/griddata.pbase) #index represents the itemid for both nodes and voltages, both of them starts from zero
Nmeas=Nmeas+1;
except psycopg2.OperationalError as e:
print('Unable to execute query!\n{0}').format(e)
sys.exit(1)
measdata1.nmeas=Nmeas
print("number of real measurements is "+str(Nmeas)+"\n")
#3rd part. we also add the pseudo measurements from the DB
profile_id=['1','2']#such profile ids, can be chosen depending on which type of profile is needed
SQLtext = ""
SQLtext += "SELECT value,accuracy FROM measurecommand.pseudo_measurements_values WHERE (profile_id="+profile_id[0]+" OR profile_id="+profile_id[1]+") "
SQLtext += "AND (60*extract(HOUR from timestamp '"+str(datetime.utcnow())+"')+extract(MINUTE from timestamp '"+str(datetime.utcnow())+"')) > (60*extract(HOUR from timestamp1)+extract(MINUTE from timestamp1))"
SQLtext += "AND (60*extract(HOUR from timestamp '"+str(datetime.utcnow())+"')+extract(MINUTE from timestamp '"+str(datetime.utcnow())+"')) < (60*extract(HOUR from timestamp2)+extract(MINUTE from timestamp2)) ORDER BY profile_id ASC"
cur.execute(SQLtext)
conn.commit()
rows1 = cur.fetchall()
ph1=1 # for the moment just on phase AN
if cur.rowcount > 0:
Lpseudo_value=float(rows1[0][0])#pseudo values for loads
Lpseudo_acc=float(rows1[0][1])
Gpseudo_value=float(rows1[1][0])#pseudo values for generators
Gpseudo_acc=float(rows1[1][1])
for x in range(griddata.nnode):
SQLtext = "SELECT nominals FROM networktopology.nodes RIGHT JOIN networktopology.injections ON networktopology.injections.node_id = networktopology.nodes.node_id WHERE (injection_type = '1' or injection_type = '2') AND networktopology.nodes.node_id = '"+str(griddata.node[x].nodeid)+"';"
cur.execute(SQLtext)
conn.commit()
rows1 = cur.fetchall()
SQLtext = "SELECT nominals FROM networktopology.nodes RIGHT JOIN networktopology.injections ON networktopology.injections.node_id = networktopology.nodes.node_id WHERE injection_type = '3' AND networktopology.nodes.node_id = '"+str(griddata.node[x].nodeid)+"' ;"
cur.execute(SQLtext)
conn.commit()
rows2 = cur.fetchall()
PLoadtot=0.0
PGentot=0.0
for row in rows1:
PLoadtot=PLoadtot+0.333*Lpseudo_value*float(row[0])
for row in rows2:
PGentot=PGentot+0.333*Gpseudo_value*float(row[0])
if x > 0:
measdata1.addmeasurementinstance(Nmeas,9,ph1,griddata.node[x].nodeid,numpy.sqrt(Lpseudo_acc**2 + Gpseudo_acc**2),(PLoadtot+PGentot)*0.97/griddata.pbase)
measdata1.addmeasurementinstance(Nmeas+1,11,ph1,griddata.node[x].nodeid,numpy.sqrt(Lpseudo_acc**2 + Gpseudo_acc**2),PLoadtot*0.2431/griddata.pbase) #reactive power of load assuming cosfi=0.97, for the generator we assume reactive power equal to zero
Nmeas=Nmeas+2;
measdata1.nmeas=Nmeas #this indicates the total number of measurements available
print("number of real + pseudo measurements is "+str(Nmeas)+"\n")
ttot2=datetime.utcnow()-start2
measdata1.time=ttot2.total_seconds()
#for x in range(measdata1.nmeas):
# print("measurement type = ",measdata1.measurementinstance[x].mtype,"; item_id = ",measdata1.measurementinstance[x].itemid,"; accuracy is = ",measdata1.measurementinstance[x].measurementaccuracy, "; value = ",measdata1.measurementinstance[x].measurementvalue)
return measdata1
def readMeasData(conn,setdsse,griddata,groupgrid): #function to read measruement data and rearrange in measdata class
tlarge=60; #if there are few measurements we extend the search in the last minute
start2=datetime.utcnow()
print("current time: "+str(start2)+"")
cur = conn.cursor()
measdata1 = measdata()
#1st part. we read the measurements that were updated before the last execution
Nmeas = 0;
ref_angle = 0;
mtype1=['1','3','5','7','9','11','13','15']
for index1,row1 in enumerate(mtype1):
mtype=mtype1[index1]
#SQLtext = selectSQL(mtype,str(datetime.utcnow() - timedelta(seconds=setdsse.treal) ),str(datetime.utcnow()),groupgrid)
SQLtext = selectSQLFLISR(mtype,str(datetime.utcnow() - timedelta(seconds=setdsse.treal) ),str(datetime.utcnow()))
#print(SQLtext)
try:
cur.execute(SQLtext)
conn.commit()
rows = cur.fetchall()
for index,row in enumerate(rows):
acc1 = float(row[3]);
if row[1]=='AN':
ph1=1
elif row[1]=='BN':
ph1=2
elif row[1]=='CN':
ph1=3
if int(mtype)==1: #V magnitude single phase
measdata1.addmeasurementinstance(Nmeas,row[4],ph1,int(row[5]),acc1,row[2]/( griddata.node[int(row[5])].vbase)) #index represents the itemid for both nodes and voltages, both of them starts from zero
elif int(mtype)==3: #V phase single phase
if int(row[5]) == 0:
ref_angle = row[2]; # reference node identified
measdata1.addmeasurementinstance(Nmeas,row[4],ph1,int(row[5]),acc1,row[2] - ref_angle)
elif int(mtype)==7: #I phase single phase
measdata1.addmeasurementinstance(Nmeas,row[4],ph1,int(row[5]),acc1,row[2] - ref_angle)
elif int(mtype)==5: #I magnitude single phase
measdata1.addmeasurementinstance(Nmeas,row[4],ph1,int(row[5]),acc1,row[2]* (griddata.node[griddata.line[int(row[5])].linefrom].vbase) / griddata.pbase) #index represents the itemid for both nodes and voltages, both of them starts from zero
elif int(mtype)==9 or int(mtype)==11 or int(mtype)==13 or int(mtype)==15:
measdata1.addmeasurementinstance(Nmeas,row[4],ph1,int(row[5]),acc1,row[2]/griddata.pbase) #index represents the itemid for both nodes and voltages, both of them starts from zero
else:
measdata1.addmeasurementinstance(Nmeas,row[4],ph1,int(row[5]),acc1,row[2]) #index represents the itemid for both nodes and voltages, both of them starts from zero
#measurementinstance(measurementid,mtype,phase,itemid,measurementaccuracy,measurementvalue):
Nmeas=Nmeas+1;
except psycopg2.OperationalError as e:
print('Unable to execute query!\n{0}').format(e)
sys.exit(1)
measdata1.nmeas = Nmeas #this indicates the total number of measurements available
#2nd part. we verify if we have enough measurements, if not we look more in the past
if (griddata.nnode/(measdata1.nmeas+1)) > setdsse.xmeasmin: #if there are few measurements we extend the search in the last minute
print("Not enough measurements in real time, looking back at "+str(setdsse.treal*setdsse.xmeaspast)+" seconds")
SQLtext0=""
SQLtext0+="INSERT INTO management.flag(ftype_id,flag_timestamp,flag_message) "
SQLtext0 += " VALUES((SELECT ftype_id FROM management.flag_index WHERE ftype_name = 'SE_FM'),'"+str(datetime.utcnow())+"','FEW measurements available for DSSE');\n"
cur.execute(SQLtext0)
conn.commit()
measdata1 = measdata()
Nmeas=0;
for index1,row1 in enumerate(mtype1):
mtype=mtype1[index1]
#SQLtext = selectSQL(mtype,str(datetime.utcnow() - timedelta(seconds=setdsse.xmeaspast) ),str(datetime.utcnow()),groupgrid)
SQLtext = selectSQLFLISR(mtype,str(datetime.utcnow() - timedelta(seconds=setdsse.xmeaspast) ),str(datetime.utcnow()))
try:
cur.execute(SQLtext)
conn.commit()
rows = cur.fetchall()
for index,row in enumerate(rows):
acc1=float(row[3]);
if row[1]=='AN':
ph1=1
elif row[1]=='BN':
ph1=2
elif row[1]=='CN':
ph1=3
if int(mtype)==1:
measdata1.addmeasurementinstance(Nmeas,row[4],ph1,int(row[5]),acc1,row[2]/griddata.node[int(row[5])].vbase ) #index represents the itemid for both nodes and voltages, both of them starts from zero
elif int(mtype)==3:
if int(row[5]) == 0 :
ref_angle = row[2]; # reference node identified
measdata1.addmeasurementinstance(Nmeas,row[4],ph1,int(row[5]),acc1,row[2] - ref_angle)
elif int(mtype)==7:
measdata1.addmeasurementinstance(Nmeas,row[4],ph1,int(row[5]),acc1,row[2] - ref_angle)
elif int(mtype)==5:
measdata1.addmeasurementinstance(Nmeas,row[4],ph1,int(row[5]),acc1,row[2]* (griddata.node[griddata.line[int(row[5])].linefrom].vbase ) / griddata.pbase) #index represents the itemid for both nodes and voltages, both of them starts from zero
elif int(mtype)==9 or int(mtype)==11 or int(mtype)==13 or int(mtype)==15:
measdata1.addmeasurementinstance(Nmeas,row[4],ph1,int(row[5]),acc1,row[2]/griddata.pbase) #index represents the itemid for both nodes and voltages, both of them starts from zero
else:
measdata1.addmeasurementinstance(Nmeas,row[4],ph1,int(row[5]),acc1,row[2]) #index represents the itemid for both nodes and voltages, both of them starts from zero
#measurementinstance(measurementid,mtype,phase,itemid,measurementaccuracy,measurementvalue):
Nmeas=Nmeas+1;
except psycopg2.OperationalError as e:
print('Unable to execute query!\n{0}').format(e)
sys.exit(1)
measdata1.nmeas=Nmeas
print("number of real measurements is "+str(Nmeas)+"\n")
#3rd part. we also add the pseudo measurements from the DB
profile_id=['1','2']#such profile ids, can be chosen depending on which type of profile is needed
SQLtext = ""
SQLtext += "SELECT value,accuracy FROM measurecommand.pseudo_measurements_values WHERE (profile_id="+profile_id[0]+" OR profile_id="+profile_id[1]+") "
SQLtext += "AND (60*extract(HOUR from timestamp '"+str(datetime.utcnow())+"')+extract(MINUTE from timestamp '"+str(datetime.utcnow())+"')) > (60*extract(HOUR from timestamp1)+extract(MINUTE from timestamp1))"
SQLtext += "AND (60*extract(HOUR from timestamp '"+str(datetime.utcnow())+"')+extract(MINUTE from timestamp '"+str(datetime.utcnow())+"')) < (60*extract(HOUR from timestamp2)+extract(MINUTE from timestamp2)) ORDER BY profile_id ASC"
cur.execute(SQLtext)
conn.commit()
rows1 = cur.fetchall()
ph1=1 # for the moment just on phase AN
if cur.rowcount > 0:
Lpseudo_value=float(rows1[0][0])#pseudo values for loads
Lpseudo_acc=float(rows1[0][1])
Gpseudo_value=float(rows1[1][0])#pseudo values for generators
Gpseudo_acc=float(rows1[1][1])
for x in range(griddata.nnode):
SQLtext = "SELECT nominals FROM networktopology.nodes RIGHT JOIN networktopology.injections ON networktopology.injections.node_id = networktopology.nodes.node_id WHERE (injection_type = '1' or injection_type = '2') AND networktopology.nodes.node_id = '"+str(griddata.node[x].nodeid)+"';"
cur.execute(SQLtext)
conn.commit()
rows1 = cur.fetchall()
SQLtext = "SELECT nominals FROM networktopology.nodes RIGHT JOIN networktopology.injections ON networktopology.injections.node_id = networktopology.nodes.node_id WHERE injection_type = '3' AND networktopology.nodes.node_id = '"+str(griddata.node[x].nodeid)+"' ;"
cur.execute(SQLtext)
conn.commit()
rows2 = cur.fetchall()
PLoadtot=0.0
PGentot=0.0
for row in rows1:
PLoadtot=PLoadtot+0.333*Lpseudo_value*float(row[0])
for row in rows2:
PGentot=PGentot+0.333*Gpseudo_value*float(row[0])
if x > 0:
measdata1.addmeasurementinstance(Nmeas,9,ph1,griddata.node[x].nodeid,numpy.sqrt(Lpseudo_acc**2 + Gpseudo_acc**2),(PLoadtot+PGentot)*0.97/griddata.pbase)
measdata1.addmeasurementinstance(Nmeas+1,11,ph1,griddata.node[x].nodeid,numpy.sqrt(Lpseudo_acc**2 + Gpseudo_acc**2),PLoadtot*0.2431/griddata.pbase) #reactive power of load assuming cosfi=0.97, for the generator we assume reactive power equal to zero
Nmeas=Nmeas+2;
measdata1.nmeas=Nmeas #this indicates the total number of measurements available
print("number of real + pseudo measurements is "+str(Nmeas)+"\n")
#4th part. we insert some virtual measurements. Basically it is the knowledge that the voltage magnitude and phase angle are normally in a certain range (nominal + 20%)
#This step will always make the system observable
for x in range(griddata.nnode):
measdata1.addmeasurementinstance(Nmeas,1,ph1,griddata.node[x].nodeid,0.2,1) #1 pu
measdata1.addmeasurementinstance(Nmeas+1,3,ph1,griddata.node[x].nodeid,0.2,0) #0 radians
Nmeas=Nmeas+2;
measdata1.nmeas=Nmeas #this indicates the total number of measurements available
print("number of total measurements is "+str(Nmeas)+"\n")
ttot2=datetime.utcnow()-start2
measdata1.time=ttot2.total_seconds()
#for x in range(measdata1.nmeas):
# print("measurement type = ",measdata1.measurementinstance[x].mtype,"; item_id = ",measdata1.measurementinstance[x].itemid,"; accuracy is = ",measdata1.measurementinstance[x].measurementaccuracy, "; value = ",measdata1.measurementinstance[x].measurementvalue)
return measdata1
def DSSEVRI1 (griddata,measdata):
# Distribution System State Estimation function
# requires grid data and measurements as input, it produces a class with the results of the estimation in output
# it is a Weighted Least Square estimator based on rectangular coordinates of voltage phasor
# measurements are converted eiher to voltage or current phasors. In this way the problem is linearized and the jacobian matrix
# does not have to be calculated at each iteration of the Netwon Rapson method.
start1 = datetime.now()
wmatrix1 = weightmatrixVRI1(griddata,measdata) #function to calculate the weight matrix
[jmatrix1,measdata] = jacobianmatrixVRI1(griddata,measdata) #function to calculate the Jacobian matrix
HW = numpy.transpose(jmatrix1).dot(wmatrix1)
G1 = HW.dot(jmatrix1); #Gain matrix
[jmatrix2,measdata2] = jacobianmatrixIRI1(griddata,measdata) #function to calculate the Jacobian matrix
HW2 = numpy.transpose(jmatrix2).dot(wmatrix1)
G2 = HW2.dot(jmatrix2); #Gain matrix
Vmagnstatus = numpy.ones((griddata.nnode,1)) #initialize the voltage magnitude and phase angle
Vphasestatus = numpy.zeros((griddata.nnode,1))
Vreal = numpy.ones((griddata.nnode,1)) #similarly we initialize the voltage real and imaginary parts
Vimag = numpy.zeros((griddata.nnode,1))
Imagn = numpy.zeros((griddata.nline,1))
Iphase = numpy.zeros((griddata.nline,1))
Imagnload = numpy.zeros((griddata.nnode,1))
Iphaseload = numpy.zeros((griddata.nnode,1))
max_delta=10; #the maximum delta is used as criteria to stop the Netwon Rapson iteration. We do here a dummy inizialitazion
iteration=0; #number of iterations of NR
while max_delta > 0.0000001 and iteration < 50: #conditions to stop NR
iteration=iteration+1;
#print("iteration number "+str(iteration)+";")
z = measurementvector1(griddata,measdata,Vreal,Vimag) #measurement vector (all measurements converted to voltages and currents)
hx = hxvector1(griddata,measdata,Vreal,Vimag)#vector with measurement quantity calculated based on the estimated state at the ith iteration
res = z - hx;#residual vector (to be minimized)
HWres = HW.dot(res)
delta = numpy.linalg.solve(G1,HWres) #the delta vector is calculated
if measdata.switchphmeas==0:
#jmatrix=numpy.delete(jmatrix, 1, 1)
Vreal[0,0]=Vreal[0,0]+delta[0,0] #it is summed to the current state. The imaginary part at the slack is by assumption always zero
Vmagnstatus[0,0]=abs(Vreal[0,0]+1j*Vimag[0,0])
n=0;
for x in range(1, griddata.nnode): #then we keep on summing delta to the other elements of the states
n=n+1;
Vreal[x,0]=Vreal[x,0]+delta[n,0]
n=n+1;
Vimag[x,0]=Vimag[x,0]+delta[n,0]
Vmagnstatus[x,0]=abs(Vreal[x,0]+1j*Vimag[x,0])
Vphasestatus[x,0]=cmath.phase(Vreal[x,0]+1j*Vimag[x,0])
else:
n=0;
for x in range(0, griddata.nnode): #then we keep on summing delta to the other elements of the states
Vreal[x,0]=Vreal[x,0]+delta[n,0]
n=n+1;
Vimag[x,0]=Vimag[x,0]+delta[n,0]
n=n+1;
Vmagnstatus[x,0]=abs(Vreal[x,0]+1j*Vimag[x,0])
Vphasestatus[x,0]=cmath.phase(Vreal[x,0]+1j*Vimag[x,0])
max_delta=max(abs(delta)) #we calculated the maximum element of delta
# print("Total number of iterations is " + str(iteration) +"")
# print("voltage magnitude is")
# print(Vmagnstatus)
# print("voltage phase angle is")
# print(Vphasestatus)
uncertaintytateV = calculatecovestimationV(G1,griddata,Vmagnstatus,Vphasestatus,measdata)#we calculate the estimated uncertainty of the state
(Ibranchireal,Ibranchiimag,Ibranchfreal,Ibranchfimag,Iloadreal,Iloadimag)=calccurrentsI(griddata,Vreal,Vimag)
# nstate=dssedata1.nstate
# dssedata1.nstate = dssedata1.nstate + 2*griddata.nnode+2*griddata.nline
# dssedata1.numnr = iteration
dssedata1=dssedata()#we write the DSSE results in a dedicated class
dssedata1.nstate=2*griddata.nnode+2*griddata.nline
dssedata1.numnr = iteration
nstate=0
for x in range(griddata.nnode):
Imagnload[x,0] = abs(Iloadreal[x,0]+1j*Iloadimag[x,0])
Iphaseload[x,0] = cmath.phase(Iloadreal[x,0]+1j*Iloadimag[x,0])
dssedata1.adddsseinstance(nstate,1,1,x,x,uncertaintytateV[2*x],Vmagnstatus[x,0]) #index represents the itemid for both nodes and voltages, both of them starts from zero
nstate=nstate+1
dssedata1.adddsseinstance(nstate,3,1,x,x,uncertaintytateV[2*x+1],Vphasestatus[x,0]) #index represents the itemid for both nodes and voltages, both of them starts from zero
nstate=nstate+1
#print("max err v.magnitude, bus "+str(x)+": "+str(300*uncertaintytateV[2*x])+" %")
#print("max err v.ph.angle, bus "+str(x)+": "+str(300*uncertaintytateV[2*x+1])+" crad")
for x in range(griddata.nline):
Imagn[x,0] = abs(Ibranchireal[x,0]+1j*Ibranchiimag[x,0])
Iphase[x,0] = cmath.phase(Ibranchireal[x,0]+1j*Ibranchiimag[x,0])
#G1inv = numpy.linalg.inv(G2) # this is covariance matrix of the estimated state in rectangular format
#G1inv = numpy.delete(G1inv, 0, 0)
#G1inv = numpy.delete(G1inv, 0, 0)
#G1inv = numpy.delete(G1inv, 0, 1)
#G1inv = numpy.delete(G1inv, 0, 1)
#var_mat=numpy.diag(G1inv); #variance output
#uncertaintytate_rect=numpy.sqrt(var_mat);
#print('uncertainty rect',uncertaintytate_rect*(4000000/(math.sqrt(3)*8660)))
uncertaintytateI = calculatecovestimationI(G2,griddata,Imagn,Iphase,measdata)#we calculate the estimated uncertainty of the state
#print('uncertainty I polar',uncertaintytateI*(4000000/(math.sqrt(3)*8660)))
for x in range(griddata.nline):
#print('line from',griddata.line[x].linefrom,'line to',griddata.line[x].lineto)
dssedata1.adddsseinstance(nstate,5,1,griddata.line[x].linefrom,griddata.line[x].lineto,uncertaintytateI[2*x],Imagn[x,0]) #the uncertainty of the branch current is yet not calculated
nstate=nstate+1
dssedata1.adddsseinstance(nstate,7,1,griddata.line[x].linefrom,griddata.line[x].lineto,uncertaintytateI[2*x+1],Iphase[x,0])
nstate=nstate+1
# print("line current magnitude is")
# print(Imagn)
# print("line current phase angle is")
# print(Iphase)
dssedata1.numnr=iteration
ttot1=datetime.now()-start1 #once we finish the estimation we calculate the required time
dssedata1.time=ttot1.total_seconds()
return dssedata1
def calccurrentsI(griddata,Vreal,Vimag): #function to calculate currents given two set of voltages
Ibranchireal=numpy.zeros((griddata.nline,1));
Ibranchiimag=numpy.zeros((griddata.nline,1));
Ibranchfreal=numpy.zeros((griddata.nline,1));
Ibranchfimag=numpy.zeros((griddata.nline,1));
Ibranchf=numpy.zeros((griddata.nline,1));
Iloadreal=numpy.zeros((griddata.nnode,1));
Iloadimag=numpy.zeros((griddata.nnode,1));
for x in range(griddata.nline):
puinvbase=(griddata.pbase/((griddata.node[griddata.line[x].linefrom].vbase)**2))**(-1);
# Ztemp = ((1/3) * (2*griddata.line[x].phaseP11[0].resistance + griddata.line[x].phaseP12[0].resistance )) +1j*((1/3) * (2*griddata.line[x].phaseP11[0].reactance + griddata.line[x].phaseP12[0].reactance ))
Ztemp = griddata.line[x].phaseP11[0].resistance +1j*griddata.line[x].phaseP11[0].reactance
Rinv = (puinvbase * (Ztemp)**(-1)).real
Xinv = (puinvbase * (Ztemp)**(-1)).imag
Gd=griddata.line[x].phaseP11[0].conductance
G0=griddata.line[x].phaseP12[0].conductance
Bd=griddata.line[x].phaseP11[0].susceptance
B0=griddata.line[x].phaseP12[0].susceptance
# Yshunt=puinvbase*((1/3)*(2*(1/(Gd+1j*Bd)) + (1/(G0+1j*B0))))**(-1)
Yshunt=(Gd+1j*Bd)*puinvbase
Ibranchi = 0.5*(Yshunt)* (Vreal[griddata.line[x].linefrom,0]+1j*Vimag[griddata.line[x].linefrom,0]) + (Rinv + 1j*Xinv)*((Vreal[griddata.line[x].linefrom,0]+1j*Vimag[griddata.line[x].linefrom,0]) - (Vreal[griddata.line[x].lineto ,0]+1j*Vimag[griddata.line[x].lineto ,0]))
Ibranchf = - 0.5*(Yshunt)* (Vreal[griddata.line[x].lineto ,0]+1j*Vimag[griddata.line[x].lineto ,0]) + (Rinv + 1j*Xinv)*((Vreal[griddata.line[x].linefrom,0]+1j*Vimag[griddata.line[x].linefrom,0]) - (Vreal[griddata.line[x].lineto ,0]+1j*Vimag[griddata.line[x].lineto ,0]))
# B = puinvbase * ((1/3) * (2*griddata.line[x].phaseP11[0].susceptance + griddata.line[x].phaseP12[0].susceptance ))
# G = puinvbase * ((1/3) * (2*griddata.line[x].phaseP11[0].conductance + griddata.line[x].phaseP12[0].conductance ))
# Ibranchi = 0.5*(G + 1j*B)* (Vreal[griddata.line[x].linefrom,0]+1j*Vimag[griddata.line[x].linefrom,0]) + (Rinv + 1j*Xinv)*((Vreal[griddata.line[x].linefrom,0]+1j*Vimag[griddata.line[x].linefrom,0]) - (Vreal[griddata.line[x].lineto ,0]+1j*Vimag[griddata.line[x].lineto ,0]))
# Ibranchf = - 0.5*(G + 1j*B)* (Vreal[griddata.line[x].lineto ,0]+1j*Vimag[griddata.line[x].lineto ,0]) + (Rinv + 1j*Xinv)*((Vreal[griddata.line[x].linefrom,0]+1j*Vimag[griddata.line[x].linefrom,0]) - (Vreal[griddata.line[x].lineto ,0]+1j*Vimag[griddata.line[x].lineto ,0]))
Ibranchireal[x,0]=Ibranchi.real;
Ibranchiimag[x,0]=Ibranchi.imag;
Ibranchfreal[x,0]=Ibranchf.real;
Ibranchfimag[x,0]=Ibranchf.imag;
Iloadreal[griddata.line[x].linefrom,0] = Iloadreal[griddata.line[x].linefrom,0] - Ibranchireal[x,0]
Iloadreal[griddata.line[x].lineto,0] = Iloadreal[griddata.line[x].lineto,0] + Ibranchfreal[x,0]
Iloadimag[griddata.line[x].linefrom,0] = Iloadimag[griddata.line[x].linefrom,0] - Ibranchiimag[x,0]
Iloadimag[griddata.line[x].lineto,0] = Iloadimag[griddata.line[x].lineto,0] + Ibranchfimag[x,0]
return Ibranchireal,Ibranchiimag,Ibranchfreal,Ibranchfimag,Iloadreal,Iloadimag
def hxvector1(griddata,measdata,Vreal,Vimag):#function to calculate the h(x) vector based on the current Netwon Raspon iteration status of voltages
hx = numpy.zeros((measdata.nmeas,1))
(Ibranchireal,Ibranchiimag,Ibranchfreal,Ibranchfimag,Iloadreal,Iloadimag)=calccurrentsI(griddata,Vreal,Vimag)
for x in range(measdata.nmeas):
if measdata.measurementinstance[x].mtype == 1 or measdata.measurementinstance[x].mtype == 2: #voltage magnitude measurements converted to real part
hx[x,0]=Vreal[measdata.measurementinstance[x].itemid,0]
if measdata.measurementinstance[x].mtype == 3 or measdata.measurementinstance[x].mtype == 4: #voltage phase angle measurements converted to imag part
hx[x,0]=Vimag[measdata.measurementinstance[x].itemid,0]
if measdata.measurementinstance[x].mtype == 5 or measdata.measurementinstance[x].mtype == 6: #current magnitude measurements converted to real part
hx[x,0]=Ibranchireal[measdata.measurementinstance[x].itemid,0];
if measdata.measurementinstance[x].mtype == 7 or measdata.measurementinstance[x].mtype == 8: #current phase angle measurements converted to imag part
hx[x,0]=Ibranchiimag[measdata.measurementinstance[x].itemid,0];
if measdata.measurementinstance[x].mtype == 9 or measdata.measurementinstance[x].mtype == 10: #active power injection measurements converted to current real part
hx[x,0]=Iloadreal[measdata.measurementinstance[x].itemid,0];
if measdata.measurementinstance[x].mtype == 11 or measdata.measurementinstance[x].mtype == 12: #reactive power injection measurements converted to current imaginary part
hx[x,0]=Iloadimag[measdata.measurementinstance[x].itemid,0];
if measdata.measurementinstance[x].mtype == 13 or measdata.measurementinstance[x].mtype == 14: #active power flow measurements converted to current real part
hx[x,0]=Ibranchireal[measdata.measurementinstance[x].itemid,0];
if measdata.measurementinstance[x].mtype == 15 or measdata.measurementinstance[x].mtype == 16: #reactive power flow measurements converted to current imaginary part
hx[x,0]=Ibranchiimag[measdata.measurementinstance[x].itemid,0];
return hx
def measurementvector1 (griddata,measdata,Vreal,Vimag):#this function restructures the input measuremenets in the z vector
z = numpy.zeros((measdata.nmeas,1))
for x in range(measdata.nmeas):
if measdata.measurementinstance[x].mtype == 1 or measdata.measurementinstance[x].mtype == 2: #voltage magnitude measurements converted to real part
for y in range(measdata.nmeas):
if measdata.measurementinstance[y].itemid == measdata.measurementinstance[x].itemid and (measdata.measurementinstance[y].mtype == 3 or measdata.measurementinstance[y].mtype == 4) : #and measdata.measurementinstance[y].measurementaccuracy == measdata.measurementinstance[x].measurementaccuracy:
ph1=measdata.measurementinstance[y].measurementvalue;
temp=measdata.measurementinstance[x].measurementvalue*cmath.exp(1j*ph1)
z[x,0]=temp.real;
if measdata.measurementinstance[x].mtype == 3 or measdata.measurementinstance[x].mtype == 4: #voltage phase angle measurements converted to imag part
for y in range(measdata.nmeas):
if (measdata.measurementinstance[y].itemid == measdata.measurementinstance[x].itemid) and (measdata.measurementinstance[y].mtype == 1 or measdata.measurementinstance[y].mtype == 2) : #and measdata.measurementinstance[y].measurementaccuracy == measdata.measurementinstance[x].measurementaccuracy:
mag1=measdata.measurementinstance[y].measurementvalue;
temp=mag1*cmath.exp(1j*measdata.measurementinstance[x].measurementvalue)
z[x,0]=temp.imag;
if measdata.measurementinstance[x].mtype == 5 or measdata.measurementinstance[x].mtype == 6: #current magnitude measurements converted to real part
for y in range(measdata.nmeas):
if (measdata.measurementinstance[y].itemid == measdata.measurementinstance[x].itemid) and (measdata.measurementinstance[y].mtype == 7 or measdata.measurementinstance[y].mtype == 8) : #and measdata.measurementinstance[y].measurementaccuracy == measdata.measurementinstance[x].measurementaccuracy:
ph1=measdata.measurementinstance[y].measurementvalue;
temp=measdata.measurementinstance[x].measurementvalue*cmath.exp(1j*ph1)
z[x,0]=temp.real;
if measdata.measurementinstance[x].mtype == 7 or measdata.measurementinstance[x].mtype == 8: #current phase angle measurements converted to imag part
for y in range(measdata.nmeas):
if (measdata.measurementinstance[y].itemid == measdata.measurementinstance[x].itemid) and (measdata.measurementinstance[y].mtype == 5 or measdata.measurementinstance[y].mtype == 6) : #and measdata.measurementinstance[y].measurementaccuracy == measdata.measurementinstance[x].measurementaccuracy:
mag1=measdata.measurementinstance[y].measurementvalue;
temp=mag1*cmath.exp(1j*measdata.measurementinstance[x].measurementvalue)
z[x,0]=temp.imag;
if measdata.measurementinstance[x].mtype == 9 or measdata.measurementinstance[x].mtype == 10: #active power injection
for y in range(measdata.nmeas):
if (measdata.measurementinstance[y].itemid == measdata.measurementinstance[x].itemid) and (measdata.measurementinstance[y].mtype == 11 or measdata.measurementinstance[y].mtype == 12):
Q1=measdata.measurementinstance[y].measurementvalue;
temp=+(measdata.measurementinstance[x].measurementvalue + 1j*Q1) / (Vreal[measdata.measurementinstance[x].itemid,0] + 1j*Vimag[measdata.measurementinstance[x].itemid,0])
z[x,0]=temp.real;
if measdata.measurementinstance[x].mtype == 11 or measdata.measurementinstance[x].mtype == 12: #reactive power injection
for y in range(measdata.nmeas):
if (measdata.measurementinstance[y].itemid == measdata.measurementinstance[x].itemid) and (measdata.measurementinstance[y].mtype == 9 or measdata.measurementinstance[y].mtype == 10):
P1=measdata.measurementinstance[y].measurementvalue;
temp=(P1 + 1j*measdata.measurementinstance[x].measurementvalue) / (Vreal[measdata.measurementinstance[x].itemid,0] + 1j*Vimag[measdata.measurementinstance[x].itemid,0])
z[x,0]=-temp.imag;# the current is taken as complex conjugate
if measdata.measurementinstance[x].mtype == 13 or measdata.measurementinstance[x].mtype == 14: #active power flow
in1 = griddata.line[measdata.measurementinstance[x].itemid].linefrom #the voltage is the one at the beginning of the line
for y in range(measdata.nmeas):
if (measdata.measurementinstance[y].itemid == measdata.measurementinstance[x].itemid) and (measdata.measurementinstance[y].mtype == 15 or measdata.measurementinstance[y].mtype == 16):
Q1=measdata.measurementinstance[y].measurementvalue;
temp= (measdata.measurementinstance[x].measurementvalue + 1j*Q1) / (Vreal[in1,0] + 1j*Vimag[in1,0])
z[x,0]=temp.real;
if measdata.measurementinstance[x].mtype == 15 or measdata.measurementinstance[x].mtype == 16: #reactive power flow
in1 = griddata.line[measdata.measurementinstance[x].itemid].linefrom #the voltage is the one at the beginning of the line
for y in range(measdata.nmeas):
if (measdata.measurementinstance[y].itemid == measdata.measurementinstance[x].itemid) and (measdata.measurementinstance[y].mtype == 13 or measdata.measurementinstance[y].mtype == 14):
P1=measdata.measurementinstance[y].measurementvalue;
temp=(P1 + 1j*measdata.measurementinstance[x].measurementvalue) / (Vreal[in1,0] + 1j*Vimag[in1,0])
z[x,0]=-temp.imag;
return z
def jacobianmatrixIRI1(griddata,measdata):#function to calculate the Jacobian of the system with current state
# in the single phase case we take only phase A. Each element of phase A, is equal to 1/3*(2*pos+zero)
R = numpy.zeros((1,griddata.nline)) #in the jacobian we always use the inverse of the series impedance
X = numpy.zeros((1,griddata.nline))
B = numpy.zeros((1,griddata.nline))
G = numpy.zeros((1,griddata.nline))
A = numpy.zeros((griddata.nline,griddata.nnode))
jmatrix = numpy.zeros((measdata.nmeas,2*(griddata.nnode)))
for m in range(griddata.nline): #logic matrix M lines X N nodes with A(m,n) = 1 if the node n is supplied by the line m
for n in range(griddata.nnode):
if griddata.line[m].lineto == n :
A[m,n] = 1;
elif griddata.line[m].linefrom == n :
for t in range(griddata.nline):
if A[t,griddata.line[m].linefrom] == 1 :
A[t,griddata.line[m].lineto] = 1;
for x in range(griddata.nline):
puinvbase=(griddata.pbase/((griddata.node[0].vbase)**2))**(-1);
pubase=(griddata.pbase/((griddata.node[0].vbase)**2));
Ztemp= griddata.line[x].phaseP11[0].resistance +1j*griddata.line[x].phaseP11[0].reactance
R[0,x] = (pubase*(Ztemp)).real
X[0,x] = (pubase*(Ztemp)).imag
for n in range(measdata.nmeas):
if measdata.measurementinstance[n].mtype == 1 or measdata.measurementinstance[n].mtype == 2: # voltage magnitude measurement
for index,row in enumerate(griddata.line): #we need to find the nodes connected
if A[index,measdata.measurementinstance[n].itemid] == 1:
jmatrix[n,2+2*index] = - R[0,index] #itemid refer to the node where the measurement is done
jmatrix[n,2+2*index+1] = + X[0,index] #itemid refer to the node where the measurement is done
jmatrix[n,0] = 1;
if measdata.measurementinstance[n].mtype == 3 or measdata.measurementinstance[n].mtype == 4: # voltage phase angle measurement
for index,row in enumerate(griddata.line): #we need to find the nodes connected
if A[index,measdata.measurementinstance[n].itemid] == 1:
jmatrix[n,2+2*index] = - X[0,index] #itemid refer to the node where the measurement is done
jmatrix[n,2+2*index+1] = - R[0,index] #itemid refer to the node where the measurement is done
jmatrix[n,1] = 1
measdata.