1 #===============================================================================
2 # PSEN SCRIPT FOR PROBABILISTIC STUDIES OF ELECTICAL NETWORKS
3 #===============================================================================
4 from openturns import *
9 from time import gmtime, strftime
11 from support_functions import *
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")
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])
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
30 for i in ([1,2,3,5]) :
31 TSoptions = getUserLaw(lines[i].split(";"))[1]
32 if TSoptions[0] == 1 :
34 f=open(TSoptions[1],"r")
38 for j in range (len(linesTS)) :
42 linesTS[j] = commaToPoint(linesTS[j])
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]"""
51 time_serie_file.append(-1)
52 time_serie_mat=zip(*time_serie_mat)
54 # Ouverture du fichier de preferences et recuperation des donnees
55 f=open("C:\B31272\Documents\PSEN\PSENdev\lib\pref.psen","r")
58 paths=lines[0].split(";")
64 print "Error in defining wind turbines characteristics"
65 WTconfig=[3.,5.,25.,1.225,0.05]
67 WTconfig.append(float(paths[i]))
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")
72 path_config_contin=lines[0].split(";")
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
82 # Creation variable nom dossier N-1
88 # Recuperation des chemins du dossier d'installation de PSSE, .SAV de l'etude et nom du rapport
91 exec_file="report.txt"
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
101 config_contingency(path_config_contin)
102 except IOError : # Si le fichier n'est pas dans un bon format on traite l'exception
104 print 'Error with contingency input file'
106 continAll = config_contingency(path_config_contin)
107 continLines = continAll[0]
108 continGroups = continAll[1]
109 continVal = continAll[2]
110 continProb = continAll[3]
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
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])
134 # Probabilistic Study: central dispersion => Monte Carlo or LHS iterations
135 montecarlosize = int(data_config[0])
137 #Extension name for the folders and files
138 day=time.strftime("%Y%m%d", gmtime())
139 hour=time.strftime("%Hh%Mm%S", gmtime())
141 #===============================================================================
142 # CHARGEMENT DE PSSE - LOADING OF PSSE
143 #===============================================================================
144 pssFolder=str(paths[3])
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']
152 _i=psspy.getdefaultint()
153 _f=psspy.getdefaultreal()
154 _s=psspy.getdefaultchar()
157 psspy.psseinit(80000)
159 # Silent execution of PSSe
160 islct=6 # 6=no output; 1=standard
161 psspy.progress_output(islct)
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))
169 #===============================================================================
170 # Fonction de wrappage - Wrapper function
171 #===============================================================================
173 # Definition des variables globales
186 global logCSVfilename
187 global logTXTfilename
194 global montecarlosize
198 ite+=1 # incrementation du compteur
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"
212 # Total initial shunt on buses
214 for i in range(len(shunt_base)) :
215 init_shunt += float(shunt_base[i][2])
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)
229 print " PSEN simulator, case number: "+str(ite)
234 for i in range (len (Xt)) :
236 if i == 0 and load_var==1 :
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
249 x[i]=float(Xt[i]) # Dans le cas d'une etude temporelle on lui donne la valeur de Xt
252 pass # Sinon on donne la valeur tiree si on etudie la variabilite de x[0]
254 x[0]=1 # Sinon on laisse la valeur de base
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
259 x[1]=-1 # Sinon on donne -1 comme marqueur
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
264 x[2]=0 # Sinon on considere qu'il n'y a pas d'eolien
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
269 x[3]=0 # Sinon on considere qu'il n'y a pas d'eolien
271 if PV_var==1 : # Si on etudie la variabilite du PV on laisse sa valeur a la va
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])
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
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])
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])
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)
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)
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])
310 for sz in range(0,nx):
311 x2.append(float(x[sz]))
315 elif x[1] < len(continLines) : # L'element tire est une ligne
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) :
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)
336 # Change the bus that is not in service
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])
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])
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])
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])
358 psspy.save(doci) #Saving .sav modifications
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)
366 if psspy.solved()==7:
367 print 'CONVERGENCE CAS '+str(ite)
371 print '==================================================================='
372 print 'NO CONVERGENCE'
373 print '==================================================================='
375 #for i in range (134) :
376 #psspy.opf_bus_indv(i,[_i,0],[_f, 0.7,_f,_f,_f])
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]
383 sizeY4=np.matrix(shunt).shape[0]
384 y=np.zeros(2*sizeY0+sizeY1+3*sizeY2+sizeY3+sizeY4)
386 rate_mat_index=Irate_num+2
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 :
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 :
408 z[6]=z[6]-z[0] # Number of lines between 90% and 100% of their limits
411 for i in range (sizeY4) :
412 final_shunt+=shunt[i][2]
413 z[7]=final_shunt-init_shunt
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])
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)
442 MyLogger(x2,y,z,logCSVfilename,logTXTfilename,ite)
443 #MyMultiLogger (x2, y, sizeY, z, ite, folder, day, fich, hour)
444 return NumericalPoint(z)
446 #===============================================================================
447 # DEFINITION DU WRAPPER - WRAPPER's DEFINITION
448 #===============================================================================
449 # Initialize size output
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]
467 class PSSEWrapperClass(OpenTURNSPythonFunction) :
469 OpenTURNSPythonFunction.__init__(self,5,8)
471 return PSSEFunction(x)
473 # Initialize the folder
474 newpath = folder+"\N"+folderN_1+day
475 if not os.path.exists(newpath): os.makedirs(newpath)
477 # Test the Num. Math. Function
478 pssefun = NumericalMathFunction(PSSEWrapperClass())
480 # Definition of the function to use
481 inputDim = pssefun.getInputDimension()
482 outputDim = pssefun.getOutputDimension()
484 # Initialization of the distribution collection:
485 #aCollection = DistributionCollection()
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
495 #Create a correlation matrix as copulas
496 corr=CorrelationMatrix(inputDim)
519 copula=Copula(NormalCopula(corr))
522 # Create the input probability distribution, args are the distributions, the correlation laws
523 inputDistribution = ComposedDistribution(collectionMarginals, copula)
525 # Create the input random vector
526 """inputRandomVector = RandomVector(inputDistribution)
528 # Create the output variable of interest
529 outputVariableOfInterest = RandomVector(pssefun, inputRandomVector)
530 outputVariableOfInterest.setDescription(pssefun.getOutputDescription())"""
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])+";")
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])+";")
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")
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;;")
608 for name in range (sizeY[0]):
609 f.write("Y:PMachine"+str(plants_base[name][0])+";")
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])+";")
617 for name in range (sizeY[1]):
618 f.write("Y:VBus"+str(buses_base[name][0])+";")
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])+";")
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])+";")
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])+";")
635 for name in range (sizeY[3]):
636 f.write("Y:Ploads "+str(loads_base[name][0])+";")
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])+";")
645 # Start the simulations
647 print "\n\n\n Starting PSEN "+str(montecarlosize)+" simulations"
649 """inputSample=inputRandomVector.getSample(montecarlosize)
650 inputSample.setDescription( ("X0","X1","X2","X3") )
651 inputSample.exportToCSVFile("InputSamples.csv")"""
654 myLHSE = LHSExperiment(inputDistribution,montecarlosize)
655 inputSample = myLHSE.generate()
657 myMCE = MonteCarloExperiment(inputDistribution,montecarlosize)
658 inputSample = myMCE.generate()
663 print 'Probabilistic'
665 outputSampleAll = pssefun(inputSample)#outputVariableOfInterest.getSample(montecarlosize)
668 for i in range (len(time_serie_mat)) :
670 RandomGenerator.SetSeed(i)
673 for j in range (len(time_serie_file)) :
674 if time_serie_file[j] == -1 :
678 Xt.append(time_serie_mat[i][j-n])
683 outputSampleAll = pssefun(inputSample)
685 outputSampleAll.add(pssefun(inputSample))
687 outputDim=outputSampleAll.getDimension()
688 outputSize=outputSampleAll.getSize()
690 outputSample=NumericalSample(0,outputDim)
691 outputSampleMissed=NumericalSample(0,outputDim)
693 for i in range (outputSize):
694 if outputSampleAll[i,5]==0 :
695 outputSampleMissed.add(outputSampleAll[i])
697 outputSample.add(outputSampleAll[i])
700 for i in range (outputDim):
701 outputDescription.append("Y"+str(i))
702 outputSample.setDescription( outputDescription )
704 # Get the empirical mean and standard deviations
705 empMeanX = inputSample.computeMean()
706 empSdX = inputSample.computeStandardDeviationPerComponent()
707 empiricalMean = outputSample.computeMean()
708 empiricalSd = outputSample.computeStandardDeviationPerComponent()
710 f=open(logCSVfilename, 'a')
713 for i in range(0,inputDim):
714 f.write("%f;" % (empMeanX[i]))
716 for i in range(0,outputDim):
717 f.write("%f;" % (empiricalMean[i]))
719 f.write("\nStandard deviation;;")
720 for i in range(0,inputDim):
721 f.write("%f;" % (empSdX[i]))
723 for i in range(0,outputDim):
724 f.write("%f;" % (empiricalSd[i]))
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))
735 print '\n\nSimulated '+str(montecarlosize)+' cases in '+ str(stop_time-start_time)+' seconds. Average '+str((stop_time-start_time)/montecarlosize)+'s per case.'
737 nMissed=int(outputSampleMissed.getSize())
739 print '\n\n Non-convergence rate is '+str(round(nMissed*100/montecarlosize,3))+' % ('+str(outputSampleMissed.getSize())+' cases on '+str(montecarlosize)+')'
741 #graphical_out(inputSample, outputSampleAll, inputDim, outputDim, montecarlosize)