Neural Networks in System Models

In this presentation, you will see how data provides new avenues and connections to physical modeling with the new support for neural networks in system modeling. Applications and examples will be explored.

System Modeler and physical modeling

Multiple domains:

SystemModel["Modelica","Children"]

Modeling based on connections:

SystemModeler["RC","Modeling"]

Differential-algebraic system of equations:

F[x[t],x'[t],t]==0

Supports real-time connection to devices and co-simulation

See a robot:

SystemModeler["RobotR3.fullRobot","Animation"]

Use data to directly create models

Interpolate data

Pass data as input

Use real-world data as input to a model:

location=$GeoLocation;
data=AirTemperatureData[location,{Yesterday,Today}]

DateListPlot[data]

Shift the data to start at time 0 and filter out date objects:

cleanData=TimeSeriesShift[data,{{0}}]["Path"];

Create model from data:

outsideTempModel=CreateDataSystemModel["OutsideTemperature",cleanData]

Connect the model with other components

Connect with model of a room with heating:

comp={"outside"∈outsideTempModel,"fromKToC"∈"Modelica.Thermal.HeatTransfer.Celsius.FromKelvin","room"∈"DocumentationExamples.Modeling.RoomHeating.Room"};
conn={"outside.y[1]"\[UndirectedEdge]"fromKToC.Kelvin","fromKToC.Celsius"\[UndirectedEdge]"room.outsideTemp"};
model=ConnectSystemModelComponents["MyModel",comp,conn]

Extract first and last dates in the original data:

{firstDate,lastDate}={data["FirstDate"],data["LastDate"]}

Simulate for as long as there is data:

simInterval=DateDifference[firstDate,lastDate];
sim=SystemModelSimulate[model,simInterval]

Plot simulation results:

SystemModelPlot[sim,{"fromKToC.Celsius","room.temperatureSensor1.T"},PlotLegends->{"Outside","Inside"}]

Use data to fit

Turn your fits into a SystemModel

Extract historical currency exchange data:

gbpTousd=TimeSeries[MapAt[UnitConvert,First[WolframAlpha["GBP to USD from January, 2022 to today",{{"Result",1},"ComputableData"}]],{All,2}]]

Plot the data:

DateListPlot[gbpTousd]

Remove units and dates:

cleanGBPData=QuantityMagnitude[MapAt[DateDifference[gbpTousd["FirstDate"],#]&,gbpTousd["Path"],{All,1}]];

Fit a linear model to the data:

fit=LinearModelFit[cleanGBPData,x,x]

Create a SystemModel from the fit:

modelFit=CreateSystemModel["GBPToUSDFit",fit]

Simulate the model and plot it together with the data:

Show[
ListPlot[cleanGBPData],
SystemModelPlot[modelFit,{"y"},280,<|"Inputs"->{"u"->Function[t,t]}|>,PlotStyle->]
]

Calibrate parameters

Calibrate parameters in a model by comparing to measurement data:

dataCal=Import["ExampleData/inertiaTestData.dat"];

Set up a criteria function for model fitting:

modelMotor=;

criteriaFun[pv:{Verbatim[Rule][_?StringQ,_?NumericQ]..},var_?StringQ]:=Module[{f,tmax=40},f=SystemModelSimulate[modelMotor,tmax,<|"ParameterValues"->pv|>][{var}][[1]];
Norm[Table[dataCal[[i,2]]-f[dataCal[[i,1]]],{i,Length[dataCal]}]]^2]

Fit parameters to the test data:

fittedpars=FindMinimum[criteriaFun[{"Resistor1.R"->p,"SpringDamper1.d"->q},"Inertia2.w"],{{p,0.1},{q,0.25}},StepMonitor:>PrintTemporary[{p,q}]]

Simulate with the fitted parameters:

simCal=SystemModelSimulate[modelMotor,{0,40},<|"ParameterValues"->({"Resistor1.R"->p,"SpringDamper1.d"->q}/. fittedpars[[2]])|>]

Show the test data and the calibrated model together:

Show[SystemModelPlot[simCal,{"Inertia2.w"},PlotStyle->Orange],ListPlot[dataCal]]

Use data to train neural networks

Basic use-case

Train a net and create a model from it

Data found at

Dua, D. and Graff, C. (2019). UCI Machine Learning Repository [http://archive.ics.uci.edu/ml]. Irvine, CA: University of California, School of Information and Computer Science, https://archive.ics.uci.edu/ml/datasets/Energy+efficiency.

With citation request:

A. Tsanas, A. Xifara: ‘Accurate quantitative estimation of energy performance of residential buildings using statistical machine learning tools’, Energy and Buildings, Vol. 49, pp. 560-567, 2012, ISSN 0378-7788, https://doi.org/10.1016/j.enbuild.2012.03.003.

Data also found at

https://www.kaggle.com/datasets/elikplim/eergy-efficiency-dataset?select=ENB2012_data.csv

Download data in the notebook directory:

descriptions = {"Relative Compactness","Surface Area","Wall Area","Roof Area","Overall Height","Orientation","Glazing Area","Glazing Area Distribution"}->{"Heating Load","Cooling Load"};

dataFile=URLDownload["https://archive.ics.uci.edu/ml/machine-learning-databases/00242/ENB2012_data.xlsx",FileNameJoin[{NotebookDirectory[],"ENB2012_data.xlsx"}]];

Import data and check no data is missing:

dataNet=Values[Normal[SemanticImport[dataFile][[;;768,;;10]]]];

Flatten[Map[NumericQ,dataNet,{2}]]//DeleteDuplicates

Split data for training and validation:

cleanNetData=Thread[dataNet[[All,;;8]]->dataNet[[All,9;;]]];

train=cleanNetData[[;;Round[768*0.8]]];
validation=cleanNetData[[Round[768*0.8]+1;;]];

Create and initialize net:

elem="RectifiedLinearUnit" ;
size=10;
net=NetInitialize[NetChain[{
LinearLayer[2size,"Input"->8],
ElementwiseLayer[elem],
LinearLayer[2size],
ElementwiseLayer[elem],
LinearLayer[size],
LinearLayer[2]}],RandomSeeding->1234,Method->"Kaiming"]

Train net:

tNet=NetTrain[net,train,All,ValidationSet->validation,TimeGoal->2*60]

Extract trained net:

fNet=tNet["TrainedNet"]

Create model from net, placing ONNX file in notebook directory:

netModel=CreateSystemModel["MyNet",fNet,GeneratedAssetLocation->FileNameJoin[{NotebookDirectory[],"MyNet.onnx"}]]

Connect the model with other components

Connect net model to room model:

testModelstr=;
ImportString[testModelstr,"MO"]

Simulate and plot:

SystemModelPlot["NetWithRoom",{"temperature.T"},PlotRange->All]

Industry use-case

Analytical Chemistry
Nanoplasmonic Avidity-Based Detection and Quantification of IgG Aggregates

Thuy (Tran^1), Erik (Martinsson^2), Sergio (Vargas^3), Ingemar ((Lundström)^4), (((Carl-Fredrik Mandenius))^5), Daniel (Aili^1)

1 Laboratory of Molecular Materials, Division of Biophysics and Bioengineering, Department of Physics, Chemistry and Biology, Linköping University, 581 83 Linköping, Sweden 
2 ArgusEye AB, Spannmålsgatan 55, 583 36 Linköping, Sweden. 
3 Wolfram MathCore AB, Teknikringen 1E SE-583 30 Linköping, Sweden
4 Sensor and Actuator Systems, Department of Physics, Chemistry and Biology, Linköping University, 581 83 Linköping, Sweden 
5 Biotechnology, Division of Biophysics and Bioengineering, Department of Physics, Chemistry and Biology, Linköping University, 581 83 Linköping, Sweden 
 
Abstract
Production of therapeutic monoclonal antibodies (mAbs) are complex processes that require extensive bioanalytical characterization to ensure high and consistent product quality. Aggregation of mAbs is common and very problematic and can result in products with altered pharmacodynamics and pharmacokinetics, and potentially increased immunogenicity. Rapid detection of aggregates, however, remains very challenging using existing analytical techniques. Here, we show a real-time and label-free fiber optical biosensor system for specific detection and quantification of immunoglobulin G (IgG) aggregates using Protein A mediated avidity effects. Compared to monomers, IgG aggregates were found to have substantially higher apparent affinity when binding to Protein A functionalized sensor chips in a specific pH range (pH3.8 - 4.0). Under these conditions, aggregates and monomers showed significantly different binding kinetics. Using principal components analysis aggregate concentrations as low as 5% at a total IgG concentration of 0.45 mg/mL could be detected.  Reliable and rapid aggregate quantification was demonstrated with a limit of detection (LOD) and limit of quantification (LOQ) of about 9 µg/mL and 30 µg/mL, respectively. Using neural network-based curve fitting, it was further possible to simultaneously quantify monomers and aggregates with high correlation between predicted and control values and small prediction errors even for aggregate concentrations lower than 30 µg/mL. Our work demonstrates a unique avidity-based biosensor approach for fast aggregate analysis that can be used for rapid at-line quality control, including lot/batch release testing. This technology can also likely be further optimized for real-time in-line monitoring of product titers and quality, facilitating process intensification.

Reinforcement Learning

https://reference.wolfram.com/language/ref/device/SimulatedCartPole.html
https://github.com/openai/gym/blob/master/gym/envs/classic_control/cartpole.py
https://www.wolfram.com/language/12/neural-network-framework/train-an-agent-in-a-reinforcement-learning-environment.html
https://community.wolfram.com/groups/-/m/t/1379860

The ball and beam model

ClearAll[model,rs,vs,thetab,omegab];

SystemModel["DocumentationExamples.Control.BallAndBeam"]

model=SystemModel["DocumentationExamples.Control.BallAndBeam.BallAndBeamDynamics"];

torqueMagnitude=0.01;
alternatedInput=TimeSeries[Table[{k,torqueMagnitude(((-1))^(Round[k/0.1]))},{k,0,3,0.1}],ResamplingMethod->{"Interpolation",InterpolationOrder->0}];
SystemModelPlot[model,{"rs","T"},1.4,<|"InitialValues"->{"rs"->0},"Inputs"->{"T"->alternatedInput}|>,PlotRange->All,Epilog->{Dashed,LightGray,InfiniteLine[{{0,1/20},{1,1/20}}],InfiniteLine[{{0,-1/20},{1,-1/20}}]}]

Model preparation and net creation

eqs=SystemModel[model,"SystemEquations"];
pars=SystemModel[model,"ParameterValues"];
cleaneqs=eqs/.{\[FormalT]->t}/.pars/.{QuantityVariable[var_,_]:>ToExpression[var]};
nlssm=NonlinearStateSpaceModel[cleaneqs,{rs[t],vs[t],thetab[t],omegab[t]},{T[t]},{rs[t],thetab[t]},t];

policyNet=NetInitialize@NetChain[{64,Tanh,3,SoftmaxLayer[]},"Input"->4,"Output"->NetDecoder[{"Class",torqueMagnitude{-1,0,1}}]]

reinforceLoss=NetGraph[
{
CrossEntropyLossLayer["Index"],
ThreadingLayer[#1*#2&]
},
{
NetPort["DiscountedReward"]->2,
{NetPort["Input"],NetPort["Action"]}->1->2
},
"Action"->NetEncoder[{"Class",torqueMagnitude{-1,0,1}}]
]

policyNet[{0.001,0,0,0}]

Training function

discountedSum[rewards_List,gamma_]:=Dot[rewards,gamma^Range[0,Length[rewards]-1]];
applyDiscount[rewards_List,gamma_]:=Table[discountedSum[Drop[rewards,i-1],gamma],{i,Length@rewards}]

collectEpisodes[nlssm_,policyNet_,totalNumberSteps_Integer,maxStepsPerEp_Integer,discountFactor_]:=Module[
{action,actions,epActions,observedState,observedStates,epObservedStates,reward,rewards,epDiscountedRewards,discountedRewards,ended,currentSteps,initialValues,sim,samplePeriod,input,finalValsrs,beamLength},

epActions={};
epObservedStates={};
epDiscountedRewards={};

currentSteps=0;

(*episode loop*)
While[currentSteps<totalNumberSteps,

actions={};
observedStates={};
rewards={};

(*start each new episode with a fresh environment*)
beamLength=1;
initialValues={RandomReal[{-beamLength/200,beamLength/200}],0,0,0};
observedState=initialValues;

(*step loop*)
Do[
AppendTo[observedStates,observedState];
action=policyNet[observedState];

samplePeriod=0.1;
input=With[{iaction=action},Function[iaction]];
sim=SystemModelSimulate[nlssm,samplePeriod,<|"InitialValues"->Thread[{rs,vs,thetab,omegab}->observedState],"Inputs"->{T->input}|>,ProgressReporting->False];

observedState=sim[{rs,vs,thetab,omegab},samplePeriod];
finalValsrs=First[observedState];
ended=Or[finalValsrs<-beamLength/20,beamLength/20<finalValsrs];
reward=Boole[!ended];

(*Echo[{finalValsrs,reward}];*)

AppendTo[rewards,reward];
AppendTo[actions,action];
If[ended,Break[]];

,
Min[maxStepsPerEp,totalNumberSteps-currentSteps]
];

(*Echo[rewards];*)

discountedRewards=applyDiscount[rewards,discountFactor];

epActions=Join[epActions,actions];
epObservedStates=Join[epObservedStates,observedStates];
epDiscountedRewards=Join[epDiscountedRewards,discountedRewards];

currentSteps+=Length@actions;
];

<|
"Action"->epActions,
"Input"->epObservedStates,
"DiscountedReward"->epDiscountedRewards
|>
];
reinforceSampler[nlssm_,net_,batchSize_]:=Block[{data},data=collectEpisodes[nlssm,net[#,"RandomSample"]&,batchSize,50,0.99];
data];

reinforceSampler[nlssm,policyNet,3]

Training ( ~ hours)

results=NetTrain[policyNet,reinforceSampler[nlssm,#Net,#BatchSize]&,All,LossFunction->reinforceLoss,Method->"RMSProp",BatchSize->1024,MaxTrainingRounds->35,TrainingProgressMeasurements-><|"Measurement"->NetPort["DiscountedReward"],"Aggregation"->"Mean","Key"->"Reward"|>]

trained=results["TrainedNet"]

netPath=FileNameJoin[{NotebookDirectory[],"RLNet3.ONNX"}];
Export[netPath,trained]

Importing trained net (instead of training)

(*trained2=NetReplacePart[trained,"Output"->None];*)
netPath=FileNameJoin[{NotebookDirectory[],"RLNet3.ONNX"}];
trained2 =Import[netPath]

Creating model from net

CreateSystemModel["RLNet",trained2,GeneratedAssetLocation->FileNameJoin[{NotebookDirectory[],"MyRLNet.onnx"}]]

Create a closed loop system with the net

bbClosedLoopSystemString=;
ImportString[bbClosedLoopSystemString,"MO"]

Testing

SystemModelPlot["BBClosedLoopSystem",{"ballAndBeam.rs","ballAndBeam.T"},3,<|"InitialValues"->{"ballAndBeam.ballAndBeamDynamics.rs"->-0.0001}|>,PlotRange->All,Epilog->{Dashed,LightGray,InfiniteLine[{{0,1/20},{1,1/20}}],InfiniteLine[{{0,-1/20},{1,-1/20}}]}]

SystemModelPlot["BBClosedLoopSystem",{"ballAndBeam.rs","ballAndBeam.T"},3,<|"InitialValues"->{"ballAndBeam.ballAndBeamDynamics.rs"->0.0001}|>,PlotRange->All,Epilog->{Dashed,LightGray,InfiniteLine[{{0,1/20},{1,1/20}}],InfiniteLine[{{0,-1/20},{1,-1/20}}]}]

System Modeler at WTC 2022

Thursday (CST times)

05:00 - 05:30 pm Cambodia Wolfram Student Stars! Mathematica, Arduino, Raspberry Pi

Friday (CST times)

09:30 - 10:30 am Meetup: Virtual Labs in Education