Examples

Source files are at https://github.com/pnnl/pecblocks/tree/master/examples/hwpv. Date files are at the Harvard Dataverse. Results from these examples are discussed in Draft Paper

To create their own examples, users will need to provide:

• A configuration file, as described in Schema

• A training dataset from laboratory tests or EMT simulations, as described in Data from EMT

Training a Model

The following script trains a model on Windows. It’s called with one argument, the root of the configuration file, e.g., call train.bat osg4, assuming that osg4_config.json has been created. The directory osg4 will be created if necessary, but not erased. The last line plots training losses to a saved PDF, but not to the screen so that a calling script may continue without user acknowledgement.

if not exist %1\ mkdir %1

python pv3_training.py %1_config.json
python loss_plot.py %1 -1

The followng script is for Linux and Mac OS X, called like ./train.sh osg4

#!/bin/bash
if [ ! -d "$1" ]; then
  echo "creating directory $1"
  mkdir $1
fi
python3 pv3_training.py $1_config.json
python3 loss_plot.py $1 -1

The following Python code supervises training a model using pecblocks.

• Line 12 disables the output of mean absolute error (MAE), because root mean square error (RMSE) is always more suitable for HWPV applications.

• Line 13 disables the plot of losses at the end of a training run. This plot waits for user input, which prevents a batch file from completing several training runs. To plot the losses afterward, invoke python loss_plot.py osg4

• Lines 15-30 open the required configuration file.

• Lines 32-38 comprise the main usage of pecblocks API

  • Line 32 initializes a model instance from the JSON configuration file

  • Line 33 loads the training data set, which may take several seconds or longer

  • Line 34 normalizes the training data to a range 0..1 on each channel, and saves the normalization factors

  • Line 35 prepares the model for training (and evaluation)

  • Line 36 trains new model coefficients, which may take several minutes or longer. Progress updates are written to the console.

  • Line 37 saves the most recent trained model coefficients. (Note: these coefficients are saved periodically during training. The set of coefficients with lowest fitting loss is also saved as needed in each epoch.)

  • Line 38 summarizes the training errors for each output channel, over all cases or events in the training data set. More detail is available in the next steps. The RMSE is most relevant to HWPV applications.

• Lines 40-59 provide summary output of the training run. This information is also availble later.

• Lines 61-73 plot the losses vs. training epoch. Use of loss_plot.py is now preferred.

# copyright 2021-2024 Battelle Memorial Institute
# supervises training of HW model for a three-phase inverter
#  arg1: relative path to the model configuration file

import numpy as np
import os
import sys
import json
import matplotlib.pyplot as plt
import pecblocks.pv3_poly as pv3_model

bWantMAE = False
bWantPlot = False

if __name__ == '__main__':
  if len(sys.argv) > 1:
    config_file = sys.argv[1]
    fp = open (config_file, 'r')
    cfg = json.load (fp)
    fp.close()
    data_path = cfg['data_path']
    model_folder = cfg['model_folder']
    model_root = cfg['model_root']
  else:
    print ('Usage: python pv3_training.py config.json')
    quit()

  print ('model_folder =', model_folder)
  print ('model_root =', model_root)
  print ('data_path =', data_path)

  model = pv3_model.pv3(training_config=config_file)
  model.loadTrainingData(data_path)
  model.applyAndSaveNormalization()
  model.initializeModelStructure()
  train_time, LOSS, VALID, SENS = model.trainModelCoefficients(bMAE=False)
  model.saveModelCoefficients()
  rmse, mae, case_rmse, case_mae = model.trainingErrors(bByCase=False)

  nlookback = 10
  recent_loss = LOSS[len(LOSS)-nlookback:]

  fp = open (os.path.join (model_folder, 'summary.txt'), 'w')
  jp = open (os.path.join (model_folder, 'normfacs.json'), 'r')
  normfacs = json.load (jp)
  jp.close()
  print ('Dataset Summary: (Mean=Offset, Range=Scale)', file=fp)
  print ('Column       Min       Max      Mean     Range', file=fp)
  for key, row in normfacs.items():
    print ('{:6s} {:9.3f} {:9.3f} {:9.3f} {:9.3f}'.format (key, row['min'], row['max'], row['offset'], row['scale']), file=fp)

  out_size = len(model.COL_Y)
  recent_loss = np.mean(recent_loss)
  valstr = ' '.join('{:.4f}'.format(rmse[j]) for j in range(out_size))
  print ('COL_Y', model.COL_Y, file=fp)
  print ('Train time: {:.2f}, Recent loss: {:.6f}, RMS Errors: {:s}'.format (train_time, recent_loss, valstr), file=fp)
  print ('COL_Y', model.COL_Y)
  print ('Train time: {:.2f}, Recent loss: {:.6f}, RMS Errors: {:s}'.format (train_time, recent_loss, valstr))
  if bWantMAE:
    valstr = ' '.join('{:.4f}'.format(mae[j]) for j in range(out_size))
    print ('                          MAE Errors: {:s}'.format (valstr), file=fp)
  fp.close()

  if bWantPlot:
    plt.figure()
    plt.title(model_root)
    plt.plot(np.log10(LOSS), label='Training Loss')
    plt.plot(np.log10(VALID), label='Validation Loss')
    if np.min(SENS) > 0.0:
      plt.plot(np.log10(SENS), label='Sensitivity Loss')
    plt.ylabel ('Log10')
    plt.xlabel ('Epoch')
    plt.legend()
    plt.grid(True)
    plt.savefig(os.path.join(model_folder, '{:s}_train_loss.pdf'.format(model_root)))
    plt.show()

Exporting a Model

The following script exports a trained model for simulation, and its accuracy metrics on Windows. It’s called with one argument, the root of the configuration file, e.g., call export.bat osg4, assuming that osg4_config.json has been created. Furthermore, the directory osg4 should have been created and populated by training the model.

python pv3_export.py %1_config.json > %1\metrics.txt
python pv3_metrics.py %1_config.json >> %1\metrics.txt

The followng script is for Linux and Mac OS X, called like ./export.sh osg4

#!/bin/bash
python3 pv3_export.py $1_config.json > $1/metrics.txt
python3 pv3_metrics.py $1_config.json >> $1/metrics.txt

The following Python code exports (post-processes) a trained model for simulation using pecblocks.

• Lines 10-27 open the required configuration file

• Lines 29-34 comprise the main usage of pecblocks API

  • Line 29 initializes a model instance from the JSON configuration file, same as for training.

  • Line 30 loads the normalization factors from training. To export the model, it’s not necessary to load the training data set.

  • Line 31 prepares the model for export

  • Line 32 loads the trained model coefficients from binary pkl files, which are not human-readable.

  • Line 33 exports the trained model coefficients in human-readable JSON file. During this process, and s-domain version of the model is prepared for simulation at variable time step. The JSON file is readable by other applications in a variety of languages, including Python, C/C++, and MATLAB.

  • Line 34 checks for stability of the exported H1 poles. If there are any warnings about unstable poles, do not use this model for simulation. Specifying the stable2ndx type for H1 in a re-trained model should avoid this issue.

# copyright 2021-2024 Battelle Memorial Institute
# exports a trained HW model and normalization coefficients to a single JSON file
#  arg1: relative path to the model configuration file

import os
import sys
import json
import pecblocks.pv3_poly as pv3_model

if __name__ == '__main__':
  if len(sys.argv) > 1:
    config_file = sys.argv[1]
    fp = open (config_file, 'r')
    cfg = json.load (fp)
    fp.close()
    model_folder = cfg['model_folder']
    model_root = cfg['model_root']
    data_path = cfg['data_path']
  else:
    print ('Usage: python pv3_export.py config.json')
    quit()
  export_path = os.path.join(model_folder,'{:s}_fhf.json'.format(model_root))
  if len(sys.argv) > 2:
    export_path = sys.argv[2]

  print ('Read Model from:', model_folder)
  print ('Export Model to:', export_path)

  model = pv3_model.pv3(training_config=config_file)
  model.loadNormalization()
  model.initializeModelStructure()
  model.loadModelCoefficients()
  model.exportModel(export_path)
  model.check_poles()

Model Metrics

The following Python code summarizes metrics of a trained model using pecblocks.

• Line 11 disables the output of mean absolute error (MAE), because root mean square error (RMSE) is always more suitable for HWPV applications.

• Lines 13-28 open the required configuration file.

• Lines 30-35 comprise the main usage of pecblocks API

  • Line 30 initializes a model instance from the JSON configuration file, same as for training or exporting.

  • Line 31 loads the training data set, which may take several seconds or longer

  • Line 32 loads and applies normalization factors established during training.

  • Line 33 prepares the model for evaluation (and training)

  • Line 34 loads the trained model coefficients from binary pkl files, which are not human-readable.

  • Line 35 summarizes the training errors for each output channel, individually by case or event, and over all cases.

• Lines 36-76 provide summary output, which includes:

  • The RMSE for each output channel and each case.

  • Identifying the case number that produced the highest RMSE for each output channel.

  • The number of cases in which the RMSE exceeded 0.05 per-unit, by output channel. This should be considered in judging whether the trained model is acceptable for use.

  • The total RMSE over all cases, for each output channel. These values are the same as output at the end of the training run. However, some cases will have higher RMSE values. Furthermore, some cases may exceed 0.05 RMSE, even when the total RMSE does not exceed 0.05 perunit.

# copyright 2021-2024 Battelle Memorial Institute
# summarizes RMSE and MAE for all training cases of an HW model for three-phase inverters
#  arg1: relative path to the model configuration file

import os
import sys
import json

import pecblocks.pv3_poly as pv3_model

bWantMAE = False

if __name__ == '__main__':
  if len(sys.argv) > 1:
    config_file = sys.argv[1]
    fp = open (config_file, 'r')
    cfg = json.load (fp)
    fp.close()
    data_path = cfg['data_path']
    model_folder = cfg['model_folder']
    model_root = cfg['model_root']
  else:
    print ('Usage: python pv3_metrics.py config.json')
    quit()

  print ('model_folder =', model_folder)
  print ('model_root =', model_root)
  print ('data_path =', data_path)

  model = pv3_model.pv3(training_config=config_file)
  model.loadTrainingData(data_path)
  model.loadAndApplyNormalization()
  model.initializeModelStructure()
  model.loadModelCoefficients()
  rmse, mae, case_rmse, case_mae = model.trainingErrors(True)
  out_size = len(model.COL_Y)
  colstr = ','.join('{:s}'.format(col) for col in model.COL_Y)
  h1str = ','.join('{:s}'.format('RMSE') for col in model.COL_Y)
  valstr = ','.join('{:.4f}'.format(rmse[j]) for j in range(out_size))
  if bWantMAE:
    h2str = ','.join('{:s}'.format('MAE') for col in model.COL_Y)
    maestr = ','.join('{:.4f}'.format(mae[j]) for j in range(out_size))
    print ('Idx,{:s},{:s}'.format(colstr, colstr))
    print ('#,{:s},{:s}'.format(h1str, h2str))
    print ('Total,{:s},{:s}'.format(valstr, maestr))
    for i in range(len(case_rmse)):
      valstr = ','.join('{:.4f}'.format(case_rmse[i][j]) for j in range(out_size))
      maestr = ','.join('{:.4f}'.format(case_mae[i][j]) for j in range(out_size))
      print ('{:d},{:s},{:s}'.format(i, valstr, maestr))
    print ('Total Error Summary')
    for j in range(out_size):
      print ('{:4s} MAE={:8.4f} RMSE={:8.4f}'.format (model.COL_Y[j], mae[j], rmse[j]))
  else:
    print ('Idx,{:s}'.format(colstr))
    print ('#,{:s}'.format(h1str))
    print ('Total,{:s}'.format(valstr))
    for i in range(len(case_rmse)):
      valstr = ','.join('{:.4f}'.format(case_rmse[i][j]) for j in range(out_size))
      print ('{:d},{:s}'.format(i, valstr))
    print ('Highest RMSE Cases')
    for j in range(out_size):
      mval = 0.0
      idx = 0
      nmax = 0
      for i in range(len(case_rmse)):
        val = case_rmse[i][j]
        if val > 0.05:
          nmax += 1
        if val > mval:
          idx = i
          mval = val
      print ('{:4s} Max RMSE={:8.4f} at Case {:d}; {:d} > 0.05'.format (model.COL_Y[j], 
                                                                        mval, idx, nmax))
    print ('Total Error Summary')
    for j in range(out_size):
      print ('{:4s} RMSE={:8.4f}'.format (model.COL_Y[j], rmse[j]))

Running a Model

Once the model has been exported, it can be run in a dynamic or EMT simulation with forward evaluation of the HWPV blocks. Two examples are described below. In addition, the models have been run in MATLAB/Simscape and ATP, with a Cigre/IEEE “real code” DLL implementation in progress.

Co-simulation using HELICS

This example runs in the z domain, with a fixed time step that must match the time step used to train the model. It runs as a co-simulation, with two Python federates communicating over the HELICS interface. One federate plays the role of the grid, while the other plays the role of a HWPV inverter model. In a more realistic use case, the grid federate might be a commercial dynamics or EMT simulator that already supports HELICS.

• See IEEE Access Paper for more information about HELICS.

• The example files and a description are at https://github.com/pnnl/pecblocks/tree/master/examples/helics

Standalone Simulation using Python

This example runs in the s domain using a Backward Euler integration of the H1 block. The time step is longer than the step used to train the model. Although not shown in this example, the Backward Euler method allows for variable time step during a simulation.

• The example files and a description are at https://github.com/pnnl/pecblocks/tree/master/examples/pi

• The file hwpv_pi.py plays the role of a grid federate, applying the test stimuli, calculating voltage from current at the PCC, and calling the HWPV model at each time step

• The file hwpv_evaluator.py implements the HWPV block model using numpy, but not pecblocks, dynonet, or torch. It can run on a Raspberry Pi or larger computer.

• The evaluator supports the z domain, the s domain with Forward Euler (less accurate and possibly unstable), and the s domain with Backward Euler (preferred).