Intelligence Layer of the Twin

Overview

With a calibrated model in place, the twin can start to reason. This day adds the services that make a twin intelligent: surrogate models that answer in milliseconds instead of seconds [1], hybrid models that let data correct the physics [2], detection of faults from the residuals between model and measurement [3], estimation of the remaining useful life [4], what-if scenarios, and a twin that plans the fan of the pack ahead of time. Each service is tested against the truth of the simulated pack.

Day at a glance

flowchart LR
  A["A calibrated core"] --> B["Surrogate models"]
  B --> C["Hybrid models"]
  A --> D["Fault detection"]
  D --> E["Remaining useful life"]
  A --> F["What-if and control"]
  E --> G["Your asset"]

A calibrated core

Every service of the day calls the calibrated model of Day 3. Two hours of simulated driving take a fraction of a second, which is fast enough for monitoring but too slow when a planner must try thousands of options, or when an engineer wants to explore a design space interactively.

The model describes the pack by its state, four numbers that change in time. They are the state of charge, abbreviated SoC, between 0 and 1, the cell temperature \(T\) in °C, the usable capacity \(Q\) in ampere-hours (Ah) and the internal resistance \(R\) in ohm. Three inputs act on the state: the current \(I\) of the load in ampere (A), the duty \(u\) of the cooling fan, from 0 for off to 1 for full speed, and the ambient temperature \(T_{\mathrm{amb}}\). The sign \(d/dt\), or a dot over a symbol, denotes a rate of change per second. The model has one equation for each state:

The first equation lets the charge fall with the current, and its factor 3600 converts ampere-hours into ampere-seconds. The second equation is an energy balance. The current heats the pack with the power \(I^2 R\) in watt (W), and the air removes the power \((hA_0 + hA_1 u)(T - T_{\mathrm{amb}})\). The difference, divided by the heat capacity \(m c_p\), is the rate at which the temperature rises. The conductances \(hA_0\) and \(hA_1\) are given in watt per kelvin (W/K) and the heat capacity in joule per kelvin (J/K). A step of one kelvin (K) equals a step of one degree Celsius.

The third equation is the ageing law. The capacity fades in proportion to the size of the current, \(|I|\), and faster when the pack is hot. Its exponential function \(\exp\) is called an Arrhenius factor, in which \(T + 273.15\) is the temperature in kelvin. The factor equals 1 at 25 °C, about 1.5 at 35 °C and about 2.2 at 45 °C. The fourth equation lets the resistance grow slowly once the capacity has fallen below its value when new, \(Q_0\). The parameter \(R_0\) is the resistance of the new pack: It is the value at which the state \(R\) starts, and it sets the speed of this growth.

The four calibrated values lie within 2 percent of the true values of the simulated pack of Days 1 to 3. Most trips of the day use the highway cycle of the course, a current that swings slowly between 34 and 46 A, \(I(t) = 40 + 6 \sin(2 \pi t / 1800)\) with the time \(t\) in seconds. The equations advance in steps of \(\Delta t\) seconds, written with the Greek capital letter delta, with the Runge-Kutta method of fourth order of Day 2. The timing of two hours uses steps of 1 s, and the services of the day use 5 to 20 s.

Two numbers summarise the thermal behaviour. With a constant current and fan duty, the temperature settles where heating and cooling balance, \(\dot{T} = 0\). This steady-state temperature \(T_{\mathrm{ss}}\), the time constant \(\tau\), the Greek letter tau, and the temperature at a time \(t\) after a start at \(T(0)\) follow from the second equation:

The time constant tells how fast the temperature approaches the steady state. After one time constant it has covered about 63 percent of the distance, and after three time constants about 95 percent.

On this day the pack on which the services are tested is itself a simulation with the calibrated parameters, to which a section adds a fault or a missing effect where it needs one. Every experiment therefore has a known right answer. Each service is also compared with a baseline, a simple alternative that costs almost nothing to build, such as a fixed alarm limit next to the cumulative sum or a fan that is always on next to the planning controller. A service earns its place only when it beats its baseline by a margin that justifies its complexity.

Python code in the Colab notebook, Section 1. Open Section 1 of the Colab notebook and run it. The variable CAL holds the calibrated parameters of the table above, and the cell prints the four of them that Day 3 determined. It then simulates two hours of the pack with simulate in steps of 1 s and prints the time that this takes, about a third of a second.

Surrogate models and their limits

A surrogate learns the input-output behaviour of the full model from a set of simulations and then answers new questions almost instantly. The surrogate of the project answers one question: how hot the pack gets during one hour of highway driving. The answer depends on three inputs, which form a query \(q\): the ambient temperature \(T_{\mathrm{amb}}\), the level of the current \(I_s\) and the fan duty \(u\). The simulation of the hour runs in steps of 5 s with the current \(I(t) = I_s\,(1 + 0.15 \sin(2 \pi t / 1800))\) and returns the highest temperature:

The planned set of simulations from which a surrogate learns is called a design of experiments. The design of the day draws \(n = 120\) queries \(q_1, \ldots, q_n\) at random, each input \(j\) uniformly between a lower bound \(\mathrm{LO}_j\) and an upper bound \(\mathrm{HI}_j\). These bounds form the training box, also called the design box. The simulated peak temperatures \(y_1, \ldots, y_n\) of these queries are the training outputs.

The baseline surrogate is a quadratic formula. Each input is first standardised, \(\tilde{q}_j = (q_j - m_j) / s_j\), with the mean \(m_j\) of input \(j\) in the training set and its standard deviation \(s_j\), the usual measure of how far values scatter around their mean. The prediction, marked with a hat, adds a constant, the three standardised inputs, their squares and their products:

The sign \(\sum\) means: Add the term for every value of the index. The ten coefficients \(w\) are chosen by least squares, the method of Day 3, with a small penalty on large coefficients that is called ridge regression.

A Gaussian process, abbreviated GP, is a good choice for a few inputs, because it also reports how uncertain each answer is [1]. It assumes no formula, only that queries close to each other have similar outputs. Closeness is measured by a kernel, a function \(\kappa(q, q')\), with the Greek letter kappa, that is large when the queries \(q\) and \(q'\) are close and falls towards zero as they move apart. The simplest kernel, which the animation of this section uses, is the squared exponential:

Here \(d\) is the distance between the two queries, in which the difference of every input \(j\) is measured in units of a length scale \(\ell_j\), a script letter l. The length scale is the distance over which the output changes noticeably. The amplitude \(\sigma_f\), with the Greek letter sigma, is the typical size of the deviations of the output from its mean. At \(d = 0\) the kernel equals \(\sigma_f^2\), and at three length scales it has fallen to about one percent of that. The surrogate of the day uses a relative with the same distance \(d\), the Matérn kernel with the smoothness 5/2:

From the kernel, the Gaussian process computes two numbers for a new query \(q\): a prediction \(\mu(q)\), with the Greek letter mu, and a standard deviation \(s(q)\) that states how uncertain the prediction is:

The bold letters are lists and tables of numbers. \(\mathbf{y}\) is the list of the \(n\) training outputs and \(\bar{y}\) their mean. \(\mathbf{k}(q)\) is the list of the kernel values \(\kappa(q, q_i)\) between the new query and the training queries, and the sign \(\top\) lays it on its side for the multiplication. \(\mathbf{K}\) is the table of the kernel values \(\kappa(q_i, q_j)\) between all pairs of training queries. \(\mathbf{I}\) is the identity table, with ones on its diagonal and zeros elsewhere, and the noise variance \(\sigma_n^2\) allows the outputs a small random error. The power \(-1\) stands for solving a system of \(n\) linear equations, which is done once during training.

In words: The prediction starts from the mean of the training outputs and adds a weighted sum of their deviations from it, in which the training queries close to \(q\) count most. The squared standard deviation, called the variance, starts from \(\sigma_f^2\), and every training query close to \(q\) takes a part of it away. Far from all training queries the prediction falls back to the mean. The standard deviation returns to \(\sigma_f\) there and never exceeds it, however wrong the prediction is.

The amplitude, the length scales and the noise variance are not set by hand. Training starts from guesses, 8 °C, 12 A and 0.4 for the three length scales, and adjusts them so that the training outputs become as probable as possible under the model, a principle called maximum likelihood.

Both surrogates are tested on 40 new random queries inside the box, with the root mean square error, abbreviated RMSE, and the largest single error:

Here \(y_i\) is the simulated peak temperature of test query \(i\), \(\hat{y}_i\) the prediction of the surrogate and \(n = 40\). The RMSE is a typical size of the error in °C, in which large errors count more than small ones: For three errors of 1, 1 and 4 °C it is \(\sqrt{(1 + 1 + 16) / 3} = 2.4\) °C. The quadratic formula has an RMSE of 0.76 °C and a largest error of 3.1 °C. The Gaussian process has an RMSE of 0.02 °C and a largest error of 0.11 °C.

The value of a surrogate is its speed, and its price is the design of experiments. The speed-up and the break-even point, the number of queries after which the surrogate has saved more time than its design cost, are

with the time \(t_{\mathrm{sim}}\) of one simulation, the time \(t_{\mathrm{sur}}\) of one prediction of the surrogate and the time \(t_{\mathrm{design}}\) of the whole design. In one run these are 32.13 ms, 0.14 ms and 3.6 s. The speed-up is \(32.13 / 0.14\), about 230: The surrogate predicts the peak temperature about 230 times faster than the simulation. Every query saves \(32.13 - 0.14 = 31.99\) ms, about 0.032 s, so a design of 3.6 s is paid back after about \(3.6 / 0.032 = 113\) queries. The unrounded times give 114. Timings change from run to run, so another run gives other values of a similar size.

Python code in the Colab notebook, Section 2. Open Section 2 of the Colab notebook and run it. The function peak_T is the simulation behind every query. The cell draws the 120 training queries Xd, simulates their outputs yd and fits the quadratic formula poly and the Gaussian process gp. It prints the RMSE and the largest error of both on 40 test queries, the time per query of the simulation and of the surrogate, and the break-even point. To see the fitted amplitude and length scales, add a cell with print(gp.kernel_).

The knowledge of a surrogate, however, ends at the edge of the range it was trained on. A test draws 25 random queries from each of three boxes, the training box and two boxes whose sides are 25 and 60 percent longer, and compares both surrogates with new simulations. Many of these queries still lie inside the training box, 14 and 5 of the 25, and the fan duty is kept between 0 and 1.

A small step outside costs the Gaussian process nothing, but in the largest box its error is eight times its error inside. Single queries show what its standard deviation is worth. With the current far outside the box, \(q = (30, 80, 0.5)\), it predicts 66.2 °C with a standard deviation of 1.0 °C, and the simulation returns 65.3 °C. Far outside on two inputs, \(q = (60, 90, 0)\), it predicts 148.1 °C with a standard deviation of 2.9 °C, and the simulation returns 156.9 °C. This error of 8.8 °C is three times the stated standard deviation.

The standard deviation grows outside the box, but it measures the distance from the training data and cannot know how the physics continues there. A surrogate must therefore refuse questions outside its box. A guard tests every input \(j\) of a query, \(\mathrm{LO}_j \le q_j \le \mathrm{HI}_j\), and passes the query to the simulation when the test fails. A real pack would have failed long before 150 °C, so the last query also lies outside the range in which the simulation itself can be trusted.

The animation shows the effect with one input and the kernel of the worked example. Its eight training points lie between 20 and 50 A, where the prediction follows the true curve within a few hundredths of a degree. At a query of 70 A the true value is 119.53 °C, the prediction 54.20 °C and the stated standard deviation 11.55 °C, so the error is 5.7 standard deviations. From 80 A on the standard deviation stays at 12.00 °C, the value of \(\sigma_f\), while the error keeps growing.

Python code in the Colab notebook, Section 3. Open Section 3 of the Colab notebook and run it. The cell prints the table above for the three boxes, and then the prediction and the standard deviation of the Gaussian process, labelled GP std, for three single queries. To compare the last of them with the simulation, add a cell with peak_T(60, 90, 0). Exercise A asks for the guard, a function in_box(x).

Hybrid models

A physics model is reliable but often incomplete; a data-driven model is flexible but knows nothing about physics. A hybrid model combines them: The physics model predicts, and a small data-driven model learns only the error that the physics leaves [2]. In the project, the physics model of the twin ignores one effect of the pack, and a hybrid model learns it from data while keeping the physics in charge.

A data-driven model, also called a machine learning (ML) model, is a flexible function whose shape is learned from examples. The pack of this experiment has a property that the model of the twin lacks: Its resistance rises with the temperature. In this pack, the resistance that produces the heat is

Every degree above 25 °C adds 1.2 percent of resistance and of heat, which makes 24 percent at 45 °C. The model without this effect is the physics model, and its prediction is written \(T_{\mathrm{phys}}\). Two learned models compete with it. The model called ML only learns the temperature \(T\) of the pack itself. The hybrid model learns the difference between the pack and the physics model and adds it to the physics prediction:

Both learned functions, \(g_{\mathrm{ML}}\) and \(g\), receive the same five inputs at every sample, collected in \(\phi\), the Greek letter phi. These inputs are the time \(t\), the ambient temperature, the fan duty, the current and \(T_{\mathrm{phys}}\). The functions differ in their training target, which is \(T\) for \(g_{\mathrm{ML}}\) and the difference \(T - T_{\mathrm{phys}}\) for \(g\).

Both are built by gradient boosting, the class GradientBoostingRegressor of scikit-learn, which adds up many small decision trees. A decision tree asks a few yes-or-no questions about the inputs, such as whether \(T_{\mathrm{phys}}\) exceeds 40 °C, and returns one constant for each combination of answers. Each new tree is fitted to the errors that the sum of the earlier trees still makes:

Here \(g_0\) is the mean of the training targets and \(\mathrm{tree}_m\) the tree number \(m\) of \(M = 100\) trees. The learning rate \(\nu = 0.1\), the Greek letter nu, lets each tree correct a tenth of the remaining error, and a tree asks at most three questions in a row. These are the default settings. A sum of constants has a consequence that decides this experiment: Beyond the training range every tree returns the same constant as at the edge of that range, so the output stays flat, however far the inputs move.

The training data are nine simulated highway trips of one hour on mild days, at 22, 26 and 30 °C with fan duties of 0, 0.5 and 1. With one sample every 10 s they hold 3240 samples. The test data are nine such trips in a heatwave, at 34, 40 and 46 °C. The test therefore asks every model about conditions outside its training data, where the twin of a pack is needed most. The table gives the RMSE against the pack.

On the mild days both learned models look excellent, and the physics model looks worst. In the heatwave the order changes. The error of the physics model roughly doubles. The ML-only model fails: Its trees saw temperatures up to 54.1 °C and never predict more than 54.0 °C, while the pack reaches 74.4 °C. The hybrid model has the smallest error, because the physics part carries the prediction into the new conditions. Only its correction is limited to the largest value that it learned, 4.7 °C, while the true difference grows to 9.0 °C.

The error in the heatwave divided by the error on the mild days is largest for the hybrid model, about 96, because its error on the mild days was tiny. A ratio of two errors says little. The number that counts is the absolute error in the conditions in which the model will be used.

Correcting the output of a physics model is one of several ways to combine physics and data. A second way gives a learned model inputs that come from physics, as Python step 4 does with the square of the current. A third way trains a neural network with a penalty for violating the physical equations [7].

Python code in the Colab notebook, Section 4. Open Section 4 of the Colab notebook and run its two cells. The function real_derivatives plays the pack with the temperature-dependent resistance, and build_dataset simulates the nine trips of the mild days, tr, and of the heatwave, te. The first cell fits the two learned models ml and hyb and prints the table above with the additional column degradation_x. The second cell draws the three predictions against the pack for one mild day and one day of the heatwave.

Detecting faults

The difference between what the sensor reads and what the twin predicts, the residual, is close to noise while the asset is healthy. A fault makes it drift. For sample number \(i\), with the reading \(y_i\) and the prediction \(\hat{y}_i\) of the twin, the residual is

A positive residual means that the pack is warmer than the twin expects. The simplest alarm is a fixed limit, which fires when a single value, the temperature or the residual, exceeds a threshold. A fixed limit reacts only to large jumps and ignores a deviation that stays below the threshold, however long it lasts. A cumulative sum, abbreviated CUSUM, the method of Page, adds up small deviations over time and raises an alarm when the sum crosses a decision limit [3]:

The allowance \(k\) is subtracted from every residual. Ordinary noise is mostly smaller than \(k\) and pulls the sum down, and the function max with 0 keeps it from falling below zero. Only a deviation that stays above \(k\) lets the sum grow, sample after sample, until it passes the decision limit \(h\). Both constants have the unit of the residual. The detector of the project uses \(k = 0.15\) °C and \(h = 4\) °C, with one sample every 10 s. The samples are numbered \(i\) in this section, because the letter \(k\) is the usual name of the allowance.

The two constants trade speed against false alarms, which are alarms without a fault. When a fault shifts the residual by a constant amount \(\delta\), the Greek letter delta, every sample adds \(\delta - k\) to the sum, and the alarm comes after about

samples. A small allowance and a small limit give early alarms and many false ones, and a large allowance makes the detector blind to every shift below it. This sum also watches only for readings above the prediction. A pack that is colder than predicted needs a second sum, in which the sign of the residual is reversed.

In the animation the fault is a slow drift of 0.036 °C per minute under noise with a standard deviation of 0.3 °C. With an allowance of 0.35 °C and a decision limit of 6, the cumulative sum raises its alarm at minute 66 and the fixed limit at minute 85, and the readout counts no false alarm. With an allowance of 0.05 °C and a decision limit of 1 the alarm comes at minute 53, but eight false alarms precede the fault.

The test of the project is a trip of 90 minutes at 35 °C with the fan at half duty, and the pack changes after 3000 s. The twin predicts the temperature with the unchanged model, so the residual shows the effect of the fault. The fixed limit of this section is a temperature of 46 °C, and a healthy trip peaks at 44.4 °C.

For the degraded fan the residual grows by about 0.035 °C per sample, and the sum crosses the decision limit 190 s after the fault. At that moment the pack is only 0.7 °C warmer than predicted, and its temperature of 43.9 °C still lies inside the range of a healthy trip. The limit of 46 °C is reached 690 s later. On the pack the cumulative sum therefore warns 11.5 minutes earlier than the fixed limit, and 12.2 minutes earlier for the blocked coolant path.

The table follows the sum through the samples after this fault. The first four residuals stay below the allowance of 0.15 °C and leave the sum at zero. From the fifth sample on, every residual adds its excess over the allowance: \(0.165 - 0.150 = 0.015\) at the fifth sample, and \(0.341 - 0.150 = 0.191\) at the tenth. The 20th sample, at 3190 s, lifts the sum from 3.867 to \(3.867 + 0.677 - 0.150 = 4.394\), above the decision limit of 4.

The second fault raises no alarm at all, and the peak temperature stays at 44.4 °C, the value of a healthy trip. The reason lies in the model of the pack. The heat comes from the state \(R\), while the parameter \(R_0\) sets only the value at which \(R\) starts and the speed of its slow growth. Raising \(R_0\) during a trip therefore leaves the heat, the temperature and the residual practically unchanged. A detector sees a fault only through its effect on the signal that it watches.

A twin should therefore state which faults its detectors cannot see, and add a second signal for them. For the resistance this signal is the voltage of the pack, which is the voltage at rest minus the drop \(I R\) across the internal resistance: At 40 A a resistance of 0.0456 ohm takes 1.8 V. A residual also depends on the model as much as on the pack, so a twin whose parameters are out of date raises alarms about its own error. Detectors that learn the normal behaviour from data follow the same idea in industrial twins [8]. Part C of the interactive lab computes a cumulative sum by hand.

Python code in the Colab notebook, Section 5. Open Section 5 of the Colab notebook and run it. The function inject_fault simulates the trip with one of the three faults, twin_residual returns the residuals, and cusum returns the sum and the index of the alarm. The cell prints the table above, in which cusum_alarm_s is the alarm of the cumulative sum and raw_threshold_s the time at which the temperature passes 46 °C. To follow the sum sample by sample, add a cell with d = inject_fault("fan degradation"); s, a = cusum(twin_residual(d)); print(s[300:320]). Exercise B asks for the warning gained.

Remaining useful life

The capacity of the pack falls slowly with use. Extrapolating the trend of the capacity to the end-of-life level, usually 80 percent of the original, gives the remaining useful life [4]. The interval of the estimate is wide early and narrows as the history grows, and a responsible twin reports the interval, not only the number. A narrow interval, however, is not yet a correct one, as the history of the pack below shows.

A quantity that tells how worn an asset is, here the capacity \(Q\), is called a health indicator. The end-of-life level, abbreviated EOL, is \(Q_{\mathrm{EOL}} = 0.8 \times 60 = 48\) Ah. The remaining useful life, abbreviated RUL, is counted in trips. The simulated history of the pack holds 180 trips of 90 minutes. Each trip has its own ambient temperature, drawn at random around 32 °C. It is an urban trip in stop-and-go traffic with a probability of 0.6 and a highway trip otherwise, and the fan follows the fixed rule of the next section.

A highway trip costs about 24 mAh of capacity and an urban trip about 12 mAh, where 1 mAh, a milliampere-hour, is a thousandth of an ampere-hour. After each trip the capacity is measured once, with a random error whose standard deviation is 0.05 Ah, about three times the loss of an average trip. The estimate fits a straight line to the measured capacities by least squares and follows it down to the end-of-life level:

Here \(n\) is the number of a trip and \(\hat{Q}(n)\) the capacity that the line gives for it. The intercept \(a\) is the capacity of the line at trip 0, and the slope \(b\) is the change of the capacity per trip, a negative number. \(n_{\mathrm{EOL}}\) is the trip at which the line reaches the end-of-life level, and \(n_{\mathrm{last}}\) is the number of the latest trip. The trips are numbered from 0, so a history of 180 trips ends with \(n_{\mathrm{last}} = 179\).

One number is not enough, because other measurements would give another line. The bootstrap estimates this uncertainty from the data themselves [9]. From the points of the history it draws as many points at random with replacement, so that some points appear several times and others not at all, fits the line to this resample and computes its RUL. This is repeated 200 times. It reports the median of the 200 values, the value in the middle when they are sorted, as the estimate. The 5th and the 95th percentile, the values below which 5 and 95 percent of them fall, bound the 90 percent interval.

The last column adds the number of the latest trip to the estimate, for example \(89 + 710 = 799\), so that estimates made at different times can be compared.

In the left panel the measured capacity falls from 60 Ah to about 57.1 Ah within 180 trips. The dashed line fitted to all of them loses 16.0 mAh per trip and reaches the dotted end-of-life level near trip 750. The right panel repeats the estimate for histories of 20 to 180 trips. Its band is the 90 percent interval, which is 165 trips wide after 20 trips and 10 trips wide after 180. If the predicted end of life stayed the same, the curve would fall by one trip for every trip that passes. It rises between 30 and 70 trips of history, and after 110 trips it falls by almost two trips for every trip that passes.

The three intervals contradict each other. The history of 90 trips places the end of life between trips 781 and 817, and the history of 180 trips between trips 745 and 755. Continuing the same simulated history to 800 trips shows the truth: The capacity falls below 48 Ah in trip 737, outside both intervals. Only the first and widest interval, trips 697 to 768, contains it.

The bootstrap is not at fault. It answers one question correctly: how much the fitted line depends on the particular points that were measured. It cannot know what the data do not show. In this history the share of the more damaging highway trips happens to be 31 percent in the first 90 trips and 40 percent afterwards, so the early line is too flat. Real cells add a second reason: Their fade often accelerates after a knee, a bend in the capacity curve that a straight line fitted before it cannot foresee. A bootstrap interval covers the uncertainty of the fit. It does not cover a model of the wrong form or a future use that differs from the past, and both caveats belong next to every remaining-useful-life number.

Python code in the Colab notebook, Section 6. Open Section 6 of the Colab notebook and run its two cells. The function fleet_history simulates the 180 trips, and rul_estimate fits the line to 200 bootstrap resamples. The first cell prints the estimate and its 90 percent interval after 40, 90 and 180 trips and tests whether the end of life predicted after 180 trips lies in the interval predicted after 90 trips. The second cell draws the figure of this section. To see the true end of life, add a cell with h = fleet_history(n_trips=800); print((h.Q_true < 48).idxmax()).

What-if scenarios and control

A twin can run the future before it happens. What-if scenarios simulate decisions that have not been taken, for example a hotter day or a heavier load. A twin in the loop goes one step further: Before every period, it simulates several fan settings over the next minutes and applies the one with the best balance of temperature, energy and ageing. In the project it beats a fan that is always on, at lower energy, and keeps the pack below its limit.

Each scenario simulates 120 trips of 90 minutes one after another and carries the capacity from trip to trip. The fan follows a fixed rule, a thermostat that starts the fan at the temperature \(T_{\mathrm{on}}\) and reaches full duty 6 °C above it:

The function clip cuts its first argument off at the two bounds, so the duty stays between 0 and 1. The baseline scenario has \(T_{\mathrm{on}} = 42\) °C, the highway cycle and an ambient temperature of 35 °C. From the capacity \(Q_N\) after \(N = 120\) trips, a scenario extrapolates the number of trips until 20 percent of the capacity \(Q_0\) are lost, and it adds up the energy of the fan in kilowatt-hours (kWh):

The sum runs over all time steps \(k\) of all trips, \(u_k\) is the fan duty in step \(k\), and \(\Delta t = 20\) s is the length of a step. The division by \(3.6 \times 10^6\) converts joule into kilowatt-hours. The 120 trips last 180 hours.

The baseline loses 5.05 percent of its capacity in 120 trips, so \(n_{\mathrm{EOL}} = 0.20 / 0.0505 \times 120\), about 476 trips. A summer that is 8 °C hotter shortens the life to 408 trips. It keeps the pack above 46 °C for 171 of the 180 hours, although the fan uses almost three times the energy. Urban duty nearly doubles the life. Starting the fan at 36 °C buys about 60 trips for about twice the fan energy, and disabling the cooling costs about a quarter of the life. Such exchange rates between energy and life are what the manager of a fleet needs, and only a model can provide them before the summer.

A twin in the loop turns such simulations into decisions, an idea known as model predictive control, abbreviated MPC. The controller of the project acts every \(\Delta t = 20\) s. For each of six candidate fan duties, 0, 0.2, 0.4, 0.6, 0.8 and 1, it simulates the next 300 s with the calibrated model, holding the duty constant, and adds up a cost:

Here \(T_j\) is the simulated temperature after step \(j\) of the 15 steps, \(Q_{15}\) the simulated capacity after the last of them, \(Q_{\mathrm{now}}\) the present capacity and \(\Delta t / 60\) the length of a step in minutes. The first term punishes every degree above 44 °C, squared so that a large excess counts much more than a small one. The second term is the energy of the fan, and the last term the capacity that the pack loses within the 300 s. The weights are \(w_{\mathrm{hot}} = 1\), \(w_{\mathrm{fan}} = 0.02\) and \(w_{\mathrm{age}} = 40000\). The controller applies the duty with the lowest cost for one step of 20 s and then plans again from the new state, which is called a receding horizon.

Four controllers drive the same trip of 90 minutes on the highway cycle at 35 °C. They are no cooling, a fan that is always on, the fixed rule with \(T_{\mathrm{on}} = 42\) °C and the twin in the loop. They are ranked by one number for the whole trip, an objective that puts three costs on one scale:

Here \(H_{44}\) is the number of hours spent above 44 °C, \(E_{\mathrm{fan}}\) the fan energy in kWh and \(\Delta Q\) the capacity lost in mAh. With these weights one hour above the limit costs as much as 0.12 kWh of fan energy or 30 mAh of capacity.

The twin in the loop has the lowest objective. It keeps the fan off while the pack is cool and starts it at about 40 °C. It never lets the temperature pass 42.6 °C, and it uses 58 percent of the energy of the fan that is always on. It beats the fixed rule by 51 percent. The fixed rule reaches full duty only at 48 °C, so the pack settles above 44 °C and stays there for more than an hour. The twin beats the always-on fan by only 7 percent: It saves \(500 \times (0.0375 - 0.022) = 7.8\) points of fan energy and pays \(2 \times (23.28 - 21.63) = 3.3\) points of additional ageing.

Two remarks keep this result in proportion. The always-on fan is a trivial baseline that needs no model and cannot be wrong about one, and a margin of 7 percent has to be weighed against the new ways to fail that a planner brings. In this test the pack is also simulated with the same model that the planner uses, so its look-ahead is exact. A real pack differs from its model, and the margin can shrink.

The weights are not facts of engineering. They state how much an hour of overheating is worth in energy and in battery life, a trade-off that belongs to the people who answer for operation and safety. Writing the weights down makes this decision visible and open to question. The planner has weights of its own, and they need not agree with those of the ranking. The weights of the planner can be converted to the units of the ranking. Its fan term charges 0.02 points for every watt and minute, and a kilowatt-hour is \(1000 \times 60 = 60000\) watt-minutes, which gives \(0.02 \times 60000 = 1200\) points per kilowatt-hour. Its ageing term charges 40000 points per ampere-hour, which is 40 points per milliampere-hour. The ranking charges 500 and 2 points for the same quantities. The planner therefore values a milliampere-hour at \(40 / 1200 = 0.033\) kilowatt-hours and the ranking at \(2 / 500 = 0.004\), so the planner protects the capacity more than the ranking asks for.

The form of a cost matters as much as its weight. The ranking counts the hours above 44 °C, and the planner squares the excess. The control room in Part A of the interactive lab uses degree-hours, the excess above 44 °C added up over time, with the temperature \(T_k\) in step \(k\) and the step \(\Delta t\) in seconds:

One degree-hour is one degree above the limit for one hour, or two degrees for half an hour. A cost that only counts whether the limit is exceeded gives a planner no reason to cool once the limit is lost within its horizon. Degree-hours and squared excesses reward every degree of cooling.

Python code in the Colab notebook, Section 7. Open Section 7 of the Colab notebook and run its cells. In the first cell the function scenario simulates 120 trips, and the cell prints the table of the five scenarios. In the second cell rollout_cost is the cost \(J(u)\) of the planner and run_controller drives one trip with one of the four controllers. The cell prints the table of the controllers with the column objective, and the third cell draws their temperatures and fan duties and prints the two margins. Exercise C asks for the margin over the fan that is always on.

Your asset

The application of the day lets the cooling of each of four assets lose half of its effect. The twin compares its own prediction with the reading once per window and watches the residuals with a cumulative sum and with a fixed limit. All four assets share one thermal model with a single state, the temperature \(T\), which advances with the explicit Euler method of Day 2:

Here \(C\) is the heat capacity in J/K, the quantity called \(m c_p\) in the model of the pack, and \(P(t)\) is the heat input in watt. The index \(k\) counts the time steps of length \(\Delta t\).

The fan follows the fixed rule of the previous section, with \(T_{\mathrm{on}}\) set 6 °C below the limit of the asset, and acts on a noisy reading. The noise has a standard deviation of one percent of the distance between the ambient temperature and the limit, 0.11 °C for the pack. After 40 percent of the simulated time the conductance \(hA_1\) of the asset drops to half, while the twin keeps the full value. At the end of every window of 30 steps the twin stores the residual, the reading minus its prediction, and then synchronises: It replaces its temperature by the reading, as the shadow of Day 1 did.

The residuals are standardised. Each residual \(r_j\) of window \(j\) is divided by \(s_0\), the standard deviation of the residuals of the first 30 percent of the windows, in which the asset is still healthy. A standardised residual \(z_j\) is therefore measured in units of its normal scatter, and the same constants serve every asset:

The cumulative sum raises its alarm when \(S_j > 5\), and the fixed limit when a single residual has \(|z_j| > 4\). For the battery pack the fault begins at minute 48. The residuals of the next four windows are 2.76, 2.94, 0.71 and 3.37 standard deviations, all below the fixed limit of 4. The cumulative sum adds each of them minus 0.5: It rises to 2.26, then to 4.71 and 4.92, and the fourth residual lifts it to \(4.92 + 3.37 - 0.5 = 7.79\), above its limit of 5. The alarm therefore comes at the end of the fourth window, 9.5 minutes after the fault. Over the rest of the run the residuals lie around 3.5 standard deviations. The fixed limit follows 22.0 minutes after the fault, when one residual happens to reach 5.1.

A detector must also stay quiet when nothing is wrong. In 80 runs without the fault, the four assets with the seeds 0 to 19, the cumulative sum raises a false alarm in only one run and the fixed limit in six. A seed is the number that starts the random number generator of a run.

Python code in the Colab notebook, Section 8. Open Section 8 of the Colab notebook. Its first line is a switch, ASSET = "battery". Replace the word by server_rack, motor or greenhouse and run the section. The function watch_asset simulates the asset and its twin and returns the standardised residuals, the cumulative sum and the two alarm times. The cell prints the delay of each alarm after the fault and draws both signals with their limits. Exercise D asks for the delay of the cumulative sum in seconds.

Going further (optional)

Large models of fields such as temperature over a whole battery can be reduced: Proper orthogonal decomposition finds the few patterns that carry most of the behaviour and simulates only those [10]. For the four-state model of the pack, two patterns already carry about 97 percent, and this section shows why such a reduction does not pay off for a model this small.

Proper orthogonal decomposition, abbreviated POD, starts from snapshots, states of the model that are stored during simulations. Twelve trips of 90 minutes, with four ambient temperatures and three fan duties, provide the snapshots, with the four states stored every 10 s. The 6492 snapshots form the rows of a table with four columns. Each column is standardised: Its mean is subtracted and the result is divided by its standard deviation, so that the state of charge counts as much as the temperature. The standardised table \(\mathbf{X}\) is split by the singular value decomposition, abbreviated SVD:

The rows of \(\mathbf{V}^{\top}\) are the modes, four fixed patterns of the four states, and the table \(\mathbf{U}\) tells how strongly each mode is present in each snapshot. The singular value \(\sigma_j\) measures how much mode \(j\) contributes over all snapshots. The singular values are sorted from the largest to the smallest, and their list is called the spectrum. The share of the total variance that the first \(m\) modes carry is called their cumulative energy, a borrowed word that has nothing to do with kilowatt-hours:

The singular values of the pack are 143.17, 68.80, 25.68 and 8.90. Their squares add up to about 25970, of which the first contributes \(143.17^2 = 20498\), so one mode carries 79 percent. Two modes carry 97.2 percent and three modes 99.7 percent.

The left panel shows the spectrum on a logarithmic axis: Each singular value is between a third and a half of the one before it. The right panel shows the cumulative energy, with a dashed line at 99.9 percent that only all four modes reach. A model reduced to two modes would follow two numbers in place of four and lose about 3 percent of the variance. Saving two numbers per step does not justify that loss, so the analysis stops at the spectrum. The same analysis pays off for models with thousands or millions of states, such as the temperature at every point of a pack. Finding that a method brings no benefit is a result as well, and it belongs in the report.

Python code in the Colab notebook, Section 9. Open Section 9 of the Colab notebook, which is optional, and run it. The cell stores the snapshots of twelve trips in the table S, standardises them and calls np.linalg.svd. It prints the four singular values and the cumulative energy of one to four modes and draws the two panels of the figure of this section.

Python step 4: Plots and the scikit-learn interface

Python code in the Colab notebook, right after Section 0 (setup). Open the notebook and run this step cell by cell: It consists of short pieces of Python code with their explanations, a quick check and three exercises. Topics: Data and a labelled scatter plot · Training and test data · Fit, predict, measure · The value of a physical feature.

The notebook continues with the hands-on sections, and the interactive lab holds three practice parts and the self-assessment. The study path, the daily task and the research assignment are on the day overview.

Review cards

Select a card to turn it over.

References

[1] Rasmussen, C. E., & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press.

[2] Karniadakis, G. E., Kevrekidis, I. G., Lu, L., Perdikaris, P., Wang, S., & Yang, L. (2021). Physics-informed machine learning. Nature Reviews Physics, 3(6), 422-440.

[3] Page, E. S. (1954). Continuous inspection schemes. Biometrika, 41(1-2), 100-115. https://doi.org/10.1093/biomet/41.1-2.100

[4] Lei, Y., Li, N., Guo, L., Li, N., Yan, T., & Lin, J. (2018). Machinery health prognostics: A systematic review from data acquisition to RUL prediction. Mechanical Systems and Signal Processing, 104, 799-834.

[5] Hunter, J. D. (2007). Matplotlib: A 2D graphics environment. Computing in Science & Engineering, 9(3), 90-95.

[6] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., & Duchesnay, E. (2011). Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12, 2825-2830.

[7] Raissi, M., Perdikaris, P., & Karniadakis, G. E. (2019). Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378, 686-707. https://doi.org/10.1016/j.jcp.2018.10.045

[8] Rodriguez-Aguilar, R., Marmolejo-Saucedo, J. A., & Kose, U. (2024). Development of a digital twin driven by a deep learning model for fault diagnosis of electro-hydrostatic actuators. Mathematics, 12(19), 3124.

[9] Efron, B., & Tibshirani, R. J. (1994). An Introduction to the Bootstrap. Chapman & Hall/CRC.

[10] Benner, P., Gugercin, S., & Willcox, K. (2015). A survey of projection-based model reduction methods for parametric dynamical systems. SIAM Review, 57(4), 483-531.