Passed
Pull Request — dev (#924)
by
unknown
01:37
created

data.datasets.sanity_checks.sanitycheck_eGon2035_electricity()   F

Complexity

Conditions 16

Size

Total Lines 230
Code Lines 97

Duplication

Lines 64
Ratio 27.83 %

Importance

Changes 0
Metric Value
eloc 97
dl 64
loc 230
rs 1.7781
c 0
b 0
f 0
cc 16
nop 0

How to fix   Long Method    Complexity   

Long Method

Small methods make your code easier to understand, in particular if combined with a good name. Besides, if your method is small, finding a good name is usually much easier.

For example, if you find yourself adding comments to a method's body, this is usually a good sign to extract the commented part to a new method, and use the comment as a starting point when coming up with a good name for this new method.

Commonly applied refactorings include:

Complexity

Complex classes like data.datasets.sanity_checks.sanitycheck_eGon2035_electricity() often do a lot of different things. To break such a class down, we need to identify a cohesive component within that class. A common approach to find such a component is to look for fields/methods that share the same prefixes, or suffixes.

Once you have determined the fields that belong together, you can apply the Extract Class refactoring. If the component makes sense as a sub-class, Extract Subclass is also a candidate, and is often faster.

1
"""
2
This module does sanity checks for both the eGon2035 and the eGon100RE scenario
3
separately where a percentage error is given to showcase difference in output
4
and input values. Please note that there are missing input technologies in the
5
supply tables.
6
Authors: @ALonso, @dana, @nailend, @nesnoj
7
"""
8
9
from sqlalchemy import Numeric
10
from sqlalchemy.sql import and_, cast, func, or_
11
import numpy as np
12
import pandas as pd
13
14
from egon.data import config, db, logger
15
from egon.data.datasets import Dataset
16
from egon.data.datasets.electricity_demand_timeseries.cts_buildings import (
17
    EgonCtsElectricityDemandBuildingShare,
18
    EgonCtsHeatDemandBuildingShare,
19
from egon.data.datasets.emobility.motorized_individual_travel.db_classes import (
20
    EgonEvCountMunicipality,
21
    EgonEvCountMvGridDistrict,
22
    EgonEvCountRegistrationDistrict,
23
    EgonEvMvGridDistrict,
24
    EgonEvPool,
25
    EgonEvTrip,
26
)
27
from egon.data.datasets.emobility.motorized_individual_travel.helpers import (
28
    DATASET_CFG,
29
    read_simbev_metadata_file,
30
)
31
from egon.data.datasets.etrago_setup import (
32
    EgonPfHvLink,
33
    EgonPfHvLinkTimeseries,
34
    EgonPfHvLoad,
35
    EgonPfHvLoadTimeseries,
36
    EgonPfHvStore,
37
    EgonPfHvStoreTimeseries,
38
)
39
from egon.data.datasets.scenario_parameters import get_sector_parameters
40
41
TESTMODE_OFF = (
42
    config.settings()["egon-data"]["--dataset-boundary"] == "Everything"
43
)
44
45
46
class SanityChecks(Dataset):
47
    def __init__(self, dependencies):
48
        super().__init__(
49
            name="SanityChecks",
50
            version="0.0.5",
51
            dependencies=dependencies,
52
            tasks={
53
                etrago_eGon2035_electricity,
54
                etrago_eGon2035_heat,
55
                residential_electricity_annual_sum,
56
                residential_electricity_hh_refinement,
57
                cts_electricity_demand_share,
58
                cts_heat_demand_share,
59
                sanitycheck_emobility_mit,
60
            },
61
        )
62
63
64
def etrago_eGon2035_electricity():
65
    """Execute basic sanity checks.
66
67
    Returns print statements as sanity checks for the electricity sector in
68
    the eGon2035 scenario.
69
70
    Parameters
71
    ----------
72
    None
73
74
    Returns
75
    -------
76
    None
77
    """
78
79
    scn = "eGon2035"
80
81
    # Section to check generator capacities
82
    logger.info(f"Sanity checks for scenario {scn}")
83
    logger.info(
84
        "For German electricity generators the following deviations between "
85
        "the inputs and outputs can be observed:"
86
    )
87
88
    carriers_electricity = [
89
        "other_non_renewable",
90
        "other_renewable",
91
        "reservoir",
92
        "run_of_river",
93
        "oil",
94
        "wind_onshore",
95
        "wind_offshore",
96
        "solar",
97
        "solar_rooftop",
98
        "biomass",
99
    ]
100
101
    for carrier in carriers_electricity:
102
103
        if carrier == "biomass":
104
            sum_output = db.select_dataframe(
105
                """SELECT scn_name, SUM(p_nom::numeric) as output_capacity_mw
106
                    FROM grid.egon_etrago_generator
107
                    WHERE bus IN (
108
                        SELECT bus_id FROM grid.egon_etrago_bus
109
                        WHERE scn_name = 'eGon2035'
110
                        AND country = 'DE')
111
                    AND carrier IN ('biomass', 'industrial_biomass_CHP',
112
                    'central_biomass_CHP')
113
                    GROUP BY (scn_name);
114
                """,
115
                warning=False,
116
            )
117
118
        else:
119
            sum_output = db.select_dataframe(
120
                f"""SELECT scn_name,
121
                 SUM(p_nom::numeric) as output_capacity_mw
122
                         FROM grid.egon_etrago_generator
123
                         WHERE scn_name = '{scn}'
124
                         AND carrier IN ('{carrier}')
125
                         AND bus IN
126
                             (SELECT bus_id
127
                               FROM grid.egon_etrago_bus
128
                               WHERE scn_name = 'eGon2035'
129
                               AND country = 'DE')
130
                         GROUP BY (scn_name);
131
                    """,
132
                warning=False,
133
            )
134
135
        sum_input = db.select_dataframe(
136
            f"""SELECT carrier, SUM(capacity::numeric) as input_capacity_mw
137
                     FROM supply.egon_scenario_capacities
138
                     WHERE carrier= '{carrier}'
139
                     AND scenario_name ='{scn}'
140
                     GROUP BY (carrier);
141
                """,
142
            warning=False,
143
        )
144
145
        if (
146
            sum_output.output_capacity_mw.sum() == 0
147
            and sum_input.input_capacity_mw.sum() == 0
148
        ):
149
            logger.info(
150
                f"No capacity for carrier '{carrier}' needed to be"
151
                f" distributed. Everything is fine"
152
            )
153
154
        elif (
155
            sum_input.input_capacity_mw.sum() > 0
156
            and sum_output.output_capacity_mw.sum() == 0
157
        ):
158
            logger.info(
159
                f"Error: Capacity for carrier '{carrier}' was not distributed "
160
                f"at all!"
161
            )
162
163
        elif (
164
            sum_output.output_capacity_mw.sum() > 0
165
            and sum_input.input_capacity_mw.sum() == 0
166
        ):
167
            logger.info(
168
                f"Error: Eventhough no input capacity was provided for carrier"
169
                f"'{carrier}' a capacity got distributed!"
170
            )
171
172
        else:
173
            sum_input["error"] = (
174
                (sum_output.output_capacity_mw - sum_input.input_capacity_mw)
175
                / sum_input.input_capacity_mw
176
            ) * 100
177
            g = sum_input["error"].values[0]
178
179
            logger.info(f"{carrier}: " + str(round(g, 2)) + " %")
180
181
    # Section to check storage units
182
183
    logger.info(f"Sanity checks for scenario {scn}")
184
    logger.info(
185
        "For German electrical storage units the following deviations between"
186
        "the inputs and outputs can be observed:"
187
    )
188
189
    carriers_electricity = ["pumped_hydro"]
190
191
    for carrier in carriers_electricity:
192
193
        sum_output = db.select_dataframe(
194
            f"""SELECT scn_name, SUM(p_nom::numeric) as output_capacity_mw
195
                         FROM grid.egon_etrago_storage
196
                         WHERE scn_name = '{scn}'
197
                         AND carrier IN ('{carrier}')
198
                         AND bus IN
199
                             (SELECT bus_id
200
                               FROM grid.egon_etrago_bus
201
                               WHERE scn_name = 'eGon2035'
202
                               AND country = 'DE')
203
                         GROUP BY (scn_name);
204
                    """,
205
            warning=False,
206
        )
207
208
        sum_input = db.select_dataframe(
209
            f"""SELECT carrier, SUM(capacity::numeric) as input_capacity_mw
210
                     FROM supply.egon_scenario_capacities
211
                     WHERE carrier= '{carrier}'
212
                     AND scenario_name ='{scn}'
213
                     GROUP BY (carrier);
214
                """,
215
            warning=False,
216
        )
217
218
        if (
219
            sum_output.output_capacity_mw.sum() == 0
220
            and sum_input.input_capacity_mw.sum() == 0
221
        ):
222
            logger.info(
223
                f"No capacity for carrier '{carrier}' needed to be "
224
                f"distributed. Everything is fine"
225
            )
226
227
        elif (
228
            sum_input.input_capacity_mw.sum() > 0
229
            and sum_output.output_capacity_mw.sum() == 0
230
        ):
231
            logger.info(
232
                f"Error: Capacity for carrier '{carrier}' was not distributed"
233
                f" at all!"
234
            )
235
236
        elif (
237
            sum_output.output_capacity_mw.sum() > 0
238
            and sum_input.input_capacity_mw.sum() == 0
239
        ):
240
            logger.info(
241
                f"Error: Eventhough no input capacity was provided for carrier"
242
                f" '{carrier}' a capacity got distributed!"
243
            )
244
245
        else:
246
            sum_input["error"] = (
247
                (sum_output.output_capacity_mw - sum_input.input_capacity_mw)
248
                / sum_input.input_capacity_mw
249
            ) * 100
250
            g = sum_input["error"].values[0]
251
252
            logger.info(f"{carrier}: " + str(round(g, 2)) + " %")
253
254
    # Section to check loads
255
256
    logger.info(
257
        "For German electricity loads the following deviations between the"
258
        " input and output can be observed:"
259
    )
260
261
    output_demand = db.select_dataframe(
262
        """SELECT a.scn_name, a.carrier,  SUM((SELECT SUM(p)
263
        FROM UNNEST(b.p_set) p))/1000000::numeric as load_twh
264
            FROM grid.egon_etrago_load a
265
            JOIN grid.egon_etrago_load_timeseries b
266
            ON (a.load_id = b.load_id)
267
            JOIN grid.egon_etrago_bus c
268
            ON (a.bus=c.bus_id)
269
            AND b.scn_name = 'eGon2035'
270
            AND a.scn_name = 'eGon2035'
271
            AND a.carrier = 'AC'
272
            AND c.scn_name= 'eGon2035'
273
            AND c.country='DE'
274
            GROUP BY (a.scn_name, a.carrier);
275
276
    """,
277
        warning=False,
278
    )["load_twh"].values[0]
279
280
    input_cts_ind = db.select_dataframe(
281
        """SELECT scenario,
282
         SUM(demand::numeric/1000000) as demand_mw_regio_cts_ind
283
            FROM demand.egon_demandregio_cts_ind
284
            WHERE scenario= 'eGon2035'
285
            AND year IN ('2035')
286
            GROUP BY (scenario);
287
288
        """,
289
        warning=False,
290
    )["demand_mw_regio_cts_ind"].values[0]
291
292
    input_hh = db.select_dataframe(
293
        """SELECT scenario, SUM(demand::numeric/1000000) as demand_mw_regio_hh
294
            FROM demand.egon_demandregio_hh
295
            WHERE scenario= 'eGon2035'
296
            AND year IN ('2035')
297
            GROUP BY (scenario);
298
        """,
299
        warning=False,
300
    )["demand_mw_regio_hh"].values[0]
301
302
    input_demand = input_hh + input_cts_ind
303
304
    e = round((output_demand - input_demand) / input_demand, 2) * 100
305
306
    logger.info(f"electricity demand: {e} %")
307
308
309
def etrago_eGon2035_heat():
310
    """Execute basic sanity checks.
311
312
    Returns print statements as sanity checks for the heat sector in
313
    the eGon2035 scenario.
314
315
    Parameters
316
    ----------
317
    None
318
319
    Returns
320
    -------
321
    None
322
    """
323
324
    # Check input and output values for the carriers "other_non_renewable",
325
    # "other_renewable", "reservoir", "run_of_river" and "oil"
326
327
    scn = "eGon2035"
328
329
    # Section to check generator capacities
330
    logger.info(f"Sanity checks for scenario {scn}")
331
    logger.info(
332
        "For German heat demands the following deviations between the inputs"
333
        " and outputs can be observed:"
334
    )
335
336
    # Sanity checks for heat demand
337
338
    output_heat_demand = db.select_dataframe(
339
        """SELECT a.scn_name,
340
          (SUM(
341
          (SELECT SUM(p) FROM UNNEST(b.p_set) p))/1000000)::numeric as load_twh
342
            FROM grid.egon_etrago_load a
343
            JOIN grid.egon_etrago_load_timeseries b
344
            ON (a.load_id = b.load_id)
345
            JOIN grid.egon_etrago_bus c
346
            ON (a.bus=c.bus_id)
347
            AND b.scn_name = 'eGon2035'
348
            AND a.scn_name = 'eGon2035'
349
            AND c.scn_name= 'eGon2035'
350
            AND c.country='DE'
351
            AND a.carrier IN ('rural_heat', 'central_heat')
352
            GROUP BY (a.scn_name);
353
        """,
354
        warning=False,
355
    )["load_twh"].values[0]
356
357
    input_heat_demand = db.select_dataframe(
358
        """SELECT scenario, SUM(demand::numeric/1000000) as demand_mw_peta_heat
359
            FROM demand.egon_peta_heat
360
            WHERE scenario= 'eGon2035'
361
            GROUP BY (scenario);
362
        """,
363
        warning=False,
364
    )["demand_mw_peta_heat"].values[0]
365
366
    e_demand = (
367
        round((output_heat_demand - input_heat_demand) / input_heat_demand, 2)
368
        * 100
369
    )
370
371
    logger.info(f"heat demand: {e_demand} %")
372
373
    # Sanity checks for heat supply
374
375
    logger.info(
376
        "For German heat supplies the following deviations between the inputs "
377
        "and outputs can be observed:"
378
    )
379
380
    # Comparison for central heat pumps
381
    heat_pump_input = db.select_dataframe(
382
        """SELECT carrier, SUM(capacity::numeric) as Urban_central_heat_pump_mw
383
            FROM supply.egon_scenario_capacities
384
            WHERE carrier= 'urban_central_heat_pump'
385
            AND scenario_name IN ('eGon2035')
386
            GROUP BY (carrier);
387
        """,
388
        warning=False,
389
    )["urban_central_heat_pump_mw"].values[0]
390
391
    heat_pump_output = db.select_dataframe(
392
        """SELECT carrier, SUM(p_nom::numeric) as Central_heat_pump_mw
393
            FROM grid.egon_etrago_link
394
            WHERE carrier= 'central_heat_pump'
395
            AND scn_name IN ('eGon2035')
396
            GROUP BY (carrier);
397
    """,
398
        warning=False,
399
    )["central_heat_pump_mw"].values[0]
400
401
    e_heat_pump = (
402
        round((heat_pump_output - heat_pump_input) / heat_pump_output, 2) * 100
403
    )
404
405
    logger.info(f"'central_heat_pump': {e_heat_pump} % ")
406
407
    # Comparison for residential heat pumps
408
409
    input_residential_heat_pump = db.select_dataframe(
410
        """SELECT carrier, SUM(capacity::numeric) as residential_heat_pump_mw
411
            FROM supply.egon_scenario_capacities
412
            WHERE carrier= 'residential_rural_heat_pump'
413
            AND scenario_name IN ('eGon2035')
414
            GROUP BY (carrier);
415
        """,
416
        warning=False,
417
    )["residential_heat_pump_mw"].values[0]
418
419
    output_residential_heat_pump = db.select_dataframe(
420
        """SELECT carrier, SUM(p_nom::numeric) as rural_heat_pump_mw
421
            FROM grid.egon_etrago_link
422
            WHERE carrier= 'rural_heat_pump'
423
            AND scn_name IN ('eGon2035')
424
            GROUP BY (carrier);
425
    """,
426
        warning=False,
427
    )["rural_heat_pump_mw"].values[0]
428
429
    e_residential_heat_pump = (
430
        round(
431
            (output_residential_heat_pump - input_residential_heat_pump)
432
            / input_residential_heat_pump,
433
            2,
434
        )
435
        * 100
436
    )
437
    logger.info(f"'residential heat pumps': {e_residential_heat_pump} %")
438
439
    # Comparison for resistive heater
440
    resistive_heater_input = db.select_dataframe(
441
        """SELECT carrier,
442
         SUM(capacity::numeric) as Urban_central_resistive_heater_MW
443
            FROM supply.egon_scenario_capacities
444
            WHERE carrier= 'urban_central_resistive_heater'
445
            AND scenario_name IN ('eGon2035')
446
            GROUP BY (carrier);
447
        """,
448
        warning=False,
449
    )["urban_central_resistive_heater_mw"].values[0]
450
451
    resistive_heater_output = db.select_dataframe(
452
        """SELECT carrier, SUM(p_nom::numeric) as central_resistive_heater_MW
453
            FROM grid.egon_etrago_link
454
            WHERE carrier= 'central_resistive_heater'
455
            AND scn_name IN ('eGon2035')
456
            GROUP BY (carrier);
457
        """,
458
        warning=False,
459
    )["central_resistive_heater_mw"].values[0]
460
461
    e_resistive_heater = (
462
        round(
463
            (resistive_heater_output - resistive_heater_input)
464
            / resistive_heater_input,
465
            2,
466
        )
467
        * 100
468
    )
469
470
    logger.info(f"'resistive heater': {e_resistive_heater} %")
471
472
    # Comparison for solar thermal collectors
473
474
    input_solar_thermal = db.select_dataframe(
475
        """SELECT carrier, SUM(capacity::numeric) as solar_thermal_collector_mw
476
            FROM supply.egon_scenario_capacities
477
            WHERE carrier= 'urban_central_solar_thermal_collector'
478
            AND scenario_name IN ('eGon2035')
479
            GROUP BY (carrier);
480
        """,
481
        warning=False,
482
    )["solar_thermal_collector_mw"].values[0]
483
484
    output_solar_thermal = db.select_dataframe(
485
        """SELECT carrier, SUM(p_nom::numeric) as solar_thermal_collector_mw
486
            FROM grid.egon_etrago_generator
487
            WHERE carrier= 'solar_thermal_collector'
488
            AND scn_name IN ('eGon2035')
489
            GROUP BY (carrier);
490
        """,
491
        warning=False,
492
    )["solar_thermal_collector_mw"].values[0]
493
494
    e_solar_thermal = (
495
        round(
496
            (output_solar_thermal - input_solar_thermal) / input_solar_thermal,
497
            2,
498
        )
499
        * 100
500
    )
501
    logger.info(f"'solar thermal collector': {e_solar_thermal} %")
502
503
    # Comparison for geothermal
504
505
    input_geo_thermal = db.select_dataframe(
506
        """SELECT carrier,
507
         SUM(capacity::numeric) as Urban_central_geo_thermal_MW
508
            FROM supply.egon_scenario_capacities
509
            WHERE carrier= 'urban_central_geo_thermal'
510
            AND scenario_name IN ('eGon2035')
511
            GROUP BY (carrier);
512
        """,
513
        warning=False,
514
    )["urban_central_geo_thermal_mw"].values[0]
515
516
    output_geo_thermal = db.select_dataframe(
517
        """SELECT carrier, SUM(p_nom::numeric) as geo_thermal_MW
518
            FROM grid.egon_etrago_generator
519
            WHERE carrier= 'geo_thermal'
520
            AND scn_name IN ('eGon2035')
521
            GROUP BY (carrier);
522
    """,
523
        warning=False,
524
    )["geo_thermal_mw"].values[0]
525
526
    e_geo_thermal = (
527
        round((output_geo_thermal - input_geo_thermal) / input_geo_thermal, 2)
528
        * 100
529
    )
530
    logger.info(f"'geothermal': {e_geo_thermal} %")
531
532
533
def residential_electricity_annual_sum(rtol=1e-5):
534
    """Sanity check for dataset electricity_demand_timeseries :
535
    Demand_Building_Assignment
536
537
    Aggregate the annual demand of all census cells at NUTS3 to compare
538
    with initial scaling parameters from DemandRegio.
539
    """
540
541
    df_nuts3_annual_sum = db.select_dataframe(
542
        sql="""
543
        SELECT dr.nuts3, dr.scenario, dr.demand_regio_sum, profiles.profile_sum
544
        FROM (
545
            SELECT scenario, SUM(demand) AS profile_sum, vg250_nuts3
546
            FROM demand.egon_demandregio_zensus_electricity AS egon,
547
             boundaries.egon_map_zensus_vg250 AS boundaries
548
            Where egon.zensus_population_id = boundaries.zensus_population_id
549
            AND sector = 'residential'
550
            GROUP BY vg250_nuts3, scenario
551
            ) AS profiles
552
        JOIN (
553
            SELECT nuts3, scenario, sum(demand) AS demand_regio_sum
554
            FROM demand.egon_demandregio_hh
555
            GROUP BY year, scenario, nuts3
556
              ) AS dr
557
        ON profiles.vg250_nuts3 = dr.nuts3 and profiles.scenario  = dr.scenario
558
        """
559
    )
560
561
    np.testing.assert_allclose(
562
        actual=df_nuts3_annual_sum["profile_sum"],
563
        desired=df_nuts3_annual_sum["demand_regio_sum"],
564
        rtol=rtol,
565
        verbose=False,
566
    )
567
568
    logger.info(
569
        "Aggregated annual residential electricity demand"
570
        " matches with DemandRegio at NUTS-3."
571
    )
572
573
574
def residential_electricity_hh_refinement(rtol=1e-5):
575
    """Sanity check for dataset electricity_demand_timeseries :
576
    Household Demands
577
578
    Check sum of aggregated household types after refinement method
579
    was applied and compare it to the original census values."""
580
581
    df_refinement = db.select_dataframe(
582
        sql="""
583
        SELECT refined.nuts3, refined.characteristics_code,
584
                refined.sum_refined::int, census.sum_census::int
585
        FROM(
586
            SELECT nuts3, characteristics_code, SUM(hh_10types) as sum_refined
587
            FROM society.egon_destatis_zensus_household_per_ha_refined
588
            GROUP BY nuts3, characteristics_code)
589
            AS refined
590
        JOIN(
591
            SELECT t.nuts3, t.characteristics_code, sum(orig) as sum_census
592
            FROM(
593
                SELECT nuts3, cell_id, characteristics_code,
594
                        sum(DISTINCT(hh_5types))as orig
595
                FROM society.egon_destatis_zensus_household_per_ha_refined
596
                GROUP BY cell_id, characteristics_code, nuts3) AS t
597
            GROUP BY t.nuts3, t.characteristics_code    ) AS census
598
        ON refined.nuts3 = census.nuts3
599
        AND refined.characteristics_code = census.characteristics_code
600
    """
601
    )
602
603
    np.testing.assert_allclose(
604
        actual=df_refinement["sum_refined"],
605
        desired=df_refinement["sum_census"],
606
        rtol=rtol,
607
        verbose=False,
608
    )
609
610
    logger.info("All Aggregated household types match at NUTS-3.")
611
612
613
def cts_electricity_demand_share(rtol=1e-5):
614
    """Sanity check for dataset electricity_demand_timeseries :
615
    CtsBuildings
616
617
    Check sum of aggregated cts electricity demand share which equals to one
618
    for every substation as the substation profile is linearly disaggregated
619
    to all buildings."""
620
621
    with db.session_scope() as session:
622
        cells_query = session.query(EgonCtsElectricityDemandBuildingShare)
623
624
    df_demand_share = pd.read_sql(
625
        cells_query.statement, cells_query.session.bind, index_col=None
626
    )
627
628
    np.testing.assert_allclose(
629
        actual=df_demand_share.groupby(["bus_id", "scenario"])[
630
            "profile_share"
631
        ].sum(),
632
        desired=1,
633
        rtol=rtol,
634
        verbose=False,
635
    )
636
637
    logger.info("The aggregated demand shares equal to one!.")
638
639
640
def cts_heat_demand_share(rtol=1e-5):
641
    """Sanity check for dataset electricity_demand_timeseries
642
    : CtsBuildings
643
644
    Check sum of aggregated cts heat demand share which equals to one
645
    for every substation as the substation profile is linearly disaggregated
646
    to all buildings."""
647
648
    with db.session_scope() as session:
649
        cells_query = session.query(EgonCtsHeatDemandBuildingShare)
650
651
    df_demand_share = pd.read_sql(
652
        cells_query.statement, cells_query.session.bind, index_col=None
653
    )
654
655
    np.testing.assert_allclose(
656
        actual=df_demand_share.groupby(["bus_id", "scenario"])[
657
            "profile_share"
658
        ].sum(),
659
        desired=1,
660
        rtol=rtol,
661
        verbose=False,
662
    )
663
664
    logger.info("The aggregated demand shares equal to one!.")
665
666
667
def sanitycheck_emobility_mit():
668
    """Execute sanity checks for eMobility: motorized individual travel
669
670
    Checks data integrity for eGon2035, eGon2035_lowflex and eGon100RE scenario
671
    using assertions:
672
      1. Allocated EV numbers and EVs allocated to grid districts
673
      2. Trip data (original inout data from simBEV)
674
      3. Model data in eTraGo PF tables (grid.egon_etrago_*)
675
676
    Parameters
677
    ----------
678
    None
679
680
    Returns
681
    -------
682
    None
683
    """
684
685
    def check_ev_allocation():
686
        # Get target number for scenario
687
        ev_count_target = scenario_variation_parameters["ev_count"]
688
        print(f"  Target count: {str(ev_count_target)}")
689
690
        # Get allocated numbers
691
        ev_counts_dict = {}
692
        with db.session_scope() as session:
693
            for table, level in zip(
694
                [
695
                    EgonEvCountMvGridDistrict,
696
                    EgonEvCountMunicipality,
697
                    EgonEvCountRegistrationDistrict,
698
                ],
699
                ["Grid District", "Municipality", "Registration District"],
700
            ):
701
                query = session.query(
702
                    func.sum(
703
                        table.bev_mini
704
                        + table.bev_medium
705
                        + table.bev_luxury
706
                        + table.phev_mini
707
                        + table.phev_medium
708
                        + table.phev_luxury
709
                    ).label("ev_count")
710
                ).filter(
711
                    table.scenario == scenario_name,
712
                    table.scenario_variation == scenario_var_name,
713
                )
714
715
                ev_counts = pd.read_sql(
716
                    query.statement, query.session.bind, index_col=None
717
                )
718
                ev_counts_dict[level] = ev_counts.iloc[0].ev_count
719
                print(
720
                    f"    Count table: Total count for level {level} "
721
                    f"(table: {table.__table__}): "
722
                    f"{str(ev_counts_dict[level])}"
723
                )
724
725
        # Compare with scenario target (only if not in testmode)
726
        if TESTMODE_OFF:
727
            for level, count in ev_counts_dict.items():
728
                np.testing.assert_allclose(
729
                    count,
730
                    ev_count_target,
731
                    rtol=0.0001,
732
                    err_msg=f"EV numbers in {level} seems to be flawed.",
733
                )
734
        else:
735
            print("    Testmode is on, skipping sanity check...")
736
737
        # Get allocated EVs in grid districts
738
        with db.session_scope() as session:
739
            query = session.query(
740
                func.count(EgonEvMvGridDistrict.egon_ev_pool_ev_id).label(
741
                    "ev_count"
742
                ),
743
            ).filter(
744
                EgonEvMvGridDistrict.scenario == scenario_name,
745
                EgonEvMvGridDistrict.scenario_variation == scenario_var_name,
746
            )
747
        ev_count_alloc = (
748
            pd.read_sql(query.statement, query.session.bind, index_col=None)
749
            .iloc[0]
750
            .ev_count
751
        )
752
        print(
753
            f"    EVs allocated to Grid Districts "
754
            f"(table: {EgonEvMvGridDistrict.__table__}) total count: "
755
            f"{str(ev_count_alloc)}"
756
        )
757
758
        # Compare with scenario target (only if not in testmode)
759
        if TESTMODE_OFF:
760
            np.testing.assert_allclose(
761
                ev_count_alloc,
762
                ev_count_target,
763
                rtol=0.0001,
764
                err_msg=(
765
                    "EV numbers allocated to Grid Districts seems to be flawed."
766
                ),
767
            )
768
        else:
769
            print("    Testmode is on, skipping sanity check...")
770
771
        return ev_count_alloc
772
773
    def check_trip_data():
774
        # Check if trips start at timestep 0 and have a max. of 35040 steps
775
        # (8760h in 15min steps)
776
        print("  Checking timeranges...")
777
        with db.session_scope() as session:
778
            query = session.query(
779
                func.count(EgonEvTrip.event_id).label("cnt")
780
            ).filter(
781
                or_(
782
                    and_(
783
                        EgonEvTrip.park_start > 0,
784
                        EgonEvTrip.simbev_event_id == 0,
785
                    ),
786
                    EgonEvTrip.park_end
787
                    > (60 / int(meta_run_config.stepsize)) * 8760,
788
                ),
789
                EgonEvTrip.scenario == scenario_name,
790
            )
791
        invalid_trips = pd.read_sql(
792
            query.statement, query.session.bind, index_col=None
793
        )
794
        np.testing.assert_equal(
795
            invalid_trips.iloc[0].cnt,
796
            0,
797
            err_msg=(
798
                f"{str(invalid_trips.iloc[0].cnt)} trips in table "
799
                f"{EgonEvTrip.__table__} have invalid timesteps."
800
            ),
801
        )
802
803
        # Check if charging demand can be covered by available charging energy
804
        # while parking
805
        print("  Compare charging demand with available power...")
806
        with db.session_scope() as session:
807
            query = session.query(
808
                func.count(EgonEvTrip.event_id).label("cnt")
809
            ).filter(
810
                func.round(
811
                    cast(
812
                        (EgonEvTrip.park_end - EgonEvTrip.park_start + 1)
813
                        * EgonEvTrip.charging_capacity_nominal
814
                        * (int(meta_run_config.stepsize) / 60),
815
                        Numeric,
816
                    ),
817
                    3,
818
                )
819
                < cast(EgonEvTrip.charging_demand, Numeric),
820
                EgonEvTrip.scenario == scenario_name,
821
            )
822
        invalid_trips = pd.read_sql(
823
            query.statement, query.session.bind, index_col=None
824
        )
825
        np.testing.assert_equal(
826
            invalid_trips.iloc[0].cnt,
827
            0,
828
            err_msg=(
829
                f"In {str(invalid_trips.iloc[0].cnt)} trips (table: "
830
                f"{EgonEvTrip.__table__}) the charging demand cannot be "
831
                f"covered by available charging power."
832
            ),
833
        )
834
835
    def check_model_data():
836
        # Check if model components were fully created
837
        print("  Check if all model components were created...")
838
        # Get MVGDs which got EV allocated
839
        with db.session_scope() as session:
840
            query = (
841
                session.query(
842
                    EgonEvMvGridDistrict.bus_id,
843
                )
844
                .filter(
845
                    EgonEvMvGridDistrict.scenario == scenario_name,
846
                    EgonEvMvGridDistrict.scenario_variation
847
                    == scenario_var_name,
848
                )
849
                .group_by(EgonEvMvGridDistrict.bus_id)
850
            )
851
        mvgds_with_ev = (
852
            pd.read_sql(query.statement, query.session.bind, index_col=None)
853
            .bus_id.sort_values()
854
            .to_list()
855
        )
856
857
        # Load model components
858
        with db.session_scope() as session:
859
            query = (
860
                session.query(
861
                    EgonPfHvLink.bus0.label("mvgd_bus_id"),
862
                    EgonPfHvLoad.bus.label("emob_bus_id"),
863
                    EgonPfHvLoad.load_id.label("load_id"),
864
                    EgonPfHvStore.store_id.label("store_id"),
865
                )
866
                .select_from(EgonPfHvLoad, EgonPfHvStore)
867
                .join(
868
                    EgonPfHvLoadTimeseries,
869
                    EgonPfHvLoadTimeseries.load_id == EgonPfHvLoad.load_id,
870
                )
871
                .join(
872
                    EgonPfHvStoreTimeseries,
873
                    EgonPfHvStoreTimeseries.store_id == EgonPfHvStore.store_id,
874
                )
875
                .filter(
876
                    EgonPfHvLoad.carrier == "land transport EV",
877
                    EgonPfHvLoad.scn_name == scenario_name,
878
                    EgonPfHvLoadTimeseries.scn_name == scenario_name,
879
                    EgonPfHvStore.carrier == "battery storage",
880
                    EgonPfHvStore.scn_name == scenario_name,
881
                    EgonPfHvStoreTimeseries.scn_name == scenario_name,
882
                    EgonPfHvLink.scn_name == scenario_name,
883
                    EgonPfHvLink.bus1 == EgonPfHvLoad.bus,
884
                    EgonPfHvLink.bus1 == EgonPfHvStore.bus,
885
                )
886
            )
887
        model_components = pd.read_sql(
888
            query.statement, query.session.bind, index_col=None
889
        )
890
891
        # Check number of buses with model components connected
892
        mvgd_buses_with_ev = model_components.loc[
893
            model_components.mvgd_bus_id.isin(mvgds_with_ev)
894
        ]
895
        np.testing.assert_equal(
896
            len(mvgds_with_ev),
897
            len(mvgd_buses_with_ev),
898
            err_msg=(
899
                f"Number of Grid Districts with connected model components "
900
                f"({str(len(mvgd_buses_with_ev))} in tables egon_etrago_*) "
901
                f"differ from number of Grid Districts that got EVs "
902
                f"allocated ({len(mvgds_with_ev)} in table "
903
                f"{EgonEvMvGridDistrict.__table__})."
904
            ),
905
        )
906
907
        # Check if all required components exist (if no id is NaN)
908
        np.testing.assert_equal(
909
            model_components.drop_duplicates().isna().any().any(),
910
            False,
911
            err_msg=(
912
                f"Some components are missing (see True values): "
913
                f"{model_components.drop_duplicates().isna().any()}"
914
            ),
915
        )
916
917
        # Get all model timeseries
918
        print("  Loading model timeseries...")
919
        # Get all model timeseries
920
        model_ts_dict = {
921
            "Load": {
922
                "carrier": "land transport EV",
923
                "table": EgonPfHvLoad,
924
                "table_ts": EgonPfHvLoadTimeseries,
925
                "column_id": "load_id",
926
                "columns_ts": ["p_set"],
927
                "ts": None,
928
            },
929
            "Link": {
930
                "carrier": "BEV charger",
931
                "table": EgonPfHvLink,
932
                "table_ts": EgonPfHvLinkTimeseries,
933
                "column_id": "link_id",
934
                "columns_ts": ["p_max_pu"],
935
                "ts": None,
936
            },
937
            "Store": {
938
                "carrier": "battery storage",
939
                "table": EgonPfHvStore,
940
                "table_ts": EgonPfHvStoreTimeseries,
941
                "column_id": "store_id",
942
                "columns_ts": ["e_min_pu", "e_max_pu"],
943
                "ts": None,
944
            },
945
        }
946
947
        with db.session_scope() as session:
948
            for node, attrs in model_ts_dict.items():
949
                print(f"    Loading {node} timeseries...")
950
                subquery = (
951
                    session.query(
952
                        getattr(attrs["table"], attrs["column_id"])
953
                    )
954
                    .filter(attrs["table"].carrier == attrs["carrier"])
955
                    .filter(attrs["table"].scn_name == scenario_name)
956
                    .subquery()
957
                )
958
959
                cols = [
960
                    getattr(attrs["table_ts"], c) for c in attrs["columns_ts"]
961
                ]
962
                query = session.query(
963
                    getattr(attrs["table_ts"], attrs["column_id"]), *cols
964
                ).filter(
965
                    getattr(attrs["table_ts"], attrs["column_id"]).in_(
966
                        subquery
967
                    ),
968
                    attrs["table_ts"].scn_name == scenario_name,
969
                )
970
                attrs["ts"] = pd.read_sql(
971
                    query.statement,
972
                    query.session.bind,
973
                    index_col=attrs["column_id"],
974
                )
975
976
        # Check if all timeseries have 8760 steps
977
        print("    Checking timeranges...")
978
        for node, attrs in model_ts_dict.items():
979
            for col in attrs["columns_ts"]:
980
                ts = attrs["ts"]
981
                invalid_ts = ts.loc[ts[col].apply(lambda _: len(_)) != 8760][
982
                    col
983
                ].apply(len)
984
                np.testing.assert_equal(
985
                    len(invalid_ts),
986
                    0,
987
                    err_msg=(
988
                        f"{str(len(invalid_ts))} rows in timeseries do not "
989
                        f"have 8760 timesteps. Table: "
990
                        f"{attrs['table_ts'].__table__}, Column: {col}, IDs: "
991
                        f"{str(list(invalid_ts.index))}"
992
                    ),
993
                )
994
995
        # Compare total energy demand in model with some approximate values
996
        # (per EV: 14,000 km/a, 0.17 kWh/km)
997
        print("  Checking energy demand in model...")
998
        total_energy_model = (
999
            model_ts_dict["Load"]["ts"].p_set.apply(lambda _: sum(_)).sum()
1000
            / 1e6
1001
        )
1002
        print(f"    Total energy amount in model: {total_energy_model} TWh")
1003
        total_energy_scenario_approx = ev_count_alloc * 14000 * 0.17 / 1e9
1004
        print(
1005
            f"    Total approximated energy amount in scenario: "
1006
            f"{total_energy_scenario_approx} TWh"
1007
        )
1008
        np.testing.assert_allclose(
1009
            total_energy_model,
1010
            total_energy_scenario_approx,
1011
            rtol=0.1,
1012
            err_msg=(
1013
                "The total energy amount in the model deviates heavily "
1014
                "from the approximated value for current scenario."
1015
            ),
1016
        )
1017
1018
        # Compare total storage capacity
1019
        print("  Checking storage capacity...")
1020
        # Load storage capacities from model
1021
        with db.session_scope() as session:
1022
            query = session.query(
1023
                func.sum(EgonPfHvStore.e_nom).label("e_nom")
1024
            ).filter(
1025
                EgonPfHvStore.scn_name == scenario_name,
1026
                EgonPfHvStore.carrier == "battery storage",
1027
            )
1028
        storage_capacity_model = (
1029
            pd.read_sql(
1030
                query.statement, query.session.bind, index_col=None
1031
            ).e_nom.sum()
1032
            / 1e3
1033
        )
1034
        print(
1035
            f"    Total storage capacity ({EgonPfHvStore.__table__}): "
1036
            f"{round(storage_capacity_model, 1)} GWh"
1037
        )
1038
1039
        # Load occurences of each EV
1040
        with db.session_scope() as session:
1041
            query = (
1042
                session.query(
1043
                    EgonEvMvGridDistrict.bus_id,
1044
                    EgonEvPool.type,
1045
                    func.count(EgonEvMvGridDistrict.egon_ev_pool_ev_id).label(
1046
                        "count"
1047
                    ),
1048
                )
1049
                .join(
1050
                    EgonEvPool,
1051
                    EgonEvPool.ev_id
1052
                    == EgonEvMvGridDistrict.egon_ev_pool_ev_id,
1053
                )
1054
                .filter(
1055
                    EgonEvMvGridDistrict.scenario == scenario_name,
1056
                    EgonEvMvGridDistrict.scenario_variation
1057
                    == scenario_var_name,
1058
                    EgonEvPool.scenario == scenario_name,
1059
                )
1060
                .group_by(EgonEvMvGridDistrict.bus_id, EgonEvPool.type)
1061
            )
1062
        count_per_ev_all = pd.read_sql(
1063
            query.statement, query.session.bind, index_col="bus_id"
1064
        )
1065
        count_per_ev_all["bat_cap"] = count_per_ev_all.type.map(
1066
            meta_tech_data.battery_capacity
1067
        )
1068
        count_per_ev_all["bat_cap_total_MWh"] = (
1069
            count_per_ev_all["count"] * count_per_ev_all.bat_cap / 1e3
1070
        )
1071
        storage_capacity_simbev = count_per_ev_all.bat_cap_total_MWh.div(
1072
            1e3
1073
        ).sum()
1074
        print(
1075
            f"    Total storage capacity (simBEV): "
1076
            f"{round(storage_capacity_simbev, 1)} GWh"
1077
        )
1078
1079
        np.testing.assert_allclose(
1080
            storage_capacity_model,
1081
            storage_capacity_simbev,
1082
            rtol=0.01,
1083
            err_msg=(
1084
                "The total storage capacity in the model deviates heavily "
1085
                "from the input data provided by simBEV for current scenario."
1086
            ),
1087
        )
1088
1089
        # Check SoC storage constraint: e_min_pu < e_max_pu for all timesteps
1090
        print("  Validating SoC constraints...")
1091
        stores_with_invalid_soc = []
1092
        for idx, row in model_ts_dict["Store"]["ts"].iterrows():
1093
            ts = row[["e_min_pu", "e_max_pu"]]
1094
            x = np.array(ts.e_min_pu) > np.array(ts.e_max_pu)
1095
            if x.any():
1096
                stores_with_invalid_soc.append(idx)
1097
1098
        np.testing.assert_equal(
1099
            len(stores_with_invalid_soc),
1100
            0,
1101
            err_msg=(
1102
                f"The store constraint e_min_pu < e_max_pu does not apply "
1103
                f"for some storages in {EgonPfHvStoreTimeseries.__table__}. "
1104
                f"Invalid store_ids: {stores_with_invalid_soc}"
1105
            ),
1106
        )
1107
1108
    def check_model_data_lowflex_eGon2035():
1109
        # TODO: Add eGon100RE_lowflex
1110
        print("")
1111
        print("SCENARIO: eGon2035_lowflex")
1112
1113
        # Compare driving load and charging load
1114
        print("  Loading eGon2035 model timeseries: driving load...")
1115
        with db.session_scope() as session:
1116
            query = (
1117
                session.query(
1118
                    EgonPfHvLoad.load_id,
1119
                    EgonPfHvLoadTimeseries.p_set,
1120
                )
1121
                .join(
1122
                    EgonPfHvLoadTimeseries,
1123
                    EgonPfHvLoadTimeseries.load_id == EgonPfHvLoad.load_id,
1124
                )
1125
                .filter(
1126
                    EgonPfHvLoad.carrier == "land transport EV",
1127
                    EgonPfHvLoad.scn_name == "eGon2035",
1128
                    EgonPfHvLoadTimeseries.scn_name == "eGon2035",
1129
                )
1130
            )
1131
        model_driving_load = pd.read_sql(
1132
            query.statement, query.session.bind, index_col=None
1133
        )
1134
        driving_load = np.array(model_driving_load.p_set.to_list()).sum(axis=0)
1135
1136
        print(
1137
            "  Loading eGon2035_lowflex model timeseries: dumb charging "
1138
            "load..."
1139
        )
1140
        with db.session_scope() as session:
1141
            query = (
1142
                session.query(
1143
                    EgonPfHvLoad.load_id,
1144
                    EgonPfHvLoadTimeseries.p_set,
1145
                )
1146
                .join(
1147
                    EgonPfHvLoadTimeseries,
1148
                    EgonPfHvLoadTimeseries.load_id == EgonPfHvLoad.load_id,
1149
                )
1150
                .filter(
1151
                    EgonPfHvLoad.carrier == "land transport EV",
1152
                    EgonPfHvLoad.scn_name == "eGon2035_lowflex",
1153
                    EgonPfHvLoadTimeseries.scn_name == "eGon2035_lowflex",
1154
                )
1155
            )
1156
        model_charging_load_lowflex = pd.read_sql(
1157
            query.statement, query.session.bind, index_col=None
1158
        )
1159
        charging_load = np.array(
1160
            model_charging_load_lowflex.p_set.to_list()
1161
        ).sum(axis=0)
1162
1163
        # Ratio of driving and charging load should be 0.9 due to charging
1164
        # efficiency
1165
        print("  Compare cumulative loads...")
1166
        print(f"    Driving load (eGon2035): {driving_load.sum() / 1e6} TWh")
1167
        print(
1168
            f"    Dumb charging load (eGon2035_lowflex): "
1169
            f"{charging_load.sum() / 1e6} TWh"
1170
        )
1171
        driving_load_theoretical = (
1172
            float(meta_run_config.eta_cp) * charging_load.sum()
1173
        )
1174
        np.testing.assert_allclose(
1175
            driving_load.sum(),
1176
            driving_load_theoretical,
1177
            rtol=0.01,
1178
            err_msg=(
1179
                f"The driving load (eGon2035) deviates by more than 1% "
1180
                f"from the theoretical driving load calculated from charging "
1181
                f"load (eGon2035_lowflex) with an efficiency of "
1182
                f"{float(meta_run_config.eta_cp)}."
1183
            ),
1184
        )
1185
1186
    print("=====================================================")
1187
    print("=== SANITY CHECKS FOR MOTORIZED INDIVIDUAL TRAVEL ===")
1188
    print("=====================================================")
1189
1190
    for scenario_name in ["eGon2035", "eGon100RE"]:
1191
        scenario_var_name = DATASET_CFG["scenario"]["variation"][scenario_name]
1192
1193
        print("")
1194
        print(f"SCENARIO: {scenario_name}, VARIATION: {scenario_var_name}")
1195
1196
        # Load scenario params for scenario and scenario variation
1197
        scenario_variation_parameters = get_sector_parameters(
1198
            "mobility", scenario=scenario_name
1199
        )["motorized_individual_travel"][scenario_var_name]
1200
1201
        # Load simBEV run config and tech data
1202
        meta_run_config = read_simbev_metadata_file(
1203
            scenario_name, "config"
1204
        ).loc["basic"]
1205
        meta_tech_data = read_simbev_metadata_file(scenario_name, "tech_data")
1206
1207
        print("")
1208
        print("Checking EV counts...")
1209
        ev_count_alloc = check_ev_allocation()
1210
1211
        print("")
1212
        print("Checking trip data...")
1213
        check_trip_data()
1214
1215
        print("")
1216
        print("Checking model data...")
1217
        check_model_data()
1218
1219
    print("")
1220
    check_model_data_lowflex_eGon2035()
1221
1222
    print("=====================================================")
1223
1224