Approximating Steam Properties: Machine Learning on IAPWS Standard

Aneet Dharmavaram Narendranath, PhD
Michigan Technological University,
Houghton, MI 49931 USA

Join the Conversation #WolframTechConf

Who am I?

Senior lecturer in mechanical engineering-engineering mechanics

ML opportunity in traditional thermodynamics pedagogy?

“Decision tree” type approach to triangulating or interpolating thermodynamic properties of steam.  
Boles, Michael A., and Yunus A. Cengel. Thermodynamics: an engineering approach. McGraw-Hill, 1989.

Abstract

Predict[] and NetTrain[] are used as thermodynamic property approximators for water and steam, by training on the IF-97 industrial standard of the International Association of Properties of Water and Steam (IAPWS IF-97).   The objective is to use the thermodynamic approximator example presented here, as an experiential learning tool to demystify machine learning and neural networks to mechanical engineering students with a traditionally non-computing coursework background.


Regions of steam per the IAPWS IF97 industrial standard.  
Wagner, W., et al. Properties of Water and Steam, IAPWS-IF97: The Industrial Standard IAPWS-IF97 for the Thermodynamic Properties and Supplementary Equations for Other Properties. Springer-Verlag: Berlin, 1998.

There exist publications on single EOS/non segmented thermal/structural physics problems: Sencan (2006, 2011), Sun (1997), Hezave (2012), Kalogirou (2000) and in empirical structural dynamics Haque (2002)

Roadmap of machine learning approximations in this demo

Graph[{"Predict" \[DirectedEdge] "full data", "Predict"\[DirectedEdge]"1 Bar", "HPNN"\[DirectedEdge]"full data", "HPNN" \[DirectedEdge] "1 Bar", "HPNN"\[DirectedEdge]"z1","HPNN"\[DirectedEdge] "z2","HPNN"\[DirectedEdge] "z3","HPNN"\[DirectedEdge] "z5", "HPNN"\[DirectedEdge] "z5","z1"\[DirectedEdge]"IAPWS","z2"\[DirectedEdge]"IAPWS","z3"\[DirectedEdge]"IAPWS","z5"\[DirectedEdge]"IAPWS", "1 Bar"\[DirectedEdge]"IAPWS","full data"\[DirectedEdge] "IAPWS"},VertexShapeFunction->"Diamond",VertexSize->Medium,VertexLabels->Automatic,ImageSize->Large,BaseStyle->{FontSize->18}]

Comparisons made

Text version of graphs

Import thermodynamic data from webpage and examine

Apply Predict/Gradient Boosted Trees (automatic ML) to approximate (T,P) -> h on training set randomly sampled from full dataset

How good/bad is the prediction?

Where in the specific enthalpy space is the most error?

Apply Predict/Gradient Boosted Trees to approximate (T,(P)) -> h on training set randomly sampled from 1 Bar Isobar

How good/bad is the prediction?

Where in the specific enthalpy space is the most error?

Heuristically create a NetTrain multilayer feedforward percepteron (neural net) to approximate (T,P) -> h on training set randomly sampled from 1 Bar Isobar

How good/bad is the prediction?

Where in the specific enthalpy space is the most error?

Discretize dataset into “thermodynamic zones” per the IAPWS standard AND then heuristically create a NetTrain neural net to approximate (T,P) -> h on training set randomly sampled in each zone

How good/bad is the prediction?

Where in the thermodynamic plane is the most error?

Outcomes have a yellow background

"Further steps" or "further investigation" has a blue background

Steam data: (T,P) -> (region,h)

(*url="http://pages.mtu.edu/~dnaneet/thermo/hptr_dense_1.txt";*)
(*url="http://pages.mtu.edu/~dnaneet/thermo/hptr.txt";*)
(*datasetDirty=Import[url, "Data"]; (*has NaN data*)*)
SetDirectory[NotebookDirectory[]];
datasetDirty = Import["hptr.txt","Data"];
(*names={'celcius','bar','region','kJPerKg'}*)
s=Dimensions[datasetDirty][[1]];
dataset=Select[datasetDirty, #[[3]]≠ 0&]; (*Select only those rows and columns with non zero region*)
s=Dimensions[dataset][[1]];

Size, shape and contents of dataset

s
TableForm@Take[dataset,5](*names={'celcius','bar','region','kJPerKg'}*)

Clustering of data: FindClusters[]

Data division

T=dataset[[2;;Length[dataset],1]]; (*temperature, C*)
P=dataset[[2;;Length[dataset],2]]; (*pressure, bar*)
r=dataset[[2;;Length[dataset],3]]/.{0->"r0",1->"r1", 2->"r2", 3->"r3", 4->"r4", 5->"r5"}; (*region, class/label*)
h=dataset[[2;;Length[dataset],4]]; (*specific enthalpy, kJ/kg*)

Clustering

subset3d=Thread[{P,T,h}];
c3d=FindClusters[subset3d,4,Method->"KMeans"];

ListPlot3D[{c3d[[1]],c3d[[2]],c3d[[3]],c3d[[4]]},ViewPoint->{-Pi,-1.95Pi,6},ImageSize->Medium,Mesh->None,PlotLegends->{"zone 1 (KMeans)","zone 2 (KMeans)","zone 3 (KMeans)","zone 4 (KMeans)"},ViewPoint->{-Pi,-1.95Pi,6},ImageSize->Large]

Clustering of data: IAPWS IF97

ZonePointCloud[zoneNum_]:= Module[{znum=zoneNum},
zoneSelection=Select[dataset, #[[3]]== znum&];
tz=Table[zoneSelection[[i]][[1]],{i,1,Length[zoneSelection],1}];
pz=Table[zoneSelection[[i]][[2]],{i,1,Length[zoneSelection],1}];
hz=Table[zoneSelection[[i]][[4]],{i,1,Length[zoneSelection],1}];
SeedRandom[1234];
Table[
{tz[[i]],pz[[i]],hz[[i]]},{i,1,Length[zoneSelection],1}]];

{ListPlot3D[{ZonePointCloud[1],ZonePointCloud[2],ZonePointCloud[3],ZonePointCloud[5]},PlotStyle->{PointSize->0.02},Mesh->None,PlotLegends->{"zone 1 (IAPWS)","zone 2 (IAPWS)","zone 3 (IAPWS)","zone 5 (IAPWS)"},ViewPoint->{-Pi,-1.95Pi,6},ImageSize->Medium],

ListPlot3D[{c3d[[1]],c3d[[2]],c3d[[3]],c3d[[4]]},PlotStyle->{PointSize->0.02},Mesh->None,PlotLegends->{"zone 1 (KMeans)","zone 2 (KMeans)","zone 3 (KMeans)","zone 5 (KMeans)"},ViewPoint->{-Pi,-1.95Pi,6},ImageSize->Medium]
}

The distribution of specific enthalpy data has some qualitative similarity between IAPWS and KMeans

Automatic Machine Learning on full data: this is baseline approximation

Gradient Boosted Trees are utilized to fit a sophisticated surface to the (T, P) -> h data

thermoState=Thread[{T,P}]; (*thermodynamic state*)
SeedRandom[123];
mldata=RandomSample[Thread[thermoState->h]]; 
s=Dimensions[mldata][[1]];
training=mldata[[1;;Ceiling[0.75*s]]];
validation=mldata[[Ceiling[0.75*s]+1;;s]];
(*pThermo0=Predict[training,Method->"DecisionTree"]*)
pThermo1=Predict[training,Method->"GradientBoostedTrees", RandomSeeding->123,TrainingProgressReporting->"ProgressIndicator"]

Break glass in case of fire

Prediction Quality

The nWorst=300 predicted examples (ground truth) are :

Off[General::noinfo];
nWorst=100;
PredictorMeasurements[pThermo1,validation];
wpe1=%["WorstPredictedExamples"->nWorst];
Take[wpe1,10]
(*Dimensions[wpe1]*)
(CopyTag[PredictorMeasurementsObject[]])

Percentage error Module

PercentageError1[nWorst_,WorstPredExamples_, PredObject_]:=Module[{nw=nWorst, wpe=WorstPredExamples, pobj=PredObject},
statesBadPred=Table[wpe[[i]][[1]],{i,1,nw,1}];
d1=Dimensions[statesBadPred][[1]];
checkWithThese=Table[Select[wpe, #[[1]]==statesBadPred[[i]]&],{i,1,d1,1}];
Flatten[
Table[{checkWithThese[[i]][[1]][[2]]-pThermo1[statesBadPred[[i]]]}*100/ checkWithThese[[i]][[1]][[2]],{i,1,d1,1}]
]
]

Percentage error between ground truth and automatic machine learning ranges between 1-100+%

Percentage error plot

ListLogPlot[Abs[PercentageError1[Dimensions[wpe1][[1]],wpe1,pThermo1]],GridLines->{None,{1,5}},GridLinesStyle->Directive[AbsoluteThickness[1.25],Blue],
Method->{"GridLinesInFront"->True},Frame->True,FrameStyle->Black,PlotStyle->{Red,PointSize->0.02},Joined->True,PlotRange->{All,{0,100}},FrameLabel->{"Prediction # x 10","%Error"},Filling->Axis]

Predict/GBDT approximator has errors close to band of phase change

satState={{81,0.5},{100,1},{150,4.76 },{200,15.55},{250,40},{300,85},{350,165}, {375, 218.14}};
Show[{
ListLogPlot[satState,PlotRange->{{0.1,500},{0.01,500}},Frame->True,FrameStyle->Black,PlotStyle->{Red,PointSize->0.035},FrameLabel->{"Temperature [C]","Pressure [bar]"},PlotLabel->"Saturation states in Red \!\(\*InterpretationBox[ButtonBox[
TooltipBox[GraphicsBox[{{GrayLevel[0], RectangleBox[{0, 0}]}, {GrayLevel[0], RectangleBox[{1, -1}]}, {RGBColor[1, 0, 0], RectangleBox[{0, -1}, {2, 1}]}},AspectRatio->1,DefaultBaseStyle->"ColorSwatchGraphics",Frame->True,FrameStyle->RGBColor[0.6666666666666666, 0., 0.],FrameTicks->None,ImageSize->Dynamic[{Automatic, 1.35 (CurrentValue["FontCapHeight"]/AbsoluteCurrentValue[Magnification])}],PlotRangePadding->None],StyleBox[
RowBox[{"RGBColor", "[", 
RowBox[{"1", ",", "0", ",", "0"}], "]"}], NumberMarks -> False]],Appearance->None,BaseStyle->{},BaselinePosition->Baseline,ButtonFunction:>With[{Typeset`box$ = EvaluationBox[]}, 
If[
Not[
AbsoluteCurrentValue["Deployed"]], SelectionMove[Typeset`box$, All, Expression]; FrontEnd`Private`$ColorSelectorInitialAlpha = 1; FrontEnd`Private`$ColorSelectorInitialColor = RGBColor[1, 0, 0]; FrontEnd`Private`$ColorSelectorUseMakeBoxes = True; MathLink`CallFrontEnd[
FrontEnd`AttachCell[Typeset`box$, 
FrontEndResource["RGBColorValueSelector"], {0, {Left, Bottom}}, {Left, Top}, "ClosingActions" -> {"SelectionDeparture", "ParentChanged", "EvaluatorQuit"}]]]],DefaultBaseStyle->{},Evaluator->Automatic,Method->"Preemptive"],RGBColor[1, 0, 0],Editable->False,Selectable->False]\) \n Worst predicted examples in Blue \!\(\*InterpretationBox[ButtonBox[
TooltipBox[GraphicsBox[{{GrayLevel[0], RectangleBox[{0, 0}]}, {GrayLevel[0], RectangleBox[{1, -1}]}, {RGBColor[0, 0, 1], RectangleBox[{0, -1}, {2, 1}]}},AspectRatio->1,DefaultBaseStyle->"ColorSwatchGraphics",Frame->True,FrameStyle->RGBColor[0., 0., 0.6666666666666666],FrameTicks->None,ImageSize->Dynamic[{Automatic, 1.35 (CurrentValue["FontCapHeight"]/AbsoluteCurrentValue[Magnification])}],PlotRangePadding->None],StyleBox[
RowBox[{"RGBColor", "[", 
RowBox[{"0", ",", "0", ",", "1"}], "]"}], NumberMarks -> False]],Appearance->None,BaseStyle->{},BaselinePosition->Baseline,ButtonFunction:>With[{Typeset`box$ = EvaluationBox[]}, 
If[
Not[
AbsoluteCurrentValue["Deployed"]], SelectionMove[Typeset`box$, All, Expression]; FrontEnd`Private`$ColorSelectorInitialAlpha = 1; FrontEnd`Private`$ColorSelectorInitialColor = RGBColor[0, 0, 1]; FrontEnd`Private`$ColorSelectorUseMakeBoxes = True; MathLink`CallFrontEnd[
FrontEnd`AttachCell[Typeset`box$, 
FrontEndResource["RGBColorValueSelector"], {0, {Left, Bottom}}, {Left, Top}, "ClosingActions" -> {"SelectionDeparture", "ParentChanged", "EvaluatorQuit"}]]]],DefaultBaseStyle->{},Evaluator->Automatic,Method->"Preemptive"],RGBColor[0, 0, 1],Editable->False,Selectable->False]\)",GridLines->{{373},{218}}],
ListLogPlot[{Take[statesBadPred,d1]},PlotStyle->{Blue,PointSize->0.3},FrameLabel->{"Temperature [C]","Pressure [bar]"},ImageSize->Large]
}]

The worst approximated values band around the saturation state.  Phase change from saturated liquid to saturated liquid vapor, may not be predicted/approximated very well.

Calculation of Power produced by a Steam Turbine

The state of steam at the entrance and the exit to a reversible, adiabatic turbine is (T1, P1) = (400 °C, 5.1 MPa) and (T2, P2) = (100 °C, 100 kPa).  If the changes in specific kinetic energy and specific potential energy are negligible in comparison to change in specific enthalpy, find the useful power generated by this turbine when operating at steady state.

Approximating specific enthalpy using “Predict” object and applying the steady state, first law of thermodynamics for an open system that is the steam turbine:

{Tin, Pin} = {400, 51};
{Tout, Pout} = {100, 1};
{
(*hin=pThermo0[{Tin, Pin}],
hout=pThermo0[{Tout, Pout}],*)
hin=pThermo1[{Tin, Pin}],
hout=pThermo1[{Tout, Pout}]}
wout0 =Quantity[hin - hout,"kJ/kg"]

Actual first law of thermodynamics for the open system that is the steam turbine:

Quantity[First[Select[mldata, #[[1]]=={Tin,Pin}&]][[2]] -First[Select[mldata, #[[1]]=={Tout,Pout}&]][[2]],"kJ/kg"] //Flatten

The predict/GBDT object overestimates the power produced by the steam turbine.   Key takeaway for students is -- validation of numerical/predictive results with empirical or classical results!

Other distributions of training and validation data should be tested to check for effect on approximation

Partition data into isobaric sample(s)

Choose data on the isobar, P = 1 Bar (atmospheric pressure).

This is an intuitive isobar with physical significance that most of us understand (we live in it, breathe it, ..., boil water at this pressure)

isobar0=Select[datasetDirty, #[[2]]≠0 && #[[2]] == 1&];
isobar1=Select[isobar0,#[[4]]>0&];
T1=isobar1[[All,1]]; P1=isobar1[[All,2]]; h1=isobar1[[All,4]];

s1=Dimensions[h1][[1]];
thermoState1=Thread[{T1,P1}]; (*thermodynamic state*)
SeedRandom[123];
mldata1=RandomSample[Thread[thermoState1->h1]];

training1=mldata1[[1;;Ceiling[0.75*s1]]]; 
validation1=mldata1[[Ceiling[0.75s1]+1;;s1]];

pThermo1Isobar1=Predict[training1,Method->"GradientBoostedTrees", RandomSeeding->123, TrainingProgressReporting->"ProgressIndicator"]

Prediction vs ground truth at phase transition point @1 bar, 100 C has an error banded mostly between 1-10% and is the same as the auto ML on full data

Off[General::noinfo];
nWorst=100;
PredictorMeasurements[pThermo1Isobar1,validation];
wpe2=%["WorstPredictedExamples"->nWorst];

{ListLogPlot[Abs[PercentageError1[Dimensions[wpe2][[1]],wpe2,pThermo1Isobar1]],GridLines->{None,{1,5}},GridLinesStyle->Directive[AbsoluteThickness[1.25],Blue],
Method->{"GridLinesInFront"->True},Frame->True,FrameStyle->Black,PlotStyle->{Red,PointSize->0.02},Joined->True,PlotRange->{All,{0,100}},FrameLabel->{"Prediction # x 10","%Error"},Filling->Axis,PlotLabel->"Auto ML on 1 Bar",ImageSize->Medium],

ListLogPlot[Abs[PercentageError1[Dimensions[wpe2][[1]],wpe2,pThermo1]],GridLines->{None,{1,5}},GridLinesStyle->Directive[AbsoluteThickness[1.25],Blue],
Method->{"GridLinesInFront"->True},Frame->True,FrameStyle->Black,PlotStyle->{Red,PointSize->0.02},Joined->True,PlotRange->{All,{0,100}},FrameLabel->{"Prediction # x 10","%Error"},Filling->Axis,PlotLabel->"Auto ML on full data",ImageSize->Medium]}

Percepteron trained on 1 bar isobar

Heuristically constructed percepteron with linear layers and Tanh nonlinear/activation layer.

net=NetChain[{
64,Tanh,32,8,1}]
s1=Dimensions[mldata1][[1]]; (*mldata1 is (T, P) -> h at 1 Bar*)
mldata2=Table[mldata1[[i,1]]->{mldata1[[i,2]]},{i,1,s1,1}]; (*mldata1 is already randomsample-d*)
training = Take[mldata2,Ceiling[0.75s1]];

trainedP1Data=NetTrain[net,training,RandomSeeding->123,TrainingProgressReporting->"ProgressIndicator"]

Break glass in case of  fire

Which prediction method is best? 
(OOTB GBDT/full data, OOTB GBDT/isobar, NN/isobar)

Percentage error module for NetTrain objects

PercentageError2[groundTruth_, predictionObject_] := 
  Module[{gt = groundTruth, po = predictionObject}, 
   d1 = Length[gt]; 
    Abs[Flatten[Table[(gt[[i,2]] - po[gt[[i]][[1]]])*
        (100/gt[[i,2]]), {i, 1, d1, 1}]]]]

Percentage error plots

{ListLogPlot[PercentageError2[mldata1,pThermo1],Frame->True,PlotRange->{All,{5,500}},FrameStyle->Black,FrameLabel->{"Temperature instance","%Error"},PlotStyle->{{Dashed,Thickness[0.01], Black}},PlotLegends->{"GBDT trained on full data \napplied to isobar"},Joined->True,ImageSize->Medium],

ListLogPlot[PercentageError2[mldata1,trainedP1Data],Frame->True,PlotRange->{All,{5,500}},FrameStyle->Black,FrameLabel->{"Temperature instance","%Error"},PlotStyle->{{Dashed,Thickness[0.01], Blue}},PlotLegends->{"NN trained on isobar \napplied to isobar"},Joined->True,ImageSize->Medium]}

The Predict[]/GBDT on full dataset, Predict[]/GBDT on isobar and HPNN on isobar, all have deficiencies.

Other distributions of training and validation data should be tested to check for effect on approximation.

More data discretization?

Percepteron trained on complete dataset

s2=Dimensions[mldata][[1]];
mldata3=Table[mldata[[i,1]]->{mldata[[i,2]]},{i,1,s2,1}];
training2=Take[mldata3,Ceiling[0.75s2]];

NNAllData=NetTrain[net,training2,RandomSeeding->123,TrainingProgressReporting->"ProgressIndicator"]

Screenshot with loss in J/kg (this takes a little long to demo live)

{ListLogPlot[PercentageError2[mldata1,NNAllData],Frame->True,PlotRange->{All,{5,500}},FrameStyle->Black,FrameLabel->{"Temperature instance","%Error"},PlotStyle->{{Dashed,Thickness[0.01], Blue}},PlotLegends->{"Percepteron trained on all data \n(includes all isobars) \napplied to 1 bar isobar"},Joined->True,ImageSize->Medium],

ListLogPlot[PercentageError2[mldata1,trainedP1Data],Frame->True,PlotRange->{All,{5,500}},FrameStyle->Black,FrameLabel->{"Temperature instance","%Error"},PlotStyle->{{Dashed,Thickness[0.01], Red}},PlotLegends->{"Percepteron trained on isobar \napplied to 1 bar isobar"},Joined->True,ImageSize->Medium]}

The percepteron trained on the full dataset and applied to 1 bar isobar is not as good as the percepteron trained on the 1 bar isobar.

Can we cluster data per the IAPWS standard and then retrain percepterons on each cluster of (P, T) -> h?

Partition of data into IAPWS zone for better approximation (error is 1-2 orders of mag lower)

Module for sequencing data into zones

ZoneData[zoneNum_]:= Module[{znum=zoneNum},
zoneSelection=Select[dataset, #[[3]]== znum&];
tz=Table[zoneSelection[[i]][[1]],{i,1,Length[zoneSelection],1}];
pz=Table[zoneSelection[[i]][[2]],{i,1,Length[zoneSelection],1}];
hz=Table[zoneSelection[[i]][[4]],{i,1,Length[zoneSelection],1}];
SeedRandom[1234];
Thread[
RandomSample[Table[
{tz[[i]],pz[[i]]}->{hz[[i]]},{i,1,Length[zoneSelection],1}]]
]
]

Zone specific training

zone1Data=ZoneData[1];
zone1Training=Take[zone1Data,Ceiling[0.7Length[zone1Data]]];
zone1Validation=Take[zone1Data,-Ceiling[0.3Length[zone1Data]]];
Length[zone1Data]

nnz1=NetTrain[net,zone1Training,RandomSeeding->123, TrainingProgressReporting->"ProgressIndicator"]

n=300;
{ListLogPlot[PercentageError2[RandomSample[zone1Validation,n],nnz1],Frame->True,PlotRange->{All,{5,500}},FrameStyle->Black,FrameLabel->{"Temperature instance","%Error"},PlotStyle->{{Dashed,Thickness[0.01], Blue}},PlotLegends->{"Percepteron trained on zone 1 \napplied to zone 1"},Joined->True,ImageSize->Medium],
ListLogPlot[PercentageError2[RandomSample[zone1Validation,n],pThermo1],Frame->True,PlotRange->{All,{5,500}},FrameStyle->Black,FrameLabel->{"Temperature instance","%Error"},PlotStyle->{{Dashed,Thickness[0.01], Red}},PlotLegends->{"GBDT trained on all data \napplied to zone 1"},Joined->True,ImageSize->Medium]}

The percepteron trained on zone 1 is superior by roughly an order of magnitude than the Predict[]/GBDT trained on all data. The error was consistently O(100) while that of the predict object was O(1000).

Partition of data based on region (zone) leads to better approximation

zone2Data=ZoneData[2];
zone2Training=Take[zone2Data,Ceiling[0.7Length[zone2Data]]];
zone2Validation=Take[zone2Data,-Ceiling[0.3Length[zone2Data]]];
Length[zone2Data]

nnz2=NetTrain[net,zone2Training,RandomSeeding->123,TrainingProgressReporting->"ProgressIndicator"]

n=300;
{ListLogPlot[PercentageError2[RandomSample[zone2Validation,n],nnz2],Frame->True,PlotRange->{All,{5,500}},FrameStyle->Black,FrameLabel->{"Temperature instance","%Error"},PlotStyle->{{Dashed,Thickness[0.01], Blue}},PlotLegends->{"Percepteron trained on zone 2 \napplied to zone 2"},Joined->True,ImageSize->Medium],
ListLogPlot[PercentageError2[RandomSample[zone2Validation,n],pThermo1],Frame->True,PlotRange->{All,{5,500}},FrameStyle->Black,FrameLabel->{"Temperature instance","%Error"},PlotStyle->{{Dashed,Thickness[0.01],Red}},PlotLegends->{"GBDT trained on all data \napplied to zone 2"},Joined->True,ImageSize->Medium]}

Length@Select[PercentageError2[RandomSample[zone2Validation,n],pThermo1],#>5&]
Length@Select[PercentageError2[RandomSample[zone2Validation,n],nnz2],#>5&]

The percepteron trained on zone 2 is superior to the Predict[]/GBDT trained on all data.  None of the approximation is greater than 5% in error.

Zone 3 is smaller datasample with less effective NN approximation

zone3Data=ZoneData[3];
z=zone3Data;
zone3Training=Take[z,Ceiling[0.7Length[z]]];
zone3Validation=Take[z,-Ceiling[0.3Length[z]]];
Length[zone3Data]

nnz3=NetTrain[net,zone3Training,RandomSeeding->123]

n=300;
{ListLogPlot[PercentageError2[RandomSample[zone3Validation,n],nnz3],Frame->True,PlotRange->{All,{5,500}},FrameStyle->Black,FrameLabel->{"Temperature instance","%Error"},PlotStyle->{{Dashed,Thickness[0.01], Blue}},PlotLegends->{"Percepteron trained on zone 3 \napplied to zone 3"},Joined->True,ImageSize->Medium],
ListLogPlot[PercentageError2[RandomSample[zone3Validation,n],pThermo1],Frame->True,PlotRange->{All,{5,500}},FrameStyle->Black,FrameLabel->{"Temperature instance","%Error"},PlotStyle->{{Dashed, Thickness[0.01],Red}},PlotLegends->{"GBDT trained on all data \napplied to zone3"},Joined->True,ImageSize->Medium]}

The percepteron trained on zone 3 is inferior to the Predict[]/GBDT model.  This could partly be because of a less dense data sample in zone 3 and/or the complexity (fold) in the data sample.

Zone 5 is smallest data sample but with more effective NN than for zone 3

zone5Data=ZoneData[5];
z=zone5Data;
zone5Training=Take[z,Ceiling[0.7Length[z]]];
zone5Validation=Take[z,-Ceiling[0.3Length[z]]];
Length[zone5Data]

nnz5=NetTrain[net,zone5Data,RandomSeeding->123]

n=300;
{ListLogPlot[PercentageError2[RandomSample[zone5Validation,n],nnz5],Frame->True,PlotRange->{All,{5,500}},FrameStyle->Black,FrameLabel->{"Temperature instance","%Error"},PlotStyle->{{Dashed, Thickness[0.01],Blue}},PlotLegends->{"percepteron trained on zone 5 \napplied to zone 5"},Joined->True,ImageSize->Medium],
ListLogPlot[PercentageError2[RandomSample[zone5Validation,n],pThermo1],Frame->True,PlotRange->{All,{5,500}},FrameStyle->Black,FrameLabel->{"Temperature instance","%Error"},PlotStyle->{{Dashed,Thickness[0.01], Red}},PlotLegends->{"GBDT trained on all data \napplied to zone 5"},Joined->True,ImageSize->Medium]}

Percepteron trained on zone 5 shows similar accuracy as the GBDT predict object.

For nearly the same size of data sample, zone 5 approximations are better than that of zone 3

Zone 3 has a “fold” in it’s surface

ListPlot3D[{ZonePointCloud[1],ZonePointCloud[2],ZonePointCloud[3],ZonePointCloud[5]},PlotStyle->{PointSize->0.02},Mesh->None,PlotLegends->{"zone 1 (IAPWS)","zone 2 (IAPWS)","zone 3 (IAPWS)","zone 5 (IAPWS)"},ViewPoint->{-Pi,-1.95Pi,6},ImageSize->Medium]

zone3Data=ZoneData[3];
z=zone3Data;
zone3Training=Take[z,Ceiling[0.7Length[z]]];
zone3Validation=Take[z,-Ceiling[0.3Length[z]]];

Figure source : Tuseeta Banerjee (2019), Applications of neural networks, Wolfram technology conference 2019

net2=NetChain[{
64,LogisticSigmoid,32,LogisticSigmoid,8,1},
"Input"->2];

nnz3Mod1=NetTrain[net2,zone3Data,RandomSeeding->123]

n=300;
SeedRandom[123];
{ListLogPlot[PercentageError2[RandomSample[zone3Validation,n],nnz3],Frame->True,PlotRange->{All,{5,20}},FrameStyle->Black,FrameLabel->{"Temperature instance","%Error"},PlotStyle->{{Dashed,Thickness[0.01], Blue}},PlotLegends->{"percepteron trained on zone 3 \napplied to zone 3"},Joined->True,ImageSize->Medium],
ListLogPlot[PercentageError2[RandomSample[zone3Validation,n],nnz3Mod1],Frame->True,PlotRange->{All,{5,20}},FrameStyle->Black,FrameLabel->{"Temperature instance","%Error"},PlotStyle->{{Dashed,Thickness[0.01], Red}},PlotLegends->{"percepteron trained on zone 3 (modified) \napplied to zone 3"},Joined->True,ImageSize->Medium]}

Zone 3 may just need greater data density for an adequate (error < 5%) neural network

For the data density available, switching the Tanh layer with a LogisticSigmoid layer helped improve the approximation

C̶o̶n̶c̶l̶u̶s̶i̶o̶n̶s̶ Summary

Water/Steam enthalpy data available from the IF-97 industrial standard is used to train ML approximators.

The Predict/GBDT and a heuristic multilayer percepteron neural network model are both trained on the data-set.  The Predict/GBDT has a better prediction than the percepteron model when training is performed on the full data-set.

The full dataset is partitioned into 4 zones per the IF-97 standard.  The percepteron is trained on these zones.  The predictions made by the zone specific percepterons (but for zone 3) are “equal” or “better” than the Predict/GBDT.

Overall: it is advantageous to partition enthalpy data into IAPWS zones and use (heuristically chosen) percpeterons.  This improved speed and accuracy

Next steps, from a teaching point of view:

ML based specific enthalpy look-up: Create a tool that locates a given (T,P) state in appropriate zone and use zone specific percepteron to evaluate simple turbines, boilers, condensers and heat exchangers.  The Mathematica ML ecosystem could be effectively used for demystification and inclusion of machine learning into undergraduate education.

Divide data into clusters using FindClusters and the “simple” KMeans algorithm.  Train percepterons on these clusters and compare with approximations made by the percepterons trained on IAPWS zones.

Neural network loss plots, stop criterion

Acknowledgements and List of References

Acknowledgements

Department chair of   @MichiganTech, Dr. William Predebon

Director of undergraduate studies MEEM, Dr.Jeffrey S. Allen

Members of the “undergraduate curriculum for Industry 4.0” workstream at the mechanical engineering-engineering mechanics department

Family and cats (who are also family)

References

Sencan A. Artificial intelligent methods for thermodynamic evaluation of ammonia–water refrigeration systems. Energy Conversion and Management 2006; 47(18-19): 3319–3332

Sencan A, Köse İİ and Selbaş R. Prediction of thermophysical properties of mixed refrigerants using artificial neural network. Energy conversion and management 2011; 52(2): 958–974

Sun DW. Thermodynamic design data and optimum design maps for absorption refrigeration systems. Applied thermal engineering 1997; 17(3): 211–221

Torrecilla JS, Rodrı́guez F, Bravo JL et al. Optimising an artificial neural network for predicting the melting point of ionic liquids. Physical Chemistry Chemical Physics 2008; 10(38): 5826– 5831

Kalogirou SA. Applications of artificial neural-networks for energy systems. Applied energy 2000; 67(1-2): 17–35

Haque ME An Artificial Neural Network Model for Preliminary Design of Reinforced Concrete Beam-Column. Proceedings of the American Society for Engineering Education, 2002 ASEE Annual Conference, 2002; 2306: 16–19

Boles, Michael A., and Yunus A. Cengel. Thermodynamics: an engineering approach. McGraw-Hill, 1989.

Wagner, W., and A. Kruse. Properties of Water and Steam, IAPWS-IF97: The Industrial Standard IAPWS-IF97 for the Thermodynamic Properties and Supplementary Equations for Other Properties. Springer-Verlag: Berlin, 1998

Appendix: Supplementary data

{Graph[{"Predict" \[DirectedEdge] "full data", "Predict" \[DirectedEdge]"1 Bar"},VertexShapeFunction->"Diamond",VertexSize->Medium,VertexLabels->Automatic,ImageSize->Small,BaseStyle->{FontSize->18}],
Graph[{"Predict" \[DirectedEdge] "1 Bar", "HPNN" \[DirectedEdge]"1 Bar"},VertexShapeFunction->"Diamond",VertexSize->Medium,VertexLabels->Automatic,ImageSize->Small,BaseStyle->{FontSize->18}],
Graph[{ "HPNN" \[DirectedEdge]"zone 1", "Predict" \[DirectedEdge]"zone 1"},VertexShapeFunction->"Diamond",VertexSize->Medium,VertexLabels->Automatic,ImageSize->Small,BaseStyle->{FontSize->18}],
Graph[{ "HPNN" \[DirectedEdge]"zone 2", "Predict" \[DirectedEdge]"zone 2"},VertexShapeFunction->"Diamond",VertexSize->Medium,VertexLabels->Automatic,ImageSize->Small,BaseStyle->{FontSize->18}],
Graph[{ "HPNN" \[DirectedEdge]"zone 3", "Predict" \[DirectedEdge]"zone 3"},VertexShapeFunction->"Diamond",VertexSize->Medium,VertexLabels->Automatic,ImageSize->Small,BaseStyle->{FontSize->18}],
Graph[{ "HPNN" \[DirectedEdge]"zone 5", "Predict" \[DirectedEdge]"zone 5"},VertexShapeFunction->"Diamond",VertexSize->Medium,VertexLabels->Automatic,ImageSize->Small,BaseStyle->{FontSize->18}],
Graph[{ "HPNN" \[DirectedEdge]"zone 3", "HPNN"\[DirectedEdge]"zone 3 \n mod"},VertexShapeFunction->"Diamond",VertexSize->Medium,VertexLabels->Automatic,ImageSize->Small,BaseStyle->{FontSize->18}]}

supHeat={{150,0.1}, {200,0.1},{300,0.1}, {400,0.1},
{250,2},{300,2},{350,2},{400,2},
{250,10},{300,10}, {350,10},{400,10},
{325,100},{375,100},{400,100},{450,100},
{250,400},{300,400},{350,400},{400,400}
};
subcooled={
{50,50},{100,50},
{50,100}, {100,100}
};

errRelief=Table[Flatten[Join[{wpe1[[i]][[1]],Abs[err[[i]]]}]],{i,1,30,1}];
ListPlot3D[errRelief,InterpolationOrder->0,Mesh->None,Filling->Bottom,ColorFunction->"TemperatureMap",PlotRange->{{0,2000},{0,1000},{0,100}},AxesLabel->{"Temperature [C]","Pressure [bar]","%Error"},PlotLegends->Automatic,PlotLabel->"Error Plot (>10%)"]//Inactivate;

Show[
ListPointPlot3D[{{99,1,100},{205,17,100},{250,40,100},{300,85,100},{325,120,100},{350,165,100}},PlotStyle->{Red,PointSize->0.03},PlotRange->{{0,2000},{0,1000},{0,100}}],ListPlot3D[{errRelief},InterpolationOrder->0,Mesh->None,Filling->Bottom,ColorFunction->"TemperatureMap",PlotRange->{{0,2000},{0,1000},{0,100}},AxesLabel->{"Temperature [C]","Pressure [bar]","%Error"},PlotLegends->Automatic,PlotLabel->"Error Plot (>10%)"]];

(*Switch off some messages*)
Off[NDSolve::mxsst]; 
Off[NDSolve::ibcinc]; 
Off[NDSolve::eerr]; 
Off[NDSolve::ndsz]; 
pbcs={ Derivative[0, 1][u][t, xMin] == Derivative[0, 1][u][t, xMax], 
      Derivative[0, 2][u][t, xMin] == Derivative[0, 2][u][t, xMax], 
      Derivative[0, 3][u][t, xMin] == Derivative[0, 3][u][t, xMax], 
      u[t, xMin] == u[t, xMax]};
met={ MaxSteps -> 100000, Method -> {"MethodOfLines", 
       "Method" -> "LSODA", "TemporalVariable" -> t, "SpatialDiscretization" -> 
        {"TensorProductGrid", "MinPoints" -> 800, "MaxPoints" -> 1200, 
         "DifferenceOrder" -> 5}}};