# Development of Liquid Capacity Measuring Algorithm Based on Data Integration from Multiple Sensors

Authors: Kiwoong Park; Si-Kyoung Lee; Hyeon Cheol Kim
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

This research proposes an algorithm using a process of integrating data from multiple sensors to measure the liquid capacity in real time regardless of the position of the liquid tank. The algorithm for measuring the capacity was created with a complementary filter using a Kalman filter in order to revise the level sensor data and IMU sensor data. The measuring precision of the proposed algorithm was assessed through repetitive experiments by varying the liquid capacity and the rotation angle of the liquid tank. The measurements of the capacity within the liquid tank were precise, even when the liquid tank was rotated. Using the proposed algorithm, one can obtain highly precise measurements, and it is affordable since an existing level sensor is used.


---

## Body

## 1. Introduction

Continuous or discontinuous measurements are used to assess the capacity of liquid contained in vehicles, tanker trucks, and industrial tanks. Continuous measurements are divided into the following types: buoyancy, weight measurement, pressure, capacitance, ultrasound, and radiation. Discontinuous measurements are divided into the following types: conductivity, heat transfer, capacitance, optical, ultrasound, and microwave [1–3]. Some sensors used for measuring liquid level are summarized in Table 1.Table 1
Comparison of sensors used for measuring liquid level.
Sensor type
Scope of measurement [m]
Precision
Temperature limits [°C]
Pressure limits [kgf/cm2]
Pressure type [9]
0.1–100
±0.05 [FS%]
70
—
Float type [10]
Customizable
6.35 mm
105
—
Displacement type [11]
0.3–3
±2 [%]
370
300
Capacitance type [12]
0.1–4
±0.1 [%]
200
102
Weight level type [13]
30.5
10 [mm]
49
—
Radiation type [14]
7
0.5 [%]
50
No limitation (noncontact)
Ultrasonic type [15]
8–11
±0.5 [%]
70
3
Radar type [16]
30
±10
By the strength of windowsThe pressure level sensor, a double-continuous measuring type, forms the core of the measurement market owing to its affordability, simplicity of its gauge, and its ability to work at high temperature and high pressure. The ultrasound- and radar-type sensors are more efficient than others in measuring the methods but these sensors are not economical [4, 5].Most demand for current level sensors comes from water treatment plants (sewage, industrial water, wastewater, and black water), a national key industry, and the sector that measures the flow rate for vehicles and vessels. The level sensor is especially adequate for measuring the flow rate of vessels and the fuel of vehicles.The existing method for measuring the vehicle fuel involves measuring poles and floating box [6]. This is structurally limited because it is immeasurable when the fuel tank of the vehicle is below 10%. In addition, measuring errors occur because the fuel gets heeled over when the vehicle turns or accelerates as shown in Figure 1.Figure 1
Disadvantages of current fuel-measuring methods for vehicles. (a) Structural limits result in inability to measure the remaining 10% of fuel. (b) Fuel-measuring error on slopes. (c) Fuel-measuring error when accelerating.
(a)
(b)
(c)This paper proposes a new fuel-measuring algorithm using a pressure level sensor that can be applied for fuel control and measurements of vehicles. The system can accurately measure the fuel capacity below 10% by improving mechanical limits of the existing method in measuring the capacity. Furthermore, the measurement errors are minimized by continuously measuring and revising the slopping or heeling of fuel when the vehicle turns or accelerates. This revised algorithm guarantees precision for laser- and ultrasound-type sensors. In addition, this algorithm is highly economical because an existing pressure level sensor [7] has been used.


## 2. Sensor Design

The proposed sensor should have a characteristic of a measuring algorithm of density and capacity for fuel measurement and a revised method in accordance with the fuel tank’s condition, acceleration, and acceleration of gravity.Figure2 demonstrates the system which this research has proposed. A pressure level sensor is used to measure capacity, density, and tilting angle, and an inertial measurement unit (IMU) is used to revise the measuring value in accordance with the rotation, acceleration, and acceleration of gravity [8]. The unique approach of this paper is that the relative value in accordance with the location of the level sensor and the rotation angle data of the IMU are both used to measure the position of a fuel tank.Figure 2
Algorithm applied in this assignment.The reasons are as follows. The IMU outputs a precise measurable value for the rotation angle and acceleration measurements, but its disadvantage is that such values are accumulated with time by an integration constant. In the case of a level sensor, it can measure the tilting angle and acceleration in accordance with the measuring value difference located structurally, but the measurements may not be accurate.The proposed algorithm of this paper guarantees precision and credibility because it uses the Kalman filter to realize a complementary filter with precision despite the error accumulation of the IMU and the imprecision and nonaccumulation of errors of the level sensor.


## 3. Sensor Measurement Algorithm for Precise Measurements of Sensors

### 3.1. Measuring Algorithm of Capacity

As demonstrated in Figure3, the capacity can be computed in accordance with the container volume and the input capacity of the liquid.Figure 3
Measuring method for capacity.The mass can be computed using the pressure applied to the bottom of a container with liquid. First, the calculation equation of the capacity is as follows:(1)
Pressure  applied  to  the  bottom  of  the  container
∝
liquid  level
×
liquid  density
.The equation to compute the force applied to the bottom in accordance with the liquid is as follows [17]:(2)
F
=
ρ
×
A
h
×
g
g
:  Gravitational  acceleration
.Therefore the pressureP applied to the unit area of the bottom side can be defined as in the equation below [2]:(3)
P
=
ρ
×
g
×
h
If  there  is  no  internal  pressure  in  the  tank
,
P
-
P
0
=
ρ
×
g
×
h
If  there  is  an  internal  pressure  in  the  tank
.In this control system, the following method has been used to measure the capacity.The sensor attached toward thex-axis and y-axis from the center of the target container is located on the opposite side of the container. The obtained measurement value is used to calculate the average value of each axis, and the capacity is computed by comparing the accumulated data from experiments: (4)
A
=
Sensor
1_out
+
Sensor
2_out
2
B
=
Sensor
3_out
+
Sensor
4_out
2
C
=
S
e
n
s
o
r
5
_
o
u
t
+
S
e
n
s
o
r
6
_
o
u
t
2
(5)
A
+
B
+
C
3
=
The  average  value  of  the  output  voltage
⟶
The  measuring  comparison  by  the  input  data
.Here the output value of each sensor does not accumulate, but the precision is low. This sensor output value changes with turning and accelerating and is used in the revised algorithm to measure the capacity.


### 3.2. Measuring Tilting Angle with Level Sensors

In Figure4, the sensor installed on the opposite side of each axis receives different pressures in accordance with the turning, and ideally the sum of the sensors on the opposite side must be the same as when the vehicle is not rotating. The tilting angle is computed based on this definition:(6)
S
e
n
s
o
r
1
_
o
u
t
+
S
e
n
s
o
r
2
_
o
u
t
A
D
C
/
θ
r
e
f
=
x
-Axis  Tilting  Value
,
S
e
n
s
o
r
3
_
o
u
t
+
S
e
n
s
o
r
4
_
o
u
t
A
D
C
/
θ
r
e
f
=
y
-Axis  Tilting  Value
0
=
z
-Axis  Tilting  Value
.Figure 4
Sensor location to measure the sensor.In the equation above, Sensor 1–Sensor 6 are digital values obtained by an ADC board used for implementation, andA
D
C
/
θ
r
e
f is the ADC value (constant) per tilting angle obtained from experiments. Sensors 5 and 6 are only used as a means to increase the accuracy because they duplicate y-axis of rotation. Therefore Sensors 5 and 6 are not used in the data analysis.


### 3.3. Measuring Algorithm of Density

Figure5 is the schematic diagram of the attached sensor in the container in order to compute the density. The density is computed with the signals outputted from the sensor attached to the container:(7)
Sensitivity  of  the  sensor
×
Pressure
=
Output  voltage
⟶
Output  Voltage
Sensitivity  of  the  sensor
=
Pressure
pa
,
ρ
×
g
×
h
2
-
h
1
=
P
2
-
P
1
⟶
ρ
=
P
1
-
P
2
g
h
2
-
h
1
,where ρ is liquid density, g is gravitational acceleration, h is height, and P is pressure.Figure 5
Liquid tank and sensor mimetic diagram to measure density.The equation above is a basic equation to compute liquid density, andP
1 and P
2 from the above equation are Sensor 7 and Sensor 4 in Figure 5, respectively [18, 19].


## 3.1. Measuring Algorithm of Capacity

As demonstrated in Figure3, the capacity can be computed in accordance with the container volume and the input capacity of the liquid.Figure 3
Measuring method for capacity.The mass can be computed using the pressure applied to the bottom of a container with liquid. First, the calculation equation of the capacity is as follows:(1)
Pressure  applied  to  the  bottom  of  the  container
∝
liquid  level
×
liquid  density
.The equation to compute the force applied to the bottom in accordance with the liquid is as follows [17]:(2)
F
=
ρ
×
A
h
×
g
g
:  Gravitational  acceleration
.Therefore the pressureP applied to the unit area of the bottom side can be defined as in the equation below [2]:(3)
P
=
ρ
×
g
×
h
If  there  is  no  internal  pressure  in  the  tank
,
P
-
P
0
=
ρ
×
g
×
h
If  there  is  an  internal  pressure  in  the  tank
.In this control system, the following method has been used to measure the capacity.The sensor attached toward thex-axis and y-axis from the center of the target container is located on the opposite side of the container. The obtained measurement value is used to calculate the average value of each axis, and the capacity is computed by comparing the accumulated data from experiments: (4)
A
=
Sensor
1_out
+
Sensor
2_out
2
B
=
Sensor
3_out
+
Sensor
4_out
2
C
=
S
e
n
s
o
r
5
_
o
u
t
+
S
e
n
s
o
r
6
_
o
u
t
2
(5)
A
+
B
+
C
3
=
The  average  value  of  the  output  voltage
⟶
The  measuring  comparison  by  the  input  data
.Here the output value of each sensor does not accumulate, but the precision is low. This sensor output value changes with turning and accelerating and is used in the revised algorithm to measure the capacity.


## 3.2. Measuring Tilting Angle with Level Sensors

In Figure4, the sensor installed on the opposite side of each axis receives different pressures in accordance with the turning, and ideally the sum of the sensors on the opposite side must be the same as when the vehicle is not rotating. The tilting angle is computed based on this definition:(6)
S
e
n
s
o
r
1
_
o
u
t
+
S
e
n
s
o
r
2
_
o
u
t
A
D
C
/
θ
r
e
f
=
x
-Axis  Tilting  Value
,
S
e
n
s
o
r
3
_
o
u
t
+
S
e
n
s
o
r
4
_
o
u
t
A
D
C
/
θ
r
e
f
=
y
-Axis  Tilting  Value
0
=
z
-Axis  Tilting  Value
.Figure 4
Sensor location to measure the sensor.In the equation above, Sensor 1–Sensor 6 are digital values obtained by an ADC board used for implementation, andA
D
C
/
θ
r
e
f is the ADC value (constant) per tilting angle obtained from experiments. Sensors 5 and 6 are only used as a means to increase the accuracy because they duplicate y-axis of rotation. Therefore Sensors 5 and 6 are not used in the data analysis.


## 3.3. Measuring Algorithm of Density

Figure5 is the schematic diagram of the attached sensor in the container in order to compute the density. The density is computed with the signals outputted from the sensor attached to the container:(7)
Sensitivity  of  the  sensor
×
Pressure
=
Output  voltage
⟶
Output  Voltage
Sensitivity  of  the  sensor
=
Pressure
pa
,
ρ
×
g
×
h
2
-
h
1
=
P
2
-
P
1
⟶
ρ
=
P
1
-
P
2
g
h
2
-
h
1
,where ρ is liquid density, g is gravitational acceleration, h is height, and P is pressure.Figure 5
Liquid tank and sensor mimetic diagram to measure density.The equation above is a basic equation to compute liquid density, andP
1 and P
2 from the above equation are Sensor 7 and Sensor 4 in Figure 5, respectively [18, 19].


## 4. Multiple Sensors Algorithm to Improve the Precision of the Result Value of the Sensor

Figure6 is a block diagram of the algorithm to accurately measure the capacity when the liquid tank has rotated by the force applied to it. An attitude Kalman filter computes accurate values, while accumulating almost no errors, by receiving the rotation angle data from the IMU and level sensor. Furthermore, a level Kalman filter precisely measures the capacity by using the obtained rotation angle data and level increment. The holistic drift reduction (HDR) filter used is excellent in avoiding drift in the obtained data from the IMU. The acceleration of gravity is distributed to the gravity value of each axis in accordance with the fuel tank conditions.Figure 6
Measuring methods of the rotation angle of the liquid tank and acting forces.


### 4.1. Revision of IMU Sensor

When the output of the IMU sensor has an angular speed ofω
~
b, ω
b with a revised drift and scale is computed using the HDR algorithm: (8)
ω
b
=
ω
~
b
-
b
g
.Here the HDR algorithm is used to compute the bias valueb
g.


### 4.2. Attitude Kalman Filter

#### 4.2.1. Initialize

The state variables to be used in the Kalman filter are indicated as a rotation matrix and as Eulerian angles, and the covariance matrix is indicated as the covariance of Eulerian angles.State variables are as follows:(9)
x
=
R
b
n
State variables for the rotation matrix
y
=
Λ
n
State variables for the Eulerian angles
.
Covariance matrix of Eulerian angles is as follows:(10)
P
=
σ
∅
∅
σ
∅
θ
σ
∅
ψ
σ
∅
θ
σ
θ
θ
σ
θ
ψ
σ
∅
φ
σ
θ
ψ
σ
ψ
ψ
.
State variables and covariance matrix are initialized as follows:(11)
x
^
k
=
0
=
I
3
×
3
,
y
^
k
=
0
=
0
0
0
T
,
P
k
=
0
=
10
8
I
3
×
3
.


#### 4.2.2. Attitude Prediction

(i) System Model. The angular speed measured by sensors has differentials of the rotation matrix and the following relation:(12)
R
˙
b
n
=
R
b
n
Ω
b
.Here, Ω
b is the skew-symmetric matrix of ω
b and is defined as follows:(13)
Ω
b
=
ω
b
×
0
-
ω
z
ω
y
ω
z
0
-
ω
x
-
ω
y
ω
x
0
.In addition, the differentials of Eulerian angles have angular speed and the following relation:(14)
Λ
˙
n
=
C
1
-
1
ω
b
.Here, C
1
-
1 is a matrix that converts the angular speed measured by the gyro sensor into Eulerian angles:(15)
C
1
-
1
=
1
sin
⁡
∅
tan
⁡
∅
cos
⁡
∅
tan
⁡
θ
0
cos
⁡
∅
-
sin
⁡
∅
0
sin
⁡
∅
cos
⁡
θ
cos
⁡
∅
cos
⁡
θ
.Ifθ
≈
±
π
/
2 in the matrix above, it is cos
⁡
θ
≈
0, which generates the denominator of matrix element to be close to 0 at times. However, such a case refers to the fact that when the sensor stands vertically, this does not occur if one considers the location and direction when attaching the sensor.(ii) Prediction. Consider(16)
X
^
k
-
=
X
^
k
-
1
1
+
Ω
b
Δ
t
=
X
^
k
-
1
R
z
ω
z
Δ
t
R
y
ω
y
Δ
t
R
x
ω
x
Δ
t
,
y
^
k
-
=
y
^
k
-
1
+
C
1
-
1
ω
b
Δ
t
.The covariance of Eulerian angles is updated as follows:(17)
P
k
-
=
P
k
-
1
+
C
1
-
1
Q
C
1
-
T
Δ
t
2
.Here, Q
=
d
i
a
g
(
σ
∅
2
,
σ
θ
2
,
σ
φ
2
) is computed as follows:(18)
σ
A
=
σ
∅
σ
θ
σ
ψ
T
=
10
+
ω
b
.


#### 4.2.3. Attitude Updates by Gravity

(i) Attitude Measurement Model of Sensors of Gravity. The matrix that indicates the attitude of sensors utilizing ∅
~ and θ
~ is revised as follows: (19)
R
b
n
⟵
R
y
θ
~
R
x
∅
~
R
b
n
.Furthermore, Eulerian angles are revised as follows:(20)
Λ
n
⟵
Λ
n
+
C
0
-
1
∅
~
θ
~
0
.Here, C
0
-
1 is as follows:(21)
C
0
-
1
=
cos
⁡
ψ
cos
⁡
θ
sin
⁡
ψ
cos
⁡
θ
0
sin
⁡
ψ
cos
⁡
ψ
0
cos
⁡
ψ
tan
⁡
θ
sin
⁡
ψ
tan
⁡
θ
1
.The attitude measurement model function h
(
·
) of the sensor utilizes Eulerian angles without changes as follows:(22)
h
Λ
n
=
Λ
n
.The attitude measurement value Z
k is computed as follows, utilizing ∅
~ and θ
~:(23)
Z
k
=
Λ
n
+
C
0
-
1
∅
~
θ
~
0
.
(ii) Attitude Update.One has(24)
x
^
k
=
R
y
K
θ
Δ
θ
R
x
K
∅
Δ
∅
x
^
k
-
,
K
^
k
=
k
^
k
-
+
K
k
z
k
-
h
y
^
k
-
.Here, K
∅ is column 1 row 1 element of K
k, and K
θ is column 2 row 2 elements of K
k:(25)
P
k
=
P
k
-
=
K
k
P
k
-
,
K
k
=
P
k
-
P
k
-
+
C
0
-
1
R
C
0
-
T
-
1
.Here, R
=
diag
⁡
(
σ
∅
2
,
σ
θ
2
,
σ
φ
2
) is computed as below because the benefit must increase when acceleration is near 1 g, and it should decrease when the size of acceleration exceeds 1 g:(26)
σ
∅
=
σ
θ
=
0.1
+
f
b
-
1
+
100
ω
b
,
σ
ψ
=
10
8
.


#### 4.2.4. Attitude Update by Level Sensors (DIFF)

The fuel tank of a vehicle has its direction of motion set to thex-axis (heading), so the attitude may be updated with the rotation angle data obtained by level sensors. First, velocity and rotation angle DIFF obtained from the IMU are brought in. Then, from the rotation matrix R
b
n, the yaw angle of the vehicle can be computed as follows:(27)
Ψ
υ
=
tan
-
1
⁡
a
21
a
11
.Here, a
i
j is the matrix “I” of the rotation matrix R
b
n and the element of column “j.”The matrix indicating the attitude of the fuel tank from the angle and direction of vehicleΨ
υ measured by level sensors is revised as follows:(28)
R
b
n
⟵
R
z
Ψ
~
R
b
n
,
Ψ
~
≅
Ψ
level
-
Ψ
.Furthermore, Eulerian angles are revised as follows:(29)
Λ
n
=
Λ
n
+
C
0
-
1
0
0
Ψ
~
.The attitude measurement model function h
(
∗
) of a fuel tank still utilizes the Eulerian angles as follows: (30)
h
Λ
n
=
Λ
n
.The attitude measurement value Z
k is computed as follows, utilizing Ψ:(31)
Z
k
=
Λ
n
+
C
0
-
1
0
0
Ψ
~
.
(i) Attitude Update. One has(32)
x
^
k
=
R
z
K
Ψ
Ψ
~
x
^
k
-
,
y
^
k
=
y
^
k
-
+
K
k
z
k
-
h
x
^
k
-
.Here, K
Ψ is a 3
×
3 element of K
k:(33)
R
k
=
P
k
-
-
K
k
P
k
-
,
K
k
=
P
k
-
P
k
-
+
C
0
-
1
R
C
0
-
T
-
1
.Here, R
=
diag
⁡
(
σ
∅
2
,
σ
θ
2
,
σ
φ
2
) is computed as follows:(34)
σ
∅
=
σ
θ
=
10
8
,
σ
ψ
=
2
σ
level
ϵ
+
V
IMU
,
ϵ
=
10
-
6
,
σ
l
e
v
e
l
=
10
.


#### 4.2.5. Attitude Measurement of Vehicles by Gravity

Since the acceleration of gravity is always toward the center of the earth,g
=
(
0,0
,
-
9.81
m
/
s
2
) is measured if another force is not applied to the acceleration sensor. The sensor attitude can be revised by comparing the acceleration measured by the acceleration sensor and the acceleration of gravity. However, this is possible if gravity is the only force applied to the acceleration sensor. Since it is not possible to separate the acceleration of gravity from other forces and measure it, the simplest method is to assess if the condition of the applied forces other than the acceleration of gravity is ‖
a
‖
≈
‖
g
‖.a
=
a
x
a
y
a
z, the acceleration measured by the acceleration sensor, includes various types of acceleration in accordance with the acceleration of gravity and the acceleration of the sensor. As an equation, this would be as follows:(35)
a
=
v
˙
+
v
×
ω
+
R
g
T
∵
R
T
=
R
-
1
=
v
˙
x
v
˙
y
v
˙
z
+
0
v
z
-
v
y
-
v
z
0
v
x
v
y
-
v
x
0
w
x
w
y
w
z
+
g
z
-
sin
⁡
θ
cos
⁡
θ
sin
⁡
ϕ
cos
⁡
θ
cos
⁡
ϕ
.


## 4.1. Revision of IMU Sensor

When the output of the IMU sensor has an angular speed ofω
~
b, ω
b with a revised drift and scale is computed using the HDR algorithm: (8)
ω
b
=
ω
~
b
-
b
g
.Here the HDR algorithm is used to compute the bias valueb
g.


## 4.2. Attitude Kalman Filter

### 4.2.1. Initialize

The state variables to be used in the Kalman filter are indicated as a rotation matrix and as Eulerian angles, and the covariance matrix is indicated as the covariance of Eulerian angles.State variables are as follows:(9)
x
=
R
b
n
State variables for the rotation matrix
y
=
Λ
n
State variables for the Eulerian angles
.
Covariance matrix of Eulerian angles is as follows:(10)
P
=
σ
∅
∅
σ
∅
θ
σ
∅
ψ
σ
∅
θ
σ
θ
θ
σ
θ
ψ
σ
∅
φ
σ
θ
ψ
σ
ψ
ψ
.
State variables and covariance matrix are initialized as follows:(11)
x
^
k
=
0
=
I
3
×
3
,
y
^
k
=
0
=
0
0
0
T
,
P
k
=
0
=
10
8
I
3
×
3
.


### 4.2.2. Attitude Prediction

(i) System Model. The angular speed measured by sensors has differentials of the rotation matrix and the following relation:(12)
R
˙
b
n
=
R
b
n
Ω
b
.Here, Ω
b is the skew-symmetric matrix of ω
b and is defined as follows:(13)
Ω
b
=
ω
b
×
0
-
ω
z
ω
y
ω
z
0
-
ω
x
-
ω
y
ω
x
0
.In addition, the differentials of Eulerian angles have angular speed and the following relation:(14)
Λ
˙
n
=
C
1
-
1
ω
b
.Here, C
1
-
1 is a matrix that converts the angular speed measured by the gyro sensor into Eulerian angles:(15)
C
1
-
1
=
1
sin
⁡
∅
tan
⁡
∅
cos
⁡
∅
tan
⁡
θ
0
cos
⁡
∅
-
sin
⁡
∅
0
sin
⁡
∅
cos
⁡
θ
cos
⁡
∅
cos
⁡
θ
.Ifθ
≈
±
π
/
2 in the matrix above, it is cos
⁡
θ
≈
0, which generates the denominator of matrix element to be close to 0 at times. However, such a case refers to the fact that when the sensor stands vertically, this does not occur if one considers the location and direction when attaching the sensor.(ii) Prediction. Consider(16)
X
^
k
-
=
X
^
k
-
1
1
+
Ω
b
Δ
t
=
X
^
k
-
1
R
z
ω
z
Δ
t
R
y
ω
y
Δ
t
R
x
ω
x
Δ
t
,
y
^
k
-
=
y
^
k
-
1
+
C
1
-
1
ω
b
Δ
t
.The covariance of Eulerian angles is updated as follows:(17)
P
k
-
=
P
k
-
1
+
C
1
-
1
Q
C
1
-
T
Δ
t
2
.Here, Q
=
d
i
a
g
(
σ
∅
2
,
σ
θ
2
,
σ
φ
2
) is computed as follows:(18)
σ
A
=
σ
∅
σ
θ
σ
ψ
T
=
10
+
ω
b
.


### 4.2.3. Attitude Updates by Gravity

(i) Attitude Measurement Model of Sensors of Gravity. The matrix that indicates the attitude of sensors utilizing ∅
~ and θ
~ is revised as follows: (19)
R
b
n
⟵
R
y
θ
~
R
x
∅
~
R
b
n
.Furthermore, Eulerian angles are revised as follows:(20)
Λ
n
⟵
Λ
n
+
C
0
-
1
∅
~
θ
~
0
.Here, C
0
-
1 is as follows:(21)
C
0
-
1
=
cos
⁡
ψ
cos
⁡
θ
sin
⁡
ψ
cos
⁡
θ
0
sin
⁡
ψ
cos
⁡
ψ
0
cos
⁡
ψ
tan
⁡
θ
sin
⁡
ψ
tan
⁡
θ
1
.The attitude measurement model function h
(
·
) of the sensor utilizes Eulerian angles without changes as follows:(22)
h
Λ
n
=
Λ
n
.The attitude measurement value Z
k is computed as follows, utilizing ∅
~ and θ
~:(23)
Z
k
=
Λ
n
+
C
0
-
1
∅
~
θ
~
0
.
(ii) Attitude Update.One has(24)
x
^
k
=
R
y
K
θ
Δ
θ
R
x
K
∅
Δ
∅
x
^
k
-
,
K
^
k
=
k
^
k
-
+
K
k
z
k
-
h
y
^
k
-
.Here, K
∅ is column 1 row 1 element of K
k, and K
θ is column 2 row 2 elements of K
k:(25)
P
k
=
P
k
-
=
K
k
P
k
-
,
K
k
=
P
k
-
P
k
-
+
C
0
-
1
R
C
0
-
T
-
1
.Here, R
=
diag
⁡
(
σ
∅
2
,
σ
θ
2
,
σ
φ
2
) is computed as below because the benefit must increase when acceleration is near 1 g, and it should decrease when the size of acceleration exceeds 1 g:(26)
σ
∅
=
σ
θ
=
0.1
+
f
b
-
1
+
100
ω
b
,
σ
ψ
=
10
8
.


### 4.2.4. Attitude Update by Level Sensors (DIFF)

The fuel tank of a vehicle has its direction of motion set to thex-axis (heading), so the attitude may be updated with the rotation angle data obtained by level sensors. First, velocity and rotation angle DIFF obtained from the IMU are brought in. Then, from the rotation matrix R
b
n, the yaw angle of the vehicle can be computed as follows:(27)
Ψ
υ
=
tan
-
1
⁡
a
21
a
11
.Here, a
i
j is the matrix “I” of the rotation matrix R
b
n and the element of column “j.”The matrix indicating the attitude of the fuel tank from the angle and direction of vehicleΨ
υ measured by level sensors is revised as follows:(28)
R
b
n
⟵
R
z
Ψ
~
R
b
n
,
Ψ
~
≅
Ψ
level
-
Ψ
.Furthermore, Eulerian angles are revised as follows:(29)
Λ
n
=
Λ
n
+
C
0
-
1
0
0
Ψ
~
.The attitude measurement model function h
(
∗
) of a fuel tank still utilizes the Eulerian angles as follows: (30)
h
Λ
n
=
Λ
n
.The attitude measurement value Z
k is computed as follows, utilizing Ψ:(31)
Z
k
=
Λ
n
+
C
0
-
1
0
0
Ψ
~
.
(i) Attitude Update. One has(32)
x
^
k
=
R
z
K
Ψ
Ψ
~
x
^
k
-
,
y
^
k
=
y
^
k
-
+
K
k
z
k
-
h
x
^
k
-
.Here, K
Ψ is a 3
×
3 element of K
k:(33)
R
k
=
P
k
-
-
K
k
P
k
-
,
K
k
=
P
k
-
P
k
-
+
C
0
-
1
R
C
0
-
T
-
1
.Here, R
=
diag
⁡
(
σ
∅
2
,
σ
θ
2
,
σ
φ
2
) is computed as follows:(34)
σ
∅
=
σ
θ
=
10
8
,
σ
ψ
=
2
σ
level
ϵ
+
V
IMU
,
ϵ
=
10
-
6
,
σ
l
e
v
e
l
=
10
.


### 4.2.5. Attitude Measurement of Vehicles by Gravity

Since the acceleration of gravity is always toward the center of the earth,g
=
(
0,0
,
-
9.81
m
/
s
2
) is measured if another force is not applied to the acceleration sensor. The sensor attitude can be revised by comparing the acceleration measured by the acceleration sensor and the acceleration of gravity. However, this is possible if gravity is the only force applied to the acceleration sensor. Since it is not possible to separate the acceleration of gravity from other forces and measure it, the simplest method is to assess if the condition of the applied forces other than the acceleration of gravity is ‖
a
‖
≈
‖
g
‖.a
=
a
x
a
y
a
z, the acceleration measured by the acceleration sensor, includes various types of acceleration in accordance with the acceleration of gravity and the acceleration of the sensor. As an equation, this would be as follows:(35)
a
=
v
˙
+
v
×
ω
+
R
g
T
∵
R
T
=
R
-
1
=
v
˙
x
v
˙
y
v
˙
z
+
0
v
z
-
v
y
-
v
z
0
v
x
v
y
-
v
x
0
w
x
w
y
w
z
+
g
z
-
sin
⁡
θ
cos
⁡
θ
sin
⁡
ϕ
cos
⁡
θ
cos
⁡
ϕ
.


## 4.2.1. Initialize

The state variables to be used in the Kalman filter are indicated as a rotation matrix and as Eulerian angles, and the covariance matrix is indicated as the covariance of Eulerian angles.State variables are as follows:(9)
x
=
R
b
n
State variables for the rotation matrix
y
=
Λ
n
State variables for the Eulerian angles
.
Covariance matrix of Eulerian angles is as follows:(10)
P
=
σ
∅
∅
σ
∅
θ
σ
∅
ψ
σ
∅
θ
σ
θ
θ
σ
θ
ψ
σ
∅
φ
σ
θ
ψ
σ
ψ
ψ
.
State variables and covariance matrix are initialized as follows:(11)
x
^
k
=
0
=
I
3
×
3
,
y
^
k
=
0
=
0
0
0
T
,
P
k
=
0
=
10
8
I
3
×
3
.


## 4.2.2. Attitude Prediction

(i) System Model. The angular speed measured by sensors has differentials of the rotation matrix and the following relation:(12)
R
˙
b
n
=
R
b
n
Ω
b
.Here, Ω
b is the skew-symmetric matrix of ω
b and is defined as follows:(13)
Ω
b
=
ω
b
×
0
-
ω
z
ω
y
ω
z
0
-
ω
x
-
ω
y
ω
x
0
.In addition, the differentials of Eulerian angles have angular speed and the following relation:(14)
Λ
˙
n
=
C
1
-
1
ω
b
.Here, C
1
-
1 is a matrix that converts the angular speed measured by the gyro sensor into Eulerian angles:(15)
C
1
-
1
=
1
sin
⁡
∅
tan
⁡
∅
cos
⁡
∅
tan
⁡
θ
0
cos
⁡
∅
-
sin
⁡
∅
0
sin
⁡
∅
cos
⁡
θ
cos
⁡
∅
cos
⁡
θ
.Ifθ
≈
±
π
/
2 in the matrix above, it is cos
⁡
θ
≈
0, which generates the denominator of matrix element to be close to 0 at times. However, such a case refers to the fact that when the sensor stands vertically, this does not occur if one considers the location and direction when attaching the sensor.(ii) Prediction. Consider(16)
X
^
k
-
=
X
^
k
-
1
1
+
Ω
b
Δ
t
=
X
^
k
-
1
R
z
ω
z
Δ
t
R
y
ω
y
Δ
t
R
x
ω
x
Δ
t
,
y
^
k
-
=
y
^
k
-
1
+
C
1
-
1
ω
b
Δ
t
.The covariance of Eulerian angles is updated as follows:(17)
P
k
-
=
P
k
-
1
+
C
1
-
1
Q
C
1
-
T
Δ
t
2
.Here, Q
=
d
i
a
g
(
σ
∅
2
,
σ
θ
2
,
σ
φ
2
) is computed as follows:(18)
σ
A
=
σ
∅
σ
θ
σ
ψ
T
=
10
+
ω
b
.


## 4.2.3. Attitude Updates by Gravity

(i) Attitude Measurement Model of Sensors of Gravity. The matrix that indicates the attitude of sensors utilizing ∅
~ and θ
~ is revised as follows: (19)
R
b
n
⟵
R
y
θ
~
R
x
∅
~
R
b
n
.Furthermore, Eulerian angles are revised as follows:(20)
Λ
n
⟵
Λ
n
+
C
0
-
1
∅
~
θ
~
0
.Here, C
0
-
1 is as follows:(21)
C
0
-
1
=
cos
⁡
ψ
cos
⁡
θ
sin
⁡
ψ
cos
⁡
θ
0
sin
⁡
ψ
cos
⁡
ψ
0
cos
⁡
ψ
tan
⁡
θ
sin
⁡
ψ
tan
⁡
θ
1
.The attitude measurement model function h
(
·
) of the sensor utilizes Eulerian angles without changes as follows:(22)
h
Λ
n
=
Λ
n
.The attitude measurement value Z
k is computed as follows, utilizing ∅
~ and θ
~:(23)
Z
k
=
Λ
n
+
C
0
-
1
∅
~
θ
~
0
.
(ii) Attitude Update.One has(24)
x
^
k
=
R
y
K
θ
Δ
θ
R
x
K
∅
Δ
∅
x
^
k
-
,
K
^
k
=
k
^
k
-
+
K
k
z
k
-
h
y
^
k
-
.Here, K
∅ is column 1 row 1 element of K
k, and K
θ is column 2 row 2 elements of K
k:(25)
P
k
=
P
k
-
=
K
k
P
k
-
,
K
k
=
P
k
-
P
k
-
+
C
0
-
1
R
C
0
-
T
-
1
.Here, R
=
diag
⁡
(
σ
∅
2
,
σ
θ
2
,
σ
φ
2
) is computed as below because the benefit must increase when acceleration is near 1 g, and it should decrease when the size of acceleration exceeds 1 g:(26)
σ
∅
=
σ
θ
=
0.1
+
f
b
-
1
+
100
ω
b
,
σ
ψ
=
10
8
.


## 4.2.4. Attitude Update by Level Sensors (DIFF)

The fuel tank of a vehicle has its direction of motion set to thex-axis (heading), so the attitude may be updated with the rotation angle data obtained by level sensors. First, velocity and rotation angle DIFF obtained from the IMU are brought in. Then, from the rotation matrix R
b
n, the yaw angle of the vehicle can be computed as follows:(27)
Ψ
υ
=
tan
-
1
⁡
a
21
a
11
.Here, a
i
j is the matrix “I” of the rotation matrix R
b
n and the element of column “j.”The matrix indicating the attitude of the fuel tank from the angle and direction of vehicleΨ
υ measured by level sensors is revised as follows:(28)
R
b
n
⟵
R
z
Ψ
~
R
b
n
,
Ψ
~
≅
Ψ
level
-
Ψ
.Furthermore, Eulerian angles are revised as follows:(29)
Λ
n
=
Λ
n
+
C
0
-
1
0
0
Ψ
~
.The attitude measurement model function h
(
∗
) of a fuel tank still utilizes the Eulerian angles as follows: (30)
h
Λ
n
=
Λ
n
.The attitude measurement value Z
k is computed as follows, utilizing Ψ:(31)
Z
k
=
Λ
n
+
C
0
-
1
0
0
Ψ
~
.
(i) Attitude Update. One has(32)
x
^
k
=
R
z
K
Ψ
Ψ
~
x
^
k
-
,
y
^
k
=
y
^
k
-
+
K
k
z
k
-
h
x
^
k
-
.Here, K
Ψ is a 3
×
3 element of K
k:(33)
R
k
=
P
k
-
-
K
k
P
k
-
,
K
k
=
P
k
-
P
k
-
+
C
0
-
1
R
C
0
-
T
-
1
.Here, R
=
diag
⁡
(
σ
∅
2
,
σ
θ
2
,
σ
φ
2
) is computed as follows:(34)
σ
∅
=
σ
θ
=
10
8
,
σ
ψ
=
2
σ
level
ϵ
+
V
IMU
,
ϵ
=
10
-
6
,
σ
l
e
v
e
l
=
10
.


## 4.2.5. Attitude Measurement of Vehicles by Gravity

Since the acceleration of gravity is always toward the center of the earth,g
=
(
0,0
,
-
9.81
m
/
s
2
) is measured if another force is not applied to the acceleration sensor. The sensor attitude can be revised by comparing the acceleration measured by the acceleration sensor and the acceleration of gravity. However, this is possible if gravity is the only force applied to the acceleration sensor. Since it is not possible to separate the acceleration of gravity from other forces and measure it, the simplest method is to assess if the condition of the applied forces other than the acceleration of gravity is ‖
a
‖
≈
‖
g
‖.a
=
a
x
a
y
a
z, the acceleration measured by the acceleration sensor, includes various types of acceleration in accordance with the acceleration of gravity and the acceleration of the sensor. As an equation, this would be as follows:(35)
a
=
v
˙
+
v
×
ω
+
R
g
T
∵
R
T
=
R
-
1
=
v
˙
x
v
˙
y
v
˙
z
+
0
v
z
-
v
y
-
v
z
0
v
x
v
y
-
v
x
0
w
x
w
y
w
z
+
g
z
-
sin
⁡
θ
cos
⁡
θ
sin
⁡
ϕ
cos
⁡
θ
cos
⁡
ϕ
.


## 5. Capacity Measuring Algorithm through Revised Rotation Angle Data and Measurement Value of Level Pressure

Capacity is computed as a digital value through the ADC board. When the acquired ADC and capacity are =S
e
n
s
o
r
a
d
c
/
(
l
/
a
d
c
), the equation to compute the capacity is as follows.


### 5.1. Level Kalman Filter

#### 5.1.1. System Model

The difference value of the level sensor measured in the coordinate system of the fuel tank is converted into a navigation system:(36)
level
˙
n
=
R
b
n
D
i
f
f
level
b
.The formula above is made into an equation of state as follows:(37)
level
n
level
˙
n
=
1
R
b
n
0
1
level
n
level
˙
n
.State variables are(38)
x
^
k
=
level
n
level
˙
n
.


#### 5.1.2. Prediction

In the formula below,G has a bias value due to noise and the level. F has 0
0
I
0. The equation depends on the general Kalman equation: (39)
x
^
k
-
=
A
x
^
k
-
1
+
G
u
k
Δ
t
,
P
k
-
=
A
P
k
-
1
A
T
+
G
Q
G
T
Δ
t
2
,
A
=
I
+
F
Δ
t
,
Q
=
diag
⁡
σ
f
x
2
,
σ
f
y
2
,
σ
f
z
2
,
σ
g
x
2
,
σ
g
y
2
,
σ
g
z
2
.


## 5.1. Level Kalman Filter

### 5.1.1. System Model

The difference value of the level sensor measured in the coordinate system of the fuel tank is converted into a navigation system:(36)
level
˙
n
=
R
b
n
D
i
f
f
level
b
.The formula above is made into an equation of state as follows:(37)
level
n
level
˙
n
=
1
R
b
n
0
1
level
n
level
˙
n
.State variables are(38)
x
^
k
=
level
n
level
˙
n
.


### 5.1.2. Prediction

In the formula below,G has a bias value due to noise and the level. F has 0
0
I
0. The equation depends on the general Kalman equation: (39)
x
^
k
-
=
A
x
^
k
-
1
+
G
u
k
Δ
t
,
P
k
-
=
A
P
k
-
1
A
T
+
G
Q
G
T
Δ
t
2
,
A
=
I
+
F
Δ
t
,
Q
=
diag
⁡
σ
f
x
2
,
σ
f
y
2
,
σ
f
z
2
,
σ
g
x
2
,
σ
g
y
2
,
σ
g
z
2
.


## 5.1.1. System Model

The difference value of the level sensor measured in the coordinate system of the fuel tank is converted into a navigation system:(36)
level
˙
n
=
R
b
n
D
i
f
f
level
b
.The formula above is made into an equation of state as follows:(37)
level
n
level
˙
n
=
1
R
b
n
0
1
level
n
level
˙
n
.State variables are(38)
x
^
k
=
level
n
level
˙
n
.


## 5.1.2. Prediction

In the formula below,G has a bias value due to noise and the level. F has 0
0
I
0. The equation depends on the general Kalman equation: (39)
x
^
k
-
=
A
x
^
k
-
1
+
G
u
k
Δ
t
,
P
k
-
=
A
P
k
-
1
A
T
+
G
Q
G
T
Δ
t
2
,
A
=
I
+
F
Δ
t
,
Q
=
diag
⁡
σ
f
x
2
,
σ
f
y
2
,
σ
f
z
2
,
σ
g
x
2
,
σ
g
y
2
,
σ
g
z
2
.


## 6. Experiment

Table2 shows the specifications of the pressure sensor applied to the system.Table 2
26 PC Series performance characteristics [7].
Min.
Type
Max.
Units
Excitation
—
10
16
Vdc
Response time
—
—
1.0
ms
Input resistance
5.5 k
7.5 k
11.5 k
Ohm
Output resistance
1.5 k
2.5 k
3.0 k
Ohm
Weight
2
GramThe display device of the algorithm built on this research has a level sensor (Figure7) to measure absolute pressure and the tilting angle; an IMU (Figure 8) for measuring the rotation angle, acceleration, and so forth; a programmable AMP (Figure 9) that can consistently program the increase rate of pressure between sensors by amplifying the values of pressure sensors in the desired offset; and a test bed (Figure 10).Figure 7
Level sensor (side view of 26 PCSPS pressure sensors used).Figure 8
IMU: analog device ADIS16365.Figure 9
AMP: analog device AD8555.Figure 10
Test bed: internally developed.Capacity Measurement of Distilled Water Using Level Sensors. In order to measure capacity using level sensors, the sensor output must have a linear output in accordance with the capacity. To assess it, the research team sequentially inserted 1 L–11 L of water (distilled water, density ≈ 1) into the test bed and assessed each sensor value. In order to avoid continuity and obtain objective data, the experiment was conducted once a day over three days. After the experiment, the team removed all distilled water from the test bed and conducted the experiment again.Figure11 demonstrates a 1 L/75 ADC_value and 0.25 V/75 ADC_value. Figure 11 indicates the difference of values in the first, second, and third measurements and each sensor value with a margin of error of ±1%. The output graph of Sensors numbers 5 and 6 differs from the output graph of Sensors numbers 1, 2, 3, and 4 owing to the installed location. Sensor number 7 reacts when it is over 11 L because it was located on the upper part of the test bed to measure density.Figure 11
Output value of sensor in accordance with water quantity (1–11 L).The level sensors used in this research are linear overall, demonstrating that capacity and density computation is possible. However, the margin of error increases with rotation. Figure12 shows the output value measured three times at a 0
°, 15
°, 30
°, and 45
° rotation, and it indicates that the value increased or decreased significantly in accordance with the rotation angle.Figure 12
Output graph of sensor in accordance with rotation. (a) Tilt toward Sensor 3; (b) tilt toward Sensor 4.
(a)
(b)When tilting toward Sensor 3, the value of Sensor 3 and Sensor 5 on the Sensor 3-axis increases, whereas the value of Sensor 4 and Sensor 6 on the opposite side decreases. When tilting toward Sensor 4, the value of Sensor 4 and Sensor 6 on the Sensor 4-axis increases, whereas the value of Sensor 0 and Sensor 2 on the opposite side decreases. Sensor 1 and Sensor 2, which are irrelevant to the tilting direction, demonstrated no change, and so the team has not included it here.The experiment results indicate that measuring capacity and density only with a level sensor generates a large margin of error when there is tilting involved. Therefore, it is inadequate to measure the capacity and density of a liquid using only level sensor.When revising the rotation angle and acceleration of gravity with the algorithm this paper proposes, the measurement results are obtained in Figure13.Figure 13
Measurement result graph using algorithm that has revised rotation angle and acceleration of gravity.The research team inserted 11 L of distilled water into the test bed, tilted it by5
°, and measured the sensor value. Figure 13 demonstrates a stable result compared with the sensor graph before the application of the algorithm. This is because the accurate estimation of the rotation angle and its revision value were designed with a Kalman filter. Thus estimating the rotation angle is important in this algorithm.The comparative results of the proposed algorithm that estimates the rotation angle are as shown in Figure14.Figure 14
Result graph of measuring the rotation angle using IMU data.Figure14 demonstrates the measurements of rotation angles using only the IMU; the measurement is carried out once every minute. As indicated in Figure 14, the margin of error continuously accumulated, demonstrating a large difference in value after 17 min from the beginning. Furthermore, when measuring the rotation angle with only level sensors (Figure 15), the margin of error does not accumulate and is drifting as with the IMU data, but the margin of error is significantly large.Figure 15
Result graph of measuring the rotation angle using level sensor data.Figure16 shows the rotation angle using the algorithm proposed by this paper. The margin of error is within ±5%, which is relatively accurate compared with existing methods.Figure 16
Result graph of measuring the rotation angle with the proposed algorithm.The results of measuring capacity and density using the algorithm this paper has proposed are as  shown in Figure17.Figure 17
Measurement graph of capacity and density using the proposed algorithm.Figure17 shows the output values of the experiment, where 1 L was inserted as a sensor value without tilting. However, the measurements were not conducted sequentially because each sensor has a slightly different minimum pressure value that computes the sensor output and because the sensor locations differ as well. For precise measurements, the research team inserted 3 L of distilled water, and when all sensors were reacting sufficiently, the team initialized the value and measured again. As demonstrated in Figure 17, accurate results from the measuring capacity could be obtained because the sensor value was close to the primary equation and each sensor had less difference in value.Figure18 is the computation of density and capacity, which is very precise because it has a margin of error of only ±0.05 L compared with the actual capacity. The density was also very precise, with a margin of error of ±0.02 kg/m3.Figure 18
Comparative graph 1 to assess the precision of capacity and density measurements. (a) Accuracy graph of suggested algorithm for capacity; (b) accuracy graph of suggested algorithm for density.
(a)
(b)Sensor number 4 and Sensor number 7 are used to measure the density. Sensor number 7 is placed on the upper part so it reacts only when the solution is more than 10 liters. Therefore the density can be measured from 10 liters. But if the sensor can be placed lower, we can measure the density by using less fluid.In addition, Figure19 demonstrates the results of the measurements with the same method as above, but tilting by 5
°.Figure 19
Comparative graph 2 to assess the precision of measurements of capacity and density.As indicated in Figure19, the margin of error of measurement increases in accordance with the tilting, and the maximum margin of error was only 0.5 L. Since it can be presumed that the larger the target capacity for measurement, the lesser the error ratio, this can be said to be a precise algorithm.As demonstrated in Figure20, compared with the results when there was no tilting, there has rarely been any difference in density, and the result was ±0.05 kg/m3.Figure 20
Comparative graph 3 to assess the precision of measurements of capacity and density.For Figure21, the program was implemented with a graphical user interface (GUI), and its resulting values were continuously obtained and observed.Figure 21
GUI to observe the measurement values in real time.


## 7. Conclusion

This research proposed an algorithm using a level sensor and IMU sensor to improve the limitations of the existing methods of measuring the liquid capacity. For this, the research team installed level sensors and IMU sensors in a liquid tank and conducted sensing in real time. The proposed algorithm utilized the strengths of each sensor and supplemented their weaknesses, which were assessed through experiments.A measurement evaluation of the liquid level sensor was conducted by increasing and decreasing the liquid in the tank. The team increased the water level by units of 1 L and analyzed the output characteristics of up to 12 L. As a result, the system demonstrated a highly linear precision that was within ±0.5%. The offset voltage of the liquid level sensor was 93.5 mV, with a sensitivity of 363.5 mV/L. The output characteristic in inclined environments resulted in a value of 1% FS. In addition, the sensitivity to density measurements was 1085.28 [kg/m3], which is consistent with the computed value of the sensor output. The reproducibility of the sensor showed minute hysteresis, but since there was repeatability, a high-precision pressure level sensor was implemented as it can ignore this repeatability.The research team believes that, with the proposed algorithm, precise liquid capacity measurements in real time would be possible even during movement or tilting of the liquid tank.



---

*Source: 4260397-2016-07-21.xml*