Daniel Varga commited on
Commit
85e0a54
1 Parent(s): 9e9b514
Files changed (1) hide show
  1. app.py +456 -0
app.py ADDED
@@ -0,0 +1,456 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ # port of
2
+ # https://colab.research.google.com/drive/1PJgcJ4ly7x5GuZy344eJeYSODo8trbM4#scrollTo=39F2u-4hvwLU
3
+
4
+ from dataclasses import dataclass
5
+ import numpy as np
6
+ import pandas as pd
7
+
8
+ import matplotlib.pyplot as plt
9
+ import matplotlib
10
+ import datetime
11
+ from scipy.interpolate import interp1d
12
+
13
+ import gradio as gr
14
+ import plotly.express as px
15
+ import plotly.graph_objects as go
16
+
17
+
18
+ #@title ### Downloading the data
19
+ # !wget "https://static.renyi.hu/ai-shared/daniel/pq/PL_44527.19-21.csv"
20
+ # !wget "https://static.renyi.hu/ai-shared/daniel/pq/pq_terheles_2021_adatok.tsv"
21
+
22
+
23
+ PATH_PREFIX = "./"
24
+
25
+ matplotlib.rcParams['figure.figsize'] = [12, 8]
26
+
27
+ START = f"2021-01-01"
28
+ END = f"2022-01-01"
29
+
30
+
31
+ def read_datasets():
32
+ #@title ### Preprocessing meteorologic data
33
+ met_data = pd.read_csv(PATH_PREFIX + 'PL_44527.19-21.csv', sep=';', skipinitialspace=True, na_values='n/a', skiprows=[0, 1, 2, 3, 4])
34
+ met_data['Time'] = met_data['Time'].astype(str)
35
+ date_time = met_data['Time'] = pd.to_datetime(met_data['Time'], format='%Y%m%d%H%M')
36
+ met_data = met_data.set_index('Time')
37
+
38
+
39
+ #@title ### Preprocessing consumption data
40
+ cons_data = pd.read_csv(PATH_PREFIX + 'pq_terheles_2021_adatok.tsv', sep='\t', skipinitialspace=True, na_values='n/a', decimal=',')
41
+ cons_data['Time'] = pd.to_datetime(cons_data['Korrigált időpont'], format='%m/%d/%y %H:%M')
42
+ cons_data = cons_data.set_index('Time')
43
+ cons_data['Consumption'] = cons_data['Hatásos teljesítmény [kW]']
44
+
45
+ # consumption data is at 14 29 44 59 minutes, we move it by 1 minute
46
+ # to sync it with production data:
47
+ cons_data.index = cons_data.index + pd.DateOffset(minutes=1)
48
+
49
+ met_2021_data = met_data[(met_data.index >= START) & (met_data.index < END)]
50
+ cons_2021_data = cons_data[(cons_data.index >= START) & (cons_data.index < END)]
51
+
52
+ return met_2021_data, cons_2021_data
53
+
54
+
55
+ @dataclass
56
+ class Parameters:
57
+ solar_cell_num: float = 114 # units
58
+ solar_efficiency: float = 0.93 * 0.96 # [dimensionless]
59
+ NOCT: float = 280 # [W]
60
+ NOCT_irradiation: float = 800 # [W/m^2]
61
+
62
+ bess_nominal_capacity: float = 330 # [Ah]
63
+ bess_charge: float = 50 # [kW]
64
+ bess_discharge: float = 60 # [kW]
65
+ voltage: float = 600 # [V]
66
+ maximal_depth_of_discharge: float = 0.75 # [dimensionless]
67
+ energy_loss: float = 0.1 # [dimensionless]
68
+
69
+ @property
70
+ def bess_capacity(self):
71
+ return self.bess_nominal_capacity * self.voltage / 1000
72
+
73
+
74
+ # mutates met_2021_data
75
+ def add_production_field(met_2021_data, parameters):
76
+ sr = met_2021_data['sr']
77
+
78
+ nop_total = sr * parameters.solar_cell_num * parameters.solar_efficiency * parameters.NOCT / parameters.NOCT_irradiation / 1e3
79
+ nop_total = nop_total.clip(0)
80
+ met_2021_data['Production'] = nop_total
81
+
82
+
83
+ def interpolate_and_join(met_2021_data, cons_2021_data):
84
+ applicable = 24*60*365 - 15 + 5
85
+
86
+ demand_f = interp1d(range(0, 365*24*60, 15), cons_2021_data['Consumption'])
87
+ #demand_f = interp1d(range(0, 6*24*60, 15), cons_2021_data['Consumption'])
88
+ demand_interp = demand_f(range(0, applicable, 5))
89
+
90
+ production_f = interp1d(range(0, 365*24*60, 10), met_2021_data['Production'])
91
+ #production_f = interp1d(range(0, 6*24*60, 10), met_2021_data['Production'])
92
+ production_interp = production_f(range(0, applicable, 5))
93
+
94
+ all_2021_datetimeindex = pd.date_range(start=START, end=END, freq='5min')[:len(production_interp)]
95
+
96
+ all_2021_data = pd.DataFrame({'Consumption': demand_interp, 'Production': production_interp})
97
+ all_2021_data = all_2021_data.set_index(all_2021_datetimeindex)
98
+ return all_2021_data
99
+
100
+
101
+ def simulator_with_solar(all_data, parameters):
102
+ demand_np = all_data['Consumption'].to_numpy()
103
+ production_np = all_data['Production'].to_numpy()
104
+ assert len(demand_np) == len(production_np)
105
+ step_in_minutes = all_data.index.freq.n
106
+ print("Simulating for", len(demand_np), "time steps. Each step is", step_in_minutes, "minutes.")
107
+ soc_series = [] # soc = state_of_charge.
108
+ # by convention, we only call end user demand, demand,
109
+ # and we only call end user consumption, consumption.
110
+ # in our simple model, demand is always satisfied, hence demand=consumption.
111
+ # BESS demand is called charge.
112
+ consumption_from_solar_series = [] # demand satisfied by solar production
113
+ consumption_from_network_series = [] # demand satisfied by network
114
+ consumption_from_bess_series = [] # demand satisfied by BESS
115
+ # the previous three must sum to demand_series.
116
+ charge_of_bess_series = [] # power taken from solar by BESS. note: power never taken from network by BESS.
117
+ discarded_production_series = [] # solar power thrown away
118
+
119
+ # 1 is not nominal but targeted (healthy) maximum charge.
120
+ # we start with an empty battery, but not emptier than what's healthy for the batteries.
121
+
122
+ #Remark from Jutka
123
+ #For the sake of simplicity 0<= soc <=1
124
+ #soc = 1 - maximal_depth_of_discharge
125
+ #and will use only maximal_depth_of_discharge percent of the real battery capacity
126
+ soc = 0
127
+ max_cap_of_battery = parameters.bess_capacity * parameters.maximal_depth_of_discharge
128
+ cap_of_battery = soc * max_cap_of_battery
129
+
130
+ time_interval = step_in_minutes / 60 # amount of time step in hours
131
+ for i, (demand, production) in enumerate(zip(demand_np, production_np)):
132
+
133
+ # these three are modified on the appropriate codepaths:
134
+ consumption_from_solar = 0
135
+ consumption_from_bess = 0
136
+ consumption_from_network = 0
137
+ charge_of_bess = 0
138
+
139
+ #Remark: If the consumption stable for ex. 10 kwh
140
+ # demand = 10
141
+ unsatisfied_demand = demand
142
+ remaining_production = production # max((production, 0))
143
+ discarded_production = 0
144
+
145
+ # crucially, we never charge the BESS from the network.
146
+ # if demand >= production:
147
+ # all goes to demand
148
+ # we try to cover the rest from BESS
149
+ # we cover the rest from network
150
+ # else:
151
+ # demand fully satisfied by production
152
+ # if exploitable production still remains:
153
+ # if is_battery_chargeable:
154
+ # charge_battery
155
+ # else:
156
+ # log discarded production
157
+
158
+ #battery_charged_enough = (soc > 1- maximal_depth_of_discharge)
159
+ is_battery_charged_enough = (soc > 0 )
160
+ is_battery_chargeable = (soc < 1.0)
161
+
162
+ if unsatisfied_demand >= remaining_production:
163
+ # all goes to demand
164
+ consumption_from_solar = remaining_production
165
+ unsatisfied_demand -= consumption_from_solar
166
+ remaining_production=0 #edited by Jutka
167
+ # we try to cover the rest from BESS
168
+
169
+ if unsatisfied_demand > 0:
170
+ if is_battery_charged_enough:
171
+ # simplifying assumption for now:
172
+ # throughput is enough to completely fulfill extra demand.
173
+ # TODO get rid of simplifying assumption.
174
+ # Remarks from Jutka
175
+ # It is a very bed assumption. The reality needs, that the BESS has limited capacity.
176
+ #
177
+ #
178
+ # cap_of_bess=soc * bess_capacity
179
+ # if cap_of_bess > unsatisfied_demand
180
+ # consumption_from_bess = unsatisfied_demand
181
+ # unsatisfied_demand = 0
182
+ # cap_of_bess -= consumption_from_bess
183
+ # soc = cap_of_bess / bess_capacity
184
+ # else: unsatisfied_demand -= cap_of_bess
185
+ # cap_of_bess = 0
186
+ # soc = 0
187
+ # if unsatisfied_demand > 0
188
+ # consumption_from_network = unsatisfied_demand
189
+ # unsatisfied_demand = 0
190
+
191
+ #Remarks: battery capacity is limited!
192
+
193
+ if cap_of_battery >= unsatisfied_demand * time_interval :
194
+
195
+ #discharge_of_bess = min ( unsatisfied_demand, bess_discharge )
196
+ #discharge = discharge_of_bess
197
+ #consumption_from_bess = discharge * time_interval
198
+ consumption_from_bess = unsatisfied_demand
199
+ #unsatisfied_demand -= consumption_from_bess
200
+ unsatisfied_demand = 0
201
+ cap_of_battery -= consumption_from_bess * time_interval
202
+ soc = cap_of_battery / max_cap_of_battery
203
+
204
+ else:
205
+ #discharge_of_bess = cap_of_battery /time_interval
206
+ #discharge = min( bess_discharge, discharge_of_bess )
207
+ consumption_from_bess = cap_of_battery / time_interval
208
+ unsatisfied_demand -= consumption_from_bess
209
+ cap_of_battery -=consumption_from_bess * time_interval
210
+ soc = cap_of_battery / max_cap_of_battery
211
+ consumption_from_network = unsatisfied_demand
212
+ unsatisfied_demand = 0
213
+ #bess_sacrifice = consumption_from_bess / (1 - energy_loss) # kW
214
+ #energy = bess_sacrifice * time_interval # kWh
215
+ #soc -= energy / bess_capacity
216
+ # print("soc after discharge", soc)
217
+ #consumption_from_network = unsatisfied_demand
218
+ #unsatisfied_demand = 0
219
+ else:
220
+ # we cover the rest from network
221
+ consumption_from_network = unsatisfied_demand
222
+ unsatisfied_demand = 0
223
+
224
+ else:
225
+ # demand fully satisfied by production
226
+
227
+ consumption_from_solar = unsatisfied_demand
228
+ remaining_production -= unsatisfied_demand
229
+ unsatisfied_demand = 0
230
+ # if exploitable production still remains:
231
+ if remaining_production > 0:
232
+ if is_battery_chargeable:
233
+ charge_of_bess = remaining_production
234
+ energy = charge_of_bess * time_interval # kWh
235
+ #Remarks: battery alowed to charge until its capacity maximum
236
+ #energy_charge = min(energy, max_cap_of_battery-cap_of_battery)
237
+ cap_of_battery += energy
238
+ #soc += energy / bess_capacity
239
+ soc = cap_of_battery / max_cap_of_battery
240
+ #print("soc after charge", soc)
241
+ else:
242
+ discarded_production = remaining_production
243
+
244
+ soc_series.append(soc)
245
+ consumption_from_solar_series.append(consumption_from_solar)
246
+ consumption_from_network_series.append(consumption_from_network)
247
+ consumption_from_bess_series.append(consumption_from_bess)
248
+ charge_of_bess_series.append(charge_of_bess)
249
+ discarded_production_series.append(discarded_production)
250
+
251
+ soc_series = np.array(soc_series)
252
+ consumption_from_solar_series = np.array(consumption_from_solar_series)
253
+ consumption_from_network_series = np.array(consumption_from_network_series)
254
+ consumption_from_bess_series = np.array(consumption_from_bess_series)
255
+ charge_of_bess_series = np.array(charge_of_bess_series)
256
+ discarded_production_series = np.array(discarded_production)
257
+
258
+ results = pd.DataFrame({'soc_series': soc_series, 'consumption_from_solar': consumption_from_solar_series,
259
+ 'consumption_from_network': consumption_from_network_series,
260
+ 'consumption_from_bess': consumption_from_bess_series,
261
+ 'charge_of_bess': charge_of_bess_series,
262
+ 'discarded_production': discarded_production_series,
263
+ 'Consumption': all_data['Consumption'],
264
+ 'Production': all_data['Production']
265
+ })
266
+ results = results.set_index(all_data.index)
267
+ return results
268
+
269
+
270
+ def visualize_simulation(results, date_range):
271
+ start_date, end_date = date_range
272
+
273
+ fig = plt.figure()
274
+ results = results.loc[start_date: end_date]
275
+
276
+ x = results.index
277
+ y = [results.consumption_from_solar, results.consumption_from_network, results.consumption_from_bess]
278
+ plt.plot(x, y[0], label='Demand served by solar', color='yellow', linewidth=0.5)
279
+ plt.plot(x, y[0]+y[1], label='Demand served by network', color='blue', linewidth=0.5)
280
+ plt.plot(x, y[0]+y[1]+y[2], label='Demand served by BESS', color='green', linewidth=0.5)
281
+ plt.fill_between(x, y[0]+y[1]+y[2], 0, color='green')
282
+ plt.fill_between(x, y[0]+y[1], 0, color='blue')
283
+ plt.fill_between(x, y[0], 0, color='yellow')
284
+
285
+ # plt.xlim(datetime.datetime.fromisoformat(start_date), datetime.datetime.fromisoformat(end_date))
286
+
287
+ plt.legend()
288
+ return fig
289
+
290
+
291
+ def plotly_visualize_simulation(results, date_range):
292
+ start_date, end_date = date_range
293
+ results = results.loc[start_date: end_date]
294
+ '''
295
+ fig = px.area(results, x=results.index, y="consumption_from_network")
296
+ return fig'''
297
+ fig = go.Figure()
298
+ fig.add_trace(go.Scatter(
299
+ x=results.index, y=results['consumption_from_network'],
300
+ hoverinfo='x+y',
301
+ mode='lines',
302
+ line=dict(width=0.5, color='blue'),
303
+ name='Network',
304
+ stackgroup='one' # define stack group
305
+ ))
306
+ fig.add_trace(go.Scatter(
307
+ x=results.index, y=results['consumption_from_solar'],
308
+ hoverinfo='x+y',
309
+ mode='lines',
310
+ line=dict(width=0.5, color='orange'),
311
+ name='Solar',
312
+ stackgroup='one'
313
+ ))
314
+ fig.add_trace(go.Scatter(
315
+ x=results.index, y=results['consumption_from_bess'],
316
+ hoverinfo='x+y',
317
+ mode='lines',
318
+ line=dict(width=0.5, color='green'),
319
+ name='BESS',
320
+ stackgroup='one'
321
+ ))
322
+ fig.update_layout(
323
+ height=400
324
+ )
325
+ return fig
326
+
327
+
328
+ def monthly_analysis(results):
329
+ consumptions = []
330
+ for month in range(1, 13):
331
+ start = f"2021-{month:02}-01"
332
+ end = f"2021-{month+1:02}-01"
333
+ if month == 12:
334
+ end = "2022-01-01"
335
+ results_in_month = results[(results.index >= start) & (results.index < end)]
336
+
337
+ total = results_in_month['Consumption'].sum()
338
+ network = results_in_month['consumption_from_network'].sum()
339
+ solar = results_in_month['consumption_from_solar'].sum()
340
+ bess = results_in_month['consumption_from_bess'].sum()
341
+ consumptions.append([network, solar, bess])
342
+
343
+ consumptions = np.array(consumptions)
344
+ step_in_minutes = results.index.freq.n
345
+ # consumption is given in kW. each tick is step_in_minutes long (5mins, in fact)
346
+ # we get consumption in kWh if we multiply sum by step_in_minutes/60
347
+ consumptions_in_mwh = consumptions * (step_in_minutes / 60) / 1000
348
+ return consumptions_in_mwh
349
+
350
+
351
+ def monthly_visualization(consumptions_in_mwh):
352
+ percentages = consumptions_in_mwh[:, :3] / consumptions_in_mwh.sum(axis=1, keepdims=True) * 100
353
+ bats = 0
354
+ nws = 0
355
+ sols = 0
356
+
357
+ print("[Mwh]")
358
+ print("==========================")
359
+ print("month\tnetwork\tsolar\tbess")
360
+ for month_minus_1 in range(12):
361
+ network, solar, bess = consumptions_in_mwh[month_minus_1]
362
+ print(f"{month_minus_1+1}\t{network:0.2f}\t{solar:0.2f}\t{bess:0.2f}")
363
+ bats += bess
364
+ nws += network
365
+ sols += solar
366
+ print(f"\t{nws:0.2f}\t{sols:0.2f}\t{bats:0.2f}")
367
+
368
+
369
+ fig, ax = plt.subplots()
370
+
371
+ ax.stackplot(range(1, 13),
372
+ percentages[:, 0], percentages[:, 1], percentages[:, 2],
373
+ labels=["hálózat", "egyenesen a naptól", "a naptól a BESS-en keresztül"])
374
+ ax.set_ylim(0, 100)
375
+ ax.legend()
376
+ plt.title('A fogyasztás hány százalékát fedezte az adott hónapban?')
377
+ plt.show()
378
+
379
+ plt.stackplot(range(1, 13),
380
+ consumptions_in_mwh[:, 0], consumptions_in_mwh[:, 1], consumptions_in_mwh[:, 2],
381
+ labels=["hálózat", "egyenesen a naptól", "a naptól a BESS-en keresztül"])
382
+ plt.legend()
383
+ plt.title('Mennyi fogyasztást fedezett az adott hónapban? [MWh]')
384
+ plt.show()
385
+
386
+
387
+ def main():
388
+ parameters = Parameters()
389
+
390
+ met_2021_data, cons_2021_data = read_datasets()
391
+
392
+ add_production_field(met_2021_data, parameters)
393
+
394
+ all_2021_data = interpolate_and_join(met_2021_data, cons_2021_data)
395
+
396
+ results = simulator_with_solar(all_2021_data, parameters)
397
+
398
+ fig = visualize_simulation(results, date_range=("2021-02-01", "2021-03-01"))
399
+ plt.show()
400
+
401
+ consumptions_in_mwh = monthly_analysis(results)
402
+ monthly_visualization(consumptions_in_mwh)
403
+
404
+
405
+ # main() ; exit()
406
+
407
+
408
+ met_2021_data, cons_2021_data = read_datasets()
409
+
410
+
411
+ def recalculate(**uiParameters):
412
+ fixed_consumption = uiParameters['fixed_consumption']
413
+ del uiParameters['fixed_consumption']
414
+
415
+ parameters = Parameters()
416
+ for k, v in uiParameters.items():
417
+ setattr(parameters, k, v)
418
+
419
+ add_production_field(met_2021_data, parameters)
420
+ all_2021_data = interpolate_and_join(met_2021_data, cons_2021_data)
421
+
422
+ if fixed_consumption:
423
+ all_2021_data['Consumption'] = 10
424
+
425
+ results = simulator_with_solar(all_2021_data, parameters)
426
+ return results
427
+
428
+
429
+ def ui_refresh(solar_cell_num, bess_nominal_capacity, fixed_consumption):
430
+ results = recalculate(solar_cell_num=solar_cell_num, bess_nominal_capacity=bess_nominal_capacity, fixed_consumption=fixed_consumption)
431
+
432
+ fig1 = plotly_visualize_simulation(results, date_range=("2021-02-01", "2021-02-07"))
433
+ fig2 = plotly_visualize_simulation(results, date_range=("2021-08-02", "2021-08-08"))
434
+
435
+ # (12, 3), the 3 indexed with (network, solar, bess):
436
+ consumptions_in_mwh = monthly_analysis(results)
437
+
438
+ network, solar, bess = consumptions_in_mwh.sum(axis=0)
439
+ html = ""
440
+ for column, column_name in zip((network, solar, bess), ("network", "solar", "BESS")):
441
+ html += f"Yearly consumption satisfied by {column_name}: {column:0.2f} MWh<br>"
442
+
443
+ return (fig1, fig2, html)
444
+
445
+
446
+ ui = gr.Interface(
447
+ ui_refresh,
448
+ inputs = [
449
+ gr.Slider(0, 2000, 114, label="Solar cell number"),
450
+ gr.Slider(0, 1000, 330, label="BESS nominal capacity"),
451
+ gr.Checkbox(value=False, label="Fixed consumption")],
452
+ outputs = ["plot", "plot", "html"],
453
+ live=True,
454
+ )
455
+
456
+ ui.launch()