Calibrate heat pump model from performance data

This example demonstrates the use of the calibrate module to obtain parameters of the heat pump model. Data is loaded from a performance file generated in the performance data generation example. Once heat pump parameters are identified, the results are verified by running the model in Dymola.

The following script can be found in IBPSA/Resources/src/fluid/heatpumps/calibration/Examples/example_calibration.py:

  1# -*- coding: utf-8 -*-
  2""" Example calibration of the heat pump model.
  3
  4This script demonstrates the use of the calibration module to obtain parameters
  5of the heat pump model. Data is loaded from a performance file generated in
  6dummy_performance_data.py. Once heat pump parameters are identified, the
  7results are verified by running the model in Dymola. A modelica record of the
  8heat pump parameters is generated.
  9
 10"""
 11from __future__ import division, print_function, absolute_import
 12
 13import numpy as np
 14import os
 15import sys
 16import datetime
 17
 18
 19def main():
 20    # Add parent directory to system path
 21    parent_dir = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
 22    sys.path.insert(1, parent_dir)
 23    # Import Heat pump and calibration module
 24    import PythonModel as hp
 25    # Change working directory to current directory
 26    os.chdir(os.path.dirname(__file__))
 27    # Author info for record generation
 28    author = 'Mr Modeler'
 29    # Make and model of the heat pump
 30    manufacturer = 'SomeManufacturer'
 31    model = 'ABC060'
 32    # Set to True if calibrating for cooling mode
 33    CoolingMode = False
 34    # File name and table name for manufacturer data in modelica
 35    tableFileName = 'manufacturerData.txt'
 36    tableName = 'ManufacturerData'
 37    # File name for performance data
 38    performanceData = 'somePerformanceData.txt'
 39
 40    # Load manufacturer data
 41    data = hp.calibrate.ManufacturerData(manufacturer, model, CoolingMode)
 42    with open(performanceData, 'r') as f:
 43        for line in f:
 44            dataPoint = line[0:-2].split('\t')
 45            EWT_Source = float(dataPoint[0])    # Entering water temperature evaporator (K)
 46            EWT_Load = float(dataPoint[1])      # Entering water temperature condenser (K)
 47            flowSource = float(dataPoint[2])    # Entering water flow rate evaporator (kg/s)
 48            flowLoad = float(dataPoint[3])      # Entering water flow rate condenser (kg/s)
 49            Capacity = float(dataPoint[4])      # Heat transfer rate on the load side (kW).
 50            HR = float(dataPoint[5])            # Heat transfer rate on the source side (kW).
 51            Power = float(dataPoint[6])         # Electrical power input to heat pump (kW)
 52            # Add data point to Data object
 53            data.add_data_point(EWT_Source, EWT_Load, flowSource,
 54                                flowLoad, Capacity, HR, Power)
 55
 56    # Data points used in calibration
 57    calData = data.calibration_data_16_points()
 58
 59    # Initialize the heat pump model
 60    P_nominal = 17.5e3
 61    COP_nominal = 4.0
 62    Q_nominal = P_nominal*COP_nominal
 63
 64    # -------------------------------------------------------------------------
 65    # Initialize all models using a value of 0. for all parameters. Parameters
 66    # will be replaced by guess values at the start of the calibration process.
 67    # -------------------------------------------------------------------------
 68    # Compressor model (Scroll)
 69    com = hp.compressors.ScrollCompressor([0., 0., 0., 0., 0., 0.])
 70    # Condenser model
 71    con = hp.heatexchangers.EvaporatorCondenser([0.])
 72    # Evaporator model
 73    eva = hp.heatexchangers.EvaporatorCondenser([0.])
 74    # Refrigerant model
 75    ref = hp.refrigerants.R410A()
 76    # Fluid model on condenser side
 77    fluCon = hp.fluids.ConstantPropertyWater()
 78    # Fluid model on evaporator side
 79    fluEva = hp.fluids.ConstantPropertyWater()
 80    # Heat pump model
 81    heaPum = hp.heatpumps.SingleStageHeatPump(com, con, eva, ref, fluCon,
 82                                              fluEva, Q_nominal, P_nominal,
 83                                              CoolingMode)
 84
 85    # Lauch the calibration of the heat pump model.
 86    optPar, optRes, gueRes = hp.calibrate.calibrate_model(heaPum, calData,
 87                                                          data, plot=True)
 88
 89    # Write the results into a record for use in Modelica
 90    write_record_scroll(author, manufacturer, model, CoolingMode,
 91                        'R410A', Q_nominal, COP_nominal,
 92                        optPar)
 93
 94    # -------------------------------------------------------------------------
 95    # Calculate heat pump performance for full dataset in Dymola using the
 96    # calibrated parameters.
 97    # -------------------------------------------------------------------------
 98    dymRes = hp.calibrate.simulate_in_dymola(heaPum, data, tableName,
 99                                             tableFileName)
100    SSE = hp.calibrate.compare_data_sets(dymRes, data, plot=True,
101                                         fname=data.name + '_dymola')
102    print('----------------------------------------------------------------\n')
103    print('Sum of square errors (dymola) : ' + str(SSE) + ' \n')
104    print('----------------------------------------------------------------\n')
105
106    # Compare the results of the Python code with the results from Dymola.
107    SSE = hp.calibrate.compare_data_sets(dymRes, optRes, plot=True,
108                                         fname='modelVerification')
109    return optRes, dymRes
110
111
112def write_record_scroll(author, manufacturer, model, CoolingMode,
113                        refrigerant, Q_nominal, COP_nominal,
114                        optPar):
115    # Evaluate current date
116    date = str(datetime.date.today())
117
118    # Extract heat pump parameters from optimized results
119    volRat = optPar[0]     # Built-in volume ratio (-)
120    v_flow = optPar[1]     # Volume flow rate at suction (m3/s)
121    leaCoe = optPar[2]     # Leakage coefficient (kg/s)
122    etaEle = optPar[3]     # Electro-mechanical efficiency (-)
123    PLos = optPar[4]       # Constant part of the power losses (W)
124    dTSup = optPar[5]      # Amplitude of superheating (K)
125    UACon = optPar[6]      # Thermal conductance of the condenser (W/K)
126    UAEva = optPar[7]      # Thermal conductance of the evaporator (W/K)
127
128    # Operation mode
129    if CoolingMode:
130        mode = 'Cooling'
131    else:
132        mode = 'Heating'
133
134    # Build string for nominal capacity. If less than 10kW, the nominal
135    # capacity is printed with the first decimal (ex. 9.1kW -> 9_1kW)
136    # otherwise the capacity is rounded (ex. 10.6kW - > 11kW).
137    if Q_nominal < 10.0e3:
138        Q_str = (str(int(Q_nominal/1.0e3)) +
139                 '_' +
140                 str(int(round(10.*(Q_nominal-np.floor(Q_nominal))/1.0e3))))
141    else:
142        Q_str = str(int(Q_nominal/1.0e3))
143
144    # Build the full name of the record, including the manufacturer, model,
145    # nominal capacity, nominal COP (keeping 2 decimals) and refrigerant type.
146    name = '_'.join([manufacturer,
147                     model,
148                     Q_str + 'kW',
149                     str(int(COP_nominal)),
150                     str(int(round(100.*(COP_nominal
151                                         -np.floor(COP_nominal))))) + 'COP',
152                     refrigerant])
153    path = name + '.mo'
154
155    # Print the record in Modelica format.
156    with open(path, 'w') as f:
157        f.write('within IBPSA.Fluid.HeatPumps.Data.ScrollWaterToWater.'
158                + mode + ';\n')
159        f.write('record ' + name + ' =\n')
160        f.write('  IBPSA.Fluid.HeatPumps.Data.ScrollWaterToWater.Generic ('
161                + '\n')
162        f.write('    volRat = ' + str(volRat) + ',\n')
163        f.write('    V_flow_nominal = ' + str(v_flow) + ',\n')
164        f.write('    leaCoe = ' + str(leaCoe) + ',\n')
165        f.write('    etaEle = ' + str(etaEle) + ',\n')
166        f.write('    PLos = ' + str(PLos) + ',\n')
167        f.write('    dTSup = ' + str(dTSup) + ',\n')
168        f.write('    UACon = ' + str(UACon) + ',\n')
169        f.write('    UAEva = ' + str(UAEva) + ')\n\n')
170
171        f.write('  annotation (\n')
172        f.write('    defaultComponentPrefixes = "parameter",\n')
173        f.write('    defaultComponentName="datHeaPum",\n')
174        f.write('    preferredView="info",\n')
175        f.write('  Documentation(info="<html>\n')
176        f.write('<p>\n')
177        f.write('Generated by ' + author + ' on ' + date + '.\n')
178        f.write('</p>\n')
179        f.write('</html>"));')
180    return
181
182# Main function
183if __name__ == "__main__":
184    main()

The scripts generates the following figures:

_images/calibration_guess_parameters.png

Fig. 1 Comparison between the performance data and the model predicted (from Python) heating capacities and power input using the initial guess parameters.

_images/calibration_final_parameters.png

Fig. 2 Comparison between the performance data and the model predicted (from Dymola) heating capacities and power input using the calibrated parameters.

The following record is generated for use in Modelica:

 1within IBPSA.Fluid.HeatPumps.Data.ScrollWaterToWater.Heating;
 2record SomeManufacturer_ABC060_70kW_4_0COP_R410A =
 3  IBPSA.Fluid.HeatPumps.Data.ScrollWaterToWater.Generic (
 4    volRat = 2.36185869323,
 5    V_flow_nominal = 0.00287086738316,
 6    leaCoe = 0.00408074950474,
 7    etaEle = 0.922319423243,
 8    PLos = 398.665383784,
 9    dTSup = 6.49606895698,
10    UACon = 7014.54967992,
11    UAEva = 49135.9514647)
12
13  annotation (
14    defaultComponentPrefixes = "parameter",
15    defaultComponentName="datHeaPum",
16    preferredView="info",
17  Documentation(info="<html>
18<p>
19Generated by Mr Modeler on 2017-01-26.
20</p>
21</html>"));