<jupyter_start><jupyter_text>DoE - Modern DoE example: Metamodels for black-box systemsThe purpose of this notebook is to show you an example of how to apply some of the modern DoE techniques available through simulation, see [the associated presentation](doe_erigrid_dtu_summerschool_2018-08-28.pdf).<jupyter_code>%matplotlib notebook
import pandas as pd
import HEMS_sim
import sobol_seq # you likely need to install this first
import numpy as np
import matplotlib.pyplot as plt
from util import block_print, enable_print
# Dictionary with basic configuration of the simulation

basic_conf = {
    'ID':'00',
    'batt_storage_capacity':20,
    'batt_charge_capacity':5,
    'pv1_scaling':1,
    'controller_change_rate':0.5,
    'climate_conditions':'sunny',
    'season':'autumn',
    'random_weather':False,
    'stochastic':False, # this time we exclude the fluctuation (deterministic simulation)
    'noise_scale':1} 
# Scenario name which determines the name of the files we will be saving with the results
scenario_name = 'without_noise'

# Do deterministic space-filling sampling with Sobol sequence
samples = 15 # more samples would be better if you have the time ;-)
z = sobol_seq.i4_sobol_generate(2, samples)
z = np.array(z) * np.array([10, 5]) + np.array([15, 2]) # Sobol samples are between 0 and 1, here we scale
# Range for factor 1 (battery storage capacity): 15 to 25
# Range for factor 2 (battery charge capacity): 2 to 7

variations = {}
for i in range(samples):
    run_id = 'run_' + str(i)
    variations[run_id] = {'ID': str(i), 'batt_storage_capacity': z[i][0], 'batt_charge_capacity': z[i][1]}
              

# Merging of the basic configuration and the variations
recipes = {key: basic_conf.copy() for key in variations}
recipes = {key: {**recipes[key],**data} for key,data in variations.items()}

# These are the parameter combinations which were selected by the Sobol sequence
plt.figure()
plt.plot(*(z.T), '*')
plt.xlabel('Battery Storage Capacity [kWh]')
plt.ylabel('Battery Charge Capacity [kW]')
# Run the simulations (delete temp files first)
for recipe in recipes:
    print("Starting run {0}".format(recipes[recipe]['ID']))
    block_print() # Print statements are blocked due to notebook message length limitations
    HEMS_sim.run_simulation(scenario_name,recipes[recipe])
    enable_print()
# Store results and get response summaries
run_store.close()
run_store = pd.HDFStore('temp_files/runs_summary_{}.h5'.format(scenario_name))
summaries = [run_store[k] for k in run_store.keys()]
run_store.close()
summaries = pd.concat(summaries, axis=0).set_index('ID')
# Get the used treatments (battery size and power)
exog = summaries[['battery storage capacity [kWh]', 'battery charge capacity[kW]']]
# Get the self consumption index as response for each treatment
endog = summaries[['Self consumption index']]
# We want to fit a metamodel to our system
# We chose Kriging (Gaussian Process) this time around. You can also choose other metamodels if you want 
# (simplest would be linear interpolation)
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C
# Kriging needs a kernel. This kernel parameterization should work, but you can also play around with it.
kernel = C(1.0, (1e-3, 1e3)) * RBF(10, (1e-2, 1e2)) 
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9)
gp.fit(exog.values, endog.values) # Fit the Kriging model to our samples and associated responses
# Next we need to evaluate our metamodel. Thus, we set up a 50x50 point evaluation grid:
x_vec, y_vec = np.meshgrid(np.linspace(15, 25, 50), np.linspace(2, 7, 50))

# x_vec and y_vec is for plotting. We reshape them to get vectors for predicting:
evalgrid = np.array([x_vec.flatten(), y_vec.flatten()]).T

# Using our metamodel for predictions (Kriging also gives us a sense of uncertainty via the sigma):
scipred, sigma = gp.predict(evalgrid, return_std=True)

# For plotting we need to reshape our prediction array:
pltsci = np.reshape(scipred, (50, 50))
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib 
# This gives us the figure in an extra window
fig = plt.figure()
ax = fig.gca(projection='3d')
# Get 3D surface plot:
surf = ax.plot_surface(x_vec, y_vec, pltsci)
# Include our original sample points:
ax.plot(exog.values[:,0], exog.values[:,1], endog.values.ravel(), 'r*')
plt.show()

# Some questions to ask yourself in the end:
# - Does the metamodel make sense to you? If not, why?
# - Did you use enough samples?<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10311750599520383, 'fraction_numerical': 0.017440592980161325, 'mean_word_length': 5.150134048257373, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '511223', 'license': 'MIT', 'n_tokens_mistral': 1534, 'n_tokens_neox': 1453, 'n_words': 521, 'path': 'DoE-exercises/4_Modern_DoE.ipynb', 'repo_name': 'ERIGrid/NA4-Summerschool_Aug18', 'size': '4587'}