# Nonlinear fluid damping of elastically mounted pitching wings in quiescent
water

Yuanhang Zhu 1<EMAIL_ADDRESS>Varghese Mathai 2  Kenneth Breuer 1  1
Center for Fluid Mechanics, School of Engineering, Brown University,
Providence, RI 02912, USA 2 Department of Physics, University of
Massachusetts, Amherst, MA 01003, USA

###### Abstract

We experimentally study the nonlinear fluid damping of a rigid but elastically
mounted pitching wing in the absence of a freestream flow. The dynamics of the
elastic mount are simulated using a cyber-physical system. We perturb the wing
and measure the fluid damping coefficient from damped oscillations over a
large range of pitching frequencies, pitching amplitudes, pivot locations and
sweep angles. A universal fluid damping scaling is proposed to incorporate all
these parameters. Flow fields obtained using particle image velocimetry are
analyzed to explain the nonlinear behaviors of the fluid damping.

##  1 Introduction

The interaction between elastically mounted pitching wings and unsteady flows
is central to many applications. With a free-stream flow, this interaction can
lead to self-sustained, flow-induced oscillations, which have been studied for
understanding classic aeroelastic behaviour (Dowell et al., 1989; Dugundji,
2008), as well as in developing oscillating foil energy harvesting devices
(Xiao & Zhu, 2014; Young et al., 2014). Without a free stream, but with
prescribed heaving or flapping (i.e. hovering), the passive flow-induced
pitching motions are used in modelling the thrust generation and maneuvering
in animal flight (Wang, 2005; Bergou et al., 2007; Shinde & Arakeri, 2013;
Kang & Shyy, 2014; Beatus & Cohen, 2015).

One of the critical parameters that govern the flow-structure interactions of
passively pitching wings is the fluid damping. According to the semi-empirical
Morison equation (Morison et al., 1950), the total fluid force exerted on a
wing submerged in unsteady viscous fluid can be divided into two parts – the
force associated with fluid inertia (i.e. the added mass force), which is in
phase with acceleration (Brennen, 1982; Corkery et al., 2019), and the force
induced by vortices in the flow (i.e. the fluid damping force), which is in
phase with velocity (Shih & Buchanan, 1971; Kang & Shyy, 2014; Su & Breuer,
2019). While the structural damping force is typically proportional to
velocity because of the constant structural damping coefficient, the fluid
damping force is expected to scale quadratically with velocity (Morison et
al., 1950; Keulegan & Carpenter, 1958), and due to this nonlinearity, the
fluid damping coefficient is usually obtained empirically as a function of the
reduced frequency, the Reynolds number, the oscillation amplitude, etc (Shih &
Buchanan, 1971). For pitching flexible wings (Alben, 2008) and heaving
membrane wings (Tzezana & Breuer, 2019), the fluid damping coefficient is
found to scale inversely with the oscillation frequency.

For elastically mounted pitching wings with a free stream, the interplay
between the fluid damping and the structural damping governs the flow-induced
oscillation. By mapping out the cycle-averaged energy transfer between the
elastic system and the ambient fluid using prescribed kinematics, Menon &
Mittal (2019) and Zhu et al. (2020) showed that the energy injected by the
negative fluid damping must be equal to the energy dissipated by the positive
structural damping in order for the flow-induced oscillations to sustain. In
other words, the total damping of the system must be zero (Dugundji, 2008).
The negative fluid damping arises primarily from the formation and shedding of
dynamic stall vortices (McCroskey, 1982; Corke & Thomas, 2015). In the absence
of a free stream, however, the fluid damping becomes positive and counteracts
the pitching motion because of the drag effect. With both the fluid damping
and the structural damping being positive, any perturbations to the system
will be damped out. However, little is known about how the fluid damping
shapes the damped oscillations, and understanding this is of critical
importance for understanding the fluid-structure interactions of elastically
mounted pitching wings under external perturbations such as gusts.

In the present study, we use laboratory experiments to characterise the fluid
damping of elastically mounted pitching wings in quiescent water, with the
elastic mount simulated using a cyber-physical system (§2). We perform ‘ring
down’ experiments to extract the fluid damping (§3.1). The effects of many
parameters are investigated, including the effects of the pitching frequency,
the pitching amplitude, the pivot location and the sweep angle (§3.2). We
propose a universal fluid damping scaling to incorporate these parameters
(§3.3), and correlate the nonlinear behaviour of the fluid damping with the
dynamics of the vortical structures measured using particle image velocimetry
(§3.4). Finally, the key findings are summarised in §4.

##  2 Experimental set-up

Figure 1(_a_) shows a schematic of the experimental set-up. We conduct all the
experiments in the Brown University free-surface water tunnel (test section
$\mathrm{width}\times\mathrm{depth}\times\mathrm{length}=0.8~{}\rm{m}\times
0.6~{}\rm{m}\times 4.0~{}\rm{m}$), with the flow speed kept at zero
($U_{\infty}=0$ m/s). A NACA 0012 wing, made of clear acrylic, is mounted
vertically in the tunnel, with an endplate on the top to skim surface waves
and eliminate wingtip vortices at the root. The wing is connected to a six-
axis force/torque transducer (ATI 9105-TIF-Delta-IP65), which measures the
fluid torque $\tau_{f}$ exerted on the wing. This $\tau_{f}$ is then fed into
the cyber-physical system (CPS). Depending on the input virtual structural
parameters, specifically the torsional stiffness $k_{v}$, damping $b_{v}$ and
inertia $I_{v}$, the CPS calculates the pitching position of the wing and
outputs the signal to the servo motor (Parker SM233AE). An optical encoder (US
Digital E3-2500) which is independent of the CPS is used to measure the
pitching position $\theta$. The CPS is operated at 4000 Hz to minimise any
phase delay between the input $\tau_{f}$ and the output $\theta$. A detailed
explanation of the CPS can be found in Zhu et al. (2020).

We use two-dimensional particle image velocimetry (PIV) to measure the flow
field around the wing. The flow is seeded using 50 $\mathrm{\mu m}$ diameter
hollow ceramic spheres and illuminated by a laser sheet at the mid-span plane.
The laser sheet is generated by a double-pulse Nd:YAG laser (532 nm, Quantel
EverGreen) with LaVision sheet optics. The transparent wing enables flow field
measurements on both sides of the wing. Due to the limitation of space beneath
the tunnel, a $45^{\circ}$ mirror is used to reflect the images into two co-
planar sCMOS cameras (LaVision). We use the DaVis software (LaVision) to
calculate (two passes at $64\times 64$ pixels, two passes at $32\times 32$
pixels, both with 50% overlap) and stitch the velocity fields from the two
cameras to form a field of view of $3.2c\times 3.2c$, where $c$ is the chord
length of the wing.

Figure 1: (_a_) A schematic of the experimental set-up. The structural
dynamics of the wing is simulated by a cyber-physical system (CPS). (_b_)
Sketches of unswept and swept wings. The leading edge (LE) and the trailing
edge (TE) are parallel. Dashed lines represent the pivot axis.

Figure 1(_b_) sketches the two types of wings we use in the present study. For
the unswept wing, a wing holder mechanism (not shown) enables the pivot axis
to be adjusted between $x/c=0$ and 1 with a step size of 0.125. For the swept
wings, the sweep angle $\Lambda$ is defined as the angle between the leading
edge and the vertical axis. Four swept wings with $\Lambda=10^{\circ}$,
$15^{\circ}$, $20^{\circ}$ and $25^{\circ}$ are used. As shown in the figure,
the pivot axis of swept wings is a vertical line passing through the mid-chord
point ($x/c=0.5$) of the mid-span plane. All the wings have a span of $s=0.3$
m and a chord length of $c=0.1$ m, which results in an aspect ratio of $AR=3$.

The governing equation of the system is

$I\ddot{\theta}+b\dot{\theta}+k\theta=\tau_{f},$ (1)

where $\theta$, $\dot{\theta}$, and $\ddot{\theta}$ are the angular position,
velocity and acceleration, respectively. $I$, $b$ and $k$ are the effective
inertia, damping and stiffness of the system. The effective inertia $I$ is the
sum of the virtual inertia $I_{v}$, which we prescribe with the CPS, and the
physical inertia $I_{p}$ of the wing (i.e. $I=I_{v}+I_{p}$). The effective
damping $b$ equals the virtual damping $b_{v}$ (i.e. $b=b_{v}$) because the
friction in the system is negligible. The effective stiffness $k$ equals the
virtual stiffness (i.e. $k=k_{v}$). $\tau_{f}$ is the nonlinear fluid torque
experienced by the wing, which can be divided into the added mass torque,
$\tau_{a}=-I_{a}\ddot{\theta}$, where $I_{a}$ is the added fluid inertia, and
the fluid damping torque, for simplicity $\tau_{b}=-b_{f}\dot{\theta}$, where
$b_{f}$ is the fluid damping coefficient (see §1). Note that $b_{f}$ is
expected to be a function of $\dot{\theta}$ (Mathai et al., 2019). Equation 1
can thus be rearranged as

$(I+I_{a})\ddot{\theta}+(b+b_{f})\dot{\theta}+k\theta=0.$ (2)

After a perturbation of amplitude $A_{0}$ is applied at time $t_{0}$, the
damped oscillations of the system can be described as

$\theta=A_{0}e^{-\gamma(t-t_{0})}\cos{[2\pi f_{p}(t-t_{0})]},$ (3)

where

$\gamma=\frac{b+b_{f}}{2(I+I_{a})}~{}~{}\mathrm{and}~{}~{}f_{p}=\frac{1}{2\pi}\sqrt{\frac{k}{I+I_{a}}-\gamma^{2}}.$
$None$

##  3 Results and Discussion

###  3.1 Extracting the fluid damping from ‘ring down’ experiments

Figure 2: (_a_) System response and amplitude decay in a typical ‘ring down’
test, where an elastically mounted unswept wing ($\Lambda=0^{\circ}$) pivots
around the mid-chord ($x/c=0.5$) at a frequency of $f_{p}=0.40$ Hz. The inset
shows the measurements of the pitching amplitude $A_{n}$ and the pitching
frequency $f_{p}$ of the $n$-th peak. The fluid damping $b_{f}$ at $A_{n}$ is
extracted by fitting an exponential curve (i.e. the red solid line) to the
adjacent three peaks. (_b_) Extracted $b_{f}$ in air and in water. The zero
value is indicated by the black dashed line. The inset compares the measured
pitching frequency $f_{p}$ in air and in water.

We conduct ‘ring down’ experiments to measure the fluid damping experienced by
elastically mounted pitching wings. In the ‘ring down’ experiment, a short-
time constant-torque impulse is applied to the CPS as the perturbation, after
which the system response and the amplitude decay of the wing is recorded and
analysed. Figure 2(_a_) shows the results from a typical ‘ring down’
experiment. In this specific case, we use an unswept wing
($\Lambda=0^{\circ}$) which pivots around the mid-chord ($x/c=0.5$) at a
frequency of $f_{p}=0.40$ Hz. We conduct the ‘ring down’ experiment twice –
once in air and once in water. The pitching amplitude of the wing decays
faster in water than in air, indicating a higher total damping in water.

To quantify this amplitude decay, the positive peaks of the system response
are identified. As shown in the inset, the amplitude of the $n$-th peak is
denoted by $A_{n}$, and the corresponding pitching frequency is measured as
$f_{p}=2/(t_{n+1}-t_{n-1})$. To measure the total damping $b+b_{f}$ at
amplitude $A_{n}$, we fit an exponential, $y=\alpha e^{-\gamma t}$, to the
three adjacent peaks, $n-1$, $n$ and $n+1$, and extract the corresponding
$\gamma$ (see equation 3). Now the only unknown in equation $None$_b_ is the
added mass, $I_{a}$. After obtaining $I_{a}$, the fluid damping, $b_{f}$, is
then calculated using equation $None$_a_ (Rao, 1995). Since $f_{p}$ and
$\gamma$ are both measured, $I_{a}$ and $b_{f}$ are also _measured_
quantities. Moreover, both $I_{a}$ and $b_{f}$ are _cycle-averaged_ , meaning
they cannot reflect the instantaneous variation of the fluid inertia and
damping. The measured fluid damping, $b_{f}$, in both air and water are
compared in figure 2(_b_). Since $\tau_{f}$ in equation 1 is negligible in air
as compared to other forces in the equation, $b_{f}$ stays near zero, which is
indicated by the good agreement between the red circles and the black dashed
line. As shown by the green squares, $b_{f}$ in water is significant because
of the existence of the fluid damping torque, $\tau_{b}$. It is also observed
that $b_{f}$ in water increases non-monotonically with $A$. This nonlinear
behaviour will be revisited later in §3.4. The inset of figure 2(_b_) shows
the measured pitching frequency, $f_{p}$, in both air and water. Due to the
combined effect of the fluid inertia and damping, we see that $f_{p}$ is
slightly lower in water than in air.

###  3.2 Frequency scaling of the fluid damping

Figure 3: (_a_) Extracted fluid damping $b_{f}$ at different pitching
frequencies for $x/c=0.5$. (_b_) A frequency scaling for the fluid damping
which collapses $b_{f}$ at different $f_{p}$ into one curve. Note that (_a_)
and (_b_) share the same legend.

We repeat the ‘ring down’ experiment for the unswept wing
($\Lambda=0^{\circ}$) pivoting at the mid-chord, $x/c=0.5$, and change the
pitching frequency by tuning the virtual inertia, $I_{v}$, and the virtual
stiffness, $k_{v}$, while keeping the virtual damping $b_{v}$ constant (Onoue
& Breuer, 2016, 2017). The extracted fluid damping, $b_{f}$, are shown in
figure 3(_a_). Note that figure 3(_a_) and (_b_) share the same legend. We
observe that $b_{f}$ increases monotonically with the pitching frequency,
$f_{p}$, and that the trend of $b_{f}$ remains consistent for all frequencies.
This observation agrees with those observed in heaving rigid plates (Keulegan
& Carpenter, 1958; Shih & Buchanan, 1971), where the fluid damping coefficient
scales inversely with the oscillation period. As we discussed earlier, $b_{f}$
derives from the fluid damping torque $\tau_{b}$, which depends strongly on
the vortex-induced forces on the wing (Kang & Shyy, 2014). Onoue & Breuer
(2016, 2017) have shown that the circulation of LEVs scales with the strength
of the feeding shear-layer velocity. In our case without a free-stream flow,
the feeding shear-layer velocity equals the leading-/trailing-edge velocity,
which is proportional to $f_{p}$. Based on this, we divide $b_{f}$ by $f_{p}$
(figure 3 _b_). It is seen that with this scaling, all of the fluid damping
curves collapse nicely.

We extend this frequency scaling to unswept wings with different pivot axes
(figure 4 _a_) and to swept wings with different sweep angles (figure 4 _b_).
For comparison, we include the previous results (figure 3 _b_) using purple
circles in both figure 4(_a_) and (_b_). Note that each symbol shape in figure
4 contains _five_ different pitching frequencies, $f_{p}=0.20$, 0.28, 0.40,
0.56 and 0.78 Hz.

For the unswept wing ($\Lambda=0^{\circ}$), we change the pivot axis from
$x/c=0$ to 1 with a step size of 0.125 (see the inset of figure 4 _a_). We
observe that $b_{f}/f_{p}$ increases as the pivot axis is moved away from the
mid-chord, $x/c=0.5$. For pivot axes that are symmetric with respect to the
mid-chord (i.e. $x/c=0.375$ & 0.625, 0.25 & 0.75, 0.125 & 0.875 and 0 & 1),
$b_{f}/f_{p}$ roughly overlap. The slight inconsistency between $b_{f}/f_{p}$
for $x/c>0.5$ and $x/c<0.5$ comes from the asymmetry of the NACA 0012 wing
geometry with respect to the mid-chord; we see that the scaled damping,
$b_{f}/f_{p}$, is always slightly higher for $x/c<0.5$. In these cases, the
damping at the trailing edge dominates due to the higher velocity and longer
moment arm, and is stronger than the cases when $x/c>0.5$, where the leading-
edge damping dominates. We will show in §3.4 that this is due to differences
in the vortex structures generated by the sharp and rounded geometries.

This frequency scaling, $b_{f}/f_{p}$, also holds for three-dimensional (3D)
swept wings (figure 4 _b_). Again, each curve includes data from five pitching
frequencies. Here, the pivot axes of swept wings are kept as a vertical line
passing through the mid-chord of the mid-span plane (see the inset of figure 4
_b_). As $\Lambda$ increases, the average pivot axes of the top and the bottom
portion of the swept wing move away from the mid-chord, leading to the
increase of the scaled damping, $b_{f}/f_{p}$, in a manner similar to that
observed for unswept wings with different pivot locations (figure 4 _a_). This
argument will be revisited in the next section.

Figure 4: (_a_) $b_{f}/f_{p}$ for an unswept wing ($\Lambda=0^{\circ}$)
pivoting at $x/c=0$ to 1 with a step size of 0.125. The pivot location for
each dataset is shown by the inset. (_b_) $b_{f}/f_{p}$ for swept wings with
$\Lambda=0^{\circ}$, $10^{\circ}$, $15^{\circ}$, $20^{\circ}$ and
$25^{\circ}$. The inset shows side views of the five swept wings and the
dashed line indicates the pivot axis. The colours of the wings correspond to
the colours of $b_{f}/f_{p}$ curves in the figure. The purple circles in (_a_)
and (_b_) are replotted from figure 3(_b_). Note that each dataset in (_a_)
and (_b_) includes _five_ different $f_{p}$.

###  3.3 Universal fluid damping scaling for unswept and swept wings

Figure 4(_a_) indicates that the pivot axis plays an important role in
determining the fluid damping of unswept wings. We extend the frequency
scaling of $b_{f}$ to take into account this effect. First, we divide the wing
into two parts, the fore part from LE to the pivot axis with a chord length of
$c_{LE}$, and the aft part from the pivot axis to TE with a chord length of
$c_{TE}$ (see the inset of figure 5 for an example when the wing pivots at
$x/c=0.5$). The Morison equation (Morison et al., 1950) indicates that the
fluid damping force $F$ scales with $0.5\rho U^{2}sc$, where $\rho$ is the
fluid density, $U\sim\dot{\theta}c$ is the characteristic velocity and $sc$ is
the wing area. We can express the total fluid damping torque as the sum of the
torque exerted on the fore and aft portions of the wing,

$\tau_{b}\sim K_{LE}F_{LE}c_{LE}+K_{TE}F_{TE}c_{TE},$ (5)

where the subscripts $LE$ and $TE$ refer to the leading- and trailing-edge
contributions, and $K_{LE}$ and $K_{TE}$ are empirical factors that account
for the subtle differences in the damping associated with the specific
geometries of the leading and trailing edges (figure 4 _a_). Since the
differences are small, $K_{LE}$ and $K_{TE}$ should be close to one, and for
consistency, their average value must equal one ($(K_{LE}+K_{TE})/2=1$).

Since the damped oscillations are observed to be near-sinusoidal (figure 2
_a_), the average angular velocity is given by $4f_{p}A$. Simplifying, we
arrive at an expression for the fluid damping:

$b_{f}\sim 2\rho f_{p}As(K_{LE}c_{LE}^{4}+K_{TE}c_{TE}^{4}),$ (6)

or, in non-dimensional form,

$B_{f}^{*}\equiv\frac{b_{f}}{2\rho
f_{p}s(K_{LE}c_{LE}^{4}+K_{TE}c_{TE}^{4})}~{}\propto A.$ (7)

Figure 5: Non-dimensional fluid damping coefficient $B_{f}^{*}$ versus
pitching amplitude $A$ for unswept wings pivoting at $x/c=0$ to 1 and swept
wings with sweep angles $\Lambda=0^{\circ}$ to $25^{\circ}$. The inset shows
the definition of the leading-edge chord $c_{LE}$ and the trailing-edge chord
$c_{TE}$, with black dashed lines indicating the pivot axes. The black dotted
line indicates the small amplitude prediction for a drag coefficient of
$C_{D}=2.8$.

For swept wings, because the pivot axis passes through $x/c=0.5$ at the mid-
span, the top half of the wing has an average pivot axis $x/c>0.5$, while the
bottom half has an average pivot axis $x/c<0.5$. Ignoring three-dimensional
effects, we approximate the swept wing by two ‘equivalent’ unswept wing
segments. We choose not to divide the wing into a large number of narrow
‘blade elements’ (Glauert, 1983), because the pivot axis of some elements near
the wing root/tip for large sweep angles may lie outside the range
$x/c=[0,1]$, where our scaling has not been tested. The inset of figure 5
shows how these two unswept wing segments are configured (rectangles with red
dotted lines). Based on the wing geometry, we see that

$\begin{split}c_{LE,top}&=c_{TE,bot}=\frac{c}{2}+\frac{s}{4}\tan{\Lambda},\\\
c_{TE,top}&=c_{LE,bot}=\frac{c}{2}-\frac{s}{4}\tan{\Lambda}.\end{split}$ (8)

Following the same analysis as for the unswept wing, and adding the fluid
damping of the top and the bottom wing segments together, we find that the
fluid damping for the full swept wing is given by

$b_{f}\sim\rho
f_{p}As(K_{LE}c_{LE,top}^{4}+K_{TE}c_{TE,top}^{4}+K_{LE}c_{LE,bot}^{4}+K_{TE}c_{TE,bot}^{4}).$
(9)

If we define an effective leading-edge chord $c_{LE}=c_{LE,top}=c_{TE,bot}$
and an effective trailing-edge chord $c_{TE}=c_{TE,top}=c_{LE,bot}$, this
scaling reduces to equation 6 with $K_{LE}$ and $K_{TE}$ cancelled out. This
cancellation results because the effective pivot axes of the top and the
bottom segments are symmetric about $x/c=0.5$ at the mid-span, which averages
out the slight differences in fluid damping experienced by the top and the
bottom segments. For the same reason, $K_{LE}$ and $K_{TE}$ also cancel out in
equation 7 for swept wings.

Figure 5 shows the non-dimensional fluid damping, $B_{f}^{*}$, as a function
of the pitching amplitude, $A$, for unswept and swept wings. Here, we have
used $K_{LE}=0.95$ and $K_{TE}=1.05$. We see that all of our measurements
collapse remarkably well under the proposed scaling, especially for
$A<1.57~{}(90^{\circ})$, despite the wide range of pitching frequencies
($f_{p}=0.20$ to 0.78 Hz), pivot axes ($x/c=0$ to 1) and sweep angles
($\Lambda=0^{\circ}$ to $25^{\circ}$) tested in the experiments. In the small-
amplitude limit ($A<0.5$), $B_{f}^{*}$ scales linearly with $A$, with a slope
that corresponds to the drag coefficient, $C_{D}$. We note that $C_{D}\approx
2.8$, which is comparable to that of an accelerated normal flat plate
(Ringuette et al., 2007). At higher pitching angles ($A>0.5$), however, the
linear approximation no longer holds and we see a decreasing slope of
$B_{f}^{*}$ as a function of $A$. This is presumably because the shed vortices
no longer follow the rotating wing and the fluid force becomes non-
perpendicular to the wing surface as $A$ increases. For
$A>1.57~{}(90^{\circ})$, the scaling works reasonably well except for the case
$\Lambda=0^{\circ}$, $x/c=0.5$, where a decreasing $B_{f}^{*}$ is observed. In
the next section, we will use insights from the velocity fields to explain
this non-monotonic behaviour.

###  3.4 Insights obtained from velocity fields

Figure 6: PIV flow field measurements for an unswept wing undergoing
prescribed sinusoidal pitching motions in quiescent water. (_a–d_) Pivot axis
(shown by green dots) $x/c=0.5$, pitching frequency $f_{p}=0.5$ Hz, pitching
amplitude
$A=0.52~{}(30^{\circ}),~{}1.05~{}(60^{\circ}),~{}1.57~{}(90^{\circ})$ and
$2.09~{}(120^{\circ})$. (_e–h_) Same as (_a–d_), except that the pivot axis is
at $x/c=0.25$. All the velocity fields are phase-averaged over 20 cycles. Only
every fifth velocity vector is shown. Spanwise vorticity $\omega$: positive
(red), counterclockwise; negative (blue), clockwise. See supplementary
materials for the full video.

To gain more insight regarding the nonlinear behaviour of $B_{f}^{*}$, we
conduct 2D PIV experiments to measure the surrounding flow fields of an
unswept wing ($\Lambda=0^{\circ}$) with a prescribed pitching motion:
$\theta=A\sin{(2\pi f_{p}t)}$. The results are shown in figure 6. The pitching
frequency is kept at $f_{p}=0.5$ Hz for all the cases and the pitching
amplitude is varied from $A=0.52~{}(30^{\circ})$ to $2.09~{}(120^{\circ})$
with a step size of $0.52~{}(30^{\circ})$. Two pivot axes are tested,
$x/c=0.5$ (figure 6 _a–d_) and $x/c=0.25$ (figure 6 _e–h_). Note that the flow
fields shown in figure 6 are _not_ sequential. Instead, all the snapshots are
taken right before $t/T=0.25$ for different pitching amplitudes, where $T$ is
the pitching period. This specific time instant is chosen because it best
reflects the difference in dynamics associated with the different pitching
amplitudes and pivot axes.

For both pivot locations (figure 6 _a–d_ : $x/c=0.5$ and _e–h_ : $x/c=0.25$),
the spanwise vorticity of the pitch-generated leading-edge vortex (LEV) and
trailing-edge vortex (TEV) increases with the pitching amplitude, $A$. This
can be explained by the increase in the feeding shear-layer velocities
associated with the higher pitching amplitudes (Onoue & Breuer, 2016). The
boundary vortices near the wing surface, which are related to the added mass
effect (Corkery et al., 2019), also become more prominent due to the increase
of the angular acceleration. When the wing pivots at $x/c=0.5$ (figure 6
_a–d_), the leading-edge velocity equals the trailing-edge velocity. As a
result, the LEV and TEV are fairly symmetric about the pivot axis, with some
subtle differences caused by the rounded and sharp edges, respectively. This
confirms the arguments given earlier for the differences between $b_{f}/f_{p}$
for $x/c>0.5$ and $x/c<0.5$ (figure 4 _a_). For $x/c=0.25$ (figure 6 _e–h_),
however, the TEV is much more prominent than the LEV because of the higher
trailing-edge velocity. Due to the low leading-edge velocity and the pitch-
induced rotational flow, the sign of the LEV even reverses and becomes
negative for $A=1.05$ to $2.09$ (figure 6 _f–h_).

For both pivot locations, due to the absence of a convective free stream, and
the existence of the pitch-induced rotational flow, the LEV and TEV (only the
TEV for $x/c=0.25$) are entrained closer to the wing surface as $A$ increases.
For $x/c=0.5$, as shown in figure 6(_c–d_), the LEV moves towards the aft
portion of the wing and the TEV moves towards the fore portion of the wing
when $A\geq 1.57~{}(90^{\circ})$. The torque generated by these two vortices,
which counteracts the wing rotation for small $A$, now assists the rotation as
the wing pitches up towards higher angular positions. This assist reduces the
fluid drag experienced by the wing and thus lowers the fluid damping. This
effect can account for the non-monotonic behaviour of $B_{f}^{*}$ for
$x/c=0.5$ (figure 5). For $x/c=0.25$ (figure 6 _g–h_), a similar scenario is
observed, in which the TEV moves towards the fore portion of the wing and gets
closer to the wing surface as $A$ increases. However, because of the existence
of a counter-rotating LEV, the TEV is not able to approach the wing surface as
closely as in the $x/c=0.5$ case. This explains why a flattening behaviour,
rather than a non-monotonic trend of $B_{f}^{*}$, is observed for $x/c=0.25$
and presumably for other pivot locations at high pitching amplitudes.

##  4 Conclusions

By utilising a cyber-physical control system to create an elastically mounted
pitching wing, we have experimentally measured the nonlinear fluid damping
associated with vortices shed from a bluff body. A theoretical scaling has
been proposed and validated, based on the Morison equation, which incorporates
the frequency, amplitude, pivot location and sweep angle. The nonlinear
behaviour of the scaled fluid damping has been correlated with the velocity
fields measured using particle image velocimetry.

One should note that our scaling may not be applicable for instantaneous fluid
damping, because the damping characterised in the present study is cycle-
averaged over near-sinusoidal oscillations. In addition, we have not
considered three-dimensional effects, which are present due to the wing tip
flows. Incorporating these may further improve the collapse of the fluid
damping coefficient, $B_{f}^{*}$ (figure 5). Lastly, in §3.4, only qualitative
analysis of the flow field has been performed thus far. In order to get more
accurate correspondence between the fluid damping and the flow dynamics,
quantitative analysis of the vortex trajectory and circulation is needed,
which will be the focus of future study.

Despite these limitations, the proposed scaling has been shown to collapse the
data over a wide range of operating conditions ($f_{p}=0.20$ to 0.78 Hz and
$A=0$ to 2.5) for both unswept ($x/c=0$ to 1) and swept wings
($\Lambda=0^{\circ}$ to $25^{\circ}$). It can be used to predict damping
associated with shed vortices, and thus benefit the future modelling of a wide
variety of flows, including unswept and swept wings in unsteady flows as well
as other bluff body geometries. The universality of this scaling reinforces
the underlying connection between swept wings and unswept wings with different
pivot locations. In addition, the results presented in this study will be of
potential value as a source of experimental data for validation and comparison
of future theoretical/computational models.

## Acknowledgments

This work is funded by Air Force Office of Scientific Research, Grant
FA9550-18-1-0322, managed by Dr. Gregg Abate.

## References

  * Alben (2008) Alben, S. 2008 Optimal flexibility of a flapping appendage in an inviscid fluid. J. Fluid Mech. 614, 355.
  * Beatus & Cohen (2015) Beatus, T. & Cohen, I. 2015 Wing-pitch modulation in maneuvering fruit flies is explained by an interplay between aerodynamics and a torsional spring. Phys. Rev. E 92 (2), 022712.
  * Bergou et al. (2007) Bergou, A. J., Xu, S. & Wang, Z. J. 2007 Passive wing pitch reversal in insect flight. J. Fluid Mech. 591, 321–337.
  * Brennen (1982) Brennen, C. E. 1982 A review of added mass and fluid inertial forces. Tech. Rep.. Naval Civil Engineering Laboratory.
  * Corke & Thomas (2015) Corke, T. C. & Thomas, F. O. 2015 Dynamic stall in pitching airfoils: aerodynamic damping and compressibility effects. Annu. Rev. Fluid Mech. 47, 479–505.
  * Corkery et al. (2019) Corkery, S. J., Babinsky, H. & Graham, W. R. 2019 Quantification of added-mass effects using particle image velocimetry data for a translating and rotating flat plate. J. Fluid Mech. 870, 492–518.
  * Dowell et al. (1989) Dowell, E. H., Curtiss, H. C., Scanlan, R. H. & Sisto, F. 1989 A modern course in aeroelasticity. Springer.
  * Dugundji (2008) Dugundji, J. 2008 Some aeroelastic and nonlinear vibration problems encountered on the journey to ithaca. AIAA J. 46 (1), 21–35.
  * Glauert (1983) Glauert, H. 1983 The elements of aerofoil and airscrew theory. Cambridge University Press.
  * Kang & Shyy (2014) Kang, C.-K. & Shyy, W. 2014 Analytical model for instantaneous lift and shape deformation of an insect-scale flapping wing in hover. J. R. Soc. Interface 11 (101), 20140933.
  * Keulegan & Carpenter (1958) Keulegan, G. H. & Carpenter, L. H. 1958 Forces on cylinders and plates in an oscillating fluid. J. Res. Natl. Bur. Stand. 60 (5), 423–440.
  * Mathai et al. (2019) Mathai, V., Loeffen, L. A. W. M., Chan, T. T. K. & Wildeman, S. 2019 Dynamics of heavy and buoyant underwater pendulums. J. Fluid Mech. 862, 348–363.
  * McCroskey (1982) McCroskey, W. J. 1982 Unsteady airfoils. Annu. Rev. Fluid Mech. 14 (1), 285–311.
  * Menon & Mittal (2019) Menon, K. & Mittal, R. 2019 Flow physics and dynamics of flow-induced pitch oscillations of an airfoil. J. Fluid Mech. 877, 582–613.
  * Morison et al. (1950) Morison, J. R., Johnson, J. W. & Schaaf, S. A. 1950 The force exerted by surface waves on piles. J. Pet. Technol. 2 (05), 149–154.
  * Onoue & Breuer (2016) Onoue, K. & Breuer, K. S. 2016 Vortex formation and shedding from a cyber-physical pitching plate. J. Fluid Mech. 793, 229–247.
  * Onoue & Breuer (2017) Onoue, K. & Breuer, K. S. 2017 A scaling for vortex formation on swept and unswept pitching wings. J. Fluid Mech. 832, 697–720.
  * Rao (1995) Rao, S. S. 1995 Mechanical Vibrations. Addison-Wesley.
  * Ringuette et al. (2007) Ringuette, M. J., Milano, M. & Gharib, M. 2007 Role of the tip vortex in the force generation of low-aspect-ratio normal flat plates. J. Fluid Mech. 581, 453–468.
  * Shih & Buchanan (1971) Shih, C. C. & Buchanan, H. J. 1971 The drag on oscillating flat plates in liquids at low reynolds numbers. J. Fluid Mech. 48 (2), 229–239.
  * Shinde & Arakeri (2013) Shinde, S. Y. & Arakeri, J. H. 2013 Jet meandering by a foil pitching in quiescent fluid. Phys. Fluids 25 (4), 041701\.
  * Su & Breuer (2019) Su, Y & Breuer, K. S. 2019 Resonant response and optimal energy harvesting of an elastically mounted pitching and heaving hydrofoil. Phys. Rev. Fluids 4 (6), 064701.
  * Tzezana & Breuer (2019) Tzezana, G. A. & Breuer, K. S. 2019 Thrust, drag and wake structure in flapping compliant membrane wings. J. Fluid Mech. 862, 871–888.
  * Wang (2005) Wang, Z. J. 2005 Dissecting insect flight. Annu. Rev. Fluid Mech. 37, 183–210.
  * Xiao & Zhu (2014) Xiao, Q. & Zhu, Q. 2014 A review on flow energy harvesters based on flapping foils. J. Fluids Struct. 46, 174–191.
  * Young et al. (2014) Young, J., Lai, J. C. S. & Platzer, M. F. 2014 A review of progress and challenges in flapping foil power generation. Prog. Aerosp. Sci. 67, 2–28.
  * Zhu et al. (2020) Zhu, Y., Su, Y. & Breuer, K. S. 2020 Nonlinear flow-induced instability of an elastically mounted pitching wing. J. Fluid Mech. 899, A35.


11institutetext: Department of Computing, Jönköping University, Sweden
11email: {ulf.johansson, tuwe.lofstrom<EMAIL_ADDRESS>
# Well-Calibrated Probabilistic Predictive Maintenance using Venn-
Abers††thanks: The authors acknowledge the Swedish Knowledge Foundation,
Jönköping University, and the industrial partners for financially supporting
the research through the AFAIR project with grant number 20200223, as part of
the research and education environment SPARK at Jönköping University, Sweden.

Ulf Johansson  Tuwe Löfström  Cecilia Sönströd

###### Abstract

When using machine learning for fault detection, a common problem is the fact
that most data sets are very unbalanced, with the minority class (a fault)
being the interesting one. In this paper, we investigate the usage of Venn-
Abers predictors, looking specifically at the effect on the minority class
predictions. A key property of Venn-Abers predictors is that they output well-
calibrated probability intervals. In the experiments, we apply Venn-Abers
calibration to decision trees, random forests and XGBoost models, showing how
both overconfident and underconfident models are corrected. In addition, the
benefit of using the valid probability intervals produced by Venn-Abers for
decision support is demonstrated. When using techniques producing opaque
underlying models, e.g., random forest and XGBoost, each prediction will
consist of not only the label, but also a valid probability interval, where
the width is an indication of the confidence in the estimate. Adding Venn-
Abers on top of a decision tree allows inspection and analysis of the model,
to understand both the underlying relationship, and finding out in which parts
of feature space that the model is accurate and/or confident.

###### Keywords:

Calibration Predictive Maintenance Venn-Abers

##  1 Introduction

Manufacturers and other companies increasingly collect data from sensors in
their factories and products. Predictive maintenance evaluates the condition
of machinery and other equipment by performing offline or online condition
monitoring. The goal of the approach is to decide when a maintenance activity
is most cost-effective but before the equipment loses performance. The
intended result is primarily a reduction in unplanned downtime costs because
of failure, but also to lower the maintenance costs. Recently, predictive
maintenance based on data-driven methods applying machine learning, has become
the most effective solution for smart manufacturing and industrial big data,
especially for performing health perception, e.g., fault diagnosis and
remaining life assessment [16].

When performing fault detection, most data sets are very unbalanced, and it is
the minority class (a fault) that is the interesting one. In fact, it could be
argued that precision must be sacrificed for increased recall, i.e., more
false alarms is acceptable, as long as most faults are in fact detected. With
this in mind, many approaches directly target increasing the number of
detected faults, using for instance cost-sensitive learning or different
sampling techniques, see e.g., [9].

In this paper, we investigate another approach, where we produce well-
calibrated probabilistic predictors. The overall purpose is to provide the
user with a correct estimation, on the instance level, of the probability that
the instance represents a fault.

For the calibration, we suggest the use of Venn-Abers predictors, which in
addition to providing calibrated probabilities, are able to output well-
calibrated probability intervals. In that setting, the width of an interval
represents the confidence in the probability estimation. In the experiments,
we first show the effects of Venn-Abers calibration, using different
underlying models and comparing the results to existing alternatives. After
that, we look specifically at the added informativeness when returning
probability intervals instead of just a single probability estimate,
differentiating between opaque and interpretable predictive models.

##  2 Background

In this section, we first discuss probabilistic prediction in general, before
briefly describing some well-known calibration techniques and providing an
introduction to Venn prediction.

###  2.1 Probabilistic Prediction

A probabilistic predictor outputs not only the predicted class label, but also
a probability distribution over the labels. Ideally, a probabilistic predictor
should be valid, i.e., the predicted probability distributions should perform
well against statistical tests based on subsequent observation of the labels.

In this paper, we focus on calibration, i.e., we want:

$p(c_{j}\mid p^{c_{j}})=p^{c_{j}},$ (1)

where $p^{c_{j}}$ is the probability estimate for class $j$. Specifically, in
order to not be misleading, predicted probabilities must be matched by
observed accuracy. As an example, if a number of predictions with the
probability estimate $0.95$ are made, these predictions should be correct in
about $95\%$ of the cases. While most predictive models are capable of
producing probability estimates, there is no guarantee that these are well-
calibrated. In fact, several models like naive Bayes [10] and decision trees
[12, 8] are notorious for being being poorly calibrated off-the-shelf. Recent
studies show that even models assumed to be generally well-calibrated like
modern (i.e., deep) neural networks [5] and traditional neural networks [6]
often are not. To rectify this, there exist a number of external methods for
calibration, where Platt scaling [11] and isotonic regression [15] are the two
most frequently used. Both these techniques, as well as Venn-Abers [13],
normally utilize a separate labeled data set for the actual calibration.

###  2.2 Platt Scaling and Isotonic Regression

Platt scaling [11] fits a sigmoid function to the scores obtained by the
model. The function is

$\hat{p}(c\mid s)=\frac{1}{1+e^{As+B}},$ (2)

where $\hat{p}(c\mid s)$ is the probability that an instance belongs to class
$c$, given its score $s$. The parameters $A$ and $B$ are found by a gradient
descent search, minimizing the loss function suggested in [11].

Isotonic regression [15], is a non-parametric method, where an isotonic, i.e.,
non-decreasing, step-wise regression function is used for generalized binning.
The function is fitted to the model scores, minimizing a mean square error
criterion, typically using the pair-adjacent violators algorithm.

###  2.3 Venn-Abers Predictors

Venn predictors [14] are multi-probabilistic predictors, i.e., they output $C$
probability distributions – where $C$ is the number of classes – with one of
them, by construction, guaranteed to be valid, as long as the data is IID. The
probabilities produced for a specific class label can be converted into a
probability interval for that label, with the size of the interval indicating
the confidence in the estimates, for details see e.g., [7].

Venn-Abers predictors are Venn predictors, i.e., they inherit the validity
guarantees. Venn-Abers predictors are, however, in the basic form restricted
to two-class problems since they assume that the underlying model is a scoring
classifier. After calibration, which is done by fitting two isotonic
regression [15] models to the calibration set, each prediction is associated
with a well-calibrated probability interval.

More formally, for each test object $x_{l+1}$ being predicted, an inductive
Venn-Abers predictor is defined in the following way.

  1. 1.

Let $\\{z_{1},\dots,z_{l}\\}$ be a training set where each instance
${z=(x,y)}$ consists of an object $x$ and a label $y$.

  2. 2.

Train a scoring classifier using a proper training set
$\\{z_{1},\dots,z_{q}\\}$ using $q<l$ instances.

  3. 3.

Produce prediction scores for the objects in the calibration set, including
the test object, $\\{s(x_{q+1}),\dots,s(x_{l}),s(x_{l+1})\\}$.

  4. 4.

Fit two isotonic calibrators, $g_{0}$ and $g_{1}$, using
$\\{(s(x_{q+1}),y_{q+1}),\dots,(s(x_{l}),y_{l}),(s(x_{l+1}),0)\\}$, and
$\\{(s(x_{q+1}),y_{q+1}),\dots,(s(x_{l}),y_{l}),(s(x_{l+1}),1)\\}$,
respectively.

  5. 5.

Finally, let the probability interval for $y_{l+1}=1$ be
$[g_{0}(s(x_{l+1})),g_{1}(s(x_{l+1}))]$.

It should be noted that when an inductive Venn-Abers predictor is applied on
top of a decision tree, all instances falling in the same leaf will obtain the
same prediction score, so after calibration every leaf will contain a specific
prediction, consisting of a label and an associated confidence (a probability
interval). This very informative model can be inspected and analyzed like any
decision tree.

##  3 Method

The data set used in the experiments was introduced in [9] with the purpose of
providing the research community with a typical predictive maintenance data
set. The data set is presented in detail in the original paper, but will be
briefly described here as well. It consists of six input features representing
quality (low/medium/high), air temperature in Kelvin, process temperature in
Kelvin, rotational speed in rounds per minute (rpm), torque in Newton meter
(Nm), and tool wear in minutes. The class labels are failure (minority class)
and no failure. Out of 10 000 instances, 339 are failures, making the problem
highly imbalanced. Failures could occur for various reasons, represented as
three hidden binary features in the data set, i.e., these should not be used
for the model learning. If one of these hidden binary features is active, then
the instance is labeled as a fault. We have transformed the quality feature
into numeric format and of course removed the three hidden features describing
the specific kinds of failures, but otherwise kept the data as is.

In this paper, we have used decision trees as the transparent models, since
they provide acceptable predictive performance while still being
interpretable. Decision trees are able to provide probability estimates for
its predictions, given as the proportion of samples of the predicted class in
a leaf. The decision tree implementation used is the DecisionTreeClassifier in
scikit-learn.

As one of the opaque models, we used standard Random forests [1], i.e.
ensembles of random trees, trained using bagging and restricted to a random
subset of the available attributes in each split. When using the
RandomForestClassifier in scikit-learn, the probability estimates of the
forest is the mean predicted class probabilities of all trees, where the class
probability of a single tree is the proportion of samples of the predicted
class in a leaf.

Extreme Gradient Boosting (XGBoost) [2], is a highly scalable implementation
of gradient boosting [4]. XGBoost includes several improvements over the
original gradient boosting algorithm, such as weighted quantile sketch for
approximate tree learning as well as enabling parallel and distributed
computation of sparse data. The same calculation to get the probability
estimate as the RandomForestClassifier is used by XGBClassifier in the xgboost
package, implementing the algorithm in Python.

Logistic Regression [3] is a regression method using the logit function to
estimate the log-odds of a binary event based on a linear combination of one
or more attributes. We include logistic regression in the experimentation
since it, by construction, should produce well-calibrated probability
estimates. The scikit-learn implementation of the algorithm was used with the
max-iter parameter set to 500.

When obtaining a single probability estimate from the Venn-Abers probability
interval $(p_{0},p_{1})$, for comparison with the other approaches, the
suggestion in [13] was used:

$p=\frac{p_{1}}{1-p_{0}+p_{1}}$ (3)

thus producing a regularized value with the estimate moved slightly towards
the neutral value $0.5$.

A number of metrics are used in the analysis. Predictive performance is
analysed using accuracy, area under the ROC curve (AUC), precision and recall.

The quality of the calibration was evaluated using the expected calibration
error (ECE). When calculating ECE, the probability estimates for class 1 are
divided into $M$, in this study ten, equally sized bins, before taking a
weighted average of the absolute differences between the fraction of correct
(foc) predictions and the mean of the prediction probabilities (mop). Here,
$n$ is the size of the data set and $B_{m}$ represents bin $m$.

$\textrm{ECE}=\sum_{m=1}^{M}\frac{\left|B_{m}\right|}{n}\Bigl{|}\textrm{foc}(B_{m})-\textrm{mop}(B_{m})\Bigr{|}$
(4)

For the evaluation, standard 10x10-fold cross-validation was used, so all
results reported are averaged over the $100$ folds.

##  4 Results

Table 1 shows the predictive performance for decision trees, random forests
and XGBoost, for Uncalibrated (Uncal) models and after calibration using Venn-
Abers (VA), Platt scaling (Platt) and isotonic regression (IsoReg). Logistic
Regression is also used as comparison. We see that both the random forest and
XGBoost, as expected, are more accurate and have higher AUC than decision
trees and logistic regression. Overall, calibration has a minor impact on
predictive performance, although the AUC decreases slightly for all three
calibrated models. The reason is the smaller training set when setting aside
instances for calibration.

Looking at precision and recall, the random forest has a precision of almost
$0.90$ and a recall of approximately $0.62$, before calibration. Once
calibrated, the precision goes down to between $0.83$ and $0.85$ and the
recall up to between $0.63$ and $0.65$. This is expected, since random forests
are generally underconfident in their probabilistic predictions. For the
decision tree, however, the calibration has the opposite effect; the precision
increases from $0.77$ to between $0.79$ and $0.81$, while the recall decreases
from $0.65$ to between $0.57$ and $0.60$. This is consistent with the fact
that decision trees tend to be overconfident in their probability estimates.
For XGBoost, calibration reduces recall, similarly to the effect seen on
decision trees, but a corresponding increase in precision is not evident.
Logistic regression achieves very low performance, especially considering
recall, and is only able to predict very few instances as positive. In the
targeted scenario, this would not be acceptable, so we exclude logistic
regression from further evaluation.

Table 1: Predictive performance Model | Cal | Acc | AUC | Prec | Rec | #PosPred
---|---|---|---|---|---|---
DT | Uncal | .981 | .904 | .768 | .645 | 2846
| VA | .981 | .902 | .790 | .595 | 2551
| Platt | .981 | .897 | .810 | .566 | 2369
| IsoReg | .981 | .897 | .796 | .588 | 2505
RF | Uncal | .985 | .966 | .895 | .616 | 2334
| VA | .984 | .962 | .831 | .648 | 2644
| Platt | .984 | .957 | .845 | .630 | 2529
| IsoReg | .983 | .959 | .834 | .634 | 2577
XGB | Uncal | .986 | .975 | .858 | .696 | 2752
| VA | .984 | .969 | .842 | .662 | 2664
| Platt | .984 | .932 | .861 | .647 | 2547
| IsoReg | .984 | .967 | .845 | .645 | 2590
LR | Uncal | .970 | .894 | .705 | .196 | 941

Turning to the calibration results, we present results for the overall
calibration, i.e., the probability estimates for the minority class 1 (i.e.,
faults), but we look specifically at the predictions that are actually for the
minority class. As explained above, the most important effect of the
calibration would be to improve the probability estimates for the minority
class predictions, i.e., having $p\geq 0.5$ and consequently predicting a
failure.

We now look at the reliability plots. Figure 1 shows an overall reliability
plot when using decision trees. The decision tree appears to be excellently
calibrated out-of-the box, with an ECE of approximately $0.012$. However, the
three calibration techniques all yield ECE:s between $0.001$ and $0.002$.
Furthermore, when looking at all bins but the leftmost, we can see that the
curves representing the calibrated values are much closer to the diagonal,
especially for Venn-Abers and isotonic regression. We can also see that the
trees are indeed overconfident in their probability estimates, as indicated by
the green line being consistently below the dotted line.

Figure 1: Reliability plot for all instances. Decision trees

Figure 2 shows the corresponding reliability plot when looking only at the
minority class predictions.

Figure 2: Reliability plot for minority class predictions. Decision trees

This is a very interesting plot, since it directly shows how calibration
affects the minority class prediction, i.e., in the investigated scenario the
predictions that would require an action. Here, it is clear that the decision
tree is very overconfident in these predictions. In fact, after calibration,
both Venn-Abers and isotonic regression result in very few predictions with a
higher confidence than 0.9. Apart from Platt scaling, which is still slightly
overconfident, calibrating the decision tree results in very good calibration
without any clear tendency to over- or underconfidence. The ECE drops from
approximately $0.09$ to between $0.02$ and $0.05$.

Turning to the random forests, Figure 3 shows the overall reliability plot.
Again, the probability estimates are well-calibrated, with an ECE below
$0.01$. Consequently, calibration has limited effect.

Figure 3: Reliability plot for all instances. Random forests

However, when looking at the corresponding reliability plot for the minority
class predictions only, in Figure 4, it is seen that the random forest is in
fact very underconfident when predicting the minority class.

Figure 4: Reliability plot for minority class predictions. Random forests

Again, calibration using Venn-Abers and isotonic regression are able to
correct this, leading to an ECE of $0.03$, compared to the original $0.15$.
Platt scaling again improves the calibration compared to the underlying model,
but to a lesser degree than the other calibrators.

Figure 5 shows the overall reliability for XGBoost. Again, the calibration is
excellent, with an ECE below $0.01$, so further calibration has a very minor
effect.

Figure 5: Reliability plot for all instances. XGBoost

When looking at the corresponding reliability plot for the minority class
predictions only, in Figure 6, it is seen that XGBoost is slightly
overconfident when predicting the minority class.

Figure 6: Reliability plot for minority class predictions. XGBoost

Even though XGBoost is fairly well calibrated on minority class predictions,
both Venn-Abers and isotonic regression are able to slightly improve
calibration.

Table 2 summarizes the calibration results. ECE is the expected calibration
error on all instances, and ECE-1 is the expected calibration error on the
minority class predictions only.

Table 2: Expected Calibration Error - ECE | ECE | ECE-1
---|---|---
Cal | DT | RF | XGB | DT | RF | XGB
Uncal | .012 | .007 | .006 | .086 | .147 | .036
VA | .001 | .003 | .005 | .016 | .026 | .031
Platt | .002 | .003 | .002 | .045 | .044 | .039
IsoReg | .001 | .001 | .002 | .022 | .024 | .032

All-in-all, the experimental results show that applying calibration is very
successful when looking at the minority class predictions. Interestingly,
calibration is able to correct both the underconfident random forests and the
overconfident techniques decision trees and XGBoost. Specifically, the
probability estimates for the minority class predictions are much closer to
the actual accuracy, i.e., they would provide a better decision support.

So far we have treated Venn-Abers as a calibration technique, only using the
regularized midpoint of the probability interval. However, for decision-
making, there are two main reasons for why it makes a lot of sense to also use
the probability intervals: i) they provide valid probability estimates of each
class label and, ii) the width of the interval conveys insight into the
uncertainty of the probability estimates.

When using opaque models, Venn-Abers will provide probability intervals per
instance for the positive class, as can be seen in Table 3. The columns are
$y=$ true target, $\hat{y}=$ prediction from the underlying model, $p_{1}=$
the probability estimate for the Failure class of the underlying model,
$\tilde{y}=$ the prediction of the Venn-Abers after calibration, and finally
$p_{1}^{L}$ and $p_{1}^{H}$ represent the probability interval for the
positive class.

Table 3: Venn-Abers predictions from a Random Forest $y$ | $\hat{y}$ | $p_{1}$ | $\tilde{y}$ | $p_{1}^{L}$ | $p_{1}^{H}$
---|---|---|---|---|---
1 | 1 | .92 | 1 | .97 | 1.0
1 | 1 | .78 | 1 | .85 | 1.0
0 | 0 | .47 | 1 (?) | .50 | .70
0 | 0 | .46 | ? | .41 | .70

The examples shown in Table 3 illustrate how the interval of a Venn-Abers
predictor provides both a valid probability interval and a level of certainty.
The third and fourth instances are interesting examples, where the
probabilities are close to $0.5$, thus indicating uncertainty about the label.
In addition, the intervals are wide, showing that the probability estimate in
itself is uncertain.

When Venn-Abers is applied to a decision tree, the leaves will, after
calibration, contain valid probability intervals, i.e., the tree can, in
addition to providing explanations of the predictions, be inspected and
analyzed to increase our understanding of the underlying relationship.
Furthermore, by looking at the widths of the probability intervals in the
leaves, we can identify different parts of feature space where the model is
confident in its estimations and not. Fig. 7 shows a stylized Venn-Abers tree
without the split criteria and the actual probability intervals. This tree is
also pruned to a maximum depth of $5$ for visualization purposes.

Figure 7: A well-calibrated and interpretable Venn tree. Colors represent the
target classes with either no failure (orange) or failure (blue). The colour
intensity corresponds to the estimated probability, while the width of the
leaves gives the sizes of the Venn-Abers probability intervals, indicating how
certain the probability estimates are.

The Venn-Abers tree consequently conveys three things:

  * •

The predicted class in the different leaves is given using different colors;
blue and orange.

  * •

The probability estimates are shown by colour intensity, with higher intensity
representing higher probabilities.

  * •

The confidence in the estimates, represented by the width of the leaves.

The right branch of the tree is more or less strongly dedicated to the
prediction of failures, with some leaves being more certain about their
probability estimates (smaller intervals and stronger color intensity),
whereas other leaves are less certain about the correct class (lower color
intensity), or indicate the model being less confident about its probability
estimate (wider leaves). Most of the central and left parts of the tree are
dedicated to the prediction of no failures with different degrees of
probability and confidence, with some leaves having high probability with high
confidence and others exhibiting a lot of uncertainty. The white and wide
leaves indicate both a probability close to 0.5 and low confidence.

Each leaf also represents a description, in the form of conjunctive conditions
defined for the input attributes of the instances in that part of the feature
space. Examples of rules describing the three numbered leaves in Fig. 7 are
listed in Fig. 8. The interval at the end of each rule is the Venn-Abers
interval for the probability.

$\displaystyle 1)$ $\displaystyle\textrm{torque [Nm] }>67.4$ $\displaystyle\&$
$\displaystyle\textrm{rotational speed [rpm] }>1239$
$\displaystyle\rightarrow$ $\displaystyle\textbf{Failure }[0.91,1.0]$
$\displaystyle 2)$ $\displaystyle\textrm{torque [Nm] }\leq 59.05$
$\displaystyle\&$ $\displaystyle\textrm{torque [Nm] }>13.05$ $\displaystyle\&$
$\displaystyle\textrm{rotational speed [rpm] }>1380.5$ $\displaystyle\&$
$\displaystyle\textrm{tool wear [min] }>201.5$ $\displaystyle\rightarrow$
$\displaystyle\textbf{No Failure }[0.62,0.97]$ $\displaystyle 3)$
$\displaystyle\textrm{torque [Nm] }\leq 59.05$ $\displaystyle\&$
$\displaystyle\textrm{torque [Nm] }>13.05$ $\displaystyle\&$
$\displaystyle\textrm{rotational speed [rpm] }>1380.5$ $\displaystyle\&$
$\displaystyle\textrm{tool wear [min] }\leq 201.5$ $\displaystyle\rightarrow$
$\displaystyle\textbf{No Failure }[0.98,1.0]$ Figure 8: Rules for the three
numbered leaves in Fig. 7

Here, the first rule is very interesting, since it with high probability and
confidence identifies faults. In a similar way, Rule 3, where a vast majority
of all instances end up, identifies no faults. Rule 2, finally, points out a
part of feature space where both the probability estimate and the confidence
are rather low. Here, it could be noted that the only difference between this
rule and Rule 3 is whether the value of the feature tool wear is higher or
lower than $201.5$, suggesting that this could be looked into.

##  5 Concluding Remarks

We have in this paper demonstrated how Venn-Abers calibration could be used in
a fault detection scenario with highly imbalanced classes. The experimental
results show that calibration using either Venn-Abers or the common
alternatives Platt scaling and isotonic regression, will lead to significantly
improved minority class predictions. Specifically, predictions from both
overconfident and underconfident models are corrected.

In the second part, we argued for using Venn-Abers not only for calibration,
but actually utilizing the added information present in the valid probability
intervals produced. When used on top of an opaque model, each predicted label
is complemented with a valid probability interval where the width indicates
the confidence in the estimate. Adding Venn-Abers to a decision tree makes it
straightforward to inspect and analyze the model in order to understand not
only the underlying relationship, but also in which parts of feature space the
model is strong and/or confident.

## References

  * [1] Breiman, L.: Random forests. Machine learning 45(1), 5–32 (2001)
  * [2] Chen, T., Guestrin, C.: Xgboost: A scalable tree boosting system. In: Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining. pp. 785–794 (2016)
  * [3] Cox, D.R.: The regression analysis of binary sequences. Journal of the Royal Statistical Society: Series B (Methodological) 20(2), 215–232 (1958)
  * [4] Friedman, J.H.: Greedy function approximation: a gradient boosting machine. Annals of statistics pp. 1189–1232 (2001)
  * [5] Guo, C., Pleiss, G., Sun, Y., Weinberger, K.Q.: On calibration of modern neural networks. In: Precup, D., Teh, Y.W. (eds.) Proceedings of the 34th International Conference on Machine Learning. Proceedings of Machine Learning Research, vol. 70, pp. 1321–1330. PMLR (2017)
  * [6] Johansson, U., Gabrielsson, P.: Are traditional neural networks well-calibrated? In: 2019 International Joint Conference on Neural Networks (IJCNN). pp. 1–8 (2019)
  * [7] Johansson, U., Löfström, T., Sundell, H., Linusson, H., Gidenstam, A., Boström, H.: Venn predictors for well-calibrated probability estimation trees. In: COPA. pp. 1–12. PMLR (2018)
  * [8] Johansson, U., Löfström, T., Boström, H.: Calibrating probability estimation trees using venn-abers predictors. In: Proceedings of the 2019 SIAM International Conference on Data Mining, SDM 2019, Calgary, Alberta, Canada, May 2-4, 2019. pp. 28–36 (2019)
  * [9] Matzka, S.: Explainable artificial intelligence for predictive maintenance applications. In: 2020 Third International Conference on Artificial Intelligence for Industries (AI4I). pp. 69–74. IEEE (2020)
  * [10] Niculescu-Mizil, A., Caruana, R.: Predicting good probabilities with supervised learning. In: ICML. pp. 625–632. ACM (2005)
  * [11] Platt, J.C.: Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods. In: Advances in Large Margin Classifiers. pp. 61–74. MIT Press (1999)
  * [12] Provost, F., Domingos, P.: Tree induction for probability-based ranking. Mach. Learn. 52(3), 199–215 (2003)
  * [13] Vovk, V., Petej, I.: Venn-abers predictors. arXiv preprint arXiv:1211.0025 (2012)
  * [14] Vovk, V., Shafer, G., Nouretdinov, I.: Self-calibrating probability forecasting. In: Advances in Neural Information Processing Systems. pp. 1133–1140 (2004)
  * [15] Zadrozny, B., Elkan, C.: Obtaining calibrated probability estimates from decision trees and naive Bayesian classifiers. In: ICML. pp. 609–616 (2001)
  * [16] Zhang, W., Yang, D., Wang, H.: Data-driven methods for predictive maintenance of industrial equipment: A survey. IEEE Systems Journal 13(3), 2213–2227 (2019)