Fault detection for a valve

This examples demonstrates how to use EstimationPy to identify faults in a valve by means of state and parameter estimation. This example has been presented in [Bonvini2014], please refer to this paper for further information.

The model

../../_images/valve.png

The considered system is a valve that regulates the water flow rate in a water distribution system. The system is described by the following equations

(1)\[\begin{split}\dot{m}(t) &= \phi(x(t)) A_v \sqrt{\rho(t)}\sqrt{\Delta p(t)}, \\ x(t) + \tau \dot{x}(t) &= u(t),\end{split}\]

where \(\dot{m}(\cdot)\) is the mass flow rate passing through the valve, \(\Delta p(\cdot)\) is the pressure difference across it, \(u(\cdot)\) is the valve opening command signal, \(x(\cdot)\) is the valve opening position, \(\tau\) is the actuator time constant, \(\phi(\cdot)\) is the power-law opening characteristic, \(A_v\) is the flow coefficient and \(\rho(\cdot)\) is the fluid density (please note that the square root of the pressure difference is regularized around zero flow in order to prevent singularities in the solution). The system has three sensors that respectively measure the pressure difference across the valve, the water temperature $T(cdot)$ and the mass flow rate passing through it. All the sensors are affected by measurement noise. In addition, the mass flow rate sensor is also affected by a thermal drift.

(2)\[\begin{split}T^{N}(t) &= T(t) + \eta_T(t) \\ \Delta p^{N}(t) &= \Delta p(t) + \eta_P(t) \\ \dot{m}^{N+D}(t) &= (1 + \lambda(T(t) - T_{ref})) \dot{m}(t) + \eta_m(t)\end{split}\]

The measurement equations are described in (2), where the superscript \(^N\) indicates a measurement affected by noise, the superscript \(^{N+D}\) indicates the presence of both noise and thermal drift, \(T_{ref}\) is the reference temperature at which the sensor has no drift, \(\lambda\) is the thermal drift coefficient and \(\eta_{T}(\cdot)\), \(\eta_{P}(\cdot)\), and \(\eta_{m}(\cdot)\) are three uniform white noises affecting respectively the temperature, pressure and mass flow rate measurements. These signals are sampled every two seconds.

Suppose during the operation, at \(t = 80 s\), the valve becomes faulty. The fault affects the ability of the valve to control its opening position. The valve opening cannot go below 20% (causing a leakage) and over 60% (it gets stuck). At \(t = 250 s\), the valve stuck position moves from 60% to 90%.

Filtering + smoothing

The fault identification procedure is asked to identify whether the valve not works as expected, that is its opening position follows the command signal. The fault identification is performed using the UKF that uses as input signals for its model the noisy pressure difference, the noisy water temperature and the command signal. The command signal is noise free because it is computed by some external controller and not measured. The UKF compares the output of its simulations with the measured mass flow rate that is affected by both noise and thermal drift. The effect of the thermal drift is visible in Figure (XX) where the green dots represent the measured mass flow rate while the green line is the actual mass flow rate passing through the valve.

The UKF and the smother estimate the value of the state variable \(x(t)\) representing the valve opening position and the parameter \(\lambda\), the thermal drift coefficient. The length of the augmented state is \(n=2\). Hence for every estimation step, the UKF performs \(1+2 \times 2=5\) simulations. The UKF has the initial conditions \(x(0) \sim N(1.0, 0.05)\) and \(\lambda \sim N(0, 0.7 \cdot 10^{-3})\), the output noise covariance matrix is \(R = [0.05]\), and the filter is parametrized by the default coefficients. Both the state and the parameter are constrained and their accepted ranges are \(x(t) \in [0, 1]\), and \(\lambda \in [-0.005, 0.025]\).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
    dir_path = os.path.dirname(__file__)
    
    # Define the path of the FMU file
    filePath = os.path.join(dir_path, "..", "..", "modelica", "FmuExamples", "Resources", "FMUs", "ValveStuck.fmu")
    
    # Path of the csv file containing the data series
    csvPath = os.path.join(dir_path, "..", "..", "modelica", "FmuExamples", "Resources", "data", "NoisyData_ValveStuck.csv")
    
    # Set the CSV file associated to the input, and its covariance
    input = m.get_input_by_name("dp")
    input.get_csv_reader().open_csv(csvPath)
    input.get_csv_reader().set_selected_column("valveStuck.dp")
    
    # Set the CSV file associated to the input, and its covariance
    input = m.get_input_by_name("cmd")
    input.get_csv_reader().open_csv(csvPath)
    input.get_csv_reader().set_selected_column("valveStuck.cmd")
    
    # Set the CSV file associated to the input, and its covariance
    input = m.get_input_by_name("T_in")
    input.get_csv_reader().open_csv(csvPath)
    input.get_csv_reader().set_selected_column("valveStuck.T_in")
    
    # Set the CSV file associated to the output, and its covariance
    output = m.get_output_by_name("m_flow")
    output.get_csv_reader().open_csv(csvPath)
    output.get_csv_reader().set_selected_column("valveStuck.m_flow")
    output.set_measured_output()
    output.set_covariance(0.05)
    
    
    #################################################################
    # Select the variable to be estimated
    m.add_variable(m.get_variable_object("command.y"))
    
    # Set initial value of parameter, and its covariance and the limits (if any)
    var = m.get_variables()[0]
    var.set_initial_value(1.0)
    var.set_covariance(0.05)
    var.set_min_value(0.0)
    var.set_constraint_low(True)
    var.set_max_value(1.00)
    var.set_constraint_high(True)
    
    #################################################################
    # Select the parameter to be estimated
    m.add_parameter(m.get_variable_object("lambda"))
    
    # Set initial value of parameter, and its covariance and the limits (if any)
    var = m.get_parameters()[0]
    var.set_initial_value(0.00)
    var.set_covariance(0.0007)
    var.set_min_value(-0.005)
    var.set_constraint_low(True)
    var.set_max_value(0.025)
    var.set_constraint_high(True)
    
    # Initialize the model for the simulation
    m.initialize_simulator()
    
    # Set models parameters
    use_cmd = m.get_variable_object("use_cmd")
    m.set_real(use_cmd, 0.0)
    Lambda = m.get_variable_object("lambda")
    m.set_real(Lambda, 0.0)
    
    #################################################################
    # Instantiate the UKF for the FMU model
    ukf_FMU = UkfFmu(m)
    
    # Start filter
    t0 = pd.to_datetime(0.0, unit = "s", utc = True)

In line 3 the model of the valve that is used for the state estimation is defined. Lines 6-18 associates inputs from CSV files to the inputs of the FMU, while lines 21-25 associate the measured output to its relative column in the CSV file.

In lines 30-39 the state variable \(x(t)\) is added as estimation variable, while lines 43-52 add \(\lambda\) to the parameters being estimated. For both state and parameter the covariance, initial values, upper and lower boundaries are defined.

In line 65 the object representing the UKF filter and smoother is instantiated using the valve FMu model, and in lines 68-71 the filter and smoothing algorithm is started.

Results

The figure below shows the temperature of the water and the pressure difference across the valve. The solid lines are the real values computed by a simulation while the dots are the noisy measurements used by the state estimation algorithm.

../../_images/Temp_Press.png

The image below shows the measured mass flow rate (green dots), the real mass flow rate (solid green line), and the estimated mass flow rate computed by the smoother. The green dots are far from the solid line because of the thermal drift of the flow sensor.

../../_images/Flow.png

The image below shows the opening position of the valve. The green sold line is the set point command \(u(t)\), while the solid red line is the actual valve opening position \(x(t)\), the unknown state of the system. The solid blue line shows the position estimated by the smoother and the blue shaded area is the confidence interval of the estimation \(\hat{x}(t)\).

../../_images/Positions.png

The image below shows how the UKF and the smoother estimate the unknown thermal drift coefficient. The solid green line represents the unknown value of the coefficient, while the red and blue line are the UKK and smoother estimation respectively. It is possible to see that the state and parameter estimation algorithm is able to identify the unknown drift coefficient.

../../_images/Drift.png

The estimated state, that is the actual position of the valve, is then converted into a probability of fault. The signal is compared with the desired position. Whenever the actual opening position is bigger than the desired value the algorithm determined that the valve is leaking, when the opening is smaller than the desired position the valve is determined to be stuck. The red colored area identifies the exact periods in which the faults were introduced in the simulation study. The black lines show the results of the identification using the algorithm.

../../_images/Probability.png