PSEN_Eficas/PSEN/PSSEWrapper.py

   1 #===============================================================================
   2 #   PSEN SCRIPT FOR PROBABILISTIC STUDIES OF ELECTICAL NETWORKS
   3 #===============================================================================
   4 from openturns import *
   5 from pylab import *
   6 from math import*
   7 import os, random
   8 import numpy as np
   9 from time import gmtime, strftime
  10 from array import *
  11 from support_functions import *
  12
  13 # Ouverture du fichier de configuration et recupecation des valeurs sous forme de listes
  14 f=open("C:\B31272\Documents\PSEN\PSENdev\lib\config.psen","r")
  15 lines=f.readlines()
  16 data_config=lines[0].split(";")
  17 data_load2 = getUserLaw(lines[1].split(";"))[0]
  18 data_wind1=getUserLaw(lines[2].split(";"))[0]
  19 data_wind2=getUserLaw(lines[3].split(";"))[0]
  20 data_corr=lines[4].split(";")
  21 data_PV=getUserLaw(lines[5].split(";"))[0]
  22 data_OPF=lines[6].split(";")
  23 Irate_num=int(lines[7].split(";")[0])
  24 f.close()
  25
  26 # Definition des variables pour les series temporelles : getUserLaw(lines[1].split(";"))[1][0] doit valoir 1
  27 # pour que le programme etudie les series temporelles
  28 time_serie_file=[]
  29 time_serie_mat=[]
  30 for i in ([1,2,3,5]) :
  31     TSoptions = getUserLaw(lines[i].split(";"))[1]
  32     if TSoptions[0] == 1 :
  33         time_serie=1
  34         f=open(TSoptions[1],"r")
  35         linesTS=f.readlines()
  36         f.close()
  37         tsm=[]
  38         for j in range (len(linesTS)) :
  39             try :
  40                 float(linesTS[j])
  41             except ValueError :
  42                 linesTS[j] = commaToPoint(linesTS[j])
  43             else :
  44                 pass
  45             tsm.append(float(linesTS[j]))
  46         time_serie_mat.append(tsm)
  47         time_serie_file.append(TSoptions[1])
  48         """time_serie_SS = TSoptions[2]
  49         time_serie_TH = TSoptions[3]"""
  50     else :
  51         time_serie_file.append(-1)
  52 time_serie_mat=zip(*time_serie_mat)
  53
  54 # Ouverture du fichier de preferences et recuperation des donnees
  55 f=open("C:\B31272\Documents\PSEN\PSENdev\lib\pref.psen","r")
  56 lines=f.readlines()
  57 f.close()
  58 paths=lines[0].split(";")
  59 WTconfig=[]
  60 for i in range (3,8):
  61     try :
  62         paths[i]
  63     except :
  64         print "Error in defining wind turbines characteristics"
  65         WTconfig=[3.,5.,25.,1.225,0.05]
  66     else :
  67         WTconfig.append(float(paths[i]))
  68
  69 # Ouverture du fichier d'analyse de N-1 et recuperation des donnees
  70 f=open("C:\B31272\Documents\PSEN\PSENdev\lib\contin.psen","r")
  71 lines=f.readlines()
  72 path_config_contin=lines[0].split(";")
  73 f.close()
  74
  75 # Definition des lois des cinq variables aleatoires
  76 N_1=int(data_config[1]) # N_1=1 to do N-1 studies
  77 wind_var1=int(data_config[3]) # To take wind1 variability into account
  78 wind_var2=int(data_config[4]) # To take wind2 variability into account
  79 PV_var=int(data_config[2]) # To take PV variability into account
  80 load_var=int(data_config[5]) # To take load variability into account
  81
  82 # Creation variable nom dossier N-1
  83 if N_1 == 1 :
  84     folderN_1 = '1_'
  85 else :
  86     folderN_1 = '_'
  87
  88 # Recuperation des chemins du dossier d'installation de PSSE, .SAV de l'etude et nom du rapport
  89 folder=paths[1]
  90 doc_base= paths[0]
  91 exec_file="report.txt"
  92
  93 # Definition des groupes de production PV, eoliennes, des intercos et des lignes en N-1
  94 ENR=config_ENR(paths[2])
  95 windTurbines1 = ENR[1] # Buses with wind turbines 1
  96 windTurbines2 = ENR[2] # Buses with wind turbines 2
  97 solarPV = ENR[0] # Buses with solar PV plant
  98 intercos=ENR[3] # Buses with interconnexions
  99 # Lines with contingency
 100 try :
 101     config_contingency(path_config_contin)
 102 except IOError :  # Si le fichier n'est pas dans un bon format on traite l'exception
 103     nb_lines=1
 104     print 'Error with contingency input file'
 105 else :
 106     continAll = config_contingency(path_config_contin)
 107     continLines = continAll[0]
 108     continGroups = continAll[1]
 109     continVal = continAll[2]
 110     continProb = continAll[3]
 111
 112
 113 # Probabilistic study information
 114 #==============================================================================
 115 # Create the marginal distributions
 116 distributionX0 = data_load2
 117 distributionX1 = getUserDefined(continVal,continProb)
 118 distributionX2 = data_wind1
 119 distributionX3 = data_wind2
 120 distributionX4 = data_PV
 121
 122 # Create the correlations between the distributions
 123 corr10=float(data_corr[0])
 124 corr20=float(data_corr[1])
 125 corr30=float(data_corr[2])
 126 corr40=float(data_corr[3])
 127 corr21=float(data_corr[4])
 128 corr31=float(data_corr[5])
 129 corr41=float(data_corr[6])
 130 corr32=float(data_corr[7])
 131 corr42=float(data_corr[8])
 132 corr43=float(data_corr[9])
 133
 134 # Probabilistic Study: central dispersion => Monte Carlo or LHS iterations
 135 montecarlosize = int(data_config[0])
 136
 137 #Extension name for the folders and files
 138 day=time.strftime("%Y%m%d", gmtime())
 139 hour=time.strftime("%Hh%Mm%S", gmtime())
 140
 141 #===============================================================================
 142 #    CHARGEMENT DE PSSE     -   LOADING OF PSSE
 143 #===============================================================================
 144 pssFolder=str(paths[3])
 145 import sys
 146 sys.path.append(pssFolder)#r"C:\Program Files\PTI\PSSE33\PSSBIN")
 147 os.environ['PATH'] = pssFolder+":"+os.environ['PATH'] #r"C:\Program Files\PTI\PSSE33\PSSBIN;"+ os.environ['PATH']
 148 os.chdir(folder)
 149 import psspy
 150 import pssarrays
 151 import redirect
 152 _i=psspy.getdefaultint()
 153 _f=psspy.getdefaultreal()
 154 _s=psspy.getdefaultchar()
 155 redirect.psse2py()
 156 #import pssdb
 157 psspy.psseinit(80000)
 158
 159 # Silent execution of PSSe
 160 islct=6 # 6=no output; 1=standard
 161 psspy.progress_output(islct)
 162
 163 # Enregistrement de l'heure de debut de simulation
 164 f=open(exec_file, 'a')
 165 start_time=time.clock()
 166 f.write("Starting time: %f;     Monte Carlo Size : %f;      " % (start_time, montecarlosize))
 167 f.close()
 168
 169 #===============================================================================
 170 #    Fonction de wrappage     -   Wrapper function
 171 #===============================================================================
 172 def PSSEFunction(x):
 173     # Definition des variables globales
 174     global TStest
 175     global Xt
 176     global sizeY0
 177     global sizeY1
 178     global sizeY2
 179     global sizeY3
 180     global sizeY4
 181     global sizeY
 182     global wind_var
 183     global PV_var
 184     global N_1
 185     global load_var
 186     global logCSVfilename
 187     global logTXTfilename
 188     global ite
 189     global folder
 190     global day
 191     global folderN_1
 192     global fich
 193     global hour
 194     global montecarlosize
 195     global WTconfig
 196     global x2
 197
 198     ite+=1 # incrementation du compteur
 199
 200     # Load data from PSSe
 201     psspy.case(doc_base) #Launching of PSSE and opening the working file
 202     all_inputs_base=read_sav(doc_base)
 203     buses_base=all_inputs_base[0]
 204     lines_base=all_inputs_base[1]
 205     transf_base=all_inputs_base[2]
 206     plants_base=all_inputs_base[3]
 207     loads_base=all_inputs_base[4]
 208     shunt_base=all_inputs_base[5]
 209     doci=folder+"\N"+folderN_1+day+"\CasNum"+str(ite)+".sav"
 210     psspy.save(doci)
 211
 212     # Total initial shunt on buses
 213     init_shunt = 0
 214     for i in range(len(shunt_base)) :
 215         init_shunt +=  float(shunt_base[i][2])
 216
 217     # Configuration de l'OPF a partir des parametres de l'utilisateur
 218     nbeOPF=5 # Nombre de lancement max de l'OPF pour atteindre la convergence de l'algorithme
 219     psspy.report_output(6,"",[0,0])
 220     psspy.produce_opf_log_file(1,r"""DETAIL""")
 221     psspy.minimize_fuel_cost(int(data_OPF[0]))
 222     psspy.minimize_adj_bus_shunts(int(data_OPF[1]))
 223     psspy.minimize_load_adjustments(int(data_OPF[2]))
 224     psspy.initial_opf_barrier_coeff(100.0)
 225     psspy.opf_fix_all_generators(1)
 226     psspy.set_opf_report_subsystem(3,1)
 227
 228
 229     print "                     PSEN simulator, case number: "+str(ite)
 230
 231     # 1. Affiche X
 232     nx = x.getSize()
 233     if TStest==1 :
 234         for i in range (len (Xt)) :
 235             if Xt[i] == -1 :
 236                 if i == 0 and load_var==1 :
 237                     pass
 238                 elif i == 1 and  N_1==1 :
 239                     x[i]=int(floor(x[i])) # Si on etudie le N-1 on arrondie la valeur tiree en floor : on obtient le numero de la ligne
 240                 elif i == 2 and wind_var1==1 :
 241                    x[i]=eol(x[i],WTconfig)
 242                 elif i == 3 and wind_var2==1 :
 243                    x[i]=eol(x[i],WTconfig)
 244                 elif i == 4 and PV_var==1 :  # Si on etudie la variabilite de l'eolien on lui donne la valeur de production de l'eolienne a partir du vent
 245                    pass
 246                 else :
 247                    x[i]=-1
 248             else :
 249                 x[i]=float(Xt[i]) # Dans le cas d'une etude temporelle on lui donne la valeur de Xt
 250     else :
 251         if load_var==1 :
 252             pass # Sinon on donne la valeur tiree si on etudie la variabilite de x[0]
 253         else :
 254             x[0]=1 # Sinon on laisse la valeur de base
 255
 256         if N_1==1 :
 257             x[1]=int(floor(x[1])) # Si on etudie le N-1 on arrondie la valeur tiree en floor : on obtient le numero de la ligne
 258         else :
 259             x[1]=-1 # Sinon on donne -1 comme marqueur
 260
 261         if wind_var1==1:
 262             x[2]=eol(x[2],WTconfig) # Si on etudie la variabilite de l'eolien on lui donne la valeur de production de l'eolienne a partir du vent
 263         else :
 264             x[2]=0 # Sinon on considere qu'il n'y a pas d'eolien
 265
 266         if wind_var2==1:
 267             x[3]=eol(x[3],WTconfig) # Si on etudie la variabilite de l'eolien on lui donne la valeur de production de l'eolienne a partir du vent
 268         else :
 269             x[3]=0 # Sinon on considere qu'il n'y a pas d'eolien
 270
 271         if PV_var==1 : # Si on etudie la variabilite du PV on laisse sa valeur a la va
 272             pass
 273         else :
 274             x[4]=0 # Sinon on considere qu'il n'y a pas de PV
 275     for i in range(0,nx):
 276         print "x[%d]=%f" % (i,x[i])
 277
 278     # 2. Fait le calcul avec PSSE
 279     #Editing some values in the PSSE .sav input file
 280     # Change the values of the different loads and treat large changes of load to help convergence
 281     if x[0] > 0.75 :
 282         for i in range(0,np.array(loads_base).shape[0]) : # On change directement toutes les charges
 283             psspy.load_chng_4(int(loads_base[i][0]),r"""1""",[1,_i,_i,_i,_i,_i],[x[0]*loads_base[i][1],x[0]*loads_base[i][2],_f,_f,_f,_f])
 284     elif x[0] > 0.4 :
 285         for i in range(0,np.array(loads_base).shape[0]) :  # On effectue un pretraitement en passant par une charge intermediaire
 286             psspy.load_chng_4(int(loads_base[i][0]),r"""1""",[1,_i,_i,_i,_i,_i],[(1+x[0])/2*loads_base[i][1],(1+x[0])/2*loads_base[i][2],_f,_f,_f,_f])
 287         psspy.fnsl([0,0,0,1,1,0,99,0]) # Load flow Newton Raphson
 288         psspy.bsys(3,0,[0.0,0.0],0,[],1,[1],0,[],0,[])
 289         psspy.set_opf_report_subsystem(3,0)
 290         psspy.nopf(0,1) # Lancement OPF
 291         for i in range(0,np.array(loads_base).shape[0]) : # On change toutes les charges
 292             psspy.load_chng_4(int(loads_base[i][0]),r"""1""",[1,_i,_i,_i,_i,_i],[x[0]*loads_base[i][1],x[0]*loads_base[i][2],_f,_f,_f,_f])
 293     else :
 294         for i in range(0,np.array(loads_base).shape[0]) : # On effectue un pretraitement en passant par une charge intermediaire
 295             psspy.load_chng_4(int(loads_base[i][0]),r"""1""",[1,_i,_i,_i,_i,_i],[0.7*loads_base[i][1],0.7*loads_base[i][2],_f,_f,_f,_f])
 296         psspy.fnsl([0,0,0,1,1,0,99,0])
 297         psspy.bsys(3,0,[0.0,0.0],0,[],1,[1],0,[],0,[])
 298         psspy.set_opf_report_subsystem(3,0)
 299         psspy.nopf(0,1)
 300         for i in range(0,np.array(loads_base).shape[0]) : # On effectue un pretraitement en passant par une charge intermediaire
 301             psspy.load_chng_4(int(loads_base[i][0]),r"""1""",[1,_i,_i,_i,_i,_i],[0.4*loads_base[i][1],0.4*loads_base[i][2],_f,_f,_f,_f])
 302         psspy.fnsl([0,0,0,1,1,0,99,0])
 303         psspy.bsys(3,0,[0.0,0.0],0,[],1,[1],0,[],0,[])
 304         psspy.set_opf_report_subsystem(3,0)
 305         psspy.nopf(0,1)
 306         for i in range(0,np.array(loads_base).shape[0]) : # On change toutes les charges
 307             psspy.load_chng_4(int(loads_base[i][0]),r"""1""",[1,_i,_i,_i,_i,_i],[x[0]*loads_base[i][1],x[0]*loads_base[i][2],_f,_f,_f,_f])
 308
 309     x2=[]
 310     for sz in range(0,nx):
 311         x2.append(float(x[sz]))
 312
 313     if x[1]<0 :
 314         pass
 315     elif x[1] < len(continLines) : # L'element tire est une ligne
 316         line_num=int(x[1])
 317         from_bus=continLines[int(line_num)][0]
 318         to_bus=continLines[int(line_num)][1]
 319         br_id=continLines[int(line_num)][2]#.replace('@','')
 320         psspy.branch_chng(from_bus,to_bus,str(br_id),[0,_i,_i,_i,_i,_i],[ _f, _f, _f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f])
 321         x2[1]='Line '+str(from_bus)+'-'+str(to_bus)+'#'+str(br_id)
 322     elif x[1] < (len(continLines)+len(continGroups)) :
 323         group_num = int(x[1])-len(continLines)
 324         bus_num = continGroups[int(group_num)][0]
 325         bus_id = continGroups[int(group_num)][1]
 326         psspy.machine_chng_2(int(bus_num),str(bus_id),[0,_i,_i,_i,_i,_i],[_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f]) # Change Interconnection disponibility
 327         psspy.opf_gendsp_indv(int(bus_num),str(bus_id),_i,0.0)
 328         x2[1]='Group '+str(bus_num)+'#'+str(bus_id)
 329     #elif x[1] < len(intercos) :
 330         #mat_num=int(x[1])
 331         #psspy.machine_chng_2(int(intercos[mat_num][0]),str(intercos[mat_num][2]),[0,_i,_i,_i,_i,_i],[_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f]) # Change Interconnection disponibility
 332         #psspy.opf_gendsp_indv(int(intercos[mat_num][0]),str(intercos[mat_num][2]),_i,0.0)
 333         #x[1]=-mat_num
 334     else :
 335         pass
 336         # Change the bus that is not in service
 337         #intercos = []
 338         #line_num=int(x[1]-len(intercos))
 339         #from_bus=lines_con[int(line_num)-1][0]
 340         #to_bus=lines_con[int(line_num)-1][1]
 341         #br_id=lines_con[int(line_num)-1][2]#.replace('@','')
 342         #psspy.branch_chng(from_bus,to_bus,str(br_id),[0,_i,_i,_i,_i,_i],[ _f, _f, _f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f])
 343         #x[1]=line_num
 344     # Change the production of the wind turbines
 345     if np.matrix(windTurbines1).shape[1]>0 :
 346         for i in range(0,np.matrix(windTurbines1).shape[0]) :
 347             psspy.machine_chng_2(windTurbines1[i][0],str(windTurbines1[i][2]),[1,_i,_i,_i,_i,_i],[x[2]*plants_base[windTurbines1[i][1]][6],_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f])
 348
 349     if np.matrix(windTurbines2).shape[1]>0 :
 350         for i in range(0,np.matrix(windTurbines2).shape[0]) :
 351             psspy.machine_chng_2(windTurbines2[i][0],str(windTurbines2[i][2]),[1,_i,_i,_i,_i,_i],[x[3]*plants_base[windTurbines2[i][1]][6],_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f])
 352
 353     # Change the production of the PV stations
 354     if np.matrix(solarPV).shape[1]>0 :
 355         for i in range(0,np.matrix(solarPV).shape[0]) :
 356             psspy.machine_chng_2(solarPV[i][0],str(solarPV[i][2]),[1,_i,_i,_i,_i,_i],[x[4]*plants_base[solarPV[i][1]][6],_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f,_f])
 357
 358     psspy.save(doci) #Saving .sav modifications
 359     ok=1
 360     while nbeOPF>=0 :
 361         #for i in (zip(*buses_base)[0]) : psspy.bus_chng_3(i,[_i,_i,_i,_i],[_f, 1.05,_f,_f,_f,_f,_f],_s)
 362         psspy.fnsl([0,0,0,1,1,1,99,0])
 363         psspy.bsys(3,0,[0.0,0.0],0,[],1,[1],0,[],0,[])
 364         psspy.set_opf_report_subsystem(3,0)
 365         psspy.nopf(0,1)
 366         if psspy.solved()==7:
 367             print 'CONVERGENCE          CAS '+str(ite)
 368             ok=1
 369             break
 370         else :
 371             print '==================================================================='
 372             print 'NO CONVERGENCE'
 373             print '==================================================================='
 374             ok=0
 375             #for i in range (134) :
 376                 #psspy.opf_bus_indv(i,[_i,0],[_f, 0.7,_f,_f,_f])
 377         nbeOPF-=1
 378     psspy.save(doci)
 379     all_inputs=read_sav(doci)
 380     buses=all_inputs[0];lines=all_inputs[1];transf=all_inputs[2];plants=all_inputs[3];loads=all_inputs[4]; shunt=all_inputs_base[5]
 381
 382     # 3. Affiche Y
 383     sizeY4=np.matrix(shunt).shape[0]
 384     y=np.zeros(2*sizeY0+sizeY1+3*sizeY2+sizeY3+sizeY4)
 385     z=np.zeros(8)
 386     rate_mat_index=Irate_num+2
 387     if ok==1 :
 388         # Creates the quantities of interest
 389         for i in range (sizeY2) :
 390             if lines [i][rate_mat_index]>100 :
 391                 z[0]+=1 # Number of lines above 100% of their limits
 392         for i in range (sizeY1):
 393             if buses[i][2]>1.06 :
 394                 z[1]+=1
 395             if buses[i][2]<0.9399 :
 396                 z[1]+=1 # Number of buses outside of their voltage limits
 397         for i in range (sizeY0) :
 398             z[2]+=float(plants[i][3]) # Total active production
 399         for i in range (sizeY3) :
 400             z[3]+=float(loads[i][1]) # Total active consumption
 401         z[4]=(z[2]-z[3])/z[2]*100 # Active power losses
 402         for i in range (sizeY2) :
 403             if lines [i][3]>z[5] :
 404                 z[5]=lines [i][rate_mat_index] # Max flow in lines
 405         for i in range (sizeY2) :
 406             if lines [i][rate_mat_index]>90 :
 407                 z[6]+=1
 408         z[6]=z[6]-z[0] # Number of lines between 90% and 100% of their limits
 409
 410         final_shunt=0
 411         for i in range (sizeY4) :
 412             final_shunt+=shunt[i][2]
 413         z[7]=final_shunt-init_shunt
 414
 415        # Creates the output vectors
 416         for Pmach in range (sizeY0):
 417             y[Pmach]=float(plants[Pmach][3])
 418         for Qmach in range (sizeY0):
 419             y[Qmach+sizeY0]=float(plants[Qmach][4])
 420         for Vbus in range (sizeY1):
 421             y[Vbus+2*sizeY0]=float(buses[Vbus][2])
 422         for Iline in range (sizeY2):
 423             y[Iline+2*sizeY0+sizeY1]=float(lines[Iline][rate_mat_index])
 424         for Pline in range (sizeY2):
 425             y[Pline+2*sizeY0+sizeY1+sizeY2]=float(lines[Pline][6])
 426         for Qline in range (sizeY2):
 427             y[Qline+2*sizeY0+sizeY1+2*sizeY2]=float(lines[Qline][7])
 428         for Pload in range (sizeY3) :
 429             y[Pload+2*sizeY0+sizeY1+3*sizeY2]=float(loads[Pload][1])
 430         for Qshunt in range (sizeY4) :
 431             y[Qshunt+2*sizeY0+sizeY1+3*sizeY2+sizeY3]=float(shunt[Qshunt][2])
 432
 433         #Ecris les sorties
 434         print "sorties:"
 435         nz = len(z)
 436         for i in range(0,nz):
 437             print "z[%d]=%f" % (i,z[i])
 438         MyLogger(x2,y,z,logCSVfilename,logTXTfilename,ite)
 439         #MyMultiLogger (x2, y, sizeY, z, ite, folder, day, fich, hour)
 440         return NumericalPoint(z)
 441     else :
 442         MyLogger(x2,y,z,logCSVfilename,logTXTfilename,ite)
 443         #MyMultiLogger (x2, y, sizeY, z, ite, folder, day, fich, hour)
 444         return NumericalPoint(z)
 445
 446 #===============================================================================
 447 #   DEFINITION DU WRAPPER -  WRAPPER's DEFINITION
 448 #===============================================================================
 449 # Initialize size output
 450 psspy.case(doc_base)
 451 all_inputs_base=read_sav(doc_base)
 452 buses_base=all_inputs_base[0]
 453 lines_base=all_inputs_base[1]
 454 trans_base=all_inputs_base[2]
 455 plants_base=all_inputs_base[3]
 456 loads_base=all_inputs_base[4]
 457 shunt_base=all_inputs_base[5]
 458 sizeY0=np.matrix(plants_base).shape[0]
 459 sizeY1=np.matrix(buses_base).shape[0]
 460 sizeY2=np.matrix(lines_base).shape[0]
 461 sizeY3=np.matrix(loads_base).shape[0]
 462 sizeY4=np.matrix(shunt_base).shape[0]
 463 sizeY=[sizeY0,sizeY1,sizeY2,sizeY3,sizeY4]
 464 sizeOutput=sizeY2
 465
 466
 467 class PSSEWrapperClass(OpenTURNSPythonFunction) :
 468   def __init__(self) :
 469      OpenTURNSPythonFunction.__init__(self,5,8)
 470   def _exec(self,x) :
 471       return PSSEFunction(x)
 472
 473 # Initialize the folder
 474 newpath = folder+"\N"+folderN_1+day
 475 if not os.path.exists(newpath): os.makedirs(newpath)
 476
 477 # Test the Num. Math. Function
 478 pssefun = NumericalMathFunction(PSSEWrapperClass())
 479
 480 # Definition of the function to use
 481 inputDim = pssefun.getInputDimension()
 482 outputDim = pssefun.getOutputDimension()
 483
 484 # Initialization of the distribution collection:
 485 #aCollection = DistributionCollection()
 486
 487 # Create a collection of the marginal distributions
 488 collectionMarginals = DistributionCollection(inputDim)
 489 collectionMarginals[0] = Distribution(distributionX0) # Load distribution
 490 collectionMarginals[1] = Distribution(distributionX1) # N-1 distribution
 491 collectionMarginals[2] = Distribution(distributionX2) # Wind 1 distribution
 492 collectionMarginals[3] = Distribution(distributionX3) # Wind 2 distribution
 493 collectionMarginals[4] = Distribution(distributionX4) # PV distribution
 494
 495 #Create a correlation matrix as copulas
 496 corr=CorrelationMatrix(inputDim)
 497
 498 corr[1,0]=corr10
 499 corr[2,0]=corr20
 500 corr[3,0]=corr30
 501 corr[4,0]=corr40
 502 corr[0,1]=corr10
 503 corr[2,1]=corr21
 504 corr[3,1]=corr31
 505 corr[4,1]=corr41
 506 corr[0,2]=corr20
 507 corr[1,2]=corr21
 508 corr[3,2]=corr32
 509 corr[4,2]=corr42
 510 corr[0,3]=corr30
 511 corr[1,3]=corr31
 512 corr[2,3]=corr32
 513 corr[4,3]=corr43
 514 corr[0,4]=corr40
 515 corr[1,4]=corr41
 516 corr[2,4]=corr42
 517 corr[3,4]=corr43
 518
 519 copula=Copula(NormalCopula(corr))
 520
 521
 522 # Create the input probability distribution, args are the distributions, the correlation laws
 523 inputDistribution = ComposedDistribution(collectionMarginals, copula)
 524
 525 # Create the input random vector
 526 """inputRandomVector = RandomVector(inputDistribution)
 527
 528 # Create the output variable of interest
 529 outputVariableOfInterest =  RandomVector(pssefun, inputRandomVector)
 530 outputVariableOfInterest.setDescription(pssefun.getOutputDescription())"""
 531
 532 #===============================================================================
 533 #   ETUDE DE DISPERSION CENTRALE    -   CENTRAL DEVIATION STUDY
 534 #===============================================================================
 535 # Initialize the logger : write the headers
 536 logCSVfilename=folder+"\N"+folderN_1+day+"\simulationDClog"+hour+".csv" # Name of the file : global variable
 537 f = open(logCSVfilename, "a")
 538 f.write("Iteration;;X:Load(pu);X:lineOff#;XProdEolienne1%Pnom;XProdEolienne2%Pnom;X:ProdPV%Pnom;;Y:NbeTransit;Y:NbeTension;Y:PProdTot;Y:PConsoTot;Y:%Losses;Y:Max%A;Y:NbeTransit_0.9-1;Y:AddedMVAR;;")
 539 # Names of the Output variables withConso the bus number
 540 for name in range (sizeY0):
 541     f.write("Y:PMachine"+str(plants_base[name][0])+";")
 542 for name in range (sizeY0):
 543     f.write("Y:QMachine"+str(plants_base[name][0])+";")
 544 for name in range (sizeY1):
 545     f.write("Y:VBus"+str(buses_base[name][0])+";")
 546 for name in range (sizeY2):
 547     f.write("Y"+str(name+1)+":%Rate "+str(lines_base[name][0])+"-"+str(lines_base[name][1])+";")
 548 for name in range (sizeY2):
 549     f.write("Y"+str(name+1)+":P "+str(lines_base[name][0])+"-"+str(lines_base[name][1])+";")
 550 for name in range (sizeY2):
 551     f.write("Y"+str(name+1)+":Q "+str(lines_base[name][0])+"-"+str(lines_base[name][1])+";")
 552 for name in range (sizeY3):
 553     f.write("Y:Load "+str(loads_base[name][0])+";")
 554 for name in range (sizeY4):
 555     f.write("Y:Shunt bus "+str(shunt_base[name][0])+";")
 556 f.write("\n")
 557 # Names of the Output variables with the bus names
 558 f.write("#;;pu;line number;%Pnom;%Pnom;%Pnom;;Nbe;Nbe;MW;MW;%;%;Nbe;MVAR;;")
 559 for name in range (sizeY0):
 560     f.write(str(plants_base[name][8])+";")
 561 for name in range (sizeY0):
 562     f.write(str(plants_base[name][8])+";")
 563 for name in range (sizeY1):
 564     f.write(str(buses_base[name][3])+";")
 565 for name in range (sizeY2):
 566     f.write(str(lines_base[name][8])+"-"+str(lines_base[name][9])+";")
 567 for name in range (sizeY2):
 568     f.write(str(lines_base[name][8])+"-"+str(lines_base[name][9])+";")
 569 for name in range (sizeY2):
 570     f.write(str(lines_base[name][8])+"-"+str(lines_base[name][9])+";")
 571 for name in range (sizeY3):
 572     f.write(str(loads_base[name][4])+";")
 573 for name in range (sizeY4):
 574     f.write(str(shunt_base[name][3])+";")
 575 f.write("\n")
 576 f.close()
 577
 578 logTXTfilename=folder+"\N"+folderN_1+day+"\simulationDClog"+hour+".txt" # Name of the file : global variable
 579 f = open(logTXTfilename, "a")
 580 f.write("Iteration\tX:Load(pu)\tX:lineOff#\tXProdEolienne1%Pnom\ttXProdEolienne2%Pnom\tX:ProdPV%Pnom\tY:NbeTransit\tY:NbeTension\tY:PProdTot\tY:PConsoTot\tY:%Losses\tY:Max%A\tY:NbeTransit_0.9-1\tY:AddedShunt\t")
 581 # Names of the Output variables withConso the bus number
 582 for name in range (sizeY0):
 583     f.write("Y:PMachine"+str(plants_base[name][0])+" - "+str(plants_base[name][8])+"\t")
 584 for name in range (sizeY0):
 585     f.write("Y:QMachine"+str(plants_base[name][0])+" - "+str(plants_base[name][8])+"\t")
 586 for name in range (sizeY1):
 587     f.write("Y:VBus"+str(buses_base[name][0])+" - "+str(buses_base[name][3])+"\t")
 588 for name in range (sizeY2):
 589     f.write("Y"+str(name+1)+":%RateA "+str(lines_base[name][0])+"-"+str(lines_base[name][1])+" - "+str(lines_base[name][8])+"-"+str(lines_base[name][9])+"\t")
 590 for name in range (sizeY2):
 591     f.write("Y"+str(name+1)+":P "+str(lines_base[name][0])+"-"+str(lines_base[name][1])+" - "+str(lines_base[name][8])+"-"+str(lines_base[name][9])+"\t")
 592 for name in range (sizeY2):
 593     f.write("Y"+str(name+1)+":Q "+str(lines_base[name][0])+"-"+str(lines_base[name][1])+" - "+str(lines_base[name][8])+"-"+str(lines_base[name][9])+"\t")
 594 for name in range (sizeY3):
 595     f.write("Y:Load "+str(loads_base[name][0])+" - "+str(loads_base[name][4])+"\t")
 596 for name in range (sizeY4):
 597     f.write("Y:Shunt "+str(shunt_base[name][0])+" - "+str(shunt_base[name][3])+"\t")
 598 f.write("\n")
 599 f.close()
 600
 601 """
 602 # Initialize the multilogger : write the headers
 603 for fich in range (np.size(sizeY,0)):
 604     multilogfilename=folder+"\N"+day+"\Y"+str(fich)+"simulationDClog"+hour+".csv"
 605     f=open(multilogfilename, 'a')
 606     f.write("Iteration;;X:Load(pu);X:lineOff#;XProdEolienne1%Pnom;XProdEolienne2%Pnom;X:ProdPV%Pnom;;Y:NbeTransit;Y:NbeTension;Y:PProdTot;Y:PConsoTot;Y:Max%A;Y:NbeTransit_0.9-1;;")
 607     if fich == 0 :
 608         for name in range (sizeY[0]):
 609             f.write("Y:PMachine"+str(plants_base[name][0])+";")
 610         f.write("\n")
 611         f.write("#;;pu;line number;%Pnom;%Pnom;%Pnom;;Nbe;Nbe;MW;MW;%;%;Nbe;;")
 612         for name in range (sizeY[0]):
 613             f.write(str(plants_base[name][8])+";")
 614         f.write("\n")
 615         f.close()
 616     elif fich == 1 :
 617         for name in range (sizeY[1]):
 618             f.write("Y:VBus"+str(buses_base[name][0])+";")
 619         f.write("\n")
 620         f.write("#;;pu;line number;%Pnom;%Pnom;%Pnom;;Nbe;Nbe;MW;MW;%;%;Nbe;;")
 621         for name in range (sizeY[1]):
 622             f.write(str(buses_base[name][3])+";")
 623         f.write("\n")
 624         f.close()
 625     elif fich == 2 :
 626         for name in range (sizeY[2]):
 627             f.write("Y"+str(name+1)+":%RateA "+str(lines_base[name][0])+"-"+str(lines_base[name][1])+";")
 628         f.write("\n")
 629         f.write("#;;pu;line number;%Pnom;%Pnom;%Pnom;;Nbe;Nbe;MW;MW;%;%;Nbe;;")
 630         for name in range (sizeY[2]):
 631             f.write(str(lines_base[name][8])+"-"+str(lines_base[name][9])+";")
 632         f.write("\n")
 633         f.close()
 634     elif fich == 3 :
 635         for name in range (sizeY[3]):
 636             f.write("Y:Ploads "+str(loads_base[name][0])+";")
 637         f.write("\n")
 638         f.write("#;;pu;line number;%Pnom;%Pnom;%Pnom;;Nbe;Nbe;MW;MW;%;%;Nbe;;")
 639         for name in range (sizeY[3]):
 640             f.write(str(loads_base[name][4])+";")
 641         f.write("\n")
 642         f.close()
 643
 644 """
 645 # Start the simulations
 646 ite=0
 647 print "\n\n\n                     Starting PSEN "+str(montecarlosize)+" simulations"
 648
 649 """inputSample=inputRandomVector.getSample(montecarlosize)
 650 inputSample.setDescription( ("X0","X1","X2","X3") )
 651 inputSample.exportToCSVFile("InputSamples.csv")"""
 652
 653 if sum(corr) == 5 :
 654     myLHSE = LHSExperiment(inputDistribution,montecarlosize)
 655     inputSample = myLHSE.generate()
 656 else :
 657     myMCE = MonteCarloExperiment(inputDistribution,montecarlosize)
 658     inputSample = myMCE.generate()
 659
 660 try :
 661     time_serie
 662 except NameError :
 663     print 'Probabilistic'
 664     TStest=0
 665     outputSampleAll = pssefun(inputSample)#outputVariableOfInterest.getSample(montecarlosize)
 666 else :
 667     TStest=1
 668     for i in range (len(time_serie_mat)) :
 669         print 'Time serie'
 670         RandomGenerator.SetSeed(i)
 671         Xt=[]
 672         n=0
 673         for j in range (len(time_serie_file)) :
 674             if time_serie_file[j] == -1 :
 675                 Xt.append(-1)
 676                 n+=1
 677             else :
 678                 Xt.append(time_serie_mat[i][j-n])
 679         Xt.insert(1,-1)
 680         try :
 681             outputSampleAll
 682         except :
 683             outputSampleAll = pssefun(inputSample)
 684         else :
 685             outputSampleAll.add(pssefun(inputSample))
 686
 687 outputDim=outputSampleAll.getDimension()
 688 outputSize=outputSampleAll.getSize()
 689
 690 outputSample=NumericalSample(0,outputDim)
 691 outputSampleMissed=NumericalSample(0,outputDim)
 692
 693 for i in range (outputSize):
 694     if outputSampleAll[i,5]==0 :
 695         outputSampleMissed.add(outputSampleAll[i])
 696     else :
 697         outputSample.add(outputSampleAll[i])
 698
 699 outputDescription=[]
 700 for i in range (outputDim):
 701     outputDescription.append("Y"+str(i))
 702 outputSample.setDescription( outputDescription )
 703
 704 # Get the empirical mean and standard deviations
 705 empMeanX = inputSample.computeMean()
 706 empSdX = inputSample.computeStandardDeviationPerComponent()
 707 empiricalMean = outputSample.computeMean()
 708 empiricalSd = outputSample.computeStandardDeviationPerComponent()
 709
 710 f=open(logCSVfilename, 'a')
 711 f.write("\n")
 712 f.write('Mean;;')
 713 for i in range(0,inputDim):
 714     f.write("%f;" % (empMeanX[i]))
 715 f.write(";")
 716 for i in range(0,outputDim):
 717     f.write("%f;" % (empiricalMean[i]))
 718 f.write(";")
 719 f.write("\nStandard deviation;;")
 720 for i in range(0,inputDim):
 721     f.write("%f;" % (empSdX[i]))
 722 f.write(";")
 723 for i in range(0,outputDim):
 724     f.write("%f;" % (empiricalSd[i]))
 725 f.write(";")
 726 f.close()
 727
 728 f=open(exec_file,'a')
 729 #stop_time=100*times()[0]
 730 stop_time=time.clock()
 731 f.write("Stop time: %f;     Duration: %f;      Time per execution: %f; " % (stop_time, stop_time-start_time, (stop_time-start_time)/montecarlosize))
 732 f.write("\n\n")
 733 f.close()
 734
 735 print '\n\nSimulated '+str(montecarlosize)+' cases in '+ str(stop_time-start_time)+' seconds. Average '+str((stop_time-start_time)/montecarlosize)+'s per case.'
 736
 737 nMissed=int(outputSampleMissed.getSize())
 738
 739 print '\n\n             Non-convergence rate is '+str(round(nMissed*100/montecarlosize,3))+' % ('+str(outputSampleMissed.getSize())+' cases on '+str(montecarlosize)+')'
 740
 741 #graphical_out(inputSample, outputSampleAll, inputDim, outputDim, montecarlosize)
 742