11_tasks_solved.ipynb
Download Notebook

Assignment #11

In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mytoolbox import solar

Exercise #07-02: Energy trading

In this exercise you will manage the distribution of energy in a household setting: A user needs electric power that they can either get from a company, a PV module on their roof, or a battery. If they generate more power than they need, they can also sell the power back to the company.

To refresh your memory from Assignment #08: Energy generated by a PV module. The power \(P\) generated by a PV module is the irradiance G (W/m2) multiplied by the area of the module A (m2) and the efficiency of the module \(\eta\), which in turn depends on temperature and a variety of other things. Let’s simplify and use \(\eta = 0.2\) for this exercise. The energy E generated by the module is then \(E = \sum P \cdot \Delta t\), where \(\Delta t\) is the sampling interval.

Read solar_Dornbirn.csv and powerusage.csv.

  1. Inspect the datetime of the new spreadsheet. Concatenate the two data frames on their DatetimeIndex using pd.concat() and subset the resulting data frame to the first week of April.
In [2]:
sdo = pd.read_csv('../08_files/solar_Dornbirn.csv', parse_dates=True, index_col='datetime')
powerusage = pd.read_csv('powerusage.csv', sep=";", parse_dates=True, index_col='datetime')
powerusage
Pload
datetime
2023-04-01 00:00:00 491.144444
2023-04-01 01:00:00 491.144444
2023-04-01 02:00:00 491.144444
2023-04-01 03:00:00 491.144444
2023-04-01 04:00:00 491.144444
... ...
2023-04-07 19:00:00 2731.144445
2023-04-07 20:00:00 2730.522222
2023-04-07 21:00:00 491.144444
2023-04-07 22:00:00 491.144444
2023-04-07 23:00:00 491.144444

168 rows × 1 columns

In [3]:
powerusage.index.min(), powerusage.index.max()
(Timestamp('2023-04-01 00:00:00'), Timestamp('2023-04-07 23:00:00'))

It looks like powerusage.csv only provides data for the first week of April anyways..

In [4]:
# Concatenate the two data frames
sdo = pd.concat([sdo, powerusage], axis=1)

# Alternatively, using DataFrame.join():
# sdo = sdo.join(powerusage)
In [5]:
# Remember the following way of subsetting on a DatetimeIndex
sdo.loc["2023-04-01":"2023-04-01", "Pload"].plot(marker='o', figsize=(12, 4))
sdo = sdo["2023-04-01":"2023-04-07"]
# However, it only works if the index is monotonically increasing,
# and no datetimes are missing/double etc
#
# Add-on question: Why did we not use that way of subsetting the data frame WX in assignment #06?
# --> Because there were multiple stations with the same datetime label.

The working plot above allows me to quickly look at Pload from one day to get a better understanding of the new data.

In [6]:
sdo
tau declination h alpha iswr_clearsky Pload
datetime
2023-04-01 00:00:00 -180.0 4.413916 -38.173684 180.000000 0.0 491.144444
2023-04-01 01:00:00 -165.0 4.413916 -36.516711 198.728610 0.0 491.144444
2023-04-01 02:00:00 -150.0 4.413916 -31.847044 215.934758 0.0 491.144444
2023-04-01 03:00:00 -135.0 4.413916 -24.861708 230.988507 0.0 491.144444
2023-04-01 04:00:00 -120.0 4.413916 -16.301545 244.108658 0.0 491.144444
... ... ... ... ... ... ...
2023-04-07 19:00:00 105.0 6.764913 -5.002620 105.662482 0.0 2731.144445
2023-04-07 20:00:00 120.0 6.764913 -14.434683 117.373406 0.0 2730.522222
2023-04-07 21:00:00 135.0 6.764913 -22.858275 130.356790 0.0 491.144444
2023-04-07 22:00:00 150.0 6.764913 -29.686049 145.142932 0.0 491.144444
2023-04-07 23:00:00 165.0 6.764913 -34.220603 161.890853 0.0 491.144444

168 rows × 6 columns

In the following I will modify sdo quite a bit. Before I do so, I create a copy of the original data frame, so I can always re-start again easily.

In [7]:
sdo_orig = sdo.copy()
  1. Compute the hourly time series of power generated by a PV module facing south at an angle of \(50^\circ\) in Dornbirn during the first week of April under idealized clear-sky conditions. Use our module solar for your calculations. The PV module has an area of 10m2. Add the time series as a new column in your data frame.
In [8]:
# Calculate irradiance at inclined PV module 
# based on the idealized clear-sky irradiance 'iswr_clearsky', which is valid for a horizontal surface
beta = 50  # degrees incline
sdo['ics50'] = solar.comp_irradiance_incline(sdo["iswr_clearsky"], 
                                             solar.comp_cos_incidenceangle(sdo['h'], sdo['alpha'], beta),
                                             sdo['h'])
# Calculate generated power
area = 10  # m2
efficieny = 0.2
sdo['Pproduced'] = sdo['ics50'] * area * efficieny

Now we have computed a whole lot of irradiances. Let’s create a quick plot to check whether the results seem meaningful:

In [9]:
fig, ax = plt.subplots()
sdo['iswr_clearsky'][0:25].plot(ax=ax, label="ics0")
sdo['ics50'][0:25].plot(ax=ax, label="ics50")
ax.legend()
<matplotlib.legend.Legend at 0x7f23bcad0390>

It looks like the results make sense in the big picture: In April, we expect a higher maximum solar radiation at a surface that is inclined toward south at an angle as opposed to a horizontal surface.

In [10]:
# Compute net power (i.e., production - usage):
sdo['Pnet'] = sdo['Pproduced'] - sdo['Pload']
In [11]:
# Define characteristics of battery
cbat_max = 5 * 1000  # (Wh)
# Initialize empty time series of hourly battery capacity, power bought, and power sold
sdo['cbat'] = np.full(sdo.shape[0], np.nan)
sdo['Pbought'] = np.full(sdo.shape[0], 0.0)
sdo['Psold'] = np.full(sdo.shape[0], 0.0)
# Define the initial energy level in the battery: let's use 50%
sdo.loc[sdo.index[0], 'cbat'] = 0.5*cbat_max
In [12]:
sdo
tau declination h alpha iswr_clearsky Pload ics50 Pproduced Pnet cbat Pbought Psold
datetime
2023-04-01 00:00:00 -180.0 4.413916 -38.173684 180.000000 0.0 491.144444 -0.0 -0.0 -491.144444 2500.0 0.0 0.0
2023-04-01 01:00:00 -165.0 4.413916 -36.516711 198.728610 0.0 491.144444 -0.0 -0.0 -491.144444 NaN 0.0 0.0
2023-04-01 02:00:00 -150.0 4.413916 -31.847044 215.934758 0.0 491.144444 -0.0 -0.0 -491.144444 NaN 0.0 0.0
2023-04-01 03:00:00 -135.0 4.413916 -24.861708 230.988507 0.0 491.144444 -0.0 -0.0 -491.144444 NaN 0.0 0.0
2023-04-01 04:00:00 -120.0 4.413916 -16.301545 244.108658 0.0 491.144444 -0.0 -0.0 -491.144444 NaN 0.0 0.0
... ... ... ... ... ... ... ... ... ... ... ... ...
2023-04-07 19:00:00 105.0 6.764913 -5.002620 105.662482 0.0 2731.144445 -0.0 -0.0 -2731.144445 NaN 0.0 0.0
2023-04-07 20:00:00 120.0 6.764913 -14.434683 117.373406 0.0 2730.522222 -0.0 -0.0 -2730.522222 NaN 0.0 0.0
2023-04-07 21:00:00 135.0 6.764913 -22.858275 130.356790 0.0 491.144444 -0.0 -0.0 -491.144444 NaN 0.0 0.0
2023-04-07 22:00:00 150.0 6.764913 -29.686049 145.142932 0.0 491.144444 -0.0 -0.0 -491.144444 NaN 0.0 0.0
2023-04-07 23:00:00 165.0 6.764913 -34.220603 161.890853 0.0 491.144444 -0.0 -0.0 -491.144444 NaN 0.0 0.0

168 rows × 12 columns

  1. Calculate and visualize the hourly (i) power bought from the electricity company, (ii) power sold to the electricity company, (iii) percentage of battery energy level (max storage capacity of 5 kWh). Also include the power generated by the PV module and the power usage into your figure. Here are more necessary decision points to run the calculations:
    • If the battery level is less than 20%, recharge it up to 20% (but not more) using power generated by PV module (if it did produce) and do not use battery for the current time step
    • Using power from the PV module or battery is preferred to buying power
    • Storing power in the battery is preferred to selling power
In [13]:
# Iterate over hourly time series
dt = 1  # [hour] knowledge of time step is required to convert power and energy
for i, idx in enumerate(sdo.index):

    # In the beginning of iteration, battery level is put forward from last time step:
    if i > 0:
        sdo.loc[idx, 'cbat'] = sdo.loc[sdo.index[i-1], 'cbat']
    # Create boolean switch that enables use of battery for current time step
    usebat = True
        
    """If the battery level is less than 20%, recharge it up to 20% (but not more) 
    using power generated by PV module (if it did produce) and do not use battery
    for the current time step""" 
    if sdo.loc[idx, 'cbat'] < 0.2*cbat_max:
        # determine how much power to send to battery:
        p2bat = min([0.2*cbat_max-sdo.loc[idx, 'cbat'], 
                    sdo.loc[idx, 'Pproduced']])
        # update battery energy content:
        sdo.loc[idx, 'cbat'] = sdo.loc[idx, 'cbat'] + p2bat*dt  # (Wh)
        # update net power 
        sdo.loc[idx, 'Pnet'] = sdo.loc[idx, 'Pproduced'] - p2bat - sdo.loc[idx, 'Pload']
        # lock battery for iteration
        usebat = False

    """Determine power sold/bought based on net power"""
    if sdo.loc[idx, 'Pnet'] >= 0:  # More power was produced than user needed
        # determine how much power to send to battery
        p2bat = min([cbat_max-sdo.loc[idx, 'cbat'],
                    sdo.loc[idx, 'Pnet']])
        # update battery energy content
        sdo.loc[idx, 'cbat'] = sdo.loc[idx, 'cbat'] + p2bat*dt  # (Wh)
        # compute power available for selling
        sdo.loc[idx, 'Psold'] = sdo.loc[idx, 'Pnet'] - p2bat
    else:  # Less power was produced than user needed
        # determine how much power to use from battery
        if usebat:
            bat2user = min([sdo.loc[idx, 'cbat'],
                       abs(sdo.loc[idx, 'Pnet'])])
        else: 
            bat2user = 0
        
        # update battery energy content
        sdo.loc[idx, 'cbat'] = sdo.loc[idx, 'cbat'] - bat2user*dt  # (Wh)
        # compute power that needs to be bought to satisfy user need
        sdo.loc[idx, 'Pbought'] = abs(sdo.loc[idx, 'Pnet']) - bat2user
        
In [14]:
# Compute percentage of battery level
sdo['pctbat'] = sdo['cbat'] / cbat_max
In [15]:
sdo
tau declination h alpha iswr_clearsky Pload ics50 Pproduced Pnet cbat Pbought Psold pctbat
datetime
2023-04-01 00:00:00 -180.0 4.413916 -38.173684 180.000000 0.0 491.144444 -0.0 -0.0 -491.144444 2008.855556 0.000000 0.0 0.401771
2023-04-01 01:00:00 -165.0 4.413916 -36.516711 198.728610 0.0 491.144444 -0.0 -0.0 -491.144444 1517.711112 0.000000 0.0 0.303542
2023-04-01 02:00:00 -150.0 4.413916 -31.847044 215.934758 0.0 491.144444 -0.0 -0.0 -491.144444 1026.566667 0.000000 0.0 0.205313
2023-04-01 03:00:00 -135.0 4.413916 -24.861708 230.988507 0.0 491.144444 -0.0 -0.0 -491.144444 535.422223 0.000000 0.0 0.107084
2023-04-01 04:00:00 -120.0 4.413916 -16.301545 244.108658 0.0 491.144444 -0.0 -0.0 -491.144444 535.422223 491.144444 0.0 0.107084
... ... ... ... ... ... ... ... ... ... ... ... ... ...
2023-04-07 19:00:00 105.0 6.764913 -5.002620 105.662482 0.0 2731.144445 -0.0 -0.0 -2731.144445 0.000000 2731.144445 0.0 0.000000
2023-04-07 20:00:00 120.0 6.764913 -14.434683 117.373406 0.0 2730.522222 -0.0 -0.0 -2730.522222 0.000000 2730.522222 0.0 0.000000
2023-04-07 21:00:00 135.0 6.764913 -22.858275 130.356790 0.0 491.144444 -0.0 -0.0 -491.144444 0.000000 491.144444 0.0 0.000000
2023-04-07 22:00:00 150.0 6.764913 -29.686049 145.142932 0.0 491.144444 -0.0 -0.0 -491.144444 0.000000 491.144444 0.0 0.000000
2023-04-07 23:00:00 165.0 6.764913 -34.220603 161.890853 0.0 491.144444 -0.0 -0.0 -491.144444 0.000000 491.144444 0.0 0.000000

168 rows × 13 columns

Visualize the results. In the following, you will find a figure that I tried to style a lot. I don’t expect you to come up with this figure, but I want you to understand how I implement the styling so that you can come back here and use some of the features as template for your future figues.

In [16]:
## Plot figure and artists
# Main (left) axis
fig, ax = plt.subplots()
l1, = ax.plot(sdo['Psold']/1000, color="green")
l2, = ax.plot(-sdo['Pbought']/1000, color="red", alpha=0.9)
a1 = ax.fill_between(sdo.index, sdo['Pproduced']/1000, color="green", alpha=0.2)
a2 = ax.fill_between(sdo.index, -sdo['Pload']/1000, color="red", alpha=0.2)

# Add second axis on the right
ax2 = ax.twinx()
l3, = ax2.plot(sdo['pctbat'], color="black")

## Style plot
ax.set_ylim(-6, 2)
ax2.set_ylim(-3, 1)
ax2.legend([l1, l2, a1, a2], 
           ['Psold', 'Pbought', 
            'Pproduced', 'Pload'], 
           ncol=4)
ax.set_ylabel('Power (kW)')
ax2.set_ylabel('Battery level (%)')
# move the right ylabel up: 
# (check out the function documentation for figure-relative coordinates)
ax2.yaxis.set_label_coords(1.1, 0.87)
# Make xticks prettier:
# (by default with this plotting approach, the date labels are plotted horizontally
# and will overlap)
# first, specify which locations should be ticked:
ax.set_xticks(sdo.index[sdo.index.hour.isin([12])])
# second, specify the labels for the ticks
# (method strftime creates strings from the datetime objects with the specified datetime formatting)
ax.set_xticklabels(sdo.index[sdo.index.hour.isin([12])].strftime('%Y-%m-%d'), rotation=45, ha='right')
# Set the desired y-axis tick locations
ax2.set_yticks(np.arange(0, 1.1, 0.5))
ax.grid(axis='x', linestyle=':')

In the following you see a figure that is barely styled at all. I expect you to create those without external help.

In [17]:
fig, ax = plt.subplots()
sdo['Psold'].plot(ax=ax, label="Psold (W)", color='green')
sdo['Pbought'].plot(ax=ax, label="Pbought (W)", color='red')
sdo['cbat'].plot(ax=ax, label="Battery energy storage (Wh)", color='black')
# sdo['Pproduced'].plot.area(ax=ax, label="Pproduced (W)", alpha=0.2, color="green")
# sdo['Pload'].plot.area(ax=ax, label="Pload (W)", alpha=0.2, color='red')
ax.legend()
<matplotlib.legend.Legend at 0x7f23bc619dd0>

Exercise#07-03: Violin plots

Run the same calculations as above, but change the angle of the PV module to \(30^\circ\). Create a violinplot of the percentage of battery level for the two scenarios (i.e., one violin per angle of PV module). Use the matplotlib axes method .violinplot(). The function takes an array or a sequence of vectors as input data.

Since we already have the code to compute all the required data above, we can either copy and paste the relevant code cells and modify the characteristics of the PV module, or we could pack the computations from above in a function and then easily run multiple computations for various PV characteristics. It’s your choice. I will go with the function in this solution.

In [18]:
def assess_power_household(df, beta=45, area=10, cbat_max=5000, cbat_init=0.5, efficieny=0.2, verbose=True):
    """
    Assess power dynamics in a household with a solar PV system and battery storage.

    Parameters
    ----------
    - df (DataFrame): DataFrame with datetime index containing columns: 
                     'Pload' (power demand),
                     'iswr_clearsky' (clear-sky irradiance),
                     'h' (elevation angle of the sun),
                     'alpha' (azimuth angle of the sun).
    - beta (float, optional): Degrees of inclination for the PV module. Default is 45 degrees.
    - area (float, optional): Area of the PV module in square meters. Default is 10 square meters.
    - cbat_max (float, optional): Maximum capacity of the battery in Watt-hours (Wh). Default is 5000 Wh.
    - cbat_init (float, optional): Initial battery level as a percentage of maximum capacity. Default is 0.5 (50%).

    Returns:
    DataFrame: Updated DataFrame with additional columns:
               'Pproduced' (power generated by the PV system),
               'Pnet' (net power: production - consumption),
               'cbat' (battery energy level in Wh),
               'Pbought' (power bought from the grid),
               'Psold' (power sold to the grid),
               'pctbat' (percentage of battery level).

    Note:
    - The function calculates the power dynamics of a household with solar PV and a battery system.
    - It computes the irradiance at an inclined PV module, generates power, and assesses the net power usage.
    - The battery is used to store excess power and provide power during periods of high demand.
    - The function returns a DataFrame (with additional columns) reflecting the power dynamics over time.

    Usage:
    ```
    import pandas as pd
    import numpy as np
    from solar_module import comp_irradiance_incline, comp_cos_incidenceangle

    # Sample data creation
    date_rng = pd.date_range('2023-01-01', '2023-01-10', freq='D')
    sdo = pd.DataFrame({'Pload': np.random.rand(len(date_rng)),
                        'iswr_clearsky': np.random.rand(len(date_rng)),
                        'h': np.random.rand(len(date_rng)),
                        'alpha': np.random.rand(len(date_rng))}, index=date_rng)

    # Assess power dynamics
    result_df = assess_power_household(sdo)
    ```
    """

    df = df.copy()  # Note this statement! What would happen if we removed this line?
    if verbose:
        print(f"Characteristics set to {beta=}, {area=}, {cbat_max=}, {cbat_init=}")

    # beta degrees incline
    # cbat_max max capacity of battery (W)
    # cbat_init initial battery level (%)
    # df dataframe with datetime index containing columns: Pload, iswr_clearsky, h (elevation angle sun), alpha (azimuth sun)
    # Calculate irradiance at inclined PV module 
    # based on the idealized clear-sky irradiance 'iswr_clearsky', which is valid for a horizontal surface
    df['ics50'] = solar.comp_irradiance_incline(df["iswr_clearsky"], 
                                                 solar.comp_cos_incidenceangle(df['h'], df['alpha'], beta),
                                                 df['h'])
    # Calculate generated power
    df['Pproduced'] = df['ics50'] * area * efficieny

    # Compute net power (i.e., production - usage):
    df['Pnet'] = df['Pproduced'] - df['Pload']

    # Initialize empty time series of hourly battery capacity, power bought, and power sold
    df['cbat'] = np.full(df.shape[0], np.nan)
    df['Pbought'] = np.full(df.shape[0], 0.0)
    df['Psold'] = np.full(df.shape[0], 0.0)
    # Define the initial energy level in the battery: let's use 50%
    df.loc[df.index[0], 'cbat'] = cbat_init*cbat_max

    # Iterate over hourly time series
    dt = 1  # time step required to convert power and energy
    for i, idx in enumerate(df.index):
    
        # In the beginning of iteration, battery level is put forward from last time step:
        if i > 0:
            df.loc[idx, 'cbat'] = df.loc[df.index[i-1], 'cbat']
        # Create boolean switch that enables use of battery for current time step
        usebat = True
            
        """If the battery level is less than 20%, recharge it up to 20% (but not more) 
        using power generated by PV module (if it did produce) and do not use battery
        for the current time step""" 
        if df.loc[idx, 'cbat'] < 0.2*cbat_max:
            # determine how much power to send to battery:
            p2bat = min([0.2*cbat_max-df.loc[idx, 'cbat'], 
                        df.loc[idx, 'Pproduced']])
            # update battery energy content:
            df.loc[idx, 'cbat'] = df.loc[idx, 'cbat'] + p2bat*dt  # (Wh)
            # update net power 
            df.loc[idx, 'Pnet'] = df.loc[idx, 'Pproduced'] - p2bat - df.loc[idx, 'Pload']
            # lock battery for iteration
            usebat = False
    
        """Determine power sold/bought based on net power"""
        if df.loc[idx, 'Pnet'] >= 0:  # More power was produced than user needed
            # determine how much power to send to battery
            p2bat = min([cbat_max-df.loc[idx, 'cbat'],
                        df.loc[idx, 'Pnet']])
            # update battery energy content
            df.loc[idx, 'cbat'] = df.loc[idx, 'cbat'] + p2bat*dt  # (Wh)
            # compute power available for selling
            df.loc[idx, 'Psold'] = df.loc[idx, 'Pnet'] - p2bat
        else:  # Less power was produced than user needed
            # determine how much power to use from battery
            if usebat:
                bat2user = min([df.loc[idx, 'cbat'],
                           abs(df.loc[idx, 'Pnet'])])
            else: 
                bat2user = 0
            
            # update battery energy content
            df.loc[idx, 'cbat'] = df.loc[idx, 'cbat'] - bat2user*dt  # (Wh)
            # compute power that needs to be bought to satisfy user need
            df.loc[idx, 'Pbought'] = abs(df.loc[idx, 'Pnet']) - bat2user

    # Compute percentage of battery level
    df['pctbat'] = df['cbat'] / cbat_max

    return df

Let me do the same for the plotting function:

In [19]:
def plot_power_household_assessment(df):
    ## Plot figure and artists
    # Main (left) axis
    fig, ax = plt.subplots()
    l1, = ax.plot(df['Psold']/1000, color="green")
    l2, = ax.plot(-df['Pbought']/1000, color="red", alpha=0.9)
    a1 = ax.fill_between(df.index, df['Pproduced']/1000, color="green", alpha=0.2)
    a2 = ax.fill_between(df.index, -df['Pload']/1000, color="red", alpha=0.2)
    
    # Add second axis on the right
    ax2 = ax.twinx()
    l3, = ax2.plot(df['pctbat'], color="black")
    
    ## Style plot
    ax.set_ylim(-6, 5)
    ax2.set_ylim(-6/5, 1)
    ax2.legend([l1, l2, a1, a2], 
               ['Psold', 'Pbought', 
                'Pproduced', 'Pload'], 
               ncol=4)
    ax.set_ylabel('Power (kW)')
    ax2.set_ylabel('Battery level (%)')
    # move the right ylabel up: 
    # (check out the function documentation for figure-relative coordinates)
    ax2.yaxis.set_label_coords(1.08, 0.87)
    # Make xticks prettier:
    # (by default with this plotting approach, the date labels are plotted horizontally
    # and will overlap)
    # first, specify which locations should be ticked:
    ax.set_xticks(df.index[df.index.hour.isin([12])])
    # second, specify the labels for the ticks
    # (method strftime creates strings from the datetime objects with the specified datetime formatting)
    ax.set_xticklabels(df.index[df.index.hour.isin([12])].strftime('%Y-%m-%d'), rotation=45, ha='right')
    # Set the desired y-axis tick locations
    ax2.set_yticks(np.arange(0, 1.1, 0.5))
    ax.grid(axis='x', linestyle=':')
    
    return fig, ax, ax2

With these two functions, I can now compute and plot any scenario I can think of. Let’s compute for \(\beta=50\) and \(area=[10, 20]\):
(I know that these are different values than I told you in the exercise. Go and use the function yourself to see the results for the actual exercise!)

In [20]:
sdo_2 = assess_power_household(sdo_orig, beta=50, area=10)
fig, ax, ax2 = plot_power_household_assessment(sdo_2)
fig.set_figwidth(12)
ax.set_title(r"area = $10m^2$")
Characteristics set to beta=50, area=10, cbat_max=5000, cbat_init=0.5
Text(0.5, 1.0, 'area = $10m^2$')

In [21]:
sdo_3 = assess_power_household(sdo_orig, beta=50, area=20)
fig, ax, ax2 = plot_power_household_assessment(sdo_3)
fig.set_figwidth(12)
ax.set_title(r"area = $20m^2$")
Characteristics set to beta=50, area=20, cbat_max=5000, cbat_init=0.5
Text(0.5, 1.0, 'area = $20m^2$')

General note:

Instead of placing the functions in this notebook here, it would be way more organized and structured to put the functions in a module. One would then have to import the module and run only few lines of code in the notebook. If the functions are specific to this project (here: exercise), I would place them in a module in the same directory. If the functions were of greater use for several (potential) future projects, I would place the module in a package like mytoolbox.

This exercise here gives you plenty of training material for the upcoming exam. Be creative and try to think of meaningful analysis tasks with this data set. For example, the logic behind the battery charging cylces could be improved so that the battery does not rest with a charge of \(<20%\) or \(>80%\) for long. And many more questions..

You also noticed that the logic of the code gets quite complicated fast. This is only a single household, though, and our rules are not even sophisticated yet. You can imagine that things get a lot more complex fast when rules need to be better and more stakeholders participate in the network. One of your follow-up lectures will look into this topic from the perspective of optimization (whereas we implemented a rule-based approach).. Have fun!

Back to the exercise

Ok, back to our violin plot exercise. Here is the task again: Run the same calculations as above, but double the area of the PV module. Create a violinplot of the percentage of battery level for the two scenarios (i.e., one violin per area of PV module). Use the matplotlib axes method .violinplot(). The function takes an array or a sequence of vectors as input data.

I will use the two scenarios from above:

In [22]:
# Version 1, as easy as possible:
fig, ax = plt.subplots()
violins = ax.violinplot([sdo_2['pctbat'], sdo_3['pctbat']])

In [23]:
# Version 2, style the violins to your liking:
fig, ax = plt.subplots()
violins = ax.violinplot([sdo_2['pctbat'], sdo_3['pctbat']], showextrema=False, showmedians=True, quantiles=np.tile(np.array([0.25, 0.75]), (2, 1)))

# Set custom colors for each violin
mycols = ['orange', 'yellow']
for i, pc in enumerate(violins['bodies']):
    pc.set_facecolor(mycols[i])
    pc.set_edgecolor('black')

# Style labels
ax.set_xticks(np.arange(1, len(violins['bodies']) + 1))
ax.set_xticklabels(["area=10", "area=20"])
ax.set_ylabel("Battery power level (%)")
ax.set_title(r"$\alpha=50^\circ$")

plt.show()