## Building a Physics-Based Model of the Refrigerator - Thermostat Modeling (Part I)

There are all kind of mathematical models: some are based on physical laws, some are derived empirically and some are in-between. There are various pros and cons in choosing a modeling approach. Physics-based models, for example, are very concise but they are notoriously hard to construct (many problems, for example, come from combining equations for different components). Statistics-based models such as the ones from machine-learning can be often automatically or semi-automatically built but they are often imprecise or incorrect. One of the goals of the refrigeration project is to determine which modeling approach is better for diagnostics.

We will start building a physics-based model by analyzing the behavior of each component in an refrigerator. Let's start with the thermostat, which as we will see is one of the most complicated components.

The thermostat is the "brains" of the refrigerator. It is responsible for switching the compressor on and off and keeping the temperature in the refrigerator close to a value set by the user. One of the criteria for evaluating the performance of a thermostat is the variance in the temperature. A good thermostat should keep the temperature almost constant. The control problem is, of course, subject to constraints as a good thermostat should not switch on and off the compress too often.

There are different types of thermostats: the main distinction is if they are mechanical or electrical. Each thermostat includes a sensor, a switch and some control mechanism (or software) that actuates the switch based on the sensor readings. Figure 1 shows the type of thermostats used in our refrigerators.

The type of the temperature sensor that is used by our refrigerators is very similar to the old mercury thermometers. As mercury is hazardous, the thermostat sensor is filled with gas. The type of gas was not identified for the purpose of this study. The gas is often nitrogen.

Modern refrigerators use semiconductor temperature sensors, micro-controllers, and relays. The problem with these thermostats is their complexity: a refrigerator controlled with such a sensor can fail due to a software bug. They also need separate power supplies for the micro-controllers and the sensors which wastes additional energy.

Our refrigerators use old-style mechanical thermostats. Part of the disassembled thermostat is shown in figure 2. The thermostat has a temperature adjusting knob, a sensing tube, and a spring-mechanism for shortening and opening the electrical terminals that start and stop the compressor. There is a single screw for calibrating the temperature. The thermostat comes factory-calibrated.

This type of thermostat is a relatively complex thermal/mechanical/electrical component. Consider, for example, the debouncing mechanism shown in figure 3. It consists of a spring-like mechanism that prevents sparkling due to small gap between the conductive plates.

Figure 4 plots data from three different refrigerators. The left plot shows at what temperature the thermostat switches the compressor on for various possible settings of the control-knob (higher knob number means colder). The right plot shows at what temperature the thermostat switches the compressor off for the full range of control-knob settings. This difference in on and off temperatures comes from the mechanical design of the thermostat.

The three thermostats shown in figure 4 behave in a very different manner. The first thermostat is the best, where the temperature of the refrigerator is linear as a function of the control-knob setting. The second thermostat has a very strange non-linear behavior where positions one and two are indistinguishable and positions four and five are also very similar. The third thermostat, when set in its coolest regime does not switch the compressor off, a behavior different from the other two refrigerators. That is the reason why the box-plots in the second row of figure 4 are missing their right-most bars (the box-plots visualize switching behavior and this thermostat does not switch). This makes the behavior of the second thermostat, set to the lowest temperature, hard or impossible to distinguish from a short-circuited thermostat.

## Bodies, Gravity, and Simulation - Part I

$\mathbf{F}_{ij} = \frac{Gm_im_j(q_j - q_i)}{||q_j - q_i||^3}$

I am writing this blog entry for my astronomy term project (not saying which one because the professor my be reading this, but after I finish the course I will remove this note and say a few good words).

Generations of mathematicians have been frustrated because they failed to find an analytic solution to the $n$-body ($n \ge 3$) problem represented by the above equation.

The $n$-body problem is easy to solve by using numerical methods. In order to use numerical integration methods one has to convert the original problem to an Ordinary Differential Equation (ODE). Here is how Newton's law of gravity looks for $n$ bodies as an ODE:

$\frac{d^2x_i(t)}{dt^2} = G\sum_{k = 1, k \ne i}^{n}{\frac{m_k(x_k(t) - x_i(t))}{||x_k(t) - x_i(t)||^3}}$

This is a second order ODE. To solve it with Python we first have to convert the second order ODE to a system of first-order ODEs. Let's substitute $\frac{dx}{dt}$ with $p$ and $x$ with $q$. We now get:

$\frac{dq}{dt} = p$

$\frac{dp}{dt} = G\sum_{k = 1, k \ne i}^{n}{\frac{m_k(q_k(t) - q_i(t))}{||q_k(t) - q_i(t)||^3}}$

Let us first solve the two-body problem in two dimensions. This, of course, can be done analytically but we want to go through the steps of writing and evolving a Python program. So, here are the four equations:

$\frac{dq_1}{dt} = p_1$

$\frac{dq_2}{dt} = p_2$

$\frac{dp_1}{dt} = \frac{Gm_1(q_2 - q_1)}{||q_2 - q_1||}$

$\frac{dp_2}{dt} = \frac{Gm_2(q_1 - q_2)}{||q_1 - q_2||}$

And this is the Python code that can solve the equations:

#!/usr/bin/env python

import numpy
import scipy
import sys
import matplotlib.pyplot

import scipy.integrate

G = 6.674e-11
m1 = 1e2
m2 = 1e2

def func(Y, t):
p1 = numpy.array([ Y[0], Y[1] ])
p2 = numpy.array([ Y[2], Y[3] ])
q1 = numpy.array([ Y[4], Y[5] ])
q2 = numpy.array([ Y[6], Y[7] ])

p1_prime = G * m2 * (q2 - q1) / \
numpy.sum((q2 - q1)**2)
p2_prime = G * m1 * (q1 - q2) / \
numpy.sum((q1 - q2)**2)
q1_prime = p1
q2_prime = p2

return [ p1_prime[0], p1_prime[1],
p2_prime[0], p2_prime[1],
q1_prime[0], q1_prime[1],
q2_prime[0], q2_prime[1] ]

def main():
t = numpy.arange(0, 20000.0, 1)
init_p1 = [ 0, 1 ]
init_p2 = [ 0, -1 ]
init_q1 = [ 0, 0 ]
init_q2 = [ 10, 0 ]
iv = init_p1 + init_p2 + init_q1 + init_q2
result = scipy.integrate.odeint(func, iv, t)

matplotlib.pyplot.scatter(result[:, 4],
result[:, 5],
s = 1,
lw = 0,
c = 'blue')
matplotlib.pyplot.scatter(result[:, 6],
result[:, 7],
s = 1,
lw = 0,
c = 'red')
matplotlib.pyplot.show()

if __name__ == '__main__':
main()


So what do we set for the masses of the objects ($m_1$ and $m_2$) and for the initial position and velocity. Let's try with very small masses ($m_1 = m_2 = 10$) and the two bodies going in different directions with the same speed. This is the result:

So, given the initial vector, the above plot looks correct. The two bodies have small mass and as they keep flying in the direction in which they were flying initially. The next check is to see if when we increase the mass of the two bodies they will start attracting each other:

Now we are going to change the initial velocity of one of the two bodies and this is the result:

Finally, we can play some more and get this really complicated looking motion:

Why is all this interesting? Simulation is a very important tool in astronomy and the $n$-body problem served as a fundamental building block of today's understanding of the formation of the Solar system. Of course, Newton's formulas work only for small speeds and we have to add relativity into the picture if we want a more accurate prediction. But before updating our equations we want to generalize this introductory example first into 2-bodies in three dimensions and then to $n$ bodies in three dimensions. This is going to happen in our next post.

## Computing Features from the Refrigeration Sensor Readings

Features in machine learning are signals containing information about a measurable property of the system being observed.

It is difficult to use raw sensor data in machine-learning. The data may be too much: the refrigerator, for example, uses an AC compressor and is connected to a 60 Hz power supply. The Nyquist-Shannon sampling theorem says that to recover the 60 Hz, one has to sample at at-least 120 Hz. In reality, most engineers over-sample in an attempt to capture higher frequencies that would preserve transients in the signal. In the refrigeration case-study, the sampling rate is 500 Hz which means that at every 2 ms there is a new data point. Even at this modest sampling frequency it may be computationally hard to apply the classification or regression formula. We are only interested in algorithms that can compute diagnosis at real-time. Summarizing the signals over time-windows solves the computational difficulty problem.

Another reason for using features is that classification and regression machine-learning formulas are typically not good enough to capture the complex physics of the devices we are trying to diagnose. Manually constructing features is bridging physics and model-based reasoning and machine learning. We can actually "hide" a fully-fledged model-based diagnosis engine as a feature. This feature will directly compute if the device is failing and then there is no work for the classifier to determine the state of the system.

Table 1: Features formulas for diagnosing the refrigerator
Description Formula
Minimum value $\min{\{x_i, x_{i - 1}, \ldots, x_{i - k}\}}$
Maximum value $\max{\{x_i, x_{i - 1}, \ldots, x_{i - k}\}}$
Mean value $\displaystyle\frac{1}{k}\sum_{j = i - k}^{i}{x_j}$
Standard deviation $\sqrt{\sum{x - \mu}}$
Derivative $x_i - x_{i - 1}, \ldots, x_{i - k + 1} - x_{i - k}$
Product $x \times y$

Table 1 shows the types of features that are used for diagnosing the refrigerator. They are all computed for a sliding window of $k$ samples. The first four features are simply common sliding-window descriptive statistics. The fifth one is the first derivative of temperature (of course, first a low-pass filter with a cut-off frequency of 1 mHz is used to smooth the original signal). The last feature is the product of two signals. Notice that after computing derivatives or products, it is necessary to again compute descriptive statistics (minimum, maximum, mean, and standard deviation) over a sliding window.

In the case of the thermostat only the present value of the thermostat at the time of diagnosis is used. This approach does not use all thermostat information. To use all information it is possible to compute a feature that shows when was the last thermostat transition. Somewhat surprisingly, this feature not only does not help with any of the classifiers we have tried, its inclusion significantly decreases the isolation accuracy.

How is the length $k$ of the sliding window determined? A more advanced approach would be to use intelligent segmentation instead of a sliding window, however segmentation implies making a step toward computing a diagnosis which imposes a bootstrapping problem. In our study we just take a range of sliding window lengths and we compute a large number of features. We also compute the product of all possible temperature pairs. This results in a large number of features: 10 681, to be precise.

Figure 1 shows the features computed from one temperature sensor for a fixed sliding window size of 30 min.

Machine learning often normalizes the feature ranges by multiplying them with suitable constants. Experiments showed that this does not help the accuracy of classification in the case of the refrigerator and we use the features without scaling.

## Common Refrigerator Test-Bed Voltage, Current, and Power Plots

Figure 1 shows 24 hours of current consumption for a nominal experiment. The thermostat is set at position 1. The refrigerator is empty and at the set-point temperature at the beginning of the experiment. The raw signal from which figure 1 is constructed contains $4.32\times10^7$ values. During this one day experiment, the thermostat switches the compressor on and off 68 times. The mean time for the compressor being on is 194.8 s with a standard deviation of 11.1 s. The mean off-time of the compressor is 1064.3 s with a standard deviation of 25.7 s. The duty-cycle of the refrigerator for this experiment is 15.5%.

Refrigerators exhibit oscillatory behavior: there is the main compressor on/off cycle and the electrical circuit is powered by a 60 Hz AC. Figure 2 shows a frequency plot of the signal shown in figure 1. The plot shows amplitudes in the frequency range 0-0.025 Hz. The largest non-DC amplitude is at $0.788\times10^{-3}$ Hz. This frequency corresponds to a period of 1269.5 s which is within 99.1% of the mean duration of the compressor on/off cycle.

Figure 3 gives us a closer look on the current signal during one compressor-on cycle. The compressor is switched-on for 172 s. In the beginning the PTC thermistor is at room temperature and has low resistance. A few seconds after that the PTC thermistor heats-up due to the current flowing through it, its resistance increases and the current that flows through the start-winding of the compressor decreases. After the start-winding of the motor is switched-off the current signal becomes a simple sinusoid of 60 Hz and constant amplitude.

If we zoom-in the current plot even further we can see the start-up of the motor as shown in figure 4. It illustrates the start-up of the compressor by running current through a special start-up winding. This particular type of refrigerator uses a Resistance Start, Induction Run (RSIR) motor. The RSIR design uses a PTC thermistor as a time delay mechanism for spinning the single-phase AC motor.

Figure 4 is important in modeling the dynamics of the PTC. The PTC is a resistor whose resistance is heavily dependent on the temperature. To emphasize the fact that it behaves similar to an ordinary (mostly linear) resistor we show the compressor voltage during the start-up of a cycle. This is done in figure 5.

The sampling of the voltage and the current is synchronized by a common clock in the power logger that we have designed and use. To verify the correct working of this power logger we compute the product of the voltage and the current signals which is the power. The resulting compressor start-up power is consistent with how this type of AC electrical motors and PTC work and is shown in figure 6.

Figure 7 shows an FFT plot of the start-up current. The 60 Hz AC mains frequency is the dominant frequency which is no surprise. The Nyquist frequency of our power logger is 250 Hz and the shape of the AC sine is well-visible. There are aliased frequencies of small amplitude at 120 Hz and 180 Hz. There are several FFT components at 80 Hz and 200 Hz which are due to the choice of the ADC (they also appear in the voltage frequency plot).

Figure 8 shows an FFT plot of the start-up voltage. The voltage is almost perfectly sinusoidal which is consistent with the working of the AC motor. In our instrumentation design we use two simultaneously clocked ADCs of the same type (LTC2440). There are small frequency components (that we can digitally filter) which are due to the ADCs. The amplitudes of these frequencies (80 Hz and 200 Hz), however, are so small that there is no need of filtering for successful diagnostic analysis. These frequencies will be filtered by the diagnostic algorithm capabilities to deal with noise and error.