-
Notifications
You must be signed in to change notification settings - Fork 2
/
test_opt_pf_30_rt_1_field_der.py
170 lines (132 loc) · 7.58 KB
/
test_opt_pf_30_rt_1_field_der.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
# -*- coding: utf-8 -*-
"""
Created on Tue Feb 5 13:58:12 2019
@author: cwhanse
2019-2-19 [email protected]: Updated with communication interfaces
"""
import pytz
from optimization.ce_api import CE_API
from optimization.feeder_class import Feeder
from forecasting.pv_system_class import PVobj
from optimization.dss_util import PFOptim
import numpy as np
import time
import os
import math
USMtn = pytz.timezone('US/Mountain')
username = "fake"
password = "fake"
# Create the Connected Energy API interface to get PV forecasts and issue DER PV setpoints
api = CE_API(username=username, password=password)
# DER devices represented by the EPRI PV simulators
pv_sizes = [684000, 93353, 5393.75, 2242.35, 44555.7, 152388, 22180.9, 18413.3, 121307, 54932.6, 27494.8, 29935.9,
59594.4, 57685.9, 280784, 129058, 58194.7, 13166.9, 52328.8, 34853.4, 22675.2, 9156.04, 1012210,
1028690, 993513, 373251, 334862, 344033, 27142.1, 23568.6, 684]
pv_names = [1, 851, 838, 869, 857, 858, 859, 872, 873, 874, 848, 849, 850, 863, 864, 865, 845, 846, 847, 854, 855, 856,
842, 843, 844, 860, 861, 862, 866, 867, 868]
# Create dictionary of all PV systems and associated forecast information
pvdict = {}
dss_to_phil_map = {}
pv_count = 0
for pv in pv_sizes:
dss_name = pv_names[pv_count]
pv_count += 1
if pv_count == 1:
pvdict['ng1'] = PVobj('ng1', dc_capacity=pv, ac_capacity=pv, lat=35.05,
lon=-106.54, alt=1657, tz=USMtn, tilt=35, azimuth=180, pf_max=0.85,
pf_min=-0.85, forecast_method='persistence')
pvdict['pvsy%d' % dss_name] = PVobj('epri%d' % pv_count, dc_capacity=pv, ac_capacity=pv, lat=35.05,
lon=-106.54, alt=1657, tz=USMtn, tilt=35, azimuth=180, pf_max=0.85,
pf_min=-0.85, surrogateid='ng1', forecast_method='persistence')
dss_to_phil_map['pvsy%d' % dss_name] = 'epri%d' % pv_count
# Setup feeder object which points to the local OpenDSS time series simulation
cwd = os.getcwd()
feeder = Feeder(filename=cwd + "\\NG_reduced_timeseries\\Master.DSS", pv=pvdict)
pvlist = feeder.pv_on_feeder
loadlist = feeder.DSS.circuit.Loads.AllNames
# stepsize hardcoded at 5m for now because that's the resolution of the PV forecasts
stepsize = '5m'
# Number of stepsize time periods that the optimization will be working over
periods = 3
# Results dataset
data_set_start_time = time.time()
results_filename = cwd + '\\Optimization %s.csv' % data_set_start_time
csvfile = open(results_filename, 'x')
csvfile.write('Time (s), Weighting 1 (W1), Weighting 2 (W2), Weighting 3 (W3), '
'Violation Threshold, Acceptance Threshold, Objective Function Threshold, '
'Power Factor DER 1, Reactive Power DER 1 (var), '
'Objective Function with Prior PFs, PF Change Bool, '
'Optimization Function, PV Forecast, Power Forecast, Reactive Power Forecast\n')
csvfile.close()
sim_duration = 5.*60.*60. # 5 hours (in seconds)
sec_per_loop = 60.0 # number of seconds between optimizations
n_iter = math.ceil(sim_duration/sec_per_loop) # number of 5-minute optimizations
starttime = time.time()
for opt_loop in range(n_iter):
print('#### Running Optimization Loop %i. Total Time: %0.2f seconds ####' % (opt_loop + 1, time.time() - starttime))
''' Update PV forecast data in the OpenDSS simulations '''
# Get the new forecast from the Connected Energy system for each of the PV systems in the feeder.pv_on_feeder
feeder.update_pv_forecasts(api=api)
# scale PV forecasts from W to p.u.
pv_forecast = {pv: feeder.pv_forecasts[pv].values[:periods] / pvdict[pv].ac_capacity for pv in feeder.pv_on_feeder}
print('pv_forecast: %s' % pv_forecast)
''' Update Load information in the OpenDSS simulations '''
# Get the new load data from the state estimator. Assume persistence and forecast the same value into the future
p_forecast, q_forecast = feeder.get_load_forecast(periods=periods)
print('p_forecast: %s' % p_forecast)
print('q_forecast: %s' % q_forecast)
''' Control and options for optimization '''
penalty = {'violation': 1.0, 'deviation': 2.0, 'power_factor': 0.05}
# Voltage violation is at ANSI Range A
# Optimization will not run if all voltages are within 0.2% of Vnom
# PFs will not change if the new PFs do not improve objective function by object%
threshold = {'violation': 0.05, 'accept': 0.002, 'object': 1.0e-7}
debug = True
swarmsize = 60
maxiter = 10
minstep = 0.001
minfunc = 1e-6
options = PFOptim(penalty=penalty, threshold=threshold, debug=debug, swarmsize=swarmsize, maxiter=maxiter,
minstep=minstep, minfunc=minfunc)
new_pf, prior_pf, opt_obj, prior_obj, change = feeder.update_power_factors(pvlist, pv_forecast,
p_forecast, q_forecast, hour=0,
sec=0, stepsize=stepsize,
numsteps=periods, options=options,
controllable_pv=[pvlist[0]])
print('The optimal power factors are %s' % new_pf)
''' Write the optimal PF values to the DER devices using the Connected Energy API '''
der = {}
new_q_values = {}
pv_names = [] # controller DER for the save file
for pv_name in new_pf.keys():
if new_pf[pv_name] >= 0: # EPRI simulator convention
excitation = "injectingQ"
# the expected reactive power based on the forecast VA level of the DERs
new_q_values[pv_name] = np.arccos(new_pf[pv_name])*pv_forecast[pv_name][0]*pvdict[pv_name].ac_capacity
else:
excitation = "absorbingQ"
# the expected reactive power based on the forecast VA level of the DERs
new_q_values[pv_name] = -np.arccos(-new_pf[pv_name])*pv_forecast[pv_name][0]*pvdict[pv_name].ac_capacity
der[dss_to_phil_map[pv_name]] = {'excitation': excitation, 'pf': new_pf[pv_name], 'forecast': None}
if dss_to_phil_map[pv_name] == 'epri1': # write to the physical device
der['ng1'] = {'excitation': excitation, 'pf': new_pf[pv_name], 'forecast': None}
pv_names.append(pv_name)
print('The optimal reactive power values are %s' % new_q_values)
print('DER dictionary to send to CE\'s API: %s' % der)
api.set_pf(der=der)
# store detailed results information in individual files for replay debugging
step_time = time.time() - starttime
# Write data from this optimization loop
results = '%s, %0.5f, %0.5f, %0.5f, %0.5f, %0.5f, %0.5f, %0.5f, %0.5f, %0.5f, ' \
'%s, %s, %s, %s, %s \n' % (step_time, penalty['violation'], penalty['deviation'], penalty['power_factor'],
threshold['violation'], threshold['accept'], threshold['object'],
new_pf[pv_names[0]], new_q_values[pv_names[0]], prior_obj, change,
opt_obj, pv_forecast, p_forecast, q_forecast)
# open/close results file each time so the latest data is stored even with a crash
csvfile = open(results_filename, 'a')
csvfile.write(results)
csvfile.close()
runtime = time.time() - starttime
wait_time = sec_per_loop - (runtime % sec_per_loop)
print('Total Run Time: %s. Waiting %0.1f seconds until recalculating optimal PFs.' % (runtime, wait_time))
time.sleep(wait_time) # Run code every sec_per_loop