Calibration and Validation of System Models

This talk will showcase the recently added functionality to calibrate, validate and visualize uncertainty in system models. Applications and examples are explored.

System Modeler

System Modeler

Multiple domains:

SystemModel["Modelica","Children"]

Modeling based on connections:

SystemModel["RC"]

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[,"Modeling"]

SystemModelPlot[,{"capacitor.v"}]

Components and Equations

An inverted pendulum in System Modeler:

SystemModelPlot[,#,<|"ParameterValues"->{"electricalMotor.J"->0},"Inputs"->{"V"->UnitBox}|>,PlotRange->All]&/@{"pendulum.pendulumDamper.phi_rel","pendulum.sliderConstraint.s"}//GraphicsRow

An inverted pendulum with equations:

eqs={
(mc+mp)x''[t]-mp l Cos[θ[t]] θ''[t]+mp l Sin[θ[t]] ((θ'[t])^2)==k i[t]-dc x'[t],
θ''[t]mp (4(l^2) /3+(a^2)/4)-mp l x''[t] Cos[θ[t]]==mp g l Sin[θ[t]]-dp θ'[t],
0==vin[t]- r i[t]-lin i'[t]-k x'[t]
};

parVals={g->9.80665,r->13,lin->1,mc->445*(0.1^3),mp->7700*π*(0.005^2)*0.61,l->0.61/2,a->0.005,k->3.7*157.48*1.6,dc->2,dp->0.01};

opssm=Sequence[{{i[t],0},{θ[t],0},{θ'[t],0},{x[t],0},{x'[t],0}},{{vin[t],0}}];

outputs={θ[t],x[t]};

nssm=NonlinearStateSpaceModel[eqs,opssm,outputs,t]//Simplify

SystemModelPlot[nssm,#,20,<|"ParameterValues"->parVals,"Inputs"->{vin->UnitBox}|>,PlotRange->All]&/@{θ,x}//GraphicsRow

Data in System Modeling

SystemModelCalibrate

SystemModelCalibrate is used to calibrate parameter values in a system model to match simulation results with real-world data.

Calibrate a Model using Data

Retrieve the data:

vdata="{{t1, v1}, {t2, v2},...}";

Calibrate the parameter value:

SystemModelCalibrate[
<|variableName->vdata|>,
"model",
"pars"
]

SystemModelCalibrate[
<|variableName->vdata|>,
{"model",parameterConstraints},
{{parameterName,initialGuess}},
"spec",
"prop"
]

?SystemModelCalibrate

Closed-Flow Heating System

Use temperature data to estimate the heat flow rate for a burner in a model of a heating system:

hmodel = ;

tdata=QuantityArray[,{"Seconds","Kelvins"}];

csmodel=SystemModelCalibrate[
<|"heater.mediums[1].T"->tdata|>,
{hmodel,1900<"burner.Q_flow"<2100},
{"burner.Q_flow"},
"CalibratedSystemModel"
]

Compute the single prediction bands and plot them together with the data and the calibrated simulation response:

csmodel["SinglePredictionBandsPlot"]//First

Plot the residuals against time:

ListPlot[Transpose[{tdata[[All,1]],First[csmodel["CalibrationDataResiduals"]]}],Frame->True,Filling->Axis,FrameLabel->{"t","residual"},PlotLabel->"Residual vs Time"]

Biochemical Reaction

Use multiple measurements of two substances’ concentrations to calibrate and validate a model of the chemical reaction:

concentrationsB=TemporalData[];

concentrationsC=TemporalData[];

Select the model, specify the reaction coefficients as calibration parameters and introduce parameter constraints:

bioModel=;

reactionCoefs={"k1","k2","k3"};

constraints=0.01<"k1"<10&&0.01<"k2"<10&&0.01<"k3"<10;

Aggregate the data points of most of the measurements to use as calibration data:

caliCsB=SortBy[Join@@concentrationsB["Paths"][[;;-2]],First];

caliCsC=SortBy[Join@@concentrationsC["Paths"][[;;-2]],First];

Calibrate the model, indicating the data for the concentration variables and the simulation interval:

csmodel=SystemModelCalibrate[
<|"B.c"->caliCsB,"C.c"->caliCsC|>,
{bioModel,constraints},
reactionCoefs,
<|"SimulationInterval"->{0,20}|>,
"CalibratedSystemModel"
]

Extract the plot of the single prediction bands:

csmodel["SinglePredictionBandsPlot"]

Compute the root mean squared errors, using the remaining data path as validation data:

csmodel["RMSE",ValidationSet-><|"B.c"->concentrationsB["Paths"][[-1]],"C.c"->concentrationsC["Paths"][[-1]]|>]

Model Simplification

Recreate the behavior of a speaker by calibrating a simpler circuit using simulation results as input data:

speaker = ;

rl = ;

Compare the behavior of both models before calibration:

sims=SystemModelSimulate[speaker];

simrl=SystemModelSimulate[rl];

SystemModelPlot[{sims,simrl},{"vs.v"}]

Retrieve simulation results as data and use them for calibration:

datav=sims["RawData","vs.v"];

csmodel=SystemModelCalibrate[
<|"vs.v"->datav|>,
{rl,0.001<l<0.5},
{l},
"CalibratedSystemModel"
]

Compare the calibrated model with the speaker:

SystemModelPlot[{sims,csmodel["CalibratedSimulationData"]},{"vs.v"}]

SystemModelUncertaintyPlot

SystemModelPlot and SystemModelManipulate

SystemModelPlot lets you plot curves:

motorModel=;

SystemModelPlot[motorModel,{"Inertia2.w"},<|"ParameterValues"->{"Inertia2.J"->2}|>]

SystemModelPlot[motorModel,{"Inertia2.w"},<|"ParameterValues"->{"Inertia2.J"->{1,2,3}}|>]

sim=SystemModelSimulateSensitivity[motorModel,10,{"Inertia2.J"},<|"ParameterValues"->{"Inertia2.J"->2}|>]

SystemModelPlot[sim,{{"Inertia2.w","Inertia2.J",0.5}}]

[motorModel,10,{"Inertia2.J"},{"Plots",{"Inertia2.w"}}]

Validate a Model visually

SystemModelUncertaintyPlot is typically used to visualize the uncertainty of key signals in a system model resulting from uncertainty in parameters, initial values or noisy inputs.

The results of multiple simulations for a given variable are aggregated as temporal data for which quantiles are computed for each time point.

SystemModelUncertaintyPlot[
motorModel,
<|"Outputs"->"Inertia2.w","ParameterValues"->{"Inertia2.J"->}|>
]

Trajectories Near Equilibrium

Plot the uncertainty in the angle and angular velocity of a simple damped pendulum when the initial conditions take values near equilibrium points:

model=CreateSystemModel["DampedPendulum",{ω'[t]==-g Sin[θ[t]]/l-b ω[t],θ'[t]==ω[t]},t,{θ∈"Units.SI.Angle"},<|"ParameterValues"->{g->10,l->1/2,b->2}|>]

When the system starts at rest at small angles, it evolves toward the stable ground position:

SystemModelUncertaintyPlot[model,<|"SimulationInterval"->5,"Outputs"->θ,"InitialValues"->{θ->CenteredInterval[0,π/4]}|>,]

When the system starts near the unstable fully inverted position, it also evolves toward the ground position:

SystemModelUncertaintyPlot[model,<|"SimulationInterval"->5,"Outputs"->θ,"InitialValues"->{θ->Interval[{3π/4,π-π/100}]}|>,]

System Tolerance

The performance of a circuit depends heavily on its components and parameter tolerances. When the resistance of one of the components in this model of a speaker follows a truncated normal distribution, the variation in the current going through the speaker can be measurable:

SystemModelUncertaintyPlot[,<|"Outputs"->"speaker.l.i","SimulationInterval"->(Quantity[40, "ms", milliseconds, "Milliseconds"]),"ParameterValues"->{"r1.R"->TruncatedDistribution[{0,∞},NormalDistribution[6,5]]}|>]

External disturbances can also affect the performance of a circuit. For instance, the input voltage of the speaker can be affected by high-frequency signals:

inputs=Table[With[{k=k},Function[{t},2.5Sin[2 π t 50]+k Sin[2 π t 2000]]],{k,0,0.6,0.05}];

Plot[Through[inputs[t]],{t,0,0.01}]

Plot the current of the speaker when affected by these high-frequency disturbances:

SystemModelUncertaintyPlot[,<|"Outputs"->"speaker.l.i","SimulationInterval"->(Quantity[10, "ms", milliseconds, "Milliseconds"]),"Inputs"->{"vi"->inputs}|>,TargetUnits->{"Milliseconds","Milliamperes"}]

Validation of Controlled Systems

Validate a controlled ball and beam system by placing the ball away from the stable position:

The system model:

model=;

The farther away the ball starts from the equilibrium point, the more the system will struggle bringing it back to the stable position:

spec=<|"SimulationInterval"->{0,1.2},"InitialValues"->{"syse.sys.ballAndBeamDynamics.rs"->Interval[{0,0.15}]}|>;

SystemModelUncertaintyPlot[model,Append[spec,"Outputs"->{"syse.sys.ballAndBeamDynamics.rs"}],PlotRange->All,PlotLabel->"r"]

The control effort peaks at the start and then later a few more times to stop the ball from overshooting the stable position:

SystemModelUncertaintyPlot[model,Append[spec,"Outputs"->{"syse.sys.T"}],PlotRange->All,PlotLabel->"Control effort"]

Controlled systems should also be validated against external disturbances. Stress the controlled ball and beam system with a white noise torque disturbance and plot the control effort as the system experiences various sampling periods of the noise:

SystemModelUncertaintyPlot[model,<|"SimulationInterval"->{0,1.2},"Outputs"->{"syse.sys.T"},"Inputs"->{"T"->WhiteNoiseProcess[0.1]}|>,SamplingPeriod->0.3,PlotRange->All,PlotLabel->"Control effort"]

Orbital Maneuver

https://en.wikipedia.org/wiki/Hohmann_transfer_orbit

System Modeler has a model already:

SystemModeler["Control.CircularOrbit","Animation"]

Create a model for the 2 boosts in the maneuver:

pulsesModel=CreateSystemModel["OrbitalManeuver.Pulses",{y[t]==Piecewise[{{f1,t<tp},{f2,tp+tn<t<2tp+tn}},0]},t,{y∈"RealOutput"},<|"ParameterValues"->{tp->1,tn->4244.4-1}|>]

SystemModelPlot[pulsesModel,{y},(Quantity[3, "h", hours, "Hours"]),<|"ParameterValues"->{f1->300,f2->200}|>,TargetUnits->{"Hours",Automatic}]

Connect the models:

model=ConnectSystemModelComponents["OrbitalManeuver.BoostModel",{"orbit"∈"Control.CircularOrbit","zero"∈"Blocks.Sources.Constant","boost"∈"Pulses"},{"orbit.fr"\[UndirectedEdge]"zero.y","orbit.fphi"\[UndirectedEdge]"boost.y"},]

Simulate the model with different magnitudes of the boosts:

sim=SystemModelSimulate[model,{"orbit.x","orbit.y","orbit.dr"},<|"ParameterValues"->{"f1"->{0,120,140,160},"f2"->{0,120,140,160}}|>];

SystemModelPlot[sim,Axes->False,Frame->True]

Plot the change in the distance to the center of force:

SystemModelPlot[sim,{"orbit.dr"}]

Create a function that basically does what we hope for the maneuver:

target[t_]=Piecewise[{{0.65 t/4320,t<4320}},0.65]

Plot[target[t],{t,0,35000}]

Create data from the function:

data=Table[{t,target[t]},{t,0,32400,32400/1500}];

Use the data to calibrate the boosts:

csmodel=SystemModelCalibrate[
<|"orbit.dr"->data|>,
{model,10<f1<200,10<f2<200},
{{f1,100},{f2,100}},
<|"CalibratedModelName"->"OrbitalManeuver.CalibratedOrbit"|>,
"CalibratedSystemModel"]

Plot the results:

Show[{SystemModelPlot[csmodel["CalibratedSimulationData"],{"orbit.dr"},TargetUnits->{"Seconds",Automatic}],ListPlot[data,PlotStyle->Darker[Green]]}]

Put together the original and the calibrated model:

comparisonModel=ConnectSystemModelComponents["OrbitalManeuver.Comparison",{"original"∈model,"calibrated"∈csmodel["CalibratedModel"]},{},]

See  the animation:

SystemModeler[comparisonModel,"Animation"]

Plot the trajectories:

SystemModelPlot[comparisonModel]

Plot the variations in the radius when the boosts are in a given interval:

uPlot=SystemModelUncertaintyPlot[csmodel["CalibratedModel"],<|"Outputs"->{"orbit.dr"},"ParameterValues"->{"f1"->CenteredInterval[csmodel["CalibratedParameters"]["f1"],10]}|>]

Create a manipulate:

man=[model,{{f1,103},90,120},{{f2,93.9},90,120},{"Plots",{"orbit.dr"}},PlotStyle->Orange];

Overlap the manipulate and the uncertainty plot:

MapAt[Show[{uPlot,#}]&,man,{1}]