Fault detection

We'll consider a problem of identifying faults that have occurred in a system based on sensor measurements of system performance.

Topic references

Samar, Sikandar, Dimitry Gorinevsky, and Stephen Boyd. "Likelihood Bounds for Constrained Estimation with Uncertainty." Decision and Control, 2005 and 2005 European Control Conference. CDC-ECC'05. 44th IEEE Conference on. IEEE, 2005. <https://web.stanford.edu/~boyd/papers/pdf/map_bounds.pdf>__

Problem statement

Each of :math:n possible faults occurs independently with probability :math:p. The vector :math:x \in \lbrace 0,1 \rbrace^{n} encodes the fault occurrences, with :math:x_i = 1 indicating that fault :math:i has occurred. System performance is measured by :math:m sensors. The sensor output is

.. math::

\begin{equation} y = Ax + v = \sum_{i=1}^n a_i x_i + v, \end{equation}

where :math:A \in \mathbf{R}^{m \times n} is the sensing matrix with column :math:a_i being the fault signature of fault :math:i, and :math:v \in \mathbf{R}^m is a noise vector where :math:v_j is Gaussian with mean 0 and variance :math:\sigma^2.

The objective is to guess :math:x (which faults have occurred) given :math:y (sensor measurements).

We are interested in the setting where :math:n > m, that is, when we have more possible faults than measurements. In this setting, we can expect a good recovery when the vector :math:x is sparse. This is the subject of compressed sensing.

Solution approach

To identify the faults, one reasonable approach is to choose :math:x \in \lbrace 0,1 \rbrace^{n} to minimize the negative log-likelihood function

.. math::

\begin{equation} \ell(x) = \frac{1}{2 \sigma^2} |Ax-y|_2^2 + \log(1/p-1)\mathbf{1}^T x + c. \end{equation}

However, this problem is nonconvex and NP-hard, due to the constraint that :math:x must be Boolean.

To make this problem tractable, we can relax the Boolean constraints and instead constrain :math:x_i \in [0,1].

The optimization problem

.. math::

\begin{array}{ll} \mbox{minimize} & |Ax-y|_2^2 + 2 \sigma^2 \log(1/p-1)\mathbf{1}^T x\ \mbox{subject to} & 0 \leq x_i \leq 1, \quad i=1, \ldots n \end{array}

is convex. We'll refer to the solution of the convex problem as the relaxed ML estimate.

By taking the relaxed ML estimate of :math:x and rounding the entries to the nearest of 0 or 1, we recover a Boolean estimate of the fault occurrences.

Example

We'll generate an example with :math:n = 2000 possible faults, :math:m = 200 measurements, and fault probability :math:p = 0.01. We'll choose :math:\sigma^2 so that the signal-to-noise ratio is 5. That is,

.. math::

\begin{equation} \sqrt{\frac{\mathbf{E}|Ax |^2_2}{\mathbf{E} |v|_2^2}} = 5. \end{equation}

.. code:: python

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(1)

n = 2000
m = 200
p = 0.01
snr = 5

sigma = np.sqrt(p*n/(snr**2))
A = np.random.randn(m,n)

x_true = (np.random.rand(n) <= p).astype(int)
v = sigma*np.random.randn(m)

y = A.dot(x_true) + v

Below, we show :math:x, :math:Ax and the noise :math:v.

.. code:: python

plt.plot(range(n),x_true)

.. parsed-literal::

[<matplotlib.lines.Line2D at 0x11ae42518>]

.. image:: fault_detection_files/fault_detection_3_1.png

.. code:: python

plt.plot(range(m), A.dot(x_true),range(m),v)
plt.legend(('Ax','v'))

.. parsed-literal::

<matplotlib.legend.Legend at 0x11aee9630>

.. image:: fault_detection_files/fault_detection_4_1.png

Recovery

We solve the relaxed maximum likelihood problem with CVXPY and then round the result to get a Boolean solution.

.. code:: python

%%time
import cvxpy as cp
x = cp.Variable(shape=n)
tau = 2*cp.log(1/p - 1)*sigma**2
obj = cp.Minimize(cp.sum_squares(A*x - y) + tau*cp.sum(x))
const = [0 <= x, x <= 1]
cp.Problem(obj,const).solve(verbose=True)
print("final objective value: {}".format(obj.value))

# relaxed ML estimate
x_rml = np.array(x.value).flatten()

# rounded solution
x_rnd = (x_rml >= .5).astype(int)

.. parsed-literal::

ECOS 2.0.4 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS

It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
 0  +7.343e+03  -3.862e+03  +5e+04  5e-01  5e-04  1e+00  1e+01    ---    ---    1  1  - |  -  - 
 1  +4.814e+02  -9.580e+02  +8e+03  1e-01  6e-05  2e-01  2e+00  0.8500  1e-02   1  2  2 |  0  0
 2  -2.079e+02  -1.428e+03  +6e+03  1e-01  4e-05  8e-01  2e+00  0.7544  7e-01   2  2  2 |  0  0
 3  -1.321e+02  -1.030e+03  +5e+03  8e-02  3e-05  7e-01  1e+00  0.3122  2e-01   2  2  2 |  0  0
 4  -2.074e+02  -8.580e+02  +4e+03  6e-02  2e-05  6e-01  9e-01  0.7839  7e-01   2  2  2 |  0  0
 5  -1.121e+02  -6.072e+02  +3e+03  5e-02  1e-05  5e-01  7e-01  0.3859  4e-01   2  3  3 |  0  0
 6  -4.898e+01  -4.060e+02  +2e+03  3e-02  8e-06  3e-01  5e-01  0.5780  5e-01   2  2  2 |  0  0
 7  +7.778e+01  -5.711e+01  +8e+02  1e-02  3e-06  1e-01  2e-01  0.9890  4e-01   2  3  2 |  0  0
 8  +1.307e+02  +6.143e+01  +4e+02  6e-03  1e-06  6e-02  1e-01  0.5528  1e-01   3  3  3 |  0  0
 9  +1.607e+02  +1.286e+02  +2e+02  3e-03  4e-07  3e-02  5e-02  0.8303  3e-01   3  3  3 |  0  0
10  +1.741e+02  +1.557e+02  +1e+02  2e-03  2e-07  2e-02  3e-02  0.6242  3e-01   3  3  3 |  0  0
11  +1.834e+02  +1.749e+02  +5e+01  8e-04  9e-08  8e-03  1e-02  0.8043  3e-01   3  3  3 |  0  0
12  +1.888e+02  +1.861e+02  +2e+01  3e-04  3e-08  2e-03  4e-03  0.9175  3e-01   3  3  2 |  0  0
13  +1.909e+02  +1.902e+02  +4e+00  7e-05  7e-09  6e-04  1e-03  0.8198  1e-01   3  3  3 |  0  0
14  +1.914e+02  +1.912e+02  +1e+00  2e-05  2e-09  2e-04  3e-04  0.8581  2e-01   3  2  3 |  0  0
15  +1.916e+02  +1.916e+02  +1e-01  2e-06  3e-10  2e-05  4e-05  0.9004  3e-02   3  3  3 |  0  0
16  +1.916e+02  +1.916e+02  +4e-02  7e-07  8e-11  7e-06  1e-05  0.8174  1e-01   3  3  3 |  0  0
17  +1.916e+02  +1.916e+02  +8e-03  1e-07  1e-11  1e-06  2e-06  0.8917  9e-02   3  2  2 |  0  0
18  +1.916e+02  +1.916e+02  +2e-03  4e-08  4e-12  4e-07  5e-07  0.8588  2e-01   3  3  3 |  0  0
19  +1.916e+02  +1.916e+02  +2e-04  3e-09  3e-13  3e-08  5e-08  0.9309  2e-02   3  2  2 |  0  0
20  +1.916e+02  +1.916e+02  +2e-05  4e-10  4e-14  4e-09  6e-09  0.8768  1e-02   4  2  2 |  0  0
21  +1.916e+02  +1.916e+02  +4e-06  6e-11  6e-15  6e-10  9e-10  0.9089  6e-02   4  2  2 |  0  0
22  +1.916e+02  +1.916e+02  +1e-06  2e-11  2e-15  2e-10  2e-10  0.8362  1e-01   2  1  1 |  0  0

OPTIMAL (within feastol=1.8e-11, reltol=5.1e-09, abstol=9.8e-07).
Runtime: 6.538894 seconds.

final objective value: 191.6347201927456
CPU times: user 6.51 s, sys: 291 ms, total: 6.8 s
Wall time: 7.5 s

Evaluation

We define a function for computing the estimation errors, and a function for plotting :math:x, the relaxed ML estimate, and the rounded solutions.

.. code:: python

import matplotlib

def errors(x_true, x, threshold=.5):
    '''Return estimation errors.
    
    Return the true number of faults, the number of false positives, and the number of false negatives.
    '''
    n = len(x_true)
    k = sum(x_true)
    false_pos = sum(np.logical_and(x_true < threshold, x >= threshold))
    false_neg = sum(np.logical_and(x_true >= threshold, x < threshold))
    return (k, false_pos, false_neg)

def plotXs(x_true, x_rml, x_rnd, filename=None):
    '''Plot true, relaxed ML, and rounded solutions.'''
    matplotlib.rcParams.update({'font.size': 14})
    xs = [x_true, x_rml, x_rnd]
    titles = ['x_true', 'x_rml', 'x_rnd']

    n = len(x_true)
    k = sum(x_true)

    fig, ax = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12, 3))

    for i,x in enumerate(xs):
            ax[i].plot(range(n), x)
            ax[i].set_title(titles[i])
            ax[i].set_ylim([0,1])
            
    if filename:
        fig.savefig(filename, bbox_inches='tight')
        
    return errors(x_true, x_rml,.5)

We see that out of 20 actual faults, the rounded solution gives perfect recovery with 0 false negatives and 0 false positives.

.. code:: python

plotXs(x_true, x_rml, x_rnd, 'fault.pdf')

.. parsed-literal::

(20, 0, 0)

.. image:: fault_detection_files/fault_detection_10_1.png