Scenario analysis - creating a heatmap¶

The following sections extend the analysis presented in the previous tutorial. Instead of examining relationships between a single system parameter and a single system variable it is also worthwhile exploring how the joint interaction of two or more system parameters affects dispatch outcomes.

Building on the preceding tutorial, in addition to augmenting the Demand Forecast for South Australia, the Demand Forecast for Victoria will also be modified. The goal of the analysis is to investigate how jointly varying these parameters would likely impact power flow over the Heywood interconnector, and also energy prices in South Australia. Results will be presented in the form of heatmaps, allowing the data to be visually represented in three dimensions. Note the same methods can be used to examine relationships between other system parameters and variables.

Imports¶

In [1]:

import os
import pickle
import itertools

import requests
import xmltodict
import numpy as np
from jsonpath_ng.ext import parse

import matplotlib.pyplot as plt
import matplotlib.ticker as ticker

# Base URL endpoint for the Dispatch API
base_url = 'http://localhost:8000/'

import os import pickle import itertools import requests import xmltodict import numpy as np from jsonpath_ng.ext import parse import matplotlib.pyplot as plt import matplotlib.ticker as ticker # Base URL endpoint for the Dispatch API base_url = 'http://localhost:8000/'

Load case file¶

In [2]:

def convert_casefile(path_to_file):
    """Load a historical NEMDE case file and convert it to a dict"""

    # Read case file contents
    with open(path_to_file, 'r') as f:
        casefile = f.read()

    # Force these nodes to always return lists
    force_list = ('Trade', 'TradeTypePriceStructure',)

    return xmltodict.parse(casefile, force_list=force_list)


def load_casefile(case_id):
    """Load case file"""
    
    # Load case file and convert to JSON
    casefile = convert_casefile(f'../../data/NEMSPDOutputs_{case_id}00.loaded')
   
    return casefile

# Load case file
case_id = '20210401001'
casefile = load_casefile(case_id=case_id)

def convert_casefile(path_to_file): """Load a historical NEMDE case file and convert it to a dict""" # Read case file contents with open(path_to_file, 'r') as f: casefile = f.read() # Force these nodes to always return lists force_list = ('Trade', 'TradeTypePriceStructure',) return xmltodict.parse(casefile, force_list=force_list) def load_casefile(case_id): """Load case file""" # Load case file and convert to JSON casefile = convert_casefile(f'../../data/NEMSPDOutputs_{case_id}00.loaded') return casefile # Load case file case_id = '20210401001' casefile = load_casefile(case_id=case_id)

Approach¶

As in the previous tutorial the case file saved to disk will be treated as our canonical 'source of truth'. This 'base' case file will always be loaded prior to performing updates.

The following two functions are used to update Demand Forecast parameters within a case file.

In [3]:

def update_region_demand_forecast(casefile, region_id, value):
    """
    Update demand forecast for a given region
    
    Parameters
    ----------
    casefile : dict
        Case information
    
    region_id : str
        Region for which Demand Forecast should be modified
    
    value : float
        New value for Demand Forecast parameter
        
    Returns
    -------
    casefile : dict
        Case file with updated Demand Forecast parameter
    """
    
    # Construct JSON path expression
    expression = (f"NEMSPDCaseFile \
                  .NemSpdInputs \
                  .PeriodCollection \
                  .Period \
                  .RegionPeriodCollection \
                  .RegionPeriod[?(@RegionID=='{region_id}')] \
                  .@DF")

    jsonpath_expr = parse(expression)

    # Update value
    jsonpath_expr.update(casefile, value)
    
    return casefile


def update_casefile(case_id, region_1, value_1, region_2, value_2):
    """
    Update Demand Forecast parameters in two regions
    
    Parameters
    ----------
    case_id : str
        Case ID for base case file
    
    region_1 : str
        First NEM region for which Demand Forecast should be updated
    
    value_1 : float
        New Demand Forecast value for region_1
    
    region_2 : str
        Second NEM region for which Demand Forecast should be updated
    
    value_2 : float
        New Demand Forecast value for region_2
    
    Returns
    -------
    casefile : dict
        Updated case file
    """
    
    # Load 'base' case file
    casefile = load_casefile(case_id=case_id)
    
    # Update parameters in both regions
    casefile = update_region_demand_forecast(casefile=casefile, region_id=region_1, 
                                            value=value_1)

    casefile = update_region_demand_forecast(casefile=casefile, region_id=region_2, 
                                            value=value_2)

    return casefile

def update_region_demand_forecast(casefile, region_id, value): """ Update demand forecast for a given region Parameters ---------- casefile : dict Case information region_id : str Region for which Demand Forecast should be modified value : float New value for Demand Forecast parameter Returns ------- casefile : dict Case file with updated Demand Forecast parameter """ # Construct JSON path expression expression = (f"NEMSPDCaseFile \ .NemSpdInputs \ .PeriodCollection \ .Period \ .RegionPeriodCollection \ .RegionPeriod[?(@RegionID=='{region_id}')] \ .@DF") jsonpath_expr = parse(expression) # Update value jsonpath_expr.update(casefile, value) return casefile def update_casefile(case_id, region_1, value_1, region_2, value_2): """ Update Demand Forecast parameters in two regions Parameters ---------- case_id : str Case ID for base case file region_1 : str First NEM region for which Demand Forecast should be updated value_1 : float New Demand Forecast value for region_1 region_2 : str Second NEM region for which Demand Forecast should be updated value_2 : float New Demand Forecast value for region_2 Returns ------- casefile : dict Updated case file """ # Load 'base' case file casefile = load_casefile(case_id=case_id) # Update parameters in both regions casefile = update_region_demand_forecast(casefile=casefile, region_id=region_1, value=value_1) casefile = update_region_demand_forecast(casefile=casefile, region_id=region_2, value=value_2) return casefile

Run scenarios¶

A modified case file is submitted to the /solve endpoint for each combination of Demand Forecast parameters, with each solution stored alongside the parameters used to obtain it. Five values are used for each region, giving 25 scenarios in total. As each request waits for its solution this takes a few minutes to complete.

In [4]:

def solve_casefile(base_url, casefile):
    """Submit a case file to the Dispatch API and return the solution"""

    url = base_url + 'solve'
    response = requests.post(url=url, json=casefile)
    response.raise_for_status()

    return response.json()


def run_scenarios(base_url, case_id, region_1, region_2, demand_forecasts):
    """Solve a modified case file for each combination of Demand Forecast parameters"""
    
    values = [i for i in itertools.product(demand_forecasts, demand_forecasts)]
    
    scenarios = []
    for value_1, value_2 in values:
        # Update the base case file and solve the model
        casefile = update_casefile(case_id=case_id,
                                   region_1=region_1, value_1=value_1,
                                   region_2=region_2, value_2=value_2)

        solution = solve_casefile(base_url=base_url, casefile=casefile)

        # Store the solution alongside the parameters used in the analysis
        parameters = {
            'region_1': region_1, 'value_1': value_1,
            'region_2': region_2, 'value_2': value_2
        }
        
        scenarios.append({'parameters': parameters, 'solution': solution})

    return scenarios

# Run scenarios with different Demand Forecast parameters
demand_forecasts = [i for i in range(-200, 201, 100)]
scenarios = run_scenarios(base_url=base_url, case_id=case_id, region_1='SA1', region_2='VIC1',
                          demand_forecasts=demand_forecasts)

def solve_casefile(base_url, casefile): """Submit a case file to the Dispatch API and return the solution""" url = base_url + 'solve' response = requests.post(url=url, json=casefile) response.raise_for_status() return response.json() def run_scenarios(base_url, case_id, region_1, region_2, demand_forecasts): """Solve a modified case file for each combination of Demand Forecast parameters""" values = [i for i in itertools.product(demand_forecasts, demand_forecasts)] scenarios = [] for value_1, value_2 in values: # Update the base case file and solve the model casefile = update_casefile(case_id=case_id, region_1=region_1, value_1=value_1, region_2=region_2, value_2=value_2) solution = solve_casefile(base_url=base_url, casefile=casefile) # Store the solution alongside the parameters used in the analysis parameters = { 'region_1': region_1, 'value_1': value_1, 'region_2': region_2, 'value_2': value_2 } scenarios.append({'parameters': parameters, 'solution': solution}) return scenarios # Run scenarios with different Demand Forecast parameters demand_forecasts = [i for i in range(-200, 201, 100)] scenarios = run_scenarios(base_url=base_url, case_id=case_id, region_1='SA1', region_2='VIC1', demand_forecasts=demand_forecasts)

Save results¶

Optionally save results. It is good practice to save results once they have been obtained. Re-analysing results simply requires loading the saved results, rather than having to re-run the scenarios.

In [5]:

def save_results(scenarios):
    """Save results to notebook directory"""
    
    output_dir = os.path.abspath(os.path.join(os.path.curdir, 'output'))
    
    with open(os.path.join(output_dir, 'scenarios.pickle'), 'wb') as f:
        pickle.dump(scenarios, f)
        
        
def load_results():
    """Load results from notebook directory"""
    
    output_dir = os.path.abspath(os.path.join(os.path.curdir, 'output'))
    
    with open(os.path.join(output_dir, 'scenarios.pickle'), 'rb') as f:
        scenarios = pickle.load(f)
        
    return scenarios


# Save / load results if required
save_results(scenarios=scenarios)
# scenarios = load_results()

def save_results(scenarios): """Save results to notebook directory""" output_dir = os.path.abspath(os.path.join(os.path.curdir, 'output')) with open(os.path.join(output_dir, 'scenarios.pickle'), 'wb') as f: pickle.dump(scenarios, f) def load_results(): """Load results from notebook directory""" output_dir = os.path.abspath(os.path.join(os.path.curdir, 'output')) with open(os.path.join(output_dir, 'scenarios.pickle'), 'rb') as f: scenarios = pickle.load(f) return scenarios # Save / load results if required save_results(scenarios=scenarios) # scenarios = load_results()

Extract results¶

The following function searches a solution dictionary for a target value corresponding to a given JSON path expression.

In [6]:

def get_target_value(solution, expression):
    """Extract value from solution dictionary given a JSON path expression"""
    
    jsonpath_expr = parse(expression)   
    values = [match.value for match in jsonpath_expr.find(solution)]
    
    # Only one value should be returned
    if len(values) != 1:
        raise Exception("Should only return 1 value, encountered", len(values))

    return values[0]

def get_target_value(solution, expression): """Extract value from solution dictionary given a JSON path expression""" jsonpath_expr = parse(expression) values = [match.value for match in jsonpath_expr.find(solution)] # Only one value should be returned if len(values) != 1: raise Exception("Should only return 1 value, encountered", len(values)) return values[0]

Functions can now be used to extract target values (energy prices in South Australia and power flow over the Heywood interconnector) from each solution.

In [7]:

def get_sa_energy_price(solution):
    """Extract South Australia's energy price from solution dictionary"""
    
    expression = "RegionSolution[?(@RegionID='SA1')].@EnergyPrice"

    return get_target_value(solution=solution, expression=expression)


def get_v_sa_flow(solution):
    """Extract V-SA (Heywood) interconnector flow from solution dictionary"""
    
    expression = "InterconnectorSolution[?(@InterconnectorID='V-SA')].@Flow"

    return get_target_value(solution=solution, expression=expression)

def get_sa_energy_price(solution): """Extract South Australia's energy price from solution dictionary""" expression = "RegionSolution[?(@RegionID='SA1')].@EnergyPrice" return get_target_value(solution=solution, expression=expression) def get_v_sa_flow(solution): """Extract V-SA (Heywood) interconnector flow from solution dictionary""" expression = "InterconnectorSolution[?(@InterconnectorID='V-SA')].@Flow" return get_target_value(solution=solution, expression=expression)

Plotting¶

A heatmap can be used to visualise the results in three dimensions. The horizontal and vertical axes will display the Demand Forecast parameters used when running the model, while colour will indicate value of the variable under investigation.

In [8]:

def get_plot_data(scenarios, func):
    """
    Organise plot data
    
    Parameters
    ----------
    scenarios : list
        Parameters used when augmenting case files, and the corresponding solutions
    
    func : function
        Function used to extract a value from each solution
    
    Returns
    -------
    x : numpy array
        x-axis values used in plot
    
    y : numpy array
        y-axis values used in plot
    
    z : numpy array
        z-axis values used in plot
    """

    # Extract x, y, z components
    x = np.array([s.get('parameters').get('value_1') for s in scenarios])
    y = np.array([s.get('parameters').get('value_2') for s in scenarios])
    z = np.array([func(s.get('solution')) for s in scenarios])

    # Add shape information to array - required when plotting
    shape = (len(set(x)), len(set(y)))

    return x.reshape(shape), y.reshape(shape), z.reshape(shape)


def format_tick_locations(ax):
    """Adjust tick spacing"""
    
    ax.xaxis.set_major_locator(ticker.MultipleLocator(100))
    ax.xaxis.set_minor_locator(ticker.MultipleLocator(50))
    ax.yaxis.set_major_locator(ticker.MultipleLocator(100))
    ax.yaxis.set_minor_locator(ticker.MultipleLocator(50))
    
    return ax


def plot_v_sa_flow(scenarios, shading=None):
    """V-SA interconnector flow as a function of SA and VIC demand forecasts"""
    
    X, Y, Z = get_plot_data(scenarios=scenarios, func=get_v_sa_flow)

    # Create heatmap and add colour bar
    fig, ax = plt.subplots()
    mesh = ax.pcolormesh(X, Y, Z, shading=shading, cmap='inferno')
    cbar = plt.colorbar(mesh)
    
    # Format labels
    ax.set_title('Heywood power flow')
    ax.set_xlabel('SA Demand Forecast (MW)')
    ax.set_ylabel('VIC Demand Forecast (MW)')
    cbar.set_label('V-SA flow (MW)')
    
    ax = format_tick_locations(ax=ax)
    ax.grid(linestyle='--')

    plt.show()
    
    
def plot_sa_energy_price(scenarios, shading=None):
    """SA energy price as a function of SA and VIC demand forecasts"""
    
    X, Y, Z = get_plot_data(scenarios=scenarios, func=get_sa_energy_price)
    
    # Create heatmap and add colour bar
    fig, ax = plt.subplots()
    mesh = ax.pcolormesh(X, Y, Z, shading=shading, cmap='inferno')
    cbar = plt.colorbar(mesh)
    
    # Format labels
    ax.set_title('SA energy price')
    ax.set_xlabel('SA Demand Forecast (MW)')
    ax.set_ylabel('VIC Demand Forecast (MW)')
    cbar.set_label('SA energy price ($/MWh)')
    
    # Format ticks and add grid lines
    ax = format_tick_locations(ax=ax)
    ax.grid(linestyle='--')

    plt.show()


# Create plots
plot_v_sa_flow(scenarios=scenarios, shading='auto')
plot_sa_energy_price(scenarios=scenarios, shading='auto')

def get_plot_data(scenarios, func): """ Organise plot data Parameters ---------- scenarios : list Parameters used when augmenting case files, and the corresponding solutions func : function Function used to extract a value from each solution Returns ------- x : numpy array x-axis values used in plot y : numpy array y-axis values used in plot z : numpy array z-axis values used in plot """ # Extract x, y, z components x = np.array([s.get('parameters').get('value_1') for s in scenarios]) y = np.array([s.get('parameters').get('value_2') for s in scenarios]) z = np.array([func(s.get('solution')) for s in scenarios]) # Add shape information to array - required when plotting shape = (len(set(x)), len(set(y))) return x.reshape(shape), y.reshape(shape), z.reshape(shape) def format_tick_locations(ax): """Adjust tick spacing""" ax.xaxis.set_major_locator(ticker.MultipleLocator(100)) ax.xaxis.set_minor_locator(ticker.MultipleLocator(50)) ax.yaxis.set_major_locator(ticker.MultipleLocator(100)) ax.yaxis.set_minor_locator(ticker.MultipleLocator(50)) return ax def plot_v_sa_flow(scenarios, shading=None): """V-SA interconnector flow as a function of SA and VIC demand forecasts""" X, Y, Z = get_plot_data(scenarios=scenarios, func=get_v_sa_flow) # Create heatmap and add colour bar fig, ax = plt.subplots() mesh = ax.pcolormesh(X, Y, Z, shading=shading, cmap='inferno') cbar = plt.colorbar(mesh) # Format labels ax.set_title('Heywood power flow') ax.set_xlabel('SA Demand Forecast (MW)') ax.set_ylabel('VIC Demand Forecast (MW)') cbar.set_label('V-SA flow (MW)') ax = format_tick_locations(ax=ax) ax.grid(linestyle='--') plt.show() def plot_sa_energy_price(scenarios, shading=None): """SA energy price as a function of SA and VIC demand forecasts""" X, Y, Z = get_plot_data(scenarios=scenarios, func=get_sa_energy_price) # Create heatmap and add colour bar fig, ax = plt.subplots() mesh = ax.pcolormesh(X, Y, Z, shading=shading, cmap='inferno') cbar = plt.colorbar(mesh) # Format labels ax.set_title('SA energy price') ax.set_xlabel('SA Demand Forecast (MW)') ax.set_ylabel('VIC Demand Forecast (MW)') cbar.set_label('SA energy price ($/MWh)') # Format ticks and add grid lines ax = format_tick_locations(ax=ax) ax.grid(linestyle='--') plt.show() # Create plots plot_v_sa_flow(scenarios=scenarios, shading='auto') plot_sa_energy_price(scenarios=scenarios, shading='auto')

The resolution of these plots can be increased by providing a more granular grid over which to explore the solution space. The trade-off is solution time - more scenarios must be solved. Interpolation can also be used to smooth the gradient, making it easier to visualise the results.

In [9]:

# Create plots with interpolation
plot_v_sa_flow(scenarios=scenarios, shading='gouraud')
plot_sa_energy_price(scenarios=scenarios, shading='gouraud')

# Create plots with interpolation plot_v_sa_flow(scenarios=scenarios, shading='gouraud') plot_sa_energy_price(scenarios=scenarios, shading='gouraud')

Summary¶

The Dispatch API allows users to explore relationships between the joint interaction of system parameters and system variables. When two system parameters are examined heatmaps can be used to visualise results in three dimensions. This type of scenario analysis can easily be extended to different system parameters and variables, demonstrating the flexibility provided by the Dispatch API when exploring the likely solution space of dispatch outcomes.